January 28, 2026
Generalized global symmetries, and in particular higher-form center symmetries of Yang–Mills (YM) theory, provide a sharp language to characterize phases and to formulate anomaly constraints. On the lattice, coupling the \(\mathbb{Z}_N\) 1-form center symmetry to a background 2-form gauge field \(B\) is equivalent to imposing ’t Hooft twisted boundary conditions, labeled by discrete fluxes \(z_{\mu\nu}\in\mathbb{Z}_N\). A crucial consequence is the fractionality of the topological charge in a background \(B\), which leads to a mixed ’t Hooft anomaly between the center symmetry and the \(2\pi\) periodicity of the \(\theta\) angle. These structures suggest that partition functions in fixed flux sectors, and their discrete Fourier transforms, are natural order parameters for confinement and its variants.
In this proceedings contribution, we summarize results from two companion papers [1], [2] and recent progress. Our main technical contribution is a “halfway-updating” hybrid Monte Carlo (HMC) algorithm that samples both link variables and discrete background flux fields, enabling direct estimation of the partition function \(Z[B]\) and the ’t Hooft partition function \(Z_{\mathrm{tH}}\). This approach provides information on (de)confinement, Higgs, Coulomb and oblique-confining phases.
We consider \(SU(N)\) lattice YM theory on a four-torus. The \(\mathbb{Z}_N\) 1-form center symmetry acts on link variables \(U_\ell\in SU(N)\) by phases determined by intersection numbers with a codimension-2 surface \(\Sigma\). Coupling to a background 2-form gauge field \(B_p\in\mathbb{Z}_N\) modifies plaquette terms schematically as \[S[U,B] \sim \sum_{p} \mathrm{Re} \mathop{\mathrm{tr}}\left(e^{-2\pi i B_p/N} U_p\right), \label{eq:coupleB}\tag{1}\] and is invariant under the combined transformations of link phases \(U\to e^{2\pi i\lambda/N} U\) and \(B\to B + d\lambda\). Globally, this background is encoded by ’t Hooft twisted boundary conditions with fluxes \(z_{\mu\nu}=\sum B_{\mu\nu}\bmod N\) [3].1
In the presence of nontrivial flux, the topological charge can become fractional. We generalize the partition function by multiplying the \(\theta\)-term \(e^{-i\theta Q}\). For example, for twisted boundary conditions one has [4] \[Q = \frac{1}{16\pi^2}\int \mathop{\mathrm{tr}}F\tilde F = -\frac{\varepsilon_{\mu\nu\rho\sigma} z_{\mu\nu} z_{\rho\sigma}}{8N} + \mathbb{Z} \in \frac{1}{N}\mathbb{Z}, \label{eq:fractionalQ}\tag{2}\] which implies that the partition function in background \(B\) is not strictly \(2\pi\) periodic in \(\theta\) but instead obeys \[Z_{\theta+2\pi}[B] = e^{-2\pi i Q[B]} Z_{\theta}[B]. \label{eq:thetaAnomaly}\tag{3}\] The explicit form (or local description) of \(Q[B]\in\frac{1}{N}\mathbb{Z}\) and the existence of the mixed ’t Hooft anomaly even on the lattice were proved in Refs. [5] based on Ref. [6].
This mixed anomaly underlies the Witten effect relation for the ’t Hooft partition function discussed below.
Let \(Z[B]\) denote the YM partition function in a fixed background \(B_{\mu\nu}\in\mathbb{Z}_N\). The ’t Hooft partition function is defined as a discrete Fourier transform over temporal flux components, \[Z_{\mathrm{tH}}[E_i;B_{ij}] \equiv \frac{1}{N^3}\sum_{B_{i4}=0}^{N-1} \exp\left(\frac{2\pi i}{N}\sum_{i=1}^3 E_i B_{i4}\right) Z[B], \label{eq:ZtH}\tag{4}\] where \(E_i\in\mathbb{Z}_N\) can be viewed as discrete “electric” fluxes and \(B_{ij}\) as fixed “magnetic” fluxes. \(Z_{\mathrm{tH}}\) is expected to detect distinct quantum phases (confinement/Higgs/Coulomb and oblique variants) and is closely related to the Wilson–’t Hooft classification of line operators [7]–[9].
Traditionally, the free energy differences between flux sectors have been estimated by reweighting [10], [11], say \(Z[z\neq0]/Z[z=0]\) in terms of the fluxes. Our goal is to measure \(Z[B]\) and hence \(Z_{\mathrm{tH}}\) directly, by constructing a Markov chain whose stationary distribution includes the discrete flux variables. To do this, for simplicity, we take gauge transformations to set \[\begin{align} B_{\mu\nu}(x) = \begin{cases} z_{\mu\nu} & \text{x_\mu=L-1 and x_\nu=L-1}\\ 0 & \text{otherwise} . \end{cases} \end{align}\] This is a representative of an equivalence class, and there are \(N^6\) patterns of different classes of \(B\).
We outline the algorithmic idea for the joint sampling of link variables \(U\) and discrete background fields \(B\). Starting from a configuration \((U,B)\), we:
Generate conjugate momenta \(\pi\) for \(U\) from a Gaussian distribution.
Evolve \((U,\pi)\) under molecular dynamics (MD) for half a trajectory length \(\tau/2\) using the Hamiltonian \(H(U,\pi,B)=\pi^2/2 + S[U,B]\), reaching \((\Check U,\Check\pi)\).
Update the discrete field by a proposal kernel \(P_F(B\to B')\) that is symmetric, \(P_F(B\to B')=P_F(B'\to B)\).
Continue the MD evolution for another half step \(\tau/2\) using \(H(\Check U,\Check\pi,B')\), reaching \((U',\pi')\).
Accept/reject the composite move \((U,B)\to(U',B')\) with the Metropolis acceptance probability \(\min\{1,\exp(-\Delta H)\}\), where \(\Delta H=H(U',\pi',B')-H(U,\pi,B)\).
The symmetry of \(P_F\) ensures detailed balance for the joint Boltzmann distribution. A uniform proposal over flux sectors—\(SU(N)/\mathbb{Z}_N\) YM theory with uniformly random \(P_F\)—is particularly effective in reducing autocorrelation by mixing boundary conditions [1].
Because \(B\) is part of the dynamical Markov chain, we can estimate the normalized weight of each flux sector by counting occurrences: \[\widehat Z[B] \equiv \frac{N_{\mathrm{conf}}(B)}{N_{\mathrm{total}}}\,. \label{eq:counting}\tag{5}\] This estimator corresponds to \(Z[B]/\sum_{B}Z[B]\) (i.e., the probability of sector \(B\) in the extended ensemble). On \(T^4\), we additionally average over Euclidean rotations of the flux tuple \[(B_{12},B_{13},B_{14},B_{23},B_{24},B_{34})\] when appropriate to improve statistics. In the \(T^4\) study reported in Ref. [2], we used \(N_{\mathrm{total}}=2590\) configurations of \(SU(2)\) YM at \(\beta=2.6\) and the box size \(L_x\times L_y\times L_z\times L_t=20\times20\times20\times20\).
Having estimated \(\widehat Z[B]\) for all relevant \(B\), we obtain \(Z_{\mathrm{tH}}\) from Eq. 4 (up to an overall normalization that cancels in ratios such as \(Z_{\mathrm{tH}}/Z_{\mathrm{tH}}[0;0]\)), without reweighting.
We focus on \(SU(2)\) lattice YM on \(T^4\) and measure \(Z_{\mathrm{tH}}[E;B]\) for all combinations of electric fluxes \(E_i\in\{0,1\}\) and magnetic fluxes \(B_{ij}\in\{0,1\}\). Figure 1 shows representative data in a confining regime. \(x\)-axis enumerates \([E_1,E_2;B_{23},B_{31}]\) with \(E_3\), \(B_{12}\) fixed as in legend. The filled symbols denote \(Z_{\mathrm{tH}}\) itself, while the unfilled ones represent the conterparts of the duality equation.
The main qualitative observation is the “light/heavy” hierarchy: \[\frac{Z_{\mathrm{tH}}[E=0;B]}{Z_{\mathrm{tH}}[E=0;B=0]}\sim 1, \qquad \frac{Z_{\mathrm{tH}}[E\neq 0;B]}{Z_{\mathrm{tH}}[E=0;B=0]}\sim 0, \label{eq:lightheavy}\tag{6}\] consistent with confinement where nontrivial electric flux costs an area-law free energy.
Figure 2 summarizes the averaged data and the observed pattern with less numerical errors. Unfortunately, while the asymptotic behavior in the large volume limit is obtained by \[\frac{Z_{\mathrm{tH}}[E=0;B]}{Z_{\mathrm{tH}}[E=0;B=0]} - 1=O(e^{-\sigma L}) \qquad \frac{Z_{\mathrm{tH}}[E\neq 0;B]}{Z_{\mathrm{tH}}[E=0;B=0]}=O(e^{-\sigma L}), \label{eq:tension}\tag{7}\] where \(\sigma\) is the (dual) string tension, it would be prohibitively expensive to determine \(\sigma\) with the present statistics.
A further diagnostic is the Witten effect at \(\theta=2\pi\), which implies a shift of electric flux by magnetic flux. The ’t Hooft partition function with the \(\theta\)-term is given by Eq. 4 by replacing \(Z[B]\) with \(Z_\theta[B]\). In our normalization this appears as \[Z_{\mathrm{tH},\theta=2\pi}[E_i;B_{jk}] = Z_{\mathrm{tH},\theta=0}[E_i + B_{jk}, B_{jk}], \label{eq:witten}\tag{8}\] illustrating oblique confinement as in Fig. 3: the \(\theta\) term effectively transmutes the “light” sector.
For thermal geometries \(T^3\times S^1\) near the deconfinement transition (small \(S^1\)), flux-sector physics becomes more delicate. So far, we assumed that in our algorithm the total partition function for \(SU(N)/\mathbb{Z}_N\) can be decomposed by \[Z_{SU(N)/\mathbb{Z}_N} = \sum_B Z[B],\qquad N_{\mathrm{total}} = \sum_B N_{\mathrm{conf}}(B).\] Is \(Z[B]\) truly a partition function? If so, configurations in \(Z[B]\) should be sufficiently thermalized in the conventional HMC for fixed \(B\). When \(P_F(B\to B')\) generates a different \(B_{i4}'\neq B_{i4}\), a thermalized configuration \(U'\) is far from the initial configuration \(U\) in \(Z[B]\). Figure 4 shows the distributions of values of plaquette and action for zero- and finite-temperature, illustrating the corresponding Boltzmann distribution. Only \(B_{34}\) can be taken as \(0\) or \(1\), otherwise \(B_{\mu\nu}=0\). For zero-temperature, the distributions for \(B_{34}=0\), \(1\) overlap each other, but for finite-temperature (when \(S^1\) becomes smaller), those move away and stay farther apart.
An important practical question is whether \(Z[B]\) computed by counting is robust against different thermalization times for different \(B\) sectors. This is tied to the “separability” of \(\sum_B Z[B]\) and to constructing an efficient proposal kernel \(P_F(B\to B')\) that respects sector-dependent thermalization scales.
Diagnosing (non-)thermalization and sector-dependent bias at finite temperature. Near the deconfinement transition on \(T^3\times S^1\), the Markov chain may sample different flux sectors with substantially different relaxation times. In that situation, the naive counting estimator in Eq. 5 can be biased if measurements are taken before the within-sector dynamics has equilibrated after a change of \(B\). To make this issue quantitative, we monitor standard observables such as the plaquette and the gauge action, and compare their distributions in each sector with those obtained from a conventional HMC run at fixed \(B\) (matched at the same bare parameters and geometry). Figure 4 illustrates these diagnostics by comparing the plaquette and action distributions between flux sectors at zero and finite temperature. As a compact diagnostic, we also track the integrated autocorrelation time \(\tau_{\mathrm{int}}\) of the plaquette within each sector, and the acceptance rate of the composite move \((U,B)\to(U',B')\) as functions of \((B,\beta,L,T)\).
A practical measurement protocol. In the thermal runs, after every accepted update that changes the temporal flux components \(B_{i4}\), we discard \(N_{\mathrm{therm}}\) subsequent trajectories (“sector re-thermalization”) before recording measurements. Equivalently, we may record measurements only when the chain has remained in the same sector for at least \(N_{\mathrm{therm}}\) trajectories. We then verify stability by checking that the sector-resolved plaquette/action histograms agree (within statistics) between the first and second halves of the Monte Carlo history, and that the inferred \(\widehat Z[B]\) is insensitive to moderate variations of \(N_{\mathrm{therm}}\). These checks provide a direct, observable-based criterion for whether sector counting can be interpreted as a reliable proxy for relative weights at finite temperature.
Tuning the flux proposal. Finally, the proposal kernel \(P_F(B\to B')\) should be tuned so that sector changes are neither too rare (leading to poor exploration) nor too frequent (leading to persistent re-thermalization transients). In practice, we adjust the proposal width (e.g.the parameter \(\alpha\) in the Gaussian proposal used below) to obtain a stable acceptance rate and manageable \(\tau_{\mathrm{int}}\), and we report these algorithmic diagnostics together with the physics observables.
As a first test, a simplified setup is that only \(B_{34}=z\in\{0,1\}\) is turned on and employed a Gaussian proposal \(P_F\propto e^{-\alpha(z-z')^2}\). We should find an empirical thermalization scale of MD time units to tune the parameter \(\alpha\). Notably, the thermalization time/acceptance ratio in MD starting from \(U[B_{34}=0]\) to \(U[B_{34}=1]\) is quite different from that from \(U[B_{34}=1]\) to \(U[B_{34}=0]\). The practical transition probability \(P_F\) in MD time units is larger than the longest thermalization time for each lattice parameter.
Based on our experience, the observation of \(Z_{\mathrm{tH}}\) is hindered by too large systematic error of order \(Z_{\mathrm{tH}}[0;0]/3\), while statistical errors are comparable with those in Fig. 2 .
Improving the finite-temperature analysis will require (i) larger statistics for finite-size scaling, (ii) systematic control of thermalization in each flux sector, and (iii) optimized \(P_F\) kernels tailored to the sector structure. It will also be interesting to compare our lattice estimates of \(Z[B]\) and \(Z_{\mathrm{tH}}\) with continuum predictions from symmetry-TFT/anomaly-inflow descriptions of center symmetry, in particular the phase factors and selection rules of twisted partition functions.
This work was partially supported by Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for Scientific Research Grant Number JP25K17402 (O.M.) and JP23K03418 (H.S.). The numerical computations in this paper were carried out on Genkai, a
supercomputer system of the Research Institute for Information Technology (RIIT), Kyushu University. Our numerical codes, which can be found in https://github.com/o-morikawa/Gaugefields.jl, is an extended version of Gaugefields.jl in JuliaQCD project [12].
O.M.acknowledges the RIKEN Special Postdoctoral Researcher Program and RIKEN FY2025 Incentive Research Projects.
In what follows, we assume that the ’t Hooft fluxes \(z\) exist, which means that the background 2-form gauge field is flat: \(dB=0\bmod N\).↩︎