# Learning nonequilibrium statistical mechanics and dynamical phase
transitions

Ying Tang<EMAIL_ADDRESS>Institute of Fundamental and Frontier
Sciences, University of Electronic Sciences and Technology of China, Chengdu
611731, China  International Academic Center of Complex Systems, Beijing
Normal University, Zhuhai 519087, China  Jing Liu  CAS Key Laboratory for
Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of
Sciences, Beijing 100190, China  School of Systems Science, Beijing Normal
University, Beijing 100875, China  Jiang Zhang  School of Systems Science,
Beijing Normal University, Beijing 100875, China  Swarma Research, Beijing
102308, China  Pan Zhang<EMAIL_ADDRESS>CAS Key Laboratory for
Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of
Sciences, Beijing 100190, China  School of Fundamental Physics and
Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou
310024, China  Hefei National Laboratory, Hefei 230088, China

###### Abstract

Nonequilibrium statistical mechanics exhibit a variety of complex phenomena
far from equilibrium. It inherits challenges of equilibrium, including
accurately describing the joint distribution of a large number of
configurations, and also poses new challenges as the distribution evolves over
time. Characterizing dynamical phase transitions as an emergent behavior
further requires tracking nonequilibrium systems under a control parameter.
While a number of methods have been proposed, such as tensor networks for one-
dimensional lattices, we lack a method for arbitrary time beyond the steady
state and for higher dimensions. Here, we develop a general computational
framework to study the time evolution of nonequilibrium systems in statistical
mechanics by leveraging variational autoregressive networks, which offer an
efficient computation on the dynamical partition function, a central quantity
for discovering the phase transition. We apply the approach to prototype
models of nonequilibrium statistical mechanics, including the kinetically
constrained models of structural glasses up to three dimensions. The approach
uncovers the active-inactive phase transition of spin flips, the dynamical
phase diagram, as well as new scaling relations. The result highlights the
potential of machine learning dynamical phase transitions in nonequilibrium
systems.

##  I Introduction

Tracking time evolution and characterizing phase transitions are fundamental
tasks in nonequilibrium statistical mechanics [1], having implications to a
wide range of fields including quantum transport [2], molecular machines [3],
quantitative biology [4, 5], and complex networks [6]. For example, the glassy
behavior has been explored in kinetically constrained models (KCM) of spin
flips on a lattice [7], where each spin can facilitate the flip of its
neighbor spins. The stochastic dynamics of spin flips give insights into the
dynamical heterogeneity of the glass transition [8]. Another example is the
voter model [9] describing consensus formation, where voters located on a
network choose the opinion based on their neighbors and can form
characteristic spatial structures.

Despite tremendous efforts, studying the dynamical phase transition remains
challenging in general. First, it requires tracking the evolving probability
distribution of microscopic configurations, which is exponentially increasing
with the system size and can be computationally prohibitive from sampling
trajectories [10]. Second, the phase transition is induced by a control
parameter, and the corresponding dynamical operator governing the time
evolution, known as the “tilted” generator [11], does not preserve the
normalization condition of the distribution. Thus, one needs to estimate the
normalization factor, i.e., the dynamical partition function, a central
quantity in nonequilibrium systems [12]. Third, studying nonequilibrium
dynamics at an arbitrary time is even more challenging, because it demands
estimating the whole spectrum of the tilted generator, rather than computing
only the largest eigenvalue in the long-time limit [13, 14] based on the large
deviation theory [11, 15]. Analytically, it is intractable except for rare
cases [16, 17, 18]; numerically, the major difficulty of studying time
evolution beyond the steady state arises from expressing the high-dimensional
probability distribution with a growing complexity when the correlation builds
up over time.

In this work, we develop a general framework to track nonequilibrium systems
over time by variational autoregressive networks (VAN), which offers an ideal
model for describing the normalized joint distribution of configurations [19].
The VAN was previously shown effective to investigate equilibrium statistical
mechanics [20, 21, 22], quantum many-body systems [23, 24, 25], chemical
reaction networks [26] and computational biology [27]. However, evaluating the
dynamical partition function and characterizing dynamical phase transitions
have not been achieved in these applications of the VAN. Here, we leverage the
VAN to propose an algorithm for tracking the evolving probability
distribution, leading to an efficient computation of the dynamical partition
function. The latter serves as the moment-generating function of dynamical
observables, which uncovers the phase transition over time (Fig. 1).

We apply the approach to representative models in nonequilibrium statistical
mechanics. We first validate the method in the voter model, by comparing with
the analytical result [6]. To demonstrate that our approach can reveal unknown
dynamical phenomena, we investigate KCMs of spin flips for glassy dynamics,
including the Fredrickson-Andersen (FA) model [28] and variants of the East
model [29], in one dimension (1D), two dimensions (2D), and three dimensions
(3D). Previously, the dynamical active-inactive phase transition in space and
time [30, 31, 32] was investigated mainly at the steady state in the long-time
limit, where KCMs have two phases with extensive or subextensive spin-flipping
activities. However, the phase transition in the full-time regime was seldom
investigated. For the 1D finite-time problem, our results agree with the
recent study by matrix product states [33]. For 2D and 3D cases, where the
phase transition at arbitrary time was not obtained either analytically or
numerically, our method uncovers the dynamical phase diagram versus time the
control parameter of the counting field, and the critical exponent of the
finite-time scaling. We also observe the emergence of characteristic spatial
structures over time, extending the steady-state result [34]. We further
discuss the applications to nonequilibrium systems with other types of
dynamics and topologies.

Figure 1: Uncovering dynamical phase transitions in nonequilibrium statistical
mechanics by machine learning. (a) The major theme of equilibrium and
nonequilibrium statistical mechanics. (b) Uncovering dynamical phase
transitions requires studying the “tilted” dynamics under $\mathbb{T}_{s}$
with various values of the control parameter $s$, where the evolved
distribution is no longer normalized (dashed light blue). As the normalization
factors, the dynamical partition function needs to be learned over time, to
discover the dynamical phase diagram. (c) A proposed algorithm for learning
the dynamical partition function by training the VAN (distributions with hat)
over time. The loss function is given by the KL-divergence between the VAN at
time $t+\delta t$ with learnable parameters $\theta_{t+\delta t}$ and the
evolved distribution learned at $t$ with parameters $\theta_{t}$ fixed. The
algorithm tracks nonequilibrium dynamics , and reveals dynamical phase
transitions, overcoming the three challenges in (b).

##  II Results

###  II.1 Nonequilibrium statistical mechanics

We consider a continuous-time discrete-state Markovian dynamics of size $N$.
Each variable has $d_{s}$ states ($d_{s}=2$ for binary spin systems), giving
total $M=d_{s}^{N}$ configurations. With configuration states
$|\textbf{x}\rangle\equiv(x_{1},x_{2},\dots,x_{N})$ forming an orthonormal
basis, the system at time $t$ is described by the probability vector
$|P_{t}\rangle=\sum_{\textbf{x}}P_{t}(\textbf{x})|\textbf{x}\rangle$. It
evolves under the stochastic master equation [35]:

$\displaystyle\frac{d}{dt}|P_{t}\rangle=\mathbb{W}|P_{t}\rangle,$ (1)

where
$\mathbb{W}=\sum_{\textbf{x},\textbf{x}^{{}^{\prime}}\neq\textbf{x}}w_{\textbf{x},\textbf{x}^{{}^{\prime}}}|\textbf{x}^{{}^{\prime}}\rangle\langle\textbf{x}|-\sum_{\textbf{x}}r_{\textbf{x}}|\textbf{x}\rangle\langle\textbf{x}|$
is the generator, $w_{\textbf{x},\textbf{x}^{{}^{\prime}}}$ is the transition
rate from $|\textbf{x}\rangle$ to $|\textbf{x}^{{}^{\prime}}\rangle$, and
$r_{\textbf{x}}=\sum_{\textbf{x}^{{}^{\prime}}\neq\textbf{x}}w_{\textbf{x},\textbf{x}^{{}^{\prime}}}$
is the escape rate from $|\textbf{x}\rangle$.

Dynamical partition function. Studying dynamical phase transitions emerging
from microscopic interactions relies on evaluating dynamical observables. In
general, a quantity of interest is the time-extensive dynamical observable
$\hat{K}$ incremented along a trajectory
$\omega_{t}=\\{\textbf{x}_{t_{0}}\rightarrow\textbf{x}_{t_{1}}\dots\rightarrow\textbf{x}_{t}\\}$,
with time-step length $\delta t$, $t=J\delta t$ and total $J$ time steps. The
probability of observing the dynamical observable with a value $K$ is obtained
by summing over all possible trajectories:
$P_{t}(K)=\sum_{\omega_{t}}p(\omega_{t})\delta[\hat{K}(\omega_{t})-K]$, where
$p(\omega_{t})$ is the probability of the trajectory $\omega_{t}$.

To extract the statistics of the dynamical observable, it is useful to
estimate its moment-generating function, i.e., the dynamical partition
function
$Z_{t}(s)=\sum_{K}P_{t}(K)e^{-sK}=\sum_{\omega_{t}}p(\omega_{t})e^{-s\hat{K}(\omega_{t})}$,
where a control parameter $s$ is introduced as the conjugate variable to the
dynamical observable. Taking derivatives of $Z_{t}(s)$ to $s$ gives moments of
the dynamical observable’s distribution. The dynamical partition function can
be evaluated by [30]:

$\displaystyle Z_{t}(s)=\langle-|e^{t\mathbb{W}_{s}}|\text{ss}\rangle,$ (2)

where $\langle-|=\sum_{\mathbf{x}}\langle\mathbf{x}|$, and $|\text{ss}\rangle$
is the steady-state probability vector under the generator $\mathbb{W}$.
Different from $\mathbb{W}$, $\mathbb{W}_{s}$ is termed as the “tilted”
generator, and its form depends on the dynamical observable $\hat{K}$ and
control parameter $s$.

A representative dynamical observable. For the spin-flip dynamics, a
characteristic dynamical observable is the “dynamical activity” [30, 36],
measuring the number of spin flips. The observable quantifies how dynamically
active the trajectories are: a more active trajectory has a higher value. The
control parameter $s$ is then termed as the “counting field”. To evaluate the
dynamical observable, the generator is modified correspondingly as the tilted
generator [30]:

$\displaystyle\mathbb{W}_{s}=\sum_{\textbf{x},\textbf{x}^{{}^{\prime}}\neq\textbf{x}}e^{-s}w_{\textbf{x},\textbf{x}^{{}^{\prime}}}|\textbf{x}^{{}^{\prime}}\rangle\langle\textbf{x}|-\sum_{\textbf{x}}r_{\textbf{x}}|\textbf{x}\rangle\langle\textbf{x}|,$
(3)

under which the probability is no longer normalized. Combining Eq. (2) and Eq.
(3), the moments can be calculated, including the average dynamical activity
per unit of time and site:

$\displaystyle k_{t}(s)=-\frac{1}{Nt}\frac{d}{ds}\ln Z_{t}(s).$ (4)

Other dynamical observables can be investigated in a similar way by using its
corresponding tilted generator.

###  II.2 Tracking nonequilibrium statistical mechanics and dynamical phase
transitions

A natural approach of studying nonequilibrium dynamics is tracking the
evolution equation, Eqs. (1) and (2). Unfortunately, the exact representation
of $|P_{t}\rangle$ requires a computational effort that is exponential in the
system size. Hence we need an efficient method to approximately represent
$|P_{t}\rangle$. Here, we consider a neural-network model, the VAN [19], as a
variational ansatz for $|P_{t}\rangle$. The VAN maps configurations to the
probability distribution $P_{t}(\textbf{x})$ in $|P_{t}\rangle$ as the product
of conditional probabilities:

$\displaystyle\hat{P}^{\theta}_{t}(\textbf{x})=\prod^{N}_{i=1}\hat{P}^{\theta}_{t}(x_{i}|x_{1},\dots,x_{i-1}),$
(5)

where the hat symbol denotes the parameterization by the neural network with
learnable parameters $\theta$.

For lattice systems, we start from an initial site and traverse the lattice in
a predetermined order to acquire all the conditional probabilities (Methods).
Each conditional probability $\hat{P}^{\theta}_{t}(x_{i}|x_{1},\dots,x_{i-1})$
is parameterized by a neural network, where the input is the sites visited
earlier $\\{x_{i^{{}^{\prime}}<i}\\}$ and the output is $x_{i}$ associated
with a probability under proper normalization [19]. The parameters of the
neural network can be shared among sites [37] to increase computational
efficiency. Since all the conditional probabilities are stored, one can
efficiently generate samples associated with normalized probabilities, which
can be used to compute quantities such as energy and entropy and to construct
loss functions to update the parameters.

The VAN is capable of expressing strongly correlated distributions [20, 22],
including equilibrium distributions in statistical mechanics [19], steady
state distributions of KCMs in nonequilibrium statistical mechanics [34] and
quantum systems [21]. Here, we find it effective to learn the time-evolving
distributions. The expressivity of the VAN and training time depends on the
architecture, as well as the depth and width of the neural network
(Supplementary Information Sect. VD). Typical neural-network architectures can
be employed, including MADE [38], PixelCNN [39, 40] or RNN [41]. For example,
we have used RNN for KCMs in 1D and 2D, PixelCNN for 2D, and MADE for 3D. In
1D, we find RNN more accurate than MADE. In 2D, RNN has a comparable accuracy
with PixelCNN but takes a longer computational time. The VAN can be further
improved by cooperating with more advanced neural network architectures and
sampling techniques.

The advantages of the VAN over the tensor network model in representing
$|P_{t}\rangle$ mainly come from its generality. The neural network ansatz can
be used in systems with various topologies, in contrast to the matrix product
states designed for one-dimensional systems. On the one hand, the VAN is
suitable for systems with arbitrary topology without modifying the structure
of VAN [19]; on the other hand, we can design the VAN to fit the topology. For
example, in 2D or 3D, one can use convolutional networks [42], as in image
recognition; for sparse networks, graph neural networks are efficient [43].

Based on the VAN representation, we evaluate Eq. (2) by applying the operator
$e^{\delta t\mathbb{W}_{s}}$ sequentially at each of the total $J$ time steps:
$Z_{t}(s)\approx\langle-|(\mathbb{T}_{s})^{J}|\text{ss}\rangle$ (Suzuki-
Trotter decomposition [44]), where the transition operator
$\mathbb{T}_{s}=(\mathbb{I}+\delta t\mathbb{W}_{s})\approx e^{\delta
t\mathbb{W}_{s}}$, $\mathbb{I}$ denotes the identity operator. Without loss of
generality, we consider the one-step evolution from a normalized probability
vector $|\hat{P}^{\theta_{j}}_{j}\rangle$ at time step $j$ ($j=0,\dots,J-1$)
with parameters $\theta_{j}$. Under the tilted generator, the evolved
probability vector $\mathbb{T}_{s}|\hat{P}^{\theta_{j}}_{j}\rangle$ becomes
unnormalized. Still, since the VAN provides a normalized probability vector,
we can use it to represent $|\hat{P}^{\theta_{j+1}}_{j+1}\rangle$ at time step
$j+1$ and approximate the normalized probability vector
$|\hat{Q}^{\theta_{j}}_{j}\rangle=\mathbb{T}_{s}|\hat{P}^{\theta_{j}}_{j}\rangle/Z_{j+1}(s)$,
by minimizing the Kullback-Leibler (KL) divergence,

$\displaystyle D_{KL}[\hat{P}^{\theta_{j+1}}_{j+1}||\hat{Q}^{\theta_{j}}_{j}]$
$\displaystyle=\mathcal{L}_{j+1}+\ln Z_{j+1}(s),$ (6)
$\displaystyle\mathcal{L}_{j+1}$
$\displaystyle=\sum_{\mathbf{x}}\hat{P}^{\theta_{j+1}}_{j+1}(\mathbf{x},s)\left[\ln\hat{P}^{\theta_{j+1}}_{j+1}(\mathbf{x},s)-\ln\mathbb{T}_{s}\hat{P}^{\theta_{j}}_{j}(\mathbf{x},s)\right].$
(7)

Minimizing the KL-divergence is equivalent to minimizing $\mathcal{L}_{j+1}$,
which plays a role analogous to the variational free energy in equilibrium
statistical mechanics. Since the VAN supports unbiased sampling in parallel to
compute $\mathcal{L}_{j+1}$, we estimate the gradients with respect to
parameters $\theta_{j+1}$ by the REINFORCE algorithm [45]:

$\begin{split}\nabla_{\theta_{j+1}}\mathcal{L}_{j+1}&=\sum_{\mathbf{x}}\hat{P}^{\theta_{j+1}}_{j+1}(\mathbf{x},s)\left\\{\nabla_{\theta_{j+1}}\ln\hat{P}^{\theta_{j+1}}_{j+1}(\mathbf{x},s)\right.\left.\cdot\left[\ln\hat{P}^{\theta_{j+1}}_{j+1}(\mathbf{x},s)-\ln\mathbb{T}_{s}\hat{P}^{\theta_{j}}_{j}(\mathbf{x},s)\right]\right\\},\end{split}$
(8)

where the summation is over the samples from the VAN.

As an essential outcome of the algorithm, the dynamical partition function is
computed as a product of the normalization constants at each time step
$Z_{t}(s)\approx\prod_{j=1}^{J}Z_{j}(s)$ (Supplementary Information Sect.
IIA). For each normalization constant, the nonnegativity of the KL-divergence
ensures that Eq. (7) provides a lower bound as:

$\ln Z_{j+1}(s)\geq-\mathcal{L}_{j+1}.$ (9)

The equality holds when the VAN faithfully learns the evolved distribution and
achieves zero KL divergence. Then, to evaluate the dynamical partition
function our algorithm starts from the steady state of the non-tilted
generator and tracks the distribution under the tilted generator (Eq. (2)).
The renormalization procedure over time points enables to extraction of the
dynamical partition function under the tilted generator, beyond the algorithm
of only tracking the evolving distribution under the non-tilted generator
[26].

###  II.3 Algorithm

The pseudocode of tracking the dynamical partition function and dynamical
phase transitions by the VAN is summarized below:

  * •

Input: The system size and dimension, the model type, the boundary condition,
time steps, values of the counting field $s$. Choose an initial distribution,
such as the steady state of the non-tilted dynamics.

  * •

Every time step $j=0,\dots,J-1$:

    1. 1.

Learn the next-step VAN $\hat{P}^{\theta_{j+1}}_{j+1}(\textbf{x})$: for every
epoch,

      1. (a)

Draw samples $\\{\textbf{x}\\}$ from
$\hat{P}^{\theta_{j+1}}_{j+1}(\textbf{x})$;

      2. (b)

Calculate relevant matrix elements of the transition operator $\mathbb{T}_{s}$
to get $\mathbb{T}_{s}\hat{P}^{\theta_{j}}_{j}(\textbf{x})$;

      3. (c)

Train the VAN by minimizing the loss function:
$\mathcal{L}_{j+1}=\mathbb{E}_{\textbf{x}\sim\hat{P}^{\theta_{j+1}}_{j+1}}[\ln\hat{P}^{\theta_{j+1}}_{j+1}(\textbf{x})-\ln\mathbb{T}_{s}\hat{P}^{\theta_{j}}_{j}(\textbf{x})]$.

    2. 2.

Calculate the variational free energy at each time step:
$\mathcal{F}_{j+1}^{\theta_{j+1}}(s)=\mathcal{L}_{j+1}$ after training.

  * •

Estimate the dynamical partition function:
$Z_{t}(s)\approx\prod_{j=1}^{J}Z_{j}(s)$, $\ln
Z_{j}(s)\approx\mathcal{F}_{j}^{\theta_{j}}(s)$; and moments of dynamical
observables.

  * •

Output: The evolved probability distributions, dynamical partition function,
and dynamical observables.

The computational complexity of tracking the dynamical phase transition
depends on the number of training steps $N_{\text{train}}$ and the number of
counting-field values $N_{s}$. Under a fixed system size, the computational
time of studying the steady-state phase transition [34] has the order of
$\mathcal{O}(N_{s}N_{\text{train}})$, and the finite-time study without the
phase transition [25] has $\mathcal{O}(N_{\text{train}}J)$ for $J$ total time
steps. Here, to uncover the phase transition at any time, the computational
complexity becomes $\mathcal{O}(N_{s}N_{\text{train}}J)$, where larger
$N_{\text{train}}$ may be required to reach high accuracy for larger system
sizes (Supplementary Fig. 4). Despite the high computational complexity, the
result under the chosen lattice sizes fulfills to uncover the finite-time
scaling of the phase transition for 2D and 3D KCMs, as demonstrated below.

Figure 2:  A schematic on kinetically constrained models and dynamical phase
diagram. (a) The rule of spin flips for the 2D SE and FA models, as
demonstrative examples. The filled (unfilled) circle is up (down) spin. The
flipping rate of each spin (red) depends on $f$, the number of up spins at
certain nearest neighbors (green). The SE model counts the left and above
neighbors, and the FA model counts all nearest neighbors. The flip-up
probability is $c$. (b) An illustrative dynamical phase diagram, with active
and inactive phases of spin flips depending on the counting field $s$ and
time. The active (inactive) phase has more (less) spin flips over time, giving
large (small) dynamical activity $k_{t}(s)$ in Eq. (4). (c) The two
architectures of the VAN used in 2D (Methods). The left panels have a
schematic of the gated recurrent unit cell (GRU) [46]. For 2D lattice systems,
a zigzag path for transmitting variables is employed, allowing the hidden
state to propagate both vertically and horizontally. The right panels are
PixelCNNs. The vanilla PixelCNN [39] encounters an issue of the “blind spot”,
i.e., the black lattice point is not conditioned on the lattice points within
the blind spot area under a $3\times 3$ causal convolution layer. This issue
is solved in the gated PixelCNN [40] by splitting the model into two parts:
vertical and horizontal stacks.

###  II.4 Applications

The present approach is generally applicable to systems in nonequilibrium
statistical mechanics. Since the computational cost is proportional to the
number of allowable state transitions at each time step in Eq. 2, the cost is
greatly reduced when compared with tracking the full distribution
exponentially proportional to the system size. This is computationally
feasible when the number of allowable transitions scales polynomially on the
system size, including the voter model [6] which validates our method of
tracking the distribution (Supplementary Fig. 1). Another representative class
of systems is KCMs modeling glasses [8]. The KCMs exhibit rich phenomena of
phase transitions (Fig. 2), where the 2D and 3D cases were unexplored either
analytically or numerically, and are our focus.

Figure 3:  Characterization of the dynamical active-inactive phase transition
of kinetically constrained models. The top, middle, and bottom panels are the
2D SE model, 2D FA model and 3D SEF model. (a) The dynamical activity
$k_{t}(s)$ denoted by color reveals the phase transition versus time $t$ and
the counting field $s$, with $L=5$ for 2D and $L=3$ for 3D. (b) The critical
point $s_{c}(N,t)$ over time, giving critical exponents $\alpha$ from
$s_{c}(N)\sim N^{-\alpha}$ at the steady state (the horizontal dotted lines
and the inset) and $\beta$ from the finite-time scaling $t^{-\beta}$ (the
black dashed line). (c) The scaled phase transition lines, with
$s_{c}(N,t)/s_{c}(N)$ and time scaled as $N^{\alpha}t^{-\beta}$, are collapsed
together, indicating the proper scaling relation. Parameters: $c=0.5$,
$L=4,5,6$ for 2D and $L=2,3,4$ for 3D.

Active-inactive phase transition over time for KCMs. The previous studies of
KCMs focused on equilibrium properties [47] or the active-inactive phase
transition in the long-time limit [30, 31, 32, 34]. For the phase transition
at arbitrary time, besides the recent result in 1D [33], the cases with
lattice dimension greater than $1$ have not been revealed. Here, the VAN
provides the phase transition over time of KCMs on 1D, 2D, and 3D lattices in
a unified way.

We consider two paradigmatic KCMs, namely, FA [28] and variants of East [29]
models, on a lattice of size $N=L$ in 1D, $N=L^{2}$ in 2D, and $N=L^{3}$ in
3D, with binary spins $x_{i}=0,1$ for $i=1,\dots,N$, $d_{s}=2$ and $2^{N}$
configurations in total. The Markovian generator is

$\displaystyle\mathbb{W}$
$\displaystyle=\sum_{i=1}^{N}f_{i}[c\sigma_{i}^{+}+(1-c)\sigma_{i}^{-}-c(1-x_{i})-(1-c)x_{i}],$
(10)

where $\sigma_{i}^{\pm}$ are the Pauli raising and lowering operators flipping
site $i$ up and down, and $c\in(0,0.5]$ controls the rate of flipping up. The
up spin number $x_{i}$ acts as the operator
$x_{i}=\sigma_{i}^{+}\sigma_{i}^{-}$, and the terms $1-x_{i}$, $x_{i}$
separately represent the escape transitions out from the down, up spins at
site $i$. The $f_{i}$ equals the number of up spins at certain directions of
the nearest neighbors (Fig. 2a): the FA model counts all directions, the 1D
East counts the left, the 2D South-East (SE) counts the left and above, and
the 3D South-East-Front (SEF) counts the left, above and back. We consider
open boundary conditions for the convenience of comparing results with the
literature. The boundary sites are up for 1D [33] and down for 2D [34] and 3D.
For the 2D FA model, the configuration with all down spins is excluded to
avoid this disconnected configuration. To access the largest ergodic element
in the configuration space, the first spin is fixed up for 2D and 3D East
models. Note that all the figures of configurations do not include boundary
sites.

The generator acts on each configuration and at each time step contributes to
the case with only one spin flip between two configurations. Each
configuration has $N$ connected configurations that transit into or out from.
The probability distribution is updated by using only the connected
configurations of batch samples in Eq. (8). This procedure reduces the
complexity from multiplying the transition matrix with the probability vector,
both of which are exponentially large, to counting an order of
$\mathcal{O}(N)$ transitions linear to the system size.

The VAN shows a high accuracy in revealing the phase transition. The obtained
dynamical partition function coincides with the numerically exact values
available for the small system sizes (Supplementary Fig. 5). Our result at
long time matches with the steady-state estimation from the variational Monte-
Carlo method [34] (Supplementary Fig. 6). For the finite-time regime which was
previously explored only for 1D, our result (Supplementary Figs. 2,3) agrees
with tensor networks [33]. For the unexplored 2D and 3D KCMs at finite time,
we obtain the dynamical activity $k_{t}(s)$ as a function of time $t$ and the
counting field $s$ with $s>0$, showing two phases with extensive or
subextensive activities (Fig. 3a). The critical point $s_{c}(N,t)$ as a
function of system size and time can be identified numerically by the peak
point of the dynamical susceptibility $\chi_{t}(s)=dk_{t}(s)/ds$. We conduct a
scaling analysis of the critical point $s_{c}(N,t)$. In the long-time regime,
the scaling of system size $s_{c}(N)\sim N^{-\alpha}$ gives the exponent
$\alpha\gtrsim 1$ for 2D and 3D (Fig. 3b). By inspecting the short-time
regime, the critical point approximately scales as $s_{c}(t)\sim t^{-1}$ for
the two models. These two regimes motivate to approximate $s_{c}(N,t)\approx
s_{c}(N)+s_{c}(t)$. We then estimate the scaling of time by fitting
$\ln[s_{c}(N,t)-s_{c}(N)]$ versus $-\ln(t)$, giving $s_{c}(t)\sim t^{-\beta}$
with $\beta\approx 1$ for 2D and 3D, which confirms the speculated value. The
finite-time scaling in 2D and 3D has a similar critical exponent as in 1D
[33], implying that the transition point evolves with time similarly in all
the dimensions for KCMs. By dividing $s_{c}(N,t)$ with $s_{c}(N)$ and time as
$N^{\alpha}t^{-\beta}$, $s_{c}(N,t)$ curves are collapsed together (Fig. 3c).

Figure 4: The emergence of the active phases. The 2D SE (a, b) and 2D FA (c,
d) models with $c=0.05$, $L=6$ and $s<0$. (a, c) The average density in the
active phase at various time points for $c=0.05$. The grey arrows point to the
chosen $\nu$ (introduced for $s<0$) for panels (b, d). (b, d) Time evolution
of the spatial profile of configurations at various time points (rows) for
different $\nu$ values (columns). For the comparison between the two models,
the left-up corner spin of both models is fixed up. The color denotes the
average number of up spins, with the same color bar in (d).

The emergence of active phases for 2D KCMs. We further analyze the emergence
of characteristic spatial structures over time, beyond the steady state [34].
For 2D SE and FA models, with a small $c=0.05$, we discover that the average
density of up spins $\langle n(t)\rangle_{\nu}$ (versus $\nu\equiv 1-e^{s}$)
of the 2D SE model gradually shows piecewise density plateaus (Fig. 4a), which
are absent in the 2D FA model (Fig. 4c). For the average spatial profile of
configurations, the 2D SE model forms structures with up-spin diagonal bands
separated by down-spin bands, and the number of bands varies over $\nu$ (Fig.
4b). The 2D FA model does not have such characteristic bands (Fig. 4d), even
when its first spin is also fixed up.

##  III Discussion

We have presented a general framework to track dynamics of nonequilibrium
statistical systems based on neural networks, and developed an efficient
algorithm to estimate the dynamical partition function. This extends
applications of the VAN from equilibrium [19] to nonequilibrium, from the
steady state [34] to arbitrary time, and from the absence of phase transitions
[25] to phase transitions. , the approach enables to reveal of the unexplored
active-inactive phase transition at finite time and scaling relations in KCMs
on 2D and 3D lattices, as well as the emergence of characteristic spatial
structures.

The error, such as quantified by the relative error of the dynamical partition
function in Eq. (16), does not keep increasing (Supplementary Figs. 4, 5),
because the time evolution of the distribution does not have a dramatic change
for all time points. The accumulation of error mainly occurs when the dynamics
change dramatically over time, e.g., when the active-inactive phase transition
occurs. Identifying these time points by trial and error helps find the most
efficient way of increasing epochs for accuracy at certain time points. An
alternative way of resolving the issue is to project the evolution into two
parts: one follows the largest eigenvalue of the tilted generator, and the
other follows from a modified tilted generator (Supplementary Sect. IIA). When
the modified generator reaches a steady value, one can stop the simulation and
extrapolate it by using the largest eigenvalue.

Capturing the active-inactive transition requires an efficient sampling of
rare inactive configurations with few up spins, which is accessible by the
importance sampling (Supplementary Sect. IIIA). This important sampling is on
the configurations from a distribution at each time point, different from the
sampling on trajectories, which may be harder to sample as the trajectory
space grows exponentially with time points. Besides, learning the multiple
probability peaks can be affected by the mode-collapse: For rugged
distributions, not all modes of the target distribution may be directly
captured by the VAN [48]. To alleviate the mode-collapse, besides the
importance sampling used here, temperature annealing [19] and variational
annealing [49] can be employed.

When the system size increases, the computational time reaches the order of
$\mathcal{{O}}(10^{2})$ hours for one value of the counting field
(Supplementary Table 1). Larger system sizes require longer computational time
to reach high accuracy (Supplementary Fig. 4). Evaluating the phase transition
needs to scan various values of the counting field. Thus, although the method
is generally applicable, computing the scaling relation for larger system
sizes meets the practical challenges of available computational resources.
However, this can be alleviated by optimizing the efficiency of tracking the
distribution, with the help of the time-dependent variational method [50] and
by paralleling multiple GPUs.

The present approach is applicable to other types of Markovian dynamics,
including stochastic reaction networks [26, 51], where phase transitions can
be analyzed after adding the control parameter. Based on generality of the
VAN, it is adaptable to other topologies, such as the voter model on graphs
[52] and epidemic spreading on networks [53], where the architecture of the
VAN can be the graph neural network [43]. It may also be generalized to the
KCMs of various dimensions in the quantum regime [54]. Another direction is to
leverage the VAN for sampling rare trajectories, with the help of the Doob
operator [55], active learning [56] and reinforcement learning [57].

##  IV Methods

###  IV.1 Variational autoregressive networks for spin-lattice systems

We use the variational autoregressive network [19] to parameterize the
probability distribution. The VAN factorizes the joint probability into a
product of conditional probabilities as Eq. (5), where $x_{i}$ denotes the
$d_{s}$-state variable of site $i$ ($d_{s}=2$ for binary spin systems). The
symbol $\theta$ represents the learnable parameters. The parameterized
distribution is automatically normalized, which is also called autoregressive
modeling in the machine learning community. Since each conditional probability
only depends on previous sites, it supports efficient ancestral sampling in
parallel.

Below, we briefly describe the architecture of the RNN and gated PixelCNN for
our problem. The setting of MADE was the same as [19].

Recurrent neural networks. For the recurrent neural network (RNN), we use a
gated recurrent unit [46] as the recurrent cell, which is capable of learning
the distribution with long-range correlations. It is more efficient than the
long-short time memory (LTSM) model and avoids the vanishing gradient problem
for vanilla recurrent neural networks.

For the 1D RNN, the conditional probability is iteratively obtained over the
one-dimensional sites. A recurrent cell processes the information from the
previous hidden state $h_{i-1}$ and the input data $x_{i-1}$ in the current
cell, generates a new hidden state $h_{i}$, gives the conditional probability
$\hat{P}^{\theta}(x_{i}|x_{1},\dots,x_{i-1})$ based on $h_{i}$, and passes on
the information of $h_{i}$ to the next cell. The dimension of the hidden
states is denoted by $d_{h}$. The GRU has a candidate hidden state
$\hat{h}_{i}$, an update gate $z_{i}$ interpolating between the previous and
candidate hidden states, and a reset gate $r_{i}$ setting the extent of
forgetting for the previous hidden state. It updates by the following gates:

$\displaystyle z_{i}$
$\displaystyle=\sigma(W_{zx}x_{i-1}+W_{zh}h_{i-1}+b_{z}),$ (11) $\displaystyle
r_{i}$ $\displaystyle=\sigma(W_{rx}x_{i-1}+W_{rh}h_{i-1}+b_{r}),$ (12)
$\displaystyle\hat{h}_{i}$
$\displaystyle=\tanh(W_{hx}x_{i-1}+W_{hh}(r_{i}\odot h_{i-1})+b_{h}),$ (13)
$\displaystyle h_{i}$ $\displaystyle=(1-z_{i})\odot
h_{i-1}+z_{i}\odot\hat{h}_{i},$ (14)

where $W$s are the weight matrices, $b$s are bias vectors, and $\sigma$ is the
sigmoid activation function, $\odot$ denotes the Hadamard product.

The conditional probability is obtained from the hidden states. The output is
acted on by a linear transform and a softmax operator
$\hat{P}^{\theta}(x_{i}|x_{1},\dots,x_{i-1})=\text{Softmax}(Wh_{i}+b)$, which
ensures the normalized condition for the output probability vector. Given an
initial hidden state $h_{0}$ and variable $x_{0}$ (chosen as a zero vector
here), the full probability is obtained by Eq. (5) with the iteratively
generated conditional probabilities. Sampling from the probability
distribution is conducted similarly: given an initial hidden state and
variable, the variable $x_{1}$ is sampled from the estimated conditional
probability, and the procedure is repeated to the last site.

For the 2D RNN, the implementation is more involved. A zigzag path [41, 34] is
used to transmit the lattice variables, both vertically and horizontally. For
the vertical and horizontal variables (hidden state $h$ and variable $x$
separately), we first concentrate the two into one and perform a linear
transform (without the bias term) to an intermediate variable with the
original dimension. They are then passed to the next GRU cell to continue the
iteration.

The gated PixelCNN. For 2D lattice systems, we find that the vanilla PixelCNN
model [39] suffers from the blind spot problem. In the worst case, the blind
spot in the receptive field only covers half of the sites above and to the
left of the current site. To circumvent the blind spot problem, we use the
gated PixelCNN [40] that combines two convolutional network stacks: the
vertical stack and the horizontal stack. The vertical stack conditions on all
the sites left and the horizontal stack conditions on all the sites above. The
activation function is also replaced by a gated activation unit:

$\mathbf{h}^{l+1}=\tanh(W_{1}^{l}\ast\mathbf{h}^{l})\odot\sigma(W_{2}^{l}\ast\mathbf{h}^{l})$
(15)

where $l$ is the number of layers, $\mathbf{h}^{l}$ is the feature map at the
$l$-th layer, $\odot$ is the Hadamard product and $\ast$ denotes the stacked
convolution operation.

###  IV.2 Details of training neural networks

When tracking the dynamics in Eq. (1), shorter Trotter time-step length
generally gives higher accuracy, with the cost of longer simulation time. For
KCMs, the range of time-step length $\delta t\in[0.01,0.1]$ is found suitable,
as also reported in the 1D case [33]. Considering both the accuracy and
efficiency, $\delta t$ here is often chosen as $\delta t=0.1$ for 1D and
$\delta t=0.05$ for 2D and 3D. The loss usually takes a number of
$>\mathcal{O}(10^{3})$ epochs to converge for the first time step of evolving
the system but requires only an order of $\mathcal{O}(10^{2})$ epochs for the
following time steps, because the probability distribution only has a small
change after each time step. The smaller number of epochs after the first time
step saves the training time for tracking the evolution.

Learning rates affect the accuracy of training. Among the tested learning
rates $10^{-5}$, $10^{-4}$, $10^{-3}$, and $10^{-2}$, we found that $10^{-3}$
typically leads to relatively lower loss values and better training. It is
possible to encounter general optimization issues such as trapping into local
minima, which may be alleviated by designing schedulers for the learning rate.

We used the Adam optimizer [58] to perform the stochastic gradient descent.
The batch size was usually set as $1000$ for each epoch. To better estimate
the variational free energy, we used the batch from the last $20\%$ of the
training epochs, where the loss converges. This average gives a more accurate
estimate by using approximately $1000\times 100\times 20\%$ (batch size,
epochs at each time point, last percent of the training epochs) batch samples.

###  IV.3 The relative error

We estimate the error of the VAN for each system separately. To quantify the
accuracy of the VAN, we calculated the relative error of the dynamical
partition functions between the VAN and the numerically exact result. The
numerically exact result of the dynamical partition function was obtained by
summing up probabilities of all possible states, which is feasible only for
systems with small sizes, e.g., approximately $L\leq 10$ for 1D, $L\leq 4$ for
2D, and $L\leq 2$ for 3D.

The relative error $e_{r}$ is defined as:

$\displaystyle e_{r}=\left|\frac{\ln Z_{t}(s)-\ln Z_{t}(s)_{\text{exact}}}{\ln
Z_{t}(s)_{\text{exact}}}\right|,$ (16)

where $\ln Z_{t}(s)_{\text{exact}}$ is the numerically exact result by summing
up probabilities of all possible states, which is feasible for systems with
small sizes. The error is shown in the inset of Supplementary Fig. 2a for 1D
and Supplementary Fig. 5 for 2D and 3D. All the figures show that the relative
error is usually smaller than $\mathcal{O}(10^{-3})$ compared with the
numerically exact result, validating the accuracy of the VAN. Based on the
accuracy and efficiency, we chose the more appropriate VAN for each dimension:
the RNN in 1D, the gated PixelCNN in 2D, and the MADE in 3D.

For larger system sizes, the numerically exact result is not accessible which
then demands the use of the present algorithm. Thus, we can only evaluate the
error from the loss function Eq. (7) based on the KL-divergence: the lower
value of the loss function implies more accurate training, with the lower
bound given by the dynamical partition function. We remark that the loss
function will not reach zero, even when the VAN is accurately learned because
the probability is not normalized under the tilted generator. Under this case,
it is not feasible to estimate the KL divergence of the two probability
distributions at consecutive time points as a quantification of accuracy.
Besides, the value of the loss function decreases with a larger number of
epochs and increases with the system size (Supplementary Fig. 4), indicating
the requirement for more epochs and longer computational time when the system
size increases.

Data availability: The authors declare that the data supporting this study are
available within the paper.

Code availability: A PyTorch implementation of the present algorithm can be
found in Supplementary Data 1 and at the GitHub repository
(https://github.com/Machine-learning-and-complex-systems/DPT).

## Acknowledgments

We thank Juan P. Garrahan, Luke Causer, and Corneel Casert for their helpful
communication. We also acknowledge Online Club Nanothermodynamica for
discussions. This work is supported by Project 11747601 (P.Z.), 12325501
(P.Z.), 12247104 (P.Z.), 12322501 (Y.T.), and 12047503 of the National Natural
Science Foundation of China. P.Z. is partially supported by the Innovation
Program for Quantum Science and Technology project 2021ZD0301900 and ZDRW-
XX-2022-3-02 of the Chinese Academy of Sciences. The high-performance
computing is supported by Dawning Information Industry Corporation Ltd.

## Author contributions

P.Z., Y.T. had the original idea for this work. Y.T., J.L. and J.Z. performed
the study, and all authors contributed to the preparation of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information The online version contains supplementary material
available at [URL will be inserted by publisher].

Correspondence and requests for materials should be addressed to Pan Zhang and
Ying Tang.

Reprints and permission information is available online at [URL will be
inserted by publisher].

## References

  * sch [2013] _Nonequilibrium Statistical Physics of Small Systems: Fluctuation Relations and Beyond_ , edited by R. Klages, W. Just, and C. Jarzynski (Wiley-VCH, New York, 2013).
  * Esposito  _et al._ [2009] M. Esposito, U. Harbola, and S. Mukamel, Nonequilibrium fluctuations, fluctuation theorems, and counting statistics in quantum systems, Rev. Mod. Phys. 81, 1665 (2009).
  * Seifert [2012] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Rep. Prog. Phys. 75, 126001 (2012).
  * Chou  _et al._ [2011] T. Chou, K. Mallick, and R. K. Zia, Non-equilibrium statistical mechanics: from a paradigmatic model to biological transport, Rep. Prog. Phys. 74, 116601 (2011).
  * Tang and Hoffmann [2022] Y. Tang and A. Hoffmann, Quantifying information of intracellular signaling: progress with machine learning, Rep. Prog. Phys. 85, 086602 (2022).
  * Krapivsky  _et al._ [2010] P. L. Krapivsky, S. Redner, and E. Ben-Naim,  _A kinetic view of statistical physics_ (Cambridge University Press, 2010).
  * Ritort and Sollich [2003] F. Ritort and P. Sollich, Glassy dynamics of kinetically constrained models, Adv. Phys. 52, 219 (2003).
  * Garrahan and Chandler [2002] J. P. Garrahan and D. Chandler, Geometrical explanation and scaling of dynamical heterogeneities in glass forming systems, Phys. Rev. Lett. 89, 035704 (2002).
  * Dornic  _et al._ [2001] I. Dornic, H. Chaté, J. Chave, and H. Hinrichsen, Critical coarsening without surface tension: The universality class of the voter model, Phys. Rev. Lett. 87, 045701 (2001).
  * Bolhuis  _et al._ [2002] P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, Transition path sampling: Throwing ropes over rough mountain passes, in the dark, Annu. Rev. Phys. Chem. 53, 291 (2002).
  * Touchette [2009] H. Touchette, The large deviation approach to statistical mechanics, Phys. Rep. 478, 1 (2009).
  * Pressé  _et al._ [2013] S. Pressé, K. Ghosh, J. Lee, and K. A. Dill, Principles of maximum entropy and maximum caliber in statistical physics, Rev. Mod. Phys. 85, 1115 (2013).
  * Ray  _et al._ [2018] U. Ray, G. K.-L. Chan, and D. T. Limmer, Exact fluctuations of nonequilibrium steady states from approximate auxiliary dynamics, Phys. Rev. Lett. 120, 210602 (2018).
  * Helms and Chan [2020] P. Helms and G. K.-L. Chan, Dynamical phase transitions in a 2d classical nonequilibrium model via 2d tensor networks, Phys. Rev. Lett. 125, 140601 (2020).
  * Ge and Qian [2012] H. Ge and H. Qian, Analytical mechanics in stochastic dynamics: Most probable path, large-deviation rate function and hamilton–jacobi equation, Int. J. Mod. Phys. B 26, 1230012 (2012).
  * Tang  _et al._ [2015] Y. Tang, R. Yuan, J. Chen, and P. Ao, Work relations connecting nonequilibrium steady states without detailed balance, Phys. Rev. E 91, 042108 (2015).
  * Tang [2021] Y. Tang, Free energy amplification by magnetic flux for driven quantum systems, Commun. Phys. 4, 1 (2021).
  * Mallick  _et al._ [2022] K. Mallick, H. Moriya, and T. Sasamoto, Exact solution of the macroscopic fluctuation theory for the symmetric exclusion process, Phys. Rev. Lett. 129, 040601 (2022).
  * Wu  _et al._ [2019] D. Wu, L. Wang, and P. Zhang, Solving statistical mechanics using variational autoregressive networks, Phys. Rev. Lett. 122, 080602 (2019).
  * Mehta  _et al._ [2019] P. Mehta, M. Bukov, C.-H. Wang, A. G. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, A high-bias, low-variance introduction to machine learning for physicists, Phys. Rep. (2019).
  * Carleo  _et al._ [2019] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Machine learning and the physical sciences, Rev. Mod. Phys. 91, 045002 (2019).
  * McNaughton  _et al._ [2020] B. McNaughton, M. V. Milošević, A. Perali, and S. Pilati, Boosting monte carlo simulations of spin glasses using autoregressive neural networks, Phys. Rev. E 101, 053312 (2020).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Westerhout  _et al._ [2020] T. Westerhout, N. Astrakhantsev, K. S. Tikhonov, M. I. Katsnelson, and A. A. Bagrov, Generalization properties of neural network approximations to frustrated magnet ground states, Nat. Commun. 11, 1593 (2020).
  * Luo  _et al._ [2022] D. Luo, Z. Chen, J. Carrasquilla, and B. K. Clark, Autoregressive neural network for simulating open quantum systems via a probabilistic formulation, Phys. Rev. Lett. 128, 090501 (2022).
  * Tang  _et al._ [2023] Y. Tang, J. Weng, and P. Zhang, Neural-network solutions to stochastic reaction networks, Nat. Mach. Intell. 5, 376 (2023).
  * Shin  _et al._ [2021] J.-E. Shin, A. J. Riesselman, A. W. Kollasch, C. McMahon, E. Simon, C. Sander, A. Manglik, A. C. Kruse, and D. S. Marks, Protein design and variant prediction using autoregressive generative models, Nat. Commun. 12, 2403 (2021).
  * Fredrickson and Andersen [1984] G. H. Fredrickson and H. C. Andersen, Kinetic ising model of the glass transition, Phys. Rev. Lett. 53, 1244 (1984).
  * Jäckle and Eisinger [1991] J. Jäckle and S. Eisinger, A hierarchically constrained kinetic ising model, Z. fur Phys. B 84, 115 (1991).
  * Garrahan  _et al._ [2007] J. P. Garrahan, R. L. Jack, V. Lecomte, E. Pitard, K. van Duijvendijk, and F. van Wijland, Dynamical first-order phase transition in kinetically constrained models of glasses, Phys. Rev. Lett. 98, 195702 (2007).
  * Nemoto  _et al._ [2017] T. Nemoto, R. L. Jack, and V. Lecomte, Finite-size scaling of a first-order dynamical phase transition: Adaptive population dynamics and an effective model, Phys. Rev. Lett. 118, 115702 (2017).
  * Bañuls and Garrahan [2019] M. C. Bañuls and J. P. Garrahan, Using matrix product states to study the dynamical large deviations of kinetically constrained models, Phys. Rev. Lett. 123, 200601 (2019).
  * Causer  _et al._ [2022] L. Causer, M. C. Bañuls, and J. P. Garrahan, Finite time large deviations via matrix product states, Phys. Rev. Lett. 128, 090605 (2022).
  * Casert  _et al._ [2021] C. Casert, T. Vieijra, S. Whitelam, and I. Tamblyn, Dynamical large deviations of two-dimensional kinetically constrained models using a neural-network state ansatz, Phys. Rev. Lett. 127, 120602 (2021).
  * Gardiner [2004] C. W. Gardiner,  _Handbook of Stochastic Methods_ , 3rd ed. (Springer-Verlag, Berlin, 2004).
  * Maes [2020] C. Maes, Frenesy: Time-symmetric dynamical activity in nonequilibria, Phys. Rep. 850, 1 (2020).
  * Uria  _et al._ [2016] B. Uria, M.-A. Côté, K. Gregor, I. Murray, and H. Larochelle, Neural autoregressive distribution estimation, J. Mach. Learn. Res. 17, 7184 (2016).
  * Germain  _et al._ [2015] M. Germain, K. Gregor, I. Murray, and H. Larochelle, Made: Masked autoencoder for distribution estimation, in  _International conference on machine learning_ (PMLR, 2015) pp. 881–889.
  * Van Oord  _et al._ [2016] A. Van Oord, N. Kalchbrenner, and K. Kavukcuoglu, Pixel recurrent neural networks, in  _International conference on machine learning_ (PMLR, 2016) pp. 1747–1756.
  * Van den Oord  _et al._ [2016] A. Van den Oord, N. Kalchbrenner, L. Espeholt, O. Vinyals, A. Graves, _et al._ , Conditional image generation with pixelcnn decoders, Advances in neural information processing systems 29 (2016).
  * Hibat-Allah  _et al._ [2020] M. Hibat-Allah, M. Ganahl, L. E. Hayward, R. G. Melko, and J. Carrasquilla, Recurrent neural network wave functions, Phys. Rev. Research 2, 023358 (2020).
  * Gu  _et al._ [2018] J. Gu, Z. Wang, J. Kuen, L. Ma, A. Shahroudy, B. Shuai, T. Liu, X. Wang, G. Wang, J. Cai, _et al._ , Recent advances in convolutional neural networks, Pattern Recognit. 77, 354 (2018).
  * Scarselli  _et al._ [2008] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini, The graph neural network model, IEEE Trans. Neural Netw. 20, 61 (2008).
  * Suzuki [1976] M. Suzuki, Generalized trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems, Commun. Math. Phys. 51, 183 (1976).
  * Williams [1992] R. J. Williams, Simple statistical gradient-following algorithms for connectionist reinforcement learning, Machine Learning 8, 229 (1992).
  * Cho  _et al._ [2014] K. Cho, B. Van Merriënboer, D. Bahdanau, and Y. Bengio, On the properties of neural machine translation: Encoder-decoder approaches, arXiv:1409.1259 (2014).
  * Berthier and Garrahan [2005] L. Berthier and J. P. Garrahan, Numerical study of a fragile three-dimensional kinetically constrained model, J. Phys. Chem. B 109, 3578 (2005).
  * Ciarella  _et al._ [2023] S. Ciarella, J. Trinquier, M. Weigt, and F. Zamponi, Machine-learning-assisted monte carlo fails at sampling computationally hard problems, Mach. Learn. Sci. Technol. (2023).
  * Hibat-Allah  _et al._ [2021] M. Hibat-Allah, E. M. Inack, R. Wiersema, R. G. Melko, and J. Carrasquilla, Variational neural annealing, Nat. Mach. Intell. 3, 952 (2021).
  * Reh  _et al._ [2021] M. Reh, M. Schmitt, and M. Gärttner, Time-dependent variational principle for open quantum systems with artificial neural networks, Phys. Rev. Lett. 127, 230501 (2021).
  * Wen  _et al._ [2023] M. Wen, E. W. C. Spotte-Smith, S. M. Blau, M. J. McDermott, A. S. Krishnapriyan, and K. A. Persson, Chemical reaction networks and opportunities for machine learning, Nat. Comput. Sci. 3, 12 (2023).
  * Nyczka  _et al._ [2012] P. Nyczka, K. Sznajd-Weron, and J. Cisło, Phase transitions in the $q$-voter model with two types of stochastic driving, Phys. Rev. E 86, 011105 (2012).
  * Biazzo  _et al._ [2022] I. Biazzo, A. Braunstein, L. Dall’Asta, and F. Mazza, A bayesian generative neural network framework for epidemic inference problems, Sci. Rep. 12, 1 (2022).
  * Pancotti  _et al._ [2020] N. Pancotti, G. Giudice, J. I. Cirac, J. P. Garrahan, and M. C. Bañuls, Quantum east model: Localization, nonthermal eigenstates, and slow dynamics, Phys. Rev. X 10, 021051 (2020).
  * Causer  _et al._ [2021] L. Causer, M. C. Bañuls, and J. P. Garrahan, Optimal sampling of dynamical large deviations via matrix product states, Phys. Rev. E 103, 062144 (2021).
  * Pickering  _et al._ [2022] E. Pickering, S. Guth, G. E. Karniadakis, and T. P. Sapsis, Discovering and forecasting extreme events via active learning in neural operators, Nat. Comput. Sci. 2, 823 (2022).
  * Rose  _et al._ [2021] D. C. Rose, J. F. Mair, and J. P. Garrahan, A reinforcement learning approach to rare trajectory sampling, New J. Phys. 23, 013013 (2021).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv:1412.6980 (2014).


This isomorphism, together with the classification of elements in
$\widehat{\Gamma^{\prime}}$, induces a bijection at the level of irreducible
unitary representations which maps
$\pi^{\prime}_{\lambda},\pi^{\prime}_{\alpha,\beta}\in\widehat{\Gamma^{\prime}}$
to the representations $\pi_{\lambda},\pi_{\alpha,\beta}\in\widehat{\Gamma}$
given by

$\pi_{\lambda}(x,y,z):=\pi^{\prime}_{\lambda}\left(x,y,z-\frac{1}{2}xy\right)\qquad\text{and}\qquad\pi_{\alpha,\beta}(x,y,z):=\pi^{\prime}_{\alpha,\beta}\left(x,y,z-\frac{1}{2}xy\right),$

respectively. More precisely, the irreducible representations
$\pi_{\lambda}(x,y,z)$ and $\pi_{\alpha,\beta}(x,y,z)$ act over a function
$\varphi\in L^{2}(\mathbb{R})$ by:

$\pi_{\lambda}(x,y,z)\varphi(u)=e^{2\pi i\lambda z}e^{2\pi i\lambda
yu}\varphi(u+x)\qquad\text{and}\qquad\pi_{\alpha,\beta}(x,y,z)=e^{2\pi
i(\alpha x+\beta y)}.$ (5.26)

Following (4.3), for the lattice $H=(\mathbb{Z}^{3},\cdot)\subset\Gamma$ we
define its (averaged) Fourier Transform $\mathcal{F}(\delta_{H})$ by the
following limit in the weak-sense:

$\mathcal{F}(\delta_{H})(\lambda):=\lim_{R\to\infty}\frac{\mathcal{F}(\delta_{H\cap
B_{R}})(\pi_{\lambda})}{R^{3}}=\lim_{R\to\infty}\frac{1}{R^{3}}\sum_{(x,y,z)\in
H\cap B_{R}}\pi_{\lambda}(x,y,z)$ (5.27)

and

$\mathcal{F}(\delta_{H})(\alpha,\beta):=\lim_{R\to\infty}\frac{\mathcal{F}(\delta_{H\cap
B_{R}})(\pi_{\alpha,\beta})}{R^{3}}=\lim_{R\to\infty}\frac{1}{R^{3}}\sum_{(x,y,z)\in
H\cap B_{R}}\pi_{\alpha,\beta}(x,y,z).$ (5.28)

More precisely, the limit (5.27) is defined in duality with functions with
compact support, namely, for every $\varphi,\psi\in C_{c}(\mathbb{R})$ there
holds

$\langle\mathcal{F}(\delta_{H})(\lambda)\varphi,\psi\rangle_{L^{2}(\mathbb{R})}=\lim_{R\to\infty}\frac{1}{R^{3}}\sum_{(x,y,z)\in
H\cap B_{R}}\int_{\mathbb{R}}e^{2\pi
i\lambda(yu+z)}\varphi(x+u)\overline{\psi(u)}du;$

whereas for $\alpha,\beta\in\mathbb{R}$ we have that

$\mathcal{F}(\delta_{H})(\alpha,\beta)=\lim_{R\to\infty}\frac{1}{R^{3}}\sum_{(x,y,z)\in
H\cap B_{R}}e^{2\pi i(\alpha x+\beta y)}.$

For the lattice $H=(\mathbb{Z}^{3},\cdot)$, consider the i.i.d random vectors
$(\xi_{p})_{p\in H}$ defined in a probability space
$(\Omega,\mathcal{T},\mathbb{P})$ with common distribution $\xi$, and the
randomly perturbed set

$H_{\xi}=\\{p\cdot\xi_{p}:\ p\in H\\}\subset\Gamma.$

With these notations, our main recovery result in the Heisenberg group reads
as follows.

######  Theorem 5.7.

Under the previous notations, assume that there is a positive number
$\varepsilon$ such that $\mathbb{E}(|\xi|^{3+\varepsilon})<\infty$. Then
almost surely, for every non-zero $\lambda\in\mathbb{R}$ and
$\alpha,\beta\in\mathbb{R}$, there holds:

$\mathcal{F}(\delta_{H_{\xi}})(\lambda)=\mathcal{F}(\delta_{H})(\lambda)\mathbb{E}[\pi_{\lambda}(\xi)]\qquad\text{and}\qquad\mathcal{F}(\delta_{H_{\xi}})(\alpha,\beta)=\mathcal{F}(\delta_{H})(\alpha,\beta)\mathbb{E}[\pi_{\alpha,\beta}(\xi)].$

Like in the Euclidean case, 5.7 will be a consequence of the approximate
recovery given by 2.8 for the finite Heisenberg groups

$\Gamma_{M,N}:=\Gamma\left(\left(\frac{1}{M}\mathbb{Z}\right)/N\mathbb{Z}\right)\quad\text{and}\quad
H_{N}:=\Gamma(\mathbb{Z}/N\mathbb{Z}).$

Observe that $|\Gamma_{M,N}|=M^{3}N^{3}$ and $|H_{N}|=N^{3}$.

Irreducible representations in $H_{n}$ and scaling reasoning:  In this case we
take $m|n$ and consider

$x,y,z\in\mathbb{Z}/n\mathbb{Z},\quad
a,b\in\mathbb{Z}/\frac{n}{m}\mathbb{Z},\quad c\in\mathbb{Z}/m\mathbb{Z}\text{
with }\mathsf{gcd}(c,m)=1,\quad f:\mathbb{Z}/m\mathbb{Z}\to\mathbb{C},$

and denoting $e(x):=e^{2\pi ix}$ we have representation formula

$\pi_{a,b,c}(x,y,z)f(j)=e\left(\frac{ax+by}{n}\right)e\left(\frac{c(yj+z)}{m}\right)f(j+x).$

Observe that by the definition of $\mathsf{gcd}$, the only case when $c=0$
corresponds to $m=1$; in this case, $\pi_{a,b,0}$ yields the $1$-dimensional
representation:

$\pi_{a,b}(x,y,z):=\pi_{a,b,0}(x,y,z)=e\left(\frac{ax+by}{n}\right).$

As we pass to the continuum we rescale $H_{n}\mapsto\Gamma_{M,N}$ with $n=MN$,
and we rescale the representations $\pi_{a,b,c}$ in order to recover continuum
representations $\pi_{\lambda},\pi_{\alpha,\beta}$ as in (5.27), (5.28).

Irreducible unitary representations in $\Gamma_{M,N}$:  To obtain the
quantitative recovery for the finite Heisenberg group $\Gamma_{M,N}$, we need
a description of its unitary irreducible representations.

Fix $m|M$, $c\in\mathbb{Z}/mM\mathbb{Z}$ with $\mathsf{gcd}(c,Mm)=1$, and
$a,b\in\frac{1}{N}\mathbb{Z}/\frac{M}{m}\mathbb{Z}$. Define the unitary
irreducible representation
$\pi_{a,b,c,m}:\Gamma_{M,N}\to\mathcal{U}\left(\ell^{2}\left(\frac{1}{M}\mathbb{Z}/m\mathbb{Z}\right)\right)$
by

$\pi_{a,b,c,m}(x,y,z)f(j):=e\left(ax+by\right)e\left(\frac{c(yj+z)}{m}\right)f(j+x),\qquad\text{where
}f\in\ell^{2}\left(\frac{1}{M}\mathbb{Z}/m\mathbb{Z}\right).$ (5.29)

Following Section 2.3, the (normalized) Discrete Fourier Transform of a
function $f\in\ell^{2}\left(\Gamma_{M,N}\right)$ associates to every
$a,b\in\frac{1}{N}\mathbb{Z}/\frac{M}{m}\mathbb{Z}$ and
$c\in\mathbb{Z}/mM\mathbb{Z}$ the operator over
$\ell^{2}\left(\frac{1}{M}\mathbb{Z}/m\mathbb{Z}\right)$ given by:

$\mathcal{F}_{M,N}(f)(a,b,c):=\mathcal{F}_{M,N}(f)(\pi_{a,b,c,m})=\frac{1}{M^{3}N^{3}}\sum_{(x,y,z)\in\Gamma_{M,N}}f(x,y,z)\pi_{a,b,c,m}(x,y,z);$
(5.30)

We point out that formula (5.30) only works in duality with continuous
functions supported on
$\ell^{2}\left(\frac{1}{M}\mathbb{Z}/m\mathbb{Z}\right)$. Since the Fourier
Transform $\mathcal{F}(\delta_{H})$ defined as in (5.27) and (5.28) acts in
duality with compactly supported functions, to establish a limit from
$\mathcal{F}_{M,N}(\delta_{H_{N}})$ to $\mathcal{F}(\delta_{H})$ we require
the discrete-to-continuum weak-limit in 4.1.

######  Lemma 5.8.

Let $\lambda$ be a non-zero real number and $k\in\mathbb{N}$ be fixed. For
every positive integer $m$, let
$\Gamma_{km,N}=\Gamma\left(\frac{1}{km}\mathbb{Z}/N\mathbb{Z}\right)$, and let
$H=(\mathbb{Z}^{3},\cdot)$ and
$H_{N}:=\Gamma(\mathbb{Z}/N\mathbb{Z})\leq\Gamma_{km,N}$. Let
$a,b\in\frac{1}{N}\mathbb{Z}/k\mathbb{Z}$ and
$c_{m}\in\mathbb{Z}/km^{2}\mathbb{Z}$ be a sequence such that
$\mathsf{gcd}(c_{m},km^{2})=1$ and $c_{m}/m\to\lambda$ when $m\to\infty$, then
the following limit holds

$\begin{split}\lim_{m\to\infty}(km)^{2}\mathcal{F}_{km,N}(\delta_{H_{N}})(\pi_{a,b,c_{m},m})=\frac{1}{N^{3}}\sum_{(x,y,z)\in
H\cap B_{N}}\pi_{a,b}(x,y,z)\pi_{\lambda}(x,y,z).\end{split}$ (5.31)

In particular, fix $k\in\mathbb{N}$, $\alpha,\beta\in[-k/2,k/2]$, and let
$a_{N},b_{N}\in\frac{1}{N}\mathbb{Z}/k\mathbb{Z}$ such that $a_{N}\to\alpha$
and $b_{N}\to\beta$ when $N\to\infty$. Then the following limits hold

$\begin{split}\lim_{N\to\infty}\lim_{m\to\infty}(km)^{2}\mathcal{F}_{km,N}(\delta_{H_{N}})(\pi_{0,0,c_{m},m})&=\mathcal{F}(\delta_{H})\left(\lambda\right)\\\
\lim_{N\to\infty}\lim_{m\to\infty}(km)^{2}\mathcal{F}_{km,N}(\delta_{H_{N}})(\pi_{a_{N},b_{N},0,m})&=\mathcal{F}(\delta_{H})\left(\alpha,\beta\right),\end{split}$
(5.32)

where all the limits for $m\to\infty$ are in the sense of 4.1 and the limit
with respect to $N$ is weak in the classical sense.

######  Lemma 5.9.

Let $\lambda\in\mathbb{R}\setminus\\{0\\}$, $\alpha,\beta\in\mathbb{R}$ and
$m_{0}\in\mathbb{N}$ large enough such that $\lambda\in[-m_{0}/2,m_{0}/2]$.
Consider $a_{N},b_{N}\in\frac{1}{N}\mathbb{Z}/k\mathbb{Z}$ and $c_{m}$ be as
in Lemma 5.8. Let $\xi^{(km,N)}$ be as in Section 4.0.1 and suppose there
exists a positive number $\varepsilon$ such that
$\mathbb{E}[|\xi|^{3+\varepsilon}]<\infty$. Then the following weak limits
hold almost surely :

$\begin{split}\lim_{N\to\infty}\lim_{m\to\infty}(km)^{2}\mathcal{F}_{km,N}\left(\delta_{H_{\xi^{(km,N)}}}\right)(\pi_{0,0,c_{m},m})&=\mathcal{F}(\delta_{H_{\xi}})(\lambda)\\\
\lim_{N\to\infty}\lim_{m\to\infty}(km)^{2}\mathcal{F}_{km,N}\left(\delta_{H_{\xi^{(km,N)}}}\right)(\pi_{a_{N},b_{N},0,m})&=\mathcal{F}(\delta_{H_{\xi}})(\alpha,\beta).\end{split}$
(5.33)

The last ingredient to prove the main recovery result on the Heisenberg group
is that the error given by 2.8 for finite Heisenberg groups converges to zero
when considering a sequence of discrete Heisenberg-torii $\Gamma_{M,N}$, for
$M,N\to\infty$.

######  Proposition 5.10.

Let $N,m,k\in\mathbb{N}$ and consider finite Heisenberg groups $\Gamma_{km,N}$
and $H_{N}$ as before. Let $\xi_{p}^{(N)}=\xi_{p},p\in H_{N}$ and
$\varepsilon>0$ satisfying the hypotheses of 2.8. Then with probability at
least $1-\frac{C}{\varepsilon^{2}|H|}$ for some suitable positive constant $C$
there holds:

$\left\|\mathcal{F}_{km,N}(\delta_{H_{\xi}})-\mathcal{F}_{km,N}(\delta_{H})\mathbb{E}(\xi(\cdot))\right\|^{2}_{L^{2}(\widehat{\Gamma}_{km,N})}\leq\frac{(\varepsilon+4C^{2})N^{7}}{(km)^{5}}.$
(5.34)

###### Proof.

From 2.8 and since for the finite Heisenberg groups $\Gamma_{km,N}$ there
holds that $\max\\{d_{\pi}:\
\pi\in\widehat{\Gamma}_{km,N}\\}=kmN=|\Gamma_{km,N}|^{1/3}$, then for
$\frac{C^{\prime}}{|H|}<\varepsilon^{2}$ we get:

$\begin{split}\left\|\mathcal{F}_{km,N}(\delta_{H_{\xi}})-\mathcal{F}_{km,N}(\delta_{H})\mathbb{E}(\xi(\cdot))\right\|^{2}_{L^{2}(\widehat{\Gamma}_{km,N})}&\leq\frac{|H|^{4}}{|\Gamma_{km,N}|^{2}}\left(\varepsilon+\frac{4C^{2}}{|H|}\right)\max_{\pi\in\widehat{\Gamma}_{km,N}}d_{\pi}\\\
&=\frac{\varepsilon N^{7}}{(km)^{5}}+\frac{4C^{2}N^{4}}{(km)^{5}}.\end{split}$
(5.35)

To conclude, it suffices to use the fact that
$\frac{N^{4}}{(km)^{5}}\leq\frac{N^{7}}{(km)^{5}}$. ∎

Hence, from (4.20), (4.21) and 5.10 we get the next result which in fact
implies (4.19).

######  Corollary 5.11.

Under the hypotheses and notations of 5.10, let $\xi_{p}^{(km,N)},p\in H_{N}$
as in (4.0.1). Then, with probability at least
$1-\frac{C}{\varepsilon^{2}|H|}$, for every
$\lambda\in\mathbb{R}\setminus\\{0\\}$,
$a,b\in\frac{1}{N}\mathbb{Z}/k\mathbb{Z}$ and every sequence
$c_{m}\in\mathbb{Z}/km^{2}\mathbb{Z}$, $\mathsf{gcd}(c_{m},km^{2})=1$ for
which $c_{m}/m\to\lambda$ when $m\to\infty$, there holds that:

$(km)^{2}\left\|\mathcal{F}_{km,N}(\delta_{H_{\xi^{(km,N)}}})(\pi)\varphi-\mathcal{F}_{km,N}(\delta_{H_{N}})\mathbb{E}(\pi(\xi^{(km,N)}))\varphi\right\|_{\ell^{2}\left(\frac{1}{km}\mathbb{Z}/m\mathbb{Z}\right)}\longrightarrow
0,$

where $\pi$ is of the form $\pi_{0,0,c_{m},m}$ or $\pi_{a,b,0,m}$ and
$\phi\in\ell^{2}\left(\frac{1}{km}\mathbb{Z}/m\mathbb{Z}\right)$.

###### Proof of 5.7.

It is a direct application of Lemma 5.8, Lemma 5.9 and Corollary 5.11.

∎

Acknowledgements: _M. P. was sponsored by the Chilean Fondecyt Regular grant
number 1210462 entitled “Rigidity, stability and uniformity for large point
configurations”, and R. V. was sponsored by the Chilean Fondecyt Postdoctoral
grant number 3210109 entitled “Geometric and Analytical aspects of Discrete
Structures”._

## References

  * [1] G. Alagic, D. Apon, D. Cooper, Q. Dang, T. Dang, J. Kelsey, J. Lichtinger, C. Miller, D. Moody, R. Peralta, Status report on the third round of the NIST post-quantum cryptography standardization process. US Department of Commerce, NIST 2022\.
  * [2] M. Baake and U. Grimm, Aperiodic Order, vol. 1: A mathematical invitation, first ed., Encyclopedia of Math. and its Appl., vol. 149, Cambridge Univ. Press, 2013.
  * [3] H. Bahouri, D. Barilari, I. Gallagher, M. Léautaud, Spectral summability for the quartic oscillator with applications to the Engel group. arXiv preprint arXiv:2206.10396.
  * [4] I. Benjamini, H. Finucane, R. Tessera, On the scaling limit of finite vertex graphs with large diameter. Combinatorica, 37(3), 333-374, 2017.
  * [5] U. Boscain, J-P. Gauthier, F. Rossi, Hypoelliptic heat kernel over 3-step nilpotent Lie groups. J. Math. Sci., 199, 614-628, 2014.
  * [6] A.H. Clifford, Representations Induced in an Invariant Subgroup. Annals of Mathematics, 38, no. 3 (1937): 533–550.
  * [7] L. Corwin, F. Greenleaf, Representations of nilpotent Lie groups and their applications: Volume 1, Part 1, Basic theory and examples. Cambridge university press, 1990.
  * [8] J. Dixmier, $C^{*}$ algebras. North-Holland Publishing Co. Translated from the French by Francis Jellett, North-Holland Mathematical Library, Vol. 15., 1977.
  * [9] V. Fischer, M. Ruzhansky, Quantization on nilpotent Lie groups. Springer Nature, 2016.
  * [10] G. Folland,A course in Abstract Harmonic Analysis. Studies in Advanced Mathematics, 1995.
  * [11] A. Hof, Diffraction by aperiodic structures at high temperatures, J. Phys. A, vol. 28, 1995.
  * [12] A.A. Kirillov, Unitary representations of nilpotent Lie groups. Russ. Math. Surv., 17 53, no. 4 (1962).
  * [13] M. A. Klatt, J. Kim, S. Torquato, Cloaking the underlying long-range order of randomly perturbed lattices, Phys. Rev. E, vol. 101, 032118, 2020.
  * [14] N. Lev, R. Peled and Y. Peres, Separating signal from noise, Proc. London Math. Soc, 110(4), 883-931, 2015.
  * [15] P. Pansu, Croissance des boules et des géodésiques fermées dans les nilvariétés., Erg. Theory and Dyn. Syst., 3(3),415-445, 1983.
  * [16] G. Pete, Probability and geometry on groups. Lecture notes for a graduate course. Available on `https://math.bme.hu/~gabor/PGG.pdf` (Last revision, June 2023) (2017).
  * [17] M. Petrache, R. Viera, Almost sure recovery in quasi-periodic structures. J. Stat. Phys., 190(2), 2023.
  * [18] T. Tao, Lecture notes of Representation Theory on finite non-abelian groups. Available on https://www.math.ucla.edu/ tao/247b.1.07w/.
  * [19] A. Terras, Fourier analysis on finite groups and applications. No. 43. Cambridge University Press, 1999.
  * [20] P. Woit, Quantum theory, groups and representations. Springer International Publishing, 2017.
  * [21] O. Yakir, Recovering the Lattice From Its Random Perturbations, Int. Math. Res. Notices,2022(8) 6243-6261, 2022.