Neural Network Quantum States analysis of the Shastry-Sutherland model
We utilize neural network quantum states (NQS) to investigate the ground
state properties of the Heisenberg model on a Shastry-Sutherland lattice using
the variational Monte Carlo method. We show that already relatively simple NQSs
can be used to approximate the ground state of this model in its different
phases and regimes. We first compare several types of NQSs with each other on
small lattices and benchmark their variational energies against the exact
diagonalization results. We argue that when precision, generality, and
computational costs are taken into account, a good choice for addressing larger
systems is a shallow restricted Boltzmann machine NQS. We then show that such
NQS can describe the main phases of the model in zero magnetic field. Moreover,
NQS based on a restricted Boltzmann machine correctly describes the intriguing
plateaus forming in magnetization of the model as a function of increasing
magnetic field.
Matěj Mezera, Jana Menšíková, Pavel Baláž, Martin Žonda
**Neural Network Quantum States analysis of the Shastry-Sutherland model**

## Abstract

**We utilize neural-network quantum states (NQSs) to investigate the ground-state properties of the Heisenberg model on a Shastry-Sutherland lattice via variational Monte Carlo method. We show that already relatively simple NQSs can be used to approximate the ground state of this model in its different phases and regimes. We first compare several types of NQSs with each other on small lattices and benchmark their variational energies against the exact diagonalization results. We argue that when precision, generality and computational costs are taken into account, a good choice for addressing larger systems is a shallow restricted Boltzmann machine NQS. We then show that such NQS can describe the main phases of the model in zero magnetic field. Moreover, NQS based on a restricted Boltzmann machine correctly describes the intriguing plateaus forming in magnetization of the model as a function of increasing magnetic field.**

###### Contents

* 1 Introduction
* 2 Shastry-Sutherland model
	* 2.1 Basic properties of the ground state
* 3 Methods
	* 3.1 Variational Monte Carlo and machine learning
	* 3.2 Neural network quantum states
* 4 Results
	* 4.1 Comparison of different NQSs architectures
* 5
* 4.2 Investigation of the ground-state phase diagrams
		* 4.2.1 Ground-state orderings
		* 4.2.2 Zero magnetic field
		* 4.2.3 Magnetization plateaus
* 5 Conclusion
* A Lattice tiles
* B Visible biases in sRBM and pRBM
* C Symmetries

## 1 Introduction

The neural network quantum states (NQSs) have recently emerged as a promising alternative to common trial states in variational Monte Carlo (VMC) studies of many-body lattice problems [1, 2, 3, 4, 5, 6, 7, 8, 9]. This research is driven by the fact that neural networks (NNs) are universal function approximators [10] as well as by the astonishing progress in the field of machine learning (ML) in general. These advancements already lead to a number of effective ML applications suitable for basic research of quantum systems and technologies [11, 12, 13, 14]. For example, even simple NQSs, such as the restricted Boltzmann machine (RBM), allow us to investigate the ground-state properties of various quantum spin models. It was already shown that RBM can outperform standard trial states in the variational search of the ground-state energies of antiferromagnetic Heisenberg model [1]. Very promising results have been obtained also for frustrated spin systems such as \(J_{1}-J_{2}\) model [15, 16, 17, 18]. Here NQSs can be trained to capture the nontrivial sign structure of the ground state and in some cases have even achieved the state-of-the-art accuracy [19] delivering cutting edge results. Nevertheless, two-dimensional frustrated quantum spin models continue to be a challenge for NQSs as well as for other methods [20]. For example, it is not clear yet how to choose an optimal neural network architecture for a particular frustrated system, how important is the role of the trial state symmetries in the learning process, or if a NQS with favorable variational energy also encodes a physically correct state.

Not all of these issues are specific to NQSs. Results of any VMC calculations are to a large extent dictated by the properties and limitations of the used trial states. An inappropriately chosen variational state, i.e., one with small overlap with the ground state, can still give a good estimate of the ground-state energy [21]. If some additional information is known about the ground state, e.g., its symmetries, one can pick a more restrictive variational state function. However, this is often not an optimal strategy if the goal is to seek for new phases or to locate a phase boundary. In principle, NQSs could be a remedy for such problems. It is reasonable to expect that a single, but expressive enough, NQS can be used to approximate distinct phases. This assumption is supported by results of Sharir et al. [22] who showed that NQSs can have even higher expressive power than matrix product states [23] and projected entangled pair states [24] as these can be efficiently
mapped to a subset of NQSs. In other words, NQSs can be effectively utilized to a larger class of quantum states than these powerful formalisms which are known primarily from their usage in Density Matrix Renormalization Group (DMRG) but are also utilized as variational states in VMC [1, 21, 25].

In practice it is not yet clear how to achieve this in a general case. Despite tremendous progress the research of frustrated quantum spin magnets is still in the stage of testing and developing NQS architectures for simple models often focusing primarily on reaching the best variational energy in particular regimes [6, 15, 16, 18]. In the present work we aim at a different target. We want to demonstrate that even shallow NQSs can be sufficient for the investigation of qualitatively different ground-state orderings including states forming only in finite magnetic field. To this goal we focus on the ground state of antiferromagnetic Heisenberg Hamiltonian on Shastry-Sutherland lattice known as the Shastry-Sutherland model (SSM) which we introduce in more detail in Sec. 2. To our knowledge, this model of frustrated quantum spin system has not been previously addressed within the NQS context, yet it seems to be an ideal testbed for our purposes.

SSM was already investigated by a number of methods including exact diagonalization (ED) techniques [26, 27, 28, 29], Quantum Monte Carlo [30], various versions of DMRG [31, 32, 33, 34, 35], perturbation theory [36] and even with quantum annealing [37]. These studies shown that SSM has a rich ground-state phase diagram. In zero magnetic field these include regions such as singlet spin dimer phase, antiferromagnetic Neel state, spin plaquette singlet phase and probably other phases. Introducing finite magnetic field further complicates the picture. Consequently, it is challenging to find a single variational function that can correctly approximate the whole ground-state phase diagram.

In addition, there are still open questions related to the ground-state phase diagram in zero as well as in the finite magnetic field even in some experimentally relevant regimes of the model. This is important because several magnetic materials have a structure topologically equivalent to SSM. The most notable examples are \(\text{SrCu}_{2}(\text{BO}_{3})_{2}\), \(\text{BaNd}_{2}\text{ZnO}_{5}\) and rare-earth tetraborides \(\text{RB}_{4}\) (R=Dy, Er, Tm, Tb, Ho) [38, 39, 40, 41, 42]. All exhibit an intriguing step-like dependence of the overall magnetization on the external magnetic field which has been found to be inherent to SSM [43, 44]. Here each plateau reflects a stable nontrivial spin ordering. The magnetic behavior of these materials is still not completely understood. This together with other open problems, e.g., the prospect of a narrow spin liquid phase in the zero magnetic field, further motivates the investigation of SSM and its generalizations [45, 46, 36, 47].

Therefore, SSM presents a model system which has the right combination of properties that are well understood and can be used to benchmark various NQSs, and of open problems that can be potentially illuminated by these variational techniques. This includes a possibility to address a rather complex behavior of a system in relation to a changing magnetic field.

The present work consists of two main parts. In the first one we explore SSM by employing a number of NQS architectures and we test them against ED results for small lattices in zero magnetic field. Here the primarily goal is to find one or few networks that are able to capture the main well understood ground-state orderings of SSM. Simultaneously, we require that these NQSs have a high chance to describe the magnetization plateaus as well. This means that the ideal network has to give a solid approximation of the ground-state orderings even when no conditions on the total magnetization are imposed. Consequently, we do not focus on getting the best possible variational energy for a particular set of parameters. Rather, we require a good approximation of the energy in distinct regimes of the model, correct description of the particular orderings, and reasonable computation complexity that allows the usage of the NQS on larger lattices. We argue that when precision, generality and computational costs are taken into account, a shallow RBM
with complex parameters is still a good choice.

In the second part we introduce a refined learning protocol for RBM NQS and test it for a broad range of model parameters and different network sizes. We then utilize it in the study of larger systems. We first investigate the zero magnetic field scenario and demonstrate that RBM is expressive enough to capture all main phases of the system. We then move to the model in finite magnetic field and show that, with the right learning strategy, RBM is able to capture the magnetization plateaus crucial for the description of real materials. This opens a possibility that NQSs could be used to investigate several open problems, such as the existence of still opaque spin liquid phase and other orderings predicted but not yet confirmed in SSM.

## 2 Shastry-Sutherland model

SSM is described by the Hamiltonian

\[\hat{H}=J\sum_{\langle i,j\rangle}\hat{\mathbf{\xi}}_{i}\cdot\hat{\mathbf{\xi}}_{j}+J^{ \prime}\sum_{\langle i,j\rangle^{\prime}}\hat{\mathbf{\xi}}_{i}\cdot\hat{\mathbf{\xi}}_ {j}-h\sum_{i}\hat{\mathbf{\xi}}_{i}^{z}, \tag{1}\]

where \(\hat{\mathbf{\xi}}_{i}=\frac{1}{2}\hat{\mathbf{\sigma}}_{i}\) is the spin-\(1/2\) operator at the \(i\)-th site with \(\hat{\mathbf{\sigma}}\) being the vector of Pauli matrices. The first term represents exchange coupling between the nearest neighbors on a square lattice (solid lines in Fig. 1). The second term is a sum over specific diagonal bonds arranged in a checkerboard pattern (dashed lines in Fig. 1). Note that these sums are interpreted in terms of nodes, i.e., there is no double-counting.

Both coupling constants are antiferromagnetic (\(J,J^{\prime}>0\)) and we set \(J^{\prime}\) as the unit of energy in the whole paper. The last term describes the influence of the external magnetic field \(h\) pointing to the \(z\)-direction.

### Basic properties of the ground state

The basic structure of the SSM ground-state phase diagram is well understood. As illustrated in Fig. 2 the SSM at \(h=0\) has at least three distinct ground-state orderings. These are the _dimer singlet_ (DS) state for (\(J^{\prime}\gg J\)), the _Neel antiferromagnetic_ (AF) ordering (\(J^{\prime}\ll J\)) and the _plaquette singlet_ (PS) state in between. The phase transition from DS to PS state is of the first-order [28]. The nature of the transition from PS to AF is still in debate. The ED study of Nakano and Sakai [29] suggests that the supposed PS phase actually consists of at least two distinct phases. In addition, some

Figure 1: (a) The Shastry-Sutherland lattice. Bonds with coupling strength \(J\) are represented by solid lines while bonds with \(J^{\prime}\) by dashed ones.
recent studies argue that there is a so-called _deconfined quantum critical point_ (DQCP), separating a line of first-order transitions and, potentially, a narrow gapless _spin liquid_ (SL) phase [34, 35, 48].

Nevertheless, even without focusing on the possible DQCP and SL phase, the three main orderings, namely DS, PS, and AF, already pose a sufficient challenge for a single variational state because of their distinctive character and symmetries.

The _DS phase_ is formed by an exactly (analytically) accessible state [49]. It was verified by numerous analytical and numerical methods that it remains to be the ground state up to \(J/J^{\prime}\approx 0.675\)[28, 29, 34]. In the limiting case of \(J\ll J^{\prime}\), the system is equivalent to an ensemble of independent spin dimers each forming a singlet ground state. The DS ground state is thus a direct product of dimer singlet states

\[\ket{\psi}_{\text{DS}}=\bigotimes_{\langle i,j\rangle^{\prime}}\frac{1}{\sqrt{ 2}}\left(\ket{\uparrow}_{i,j}-\ket{\downarrow\uparrow}_{i,j}\right)\,. \tag{2}\]

As such, it is antisymmetric with respect to the exchange of two intradimer spins and symmetric with respect to transformations rearranging only the spin pairs without swapping the intradimer spins. The energy of the dimer ground state is

\[E_{\text{DS}}=-\frac{3}{8}J^{\prime}N\,, \tag{3}\]

where \(N\) is the number of lattice sites (twice the number of dimers).

The _PS phase_ can be understood as weakly coupled plaquette singlet states illustrated in Fig. 2. The plaquette singlet is a ground state of an isolated 4-spin Heisenberg cluster with four bonds arranged in a cycle [28]. The pattern of the plaquette singlets in Fig. 2 indicates that the PS state is two-fold degenerate.

Figure 2: Illustration of the SSM phase diagram of the SSM for small \(h\) based on the results from Ref. [35]. There is a first-order transition at \(J/J^{\prime}\approx 0.675\) between DS and PS phases. The gray squares in PS depict the plaquette singlets. The nature of the transition between the PS and AF phases remains unresolved. It is not clear whether there is a narrow spin liquid phase, a DQCP or just a second order transition in the region labeled with a question mark.
It is important to stress again that the relevant range \(J/J^{\prime}\) discussed here (\(0.675\lesssim J/J^{\prime}\)\(\lesssim 0.82\)) might be much more complex. As we already stated, it has been argued that at \(J/J^{\prime}\approx 0.70\) the PS phase splits into two distinct regions with quantitatively different behaviors [29, 34, 35, 48]. For the sake of simplicity we omit this possibility in most of our discussion. Nevertheless, this might be important for future more detailed studies.

The _AF phase_ stabilizes when \(J/J^{\prime}\gtrsim 0.82\). When \(J^{\prime}\) becomes negligible, the ground state of SSM is approaching the ground state of the antiferromagnetic Heisenberg model with only nearest neighbor bonds on a square lattice. Although this state is not analytically accessible, it can be explored by Monte Carlo (MC) simulations [26]. Using the first-order correction to these quantum MC results, the energy of the SSM in the AF phase was estimated [26] to be

\[E_{\text{AF}}=(0.102J^{\prime}-0.669J)N \tag{4}\]

where \(N\) is assumed to be large.

A more detailed discussion of the symmetries of these three states is postponed to the Appendix. C. Note that the three main phases DS, PS and AF are reasonably understood and simultaneously they differ qualitatively. This is one of several qualities of the model that make SSM a suitable testbed for NQSs.

So far we have discussed the \(h=0\) case. When we introduce finite magnetic field to the DS phase (Eq. (2)) some dimers can morph into triplet states. These triplets are formed in repeating patterns, e.g., checkerboard, stripes, or more complex configurations (for illustration see Fig. 3), giving rise to stable plateaus of constant magnetization in increasing magnetic field.

Because each plateau signals a distinct stable ordering it also presents a challenge for the NQSs. Particularly so, because a finite magnetic field does not allow for a simple restriction of the Hilbert space to its zero magnetization part. This restriction was heavily utilized in previous NQS investigations of quantum spin models. Note, that it is mostly these plateaus which make SSM interesting experimentally. Good examples are SrCu\({}_{2}\)(BO\({}_{3}\))\({}_{2}\), BaNd\({}_{2}\)ZnO\({}_{5}\), CaCo\({}_{2}\)Al\({}_{8}\) and rare-earth tetraborides RB\({}_{4}\) (R=Dy, Er, Tm, Tb, Ho)) [38, 39, 40, 41, 42] which all exhibit the intriguing step-like dependence of the overall magnetization on the external magnetic field or show magnetic frustration and can be modeled by SSM or its generalizations.

Figure 3: A simplified illustration of magnetization as a function of external magnetic field \(h\) and coupling constant \(J\) inspired by Ref. [50]. A more detailed illustration would contain additional steps (e.g., supersolid phase), however, their actual position and width is not clear yet. The arrangements of singlets and triplets are displayed for some of the plateaus (namely the \(m^{z}=1,1/2,1/3\) and \(1/4\) plateaus are shown).
## 3 Methods

### Variational Monte Carlo and machine learning

VMC is a standard method that allows us to stochastically evaluate the expectation values of quantum operators without the need to probe the full Hilbert space. In short, suppose we have a Hamiltonian operator \(\hat{H}\) and a trial wave function \(|\psi_{\mathbf{\theta}}\rangle\) depending continuously on a vector of parameters \(\mathbf{\theta}\). We are seeking a ground state of \(\hat{H}\) or its approximation in a variational way. The goal is to minimize the variational energy

\[E_{\mathbf{\theta}}=\langle\hat{H}\rangle_{\mathbf{\theta}}:=\frac{\langle\psi_{\mathbf{ \theta}}|\hat{H}|\psi_{\mathbf{\theta}}\rangle}{\langle\psi_{\mathbf{\theta}}|\psi_{ \mathbf{\theta}}\rangle}\geq E_{0} \tag{5}\]

with respect to the vector of parameters \(\mathbf{\theta}\), where \(E_{0}\) is the true ground-state energy. We utilize a fixed orthonormal basis \(\{|\mathbf{\sigma}^{z}\rangle\}\) of the \(z\)-projected \(\frac{1}{2}\)-spins and use the notation

\[|\psi_{\mathbf{\theta}}\rangle=\sum_{\mathbf{\sigma}^{z}}\psi_{\mathbf{\theta}}(\mathbf{ \sigma}^{z})\,|\mathbf{\sigma}^{z}\rangle\,,\text{ where }\,\langle\mathbf{\sigma}^{z}|\psi_{\mathbf{ \theta}}\rangle\equiv\psi_{\mathbf{\theta}}(\mathbf{\sigma}^{z}) \tag{6}\]

as typical in NQS studies [51]. The variational energy in Eq. (5) is, in the jargon of ML, a _loss function_. Using this loss function the parameters \(\mathbf{\theta}\) are optimized to obtain the lowest-energy state that the chosen variational function can represent. In practice, we use the VMC implementation from the NQS toolbox NetKet [9, 51].

In general, the form of the trial wave function \(\psi_{\mathbf{\theta}}(\mathbf{\sigma}^{z})\) restricts the optimization process to a subset of the Hilbert space. An improper choice of the ansatz can bias the approximation towards a wrong phase, or even can make approaching the correct state impossible. Clearly, this is where one can expect that NQSs could outperform standard variational states due to their high expressiveness.

### Neural network quantum states

In the first part of our work we explore several NQS architectures [6, 9]. We chose these particular networks due to their successful application in previous studies of other Heisenberg models.

_Restricted Boltzmann machine (RBM)_ is a generative artificial NN constituted of a visible layer with \(N\) nodes (one for each lattice site) fully connected with a single hidden layer with \(M=\alpha N\) nodes (hidden degrees of freedom) where \(\alpha\) is the _hidden layer density_[1]. It can be used to define a NQS

\[\log\psi_{\mathbf{\theta}}(\mathbf{\sigma}^{z})=\sum_{i}\sigma_{i}^{z}a_{i}+\sum_{j} \log\left[2\cosh(\sum_{i}W_{ij}\sigma_{i}^{z}+b_{j})\right], \tag{7}\]

where vector \(\mathbf{\theta}\) contains the variation network parameters \(\mathbf{\theta}=\{\mathbf{a},\mathbf{b},\mathbf{W}\}\). This NQS can be interpreted as a one-layered fully-connected neural network with \(\log\cosh\) activation function followed by summation of outputs and additional summation of visible biases [1]. Note that complex-valued parameters are necessary in order to represent generally complex-valued wave function outputs.

The size of the visible layer \(N\) is fixed by the size of the investigated spin system. The expressive power of RBM can be, however, modified by changing \(\alpha\). According to the universal approximation theorem [52], the RBM is theoretically able to express any wave function to any desired degree of
accuracy when \(\alpha\) can be arbitrarily large. In practice, we aim for a reasonably small \(\alpha\) to restrict the total number of variational parameters of RBM, which is \(MN+M+N=\mathcal{O}(\alpha N^{2})\).

_Modulus-phase split real-valued RBM (rRBM):_ One can avoid complex parameters, which in general make the learning process harder, by introducing two independent real-valued NNs [17, 53] to represent the modulus \(A(\mathbf{\sigma}^{z})\) and the phase \(\Phi(\mathbf{\sigma}^{z})\) of the wave function separately

\[\log\psi_{\mathbf{\theta}}(\mathbf{\sigma}^{z})=A(\mathbf{\sigma}^{z})+\mathrm{i}\Phi(\bm {\sigma}^{z}). \tag{8}\]

Unlike in the Ref. [53] where rRBM architecture proved to be advantageous in the investigation of transverse-field Ising model, we have experienced that for SSM, the rRBM shows worse results than complex-valued RBM. This is in accord with the recent study of other frustrated system, namely \(J_{1}-J_{2}\) model [18]. Consequently, we discuss the results of this network only briefly in chapter 4.1 and focus predominately on complex valued architectures.

_Symmetric variant of RBM (sRBM):_ Carleo and Troyer [1] used translational symmetries to reduce the number of variational parameters in RBM. They replaced the fully-connected layer with a convolutional layer and set the visible biases to the constant value \(a^{f}\) across each convolutional filter \(f\). The resulting expression for its output is

\[\log\psi_{\mathbf{\theta}}(\mathbf{\sigma}^{z})=\sum_{f=1}^{F}\sum_{g\in G}\left\{a^{ f}\underbrace{\sum_{i=1}^{N}T_{g}(\mathbf{\sigma}^{z})_{i}}_{\sum\limits_{i}^{N} \sigma_{i}^{z}=m^{z}}\right. \tag{9}\]

Here \(\mathbf{T}_{g}\) denotes a symmetry transformation of a spin configuration according to an element \(g\) from the symmetry group \(G\) of order \(|G|\). Index \(f\) denotes different feature filters. The number of these filters \(F\) determines the network size \(M=F|G|\). Resulting sRBM has fewer variational parameters than the RBM by a factor of \(|G|\). We can view this approach as binding the values of some of the \(\mathcal{O}(\alpha N^{2})\) parameter making the total asymptotic number of parameters \(\mathcal{O}(\alpha N)\). Carleo and Troyer [1] also showed that this approach significantly improves convergence and accuracy of the ground states of the antiferromagnetic Heisenberg model on a square lattice. However, this approach suffers from two crucial disadvantages in more general circumstances. First drawback is that visible biases are inherently constant for each filter \(f\) which significantly lowers the expressiveness of the network as discussed later in this section. As we show in Appendix B, the sRBM architecture cannot be modified to ease this condition while preserving symmetries. The second drawback is that sRBM is not applicable if the ground state does not transform under the trivial irreducible representation (irrep) of a given symmetry group.

To illustrate the problem let us consider a single spin dimer (i.e., a single bond of SSM with \(J=0\), \(J^{\prime}=1\) and \(h=0\)). Its ground state is a singlet \(|\psi_{0}\rangle=(|\uparrow\downarrow\rangle-|\uparrow\downarrow\rangle)/ \sqrt{2}\). The symmetry group of the single dimer Hamiltonian contains just two operations - an identity and a swap of both spins \(G=\{g_{12},g_{21}\}\). If we apply the swap operation to the ground state, we obtain \(\hat{\mathbf{T}}_{g_{21}}\ket{\psi_{0}}=(|\downarrow\uparrow\rangle-|\downarrow \uparrow\rangle)/\sqrt{2}=-|\psi_{0}\rangle\). Although this state is a multiple of the ground state, we see that it does not transform under the trivial irrep because one of its characters is \(\chi_{g_{21}}=-1\). Because sRBM represents only states where \(\hat{\mathbf{T}}_{g}\ket{\psi}=|\psi\rangle\); \(\forall g\in G\) this symmetry should not be used in sRBM. Note, that we do not strictly follow this rule and sometimes use all available lattice symmetries. The reason is that this leads to NQS with a small number of parameters that is easy to optimize. The resulting variational energy can be then compared with the energy obtained with RBM with the same \(\alpha\) to check how well is the full network optimized, i.e., if it leads to lower
energy than sRBM. If not, this signals that the variational energy of RBM can be lowered by a better learning.

_Projected RMB (pRBM)_: Recently, Nomura [54] introduced an alternative way to symmetrize RBM (or any other NN) by using quantum-number projection (also called _incomplete symmetrization operator_)

\[\psi_{\mathbf{\theta}}^{G}(\mathbf{\sigma}^{z})=\sum_{g\in G}\chi_{g^{-1}} \psi_{\theta}(\mathbf{T}_{g}(\mathbf{\sigma}^{z}))\,, \tag{10}\]

where \(g\) is an element of the given symmetry group \(G\) and \(\chi_{g}\) is its character from the irrep in question. The wave function on the right-hand side may be arbitrary and the function on the left-hand side satisfies the desired transformation property \(\psi_{\mathbf{\theta}}^{G}(\mathbf{T}_{g}(\mathbf{\sigma}^{z}))=\chi_{g}\psi_{\mathbf{\theta}}^ {G}(\mathbf{\sigma}^{z})\) in case of one-dimensional representation or \(\psi_{\mathbf{\theta},a}^{G}(\mathbf{T}_{g}(\mathbf{\sigma}^{z}))=\sum_{b=1}^{d}D(g)_{b}^{ a}\psi_{\mathbf{\theta},b}^{G}(\mathbf{\sigma}^{z})\) for more-dimensional irreps where functions \(\psi_{\mathbf{\theta},a}^{G}\) form a basis of \(d\)-dimensional irrep \(D(g)_{b}^{a}\). Unfortunately, pRBM makes the learning process of NN much more expensive than sRBM. The computational time increases by a factor of \(|G|\) producing a computational cost \(\mathcal{O}(\alpha N^{2}|G|)\). On the other hand, pRBM implementation does not suffer from the problems mentioned for sRBM and it can be generalized by setting mutually independent visible biases (see Appendix B).

_Group-convolutional NN (GCNN)_: The group equivariant convolutional NNs represent a promising class of NNs built inherently on symmetries. They were proposed by Cohen and Ni [55] as a natural extension of the well-known convolutional neural networks. While convolutional networks preserve invariance under translations, GCNN are equivariant under the action of an arbitrary group \(G\) (which may contain a subgroup of translations). Roth and MacDonald [56] further improved GCNNs so that they can transform under an arbitrary irreducible representation of \(G\), which is more suitable for NQSs for SSM. GCNN can be composed of any number of hidden layers. The first and subsequent layers are given by

\[\mathbf{f}_{g}^{1}=\mathbf{f}\left(\sum_{i=1}^{N}W_{g^{-1}i}^{0}\sigma_{ i}^{z}+\mathbf{b}^{0}\right),\ \ \mathbf{f}_{g}^{k+1}=\mathbf{f}\left(\sum_{h\in G}W_{g^{-1}h}^{k}\mathbf{f}_{h}^{k}+\mathbf{b}^ {k}\right), \tag{11}\]

where \(\mathbf{f}\) is a nonlinear activation function (the output is typically a vector since GCNN can have multiple parallel feature filters) and \(\mathbf{f}_{g}^{1}\) is a 1st-layer feature vector corresponding to group element \(g\). The result of the last layer \(\mathbf{f}_{g}^{K}=\mathbf{f}_{g}^{(J)K}\), where \((j)\) denotes the individual features of the layer, is then projected in a similar fashion as pRBM

\[\psi(\mathbf{\sigma}^{z})=\sum_{g\in G}\sum_{j}\chi_{g^{-1}}\exp\Bigl{(} \mathbf{f}_{g}^{(j)K}\Bigr{)}\,. \tag{12}\]

The main advantage over symmetrizing an arbitrary deep network by formula (10) is that we do not need to evaluate the forward pass of the non-symmetric wavefunction \(|G|\) times. This is achieved thanks to the fact that each layer of the GCNN fulfills _equivariance_. GCNN with \(K\) layers and a typical number of feature filters \(F\) in each layer has \(\mathcal{O}(FN+KF^{2}|G|)\) parameters.

_Jastrow network:_ As a baseline, we also use a Jastrow network based on the standard Jastrow ansatz [57, 58]

\[\psi_{\mathbf{\theta}}=\exp\left(\sum_{i,j}\mathbf{\sigma}_{i}^{z}W_{i,j} \mathbf{\sigma}_{j}^{z}\right)\,, \tag{13}\]
where the variational parameters \(\mathbf{\theta}=\left\{W_{i,j}\right\}\) form a matrix of size \(N\times N\). The Jastrow ansatz is physically motivated by two-body interactions and assigns trainable parameters \(W_{i,j}\) to pairwise spin correlations. The number of its parameters scales as \(\mathcal{O}(N^{2})\).

The complicated sign structure of the complex phases of the basis coefficients forming the ground-state wavefunction presents a major challenge in optimizing parameters of a variational function of a frustrated spin system. In case of Heisenberg model on a bipartite lattice consisting of sublattices A and B (i.e., SSM with \(J^{\prime}=0\)) this can be solved using the Marshal sign rule (MSR) [59]. The MSR states that the sign of \(\psi(\sigma^{z})\) is given by \((-1)^{N_{A}^{\dagger}(\sigma^{z})}\) where \(N_{A}^{\dagger}(\sigma^{z})\) is the total number of up-spins on a sublattice A. Because this alternates with a spin-flip, it can be difficult for NN to learn the correct signs. However, it is possible to circumvent this problem in two analogous ways.

If the sign structure is dictated by MSR, the Hamiltonian can be gauge transformed by changing the signs of some terms to make all wave-function coefficients positive in the transformed basis. In particular, we change \(\sigma^{x}\rightarrow-\sigma^{x}\) and \(\sigma^{y}\rightarrow-\sigma^{y}\); for \(\forall\sigma\in\text{A}\). The same result can be also obtained by setting the visible biases to \(a_{i}=i\pi/2\) for \(i\in A\) and \(a_{i}=0\) for \(i\in B\) as this exactly reconstructs the Marshall sign factor (up to an overall constant factor). In another words, the biases can be set to play the role of Marshall basis. What is important here is that in general case the simple Marshall sign rule is not always applicable. Especially problematic are systems with strong frustration [17, 18, 60]. The advantage of using the visible biases instead, is that their setting does not have to be known ahead as it can be, despite possible technical difficulties, learned. It is, therefore, beneficial to include the visible biases whenever allowed by the architecture. An additional bonus is that free visible biases also allow to overcome an improper initialization of weights.

## 4 Results

### Comparison of different NQSs architectures

It is not feasible to apply all NQSs introduced above to investigate the ground-state phase diagram of SSM at large lattices. Therefore, in the first part of our investigation, we benchmark these NQSs against the exact results obtained by Lanczos ED method. The aim is to identify a network that is both expressive enough to cover various phases and computationally tractable even for large lattices. We focus on a regular lattice with \(N=4\times 4=16\) points and irregular lattice with \(N=20\) (see Appendix A). Throughout this paper, we apply periodic boundary conditions for all lattices used. The irregular \(N=20\) is considered because \(N=16\) lattice has some undesirable properties, e.g., some extra symmetries with trivial irrep which favor symmetric networks. It also suffers by stronger finite size effects and does not exhibit the PS phase. On the other hand it is regular and easy to calculate.

We initially focus on just the two cases represented by \(J/J^{\prime}=0.2\) (DF phase) and \(J/J^{\prime}=0.9\) (AF phase) and investigate the model with and without MSR. Because the goal here is to compare different networks we estimate the accuracy of each architecture by comparing the average energy of the last 50 learning iterations \(E_{50}\) with the exact result \(E_{\text{ex}}\). The same computational protocol is used for each architecture. In particular, we use 2000 MC samples and 1000 training iterations for three values of fixed learning rates \((0.2,0.05,0.01)\). Each particular combination of architecture and basis (MSR or direct) was computed four times for each learning rate (yielding 12 independent
runs for each case of interest). This is to eliminate occasional events when NN gets stuck in a local energy minimum too far from the ground state. The zero magnetization was not implicitly assumed (i.e., we use local single-spin-flip Metropolis updates in VMC). We summarize our results in Table 1 where the values are \(\min\limits_{i}\frac{\left|E_{50}^{i}-E_{\mathrm{ex}}\right|}{E_{\mathrm{ex}}}\), with \(i\) enumerating the twelve independent runs.

There are several results in Table 1 which were important for our decision on which network should be used in the detailed study of the phase diagram at larger lattices. Starting with RBM one can see that networks with \(\alpha=2\) (560 parameters for \(N=16\) and 860 for \(N=20\)) and 4 (3380 parameters for \(N=20\)) and 8 (4398 parameters for \(N=16\)) show similar precision where the significantly larger networks are notably better (approximately three times) only in the AF phase. For general case this favors the computationally less demanding network with \(\alpha=2\). Interesting is also the comparison with the Jastrow network. Both architectures have comparable precision in the DS phase for \(N=16\), however, RBM in the AF phase as well as in the DS phase for \(N=20\)

\begin{table}
\begin{tabular}{l r|r r|r r} \hline \multicolumn{2}{c|}{\(N=4\times 4\)} & \multicolumn{2}{c|}{\(J/J^{\prime}=0.2\) (DS)} & \multicolumn{2}{c}{\(J/J^{\prime}=0.9\) (AF)} \\ architecture & params & direct & MSR & direct & MSR \\ \hline Jastrow & 256 & 5.8\(\times 10^{-5}\) & 1.1\(\times 10^{-5}\) & 6.2\(\times 10^{-2}\) & 2.1\(\times 10^{-2}\) \\ RBM (\(\alpha=2\)) & 560 & 1.9\(\times 10^{-5}\) & 1.6\(\times 10^{-5}\) & 4.7\(\times 10^{-3}\) & 5.1\(\times 10^{-3}\) \\ RBM (\(\alpha=16\)) & 4368 & 2.8\(\times 10^{-5}\) & 1.4\(\times 10^{-5}\) & 1.6\(\times 10^{-3}\) & 1.1\(\times 10^{-3}\) \\ rRBM (\(\alpha=2\)) & 1088 & 2.3\(\times 10^{-4}\) & 2.1\(\times 10^{-4}\) & 8.3\(\times 10^{-3}\) & 7.9\(\times 10^{-3}\) \\ rRBM (\(\alpha=8\)) & 4352 & 2.3\(\times 10^{-4}\) & 2.2\(\times 10^{-4}\) & 9.3\(\times 10^{-3}\) & 7.8\(\times 10^{-3}\) \\ sRBM (\(\alpha=4\)) & 18 & 0.0 & 0.0 & 4.8\(\times 10^{-3}\) & 2.5\(\times 10^{-3}\) \\ sRBM (\(\alpha=16\)) & 69 & 0.0 & 0.0 & 8.9\(\times 10^{-4}\) & 3.8\(\times 10^{-4}\) \\ sRBM (\(\alpha=128\)) & 545 & 0.0 & 0.0 & 1.3\(\times 10^{-3}\) & 8.5\(\times 10^{-5}\) \\ pRBM (\(\alpha=0.5\)) & 136 & 7.9\(\times 10^{-5}\) & 1.2\(\times 10^{-4}\) & 1.5\(\times 10^{-3}\) & 8.5\(\times 10^{-4}\) \\ pRBM (\(\alpha=2\)) & 544 & 9.1\(\times 10^{-6}\) & 2.2\(\times 10^{-5}\) & 7.1\(\times 10^{-5}\) & 2.0\(\times 10^{-5}\) \\ GCNN & 2188 & 6.2\(\times 10^{-6}\) & 3.6\(\times 10^{-6}\) & 3.2\(\times 10^{-5}\) & 3.6\(\times 10^{-5}\) \\ GCNNt & 268 & 4.2\(\times 10^{-7}\) & 5.1\(\times 10^{-7}\) & 4.9\(\times 10^{-3}\) & 4.7\(\times 10^{-3}\) \\ \hline \multicolumn{2}{c|}{\(N=20\)} & \multicolumn{2}{c|}{\(J/J^{\prime}=0.2\) (DS)} & \multicolumn{2}{c}{\(J/J^{\prime}=0.9\) (AF)} \\ \hline Jastrow & 400 & 1.1\(\times 10^{-5}\) & 1.0\(\times 10^{-3}\) & 1.4\(\times 10^{-1}\) & 3.0\(\times 10^{-2}\) \\ RBM (\(\alpha=2\)) & 860 & 2.2\(\times 10^{-5}\) & 1.4\(\times 10^{-5}\) & 6.6\(\times 10^{-3}\) & 6.2\(\times 10^{-3}\) \\ RBM (\(\alpha=4\)) & 3380 & 1.2\(\times 10^{-5}\) & 1.7\(\times 10^{-5}\) & 2.2\(\times 10^{-3}\) & 2.1\(\times 10^{-3}\) \\ sRBM (\(\alpha=4\)) & 85 & 1.2\(\times 10^{-1}\) & 1.5\(\times 10^{-1}\) & 5.0\(\times 10^{-2}\) & 1.4\(\times 10^{-3}\) \\ pRBM (for AF phase) & 336 & 2.3\(\times 10^{-1}\) & 2.3\(\times 10^{-1}\) & 3.5\(\times 10^{-3}\) & 3.0\(\times 10^{-3}\) \\ pRBM (for DS phase) & 336 & 7.1\(\times 10^{-4}\) & 6.8\(\times 10^{-5}\) & 4.4\(\times 10^{-2}\) & 4.7\(\times 10^{-2}\) \\ \hline \end{tabular}
\end{table}
Table 1: Comparison of the precision of NQS variational results on lattices \(N=16\) and \(N=20\). The listed values were calculated as \(\min\limits_{i}\frac{\left|E_{50}^{i}-E_{\mathrm{ex}}\right|}{E_{\mathrm{ex}}}\), where \(E_{50}^{i}\) is the average energy of the last 50 iterations of the \(i\)-th run. A number of variational parameters is also shown for each architecture. The difference between GCNN and GCNNt is that for GCNN we used the correct characters for the expected ground state, GCNNt utilized only the translation symmetry. Error 0.0 here means relative error less than \(10^{-7}\) which we consider as “numerical precision” due to the standard MC errors which are typically larger even for \(L=16\).
with MSR is by one or even two orders of magnitude more precise than the Jastrow ansatz.

For \(N=16\) even better seems to be the sRBM architecture. In its implementation we have utilized full automorphism group of the finite lattice. Despite the resulting small number of variational parameters it shows excellent precision. Actually, significantly increasing \(\alpha\) is not that beneficial (compare \(\alpha=4\) and \(\alpha=128\) case). In the DS phase the usage of the symmetries allowed sRBM with \(\alpha=4\) to find the ground-state energies within the numerical precision (hence the zero error). Because sRBM is just a restricted RBM, this already suggests that the learning protocol for RBM can be improved which we demonstrate in the next section. However, it is important to stress that the excellent results are a consequence of the special symmetries of the \(N=16\) lattice. Both DS and AF states transform under the trivial irreducible representation and automorphism group is therefore applicable. This is not true for DS ground state in different tiles, including regular ones such as \(N=6\times 6\) (for a more detailed discussion of the symmetries see Appendix C). This is illustrated in the second part of Table 1 where sRBM with \(\alpha=4\) gives very bad results in the DS phase of \(N=20\) due to the improper symmetries. In short, using inappropriate symmetries in sRBM for states that do not transform under trivial irrep can make the variational energy significantly worse than for simple RBM. For \(N=20\) sRBM also fails in the AF phase but only when the direct basis is used. This implies that sRBM has troubles to learn the correct sign structure of the state for larger lattices, which can be attributed to the fixed visible biases.

Both remaining architectures, namely pRBM and GCNN, show excellent accuracy for \(N=16\). They clearly outperform all other networks in the AF phase. However, results at \(N=20\) are less convincing especially when one takes into account that these networks are more computationally demanding than RBM even for cases when RBM contains more parameters. In addition, the reached precision required a usage of correct symmetries of the expected state, i.e., the proper line form Table 2 in Appendix C. If one uses an improper one, i.e., if different state is expected, as illustrated by the last two lines in Table 1, the precision can drop by several orders of magnitude. Similarly, the precision decreases significantly for both GCNN and pRBM when we use only the group of translations instead of the full symmetry group, as illustrated by GCNN in Table 1. Note that for this case the precision in the AF phase drops to the level of a simple RBM with \(\alpha=2\). The network is much better in the DS phase, but as we will discuss in the next chapter, already RBM with \(\alpha=2\) and modified learning protocol can reach the numerical precision in this phase as well. Although we can not exclude that much better results could be obtained for the symmetrized pRBM and GCNN networks with a different learning protocol, considering their much higher computational demands and the necessity to a priori identify the correct irrep symmetries for each lattice type to make the learning efficient the presented results favor RBM for the investigation of larger clusters.

The last question that must be addressed here is if it is worth to use MSR. The Table 1 shows several cases where MSR is favorable in the AF phase (e.g., for sRBM and \(N=16\)), but this is not a general rule. In addition, its usage comes with a price as well. We have noticed that MSR basis seems to strongly favor the AF ordering even for \(J/J^{\prime}\) where PS is already the ground state in exact results. We will discuss this briefly when addressing larger lattices.

To wrap it up, in general, the usage of MSR basis does not lead to significantly better results. With some exceptions, here presented networks are able to approximate the ground-state energy quite well even without MSR. Therefore, we will mostly omit the MSR from further discussion. Furthermore, if the symmetry of the ground state is known, it is worth using this information in building the NN. If not, then the usage of just translations does not lead to a significant improvement of the precision. Fortunately, the complex valued RBM with visible biases can give a very good approximation of the ground-state energy without any restrictions. Its clear advantage is that
no preliminary information about the ground-state properties is needed. As such, it is suitable for problems where the character of the ground state or position of the phase boundary is unknown. In addition, the precision of RBM for SSM can be significantly improved using a different learning strategy discussed in the following section.

### Investigation of the ground-state phase diagrams

Focusing only on RBM allowed us to test several learning strategies and to use more precise MC calculations. What follows is a description of the best learning protocol that we have found, which we used to produce all results discussed bellow. It proved to be beneficial to use a more precise MC calculations already during the training. We typically generate 4000 - 12 000 MC samples at every sampling step. It was also more advantageous to run 10 - 30 independent learnings (with random initial variational parameters) with shorter learning times than to use few runs with a lot of learning iterations. We used approximately 2000 convergence iterations in each run. During the learning we have been lowering the learning rate \(\eta\) by several discrete steps. Typically we started with \(\eta\) = 0.08 (\(\approx\)200 iterations), then changed it to \(\eta\) = 0.04 (\(\approx\)1600 iterations), followed by \(\eta\) = 0.02 (\(\approx\)100 iterations), \(\eta\) = 0.01 (\(\approx\)100 iterations ) and \(\eta\) = 0.001 (\(\approx\)50 iterations). The trained RBM was then used to calculate expectation values of energy and order parameters, introduced in the next section, where we used 12 - 60 thousand states. Consequently, the Monte Carlo error bars in all presented figures are for small lattices negligible. The relevant absolute error comes from the learning process or limitations of the used NQS. The state with the lowest energy (evaluated more precisely after training) from all independent runs was kept as the final result in the following discussion. Due to the stochastic fluctuations in the learned parameters it was for some cases advantageous to refine the results by again rerunning the final state multiple times with high number of MC samples but small number (5-10) of iterations and small learning rate (\(\eta\)\(\leq\) 0.001) keeping again the result with the lowest energy. Also, as discussed bellow, we have utilized transfer learning in some problematic regimes.

#### 4.2.1 Ground-state orderings

As already discussed, a good agreement of the variational energy with the exact one does not guarantee that the variational state correctly captures the character of the exact ground state, i.e., that it reflects the correct phase. To examine this, and with the aim to see if RMB NQS can correctly capture the transitions between the phases, we calculate the order parameters for the three main expected orderings. They are constructed to be large (close to one) whenever the state is in the respective phase and small in other domains.

In particular, we define the order parameter for the DS phase as

\[\mathcal{P}_{\text{DS}}=-\frac{4}{3N\!\!\sum_{\langle i,j\rangle^{\prime}}^{ L}\,\left\langle\hat{\mathbf{s}}_{i}\cdot\hat{\mathbf{s}}_{j}\right\rangle, \tag{14}\]

which reflects the fact that operator \(\hat{\mathbf{s}}_{1}\cdot\hat{\mathbf{s}}_{2}\) has for isolated dimer the expectation value \(-\frac{3}{4}\hbar^{2}\) (singlet state). Therefore, \(\mathcal{P}_{\text{DS}}\) is one in the DS phase and strictly lower in other phases.

For the PS order parameter we use a definition based on order parameter from Ref. [35]

\[\mathcal{P}_{\text{PS}}=\frac{4}{N}\left|\left\langle\sum_{\mathbf{r}\in\text{ signlet}}\hat{Q}_{\mathbf{r}}-\sum_{\mathbf{r}\,\text{empty}}\hat{Q}_{\mathbf{r}}\right\rangle \right|, \tag{15}\]
where the order parameter is given by the difference \(\hat{Q}_{r}=\frac{1}{2}\left(\hat{P}_{r}+\hat{P}_{r}^{-1}\right)\), with \(\hat{P}_{r}\) being the permutation operator. This operator performs a cyclic permutation of four spins on a plaquette (a square on the lattice without the diagonal \(J^{\prime}\) bond) at position \(r\). The first sum in Eq. 15 runs over all singlet squares (see Fig. 0(b)) and the second sum is over all empty squares. The meaning of this construction can be understood looking at the the Fig. 0(b). Operator \(\hat{Q}_{r}\) gives a large mean value in the plaquette singlet (gray square) and a value close to zero in the empty square in between four plaquette singlets. In practice, we don't know which set of squares will become singlets, as the state is degenerated, therefore, we use the absolute value.

For the AF phase we employ the standard structure factor

\[\mathcal{P}_{\text{AF}}=\frac{1}{N^{2}}\sum_{ij}\text{e}^{\text{i}\mathbf{q}\cdot \mathbf{r}_{ij}}\left(\hat{\mathbf{S}}_{i}\cdot\hat{\mathbf{S}}_{j}\right), \tag{16}\]

where \(\mathbf{r}_{ij}\) denotes the difference in discrete coordinates of spin \(i\) and \(j\), and we take \(\mathbf{q}=(\pi,\pi)\) which measures the antiferromagnetic checkerboard ordering. Finally, in the case of finite magnetic field we use the normalized magnetization in the \(z\)-direction

\[\mathcal{M}=\frac{2}{N}\sum_{i}\left\langle\hat{\mathbf{S}}_{i}^{z}\right\rangle, \tag{17}\]

to identify the expected plateaus in the magnetization. These expectation values are calculated using VMC for trained RBM NQS.

#### 4.2.2 Zero magnetic field

We first investigate the phases of SSM in zero magnetic field. Unlike in procedure used to compare different network architectures, here we restrict the Hilbert space by the condition \(\mathcal{M}=0\). Before moving to larger lattices we test the RBM for \(N=20\) in a wide range of \(J/J^{\prime}\). We use the irregular lattice \(N=20\) because it shows an onset of the PS ordering (see black dashed line in Fig. 4(a)) not present for smaller regular lattices. We also readdress the role of the parameter \(\alpha\) within the new learning protocol, but start our discussion with the \(\alpha=2\) case.

As it is clear from the comparison of the ground-state energies in panels Fig. 4(b) and Fig. 4(c) the RBM variational energy is in a very good agreement with ED. The updated learning protocol ensures that the relative error in the \(J/J^{\prime}<0.68\) region, i.e., for the DS phase is on the order of the numerical precision already for \(\alpha=2\) despite not using any symmetries except the condition \(\mathcal{M}=0\). The largest error is in the vicinity of the expected first-order phase transition from DS to PS phases, but only from the side of expected PS phase. Nevertheless, even here the largest observed relative error in energy was approximately 1% for \(\alpha=2\).

Considering the focus of our study, even more important than the error in energy is the nature of optimized variational states. Panel (a) in Fig. 4 shows that already a shallow network, i.e., RBM with complex parameters and \(\alpha=2\) is expressive enough to correctly capture the formation of the distinct DS (blue diamonds) and AF ordering (red crosses), as well as the onset of the PS phase (black circles). The agreement is far from perfect though. In agreement with the results for the energy, the largest differences in order parameter values between RBM and ED are in the right vicinity of the expected phase transition. Here an error of 1% and less in the estimation of the ground-state energy translates in an error of tens of percents in the order parameters. Still, even here the RBM gives a correct qualitative picture. The position of the abrupt change of phase matches the exact result and there is a clear onset of the PS ordering. With increasing \(J/J^{\prime}\) the RBM results again align with the exact ones.
This benchmark shows that RBM with \(\alpha=2\) can easily capture the correct state in the DS phase, but it gives worse results above the critical \(J/J^{\prime}\approx 0.68\). What is not clear, is if the relative errors in panel (c) represent some inherent limitation of the RBM with small \(\alpha\), e.g., a difficulty to set the correct sign structure of the frustrated state, or are related to the learning process. Gradually increasing \(\alpha\) from 2 (blue circles) to 4 (red pluses), 8 (green pluses) and 16 (black diamonds with yellow cores) in the problematic region lowers the relative error in energy. However, this

Figure 4: Comparison of exact (lines) and various RBM variational results (symbols) at irregular lattice \(N=20\). (a) Evolution of the order parameters. Here blue solid line (ED), pure blue diamonds (RBM with \(\alpha=2\)) and blue diamonds with red edge (RBM with \(\alpha=16\)) show the DS order parameter; black dashed line (ED), black circles (RBM with \(\alpha=2\)) and black circles with yellow edge (RBM with \(\alpha=16\)) show the PS order parameter; and red dot-dashed line (ED), red crosses (RBM with \(\alpha=2\)) and red crosses with blue edge (RBM with \(\alpha=16\)) show the AF order parameter. Results of symmetric variants of RBM are not shown, as they were comparable to presented results for \(J/J^{\prime}>0.68\) and way off the exact results for \(J/J^{\prime}\leq 0.68\). (b) The exact (red line) and RBM \(\alpha=2,16\) ground-state energies. (c) Relative error in the ground-state energy for of the RBM with \(\alpha=2\) (blue circles), 4 (red pluses), 8 (green crosses) and 16 (black-yellow diamonds). Note that the relative error in the DS phase for RBM \(\alpha=2\) is at the level of numerical precision.
significant improvement in energy leads only to a small improvement for the order parameters near the critical point. This is shown in panel (a) where results calculated with RBM with \(\alpha=16\) are marked by the same symbols as for \(\alpha=2\) but highlighted via differently colored edges.

Using a symmetric NQS symmetries did not significantly improve the results. We have tested the sRBM architecture with \(\alpha=4\) in direct as well as MSR basis using the same protocol as for RBM. The sRBM results have been comparable to RBM for \(J/J^{\prime}>0.68\) and much worse than RBM results bellow this critical value. This suggests that the issue is not entirely due to the insufficient learning. On the other hand, the learning was the most difficult in the vicinity of the observed discontinuity. A significant fraction (often more than a half) of the independent runs for \(0.69\leq J/J^{\prime}\leq 0.72\) ended either in the wrong (DS) phase or even in a state with an energy much higher than the real ground state. This was not true for the rest of the \(J/J^{\prime}\) interval where most of the independent runs with the same \(\alpha\) showed very similar variational energies. In addition, the relative errors for all investigated RBM variants (including ones not presented here) follow the same pattern. They are maximal just above the critical point and then, if we neglect some noise, monotonously decrease with increasing \(J/J^{\prime}\). Yet, increasing \(\alpha\) significantly lowers the variational energy even for \(J/J^{\prime}>0.74\). This again suggests that the problem is indeed the small \(\alpha\). Ultimately, both statement seem to be correct. A significantly larger \(\alpha\) than \(\alpha=16\) is needed to capture the critical region together with a high precision learning, i.e., many independent runs.

After testing the RBM on small lattices and understanding its strength and limitations we can now approach larger ones. We focus on \(\alpha=2\) as the increased precision of variational energy obtained with larger \(\alpha\)'s does not significantly improve the estimates of the order parameters. Although we can not easily compare the VMC results with exact diagonalization for larger lattices, we can use the exact asymptotic results for the energy in DS Eq. 3 and AF phase Eq. 4 to guide us.

Fig. 5 shows the evolution of the order parameters and energy for \(N=20\), \(36\), \(64\) as well as several points calculated for \(N=100\). The results for \(N=20\), \(36\) and \(64\) are in very good agreement with the exact result in the assumed DS phase and are in between the \(N=20\) exact energy and the large \(N\) asymptotic in the supposed \(AF\) phase (up to one point discussed later). The advantage of the SSM model is that we can see, that \(N=100\) results are not converged enough even in the DS phase, which is much easier to access, with relative error of \(\approx 2\%\). Note, that this does not signal problems with the chosen NQS as it mostly reflects our limited computational resources. Because of the size of the lattice, we used less iterations in the learning process (\(\approx\)800) and only several independent runs, which proved to be enough only for a qualitative study confirming the general picture discussed below.

Leaving \(N=100\) aside Fig. 5 illustrates the usability of RBM for larger clusters. Although a much more thorough finite-size analysis would be necessary to asses the phase boundaries, the presented results confirm the overall picture of DS and AF phase separated by a narrow PS or at least its indication. There is, however, at least one issue. The point of the discontinuous phase transition from DS to PS phase should be \(J/J^{\prime}\simeq 0.675\), but our results at larger lattices push it to \(J/J^{\prime}\simeq 0.7\). Besides finite size effects, this could be also related to two technical problems. First is the already discussed difficulty to train the NQS in the vicinity of the discontinuous phase transition. The second is the tendency of the direct base to prefer DS over AF ordering. Both these issues can be seen in the panel (e) (inset of panel (d)) with the detail of the \(N=64\) results. Here green diamonds show the RBM data, red empty diamonds are sRBM data with MSR basis and blue stars are the sRBM data for direct basis all with \(\alpha=2\). Clearly, all these networks show (different) problems around the expected point of the phase-transition. For \(J/J^{\prime}=0.7\) and \(0.72\) sRBM with MSR gives lower energy than RBM and even lower than the energy of the DS ordering. The sharp
transition must be therefore placed bellow \(J/J^{\prime}=0.7\). On the other hand, sRBM with MSR can not correctly capture the onset of the DS ordering. The sRBM network with direct basis illustrates the opposite problem. It overestimates the stability of the DS ordering.

Investigation of sRBM showed that the RBM results at \(J/J^{\prime}=0.7\) are not fully converged yet. Because we have not been able to solve this problem using the direct approach, we utilized transfer learning. We used the RBM parameters trained for \(J/J^{\prime}=0.74\) as a starting point to train the network at \(J/J^{\prime}=0.72\), then used these results as a starting point for \(J/J^{\prime}=0.70\) and finally these results for \(0.69\). That way we got smaller variational energies for \(J/J^{\prime}=0.72\) and \(0.70\) than in the direct approach or the sRBM results and the \(J/J^{\prime}=0.72\) result dropped even bellow the DS energy. Interestingly, this also lead to an observable change in the order parameters (black crosses in all panels). In contrast to \(N=20\) case, especially sensitive is the PS ordering as seen from the comparison of black crosses and green diamonds in panel (d). Even if if is suggested by the order parameters, the transfer learning technique have not actually pushed the point of the expected phase transition bellow \(J/J^{\prime}=0.69\). The reason is, that the obtained energy is at this point larger than DS energy already reproduced by the direct approach. This shows, that although useful the transfer learning has to be used with care.

Figure 5: Evolution of order parameters for DS (a), PS (b), AF (c) and variational energy (d) as a function of \(J/J^{\prime}\) for \(h=0\) and various lattice sizes. All results in panels (a)-(d) have been obtained using RBM NQS with \(\alpha=2\) and VMC with exchange updates (simultaneous flip of two opposite spins in the basis state) for the Hilbert subspace restricted to \(\mathcal{M}=0\). The black dashed lines in panels (d) and (e) show the asymptotic energies for DS (horizontal) and AF phase (tilted). The black crosses represent \(N=64\) results for which we have utilized transfer learning. The inset (e) shows the detail of the variational energy for \(N=64\) in the vicinity of the phase transition calculated using RBM (green diamonds), sRBM in direct base (blue stars), sRBM with MSR (red empty diamonds) and three points calculated with RBM utilizing transfer learning (black crosses).
Taking into account these difficulties as well as the fact that our RBM results underestimate the PS order parameter even for \(N=20\) and large \(\alpha\), the investigation of possible SL phase and related DQCP, predicted in this problematic region, is currently beyond our reach. Nevertheless, here demonstrated expressiveness of a simple RBM with \(\alpha=2\) suggests that the problem can be indeed attacked by a more expressive or specialized networks. A good candidate might be a composed GCNN that would combine networks for different characters of the symmetry group for particular lattice size and boundary conditions.

#### 4.2.3 Magnetization plateaus

Historically, the most intriguing property of the SSM is its ability to describe fractional plateaus in magnetization as a function of external magnetic field which are observed also in real materials. To address this problem via VMC one has to drop the restriction of fixed \(\mathcal{M}=0\). In addition to significantly enlarging the Hilbert space, this also makes the optimization (learning) process a harder task. Moreover, each plateau represents a different ordering and, therefore, a challenge to NQS. Yet, as we demonstrate here already a simple RBM NQS with \(\alpha=2\) is sufficiently expressive to capture the main plateaus.

We focus on the case \(J/J^{\prime}=0.45\) which is inside of the DS phase (at \(h=0\)) where several broad plateaus are expected to form. The most stable ones, if allowed by the lattice size, should be the \(\mathcal{M}=1/2\) and \(1/3\) plateaus [27, 36]. We start the discussion by benchmarking the RBM NQS results (blue filled diamonds in all panels of Fig. 6) against the ED results for \(N=20\) lattice (blue solid lines). Clearly, the variational energy in panel (d) is in a very good agreement with the exact one. The relative error plotted in panel (c) is way bellow 1% in the whole range of \(h\). In addition, it shows a structure which can be understood by comparing the profile of the relative

Figure 6: Comparison of ED (blue solid lines) and RBM with \(\alpha=2\) (symbols) results for \(N=20\) and \(J/J^{\prime}=0.45\). Panels (a) and (b) show the magnetization and dimer state order parameter as functions of magnetic field. Panel (c) presents the relative error of the variational energy with respect to the ED result where blue dotted lines are just a guide to the eyes. Panel (d) shows the evolution of the normalized energy on external magnetic field. Blue filled diamonds represent the direct approach, empty red diamonds were obtained by utilizing the transfer learning discussed in the main text and the empty red squares by fixing \(\mathcal{M}N\) to integer values from the vicinity of the direct approach.
error dependence on \(h\) with the normalized magnetization plotted in panel (a) and the DS order parameter in panel (b). Panel (a) shows that RBM NQS with \(\alpha=2\) is able to capture all main steps in the magnetization observed in the ED curve. The most stable are \(\mathcal{M}=0\), \(1/2\) and \(1\) followed by plateaus \(1/5\) and \(3/10\) forming in the range \(0.7\lesssim h/J^{\prime}\lesssim 1.2\).

The stability of these plateaus is also reflected in the relative error. Despite not using any restriction on \(\mathcal{M}\), the relative error for \(h/J^{\prime}<0.7\), where \(\mathcal{M}=0\), is negligible. In this region the system stays in the DS ordering as revealed by panel (b). Similar situation is for \(h/J^{\prime}\geq 2.1\). Here the state is fully polarized (\(\mathcal{M}=1\)) and, therefore, easy to reproduce with variational techniques. Another regions with very small errors in the variational energy are the central parts of the above discussed stable plateaus as best illustrated by the \(1/2\) one. Here RBM NQS gives relative error below \(0.1\%\). Consequently, the regions with the largest errors are related to the transitions between the stable plateaus. Here we also observe the largest deviations of the NQS magnetization (and \(\mathcal{P}_{\text{DS}}\)) from the ED results. These problematic regions can be divided into two types. Into the first one belong the step edges, i.e., the abrupt changes of the magnetization for \(\mathcal{M}\leq 1/2\). The related convergence problems are similar to the difficulties to correctly capture the precise position of the discontinuous phase transition discussed for \(h=0\) and \(J/J^{\prime}\approx 0.69\). As such they can be also treated by the transfer learning. The red hollow diamonds in Fig. 6 were obtained by approaching the step edges from left and right using the RBM parameters learned in the centers of the neighboring plateaus as the initial input. Transfer learning clearly suppresses the errors and gives the correct value of \(\mathcal{M}\) even very close to the discontinuities.

The second problematic region is at large magnetic field where the \(\mathcal{M}=1/2\) plateau transits into \(\mathcal{M}=1\). Here, only one additional sudden step up of \(\mathcal{M}\) from \(1/2\) is expected in the thermodynamic limit, which is followed by continuous increase of \(\mathcal{M}\) to \(1\) as \(h\) raises. Finite \(N=20\) lattice shows in this region a number of very narrow transient steps. This makes this region unsuited for transfer learning, unless a much more refined grid of \(h\)'s is applied. On the other hand, the small lattice allowed us to test the actual expressiveness of RBM by fixing \(\mathcal{M}N\) to integer values taken from the vicinity on the direct RBM results for \(\mathcal{M}N\). The results with the lowest energies are depicted by the red empty squares and they reproduce both \(\mathcal{M}\) a DS of the exact study. This proves, that with correct learning strategy, RBM with \(\alpha=2\) is sufficient for capturing this rather complex evolution of the SSM ground state in the increasing magnetic field.

The stability of the magnetization plateaus must be confirmed on large lattices, because the magnetization could be always discrete on finite clusters yet continuous in thermodynamic limit. Moreover, lattice \(N=20\) is not divisible by three, thus can not hold the important \(1/3\) plateau. To show that RBM NQS can really capture these features we address larger clusters. Fig. 7 presents, in addition to the exact (solid red line) and RBM (red diamonds) results for \(N=20\), RBM results for \(N=36\) (blue squares) and \(N=64\) (yellow triangles). We stress here, that these results were obtained with the direct approach. We have not used the transfer learning and fixed \(\mathcal{M}\) to avoid the possibility, that we that way introduce a bias towards seemingly stable plateaus. Still, the results for \(N=36\) show stable flat steps in the magnetization which holds both the \(1/2\) and the \(1/3\) plateaus. Although the results for \(N=64\) are less stable, they confirm the \(1/2\) plateau and clearly signal formation of additional two plateaus for \(h/J^{\prime}<1.2\). These are very encouraging results as they again show that already a simple RBM with small number of parameters is expressive enough to correctly capture the complicated magnetization dependence reflecting the underlying complex ordering of the quantum spins.
## 5 Conclusion

We have investigated the ground-state properties of the Shastry-Sutherland model via variational Monte Carlo with NQS variational functions. Our main goal was to show that a single and relatively simple NQS architecture can be used to approximate a broad range of regimes of this model. We have first tested and benchmarked several NQS architectures known from the literature to be suitable for different variations of the Heisenberg model. We discussed the role, advantages and disadvantages of the NQSs incorporating lattice symmetries and biases on the visible layer. We conclude that when precision, generality and computational costs are taken into account, a good choice for addressing larger SSM lattices without as well as with external magnetic field is a restricted Boltzmann machine NQS with complex parameters.

Focusing on RBM NQS allowed us to refine the learning strategy, where we realized that if a more precise MC sampling is used then it is advantageous to run several (tens) short independent optimizations instead of few long learnings. Using this strategy, we have shown that already RBM NQS with \(\alpha=2\) can well approximate the DS and AF phases and shows the onset of the PS ordering. It also correctly identifies the point of the discontinuous change of the DS regime to PS/AF. However, in its vicinity the largest deviations from the exact results are observed. Here, the variational energy can be significantly lowered by increasing \(\alpha\), but this leads only to a small improvement of the estimation of order parameters. Consequently, we used RBM NQS with \(\alpha=2\) to address larger lattices and we confirmed the ability of RBM NQS to describe the main orderings of the SSM in zero magnetic field.

A gradual increase of the magnetic field in SSM leads to formation of stable plateaus in the magnetization each reflecting a differed ground state ordering. We have shown that RBM with \(\alpha=2\) can capture the relevant plateaus which form for here studied lattice sizes. Transfer learning can then be utilized to refine the results.

Figure 7: Comparison of the exact (red solid line) magnetization (a) and variational energy (b) results with VMC calculations utilizing RBM NQS with \(\alpha=2\) for lattices \(N=20,36\) and \(N=64\) as functions of external magnetic field.
To wrap it up, we have demonstrated that SSM is a good system for benchmarking NQSs and that a simple RBM NQS can be used to address its ground state in a broad range of regimes. This opens a possibility that NQSs could be used to address some unresolved questions related to the SSM, e.g., the existence of the spin-liquid phase and DQCP or the position, or to precisely capture the size and character of additional steps in the magnetization for larger lattices. We, however, leave this for future more focused studies.

## Acknowledgements

We thank Artur Slobodeniuk for the helpful discussions and Alberto Marmodoro for helping us to access additional computational resources.

Funding informationThis research was supported by the project e-INFRA CZ (ID:90140) of the Czech Ministry of Education, Youth and Sports (M.M., J.M.), BB. acknowledges a support from the Czech Science Foundation via Project No. 19-28594X, and M.Z. acknowledges a support from the Czech Science Foundation via Project No. 23-05263K.

## Appendix A Lattice tiles

To benchmark various architectures we utilize ED. We use the Lanczos algorithm implemented in SciPy library [61]. The only square (regular) systems tractable by this implementation, without extensive utilization of expected state symmetries, are of the size \(2\times 2\) and \(4\times 4\). Therefore we also constructed the so-called tilted regular-square clusters. They are depicted in figure 8 and each of them can be thought of as a single repeating building block of the infinite lattice. Clusters of sizes \(N=4,8,16,20\) are accessible via ED and used to benchmark our NQS implementations. However, in the text we are discussing only results for \(N\geq 16\)). We always use periodic boundary conditions.

Figure 8: Shapes of tilted tiles of sizes \(N=4,8,16,20\) used with periodic boundary conditions.
Visible biases in sRBM and pRBM

Here we show by contradiction that allowing uneven biases for sRBM is equivalent to constant biases when we enforce enough symmetries. Let us suppose visible biases are kept non-constant \(\left(a^{f}\to a_{i}^{f}\right)\) in the expression 9. We further suppose the condition that \(\forall i,j\ \exists g:g\sigma_{i}=\sigma_{j}\). This condition holds for every tile of SSM.

It follows that \(\sum\limits_{g\in G}T_{g}(\mathbf{\sigma}^{z})_{i}=C\sum\limits_{j=1}^{N}\sigma_{i }^{z}=Cm^{z}\), where \(C\) is the number of unique \(g\) fulfilling the condition. The first term in (9), after the generalization \(a^{f}\to a_{i}^{f}\), can be rewritten as

\[\sum\limits_{f=1}^{F}\sum\limits_{g\in G}\sum\limits_{i=1}^{N}a_{i}^{f}\,T_{g} (\mathbf{\sigma}^{z})_{i}=\sum\limits_{f=1}^{F}\sum\limits_{i=1}^{N}a_{i}^{f}\, \sum\limits_{g\in G}T_{g}(\mathbf{\sigma}^{z})_{i}=Cm^{z}\sum\limits_{f=1}^{F}\sum \limits_{i=1}^{N}a_{i}^{f}=Cm^{z}\sum\limits_{f=1}^{F}a^{f}. \tag{18}\]

Thus non-constant biases can be replaced by a constant value without loss of generality. Therefore, visible biases cannot be built into sRBM as independent variational parameters.

On the other hand, pRBM is not limited in this way. This can be clearly seen after rewriting both ansatze into similar forms. First the sRBM

\[\log\psi_{\mathbf{\theta}}(\mathbf{\sigma}^{z})=\log\prod\limits_{g\in G}\exp\sum \limits_{f=1}^{F}\left\{a^{f}\,\sum\limits_{i=1}^{N}T_{g}(\mathbf{\sigma}^{z})_{i} +\log\left[2\cosh\left(\sum\limits_{i=1}^{N}w_{i}^{f}\,T_{g}(\mathbf{\sigma}^{z})_ {i}+b^{f}\right)\right]\right\},\]

and then pRBM

\[\log\psi_{\mathbf{\theta}}^{G}(\mathbf{\sigma}^{z})=\log\sum\limits_{g\in G}\chi_{g^{- 1}}\exp\left\{\sum\limits_{i=1}^{N}a_{i}\,T_{g}(\mathbf{\sigma}^{z})_{i}+\sum \limits_{j=1}^{M}\log\left[2\cosh\left(\sum\limits_{i=1}^{N}W_{ij}\,T_{g}(\bm {\sigma}^{z})_{i}+b_{j}\right)\right]\right\}.\]

The sum (rather than product) of exponentials makes it impossible to use an analogous reduction of visible biases as discussed above for the sRBM. Note, that the usage of visible biases typically does not lead to a significant increase of parameters (\(+N\)). Yet, they usually improve the convergence of the learning process for frustrated systems, because they help to set the correct sign structure of the approximated state. Therefore, it is beneficial to include the visible biases to the NQS parameters whenever possible.

## Appendix C Symmetries

An infinite Shastry-Sutherland lattice has a \(p4g\) wallpaper group symmetry whose point group is \(C_{4v}\)[62]. The character table of \(C_{4v}\) is shown in Tab. 2. Each eigenstate of the SSM at infinite lattice must transform following one of the rows in the character table which, however, disregards the translations or glide reflections.

For finite lattices investigated in the paper the table and the number of additional translations depends on the system size and shape (note that we are using also irregular lattices). Different small clusters can have different character tables with varying number of irreducible representations (irreps) [63, 64]. A detailed analysis of each lattice goes beyond the scope of our paper. In practical implementations, we used the automorphisms of the graph via routines implemented in
NetKet [9, 51] and a particular line from its character table 2. For illustrative purposes, it is still useful to discuss the irreps of individual phases of SSM on the infinite lattice.

_DS_, described by Eq. 2, changes sign when we swap the spins in a dimer. More generally, the parity of the permutation determines the sign change. Consider a \(L\times L\) square lattice, where \(L\) is even, and a reflection symmetry along its diagonal axis (\(\sigma_{v}\)) within the squares containing the \(J^{\prime}\)-bonds. This axis cuts perpendicularly through \(L/2\)\(J^{\prime}\)-bonds. For each of these bonds a sign change occurs during the reflection while the sign of other dimers does not change. A similar argument can be constructed for the \(C_{4}\) rotation. This has an important implication even for finite lattices. Namely, for regular lattices the DS ground state transforms under the trivial irrep \(A_{1}\) if \(L\) is divisible by \(4\), and under the antisymmetric irrep (corresponding to \(B_{2}\)) otherwise. This has some important consequences for the use of symmetries of some finite lattices, as discussed in the main text.

_PS_ is twofold degenerate. Leaving out the translations this means that it transforms under irrep E, which is the only irrep of dimension \(2\).

Similar analysis of _AF_ state for finite lattices is rather complicated [63, 64]. If needed, we have assumed that AF transforms under trivial irrep \(A_{1}\) (with as well as without the application of MSR).
