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 ArpackThen we define the Hilbert space with N sites and parity conservation.
N = 12
H = hilbert_space(1:N, ParityConservation())4096-dimensional SymmetricFockHilbertSpace:
12 fermions: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
FockSymmetry(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 potential10-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.0We 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]+f†[4]*f[3]+2.6*f[4]*f[5] + ...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.
import FermionicHilbertSpaces: indices, sector, quantumnumbers
(Eo, o), (Ee, e) = map(quantumnumbers(H)) do parity
Hsec = sector(parity, H)
ham = matrix_representation(hsym, Hsec)
vals, vecs = eigs(ham; nev=1, which=:SR)
inds = indices(Hsec, H)
ground_state = zeros(eltype(vecs), dim(H))
ground_state[inds] = vecs[:, 1]
(; energy=first(vals), ground_state)
end2-element Vector{@NamedTuple{energy::Float64, ground_state::Vector{Float64}}}:
(energy = -30.884311784513162, ground_state = [-0.01719924320326297, 0.018548790166372446, -0.016890092086739035, -0.01619378979334392, 0.01714560172625711, 0.01797227143022799, -0.01920782861443476, 0.017876658689792272, -0.017218005481788652, -0.017685368150297053 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
(energy = -30.884331431867, ground_state = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.020300746748632702, -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. Then, we can construct the ground state Majorana operators as γ = o * e' + hc and γ̃ = 1im * o * e' + hc but that takes a lot of memory for large systems. We can use LowRankMatrices.jl to avoid this
using LowRankMatrices
γ = LowRankMatrix(o, e) + hc
γ̃ = 1im * LowRankMatrix(o, e) + hc
δρ = LowRankMatrix(o, o) - LowRankMatrix(e, e)4096×4096 LowRankMatrices.LowRankMatrix{Float64}:
0.000295814 -0.000319025 0.000290497 … 0.0 0.0
-0.000319025 0.000344058 -0.000313291 0.0 0.0
0.000290497 -0.000313291 0.000285275 0.0 0.0
0.000278521 -0.000300375 0.000273515 0.0 0.0
-0.000294891 0.00031803 -0.000289591 0.0 0.0
-0.000309109 0.000333364 -0.000303553 … 0.0 0.0
0.00033036 -0.000356282 0.000324422 0.0 0.0
-0.000307465 0.00033159 -0.000301938 0.0 0.0
0.000296137 -0.000319373 0.000290814 0.0 0.0
0.000304175 -0.000328042 0.000298707 0.0 0.0
⋮ ⋱ ⋮
0.0 0.0 0.0 0.000303988 0.000295419
0.0 0.0 0.0 -0.000260864 -0.000253511
0.0 0.0 0.0 -0.0002681 -0.000260544
0.0 0.0 0.0 … 0.00034381 0.000334119
0.0 0.0 0.0 -0.000303638 -0.00029508
0.0 0.0 0.0 -0.000286919 -0.000278832
0.0 0.0 0.0 0.000320621 0.000311584
0.0 0.0 0.0 -0.000257955 -0.000250684
0.0 0.0 0.0 … -0.000250684 -0.000243618Now we can compute the reduction of the Majorana operators to each mode.
Hmodes = [hilbert_space(i:i) for i in 1:N]
γR = [partial_trace(γ, H => Hmode) for Hmode in Hmodes];
γ̃R = [partial_trace(γ̃, H => Hmode) for Hmode in Hmodes];
δρR = [partial_trace(δρ, H => Hmode) for Hmode in Hmodes];
γ_reductions = [norm(svdvals(γ), 1) for γ in γR]
γ̃_reductions = [norm(svdvals(γ̃), 1) for γ̃ in γ̃R]
LD = [norm(svdvals(δρ), 1) for δρ in δρR]12-element Vector{Float64}:
6.532101248790667e-15
3.10775710721245e-15
1.6837659738699884e-15
8.996709627284716e-16
1.1753835751915354e-15
1.0685896612017132e-15
3.649424512586208e-16
3.3837949803272593e-16
3.5843723822370777e-16
1.6894580352755373e-15
3.0939335295132597e-15
3.6660670159338604e-15We 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.