Interacting Kitaev Chain Tutorial

In this example, we study the interacting Kitaev chain. We show how the Hamiltonian can be constructed symbolically or by using matrix representations of the fermionic operators, and how to restrict the Hilbert space to a subspace of a given parity. We solve for the ground states and characterize the locality of the many-body Majoranas.

We start by importing the necessary packages.

using FermionicHilbertSpaces, LinearAlgebra, Plots
using Arpack

Then we define the Hilbert space with N sites and parity conservation.

N = 12
H = hilbert_space(1:N, ParityConservation())

4096⨯4096 SymmetricFockHilbertSpace:
modes: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
ParityConservation([-1, 1])

Symbolic fermions can be defined using the @fermions macro,

@fermions f

(FermionicHilbertSpaces.SymbolicFermionBasis(:f, 0x39415815c723b896),)

Let's define the interacting Kitaev chain Hamiltonian. It is a function of the fermions f and parameters N, μ, t, Δ, and U, representing the number of sites, chemical potential, hopping amplitude, pairing amplitude, and interaction strength, respectively. Note the use of the Hermitian conjugate hc, which simplifies the expression for the Hamiltonian.

kitaev_chain(f, N, μ, t, Δ, U) = sum(t * f[i]' * f[i+1] + hc for i in 1:N-1) +
                                 sum(Δ * f[i] * f[i+1] + hc for i in 1:N-1) +
                                 sum(U * f[i]' * f[i] * f[i+1]' * f[i+1] for i in 1:N-1) +
                                 sum(μ[i] * f[i]' * f[i] for i in 1:N)

kitaev_chain (generic function with 1 method)

Define parameters close to the sweet spot with perfectly localized Majoranas.

U = 4.0
t = 1.0
δΔ = 0.4
Δ = t + U / 2 - δΔ # slightly detuned from the sweet spot
μ = fill(-U / 2, N) # edge chemical potential
μ[2:N-1] .= -U # bulk chemical potential

10-element view(::Vector{Float64}, 2:11) with eltype Float64:
 -4.0
 -4.0
 -4.0
 -4.0
 -4.0
 -4.0
 -4.0
 -4.0
 -4.0
 -4.0

We can now construct the Hamiltonian using symbolic fermions for a symbolic representation

hsym = kitaev_chain(f, N, μ, t, Δ, U)

Sum with 67 terms: 
 - 4.0*f†[8]*f†[9]*f[8]*f[9] + 2.6*f[4]*f[5] + f†[4]*f[3] + ...

To convert the symbolic Hamiltonian to a matrix representation, we can use the matrix_representation function.

matrix_representation(hsym, H)

4096×4096 SparseArrays.SparseMatrixCSC{Float64, Int64} with 49151 stored entries:
⎡⢿⣷⣷⡈⠳⠄⠀⠘⢦⡀⠀⠀⠀⠀⠀⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⡙⠻⣿⣿⣴⢦⡀⠀⠀⠙⠂⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠙⠆⠰⣟⠿⣧⣝⣦⠠⣄⣄⠀⠀⠀⠀⠀⠀⠀⠈⠳⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⣀⠀⠀⠈⠳⣽⣿⣿⣦⢈⠉⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠈⠳⣄⠀⠀⢦⡈⢛⣿⣿⡅⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠈⠀⠀⠙⢧⡀⠁⠉⣿⣿⣷⡈⠻⠅⠀⠘⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⡙⠻⣿⣿⣶⢦⡀⠀⠀⠙⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⢤⡀⠀⠀⠀⠀⠀⠀⠀⠙⠟⠆⠸⣟⠿⣧⣽⣦⠠⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⣀⠀⠀⠈⠳⣿⣿⣿⣦⣈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠈⠳⣄⠀⠀⢦⡈⢻⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣿⣿⣧⡈⠳⠀⠀⠙⢦⡀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡉⠻⣿⣿⣿⢦⡀⠀⠀⠉⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠂⠻⣟⢻⣶⣽⡆⠰⣴⣄⠀⠀⠀⠀⠀⠀⠀⠈⠓⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣄⠀⠀⠈⠳⠿⣿⣿⣦⣌⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⡄⠀⢐⣦⡈⢿⣿⣿⣀⢀⠈⢳⣄⠀⠀⡀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⢘⣿⣿⣥⡈⠳⠀⠀⠙⢦⡀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⣀⡁⠻⣿⣿⣟⢦⡀⠀⠀⠉⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢦⡀⠀⠀⠀⠀⠀⠀⠀⠙⠙⠂⠻⣝⢻⣶⣽⠆⠰⣄⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠠⣄⠀⠀⠈⠳⠟⣿⣿⣦⣌⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠀⠀⠀⠀⠈⠳⡄⠀⠐⢦⡈⢿⢿⣷⎦

Now, let's diagonalize the system. Since parity is conserved, we can work in the even and odd parity sectors separately. To do this, we can create two new Hilbert spaces for the even and odd sectors

Heven = hilbert_space(1:N, ParityConservation(1))
Hodd = hilbert_space(1:N, ParityConservation(-1))
heven = matrix_representation(hsym, Heven)
hodd = matrix_representation(hsym, Hodd)

2048×2048 SparseArrays.SparseMatrixCSC{Float64, Int64} with 24576 stored entries:
⎡⢿⣷⣷⡈⠳⠄⠀⠘⢦⡀⠀⠀⠀⠀⠀⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⡙⠻⣿⣿⣴⢦⡀⠀⠀⠙⠂⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⎥
⎢⠙⠆⠰⣟⠿⣧⣝⣦⠠⣄⣄⠀⠀⠀⠀⠀⠀⠀⠈⠳⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⎥
⎢⣀⠀⠀⠈⠳⣽⣿⣿⣦⢈⠉⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⎥
⎢⠈⠳⣄⠀⠀⢦⡈⢛⣿⣿⡅⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⎥
⎢⠀⠀⠈⠀⠀⠙⢧⡀⠁⠉⣿⣿⣷⡈⠻⠅⠀⠘⢦⡀⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⡙⠻⣿⣿⣶⢦⡀⠀⠀⠙⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⢤⡀⠀⠀⠀⠀⠀⠀⠀⠙⠟⠆⠸⣟⠿⣧⣽⣦⠠⣄⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⣀⠀⠀⠈⠳⣿⣿⣿⣦⣈⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠈⠳⣄⠀⠀⢦⡈⢻⣿⣿⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⣿⣿⣧⡈⠳⠀⠀⠙⢦⡀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⡉⠻⣿⣿⣿⢦⡀⠀⠀⠉⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠙⠂⠻⣟⢻⣶⣽⡆⠰⣴⣄⠀⠀⠀⠀⠀⠀⠀⠈⠓⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⣄⠀⠀⠈⠳⠿⣿⣿⣦⣌⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠈⠳⡄⠀⢐⣦⡈⢿⣿⣿⣀⢀⠈⢳⣄⠀⠀⡀⠀⠀⎥
⎢⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⢘⣿⣿⣥⡈⠳⠀⠀⠙⢦⡀⎥
⎢⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⣀⡁⠻⣿⣿⣟⢦⡀⠀⠀⠉⎥
⎢⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢦⡀⠀⠀⠀⠀⠀⠀⠀⠙⠙⠂⠻⣝⢻⣶⣽⠆⠰⣄⎥
⎢⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠠⣄⠀⠀⠈⠳⠟⣿⣿⣦⣌⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠀⠀⠀⠀⠈⠳⡄⠀⠐⢦⡈⢿⢿⣷⎦

and then diagonalize each sector separately.

oddeigs = eigs(hodd; nev=1, which=:SR)
eveneigs = eigs(heven; nev=1, which=:SR)
oddvec = oddeigs[2][:, 1] # odd ground state
evenvec = eveneigs[2][:, 1] # even ground state

2048-element Vector{Float64}:
  0.015608274175071248
  0.016060979381801645
 -0.01996274836346006
  0.01786437579224462
  0.018905344284311718
 -0.021406546465997184
  0.016692660286701376
  0.01624211725860512
 -0.01892709491454756
  0.02030074674863281
  ⋮
 -0.018927094914547606
  0.016242117258605148
  0.016692660286701418
 -0.021406546465997045
  0.018905344284311825
  0.017864375792244647
 -0.01996274836346016
  0.016060979381801885
  0.015608274175071112

The ground states are almost degenerate, as expected.

first(oddeigs[1]) - first(eveneigs[1])

1.964735383808147e-5

Now, we construct the ground state Majoranas. First, we need to embed the lowest energy odd and even states into the full Hilbert space.

function extend(v, p)
    mapreduce(H -> H == first(p) ? v : zeros(size(H, 1)), vcat, last(p))
end

Hsum = (Hodd, Heven)
o = extend(oddvec, Hodd => Hsum) # odd ground state in the full Hilbert space
e = extend(evenvec, Heven => Hsum) # even ground state in the full Hilbert space

4096-element Vector{Float64}:
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  ⋮
 -0.018927094914547606
  0.016242117258605148
  0.016692660286701418
 -0.021406546465997045
  0.018905344284311825
  0.017864375792244647
 -0.01996274836346016
  0.016060979381801885
  0.015608274175071112

Then, we can construct the ground state Majorana operators as γ = o * e' + hc γ̃ = 1im * o * e' + hc but that takes a lot of memory for large systems. Let's make a quick struct representing a rank-1 matrix and use that instead.

struct Rank1Matrix{T} <: AbstractMatrix{T}
    vec1::Vector{T}
    vec2::Vector{T}
end
Base.getindex(m::Rank1Matrix, i::Int, j::Int) = m.vec1[i] * conj(m.vec2[j])
Base.size(m::Rank1Matrix) = (length(m.vec1), length(m.vec2))

oe = Rank1Matrix(o, e)
ee = Rank1Matrix(e, e)
oo = Rank1Matrix(o, o)

4096×4096 Main.Rank1Matrix{Float64}:
  0.000295814  -0.000319025  …  -0.0  -0.0  -0.0  -0.0  -0.0  -0.0
 -0.000319025   0.000344058      0.0   0.0   0.0   0.0   0.0   0.0
  0.000290497  -0.000313291     -0.0  -0.0  -0.0  -0.0  -0.0  -0.0
  0.000278521  -0.000300375     -0.0  -0.0  -0.0  -0.0  -0.0  -0.0
 -0.000294891   0.00031803       0.0   0.0   0.0   0.0   0.0   0.0
 -0.000309109   0.000333364  …   0.0   0.0   0.0   0.0   0.0   0.0
  0.00033036   -0.000356282     -0.0  -0.0  -0.0  -0.0  -0.0  -0.0
 -0.000307465   0.00033159       0.0   0.0   0.0   0.0   0.0   0.0
  0.000296137  -0.000319373     -0.0  -0.0  -0.0  -0.0  -0.0  -0.0
  0.000304175  -0.000328042     -0.0  -0.0  -0.0  -0.0  -0.0  -0.0
  ⋮                          ⋱   ⋮                             ⋮
 -0.0           0.0              0.0   0.0   0.0   0.0   0.0   0.0
 -0.0           0.0              0.0   0.0   0.0   0.0   0.0   0.0
 -0.0           0.0              0.0   0.0   0.0   0.0   0.0   0.0
 -0.0           0.0          …   0.0   0.0   0.0   0.0   0.0   0.0
 -0.0           0.0              0.0   0.0   0.0   0.0   0.0   0.0
 -0.0           0.0              0.0   0.0   0.0   0.0   0.0   0.0
 -0.0           0.0              0.0   0.0   0.0   0.0   0.0   0.0
 -0.0           0.0              0.0   0.0   0.0   0.0   0.0   0.0
 -0.0           0.0          …   0.0   0.0   0.0   0.0   0.0   0.0

Now we can compute the reduction of the Majorana operators to each mode.

Hmodes = [hilbert_space(i:i) for i in 1:N]
eoR = [partial_trace(oe, H => Hmode) for Hmode in Hmodes];
eeR = [partial_trace(ee, H => Hmode) for Hmode in Hmodes];
ooR = [partial_trace(oo, H => Hmode) for Hmode in Hmodes];
γ_reductions = [norm(svdvals(eoR + hc), 1) for eoR in eoR]
γ̃_reductions = [norm(svdvals(1im * eoR + hc), 1) for eoR in eoR]
LD = [norm(svdvals(ooR - eeR), 1) for (ooR, eeR) in zip(ooR, eeR)]

12-element Vector{Float64}:
 7.66053886991358e-15
 2.7755575615628914e-15
 8.881784197001252e-16
 5.551115123125783e-16
 1.942890293094024e-15
 4.996003610813204e-16
 1.1102230246251565e-15
 6.661338147750939e-16
 7.771561172376096e-16
 1.7208456881689926e-15
 2.275957200481571e-15
 4.9960036108132044e-15

We can plot the reductions to visualize the localization of the Majorana modes.

#
lw = 4
legendfontsize = 15
marker = true
markerstrokewidth = 2
plot(xlabel="Site", title="Majorana quality measures"; frame=:box, size=(500, 300), xticks=1:3:N, yscale=:identity, legendfontsize, ylims=(-1e-1, 2), legendposition=:top, labelfontsize=15)
plot!(1:N, γ_reductions; label="‖γₙ‖", lw, legendfontsize, marker, markerstrokewidth)
plot!(1:N, γ̃_reductions; label="‖γ̃ₙ‖", lw, legendfontsize, marker, markerstrokewidth)
plot!(1:N, LD; label="‖(iγγ̃)ₙ‖", lw, legendfontsize, marker, markerstrokewidth)

This page was generated using Literate.jl.