# Highly resolved spectral functions of two-dimensional systems with neural
quantum states

Tiago Mendes-Santos  Theoretical Physics III, Center for Electronic
Correlations and Magnetism, Institute of Physics, University of Augsburg,
86135 Augsburg, Germany  Markus Schmitt  Forschungszentrum Jülich GmbH, Peter
Grünberg Institute, Quantum Control (PGI-8), 52425 Jülich, Germany  Markus
Heyl  Theoretical Physics III, Center for Electronic Correlations and
Magnetism, Institute of Physics, University of Augsburg, 86135 Augsburg,
Germany

###### Abstract

Spectral functions are central to link experimental probes to theoretical
models in condensed matter physics. However, performing exact numerical
calculations for interacting quantum matter has remained a key challenge
especially beyond one spatial dimension. In this work, we develop a versatile
approach using neural quantum states to obtain spectral properties based on
simulations of the dynamics of excitations initially localized in real or
momentum space. We apply this approach to compute the dynamical structure
factor in the vicinity of quantum critical points (QCPs) of different two-
dimensional quantum Ising models, including one that describes the complex
density wave orders of Rydberg atom arrays. When combined with deep network
architectures we find that our method reliably describes dynamical structure
factors of arrays with up to $24\times 24$ spins, including the diverging time
scales at critical points. Our approach is broadly applicable to interacting
quantum lattice models in two dimensions and consequently opens up a route to
compute spectral properties of correlated quantum matter in yet inaccessible
regimes.

Introduction. Spectral functions are key tools to characterize and probe
quantum many-body phases and their transitions. In addition, they serve as a
common framework to connect theoretical descriptions with experimental probes
such as photoemission or inelastic neutron scattering. In this context, a
regime of particular interest is two-dimensional (2D) interacting quantum
matter, where experimental probes can indicate the occurrence of prominent
properties, such as exotic fractionalized quasiparticles in candidate
materials realizing 2D spin-liquid phases [1] or universal features associated
to quantum critical points [2].

Figure 1:  Panel (a) shows a schematic picture of the NQS architecture used in
this work. Besides the convolutional neural network
$\psi_{\mathbf{\theta}}(\mathbf{s})$ that forms the variational part, the
initial excitations in real or momentum space ($\mathbf{r}$ or $\mathbf{q}$)
are captured by directly incorporating the corresponding operator matrix
elements
$f_{\mathbf{r}/\mathbf{q}}(\mathbf{s})=\braket{\mathbf{s}}{\hat{\sigma}^{z}_{\mathbf{r}/\mathbf{q}}}{\mathbf{s}}$
in the ansatz. The logarithmic wave function coefficients of the state
$\ket{\phi^{z}_{\boldsymbol{r}/\boldsymbol{q}}}$ then take the form
$\log\braket{\mathbf{s}}{\phi^{z}_{\boldsymbol{r}/\boldsymbol{q}}}=\log\psi_{\mathbf{\theta}}(\mathbf{s})+\log
f_{\mathbf{r}/\mathbf{q}}(\mathbf{s})$ as indicated by the computational
graph. Panels (b1,c1) show some exemplary results of real-space dynamical
correlators and the corresponding DSF close to the quantum critical point of
the two-dimensional quantum Ising model. Analogously, panels (b2,c2) show
results obtained directly in momentum space for larger system sizes. We
consider the path $X=(\pi,0)\to M=(\pi,\pi)\to\Gamma=(0,0)\to X\to
S=(\pi/2,\pi/2)$ of the Brilloin Zone.

At the theoretical level, accessing and describing spectral functions is,
thus, of great interest in strongly interacting solid-state materials. But
addressing dynamical properties of correlated matter in a controlled manner
poses, at the same time, substantial challenges. Quantum Monte Carlo is poised
by a sign problem [3] and the applicability of dynamical mean field theory [4]
is limited in low dimensions. Tensor network approaches, which render the
treatment of weakly entangled states feasible, can be used to obtain
numerically exact results for one-dimensional systems [5, 6, 7, 8, 9, 10].
While extensions to higher dimensions exist [11, 12, 13, 14, 15], the growth
of entanglement in time together with the two-dimensional lattice structure
that increases the complexity of tensor contractions remains as a challenge
for tensor network methods. In addition, variational methods for capturing
excitations based on Gutzwiller-projected mean-field states are restricted to
specific cases due to their built-in bias [16, 17]. Finally, programmable
quantum simulation could emerge as a new route [18, 19, 20], but it is still
in its infancy at this point.

Recently, the idea to combine the variational Monte Carlo (VMC) framework with
neural quantum states (NQSs) [21] has been shown to be very fruitful for
investigations of correlated matter, including the simulation of ground states
of frustrated Hamiltonians [22, 23, 24, 25, 26, 27] and the dynamics of two-
dimensional systems [28, 29, 30, 31, 32, 33]. For spectral functions, first
attempts proposed NQS-based algorithms built directly in the frequency domain
[34, 35, 36] or a method simulating the response to an initial time-dependent
perturbation of the system [37]. However, it remains desirable to enhance the
resolution and reachable system sizes over what has been achieved so far in
order to address open physical questions.

In this work, we introduce an alternative versatile scheme for the simulation
of spectral functions based on the direct encoding of local excitations in the
neural network architecture – the specNQS, see Fig. 1(a). The spectral
information is then extracted from dynamical correlation functions that are
obtained by real-time evolution. When combined with convolutional neural
networks, we demonstrate that our scheme allows us to access dynamical
properties beyond what has been feasible with other state-of-the-art
approaches. As a benchmark, we simulate the dynamical structure factor (DSF)
of the 2D quantum Ising model (QIM), and we showcase that the specNQS reliably
describes spectral features associated with a diverging correlation length for
system sizes up to $24\times 24$ sites. Furthermore, we contribute to the
characterization of quantum phase transitions in experimentally realized long-
range interacting Rydberg atom arrays [38, 39] by revealing spectral
properties close to phase boundaries, the nature of which is under ongoing
debate [40, 41].

Method. In the following, we will be interested in computing the DSF

$\displaystyle S^{zz}(\boldsymbol{q},\omega)$
$\displaystyle=\frac{1}{N_{s}}\int_{-\infty}^{\infty}dte^{i\omega
t}\big{\langle}\hat{\sigma}^{z}_{-\boldsymbol{q}}(t)\hat{\sigma}^{z}_{\boldsymbol{q}}\big{\rangle}.$
(1)

Here,
$\langle\hat{A}(t)\hat{B}\rangle=\braket{\psi_{0}}{e^{i\hat{H}t}\hat{A}e^{-i\hat{H}t}\hat{B}}{\psi_{0}}$
denotes the dynamical correlation function in the ground state
$\ket{\psi_{0}}$ of a given Hamiltonian $\hat{H}$,
$\hat{H}\ket{\psi_{0}}=E_{0}\ket{\psi_{0}}$ with $E_{0}$ being the ground-
state energy.
$\hat{\sigma}^{z}_{\boldsymbol{q}}=\sum_{\boldsymbol{r}}e^{-i\boldsymbol{q}\cdot\boldsymbol{r}}\hat{\sigma}^{z}_{\boldsymbol{r}}$
is the spin operator in momentum space, where $\hat{\sigma}^{z}_{\mathbf{r}}$
denotes the Pauli-$z$ operator at lattice site $\boldsymbol{r}$ and $N_{s}$ is
the number of lattice sites.

The central idea of our approach is to obtain the DSF from a variational
representation of suitably time-evolved wave functions. NQSs constitute a
versatile family of variational wave functions relying on the proven
representational power of artificial neural networks (ANNs). In particular,
any function can be accurately approximated by an ANN in the limit of large
network sizes [42, 43, 44, 45]. This means that the accuracy of the proposed
approach can be asserted self-consistently by checking the convergence with
increasing network size.

As a first step to access the DSF, we compute the NQS representation of the
ground state

$\ket{\psi_{0}}=\sum_{\boldsymbol{s}}\psi_{\boldsymbol{\theta}}^{(0)}(\boldsymbol{s})\ket{\boldsymbol{s}},$
(2)

where $\boldsymbol{s}=(s_{1},s_{2},...,s_{N_{s}})$ labels the Pauli-Z basis of
spin configurations. The variational ansatz
$\psi_{\boldsymbol{\theta}}^{(0)}(\boldsymbol{s})$ is parameterized by
$\boldsymbol{\theta}=(\theta_{1},\theta_{2},...,\theta_{M})$ and it takes the
form of an ANN. The ground state is then obtained by optimizing
$\boldsymbol{\theta}$ to minimize the energy expectation value
$\mathcal{E}(\boldsymbol{\theta})=\braket{\psi(\boldsymbol{\theta})}{\hat{H}}{\psi(\boldsymbol{\theta})}/\braket{\psi(\boldsymbol{\theta})}{\psi(\boldsymbol{\theta})}$.
For the results presented throughout this manuscript, we employed the
Stochastic Reconfiguration algorithm to find the ground states [46, 47].

Our approach to access the dynamics relies on computing the time-evolved wave
functions following an excitation,

$\displaystyle\ket{\phi_{\boldsymbol{r}/\boldsymbol{q}}^{\alpha}(t)}=e^{-i\hat{H}t}\hat{\sigma}_{\boldsymbol{r}/\boldsymbol{q}}^{\alpha}\ket{\psi_{0}}\equiv
e^{-i\hat{H}t}\ket{\phi_{\boldsymbol{r}/\boldsymbol{q}}^{\alpha}(0)}\ .$ (3)

Here, $\hat{\sigma}_{\boldsymbol{r}/\boldsymbol{q}}^{\alpha}$ is an operator
in either position ($\boldsymbol{r}$) or momentum space ($\boldsymbol{q}$). In
our numerical approach we rely in either case on an exact representation of
$\ket{\phi_{\boldsymbol{r}/\boldsymbol{q}}^{\alpha}(0)}$, which will be
discussed in the following paragraphs. Let us first discuss the time-evolution
algorithm, assuming that the representation of the initial state is given. We
employ a time-dependent variational principle (TDVP) that is based on the
minimization of the Fubini-Study distance between the variational time-evolved
state $\ket{\psi_{\boldsymbol{\theta}(t+\delta t)}}$ and the exact one
$e^{-i\hat{H}\delta t}\ket{\psi(t)}$, where $\delta t$ is an infinitesimal
time interval. It yields an ordinary non-linear differential equation
prescribing the optimal evolution of the variational parameters [21, 28],

$\displaystyle S_{k,k^{\prime}}\dot{\theta}_{k^{\prime}}=-iF_{k}\ ,$ (4)

where $\dot{\theta}_{k^{\prime}}$ is the time derivative,
$S_{k,k^{\prime}}=\braket{\partial_{\theta_{k}}\psi_{\boldsymbol{\theta}}}{\partial_{\theta_{k^{\prime}}}\psi_{\boldsymbol{\theta}}}-\braket{\partial_{\theta_{k}}\psi_{\boldsymbol{\theta}}}{\psi_{\boldsymbol{\theta}}}\braket{\psi_{\boldsymbol{\theta}}}{\partial_{\theta_{k^{\prime}}}\psi_{\boldsymbol{\theta}}}$
is the quantum metric tensor and
$F_{k}=\braket{\partial_{\theta_{k}}\psi_{\boldsymbol{\theta}}}{\hat{H}}{\psi_{\boldsymbol{\theta}}}-\braket{\partial_{\theta_{k}}\psi_{\boldsymbol{\theta}}}{\psi_{\boldsymbol{\theta}}}\braket{\psi_{\boldsymbol{\theta}}}{\hat{H}}{\psi_{\boldsymbol{\theta}}}$.
Both $S_{k,k^{\prime}}$ and $F_{k}$ can be estimated efficiently via Monte-
Carlo sampling of the Born distribution $|\psi_{\boldsymbol{\theta}}(s)|^{2}$.

Upon integration, the TDVP in Eq. (4) in general yields the time-evolved state
up to a global phase and normalization [48]. The phase is irrelevant for
equal-time correlation functions. For computing the dynamical correlation
functions, however, it becomes important, because we are interested in
evaluating the overlap of two time-evolved states. To keep track of relative
changes in phase between such time-evolved states we consider the
(logarithmic) prefactor $\theta_{0}$ as an additional variational parameter:
$\ket{\psi_{\theta_{0},\boldsymbol{\theta}}}=e^{\theta_{0}}\ket{\psi_{\boldsymbol{\theta}}}$.
By using the TDVP, one can establish the following equation of motion for
$\theta_{0}$ [49]

$\dot{\theta}_{0}=-i\braket{\hat{H}}-\dot{\theta}_{k}\braket{\psi_{\boldsymbol{\theta}}}{\partial_{\theta_{k}}\psi_{\boldsymbol{\theta}}}.$
(5)

Equation (4) for the other parameters $\boldsymbol{\theta}$ remains unchanged.
Thus, for each time step, we first obtain $\dot{\theta}_{k}$ and then use the
result to solve Eq. (5) for the evolution of $\theta_{0}$.

Momentum-space scheme. First, we discuss how to access the dynamical spin
structure factor directly in momentum space,

$S^{zz}(\boldsymbol{q},t)=e^{iE_{0}t}\bra{\psi_{0}}\hat{\sigma}_{\boldsymbol{q}}^{z}e^{-i\hat{H}t/2}e^{-i\hat{H}t/2}\hat{\sigma}_{\boldsymbol{q}}^{z}\ket{\psi_{0}}\
.$ (6)

To simplify the discussion we focus on the $z$ component, but our approach can
be straightforwardly generalized for $\alpha=x,y$ by choosing the
computational basis accordingly. The central idea is that the action of the
operator $\hat{\sigma}_{\boldsymbol{q}}^{z}$ on the initial state
$\ket{\psi_{0}}$ can be captured explicitly and efficiently by modifying the
individual wave function coefficients with corresponding prefactors.
Concretely, the excitation $\sigma_{\boldsymbol{q}}^{z}\ket{\psi_{0}}$ is
encoded on the NQS ansatz by adding a configuration-dependent factor
$f(\boldsymbol{q},\boldsymbol{s})=\sum_{j}e^{-i\boldsymbol{r}_{j}\boldsymbol{q}}e^{i\pi/2(s_{j}-1)}$
on the initial-time quantum state

$\ket{\phi_{\boldsymbol{q}}^{z}(0)}=\hat{\sigma}_{\boldsymbol{q}}^{z}\ket{\psi_{0}}=\sum_{\boldsymbol{s}}f(\boldsymbol{q},\boldsymbol{s})\psi_{\boldsymbol{\theta}}(\boldsymbol{s})\ket{\boldsymbol{s}}\
.$ (7)

We then obtain $S^{zz}(\boldsymbol{q},t)$ by performing a two-sided time
evolution with the TDVP approach, followed by evaluating the overlap of the
two time evolved states,

$S^{zz}(\boldsymbol{q},t)=e^{iE_{0}t}\braket{\phi^{z}_{\boldsymbol{q}}(-t/2)}{\phi^{z}_{\boldsymbol{q}}(t/2)}.$
(8)

This is a the central object to calculate the dynamical structure factor by
means of our NQS approach. Notice that this implies that for the dynamical
correlation function up to time $t$ numerical integration is only required up
to time $t/2$. Here, it is important to account for the global factor
associated to each state of the overlap as discussed above, see Eq. (5). Refs
[50, 51] presents further details about the calculation of the overlap.

Real-space scheme. Second, we discuss a strategy to obtain dynamical
correlations in real space. Our scheme is based on a many-body Ramsey protocol
[18], that is used to simulate the retarded Green’s function (GF)

$G^{\alpha\alpha}_{\boldsymbol{r},\boldsymbol{r}^{\prime}}(t)=-\frac{i}{2}\bra{\psi_{0}}[\hat{\sigma}_{\boldsymbol{r}}^{\alpha}(t),\hat{\sigma}_{\boldsymbol{r}^{\prime}}^{\alpha}(0)]\ket{\psi_{0}},$
(9)

where $\boldsymbol{r}$, $\boldsymbol{r}^{\prime}$ are the sites of a lattice,
and $[\hat{A},\hat{B}]$ denotes the commutator. In particular, to access the
longitudinal $G^{zz}_{\boldsymbol{r},\boldsymbol{r}^{\prime}}(t)$ component,
we start the protocol with the following quantum state
$\ket{\phi_{\boldsymbol{r}}^{z}(0)}=e^{i\frac{\pi}{4}\hat{\sigma}^{z}_{\boldsymbol{r}}}\ket{\psi_{0}}$,
where the local pertubation is represented by the $\frac{\pi}{4}$-rotation of
a spin at site $\boldsymbol{r}$. Further, we obtain the time-evolved state
$\ket{\phi_{\boldsymbol{r}}^{z}(t)}=e^{-i\hat{H}t}\ket{\phi_{\boldsymbol{r}}^{z}(0)}$.
Following Ref. [18], the result of a local measurement of
$\hat{\sigma}_{\boldsymbol{r}}^{z}$ at a time $t$ depends on the GF, i.e.,

$\displaystyle\bra{\phi_{\boldsymbol{r}^{\prime}}^{z}(t)}\hat{\sigma}_{\boldsymbol{r}}^{z}\ket{\phi_{\boldsymbol{r}^{\prime}}^{z}(t)}=$
$\displaystyle\bra{\psi_{0}}\hat{\sigma}_{\boldsymbol{r}}^{z}(t)\ket{\psi_{0}}+G^{zz}_{\boldsymbol{r},\boldsymbol{r}^{\prime}}(t)+\bra{\psi_{0}}\hat{\sigma}_{\boldsymbol{r}}^{z}\hat{\sigma}_{\boldsymbol{r}^{\prime}}^{z}(t)\hat{\sigma}_{\boldsymbol{r}}^{z}\ket{\psi_{0}}.$
(10)

The $G^{zz}_{\boldsymbol{r},\boldsymbol{r}^{\prime}}(t)$ is obtained by
reconstructing the terms of the Eq. (10). The first term on the right hand
side is accessed from the ground state. For the remaining contributions we
time-evolve the initial states
$e^{i\frac{\pi}{4}\hat{\sigma}^{z}_{\boldsymbol{r}}}\ket{\psi_{0}}$ and
$\hat{\sigma}_{\boldsymbol{r}}^{z}\ket{\psi_{0}}$ after incorporating the
operator action into the variational ansatz in analogy to Eq. (7).

One central difference between the momentum- and the real-space approach is
that the latter does not require the calculation of state overlaps. Moreover,
they differ in the way translational symmetry can be exploited for the
efficiency of simulations. In the real space approach all correlation
functions in Eq. (10) depend only on relative positions
$\boldsymbol{r}^{\prime}-\boldsymbol{r}$, which means that all the momentum
points of $S^{zz}(\boldsymbol{q},\omega)$ can be obtained from the two time-
evolved states $\ket{\phi_{\boldsymbol{r}=0}^{z}(t)}$ and
$e^{-i\hat{H}t}\hat{\sigma}_{\boldsymbol{r}=0}^{z}\ket{\psi_{0}}$. However,
with this approach, translational symmetry cannot be built into the
variational ansatz to enhance efficiency. By contrast, the time-evolved states
in Eq. (8) for the momentum-space approach preserve translational symmetry,
which can be exploited to introduce beneficial bias through built-in
invariance of the NQS. This comes, however, at the cost of individual
simulations required for each point in momentum space.

Finally, it is worth mentioning that we compute
$S^{zz}(\boldsymbol{q},\omega)$ by performing a Fourier transform with a
Gaussian envelope to avoid the finite-time effects of our simulation (our
simulations are performed up to a time $t_{max}$) [50].

Neural quantum states architectures. In conjunction with the momentum-space
scheme, we employ convolutional neural networks (CNNs) [52, 28] as the
variational part of the NQS architecture, which allows us to exploit the
translational symmetry, that is also conserved after applying the operators in
momentum space. The hyperparameters of the CNNs are the total number of layers
$l$, the number of channels in each layer $k$, $\alpha_{k}$, and the linear
size of the square filter, $F$; in the following, we characterize the CNN
architecture with the tuples
$\boldsymbol{\alpha}=(\alpha_{1},...,\alpha_{l};F)$. Meanwhile, to implement
the real-space scheme, we use Restricted Boltzmann Machines (RBMs) [21]
composed of a single fully-connected hidden layer with $M$ nodes, where $M$ is
a hyperparameter; further details about the NQSs are discussed in the SM [50].

Figure 2: 2D quantum Ising model, vicinity of the quantum critical point.
Panel (a) shows the finite-size scaling of the gap
$\Delta(\boldsymbol{q}=\Gamma)$; the star points are series-expansion results
up to forth order, while the triangle points are exact-diagonalization results
for $L=4$. In panel (b), we show the corresponding DSF
$S(\boldsymbol{q}=\Gamma,\omega)$ at QCP for different values of $L$. In panel
(c) we present $S(\boldsymbol{q},\omega)$ at the ferromagnetic side of the
quantum phase transition.

Figure 3: Two-dimensional Rydberg atom arrays. In panel (a), we show a
schematic phase diagram of the RyM and the color plot of the spin correlation
function,
$C^{z}(\boldsymbol{r})=\left<\left<\sigma^{z}_{\boldsymbol{0}}\sigma^{z}_{\boldsymbol{r}}\right>\right>$,
characteristic of CbDW and StDW phases; we consider (a1) $R_{b}=1.2$,
$\delta=1.05$ and (a2) $R_{b}=1.45$, $\delta=2.00$. In panels (b1) and (b2),
we show the corresponding DSF for $L=12$ along the following path of the
Brillouin zone: $S=(\pi/2,\pi/2)\to M=(\pi,\pi)\to X=(\pi,0)$. The circle
(red) points mark the peak of the DSF for each value of $\boldsymbol{q}$.
Finally, considering the same set of parameters $R_{b}$ and $\delta$, panels
(c1.1,c1.2) and (c2.1,c2.2) show the respective finite-size scaling of the DSF
and the dynamical correlation at the $\boldsymbol{q}=\boldsymbol{M}$ point;
the dashed lines are exact results for $L=4$.

Results I: Two-dimensional Quantum Ising model. To benchmark our approach, we
consider the paradigmatic 2D QIM on a square lattice:

$\hat{H}=-\sum_{\left<i,j\right>}\hat{\sigma}^{z}_{i}\hat{\sigma}^{z}_{j}-g\sum_{i}\hat{\sigma}^{x}_{i}.$
(11)

This model describes a second-order QCP at $g_{c}\approx 3.044$, separating a
ferromagnetic from a paramagnetic phase. In the vicinity of $g_{c}$, the major
contributions to $S^{zz}(\boldsymbol{q},\omega)$ come from low-energy
quasiparticle excitations whose frequency and spectral weight are expected to
scale as $\omega\sim L^{-z}$ and $\tilde{S}\sim L^{(1-\eta)}$, respectively
[2]; here, $z=1$ and $\eta\approx 0.04$ are critical exponents of the 3D-Ising
universality class [53].

We start discussing the DSF obtained with the real-space scheme shown in Fig.
1(c1). The plots show the DSF along a path through the first Brillouin zone
for a lattice of size $N_{s}=10\times 10$ at transverse field $g=3.1$. The
main feature of $S^{zz}(\boldsymbol{q},\omega)$ is a minimum at the $\Gamma$
point, corresponding to the low-energy gap. We compare our results to series-
expansion results up to fourth order in $g^{-1}$ [54, 55] and find very good
agreement almost everywhere in the Brillouin zone. The most notable deviation
appears at the gap closing point when approaching the critical value $g_{c}$,
which we will discuss in more detail below. Figure 1 (c2) shows similar
results obtained with the momentum scheme in the vicinity of the $\Gamma$
point. The fact that we can in this case use more efficient NQS architectures
with built-in translational invariance allows us to simulate substantially
larger system sizes up to $N_{s}=20\times 20$ lattice sites and we again find
good agreement with the series expansion.

Let us now focus on the behavior of the DSF in the vicinity of the QCP. As a
demanding benchmark for the accuracy of the NQS approach, we investigate the
finite-size scaling of the spectral gap extracted from the DSF,
$\Delta(\boldsymbol{q}=\Gamma)=\text{argmax}_{\omega}S^{zz}(\boldsymbol{q}=\Gamma,\omega)$.
The expected universal scaling behavior $\Delta(\boldsymbol{q}=\Gamma)\sim
L^{-z}$, with $L$ the linear system size, can be compellingly confirmed by our
finite-size simulations up to $N_{s}=24\times 24$; we obtain $z=1.00\pm 0.05$
by fitting such universal scaling, see SM [50]. We emphasize that our results
show that with our NQS approach such large system sizes required to extract
this universal behavior have now become within reach, while simultaneously
also the diverging time scale associated with the closing of the spectral gap
can be captured.

In Fig. 2(c) we moreover show the frequency-resolved DSF at a few momentum
points at $g=2.9$, i.e. on the ferromagnetic side of the QCP. These cuts
reveal a double peak structure of $S^{zz}(\boldsymbol{q},\omega)$, in
particular an excitation at very low energy for $\boldsymbol{q}=\Gamma$. This
feature is a signature of the $\mathbb{Z}_{2}$ symmetry-broken phase, where
the degeneracy of the ground state is lifted to an exponentially small energy
gap due to the finite system size.

Results II: Rydberg atom arrays. Now, to take our approach to a next level,
beyond the QIM benchmark, we consider a long-range interacting model
describing Rydberg atoms arrays on a square lattice (RyM). The RyM is defined
as

$\hat{H}=-\frac{\Omega}{2}\sum_{i}\hat{\sigma}^{x}_{i}-\delta\sum_{i}\hat{n}_{i}+\sum_{i,j}\left(\frac{R_{b}}{r_{i,j}}\right)^{6}\hat{n}_{i}\hat{n}_{j},$
(12)

where $\hat{n}_{i}\equiv(\hat{\sigma}^{z}_{i}+1)/2$. The parameter $\Omega$
represents the Rabi frequency, while $\delta$ denotes the detuning. The term
$\left(R_{b}/r_{i,j}\right)^{6}$ represents the interaction between atoms in
Rydberg states, $r_{i,j}$ being the distance between atoms at sites $i$ and
$j$, and $R_{b}$ the so-called Rydberg blockade radius 111Here, we consider a
finite cutoff for the interaction potential such that the interactions are set
to $0$ for $r_{i,j}>r_{0}=3a$..

Recently, quantum simulation experiments [38, 39] have motivated theoretical
efforts to understand ground-state properties of the square-lattice RyM [57,
40, 41]. The interplay between the parameters $\delta$ and $R_{b}$ leads to
density-wave phases and related QCPs. As an example, we present a schematic
phase diagram for RyM in Fig. 3(a), indicating the emergence of checkerboard
(CbDW) and striated (StDW) density-wave. Ground state searches and the
estimations for the disordered-CbDW and disordered-StDW transitions are
discussed in the SM [50], which is consistent with the literature [41]. Here,
we focus in computing spectral functions in the vicinity of QCPs.

Near the disordered-CbDW QCP, $S^{zz}(\boldsymbol{q},\omega)$ is characterized
by a dominant low energy mode at $\boldsymbol{q}=\boldsymbol{M}$, see Fig. 3
(b1). The finite-size scaling of the DSF (we obtain $z=0.9\pm 0.2$ by fitting
the size scaling of the the lowest-energy gap, see SM [50]) and the dynamical
correlator $S^{zz}(\boldsymbol{M},t)$ is consistent with a second-order QCP in
the same universality class as the previously discussed 2D QIM; see Fig. 3
(c1.1) and (c1.2).

Close to the disordered-StDW QCP, however, $S^{zz}(\boldsymbol{q},\omega)$
exhibits qualitatively different behavior. It is characterized by two dominant
low energy modes, occurring at $\boldsymbol{M}$ and $\boldsymbol{X}$; see Fig.
3 (b2). In addition, the spectral weight associated with the lower-energy peak
of $S^{zz}(\boldsymbol{q}=\boldsymbol{M},\omega)$ decreases with system size,
while we observe a spectral weight transfer to higher energies for larger
values of $L$; see Fig. 3 (c2.1). These results are consistent with the
prediction of a first-order QCP [41]. We note, however, that for $L=12$ our
simulations cannot resolve a lower energy peak. For first-order phase
transitions, the lowest energy gap is expected to vanish exponentially as $L$
increases [58]. This implies that to resolve the spectral gap accurately, we
have to perform simulations up to an exponentially long time, $t_{max}$, and
with a Fourier transform done with a Gaussian broadening factor scaling as
$1/t_{max}$ [50]. In this regime, the phase transition is better evidenced by
ground-state static properties related to the DSF [41].

Discussion and conclusions. In summary, we have proposed a cutting-edge method
to simulate spectral properties of 2D quantum many-body systems, which relies
on the representational power of NQSs to access spectral functions during the
real-time dynamics of local excitations. We demonstrated that this scheme
allows us to reliably perform finite-size scaling of dynamical properties near
2D quantum critical points for unprecedented system sizes. A promising future
direction is to characterize spectral features of 2D models stabilizing
quantum spin liquid phases. Of particular current interest are spin liquids
with Ising-like interactions, such as the 2D RyM in a ruby lattice, which has
recently been proposed as a way to realize $\mathbb{Z}_{2}$ topological order
states in programmable quantum simulators [59, 60], and also Kitaev-type spin
models [61], where highly-accurate simulations of DSF are essential to
characterize exotic fractionalized excitations [62, 63].

###### Acknowledgements.

### Acknowledgements. -

We thank Matteo Rizzi, Michael Knapp and Ao Chen for fruitful discussions.
This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No. 853443). MS was supported through the Helmholtz
Initiative and Networking Fund. The authors gratefully acknowledge the Gauss
Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project
by providing computing time through the John von Neumann Institute for
Computing (NIC) on the GCS Supercomputer JUWELS [64] at Jülich Supercomputing
Centre (JSC). We used the jVMC codebase [65] that is built on the JAX library
[66] to implement our approach. We use the quspin package [67, 68] to obtain
the exact diagonalization results. The data shown in the figures are available
on Zenodo [69].

## References

  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Sachdev [2011] S. Sachdev,  _Quantum Phase Transitions_, 2nd ed. (Cambridge University Press, 2011).
  * Troyer and Wiese [2005] M. Troyer and U.-J. Wiese, Computational complexity and fundamental limitations to fermionic quantum monte carlo simulations, Phys. Rev. Lett. 94, 170201 (2005).
  * Aoki  _et al._ [2014] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Nonequilibrium dynamical mean-field theory and its applications, Rev. Mod. Phys. 86, 779 (2014).
  * Hallberg [1995] K. A. Hallberg, Density-matrix algorithm for the calculation of dynamical properties of low-dimensional systems, Phys. Rev. B 52, R9827 (1995).
  * Kühner and White [1999] T. D. Kühner and S. R. White, Dynamical correlation functions using the density matrix renormalization group, Phys. Rev. B 60, 335 (1999).
  * Jeckelmann [2002] E. Jeckelmann, Dynamical density-matrix renormalization-group method, Phys. Rev. B 66, 045114 (2002).
  * Benthien  _et al._ [2004] H. Benthien, F. Gebhard, and E. Jeckelmann, Spectral function of the one-dimensional hubbard model away from half filling, Phys. Rev. Lett. 92, 256401 (2004).
  * White and Feiguin [2004] S. R. White and A. E. Feiguin, Real-time evolution using the density matrix renormalization group, Phys. Rev. Lett. 93, 076401 (2004).
  * Paeckel  _et al._ [2019] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Time-evolution methods for matrix-product states, Annals of Physics 411, 167998 (2019).
  * Zaletel  _et al._ [2015] M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore, and F. Pollmann, Time-evolving a matrix product state with long-ranged interactions, Phys. Rev. B 91, 165112 (2015).
  * Gohlke  _et al._ [2017] M. Gohlke, R. Verresen, R. Moessner, and F. Pollmann, Dynamics of the kitaev-heisenberg model, Phys. Rev. Lett. 119, 157203 (2017).
  * Verresen  _et al._ [2018] R. Verresen, F. Pollmann, and R. Moessner, Quantum dynamics of the square-lattice heisenberg model, Phys. Rev. B 98, 155102 (2018).
  * Verresen  _et al._ [2019] R. Verresen, R. Moessner, and F. Pollmann, Avoided quasiparticle decay from strong quantum interactions, Nature Physics 15, 750–753 (2019).
  * Van Damme and Vanderstraeten [2022] M. Van Damme and L. Vanderstraeten, Momentum-resolved time evolution with matrix product states, Phys. Rev. B 105, 205130 (2022).
  * Ferrari  _et al._ [2018] F. Ferrari, A. Parola, S. Sorella, and F. Becca, Dynamical structure factor of the ${J}_{1}-{J}_{2}$ heisenberg model in one dimension: The variational monte carlo approach, Phys. Rev. B 97, 235103 (2018).
  * Ferrari and Becca [2019] F. Ferrari and F. Becca, Dynamical structure factor of the ${J}_{1}-{J}_{2}$ heisenberg model on the triangular lattice: Magnons, spinons, and gauge fields, Phys. Rev. X 9, 031026 (2019).
  * Knap  _et al._ [2013] M. Knap, A. Kantian, T. Giamarchi, I. Bloch, M. D. Lukin, and E. Demler, Probing real-space and time-resolved correlation functions with many-body ramsey interferometry, Phys. Rev. Lett. 111, 147205 (2013).
  * Baez  _et al._ [2020] M. L. Baez, M. Goihl, J. Haferkamp, J. Bermejo-Vega, M. Gluza, and J. Eisert, Dynamical structure factors of dynamical quantum simulators, Proceedings of the National Academy of Sciences 117, 26123 (2020).
  * Sun  _et al._ [2023] J. Sun, L. Vilchez-Estevez, V. Vedral, A. T. Boothroyd, and M. S. Kim, Probing spectral features of quantum many-body systems with quantum simulators, arXiv e-prints , arXiv:2305.07649 (2023), arXiv:2305.07649 [quant-ph] .
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Nomura and Imada [2021] Y. Nomura and M. Imada, Dirac-type nodal spin liquid revealed by refined quantum many-body solver using neural-network wave function, correlation ratio, and level spectroscopy, Phys. Rev. X 11, 031034 (2021).
  * Nomura [2021] Y. Nomura, Helping restricted boltzmann machines with quantum-state representation by restoring symmetry, Journal of Physics: Condensed Matter 33, 174003 (2021).
  * Astrakhantsev  _et al._ [2021] N. Astrakhantsev, T. Westerhout, A. Tiwari, K. Choo, A. Chen, M. H. Fischer, G. Carleo, and T. Neupert, Broken-symmetry ground states of the heisenberg model on the pyrochlore lattice, Phys. Rev. X 11, 041021 (2021).
  * Roth  _et al._ [2022] C. Roth, A. Szabó, and A. MacDonald, High-accuracy variational monte carlo for frustrated magnets with deep neural networks, arXiv:2211.07749 (2022).
  * Reh  _et al._ [2023] M. Reh, M. Schmitt, and M. Gärttner, Optimizing design choices for neural quantum states (2023).
  * Chen and Heyl [2023] A. Chen and M. Heyl, Efficient optimization of deep neural quantum states toward machine precision (2023).
  * Schmitt and Heyl [2020] M. Schmitt and M. Heyl, Quantum many-body dynamics in two dimensions with artificial neural networks, Phys. Rev. Lett. 125, 100503 (2020).
  * 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).
  * Fabiani  _et al._ [2021] G. Fabiani, M. D. Bouman, and J. H. Mentink, Supermagnonic propagation in two-dimensional antiferromagnets, Phys. Rev. Lett. 127, 097202 (2021).
  * Schmitt  _et al._ [2022] M. Schmitt, M. M. Rams, J. Dziarmaga, M. Heyl, and W. H. Zurek, Quantum phase transition dynamics in the two-dimensional transverse-field ising model, Science Advances 8, 10.1126/sciadv.abl6850 (2022).
  * Donatella  _et al._ [2022] K. Donatella, Z. Denis, A. L. Boité, and C. Ciuti, Dynamics with autoregressive neural quantum states: application to critical quench dynamics, arXiv:2209.03241 (2022).
  * Mendes-Santos  _et al._ [2023a] T. Mendes-Santos, M. Schmitt, A. Angelone, A. Rodriguez, P. Scholl, H. J. Williams, D. Barredo, T. Lahaye, A. Browaeys, M. Heyl, and M. Dalmonte, Wave function network description and kolmogorov complexity of quantum many-body systems, arXiv:2301.13216 (2023a).
  * Choo  _et al._ [2018] K. Choo, G. Carleo, N. Regnault, and T. Neupert, Symmetries and many-body excitations with neural-network quantum states, Phys. Rev. Lett. 121, 167204 (2018).
  * Hendry and Feiguin [2019] D. Hendry and A. E. Feiguin, Machine learning approach to dynamical properties of quantum many-body systems, Phys. Rev. B 100, 245123 (2019).
  * Hendry  _et al._ [2021] D. Hendry, H. Chen, P. Weinberg, and A. E. Feiguin, Chebyshev expansion of spectral functions using restricted boltzmann machines, Phys. Rev. B 104, 205130 (2021).
  * Fabiani and Mentink [2019] G. Fabiani and J. H. Mentink, Investigating ultrafast quantum magnetism with machine learning, SciPost Phys. 7, 4 (2019).
  * Ebadi  _et al._ [2021] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum phases of matter on a 256-atom programmable quantum simulator, Nature 595, 227–232 (2021).
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Quantum simulation of 2d antiferromagnets with hundreds of rydberg atoms, Nature 595, 233–238 (2021).
  * Felser  _et al._ [2021] T. Felser, S. Notarnicola, and S. Montangero, Efficient tensor network ansatz for high-dimensional quantum many-body problems, Phys. Rev. Lett. 126, 170603 (2021).
  * Kalinowski  _et al._ [2022] M. Kalinowski, R. Samajdar, R. G. Melko, M. D. Lukin, S. Sachdev, and S. Choi, Bulk and boundary quantum phase transitions in a square rydberg atom array, Phys. Rev. B 105, 174417 (2022).
  * Cybenko [1989] G. Cybenko, Approximation by superpositions of a sigmoidal function, Math. Control. Signals, Syst. 2, 303 (1989).
  * Hornik [1991] K. Hornik, Approximation capabilities of multilayer feedforward networks, Neural Netw 4, 251 (1991).
  * Kim and Adalı [2003] T. Kim and T. Adalı, Approximation by fully complex multilayer perceptrons, Neural Comput. 15, 1641 (2003).
  * Roux and Bengio [2008] N. L. Roux and Y. Bengio, Representational power of restricted boltzmann machines and deep belief networks, Neural Comput. 20, 1631 (2008).
  * Sorella [1998] S. Sorella, Green function monte carlo with stochastic reconfiguration, Phys. Rev. Lett. 80, 4558 (1998).
  * Sorella [2001] S. Sorella, Generalized lanczos algorithm for variational quantum monte carlo, Phys. Rev. B 64, 024512 (2001).
  * Hackl  _et al._ [2020] L. Hackl, T. Guaita, T. Shi, J. Haegeman, E. Demler, and J. I. Cirac, Geometry of variational methods: dynamics of closed quantum systems, SciPost Phys. 9, 48 (2020).
  * Carleo [2011] G. Carleo, Spectral and dynamical properties of strongly correlated systems (2011).
  * [50] Supplemental material shows (i) further details for the implementation of the real-space and momentum-space scheme with tVMC, (ii) details about the NQS architecture used here, and (iii) convergence checks for the results with network sizes.
  * Wu  _et al._ [2020] Y. Wu, L.-M. Duan, and D.-L. Deng, Artificial neural network based computation for out-of-time-ordered correlators, Phys. Rev. B 101, 214308 (2020).
  * Sharir  _et al._ [2020] O. Sharir, Y. Levine, N. Wies, G. Carleo, and A. Shashua, Deep autoregressive models for the efficient variational simulation of many-body quantum systems, Phys. Rev. Lett. 124, 020503 (2020).
  * Pelissetto and Vicari [2002] A. Pelissetto and E. Vicari, Critical phenomena and renormalization-group theory, Physics Reports 368, 549–727 (2002).
  * Hamer  _et al._ [2006a] C. J. Hamer, J. Oitmaa, and W. Zheng, One-particle dispersion and spectral weights in the transverse ising model, Phys. Rev. B 74, 174428 (2006a).
  * Hamer  _et al._ [2006b] C. J. Hamer, J. Oitmaa, Z. Weihong, and R. H. McKenzie, Critical behavior of one-particle spectral weights in the transverse ising model, Phys. Rev. B 74, 060402 (2006b).
  * Note [1] Here, we consider a finite cutoff for the interaction potential such that the interactions are set to $0$ for $r_{i,j}>r_{0}=3a$.
  * Samajdar  _et al._ [2020] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Complex density wave orders and quantum phase transitions in a model of square-lattice rydberg atom arrays, Phys. Rev. Lett. 124, 103601 (2020).
  * Campostrini  _et al._ [2014] M. Campostrini, J. Nespolo, A. Pelissetto, and E. Vicari, Finite-size scaling at first-order quantum transitions, Phys. Rev. Lett. 113, 070402 (2014).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021).
  * Verresen  _et al._ [2021] R. Verresen, M. D. Lukin, and A. Vishwanath, Prediction of toric code topological order from rydberg blockade, Phys. Rev. X 11, 031005 (2021).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Annals of Physics 321, 2 (2006), january Special Issue.
  * Becker and Wessel [2018] J. Becker and S. Wessel, Diagnosing fractionalization from the spin dynamics of ${Z}_{2}$ spin liquids on the kagome lattice by quantum monte carlo simulations, Phys. Rev. Lett. 121, 077202 (2018).
  * Sun  _et al._ [2018] G.-Y. Sun, Y.-C. Wang, C. Fang, Y. Qi, M. Cheng, and Z. Y. Meng, Dynamical signature of symmetry fractionalization in frustrated magnets, Phys. Rev. Lett. 121, 077201 (2018).
  * Jülich Supercomputing Centre [2019] Jülich Supercomputing Centre, JUWELS: Modular Tier-0/1 Supercomputer at the Jülich Supercomputing Centre, J. Large-Scale Res. Facilities 5 (2019).
  * Schmitt and Reh [2022] M. Schmitt and M. Reh, jVMC: Versatile and performant variational Monte Carlo leveraging automated differentiation and GPU acceleration, SciPost Phys. Codebases , 2 (2022).
  * Bradbury  _et al._ [2018] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs (2018).
  * Weinberg and Bukov [2017] P. Weinberg and M. Bukov, QuSpin: a Python package for dynamics and exact diagonalisation of quantum many body systems part I: spin chains, SciPost Phys. 2, 003 (2017).
  * Weinberg and Bukov [2019] P. Weinberg and M. Bukov, QuSpin: a Python package for dynamics and exact diagonalisation of quantum many body systems. Part II: bosons, fermions and higher spins, SciPost Phys. 7, 020 (2019).
  * Mendes-Santos  _et al._ [2023b] T. Mendes-Santos, M. Schmitt, and M. Heyl, Highly resolved spectral functions of two- dimensional systems with neural quantum states: data (2023b).

##  I Further details on the simulation of dynamical correlations with tVMC

In this section, we discuss further important details to simulate the
dynamical correlations $S^{zz}(\boldsymbol{q},t)$ and
$G^{zz}_{\boldsymbol{r},\boldsymbol{r}^{\prime}}(t)$.

Momentum-space scheme. As described in the main text, we obtain
$S^{zz}(\boldsymbol{q},t)$ by (i) performing a two-sided time evolution
followed by (ii) the evaluation of the overlap of two time-evolved states; see
Eq. (8).

First, to time evolve the states $\bra{\phi^{z}_{\boldsymbol{q}}(-\delta t)}$
and $\ket{\phi^{z}_{\boldsymbol{q}}(\delta t)}$, we solve the TDVP equation
[Eq. (4)] with a second-order integration method. Specifically, we implement
the adaptative-time-step scheme to describe the time evolution of the ket-
state $\ket{\phi^{z}_{\boldsymbol{q}}(\delta t)}$. Such an approach allows us
to estimate errors based on varying time step sizes and adjust $\delta t$
during the simulation. The adaptative-time-step scheme is implemented with the
Heun method; we refer the reader to ref. [28] for further details. The time
evolution of the bra-state $\bra{\phi^{z}_{\boldsymbol{q}}(-\delta t)}$ is
then performed by using the Heun method with the time step $\delta t$ defined
by the ket-state time evolution.

Second, to compute the overlap between the two unnormalized states
$\bra{\phi^{z}_{\boldsymbol{q}}(-t/2)}$ and
$\ket{\phi^{z}_{\boldsymbol{q}}(t/2)}$ we employ Monte Carlo sampling to
obtain the normalized overlap (to simplify the notation, let us define
$\bra{\phi^{z}_{\boldsymbol{q}}(-t/2)}=\bra{\phi_{1}}$ and
$\ket{\phi^{z}_{\boldsymbol{q}}(t/2)}=\ket{\phi_{2}}$)

$O=\frac{\braket{\phi_{1}}{\phi_{2}}}{\sqrt{\braket{\phi_{1}}{\phi_{1}}\braket{\phi_{2}}{\phi_{2}}}}.$
(13)

The idea is that the quantity
$p_{1}=\braket{\phi_{1}}{\phi_{2}}/\braket{\phi_{1}}{\phi_{1}}$, can be
estimated via Monte Carlo sampling, i.e.,

$p_{1}=\sum_{\boldsymbol{s}}\frac{|\phi_{1}(\boldsymbol{s})|^{2}}{\braket{\phi_{1}}{\phi_{1}}}\frac{\phi_{2}(\boldsymbol{s})}{\phi_{1}(\boldsymbol{s})}\approx\frac{1}{N_{MC}}\sum_{i}\frac{\phi_{2}(\boldsymbol{s}_{i})}{\phi_{1}(\boldsymbol{s}_{i})},$
(14)

where $\boldsymbol{s}_{1},...,\boldsymbol{s}_{N_{MC}}$ represent a set of
$N_{MC}$ MC samples generated by the Born distribution
$|\phi_{1}(\boldsymbol{s})|^{2}$. Analogously,
$p_{2}=\braket{\phi_{2}}{\phi_{1}}/\braket{\phi_{2}}{\phi_{2}}$ is estimated
via MC sampling of $|\phi_{2}(\boldsymbol{s})|^{2}$. By combining the two
results, we obtain $O=\sqrt{p_{1}p_{2}^{*}}$ [51].

Finally, by computing
$O_{\boldsymbol{q}}=\braket{\phi^{z}_{\boldsymbol{q}}(-t/2)}{\phi^{z}_{\boldsymbol{q}}(t/2)}/N_{\boldsymbol{q}}$,
where the normalization constant

$N_{\boldsymbol{q}}=\sqrt{\braket{\phi^{z}_{\boldsymbol{q}}(-t/2)}{\phi^{z}_{\boldsymbol{q}}(-t/2)}\braket{\phi^{z}_{\boldsymbol{q}}(t/2)}{\phi^{z}_{\boldsymbol{q}}(t/2)}}$
(15)

is simply related to the ground-state static structure fractor,
$N_{\boldsymbol{q}}=N_{s}\sum_{\boldsymbol{r}}e^{i\boldsymbol{r}\boldsymbol{q}}\left<S^{z}_{\boldsymbol{r}}S^{z}_{\boldsymbol{0}}\right>$,
we obtain the dynamical correlation in momentum space

$S^{zz}(\boldsymbol{q},t)=e^{iE_{0}t}N_{\boldsymbol{q}}O_{\boldsymbol{q}}.$
(16)

Real-space scheme. We use Eq. (10) to compute the longitudinal component of
the Green function $G^{zz}_{\boldsymbol{r},\boldsymbol{r}^{\prime}}(t)$, for a
deduction of such expression see refs. [18, 19]. The expectation value
$\bra{\phi_{\boldsymbol{r}^{\prime}}^{z}(t)}\hat{S}_{\boldsymbol{r}}^{z}\ket{\phi_{\boldsymbol{r}^{\prime}}^{z}(t)}$
is obtained by simulating the time evolution of the local excitation
$\ket{\phi_{\boldsymbol{r}^{\prime}}^{z}(t)}=e^{-i\hat{H}t}e^{i\frac{\pi}{4}\hat{S}_{\boldsymbol{r^{\prime}}}^{z}}\ket{\psi_{0}}$
with the tVMC, where the local excitation is encoded on the NQS ansatz by
adding the configuration-dependent factor
$f(\boldsymbol{r^{\prime}},\boldsymbol{s})=1/\sqrt{2}\left(1-is_{\boldsymbol{r^{\prime}}}\right)$
on the initial-time quantum state

$\ket{\phi_{\boldsymbol{r}^{\prime}}^{z}(0)}=\sum_{\boldsymbol{s}}f(\boldsymbol{r^{\prime}},\boldsymbol{s})\psi_{\boldsymbol{\theta}}(\boldsymbol{s})\ket{\boldsymbol{s}}.$
(17)

Similarly, the term
$\bra{\psi_{0}}\hat{S}_{\boldsymbol{r}}^{z}\hat{S}_{\boldsymbol{r}^{\prime}}^{z}(t)\hat{S}_{\boldsymbol{r}}^{z}\ket{\psi_{0}}$
is obtained by simulating the time evolution of the local excitation
$\hat{S}_{\boldsymbol{r}}^{z}\ket{\psi_{0}}$, which is encoded on the NQS
ansatz by adding the factor
$f^{\prime}(\boldsymbol{r},\boldsymbol{s})=s_{\boldsymbol{r}}$.

Fourier transform Let us now discuss the details of the Fourier transform (FT)
used to obtain the dynamical structure factor in the frequency-momentum
domain: $S^{zz}(\boldsymbol{q},\omega)$.

In the momentum-space scheme, we perform a FT of $S^{zz}(\boldsymbol{q},t)$ to
the frequency domain

$S^{zz}(\boldsymbol{q},\omega)=\frac{1}{2\pi}\int_{-\infty}^{\infty}dte^{i\omega
t}e^{-\gamma^{2}t^{2}}S^{zz}(\boldsymbol{q},t),$ (18)

where we apply a Gaussian envelope to avoid the finite-time effects of our
simulation; the parameter $\gamma=2/t_{max}$ controls the broadening factor of
the Gaussian peaks, with $t_{max}$ being the maximum times reliably obtained
during our simulations. For the results shown in the main text, we typically
performed simulations with $t_{max}\approx 20J$ and $t_{max}\approx 40\Omega$
for the 2D QIM and RyM, respectively. Furthermore, we employ the trapezoidal
rule to perform the numerical integration of Eq. (18).

In the real-space scheme, we calculate the Green function,
$G^{z}_{\boldsymbol{r},\boldsymbol{r}^{\prime}}(t)$, for all pairs of spins
($\boldsymbol{r},\boldsymbol{r}^{\prime}$), and then perform a FT in real
space and time to obtain $G^{z}(\boldsymbol{q},\omega)$ (the FT to the
frequency domain is performed as described above). Finally, we relate the
retarded Green function to the DSF via the fluctuation-dissipation theorem

$S^{zz}(\boldsymbol{q},\omega)=-\frac{1}{\pi}Im[G^{z}(\boldsymbol{q},\omega)].$
(19)

##  II Network architectures and hyperparameters

As mentioned in the main text, the real-space-scheme is implemented with
restricted Boltzmann machine (RBM) NQSs [21], which are architectures composed
by a single fully-connected (dense) hidden layer. Concretely, RBMs are defined
as

$\psi_{\boldsymbol{\theta}}(\boldsymbol{s})=\prod_{i=1}^{M}\cosh\left(\sum_{j=1}^{N}W_{i,j}s_{j}+b_{i}\right),$
(20)

where the set of variational parameters
$\boldsymbol{\theta}=(W,\boldsymbol{b})$ are the complex weight matrix, $W$,
and the hidden bias, $\boldsymbol{b}$. The hyperparameter $M$ represents the
number of nodes in the single hidden layer, which sets the total number of
variational parameters, i.e., $N_{p}=2(NM+M)$.

For the results shown in Fig. 1 (b1) and (c1) we consider $M=40$.

The momentum-space scheme is implemented with the convolutional neural network
ansatz [28]. Specifically, we consider architectures composed of a set of
$L_{l}$ stacked convolutional layers. Each output neuron of a convolutional
layer is described by

$a^{(l)}_{c,k}(s)=f_{l}\left(\sum_{c^{\prime}=1}^{\alpha_{(l-1)}}\sum_{t=1}^{F^{(l)}}\tilde{W}^{(c,l)}_{k,c^{\prime},t}a^{(l-1)}_{c^{\prime},T_{k}(t)}+b_{c}^{(l)}\right),$
(21)

where $a^{(l)}_{c,k}$ represent a value at a channel $c$ of an output layer
$l$; the index $c=1,...,\alpha_{l}$ runs over the total number of output
channels. The value of a neuron, $a^{(l)}_{c,k}$, is defined by an affine map
that resembles a convolution of the previous layer with a filter along an
orbit generated by translations. Specifically, the individual neurons
$a^{(l)}_{c,k}$ are coupled with identical weights $\tilde{W}$ to the neurons
of the previous layer, which are transformed by translations, $T_{k}(t)$, by
some fixed stride. In particular, we consider a $2D$ filter with size $F\times
F$. Such map can also include a bias $b_{m}^{(l)}$ for each channel. Finally,
$f_{l}(...)$ represents an activation function.

Starting with the initial layer, which is the spin configuration
$a^{l=0}_{1,j}(s)=s_{j}$, and considering the recursive relation for the
$L_{l}$ convolutional layers, we can define the logarithm of the wave function
in terms of the last convolutional layer

$\ln\psi_{\boldsymbol{\theta}}=\sum_{c,k}a^{(L_{l})}_{c,k}(s).$ (22)

The total number of variational parameters a CNN NQS is given by

$N_{p}^{CNN}=2\left(F^{2}\sum_{l=1}^{L_{l}}\alpha_{l-1}\alpha_{l}+N_{b}\right),$
(23)

where $\alpha_{0}=1$ and $N_{b}$ is the number of biases. In this work, we
consider CNN NQSs with the following structure
$\boldsymbol{\alpha}=(\alpha_{1},...,\alpha_{L_{l}};F)$, where $\alpha_{k}$ is
the number of channels in each layer $k$. Furthermore, we use a stride
$\boldsymbol{s}=[1,1]$. For the results present in the main text we adopt the
following hyperparameters:

Figure 4: Additional results for the 2D QIM. In panel (a) we show a comparison
of the imaginary part of dynamical correlation $S^{zz}_{im}(\boldsymbol{q},t)$
with exat results for $g=3.1$, while in panel (b) we check the convergence of
the results for $L=12$ by increasing the linear size of CNN filter $F$; in
this case we consiter $\boldsymbol{\alpha}=(4,3;F)$ and $g=3.044$.

Figure 5: Additional results for the 2D RyM. The upper pannel show
$S^{zz}_{im}(\boldsymbol{q}=M,t)$ for (a) $L=4$, (b) $L=8$ and (c) $L=12$ and
$R_{b}=1.45$ and $\delta=-2.00$. In the lower panels we show the correspond
Fourier transform to frequency $S^{zz}_{im}(\boldsymbol{q}=M,\omega)$.

  * •

Results for the 2D QIM (Figs. 1 and 2). For $L<20$, we consider
$\boldsymbol{\alpha}=(4,3;F=L/2)$, while for $L=20$,
$\boldsymbol{\alpha}=(4,3;F=8)$, and $L=24$, $\boldsymbol{\alpha}=(4,3;F=10)$.
We employ a sixth-degree polynomial in the first layer and a fifth-degree
polynomial in subsequent layers.

  * •

Results for the 2D RyM (Fig. 3). We consider: $\boldsymbol{\alpha}=(6,5;F=4)$
for $L=12$ and $\boldsymbol{\alpha}=(4,3;F=4)$ for $L=8$ and $4$ results. We
use fifth-degree polynomial activation functions in all the CNN layers.

We perform all the simulations considering $N_{MC}=4\times 10^{4}$ Monte
Carlos samples.

##  III Convergence checks and bechmarks

We now provide additional results to check the convergence of the results
presented in Figs. 2 and 3 of the main text. Furthermore, we perform further
benchmarks with exact results and theretical descriptions of QCPs.

For the smaller system size considered here, $L=4$, we make a direct
comparison with exact results. As can be seen from Fig. 4 (a),
$S^{zz}_{im}(\boldsymbol{q}=\Gamma,t)$ exhibit an almost perfect agreement
with ED results. To be more quantitative in the comparison, we present the
values of the lowest energy gap, $\Delta(\boldsymbol{q}=\Gamma)$ computed with
NQS and ED on table 1. In the exact case, we directly compute the gap from
$\Delta(\boldsymbol{q}=\Gamma)=E_{1}-E_{0}$, where $E_{0}$ and $E_{1}$ are the
ground-state and the first-excited-state energies, respectively. Using the NQS
approach, we extract the gap from the DSF, i.e.,
$\Delta(\boldsymbol{q}=\Gamma)=\text{argmax}_{\omega}S^{zz}(\boldsymbol{q}=\Gamma,\omega)$.
The discrepancy between such results are
$\epsilon=|\Delta-\Delta^{EX}|/\Delta^{EX}<0.2\%$ for all the values $g$
considered in Fig. 2, which is within the statistical error bars for
$N_{MC}=4\times 10^{4}$ samples used here.

For $L>4$, a direct comparison with exact results is not possible.
Nevertheless, an important aspect of the neural quantum state approach is that
we can test the accuracy of the results by comparing NQS with different sizes
or architectures. For this, we provide additional results for
$S^{zz}_{im}(\boldsymbol{q}=\Gamma,t)$ and different
$\boldsymbol{\alpha}=(4,3;F)$. By focusing on the challenging regime
$g=3.044$, we observe the convergence of the results for the larger values of
$F$ considered; see Fig. 4 (b).

Spectral gap (QIM) - L = 4
---
– |  Exact |  NQS |  $\epsilon$
$g=3.5$ |  $2.1124$ |  $2.1138$ |  $0.0007$
$g=3.1$ |  $1.2159$ |  $1.2178$ |  $0.0016$
$g=3.044$ |  $1.0968$ |  $1.0970$ |  $0.0001$
Table 1: Spectral gap of the 2D QIM.

We now consider the convergence of our results for the Rydberg-atoms-array
model (RyM). In particular, we focus on a regime of parameters in which the
system is near the disordered-StDW transition (i.e., $R_{b}=1.45$ and
$\delta=-2.00$). In this case, we observe that
$S^{zz}_{im}(\boldsymbol{q}=M,\omega)$ exhibit a non-trivial size scaling, as
discussed in Fig. 3 (c2) of the main text the spectral weight associated with
the lower-energy peak of $S^{zz}(\boldsymbol{q}=\boldsymbol{M},\omega)$
decreases with system size, while we observe a spectral weight transfer to
higher energies for larger values of $L$; see Fig. 3 (c2.1). As can be seen in
Fig. 5, results for $S^{zz}_{im}(\boldsymbol{q}=M,t)$ for different
$\boldsymbol{\alpha}$ converge for $t\lesssim 30\Omega$. By chosing
$t_{max}=30\Omega$, we show in the lower panels of 5 the results in frequency,
confirming the non-trivial scaling of $S^{zz}_{im}(\boldsymbol{q}=M,\omega)$.

Figure 6: Data for the finite-size dependence of the spectral gap $\Delta$ as
obtained through the NQS simulations close to both the QIM and the disordered-
CbDW quantum critical points. We further include the fitted values of the
critical exponent $z$ obtained by assuming that $\Delta\sim L^{-z}$ and the
standard errors of the coefficient of the linear fitting.

Finally, as a more demanding benchmark of our results, we perform a
quantitative verification of the size-scaling of the low-energy spectral gap,
$\Delta\sim L^{-z}$, at the QCPs considered in this work. Fig. 6 displays our
data for $\Delta$ as a function of system size. Further, we also show fitted
values of $z$ for both the QIM (for $g=3.044$) and the disordered-CbDW (for
$R_{b}=1.2$ and $\delta=1.05$) transitions obtained by assuming a functional
dependence $\Delta\sim L^{-z}$ of the gap. Our estimation of $z$ is consistent
with the theoretical expectations for the 2D QIM, $z=1$. Going beyond a
benchmark, we also estimate the dynamical exponent for the disordered-CbDW [in
this case, we consider $\Delta(\boldsymbol{q}=\boldsymbol{M})$] the obtained
value is consistent with the Ising universality class [57, 41].

##  IV Ground state results for the 2D RyM

As a complement to the results of the 2D RyM, we discuss the onset of the CbDW
and StDW phases by considering the ground-state (GS) correlations

$C^{z}(\boldsymbol{r})=\left<\sigma^{z}_{\boldsymbol{0}}\sigma^{z}_{\boldsymbol{r}}\right>-\left<\sigma^{z}_{\boldsymbol{0}}\right>\left<\sigma^{z}_{\boldsymbol{r}}\right>.$
(24)

Particularly, we show in Figs. 7 (a) and (b) the behavior of the static
structure factor,
$S_{\boldsymbol{q}}=1/N_{s}\sum_{\boldsymbol{r}}e^{i\boldsymbol{r}\boldsymbol{q}}C^{z}(\boldsymbol{r})$,
along the lines $R_{b}=1.2$ and $R_{b}=1.45$, respectively. For $R_{b}=1.2$,
we note an increase of the $S_{\boldsymbol{q=M}}$ for $\delta\approx 1.05$,
indicating the emergence of a phase with CbDW pattern of correlations, while
for $R_{b}=1.45$ we note an increase of both $S_{\boldsymbol{q=M}}$ as
$S_{\boldsymbol{q=X}}$ for $\delta\approx 2.1$, indicating the emergence of
the StDW phase. Those results are in agreement with previous simulations of GS
properties the 2D RyM [57, 41].

Figure 7: Ground-state results for the 2D RyM. Static structure factor
$S_{\boldsymbol{q}}$ as function of $\delta$ for (a) $R_{b}=1.2$ and (b)
$R_{b}=1.45$, and $L=8$.