Skip to content

Heisenberg, U(1)

The ground state of the spin-1/2 Heisenberg chain on \(N = 20\) sites, open boundaries:

\[ H = \sum_{i=0}^{N-2}\mathbf{S}_i\cdot\mathbf{S}_{i+1} = \sum_{i=0}^{N-2}\left[S^z_iS^z_{i+1} + \tfrac12\left(S^+_iS^-_{i+1} + S^-_iS^+_{i+1}\right)\right], \qquad E_0 = \min_{\psi}\frac{\langle\psi\vert H\vert\psi\rangle}{\langle\psi\vert\psi\rangle}. \]

Representation. spin_half() grades the physical doublet by U(1) with charge \(2S^z\), so \(S^\pm\) are rank-3 tensors carrying \(\mp 2\) on a D=1 charge leg. tenet.models.heisenberg(N) writes the three summands above as three MPO.from_terms entries, coefficient for coefficient, and lets the builder derive the graded MPO bond from the operators' own charges; Heisenberg, U(1) walkthrough spells that term list out. The state is an MPS, \(\Psi_{s_1\cdots s_N} = A^{s_1}\cdots A^{s_N}\), seeded as the Néel product state, whose D=1 bond-0 leg carries charge \(0\) — that is \(S^z_{\mathrm{tot}} = 0\), structurally.

Approximation. dmrg_ sweeps two-site updates, each a Lanczos ground eigenpair of the two-site effective Hamiltonian followed by a truncated SVD back onto the chosen bond. The schedule ramps \(\chi\) — three sweeps at 16 with noise \(10^{-4}\), three at 32 with \(10^{-5}\), then 64 noiseless — because noise repopulates bond sectors a D=1 product seed left empty, which a rescaling update can never reach on its own. The energy is variational, so it approaches \(E_0\) from above.

Checks.

  • \(E\) matches the recorded \(N = 20\) exact-diagonalization energy to \(10^{-10}\).
  • The bond energies \(\langle\mathbf{S}_i\cdot\mathbf{S}_{i+1}\rangle\), measured with a single invariant \(\mathbf{S}\cdot\mathbf{S}\) operator rather than the three-term sum the MPO was built from, sum to out.energy — the measurement and the sweep agree.
  • \(\langle S^z_n\rangle\) is float noise on every site: the sector is fixed by the boundary leg and the tensors' invariance, not by a penalty term.

Step by step in the DMRG tutorial; the semantics of every call in DMRG.

Source

"""U(1) Heisenberg ground state through ``tenet.network``, end to end.

Run it standalone::

    uv run python examples/heisenberg.py

What ``examples/toy_codes/dmrg.py`` writes out by hand -- the 5x5 ``W`` matrix, its
channel table, the reachable bond spaces -- this file never mentions: the Hamiltonian is
one call to ``tenet.models.heisenberg``, which writes the term list and lets
``MPO.from_terms`` derive the graded MPO bond, and the Neel product state seeds the
``S^z_tot = 0`` sector by its own charges. ``examples/heisenberg_walkthrough.py``
is the middle of the three: the same library calls, with the ``W`` and the bond spaces
spelled out beside them. The schedule below ramps ``chi``
with noise because writing one is what a user does; on this chain it buys nothing -- the
flat ``chi=64`` run reaches the same energy, since a degeneracy-1 U(1) seed already
ramps for free (see docs/tutorials/dmrg.md, "Schedules and noise", for when it pays).
"""

from tenet.models import heisenberg, spin_half
from tenet.network import MPS, Sweep, dmrg_, expectation_1site, expectation_2site, local_op
from tenet.symmetry import U1Sector

# The standard spin-1/2 site, graded by U(1): charge t = 2 S^z, so the doublet is
# {-1, +1}. Sz/S+/S- arrive as local_op's rank-3 charge-leg operators and S.S as the
# matrix of the invariant two-site term; no spin matrix is written out in this file.
SITE = spin_half()
PHYS = SITE.phys
SZ = SITE.matrices["Sz"]


def main(n_sites: int = 20, chi: int = 64):
    """Ground state at the defaults CI runs; returns the DMRG_out and the bond profile."""
    # Neel: up, down, up, ... Its charges sum to zero, and DMRG never leaves the sector
    # its seed is in, so this one product state is the whole S^z_tot = 0 input.
    psi = MPS.product(PHYS, [U1Sector(1 if n % 2 else -1) for n in range(n_sites)])
    # Staged chi: early sweeps at a small bond are cheap and move the state most, and
    # they hand the later, expensive sweeps a starting point already near the minimum.
    # Noise repopulates bond charges the product seed left empty -- without it a sector
    # that starts at zero weight can never be reached, since the update only rescales
    # what is already there. It decays because by then the missing sectors are found
    # and further noise would only be energy the last, noiseless sweep has to undo.
    schedule = [Sweep(16, noise=1e-4)] * 3 + [Sweep(32, noise=1e-5)] * 3 + [Sweep(chi)]
    # H = sum_i S_i . S_{i+1}, J = 1, open boundaries. Under U(1) that term list is
    # S^z S^z with the transverse half split into S^+ S^- and S^- S^+, since raising and
    # lowering are separate operators here; examples/heisenberg_walkthrough.py writes it
    # out, and this file asks the model function for it.
    out = dmrg_(psi, heisenberg(n_sites), schedule=schedule)
    # legs[0] of site n is its left virtual bond, so this is the cut through the middle
    # of the chain -- the one that carries the most entanglement and the largest bond.
    mid = out.psi[n_sites // 2].legs[0].space
    print(f"N={n_sites}  ~{out.sweeps} sweeps  E = {out.energy:.15f}  mid bond: {mid.dim} states")

    # S.S as a single two-site matrix, so the bond energy is one expectation value
    # rather than the three-term sum the MPO was built from.
    ss = local_op(SITE.matrices["S.S"], phys=PHYS)
    profile = [expectation_2site(out.psi, ss, n) for n in range(n_sites - 1)]
    print("bond energies:", " ".join(f"{e:+.4f}" for e in profile))
    # H is the sum of those bonds, so the profile summing to out.energy checks the
    # variational energy against a route that never touches the environment caches.
    print(f"sum of bond energies = {sum(profile):.15f}  vs  out.energy = {out.energy:.15f}")

    # <S^z_n> = 0 site by site is not convergence but symmetry: a state living in one
    # U(1) sector has no local magnetisation to round-off, however far from the minimum.
    op_sz = local_op(SZ, phys=PHYS)
    max_sz = max(abs(expectation_1site(out.psi, op_sz, n)) for n in range(n_sites))
    print(f"max_n |<S^z_n>| = {max_sz:.1e}")
    return out, profile


if __name__ == "__main__":
    main()

Output

Produced by heisenberg.main() at its defaults — exactly python examples/heisenberg.py — as run by tests/test_examples.py.

N=20  ~9 sweeps  E = -8.682473334397722  mid bond: 64 states
bond energies: -0.6534 -0.2943 -0.5664 -0.3370 -0.5401 -0.3540 -0.5286 -0.3616 -0.5239 -0.3638 -0.5239 -0.3616 -0.5286 -0.3540 -0.5401 -0.3370 -0.5664 -0.2943 -0.6534
sum of bond energies = -8.682473334397697  vs  out.energy = -8.682473334397722
max_n |<S^z_n>| = 3.9e-13