torch_vn_algebra.su2¶
Source: torch_vn_algebra/su2.py. Signatures and docstrings as in the source.
Module documentation¶
SU(2) symmetry of spin-1/2 chains: spin multiplets, highest-weight vectors and the Wigner-Eckart theorem.
Mathematics
-----------
The total spin S = sum_i S_i of L spins 1/2 (hbar = 1) acts on H = (C^2)^{(x) L}. By complete reducibility and
Schur's lemma (in the form of Schur-Weyl duality),
H = (+)_S C^{m_S} (x) V_S, dim V_S = 2S + 1, S = L/2, L/2 - 1, ..., >= 0,
with V_S the spin-S irreducible representation and m_S its multiplicity. An operator commutes with SU(2) if and
only if it acts as A_S (x) 1_{2S+1} on every summand, so the commutant of SU(2) is the Type I algebra
SU(2)' = (+)_S M_{m_S}(C) (block S of size m_S, one block per total spin),
and an SU(2)-invariant Hamiltonian, e.g. the Heisenberg model H = sum_<ij> J_ij S_i . S_j, is H = (+)_S H_S (x) 1:
the spectrum of H is the union of the spectra of the H_S, each level 2S+1 times degenerate. Compared with the
U(1) sectors N (number of up spins, dimension binom(L, N)) of `chains.SpinChain`, every block shrinks from
binom(L, L/2 + S) to m_S, and the 2S+1 copies of each multiplet are computed once.
Highest-weight vectors. The multiplicity space C^{m_S} is realised inside the U(1) sector N = L/2 + S
(S^z_total = N - L/2 = S) as the space of highest-weight vectors
H_S^hw = ker S^+ cap {S^z = S}, dim H_S^hw = m_S,
because a spin-S' multiplet has a state with S^z = S iff S' >= S, and that state is annihilated by S^+ iff
S' = S. With an orthonormal basis B_S (a binom(L, N) x m_S isometry) the block of an invariant operator A is
A_S = B_S^* A_N B_S, where A_N is the U(1) block of A in sector N: A_N maps H_S^hw to itself because A commutes
with S^+ and S^z. The map A -> (A_S)_S is a *-isomorphism of SU(2)' onto (+)_S M_{m_S}.
Multiplicities. Every multiplet with S' >= |M| contributes exactly one state to S^z = M, so the U(1) sector
dimensions are d_{L/2+M} = sum_{S' >= |M|} m_S'. Taking differences,
m_S = d_{L/2+S} - d_{L/2+S+1} = binom(L, L/2 - S) - binom(L, L/2 - S - 1), sum_S (2S+1) m_S = 2^L.
(Equivalently: S^+ maps sector N onto sector N + 1 for N >= L/2, and H_S^hw is its kernel.) m_S is also the
number of standard Young tableaux of the two-row shape (L/2 + S, L/2 - S).
Construction of B_S (Yamanouchi basis). Adding one site at a time, a spin-S multiplet of n sites arises from a
spin-(S -+ 1/2) multiplet |S', M'> of the first n - 1 sites coupled with the new spin 1/2 by Clebsch-Gordan
coefficients (Condon-Shortley phases, old (x) new):
S' = S - 1/2: |S, S> = |S', S'> (x) |up>
S' = S + 1/2: |S, S> = sqrt(2S'/(2S'+1)) |S', S'> (x) |down> - sqrt(1/(2S'+1)) |S', S'-1> (x) |up>,
|S', S'-1> = S^- |S', S'> / sqrt(2S').
Each basis vector is labelled by its path of intermediate spins (S_1, ..., S_L) = (1/2, ..., S), S_n >= 0,
|S_n - S_{n-1}| = 1/2 (a Yamanouchi symbol / standard Young tableau). Different paths give orthogonal vectors and
each step is an isometry, so B_S is orthonormal by construction; no null space, SVD or Gram-Schmidt is needed and
the cost is O(L * sum_S m_S binom(L, L/2+S)) for all S together (sparse lowering operators only). The blocks
B_S^* H B_S then cost O(m_S^2 binom(L, L/2+S)) dense flops. The basis of the n-site sector N is the
concatenation [(sector N-1 of n-1 sites) (x) up; (sector N of n-1 sites) (x) down], which is exactly the order of
the iterated fused product used by SpinChain (labels in decreasing order), so B_S is expressed directly in the
U(1) basis of `SpinChain(L)`.
Wigner-Eckart theorem. The site spins form a rank-1 tensor operator with spherical components
T_{+1} = -S^+_i / sqrt 2, T_0 = S^z_i, T_{-1} = S^-_i / sqrt 2. For multiplets |S, M, a> (a = Yamanouchi label),
<S', M', a | T_q | S, M, b> = <S M; 1 q | S' M'> <S', a || S_i || S, b> / sqrt(2S' + 1)
(Edmonds' / Racah's normalisation of the reduced matrix element, for which
sum_{M, M', q} |<S' M' a|T_q|S M b>|^2 = |<S' a||S_i||S b>|^2; Sakurai divides by sqrt(2S+1) instead). Only
S' - S in {-1, 0, 1} occur. The reduced elements are computed from highest-weight vectors only, choosing M = S,
q = S' - S, M' = S'. Scalar products of two rank-1 operators then follow without any M sum on the U(1) level:
<S a| S_i . S_j |S b> = sum_{S''} c(S, S'') sum_c <S a||S_i||S'' c> <S'' c||S_j||S b>,
c(S, S'') = (-1)^{S - S''} / (2S + 1),
with S_i . S_j = sum_q (-1)^q T^{(i)}_q T^{(j)}_{-q}; `spin_correlation` evaluates c(S, S'') from Clebsch-Gordan
coefficients.
Conventions
-----------
* Basis and sectors as in `chains.SpinChain`: bit i of a label is site i, 1 = up, S^z_i = bit - 1/2, sector N =
number of up spins, S^z_total = N - L/2; the highest-weight vectors of spin S live in sector N = L/2 + S.
* Sectors of the algebra are ordered by increasing S (S_min = 0 or 1/2, ..., L/2) and carry the charges 2S
(integers, exact for half-integer spins); `spins` gives the S values. Functions accept spins as int, float
(multiples of 1/2) or fractions.Fraction.
* H = sum_bonds J_ij S_i . S_j + h sum_i S^z_i, J > 0 antiferromagnetic (the convention of SpinChain.xxz with
Delta = 1). A uniform field commutes with S^2 and S^z but splits the multiplets: block S of `heisenberg` is H
restricted to the highest-weight vectors (M = S), i.e. H_S^(h=0) + h S 1, and the multiplet levels are
E + h (M - S), M = -S..S (`full_spectrum` adds them).
References: E. Lieb and D. Mattis, J. Math. Phys. 3, 749 (1962) (ordering of energy levels); A. R. Edmonds,
Angular Momentum in Quantum Mechanics (Princeton, 1957), ch. 3-5; J. J. Sakurai, Modern Quantum Mechanics, sec.
3.8-3.11; G. Racah, Phys. Rev. 62, 438 (1942) (closed formula for the Clebsch-Gordan coefficients);
H. Weyl / I. Schur (Schur-Weyl duality); T. Yamanouchi, Proc. Phys.-Math. Soc. Japan 19, 436 (1937).
chain = SU2Chain(L=12, boundary='periodic')
chain.multiplicities # [132, 297, 275, 154, 54, 11, 1] for S = 0..6
H = chain.heisenberg(J=1.0, J2=0.2) # Operator of chain.algebra = (+)_S M_{m_S}
E = chain.full_spectrum(H) # all 2^L levels, sorted
gs = chain.ground_state(H) # energy, spin S, S^z and highest-weight coefficients
psi = chain.embed(gs.spin[0], gs.coefficients[0, :chain.multiplicity(gs.spin[0])]) # SpinChain vector
Limits: the bases of all spins together hold sum_S m_S binom(L, L/2+S) numbers in float64 (L = 14: 7.2e6,
55 MB; L = 16: 1.0e8, 0.75 GB; L = 18: 1.4e9, 10.6 GB), and the largest block has m_S = 1001 (L = 14), 3640
(L = 16), 13260 (L = 18) rows, at S = 1 or 2, against binom(L, L/2) = 3432, 12870, 48620 for U(1). The
padded algebra of all spins holds C m_max^2 numbers per operator (L = 16: 1.2e8, 0.9 GB real, 1.8 GB complex);
`heisenberg(..., spin=S)` builds one block in `spin_algebra(S)` without padding, and the bases of a single spin
are built by a pruned recursion (L = 18, S = 0: 4862 x 48620, 1.9 GB). Dense methods are practical up to
L = 16 on a workstation (L = 14: about 1 s for the full spectrum on one CPU core, against 10 s for the dense
U(1) sectors); `cost.check_memory` guards the large allocations.
clebsch_gordan_exact¶
def clebsch_gordan_exact(j1: Spin, m1: Spin, j2: Spin, m2: Spin, J: Spin, M: Spin) -> Tuple[int, Fraction]
Clebsch-Gordan coefficient <j1 m1; j2 m2 | J M> in exact arithmetic, returned as (sign, square) with sign in
{-1, 0, 1} and square a Fraction, so that the coefficient is sign * sqrt(square). Racah's formula with the
Condon-Shortley phase convention (<j1 j1; j2 (J - j1) | J J> > 0):
<j1 m1 j2 m2|J M> = delta_{M, m1+m2} sqrt[(2J+1) (J+j1-j2)! (J-j1+j2)! (j1+j2-J)! / (j1+j2+J+1)!]
* sqrt[(J+M)! (J-M)! (j1-m1)! (j1+m1)! (j2-m2)! (j2+m2)!]
* sum_k (-1)^k / [k! (j1+j2-J-k)! (j1-m1-k)! (j2+m2-k)! (J-j2+m1+k)! (J-j1-m2+k)!].
Arguments are integers or half-integers (int, float or Fraction); coefficients violating the selection rules
(triangle, |m| <= j, j - m integer, M = m1 + m2) are 0.
clebsch_gordan¶
def clebsch_gordan(j1: Spin, m1: Spin, j2: Spin, m2: Spin, J: Spin, M: Spin) -> float
<j1 m1; j2 m2 | J M> as a float64 (exact rational arithmetic, one rounding at the end); see
`clebsch_gordan_exact`.
wigner_eckart¶
def wigner_eckart(reduced, j_out: Spin, m_out: Spin, k: Spin, q: Spin, j_in: Spin, m_in: Spin)
Matrix element <j_out m_out | T^k_q | j_in m_in> = <j_in m_in; k q | j_out m_out> <j_out||T^k||j_in> /
sqrt(2 j_out + 1) from a reduced matrix element (number or tensor) in Edmonds' normalisation.
multiplicity¶
def multiplicity(L: int, S: Spin) -> int
m_S = binom(L, L/2 - S) - binom(L, L/2 - S - 1): the number of spin-S multiplets of L spins 1/2
(0 if S is not one of L/2, L/2 - 1, ..., >= 0).
two_spins¶
def two_spins(L: int) -> List[int]
The values 2S of the total spins of L spins 1/2, increasing: L mod 2, L mod 2 + 2, ..., L.
SU2GroundState¶
class SU2GroundState()
SU2Chain¶
class SU2Chain(L: int, boundary: str='open', complex_valued: bool=True, precision: str='double', device=None)
Spin-1/2 chain of L sites with its SU(2) symmetry: the algebra (+)_S M_{m_S} of SU(2)-invariant operators,
orthonormal highest-weight bases B_S in the U(1) sectors of `SpinChain`, invariant Heisenberg Hamiltonians
(nearest and next-nearest neighbours or arbitrary bonds), the full spectrum with its degeneracies, reduced
matrix elements of the site spins (Wigner-Eckart) and spin correlations.
The U(1) chain with the same basis conventions is `self.u1` (a SpinChain); `embed`, `multiplet_basis`,
`restrict` and `lift` translate between the two descriptions.
SU2Chain.algebra¶
@property
def algebra(self) -> TypeIAlgebra
(+)_S M_{m_S}, sectors ordered by increasing S, charges 2S.
SU2Chain.multiplicity¶
def multiplicity(self, S: Spin) -> int
SU2Chain.sector_index¶
def sector_index(self, S: Spin) -> int
Position of spin S among the sectors of `algebra`.
SU2Chain.u1_sector¶
def u1_sector(self, S: Spin, M: Optional[Spin]=None) -> int
U(1) sector N = L/2 + M (default M = S) holding the states |S, M, a>.
SU2Chain.bonds¶
@property
def bonds(self) -> List[Tuple[int, int]]
Nearest-neighbour bonds, as SpinChain.bonds (periodic L = 2: one bond).
SU2Chain.next_nearest_bonds¶
@property
def next_nearest_bonds(self) -> List[Tuple[int, int]]
Distinct pairs {i, i+2}: open chains i = 0..L-3; periodic chains i + 2 mod L, every unordered pair once
(L = 4: the two diagonals once each; L = 3: the pairs coincide with the nearest-neighbour bonds).
SU2Chain.build_bases¶
def build_bases(self, spins: Optional[Sequence[Spin]]=None) -> 'SU2Chain'
Build (and cache) the highest-weight bases of the given spins (default: all) in one recursion; otherwise
they are built on first use. Returns self.
SU2Chain.lower¶
def lower(self, X: torch.Tensor, N: int) -> torch.Tensor
S^-_total on vectors of the U(1) sector N of `u1`: X of shape (binom(L, N),) or (binom(L, N), r).
SU2Chain.raise_¶
def raise_(self, X: torch.Tensor, N: int) -> torch.Tensor
S^+_total on vectors of the U(1) sector N of `u1` (sector N -> N + 1).
SU2Chain.hw_basis¶
def hw_basis(self, S: Spin, dtype: Optional[torch.dtype]=None) -> torch.Tensor
Orthonormal basis of the highest-weight vectors of spin S, as the columns of a real (binom(L, N), m_S)
matrix in the basis of the U(1) sector N = L/2 + S of `u1` (Yamanouchi order, see `paths`). Computed once
on first use (float64) and cached; dtype: cast (cached), default float64. The cached tensor is returned:
clone it before modifying it in place.
SU2Chain.paths¶
def paths(self, S: Spin) -> List[Tuple[int, ...]]
Yamanouchi labels of the columns of hw_basis(S): (2 S_1, ..., 2 S_L) with S_n the spin of sites 0..n-1.
SU2Chain.multiplet_basis¶
def multiplet_basis(self, S: Spin, M: Spin) -> torch.Tensor
The states |S, M, a> (a = 1..m_S) as columns in the U(1) sector N = L/2 + M of `u1`:
(S^-)^{S-M} B_S, normalised by prod sqrt((S + M')(S - M' + 1)) (Condon-Shortley); float64 (for M = S the
cached hw_basis(S) itself).
SU2Chain.embed¶
def embed(self, S: Spin, coefficients: torch.Tensor, M: Optional[Spin]=None) -> torch.Tensor
The vector sum_a c_a |S, M, a> (default M = S, the highest weight) in the U(1) sector N = L/2 + M of `u1`,
for coefficients of shape (m_S,) or (batch, m_S) (e.g. eigenvectors of block S); dtype of `u1`.
SU2Chain.multiplet¶
def multiplet(self, S: Spin, coefficients: torch.Tensor) -> List[Tuple[int, torch.Tensor]]
All 2S+1 states of the multiplet with highest-weight coordinates `coefficients`: [(N, vector)] for
M = S, S-1, ..., -S (vectors in the U(1) sector N of `u1`).
SU2Chain.spin_algebra¶
def spin_algebra(self, S: Spin) -> TypeIAlgebra
The single block M_{m_S} (no padding, charge 2S), cached; for large L, where the padded algebra of all
spins is mostly padding.
SU2Chain.heisenberg¶
def heisenberg(self, J=1.0, h=0.0, J2=0.0, bonds: Optional[Sequence[Tuple[int, int]]]=None, spin: Optional[Spin]=None) -> Operator
SU(2)-invariant Heisenberg Hamiltonian in the algebra (+)_S M_{m_S}:
H = J sum_<ij> S_i . S_j + J2 sum_<<ij>> S_i . S_j + h sum_i S^z_i
(nearest-neighbour bonds `bonds`, next-nearest `next_nearest_bonds`; J2 != 0 gives the frustrated J1-J2
chain). With an explicit `bonds` list [(i, j), ...] the exchange part is sum_b J_b S_i . S_j with J a number,
a tensor (n_bonds,) or (batch, n_bonds) (one Hamiltonian per row, e.g. random couplings), and J2 must be 0.
J, J2, h may be numbers or tensors (batch,). Block S is H restricted to the highest-weight vectors,
H_S^(h=0) + h S 1 (the field is recorded as `su2_field` on the result for `full_spectrum`).
spin=None: operator of `algebra` (all spins); spin=S: the block of spin S only, an operator of
`spin_algebra(S)`. Cost: building the bases once, then O(m_S^2 binom(L, L/2+S)) per term and spin.
SU2Chain.restrict¶
def restrict(self, A: Operator) -> Operator
The SU(2)-invariant operator A of the full U(1) algebra `u1.algebra` (e.g. SpinChain.xxz with Delta = 1)
as an element of (+)_S M_{m_S}: block S = B_S^* A_N B_S with N = L/2 + S (A restricted to the highest-weight
vectors). For an A that is not SU(2)-invariant this is a compression, not a homomorphism; for an A that
commutes with S^2 and S^z (e.g. with a uniform field) it is the restriction to S^z = S.
SU2Chain.lift¶
def lift(self, A: Operator) -> Operator
The inverse of `restrict` on SU(2)' : the element (+)_S A_S of (+)_S M_{m_S} as the operator
(+)_S A_S (x) 1_{2S+1} of the U(1) algebra `u1.algebra`, block N = sum_{S >= |M|} X_{S,M} A_S X_{S,M}^*
with X_{S,M} = multiplet_basis(S, M). Memory: one dense block per U(1) sector (sum_N binom(L, N)^2).
SU2Chain.casimir¶
def casimir(self) -> Operator
S^2 = (+)_S S(S+1) 1 (central).
SU2Chain.full_spectrum¶
def full_spectrum(self, H: Operator, h=None, return_quantum_numbers: bool=False)
All 2^L energies of H (an operator of `algebra` with blocks H restricted to highest-weight vectors, e.g.
from `heisenberg`), sorted, shape (batch, 2^L): every eigenvalue E of block S gives the 2S+1 levels
E + h (M - S), M = -S..S. h: the uniform field (number or (batch,)); default the field recorded by
`heisenberg`, otherwise 0 (pass it explicitly for operators derived by arithmetic).
return_quantum_numbers: also return S and M of every level (tensors of the same shape).
SU2Chain.lowest_energies¶
def lowest_energies(self, H: Operator, h=None) -> torch.Tensor
(batch, C): the lowest energy among the states of total spin S, for every S (with a field h the
multiplet member M = -S for h > 0).
SU2Chain.ground_state¶
def ground_state(self, H: Operator, h=None) -> SU2GroundState
Lowest level of H (an operator of `algebra`): energy, total spin S (the smallest S on exact ties), S^z,
and the eigenvector of block S in the highest-weight basis, zero padded to k_max (pass
`coefficients[b, :m_S]` to `embed(S, ..., M=sz)` for the vector in the U(1) sector).
SU2Chain.expectation¶
def expectation(self, A: Operator, S: Spin, coefficients: torch.Tensor) -> torch.Tensor
<psi|A|psi> for psi = sum_a c_a |S, M, a> (any M: A is SU(2)-invariant), c of shape (m_S,) or
(batch, m_S); A an operator of `algebra` or of `spin_algebra(S)`. Returns (batch,).
SU2Chain.reduced_matrix_elements¶
def reduced_matrix_elements(self, sites: Optional[Sequence[int]]=None) -> Dict[Tuple[int, int], torch.Tensor]
Reduced matrix elements <S', a || S_i || S, b> (Edmonds' normalisation, module docstring) of the site
spins, from highest-weight vectors only: <S' S' a | T_q | S S b> with q = S' - S, divided by
<S S; 1 q | S' S'> / sqrt(2S'+1). Returns {(2S', 2S): tensor (n_sites, m_S', m_S)} for all S' - S in
{-1, 0, 1} (zero blocks for S = S' = 0), real dtype of the chain. sites: default all, in order.
Hermiticity: R_{S,S'} = (-1)^{S'-S} R_{S',S}^T.
SU2Chain.spin_correlation¶
def spin_correlation(self, pairs: Sequence[Tuple[int, int]], method: str='wigner_eckart') -> Operator
The SU(2)-invariant operators S_i . S_j for pairs [(i, j), ...] as a batch of operators of `algebra`
(batch index = pair). method='wigner_eckart': from the reduced matrix elements of the site spins,
block S = sum_{S''} c(S, S'') R^i_{S,S''} R^j_{S'',S}; method='direct': B_S^T (S_i . S_j) B_S with the
sparse bond operator. Expectation values in a multiplet: `expectation(C, S, coefficients)`.