Long-Range Entangled-Plaquette States for Critical and Frustrated Quantum Systems on a Lattice

TL;DR

This paper introduces a novel variational wave function called long-range entangled-plaquette state, capable of accurately modeling critical and frustrated quantum lattice systems by capturing long-range correlations and complex sign structures.

Contribution

The authors develop and validate a flexible, scale-free variational Ansatz that efficiently describes long-range entanglement and frustration effects in quantum lattice models.

Findings

Accurately reproduces ground state energies of XX and Heisenberg chains.

Captures long-range properties like Luttinger exponent and central charge.

Successfully models incommensurate correlations in frustrated systems.

Abstract

We explore a variational Ansatz for lattice quantum systems -- named long-range entangled-plaquette state -- in which pairs of clusters of adjacent lattice sites are correlated at any distance. The explicit scale-free structure of correlations built in this wave function makes it fit to reproduce critical states with long-range entanglement. The use of complex weights in the Ansatz allows for an efficient optimization of non positive definite states in a fully variational fashion, namely without any additional bias (arising \emph{e.g.} from pre-imposed sign structures) beyond that imposed by the parametrization of the state coefficients. These two features render the Ansatz particularly appropriate for the study of quantum phase transitions in frustrated systems. Moreover, the Ansatz can be systematically improved by increasing the long range plaquette size, as well as by the inclusion of even larger adjacent-site plaquettes. We validate our Ansatz in the case of the XX and Heisenberg chain, and further apply it to the case of a simple, yet paradigmatic model of frustration, namely the $J_1-J_2$ antiferromagnetic Heisenberg chain. For this model we provide clear evidence that our trial wave function faithfully describes both the short-range physics (particularly in terms of ground state energy) and the long-range one expressed by the Luttinger exponent, and the central charge of the related conformal field theory, which govern the decay of correlations and the scaling of the entanglement entropy, respectively. Finally we successfully reproduce the incommensurate correlations developing in the system at strong frustration, as a result of the flexible representation of sign (phase) structures via complex weights.

Figures

Equations

$$|\Psi ⟩=\sum _{\mathbf{\sigma }}\psi (\mathbf{\sigma })|\mathbf{\sigma }⟩$$

$${\psi }_{\mathrm{EPS}}(\mathbf{\sigma })=\prod _p{C}_p({\mathbf{\sigma }}_p)$$

$${\psi }_{n\mathrm{LR}-n^{\prime }\mathrm{A}-\mathrm{EPS}}(\mathbf{\sigma })={\psi }_{n\mathrm{LR}-\mathrm{EPS}}(\mathbf{\sigma }){\psi }_{n^{\prime }\mathrm{A}-\mathrm{EPS}}(\mathbf{\sigma }).$$

$${\psi }_{2LR-\mathrm{EPS}}(\mathbf{\sigma })$$

$$C_{ij}({\sigma }_i,{\sigma }_j)$$

$$⟨\mathcal{H}⟩=⟨E(\mathbf{\sigma };\mathbf{C})⟩=\sum _{\mathbf{\sigma }}E(\mathbf{\sigma };\mathbf{C})P(\mathbf{\sigma };\mathbf{C})$$

$$P(\mathbf{\sigma };\mathbf{C})=\frac{|\psi (\mathbf{\sigma };\mathbf{C}){|}^2}{{\sum }_{{\mathbf{\sigma }}^{\prime }}|\psi ({\mathbf{\sigma }}^{\prime };\mathbf{C}){|}^2}$$

$$E(\mathbf{\sigma };\mathbf{C})=\sum _{{\mathbf{\sigma }}^{\prime }}⟨\mathbf{\sigma }|\mathcal{H}|{\mathbf{\sigma }}^{\prime }⟩\frac{\psi ({\mathbf{\sigma }}^{\prime };\mathbf{C})}{\psi (\mathbf{\sigma };\mathbf{C})}$$

$$L_l(\mathbf{\sigma };\mathbf{C})=\frac{1}{\psi (\mathbf{\sigma };\mathbf{C})}\frac{\partial \psi (\mathbf{\sigma };\mathbf{C})}{\partial C_l},$$

$$g_l=⟨L_l^{\ast }E⟩-⟨L_l^{\ast }⟩⟨E⟩,$$

$$S_{l,m}=⟨L_l^{\ast }{L}_m⟩-⟨L_l^{\ast }⟩⟨L_m⟩,$$

$$\dot{\mathbf{C}}=-{\mathbf{S}}^{-1}\mathbf{g}.$$

$$L_{C_p({\bar{\mathbf{\sigma }}}_p)}=\frac{{\delta }_{{\mathbf{\sigma }}_p,{\bar{\mathbf{\sigma }}}_p}}{C_p({\mathbf{\sigma }}_p)}.$$

$$C^{\alpha \alpha }(r)=\frac{1}{N}\sum _i⟨S_i^{\alpha }{S}_{i+r}^{\alpha }⟩$$

$${\mathcal{S}}^{\alpha \alpha }(k)=\frac{1}{N}\sum _r{C}^{\alpha \alpha }(r)e^{ikr};$$

$${\mathcal{D}}_N=[D(N/2)-D(N/2-1)]/2$$

$$D(r)=\frac{1}{N}\sum _i\left[⟨({\mathbf{S}}_i\cdot {\mathbf{S}}_{i+1})({\mathbf{S}}_{i+r}\cdot {\mathbf{S}}_{i+r+1})⟩-⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+1}{⟩}^2\right].$$

$$\mathrm{Tr}{\rho }_A^2=⟨{\mathrm{SWAP}}_A{⟩}_2={⟨\frac{\psi ({\mathbf{\sigma }}_A^{\prime },{\mathbf{\sigma }}_B)\psi ({\mathbf{\sigma }}_A,{\mathbf{\sigma }}_B^{\prime })}{\psi ({\mathbf{\sigma }}_A,{\mathbf{\sigma }}_B)\psi ({\mathbf{\sigma }}_A^{\prime },{\mathbf{\sigma }}_B^{\prime })}⟩}_2$$

$$⟨...{⟩}_2=\sum _{{\mathbf{\sigma }}_A,{\mathbf{\sigma }}_B}\sum _{{\mathbf{\sigma }}_A^{\prime },{\mathbf{\sigma }}_B^{\prime }}|\psi ({\mathbf{\sigma }}_A,{\mathbf{\sigma }}_B){|}^2|\psi ({\mathbf{\sigma }}_A^{\prime },{\mathbf{\sigma }}_B^{\prime }){|}^2(...).$$

$$\mathcal{H}=-J\sum _i\left(S_i^x{S}_{i+1}^x+S_i^y{S}_{i+1}^y\right)$$

$${\psi }_{XX}(\mathbf{\sigma })=\mathcal{N}\prod _{i<j}\sin {\left[\frac{\pi }{N}(j-i)\right]}^{(2{\sigma }_i+1)(2{\sigma }_j+1)/4}$$

$$R_2(l)=\frac{c}{4}\log (l)+c_1+...$$

$$\mathcal{H}=J_1\sum _i{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+1}+J_2\sum _i{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+2},$$

$$⟨S_i^{+}{S}_{i+r}^{-}+S_i^{-}{S}_{i+r}^{+}⟩=\sum _{\mathbf{\sigma }}{\psi }^{\ast }(\mathbf{\sigma })\psi ({\mathbf{\sigma }}_{i,i+r})$$

$$|{\Psi }_{\mathrm{MF}}⟩={\otimes }_j\left(|{\uparrow }_j⟩+e^{iQ{r}_j}|{\downarrow }_j⟩\right)$$

$$R_2(l;N)=\frac{c}{4}\log d(l|N)+c_1+...$$

$$c_{N^{\prime },N}=\frac{4R_2^{N^{\prime }-N}}{\log (N^{\prime }/N)}$$

$${\theta }_{ij}({\sigma }_i,{\sigma }_j)=\frac{\pi }{N_i}\frac{{\sigma }_i+1}{2}$$

$${\theta }_{ij}({\sigma }_i,{\sigma }_j)=\pi N_{A,\uparrow }^{(ij)}=\frac{\pi }{2}[({\sigma }_i+1){\delta }_{i,A}+({\sigma }_j+1){\delta }_{j,A}]$$

$$\mathrm{sgn}[\psi (\mathbf{\sigma })]$$

$$\mathcal{M}\!=\!\frac{1}{N_{C}}\!\!\sum_{\bm{\sigma}}\frac{1}{2}\Bigg{|}{\rm sgn}\left[{\rm Re}\left(\frac{\psi(\bm{\sigma})}{\psi(\bm{\sigma}_{\rm{REF}})}\right)\right]\!+\!\frac{(-1)^{N_{A,\uparrow}(\bm{\sigma})}}{(-1)^{N_{A,\uparrow}(\bm{\sigma}_{\rm{REF}})}}\Bigg{|}$$

Full text

I Introduction

Knowing the equilibrium state of quantum many-body systems is one of the central problems of modern theoretical physics. Similar to the classical case, this problem can be generally cast in the form of the evaluation of a statistical sum over a number of configurations growing exponentially with the number of degrees of freedom, and addressed with a stochastic approach, defining the general strategy of quantum Monte Carlo (QMC) techniques Sandvik (2010); Becca and Sorella (2017). In contrast to the classical case, however, the statistical sum may have weights which are not positive definite (or may even be complex) leading to the so called sign (or phase) problem Chandrasekharan and Wiese (1999); the latter imposes the price of an exponentially large statistics (in the system size or in the inverse temperature) for reliable results to be obtained. Such a fundamental limitation currently impairs significantly our understanding of strongly correlated fermionic systems (from models of electrons in solids to models of elementary particles in relativistic quantum field theory), frustrated quantum magnetism, and bosonic quantum particles in gauge fields, to cite a few examples.

In this context, a general alternative strategy is offered by the variational approach Becca and Sorella (2017), mostly focusing on the ground state of the many-body problem of interest, and consisting in a chosen parametrization (Ansatz) of the ground-state wave function in terms of a reduced number of parameters (polynomial in the system size). In the following, we shall specialize our discussion to $S=1/2$ quantum spin models on a lattice, and to the corresponding bosonic or fermionic Hamiltonian that they can be mapped onto. Hence, the general form of the wave function reads:

[TABLE]

where, for a lattice of $N$ sites, $|\bm{\sigma}\rangle=|\sigma_{1},...,\sigma_{N}\rangle$ is the eigenvector of e.g. the $S_{i}^{z}$ operator at each site $i$ ($\sigma_{i}=\pm 1/2$). Providing a variational wave function ultimately means choosing a suitable form of the wave function coefficients (or weights) $\psi(\bm{\sigma})\rightarrow\psi(\bm{\sigma};\bm{C})$ depending on a set of adjustable parameters $\bm{C}=\{C_{l}\}$, where the index $l$ is here used to enumerate the elements of $\bm{C}$. Over the last decades many variational Ansätze have been formulated to describe the ground-state physics of lattice spin models escaping the reach of unbiased QMC approaches. A few examples are matrix product states (MPS) Schollwöck (2005, 2011), various tensor network states (TNS) Orús (2014); Verstraete et al. (2008), resonating valence bond (RVB) states Liang et al. (1988); Becca et al. (2010), neural network quantum states (NNQS) Carleo and Troyer (2017), as well as entangled-plaquette states (EPS) Mezzacapo et al. (2009), which are also known as correlated-product states (CPS) Changlani et al. (2009), and constitute the focus of this work. The crucial aspect for the success of a variational Ansatz is its ability to reproduce the entanglement and correlations expected for the ground state of the Hamiltonian of interest. The entanglement content and correlation properties of MPS and different TNS have been extensively characterized Schollwöck (2011) and recognized as one of their main limitations, while other variational states are more flexible on the entanglement content Gao and Duan (2017). An important bias, common to any variational form, is offered by the specific parametrization of the wave function coefficients $\psi(\bm{\sigma})$, and in particular of their sign structure, which generally becomes non-trivial when dealing with the ground state of a frustrated spin model Becca and Sorella (2017). In the case of RVB states the sign structure of the coefficients is inherited from that of a fermionic state projected onto the spin Hilbert space to provide the Ansatz Becca and Sorella (2017); only a few studies have dealt so far with the sign structure of NNQS, either imposing it a priori Carleo and Troyer (2017); Choo et al. (2019) or modeling it using an auxiliary neural network Cai and Liu (2018).

In this work we introduce the long-range entangled-plaquette states (LR-EPS) which offer a very flexible parametrization of highly entangled quantum spin wave functions. The major strength of EPS in general relies on the possibility of explicitly correlating different groups of lattice sites into plaquettes, expressing $\psi(\bm{\sigma})$ as a product of plaquette coefficients $C_{p}(\bm{\sigma}_{p})$, where $p$ is the plaquette index and $\bm{\sigma}_{p}$ is the plaquette configuration. Designing the EPS Ansatz by considering overlapping (i.e. entangled) plaquettes constitutes a fundamental aspect which will be discussed in the following section. LR-EPS are a natural generalization of the entangled-plaquette wave functions based on plaquettes of adjacent sites (hereafter referred as A-EPS) mostly adopted in previous works on the subject Mezzacapo et al. (2009); Mezzacapo and Cirac (2010); Mezzacapo (2011); Al-Assam et al. (2011); Mezzacapo and Boninsegni (2012); Mezzacapo (2012); Durič et al. (2014, 2016); Mezzacapo et al. (2016). The idea of the LR-EPS Ansatz is to consider plaquettes directly correlating clusters of adjacent lattice sites at arbitrary distances. This results in a potentially scale-invariant wave function which can accurately capture long-range entanglement and effectively describe critical ground states (associated with quantum critical points, gapless ordered phases or extended critical ones). Moreover the use of complex plaquette coefficients $C_{p}(\bm{\sigma}_{p})$ allows one to reproduce various sign/phase structures, turning the latter into a feature that can be fully optimized variationally. Hence, the LR-EPS Ansatz marries the relative simplicity of the EPS concept with a high degree of flexibility in coding the correlations and the sign structure of a quantum state. Moreover, as any EPS wave function, it can be systematically improved by increasing the size of the (long-range) plaquettes, as well as by combining adjacent-site and long-range plaquettes. We first validate our Ansatz in the case of the $S=1/2$ XX chain. Remarkably the ground state of this model, along with that of other exactly solvable ones Sutherland (2004), is a LR-EPS, whose form our optimization algorithm is able to accurately reconstruct without any bias except for the constraint of translational invariance. This is proven by our variational results being able to reproduce the exact correlations and entanglement structure of the $S=1/2$ XX chain with excellent precision. We then move on to the frustrated $J_{1}-J_{2}$ chain, for which our Ansatz is shown to successfully reproduce the most challenging traits of the physics, namely the quantum phase transition from gapless spin liquid to valence-bond crystal; and the appearance of incommensurate correlations as the degree of frustration is increased. In particular we provide extensive results for the evolution of entanglement properties across the above cited transition.

The paper is structured as follows. Sec. II briefly recalls the basics of the EPS wave function, presents in detail the form of the LR-EPS Ansatz, its parametrization of the sign structure, the optimization strategy, and the observables relevant to this study; Sec. III shows a validation of the Ansatz in the case of the exactly solvable XX chain while Sec. IV, discusses our results for the $J_{1}-J_{2}$ chain; conclusions and perspectives are presented in Sec. V.

II Long-range entangled-plaquette states: form, optimization and observables

II.1 Structure of the Ansatz

The EPS Ansatz amounts to parametrizing the wavefunction weights as

[TABLE]

where $C_{p}(\bm{\sigma}_{p})$ are coefficients associated with plaquettes indexed with $p$ and of size $n_{p}$, and $\bm{\sigma}_{p}=(\sigma_{1,p},...,\sigma_{n_{p},p})$ is the plaquette configuration in the computational basis ($j,p$ is the index of the $j$-th site of the $p$-th plaquette). In general terms, the plaquettes may have a completely arbitrary geometry, although their size $n_{p}$ is necessarily limited to $\sim{\cal O}(10)$, due to the exponential growth of the parameter space with $n_{p}$. Assuming for simplicity that all plaquettes have the same size $n_{p}=n$ and that the wave function contains $M$ of them, the number of variational parameters $C_{p}(\bm{\sigma}_{p})$ is $M\times 2^{n}$. If non overlapping (i.e., disentangled) plaquettes with $n>1$ are considered, one recovers a cluster mean-field Ansatz in which inter-plaquette correlations and entanglement are fully neglected (see, e.g., Fig. 1(a) for plaquettes of 2 adjacent sites). The crucial aspect of EPS is that plaquettes can be overlapped (i.e., entangled), introducing entanglement and correlations over distances larger than the plaquette size.

Most of the literature on EPS Mezzacapo et al. (2009); Mezzacapo and Cirac (2010); Al-Assam et al. (2011); Mezzacapo and Boninsegni (2012); Mezzacapo (2012); Durič et al. (2014, 2016); Stojevic et al. (2016) has considered Ansätze based on $M\leq N$ overlapping plaquettes of $n$ adjacent sites ($n$A-EPS), and improved them systematically by increasing $n$ (Fig. 1(b,c)). While the $n$A-EPS Ansatz is asymptotically exact in the large $n$ limit, and any finite-$n$ A-EPS is a valid variational choice, convergence of the results in the plaquette size is often not achievable, as a result of the mentioned exponential cost of increasing $n$. This aspect becomes particularly serious in the vicinity of a quantum phase transition. For instance, as numerically found in Ref. Stojevic et al. (2016), an $n$A-EPS Ansatz with partially overlapping uniform plaquettes generally exhibits a finite correlation length (proportional to $n$), preventing the correct description of the critical regime at fixed $n$. In general, the extrapolation of the results to the large-$n$ limit is necessary to correctly determine the ground state phase boundaries of a given model of interest Mezzacapo (2012); Mezzacapo and Boninsegni (2012).

An extremely simple alternative to the above scheme consists in fixing the size of the plaquettes while playing with their geometry, which can be promoted to include (clusters of) sites at arbitrary distances Changlani et al. (2009); Al-Assam et al. (2011). In this way correlations can be established directly among arbitrarily distant sites. This is the strategy underlying the long-range EPS Ansatz that we consider here in order to tackle systems that develop a diverging correlation length in their ground state. In particular, we shall focus on plaquettes with an even number of sites, $n$, and establish their geometry as composed of two clusters of $n/2$ adjacent sites at an arbitrary distance from each other: this ultimately defines the $n$LR-EPS wave function as illustrated in Fig 1(d,e). The $n$LR-EPS clearly extends the $n$A-EPS, explicitly incorporating dominant $n$-site correlations at all length scales. Such an improvement occurs at a computational price now independent of $n$, and polynomial in the system size.

In general, the $n$LR-EPS Ansatz is systematically improvable upon increasing $n$, similarly to the A-EPS Ansatz. However, unlike in the A-EPS Ansatz, this improvement procedure is not a priori crucial, since the long-range physics of a given Hamiltonian can be well captured already with $n=2$ (see Secs. III and IV). Another strategy for the improvement of the Ansatz that we shall pursue here is to combine both adjacent-site as well as long-range plaquettes, to give rise to $n$LR-$n^{\prime}$A-EPS (where $n$ ($n^{\prime}$) is the size of the long-range (adjacent-site) plaquettes, with $n^{\prime}>n$) exemplified in Fig. 1(f). The wave-function coefficients of this Ansatz are defined via:

[TABLE]

The latter strategy allows one in principle to faithfully describe $n$ point correlations at all distances, also improving on local properties such as the energy of short-range interacting models.

In the particular case of the 2LR-EPS Ansatz and $S=1/2$ spins, for a $N-$site lattice $M=N(N-1)/2$ plaquettes are formed by each pair of sites $i$ and $j$, i.e., $p=(ij)$ in Eq. (2), and the wave function coefficients can be explicitly rewritten as

[TABLE]

establishing a link with what goes under the name of spin-Jastrow Ansatz in the previous literature Pang (1991). Therefore, our approach generalizes systematically spin-Jastrow states, both in terms of the size of the plaquettes that are explicitly correlated in the form of the Ansatz; as well as via the use of complex plaquette coefficients, [see detail in Sec. II.3]. It goes without saying that the EPS Ansatz applies to any lattice geometry and to any number of spatial dimensions without requiring particular modifications. The explicit correlations introduced within all pairs of clusters allows one to describe long-range entanglement in any such situation.

II.2 Variational energy minimization

Given the Hamiltonian ${\cal H}$, the variational optimum is searched via imaginary-time evolution projected onto the space of states compatible with the variational Ansatz via the time-dependent variational principle Kramer and Saraceno (1981) - this approach is equivalent to the so-called stochastic reconfiguration scheme Becca and Sorella (2017); Carleo and Troyer (2017). The variational energy to be minimized reads

[TABLE]

where

[TABLE]

and $E(\bm{\sigma};\bm{C})$ is the energy estimator

[TABLE]

Introducing the logarithmic derivatives with respect to the $l$-th variational parameters of the wave function coefficients

[TABLE]

the gradient

[TABLE]

as well as the covariance matrix

[TABLE]

the projected imaginary-time dynamics of the variational Ansatz is then described by the equation

[TABLE]

The simple form of the EPS Ansatz in Eq. (2) allows for a straightforward calculation of the $L_{l}$’s, whose expression reduces, for $\bm{C}=\{C_{p}(\bm{\sigma}_{p})\}$, to

[TABLE]

Statistical sums over the configurations $\bm{\sigma}$, contained in Eqs. (5), (9) and (10), are sampled via a Monte Carlo scheme which makes use of Metropolis updates based on the exchange of spins with opposite $\sigma_{z}$. Hence, the total spin along $z$, set to [math] in the initial state, is conserved along the simulation. Further details on optimization strategies will be described in the next paragraph.

II.3 Complex wave function coefficients

In order to search for the ground state of frustrated magnetic models with a real-valued Hamiltonian, it is crucial to be able to account for weights $\psi(\bm{\sigma})$ of the wave function with positive as well negative sign Becca and Sorella (2017). The function ${\rm sgn}[\psi(\bm{\sigma})]$ defines the so-called sign structure of the state, and it is in general unknown. In fact, one can argue that knowing it a priori would result in being able to solve the ground state problem by using QMC techniques based on ground-state projection (such as Green’s function Monte Carlo Becca and Sorella (2017)) with fixed-node constraints. Indeed, the sign structure of the ground state is known to be trivial (all positive weights) when the off-diagonal elements of the Hamiltonian matrix $\langle\bm{\sigma}^{\prime}|{\cal H}|\bm{\sigma}\rangle$ are semi-negative definite (Perron-Frobenius theorem), or can be made so by an unitary transformation. This is the case for instance of antiferromagnetic Hamiltonians defined on a bipartite lattice, namely containing only antiferromagnetic interactions between sites belonging to different sublattices ($A$ and $B$): in that case the unitary transformation (amounting to a $\pi$-angle spin rotation on all $A$ spins) leads to a sign structure determined by the Marshall sign rule Auerbach (2006): ${\rm sgn}[\psi(\bm{\sigma})]=(-1)^{N_{A,\uparrow}(\bm{\sigma})}$ (where $N_{A,\uparrow}(\bm{\sigma})$ is, for the generic configuration $\bm{\sigma}$, the number of $\uparrow$ spins on the $A$ sublattice).

In order to tackle generic systems, in which the ground-state signs are completely unknown, an ideal variational Ansatz should therefore be able to reproduce a wide variety of possible sign structures, and to do so in a continuous manner, so that the sign structure can be a subject of variational optimization. In the EPS Ansatz the sign of $\psi(\bm{\sigma})$ is simply given by the product of the signs of the plaquette coefficients $C_{p}(\bm{\sigma}_{p})$, introducing a fundamental bias on the possible sign structures that the Ansatz can realize. Nonetheless it is a simple exercise to prove that the Marshall sign structure can indeed be reproduced by the EPS Ansatz, as discussed in Appendix A. From the point of view of the numerical optimization, an EPS Ansatz with real plaquette parameters (obeying the dynamics governed by the equation Eq. (11)) should in principle lead to changing the signs of the wave function weights, so as to explore non-trivial sign structures in a purely variational way. However, in practice, a real coefficient $C_{p}(\bm{\sigma}_{p})$ approaching zero (necessary for a sign change) entails that the weight $\psi(\bm{\sigma})$ of all the configurations $\bm{\sigma}$ compatible with $\bm{\sigma}_{p}$ on the $p$-th plaquette also approaches zero, due to the multiplicative structure of the EPS Ansatz. This in turn makes the appearance of the $\bm{\sigma}_{p}$ configuration in the update very rare, so that an enormous statistics has to be accumulated in order to properly sample the corresponding gradient (Eq. (9)) and covariance matrix (Eq. (10)). A very simple strategy to circumvent this issue is to extend the weights to the complex plane by considering complex-valued $C_{p}(\bm{\sigma}_{p})$ so that sign changes (or $\pi$ phase shifts) can be made without ever crossing the origin: this allows then for a full variational optimization of the phase structure of the Ansatz, which is the strategy that we pursue in this study. The use of complex $C_{p}(\bm{\sigma}_{p})$ should be seen as an extension of the variational space in order to be able to recover the correct variational optimum which, for a real-valued Hamiltonian matrix, should also be real up to a global phase. Our variational optimization procedure indeed leads to states consistent with this scenario.

From a technical point of view, it is convenient to parametrize the plaquette coefficients in their polar decomposition $C_{p}(\bm{\sigma}_{p})=A_{p}(\bm{\sigma}_{p})e^{i\theta_{p}(\bm{\sigma}_{p})}$ This allows, for example, to optimize the amplitude $A$ and phase $\theta$ as independent variables. Also, one can optimize the phases at first while keeping the amplitudes fixed and equal to one: in doing so one produces a Monte Carlo dynamics for the sampling of the statistical sums of Eqs. (5), (9) and (10), where all configurations are equally probable. After a possible pre-optimization of the phases (achieved when the variational energy ceases to decrease), the amplitudes are left free to vary (i.e., to depart from their initially unit value) and optimized alongside with the phases. In this way any nodal and sign structure compatible with the Ansatz may emerge from the optimization in an unbiased fashion. Moreover, the introduction of spatial symmetries in the plaquette parameters can be done separately for the amplitude and phase variables, namely the functions $A_{p}$ and $\theta_{p}$ can be made to depend on the plaquette $p$ index in a different way. Explicit examples will be provided in Sec. IV.

II.4 Correlation functions and entanglement entropies

The main focus of our present work is on the ability of the nLR-EPS Ansatz to correctly capture the correlation and entanglement properties of quantum spin states with small (or even minimal, i.e., $n=2$) plaquette size. All of the findings presented in this study concern one-dimensional $S=1/2$ quantum spin models, described by spin operators $S_{i}^{\alpha}$, where $\alpha=x,y,z$ and the index $i$ runs on the lattice sites of a linear chain. In the following we shall present results for the spin-spin correlation function

[TABLE]

and the related structure factor

[TABLE]

as well as for the dimer order paremeter defined as in Ref. Sandvik, 2010 as

[TABLE]

where

[TABLE]

is the dimer-dimer correlation function. Moreover we will concentrate on the 2-Rényi entanglement entropy $R_{2}(A)=-\log({\rm Tr}\rho^{2}_{A})$ where $\rho_{A}={\rm Tr}_{B}|\Psi\rangle\langle\Psi|$ is the reduced density matrix describing the subsystem $A$ after having traced out the degrees of freedom of its complement $B$. The purity ${\rm Tr}\rho^{2}_{A}$ can be conveniently calculated as the expectation value of the SWAP operator of the configuration of the $A$ subsystem between two replicas of the whole system Hastings et al. (2010). Denoting with $\bm{\sigma}_{A}$ and $\bm{\sigma}_{B}$ the configurations of the subsystems $A$ and $B$ in a state $\bm{\sigma}=(\bm{\sigma}_{A},\bm{\sigma}_{B})$ of the computational basis, one has

[TABLE]

where the two-replica statistical average $\langle...\rangle_{2}$ is defined as

[TABLE]

II.5 Alternative variational Ansätze for quantum spin models and comparison with LR-EPS

In this section we briefly review some of the most popular variational states for lattice spin models, and contrast their properties with those of our LR-EPS Ansatz.

The most successful example of a variational Ansatz is represented by matrix-product states (MPS), which represent the variational Ansatz optimized by the density-matrix renormalization group (DMRG) algorithm as well as other related techniques Schollwöck (2005). The MPS capture with impressive precision the physics of one-dimensional quantum systems with entanglement entropies obeying an area-law scaling with subsystem size, including possible logarithmic corrections. Indeed the maximum sub-system entanglement entropy that the Ansatz can allow for is given by $\log D$ where $D$ is the linear dimension of the matrices composing the Ansatz and allows for its systematic improvement. In this framework, only logarithmic scalings of the entanglement entropy with subsystem size are tolerable in order to achieve a polynomial scaling of the number of variational parameters with system size.

The direct application of MPS to models in higher spatial dimension is in principle problematic. In fact, in this case, the $D$ parameter scales exponentially with the subsystem size even in the case of area-law states Stoudenmire and White (2012). This issue can be circumvented by generalizing MPS to tensor-network states (TNS) Orús (2014), which are suited to study the physics of area-law states in two and higher dimensions. Like MPS, TNS can be systematically improved by enlarging the number of parameters with a concomitant polynomial scaling in the computational cost; nonetheless TNS generally exhibit a finite correlation length Corboz et al. (2018); in several formulations (such as the celebrated projected entangled-pair states Verstraete et al. (2008)) they do not allow for an efficient exact calculation of the wave function coefficients Glasser et al. (2018) nor for an efficient representation of quantum states with faster than area-law entanglement scaling. Many Ansätze offer valuable alternatives to MPS/TNS with rather complementary properties. A famous example is given by resonating valence-bond (RVB) states Liang et al. (1988), efficiently parametrized as projected Bardeen-Cooper-Schrieffer (pBCS) states Anderson (1987). Their sign structure descends from that of a fermionic determinant which provides the ground state of a BCS-like Hamiltonian, and as such it can be highly non-trivial. As a consequence the pBCS states have been successfully applied to several frustrated models of quantum magnetism Becca and Sorella (2017); nonetheless they cannot be systematically improved in their sign structure – although one may argue that they could be combined e.g. with complex-valued Jastrow factors altering both the amplitude and sign structure of the coefficients. More recently the ability of neural networks to reproduce a function of many variables has been exploited to parametrize in this form the weights $\psi(\bm{\sigma};\bm{C})$, defining the Ansatz called neural-network quantum states (NNQS) Carleo and Troyer (2017). The latter proves to be very effective both for unfrustrated Carleo and Troyer (2017) as well as frustrated models of magnetism, with a fixed sign structure Choo et al. (2019) or by using complex coefficients Glasser et al. (2018). Ref. Cai and Liu (2018) recently explored the possibility of parametrizing the signs by using a dedicated neural network, but without complexification of the coefficients. NNQS as formulated in Ref. Carleo and Troyer (2017) can be systematically improved by increasing the depth of the neural network, although networks beyond single-layer ones do not allow for an efficiently calculable form of the coefficients (see Ref. Choo et al. (2019) for a multi-layer convolutional network representation). LR-EPS with complex coefficients appear as a valuable alternative to all the above variational schemes, because of the relative simplicity of their structure and the flexibility of their formulation, with correlations and long-range entanglement explicitly built in. EPS come with a large variety of improvement schemes (e.g., extension of the plaquette sizes, combination of long-range and adjacent-site plaquettes) as well as with the possibility to achieve a full variational optimization of the sign structure (within the structures compatible with the Ansatz) by generalizing it to the phase structure of a complex wave function, as described above.

III Validation of the Ansatz: XX chain

A first validation of our approach comes from the case of the $S=1/2$ XX chain, with Hamiltonian

[TABLE]

which is exactly solvable by mapping it to free fermions via the Jordan-Wigner transformation Lieb et al. (1961). Here and throughout the rest of the paper we consider chains with periodic boundary conditions. This Hamiltonian is not frustrated and satisfies the Perron-Frobenius theorem, so that its ground state wave function has all positive-defined coefficients. One can therefore search for its ground state in the form of a LR-EPS with real coefficients. In the context of LR-EPS, the case of the XX chain is in fact extremely special, because the exact coefficients of the fermionic ground state, given by the Slater determinant, can be recast, in the form of a Vandermonde determinant Sutherland (2004), as

[TABLE]

where ${\cal N}$ is a normalization factor. It is immediate to recognize that the above espression is fully compatible with that of the 2LR-EPS Ansatz, Eq. (4). In fact the 2LR-EPS form encompasses an entire family of many-body wavefunctions, such as the exact ground state of one-dimensional bosons or fermions with inverse-square interaction potential Sutherland (1971), the exact ground state of the Haldane-Shastry model Haldane (1988); Shastry (1988), or the Laughlin wavefunction Changlani et al. (2009), to cite a few relevant examples.

In the specific case of the XX chain, the ground-state physics features a critical Luttinger-liquid phase with algebraically decaying correlations $C^{xx(yy)}\sim r^{-1/(2K)}$ with Luttinger liquid exponent $K=1$ Giamarchi (2003), and a logarithmic scaling of the entanglement entropy of a subsystem of linear size $l$:

[TABLE]

with $c_{1}$ a non universal constant and $c=1$ giving the central charge of the related conformal field theory Calabrese and Cardy (2004) (the missing terms are subdominant corrections). These traits are in fact universal to many one-dimensional systems admitting a description as a Luttinger liquid Giamarchi (2003), whose scale-free, critical nature is nicely captured by the spatial structure of correlations inscribed in the LR-EPS Ansatz. Indeed Jastrow wave functions (to which a real-valued 2LR-EPS can be specialized) have already been assessed in the past as very efficient descriptions of Luttinger-liquid phases of both interacting bosons and fermions Capello et al. (2005, 2008).

Here we use the simplest form of the Ansatz, namely a 2LR-EPS where we enforce translational symmetry, so that only ${\cal O}(N)$ variational parameters are required. A very high accuracy can be obtained for both the off-diagonal spin-spin correlation function $C^{xx}(r)$ [Fig. 2] and the 2-Rényi entropy [Fig. 3]. In particular it is remarkable to see that, for the energy, the error relative to the exact results is less than $10^{-8}$; in fact slightly lower than the statistical error that we have on the energy estimator itself. This means that, for all purposes, our optimization algorithm finds the absolute minimum of the variational energy, and that our accuracy on the reconstruction of the ground state is only limited by the statistical uncertainty. The relative accuracy of the LR-EPS predictions remains rather good ($\lesssim 10^{-3}$ for correlations, $\lesssim 10^{-4}$ for the entanglement entropy) when looking at long-distances, without any significant degradation of our estimates being observed for increasing $r$, or $l$.

IV $\bm{J_{1}-J_{2}}$ chain

IV.1 Model and phase diagram

A more stringent test of our Ansatz is provided by the frustrated antiferromagnetic $J_{1}-J_{2}$ chain whose Hamiltonian is

[TABLE]

This model features a very rich phase diagram upon changing the degree of frustration $\alpha=J_{2}/J_{1}$. At low frustration, its ground state realizes a gapless phase (described as a Luttinger liquid) continuously connected to its unfrustrated limit $J_{2}=0$. When $\alpha=\alpha_{c}=0.241167....$ Eggert (1996) the system undergoes a quantum phase transition to a gapped, valence-bond crystal (VBC) with spontaneous dimerization. Beyond the exactly solvable Majumdar-Ghosh (MG) point $\alpha=1/2$ Majumdar and Ghosh (1969), correlations start developing an incommensurate structure Bursill et al. (1995); White and Affleck (1996); Gehring et al. (1997); Deschner and Sørensen (2013) with a pitch vector that evolves continuously towards $\pi/2$ as $\alpha$ increases. After its seminal numerical investigations via DMRG Bursill et al. (1995); White and Affleck (1996); Gehring et al. (1997), Eq. (22) has been successfully investigated in the more recent past with RVB states Ferrari et al. (2018), as well as with NNQS Cai and Liu (2018).

In order to capture the rich phenomenology correctly with a variational approach, it is crucial to use a wave function with a flexible, continuously adjustable sign structure. This is clearly evidenced by the fact that correlation functions can exhibit different sign patterns, going from the staggered one (related to a pitch vector $\phi=\pi$) for small $\alpha$ to incommensurate sign patterns for $\alpha>0.5$. Indeed the sign pattern of off-diagonal correlations $C^{xx(yy)}(r)$ is related to the mutual signs of coefficients corresponding to configurations connected by two spin flips at distance $r$. Indeed,

[TABLE]

where $\bm{\sigma}_{i,i+r}$ corresponds to the $\bm{\sigma}$ configuration with flipped $i$-th and $(i+r)$-th spins. Therefore, off-diagonal correlations with a non-trivial sign pattern are consequence of the sign structure of the wave function coefficients. A closer look at the ground-state sign structure of model (22) has been offered in Ref. Zeng and Parkinson (1995), pointing out that a Marshall sign for a two-sublattice structure $ABAB...$ appears at weak frustration $\alpha\ll 1$, while a Marshall-like sign with sublattice structure $ABBA...$ emerges in the opposite limit $\alpha\gg 1$. For intermediate values $\alpha\sim{\cal O}(1)$ a definite sign structure could not be identified, and most likely it evolves with $\alpha$ along with the incommensurate correlations.

It is instructive to remind the reader that the mean-field (MF) solution of the problem gives a wavefunction

[TABLE]

where $Q=\pi$ for $\alpha<1/4$ and $Q=\cos^{-1}[-1/(4\alpha)]$ for $\alpha>1/4$. Therefore dressing the MF wavefunction with a real correlation term (e.g. of the Jastrow type) would already produce non-trivial incommensurate correlations. Nonetheless the exact solution of the ground-state problem of a time-reversal invariant Hamiltonian (such as the one under investigation) is given by a state $|\Psi_{0}\rangle$ that can only possess a time-reversal invariant total momentum – namely, under a translation $T_{\delta}$ of a distance $\delta$, $T_{\delta}|\Psi_{0}\rangle=e^{iP\delta}|\Psi_{0}\rangle$ with $P=0$ or $\pi$. On the other hand, projecting the above mean-field state onto the sector at zero total magnetization ${\cal P}_{S^{z}=0}|\Psi_{\rm MF}\rangle$ produces a state which has a momentum $P=NQ/2$ (because upon translation of $\delta$ sites the $N/2$ $\downarrow$ spins produce a phase factor $\exp(iNQ\delta/2)$). This situation, however, results in a ground-state wave function leading to both incommensurate correlations and an unphysical breaking of the time-reversal symmetry (i.e., the ground state remains complex-valued). Conversely, since our optimized Ansatz is essentially real (see discussion in Sec. IV.2, IV.3.2, and Appendix B) the emergence of incommensurate correlations discussed in this work cannot be due to the above scenario.

IV.2 Validation of the Ansatz: comparison with exact results

In light of the above observations, the 1$d$ $J_{1}-J_{2}$ model offers a rather challenging testbed for our variational Ansatz. In the following we shall check the accuracy of our results against exact ones by comparing ground-state energies and correlation functions. First of all, in the case of the unfrustrated Heisenberg chain ($J_{2}=0$) the optimization of complex-valued LR-EPS reproduces faithfully the physics of the system: when, for instance, the calculation starts from a state with random coefficients of unit norm (i.e., $C_{p}(\bm{\sigma}_{p})=A_{p}(\bm{\sigma}_{p})e^{i\theta_{p}(\bm{\sigma}_{p})}$ with $\theta_{p}(\bm{\sigma}_{p})$, and $A_{p}(\bm{\sigma}_{p})$ initially chosen as a random phase between [math] and $2\pi$, and $1$, respectively), the optimization algorithm is capable of finding efficiently the Marshall sign structure, as witnessed by the correct sign pattern reconstructed for the correlation function $C(r)=\sum_{\alpha}C^{\alpha\alpha}(r)$ in Fig. 4. This shows that the variational search of the sign structure, when compatible with the LR-EPS multiplicative form, can indeed be efficiently performed. Moreover, our estimates are in very good agreement with numerically exact QMC results based on the stochastic series expansion Sandvik (2010). We note that, by means of a $n$A-EPS wave function explicitly including the Marshall signs, one may achieve for the present model and system size an accuracy on the ground state energy similar to that obtained with the 2LR-EPS, when $n$ is as large as $12-14$ sites (i.e., via a considerably larger number of variational parameters than in the 2LR-EPS case).

Results at finite frustration compared with exact diagonalization ones are reported in Fig. 5 and 6. In paricular, Fig. 5 (upper panel) shows the $\alpha$ dependence of the variational energy of different LR-EPS states for a $N=16$ chain with periodic boundary conditions. The simplest Ansatz that we test is the 2LR-EPS one, where we do not impose any symmetry, namely we optimize $2N(N-1)$ independent complex coefficients. We observe that the 2LR-EPS Ansatz remains very accurate (with relative errors in the $0.1\%$ range) up to $\alpha\simeq 0.5$. For larger values of $\alpha$ the precision degrades, with relative errors rising to the $1\%$ level) – yet the Ansatz can be systematically improved by moving to the 4LR-EPS or to the mixed long-range/adjacent-site 2LR-8A-EPS Ansatz (i.e., by increasing, if no symmetries are considered, the number of variational parameter by a factor of about $3.5$ or $9.5$, respectively). Nonetheless the minimal 2LR-EPS wave function is already capable of reproducing the non-trivial and $\alpha$-dependent sign pattern of the correlations as shown in Fig. 6. In our view, this offers valuable evidence of the fact that our simplest wavefunction can describe the different physical regimes of the model in question, without the need to further enrich its parametrization. Therefore we shall focus on this wavefunction for the study of larger system sizes, further reducing its variational parameters by imposing symmetries (see below).

A further significant test of the quality of our results is offered by the real vs. complex nature of the wavefunction coefficients. Randomly choosing pairs of basis states $|\bm{\sigma}\rangle$ and $|\bm{\sigma}^{\prime}\rangle$ we find that the ratio of their coefficients $\psi_{\rm 2LR-EPS}(\bm{\sigma})/\psi_{\rm 2LR-EPS}(\bm{\sigma}^{\prime})$ is a real number (whenever both coefficients are larger than machine precision in modulus). This is a direct proof that the variational optimization of the complex Ansatz is able to eventually align the phases of all the wavefunction coefficients (modulo $\pi$), returning a real valued wavefunction up to an irrelevant global phase. This offers additional evidence that the strategy of producing non-trivial sign structure by optimizing a complex wavefuction is successful.

IV.3 Correlations and entanglement on larger lattices

Having validated our Ansatz for the frustrated $J_{1}-J_{2}$ chain for a small lattice size, we move on to studying the physics of the same chain for larger lattices (up to $N=80$). Our goal is to show that the LR-EPS Ansatz is a very good tool to study fundamental features of the correlations and the entanglement of the ground state of the system, and their evolution across the phase diagram of the model. Throughout this section we specialize our attention to a 2LR-EPS Ansatz with coefficients $C_{ij}(\sigma_{i},\sigma_{j})=A_{ij}(\sigma_{i},\sigma_{j})e^{i\theta_{ij}(\sigma_{i},\sigma_{j})}$ where, without loss of generality, we consider $i<j$. We parametrize the spatial dependence of both amplitudes $A_{ij}$ and phases $\theta_{ij}$ in terms of the coordinate $i$ of the first site and of the distance between the sites $d=j-i$. In order to reduce the number of variational parameters, we chose the dependence on $(i,d)$ to be periodic of period $(p_{i},p_{d})$. All the results for lattices of size $N>16$ and in particular those in the following sections are obtained via wave functions with periods $(p_{i},p_{d})=(2,N)$ for the amplitudes (allowing to reproduce correlations at all distances, and a possible spontaneous dimerization of the lattice) and $(N,2)$ for the phases (allowing to describe relevant sign patterns with a number of phase parameters linear in $N$, as discussed in Appendix A). Further insight into the accuracy of the 2LR-EPS Ansatz on large lattices can be obtained by comparing our estimated ground state energies for the $J_{1}-J_{2}$ model with DMRG results. As an example, here we focus on a chain of $80$ sites and $\alpha=0.7$, i.e., in the parameter region where the Ansatz is quantitatively less accurate (see Fig. 5). In the mentioned case we find that the 2LR-EPS Ansatz yields a ground state energy characterized by an error relative to the DMRG estimate Vodola of approximately 0.76%. When the 2LR-EPS Ansatz is improved by means of the 10A-EPS one resulting in the 2LR-10A-EPS Ansatz such a relative error decreases to approximately 0.37%.

IV.3.1 Spin-spin correlations in the gapless phase

Throughout the gapless phase $\alpha\leq\alpha_{c}$, $2$-point correlations exhibit a power-law decay as $d(r|N)^{-1}$ [with the chord length $d(r|N)=(N/\pi)\sin(\pi r/N)$]. The latter is the only one compatible with SU(2) symmetry for a gapless Luttinger liquid. Indeed, when mapping the spin model to a hardcore boson chain, the predictions of Luttinger liquid theory for the decay of the off-diagonal correlation function $\langle S^{+}_{i}S^{-}_{i+r}\rangle=\langle b_{i}^{\dagger}b_{i+r}\rangle\sim d(r|N)^{-1/(2K)}$, and of the diagonal correlation function $\langle S^{z}_{i}S^{z}_{i+r}\rangle=\langle(n_{i}-1/2)(n_{i+r}-1/2)\rangle\sim d(r|N)^{-2K}$ Cazalilla (2004) must coincide in the presence of SU(2) symmetry, hence the value $K=1/2$ for the Luttinger exponent. A more detailed analysis Eggert (1996) points out the existence of a multiplicative logarithmic correction to the power-law decay, in the form $\sqrt{\log(r/r_{0})\lambda_{0}}$, where $r_{0}$ and $\lambda_{0}$ are coefficients continuously depending on $\alpha$. Both coefficients vanish as one approaches the critical point with the constraint $-\lambda_{0}\ln r_{0}\to 1$, so that the logarithmic correction disappears in the same limit.

Fig. 7 shows the absolute value of the spin-spin correlation multiplied by the chord length for $\alpha=0$, and $\alpha=\alpha_{c}$, as a function of the square root of the logarithm of the chord length. The expected logarithmic correction is evident for $\alpha=0$, and essentially absent at the transition point: hence the simple 2LR-EPS Ansatz is fully capable of capturing this subtle aspect. Furthermore, a power-law fit of the correlations as a function of the chord length for $\alpha=\alpha_{c}$ (see inset of Fig. 7) leads to a value of the Luttinger parameter in agreement with the theoretical expectation. Indeed, our estimated value of $K$ is $0.497(5)$. For $\alpha>\alpha_{c}$, correlations turn to an exponential decay, albeit with an exponentially divergent correlation length for $\alpha\to\alpha_{c}^{+}$ White and Affleck (1996), introducing very significant finite-size effects.

IV.3.2 Incommensurability in the correlations for $\alpha>1/2$

The incommensurability developing in the correlations for $\alpha>1/2$ is best seen in the structure factor ${\cal S}(k)=\sum_{\alpha}{\cal S}^{\alpha\alpha}(k)$ plotted in Fig. 8. There, one clearly observes that the peak at $\pi$, characteristic of Néel-like correlations dominant up to $\alpha=1/2$, splits into two peaks for $\alpha>1/2$, and the twin peaks move continuously towards the values $\pi/2$ and $3\pi/2$ as $\alpha$ increases. In Fig. 9 we compare the pitch vector of correlations $\phi$ namely the position of the left peak in Fig. 8) obtained in this work with previous DMRG results from Refs. White and Affleck (1996); Gehring et al. (1997) as well as with the MF prediction. The agreement between our estimates and the DMRG ones is acceptable, taking into account that the DMRG results are obtained with open boundary conditions and different system sizes than ours. In general, we consider the successful description of the development of incommensurate helimagnetism at short range as a significant achievement of our Ansatz, considering that it is not accompanied by an obvious breaking of time-reversal invariance. Testing that the phases of all coefficients in the optimized wavefunction are equal modulo $\pi$ is prohibitive when considering the system sizes of interest in this section. Therefore we opt for an alternative test based on spin currents ${\cal J}_{ij}=i\langle S^{+}_{i}S^{-}_{j}-S^{-}_{i}S^{+}_{j}\rangle$, whose average value should be zero on a time-reversal invariant wavefunction, and finite otherwise. Our optimized wavefunctions give a null value of the currents (within the statistical error bar) for all values of $\alpha$, strongly suggesting the real-valuedness of all coefficients. This test is highly non-trivial, as it clearly indicates that the variational optimization exploits the complex nature of coefficients of our Ansatz to introduce incommensurate correlations beyond the reach of a MF-based Ansatz (see Fig. 9 and previous discussion in Sec. IV.1).

IV.3.3 VBC order

The appearance of VBC order for $\alpha>\alpha_{c}$ can be efficiently captured by the dimer order parameter $\mathcal{D}_{N}$, defined in Eq. (15).

We show ${\cal D}_{N}$ for various system sizes as a function of $\alpha$ in Fig. 10, observing that it marks quite clearly the VBC transition. In particular the dimer-dimer correlation function exhibits a decay as $D(r)\sim r^{-1}$ at long distances in the gapless phase, so that we expect that ${\cal D}_{N}\sim N^{-1}$ in this phase: this is indeed observed in the inset of Fig. 10 for $\alpha=\alpha_{c}$, whereas for $\alpha=0.3>\alpha_{c}$ one observes that ${\cal D}_{N}$ extrapolates to a small but finite value.

IV.3.4 Scaling of the entanglement entropy

We conclude our analysis of the 2LR-EPS results for the $J_{1}-J_{2}$ chain with a study of the Rényi entanglement entropy. Considering a total system of size $N$, in the gapless phase $\alpha<\alpha_{c}$ the entanglement entropy of a subsystem of linear size $l$ is expected to scale as

[TABLE]

where again $c=1$, and $c_{1}$ is a non-universal constant Calabrese and Cardy (2004) - extending Eq. (21) to account for finite-size effects. On the other hand in the gapped VBC phase the entanglement entropy saturates to a constant (which takes a value $\approx(c/4)\log(\xi)$ close to the critical point), exhibiting therefore an area law. In particular the correlation length is minimal at the MG point, where $R_{2}(l)\to\log 2$. Fig. 11 shows that all these features are very well captured by the 2LR-EPS Ansatz, and in particular the universality of the central charge throughout the gapless phase.

In order to further contrast the gapless and gapped regime using the entanglement entropy, we concentrate on the difference $R^{N^{\prime}-N}_{2}=R_{2}(N^{\prime}/2)-R_{2}(N/2)\simeq(c/4)\log(N^{\prime}/N)$, which allows one to obtain a finite-size estimate of the central charge as

[TABLE]

while eliminating the size-independent subleading corrections to the dominant logarithmic scaling (in the gapless phase, see Eq. (25)) Läuchli and Kollath (2008). Fig. 12 shows $c_{80,40}$ as a function of the frustration parameter. The transition between the gapless and the gapped phase is clearly signalled by the drop of the central charge estimator from values close to unity to vanishing values. On the other hand the estimator increases again for large $\alpha$, consistent with the fact that, in the limit $\alpha\to\infty$, the $J_{1}-J_{2}$ model reproduces two decoupled Heisenberg chains with $c=1$.

V Conclusions

We have discussed a very flexible variational Ansatz (the long-range entangled-plaquette state or LR-EPS) for strongly interacting quantum lattice models, in which explicit quantum correlations within (clusters of) sites are introduced coupling them at all distances into overlapping plaquettes. This Ansatz generalizes both the well-known Jastrow Ansatz (including two-site correlations) Pang (1991) as well as the adjacent-site EPS Ansatz, and offers the possibility of combining both forms of overlapping plaquettes in the same wave function. Moreover, the use of complex coefficients opens the possibility of reproducing non-trivial sign (phase) structures emerging naturally from the variational optimization without any explicit bias. We have demonstrated the effectiveness of this Ansatz in the case of the frustrated $J_{1}-J_{2}$ quantum spin chain: there we show that it captures both the universal long-wavelength features of the gapless phase of the model (central charge and Luttinger exponent of the corresponding field theory); as well as the appearance of incommensurate short-range correlations at strong frustration.

As it currently stands, our Ansatz can be systematically improved by increasing the size of the plaquettes (either long-range ones or adjacent-site ones) – this improvement strategy is limited, as it introduces a concomitant exponential growth of the variational parameters. However, our findings show that improving the Ansatz is not crucial to capture the long-range physics, which is already accounted for by using minimal (2-site) long-range plaquettes. Further approaches making the Ansatz systematically improvable with a polynomial growth of the parameters are currently under investigation. Additionally, the success of our wave function in reconstructing non-trivial sign structures variationally suggests its potential application to other models of frustrated quantum magnetism, especially focusing on quantum critical points for which a scale-free variational form is particularly well suited. It would also be very tempting to apply this Ansatz to lattice models of strongly correlated fermions, whose nodal surface becomes the object of variational optimization when using complex coefficients.

Acknowledgements

We thank D. Vodola, A. Angelone and G. Misguich for insightful discussions. TR, and FM acknowledge support from the Fédération de Recherche André-Marie Ampère (FRAMA). Computer calculations have been performed on the PSMN cluster of the ENS de Lyon.

Appendix A EPS Ansatz and Marshall-like signs

In the following we show that our simplest Ansatz (namely the 2LR-EPS one) is fully compatible with Marshall-like sign structures, namely sign structures depending uniquely on the parity of the number of $\uparrow$-spins on one of two (arbitrary) sublattices $A$ and $B$ of equal size in which the lattice has been decomposed. A Marshall sign rule with $ABAB...$ sublattice structure is exact in the $J_{2}\to 0$ limit. Yet a Marshall-like sign with a sublattice structure of the kind $ABBA...$ – or with a longer period – is also reproducible with the LR-EPS Ansatz; such structures are indeed relevant for the $J_{1}-J_{2}$ chain under investigation for sufficiently large $J_{2}$, as pointed out in Zeng and Parkinson (1995). It is also important to stress that our Ansatz is compatible with sign structures that are richer than Marshall-like ones.

Marshall-like signs can be represented within the 2LR-EPS Ansatz in several different ways : as the sign/phase of each weight in the wave function results from the sum of the phases of the plaquette coefficients – similar to the relationship between a vector potential and the flux of the resulting magnetic field – different parametrizations of Marshall-like signs can be viewed as different gauge choices.

The simplest scheme representing Marshall-like signs within a 2LR-EPS amounts to considering that the phase of each plaquette, $\theta_{ij}(\sigma_{i},\sigma_{j})$, only depends on the configuration of one spin, and it reads

[TABLE]

if site $i$ belongs to sublattice $A$ or $\theta_{ij}(\sigma_{i},\sigma_{j})=0$ otherwise; here and in the following of this Appendix $\sigma_{i}=\pm 1$ is two times the eigenvalue of the $S_{i}^{z}$ operator, and $N_{i}$ is the number of $(ij)$ plaquettes for which site $i$ is chosen as reference site. This requires therefore the phase of the plaquette to depend on the absolute position of at least one of its sites, leading to ${\sim O}(N)$ parameters. Such a requirement is fully compatible with the choice of modularity of the phases that we made in the study of our largest lattices.

Alternatively Marshall-like signs can be enforced at the level of each pair:

[TABLE]

where $N^{(ij)}_{A,\uparrow}=\frac{\pi}{2}[(\sigma_{i}+1)\delta_{i,A}+(\sigma_{j}+1)\delta_{j,A}]$ is the number of $\uparrow$ spins on the sites of the $(ij)$ plaquette belonging to sublattice $A$ ($\delta_{i,A}=1$ if $i$ belongs to $A$ and 0 otherwise). Then one can easily show that, up to a global phase factor:

[TABLE]

which is precisely the Marshall sign (up to the constant phase factor $e^{i\pi N^{2}/2}$). Here we have used the fact that $N$, being composed of two sublattices, is even, so that $\pi N_{A,\uparrow}N$ is an integer multiple of $2\pi$. We then observe that, in order to reproduce the Marshall-like sign, the phase of the coefficient $e^{i\theta_{ij}}$ must be able to distinguish between $AA$ plaquettes (both sites belonging to $A$) and $BB$ plaquettes (both sides belonging to $B$); and between $AB$ and $BA$ plaquettes. For arbitrary geometries of the $A$ and $B$ sublattices, this requires the phases $\theta_{ij}$ to depend on the absolute positions of both sites, increasing substantially the number of parameters to ${\sim O}(N^{2})$. We have therefore not pursued this parametrization for the largest lattices studied here.

Appendix B Emergence of the sign structure and of real wavefunction coefficients during the optimization

In order to quantitatively exemplify how the 2LR-EPS wave function is able to reproduce the correct ground state sign structure (when the latter is compatible with the parametrization of the Ansatz) we discuss in detail the case of the $J_{1}-J_{2}$ model for $\alpha=0$. Indeed, in this situation, the signs of the coefficients of the exact ground state are known, as they are dictated exactly by the Marshall sign rule, and they can be exactly reproduced by optimizing our Ansatz.

To elucidate how such a sign structure emerges along the optimization procedure we closely examine the evolution with the number of optimization steps of the fraction of wave function coefficients possessing the expected Marshall sign. Hence we define the quantity

[TABLE]

where $N_{C}$ is the total number of ${\bm{\sigma}}$ configurations, $\bm{\sigma}_{\rm{REF}}$ is a reference configuration (here taken as the Néel state $\uparrow\downarrow\uparrow\downarrow...$) , and the term in the sum is 1 (0) if the sign of the coefficient of the generic configuration $\bm{\sigma}$ follows (breaks) the Marshall sign rule.

Fig. 13 clearly shows how the desired sign structure is fully reconstructed (i.e., $\mathcal{M}$ reaches 1) by our optimization procedure within a few thousand steps.

Another important point of our work is that our optimized 2LR-EPS Ansatz is essentially real (up to a global phase factor), recovering a fundamental property of the ground state of the $J_{1}-J_{2}$ chain. This aspect can be elucidated by analyzing the phase relationship between e.g. the wave-function coefficient of all allowed configurations and that of a reference one (here taken again to be the Néel state). If the state is real up to a global phase, then the following quantity

[TABLE]

is vanishing. The evolution of $I(\bm{\sigma},\bm{\sigma}_{\rm{REF}})$ along the optimization procedure is shown in Fig. 14 for the same system parameters of Fig. 13, and for 3 randomly chosen wave function coefficients. $I(\bm{\sigma},\bm{\sigma}_{\rm{REF}})$ is seen to approach a null value (i.e. $<10^{-15}$) after a few thousand optimization steps. The same behavior can be verified for all the wave function coefficients, and this despite the complex form of the initial state.