MPSKit.jl

Appearance

Quasiparticle excitations ​

The previous tutorials ended with a ground state: the lowest-energy state of the transverse-field Ising model, first on a finite chain and then directly in the thermodynamic limit. The natural next question is what lies above it: how much energy does it cost to excite the system? For a translation-invariant chain the answer is organized by momentum — for each momentum there is a lowest excitation energy , and the resulting curve is the dispersion relation of the model. Its minimum over all momenta is the energy gap, one of the most basic characterizations of a quantum phase. <!– REVIEW: physics framing — "lowest excitation energy at each momentum defines the dispersion relation, whose minimum is the gap" is standard for translation-invariant systems, but please confirm this phrasing is acceptable as stated. –>

In this tutorial we compute the dispersion relation of the infinite transverse-field Ising chain with MPSKit's quasiparticle ansatz, and finish with a plot of across the Brillouin zone — compared against the exact solution.

Loading the packages ​

As in the previous tutorials, every code block on this page shares one Julia session, so we load the packages once.

julia
using MPSKit, MPSKitModels, TensorKit
using Plots

1. Find the ground state ​

Excitations are computed on top of a ground state, so the first step is the calculation you already know from The thermodynamic limit: build the infinite Hamiltonian, make a random InfiniteMPS, and converge it with VUMPS.

This time we set the field to g = 2.0, deep in the paramagnetic phase, where the model is gapped: the lowest excitation costs a finite amount of energy, which is exactly what we want to measure. <!– REVIEW: physics claim — TFIM at g = 2 is in the gapped paramagnetic (disordered) phase; the transition is at g = 1. –>

julia
g = 2.0
H = transverse_field_ising(; g)
ψ₀ = InfiniteMPS(ℂ^2, ℂ^12)
ψ, envs, ϵ = find_groundstate(ψ₀, H, VUMPS(; verbosity = 0))
(InfiniteMPS{TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}}(TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}[TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.42853030561648936 + 0.5522043709784232im, 0.10926092748794423 + 0.06114540663145021im, 0.08250101819489453 + 0.035278113354461724im, -0.008242947358678145 - 0.015477155056059275im, 0.08440547845492 - 0.0228118197023556im, -0.00974599566789386 + 0.03516452930579776im, -0.012817885896088147 + 0.012421751571305916im, 0.003130631075264648 + 0.013657614737757562im, -0.0272845994967541 - 0.02612199749235388im, 0.004832797745529021 + 0.07612467620252321im  …  -0.2307570649523536 + 0.21910570033645713im, 0.18772181535510304 - 0.14957166774020614im, -0.0022240894209827424 + 0.06421715959182803im, 0.15823532554643238 - 0.05170274861372345im, 0.31083324532694473 + 0.01652434564723874im, -0.18103630085908673 - 0.10714486043681204im, -0.06485373718423737 - 0.0022762571505187492im, 0.08500441014450942 - 0.17757649034401699im, 0.27831501262449904 + 0.06249323233732048im, 0.1959255524158215 - 0.10098255810523535im], (ℂ^12 ⊗ ℂ^2) ← ℂ^12)], TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}[TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.49614524618256417 + 0.5018139927394476im, -0.170081319122884 + 0.6587224367233614im, 0.1489259403571852 + 0.0720259485949183im, 0.020001260475815456 + 0.09829161153807245im, -0.007590846281562984 - 0.03308858138409474im, 0.014892999248544859 + 0.01788377352651017im, -0.0047652959473761065 - 0.0003475471065039125im, -0.004937742712259875 + 0.002149193746261737im, 0.00042808942785697537 - 0.001937273180872961im, 0.0003038668191863187 + 0.0003665307013296668im  …  -1.2367569523513612e-6 + 1.3606925174788282e-6im, 1.1467553595122055e-5 - 3.1207376352601783e-6im, 7.093545839957616e-6 + 2.9694657758568123e-5im, 0.0007391569258114913 - 0.0006012588449221609im, 0.003941349191100615 + 0.0005119971772077444im, -0.005200601887047545 + 4.66366825065194e-5im, 0.023487313140735454 + 0.012445087391549276im, 0.023524833674635343 + 0.00942971592356874im, 0.038114378800696826 + 0.026103774077376445im, 0.0005121980094439706 + 0.016008466990638487im], (ℂ^12 ⊗ ℂ^2) ← ℂ^12)], TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}[TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.9628982101687036 + 0.0im, 0.08230732034273702 - 0.0815189668597753im, 0.0998235446626757 - 0.017714404568099786im, -0.06021263833274354 - 0.05450603141549537im, 0.05476868170963711 - 0.1065782741875344im, 0.0263458067077055 + 0.005546664834898234im, 0.03338737771423729 - 0.016471924524305497im, 0.02027514772576903 + 0.00865347918805032im, -0.024989887462760246 - 0.021082633739052847im, 0.06264065915846681 + 0.006125625720914012im  …  0.0 - 0.0im, 0.0 - 0.0im, 0.0 - 0.0im, 0.0 - 0.0im, 0.0 - 0.0im, 0.0 - 0.0im, 0.0 - 0.0im, 0.0 - 0.0im, 0.0 - 0.0im, 9.86604222429129e-9 + 0.0im], ℂ^12 ← ℂ^12)], TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}[TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.4777373695329012 + 0.4831957954464246im, 0.06042102374680439 + 0.08344029759559829im, 0.03746731797881165 + 0.02799127031178028im, -0.0035975307994301443 - 0.06979317652458972im, 0.079178164379453 - 0.029765446613518383im, 0.0002467569220941605 + 0.016075265283805304im, 0.023386775609607856 + 0.007909693479893114im, 0.003989140656024261 + 0.014986593971884975im, 0.001292082787836938 - 0.029435943635811243im, 0.030472194789236975 + 0.03362860199253573im  …  -2.2766589463734482e-9 + 2.1617060911024003e-9im, 1.8520713567140596e-9 - 1.4756803894825411e-9im, -2.1942960138015302e-11 + 6.335692080570277e-10im, 1.56115640321558e-9 - 5.101015009349134e-10im, 3.06669392310913e-9 + 1.630298918844414e-10im, -1.786111788405251e-9 - 1.0570957171853847e-9im, -6.398497094627759e-10 - 2.2457649160362952e-11im, 8.386570997367048e-10 - 1.751977151775526e-9im, 2.7458676662074707e-9 + 6.165608689724497e-10im, 1.933009772952091e-9 - 9.962981821832005e-10im], (ℂ^12 ⊗ ℂ^2) ← ℂ^12)]), MPSKit.InfiniteEnvironments{BlockTensorKit.BlockTensorMap{TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}, ComplexF64, TensorKit.ComplexSpace, 2, 1, 3}}(BlockTensorKit.BlockTensorMap{TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}, ComplexF64, TensorKit.ComplexSpace, 2, 1, 3}[BlockTensorKit.BlockTensorMap{TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}, ComplexF64, TensorKit.ComplexSpace, 2, 1, 3}(TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}[TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im], (ℂ^12 ⊗ (ℂ^1)') ← ℂ^12) TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[-0.053311963206765796 - 1.1191852949027205e-17im, -0.358251736527705 - 0.5843617581236721im, 0.05078770391717339 + 0.0020020264803069193im, 0.03140320407758463 + 0.020204083934461416im, 0.02113562929331718 - 0.033110409037399494im, 0.02784276176987425 - 0.06347295163711975im, 0.00564235478392521 + 0.05782980197192009im, -0.016792117320959185 - 0.02528540597389974im, -0.01667518844875291 + 0.05534914006385503im, -0.02478566590881191 + 0.042620129699132585im  …  -0.06008823568456112 - 0.05847307940714678im, 0.006625116333653801 + 0.06682149600698672im, -0.23890155318621956 + 0.1443530045323255im, 0.44220795633890164 + 0.08268604374140834im, 0.00945009561038055 - 0.0018612138469955854im, -0.1952124146421941 - 0.03205691807310841im, -0.0828407759216016 - 0.1515646899933442im, -0.05056911193898503 + 0.06208179263375756im, 0.08946518961360396 - 0.027054575206560048im, 0.02288913634675515 - 9.878828038172291e-18im], (ℂ^12 ⊗ (ℂ^1)') ← ℂ^12) TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.267700208657925 - 1.0862600033481454e-16im, -0.41141236681441395 + 0.3407360236289823im, -0.5649383573911897 + 0.12597614820295006im, 0.36113011909190956 + 0.42225501540357396im, -0.34323551012734654 + 0.8371249542015499im, 0.18194027158588746 + 0.015402135405870114im, -0.26250139983595067 + 0.06416611310552961im, -0.08074822927167143 - 0.12348610323437025im, 0.21011820884597268 + 0.00016292666024734576im, -0.38981164353391345 - 0.008608589134343727im  …  0.27150008412747334 - 0.4135408423468176im, 0.219998052300909 + 0.030603320744225608im, 0.2967517780387321 + 0.3838582879731897im, -0.12776117715889937 - 0.08949920535578022im, -1.055255639576072 + 0.11040592205664812im, 0.9016235115684244 - 0.22216483950360272im, -1.111413814095592 - 0.22011335619445244im, -0.1308040336922776 - 0.509005125619096im, 0.0773651026469476 - 0.3530951201319407im, 5.800809984097201 - 1.3569633453373572e-15im], (ℂ^12 ⊗ (ℂ^1)') ← ℂ^12);;;], (⊞(ℂ^12) ⊗ ((ℂ^1)' ⊞ (ℂ^1)' ⊞ (ℂ^1)')) ← ⊞(ℂ^12))], BlockTensorKit.BlockTensorMap{TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}, ComplexF64, TensorKit.ComplexSpace, 2, 1, 3}[BlockTensorKit.BlockTensorMap{TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}, ComplexF64, TensorKit.ComplexSpace, 2, 1, 3}(TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}[TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[-0.06264671482047657 + 8.315306845948668e-18im, -0.033113190965104074 + 0.04467149826100846im, 0.010857033704262813 + 0.01149260323792192im, -0.013469338889962918 - 0.009542057769389082im, 0.000491101403614739 - 0.002098809637669323im, -0.0024443019383136888 + 0.00017752862102483866im, 0.00021425207644883436 + 0.00036381594968610414im, 0.00017786249741852102 - 0.00039291949060648847im, 2.5085002829184418e-5 - 1.614708106089312e-6im, 6.319009041193194e-5 - 3.707316572589322e-6im  …  0.005530159999667325 + 0.0005167319432060265im, -0.010347102020289023 + 0.031024865883165068im, -0.005624252348371371 - 0.0532462676982797im, 0.009767989048132943 + 0.14637699533512355im, -0.28496269776275196 + 0.3211728450551784im, 0.27367849307393827 - 0.16098724635210845im, -0.8199992099488137 + 0.8062759386213453im, -0.5452799881323059 + 0.08038050283528242im, -0.9091909977331638 - 0.14119788284416532im, 6.8583495153436775 - 5.622561183116999e-16im], (ℂ^12 ⊗ ℂ^1) ← ℂ^12) TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.016604552058207198 + 4.512919840608712e-17im, 0.7204228086902094 + 1.203000676289737im, 0.17275778175423784 - 0.05286436391704251im, 0.014607826360677539 + 0.006407427409276969im, -0.02018806587791769 - 0.013646488347112145im, 0.0006610098064023305 + 0.00267779961031096im, 0.0002131597738625616 + 0.0029036702657698815im, -0.0001234188248281301 - 0.00021717570621613937im, -0.00012885257912906162 - 0.0003698679552948068im, 0.0002685352096943389 - 0.00015237950492903352im  …  0.0024940907629097267 + 0.002941784752325173im, 0.0014670123563556025 + 0.001229125204971041im, -0.006057584442412711 + 0.011561307785657903im, -0.017512890821162943 - 0.027301405184900023im, 0.08104028854900334 - 0.11168649540673306im, -0.004224873642330732 + 0.07631637726442288im, 0.2130879535965599 - 0.015743194268480386im, 0.32867132809225136 + 0.017819655247381234im, 0.36686221226300025 + 0.1759602271281763im, 0.5508550456934962 - 1.521040494344052e-17im], (ℂ^12 ⊗ ℂ^1) ← ℂ^12) TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im], (ℂ^12 ⊗ ℂ^1) ← ℂ^12);;;], (⊞(ℂ^12) ⊗ (ℂ^1 ⊞ ℂ^1 ⊞ ℂ^1)) ← ⊞(ℂ^12))]), 9.815144034293853e-11)

We keep all three return values this time: the optimized state ψ and the environments envs both feed directly into the excitation calculation below, so nothing has to be recomputed.

2. One excitation at one momentum ​

The quasiparticle ansatz builds an excited state directly on top of the uniform ground state. The idea is simple to picture: take the converged ground state and perturb it locally, replacing the tensor at one site with a new one that we get to optimize. Because the chain is infinite and translation invariant, we do not place this perturbation at any particular site; instead we superpose it across all sites with a plane-wave phase, which gives the excitation a definite momentum . Optimizing the perturbation then yields the lowest excited state at that momentum. <!– REVIEW: quasiparticle-ansatz framing — "local perturbation of the uniform ground state, momentum-superposed over all sites, then variationally optimized" — please confirm this hand-held description does not oversimplify in a misleading way. –>

The call is excitations with the QuasiparticleAnsatz algorithm, a momentum (a real number, in radians per site), and the ground state with its environments. Let us ask for the excitation at the edge of the Brillouin zone,  :

julia
E, ϕ = excitations(H, QuasiparticleAnsatz(), π, ψ, envs)
E
1-element Vector{Float64}:
 6.000000000013939

Two things to note about the return values:

  • E is a vector of excitation energies, of length num — the keyword controlling how many excitations to compute at this momentum, which defaults to num = 1, so here it has a single entry.

  • The entries of E are energies above the ground state — gaps at this momentum, not total energies. The ground-state energy is subtracted internally, so you can read them off directly.

The second return value ϕ holds the corresponding quasiparticle states, which can be used for further post-processing; we will not need them in this tutorial.

3. The full dispersion ​

To trace out the whole dispersion relation we simply pass a range of momenta instead of a single number. By symmetry it is enough to scan from to , and we use 16 points to keep the runtime modest. <!– REVIEW: physics claim — restricting the scan to [0, π] uses the k → −k symmetry of the TFIM dispersion. –>

julia
momenta = range(0, π, 16)
Es, ϕs = excitations(H, QuasiparticleAnsatz(), momenta, ψ, envs; verbosity = 0)
size(Es)
(16, 1)

With a range of momenta the energies come back as a matrix of size (length(momenta), num) — here (16, 1), one row per momentum and one column because we kept the default num = 1. We pass verbosity = 0 to silence the progress line this method otherwise prints for every momentum. The momenta are independent of one another, so MPSKit works on them in parallel by default.

4. Plot the dispersion ​

Now for the payoff. This particular model is exactly solvable, so we can plot our numerical dispersion right on top of the known answer:

<!– REVIEW: exact dispersion — the single-particle dispersion of the TFIM in the Pauli-matrix convention H = -(Σ σᶻσᶻ + g Σ σˣ) used by transverse_field_ising, obtained from its free-fermion (Jordan-Wigner) solution; at k = 0 it gives the gap 2(g − 1) for g > 1. Please verify formula and convention. –>

For this Hermitian problem the computed energies come back as real numbers; the real.(...) below is a harmless safeguard for the general case, where the eigenvalue solver may return a complex number type with numerically vanishing imaginary parts.

julia
k_exact = range(0, π, 200)
ΔE_exact = @. 2 * sqrt(1 + g^2 - 2g * cos(k_exact))
plot(k_exact, ΔE_exact; label = "exact", xlabel = "momentum k", ylabel = "ΔE(k)", title = "TFIM dispersion (g = $g)")
scatter!(momenta, real.(Es); label = "quasiparticle ansatz (D = 12)")

The 16 computed points fall right on the exact curve. The dispersion rises monotonically from   to  , so its minimum — the gap — sits at zero momentum, where the exact value is   . Our first matrix entry is precisely that point, so we can close with a numerical check:

julia
real(Es[1, 1]), 2 * (g - 1)
(2.00000000006439, 2.0)

A ground state at bond dimension 12 plus a variational quasiparticle on top reproduces the exact gap of the model — that is the quasiparticle ansatz working as intended. <!– REVIEW: convergence claim — please confirm that D = 12 is ample for the TFIM quasiparticle gap at g = 2 to the displayed precision, so promising agreement here is safe. –>

Where to go next ​

You have computed a full dispersion relation on top of an infinite ground state and read off the energy gap.

The excitations entry point can do considerably more than what we used here: it can target excitations carrying a nontrivial symmetry charge, build topological (domain-wall) excitations that interpolate between two different ground states, and compute excited states of finite chains, where momentum is no longer a good quantum number and different algorithms take over. All of these are recipes in Excited states. For what each excitation algorithm actually does and when to choose it, see the library reference Excitations.

<!– NOTE(other pages, do not edit here): tutorials/thermodynamic_limit.md line 134 has "TODO(link): excitations tutorial" which can now point to (@ref tutorial_excitations); the orchestrator wires that separately. This page also needs a nav entry in docs/make.jl. –>