"""U(1) lattice gauge-theory (QED) model helpers."""
import logging
import time
import numpy as np
from numba import typed
from scipy.sparse import identity
from edlgt.entanglement import flux_block_projector
from edlgt.modeling import (
LocalTerm,
NBodyTerm,
PlaquetteTerm,
TwoBodyTerm,
check_link_symmetry,
enumerate_block_classes,
get_origin_surfaces,
get_site_label,
outward_link_ops,
staggered_mask,
tile_plaquette_blocks,
zig_zag,
)
from edlgt.operators import (
QED_dressed_site_operators,
QED_gauge_integrated_operators,
QED_gauge_invariant_operators,
QED_gauge_invariant_states,
QED_plq_supersite_operators,
)
from .quantum_model import (
QuantumModel,
format_loc_dims,
validate_staggered_matter_lattice,
)
logger = logging.getLogger(__name__)
__all__ = ["QED_Model"]
[docs]
class QED_Model(QuantumModel):
"""QED lattice gauge model with optional matter fields."""
def __init__(
self,
spin,
pure_theory,
bg_list=None,
plaq_basis=False,
link_symmetries=True,
get_only_bulk=False,
block_projectors=None,
**kwargs,
):
"""Initialize the QED model and construct its symmetry sector.
Parameters
----------
spin : float or str = integrated
Gauge-link spin representation.
pure_theory : bool
If True, build the pure-gauge theory (no matter fields).
bg_list : list, optional
Optional background-charge configuration used during local-basis
projection.
get_only_bulk : bool, optional
Restrict gauge-invariant local states to bulk-compatible ones when
supported by the operator factory.
block_projectors : dict, optional
Plaquette-block-LBO projectors, required when ``plaq_basis=True``.
``lvals`` is then the coarse super-lattice (one 2x2 block per site).
Keyed by position class: ``<class>`` is the dense ``(D, k)`` block
projector and ``<class>_embedding`` the sparse ``(prod_dim, D)``
subsystem embedding, both from
:func:`extract_plaquette_block_projectors`.
**kwargs
Arguments forwarded to :class:`~edlgt.models.quantum_model.QuantumModel`.
"""
# Initialize base class with the common parameters
super().__init__(**kwargs)
if not pure_theory:
validate_staggered_matter_lattice(self.lvals, "QED", self.has_obc)
# Build local operator dictionaries in complex mode, then optionally
# switch to real mode at Hamiltonian-assembly time when allowed.
self.configure_dtype_mode(dtype_mode="complex", auto_mode="complex")
# Check Background charge list
if bg_list is None:
self.bg_list = None
else:
static_charges = np.asarray(bg_list, dtype=float)
if static_charges.shape != (self.n_sites,):
msg = f"len(bg_list) must match n_sites = {self.n_sites}: got {static_charges.shape}"
raise ValueError(msg)
self.bg_list = static_charges.tolist()
# Check spin truncation: it can be integer or "integrated"
if isinstance(spin, str):
if spin != "integrated":
raise ValueError(f"spin must be integer or 'integrated': got {spin}")
else:
if not np.all([pure_theory == False, self.dim == 1]):
msg = "QED with integrated gauge fields is only in 1D with matter"
raise ValueError(msg)
self.spin = spin
self.pure_theory = pure_theory
self.plaq_basis = False
self.background = 0
self.staggered_basis = False
logger.info("----------------------------------------------------")
logger.info("(%s+1)D QED MODEL with gauge fields integrated", self.dim)
if self.bg_list is not None:
logger.info("bg_list=%s", self.bg_list)
ops = QED_gauge_integrated_operators()
self.project_operators(ops)
elif np.isscalar(spin):
self.spin = spin
self.pure_theory = pure_theory
self.plaq_basis = plaq_basis
self.link_symmetries = link_symmetries
self.staggered_basis = False if self.pure_theory else True
pure_label = "pure" if self.pure_theory else "with matter"
self.background = 0
if self.bg_list is not None:
self.background = int(max(np.abs(self.bg_list)))
if self.background != 0:
logger.info("bg_list=%s", self.bg_list)
else:
self.bg_list = None
if not self.plaq_basis:
logger.info("----------------------------------------------------")
msg = f"({self.dim}+1)D QED MODEL {pure_label} j={spin}"
logger.info(msg)
# -------------------------------------------------------------------------------
# Acquire gauge invariant basis and states
self.gauge_basis, self.gauge_states = QED_gauge_invariant_states(
self.spin,
self.pure_theory,
lattice_dim=self.dim,
background=self.background,
get_only_bulk=get_only_bulk,
)
# -------------------------------------------------------------------------------
# Acquire operators
ops = QED_dressed_site_operators(
self.spin, self.pure_theory, self.dim, background=self.background
)
# Initialize the operators, local dimension and lattice labels
self.project_operators(ops, self.bg_list)
else:
logger.info("----------------------------------------------------")
logger.info(
"(%s+1)D QED PLAQUETTE MODEL %s j=%s", self.dim, pure_label, spin
)
# PLAQUETTE-BLOCK LBO BASIS
# self.lvals is the coarse super-lattice; each site is a 2x2 block
# of the original 2*self.lvals dressed lattice, reduced by the
# block_projectors. Fills self.ops / self.loc_dims and stores the
# boundary-link generators for the symmetry sector below.
self.block_projectors = block_projectors
self._setup_plaquette_basis()
else:
raise ValueError(f"spin must be integer or 'integrated': got {spin}")
# -------------------------------------------------------------------------------
# GLOBAL SYMMETRIES
if self.pure_theory:
global_ops = None
global_sectors = None
else:
global_ops = [self.ops["N"]]
global_sectors = [int(self.n_sites / 2)]
# -------------------------------------------------------------------------------
# LINK SYMMETRIES (only present in the truncated version of gauge fields)
if self.spin == "integrated":
link_ops = None
link_sectors = None
else:
if not self.plaq_basis:
if self.link_symmetries:
dirs = self.directions
link_ops = [
[self.ops[f"E_p{d}"], -self.ops[f"E_m{d}"]] for d in dirs
]
link_sectors = [0 for _ in self.directions]
else:
link_ops = None
link_sectors = None
else:
# Inter-super-site link continuity (set by _setup_plaquette_basis):
# one combined generator per direction packing both shared rungs.
link_ops = self._plaq_link_ops
link_sectors = self._plaq_link_sectors
# ---------------------------------------------------------------------------
# ELECTRIC-FLUX “NBODY” SYMMETRIES
# only in the pure (no-matter) theory, more than 1D, *and* PBC
# Constrain, for each cartesian direction, the corresponding
# Electric flux on the face/line through the origin:
# 3D: (Ex → yz-face at x=0) (Ey → xz-face at y=0) (Ez → xy-face at z=0)
# 2D: (Ex → y-axis at x=0) (Ey → x-axis at y=0)
if self.pure_theory and not any(self.has_obc):
logger.info("fixing surface electric fluxes")
# one flux‐constraint per cartesian direction
nbody_sectors = np.zeros(self.dim, dtype=float)
nbody_ops = []
nbody_sites_list = typed.List()
surfaces = get_origin_surfaces(self.lvals)
nbody_sym_type = "U"
if self.dim == 2:
# in 2D we have two lines through (0,0):
line_of = {"x": "y", "y": "x"}
for direction in self.directions:
sites = np.array(surfaces[line_of[direction]][1], dtype=np.uint8)
nbody_sites_list.append(sites)
nbody_ops.append(self.ops[f"E_p{direction}"])
elif self.dim == 3:
# in 3D we have three faces through (0,0,0):
face_of = {"x": "yz", "y": "xz", "z": "xy"}
for direction in self.directions:
sites = np.array(surfaces[face_of[direction]][1], dtype=np.uint8)
nbody_sites_list.append(sites)
logger.debug(
"%s sites: %s %s",
direction,
sites,
surfaces[face_of[direction]][0],
)
nbody_ops.append(self.ops[f"E_p{direction}"])
else:
# no electric‐flux constraint in 1D, or in OBC or with matter
nbody_sectors = None
nbody_ops = None
nbody_sites_list = None
nbody_sym_type = None
# GET SYMMETRY SECTOR
self.get_abelian_symmetry_sector(
global_ops=global_ops,
global_sectors=global_sectors,
link_ops=link_ops,
link_sectors=link_sectors,
nbody_ops=nbody_ops,
nbody_sectors=nbody_sectors,
nbody_sites_list=nbody_sites_list,
nbody_sym_type=nbody_sym_type,
)
# DEFAULT PARAMS
self.default_params()
def _get_auto_dtype_mode(self, m=None, theta=0.0):
"""Heuristic dtype mode for QED Hamiltonian assembly."""
if self.plaq_basis:
# LBO projectors (RDM eigenvectors) make the super-site ops complex in general.
return "complex"
if not self.pure_theory:
return "complex"
if np.abs(theta) <= 1e-12:
return "real"
return "complex"
[docs]
def build_Hamiltonian(self, g, m=None, theta=0.0, dtype_mode="auto"):
"""Dispatch to the appropriate QED Hamiltonian builder.
Parameters
----------
dtype_mode : str or bool, optional
"auto", "real", "complex", or legacy bool flag.
"""
resolved_mode = self.configure_dtype_mode(
dtype_mode=dtype_mode,
auto_mode=self._get_auto_dtype_mode(m=m, theta=theta),
)
if not self.plaq_basis:
if self.spin == "integrated":
return self.build_integrated_Hamiltonian(
g, m, theta, dtype_mode=resolved_mode
)
else:
return self.build_truncated_Hamiltonian(
g, m, theta, dtype_mode=resolved_mode
)
else:
return self.build_plaquette_Hamiltonian(g, dtype_mode=resolved_mode)
[docs]
def build_truncated_Hamiltonian(self, g, m=None, theta=0.0, dtype_mode="auto"):
"""Assemble the QED Hamiltonian.
Parameters
----------
g : float
Gauge coupling.
m : float, optional
Bare fermion mass (used only when matter is present).
theta : float, optional
Topological-angle coupling parameter.
dtype_mode : str or bool, optional
"auto", "real", "complex", or legacy bool flag.
"""
self.configure_dtype_mode(
dtype_mode=dtype_mode,
auto_mode=self._get_auto_dtype_mode(m=m, theta=theta),
)
logger.info("----------------------------------------------------")
logger.info("BUILDING truncated HAMILTONIAN")
# Hamiltonian Coefficients
self.QED_Hamiltonian_couplings(g, m, theta)
h_terms = {}
# -------------------------------------------------------------------------------
# ELECTRIC ENERGY
op_name = "E2"
h_terms["E2"] = LocalTerm(self.ops[op_name], op_name, **self.def_params)
self.H.add_term(h_terms["E2"].get_Hamiltonian(strength=self.coeffs["E"]))
# -------------------------------------------------------------------------------
# PLAQUETTE TERM: MAGNETIC INTERACTION
if self.dim > 1:
op_names_list = ["C_px,py", "C_py,mx", "C_my,px", "C_mx,my"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["plaq_xy"] = PlaquetteTerm(
["x", "y"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["plaq_xy"].get_Hamiltonian(
strength=self.coeffs["B"], add_dagger=True
)
)
if self.dim == 3:
# XZ Plane
op_names_list = ["C_px,pz", "C_pz,mx", "C_mz,px", "C_mx,mz"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["plaq_xz"] = PlaquetteTerm(
["x", "z"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["plaq_xz"].get_Hamiltonian(
strength=self.coeffs["B"], add_dagger=True
)
)
# YZ Plane
op_names_list = ["C_py,pz", "C_pz,my", "C_mz,py", "C_my,mz"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["plaq_yz"] = PlaquetteTerm(
["y", "z"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["plaq_yz"].get_Hamiltonian(
strength=self.coeffs["B"], add_dagger=True
)
)
# -------------------------------------------------------------------------------
# TOPOLOGICAL TERM
# -------------------------------------------------------------------------------
if self.dim == 2 and np.abs(self.coeffs["theta"]) > 10e-10:
logger.info("Adding topological term")
op_names_list = ["C_px,py", "C_py,mx", "C_my,px", "C_mx,my"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["plaq_xy"] = PlaquetteTerm(
["x", "y"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["plaq_xy"].get_Hamiltonian(
strength=self.coeffs["theta"], add_dagger=True
)
)
if self.dim == 3 and np.abs(self.coeffs["theta"]) > 10e-10:
logger.info("Adding topological term")
# XY Plane
op_names_list = ["EzC_px,py", "C_py,mx", "C_my,px", "C_mx,my"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["Ez_Bxy"] = PlaquetteTerm(
["x", "y"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["Ez_Bxy"].get_Hamiltonian(
strength=self.coeffs["theta"], add_dagger=True
)
)
# XZ Plane
op_names_list = ["EyC_px,pz", "C_pz,mx", "C_mz,px", "C_mx,mz"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["Ey_Bxz"] = PlaquetteTerm(
["x", "z"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["Ey_Bxz"].get_Hamiltonian(
strength=self.coeffs["theta"], add_dagger=True
)
)
# YZ Plane
op_names_list = ["ExC_py,pz", "C_pz,my", "C_mz,py", "C_my,mz"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["Ex_Byz"] = PlaquetteTerm(
["y", "z"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["Ex_Byz"].get_Hamiltonian(
strength=self.coeffs["theta"], add_dagger=True
)
)
# -------------------------------------------------------------------------------
if not self.pure_theory:
# ---------------------------------------------------------------------------
# STAGGERED MASS TERM
for stag_label in ["even", "odd"]:
h_terms[f"N_{stag_label}"] = LocalTerm(
operator=self.ops["N"], op_name="N", **self.def_params
)
self.H.add_term(
h_terms[f"N_{stag_label}"].get_Hamiltonian(
self.coeffs[f"m_{stag_label}"],
staggered_mask(self.lvals, stag_label),
)
)
# ---------------------------------------------------------------------------
# HOPPING
for d in self.directions:
for stag_label in ["even", "odd"]:
# Define the list of the 2 non trivial operators
op_names_list = [f"Q_p{d}_dag", f"Q_m{d}"]
op_list = [self.ops[op] for op in op_names_list]
# Define the Hamiltonian term
h_terms[f"{d}_hop_{stag_label}"] = TwoBodyTerm(
d, op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms[f"{d}_hop_{stag_label}"].get_Hamiltonian(
strength=self.coeffs[f"t{d}_{stag_label}"],
add_dagger=True,
mask=staggered_mask(self.lvals, stag_label),
)
)
# ---------------------------------------------------------------------------
# APPLY LINK PENALTIES IN CASE LINK SYMMETRIES ARE NOT APPLIED
if not self.link_symmetries:
self.H.add_term(h_terms["E2"].get_Hamiltonian(strength=100))
for d in self.directions:
# Define the list of the 2 non trivial operators
op_names_list = [f"E_p{d}", f"E_m{d}"]
op_list = [self.ops[op] for op in op_names_list]
# Define the Hamiltonian term
h_terms["penalty"] = TwoBodyTerm(
d, op_list, op_names_list, **self.def_params
)
self.H.add_term(h_terms["penalty"].get_Hamiltonian(strength=-100))
self.H.build(format=self.ham_format)
[docs]
def build_integrated_Hamiltonian(self, g, m, theta=0.0, dtype_mode="auto"):
"""Assemble the integrated-gauge 1D QED Hamiltonian.
Parameters
----------
dtype_mode : str or bool, optional
"auto", "real", "complex", or legacy bool flag.
"""
self.configure_dtype_mode(dtype_mode=dtype_mode, auto_mode="complex")
# Hamiltonian Coefficients
self.QED_Hamiltonian_couplings(g, m, theta)
logger.info("BUILDING integrated HAMILTONIAN")
if self.dim != 1 or self.pure_theory:
msg = "Integrated QED Hamiltonian is defined only in 1D with matter."
raise ValueError(msg)
h_terms = {}
n_links = self.n_sites - 1
# ---------------------------------------------------------------------------
# STAGGERED MASS TERM
h_terms["N"] = LocalTerm(operator=self.ops["N"], op_name="N", **self.def_params)
for _, stag_label in enumerate(["even", "odd"]):
mask = staggered_mask(self.lvals, stag_label)
mass_coeff = self.coeffs[f"m_{stag_label}"]
self.H.add_term(h_terms["N"].get_Hamiltonian(mass_coeff, mask))
# ---------------------------------------------------------------------------
# HOPPING TERM
op_names_list = ["Q_dag", "Q"]
op_list = [self.ops[op] for op in op_names_list]
# Define the Hamiltonian term
h_terms["hop"] = TwoBodyTerm("x", op_list, op_names_list, **self.def_params)
self.H.add_term(
h_terms["hop"].get_Hamiltonian(strength=self.coeffs["t"], add_dagger=True)
)
# ---------------------------------------------------------------------------
# ELECTRIC ENERGY: long-range Hamiltonian
# Eq. structure:
# H_E = g^2/2 * sum_{n=0}^{L-2}[sum_{k=0}^{n}(q_k + (-Sz_k + (-1)^k)/2)]^2
# = sum_i h_i * Sz_i + sum_{i<j} J_ij * Sz_i * Sz_j + const
static_charges = self._integrated_static_charges()
stagger_sign = np.ones(self.n_sites, dtype=float)
for ii in range(self.n_sites):
if ii % 2 == 1:
stagger_sign[ii] = -1.0
charge_offsets = static_charges + 0.5 * stagger_sign
# A[n] constant factors
prefix_charge = np.cumsum(charge_offsets)
# ---------------------------------------------------------------------------
# One-body coefficients h_i = -E * sum_{n=i}^{L-2} prefix_charge[n]
linear_sz_coeff = np.zeros(self.n_sites, dtype=float)
for kk in range(n_links):
linear_sz_coeff[kk] = -self.coeffs["E"] * np.sum(prefix_charge[kk:n_links])
# Build and cache single-site masks once
site_masks = [self._single_site_mask(ii) for ii in range(self.n_sites)]
# One-body Sz part
h_terms["Sz"] = LocalTerm(
operator=self.ops["Sz"], op_name="Sz", **self.def_params
)
for ii in range(self.n_sites):
coeff_ii = linear_sz_coeff[ii]
if np.abs(coeff_ii) < 1e-14:
continue
self.H.add_term(
h_terms["Sz"].get_Hamiltonian(strength=coeff_ii, mask=site_masks[ii])
)
# ---------------------------------------------------------------------------
# Two-body Sz_i Sz_j part with J_ij = g^2/2 * 0.5 * (L-1-j), i<j
szsz_terms = {}
op_list = [self.ops["Sz"], self.ops["Sz"]]
op_names_list = ["Sz", "Sz"]
for dist in range(1, self.n_sites):
szsz_terms[dist] = NBodyTerm(
op_list, op_names_list, distances=[(dist,)], **self.def_params
)
for jj in range(1, n_links):
right_weight = 0.5 * self.coeffs["E"] * (self.n_sites - 1 - jj)
if np.abs(right_weight) > 1e-14:
for ii in range(jj):
dist = jj - ii
self.H.add_term(
szsz_terms[dist].get_Hamiltonian(
strength=right_weight, mask=site_masks[ii]
)
)
# ---------------------------------------------------------------------------
# Constant electric energy term:
# E * sum_n [ A_n^2 + (1/4) * sum_{k<=n} (Sz_k)^2 ], with A_n as prefix_charge
# and (Sz_k)^2 = I for Pauli operators.
if n_links > 0:
sum_sz2_counts = 0.25 * 0.5 * n_links * (n_links + 1)
electric_constant = np.sum(prefix_charge[:n_links] ** 2) + sum_sz2_counts
electric_constant *= self.coeffs["E"]
else:
electric_constant = 0.0
if np.abs(electric_constant) > 1e-14:
hdim = self.H.shape[0]
self.H.add_term(float(electric_constant) * identity(hdim, format="csc"))
self.H.build(format=self.ham_format)
def _integrated_static_charges(self) -> np.ndarray:
"""Return static charges q_n for integrated-QED sectors."""
if self.bg_list is None:
return np.zeros(self.n_sites, dtype=float)
static_charges = np.asarray(self.bg_list, dtype=float)
if static_charges.shape != (self.n_sites,):
logger.info("Integrated QED static charges must have one value per site")
msg = f"expected shape ({self.n_sites},), got {static_charges.shape}."
raise ValueError(msg)
return static_charges
[docs]
def reconstruct_integrated_E2_from_N(
self,
density_obs_name: str = "N",
density_corr_obs_name: str = "N_N",
state_index: int | None = None,
dynamics: bool = False,
compute_density_corr: bool = True,
print_values: bool = True,
) -> np.ndarray:
"""Reconstruct link-resolved <E^2> in integrated 1D QED from matter density.
Parameters
----------
density_obs_name : str, optional
Key in self.res containing measured site-resolved <N_k>.
density_corr_obs_name : str, optional
Key in self.res containing measured two-point correlator
<N_k N_l>.
state_index : int or None, optional
Eigenstate index used to compute <N_k N_l> on the fly when
density_corr_obs_name is missing and compute_density_corr=True.
If None and only one eigenstate is available, index 0 is used.
dynamics : bool, optional
If True, interpret state_index as a time index and compute
missing correlators on self.H.psi_time[state_index] instead of
self.H.Npsi[state_index].
compute_density_corr : bool, optional
If True, compute <N_k N_l> when not already present in
self.res. If False, missing correlations raise KeyError.
print_values : bool, optional
If True, print link-resolved reconstructed values and their
average using the same style as local observable measurements.
Returns
-------
numpy.ndarray
Link-resolved reconstructed Casimir values <E_{k,k+1}^2> with
shape (n_sites - 1,).
Notes
-----
The reconstruction uses
E_n = sum_{k=0}^n [ q_k + N_k + ((-1)^k-1)/2 ].
Therefore:
<E_n^2> = B_n^2 + 2 B_n sum_{k<=n}<N_k> + sum_{k,l<=n}<N_k N_l>,
where B_n = sum_{k=0}^n [q_k + ((-1)^k-1)/2].
"""
if self.spin != "integrated" or self.dim != 1 or self.pure_theory:
msg = "Integrated E2 reconstruction is defined only for 1D QED with matter"
raise ValueError(msg)
if density_obs_name not in self.res:
msg = f"Missing '{density_obs_name}' in self.res. Measure local density first."
raise KeyError(msg)
density_values = np.asarray(self.res[density_obs_name], dtype=float)
if density_values.shape != (self.n_sites,):
msg = f"Expected {density_obs_name}.shape=({self.n_sites},) got {density_values.shape}."
raise ValueError(msg)
n_links = self.n_sites - 1
if n_links <= 0:
link_casimir = np.zeros(0, dtype=float)
self.res["E2"] = link_casimir
self.res["E2_avg"] = 0.0
return link_casimir
# Acquire or compute density correlator <N_k N_l>.
#
# In dynamics, or when a specific eigenstate/time-step index is
# requested, reusing a previously cached correlator would silently mix
# data from different states. Recompute in those cases and overwrite the
# cache with the current state's correlator.
reuse_cached_density_corr = (
density_corr_obs_name in self.res
and not compute_density_corr
and not dynamics
and state_index is None
)
if reuse_cached_density_corr:
density_corr = np.asarray(self.res[density_corr_obs_name], dtype=float)
else:
if not compute_density_corr and density_corr_obs_name not in self.res:
raise KeyError(
f"Missing '{density_corr_obs_name}' in self.res. "
"Measure ['N','N_N'] or set compute_density_corr=True."
)
if not compute_density_corr:
density_corr = np.asarray(self.res[density_corr_obs_name], dtype=float)
else:
if not hasattr(self, "H"):
raise ValueError(
"Cannot compute density correlator: Hamiltonian not available."
)
if dynamics:
if not hasattr(self.H, "psi_time"):
raise ValueError(
"Cannot compute density correlator in dynamics mode: "
"time-evolved states (H.psi_time) are not available."
)
if state_index is None:
if len(self.H.psi_time) == 1:
state_index = 0
else:
raise ValueError(
"state_index is required to compute density correlator "
"when multiple time steps are available."
)
target_state = self.H.psi_time[state_index]
else:
if not hasattr(self.H, "Npsi"):
raise ValueError(
"Cannot compute density correlator: Hamiltonian eigenstates not available."
)
if state_index is None:
if len(self.H.Npsi) == 1:
state_index = 0
else:
raise ValueError(
"state_index is required to compute density correlator "
"when multiple eigenstates are available."
)
target_state = self.H.Npsi[state_index]
op_list = [self.ops["N"], self.ops["N"]]
op_names_list = ["N", "N"]
density_corr_term = TwoBodyTerm(
"x", op_list=op_list, op_names_list=op_names_list, **self.def_params
)
density_corr_term.get_expval(target_state)
density_corr = np.asarray(density_corr_term.corr, dtype=float)
self.res[density_corr_obs_name] = density_corr
if density_corr.shape != (self.n_sites, self.n_sites):
raise ValueError(
f"Expected '{density_corr_obs_name}' shape ({self.n_sites}, {self.n_sites}), "
f"got {density_corr.shape}."
)
density_corr = np.array(density_corr, dtype=float, copy=True)
# For fermions N_k^2 = N_k, so force the diagonal explicitly.
np.fill_diagonal(density_corr, density_values)
static_charges = self._integrated_static_charges()
site_indices = np.arange(self.n_sites, dtype=float)
staggered_offsets = 0.5 * ((-1.0) ** site_indices - 1.0)
charge_offsets = static_charges + staggered_offsets
prefix_charge = np.cumsum(charge_offsets)
prefix_density = np.cumsum(density_values)
# Prefix sums of <N_k N_l>: top-left block sum at link n is [n,n].
density_corr_prefix = density_corr.cumsum(axis=0).cumsum(axis=1)
prefix_density_corr = np.diag(density_corr_prefix)
link_casimir = (
prefix_charge[:n_links] ** 2
+ 2.0 * prefix_charge[:n_links] * prefix_density[:n_links]
+ prefix_density_corr[:n_links]
)
self.res["E2"] = link_casimir
self.res["E2_avg"] = float(np.mean(link_casimir))
if print_values:
logger.info("----------------------------------------------------")
logger.info("E2 ")
for ii in range(n_links):
logger.info("%s %.16f", (ii,), link_casimir[ii])
logger.info("%.16f", self.res["E2_avg"])
return link_casimir
[docs]
def get_fermionic_string_correlator(
self,
state=None,
state_index=None,
dynamics: bool = False,
print_values: bool = False,
):
"""Measure the gauge-invariant fermionic string correlator matrix.
Notes
-----
For the integrated 1D QED model, the off-diagonal entries are built as
<Q_dag(i) U(i+1) ... U(j-1) Q(j)> for i < j,
with Q_dag = Sm, Q = Sp, and U = P_psi = Sz.
For the truncated dressed-site model, the same gauge-invariant string
is represented by <Q_px_dag(i) U(i+1) ... U(j-1) Q_mx(j)>.
In both cases, the diagonal is the local density <N(i)>.
"""
if self.dim != 1 or self.pure_theory:
msg = "Fermionic string correlator is implemented only for 1D QED with matter."
raise NotImplementedError(msg)
if self.spin == "integrated":
required_ops = {"U", "Q_dag", "Q", "N"}
left_op_name = "Q_dag"
right_op_name = "Q"
else:
required_ops = {"U", "Q_px_dag", "Q_mx", "N"}
left_op_name = "Q_px_dag"
right_op_name = "Q_mx"
string_op_name = "U"
if not required_ops.issubset(self.ops):
raise NotImplementedError(
"Fermionic string correlator requires the appropriate 1D QED "
"operators for the selected representation."
)
if not all(self.has_obc):
msg = "Fermionic string correlator is currently implemented only with OBC."
raise NotImplementedError(msg)
if not hasattr(self, "def_params"):
self.default_params()
target_state = self._get_measurement_state(
state=state, state_index=state_index, dynamics=dynamics
)
corr = np.zeros((self.n_sites, self.n_sites), dtype=np.complex128)
density_term = LocalTerm(self.ops["N"], "N", **self.def_params)
density_term.get_expval(target_state, print_values=False)
np.fill_diagonal(corr, density_term.obs.astype(np.complex128))
for dist in range(1, self.n_sites):
op_list = [self.ops[left_op_name]]
op_list += [self.ops[string_op_name] for _ in range(dist - 1)]
op_list += [self.ops[right_op_name]]
op_names_list = [left_op_name]
op_names_list += [string_op_name for _ in range(dist - 1)]
op_names_list += [right_op_name]
distances = [(step,) for step in range(1, dist + 1)]
corr_term = NBodyTerm(
op_list=op_list,
op_names_list=op_names_list,
distances=distances,
**self.def_params,
)
corr_term.get_expval(target_state, component="real", print_values=False)
real_values = np.array(corr_term.obs, dtype=float, copy=True)
corr_term.get_expval(target_state, component="imag", print_values=False)
imag_values = np.array(corr_term.obs, dtype=float, copy=True)
string_values = real_values + 1j * imag_values
for ii, value in enumerate(string_values):
jj = ii + dist
corr[ii, jj] = value
corr[jj, ii] = np.conjugate(value)
self.res["fermionic_string_correlator"] = corr
if print_values:
logger.info("----------------------------------------------------")
logger.info("FERMIONIC STRING CORRELATOR")
for ii in range(self.n_sites):
msg = np.array2string(corr[ii], precision=8, suppress_small=False)
logger.info(msg)
return corr
[docs]
def measure_fermionic_nongaussianity(
self,
state=None,
state_index=None,
dynamics: bool = False,
print_value: bool = True,
eig_tol: float = 1e-10,
):
"""Measure the fermionic non-Gaussianity from the string correlator."""
corr = self.get_fermionic_string_correlator(
state=state,
state_index=state_index,
dynamics=dynamics,
print_values=False,
)
herm_corr = 0.5 * (corr + corr.conj().T)
eigvals = np.linalg.eigvalsh(herm_corr)
eigvals = np.real_if_close(eigvals, tol=1000).astype(float)
min_eval = float(np.min(eigvals))
max_eval = float(np.max(eigvals))
if min_eval < -eig_tol or max_eval > 1 + eig_tol:
logger.warning(
"Fermionic string correlator has eigenvalues outside [0, 1]: "
"min=%.6e, max=%.6e. Clipping for NG.",
min_eval,
max_eval,
)
clipped = np.clip(eigvals, 0.0, 1.0)
entropy_vals = np.zeros_like(clipped)
valid = np.logical_and(clipped > eig_tol, clipped < 1.0 - eig_tol)
vv = clipped[valid]
entropy_vals[valid] = -(vv * np.log(vv) + (1.0 - vv) * np.log(1.0 - vv))
ng_value = float(np.sum(entropy_vals))
self.res["fermionic_string_correlator_eigvals"] = eigvals
self.res["fermionic_nongaussianity"] = ng_value
if print_value:
logger.info("----------------------------------------------------")
logger.info("FERMIONIC NON-GAUSSIANITY")
logger.info("NG = %.16f", ng_value)
return ng_value, eigvals, corr
[docs]
def boundary_flux_signatures(
self, keep_sites, block_coords, unique_subsys_configs, axes=("x", "y")
):
"""Compute the boundary electric-flux signature of every block RDM row.
Parameters
----------
keep_sites : list[int]
The four block site indices, ascending (column order of
``unique_subsys_configs``).
block_coords : list[tuple]
Coordinates of ``keep_sites`` (same order).
unique_subsys_configs : numpy.ndarray
``(D, 4)`` table of the distinct 4-site sub-configurations appearing
in the sector; row ``r`` is the basis label of RDM row ``r``, column
``j`` the local-state label at ``keep_sites[j]``.
axes : tuple[str, str], optional
Plaquette-plane axes.
Returns
-------
numpy.ndarray
``(D, 8)`` real array: the eight outward-link electric fields (two
per site) of each RDM basis row. Rows with identical signatures
share a boundary flux (one superselection block).
Notes
-----
In the gauge-invariant (electric) basis the link operators
``E_m{d}``/``E_p{d}`` are diagonal, so the electric field a link carries
for a given local-state label is just that operator's diagonal entry at
that label. Columns are laid out leg-major (site 0's two links, then
site 1's, ...).
"""
n_states = unique_subsys_configs.shape[0]
# outward_link_ops returns one link name per plane axis, so each of the
# len(keep_sites) sites contributes len(axes) outward links.
n_links = len(axes) * len(keep_sites)
flux_signatures = np.zeros((n_states, n_links))
for row in range(n_states):
# local-state label of each block site in this RDM row
sub_config = unique_subsys_configs[row]
for site_column, site in enumerate(keep_sites):
outward_ops_names = outward_link_ops(
block_coords[site_column], block_coords, axes
)
for op_idx, op_name in enumerate(outward_ops_names):
if op_name in self.ops:
# self.ops[op_name][site] is the (d, d) operator at this
# site; its diagonal holds the electric field per state.
field_per_state = np.real(
np.diagonal(self.ops[op_name][site])
)
else:
# Link points off the lattice (trivial OBC border):
# zero flux for every local state of this site.
field_per_state = np.zeros(int(self.loc_dims[site]))
# leg-major column: site 0's links first, then site 1's, ...
col = len(outward_ops_names) * site_column + op_idx
local_state = int(sub_config[site_column])
flux_signatures[row, col] = field_per_state[local_state]
return flux_signatures
def _setup_plaquette_basis(self, axes=("x", "y")):
"""Build the plaquette-block-LBO super-site operators and link generators.
``self.lvals`` is the coarse super-lattice: each site is a 2x2 block of
the original ``2 * self.lvals`` dressed lattice. For every position class
we reconstruct the block's four leg labels, build the gauge-invariant
super-site operators (:func:`QED_plq_supersite_operators`), and project
them in two stages -- the sparse embedding
``self.block_projectors[<class>_embedding]`` (prod -> D-dim subset) then
the dense block projector ``self.block_projectors[<class>]`` (D -> k). The
result fills
- ``self.ops`` -- per-super-site operators padded to ``max_loc_dim``,
- ``self.loc_dims`` -- the kept dimension of each super-site,
- ``self._plaq_link_ops`` / ``self._plaq_link_sectors`` -- the boundary
generators that impose inter-super-site link continuity in
:meth:`get_abelian_symmetry_sector`.
The two shared rungs of each super-bond are packed into a single combined
generator ``flux1 * (2*spin + 1) + flux2`` per direction, so one additive
U(1) constraint per direction matches both rungs at once -- fitting the
one-pair-per-direction contract of the link-sector builder.
"""
if self.dim != 2:
raise ValueError("Plaquette-block basis is only defined in 2D")
if self.block_projectors is None:
raise ValueError("plaq_basis=True requires block_projectors")
# Plaquette-block LBO is implemented for fully-OBC pure-gauge lattices with
# hard link-symmetry constraints. Reject the untested/unsupported combos
# up front with actionable errors rather than a confusing crash later
# (outward_link_ops and the boundary-flux generators assume open borders;
# the projectors are extracted at background=0; the inter-super-site
# continuity is a hard sector constraint, with no penalty mode).
if any(not obc for obc in self.has_obc):
raise NotImplementedError(
f"plaq_basis is only implemented for fully-OBC lattices, "
f"got has_obc={self.has_obc}."
)
if self.bg_list is not None:
raise NotImplementedError(
"plaq_basis does not support background charges (bg_list)."
)
if not self.link_symmetries:
raise NotImplementedError(
"plaq_basis requires link_symmetries=True (no penalty mode)."
)
# Operators the plaquette Hamiltonian needs (intra + inter-block seams).
ham_op_names = [
"E2_plq",
"E2_plq_px",
"E2_plq_py",
"B2_plq",
"B2_plq_px",
"B2_plq_mx",
"B2_plq_py",
"B2_plq_my",
"B2_plq_px_py",
"B2_plq_mx_py",
"B2_plq_px_my",
"B2_plq_mx_my",
]
link_op_names = [
f"E_plq_link_{sign}{d}" for d in self.directions for sign in "pm"
]
# injective packing base for the two rungs
flux_pack = 2 * self.spin + 1
# ---------------------------------------------------------------------------
# Tile the original lattice into disjoint 2x2 blocks (one per super-site),
# ordered to match the coarse-lattice zig-zag (== self.lvals site order).
original_lvals = [2 * length for length in self.lvals]
blocks, _, _, _ = tile_plaquette_blocks(original_lvals, self.has_obc, axes)
if len(blocks) != self.n_sites:
msg = f"super-lattice {self.lvals} expects {self.n_sites} blocks, tiling gave {len(blocks)}"
raise ValueError(msg)
# ---------------------------------------------------------------------------
# The dressed-site gauge-invariant operators are the SAME for every
# super-site (they depend only on spin/theory/dim, not on the block's
# position), so build them ONCE here and feed them to each per-class
# assembly below. Doing this avoids rebuilding -- and re-logging -- the
# whole GI basis once per position class.
class_labels = sorted({block["label"] for block in blocks})
# Fail fast if the projector dict is missing entries: each position class
# the tiling uses needs both its (D, k) projector "<label>" and its sparse
# embedding "<label>_embedding". The usual cause is calling
# extract_plaquette_block_projectors with the default return_full_space=False
# (no embeddings) or with only_classes/an instance-A lattice omitting a class.
missing = [
key
for label in class_labels
for key in (label, f"{label}_embedding")
if key not in self.block_projectors
]
if missing:
raise KeyError(
f"block_projectors is missing {missing}. Extract them with "
"extract_plaquette_block_projectors(..., return_full_space=True) "
f"for every position class the tiling uses ({class_labels})."
)
logger.info("----------------------------------------------------")
logger.info(
"PLAQUETTE BASIS: %s dressed lattice -> %s super-lattice",
original_lvals,
list(self.lvals),
)
logger.info(
"%d super-sites in %d position classes: %s",
self.n_sites,
len(class_labels),
class_labels,
)
dressed_gi_ops = QED_gauge_invariant_operators(
self.spin, self.pure_theory, self.dim
)
# ---------------------------------------------------------------------------
# Build + project the super-site operators ONCE per position class: every
# super-site of a class (e.g. all "px" right-edge blocks) shares the same
# four leg labels and the same projector, so we never rebuild per site.
projected_by_class = {}
for label in class_labels:
# A representative block of this class; its four sites' coordinates fix
# the per-leg gauge-invariant bases.
representative = next(b for b in blocks if b["label"] == label)
# Each leg's gauge-invariant label ("site", "site_px", "site_px_py", ...)
# read from its position in the ORIGINAL dressed lattice.
leg_labels = []
for coords in representative["coords"]:
site_label = get_site_label(
coords, original_lvals, self.has_obc, all_sites_equal=False
)
leg_labels.append(site_label)
# The (gauge-invariant-leg-product -> kept-block) isometry for this
# class, learned from the dressed RDM. Shape: (prod(leg dims), k).
# Two-stage projection inputs: the SPARSE embedding S (prod_dim x D)
# maps a GI-leg-product operator into the small sector-reachable subset
# basis; the dense block projector P (D x k) then truncates it.
embedding = self.block_projectors[f"{label}_embedding"] # sparse (prod, D)
block_proj = self.block_projectors[label] # dense (D, k)
logger.info("%s: legs: %s", label, leg_labels)
logger.info(
"full-space: %d subset: %d kept: %d",
embedding.shape[0],
embedding.shape[1],
block_proj.shape[1],
)
# Super-site operators on the GI-leg product space for these legs,
# reusing the shared dressed-site GI operators built above.
supersite_ops = QED_plq_supersite_operators(
self.spin,
self.pure_theory,
self.dim,
leg_labels=leg_labels,
gi_ops=dressed_gi_ops,
)
def project(name, embed=embedding, proj=block_proj, ops=supersite_ops):
# Two cheap stages (defaults bind this class's S / P / ops):
# 1) S^dag O S reduces O to the (D, D) subset basis via the sparse
# embedding; 2) P^dag (.) P truncates to the (k, k) kept basis.
# Never forms a dense (prod_dim, k) matrix -- ~30x faster than that.
o_subset = (embed.conj().T @ ops[name] @ embed).toarray()
return proj.conj().T @ o_subset @ proj
# The Hamiltonian operators, reduced to this class's kept block basis.
class_ops = {name: project(name) for name in ham_op_names}
# Boundary-link generators. Each super-bond shares TWO links (rungs 1,2),
# but get_abelian_symmetry_sector accepts only ONE generator per
# direction. So PACK the two rung fluxes injectively into one value,
# f1*(2s+1) + f2 (with f1,f2 in [-s, s]): then "combined_p(i) ==
# combined_m(j)" holds iff BOTH rungs match. The packed generator is a
# sum/scaling of diagonal operators, hence still diagonal in the kept
# basis (required by get_abelian_symmetry_sector's np.diagonal read).
for d in self.directions:
for sign in "pm":
rung1 = project(f"E_plq_{sign}{d}1") # lower shared link
rung2 = project(f"E_plq_{sign}{d}2") # upper shared link
# The packed flux is an exact integer; round away the ~1e-15
# projection noise so get_abelian_symmetry_sector (which
# compares link charges) keeps every continuous configuration.
class_ops[f"E_plq_link_{sign}{d}"] = np.round(
flux_pack * rung1 + rung2
)
projected_by_class[label] = class_ops
# ---------------------------------------------------------------------------
# Per-super-site kept dimension (k may differ by class under OBC).
self.loc_dims = np.array(
[self.block_projectors[b["label"]].shape[1] for b in blocks],
dtype=int,
)
# Show the kept dimension of each super-site arranged on the coarse
# super-lattice geometry (same grid layout as the dressed-site model's
# "LOCAL DIMENSION per SITE" print).
logger.info("====================================================")
logger.info("LOCAL DIMENSION per SUPER-SITE")
format_loc_dims(self.loc_dims, self.lvals)
# Per-super-site class labels (informational). No per-site gauge_basis is
# needed: each super-site already carries its own reduced operators.
self.lattice_labels = np.array([block["label"] for block in blocks])
self.gauge_basis = None
self.gauge_states = None
# Stack each operator as (n_super_sites, max_k, max_k), zero-padding every
# super-site's (k_s, k_s) block to the common max_k so the Numba sector
# kernels see a uniform array (only states < k_s ever appear in
# sector_configs, so the padding is never addressed).
max_loc_dim = int(np.max(self.loc_dims))
self.ops = {}
for name in ham_op_names + link_op_names:
stack = np.zeros((self.n_sites, max_loc_dim, max_loc_dim), dtype=complex)
for site_idx, block in enumerate(blocks):
loc_dim = self.loc_dims[site_idx]
stack[site_idx, :loc_dim, :loc_dim] = projected_by_class[
block["label"]
][name]
self.ops[name] = stack
# ---------------------------------------------------------------------------
# Link generators for get_abelian_symmetry_sector: one [left, right] pair
# per direction. The "- combined_m" makes the constraint a difference, so
# the sector keeps only configs with combined_p(i) == combined_m(i+d) on
# every super-bond -- i.e. the boundary links between super-sites are
# continuous. This is the external link symmetry, imposed as a hard
# constraint on the basis (not approximately).
self._plaq_link_ops = [
[self.ops[f"E_plq_link_p{d}"], -self.ops[f"E_plq_link_m{d}"]]
for d in self.directions
]
self._plaq_link_sectors = [0 for _ in self.directions]
[docs]
def build_plaquette_Hamiltonian(self, g, dtype_mode="auto"):
"""Assemble the plaquette-basis QED Hamiltonian.
Parameters
----------
dtype_mode : str or bool, optional
"auto", "real", "complex", or legacy bool flag.
"""
self.configure_dtype_mode(
dtype_mode=dtype_mode,
auto_mode=self._get_auto_dtype_mode(m=None, theta=0.0),
)
logger.info("----------------------------------------------------")
logger.info("BUILDING plaquette-SUPERSITE HAMILTONIAN")
# Hamiltonian Coefficients. E=g/2 matches the dressed normalization: the
# dressed "E2" is already the per-link energy, and the super-site electric
# operators likewise count each link once (see _setup_plaquette_basis).
self.coeffs = {"E": g / 2, "B": -1 / (2 * g)}
h_terms = {}
assert self.dim == 2, "Plaquette Hamiltonian only defined for dim >=2"
# -------------------------------------------------------------------------------
# ELECTRIC ENERGY (interior + outgoing-boundary links of each super-site)
for op_name in ["E2_plq", "E2_plq_px", "E2_plq_py"]:
h_terms[op_name] = LocalTerm(self.ops[op_name], op_name, **self.def_params)
self.H.add_term(h_terms[op_name].get_Hamiltonian(strength=self.coeffs["E"]))
# -------------------------------------------------------------------------------
# MAGNETIC ENERGY
# Intra-block plaquette loop (B2_plq already includes its h.c.).
h_terms["B2_plq"] = LocalTerm(self.ops["B2_plq"], "B2_plq", **self.def_params)
self.H.add_term(h_terms["B2_plq"].get_Hamiltonian(strength=self.coeffs["B"]))
for d in self.directions:
# Define the list of the 2 non trivial operators
op_names_list = [f"B2_plq_p{d}", f"B2_plq_m{d}"]
op_list = [self.ops[op] for op in op_names_list]
# Define the Hamiltonian term
h_terms[f"plq_{d}"] = TwoBodyTerm(
d, op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms[f"plq_{d}"].get_Hamiltonian(
strength=self.coeffs["B"], add_dagger=True
)
)
# Plaquette operator between sites:
op_names_list = ["B2_plq_px_py", "B2_plq_mx_py", "B2_plq_px_my", "B2_plq_mx_my"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["B2_plaq_xy"] = PlaquetteTerm(
["x", "y"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["B2_plaq_xy"].get_Hamiltonian(
strength=self.coeffs["B"], add_dagger=True
)
)
self.H.build(format=self.ham_format)
[docs]
def check_symmetries(self):
"""Check link-symmetry constraints on measured electric fields."""
# CHECK LINK SYMMETRIES
for ax in self.directions:
check_link_symmetry(
ax,
self.obs_list[f"E_p{ax}"],
self.obs_list[f"E_m{ax}"],
value=0,
sign=1,
)
[docs]
def QED_Hamiltonian_couplings(self, g, m=None, theta=0.0, magnetic_basis=False):
"""Set QED Hamiltonian couplings from physical parameters.
Parameters
----------
g : float
Gauge coupling.
m : float, optional
Bare mass parameter (used only with matter).
theta : float, optional
Topological-angle parameter.
magnetic_basis : bool, optional
If True, use the alternative magnetic-basis normalization.
Returns
-------
None
Stores couplings in self.coeffs.
"""
if not magnetic_basis:
if self.dim == 1:
E = g / 2
B = -1 / (2 * g)
elif self.dim == 2:
E = g / 2
B = -1 / (2 * g)
else:
E = g / 2
B = -1 / (2 * g)
# DICTIONARY WITH MODEL COEFFICIENTS
coeffs = {
"g": g,
"E": E, # ELECTRIC FIELD coupling
"B": B, # MAGNETIC FIELD coupling
"theta": -theta * g, # THETA TERM coupling
}
if not self.pure_theory and m is not None:
t = 1 / 2
coeffs |= {
"t": -complex(0, t), # integrated-gauge hopping coefficient
"tx_even": -complex(0, t), # x HOPPING (EVEN SITES)
"tx_odd": -complex(0, t), # x HOPPING (ODD SITES)
"ty_even": -t, # y HOPPING (EVEN SITES)
"ty_odd": t, # y HOPPING (ODD SITES)
"tz_even": -complex(0, t), # z HOPPING (EVEN SITES)
"tz_odd": complex(0, t), # z HOPPING (ODD SITES)
"m": m,
"m_odd": -m, # EFFECTIVE MASS for ODD SITES
"m_even": m, # EFFECTIVE MASS for EVEN SITES
}
else:
# DICTIONARY WITH MODEL COEFFICIENTS
coeffs = {
"g": g,
"E": -g,
"B": -0.5 / g,
}
self.coeffs = coeffs
[docs]
def overlap_QMB_state(self, name):
"""Return predefined benchmark basis configurations for selected labels.
Parameters
----------
name : str
Label of a predefined reference configuration.
Returns
-------
numpy.ndarray
Configuration in the model symmetry-sector basis.
"""
# ---------------------------------------------------------------------------
# 1D QED with matter: simple benchmark labels
# V : vacuum
# meson_center : local meson in the central (even, odd) pair
# Backward-compatible alias:
# PV -> V
if self.dim == 1 and not self.pure_theory:
name_norm = name.strip().lower().replace(" ", "_")
if name_norm == "pv":
name_norm = "v"
if name_norm in ("v", "meson_center"):
# Integrated case (matter-only local basis: 0/1 occupancies)
if self.spin == "integrated":
config_state = np.zeros(self.n_sites, dtype=np.int32)
# Integrated occupation convention:
# basis index 0 -> N=0 (empty), index 1 -> N=1 (occupied).
# Vacuum occupancy (even empty, odd occupied) is [0,1,0,1,...].
for ii in range(self.n_sites):
config_state[ii] = 0 if (ii % 2 == 0) else 1
if name_norm == "meson_center":
even_site = self.n_sites // 2
if even_site % 2 == 1:
even_site -= 1
even_site = max(0, min(even_site, self.n_sites - 2))
odd_site = even_site + 1
# Meson on central (even, odd) pair relative to V.
config_state[even_site] = 1
config_state[odd_site] = 0
return config_state
# Truncated case (hard-coded local-state labels for spin=1)
if np.isscalar(self.spin) and int(self.spin) == 1:
config_state = np.zeros(self.n_sites, dtype=np.int32)
if all(self.has_obc):
# OBC:
# left edge = 0, right edge = 1
# core: even -> 1, odd -> 3
if self.n_sites >= 1:
config_state[0] = 0
for ii in range(1, self.n_sites - 1):
config_state[ii] = 1 if (ii % 2 == 0) else 3
if self.n_sites >= 2:
config_state[self.n_sites - 1] = 1
else:
# PBC: alternance 1,3
for ii in range(self.n_sites):
config_state[ii] = 1 if (ii % 2 == 0) else 3
if name_norm == "meson_center":
even_site = self.n_sites // 2
if even_site % 2 == 1:
even_site -= 1
even_site = max(0, min(even_site, self.n_sites - 2))
odd_site = even_site + 1
config_state[even_site] = 4
config_state[odd_site] = 1
return config_state
raise NotImplementedError(
"1D labels 'V'/'meson_center' are implemented for "
"spin='integrated' or truncated spin=1."
)
# MINIMAL STRING CONFIGURATIONS IN 3D QED WITH BACKGROUND charges
if len(self.lvals) == 3 and self.bg_list == [-1, 0, 0, 0, 0, 0, 0, 1]:
if self.spin == 1:
if name == "S1":
config_state = np.array([4, 9, 9, 3, 6, 0, 3, 10])
elif name == "S2":
config_state = np.array([3, 9, 6, 0, 9, 3, 3, 11])
elif name == "S3":
config_state = np.array([1, 7, 9, 3, 9, 2, 3, 10])
elif name == "S4":
config_state = np.array([4, 9, 9, 3, 7, 3, 0, 8])
elif name == "S5":
config_state = np.array([1, 6, 9, 2, 9, 3, 3, 11])
elif name == "S6":
config_state = np.array([3, 9, 7, 3, 9, 3, 2, 8])
else:
raise ValueError(f"Unknown String name: {name}")
elif self.spin == 2:
if name == "S1":
config_state = np.array([11, 27, 27, 9, 22, 4, 9, 29])
elif name == "S2":
config_state = np.array([10, 27, 22, 4, 27, 9, 9, 30])
elif name == "S3":
config_state = np.array([10, 27, 23, 9, 27, 9, 8, 25])
elif name == "S4":
config_state = np.array([6, 22, 27, 8, 27, 9, 9, 30])
elif name == "S5":
config_state = np.array([6, 23, 27, 9, 27, 8, 9, 29])
elif name == "S6":
config_state = np.array([11, 27, 27, 9, 23, 9, 4, 25])
else:
raise ValueError(f"Unknown String name: {name}")
elif self.spin == 3:
if name == "S1":
config_state = np.array([20, 54, 47, 11, 54, 18, 18, 58])
elif name == "S2":
config_state = np.array([21, 54, 54, 18, 47, 11, 18, 57])
elif name == "S3":
config_state = np.array([21, 54, 54, 18, 48, 18, 11, 51])
elif name == "S4":
config_state = np.array([20, 54, 48, 18, 54, 18, 17, 51])
elif name == "S5":
config_state = np.array([14, 48, 54, 18, 54, 17, 18, 57])
elif name == "S6":
config_state = np.array([14, 47, 54, 17, 54, 18, 18, 58])
else:
raise ValueError(f"Unknown String name: {name}")
return config_state
else:
raise NotImplementedError("Only 3D QED states are supported")
[docs]
def print_state_config(self, config, amplitude=None):
"""Log a readable per-site decomposition of a QED basis configuration."""
if self.spin == "integrated":
raise NotImplementedError("It does not work with spin='integrated'")
logger.info("----------------------------------------------------")
msg = f"CONFIG {config}"
if amplitude is not None:
msg += f" |psi|^2={np.abs(amplitude)**2:.8f}"
logger.info(msg)
logger.info("----------------------------------------------------")
# Choose width (sign + digits).
max_abs = 1
max_abs = max(max_abs, int(abs(getattr(self, "spin", 1))))
max_abs = max(max_abs, int(abs(getattr(self, "background", 0))))
entry_width = len(str(max_abs)) + 2
logger.info("%s%s", " " * 19, self._format_header(entry_width))
for ii, cfg_idx in enumerate(config):
# Given the 1D point on the d-dimensional lattice, get the corresponding coords
coords = zig_zag(self.lvals, ii)
# Get the corresponding local basis state from the gauge sector basis
loc_basis_state = self.gauge_states_per_site[ii][cfg_idx]
state_str = (
"["
+ " ".join(f"{val:>{entry_width}d}" for val in loc_basis_state)
+ " ]"
)
logger.info("SITE %2d in %s, state %3d: %s", ii, coords, cfg_idx, state_str)
def _local_state_labels(self) -> list[str]:
"""Return column labels used by :meth:`print_state_config`."""
labels: list[str] = []
# Background first (only if present in your effective gauge states)
if getattr(self, "background", 0) > 0:
labels.append("bg")
# Matter occupation (only if not pure theory)
if not getattr(self, "pure_theory", True):
labels.append("M")
# Links: -x,-y,-z,+x,+y,+z up to dim
dim = int(getattr(self, "dim", 0))
dirs = "xyz"[:dim]
labels += [f"-{d}" for d in dirs] + [f"+{d}" for d in dirs]
return labels
def _format_header(self, entry_width: int) -> str:
"""Format the header row for :meth:`print_state_config`."""
labels = self._local_state_labels()
header = "[" + " ".join(f"{lab:>{entry_width}s}" for lab in labels) + " ]"
return header