1 Introduction↩︎

Standard lattice simulations of four-dimensional Yang–Mills theories sample gauge fields \(U\) from the target Boltzmann weight \(p(U)\propto \mathrm{e}^{-S[U]}\) using Markov Chain Monte Carlo (MCMC) methods. As the continuum limit is approached, autocorrelations grow because of critical slowing down, and the slowest modes are typically topological. This leads to a severe loss of ergodicity, known as topological freezing [1]–[3]. A robust mechanism to accelerate topology sampling is the use of open boundary conditions (OBC) in time, which remove topological barriers and turn the evolution of topological modes into a diffusion process [4]–[6]. However, OBC break translation invariance and introduce unphysical boundary effects: the key algorithmic question is therefore whether it is possible to exploit the fast Monte Carlo dynamics of topological fluctuations of OBC while computing observables with periodic boundary conditions (PBC).

A successful answer to this problem is provided by Parallel Tempering on Boundary Conditions (PTBC), where replicas with different boundary conditions interpolating between OBC and PBC are simulated in parallel and allowed to swap configurations [7]–[16]. In this way, fast topological fluctuations generated with OBC are transferred to the physical PBC ensemble, while maintaining exactness. Our approach [17] shares the same philosophy, but follows a different, flow-based route.

In this conference proceeding we focus on a flow-based strategy based on Stochastic Normalizing Flows (SNFs) [17]–[22], a combination of Normalizing Flows [23], [24] and Non-Equilibrium Markov Chain Monte Carlo (NE-MCMC) calculations based on Jarzynski and Crooks identities [25]–[32]. Crucially, SNFs preserve the favorable scaling of NE-MCMC with the number of degrees of freedom undergoing a transformation, while significantly reducing the dissipated work compared to purely stochastic protocols, translating into a substantial gain in efficiency. In particular, we exploit the NE-MCMC structure to design a controlled non-equilibrium interpolation between OBC and PBC, already explored in Ref. [32] in \(2d\) \(\mathrm{CP}^{N-1}\) models and in Ref. [17] for \(4d\) \(\mathrm{SU}(3)\) Yang–Mills theory, that preserves exactness through reweighting and improves the sampling of topological sectors near the continuum limit.

More broadly, this work is part of a recent effort to go beyond purely local Monte Carlo updates by constructing learned transformations between ensembles in lattice field theory simulations. Normalizing flows have been widely investigated in lattice gauge theory [33]–[43] and realize exact, invertible maps between probability distributions, while related approaches based on diffusion models [44]–[50], generate configurations by reversing a stochastic noise process.

2 Stochastic Normalizing Flows↩︎

We consider the problem of transporting gauge configurations from a prior ensemble (that is easier to sample from) to a physical target one. Concretely, we introduce a prior distribution \(q_0(U_0)\propto \mathrm{e}^{-S_0[U_0]}\) and a target distribution \(p(U)\propto \mathrm{e}^{-S[U]}\) (here, OBC and PBC in \(4d\) \(\mathrm{SU}(3)\) gauge theory). The transport is implemented through a non-equilibrium evolution, defined as a sequence \(\mathcal{U}=[U_0,\dots,U_{n_{\mathrm{step}}}]\equiv[U_0,\dots,U]\) generated by a protocol \(\lambda(n)\) via Markov updates \(P_{\lambda(n)}\) with equilibrium weight \(\propto \mathrm{e}^{-S_{\lambda(n)}}\). A Jarzynski-based reweighting yields an unbiased estimator for target expectation values: \[\langle \mathcal{O} \rangle_p = \frac{\langle \mathcal{O}(U)\, \mathrm{e}^{-W(\mathcal{U})} \rangle_{\mathrm{f}}}{\langle \mathrm{e}^{-W(\mathcal{U})} \rangle_{\mathrm{f}}}\] where \(\langle \cdot \rangle_{\mathrm{f}}\) is an average over forward trajectories and the work \(W\) is \[W(\mathcal{U})=\sum_{n=0}^{n_{\mathrm{step}}-1}\Bigl(S_{\lambda(n+1)}[U_n]-S_{\lambda(n)}[U_n]\Bigr)\,.\] In practice, the exponential reweighting becomes inefficient if the dissipated work \(W_{\mathrm{d}}=W-\Delta F\), with \(\Delta F\) being the free energy difference between prior and target, is large.

SNFs enhance these trajectories by inserting parametric, deterministic, invertible layers \(g_{\rho(n)}\) between stochastic updates: \[U_0 \xrightarrow{\, g_{\rho(1)} \,} g_{\rho(1)}(U_0) \xrightarrow{\, P_{\lambda(1)} \,} U_1 \xrightarrow{\, g_{\rho(2)} \,} \cdots \xrightarrow{\, P_{\lambda(n_{\mathrm{step}})} \,} U_{n_{\mathrm{step}}}\equiv U\,.\] The correct reweighting uses a variational work including Jacobians [18], [19], [51]: \[\begin{align} W^{(\rho)}(\mathcal{U}) &= S[U]-S_0[U_0]-Q^{(\rho)}(\mathcal{U}) -\sum_{n=0}^{n_{\mathrm{step}}-1}\log\left|\det J_{g_{\rho(n+1)}}(U_n)\right|\,, \end{align}\] where \(Q^{(\rho)}(\mathcal{U})\) is a generalized pseudo-heat term (including the effect of inserting \(g_{\rho}\) layers).

A key quantity to control the performances of non-equilibrium samplers is the reverse Kullback–Leibler (KL) divergence between forward and reverse path probability densities: \[\tilde{D}_{\mathrm{KL}}\bigl(q_0\,\mathcal{P}_{\mathrm{f}}\,\|\, p\,\mathcal{P}_{\mathrm{r}}\bigr)=\langle W_{\mathrm{d}}\rangle_{\mathrm{f}} \ge 0\,,\] so making trajectories more reversible (smaller \(\langle W_{\mathrm{d}}\rangle\)) directly stabilizes reweighting. A commonly used empirical proxy for the efficiency of the reweighting estimator, expressed as an effective sample size, is \[\hat{\mathrm{ESS}}\equiv \frac{\langle \mathrm{e}^{-W}\rangle_{\mathrm{f}}^2}{\langle \mathrm{e}^{-2W}\rangle_{\mathrm{f}}} = \frac{1}{\langle \mathrm{e}^{-2W_{\mathrm{d}}}\rangle_{\mathrm{f}}}\,,\] showing that the effective sample size is controlled by fluctuations of the dissipated work, and that large positive values of \(W_{\mathrm{d}}\) exponentially suppress \(\hat{\mathrm{ESS}}\).

2.1 Defect gauge-equivariant coupling layers↩︎

To construct efficient layers for gauge theories, we use gauge equivariant updates based on masked stout smearing [37], [52], [53]. At layer \(n\), a link is updated as \[U'_\mu(x)=\exp\,\bigl(i\,Q^{(n)}_\mu(x)\bigr)\,U_\mu(x)\,,\] with \(Q^{(n)}_\mu(x)\) traceless Hermitian, built from staples through \[\require{physics} \Omega^{(n)}_\mu(x)=C^{(n)}_\mu(x)\,U^\dagger_\mu(x)\,, \qquad Q^{(n)}_\mu(x)=\frac{i}{2}\bigl(\Omega^{\dagger}-\Omega\bigr)-\frac{i}{2N}\Tr\bigl(\Omega^{\dagger}-\Omega\bigr)\,.\] The staple sum is \[\begin{align} C^{(n)}_\mu(x)=\sum_{\nu\neq\mu}\Bigl[ &\rho^{+}_{\mu\nu}(n,x)\,U_\nu(x)U_\mu(x+\hat{\nu})U^\dagger_\nu(x+\hat{\mu})\\ +&\rho^{-}_{\mu\nu}(n,x)\,U^\dagger_\nu(x-\hat{\nu})U_\mu(x-\hat{\nu})U_\nu(x-\hat{\nu}+\hat{\mu}) \Bigr]\,. \end{align}\] To make the Jacobian determinant tractable, we use an even–odd masking schedule, which ensures a (block-)triangular Jacobian and thus allows an efficient computation of the determinant as the product of the diagonal blocks.

In boundary-condition evolutions the action is modified only in a localized region (the defect), thus, following the work in Ref. [54], we restrict the deterministic layer support to the defect and its immediate neighborhood. Refer to Fig. 1 for a schematic representation of the defect coupling layer. This particular coupling layer reduces the number of trainable parameters and focuses the flow capacity where the mismatch to PBC is largest. Global propagation of the defect information is still guaranteed by the interleaved stochastic gauge update, which acts on the full lattice. In this work, we use the standard heatbath plus 4 overrelaxation steps as MCMC update.

Defect SNFs training is done by minimizing the average dissipated work, \[\mathcal{L}(\rho)=\langle W_{\mathrm{d}}^{(\rho)} \rangle_{\mathrm{f}},\] with respect to the parameters \(\rho\) of the layer. This procedure favors reversible trajectories and improves both variance and stability of reweighting. Furthermore, as shown in Ref. [17], [21], in practice, one can train at fixed target coupling and small \(n_{\mathrm{step}}\) (e.g.\(n_{\mathrm{step}}=8,16\)) and observe smooth profiles as functions of \(n/n_{\mathrm{step}}\). This enables a simple transfer procedure: group parameters into a small number of geometric classes (suggested by symmetries of the defect), and interpolate with splines a rescaled profile \(\rho_{\mathrm{class}}(n)\,n_{\mathrm{step}}\) as a function of \(n/n_{\mathrm{step}}\). The resulting fit is then used to instantiate larger \(n_{\mathrm{step}}\) flows at negligible additional training cost. See Ref. [17], [21] for further details on the training procedure.

3 Scaling and numerical performance↩︎

For boundary-condition flows, the relevant size parameter is the number of degrees of freedom affected by the defect, \(n_{\mathrm{dof}}\propto (L_d/a)^3\). A robust empirical scaling observed for purely stochastic non-equilibrium flows is \[\langle W_{\mathrm{d}}\rangle_{\mathrm{f}} \propto \frac{n_{\mathrm{dof}}}{n_{\mathrm{step}}}\,,\] and \[\hat{\mathrm{ESS}}\approx \exp\,\Bigl(-k'\,\frac{n_{\mathrm{dof}}}{n_{\mathrm{step}}}\Bigr)\,,\] so controlling \(\hat{\mathrm{ESS}}\) at fixed defect size requires \(n_{\mathrm{step}}\propto n_{\mathrm{dof}}\). These scaling relations are illustrated in Fig. 2, where we compare purely stochastic NE-MCMC and defect SNFs at \(\beta=6.0\) on a \(L/a=16\) lattice for several defect sizes. When the performance metrics are plotted against the scaling variable \(n_{\mathrm{step}}/n_{\mathrm{dof}}\), data corresponding to different \(L_d/a\) values collapse onto a common curve, showing that the efficiency is primarily controlled by the ratio between the number of non-equilibrium steps and the number of degrees of freedom touched by the defect. Increasing \(n_{\mathrm{step}}/n_{\mathrm{dof}}\) leads to smaller \(\langle W_{\mathrm{d}}\rangle_{\mathrm f}\) (left panel) and, consistently, to larger \(\hat{\mathrm{ESS}}\) (right panel), i.e., to a better-behaved reweighting estimator. At fixed \(n_{\mathrm{step}}/n_{\mathrm{dof}}\), defect SNFs systematically yield more reversible trajectories than NE-MCMC, resulting in a higher \(\hat{\mathrm{ESS}}\); equivalently, for a fixed target \(\hat{\mathrm{ESS}}\) the SNF reaches the same estimator quality with fewer non-equilibrium steps, corresponding to an overall speedup of about a factor \(\sim 3\) in this setup.

As a physics validation of the method, Fig. 3 shows the topological susceptibility in lattice units \(a^{4}\chi_{_{\scriptscriptstyle{\rm L}}}\) obtained from the reweighted PBC ensemble as a function of the proxy effective sample size \(\hat{\mathrm{ESS}}\). Results at \(\beta=6.4\) on a \(30^{4}\) lattice and at \(\beta=6.5\) on a \(34^{4}\) lattice are consistent across different flow setups and defect sizes. Moreover, our determinations agree with high-statistics reference computations, shown as horizontal bands, providing a non-trivial check that the reweigthing procedure suggested by Jarzynski equality correctly removes the defect-induced boundary artifacts and reproduces the theory with PBC.

4 Conclusions and future outlooks↩︎

We have introduced a defect-based Stochastic Normalizing Flow strategy to exploit the fast topological dynamics of open boundaries while recovering expectation values of the physical periodic theory through Jarzynski equality. By interleaving global non-equilibrium MCMC steps with localized, gauge-equivariant defect coupling layers, the SNF produces more reversible trajectories than purely stochastic NE-MCMC, reducing dissipation and increasing \(\hat{\mathrm{ESS}}\) at fixed scaling variable \(n_{\mathrm{step}}/n_{\mathrm{dof}}\); this yields a tangible speedup at fixed estimator quality.

Future developments include several clear directions. On the machine learning side, richer gauge-equivariant layers and multiscale architectures [56]–[58] can improve the gain of SNFs on NE-MCMC. On the algorithmic side, improving the non-equilibrium evolution itself is a natural direction, in particular by optimizing the schedule for the boundary-interpolation parameter (beyond a linear protocol) to reduce the dissipated work at fixed computational cost. Finally, the same framework can be extended to more challenging settings (including dynamical-fermion simulations), with the potential to enable controlled topology sampling closer to the continuum limit.

Acknowledgments↩︎

We thank M. Caselle, G. Kanwar and M. Panero for insightful and helpful discussions. C. B. acknowledges support by the Spanish Research Agency (Agencia Estatal de Investigación) through the grant IFT Centro de Excelencia Severo Ochoa CEX2020-001007-S and, partially, by grant PID2021-127526NB-I00, both funded by MCIN/AEI/10.13039/ 501100011033. A. B., A. N., D. P. and L. V. acknowledge support by the Simons Foundation grant 994300 (Simons Collaboration on Confinement and QCD Strings). A. B. was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) as part of the CRC 1639 NuMeriQS – project no. 511713970 and under Germany’s Excellence Strategy – Cluster of Excellence Matter and Light for Quantum Computing (ML4Q) EXC 2004/1 – 390534769. A. N. acknowledges support from the European Union - Next Generation EU, Mission 4 Component 1, CUP D53D23002970006, under the Italian PRIN “Progetti di Ricerca di Rilevante Interesse Nazionale – Bando 2022” prot. 2022ZTPK4E. A. B., A. N., D. P. and L. V. acknowledge support from the SFT Scientific Initiative of INFN. The work of D. V. is supported by STFC under Consolidated Grant No. ST/X000680/1. We acknowledge EuroHPC Joint Undertaking for awarding the project ID EHPC-DEV-2024D11-010 access to the LEONARDO Supercomputer hosted by the Consorzio Interuniversitario per il Calcolo Automatico dell’Italia Nord Orientale (CINECA), Italy. This work was partially carried out using the computational facilities of the “Lovelace” High Performance Computing Centre, University of Plymouth.

References↩︎