Skip to content

The thermodynamic limit

In Your first ground state we put the transverse-field Ising model on a finite chain of L = 16 sites. That is a perfectly good calculation, but it carries two prices: the open ends of the chain are physically different from its middle (boundary effects), and every quantity we measured still depends on the length L (finite-size effects). To read off the true physics of the model we would have to repeat the calculation at several lengths and extrapolate L → ∞.

MPSKit lets you skip that extrapolation and work directly at L = ∞. The trick is translation invariance: instead of storing one tensor per site, we store a single tensor and imagine it repeated forever along the chain — an InfiniteMPS. There are no ends, so there are no boundary effects, and there is no L to extrapolate. Best of all, as you are about to see, the code barely changes: the same model, the same workflow, two edits.

Infinite states are always normalized

An InfiniteMPS is normalized to 1 by construction, and you cannot choose otherwise. Any other normalization would make expectation values either blow up or vanish as the (infinite) chain length is taken to infinity, so per-site quantities are the only ones that make sense here.

Loading the packages

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

julia
using MPSKit, MPSKitModels, TensorKit
using TensorKitTensors.SpinOperators: σˣ, σᶻ
using Plots

1. Build the Hamiltonian and initial state

Here are the only two lines that differ from the finite tutorial.

For the Hamiltonian, we drop the lattice argument. Where the finite version wrote transverse_field_ising(FiniteChain(L); g = 0.5), we simply omit FiniteChain(L): with no lattice, transverse_field_ising builds the Hamiltonian for the infinite chain.

julia
H = transverse_field_ising(; g = 0.5)
1-site InfiniteMPOHamiltonian(ComplexF64, TensorKit.ComplexSpace) with maximal dimension 3:
| ⋮
| (ℂ^1 ⊞ ℂ^1 ⊞ ℂ^1)
┼─[1]─ ℂ^2
│ (ℂ^1 ⊞ ℂ^1 ⊞ ℂ^1)
| ⋮

For the state, we swap FiniteMPS for InfiniteMPS. There is no length to pass, so the constructor takes just the physical and virtual spaces — the physical space ℂ^2 of a spin-1/2 site and the bond space ℂ^D whose dimension D is again the accuracy knob.

julia
D = 4
ψ₀ = InfiniteMPS(ℂ^2, ℂ^D)
1-site InfiniteMPS(ComplexF64, TensorKit.ComplexSpace) with maximal dimension 4:
| ⋮
| ℂ^4
├─[1]─ ℂ^2
│ ℂ^4
| ⋮

That is the whole difference. The bond dimension means exactly what it did on the finite chain (see Controlling bond dimension), and ℂ^2/ℂ^D are the same physical/virtual spaces.

InfiniteMPS also accepts bare integers

Unlike FiniteMPS, the infinite constructor happily takes plain integers: InfiniteMPS(2, D) is equivalent to InfiniteMPS(ℂ^2, ℂ^D). We stick with the explicit spaces to match the rest of the documentation.

2. Find the ground state

We optimize with VUMPS, the infinite-chain workhorse, passing it explicitly so it is visible.

julia
ψ, envs, ϵ = find_groundstate(ψ₀, H, VUMPS())
(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.2697844904199237 - 0.17140005961620117im, 0.10463202821559439 - 0.5269242785940307im, -0.2954251486026092 + 0.06599635300923569im, 0.0008790087492598284 - 0.12138258713104633im, -0.32906639101740537 - 0.3472276612139481im, -0.11239447210794645 + 0.3777440045021627im, -0.1863125150932261 - 0.28818343022431125im, 0.008804752124413145 - 0.02934100094418994im, -0.3742078798456634 + 0.08839216894519714im, -0.2829183673748616 - 0.23444753019145803im  …  -0.3129935400521174 - 0.0008884787301747521im, 0.07883928835301039 - 0.11476272773408097im, -0.4163753987127318 - 0.2906256682271941im, -0.4177900349208379 + 0.5486157258531549im, -0.1320421442876951 - 0.2636457001745064im, -0.08247838700885002 - 0.029943155814042936im, -0.0820678844609076 + 0.001996340858228862im, -0.2996604670374363 + 0.25658979193161174im, 0.07764822063557134 - 0.014959763638265774im, -0.05049387750431271 + 0.028645211864645984im], (ℂ^4 ⊗ ℂ^2) ← ℂ^4)], TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}[TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[-0.7815558014465971 - 0.6080114756072706im, -0.029417170634703442 - 0.13455231261506692im, 0.019291766247287224 - 0.0037914378955690013im, -0.010061899862282432 - 0.004641867738797722im, -0.11205051053040363 - 0.08049231455802044im, 0.23899314810465894 + 0.9251151741224246im, -0.24712618660084068 + 0.06617098309088693im, -0.012986174019472108 + 0.0008932374090170891im, -0.009664107452328135 + 0.005316457112292908im, -0.016629695591154903 - 0.011100660065953474im  …  -0.18942301916780252 - 0.15435211003675078im, -0.01848918093123513 - 0.023444581412601057im, -4.220867134365101e-6 - 2.9461210609836875e-6im, -0.0001497623276700045 + 0.0003389654395905836im, 0.0038131737282211327 + 0.001835782777705243im, -0.04143879619130808 - 0.03297620964336972im, -8.319358861710894e-7 + 2.0237241545824245e-8im, -9.610345545065498e-5 + 0.00011164354689714369im, 0.001491280923167713 + 0.001457365551621201im, -0.0016450964011766405 - 0.0011717211348618202im], (ℂ^4 ⊗ ℂ^2) ← ℂ^4)], TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}[TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.6840652403222441 + 0.0im, 0.3004421854013194 + 0.40961331944906654im, 0.42049879396602985 - 0.24508760027611232im, 0.1331909894087692 + 0.13772779620490025im, 0.0 + 0.0im, 0.018935350350719443 + 0.0im, -0.005871887149844173 + 0.0021273980687485515im, 0.0020929754125719267 - 0.0028475238189079196im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0003731826019313016 + 0.0im, -0.00018051634330144305 + 4.989304240902033e-5im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 6.93448388024926e-6 + 0.0im], ℂ^4 ← ℂ^4)], TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}[TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[-0.5346351571418106 - 0.41591951617996975im, 0.013680241354629037 - 0.5053557578256118im, -0.47719316543461043 - 0.06339237728717219im, -0.020804183286872904 - 0.18881984198789417im, -0.07664985941421078 - 0.05506199450222615im, 0.003831442887944153 - 0.05256328852136812im, -0.07030842498038352 - 0.011283774941263397im, -0.0006623761980322218 - 0.024921884890387727im, -0.006610879986876706 + 0.003636803512183537im, -0.0053960863206497945 - 0.0025714540271239107im  …  -0.0001300741283957023 + 6.243022989501925e-6im, 3.710722411779995e-5 - 5.0517675411157434e-5im, -2.8873484905057975e-6 - 2.015339011508147e-6im, -2.8971582624873256e-6 + 3.8043669073799496e-6im, -9.156441210765685e-7 - 1.828246857957144e-6im, -5.719450451818305e-7 - 2.0764033131627266e-7im, -5.690984218803225e-7 + 1.3843593500871019e-8im, -2.077990678219067e-6 + 1.7793177759862734e-6im, 5.384503343274074e-7 - 1.0373823980189303e-7im, -3.501489796049372e-7 + 1.986397599217124e-7im], (ℂ^4 ⊗ ℂ^2) ← ℂ^4)]), 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, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im], (ℂ^4 ⊗ (ℂ^1)') ← ℂ^4) TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.0040662167527948205 + 0.0im, -0.3611326403100936 - 0.45716179411714275im, -0.1899640619846278 + 0.26385729471843455im, 0.20573595479808188 - 0.1119170389592164im, -0.3611326403100936 + 0.45716179411714275im, 0.1781842768169979 - 1.3877787807814457e-17im, -0.016185927514883222 + 0.3265829071056111im, -0.14000629428877892 - 0.048025796100258866im, -0.1899640619846278 - 0.26385729471843455im, -0.016185927514883222 - 0.3265829071056111im, 0.28161571213005654 + 5.551115123125783e-17im, -0.45009938587138276 + 0.07394377768522864im, 0.20573595479808188 + 0.11191703895921642im, -0.14000629428877892 + 0.048025796100258866im, -0.45009938587138276 - 0.07394377768522861im, -0.4638662056954297 - 4.336808689942018e-18im], (ℂ^4 ⊗ (ℂ^1)') ← ℂ^4) TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[1.0650569897299038 - 1.7313892412859682e-16im, -0.23076849013218334 - 0.22569584796883507im, -1.0348917761860081 + 0.22014524862243162im, -0.8270030564931024 - 0.157427840918757im, -0.23076849013218337 + 0.22569584796883507im, 1.5051779105110117 - 1.5792321966845613e-16im, 0.11778679475306703 + 0.714117005549583im, -0.10965809806865406 - 0.16255920130414492im, -1.0348917761860084 - 0.22014524862243134im, 0.117786794753067 - 0.714117005549583im, 1.9298311464425932 - 2.1551174382066316e-16im, 0.5451293690106087 - 0.24775707319848353im, -0.8270030564931025 + 0.15742784091875728im, -0.10965809806865415 + 0.16255920130414514im, 0.5451293690106087 + 0.24775707319848353im, 2.9640480439962276 - 1.3923097030678674e-16im], (ℂ^4 ⊗ (ℂ^1)') ← ℂ^4);;;], (⊞(ℂ^4) ⊗ ((ℂ^1)' ⊞ (ℂ^1)' ⊞ (ℂ^1)')) ← ⊞(ℂ^4))], 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.004049804289442349 - 1.29892721929335e-19im, 0.18597506877618417 + 0.13020403789322313im, -0.01845296967358574 + 0.02641992183396776im, -0.00032854188983181745 - 0.00014545452558323348im, 0.1859750687761842 - 0.13020403789322316im, 1.856895425774736 - 8.757160474804445e-17im, 0.00514630741787124 + 0.49079465902207065im, 0.007948634841435058 - 0.040048523542333334im, -0.018452969673585718 - 0.026419921833967755im, 0.005146307417871347 - 0.49079465902207076im, 1.8838924910816743 - 2.7669505421248375e-16im, 0.25595337229095433 + 0.0467534580179301im, -0.0003285418898318809 + 0.00014545452558324817im, 0.00794863484143505 + 0.04004852354233332im, 0.2559533722909544 - 0.04675345801792995im, 3.727376152629918 - 5.477185098852882e-16im], (ℂ^4 ⊗ ℂ^1) ← ℂ^4) TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[1.3594243903579712 + 1.1102230246251565e-16im, 0.2967089363036857 + 0.24740266509009773im, -0.032838421165208576 + 0.04001860244189013im, 0.0008563982787859744 - 0.0005018526714359178im, 0.2967089363036857 - 0.24740266509009773im, -1.3594118166942744 + 0.0im, 0.00012276279788626943 + 0.0008974408593880878im, -0.011858405540339043 + 0.050730457380942im, -0.032838421165208576 - 0.04001860244189013im, 0.00012276279788625555 - 0.0008974408593880878im, -1.3725288236140143 + 0.0im, -0.32873428598895915 - 0.07333148669710672im, 0.0008563982787859735 + 0.0005018526714359178im, -0.011858405540339047 - 0.05073045738094201im, -0.32873428598895915 + 0.07333148669710672im, 1.3725162499236192 + 0.0im], (ℂ^4 ⊗ ℂ^1) ← ℂ^4) 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, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im], (ℂ^4 ⊗ ℂ^1) ← ℂ^4);;;], (⊞(ℂ^4) ⊗ (ℂ^1 ⊞ ℂ^1 ⊞ ℂ^1)) ← ⊞(ℂ^4))]), 8.598666533245746e-11)

The lines printed above are VUMPS's per-iteration convergence log, shown at the default verbosity. VUMPS (the variational uniform matrix product state algorithm) optimizes the single repeated tensor directly in the thermodynamic limit, iterating until it reaches a fixed point.

The return value has the same shape as on the finite chain: the optimized state ψ, the reusable envs, and a convergence-error measure ϵ.

The algorithm is optional here too

Just as find_groundstate(ψ₀, H) selected DMRG for a finite input, calling it with no algorithm on an infinite input selects VUMPS automatically. VUMPS accepts the familiar keywords tol (default 1e-10), maxiter (default 200), and verbosity (default 3); we use verbosity = 0 later to silence the log inside a loop. Note there is no find_groundstate! for infinite states — VUMPS returns a fresh state and leaves ψ₀ untouched.

3. Measure observables

For the default single-site unit cell used here, expectation_value(ψ, H) returns the energy of that one-site unit cell, which is exactly the energy per site:

julia
E = expectation_value(ψ, H)
-1.0635444098648656 + 5.54963149165599e-17im

The magnetization is the local order parameter . Because the state is translation-invariant, every site is identical, so we measure it at site 1 of the unit cell:

julia
expectation_value(ψ, 1 => σᶻ())
0.9646786346216887 + 7.917296201993681e-17im

So far these are the same quantities we computed on the finite chain. The infinite setting also unlocks an observable with no finite-chain analogue: the correlation_length, extracted from the transfer-matrix spectrum of the uniform state.

julia
correlation_length(ψ)
0.5833500835890653

The correlation length tells us how far apart two spins can still "feel" each other; it is measured in units of the lattice spacing. It grows as we approach the critical point g = 1, where correlations become long-ranged. We can see this by optimizing a second state right at criticality and comparing:

julia
H_crit = transverse_field_ising(; g = 1.0)
ψ_crit, = find_groundstate(ψ₀, H_crit, VUMPS(; verbosity = 0))
correlation_length(ψ_crit)
7.454602404710531

At a genuine critical point the correlation length diverges, but a finite bond dimension D can only capture correlations out to a finite range, so what we measure is large but capped rather than infinite.

4. Magnetization across the transition

As on the finite chain, we finish by sweeping the field g and recording the magnetization. The structure mirrors the finite sweep exactly — only InfiniteMPS and VUMPS have changed.

julia
g_values = 0.1:0.1:2
M = map(g_values) do g
    Hg = transverse_field_ising(; g = g)
    ψg, = find_groundstate(ψ₀, Hg, VUMPS(; verbosity = 0))
    return abs(expectation_value(ψg, 1 => σᶻ()))
end
scatter(g_values, M; xlabel = "g", ylabel = "M", label = "D = $D", title = "TFIM magnetization (L = ∞)")

Compare this with the finite-chain sweep of the previous tutorial, where the magnetization dropped to zero well before g = 1, at a point set by the algorithm rather than by the physics. The infinite curve instead tracks the transition itself: the magnetization stays on its ordered branch all the way up to the critical point and collapses to zero right at g = 1. What little smearing remains around the critical point is a finite-bond-dimension effect, and it shrinks as D grows.

We still take the absolute value of the magnetization, but for a subtly different reason than on the finite chain. On the finite chain the nonzero magnetization was an artifact of the algorithm: the exact ground state there is symmetric, and DMRG landed on a symmetry-broken state only because it carries less entanglement. In the thermodynamic limit the symmetry breaking is genuine — the two oppositely magnetized states become true ground states — and an infinite MPS at finite bond dimension settles into one of them on the ordered side, landing on a definite nonzero magnetization of either sign; abs again puts both branches onto a single order-parameter curve.

Where to go next

You have now run the same TFIM calculation twice — once at finite size, once directly at L = ∞ — and seen how little the code had to change.

From here you can go beyond ground states. A natural next step is to compute the excitations above this infinite ground state (the model's quasiparticle spectrum), or to exploit the symmetries of the model to make the calculation cheaper and more accurate.

To go deeper on the individual steps used here, see Constructing states, Controlling bond dimension, and Entanglement entropy and spectrum; the algorithm reference is Ground-state algorithms.