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.


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 manybody 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]29]. Unlike many other Monte Carlo methods, neural network quantum states, as a variant of variational Monte Carlo, can be applied to frustrated spin models. Moreover, their use has been extended to fermionic [22,30] and bosonic [31][32][33] systems, as well as to molecules [22,23] 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 ex-FIG. 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 × 3 × 3. (b) A neural network takes a configuration of the spins, encoded in the vector ⃗ σ = (σ1, ..., σ N ), and outputs the corresponding value of the wave function, ψ⃗ θ (⃗ σ) = ⟨⃗ σ|ψ⟩, which depends on the variational parameters in ⃗ θ. (c) From the fluctuations of the magnetization, we extract the zeros of the moment generating function of the magnetization (shown as ×) and investigate their motion in the complex plane as we increase the system size (indicated with numbers). Above the critical field, h > hc, the zeros remain complex in the thermodynamic limit, signaling that the system is in the paramagnetic phase (PM). For h < hc, the system is in the ferromagnetic phase (FM) with finite magnetization and the zeros approach the origin in the thermodynamic limit. (d) Phase diagram of the transverse-field Ising model. ample an external magnetic field or the inverse temperature [34][35][36][37]. 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 axis in the thermodynamic limit of large system sizes, giving rise to a non-analytic behavior of the free energy density [38][39][40][41][42][43][44][45][46][47][48][49]. Lee-Yang zeros are not just a theoretical concept, but they can also be determined experimentally [50][51][52][53][54]. In recent years, applications of Lee-Yang theory have been expanded to dynamical quantum phase transitions in quantum many-body systems after a quench [55][56][57] and to quantum phase transitions in systems at zero temperature [58,59].
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 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.

A. Transverse-field Ising model
We consider the transverse-field Ising model on a lattice of spin-1/2 sites as described by the Hamiltonian 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 in the thermodynamic limit 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 [60].

B. 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 where the vector ⃗ θ 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 timeevolution of quantum systems [8,61,62], quantum state tomography [63][64][65], as well as finite-temperature equilibrium physics [66][67][68]. 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 where we have defined the probability and the local spin Hamiltonian 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 [69]. 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 sinceM z is diagonal in the computational basis, such that M n z (⃗ σ) = (⟨⃗ σ|M z |⃗ σ⟩) n = ⟨⃗ σ|M n z |⃗ σ⟩. Additional details of these calculations are provided in Appendix A.

C. 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 [34][35][36][37]. 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 [38,41,42,[45][46][47], atomic physics [39] and particle physics [40,43,44,48,[70][71][72][73].
Recently, it has been extended to the zeros of the moment generating function that describes the fluctuations of the order parameter [58,59] and thereby allows for the detection of quantum phase transitions. In Ref. [58], the method was implemented for one-dimensional chains, in particular the transverse-field Ising chain and an anisotropic Heisenberg model. In addition, preliminary results for a J 1 -J 2 Heisenberg model in two dimensions were presented. It was observed that the Lee-Yang approach made it possible to determine the critical points using rather short chains. Following this approach, we now define the moment generating function whereM 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 |ψ Lee-Yang theory, and the cumulant generating function, Θ(s) = ln χ(s), becomes the corresponding free energy. The moments and cumulants of the magnetization are given by derivatives with respect to the counting field as and ⟨⟨M n z ⟩⟩ = ∂ n s Θ(s)| s=0 .
Importantly, away from a phase transition, the cumulants are expected to grow linearly with the system size, such that the normalized cumulants ⟨⟨M n z ⟩⟩/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 [41,74]. This nonanalytic behavior emerges in the thermodynamic limit, 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. [41,53,54,58,59]. 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 [41,53,58,59]. 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,Û †ĤÛ =Ĥ, with respect to the unitary operatorÛ = ∏ iσ x i that flips all spins, implies that all odd cumulants vanish, and in this model the zeros are purely imaginary [58,59]. In that case, the zeros that are closest to s = 0 can be approximated as [59] Im(s 0 ) ≃ for large enough cumulant orders, n ≫ 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 [58,59] Im(s 0 ) ≃ Im(s 0,c ) + αL −γ (11) to predict the convergence point, Im(s 0,c ), in the thermodynamic limit, where L → ∞ 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
A. 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 ≃ 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 × L = 10 × 10, while in three dimensions, the biggest lattice is of size L×L×L = 4×4×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 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 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 twodimensional square lattice, and h c = 5.16J for the threedimensional cubic lattice. These values are all within less than 1% difference from other numerical results [60]. 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 folddegenerate, and the system will exhibit an abrupt change if a small magnetic field is applied in the z-direction.

B. 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 × 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 × 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 I, 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 max , and the number of sites in a unit cell, N cell . Those parameters control the maximum number of spins in the lattice that we have considered, N 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.

A. Dimensionality and lattice geometry
The importance of the lattice geometry and the dimension of the system can be understood from the results in Table I. 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.

B. 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 MF c = CJ [75]. The physical picture here is that the spins experience a competition between two opposing effects. On the one hand, the external magnetic field tends to align the spins in the x-direction. On the other hand, the coupling between them tends to align them along the z-axis. The phase transition then occurs when the two effects are equally strong. Since the coupling between neighboring spins is proportional to the coordination number, so is the critical field. Still, this mean-field approximation does not fully capture the actual physics, since it ignores the effects of the lattice geometry. We show the 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 threedimensional lattices are much closer.
To further support these observations, we show in Fig. 5b the ratio of the critical fields over the meanfield 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 meanfield 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. We also show the maximum linear dimension, Lmax, and the number of sites per unit cell, N cell , which give the maximum number of sites that we have used as Nmax = N cell × L D max . The two last columns contain our predictions of the critical field as well as the results from Ref. [60].

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 frustrated 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. The application of neural network quantum states to fermionic models is currently being developed, which in the future may provide a better understanding of interacting fermionic systems using Lee-Yang theory. All of our calculations were implemented in Netket 3.3 [69,76]. 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 α = 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 [15,77] 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 energy, so that the network was as similar as possible to a ground state of the Hamiltonian.
We evaluated the moments of the magnetization using regular sampling with an unbiased Markov chain, since For the two-and three-dimensional lattices with up to N max = 128 sites, we took 100 × 1024 × 128 ≃ 13 million samples. For the one-dimensional lattice, we took up to 1000 × 1024 × 128 ≃ 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.

Appendix B: Error analysis
The position of the zeros can be slightly imprecise for several reasons. For example, the ground state is never completely accurate, which may lead to errors in the zeros. The Monte Carlo sampling itself may also lead to inaccuracies. In addition, one has to ensure that only the closest zeros contribute to the cumulants by using a sufficiently high cumulant order. The errors from the Monte Carlo sampling are statistical in nature and can easily be quantified. Regarding the ground state, we make sure that it has converged so that the zeros obtained with the cumulant method remained unchanged. In Fig. 6(a) we show the relative error in the ground state energy for the square lattice, and we see that they are small for all field strengths. In Fig. 6(b), we show the extracted zero for different relative errors, and we see how the position of the zero converges as the error is reduced. These results were obtained by running the algorithm for finding the ground state ten times, taking ten snapshots each time at different stages of the training. These 100 data points have different errors in the energy as shown in Fig. 6(b) together with the extracted zero. As the error is reduced, we see a clear convergence of the zero. In particular, for ∆E/|⟨E⟩| < 0.001, the position of the zero remains the same as the error is further reduced.