Lee-Yang theory of quantum phase transitions with neural network quantum
  states
Predicting the phase diagram of interacting quantum many-body systems is a
central problem in condensed matter physics and related fields. A variety of
quantum many-body systems, ranging from unconventional superconductors to spin
liquids, exhibit complex competing phases whose theoretical description has
been the focus of intense efforts. Here, we show that neural network quantum
states can be combined with a Lee-Yang theory of quantum phase transitions to
predict the critical points of strongly-correlated spin lattices. Specifically,
we implement our approach for quantum phase transitions in the transverse-field
Ising model on different lattice geometries in one, two, and three dimensions.
We show that the Lee-Yang theory combined with neural network quantum states
yields predictions of the critical field, which are consistent with large-scale
quantum many-body methods. As such, our results provide a starting point for
determining the phase diagram of more complex quantum many-body systems,
including frustrated Heisenberg and Hubbard models.
Pascal M. Vecsei, Christian Flindt, Jose L. Lado
# Lee-Yang theory of quantum phase transitions with neural network quantum states

###### Abstract

Predicting the phase diagram of interacting quantum many-body systems is a central problem in condensed matter physics and related fields. A variety of quantum many-body systems, ranging from unconventional superconductors to spin liquids, exhibit complex competing phases whose theoretical description has been the focus of intense efforts. Here, we show that neural network quantum states can be combined with a Lee-Yang theory of quantum phase transitions to predict the critical points of strongly-correlated spin lattices. Specifically, we implement our approach for quantum phase transitions in the transverse-field Ising model on different lattice geometries in one, two, and three dimensions. We show that the Lee-Yang theory combined with neural network quantum states yields predictions of the critical field, which are consistent with large-scale quantum many-body methods. As such, our results provide a starting point for determining the phase diagram of more complex quantum many-body systems, including frustrated Heisenberg and Hubbard models.

## I Introduction

Solving a generic family of quantum many-body problems and ultimately predicting their phase diagram is a challenging task [1; 2]. The exponential growth of the Hilbert space with the system size, especially for high dimensional systems, makes most realistic models intractable in practice. Some problems, such as the transverse-field Ising model in one dimension, can be solved analytically [3]. However, more generally, obtaining the phase diagram of an interacting quantum many-body system is a critical open problem. To this end, several numerical tools have been developed, including Monte Carlo simulations [4], and tensor-network algorithms [5]. Nevertheless, despite considerable progress, the phase diagram of many quantum systems in two and three dimensions remain unknown [6; 7].

Neural network quantum states are a recently developed class of variational states [8] that have shown great potential for parametrizing and finding the ground state of interacting quantum many-body systems [9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26]. Neural network quantum states represent the wave function of a quantum many-body system as a neural network. Specifically, the neural network is a parametrized function that takes the configuration of a many-body system as the input and outputs the corresponding amplitude and phase of the wave function. By optimizing the parameters of the neural network, so that the energy is minimized, an accurate approximation of the ground state can be found. Neural network quantum states exploit the fact that neural networks can faithfully represent many complex functions [27], including a variety of quantum many-body wave functions. They have already been applied to find the wave functions of several spin models [9; 10; 11; 12; 13; 14; 28], including the \(J_{1}-J_{2}\) Heisenberg model [15; 16; 17; 18; 19; 20; 21]. Moreover, their use has been extended to fermionic [22; 29] and bosonic [30; 31; 32] systems, as well as to molecules [23; 22] and nuclei [24; 25; 26].

In the context of critical behavior, a rigorous foundation of phase transitions was established by Lee and Yang, who considered the zeros of the partition function in the complex plane of the control parameters, for example an external magnetic field or the inverse temperature [33; 34; 35; 36]. This approach relies on the fact that for systems of finite size, the partition function zeros are all complex. However, if a system exhibits a phase transition, the zeros will approach the critical value on the real

Figure 1: Neural network approach to quantum phase transitions. (a) Cubic Ising lattice of interacting spins in a transverse magnetic field, here a system of size \(3\times 3\times 3\). (b) A neural network takes a configuration of the spins, encoded in the vector \(\vec{\sigma}=(\sigma_{1},...,\sigma_{N})\), and outputs the corresponding value of the wave function, \(\psi_{\vec{\delta}}(\vec{\sigma})=(\vec{\sigma}|\psi)\), which depends on the variational parameters in \(\vec{\theta}\). (c) From the fluctuations of the magnetization, we extract the zeros of the moment generating function of the magnetization and investigate their motion in the complex plane as we increase the system size. (d) Above the critical field, \(h>h_{c}\), the zeros remain complex in the thermodynamic limit, and the system is in the paramagnetic phase (PM). At \(h=h_{c}\), the zeros reach the real-axis, signaling a quantum phase transition. For \(h<h_{c}\), the system is in the ferromagnetic phase (FM) with finite magnetization.
axis in the thermodynamic limit of large system sizes, giving rise to a non-analytic behavior of the free energy density [37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48]. Lee-Yang zeros are not just a theoretical concept, but they can also be determined experimentally [49; 50; 51; 52; 53]. In recent years, applications of Lee-Yang theory have been expanded to dynamical quantum phase transitions in quantum many-body systems after a quench [54; 55; 56] and to quantum phase transitions in systems at zero temperature [57; 58].

Here, we combine neural network quantum states with a Lee-Yang theory of quantum phase transitions to predict the critical behavior of interacting spin lattices in one, two, and three dimensions. As illustrated in Fig. 1(a), we consider the transverse-field Ising model in different dimensions and lattice geometries. We then find the ground state of the system as well as the fluctuations of the magnetization using neural network quantum states, Fig. 1(b). From these fluctuations, we determine the complex zeros of the moment generating function of the magnetization and follow their motion as the system size is increased. As illustrated in Fig. 1(c), the zeros remain complex in the thermodynamic limit in case there is no phase transition. On the other hand, if the magnetic field is tuned to its critical value, the zeros of the moment generating function will reach the real axis, signaling a phase transition. Thus, by investigating the positions of the zeros for different magnetic fields, we can map out the phase diagram of the system, Fig. 1(d).

Our manuscript is organized as follows: In Sec. II, we describe the methods that we use throughout this work. In particular, we introduce the transverse-field Ising model, we discuss our calculations of the magnetization cumulants in the ground state using neural network quantum states, and we provide the details of the Lee-Yang theory that we use to predict the critical magnetic field for a given lattice geometry. In Sec. III, we present the results of our calculations. As examples, we first discuss our procedure for the transverse-field Ising model on a one-dimensional chain, a two-dimensional square lattice, and a cubic lattice in three dimensions. We then provide predictions of the critical fields for several other lattice geometries. In Sec. IV, we discuss our results and the role of the coordination number and dimensionality of a given lattice. We also compare our predictions with mean-field theory, which becomes increasingly accurate in higher dimensions. Finally, in Sec. V, we summarize our conclusions. Technical details of our neural network calculations are provided in Appendix A.

## II Methods

### Transverse-field Ising model

We consider the transverse-field Ising model on a lattice of spin-\(1/2\) sites as described by the Hamiltonian

\[\hat{\mathcal{H}}=-J\sum_{\{i,j\}}\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}-h \sum_{i}\hat{\sigma}_{i}^{x}. \tag{1}\]

Here, the first sum runs over all nearest neighbors, denoted by \(\{i,j\}\), the coupling between them is \(J\), and \(h\) is the transverse magnetic field. The one-dimensional version of this model can be solved analytically and it is known to exhibit a continuous phase transition at the critical field \(h_{c}=J\)[3]. Above the critical field, the system is in a paramagnetic phase with vanishing magnetization. Below it, the system exhibits spontaneous symmetry-breaking and enters a ferromagnetic phase with a non-vanishing magnetization. In the following we will investigate the model in different dimensions and geometries. The two-dimensional systems we consider are square, honeycomb, Kagome, and triangular lattices. In three dimensions, we consider cubic, face-centred cubic, body-centred cubic, and diamond lattices. In all of these cases, we impose periodic boundary conditions, and we compare our predictions with earlier results based on large-scale quantum Monte Carlo simulations [59].

### Neural network quantum states

To find the ground state of the system together with the moments and cumulants of the magnetization, we use neural network quantum states. The neural network quantum states are variational states of the form

\[\psi_{\vec{\theta}}(\vec{\sigma})=\langle\vec{\sigma}|\psi_{\vec{\theta}}\rangle, \tag{2}\]

where the vector \(\vec{\theta}\) contains the variational parameters that we need to determine to minimize the energy and thereby find the ground state. The neural network provides a compressed algorithmic representation of the coefficients of the wavefunction, and it takes a spin configuration in the computational basis as the input, and outputs the wave function in response. The energy is minimized using stochastic reconfiguration, which is an approximate imaginary time-evolution within the variational space of the neural network. Neural network state methodologies have been extended to the time-evolution of quantum systems [60; 8; 61], quantum state tomography [62; 63; 64], as well as finite-temperature equilibrium physics [65; 66; 67]. Importantly, while many other approaches are not able to exploit the computational power of massive parallel computing, neural network quantum states can be implemented with modern graphics processing units.

The energy is evaluated by sampling over the wave function as

\[\langle\hat{\mathcal{H}}\rangle=\frac{\sum_{\vec{\sigma}\vec{\sigma}^{\prime}} \psi^{*}(\vec{\sigma})\langle\vec{\sigma}|\hat{\mathcal{H}}|\vec{\sigma}^{ \prime}\rangle\psi(\vec{\sigma}^{\prime})}{\sum_{\vec{\sigma}^{\prime}}|\psi( \vec{\sigma}^{\prime})|^{2}}=\sum_{\vec{\sigma}}P_{\psi}(\vec{\sigma}) \mathcal{H}_{\text{loc}}(\vec{\sigma}), \tag{3}\]

where we have defined the probability

\[P_{\psi}(\vec{\sigma})=\frac{|\psi(\vec{\sigma})|^{2}}{\sum_{\vec{\sigma}^{ \prime}}|\psi(\vec{\sigma}^{\prime})|^{2}} \tag{4}\]
and the local spin Hamiltonian

\[\mathcal{H}_{\rm loc}(\vec{\sigma})=\sum_{\vec{\sigma}}\langle\hat{\sigma}|\hat{ \mathcal{H}}|\hat{\sigma}^{\prime}\rangle\frac{\psi(\vec{\sigma}^{\prime})}{ \psi(\vec{\sigma})}. \tag{5}\]

Since Eq. (3) is just an average with respect to a normalized probability distribution, Markov-chain Monte Carlo can be used for evaluating the energy and the gradients [68]. It is worth noting that the spin Hamiltonian in Eq. (5) is given by only a few terms in the sum, since only nearest neighbors are coupled. We will also need the expectation value of the total magnetization and its moments, which we express as

\[(\hat{M}_{z}^{n})=\sum_{\vec{\sigma}}P_{\psi}(\vec{\sigma})M_{z}^{n}(\vec{ \sigma}), \tag{6}\]

since \(\hat{M}_{z}\) is diagonal in the computational basis, such that \(M_{z}^{n}(\vec{\sigma})=(\langle\vec{\sigma}|\hat{M}_{z}|\vec{\sigma}\rangle)^ {n}=\langle\vec{\sigma}|\hat{M}_{z}^{n}|\vec{\sigma}\rangle\). Additional details of these calculations are provided in Appendix A.

### Lee-Yang theory

The classical Lee-Yang theory of phase transitions considers the zeros of the partition function in the complex plane of the control parameter, for instance magnetic field or inverse temperature [33, 34, 35, 36]. For finite systems, the partition function zeros are situated away from the real axis. However, in case of a phase transition, they will approach the critical value on the real axis in the thermodynamic limit. One may thereby predict the occurrence of a phase transition by investigating the position of the zeros as the system size is increased. The Lee-Yang theory of phase transitions has found applications in condensed matter physics [37, 40, 41, 44, 45, 46], atomic physics [38] and particle physics [47, 48, 49, 42, 43, 69, 70, 71, 72]. Recently, it has been extended to the zeros of the moment generating function that describes the fluctuations of the order parameter [57, 58] and thereby allows for the detection of quantum phase transitions. Following this approach, we define the moment generating function

\[\chi(s)=\langle e^{s\hat{M}_{z}}\rangle=\frac{1}{g}\sum_{k=1}^{g}\langle\psi_{ k}^{(0)}|e^{s\hat{M}_{z}}|\psi_{k}^{(0)}\rangle, \tag{7}\]

where \(\hat{M}_{z}\) is the total magnetization, and \(s\) is referred to as the counting field. Here, we have included the possibility that the system may have \(g\) degenerate and normalized ground states that we denote by \(|\psi_{k}^{(0)}\rangle\), \(k=1,\ldots,g\). Within this framework, the moment generating function plays the role of the partition function in the classical Lee-Yang theory, and the cumulant generating function, \(\Theta(s)=\ln\chi(s)\), becomes the corresponding free energy. The moments and cumulants of the magnetization are given by derivatives with respect to the counting field as

\[\langle\hat{M}_{z}^{n}\rangle=\partial_{s}^{n}\chi(s)|_{s=0} \tag{8}\]

and

\[\langle\!\langle\hat{M}_{z}^{n}\rangle\!\rangle=\partial_{s}^{n}\Theta(s)|_{s =0}. \tag{9}\]

Importantly, away from a phase transition, the cumulants are expected to grow linearly with the system size, such that the normalized cumulants \(\langle\!\langle\hat{M}_{z}^{n}\rangle\!\rangle/N\) converge to finite values as the number of spins \(N\) approaches infinity. By contrast, at a phase transition, a different scaling behavior is expected due to as non-analytic behavior of the cumulant generating function at \(s=0\)[40, 73]. This non-analytic behavior emerges in the thermodynamic limit,

Figure 2: Extraction of zeros from the cumulants of the magnetization. (a) Extracted zeros for a linear Ising chain in different magnetic fields, \(h=0.6,0.7,0.8,0.9,0.95,1.0,1.05,1.1,1.15,1.2J\) (starting from the lower curve), as a function of the inverse system size, \(1/L\). The solid lines are the finite-size scaling ansatz in Eq. (11), which allows us to determine the value in the thermodynamic limit, where \(1/L\) approaches zero. (b) Similar results for a two-dimensional square lattice with the following values of the magnetic field, \(h=0.5,1.0,1.5,2.0,2.5,2.9,3.0,3.1,3.2,3.3,3.4,3.5J\) (starting from the lower curve). (c) Results for a cubic lattice in three dimensions with \(h=0.0,1.0,2.0,4.0,5.0,5.2,5.4,5.6,5.8,6.0J\) (starting from the lower curve).
if the complex zeros of the moment generating function approach \(s=0\).

To determine the position of the zeros that are closest to \(s=0\), we use the cumulant method that was developed in Refs. [52; 53; 57; 58; 40]. In this approach, the zeros of the moment generating function can be determined from the high cumulants of the order parameter. By doing so for different system sizes, we can then find the convergence points in the thermodynamic limit using finite-size scaling [57; 58; 40; 52]. The cumulant method allows us to express the zeros in terms of the high cumulants of the magnetization. Moreover, for the transverse-field Ising model, the symmetry, \(\hat{U}^{\dagger}\hat{H}\hat{U}=\hat{H}\), with respect to the unitary operator \(\hat{U}=\prod_{i}\hat{\sigma}_{i}^{x}\) that flips all spins, implies that all odd cumulants vanish, and in this model the zeros are purely imaginary [57; 58]. In that case, the zeros that are closest to \(s=0\) can be approximated as [58]

\[\mathrm{Im}(s_{0})\simeq\sqrt{2n(2n+1)|(\hat{M}_{z}^{2n})/\!(\hat{M}_{z}^{2n+2 })|} \tag{10}\]

for large enough cumulant orders, \(n\gg 1\). Thus, in the following, we find the zeros from the high magnetization cumulants, which we calculate using neural network quantum states, and we ensure that the results from Eq. (10) are unchanged if we increase the cumulant order. We then use the scaling ansatz [57; 58]

\[\mathrm{Im}(s_{0})\simeq\mathrm{Im}(s_{0,c})+\alpha L^{-\gamma} \tag{11}\]

to predict the convergence point, \(\mathrm{Im}(s_{0,c})\), in the thermodynamic limit, where \(L\to\infty\) is the linear system size. We carry out this procedure for different magnetic fields to find the critical field, where the zeros reach \(s=0\), and the system exhibits a phase transition.

## III Results

### Extracted zeros

Figure 2 shows zeros obtained for the transverse-field Ising model in one (chain), two (square), and three (cube) dimensions. In each case, we have determined the zeros from Eq. (10) using magnetization cumulants of up to order \(n=10\) for a fixed magnetic field and a given system size. We then obtain the imaginary part of the zeros, and using the finite-size scaling ansatz from Eq. (11), we find the convergence point in the thermodynamic limit as illustrated in the figure. As an example, we see in Fig. 2a how the zeros eventually reach \(s=0\) as we decrease the magnetic field from above to \(h\simeq J\), where the system exhibits a quantum phase transition. In Figs. 2b and 2c, we show similar results for the two-dimensional square lattice and for the three-dimensional cubic lattice. For increased dimensionality, we observe that the quantum phase transitions occurs at higher magnetic fields, as expected for an increasing number of nearest neighbors. In one dimension, we use chains of up to a length of \(L=100\). For the two-dimensional square lattices, we consider systems of sizes up to \(L\times L=10\times 10\), while in three dimensions, the biggest lattice is of size \(L\times L\times L=4\times 4\times 4\). The figure includes small error bars that represent sampling errors in the neural network quantum states. We note that additional errors could potentially arise from small inaccuracies in the variational ground state.

The results for the three different geometries are combined in Fig. 3, where we show the extracted convergence points as a function of the transverse magnetic field. The extrapolation is performed by a constrained minimization of \(\mathrm{Im}(s_{0,c})\), imposing that the imaginary part is not negative. At large magnetic fields, the systems are in the paramagnetic phase with the spins mostly pointing

Figure 3: Convergence points of the zeros in the thermodynamic limit. (a) Convergence points for a linear Ising chain as a function of the magnetic field. A quantum phase transition occurs at \(h_{c}=1.00J\), where the curve exhibits a kink, and the zeros reach the real-axis. Above the critical field, the system is in the paramagnetic phase, while it is in the ferromagnetic phase below it. (b,c) Similar results for the two-dimensional square lattice (b) and the cubic lattice in three dimensions (c).
along the direction of the field. In that case, the zeros of the moment generating function do not converge to \(s=0\) in the thermodynamic limit. By contrast, as the magnetic field is lowered, the zeros eventually reach \(s=0\), signaling a quantum phase transition. Based on our calculations, we estimate the critical fields to be \(h_{c}=1.00J\) for the one-dimensional chain, \(h_{c}=3.05J\) for the two-dimensional square lattice, and \(h_{c}=5.16J\) for the three-dimensional cubic lattice. These values are all within less than \(1\%\) difference from other numerical results [59]. Below the critical field, the zeros also reach \(s=0\), since the system is in the ferromagnetic phase with spontaneous magnetization. In that case, the ground state is two fold-degenerate, and the system will exhibit an abrupt change if a small magnetic field is applied in the \(z\)-direction.

### Critical magnetic fields

We have considered other geometries in two and three dimensions as illustrated in Fig. 4, where we show results for a honeycomb lattice, a Kagome lattice, and a diamond lattice. The honeycomb lattice has two sites per unit cell, and we restrict ourselves to a linear dimension of \(L=8\), which corresponds to \(2\times L^{2}=128\) sites. Similarly, for the Kagome lattice, we go up to \(L=6\), while for the diamond lattice, we consider systems of linear size up to \(L=4\), which corresponds to \(2\times L^{3}=128\) sites. The results in Fig. 4 are qualitatively similar to those in Fig. 3, but with different critical fields. In particular, we find \(h_{c}=2.14J\) for the honeycomb lattice, \(h_{c}=2.95J\) for the Kagome lattice, and \(h_{c}=3.20J\) for the diamond lattice.

The predictions of the critical fields are summarized in Table 1, where we also show results for triangular lattices in two dimensions and face-centred cubic (FCC) and body-centred cubic (BCC) lattices in three dimensions. The results are ordered according to the dimension \(D\) as well as the number of nearest neighbors, the coordination number \(C\). In addition, we indicate the maximum linear dimension that we have used, \(L_{\rm max}\), and the number of sites in a unit cell, \(N_{\rm cell}\). Those parameters control the maximum number of spins in the lattice that we have considered, \(N_{\rm max}\). The last column contains the critical magnetic fields that we predict with the combination of Lee-Yang theory and neural network quantum states. We note that our methodology provides accurate predictions even with a rather low number of lattice sites.

## IV Discussion

### Dimensionality and lattice geometry

The importance of the lattice geometry and the dimension of the system can be understood from the results in Table 1. The chain and the honeycomb lattice, which have the lowest coordination numbers, also have the lowest critical fields. The coordination numbers are larger for the Kagome and the square lattices, where each spin has four nearest neighbors, as well as for the triangular lattice with six nearest neighbors, and we see that the critical fields increase accordingly. For the lattices in three dimensions, the coordination numbers and the critical fields are even larger. Despite this general behavior, we also see that lattices with the same dimension and coordination number (the square and Kagome lattices) still have different critical fields, which are directly related to their specific lattice geometries.

Figure 4: Convergence points of the zeros in the thermodynamic limit. (a) Convergence points for honeycomb lattice as a function of the magnetic field. A quantum phase transition occurs at \(h_{c}\approx 2.14J\), where the curve exhibits a kink, and the zeros reach the real-axis. Above the critical field, the system is in the paramagnetic phase, while it is in the ferromagnetic phase below it. (b,c) Similar results for the Kagome lattice (b) and the diamond lattice (c).
### Mean-field approximation

To better understand the role of the coordination number, we show in Fig. 5 the critical fields as a function of the coordination number. In Fig. 5a, we see the clear trend that the critical fields increase with the coordination number. Indeed, within a simple mean-field approximation, we would expect that the critical field is directly related to the coordination number as \(h_{c}^{\rm MF}=CJ\)[74]. We show this mean-field approximation with a dashed line in the figure and find good qualitative agreement with our predictions. We also see that our results come closer to the mean-field approximation as the dimension of the system is increased. In particular, it is clear that the critical field for the one-dimensional chain is furthest away from the mean-field approximation, while the results for the three-dimensional lattices are much closer.

To further support these observations, we show in Fig. 5b the ratio of the critical fields over the mean-field approximation. This ratio allows us to characterize how the relative deviations from the mean-field prediction decrease for larger coordination numbers. Still, we see that the critical fields are all smaller than the mean-field approximation, which ignores quantum fluctuations. The results for the critical fields in three dimensions are closer to the mean-field approximation as compared with one and two dimensions. This observation is in line with the expectation that mean-field theory becomes more accurate in higher dimensions.

## V Conclusions

We have combined a Lee-Yang theory of quantum phase transitions with neural network quantum states to predict the critical field of the transverse-field Ising model in different dimensions and lattice geometries. Specifically, we have used neural network quantum states to find the ground state of the interacting spin system, which further makes it possible to extract the cumulants of the magnetization. From these cumulants, we determine the complex zeros of the moment-generating function, which reach the real-axis in the thermodynamic limit if the system exhibits a phase transition. Our method works with rather small systems, which in turn allows us to treat lattices in two and three dimensions. Our predictions agree well with results that were obtained using large-scale quantum many-body methods. We have also analyzed the differences between our predictions and a simple mean-field approximation, which becomes increasingly accurate for higher coordination numbers and dimensions. Thanks to the flexibility of neural network quantum states, the method can potentially treat frus

Figure 5: Comparison with mean-field theory. (a) The critical fields are shown as functions of the coordination number, \(C\). The dashed line is a simple mean-field approximation that directly links the critical field to the coordination number as \(h_{c}^{\rm MF}=CJ\). (b) The ratio of the critical fields over the mean-field approximation as functions of the coordination number, \(C\). For large coordination numbers and dimensions, the critical fields approach the mean-field approximation indicated with a dashed line.

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline Lattice & \(D\) & \(C\) & \(L_{\rm max}\) & \(N_{\rm cell}\) & \(N_{\rm max}\) & \(h_{c}/J\) \\ \hline Chain & 1 & 2 & 60 & 1 & 60 & 1.00 \\ Honeycomb & 2 & 3 & 8 & 2 & 128 & 2.14 \\ Kagome & 2 & 4 & 6 & 3 & 108 & 2.95 \\ Square & 2 & 4 & 10 & 1 & 100 & 3.05 \\ Triangular & 2 & 6 & 10 & 1 & 100 & 4.78 \\ Diamond & 3 & 4 & 4 & 2 & 128 & 3.20 \\ Cubic & 3 & 6 & 4 & 1 & 64 & 5.16 \\ BCC & 3 & 8 & 4 & 1 & 64 & 7.10 \\ FCC & 3 & 12 & 4 & 1 & 64 & 10.8 \\ \hline \end{tabular}
\end{table}
Table 1: Summary of critical fields. For each lattice, we indicate the dimension, \(D\), and the coordination number, \(C\). We also show the maximum linear dimension, \(L_{\rm max}\), and the number of sites per unit cell, \(N_{\rm cell}\), which give the maximum number of sites that we have used as \(N_{\rm max}=N_{\rm cell}\times L_{\rm max}^{D}\). The last column contains our predictions of the critical field.
trated problems, in stark contrast to quantum Monte Carlo approaches that suffer from sign-problems. Our results show that the combination of Lee-Yang theories of phase transitions with neural network quantum states provides a viable way forward to predict the phase behavior of complex quantum many-body systems such as Heisenberg models and fermionic Hubbard models.

###### Acknowledgements.

 We acknowledge the computational resources provided by the Aalto Science-IT project and the support from the Finnish National Agency for Education (Opetushallitus), the Academy of Finland through grants (Grants No. 331342 and No. 336243) and the Finnish Centre of Excellence in Quantum Technology (Projects No. 312057 and No. 312299), and from the Jane andatos Erkko Foundation.

## Appendix A Details of calculations

All of our calculations were implemented in Netket 3.3 [68, 75]. In one dimension, we found that a restricted Boltzmann machine works well, while in two dimensions, a group convolutional neural network functions better. In three dimensions, we used a simple and shallow symmetric architecture with real weights, which is sufficient, since the transverse-field Ising model is stoquastic.

In one dimension, we used a simple real restricted Boltzmann machine with a number of hidden units per visible unit of \(\alpha=20\). For each training iteration, 8192 samples were used, taken from 128 parallel chains. The network was trained for 3000 iterations with a learning rate of 0.02, and then for further 1000 iterations with a learning rate of 0.01. Stochastic reconfiguration with a diagonal shift of 0.01 was used.

In two dimensions, we used a group convolutional neural network [76, 15] defined over the group of all translations with four layers of feature dimension 8 each and complex parameters. We used 32 parallel Markov chains constructed using a Metropolis algorithm with local updates, and we took 1024 samples per iteration step. Stochastic reconfiguration with a diagonal shift of 0.01 was used, and the network was trained with a learning rate of 0.01 for 2000 iterations. If necessary, we trained the network multiple times and chose the network with the lowest variance of the energy.

In three dimensions, we applied a dense symmetric layer with real weights and 40 features to the input, and we then activated it with the ReLu function, which was then summed over to obtain the wave function. We used a local Metropolis update Markov-chain with 128 parallel chains and 8192 samples per training step. A learning rate of 0.002 and stochastic reconfiguration with a diagonal shift of 0.01 were applied. We then trained the network for 2000 iterations. If necessary, we ran this training multiple times for the same configuration (system size and magnetic field), and we chose the network parameters that resulted in the lowest variance of the ground state energy, so that the network was as similar as possible to an eigenstate of the Hamiltonian.

We evaluated the moments of the magnetization using regular sampling with an unbiased Markov chain, since

\[\langle\hat{M}_{z}^{n}\rangle=\sum_{\bar{\sigma}}P_{\psi}(\bar{\sigma})M_{z}^ {n}=\sum_{\bar{\sigma}}P_{\sigma}(\bar{\sigma})\left[\sum_{i}\sigma_{i}\right] ^{n}. \tag{10}\]

For the two- and three-dimensional lattices with up to \(N_{\text{max}}=128\) sites, we took \(100\times 1024\times 128\simeq 13\) million samples. For the one-dimensional lattice, we took up to \(1000\times 1024\times 128\simeq 270\) million samples. For the sampling, we used 128 parallel chains and discarded the first 64 entries. From the moments, we then obtained the cumulants using a standard recursion relation between them.
