"""SU(2) singlet construction utilities and helper data structures.
This module is the user-facing surface for building local SU(2)
gauge-invariant (singlet) states. The flow when a caller asks for the
singlets of a chain of spin irreps ``[j_1, j_2, ..., j_N]`` is:
1. :func:`get_SU2_singlets` converts the spin list to doubled integers
and dispatches to the numba enumerator
:func:`._su2_singlet_enum._enumerate_singlet_paths_doubled`, which
returns every chain-of-couplings path that lands on a singlet
(``J_N = 0``), as flat ``(m_paths, J_paths, cg_products)`` arrays.
2. The output is sorted by the intermediate-J chain
``(J_2, J_3, ..., J_{N-1})``; rows sharing the same chain form one
multiplet. Each multiplet is wrapped into an :class:`SU2_singlet`
object: sympy-typed ``J_config`` / ``M_configs`` for downstream
compatibility, ``CG_values`` kept as floats.
3. :func:`SU2_singlet_canonical_vector` turns one :class:`SU2_singlet`
into a normalized vector in the local Hilbert basis. The per-leg
``(j, m) → index`` lookups built from :func:`get_spin_Hilbert_spaces`
are memoized by :func:`_cached_leg_lookup`, so repeated calls with
the same leg layout (as in the outer loop of
:func:`SU2_gauge_invariant_states`) skip the layout reconstruction.
Matter sector
-------------
Each matter site carries one fermionic doublet (red + green) per flavor.
A single matter leg has four states ``|n_r, n_g⟩``:
|0,0⟩ vacuum (n = 0, J = 0, m = 0)
|r,0⟩ (n = 1, J = 1/2, m = +1/2)
|0,g⟩ (n = 1, J = 1/2, m = -1/2)
|r,g⟩ pair (n = 2, J = 0, m = 0)
The vacuum and the pair share the same SU(2) labels, so we tag each
matter slot with its fermion occupation ``n ∈ {0, 1, 2}``. The
:class:`SU2_singlet` records this in its ``occupations`` field (one
entry per flavor); the per-leg lookup in :func:`_cached_leg_lookup`
uses ``(J, m, n)`` keys on matter legs (and plain ``(J, m)`` on gauge
and background legs).
For ``N_f > 1`` flavors, each flavor contributes one independent matter
leg to the chain. Flavors transform under the same local SU(2), so the
singlet-formation step automatically explores chains that combine
different flavors into local color singlets — the enumerator is
flavor-blind.
The Clebsch-Gordan kernel lives in :mod:`._clebsch_gordan` and the path
enumerator in :mod:`._su2_singlet_enum`; both work in pure integer
arithmetic on doubled spins / doubled magnetic numbers.
"""
import logging
from functools import lru_cache
import numpy as np
from sympy import S
from edlgt.tools import get_time, validate_parameters
from ._su2_singlet_enum import _enumerate_singlet_paths_doubled, spins_to_doubled
from .spin_operators import m_values, spin_space
logger = logging.getLogger(__name__)
__all__ = [
"SU2_singlet",
"get_SU2_singlets",
"SU2_singlet_canonical_vector",
"get_spin_Hilbert_spaces",
]
S0 = S(0)
S12 = S(1) / 2
# Conventional per-state fermion occupation of a matter leg. Keep in sync
# with the ``(j, m)`` layout in :func:`get_spin_Hilbert_spaces`.
_MATTER_LEG_OCCUPATIONS = (0, 1, 1, 2)
def _default_n_flavors(pure_theory, n_flavors):
"""Resolve the default ``n_flavors`` (``0`` for pure, ``1`` otherwise)."""
if n_flavors is not None:
return int(n_flavors)
return 0 if pure_theory else 1
# ----------------------------------------------------------------------------------
# SU2_singlet — central data type
# ----------------------------------------------------------------------------------
[docs]
class SU2_singlet:
"""Representation of one SU(2)-singlet multiplet in coupled-spin form."""
def __init__(
self,
J_config,
M_configs,
CG_values,
pure_theory=True,
occupations=None,
background=0,
):
"""Initialize an SU(2)-singlet descriptor.
Parameters
----------
J_config : list
Total-spin labels for the constituent degrees of freedom.
Matter slots carry the physical total spin (``S(0)`` for
vacuum / pair, ``S(1)/2`` for half-occupied).
M_configs : list
Allowed sets of magnetic quantum numbers producing the singlet.
CG_values : list
Lists of intermediate Clebsch-Gordan coefficients associated with
each ``M`` configuration.
pure_theory : bool, optional
If ``False``, matter degrees of freedom are included in the
leading entries of the configuration.
occupations : tuple of int, optional
Per-flavor fermion occupation ``n_f ∈ {0, 1, 2}`` (vacuum /
half / pair). Required when ``pure_theory=False`` and the
length determines the number of matter flavors.
background : int, optional
Background-charge sector information used for labeling.
Raises
------
TypeError
If inputs are not in the expected formats.
ValueError
If configuration lengths are inconsistent or ``occupations``
is missing / out of range for a matter theory.
"""
# CHECK ON TYPES
validate_parameters(spin_list=J_config, pure_theory=pure_theory)
n_spins = len(J_config)
if not isinstance(M_configs, list):
raise TypeError(f"M_configs must be a list of lists, not {type(M_configs)}")
if not isinstance(CG_values, list):
raise TypeError(f"CG_values must be a list of lists, not {type(CG_values)}")
if len(M_configs) != len(CG_values):
raise ValueError(f"M_configs and CG_values must have the same # of entries")
# Check each configuration of M values
for ii, conf in enumerate(M_configs):
validate_parameters(sz_list=conf)
if len(conf) != n_spins:
raise ValueError(f"{ii} M-config has {len(conf)} vals, not {n_spins}")
if len(CG_values[ii]) != (n_spins - 1):
n_CGs = len(CG_values[ii])
raise ValueError(f"{ii} M config has {n_CGs} CGs, not {n_spins-1}")
# Matter-sector occupations: one entry per flavor.
if pure_theory:
if occupations not in (None, ()):
raise ValueError(
"occupations must be empty / None when pure_theory=True"
)
occupations_tuple: tuple[int, ...] = ()
else:
if occupations is None:
raise ValueError(
"occupations is required when pure_theory=False (one entry per flavor)"
)
occupations_tuple = tuple(int(n) for n in occupations)
for n in occupations_tuple:
if n not in (0, 1, 2):
raise ValueError(f"occupation must be in {{0, 1, 2}}, got {n}")
# ----------------------------------------------------------------------------------
# Shallow copies suffice: J_config / M_configs entries are immutable
# sympy / python scalars; the singlet keeps its own list view.
self.n_spins = n_spins
self.J_config = list(J_config)
self.pure_theory = pure_theory
self.occupations = occupations_tuple
self.n_flavors = len(occupations_tuple)
# M-configs and the product of intermediate CGs for each path
self.M_configs = [list(m) for m in M_configs]
self.CG_values = [np.prod(CG_set) for CG_set in CG_values]
self.JM_configs = [self.J_config + m for m in self.M_configs]
[docs]
def display_singlet(self, msg: str | None = None):
"""Log the singlet components and their Clebsch-Gordan coefficients."""
if msg is not None:
msg_len = len(msg)
line = "=" * int(np.ceil((52 - 2 - msg_len) / 2))
logger.info(line + " " + msg + " " + line)
else:
logger.info("====================================================")
logger.info("J: %s", self.J_config)
if self.occupations:
logger.info("n: %s", list(self.occupations))
for m, CG in zip(self.M_configs, self.CG_values):
logger.info("M: %s CG:%s", m, float(CG))
# ----------------------------------------------------------------------------------
# Main entry points
# ----------------------------------------------------------------------------------
[docs]
@get_time
def get_SU2_singlets(
spin_list,
pure_theory=True,
occupations=None,
background=0,
):
"""Enumerate all SU(2) singlets compatible with a list of spin irreps.
Parameters
----------
spin_list : list
Spin irreps to be coupled. Matter slots (if any) come first and
must carry the physical spin associated with each flavor's
occupation (``S(0)`` for vacuum / pair, ``S(1)/2`` for half).
pure_theory : bool, optional
If ``True``, ``spin_list`` is treated as gauge-only and
``occupations`` must be empty / ``None``.
occupations : tuple of int, optional
Per-flavor fermion occupations ``n_f ∈ {0, 1, 2}``. Required when
``pure_theory=False``.
background : int, optional
Background-charge sector information.
Returns
-------
list or None
List of :class:`SU2_singlet` objects, or ``None`` if no singlet exists.
"""
validate_parameters(spin_list=spin_list, pure_theory=pure_theory)
n_spins = len(spin_list)
if n_spins < 2:
raise ValueError("2 is the minimum number of spins to form a singlet")
if not pure_theory and occupations is None:
raise ValueError(
"occupations is required when pure_theory=False (one entry per flavor)"
)
# ----------------------------------------------------------------------------------
# Hot path: int-encoded enumerator + sympy wrapper.
spins_d = spins_to_doubled(spin_list)
m_paths, J_paths, cg_products = _enumerate_singlet_paths_doubled(spins_d)
if m_paths.shape[0] == 0:
return None
# Sort rows by the intermediate-J chain (J_2d, J_3d, ..., J_Nd) so that
# paths belonging to the same multiplet land in a contiguous slice.
# ``np.lexsort`` treats the last key as primary, hence the reversed
# column order below.
sort_keys = [J_paths[:, k] for k in range(J_paths.shape[1] - 1, -1, -1)]
order = np.lexsort(sort_keys)
m_paths = m_paths[order]
J_paths = J_paths[order]
cg_products = cg_products[order]
# Multiplet boundaries: rows where the J chain changes.
boundaries = [0]
for i in range(1, m_paths.shape[0]):
if not np.array_equal(J_paths[i], J_paths[i - 1]):
boundaries.append(i)
boundaries.append(m_paths.shape[0])
# Build one SU2_singlet per multiplet. ``SU2_singlet`` expects
# ``CG_values`` as per-path lists of length ``n_spins - 1`` whose product
# gives the path's CG — we already have the product, so we pad with ones
# to match the shape contract.
n_pad = n_spins - 2
singlets = []
for grp_start, grp_end in zip(boundaries[:-1], boundaries[1:]):
M_configs = []
CG_values = []
for row in range(grp_start, grp_end):
m_config = [S(int(m_paths[row, k])) / 2 for k in range(n_spins)]
M_configs.append(m_config)
cg_list = [float(cg_products[row])]
if n_pad > 0:
cg_list.extend([1.0] * n_pad)
CG_values.append(cg_list)
singlets.append(
SU2_singlet(
J_config=spin_list,
M_configs=M_configs,
CG_values=CG_values,
pure_theory=pure_theory,
occupations=occupations,
background=background,
)
)
return singlets
[docs]
@get_time
def SU2_singlet_canonical_vector(spin_list, singlet, background=False):
"""Construct the canonical basis vector of a specific SU(2) singlet.
Parameters
----------
spin_list : list
Maximum spin irreps for each (gauge) degree of freedom.
singlet : SU2_singlet
Singlet descriptor returned by :func:`get_SU2_singlets`.
background : bool or int, optional
Background-charge flag/sector forwarded to
:func:`get_spin_Hilbert_spaces`.
Returns
-------
numpy.ndarray
Normalized canonical state vector for the requested singlet.
"""
validate_parameters(spin_list=spin_list)
if not isinstance(singlet, SU2_singlet):
raise TypeError(f"singlet is not {SU2_singlet}, not {type(singlet)}")
# Per-leg index dicts and strides for the local Hilbert space; cached.
leg_index_maps, strides, len_basis, matter_start, n_matter_legs = _cached_leg_lookup(
_spin_list_cache_key(spin_list),
singlet.pure_theory,
background,
singlet.n_flavors,
)
state = np.zeros(len_basis)
n_spins = singlet.n_spins
j_config = singlet.J_config
occupations = singlet.occupations
n_legs = len(strides)
matter_end = matter_start + n_matter_legs
# O(K) loop over the K paths contributing to this singlet: each path's
# (J_config, M_config) maps to a unique canonical-basis index. Matter
# legs additionally key on the per-flavor occupation so the vacuum and
# pair states (both J = 0, m = 0) are disambiguated.
for jm_config, cg in zip(singlet.JM_configs, singlet.CG_values):
m_config = jm_config[n_spins:]
idx = 0
for k in range(n_legs):
if matter_start <= k < matter_end:
n_f = occupations[k - matter_start]
idx += strides[k] * leg_index_maps[k][(j_config[k], m_config[k], n_f)]
else:
idx += strides[k] * leg_index_maps[k][(j_config[k], m_config[k])]
state[idx] = cg
# Sanity check: a properly built singlet vector is unit-normalized.
state_norm = float(np.sum(state**2))
if np.abs(state_norm - 1) > 1e-10:
raise ValueError(f"The state is not normalized: norm {state_norm}")
return state
# ----------------------------------------------------------------------------------
# Local Hilbert-space construction (used by both entry points and by callers
# that need the per-leg (j, m) lists directly).
# ----------------------------------------------------------------------------------
[docs]
@get_time
def get_spin_Hilbert_spaces(
max_spin_irrep_list,
pure_theory,
background=0,
n_flavors=None,
):
"""Build local spin Hilbert spaces used in singlet construction.
Parameters
----------
max_spin_irrep_list : list
Maximum spin irrep kept for each gauge degree of freedom.
pure_theory : bool
If ``False``, prepend the matter Hilbert space(s) used by the
SU(2) dressed-site construction.
background : int, optional
If nonzero, prepend the background-charge Hilbert space.
n_flavors : int, optional
Number of matter flavors. Each flavor contributes one independent
matter leg (4 fermion-occupation states). Defaults to ``1`` when
``pure_theory=False`` and ``0`` otherwise.
Returns
-------
tuple
``(j_list, m_list)`` with per-degree-of-freedom lists of spin irreps and
magnetic quantum numbers. Matter legs carry the physical
``J ∈ {0, 1/2, 1/2, 0}`` and ``m ∈ {0, +1/2, -1/2, 0}``; the
vacuum and pair states are distinguished by an occupation tag
carried separately on :class:`SU2_singlet`.
"""
validate_parameters(spin_list=max_spin_irrep_list, pure_theory=pure_theory)
n_flavors = _default_n_flavors(pure_theory, n_flavors)
spin_dof = []
m_list = []
j_list = []
# For each single spin particle in the list
for max_irrep in max_spin_irrep_list:
# Create an array with all the spin irreps up to the max one
spin_irreps = np.arange(S0, spin_space(max_irrep), 1) / 2
# For each spin dof save the list of allowed irreps
spin_dof.append(spin_irreps)
# run over each spin dof
for spin_irreps in spin_dof:
m_set = []
j_set = []
# run over the irreps of each spin dof
for irrep in spin_irreps:
# save the values of z-component of the irrep
m_set += list(m_values(irrep))
# save the irrep for each z-component
j_set += [irrep for i in m_values(irrep)]
# At the end, or each spin dof, len(m_set) = len(j_set) = direct sum of each irrep hilbert space
m_list.append(m_set)
j_list.append(j_set)
if not pure_theory and n_flavors > 0:
# Each flavor's matter leg has 4 fermion-occupation states. The
# (J, m) labels follow the (0, +1/2, -1/2, 0) layout; the vacuum
# and pair states have the same SU(2) numbers and are
# disambiguated by the per-flavor occupation tag stored on
# :class:`SU2_singlet`.
matter_j = [S0, S12, S12, S0]
matter_m = [S0, S12, -S12, S0]
for _ in range(n_flavors):
j_list.insert(0, list(matter_j))
m_list.insert(0, list(matter_m))
if background != 0:
m_set_bg = []
j_set_bg = []
for irrep in np.arange(S0, spin_space(background), 1) / 2:
# save the values of z-component of the irrep
m_set_bg += list(m_values(irrep))
# save the irrep for each z-component
j_set_bg += [irrep for i in m_values(irrep)]
# add the Hilbert space of the background charge
j_list.insert(0, j_set_bg)
m_list.insert(0, m_set_bg)
return j_list, m_list
# ----------------------------------------------------------------------------------
# Internal helpers for the canonical-vector lookup cache
# ----------------------------------------------------------------------------------
@lru_cache(maxsize=128)
def _cached_leg_lookup(spin_list_key, pure_theory, background, n_flavors):
"""Memoized leg layout for the canonical-vector construction.
Returns ``(leg_index_maps, strides, len_basis, matter_start, n_matter_legs)``.
The per-leg ``leg_index_maps[k]`` is a dict keyed by ``(j, m, n)`` on
matter legs and ``(j, m)`` on gauge / background legs. ``matter_start``
is the index of the first matter leg in the global leg list, and
``n_matter_legs`` the count of consecutive matter legs (zero for a
pure theory). Callers select the right key form based on whether the
leg index ``k`` falls in ``[matter_start, matter_start + n_matter_legs)``.
"""
# Rebuild a sympy spin_list from the doubled-int key.
spin_list = [S(j2) / S(2) if (j2 % 2) else S(j2 // 2) for j2 in spin_list_key]
j_list, m_list = get_spin_Hilbert_spaces(
spin_list, pure_theory, background, n_flavors=n_flavors
)
n_matter_legs = 0 if pure_theory else int(n_flavors)
matter_start = 1 if background != 0 else 0
leg_index_maps = []
leg_dims = []
for k, (j_leg, m_leg) in enumerate(zip(j_list, m_list)):
leg_dims.append(len(m_leg))
if matter_start <= k < matter_start + n_matter_legs:
leg_index_maps.append(
{
(j, m, n): i
for i, (j, m, n) in enumerate(
zip(j_leg, m_leg, _MATTER_LEG_OCCUPATIONS)
)
}
)
else:
leg_index_maps.append(
{(j, m): i for i, (j, m) in enumerate(zip(j_leg, m_leg))}
)
# Row-major strides matching product(*m_list): the last leg varies fastest.
n_legs = len(leg_dims)
strides = [1] * n_legs
for k in range(n_legs - 2, -1, -1):
strides[k] = strides[k + 1] * leg_dims[k + 1]
len_basis = strides[0] * leg_dims[0] if n_legs else 1
return tuple(leg_index_maps), tuple(strides), len_basis, matter_start, n_matter_legs
def _spin_list_cache_key(spin_list):
"""Hashable doubled-int tuple ``(2j_1, 2j_2, ...)`` for ``spin_list``."""
key = []
for s in spin_list:
v = 2.0 * float(s)
r = int(round(v))
if abs(v - r) > 1e-9:
raise ValueError(f"spin {s!r} is not a half-integer (2s = {v})")
key.append(r)
return tuple(key)