using Markdown

Hubbard chain at half filling

The Hubbard model is a model of interacting fermions on a lattice, which is often used as a somewhat realistic model for electrons in a solid. The Hamiltonian consists of two terms that describe competing forces of each electron: a kinetic term that allows electrons to hop between neighboring sites, and a potential term reflecting on-site interactions between electrons. Often, a third term is included which serves as a chemical potential to control the number of electrons in the system.

\[H = -t \sum_{\langle i, j \rangle, \sigma} c^{\dagger}_{i,\sigma} c_{j,\sigma} + U \sum_i n_{i,\uparrow} n_{i,\downarrow} - \mu \sum_{i,\sigma} n_{i,\sigma}\]

At half-filling, the system exhibits particle-hole symmetry, which can be made explicit by rewriting the Hamiltonian slightly. First, we fix the overall energy scale by setting t = 1, and then shift the total energy by adding a constant U / 4, as well as shifting the chemical potential to N U / 2. This results in the following Hamiltonian:

\[H = - \sum_{\langle i, j \rangle, \sigma} c^{\dagger}_{i,\sigma} c_{j,\sigma} + U / 4 \sum_i (1 - 2 n_{i,\uparrow}) (1 - 2 n_{i,\downarrow}) - \mu \sum_{i,\sigma} n_{i,\sigma}\]

Finally, setting \mu = 0 and defining u = U / 4 we obtain the Hubbard model at half-filling.

\[H = - \sum_{\langle i, j \rangle, \sigma} c^{\dagger}_{i,\sigma} c_{j,\sigma} + u \sum_i (1 - 2 n_{i,\uparrow}) (1 - 2 n_{i,\downarrow})\]

using TensorKit
using MPSKit
using MPSKitModels
using SpecialFunctions: besselj0, besselj1
using QuadGK: quadgk
using Plots
using Interpolations
using Optim

const t = 1.0
const mu = 0.0
const U = 3.0

3.0

For this case, the groundstate energy has an analytic solution, which can be used to benchmark the numerical results. It follows from Eq. (6.82) in .

\[e(u) = - u - 4 \int_0^{\infty} \frac{d\omega}{\omega} \frac{J_0(\omega) J_1(\omega)}{1 + \exp(2u \omega)}\]

We can easily verify this by comparing the numerical results to the analytic solution.

function hubbard_energy(u; rtol=1e-12)
    integrandum(ω) = besselj0(ω) * besselj1(ω) / (1 + exp(2u * ω)) / ω
    int, err = quadgk(integrandum, 0, Inf; rtol=rtol)
    return -u - 4 * int
end

function compute_groundstate(psi, H;
                             svalue=1e-3,
                             expansionfactor=(1 / 5),
                             expansioniter=20)
    verbosity = 1
    psi, = find_groundstate(psi, H; tol=svalue * 10, verbosity)
    for _ in 1:expansioniter
        D = maximum(x -> dim(left_virtualspace(psi, x)), 1:length(psi))
        D′ = max(2, round(Int, D * expansionfactor))
        trscheme = truncbelow(svalue / 10) & truncdim(D′)
        psi′, = changebonds(psi, H, OptimalExpand(; trscheme=trscheme))
        all(left_virtualspace.(Ref(psi), 1:length(psi)) .==
            left_virtualspace.(Ref(psi′), 1:length(psi))) && break
        psi, = find_groundstate(psi′, H, VUMPS(; tol=svalue / 5, verbosity))
    end

    # convergence steps
    psi, = changebonds(psi, H, SvdCut(; trscheme=truncbelow(svalue)))
    psi, = find_groundstate(psi, H,
                            VUMPS(; tol=svalue, verbosity) &
                            GradientGrassmann(; tol=svalue / 100, verbosity))

    return psi
end

H = hubbard_model(InfiniteChain(2); U, t, mu=U / 2)
psi = InfiniteMPS(H.data.pspaces, H.data.pspaces)
psi = compute_groundstate(psi, H)
E = real(expectation_value(psi, H)) / 2
@info """
Groundstate energy:
    * numerical: $E
    * analytic: $(hubbard_energy(U / 4) - U / 4)
"""

┌ Warning: ignoring imaginary component 7.807755346263268e-6 from total weight 1.519432279920462: operator might not be hermitian?
│   α = 0.408644890723205 + 7.807755346263268e-6im
│   β₁ = 0.9383819653579948
│   β₂ = 1.1229973702462615
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 8.405210114706407e-6 from total weight 2.1738960127263907: operator might not be hermitian?
│   α = 0.27297674897930924 + 8.405210114706407e-6im
│   β₁ = 1.1229973702462615
│   β₂ = 1.8412453598077907
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.0370915492496913e-5 from total weight 2.827735337332976: operator might not be hermitian?
│   α = 1.1648009191635844 + 1.0370915492496913e-5im
│   β₁ = 1.8412453598077907
│   β₂ = 1.8025375118415796
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 9.968737032850505e-6 from total weight 2.317491673659234: operator might not be hermitian?
│   α = 0.7920195374082283 + 9.968737032850505e-6im
│   β₁ = 1.8025375118415796
│   β₂ = 1.2224284143245874
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 7.730765716737475e-6 from total weight 1.931808439300787: operator might not be hermitian?
│   α = 0.4590455194221491 + 7.730765716737475e-6im
│   β₁ = 1.2224284143245874
│   β₂ = 1.4236677382887843
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 9.461798323513737e-6 from total weight 1.885442422423453: operator might not be hermitian?
│   α = 0.7494644422442023 + 9.461798323513737e-6im
│   β₁ = 1.4236677382887843
│   β₂ = 0.9830393425248946
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.4794589075523001e-5 from total weight 1.7036029158608066: operator might not be hermitian?
│   α = 0.8543239115625937 + 1.4794589075523001e-5im
│   β₁ = 0.9830393425248946
│   β₂ = 1.0981926970669473
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.4779092975814034e-5 from total weight 1.801371752899485: operator might not be hermitian?
│   α = 0.49812454899961545 + 1.4779092975814034e-5im
│   β₁ = 1.0981926970669473
│   β₂ = 1.3382021243887245
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.5082805806504451e-5 from total weight 2.3893902682891692: operator might not be hermitian?
│   α = 1.2559157134638745 + 1.5082805806504451e-5im
│   β₁ = 1.3382021243887245
│   β₂ = 1.5300577273173481
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.4763764194347662e-5 from total weight 2.1971628324934835: operator might not be hermitian?
│   α = 0.5374925927117064 + 1.4763764194347662e-5im
│   β₁ = 1.5300577273173481
│   β₂ = 1.482413429556458
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.6812044682155275e-5 from total weight 2.518373809924349: operator might not be hermitian?
│   α = 1.6677501130681405 + 1.6812044682155275e-5im
│   β₁ = 1.482413429556458
│   β₂ = 1.1675900952226779
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.1656237645502748e-5 from total weight 1.836376541200128: operator might not be hermitian?
│   α = 0.8424638093750539 + 1.1656237645502748e-5im
│   β₁ = 1.1675900952226779
│   β₂ = 1.1398538943063923
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.017848068659774e-5 from total weight 1.5505038790334096: operator might not be hermitian?
│   α = 0.5675363971143528 + 1.017848068659774e-5im
│   β₁ = 1.1398538943063923
│   β₂ = 0.8847021060103001
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 8.048645567207413e-6 from total weight 1.50620841149816: operator might not be hermitian?
│   α = 0.7404925865623019 + 8.048645567207413e-6im
│   β₁ = 0.8847021060103001
│   β₂ = 0.9683164212540453
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.2039057236130912e-5 from total weight 1.6744528808996058: operator might not be hermitian?
│   α = 0.9406358537567846 + 1.2039057236130912e-5im
│   β₁ = 0.9683164212540453
│   β₂ = 0.9906361335854118
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.416432382847499e-5 from total weight 1.8325177435853954: operator might not be hermitian?
│   α = 1.5919932470890588 + 1.416432382847499e-5im
│   β₁ = 0.646450554019812
│   β₂ = 0.6370089974129052
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.4532453500901155e-5 from total weight 1.6446897126802875: operator might not be hermitian?
│   α = 1.2841210428311594 + 1.4532453500901155e-5im
│   β₁ = 0.6370089974129052
│   β₂ = 0.8063851036312584
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.3416927348061658e-5 from total weight 1.521631679409698: operator might not be hermitian?
│   α = 1.1434117013792926 + 1.3416927348061658e-5im
│   β₁ = 0.8063851036312584
│   β₂ = 0.5980933985538383
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.2158232090938442e-5 from total weight 1.5502040087386693: operator might not be hermitian?
│   α = 1.2769345460168038 + 1.2158232090938442e-5im
│   β₁ = 0.5980933985538383
│   β₂ = 0.6440923228519313
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.0742972417238095e-5 from total weight 1.1117183071326222: operator might not be hermitian?
│   α = 0.7297689321905663 + 1.0742972417238095e-5im
│   β₁ = 0.6440923228519313
│   β₂ = 0.5371219410441697
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.2433513098084048e-5 from total weight 1.9525016713436172: operator might not be hermitian?
│   α = 1.8430983693388352 + 1.2433513098084048e-5im
│   β₁ = 0.5371219410441697
│   β₂ = 0.35602134463329077
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 6.78403706365714e-6 from total weight 1.66019539364691: operator might not be hermitian?
│   α = 1.6215725537902703 + 6.78403706365714e-6im
│   β₁ = 0.35602134463329077
│   β₂ = 2.4874957371557414e-21
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 4.420632150038684e-6 from total weight 1.5700523001952702: operator might not be hermitian?
│   α = 0.5239626218600479 + 4.420632150038684e-6im
│   β₁ = 0.8717165441447274
│   β₂ = 1.1960926648411632
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 3.601080797851694e-6 from total weight 2.0345126979997414: operator might not be hermitian?
│   α = 0.426288757350352 + 3.601080797851694e-6im
│   β₁ = 1.1960926648411632
│   β₂ = 1.589616982414027
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 4.19976496114291e-6 from total weight 2.0484867724942677: operator might not be hermitian?
│   α = 0.5775202889239899 + 4.19976496114291e-6im
│   β₁ = 1.589616982414027
│   β₂ = 1.1558054430432567
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 3.640311513675165e-6 from total weight 1.6360531190605416: operator might not be hermitian?
│   α = 0.07558256739995857 + 3.640311513675165e-6im
│   β₁ = 1.1558054430432567
│   β₂ = 1.155452665283735
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 2.5033396548745923e-6 from total weight 1.871835891862734: operator might not be hermitian?
│   α = 0.0559770634887897 + 2.5033396548745923e-6im
│   β₁ = 1.155452665283735
│   β₂ = 1.4715859854969997
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 2.6430412059241792e-6 from total weight 1.9064423191724782: operator might not be hermitian?
│   α = 0.22488899189261616 + 2.6430412059241792e-6im
│   β₁ = 1.4715859854969997
│   β₂ = 1.1909584144457346
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 3.0997451073364046e-6 from total weight 1.6098273943138954: operator might not be hermitian?
│   α = -0.18729098151608492 + 3.0997451073364046e-6im
│   β₁ = 1.1909584144457346
│   β₂ = 1.0668103780792073
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 3.8430909460249185e-6 from total weight 1.4501419373064255: operator might not be hermitian?
│   α = -0.1756932876023117 + 3.8430909460249185e-6im
│   β₁ = 1.0668103780792073
│   β₂ = 0.9664156063692543
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 4.737672034002549e-6 from total weight 1.6650488951358715: operator might not be hermitian?
│   α = 0.3642826605632644 + 4.737672034002549e-6im
│   β₁ = 0.9664156063692543
│   β₂ = 1.3060347783078508
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 4.0463279378337436e-6 from total weight 1.809699994321244: operator might not be hermitian?
│   α = 0.4792707902891812 + 4.0463279378337436e-6im
│   β₁ = 1.3060347783078508
│   β₂ = 1.1574051740232911
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 5.543573193336759e-6 from total weight 1.4866339796856036: operator might not be hermitian?
│   α = -0.07750469979704211 + 5.543573193336759e-6im
│   β₁ = 1.1574051740232911
│   β₂ = 0.9297778628138067
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 4.437985246134174e-6 from total weight 1.2526483933898018: operator might not be hermitian?
│   α = 0.14550455114158256 + 4.437985246134174e-6im
│   β₁ = 0.9297778628138067
│   β₂ = 0.826722171506748
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 3.479709569688405e-6 from total weight 1.2697888613350266: operator might not be hermitian?
│   α = 0.5191679357429178 + 3.479709569688405e-6im
│   β₁ = 0.826722171506748
│   β₂ = 0.8120091489595418
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 4.309653650808043e-6 from total weight 1.1619455195047197: operator might not be hermitian?
│   α = -0.03585372914963633 + 4.309653650808043e-6im
│   β₁ = 0.8120091489595418
│   β₂ = 0.8303451345016631
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 3.3876945127499225e-6 from total weight 1.5863366135243355: operator might not be hermitian?
│   α = 0.10179881851551338 + 3.3876945127499225e-6im
│   β₁ = 0.8303451345016631
│   β₂ = 1.3478233599231915
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 2.4515236037372934e-6 from total weight 1.6563602496932375: operator might not be hermitian?
│   α = 0.2942178379823914 + 2.4515236037372934e-6im
│   β₁ = 1.3478233599231915
│   β₂ = 0.9166991496757051
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 4.820295013530096e-6 from total weight 1.380724305754562: operator might not be hermitian?
│   α = 0.05175251300778737 + 4.820295013530096e-6im
│   β₁ = 0.9166991496757051
│   β₂ = 1.0312050983482381
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 4.880373024380991e-6 from total weight 1.4083502233166574: operator might not be hermitian?
│   α = 0.4869437492653818 + 4.880373024380991e-6im
│   β₁ = 1.0312050983482381
│   β₂ = 0.8264092095833846
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 2.4708468538414685e-6 from total weight 1.3679952162146665: operator might not be hermitian?
│   α = 0.3365703089542389 + 2.4708468538414685e-6im
│   β₁ = 0.8264092095833846
│   β₂ = 1.0369084612569812
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 4.845211715443637e-6 from total weight 1.9535857899234355: operator might not be hermitian?
│   α = 1.648360457701184 + 4.845211715443637e-6im
│   β₁ = 0.697422390792335
│   β₂ = 0.782947794476709
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 2.3929177318110506e-6 from total weight 1.6460071714797488: operator might not be hermitian?
│   α = 1.0014749966958667 + 2.3929177318110506e-6im
│   β₁ = 0.782947794476709
│   β₂ = 1.0456482155457962
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 3.4887348726309053e-6 from total weight 2.297126647611539: operator might not be hermitian?
│   α = 1.9544737849469627 + 3.4887348726309053e-6im
│   β₁ = 1.0456482155457962
│   β₂ = 0.602862230062428
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 3.893633891298431e-6 from total weight 1.51990288741988: operator might not be hermitian?
│   α = 1.286400942835462 + 3.893633891298431e-6im
│   β₁ = 0.602862230062428
│   β₂ = 0.5402171165451158
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 4.977927394594905e-6 from total weight 1.3906451753871552: operator might not be hermitian?
│   α = 1.2402536191060138 + 4.977927394594905e-6im
│   β₁ = 0.5402171165451158
│   β₂ = 0.3222272972435882
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 5.626639152545021e-6 from total weight 1.5689756337692407: operator might not be hermitian?
│   α = 1.4337246333918294 + 5.626639152545021e-6im
│   β₁ = 0.3222272972435882
│   β₂ = 0.5498070423762209
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 2.929774392332547e-6 from total weight 1.435283449174042: operator might not be hermitian?
│   α = 1.3258019443408067 + 2.929774392332547e-6im
│   β₁ = 0.5498070423762209
│   β₂ = 1.8809234760612477e-21
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
[ Info: CG: converged after 74 iterations: f = -3.630108220530, ‖∇f‖ = 8.8587e-06
┌ Info: Groundstate energy:
│     * numerical: -1.8150541102648536
└     * analytic: -2.190038374277775

Symmetries

The Hubbard model has a rich symmetry structure, which can be exploited to speed up simulations. Apart from the fermionic parity, the model also has a $U(1)$ particle number symmetry, along with a $SU(2)$ spin symmetry. Explicitly imposing these symmetries on the tensors can greatly reduce the computational cost of the simulation.

Naively imposing these symmetries however, is not compatible with our desire to work at half-filling. By construction, imposing symmetries restricts the optimization procedure to a single symmetry sector, which is the trivial sector. In order to work at half-filling, we need to effectively inject one particle per site. In MPSKit, this is achieved by the add_physical_charge function, which shifts the physical spaces of the tensors to the desired charge sector.

H_u1_su2 = hubbard_model(ComplexF64, U1Irrep, SU2Irrep, InfiniteChain(2); U, t, mu=U / 2);
charges = fill(FermionParity(1) ⊠ U1Irrep(1) ⊠ SU2Irrep(0), 2);
H_u1_su2 = MPSKit.add_physical_charge(H_u1_su2, dual.(charges));

pspaces = H_u1_su2.data.pspaces
vspaces = [oneunit(eltype(pspaces)), first(pspaces)]
psi = InfiniteMPS(pspaces, vspaces)
psi = compute_groundstate(psi, H_u1_su2; svalue=1e-3, expansionfactor=1 / 3,
                          expansioniter=20)
E = real(expectation_value(psi, H_u1_su2)) / 2
@info """
Groundstate energy:
    * numerical: $E
    * analytic: $(hubbard_energy(U / 4) - U / 4)
"""

[ Info: CG: converged after 51 iterations: f = -4.379995918984, ‖∇f‖ = 9.4208e-06
┌ Info: Groundstate energy:
│     * numerical: -2.189997959491911
└     * analytic: -2.190038374277775

Excitations

Because of the integrability, it is known that the Hubbard model has a rich excitation spectrum. The elementary excitations are known as spinons and holons, which are domain walls in the spin and charge sectors, respectively. The fact that the spin and charge sectors are separate is a phenomenon known as spin-charge separation.

The domain walls can be constructed by noticing that there are two equivalent groundstates, which differ by a translation over a single site. In other words, the groundstates are $\psi_{AB}` and$\psi_{BA}$, where$A$and$B`` are the two sites. These excitations can be constructed as follows:

alg = QuasiparticleAnsatz(; tol=1e-3)
momenta = range(-π, π; length=33)
psi_AB = psi
envs_AB = environments(psi_AB, H_u1_su2);
psi_BA = circshift(psi, 1)
envs_BA = environments(psi_BA, H_u1_su2);

spinon_charge = FermionParity(0) ⊠ U1Irrep(0) ⊠ SU2Irrep(1 // 2)
E_spinon, ϕ_spinon = excitations(H_u1_su2, alg, momenta,
                                 psi_AB, envs_AB, psi_BA, envs_BA;
                                 sector=spinon_charge, num=1);

holon_charge = FermionParity(1) ⊠ U1Irrep(1) ⊠ SU2Irrep(0)
E_holon, ϕ_holon = excitations(H_u1_su2, alg, momenta,
                               psi_AB, envs_AB, psi_BA, envs_BA;
                               sector=holon_charge, num=1);

[ Info: Found excitations for momentum = 1.5707963267948966
[ Info: Found excitations for momentum = -1.5707963267948966
[ Info: Found excitations for momentum = -2.945243112740431
[ Info: Found excitations for momentum = -3.141592653589793
[ Info: Found excitations for momentum = -0.19634954084936207
[ Info: Found excitations for momentum = 3.141592653589793
[ Info: Found excitations for momentum = 2.748893571891069
[ Info: Found excitations for momentum = -0.39269908169872414
[ Info: Found excitations for momentum = -2.748893571891069
[ Info: Found excitations for momentum = 2.945243112740431
[ Info: Found excitations for momentum = -2.552544031041707
[ Info: Found excitations for momentum = 1.9634954084936207
[ Info: Found excitations for momentum = 0.0
[ Info: Found excitations for momentum = 0.19634954084936207
[ Info: Found excitations for momentum = 2.552544031041707
[ Info: Found excitations for momentum = -1.7671458676442586
[ Info: Found excitations for momentum = -1.3744467859455345
[ Info: Found excitations for momentum = 0.39269908169872414
[ Info: Found excitations for momentum = -2.356194490192345
[ Info: Found excitations for momentum = 0.5890486225480862
[ Info: Found excitations for momentum = -0.5890486225480862
[ Info: Found excitations for momentum = 1.7671458676442586
[ Info: Found excitations for momentum = 0.7853981633974483
[ Info: Found excitations for momentum = -0.9817477042468103
[ Info: Found excitations for momentum = 2.356194490192345
[ Info: Found excitations for momentum = 1.3744467859455345
[ Info: Found excitations for momentum = -0.7853981633974483
[ Info: Found excitations for momentum = -1.9634954084936207
[ Info: Found excitations for momentum = -1.1780972450961724
[ Info: Found excitations for momentum = 0.9817477042468103
[ Info: Found excitations for momentum = 2.1598449493429825
[ Info: Found excitations for momentum = -2.1598449493429825
[ Info: Found excitations for momentum = 1.1780972450961724
[ Info: Found excitations for momentum = 0.0
[ Info: Found excitations for momentum = -3.141592653589793
[ Info: Found excitations for momentum = 3.141592653589793
[ Info: Found excitations for momentum = -1.3744467859455345
[ Info: Found excitations for momentum = 1.3744467859455345
[ Info: Found excitations for momentum = 2.1598449493429825
[ Info: Found excitations for momentum = -1.7671458676442586
[ Info: Found excitations for momentum = 1.7671458676442586
[ Info: Found excitations for momentum = -1.5707963267948966
[ Info: Found excitations for momentum = 1.5707963267948966
[ Info: Found excitations for momentum = 1.1780972450961724
[ Info: Found excitations for momentum = -0.39269908169872414
[ Info: Found excitations for momentum = 0.7853981633974483
[ Info: Found excitations for momentum = -2.1598449493429825
[ Info: Found excitations for momentum = -2.552544031041707
[ Info: Found excitations for momentum = 1.9634954084936207
[ Info: Found excitations for momentum = 0.9817477042468103
[ Info: Found excitations for momentum = -1.1780972450961724
[ Info: Found excitations for momentum = -0.9817477042468103
[ Info: Found excitations for momentum = 2.748893571891069
[ Info: Found excitations for momentum = -2.748893571891069
[ Info: Found excitations for momentum = -1.9634954084936207
[ Info: Found excitations for momentum = 2.552544031041707
[ Info: Found excitations for momentum = 0.5890486225480862
[ Info: Found excitations for momentum = -2.356194490192345
[ Info: Found excitations for momentum = -0.7853981633974483
[ Info: Found excitations for momentum = -2.945243112740431
[ Info: Found excitations for momentum = 2.945243112740431
[ Info: Found excitations for momentum = -0.5890486225480862
[ Info: Found excitations for momentum = -0.19634954084936207
[ Info: Found excitations for momentum = 0.19634954084936207
[ Info: Found excitations for momentum = 0.39269908169872414
[ Info: Found excitations for momentum = 2.356194490192345

Again, we can compare the numerical results to the analytic solution. Here, the formulae for the excitation energies are expressed in terms of dressed momenta:

function spinon_momentum(Λ, u; rtol=1e-12)
    integrandum(ω) = besselj0(ω) * sin(ω * Λ) / ω / cosh(ω * u)
    return π / 2 - quadgk(integrandum, 0, Inf; rtol=rtol)[1]
end
function spinon_energy(Λ, u; rtol=1e-12)
    integrandum(ω) = besselj1(ω) * cos(ω * Λ) / ω / cosh(ω * u)
    return 2 * quadgk(integrandum, 0, Inf; rtol=rtol)[1]
end

function holon_momentum(k, u; rtol=1e-12)
    integrandum(ω) = besselj0(ω) * sin(ω * sin(k)) / ω / (1 + exp(2u * abs(ω)))
    return π / 2 - k - 2 * quadgk(integrandum, 0, Inf; rtol=rtol)[1]
end
function holon_energy(k, u; rtol=1e-12)
    integrandum(ω) = besselj1(ω) * cos(ω * sin(k)) * exp(-ω * u) / ω / cosh(ω * u)
    return 2 * cos(k) + 2u + 2 * quadgk(integrandum, 0, Inf; rtol=rtol)[1]
end

Λs = range(-10, 10; length=51)
P_spinon_analytic = rem2pi.(spinon_momentum.(Λs, U / 4), RoundNearest)
E_spinon_analytic = spinon_energy.(Λs, U / 4)
I_spinon = sortperm(P_spinon_analytic)
P_spinon_analytic = P_spinon_analytic[I_spinon]
E_spinon_analytic = E_spinon_analytic[I_spinon]
P_spinon_analytic = [reverse(-P_spinon_analytic); P_spinon_analytic]
E_spinon_analytic = [reverse(E_spinon_analytic); E_spinon_analytic];

ks = range(0, 2π; length=51)
P_holon_analytic = rem2pi.(holon_momentum.(ks, U / 4), RoundNearest)
E_holon_analytic = holon_energy.(ks, U / 4)
I_holon = sortperm(P_holon_analytic)
P_holon_analytic = P_holon_analytic[I_holon]
E_holon_analytic = E_holon_analytic[I_holon];

p = let p_excitations = plot(; xaxis="momentum", yaxis="energy")
    scatter!(p_excitations, momenta, real(E_spinon); label="spinon")
    plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label="spinon (analytic)")

    scatter!(p_excitations, momenta, real(E_holon); label="holon")
    plot!(p_excitations, P_holon_analytic, E_holon_analytic; label="holon (analytic)")

    p_excitations
end

The plot shows some discrepancies between the numerical and analytic results. First and foremost, we must realize that in the thermodynamic limit, the momentum of a domain wall is actually not well-defined. Concretely, only the difference in momentum between the two groundstates is well-defined, as we can always shift the momentum by multiplying one of the groundstates by a phase. Here, we can fix this shift by realizing that our choice of shifting the groundstates by a single site, differs from the formula by a factor $\pi/2$.

momenta_shifted = rem2pi.(momenta .- π / 2, RoundNearest)
p = let p_excitations = plot(; xaxis="momentum", yaxis="energy", xlims=(-π, π))
    scatter!(p_excitations, momenta_shifted, real(E_spinon); label="spinon")
    plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label="spinon (analytic)")

    scatter!(p_excitations, momenta_shifted, real(E_holon); label="holon")
    plot!(p_excitations, P_holon_analytic, E_holon_analytic; label="holon (analytic)")

    p_excitations
end

The second discrepancy is that while the spinon dispersion is well-reproduced, the holon dispersion is not. This is due to the fact that the excitation ansatz captures the lowest-energy excitation, and not the elementary single-particle excitation. To make this explicit, we can consider the scattering states comprising of a holon and two spinons. If these are truly scattering states, the energy of the scattering state should be the sum of the energies of the individual excitations, and the momentum is the sum of the momenta. Thus, we can find the lowest-energy scattering states by minimizing the energy over the combination of momenta for the constituent elementary excitations.

holon_dispersion_itp = linear_interpolation(P_holon_analytic, E_holon_analytic;
                                            extrapolation_bc=Line())
spinon_dispersion_itp = linear_interpolation(P_spinon_analytic, E_spinon_analytic;
                                             extrapolation_bc=Line())
function scattering_energy(p1, p2, p3)
    p1, p2, p3 = rem2pi.((p1, p2, p3), RoundNearest)
    return holon_dispersion_itp(p1) + spinon_dispersion_itp(p2) + spinon_dispersion_itp(p3)
end;

E_scattering_min = map(momenta_shifted) do p
    e = Inf
    for i in 1:10 # repeat for stability
        res = optimize((rand(2) .* (2π) .- π)) do (p₁, p₂)
            p₃ = p - p₁ - p₂
            return scattering_energy(p₁, p₂, p₃)
        end

        e = min(Optim.minimum(res), e)
    end
    return e
end
E_scattering_max = map(momenta_shifted) do p
    e = -Inf
    for i in 1:10 # repeat for stability
        res = optimize((rand(Float64, 2) .* (2π) .- π)) do (p₁, p₂)
            p₃ = p - p₁ - p₂
            return -scattering_energy(p₁, p₂, p₃)
        end

        e = max(-Optim.minimum(res), e)
    end
    return e
end;

p = let p_excitations = plot(; xaxis="momentum", yaxis="energy", xlims=(-π, π),
                             ylims=(-0.1, 5))
    scatter!(p_excitations, momenta_shifted, real(E_spinon); label="spinon")
    plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label="spinon (analytic)")

    scatter!(p_excitations, momenta_shifted, real(E_holon); label="holon")
    plot!(p_excitations, P_holon_analytic, E_holon_analytic; label="holon (analytic)")

    I = sortperm(momenta_shifted)
    plot!(p_excitations, momenta_shifted[I], E_scattering_min[I]; label="scattering states",
          fillrange=E_scattering_max[I], fillalpha=0.3, fillstyle=:x)

    p_excitations
end

This page was generated using Literate.jl.