"""Operator factories and gauge-invariant local bases for SU(2) lattice models."""
import logging
import os
from concurrent.futures import ThreadPoolExecutor
from itertools import product
import numpy as np
from numpy.linalg import matrix_rank
from scipy.sparse import csr_matrix, identity, isspmatrix, kron
from scipy.sparse.linalg import norm
from sympy import S
from edlgt.modeling import LGT_border_configs, get_lattice_borders_labels
from edlgt.modeling import qmb_operator as qmb_op
from edlgt.tools import validate_parameters
from .bose_fermi_operators import fermi_operators as SU2_matter_operators
from .spin_operators import (
SU2_generators,
get_Pauli_operators,
get_spin_operators,
m_values,
spin_space,
)
from .SU2_rishons import SU2_Rishon, SU2_Rishon_gen
from .SU2_singlets import (
SU2_singlet_canonical_vector,
get_spin_Hilbert_spaces,
get_SU2_singlets,
)
logger = logging.getLogger(__name__)
__all__ = [
"SU2_dressed_site_operators",
"SU2_gen_dressed_site_operators",
"SU2_check_gauss_law",
"SU2_gauge_invariant_states",
"SU2_gauge_invariant_operators",
"SU2_gauge_integrated_operators",
]
[docs]
def SU2_dressed_site_operators(
spin, pure_theory, lattice_dim, background=0, n_flavors=None
):
"""Build SU(2) dressed-site operators using the hardcoded rishon construction.
Parameters
----------
spin : float
Maximum rishon spin representation.
pure_theory : bool
If ``True``, exclude matter fields.
lattice_dim : int
Number of spatial lattice dimensions (1, 2, or 3).
background : int, optional
Maximum background-charge irrep included at each site.
n_flavors : int, optional
Number of matter flavors. Defaults to ``1`` for matter theories.
Multi-flavor operators are not yet implemented; passing
``n_flavors > 1`` raises :class:`NotImplementedError`.
Returns
-------
dict
Dictionary of dressed-site operators used by the SU(2) model builders.
"""
validate_parameters(
spin_list=[spin], pure_theory=pure_theory, lattice_dim=lattice_dim
)
n_flavors_eff = 0 if pure_theory else (1 if n_flavors is None else int(n_flavors))
if n_flavors_eff > 1:
raise NotImplementedError(
"Flavored matter operators (n_flavors > 1) are not yet supported "
"in SU2_dressed_site_operators. The flavored gauge-invariant basis "
"is available via SU2_gauge_invariant_states(..., n_flavors=N); the "
"per-flavor fermionic operators and the corresponding Hamiltonian "
"terms are tracked in docs/su2_singlets_optimization/follow_ups.md "
"(F2). Use n_flavors=1 (or omit) for the single-flavor build."
)
# Lattice directions
dimensions = "xyz"[:lattice_dim]
# Get SU2 rishon operators
in_ops = SU2_Rishon(spin).ops
in_ops |= SU2_generators(spin)
if background > 0:
in_ops["T2_bg"] = SU2_generators(background)["T2"]
if not pure_theory:
in_ops |= SU2_matter_operators(has_spin=True, colors=True)
in_ops |= SU2_generators(1 / 2, matter=True)
for op in in_ops.keys():
in_ops[op] = csr_matrix(in_ops[op])
# Dictionary for dressed site operators
ops = {}
if lattice_dim == 1:
# T generators for electric term
for op in ["T2", "P", "Tx", "Ty", "Tz"]:
ops[f"{op}_mx"] = qmb_op(in_ops, [op, "Iz"])
ops[f"{op}_px"] = qmb_op(in_ops, ["Iz", op])
# Local scalar transporter obtained after contracting the intermediate
# color index inside a dressed site.
ops["W"] = 0
for col in "rg":
ops["W"] += qmb_op(in_ops, [f"Z{col}_dag_P", f"Z{col}"])
if not pure_theory:
# Update Electric operators
for op in ops.keys():
ops[op] = kron(in_ops["ID_psi"], ops[op])
# Add Hopping operators
for side in "pm":
ops[f"Q{side}x_dag"] = 0
for col in "rg":
ops["Qpx_dag"] += qmb_op(in_ops, [f"psi_{col}_dag_P", "P", f"Z{col}"])
ops["Qmx_dag"] += qmb_op(in_ops, [f"psi_{col}_dag_P", f"Z{col}", "Iz"])
# add their dagger operators
Qs = {}
for op in ops:
dag_op = op.replace("_dag", "")
Qs[dag_op] = csr_matrix(ops[op].conj().transpose())
ops |= Qs
elif lattice_dim == 2:
# T generators for electric term
for op in ["T2", "P", "Tx", "Ty", "Tz"]:
ops[f"{op}_mx"] = qmb_op(in_ops, [op, "Iz", "Iz", "Iz"])
ops[f"{op}_my"] = qmb_op(in_ops, ["Iz", op, "Iz", "Iz"])
ops[f"{op}_px"] = qmb_op(in_ops, ["Iz", "Iz", op, "Iz"])
ops[f"{op}_py"] = qmb_op(in_ops, ["Iz", "Iz", "Iz", op])
# Corner Operators
for corner in ["px,py", "py,mx", "mx,my", "my,px"]:
ops[f"C_{corner}"] = 0
for col in ["r", "g"]:
ops["C_px,py"] += qmb_op(in_ops, ["Iz", "Iz", f"Z{col}_P", f"Z{col}_dag"])
ops["C_py,mx"] += qmb_op(in_ops, [f"P_Z{col}_dag", "P", "P", f"Z{col}"])
ops["C_mx,my"] += qmb_op(in_ops, [f"Z{col}_P", f"Z{col}_dag", "Iz", "Iz"])
ops["C_my,px"] += qmb_op(in_ops, ["Iz", f"Z{col}_P", f"Z{col}_dag", "Iz"])
if not pure_theory:
# Update Electric and Corner operators
for op in ops.keys():
ops[op] = kron(in_ops["ID_psi"], ops[op])
# Add Hopping operators
for side in "pm":
for ax in dimensions:
ops[f"Q{side}{ax}_dag"] = 0
for col in "rg":
# ----------------------------------------------------------------------
op_list = [f"psi_{col}_dag_P", f"Z{col}", "Iz", "Iz", "Iz"]
ops["Qmx_dag"] += qmb_op(in_ops, op_list)
# ----------------------------------------------------------------------
op_list = [f"psi_{col}_dag_P", "P", f"Z{col}", "Iz", "Iz"]
ops["Qmy_dag"] += qmb_op(in_ops, op_list)
# ----------------------------------------------------------------------
op_list = [f"psi_{col}_dag_P", "P", "P", f"Z{col}", "Iz"]
ops["Qpx_dag"] += qmb_op(in_ops, op_list)
# ----------------------------------------------------------------------
op_list = [f"psi_{col}_dag_P", "P", "P", "P", f"Z{col}"]
ops["Qpy_dag"] += qmb_op(in_ops, op_list)
# add their dagger operators
Qs = {}
for op in ops:
dag_op = op.replace("_dag", "")
Qs[dag_op] = csr_matrix(ops[op].conj().transpose())
ops |= Qs
elif lattice_dim == 3:
# T generators for electric term
for op in ["T2", "P", "Tx", "Ty", "Tz"]:
ops[f"{op}_mx"] = qmb_op(in_ops, [op, "Iz", "Iz", "Iz", "Iz", "Iz"])
ops[f"{op}_my"] = qmb_op(in_ops, ["Iz", op, "Iz", "Iz", "Iz", "Iz"])
ops[f"{op}_mz"] = qmb_op(in_ops, ["Iz", "Iz", op, "Iz", "Iz", "Iz"])
ops[f"{op}_px"] = qmb_op(in_ops, ["Iz", "Iz", "Iz", op, "Iz", "Iz"])
ops[f"{op}_py"] = qmb_op(in_ops, ["Iz", "Iz", "Iz", "Iz", op, "Iz"])
ops[f"{op}_pz"] = qmb_op(in_ops, ["Iz", "Iz", "Iz", "Iz", "Iz", op])
# Corner Operators
corner_list = []
for pdir in ["xy", "xz", "yz"]:
# DEFINE THE LIST OF CORNER OPERATORS
corner_list += [
f"p{pdir[0]},p{pdir[1]}",
f"p{pdir[1]},m{pdir[0]}",
f"m{pdir[1]},p{pdir[0]}",
f"m{pdir[0]},m{pdir[1]}",
]
for corner in corner_list:
ops[f"C_{corner}"] = 0
for col in "rg":
# --------------------------------------------------------------------------
# XY Plane
op_list = ["Iz", "Iz", "Iz", f"Z{col}_P", f"Z{col}_dag", "Iz"]
ops["C_px,py"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = [f"P_Z{col}_dag", "P", "P", "P", f"Z{col}", "Iz"]
ops["C_py,mx"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = [f"Z{col}_P", f"Z{col}_dag", "Iz", "Iz", "Iz", "Iz"]
ops["C_mx,my"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = ["Iz", f"Z{col}_P", "P", f"Z{col}_dag", "Iz", "Iz"]
ops["C_my,px"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
# XZ Plane
op_list = ["Iz", "Iz", "Iz", f"Z{col}_P", "P", f"Z{col}_dag"]
ops["C_px,pz"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = [f"P_Z{col}_dag", "P", "P", "P", "P", f"Z{col}"]
ops["C_pz,mx"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = [f"Z{col}_P", "P", f"Z{col}_dag", "Iz", "Iz", "Iz"]
ops["C_mx,mz"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = ["Iz", "Iz", f"Z{col}_P", f"Z{col}_dag", "Iz", "Iz"]
ops["C_mz,px"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
# YZ Plane
op_list = ["Iz", "Iz", "Iz", "Iz", f"Z{col}_P", f"Z{col}_dag"]
ops["C_py,pz"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = ["Iz", f"P_Z{col}_dag", "P", "P", "P", f"Z{col}"]
ops["C_pz,my"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = ["Iz", f"Z{col}_P", f"Z{col}_dag", "Iz", "Iz", "Iz"]
ops["C_my,mz"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = ["Iz", "Iz", f"Z{col}_P", "P", f"Z{col}_dag", "Iz"]
ops["C_mz,py"] += qmb_op(in_ops, op_list)
# THETA TERM
ops["EzC_px,py"] = 0
ops["ExC_py,pz"] = 0
ops["EyC_px,pz"] = 0
sigma_ops = get_spin_operators(0.5)
for nu in "xyz":
Ez = ops[f"T{nu}_pz"] + ops[f"T{nu}_mz"]
Ey = ops[f"T{nu}_py"] + ops[f"T{nu}_my"]
Ex = ops[f"T{nu}_px"] + ops[f"T{nu}_mx"]
for ii, c1 in enumerate("rg"):
for jj, c2 in enumerate("rg"):
factor = sigma_ops[f"S{nu}"].todense()[jj, ii]
# ------------------------------------------------------------------
op_list = ["Iz", "Iz", "Iz", f"Z{c1}_P", f"Z{c2}_dag", "Iz"]
ops["EzC_px,py"] += 1j * qmb_op(in_ops, op_list) @ (factor * Ez)
# ------------------------------------------------------------------
op_list = ["Iz", "Iz", "Iz", f"Z{c1}_P", "P", f"Z{c2}_dag"]
ops["EyC_px,pz"] += -1j * qmb_op(in_ops, op_list) @ (factor * Ey)
# ------------------------------------------------------------------
op_list = ["Iz", "Iz", "Iz", "Iz", f"Z{c1}_P", f"Z{c2}_dag"]
ops["ExC_py,pz"] += 1j * qmb_op(in_ops, op_list) @ (factor * Ex)
# ------------------------------------------------------------------
if not pure_theory:
# Update Electric and Corner operators
for op in ops.keys():
ops[op] = kron(in_ops["ID_psi"], ops[op])
# Add Hopping operators
for side in "pm":
for ax in dimensions:
ops[f"Q{side}{ax}_dag"] = 0
for col in "rg":
# ----------------------------------------------------------------------
op_list = [f"psi_{col}_dag_P", f"Z{col}", "Iz", "Iz", "Iz", "Iz", "Iz"]
ops["Qmx_dag"] += qmb_op(in_ops, op_list)
# ----------------------------------------------------------------------
op_list = [f"psi_{col}_dag_P", "P", f"Z{col}", "Iz", "Iz", "Iz", "Iz"]
ops["Qmy_dag"] += qmb_op(in_ops, op_list)
# ----------------------------------------------------------------------
op_list = [f"psi_{col}_dag_P", "P", "P", f"Z{col}", "Iz", "Iz", "Iz"]
ops["Qmz_dag"] += qmb_op(in_ops, op_list)
# ----------------------------------------------------------------------
op_list = [f"psi_{col}_dag_P", "P", "P", "P", f"Z{col}", "Iz", "Iz"]
ops["Qpx_dag"] += qmb_op(in_ops, op_list)
# ----------------------------------------------------------------------
op_list = [f"psi_{col}_dag_P", "P", "P", "P", "P", f"Z{col}", "Iz"]
ops["Qpy_dag"] += qmb_op(in_ops, op_list)
# ----------------------------------------------------------------------
op_list = [f"psi_{col}_dag_P", "P", "P", "P", "P", "P", f"Z{col}"]
ops["Qpz_dag"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
# add their dagger operators
Qs = {}
for op in ops:
dag_op = op.replace("_dag", "")
Qs[dag_op] = csr_matrix(ops[op].conj().transpose())
ops |= Qs
# -----------------------------------------------------------------------------
if not pure_theory:
# Psi NUMBER OPERATORS
for label in ["r", "g", "tot", "single", "pair", "zero"]:
op_list = [f"N_{label}"] + ["Iz" for _ in range(2 * lattice_dim)]
ops[f"N_{label}"] = qmb_op(in_ops, op_list)
ops["N-1"] = ops["N_tot"] - identity(ops["N_tot"].shape[0])
# -----------------------------------------------------------------------------
# CASIMIR/ELECTRIC OPERATOR
ops[f"E_square"] = 0
for s in "pm":
for d in dimensions:
ops[f"E_square"] += 0.5 * ops[f"T2_{s}{d}"]
# -----------------------------------------------------------------------------
# BACKGROUND FIELD OPERATORS
if background > 0:
bg_len = 0
j_bg_set = np.arange(0, spin_space(background), 1) / 2
for irrep in j_bg_set:
bg_len += len(m_values(irrep))
for op in ops.keys():
ops[op] = kron(identity(bg_len), ops[op])
if pure_theory:
id_list = ["Iz" for _ in range(2 * lattice_dim)]
else:
id_list = ["ID_psi"] + ["Iz" for _ in range(2 * lattice_dim)]
ops["bg"] = qmb_op(in_ops, ["T2_bg"] + id_list)
# GAUSS LAW VIOLATING OPERATORS
if lattice_dim == 1:
# Add background hopping operators
for side in "pm":
for ax in dimensions:
ops[f"V{side}{ax}_dag"] = 0
for col in "rg":
# ----------------------------------------------------------------------
op_list = [f"Z{col}_dag", "P_psi", f"Z{col}", "Iz"]
ops["Vmx_dag"] += qmb_op(in_ops, op_list)
# ----------------------------------------------------------------------
op_list = [f"Z{col}_dag", "P_psi", "P", f"Z{col}"]
ops["Vpx_dag"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
# add their dagger operators
ops["Vmx"] = csr_matrix(ops["Vmx_dag"].conj().transpose())
ops["Vpx"] = csr_matrix(ops["Vpx_dag"].conj().transpose())
return ops
[docs]
def SU2_gen_dressed_site_operators(
spin, pure_theory, lattice_dim, background=0, n_flavors=None
):
"""Build SU(2) dressed-site operators using the generalized rishon construction.
Parameters
----------
spin : float
Maximum rishon spin representation.
pure_theory : bool
If ``True``, exclude matter fields.
lattice_dim : int
Number of spatial lattice dimensions (1, 2, or 3).
background : int, optional
Maximum background-charge irrep included at each site.
n_flavors : int, optional
Number of matter flavors. Defaults to ``1`` for matter theories.
Multi-flavor operators are not yet implemented; passing
``n_flavors > 1`` raises :class:`NotImplementedError`.
Returns
-------
dict
Dictionary of dressed-site operators for the generalized SU(2) model.
"""
n_flavors_eff = 0 if pure_theory else (1 if n_flavors is None else int(n_flavors))
if n_flavors_eff > 1:
raise NotImplementedError(
"Flavored matter operators (n_flavors > 1) are not yet supported "
"in SU2_gen_dressed_site_operators. See "
"docs/su2_singlets_optimization/follow_ups.md (F2)."
)
validate_parameters(
spin_list=[spin], pure_theory=pure_theory, lattice_dim=lattice_dim
)
# Lattice directions
dimensions = "xyz"[:lattice_dim]
# Get SU2 rishon operators
in_ops = SU2_Rishon_gen(spin).ops
in_ops |= SU2_generators(spin)
if background > 0:
in_ops["T2_bg"] = SU2_generators(background)["T2"]
# Add SU2 matter operators / generators
if not pure_theory:
in_ops |= SU2_matter_operators(has_spin=True, colors=True)
in_ops |= SU2_generators(1 / 2, matter=True)
for op in in_ops.keys():
in_ops[op] = csr_matrix(in_ops[op])
# Dictionary for dressed site operators
ops = {}
if lattice_dim == 1:
# T generators for electric term
for op in ["T2", "P"]:
ops[f"{op}_mx"] = qmb_op(in_ops, [op, "Iz"])
ops[f"{op}_px"] = qmb_op(in_ops, ["Iz", op])
if not pure_theory:
# Update Electric operators
for op in ops.keys():
ops[op] = kron(in_ops["ID_psi"], ops[op])
# Hopping operators
# Q1dag = psi_r_dag Zg_dag - psi_g_dag Zr_dag
# ------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "Zg_dag", "Iz"]
op_list2 = ["psi_g_dag_P", "Zr_dag", "Iz"]
ops["Q1_mx_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "Zg_dag"]
op_list2 = ["psi_g_dag_P", "P", "Zr_dag"]
ops["Q1_px_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
# Hopping operators Q2dag = psi_r_dag Zr + psi_g_dag Zg
op_list1 = ["psi_r_dag_P", "Zr", "Iz"]
op_list2 = ["psi_g_dag_P", "Zg", "Iz"]
ops["Q2_mx_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "Zr"]
op_list2 = ["psi_g_dag_P", "P", "Zg"]
ops["Q2_px_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
# The reduced dressed-site fermionic mode used for the NG
# covariance matrix is the singlet combination Q = Q1 + Q2.
# For j=1/2 this reproduces exactly the hardcoded Q operators after
# projection to the local gauge-invariant basis.
ops["Qmx_dag"] = ops["Q1_mx_dag"] + ops["Q2_mx_dag"]
ops["Qpx_dag"] = ops["Q1_px_dag"] + ops["Q2_px_dag"]
# ------------------------------------------------------------------
# add their dagger operators
Qs = {}
for op in ops:
dag_op = op.replace("_dag", "")
Qs[dag_op] = csr_matrix(ops[op].conj().transpose())
ops |= Qs
# In the generic 1D rishon decomposition, the local link matrix is
# resolved into two half-link channels denoted A and B. After
# contracting the intermediate color index inside a dressed site,
# the effective scalar transporter entering the reduced LxL
# covariance matrix is
#
# W = W_BA + W_AB - W_BB - W_AA
#
# with W_XY = sum_c ZX_c^dag P ZY_c.
#
# We build it explicitly on the three dressed-site factors
# [matter, minus-rishon, plus-rishon]. This ordering matches the
# parity placement used by the Q endpoint operators and reproduces
# exactly the hardcoded j=1/2 transporter after GI projection.
ops["W"] = 0
for col in "rg":
ops["W"] += qmb_op(in_ops, ["ID_psi", f"ZB_{col}_dag_P", f"ZA_{col}"])
ops["W"] += qmb_op(in_ops, ["ID_psi", f"ZA_{col}_dag_P", f"ZB_{col}"])
ops["W"] -= qmb_op(in_ops, ["ID_psi", f"ZB_{col}_dag_P", f"ZB_{col}"])
ops["W"] -= qmb_op(in_ops, ["ID_psi", f"ZA_{col}_dag_P", f"ZA_{col}"])
elif lattice_dim == 2:
# T generators for electric term
for op in ["T2", "P"]:
ops[f"{op}_mx"] = qmb_op(in_ops, [op, "Iz", "Iz", "Iz"])
ops[f"{op}_my"] = qmb_op(in_ops, ["Iz", op, "Iz", "Iz"])
ops[f"{op}_px"] = qmb_op(in_ops, ["Iz", "Iz", op, "Iz"])
ops[f"{op}_py"] = qmb_op(in_ops, ["Iz", "Iz", "Iz", op])
# Corner Operators
for l1, l2 in product(["A", "B"], repeat=2):
for corner in ["px,py", "py,mx", "mx,my", "my,px"]:
ops[f"C{l1}{l2}_{corner}"] = 0
for s in ["r", "g"]:
# ------------------------------------------------------------------
op_list = ["Iz", "Iz", f"Z{l1}_{s}_P", f"Z{l2}_{s}_dag"]
ops[f"C{l1}{l2}_px,py"] += qmb_op(in_ops, op_list)
# ------------------------------------------------------------------
op_list = [f"P_Z{l2}_{s}_dag", "P", "P", f"Z{l1}_{s}"]
ops[f"C{l1}{l2}_py,mx"] += qmb_op(in_ops, op_list)
# ------------------------------------------------------------------
op_list = [f"Z{l1}_{s}_P", f"Z{l2}_{s}_dag", "Iz", "Iz"]
ops[f"C{l1}{l2}_mx,my"] += qmb_op(in_ops, op_list)
# ------------------------------------------------------------------
op_list = ["Iz", f"Z{l1}_{s}_P", f"Z{l2}_{s}_dag", "Iz"]
ops[f"C{l1}{l2}_my,px"] += qmb_op(in_ops, op_list)
if not pure_theory:
# Update Electric and Corner operators
for op in ops.keys():
ops[op] = kron(in_ops["ID_psi"], ops[op])
# Hopping operators
# Q1dag = psi_r_dag Zg_dag - psi_g_dag Zr_dag
# ------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "Zg_dag", "Iz", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "Zr_dag", "Iz", "Iz", "Iz"]
ops["Q1_mx_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "Zg_dag", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "Zr_dag", "Iz", "Iz"]
ops["Q1_my_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "Zg_dag", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "P", "Zr_dag", "Iz"]
ops["Q1_px_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "P", "Zg_dag"]
op_list2 = ["psi_g_dag_P", "P", "P", "P", "Zr_dag"]
ops["Q1_py_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
# Hopping operators Q2dag = psi_r_dag Zr + psi_g_dag Zg
op_list1 = ["psi_r_dag_P", "Zr", "Iz", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "Zg", "Iz", "Iz", "Iz"]
ops["Q2_mx_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "Zr", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "Zg", "Iz", "Iz"]
ops["Q2_my_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "Zr", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "P", "Zg", "Iz"]
ops["Q2_px_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "P", "Zr"]
op_list2 = ["psi_g_dag_P", "P", "P", "P", "Zg"]
ops["Q2_py_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# ------------------------------------------------------------------
# add their dagger operators
Qs = {}
for op in ops:
dag_op = op.replace("_dag", "")
Qs[dag_op] = csr_matrix(ops[op].conj().transpose())
ops |= Qs
elif lattice_dim == 3:
# T generators for electric term
for op in ["T2", "P"]:
ops[f"{op}_mx"] = qmb_op(in_ops, [op, "Iz", "Iz", "Iz", "Iz", "Iz"])
ops[f"{op}_my"] = qmb_op(in_ops, ["Iz", op, "Iz", "Iz", "Iz", "Iz"])
ops[f"{op}_mz"] = qmb_op(in_ops, ["Iz", "Iz", op, "Iz", "Iz", "Iz"])
ops[f"{op}_px"] = qmb_op(in_ops, ["Iz", "Iz", "Iz", op, "Iz", "Iz"])
ops[f"{op}_py"] = qmb_op(in_ops, ["Iz", "Iz", "Iz", "Iz", op, "Iz"])
ops[f"{op}_pz"] = qmb_op(in_ops, ["Iz", "Iz", "Iz", "Iz", "Iz", op])
# Corner Operators
corner_list = []
for pdir in ["xy", "xz", "yz"]:
# DEFINE THE LIST OF CORNER OPERATORS
corner_list += [
f"p{pdir[0]},p{pdir[1]}",
f"p{pdir[1]},m{pdir[0]}",
f"m{pdir[1]},p{pdir[0]}",
f"m{pdir[0]},m{pdir[1]}",
]
for l1, l2 in product(["A", "B"], repeat=2):
for corner in corner_list:
ops[f"C{l1}{l2}_{corner}"] = 0
for s in ["r", "g"]:
# XY Plane
# --------------------------------------------------------------------------
op_list = ["Iz", "Iz", "Iz", f"Z{l1}_{s}_P", f"Z{l2}_{s}_dag", "Iz"]
ops[f"C{l1}{l2}_px,py"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = [f"P_Z{l2}_{s}_dag", "P", "P", "P", f"Z{l1}_{s}", "Iz"]
ops[f"C{l1}{l2}_py,mx"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = [f"Z{l1}_{s}_P", f"Z{l2}_{s}_dag", "Iz", "Iz", "Iz", "Iz"]
ops[f"C{l1}{l2}_mx,my"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = ["Iz", f"Z{l1}_{s}_P", "P", f"Z{l2}_{s}_dag", "Iz", "Iz"]
ops[f"C{l1}{l2}_my,px"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
# XZ Plane
op_list = ["Iz", "Iz", "Iz", f"Z{l1}_{s}_P", "P", f"Z{l2}_{s}_dag"]
ops[f"C{l1}{l2}_px,pz"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = [f"P_Z{l2}_{s}_dag", "P", "P", "P", "P", f"Z{l1}_{s}"]
ops[f"C{l1}{l2}_pz,mx"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = [f"Z{l1}_{s}_P", "P", f"Z{l2}_{s}_dag", "Iz", "Iz", "Iz"]
ops[f"C{l1}{l2}_mx,mz"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = ["Iz", "Iz", f"Z{l1}_{s}_P", f"Z{l2}_{s}_dag", "Iz", "Iz"]
ops[f"C{l1}{l2}_mz,px"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
# YZ Plane
op_list = ["Iz", "Iz", "Iz", "Iz", f"Z{l1}_{s}_P", f"Z{l2}_{s}_dag"]
ops[f"C{l1}{l2}_py,pz"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = ["Iz", f"P_Z{l2}_{s}_dag", "P", "P", "P", f"Z{l1}_{s}"]
ops[f"C{l1}{l2}_pz,my"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = ["Iz", f"Z{l1}_{s}_P", f"Z{l2}_{s}_dag", "Iz", "Iz", "Iz"]
ops[f"C{l1}{l2}_my,mz"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
op_list = ["Iz", "Iz", f"Z{l1}_{s}_P", "P", f"Z{l2}_{s}_dag", "Iz"]
ops[f"C{l1}{l2}_mz,py"] += qmb_op(in_ops, op_list)
# --------------------------------------------------------------------------
if not pure_theory:
# Update Electric and Corner operators
for op in ops.keys():
ops[op] = kron(in_ops["ID_psi"], ops[op])
# Add Hopping operators
# Q1dag = psi_r_dag Zg_dag - psi_g_dag Zr_dag
op_list1 = ["psi_r_dag_P", "Zg_dag", "Iz", "Iz", "Iz", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "Zr_dag", "Iz", "Iz", "Iz", "Iz", "Iz"]
ops["Q1_mx_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# --------------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "Zg_dag", "Iz", "Iz", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "Zr_dag", "Iz", "Iz", "Iz", "Iz"]
ops["Q1_my_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# --------------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "Zg_dag", "Iz", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "P", "Zr_dag", "Iz", "Iz", "Iz"]
ops["Q1_mz_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# --------------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "P", "Zg_dag", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "P", "P", "Zr_dag", "Iz", "Iz"]
ops["Q1_px_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# --------------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "P", "P", "Zg_dag", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "P", "P", "P", "Zr_dag", "Iz"]
ops["Q1_py_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# --------------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "P", "P", "P", "Zg_dag"]
op_list2 = ["psi_g_dag_P", "P", "P", "P", "P", "P", "Zr_dag"]
ops["Q1_pz_dag"] = qmb_op(in_ops, op_list1) - qmb_op(in_ops, op_list2)
# Hopping operators Q2dag = psi_r_dag Zr + psi_g_dag Zg
# --------------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "Zr", "Iz", "Iz", "Iz", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "Zg", "Iz", "Iz", "Iz", "Iz", "Iz"]
ops["Q2_mx_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# --------------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "Zr", "Iz", "Iz", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "Zg", "Iz", "Iz", "Iz", "Iz"]
ops["Q2_my_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# --------------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "Zr", "Iz", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "P", "Zg", "Iz", "Iz", "Iz"]
ops["Q2_mz_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# --------------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "P", "Zr", "Iz", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "P", "P", "Zg", "Iz", "Iz"]
ops["Q2_px_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# --------------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "P", "P", "Zr", "Iz"]
op_list2 = ["psi_g_dag_P", "P", "P", "P", "P", "Zg", "Iz"]
ops["Q2_py_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# --------------------------------------------------------------------------
op_list1 = ["psi_r_dag_P", "P", "P", "P", "P", "P", "Zr"]
op_list2 = ["psi_g_dag_P", "P", "P", "P", "P", "P", "Zg"]
ops["Q2_pz_dag"] = qmb_op(in_ops, op_list1) + qmb_op(in_ops, op_list2)
# --------------------------------------------------------------------------
# add their dagger operators
Qs = {}
for op in ops:
dag_op = op.replace("_dag", "")
Qs[dag_op] = csr_matrix(ops[op].conj().transpose())
ops |= Qs
# -----------------------------------------------------------------------------
if not pure_theory:
# Psi NUMBER OPERATORS
for label in ["r", "g", "tot", "single", "pair", "zero"]:
ops[f"N_{label}"] = qmb_op(
in_ops, [f"N_{label}"] + ["Iz" for i in range(2 * lattice_dim)]
)
# CASIMIR/ELECTRIC OPERATOR
ops[f"E_square"] = 0
for s in "pm":
for d in dimensions:
ops[f"E_square"] += 0.5 * ops[f"T2_{s}{d}"]
if background > 0:
bg_len = 0
j_bg_set = np.arange(0, spin_space(background), 1) / 2
for irrep in j_bg_set:
bg_len += len(m_values(irrep))
for op in ops.keys():
ops[op] = kron(identity(bg_len), ops[op])
if pure_theory:
id_list = ["Iz" for _ in range(2 * lattice_dim)]
else:
id_list = ["ID_psi"] + ["Iz" for _ in range(2 * lattice_dim)]
ops["bg"] = qmb_op(in_ops, ["T2_bg"] + id_list)
return ops
def SU2_gauge_integrated_operators():
# We assume an internal ordering like [red , green]
in_ops = get_Pauli_operators()
in_ops["Ir"] = identity(2, dtype=float)
in_ops["Ig"] = identity(2, dtype=float)
ops = {}
ops["I"] = qmb_op(in_ops, ["Ir", "Ig"])
# Hoppping operators
ops["Q_r"] = qmb_op(in_ops, ["Sp", "Ig"])
ops["Q_r_dag"] = qmb_op(in_ops, ["Sm", "Sz"])
ops["Q_g"] = qmb_op(in_ops, ["Sz", "Sp"])
ops["Q_g_dag"] = qmb_op(in_ops, ["Ir", "Sm"])
# Mass operators
ops["Sm_r"] = qmb_op(in_ops, ["Sm", "Ig"])
ops["Sp_r"] = qmb_op(in_ops, ["Sp", "Ig"])
ops["Sm_g"] = qmb_op(in_ops, ["Ir", "Sm"])
ops["Sp_g"] = qmb_op(in_ops, ["Ir", "Sp"])
ops["N_r"] = ops["Sm_r"] @ ops["Sp_r"]
ops["N_g"] = ops["Sm_g"] @ ops["Sp_g"]
ops["N_tot"] = ops["N_r"] + ops["N_g"]
# Diagonal electric operators
ops["Sz_rg"] = qmb_op(in_ops, ["Sz", "Sz"])
ops["inv_Sz"] = ops["I"] - ops["Sz_rg"]
# Off-diagonal QxQy electric operators
ops["Sp_r_Sm_g"] = qmb_op(in_ops, ["Sp", "Sm"])
ops["Sm_r_Sp_g"] = qmb_op(in_ops, ["Sm", "Sp"])
# Off-diagonal Qz electric operators
ops["Sz_r"] = qmb_op(in_ops, ["Sz", "Ig"])
ops["Sz_g"] = qmb_op(in_ops, ["Ir", "Sz"])
ops["Delta_Z"] = ops["Sz_g"] - ops["Sz_r"]
# Extra operators
ops["N_pair"] = ops["N_r"] @ ops["N_g"]
ops["N_single"] = ops["N_tot"] - 2 * ops["N_pair"]
ops["N_zero"] = ops["I"] - ops["N_single"] - ops["N_pair"]
return ops
[docs]
def SU2_gauge_invariant_states(
s_max, pure_theory, lattice_dim, background=0, n_flavors=None
):
"""Construct local SU(2) gauge-invariant basis states and basis matrices.
Parameters
----------
s_max : float
Maximum spin irrep used for the rishon Hilbert spaces.
pure_theory : bool
If ``True``, exclude matter fields.
lattice_dim : int
Number of spatial lattice dimensions.
background : int, optional
Maximum background-charge irrep included at the site.
n_flavors : int, optional
Number of matter flavors. Each flavor contributes one independent
SU(2) doublet on the site (4 fermion-occupation states). Defaults
to ``1`` for matter theories and is ignored when ``pure_theory=True``.
Returns
-------
tuple
``(gauge_basis, gauge_states)`` dictionaries for bulk and border site
classes. When ``background != 0``, exact background-resolved entries
keyed by ``(bg, label)`` are added alongside the legacy aggregate
``label`` keys.
"""
validate_parameters(
spin_list=[s_max], pure_theory=pure_theory, lattice_dim=lattice_dim
)
n_flavors = 0 if pure_theory else (1 if n_flavors is None else int(n_flavors))
spin_list = [S(s_max) for i in range(2 * lattice_dim)]
# Rishon-only outer product. Matter / background sectors are iterated
# explicitly below so that matter slots can be tagged with their
# per-flavor fermion occupation independently of the SU(2) labels.
rishon_spins = []
for s in spin_list:
tmp = np.arange(S(0), spin_space(s), 1)
rishon_spins.append(tmp / 2)
if background != 0:
bg_tmp = np.arange(S(0), spin_space(background), 1)
background_values = [S(irrep) for irrep in (bg_tmp / 2)]
else:
background_values = [S(0)]
# Total size of the full dressed-site basis (before GI projection).
j_list, _ = get_spin_Hilbert_spaces(
spin_list, pure_theory, background, n_flavors=n_flavors
)
site_dim = int(np.prod([len(j_set) for j_set in j_list], dtype=np.int64))
# Set rows and col counters list for the basis
gauge_states = {}
gauge_basis = {}
def register_basis_key(key):
gauge_states[key] = []
gauge_basis[key] = []
register_basis_key("site")
if background != 0:
for bg_value in background_values:
register_basis_key((bg_value, "site"))
# List of borders/corners of the lattice
borders = get_lattice_borders_labels(lattice_dim)
for label in borders:
border_key = f"site_{label}"
register_basis_key(border_key)
if background != 0:
for bg_value in background_values:
register_basis_key((bg_value, border_key))
# Outer iteration unit: one ``(occupations, bg_value, rishon_config)``
# triple per row. For pure theory we run a single dummy occupation
# tuple ``()`` and a single dummy bg sector. For a matter theory with
# ``N_f`` flavors, ``occupations`` cycles over ``product([0,1,2], …)``
# — the 3^N_f matter sectors — and feeds the right matter J's into
# the spin chain.
if pure_theory:
occupation_iter = [()]
else:
occupation_iter = list(product([0, 1, 2], repeat=n_flavors))
def _matter_js(occupations):
"""Physical SU(2) J for each matter flavor given its occupation."""
return [S(1) / 2 if n == 1 else S(0) for n in occupations]
def _process_row(args):
"""Per-row heavy compute: enumerate singlets and build canonical vectors.
Returns ``None`` when the row fails the cheap pre-checks (parity /
single-largest spin) or yields no singlets. Dict bookkeeping is
kept in the main thread below, so the output column ordering is
fully deterministic regardless of the number of worker threads.
"""
occupations, bg_value, rishon_cfg = args
matter_js = _matter_js(occupations)
spins_config = list(rishon_cfg)
if background != 0:
spins_config = [bg_value] + matter_js + spins_config
else:
spins_config = matter_js + spins_config
# 1) parity: total 2*j must be even
j2 = [int(2 * dd) for dd in spins_config]
total = sum(j2)
if total & 1:
return None
# 2) the largest spin cannot exceed the sum of the others
max_j2 = max(j2)
if max_j2 > total - max_j2:
return None
singlets_local = get_SU2_singlets(
spins_config,
pure_theory,
occupations if not pure_theory else None,
background,
)
if singlets_local is None:
return None
vectors = [
SU2_singlet_canonical_vector(spin_list, s, background)
for s in singlets_local
]
return spins_config, bg_value, singlets_local, vectors
# Thread-pool the per-row compute. The numba kernels are ``nogil=True``
# so they release the GIL while running. Honor ``NUMBA_NUM_THREADS`` for
# consistency with the validation suite ``--threads`` knob.
# NOTE: with the current leg-lookup cache the GIL-held Python work per
# row dominates; threading does not yet deliver a wall-clock speedup.
# See docs/su2_singlets_optimization/MILESTONES.md for the planned
# follow-up (defer SU2_singlet construction past the parallel section).
env_threads = os.environ.get("NUMBA_NUM_THREADS")
n_workers = int(env_threads) if env_threads else (os.cpu_count() or 1)
def _row_iter():
# Nesting order must match the legacy ``product(*spins)`` with
# ``spins = [bg, matter*N_f, rishons]``: bg slowest, then the matter
# occupation tuples, then rishons (innermost). Downstream code
# (e.g. ``DFL_Model.get_background_charges_configs``) relies on
# the resulting column ordering of the gauge basis.
for bg_value in background_values:
for occupations in occupation_iter:
for rishon_cfg in product(*rishon_spins):
yield occupations, bg_value, rishon_cfg
with ThreadPoolExecutor(max_workers=n_workers) as ex:
row_results = ex.map(_process_row, _row_iter())
for result in row_results:
if result is None:
continue
spins_config, bg_value, singlets_local, vectors = result
# Slot layout of ``spins_config``: [bg?, matter*N_f, rishons*2D].
# ``LGT_border_configs`` only looks at the rishon part — strip the
# bg and all matter slots before calling it and use
# ``pure_theory=True`` so it consumes the slice as-is.
n_skip = (1 if background != 0 else 0) + n_flavors
for s, singlet_state in zip(singlets_local, vectors):
gauge_states["site"].append(s)
gauge_basis["site"].append(singlet_state)
if background != 0:
exact_key = (bg_value, "site")
gauge_states[exact_key].append(s)
gauge_basis[exact_key].append(singlet_state)
spin_sizes = [spin_space(ss) for ss in spins_config[n_skip:]]
label = LGT_border_configs(
config=spin_sizes, offset=1, pure_theory=True
)
if label:
for ll in label:
border_key = f"site_{ll}"
gauge_states[border_key].append(s)
gauge_basis[border_key].append(singlet_state)
if background != 0:
exact_key = (bg_value, border_key)
gauge_states[exact_key].append(s)
gauge_basis[exact_key].append(singlet_state)
# Build the basis combining the states into a matrix
for label in list(gauge_basis.keys()):
if gauge_basis[label]:
gauge_basis[label] = csr_matrix(np.column_stack(tuple(gauge_basis[label])))
else:
gauge_basis[label] = csr_matrix((site_dim, 0))
return gauge_basis, gauge_states
[docs]
def SU2_gauge_invariant_operators(
spin,
pure_theory,
lattice_dim,
background=0,
use_generic_operators=False,
):
"""Project truncated SU(2) dressed-site operators onto GI local bases.
Parameters
----------
spin : float
Gauge-link spin representation.
pure_theory : bool
If ``True``, exclude matter fields.
lattice_dim : int
Number of spatial dimensions.
background : scalar, optional
Maximum background irrep included in the local basis.
use_generic_operators : bool, optional
If ``True``, force the generalized truncated operator construction.
Otherwise the hardcoded construction is used for ``spin < 1`` and the
generalized one for higher truncations.
Returns
-------
dict
Dictionary keyed by ``(label, op_name)``. When ``background == 0``,
``label`` is a string such as ``"site"`` or ``"site_mx"``. When
``background != 0``, ``label`` is an exact basis label of the form
``(bg, geom_label)``.
"""
if spin < 1 and not use_generic_operators:
in_ops = SU2_dressed_site_operators(
spin, pure_theory, lattice_dim, background=background
)
else:
in_ops = SU2_gen_dressed_site_operators(
spin, pure_theory, lattice_dim, background=background
)
gauge_basis, _ = SU2_gauge_invariant_states(
spin, pure_theory, lattice_dim, background
)
use_exact_bg_labels = background != 0
ops = {}
for label, basis in gauge_basis.items():
if use_exact_bg_labels != isinstance(label, tuple):
continue
if basis.shape[1] == 0:
continue
for op_name, op in in_ops.items():
ops[(label, op_name)] = basis.transpose() @ op @ basis
return ops
[docs]
def SU2_check_gauss_law(
gauge_basis: csr_matrix, gauss_law_op: csr_matrix, threshold=1e-14
):
"""Validate an SU(2) gauge-invariant basis against a Gauss-law operator.
Parameters
----------
gauge_basis : scipy.sparse.csr_matrix
Basis matrix whose columns span the candidate gauge-invariant subspace.
gauss_law_op : scipy.sparse.csr_matrix
Local Gauss-law operator.
threshold : float, optional
Numerical tolerance used in the checks.
Returns
-------
None
Raises
------
TypeError
If inputs are not sparse matrices.
ValueError
If the basis is not isometric/projective or does not lie in the
Gauss-law kernel.
"""
if not isspmatrix(gauge_basis):
raise TypeError(f"gauge_basis should be csr_matrix, not {type(gauge_basis)}")
if not isspmatrix(gauss_law_op):
raise TypeError(f"gauss_law_op shoul be csr_matrix, not {type(gauss_law_op)}")
# This functions performs some checks on the SU2 gauge invariant basis
logger.info("CHECK GAUSS LAW")
# True and the Effective dimensions of the gauge invariant dressed site
site_dim = gauge_basis.shape[0]
eff_site_dim = gauge_basis.shape[1]
# Check that the Matrix Basis behave like an isometry
norm_isometry = norm(gauge_basis.transpose() * gauge_basis - identity(eff_site_dim))
if norm_isometry > threshold:
raise ValueError(f"Basis must be Isometry: B^T*B=1; got norm {norm_isometry}")
# Check that B*B^T is a Projector
Proj = gauge_basis * gauge_basis.transpose()
norm_Proj = norm(Proj - Proj**2)
if norm_Proj > threshold:
raise ValueError(f"P=B*B^T: expected P-P**2=0: obtained norm {norm_Proj}")
# Check that the basis is the one with ALL the states satisfying Gauss law
norma_kernel = norm(gauss_law_op * gauge_basis)
if norma_kernel > threshold:
raise ValueError(f"Gauss Law Kernel with norm {norma_kernel}; expected 0")
GL_rank = matrix_rank(gauss_law_op.todense())
if site_dim - GL_rank != eff_site_dim:
logger.info("Large dimension %s", site_dim)
logger.info("Effective dimension %s", eff_site_dim)
logger.info(GL_rank)
logger.info("Some gauge basis states are missing")
logger.info("GAUSS LAW SATISFIED")