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, PlotsThen we define the Hilbert space with N sites and parity conservation.
N = 10
H = hilbert_space(1:N, ParityConservation())1024⨯1024 SymmetricFockHilbertSpace:
modes: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
ParityConservation([-1, 1])Fermions are defined either as symbolic fermions f using the @fermions macro,
@fermions fsym(FermionicHilbertSpaces.SymbolicFermionBasis(:fsym, 0xaab5e7db942062a0),)or as sparse matrix representations fmat using the fermions function, taking H as argument.
fmat = fermions(H)OrderedCollections.OrderedDict{Int64, SparseArrays.SparseMatrixCSC{Int64, Int64}} with 10 entries:
1 => sparse([513, 516, 518, 519, 522, 523, 525, 528, 530, 531 … 494, 495, …
2 => sparse([513, 515, 517, 519, 521, 523, 525, 527, 529, 531 … 493, 495, …
3 => sparse([513, 514, 517, 518, 521, 522, 525, 526, 529, 530 … 493, 494, …
4 => sparse([513, 514, 515, 516, 521, 522, 523, 524, 529, 530 … 491, 492, …
5 => sparse([513, 514, 515, 516, 517, 518, 519, 520, 529, 530 … 487, 488, …
6 => sparse([513, 514, 515, 516, 517, 518, 519, 520, 521, 522 … 487, 488, …
7 => sparse([513, 514, 515, 516, 517, 518, 519, 520, 521, 522 … 471, 472, …
8 => sparse([513, 514, 515, 516, 517, 518, 519, 520, 521, 522 … 439, 440, …
9 => sparse([513, 514, 515, 516, 517, 518, 519, 520, 521, 522 … 375, 376, …
10 => sparse([513, 514, 515, 516, 517, 518, 519, 520, 521, 522 … 247, 248, …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. Here, f can be either symbolic (fsym) or a matrix representation (fmat). 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
params = (μ, t, Δ, U)([-2.0, -4.0, -4.0, -4.0, -4.0, -4.0, -4.0, -4.0, -4.0, -2.0], 1.0, 2.6, 4.0)We can now construct the Hamiltonian using symbolic fermions for a symbolic representation
hsym = kitaev_chain(fsym, N, μ, t, Δ, U)Sum with 55 terms:
- 2.6*fsym†[1]*fsym†[2] - 4.0*fsym†[1]*fsym†[2]*fsym[1]*fsym[2] - 2.0*fsym†[1]*fsym[1] + ...or using matrix representations for the matrix representation.
hmat = kitaev_chain(fmat, N, μ, t, Δ, U)1024×1024 SparseArrays.SparseMatrixCSC{Float64, Int64} with 10238 stored entries:
⎡⢿⣷⣷⡈⠳⠄⠀⠘⢦⡀⠀⠀⠀⠀⠀⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⡙⠻⣿⣿⣴⠦⡀⠀⠀⠙⠂⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠙⠆⠰⡟⠿⣧⣌⣦⠠⣄⣄⠀⠀⠀⠀⠀⠀⠀⠈⠳⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⣀⠀⠀⠈⠢⣽⣿⣿⣦⢈⠉⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠈⠳⣄⠀⠀⢦⡈⢛⣿⣿⡅⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠈⠀⠀⠙⢧⡀⠁⠉⣿⣿⣷⡈⠻⠅⠀⠘⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⡙⠻⣿⣿⣶⠦⡀⠀⠀⠙⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⢤⡀⠀⠀⠀⠀⠀⠀⠀⠙⠟⠆⠸⡟⠿⣧⣬⣦⠠⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⣀⠀⠀⠈⠢⣿⣿⣿⣦⣈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠈⠳⣄⠀⠀⢦⡈⢻⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣿⣿⣧⡈⠳⠀⠀⠙⢦⡀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡉⠻⣿⣿⣿⠢⡀⠀⠀⠉⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠂⠻⡛⢻⣶⣼⡆⠰⣴⣄⠀⠀⠀⠀⠀⠀⠀⠈⠓⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣄⠀⠀⠈⠲⠿⣿⣿⣦⣌⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⡄⠀⢐⣦⡈⢿⣿⣿⣀⢀⠈⢳⣄⠀⠀⡀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⢘⣿⣿⣥⡈⠳⠀⠀⠙⢦⡀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⣀⡁⠻⣿⣿⣟⠢⡀⠀⠀⠉⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢦⡀⠀⠀⠀⠀⠀⠀⠀⠙⠙⠂⠻⡙⢻⣶⣼⠆⠰⣄⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠠⣄⠀⠀⠈⠲⠟⣿⣿⣦⣌⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠀⠀⠀⠀⠈⠳⡄⠀⠐⢦⡈⢿⢿⣷⎦To convert the symbolic Hamiltonian to a matrix representation, we can use the matrix_representation function.
matrix_representation(hsym, H) ≈ hmat # truetrueNow, 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)512×512 SparseArrays.SparseMatrixCSC{Float64, Int64} with 5120 stored entries:
⎡⢿⣷⣷⡈⠳⠄⠀⠘⢦⡀⠀⠀⠀⠀⠀⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⡙⠻⣿⣿⣴⠦⡀⠀⠀⠙⠂⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⎥
⎢⠙⠆⠰⡟⠿⣧⣜⣦⠠⣄⣄⠀⠀⠀⠀⠀⠀⠀⠈⠳⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⎥
⎢⣀⠀⠀⠈⠲⣽⣿⣿⣦⢈⠉⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⎥
⎢⠈⠳⣄⠀⠀⢦⡈⢛⣿⣿⡅⠀⠀⠓⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⎥
⎢⠀⠀⠈⠀⠀⠙⢧⡀⠁⠉⣿⣿⣷⡈⠳⠅⠀⠘⢦⡀⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠙⢤⠀⡙⠻⣿⣿⣶⠦⡀⠀⠀⠙⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⢤⡀⠀⠀⠀⠀⠀⠀⠀⠙⠝⠆⠸⡟⠿⣧⣭⡆⠠⣄⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⣀⠀⠀⠈⠣⠿⣿⢟⣦⣈⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠈⠳⣄⠀⠀⢦⡈⢻⣿⣿⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⣿⣿⣧⡈⠳⠀⠀⠙⢦⡀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⡉⠻⣵⣿⣶⢢⡀⠀⠀⠉⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠙⠂⠸⣛⢻⣶⣼⡆⠰⣔⣄⠀⠀⠀⠀⠀⠀⠀⠈⠓⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⣄⠀⠀⠈⠲⠿⣿⣿⣦⣌⠀⠓⣄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠈⠳⡄⠀⢐⢦⡈⢿⣿⣿⣀⢀⠈⢳⣄⠀⠀⡀⠀⠀⎥
⎢⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢤⠀⠀⢘⣿⣿⣥⡈⠳⠀⠀⠙⢦⡀⎥
⎢⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⣀⡁⠻⣿⣿⣟⠦⡀⠀⠀⠉⎥
⎢⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢦⡀⠀⠀⠀⠀⠀⠀⠀⠙⠙⠂⠻⡝⢻⣶⣼⠆⠰⣄⎥
⎢⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠠⣄⠀⠀⠈⠲⠟⣿⣿⣦⣌⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠀⠀⠀⠀⠈⠳⡄⠀⠐⢦⡈⢿⢿⣷⎦and then diagonalize each sector separately.
oddvec = eigvecs(Matrix(hodd))[:, 1] # odd ground state
evenvec = eigvecs(Matrix(heven))[:, 1] # even ground state512-element Vector{Float64}:
0.03288119401730312
0.033835134895431766
-0.04205167644775886
0.037632282939188876
0.03982007966438107
-0.04508591792174278
0.03516280415418946
0.034213395560383386
-0.03989305953152109
0.042787162289959324
⋮
-0.03989305953152123
0.03421339556038334
0.03516280415418962
-0.04508591792174299
0.039820079664381336
0.03763228293918863
-0.042051676447758896
0.03383513489543187
0.03288119401730323The ground states are almost degenerate, as expected.
first(eigvals(Matrix(hodd))) - first(eigvals(Matrix(heven))) # ≈ 10^-4-0.00014249880423022887Now, we construct the ground state Majoranas. First, we need to embed the lowest energy odd and even states into the full Hilbert space. We can do this by defining a DirectSum type that holds the spaces of the odd and even sectors.
struct DirectSum{HS}
spaces::HS
endThen, we can extend the odd and even vectors to the full Hilbert space.
function extend(v, p=Pair{<:AbstractHilbertSpace,<:DirectSum})
mapreduce(H -> H == first(p) ? v : zeros(size(H, 1)), vcat, last(p).spaces)
end
Hsum = DirectSum((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 space1024-element Vector{Float64}:
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
⋮
-0.03989305953152123
0.03421339556038334
0.03516280415418962
-0.04508591792174299
0.039820079664381336
0.03763228293918863
-0.042051676447758896
0.03383513489543187
0.03288119401730323Then, we can construct the ground state Majorana operators as
γ = o * e' + hc
γ̃ = 1im * o * e' + hc1024×1024 Matrix{ComplexF64}:
-0.0+0.0im 0.0-0.0im … 0.0-0.00122576im 0.0-0.00119121im
0.0+0.0im 0.0+0.0im 0.0+0.00132197im 0.0+0.0012847im
-0.0+0.0im 0.0-0.0im 0.0-0.00120366im 0.0-0.00116973im
-0.0+0.0im 0.0-0.0im 0.0-0.00115405im 0.0-0.00112152im
0.0+0.0im 0.0+0.0im 0.0+0.00122173im 0.0+0.00118728im
0.0+0.0im 0.0+0.0im … 0.0+0.0012806im 0.0+0.0012445im
-0.0+0.0im 0.0-0.0im 0.0-0.00136883im 0.0-0.00133024im
0.0+0.0im 0.0+0.0im 0.0+0.00127393im 0.0+0.00123801im
-0.0+0.0im 0.0-0.0im 0.0-0.00122773im 0.0-0.00119311im
-0.0+0.0im 0.0-0.0im 0.0-0.00126104im 0.0-0.00122549im
⋮ ⋱
0.0-0.00144523im 0.0+0.00155866im … 0.0+0.0im 0.0+0.0im
0.0+0.00123947im 0.0-0.00133675im 0.0+0.0im 0.0+0.0im
0.0+0.00127386im 0.0-0.00137384im 0.0+0.0im 0.0+0.0im
0.0-0.00163335im 0.0+0.00176155im 0.0+0.0im 0.0+0.0im
0.0+0.00144258im 0.0-0.0015558im 0.0+0.0im 0.0+0.0im
0.0+0.00136333im 0.0-0.00147033im … 0.0+0.0im 0.0+0.0im
0.0-0.00152343im 0.0+0.001643im 0.0+0.0im 0.0+0.0im
0.0+0.00122576im 0.0-0.00132197im 0.0+0.0im 0.0+0.0im
0.0+0.00119121im 0.0-0.0012847im 0.0+0.0im 0.0+0.0imFinally, we can check the locality of the Majorana operators. This is done by tracing down the Majorana operators to each mode. So let's first define the Hilbert space for each mode.
Hmodes = [hilbert_space(i:i) for i in 1:N]10-element Vector{SimpleFockHilbertSpace{Int64}}:
2⨯2 SimpleFockHilbertSpace:
modes: [1]
2⨯2 SimpleFockHilbertSpace:
modes: [2]
2⨯2 SimpleFockHilbertSpace:
modes: [3]
2⨯2 SimpleFockHilbertSpace:
modes: [4]
2⨯2 SimpleFockHilbertSpace:
modes: [5]
2⨯2 SimpleFockHilbertSpace:
modes: [6]
2⨯2 SimpleFockHilbertSpace:
modes: [7]
2⨯2 SimpleFockHilbertSpace:
modes: [8]
2⨯2 SimpleFockHilbertSpace:
modes: [9]
2⨯2 SimpleFockHilbertSpace:
modes: [10]Now we can compute the reduction of the Majorana operators to each mode.
γ_reductions = [norm(partial_trace(γ, H => Hmode)) for Hmode in Hmodes]
γ̃_reductions = [norm(partial_trace(γ̃, H => Hmode)) for Hmode in Hmodes]10-element Vector{Float64}:
1.401663147878615
3.7964344997080343e-16
0.17777260969374328
2.082212459856022e-16
0.024489122150747833
6.470498218403839e-17
0.00337629285483671
3.986562883376773e-17
0.0004621281414238482
1.5608926981836748e-16We can plot the reductions to visualize the localization of the Majorana modes.
plot(xlabel="Site", ylabel="||γᵢ|| / √2", title="Majorana Locality", frame=:box, size=(500, 300))
plot!(1:N, γ_reductions / sqrt(2), label="γ", lw=2)
plot!(1:N, γ̃_reductions / sqrt(2), label="γ̃", lw=2)
This page was generated using Literate.jl.