Neural Gutzwiller-projected variational wave functions

TL;DR

This paper introduces neural network-augmented Gutzwiller-projected variational wave functions to improve the accuracy of ground state approximations in complex quantum spin models, especially in frustrated regimes.

Contribution

It combines neural networks with fermionic Gutzwiller-projected wave functions to enhance their effectiveness in describing both ordered and disordered quantum spin phases.

Findings

Neural augmentation improves accuracy in frustrated spin regimes.

RBM states excel for simple sign structures like Nel states.

Significant improvements are observed in highly frustrated regimes.

Abstract

Variational wave functions have enabled exceptional scientific breakthroughs related to the understanding of novel phases of matter. Examples include the Bardeen-Cooper-Schrieffer theory of superconductivity, the description of the fractional quantum Hall effect through the Laughlin state, and Feynman's variational understanding of large-scale quantum effects in liquid Helium. More recently, Gutzwiller-projected wave functions, typically constructed from fermionic degrees of freedom, have been employed to examine quantum spin models in the presence of competing interactions, where exotic phases with no spontaneous symmetry breaking and fractional excitations may exist. In this work, we investigate the aforementioned fermionic wave functions supplemented with neural networks, specifically with the so-called restricted Boltzmann machine (RBM), to boost their accuracy and obtain reliable approximations to the ground state of generic spin models. In particular, we apply our neural augmented fermionic construction to the description of both magnetically ordered and disordered phases of increasing complexity, including cases where the ground state displays a non-trivial sign structure. Even though the RBM state is by far more effective for N\'eel states endowed with a particularly simple sign structure, it provides a significant improvement over the original fermionic state in highly frustrated regimes where a complex sign structure is anticipated, thus marking the path to an understanding of strongly-correlated spin models on the lattice via neural Gutzwiller-projected variational wave functions.

Equations

$$\mathcal{H}=J_1\sum _{⟨i,j⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+J_2\sum _{⟨⟨i,j⟩⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$\mathcal{H}=J\sum _{⟨i,j⟩}(S_i^x{S}_j^x+S_i^y{S}_j^y).$$

$$|{\Psi }_{\mathcal{C}}⟩=\hat{\mathcal{C}}|{\Psi }_f⟩=\sum _{\sigma }\mathcal{C}(\sigma )⟨\sigma |{\Psi }_f⟩|\sigma ⟩.$$

$${\mathbf{S}}_i=\frac{1}{2}\sum _{\alpha ,\beta }{f}_{i,\alpha }^{†}{\mathbf{\sigma }}_{\alpha ,\beta }{f}_{i,\beta }.$$

$$|{\Psi }_f⟩={\mathcal{P}}_{S_{tot}^z=0}{\mathcal{P}}_G|{\Phi }_0⟩.$$

$${\mathcal{H}}_0={\mathcal{H}}_{\mathrm{BCS}}+{\mathcal{H}}_{\mathrm{AF}},$$

$${\mathcal{H}}_{\mathrm{BCS}}=\sum _{i,j,\sigma }{t}_{i,j}{f}_{i,\sigma }^{†}{f}_{j,\sigma }^{}+\sum _{i,j}{\Delta }_{i,j}{f}_{i,\uparrow }^{†}{f}_{j,\downarrow }^{†}+h.c..$$

$${\mathcal{H}}_{\mathrm{AF}}={\Delta }_{\mathrm{AF}}\sum _i\left(e^{iQ\cdot R_i}{f}_{i,\uparrow }^{†}{f}_{i,\downarrow }^{}+e^{-iQ\cdot R_i}{f}_{i,\downarrow }^{†}{f}_{i,\uparrow }^{}\right).$$

$${\mathcal{C}}_{\mathrm{Jastrow}}(\sigma )=\exp \left(\frac{1}{2}\sum _{i,j}{v}_{i,j}{\sigma }_i^z{\sigma }_j^z\right).$$

$$E_{\mathrm{RBM}}=\sum _{i=1}^N\sum _{\alpha =1}^{N_{\alpha }}{h}^{\alpha }{W}_i^{\alpha }{\sigma }_i^z+\sum _{\alpha =1}^{N_{\alpha }}{b}^{\alpha }{h}^{\alpha }+\sum _{i=1}^N{a}_i{\sigma }_i^z.$$

$${\mathcal{C}}_{\mathrm{RBM}}(\sigma )=\exp \left[\sum _{\alpha }\log \cosh \left(b^{\alpha }+\sum _i{W}_i^{\alpha }{\sigma }_i^z\right)\right].$$

$${\mathcal{C}}_{\mathrm{tRBM}}(\sigma )=\prod _R{C}_{\mathrm{RBM}}[{\mathcal{T}}_R(\sigma )]$$

$$=\exp \left[\sum _R\sum _{\alpha }\log \cosh \left(b^{\alpha }+\sum _i{W}_i^{\alpha }{\sigma }_{i{+}_R}^z\right)\right].$$

$${\mathcal{C}}_{\mathrm{sRBM}}(\sigma )=\prod _{\Sigma }\prod _R{C}_{\mathrm{RBM}}[\Sigma {\mathcal{T}}_R(\sigma )]$$

$$=\exp \left[\sum _{\Sigma }\sum _R\sum _{\alpha }\log \cosh \left(b^{\alpha }+\sum _i{W}_i^{\alpha }{\sigma }_{\Sigma (i{+}_R)}^z\right)\right],$$

$$E_{\mathcal{C}}=⟨\mathcal{H}{⟩}_{\mathcal{C}}=\frac{⟨{\Psi }_{\mathcal{C}}|\mathcal{H}|{\Psi }_{\mathcal{C}}⟩}{⟨{\Psi }_{\mathcal{C}}|{\Psi }_{\mathcal{C}}⟩},$$

$${\mathcal{O}}_k(\sigma )=\frac{\partial \log \left[⟨\sigma |{\Psi }_{\mathcal{C}}⟩\right]}{\partial {\omega }_k}=\frac{1}{⟨\sigma |{\Psi }_{\mathcal{C}}⟩}\frac{\partial ⟨\sigma |{\Psi }_{\mathcal{C}}⟩}{\partial {\omega }_k}.$$

$$⟨\sigma |{\hat{E}}_{\mathrm{loc}}|{\sigma }^{\prime }⟩=E_{\mathrm{loc}}(\sigma ){\delta }_{\sigma ,{\sigma }^{\prime }}=\frac{⟨\sigma |\mathcal{H}|{\Psi }_{\mathcal{C}}⟩}{⟨\sigma |{\Psi }_{\mathcal{C}}⟩}{\delta }_{\sigma ,{\sigma }^{\prime }}.$$

$${\omega }_k^{\prime }={\omega }_k+\eta \sum _{k^{\prime }}{\mathcal{S}}_{k,k^{\prime }}^{-1}{f}_{k^{\prime }}.$$

$${\mathcal{S}}_{k,k^{\prime }}=\Re \left[⟨{\hat{\mathcal{O}}}_k^{†}{\hat{\mathcal{O}}}_{k^{\prime }}^{}{⟩}_{\mathcal{C}}-⟨{\hat{\mathcal{O}}}_k^{†}{⟩}_{\mathcal{C}}⟨{\hat{\mathcal{O}}}_{k^{\prime }}^{}{⟩}_{\mathcal{C}}\right],$$

$$f_k=-\frac{\partial E_{\mathcal{C}}}{\partial {\omega }_k}=-2\Re \left[⟨\mathcal{H}{\hat{\mathcal{O}}}_k{⟩}_{\mathcal{C}}-⟨\mathcal{H}{⟩}_{\mathcal{C}}⟨{\hat{\mathcal{O}}}_k{⟩}_{\mathcal{C}}\right].$$

$$\Delta E=\left|\frac{E_{\mathcal{C}}-E_0}{E_0}\right|,$$

$${\mathcal{S}}_R^{xy}\equiv \frac{1}{2N}\sum _{i=1}^N\left(S_i^x{S}_{i+R}^x+S_i^y{S}_{i+R}^y\right),$$

$${\mathcal{S}}_R^z\equiv \frac{1}{N}\sum _{i=1}^N{S}_i^z{S}_{i+R}^z.$$

$$\Delta {\mathcal{S}}_R^{\alpha }=\frac{⟨{\mathcal{S}}_R^{\alpha }{⟩}_{\mathcal{C}}-⟨{\mathcal{S}}_R{⟩}_0}{⟨{\mathcal{S}}_R{⟩}_0},$$

Full text

I Introduction

Fractionalization refers to a set of phenomena in quantum many-body physics where a collection of strongly interacting microscopic degrees of freedom, such as electrons or spins on a lattice, break down into multiple quasiparticles. One of the most prominent examples is the fractional quantum Hall effect Tsui et al. (1982), where the constituent electrons decay into quasiparticles carrying fractions of the electron charge. Another prime example of fractionalization arises in the study of spin liquids Balents (2010); Savary and Balents (2016), which are highly-entangled states of matter that host a wide array of exotic phenomena such as emergent gauge structures, non-local excitations, and fractional quasiparticles. These quasiparticles, often termed “spinons”, behave as fractions of ordinary magnons and spin waves.

Traditionally, gauge theories have played a central role in the study of quantum spin liquids since their particular structure can capture the highly peculiar non-local excitations and entanglement properties present in these exotic phases Balents (2010); Savary and Balents (2016); Wen (2004). To construct such gauge theories, physicists rely on representations of the original degrees of freedom in terms of pairs of bosons or fermions, whose local Hilbert space is constrained to a subspace corresponding to the original Hilbert space of the spin system. This reformulation results in a rich theory of quasiparticles coupled to gauge fields with the potential to shed light onto the properties of the original system.

Since a gauge theoretical reformulation of a quantum spin system is formally exact, an explanation of all physically conceivable phases in the original models is always possible. Whereas conventional phases of matter, such as Néel antiferromagnets and valence bond solids, are associated with the confining phases of the gauge theory, fractionalization is associated with the deconfined phases where the spinons materialize as legitimate excitations of the system. Thus, the hope is that this reformulation may help us understand exotic phases which naturally realize the quasiparticles and the emerging gauge fields of the theory. Originally suggested by Baskaran and Anderson in the context of high temperature superconductivity Baskaran and Anderson (1988), fractionalization and emergent gauge structures have been demonstrated in several exactly solvable systems such as Kitaev’s toric code Kitaev (2003) and honeycomb model Hermanns et al. (2018), string-net condensates Levin and Wen (2005), as well confirmed in microscopic models Cincio and Vidal (2013); Isakov et al. (2011); Hickey and Trebst (2019), variational wave functions Zhang et al. (2012), and also experiments Tsui et al. (1982).

In practice, however, gauge theories of quantum spin systems are typically intractable. In such cases, the generic procedure, also known as the “parton” construction, is to rewrite the original model in terms of canonical bosons or fermions, perform a mean-field decoupling of the gauge theory, and then study the implications of the mean-field theory and its stability against fluctuations Wen (2004). Alternatively, a numerical approach can be introduced, where wave functions from a mean-field treatment of the gauge theory are numerically projected back to the original Hilbert space through the variational Monte Carlo method (VMC) Becca and Sorella (2017).

The accuracy and degree of applicability of these states, here referred to as Gutzwiller projected wave functions, have been demonstrated in several prototypical models of frustrated magnetism Sheng et al. (2009); Iqbal et al. (2011, 2013); Gong et al. (2013); Hu et al. (2015); Iqbal et al. (2016). The main advantage of this strategy is that the states it produces can lead to very accurate estimates of the ground-state energy but their most important characteristic is the interpretability, as these states often provide us with a physically transparent description of the phases of the model under consideration. On the other hand, the main disadvantage of this approach is the lack of systematic ways of improving the quality of the approximations the construction introduces. While the parton construction remains a powerful approach to quantum spin liquids, it is clear that the Gutzwiller projector only partially reintroduces the missing gauge fluctuations of the mean-field treatment, which may not be sufficient to capture important long wavelength properties of the system. This effect is particularly important in special cases, e.g., in $U(1)$ spin liquids, since Gutzwiller projected states do not account for crucial spatial gauge fluctuations Tay and Motrunich (2011); Liu et al. (2014).

Here, we explore the possibility of improving the accuracy of the parton construction through neural networks. In particular, because of their extensive use in condensed matter physics and quantum information Carleo and Troyer (2017); Chen et al. (2018); Torlai et al. (2018); Torlai and Melko (2018); Deng et al. (2017a); Gao and Duan (2017); Deng et al. (2017b); Nomura et al. (2017); Clark (2018), we consider restricted Boltzmann machines (RBM), which, due to the nature of their non-local structure, have been shown to represent some highly-entangled many-body states using a relatively small number of parameters Deng et al. (2017a); Gao and Duan (2017); Chen et al. (2018); Glasser et al. (2018). The family of states we consider can be thought of as quantum mechanical version of the so-called product of experts idea used in probabilistic models in machine learning Hinton (2002), where the core strategy is to combine several simpler wave functions (the “experts”) by multiplying their amplitudes in a certain computational basis Becca and Sorella (2017); Toulouse and Umrigar (2008); Carrasquilla et al. (2013); Nomura et al. (2017); Glasser et al. (2018). Concretely, the family of states we consider is written as a product of a Gutzwiller-projected fermionic state and a complex-valued RBM. By virtue of the high representation power of the RBM Le Roux and Bengio (2008); Gao and Duan (2017), the expectation is that the RBM may help us improve the accuracy of the Gutzwiller projected states and in part alleviate the problem of the missing gauge fluctuations in these family of states.

The paper is organized as follow: In Sec. II, we introduce the spin models that are considered in this study; in Sec. III, we define the variational wave functions; in Sec. IV, we briefly describe the numerical optimization of the wave functions; in Sec. V, we discuss the results; finally, in Sec. VI, we draw our conclusions.

II Spin models

We consider two-dimensional (2D) quantum lattice models that are traditionally studied in condensed matter physics. These models host a wide array of effects that are relevant to the understanding of physical phenomena such as low-temperature experiments with superfluid helium Daunt and Smith (1954), Fe-based superconducting materials Kamihara et al. (2008), triangular lattice compounds Paddison et al. (2017), fractional excitations in the quantum antiferromagnets Dalla Piazza et al. (2015), among many others. We begin by introducing the $J_{1}-J_{2}$ Heisenberg model, an archetypical model of frustration in quantum antiferromagnets. Its Hamiltonian is given by

[TABLE]

where $\mathbf{S}_{i}=(S^{x}_{i},S^{y}_{i},S^{z}_{i})$ are spin-1/2 degrees of freedom sitting on the sites of a square lattice and $J_{1},J_{2}\geq 0$ characterize the antiferromagnetic couplings between the magnetic degrees of freedom. Here, $\langle\cdots\rangle$ and $\langle\langle\cdots\rangle\rangle$ restrict the summations over first- and second-nearest neighboring sites, respectively. Throughout this paper, we focus on two different values of the frustrating ratio, namely $J_{2}/J_{1}=0$ and $J_{2}/J_{1}=0.5$, to assess the accuracy of the aforementioned variational wave functions in two different regimes. For $J_{2}=0$ the model reduces to the (unfrustrated) Heisenberg model, whose ground state displays long-range Néel magnetic order Sandvik (1997); Calandra Buonaura and Sorella (1998), At finite $J_{2}$, the energetic competition between the antiferromagnetic couplings at first- and second-nearest distances introduces frustration, resulting in a reduction of the magnetic order. In particular, the system exhibits a high degree of frustration around $J_{2}/J_{1}\approx 0.5$, where the nature of the ground state of the model is a long standing problem Read and Sachdev (1989); Schulz et al. (1996). Different scenarios have been suggested, such as the onset of a valence-bond solid phase, with columnar Poilblanc and Mambrini (2017); Haghshenas and Sheng (2018) or plaquette Mambrini et al. (2006); Gong et al. (2014) order, and/or the existence of a spin-liquid phase Hu et al. (2013); Morita et al. (2015); Wang and Sandvik (2018).

Due to its relevance to physical situations ranging from low-temperature experiments with superfluid helium Daunt and Smith (1954), magnetic insulators Giamarchi et al. (2008), and ultracold gases in optical lattices Bloch et al. (2008); Carrasquilla and Rigol (2012), we also consider the $XY$ model on the square lattice, whose Hamiltonian is given by

[TABLE]

In this case, the ground state exhibits Néel order with the spins aligned in the $xy$ plane. A fundamental difference between the XY and the unfrustrated Heinseberg models comes from their symmetries: while the Heisenberg exchange interaction is invariant under global spin rotations around any axis (i.e., it exhibits $SU(2)$ symmetry), the $XY$ coupling is invariant only under global rotations around $S_{z}$ (leading to a $U(1)$ symmetry).

Finally, we consider the Heisenberg model on the triangular lattice, a model of geometric frustration relevant to the understanding of organic Mott insulators displaying spin liquid behavior Shimizu et al. (2003). Its Hamiltonian is given by Eq. (1) with $J_{2}=0$. Due to frustration, this system exhibits strong quantum fluctuations. Even though this model was initially proposed as a candidate host for the resonating valence-bond spin liquid Anderson (1973), which motivated the entire field of spin liquids, numerical evidence Capriotti et al. (1999); White and Chernyshev (2007) suggested that the ground state of the system is magnetically ordered with a $120^{\circ}$ pattern, which is the accepted consensus.

III Construction of the wave functions

In general, the variational Ansätze employed in this work can be written as a product of a Gutzwiller-projected fermionic wave function, $|\Psi_{f}\rangle$, and a many-body correlator, $\hat{\mathcal{C}}$:

[TABLE]

In the above formula we have inserted a resolution of the identity ($\hat{\mathds{1}}=\sum_{\sigma}|\sigma\rangle\langle\sigma|$) and exploited the fact that $\hat{\mathcal{C}}$ is diagonal in the many-body computational basis $\{|\sigma\rangle\}$.

In order to construct the fermionic part of the Ansatz we turn to a parton construction, representing the spin degrees of freedom in terms of Abrikosov fermions Wen (2004):

[TABLE]

Here $\bm{\sigma}=(\sigma^{x},\sigma^{y},\sigma^{z})$ is the vector of Pauli matrices and the operator $f^{\dagger}_{i,\alpha}$ ($f^{\phantom{\dagger}}_{i,\alpha}$) creates (annihilates) a fermion with spin $\alpha$ at site $i$. Within the Abrikosov picture, the commutation relations among spins are a consequence of the fermionic anticommutation relations. However, the fermionic representation of spins enlarges the Hilbert space of the system, allowing each lattice site to be in four different states, namely $|0\rangle$, $|\uparrow\rangle$, $|\downarrow\rangle$, and $|\uparrow\downarrow\rangle$. Therefore, the mapping of Eq. (4) is a faithful representation of spin degrees of freedom only if the fermions are restricted to the subspace of configurations with one fermion per site. This constraint can be enforced by the Gutzwiller projection operator $\mathcal{P}_{G}=\prod_{i}(n_{i,\uparrow}-n_{i,\downarrow})^{2}$ (where $n_{i,\sigma}=f^{\dagger}_{i,\sigma}f^{\phantom{\dagger}}_{i,\sigma}$). The wave function $|\Psi_{f}\rangle$ is obtained by Gutzwiller-projecting a certain fermionic state $|\Phi_{0}\rangle$ (see below):

[TABLE]

In addition to the Gutzwiller projector, here we consider a second projector, $\mathcal{P}_{S^{z}_{tot}=0}$, which restricts the wave function to the subspace of configurations with ${S^{z}_{tot}=\sum_{i}S^{z}_{i}=0}$. We emphasize the fact that both projections are treated exactly by sampling only the fermionic configurations which satisfy the desired constraints. Finally, $|\Phi_{0}\rangle$ is obtained by computing the ground state of a quadratic (i.e., a mean-field) Hamiltonian ($\mathcal{H}_{0}$) of Abrikosov fermions,

[TABLE]

which we have split into two parts for clarity. The first part is a Bardeen-Cooper-Schrieffer (BCS) Hamiltonian containing a (complex) hopping term ($t_{i,j}$) and a singlet pairing term ($\Delta_{i,j}=\Delta_{j,i}$):

[TABLE]

The above Hamiltonian is invariant under global spin rotations and, therefore, the spin wave functions obtained by the Gutzwiller projection of its ground states are $SU(2)$ symmetric. Ansätze of these form are particularly suited to describe magnetically disordered phases of matter, such as quantum spin liquids and valence-bond solids Savary and Balents (2016). The Gutzwiller-projected BCS wave functions fall back into the class of resonating valence-bond states, first introduced by Anderson in the context of high-Tc superconductivity Anderson (1987). The second term of Eq. (6), ${\cal H}_{\rm AF}$, contains a magnetic field in the $xy$-plane which induces magnetic order in the variational Ansatz:

[TABLE]

The vector $Q$ determines the periodicity of the magnetic order: the Néel phase on the square lattice is obtained by setting $Q=(\pi,\pi)$, while the $120^{\circ}$ phase on the triangular lattice corresponds to $Q=(\frac{4\pi}{3},0)$ [or, equivalently, $Q=(\frac{2\pi}{3},\frac{2\pi}{\sqrt{3}})$]. Typically, the wave function obtained by Gutzwiller-projecting the ground state of ${\cal H}_{\rm AF}$ overestimates the magnetic order parameter Becca et al. (2011). To improve this construction, further quantum fluctuations can be added by including hopping terms in the BCS Hamiltonian and by applying a two-body Jastrow factor correlator to the fermionic wave function [as in Eq. (3)]:

[TABLE]

The Jastrow factor has the effect of adding correlations which are perpendicular to the in-plane magnetic field $\Delta_{\rm{AF}}$. In this work, we consider Jastrow factors with translationally invariant long-range pseudopotentials ${v_{i,j}=v(|R_{i}-R_{j}|)\in\mathbb{R}}$. We emphasize the fact that all the coupling constants defining ${\cal H}_{0}$ (i.e. $t_{i,j},\Delta_{i,j},\Delta_{\rm{AF}}$) and the Jastrow pseudopotentials play the role of variational parameters, which are optimized in order to find the best variational energy. The optimal parametrization of the auxiliary Hamiltonian ${\cal H}_{0}$ depends on the model under investigation.

The central aim of this work is to explore the possibility of improving Gutzwiller-projected fermionic wave functions by applying a stronger many-body correlator than the two-body Jastrow factor. For this purpose, a neural network is employed in the form of a restricted Boltzmann machine. This network is defined by introducing a set of auxiliary Ising variables, $\{h^{\alpha}\}_{\alpha=1,...,N_{\alpha}}$, which form the so-called hidden layer. These variables are coupled to the z-components of the spins of the lattice (dubbed as visible layer, $\{\sigma_{i}^{z}\}_{i=1,...,N}$) through a classical energy functional of the form:

[TABLE]

The RBM correlator is then obtained by computing the Boltzmann factor $e^{E_{\mathrm{RBM}}}$ and taking its trace over the hidden variables degrees of freedom. This operation can be performed exactly due to the particular form of the classical energy functional of Eq. (10), which only contains interactions between variables belonging to the two different layers (i.e., no intralayer couplings are allowed). Since we consider translationally invariant states with conserved magnetization $S^{z}_{tot}=0$, we can set $a_{i}=0$ and obtain the final form of the RBM correlator Carleo and Troyer (2017):

[TABLE]

The parameters of $\hat{\mathcal{C}}_{\mathrm{RBM}}$ are called biases ($b^{\alpha}$) and weights ($W^{\alpha}_{i}$), and are (in general) assumed to be complex numbers in this work. The complex parametrization of the RBM allows the correlator to change both the amplitudes and the phases of the fermionic wave function to which it is applied. The expression of Eq. (11) can be regarded as a sort of many-body Jastrow factor, since a series expansion of the $\log\cosh(\dots)$ function contains the n-body terms of the $\sigma^{z}$ variables. Unlike the Jastrow factor, the RBM correlator not only breaks the $SU(2)$ symmetry of spin, but also the $Z_{2}$ symmetry $\sigma^{z}\mapsto-\sigma^{z}$. This happens if the biases $b^{\alpha}$ are nonzero, since the aforementioned expansion can contain products of odd numbers of spins.

An important question to address is the implementation of lattice symmetries in the correlator $\hat{\mathcal{C}}$. As already pointed out, a symmetric two-body Jastrow factor can be obtained by simply taking a symmetric pseudopotential, e.g., $v_{i,j}=v(|R_{i}-R_{j}|)$. This procedure, however, cannot be applied to the RBM correlator, since its parameters depend upon the index $\alpha$ (labelling the hidden units), which does not have any physical meaning. Therefore, the most straightforward way of including symmetries is to implement them a posteriori. Concretely, if we want to enforce translational symmetry, we symmetrize the RBM correlator through a product over all possibile Bravais lattice translations $\{\mathcal{T}_{R}\}$:

[TABLE]

The above expression is translationally invariant with momentum $K=(0,0)$. In addition, the point group symmetries $\{\Sigma\}$ of the lattice can be implemented on top of the translationally invariant correlator. The procedure is the same as the one employed in Eq. (III):

[TABLE]

where $\Sigma(j)$ indicates the position of the site obtained by applying the symmetry $\Sigma$ to the site $j$. The quantum numbers of the above correlators, which are associated to the different point group symmetries, are all zero by construction. A more general strategy to implement symmetries with the desired quantum numbers is outlined in Ref. Choo et al., 2018.

The main advantage of the RBM correlator with respect to the Jastrow factor comes from the fact that its accuracy can be, in principle, systematically improved by increasing the number of hidden variables $N_{\alpha}$. This is due to the fact that a RBM is a universal function approximator Le Roux and Bengio (2008), which means that it can approximate any function with arbitrary accuracy if the number of variables in the hidden layer is allowed to grow arbitarily large. In addition, the non-local structure of this neural network makes it capable of capturing highly-entangled phases of matter. On the other hand, the disadvantage of employing the RBM correlator in relation to the simpler two-body Jastrow factor mainly resides in its higher computational cost. Crucial for an efficient variational Monte Carlo, computing ratios of translationally invariant RBMs has a cost which scales linearly with the number of sites and the number of hidden units ($O(N\times N_{\alpha})$), while computing ratios of two-body Jastrow factors simply scales as $O(1)$ Becca and Sorella (2017). Finally, another disadvantage of the RBM is the lack of straightforward physical interpretability of its variational parameters, which are associated to many-body spin-spin correlations at all distances. Instead, the pseudopotential $v_{i,j}$ of the Jastrow factor clearly accounts for the two-body correlation of spins $i$ and $j$, and typically shows a clear physical behavior Capello et al. (2005), decaying with the distance $|R_{i}-R_{j}|$.

IV Optimization of the variational wave function

In order to minimize the variational energy of the wave function (3), namely

[TABLE]

we employ the stochastic reconfiguration (SR) technique Sorella (2005); Becca and Sorella (2017), which is briefly summarized below.

Let us denote by $\omega=\{\omega_{k}\}$ the set of all the variational parameters of the RBM-fermionic state, which is formed by the couplings included in the auxiliary Hamiltonian $\mathcal{H}_{0}$ and by the weights and biases of the RBM correlator. For each parameter $\omega_{k}$, we can define a corresponding operator $\hat{\mathcal{O}}_{k}$ that is diagonal in the basis of spin configurations, i.e., ${\langle\sigma|\hat{\mathcal{O}}_{k}|\sigma^{\prime}\rangle=\mathcal{O}_{k}(\sigma)\delta_{\sigma,\sigma^{\prime}}}$, and yields the logarithmic derivative of the amplitudes of $|\Psi_{\mathcal{C}}\rangle$:

[TABLE]

Let us also introduce another diagonal operator, the so-called local energy $\hat{E}_{\rm loc}$, whose matrix elements read

[TABLE]

At each step of the SR algorithm, the parameters of the variational wave functions are updated according to

[TABLE]

Here $\eta$ is an arbitrary hyperparameter dubbed as the learning rate, $\mathcal{S}^{-1}$ is the inverse of the covariance matrix

[TABLE]

and $f_{k}$ are the forces Becca and Sorella (2017)

[TABLE]

Both the elements of the covariance matrix and the forces are computed stochastically by Monte Carlo sampling (through the Metropolis algorithm Metropolis et al. (1953)), and then employed to update the variational parameters as in Eq. (17). We note that a redundant parametrization of the variational wave function causes the covariance matrix to be non-invertible. This situation is commonly encountered when neural network correlators with a large number of variational parameters are considered. To avoid numerical instabilities connected to the singularity of the covariance matrix, we apply an explicit regularization of the form $\mathcal{S}_{k,k^{\prime}}\mapsto(1+\epsilon\delta_{k,k^{\prime}})\mathcal{S}_{k,k^{\prime}}$, where $\epsilon$ is an arbitary small parameter Becca and Sorella (2017).

For most of the variational results presented in this work, we performed few distinct optimizations starting from different set of initial parameters, and we selected as optimal wave function the one which provided the lowest variational energy. In general, we used $\approx 2\times 10^{4}$ up to $\approx 2\times 10^{5}$ Monte Carlo samples to evaluate the forces and the covariance matrix of the SR method. We typically chose a learning rate $\eta$ in the range $[0.02,0.05]$ and a regularization parameter $\epsilon=10^{-3}$. For what concerns the latter, in the case of variational wave functions with a large number of variational parameters (e.g. $N_{\alpha}=8$ on a $10\times 10$ lattice), we observed that starting the optimization with a larger value of $\epsilon$, e.g. $10^{-2}$, helps reducing some numerical instabilities and yields a smoother decrease of the variational energy; then, after few hundreds SR steps, the value of $\epsilon$ can be safely changed to $10^{-3}$ for the rest of the optimization.

V Results

V.1 The $J_{1}-J_{2}$ model on the square lattice

We first discuss the variational Monte Carlo results for the $J_{1}-J_{2}$ model on a $6\times 6$ square lattice. For $J_{2}=0$, the auxiliary fermionic Hamiltonian $\mathcal{H}_{0}$ which defines the variational Ansatz contains a Néel magnetic field (with $Q=(\pi,\pi)$) in Eq. (8) and a complex hopping term in Eq. (7), which induces a staggered flux through the square plaquettes Affleck and Marston (1988). This fermionic wave function possesses the sign structure of the exact ground state of the Heisenberg model, i.e., it follows the so-called Marshall-Peierls sign rule Marshall and Peierls (1955); Becca et al. (2011). Therefore, it is sufficient to use an RBM correlator with real parameters, so that $C_{\mathrm{RBM}}(\sigma)\geq 0$. On the contrary, in the frustrated regime, the sign structure of the exact ground state is unknown. For $J_{2}/J_{1}=0.5$ we combine a complex-valued RBM correlator and a spin liquid fermionic wave function. The latter is obtained by considering an auxiliary Hamiltonian with a uniform real hopping and a $d_{x^{2}-y^{2}}$ pairing at first neighbors, and a $d_{xy}$ pairing at fifth neighbors Hu et al. (2013). This wave function corresponds to a $Z_{2}$ spin liquid state, which satisfies all the symmetries of the model after the Gutzwiller projection Wen (2002). Once more, we point out that all the parameters in the variational state, i.e., the weights and biases of the RBM and the couplings included in $\mathcal{H}_{0}$, are fully optimized through the stochastic reconfiguration technique Sorella (2005); Becca and Sorella (2017).

One of the purposes of this work is evaluating the accuracy gain provided by the use of the RBM correlator. For magnetically ordered phases (e.g., $J_{2}=0$), the best fermionic state breaks the $SU(2)$ symmetry [since $\Delta_{\rm AF}$ is finite in Eq. (8)]; in this case, the reference state also contains the Jastrow factor, since it plays an important role to include quantum fluctuations. Instead, within the non-magnetic phase (e.g., $J_{2}/J_{1}=0.5$), the best fermionic state does not break the $SU(2)$ symmetry (i.e., $\Delta_{\rm AF}=0$) and the Jastrow factor typically gives a negligible contribution. In this case, we choose to compare the accuracy of the combined RBM-fermionic wave function to the pure Gutzwiller-projected BCS state. In both cases, the RBM correlator is expected to yield an improvement of the accuracy when all the parameters are properly optimized, since it can be shown that the Jastrow factor can be represented exactly as a RBM Nomura et al. (2017); Glasser et al. (2018). The underlying question we address is understanding to which extent the application of the RBM correlator improves the accuracy and the physical content of the variational wave function. For the $J_{1}-J_{2}$ model, we consider both the translationally invariant [$\hat{\mathcal{C}}_{\mathrm{tRBM}}$, Eq. (III)] and the fully symmetric [$\hat{\mathcal{C}}_{\mathrm{sRBM}}$, Eq. (III)] RBM correlators.

In Fig. 1 we report the relative error of the variational energy of the $J_{1}-J_{2}$ model on the $6\times 6$ lattice with respect to the exact value, obtained by Lanczos diagonalization. This quantity is defined as

[TABLE]

where $E_{\mathcal{C}}$ is the energy of a given variational Ansatz and $E_{0}$ is the exact ground-state energy. The results clearly show that the variational wave function is overall more accurate in the unfrustrated regime with respect to the frustrated one. In particular, at $J_{2}=0$ the inclusion of the RBM provides a large energy gain with respect to the fermionic reference state ($\Delta E\approx 0.4\%$): the relative error of the energy improves by a factor $\approx 8$ when applying the translationally invariant RBM ($\Delta E\approx 0.05\%$ for $N_{\alpha}=12$) and a factor of $\approx 20$ when applying the fully symmetric RBM ($\Delta E\approx 0.02\%$ for $N_{\alpha}=12$). In general, for each value of $N_{\alpha}$ we observe that the inclusion of the point group symmetries in the RBM correlator halves the relative error of the variational energy. By contrast, at $J_{2}/J_{1}=0.5$ the accuracy gain is considerably smaller: the relative error of the energy changes from $\Delta E\approx 0.6\%$ to $\Delta E\approx 0.3\%$ by applying a translationally invariant RBM correlator, and the addition of point group symmetries is far less effective than what we observe for the unfrustrated case. In Fig. 1 the variational energies obtained with the RBM-fermionic wave functions are compared to the ones of Ref. Choo et al., 2019, in which a convolutional neural network (CNN) quantum state is employed: the CNN wave function is more accurate in the unfrustrated phase, while the RBM-fermionic Ansatz gives better energies in the frustrated regime.

To further elucidate the ability of the RBM-fermionic wave functions in capturing the ground state properties of the model, we compute the spin-spin correlations at different distances on the lattice. Whereas the exact ground state wave function of the $J_{1}-J_{2}$ model on finite cluster possesses all the symmetries of its Hamiltonian, most of the variational Ansätze used in our calculations break the spin $SU(2)$ symmetry of the model due to the presence of the RBM correlator, the Jastrow factor, or the magnetic field $\Delta_{\rm AF}$. Therefore, to investigate the spin symmetry properties of the resulting states, we separate the computation of the spin-spin correlations at distance $R$ into the in-plane part,

[TABLE]

and the out-of-plane part,

[TABLE]

We then compare the variational estimates of the two contributions separately with the exact value ${\langle\mathcal{S}_{R}\rangle_{0}\equiv 1/3\langle\mathbf{S}_{0}\cdot\mathbf{S}_{R}\rangle_{0}=\langle S^{a}_{0}S^{a}_{R}\rangle_{0}}$, $a=x,y,z$. Here $\langle\dots\rangle_{0}$ indicates the expectation value over the exact ground state $|\Psi_{0}\rangle$, obtained by Lanczos diagonalization. Thus, the relative error of the spin-spin correlations at distance $R$ is computed as

[TABLE]

where $\langle\dots\rangle_{\mathcal{C}}$ indicates the expectation value over a given variational wave function, and $\alpha=xy,z$ refer to the in-plane, Eq. (21), and out-of plane, Eq. (22), estimators.

In Fig. 2, we show the in-plane and out-of-plane variational correlations as a function of distance, in comparison to the exact value. The results are obtained with the aforementioned RBM-fermionic wave function $\hat{\mathcal{C}}_{\mathrm{sRBM}}|\Psi_{f}\rangle$ ($N_{\alpha}=1$) for $J_{2}=0$. Due to the presence of the antiferromagnetic parameter $\Delta_{\rm AF}$, which induces magnetic ordering in the $xy$ plane, the in-plane correlations overestimate the exact ones in absolute value, while the out-of-plane correlations are underestimated. This tendency is observed for the spin-spin correlations at any distance and for any value of $N_{\alpha}$. The results obtained by increasing the number of hidden units are shown in Fig. 3, where the relative error of the variational estimates of the correlations with respect to the exact value is reported for some selected distances $R$. We observe a systematic improvement of the accuracy when the number of hidden units is increased, with both the in-plane and the out-of-plane correlations approaching the exact value. The fact that these two terms tend to get closer to each other when $N_{\alpha}$ is increased indicates that the RBM correlator tries to restore the anticipated spin $SU(2)$ symmetry in the wave function. As expected, since the optimal wave function is computed by minimizing the ground state energy, the most accurate values for the correlations are obtained at first-neighbors. We note also that, at any distance, the RBM correlator systematically provides a more accurate estimation of the spin-spin correlations than the simple Jastrow factor.

Before moving to the frustrated case, we explain the motivation of our choice to consider the accuracy of $\mathcal{S}^{xy}_{R}$ and $\mathcal{S}^{z}_{R}$ separately, instead of using the isotropic spin-spin correlation, $(\mathcal{S}^{z}_{R}+2\mathcal{S}^{xy}_{R})/3$. When performing the optimization of the variational parameters of the wave function, we observed that the final set of weights and biases of the RBM can depend on their initial values. Indeed, the parametrization of the RBM displays a considerable degree of redundancy, as indicated by the presence of zero eigenvalues of the covariance matrix of Eq. (18) Sorella (2005); Becca and Sorella (2017). Therefore, the optimization procedure can end into different local minima which can have equivalent energies and slightly different in-plane and out-of-plane spin-spin correlations. However, this small difference, which is typically negligible in the unfrustrated regime, can be enhanced or reduced in a random fashion by cancellation of error when the isotropic correlation is computed. Therefore, we argue that a study of the separate components of the correlation function, rather than the isotropic counterpart, provides a better characterization of the symmetry properties of our Ansatz.

At $J_{2}/J_{1}=0.5$, the situation is considerably different from the unfrustrated case, as demostrated in Fig. 4. The accuracy of the variational correlations does not show a systematic improvement with the number of hidden units. Indeed, even if in general the results are more accurate for $N_{\alpha}=12$ than $N_{\alpha}=1$, the behavior of the relative error is not as smooth as the respective one observed in Fig. 3 for $J_{2}=0$. Moreover, in some cases the relative error obtained by applying the RBM correlator is larger than the one obtained from the use of the simple Jastrow factor. A regular improvement is observed only at first and second neighbors, which are the correlations contributing to the value of the energy. Most importantly, the role of the RBM correlator regarding the expected $SU(2)$ symmetry is not clear; we find that the out-of-plane correlations display a better improvement than the in-plane ones when $N_{\alpha}$ is increased. We argue that this irregular behaviour of the accuracy of the spin-spin correlations is a consequence of the presence of several local minima in the optimization of the variational parameters of the RBM correlator, which lead to states with markedly different energies and correlation functions. While the effect is already present in the unfrustrated Heinseberg model, it is enhanced in the highly-frustrated regime of the $J_{1}-J_{2}$ model.

To investigate this issue, we performed $40$ distinct optimizations of the $\hat{\mathcal{C}}_{\mathrm{tRBM}}|\Psi_{f}\rangle$ Ansatz with $N_{\alpha}=4$ at $J_{2}/J_{1}=0.5$, choosing different values of the RBM parameters as starting point. The relative error of the spin-spin correlations obtained by the $40$ resulting wave functions is plotted as a function of the relative error of the variational energy in Fig. 5. We observe that the accuracy of the results show considerably large fluctuations. In particular, while the out-of-plane correlations seem to be more accurate when the variational energy is lower, an opposite effect is observed for the in-plane terms at some distances. The RBM correlator, which is a function of the $\sigma_{z}$ degrees of freedom, tends to “sacrifice” the accuracy of the in-plane correlations for the sake of improving the variational energy. This numerical experiment suggests that the complexity of the optimization landscape of RBM in the presence of frustration is significantly different from the unfrustrated case: while in the unfrustrated case all the local minima display similar energies and correlation functions, the highly frustrated regime exhibits a wide array of local minima with similar variational energy but strinkingly different correlation functions. We speculate that these minima are due to the presence of a possible glassy phase in the optimization induced by frustration Baity-Jesi et al. (2018).

We conclude our analysis of the $J_{1}-J_{2}$ model by presenting the variational energies for the $10\times 10$ lattice. The wave functions we employ have the same form of the ones used for the $6\times 6$ lattice. However, while the number of parameters entering $\mathcal{H}_{0}$ is the same, the number of weights of the RBM increases linearly with the size. For the Heinseberg model ($J_{2}=0$), in Fig. 6 we show the relative error of the variational energy with respect to the exact one, computed by quantum Monte Carlo Sandvik (1997); Calandra Buonaura and Sorella (1998). The relative error of the RBM-fermionic wave function is of the same order of magnitude of the one obtained on the $6\times 6$ lattice and shows a remarkable energy gain with respect to the Jastrow-fermionic state. However, the accuracy gain provided by the inclusion of the point group symmetries is slightly smaller than the one observed for the $6\times 6$ lattice. We note that the variational energy of the CNN quantum state of Ref. Choo et al., 2019 is lower than the best RBM-fermionic energy found in this work. Here, we emphasize the fact that the CNN state employs a larger number of variational parameters (3838 complex numbers Choo et al. (2019)) than the RBM-fermionic wave function with $N_{\alpha}=8$ (810 real numbers). However, the local structure of the CNN is advantageous in the process of optimization, since the optimal parameters obtained for a smaller lattice can be employed as a starting point for the optimization of the wave function on a larger lattice Choo et al. (2019). The same procedure cannot be applied in the case of the RBM correlator, due to its highly nonlocal structure, which implies that the optimization of the parameters of this network necessarily becomes harder when the size of the system increases.

In Fig. 7 we compare our variational energies in the frustrated phase, $J_{2}/J_{1}=0.5$, with several different results from literature. Here, at variance with the unfrustrated case, the variational energies obtained by using the RBM-fermionic wave function are better than the ones of the CNN of Ref. Choo et al., 2019, and are very close to the best density matrix renormalization group (DMRG) estimates of Ref. Gong et al., 2014. However, a considerably lower variational energy is obtained in Ref. Hu et al., 2013, where a fermionic wave function, defined by a BCS Hamiltonian $\mathcal{H}_{0}$ which contains two additional $d_{x^{2}-y^{2}}$ pairings (at fourth- and sixth-neighbors) with respect to the one employed in this work, is improved by the application of few Lanczos steps. The relative energy gain provided by two Lanczos steps, which require the addition of only two variational parameter, is remarkably larger than the improvement which is obtained by the application of the RBM correlator, which contains more than $2000$ parameters for $N_{\alpha}=12$. In general, as already observed for the $6\times 6$ lattice, in the frustrated regime the RBM correlator yields a much smaller energy gain with respect to the unfrustrated case.

In summary, these results suggest that the RBM correlator provides a systematic way of improving the description of magnetically ordered phases beyond the Jastrow factor, where the RBM effectivetively induces out-of-plane fluctuations that counterbalance the in-plane magnetic order induced by $\Delta_{\rm AF}$. In the frustrated regime, even though the application of the RBM leads to better variational energies, it generally does not improve the description of the correlation functions beyond nearest neighbors. Apart from their numerically expensive training procedure, the RBM’s energetic enhancement comes at the price of breaking of the $SU(2)$ symmetry of the fermionic wave function. The symmetry breaking we observe is especially evident in the correlation functions beyond nearest neighbors, which do not directly affect the variational energy during the optimization procedure. Fortunately, all these aspects can be addressed in the future using a recently introduced parametrization of the RBM correlator, which by construction satisfies the $SU(2)$ symmetry of the models considered in our study Vieijra et al. (2019).

V.2 The $XY$ model on the square lattice

As already pointed out, one of the main drawbacks of the introduction of the RBM correlator is the breaking of the spin $SU(2)$ symmetry of the wave function. For this reason, we evaluate the accuracy of the RBM-fermionic construction for a model whose Hamiltonian has lower symmetry, i.e., $XY$ model of Eq. (2). The exact ground state of the $XY$ model has the same sign structure of the one of the Heinseberg model, i.e. it follows the Marshall-Peierls rule Marshall and Peierls (1955). Thus, we employ an analogous RBM-fermionic wave function like the one used for the Heinseberg model (with real weights and biases). The variational results are reported in Fig. 8 for the translationally invariant correlator $\hat{\mathcal{C}}_{\mathrm{tRBM}}$. The relative error of the variational energy with respect to the exact one is at least a factor of $10$ smaller than the one obtained for the Heisenberg model and the accuracy gain provided by the RBM correlator is remarkable ($\Delta E\approx 0.0004\%$ for $N_{\alpha}=12$). The higher accuracy of the wave function is related to the fact that here the symmetry of the variational Ansatz is consistent with the spin symmetry of the model.

V.3 The Heisenberg model on the triangular lattice

Our previous results suggest the idea that the application of the RBM correlator is more effective for magnetically ordered phases, as exemplified in the Heisenberg and $XY$ model, rather than for non-magnetic ones, as in the frustrated region of the $J_{1}-J_{2}$ model. However, in the Néel phases considered above, the exact sign structure of the wave function is particularly simple and exactly captured by the fermionic part of the variational Ansatz. To try to disentangle whether the successes observed in our simulations are related to the special structure of the sign or to the presence of magnetic order, we consider a model whose ground state is magnetically ordered but displays a non-trivial sign structure: the Heisenberg model on the triangular lattice. In this case, the fermionic part of the Ansatz is constructed via an auxiliary Hamiltonian $\mathcal{H}_{0}$ which features a magnetic field $\Delta_{\mathrm{AF}}$ with pitch vector $Q=(\frac{4\pi}{3},0)$ [or, equivalently, $Q=(\frac{2\pi}{3},\frac{2\pi}{\sqrt{3}})$] and a real nearest-neighbor hopping $t$. The sign structure of the hopping generates an alternation of [math] and $\pi$ fluxes threading the triangular plaquettes Iqbal et al. (2016). While in square lattice a real parametrization on top of the fermionic state already gives an accurate representation for both signs and amplitudes of the exact ground state, we anticipate that the triangular lattice Heisenberg model requires a complex-valued correlator to approximate the unknown sign structure of the wave function induced by the geometric frustration of the problem. Thus, on top of the fermionic state, we apply a translationally invariant RBM correlator with complex parameters ($\hat{\mathcal{C}}_{\mathrm{tRBM}}$).

In Fig. 9, we compare the accuracy of the RBM-fermionic wave function for the Heisenberg model on the $6\times 6$ square and triangular lattices. The energy gain provided by the application of the RBM correlator is considerably larger in the case of the square lattice, where the relative error of the energy decreases of a factor $\approx 8$, with respect to the case of the triangular lattice, where it decreases of a factor $\approx 1.5$ (from $\Delta E\approx 2.2\%$ to $\Delta E\approx 1.4\%$). Overall, the variational energy is more accurate on the square lattice than on the triangular lattice, and the relative errors differ by an order of magnitude.

In Fig. 9, we also compare the results of the aforementioned RBM-fermionic wave functions to the ones obtained by simpler Ansätze, which are constructed by setting the hopping terms to zero and considering an auxiliary fermionic Hamiltonian with only magnetic field ($\mathcal{H}_{0}=\mathcal{H}_{\mathrm{AF}}$). In this way, the fermionic degrees of freedom are localized and $|\Psi_{f}\rangle$ reduces to a (projected) product state $|\Psi_{Q}\rangle=\mathcal{P}_{S^{z}_{tot}=0}\prod_{i=1}^{N}(|\uparrow\rangle_{i}+e^{iQR_{i}}|\downarrow\rangle_{i})$, which displays “classical” order in the $xy$-plane. This wave function can be employed as a reference state for the application of a Jastrow factor Carrasquilla et al. (2013) or the RBM correlator. Here, we apply a translationally invariant RBM correlator to $|\Psi_{Q}\rangle$. We observe that the energy gain provided by the presence of the hopping term in the fermionic Ansatz is remarkably large when the simple Jastrow factor is applied to the reference state. However, when the RBM correlator is employed, the contribution of the hopping term becomes less important, decreasing considerably with the number of hidden units, which suggests that the RBM replaces the effect of the fermionic hopping term in the state.

Finally, in Fig. 10 we present the variational energies obtained for the Heisenberg model on the $12\times 12$ triangular lattice. For a large enough $N_{\alpha}$, the variational energy of the RBM-fermionic wave function is more accurate than the Green’s function Monte Carlo (GFMC) result of Ref. Capriotti et al., 1999. To summarize our numerical experiments, we surmise that the high accuracy of the results for the unfrustrated square lattice Heisenberg model is due to having an exact representation of the sign structure built in the Ansatz, which in turn, alleviates the energy optimization problem. On the other hand, the lower accuracy in the triangular lattice model is presumably due to a combination of the approximate nature of the sign structure imposed by our Ansatz and the energy optimization problem in the presence of complex numbers. Since magnetically ordered states are pervasive in frustrated magnetism, it remains an important issue to unequivocally establish whether the origin of the high accuracy of the results for the Heisenberg model on the square lattice is only due to the absence of frustration or to the fact that the sign structure is exactly known in the unfrustrated case.

VI Conclusion

Inspired by Hinton’s product of experts idea Hinton (2002), we have studied a neural augmentation of the parton construction that combines a family of physically motivated Gutzwiller-projected variational states with a complex-valued RBM correlator. We focused our attention on various prototypical spin models traditionally used in condensed matter physics which exhibit a wide array of conventional ground states with magnetic order, as well as more exotic ones where a spin-liquid behavior has been anticipated.

In agreement with previous results based upon neural networks alone Carleo and Troyer (2017); Choo et al. (2019), our calculations showed that RBMs are very effective in the unfrustrated Heisenberg and $XY$ models on the square lattice, where the knowledge of the exact sign structure of the ground state allows us to use a real parametrization of the RBM. Here, a few hidden units in the neural network are sufficient to reach a striking accuracy, which for the Heisenberg case is comparable with the best variational wave functions defined within the bosonic resonating valence-bond picture Liang et al. (1988). Moreover, we emphasize the remarkable ability of the RBM to systematically recover the spin $SU(2)$ symmetry upon increasing the number of hidden units. This is clear not only by looking at first-neighbor correlations, but also at further distances: both out-of-plane correlations (i.e., along $z$) and in-plane ones (i.e., in the $xy$ plane) exhibit a clear tendendcy to converge toward the exact (isotropic) result upon increasing $N_{\alpha}$.

In the highly frustrated regime, the exact ground-state sign is not known a priori, even when the system displays magnetic order, like in the Heisenberg model on the triangular lattice. Here, we find a substantial energy gain with respect to the original parton wave function, although the variational procedure does not yield the same accuray as in the square lattice, even for a relatively large number of hidden units. A similar effect is also observed in the highly-frustated regime of the $J_{1}-J_{2}$ Heisenberg model, where for $J_{2}/J_{1}=0.5$ our neural Gutzwiller-projected wave functions exhibit accuracies beyond recent neural network calculations based on CNNs Choo et al. (2019), though with a lower accuracy than in the $J_{2}=0$ case. Based on the numerical experiments presented in Fig. 5, we surmise that the application of neural network variational states to frustrated quantum spin systems requires an extensive investigation of the intricate relation between the representation power of neural networks to capture highly-entangled states of matter with a complicated sign structure and the complexity of the optimization landscape of the problem induced by frustration. Motivated by Ref. Day et al., 2019, the disentangling of these factors could be approached through a clustering analysis of the trained RBM parameters and their associated spin-spin correlation functions, which may shed light onto the interplay between frustration, sign structure and entanglement Grover and Fisher (2015), and the rough optimization landscape of the problem Baity-Jesi et al. (2018).

At variance with the unfrustrated limit, the $SU(2)$ symmetry is broken and hardly recovered when increasing $N_{\alpha}$ in the frustrated case: the out-of-plane correlations are clearly more accurate than in-plane ones, highligthing the “asymmetry” of the RBM form, which is defined in terms of the $z$-component of the local spin. In this regard, our calculations provide a first investigation into the inductive bias of the RBM applied to ground states of frustrated systems, where we observe a strong tendency to break $SU(2)$ symmetry in our setting. Thus, the recent proposal to generalize the RBM to fulfill the $SU(2)$ symmetry Vieijra et al. (2019) may pave the way for studying fully symmetric neural Gutzwiller-projected wave functions, which will be especially relevant in the understanding of models for spin liquid phases and other states exhibiting fractionalized excitations and gauge structures.

Acknowledgements.

We would like to thank K. Choo, G. Carleo, G. Torlai, and R. Melko for useful discussions. F.F. acknowledges the kind hospitality and financial support of the Vector Institute for Artificial Intelligence in Toronto, where this project started. We also thank the Kavli Institute for Theoretical Physics (KITP) in Santa Barbara and the program “Machine Learning for Quantum Many-Body Physics”. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. J.C. acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Canada CIFAR AI chair program.