MPSKit.jl

Appearance

Statistical mechanics ​

This page collects recipes for the boundary-MPS (transfer-matrix) approach to two-dimensional classical statistical mechanics. The partition function of a classical lattice model is written as a contraction of a two-dimensional tensor network, one row of which is an InfiniteMPO — the transfer matrix. Contracting the network in the thermodynamic limit amounts to finding the leading eigenvector of that transfer matrix, which MPSKit approximates by an InfiniteMPS via leading_boundary.

The examples on this page use MPSKit.jl, MPSKitModels.jl, and TensorKit.jl. See Installation for how to add these packages to your environment.

julia
using MPSKit, MPSKitModels, TensorKit

For the structure of MPOs see Operators and Hamiltonians. For a full worked case using anyonic symmetries, see the gallery example The Hard Hexagon model.

1. Build the transfer matrix ​

MPSKitModels ships ready-made transfer-matrix MPOs for several classical models. The two-dimensional classical Ising model is provided by classical_ising; other options include sixvertex and hard_hexagon.

julia
mpo = classical_ising()
1-site InfiniteMPO(ComplexF64, TensorKit.ComplexSpace) with maximal dimension 2:
| ⋮
| ℂ^2
┼─[1]─ ℂ^2
│ ℂ^2
| ⋮

The returned object is an InfiniteMPO with a single-site unit cell: one tensor whose four legs are the two horizontal (virtual) and two vertical (physical) bonds of the Boltzmann-weight tensor. Its physical space is read off with physicalspace:

julia
P = physicalspace(mpo, 1)
ℂ^2

By default classical_ising uses the inverse temperature beta = log(1 + sqrt(2)) / 2.

2. Find the leading boundary MPS ​

leading_boundary approximates the dominant eigenvector of the transfer matrix by an InfiniteMPS. Supply an initial guess with a chosen bond dimension and an optimization algorithm — VUMPS is the usual choice.

The transfer matrix of a classical model is generally not Hermitian, so pass a non-Hermitian eigensolver to VUMPS:

julia
alg = VUMPS(;
    verbosity = 0,
    alg_eigsolve = MPSKit.Defaults.alg_eigsolve(; ishermitian = false),
)

ψ₀ = InfiniteMPS([P], [ℂ^16])          # initial guess, bond dimension 16
ψ, envs, ϵ = leading_boundary(ψ₀, mpo, alg)
ϵ                                       # final convergence error
9.57320793169616e-11

leading_boundary returns a triple (ψ, environments, ϵ): the converged boundary MPS, its environment manager, and the final convergence error. Reuse the returned envs in subsequent calls to avoid recomputing environments.

Note

The initial guess sets the bond dimension D of the boundary MPS. A larger D gives a better approximation of the leading eigenvector; near a critical point the accessible correlation length grows with D (see Controlling bond dimension for growing D on the fly).

3. Free energy and partition function per site ​

The expectation value of the transfer matrix in the converged boundary MPS is the partition function per site   in the thermodynamic limit:

julia
Λ = expectation_value(ψ, mpo)
2.5337371932371187 - 4.4375469636638844e-16im

For a Hermitian-normalised model the imaginary part is zero up to floating-point noise. The free energy per site follows from   :

julia
β = log(1 + sqrt(2)) / 2
f = -1 / β * log(real(Λ))
-2.109651067364864

4. Correlation length and entanglement entropy ​

The boundary MPS encodes the correlations of the two-dimensional system. correlation_length returns the (largest) correlation length of the transfer matrix, and entropy the entanglement entropy of the boundary MPS across a virtual bond:

julia
ξ = maximum(values(correlation_length(ψ)))
220.73426283874608
julia
S = real(first(entropy(ψ)))
0.4783482587753837

At a critical point the true correlation length diverges; the finite bond dimension of the boundary MPS cuts it off at a finite value that grows with D. This finite-entanglement scaling of S against log(ξ) is exactly what the gallery example The Hard Hexagon model exploits to extract a central charge.

5. Symmetric transfer matrices ​

When the classical model has a global symmetry, the transfer matrix can be built from symmetric tensors, which makes the boundary computation cheaper and more stable. classical_ising accepts a symmetry type; the spin-flip symmetry gives a Z2Irrep-graded MPO:

julia
mpo_z2 = classical_ising(Z2Irrep)
P_z2 = physicalspace(mpo_z2, 1)
Rep[ℤ₂](…) of dim 2:
 0 => 1
 1 => 1

The workflow is identical — only the virtual space of the initial guess is now a graded space:

julia
V_z2 = Z2Space(0 => 8, 1 => 8)
ψ_z2, = leading_boundary(InfiniteMPS([P_z2], [V_z2]), mpo_z2, alg)
real(expectation_value(ψ_z2, mpo_z2))
2.5337369771436005

For anyonic (non-invertible) symmetries the same recipe applies with a Vect[FibonacciAnyon] virtual space — this is the case worked out in The Hard Hexagon model. See Symmetries for how to choose graded spaces.

6. Multi-row unit cells ​

If the transfer matrix has a unit cell spanning several rows, the boundary object becomes a MultilineMPS and the operator a MultilineMPO. repeat stacks copies of a single-row MPO into a multi-row transfer matrix:

julia
mmpo = repeat(mpo, 2, 1)     # 2 rows × 1 column
MultilineMPO{InfiniteMPO{TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 2, Vector{ComplexF64}}}}(InfiniteMPO{TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 2, Vector{ComplexF64}}}[[TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 2, Vector{ComplexF64}}(ComplexF64[2.2071067811865452 + 0.0im, 0.49999999999999933 + 0.0im, 0.49999999999999933 + 0.0im, 0.20710678118654724 + 0.0im, 0.49999999999999933 + 0.0im, 0.2071067811865472 + 0.0im, 0.2071067811865472 + 0.0im, 0.49999999999999933 + 0.0im, 0.49999999999999933 + 0.0im, 0.2071067811865472 + 0.0im, 0.2071067811865472 + 0.0im, 0.4999999999999994 + 0.0im, 0.20710678118654724 + 0.0im, 0.4999999999999994 + 0.0im, 0.4999999999999994 + 0.0im, 2.2071067811865452 + 0.0im], (ℂ^2 ⊗ ℂ^2) ← (ℂ^2 ⊗ ℂ^2))], [TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 2, Vector{ComplexF64}}(ComplexF64[2.2071067811865452 + 0.0im, 0.49999999999999933 + 0.0im, 0.49999999999999933 + 0.0im, 0.20710678118654724 + 0.0im, 0.49999999999999933 + 0.0im, 0.2071067811865472 + 0.0im, 0.2071067811865472 + 0.0im, 0.49999999999999933 + 0.0im, 0.49999999999999933 + 0.0im, 0.2071067811865472 + 0.0im, 0.2071067811865472 + 0.0im, 0.4999999999999994 + 0.0im, 0.20710678118654724 + 0.0im, 0.4999999999999994 + 0.0im, 0.4999999999999994 + 0.0im, 2.2071067811865452 + 0.0im], (ℂ^2 ⊗ ℂ^2) ← (ℂ^2 ⊗ ℂ^2))]])

Build a MultilineMPS with one InfiniteMPS per row and call leading_boundary exactly as before; the returned environments are a MultilineEnvironments:

julia
mψ = MultilineMPS([InfiniteMPS([P], [ℂ^12]), InfiniteMPS([P], [ℂ^12])])
mψ, menvs, = leading_boundary(mψ, mmpo, alg)
real(expectation_value(mψ, mmpo))
6.419823378983467

The expectation_value of a multi-row transfer matrix is the product of the per-site weights over the rows of the unit cell; divide log of it by the number of rows to recover the per-site free energy.

See also ​