Computing observables

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

This page collects recipes for extracting physical quantities from an MPS: local and multi-site expectation values, the energy of a Hamiltonian, two-point correlators, and the energy variance as a convergence diagnostic. All examples share a single namespace and build on state and operator objects you would have in hand after a ground-state calculation.

using MPSKit, TensorKit
using TensorKitTensors.SpinOperators: σˣ, σᶻ

For building MPS objects see Constructing states. For controlling the bond dimension during optimization see Controlling bond dimension. The reference page for ground-state algorithms is Ground-state algorithms.

---

Setup: state and operators

The examples below use a spin-1/2 FiniteMPS together with the Pauli operators from TensorKitTensors.jl. These are ComplexF64 TensorMaps, matching the default element type of the state.

L = 8
ψ = FiniteMPS(L, ℂ^2, ℂ^8)   # random finite MPS, bond dim ≤ 8

# single-site Pauli operators
X = σˣ()
Z = σᶻ()

2←2 TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}:
 codomain: ⊗(ℂ^2)
 domain: ⊗(ℂ^2)
 blocks: 
 * Trivial() => 2×2 reshape(view(::Vector{ComplexF64}, 1:4), 2, 2) with eltype ComplexF64:
 1.0+0.0im   0.0+0.0im
 0.0+0.0im  -1.0+0.0im

The finite TFIM Hamiltonian used in recipes 3 and 5 is built from these:

lattice = fill(ℂ^2, L)
H = FiniteMPOHamiltonian(lattice, (i, i + 1) => -(X ⊗ X) for i in 1:(L - 1)) +
    FiniteMPOHamiltonian(lattice, (i,) => -0.5 * Z for i in 1:L)

8-site FiniteMPOHamiltonian(ComplexF64, TensorKit.ComplexSpace) with maximal dimension 3:
┬─[8]─ ℂ^2
│ (ℂ^1 ⊞ ℂ^1 ⊞ ℂ^1)
┼─[7]─ ℂ^2
│ (ℂ^1 ⊞ ℂ^1 ⊞ ℂ^1)
┼─[6]─ ℂ^2
│ (ℂ^1 ⊞ ℂ^1 ⊞ ℂ^1)
┼─[5]─ ℂ^2
│ (ℂ^1 ⊞ ℂ^1 ⊞ ℂ^1)
┼─[4]─ ℂ^2
│ (ℂ^1 ⊞ ℂ^1 ⊞ ℂ^1)
┼─[3]─ ℂ^2
│ (ℂ^1 ⊞ ℂ^1 ⊞ ℂ^1)
┼─[2]─ ℂ^2
│ (ℂ^1 ⊞ ℂ^1 ⊞ ℂ^1)
┴─[1]─ ℂ^2

---

1. Local (one-site) expectation value

Use expectation_value(ψ, i => O) to evaluate ⟨ψ|Oᵢ|ψ⟩ at a single site i. The pair i => O identifies the site and the single-site operator.

expectation_value(ψ, 4 => Z)   # ⟨Z⟩ at site 4

0.02669478621326907 + 1.50038454832517e-18im

To compute a local observable at every site, broadcast over the indices:

[expectation_value(ψ, i => Z) for i in 1:L]

8-element Vector{ComplexF64}:
   -0.2698358582460263 - 1.5166165747195903e-17im
   -0.0658453463520652 - 3.700847778157179e-18im
 -0.008512944774245414 - 4.784713650815622e-19im
   0.02669478621326907 + 1.50038454832517e-18im
   0.09317294153598002 + 5.2367994516135655e-18im
   0.26396763439492227 + 1.483634132672786e-17im
  -0.10298039572249759 - 5.788028916509736e-18im
   0.24558878599380093 + 1.380335533700128e-17im

Note

The state ψ must be normalised for the expectation value to be meaningful. A freshly constructed FiniteMPS is normalised by default; if you modified the tensors by hand, call normalize!(ψ) first.

---

2. Multi-site (contiguous) expectation value

For a product of operators on a contiguous range of sites, pass a tuple of indices together with a multi-site operator formed by taking tensor products ⊗:

# ⟨X₂ X₃⟩ — two-site operator on sites 2 and 3
expectation_value(ψ, (2, 3) => X ⊗ X)

0.9526406959270952 - 1.9678349875681903e-18im

The operator X ⊗ X is a {2,2} TensorMap (two incoming, two outgoing legs) matching the two-site index tuple (2, 3). The tuple must be contiguous; arbitrary non-adjacent index sets are not supported by this form.

# ⟨Z₁ Z₂ Z₃⟩ — three-site operator
expectation_value(ψ, (1, 2, 3) => Z ⊗ Z ⊗ Z)

0.00904802640992798 + 2.243269186357736e-18im

---

3. Energy (full-MPO expectation value)

When the operator is an AbstractMPO (e.g. a Hamiltonian), pass it directly without an index argument. MPSKit evaluates the full contraction ⟨ψ|H|ψ⟩:

E = expectation_value(ψ, H)

-6.715460830564508 - 4.00259619932691e-16im

The result is a scalar; for a Hermitian H and a normalised ψ its imaginary part is zero up to floating-point noise.

The same form works for InfiniteMPS with an InfiniteMPOHamiltonian, where the returned value is the energy per unit cell.

Note

The full-MPO form automatically computes and caches the environments. If you already have environments from a prior find_groundstate call you can pass them as a trailing argument to avoid recomputation, but this is optional; omitting them is always safe and correct.

---

4. Two-point correlators

correlator computes ⟨O₁ᵢ O₂ⱼ⟩ for two sites with i < j. The recommended call uses a single two-site operator O₁₂:

# ⟨Z₂ Zⱼ⟩ for a single target site j = 6
correlator(ψ, Z ⊗ Z, 2, 6)

-0.01742418741370863 + 1.3877787807814457e-17im

Warning

i must be strictly less than j. Calling correlator(ψ, O₁₂, i, j) with i ≥ j will throw an error.

Correlation profile over a range

Pass a range as j to obtain a vector of correlators — one entry per target site. This is the efficient route for a full correlation profile:

# ⟨Z₂ Zⱼ⟩ for j = 3, 4, …, L
corr = correlator(ψ, Z ⊗ Z, 2, 3:L)

6-element Vector{ComplexF64}:
  -0.025076420488843384 - 1.0408340855860843e-17im
 -0.0010918344995656793 - 1.0842021724855044e-17im
  -0.005992104279236707 + 8.239936510889834e-18im
   -0.01742418741370863 + 1.3877787807814457e-17im
  0.0068138574049169285 + 1.3010426069826053e-18im
  -0.016169918384384016 + 5.204170427930421e-18im

The result is a Vector whose k-th element corresponds to j = 3 + k - 1.

A common pattern is to normalise the correlator by ⟨Z⟩² to extract the connected part:

z_mean = expectation_value(ψ, 2 => Z)
connected = [c - z_mean * expectation_value(ψ, j => Z) for (j, c) in zip(3:L, corr)]

6-element Vector{ComplexF64}:
  -0.025636958285979578 - 1.0471351081367524e-17im
  0.0006658929444413575 - 1.0644435044363528e-17im
 0.00014290032684061843 + 8.929574258225434e-18im
  -4.314710123965587e-5 + 1.5831595874326166e-17im
  3.3077581096341464e-5 + 5.388130695759009e-19im
   9.602895610257012e-7 + 7.021943853901375e-18im

This subtracts the disconnected part to leave the connected correlator  .

---

5. Energy variance as a convergence check

variance returns ⟨H²⟩ − ⟨H⟩², which is zero if and only if ψ is an exact eigenstate of H. Use it as a quantitative convergence diagnostic after a ground-state search:

var_E = variance(ψ, H)

2.7977337579228774

A smaller variance indicates that ψ is closer to a true eigenstate.

After running a ground-state algorithm the variance should have dropped significantly compared to the random starting state above:

ψ_gs, envs, _ = find_groundstate(ψ, H, DMRG(; maxiter = 10))
variance(ψ_gs, H)

1.5120349416974932e-11

Note

The variance function also accepts an optional pre-computed envs argument. Pass the environments returned by find_groundstate to skip recomputation:

variance(ψ_gs, H, envs)