Parent Hamiltonian Reconstruction of Jastrow-Gutzwiller Wavefunctions

Variational wave functions have been a successful tool to investigate the properties of quantum spin liquids. Finding their parent Hamiltonians is of primary interest for the experimental simulation of these strongly correlated phases, and for gathering additional insights on their stability. In this work, we systematically reconstruct approximate spin-chain parent Hamiltonians for Jastrow-Gutzwiller wave functions, which share several features with quantum spin liquid wave-functions in two dimensions. Firstly, we determine the different phases encoded in the parameter space through their correlation functions and entanglement content. Secondly, we apply a recently proposed entanglement-guided method to reconstruct parent Hamiltonians to these states, which constrains the search to operators describing relativistic low-energy field theories - as expected for deconfined phases of gauge theories relevant to quantum spin liquids. The quality of the results is discussed using different quantities and comparing to exactly known parent Hamiltonians at specific points in parameter space. Our findings provide guiding principles for experimental Hamiltonian engineering of this class of states.


Introduction
Variational wave functions play a key role in the understanding of quantum phases of matter [1][2][3][4][5][6][7][8]. A paradigmatic example is Laughlin wave functions [5], which can be formulated as parametric Jastrow states reproducing several key features of certain fractional quantum Hall effects [9]. Shortly after this, resonating valence bond (RVB) states have been employed as effective descriptions of high-temperature superconductors [6,7,10], and later on, have been linked to fractional quantum Hall physics in Ref. [8]. These early successes boosted variational wave functions as theoretical tools to provide simple pictures for a variety of quantum phases, including topological matter, low-dimensional systems, and tensor networks [11][12][13][14][15].
Perhaps, among these applications, one of the most fruitful has been in the field of quantum spin liquids [16][17][18][19][20][21][22]. These are quantum phases characterized by strong correlations and longrange entanglement among arbitrary far subregions of the system [23], and for these reasons, semi-classical pictures fail in describing the phenomena involved. Variational wave functions have been used to distill generic properties such as correlation functions and entanglement [14].
Interestingly, despite the conceptual simplicity of Jastrow wave functions, it is often challenging to find the corresponding parent Hamiltonians -that is, the Hamiltonians supporting these wave functions as ground states. The major obstruction is that, given a Hamiltonian on a lattice (possibly with frustration terms), quantum fluctuations may cooperate and induce an ordered ground state. This phenomenon is typically referred to as "order-by-disorder" [13]. This problem is of primary importance also due to the latest experimental breakthrough in quantum engineering of synthetic systems [24][25][26][27][28]. In fact, the high degree of interaction tunability of these platforms offers new perspectives and possibilities in otherwise hardly achievable phases of matter, including spin liquids, once parent Hamiltonians are (approximately) identified.
Most of the works in parent Hamiltonian construction studied specific variational states using insightful analytic manipulations [29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44]. Very recently, a series of novel techniques based on systematic approaches have been considered in Ref. [45][46][47][48][49][50]. Indeed, the authors of the latter works introduced new efficient computational algorithms, which remarkably scale polynomially in the system size when restricting the search to local Hamiltonians that have a given initial state as the input eigenstate. To benchmark their techniques, they considered the ground state of some a priori known Hamiltonian as input and checked if the output reconstructed operator coincided with that Hamiltonian. So far, however, there have been no applications of such methods to generic spin liquid variational wave functions, whose parent Hamiltonians are still undetermined.
The present work is the first step in this direction. For concreteness, here we study the class of 1D Jastrow-Gutzwiller variational wave functions [30,51]. These states share two key features with their two-dimensional cousins employed as effective descriptions of quantum spin liquids: they describe extensive superpositions over some (spatially local) state basis, and they have in general as weights analytic functions of the space coordinates. Despite their common appearance, their parent Hamiltonians are not known except for a few fine-tuned cases, amenable to exact solutions. We use an entanglement-guided algorithm presented in Ref. [50] to search local parent Hamiltonians for these states. This method relies on the Bisognano Wichmann theorem [52,53], a quantum field theory result that links systematically the local Hamiltonian density to its ground state reduced density matrix. Its advantage with respect to the other above-mentioned techniques resides in certifying the input state as the ground state of the reconstructed parent Hamiltonian. Indeed, although the methods in Ref. [45][46][47][48][49] are of broader applicability (for instance, they allow for extensions to timedependent problems), they typically certify the ansatz state to be a generic eigenstate, and not the ground state, of the output operator. The main disadvantage is that the method is not applicable in case the wave function cannot be cast as the ground state of Hamiltonian operator supporting low-energy relativistic excitations.
Since the Bisognano-Wichmann technique requires the input state to exhibits relativistic low lying physics, we first investigate the entanglement and correlation properties of these wave functions, identifying a region where the algorithm is expected to perform better. In this regime, we obtain local approximate parent Hamiltonian searching through different algebras of local operators. To check our results, we computed the relative entropy, the correlation functions and the overlap between their ground state and the Jastrow-Gutzwiller wave functions, obtaining fidelities ranging between 95% to over 99%. In addition, we computed the relative error between the ground state energy and the Jastrow-Gutzwiller variational energy of the reconstructed Hamiltonian. In all the considered cases, the relative error is less than 1%, even in the extrapolated thermodynamic limit. We perform systematic searches by increasing both system sizes and interaction range. These results suggest that the exact, yet unknown, parent Hamiltonians of these states exhibit long-range features.
In addition, the method allows us to perform direct parent Hamiltonian searches utilizing simple long-range interactions in the form of monotonous power-law potentials. We find that, while considerably improving the parent Hamiltonian search, such simple long-range interactions are not always sufficiently rich to capture the (unapparent) complexity of Jastrow-Gutzwiller wave functions. These results indicate that the search for exact -albeit long-ranged -parent Hamiltonians for 2D Jastrow-Gutzwiller might be particularly challenging, a fact which is compatible with the scarcity of exact results in this context (with some notable exceptions, see Ref. [29,40]).
The remaining of this paper is structured as follows. In Section 2 we introduce the Jastrow-Gutzwiller states and discuss their physical content through participation spectrum, entanglement entropy and correlation functions. In Section 3 we summarize the Bisognano-Wichmann Ansatz method which we employ in Section 4 to reconstruct various parent Hamiltonians for the above-considered states. The last section is for conclusions and outlooks.

Model wave functions
The Jastrow-Gutzwiller (JG) wave functions are paradigmatic states appearing in several contexts, from integrability to topology (e.g. Laughlin states), to quantum spin liquids. They are characterized by an extensive superposition of spatially local states, and the local weights of the wave functions are captured by polynomials. Throughout this paper, we investigate the one-dimensional case defined on a periodic chain Λ of length L. This setting permits the understanding of finite-volume effects in a systematic manner, as well as enables comparison to exact results.
Let us introduce the wave functions of interest, through the variables n i ∈ {0, 1} defined at each site i ∈ Λ. In the basis {|n 1 n 2 . . . n L }, these states read: Here the sum is over combinations P N {n} constrained by i n i = N . Pictorially, the {n i } variables are occupation numbers of hard-core bosons living on the lattice. The real parameter α and the filling fraction ν = N/L control the properties of the states. For specific combined values of ν and α, conformal field theory calculations have been used to derive exact results pertaining the parent Hamiltonians of these states [42-44, 54, 55]. Throughout this paper, we will consider exclusively the half-filling case ν = 1/2 and L even; the main motivation being that, in spin language, this regime captures both paramagnetic and antiferromagnetic phenomenology. Within this setting, exact results are available only for α ∈ {0, 1, 2}. In Ref. [56], it was proven that α = 0 corresponds to the XXZ chain at ∆ = −1, while the state at α = 2 is the ground state of the Haldane-Shastry Hamiltonian [30,31]. The case α = 1 corresponds to a (symmetrized) Slater determinant, and its parent Hamiltonian is a free fermionic one (up to boundary contributions).

Participation spectrum
To obtain insights for generic values of α, it is instructive to rephrase Eq. (1) in the language of participation spectroscopy [58][59][60][61][62]. This consists of rewriting the wave functions Eq. (1) in a pseudo-energy fashion: In the last equality, we defined the function H α [{n}]: The functional coefficient E 0 [{n}] is an energy constant, while V (i, j) is a logarithmic interaction between occupied particles mediated by chord distances. Thus, we recognize H α [{n}] to be a 2D Coulomb gas (classical) Hamiltonian constrained in a 1D circular lattice [54,57]. Analogously, the wave function normalization Z α is a classical partition function: The parameter α plays the role of temperature and controls the leading weights in the JG states. The modulus squared coefficients in Eq. (1): are Boltzmann weights with classical Hamiltonian 2H α and partition function Z 2 α . The pseudo-energies of 2H α are collectively named participation spectrum and denoted ε({n}).
The ground state ε min determines the larger weights in the sum Eq. (1). For α > 0, the Hamiltonian favors repulsion among particles, constrained by the half-filling condition. Thus, the most probable configurations come from alternating occupation numbers. At negative temperature α < 0 the dominant coefficients are those maximizing the number of occupied nearest neighboring sites. For both cases, such configurations are not unique but degenerate, and for large values of α these states are expected to be the most relevant contributions to the Jastrow-Gutzwiller wave functions. Consequently the JG state are captured by the coherent superposition of these degenerate configurations, which leads, for α 1 and α −1, respectively to antiferromagnetic and ferromagnetic Greenberger-Horne-Zeilinger (GHZ) states [63]: The former state is usually dubbed Néel/anti-Néel state and corresponds to a global Schrödinger cat state. Apart from these extreme limit, at intermediate values of α the system exhibits competing weights, which render analytical arguments demanding if not impossible.
To test this heuristic argument, we consider the gap G = ε min − ε 1 st between the ground state energy of Eq. (4) and its first excited energy, which we refer to as participation gap. Let us discuss the case α > 0. It is convenient to introduce the number of ferromagnetic domain walls as the number of consecutive occupied/unoccupied sites N dws . For example N dws (|010101 ) = 0, while N dws (|011001 ) = 2. The Néel and anti-Néel states, i.e. the most The participation gap G increases linearly with α, with a coefficient that saturates to a constant g ∞ already at modest system sizes. (d) Pseudo-energy differences between two domain walls as a function of the domain-wall separation r. This is a measure of the confining potential between domain walls. The black solid line is the Luttinger liquid prediction [60] with Luttinger parameter K = 1/α. The fit describe extremely well our data, for 4 ≤ L ≤ 64.
probable states, are the only ones with N dws = 0, and all other pseudo-energy excitations can be easily labelled with this number. In Fig. 1 we present the participation spectrum of the JG states for α = 2, 6 and L = 16. The gap G between the most probable and the second most probable state increases linearly with α, with an exactly computable L-dependent constant g L . This saturates a thermodynamic value 1 g ∞ already for modest system sizes.
It is important to emphasize one aspect that is relevant in determining the system properties in the thermodynamic limit. The ground state pseudo-energy with alternating occupied sites is doubly degenerate for every system size. Instead, although the configurations with domain walls are exponentially suppressed in α, their degeneracy scales linearly with system size. In particular at L ∼ exp (cα) for some constant c, we expect a competing and non-trivial behavior between the Néel sector and the first excited sector. This has potentially relevant consequences, which are difficult to predict with the present study. In particular is unclear what effect this pseudo-energy thermodynamics have on quantum observables.
At a practical level, our results are consistent with the intuition above, that the Néel state predominately contributes for large α. In order clearly see the effects of the aforementioned thermodynamic competition for α = 6, we would have needed around L ∼ 10 4 sites. The large 1 We get an analytic expression for the constant: where the ellipsis indicate further computable digits.
gap for any computable finite L considered, renders these excited sectors negligible. The results for α < 0 are analogous to the latter, whereas the most probable configurations are the ferromagnetic ones and the excited pseudo-energy states are obtained as functions of antiferromagnetic domain walls, i.e. number of alternating occupied/unoccupied sites. However the most probable states there are L-degenerate: in the thermodynamics of the Coulomb gas this implies the low-lying pseudo-energy excitation are negligible even at small negative values of α.
Finally, from the substructure of the N dws = 2 sector we can extract how these domain walls interact. In particular, the pseudo-energy difference ∆ε 2dw = ε 2dw (r) − ε 2dw (2) between domains separated by a distance r and those close together (r = 2) has been used for local antiferromagnetic quantum Hamiltonian systems to distinguish between critical and symmetry broken phases of matter. In the former case, the domain walls are logarithmically confined with the separation distance; instead, in the latter this confining is linear. Moreover, the prefactor of this potential for 1D Luttinger liquids [66,67] is related to the Luttinger parameter. This has been tested in Ref. [60], where its authors analyze the XXZ chain.
Because of the explicit form of the classical Hamiltonian density Eq. (6), the interaction between two domain walls is expected to be logarithmic with their separation distance ( Fig. 1, panel (d)). By analogy with the XXZ phenomenology, one is tempted to conclude the JG states are gapless. If furthermore one assumes these states are representatives of Luttinger liquids, the fitted pre-factor suggests a Luttinger parameter K = 1/α. The latter statement has been recently conjectured [64]. This hypothesis is supported by CFT arguments [56] and from studies on the Resta polarization [65]. Here the authors estimate α c = 4 as critical value separating a conducting Luttinger phase to an insulating Néel ordered phase. Our data do not exhibits any transition point in the participation gap, nor a clear distinction between a gapped and a gapless phase. As remarked earlier, this may be due to a finite size effect, which we are not able to resolve at computationally affordable system sizes. In fact it is possible that the N dws = 2 domain walls sector results as decoupled for physical observables of the system, after a critical value of α. At present, however, the consequences of the participation spectroscopy to physical observables are unclear, and further studies are needed in this direction. In the next two subsections we improve our understanding of the Jastrow-Gutzwiller wave functions by numerically studying the entanglement entropy and the correlation functions of the Jastrow-Gutzwiller wave functions. We focus on these properties among others because they serve in the reconstruction technique and its quality checks. The considered system sizes suggests the existence of a critical phase between a Néel and a ferromagnetic GHZ regimes. Using finite size scaling we can bound the former in the interval α ∈ (0, 4.3).

Entanglement entropy
In this subsection, we discuss the entanglement entropy properties of the JG states (for related studies of Rényi entropies in 2D, see Ref. [14]). Entanglement is a fundamental quantity measuring quantum correlations among subregions of the system [68][69][70][71][72][73]. For pure states, this is determined by the spectrum of the reduced density matrix [74,75]. This operator is defined by giving a bipartition of the chain Λ = A ∪Ā and a state |Φ : Here ρ JG is the half-system reduced density matrix of the JG state. The results in green line are obtained through ED using symmetry restrictions. The red lines are the ferromagnetic GHZ predictions for the corresponding system sizes, while the black one is the Néel/anti-Néel cat state entanglement entropy.
Given its spectrum σ(ρ A ), we define the von Neumann entropy by: This function is a bona fide measure of entanglement for pure states when the Hilbert space factorizes in a tensor product form, H = H A ⊗ HĀ, and for this reason is usually referred to as entanglement entropy [76,77]. Fixing A = {1, 2, . . . , L/2}, we compute through exact diagonalization (ED) the von Neumann entropy for the state Eq. (1). We check the GHZ limits by comparing with the analytic calculations for the states in Eq. (9): The agreement is shown in Fig. 2. We isolate an intermediate region between the GHZ regimes by introducing the function: We plot this function in Fig. 3. Within this interval,S vN is logarithmic, with a pre-factor close to 1/3. This is consistent with exact solutions, where the systems display a critical regime. For instance, at α = 1 the system is a linear combination of Slater determinant. At this point the JG state correspond to a free fermion gas and the entanglement entropy can be computed analytically [78,79]: Here c is the central charge (c = 1 for free fermions) and the sub-leading term is a constant. The same scaling holds at α = 2, since the Haldane-Shastry Hamiltonian share the same universality class of the Heisenberg antiferromagnet [30,56]. By continuity, we argue the same critical behavior extends to the whole intermediate region. This is in line with the Luttinger liquid conjecture (see Sec. 2.2). Since the latter is of interest for the subsequent analysis, we estimate its bounding transition points. From Fig. 3 is clear that there is a transition in parameter space at α = 0.
We perform finite-size scaling on our data to estimate the critical value α c of the JG wave functions separating a critical phase with respect to a Néel ordered state. This is a phenomenological finite-size scaling procedure, since it is inherently related to a parameter characterizing the variational wave functions, and not associated to a coupling term in a Hamiltonian. Nevertheless, it is useful to bound the region of validity of the reconstruction method (Sec. 3), which relies on relativistic invariance. We consider the scaled entanglement entropyS(α) as an order parameter, as well as its derivative: which is roughly a susceptibility. We choose to consider both these quantities since the scaling we have is very mild with system size. From Fig. 3, introducing t = log(L) and α = (α − α c )/α c , we use the following simplified scaling ansatz: To perform the finite size scaling we vary the exponents ν, γ and the critical value α c over a suitable range of parameters. The fit is the best over different degrees of polynomials, test with a least-square method against the data [104]. By requiring the exponents to obey scaling relations γ = β − 1/ν we are able to reduce the fitting regime. We estimate the transition at α c = 4.3±0.1 with ν = 2.1±0.2 and β = −0.15±0.3. Value and error bars are the average and standard deviations of the best fits varying the range of system sizes considered. In Fig. 4 we plot both the order parameters of interest and the optimal data collapse. While the quality of the collapses is generically good, the modest system size are not able to resolve more efficiently the exponent landscape, which results quite flat. We believe a more systematic analysis is needed to better characterize the entanglement entropy and its phase transition for the JG wave functions. This would be a useful test also for the Luttinger liquid conjecture in Ref. [65], where it is argued the transition is around the value α conj c = 4. In this paper we choose to follow a more restrictive and cautious approach, focusing on subintervals of α ∈ (0, 4) in the rest of the paper. A concluding remark, which will be useful later, is about the α = 0 point. As previously discussed in the context of participation spectrum, this point is peculiar since the JG state is in an equal-weight combinatorial superposition. Its exact entanglement entropy can be computed [80]: We see that the pre-factor is different from the one in Eq. (15), signal that the state is not representative of the same phase. One can see this by investigating the properties of the exact parent Hamiltonian at α = 0: the XXZ chain at the ferromagnetic transition [56,80]. This Hamiltonian has a gapless quadratic spectrum, thus it breaks relativistic invariance due to a different dynamical exponent 2 z = 2. This observation will be important when trying to reconstruct local Hamiltonians using a relativistic ansatz. Indeed, as we shall comment in Section 4, for α = 0 the algorithm will not be able to return a correct parent Hamiltonian, as expected.

Correlation functions
To further characterize and resolve the Jastrow-Gutzwiller states, we compute the one-body and two-body spin correlation functions {σ z , σ + , σ − }. Their scaling properties resolve the nature of the state being critical or not. Due to the binary nature of the n i variables, for notational convenience we introduce the unary-not operator F ij acting on the site i, j, whose action on basis state is defined by logical negation on n i and n j . Since the system exhibits a U (1) symmetry related to number conservation, we compute only U (1) invariant correlation functions. Recalling σ z = 2n − 1 with n the number operator we have: At half-filling the first one is identically zero. The latter ones can be easily implemented numerically. The correlation length can be extrapolated through finite size scaling of the connected correlation function σ z i σ z i+L/2 c : Here a is a constant, while γ characterize the algebraic decay. In all the above equations, we exploited periodic boundary conditions. Let us stress that the definition Eq. (22) is meaningful only when the cluster decomposition principle holds. This requires the connected correlation function to decay to zero with the distance between the spins. This definition is used throughout in the literature of critical phenomena, where the phase is defined through the ground state manifold of specific Hamiltonians [82]. In this context, symmetry broken phases at finite system size manifest themselves as a coherent superposition of the ground states in the different symmetry sectors (GHZ states) [23]. The latter are a remarkable example of states which do not respect the cluster decomposition. To avoid odd/even effects, we present only L multiples of four.
Having the above remark in mind, we consider the definition Eq. (22) to also characterize the parameter space of the JG wave functions. Here we first check the system is fulfilling the cluster decomposition principle condition. When this is not the case, we expect the JG state to be representative of a finite size symmetry broken phase. Within this setting, if the parameter ξ is finite, the exponential behavior dominates on the algebraic one and the system is gapped, while if ξ → ∞ the system behaves as critical.
In Fig. 5, we show the results of our fitting procedure, plotting the inverse correlation length versus 1/L. For the chain lengths considered, the thermodynamic limit is difficult to estimate since at finite size the inverse correlation length 1/ξ L can be trusted upon the value 1/L. However, all values α < 4.0 are compatible with an infinite correlation length.
For large positive values and negative values of α, the cluster decomposition principle fails. The corresponding GHZ states (introduced in Sec. 2.2), representatives of symmetry broken phases, are confirmed to reproduce the correlation functions of the JG wave functions. A detailed discussion is given in Appendix A.

Entanglement guided search for parent Hamiltonians
In this section we summarize the scheme we employ to reconstruct parent Hamiltonians [50]. As previously remarked, this method requires additional conditions to work. This in contrast to other techniques [45,46] based on the quantum covariance matrix (QCM). The latter are simpler to implement since are based on requiring the input state to satisfy the zero energy variance condition. Thus, those methods generically guarantee that the input state is an eigenstate (not the ground state) of the parent Hamiltonian. Here comes the reason we have chosen to use the Bisognano-Wichmann Ansatz (BWA) scheme: the additional physical constraints guarantee the parent Hamiltonian of the input state as the ground state. This condition is at the core of eventual simulation protocols, since excited state are less robust in analogue experiments. Nevertheless, the relativistic requirement can be applied only to a narrow number of settings: for example if non-translational system are considered, such as disordered systems, BWA fails while QCM still gives meaningful results [103], provided a a fortiori analysis is done on the parent Hamiltonian space and their spectra. The method we adopt is based on the Bisognano Wichmann (BW) theorem, which for convenience we recap in the first subsection. Then, we introduce the common ingredients shared with other aformentioned techniques [45][46][47]49,50]. We conclude this section by presenting the algorithm and our chosen implementation.

Bisognano-Wichmann theorem and lattice models
By definition, reduced density matrices are positive operators with bounded spectrum σ(ρ A ) ⊂ [0, 1]. Consequently, it is always possible to find a lower bounded operator K A such that ρ A ∼ exp(−K A ). This object is usually referred to as entanglement or modular Hamiltonian, and in general is highly non-local, being the logarithm of the non-local operator ρ A .
Remarkably, Bisognano and Wichmann proved that the entanglement Hamiltonian acquire a local density when considering the ground state of a relativistic quantum field theory partitioned into two half-spaces [52,53,77,81]. Moreover, the density of this modular operator is proportional to the one of the theory Hamiltonian. The statement is the following.
Theorem (Bisognano Wichmann) Given a local relativistic QFT in d + 1 spacetime dimensions, described by an Hamiltonian H = d d xH(x) the half-space reduced density matrix of the vacuum |Ω is: Here A and B are respectively the manifolds A = {x ∈ R d : x 1 ≥ 0} and its complementary, while v is the sound velocity of the relativistic excitations. Sometimes, the pre-factor β ≡ 2π/v is dubbed entanglement temperature due to the analogy with respect to thermal density matrices.
More recently, this result has been revisited in the context of holography and many-body physics [83][84][85][86][87][88][89][90][91][92][93][94]. In particular, the theorem has been extended for theories with conformal invariance [83,85,86]. Given the subsystem A = {x ∈ R d |0 ≤ r ≤ R, r = ||x||}, its entanglement Hamiltonian reads: Interestingly, when considering lattice systems exhibiting relativistic low-lying excitations, the discretisation of Eq. (23) and Eq. (25) gives a fine approximation of their reduced density matrices [78,93,[95][96][97][98][99][100], with even exact results for specific models [101,102]. Moreover, the discrepancies due to the lattice structure disappear in the thermodynamic limit. This motivates the core idea behind the BWA method: to find optimal BW entanglement Hamiltonian describing the reduced density matrix of state of interest, in our case the Jastrow-Gutzwiller wave functions. For concreteness, in the remaining of this paper we make use of the discrete version of Eq. (25) in 1D system of size L and A = {1, 2, . . . , L/2}: Here r label the sites, h r is the lattice density of the Hamiltonian H, while K A the corresponding modular operator. Conventionally, we chose to absorb the entanglement temperature in the Hamiltonian density couplings h r .

Basis of local operators
To quantitatively describe the theory and entanglement Hamiltonians on the lattice we introduce the basis of local operators. As previously mentioned, these fully characterize the operator space of the parent Hamiltonian search. We say an operator is k-local if either (1) it has finite domain k-nearby few body operators, or (2) it is written as a linear combination of the latter. Furthermore, we require k to be constant for any finite system size L we consider. If these conditions are not fulfilled, we say the operator is non-local.
We define a basis of k-local operators as the set of matrices {O µ,r } µ∈I,r∈Γ . Here I is a set of internal indices, while Γ ⊂ Λ is a set of sub-lattice ones. Depending on the values of I and Γ, these basis span different vector spaces of local operators, whose generic element is: The dimension of these spaces is thus given by the combined cardinality of the label sets D = |I||Γ|. Before moving on, we clarify the above notation through few examples. Let us first consider the Pauli algebra at each site r ∈ Γ = Λ: The generic linear combination is: We see the total dimension is D = 4L in this case. A less trivial example is the two-body nearest neighboring interactions: Here α covers, in addition to the elements in Eq. (31), the following two-body operators at each site r: O 4,r = σ x r σ x r+1 , O 5,r = σ x r σ y r+1 , . . . , O 10,r = σ z r σ y r+1 , O 11,r = σ z r σ z r+1 .
The linear space has dimension D = 12L. Imposing symmetries one can reduce the dimension D of the operator space, in the same fashion symmetry constraints can be used to block diagonalize observables. For example, imposing U (1) and translational symmetry, a possible operator basis is the following: Here, the index α takes three values (D = 3) and the Hamiltonian is: In the second step of the above equation, we wrote the operators h α in terms of Eq. (31). Thus, the freedom of choosing the operator basis enables us to specify the required symmetries of the parent Hamiltonian, and it allows a reduction of complexity (for translational invariant systems, D ∼ O(1) in system size). Motivated by the symmetries of the JG states, we will consider the following basis for k ≥ 2: Varying the value of k we consider an increasing number of nearest-neighboring hopping and exchange operators. Finally, since the physics of the JG state at α = 2 is captured by a long range model, we shall consider the basis of non-local operators: These basis are both U (1) and translationally invariant, thus exhibits coefficients w α not depending on lattice sites. In literature, non-translational invariant basis have been employed in the reconstruction of disorder system Hamiltonians [45,48,103], or to enlarge the set of Hamiltonians having the input state as an eigenstate [46].

Parent Hamiltonian reconstruction method
We are now in position to present the BWA scheme. Let ρ input A be the half-system reduced density matrix of the the input state. We want to find optimal coefficients w α in Eq. (34) such that: This optimization can be implemented using any estimator of distance between ρ input A and the model reduced density matrix ρ BW A ({w α }). For example one can use the Kullback-Leibler divergence between the participation spectra of the reduced density matrices [60]. This estimator has the advantage of being easy to implement even for larger spacetime dimensions, but has the drawback of leading in general to a non-convex optimization. Such obstacle can be anyway surpassed using stochastic optimization algorithms. Instead, for the class of models described by the basis in Eq. (35) and Eq. (36), it can be proven that any convex estimator acting on the space of density matrices leads to a convex optimization problem (with a unique solution). Among these, we have found particularly useful for numerical implementations the relative entropy, which we adopt in the remaining of this paper. Given two density operators ρ and σ, it is defined as: This function quantifies the distance between between ρ and σ, it is non-negative S(ρ|σ) ≥ 0 (with the equality holding only if ρ = σ) and it is jointly convex. In particular, its restriction to a single argument is a convex function. As already stated, the relative entropy leads to a convex optimization admitting, up to numerical precision, a unique solution [50]: The relative entropy value express a "distance" in the reduced density matrix manifold, and quantify the difference between the initial wave function and the closer one fulfilling the BW theorem.
We implement a gradient descent on the relative entropy. Introducing the notation ∂ α = ∂/∂w α and: the gradient of the relative entropy reads We remark that the actual input needed are just the expectation values over the ground state and over the "thermal" BW density matrix. The former can be sometimes computed analytically, as in the JG states (see Section. 2), while the latter can be implemented with different numerical methods, including quantum Monte Carlo when no sign problem is present.

Reconstruction of Jastrow-Gutzwiller parent Hamiltonians
In this section, we apply the entanglement based reconstruction technique to JG wave functions, considering different choices for the operator basis. We quantify the quality of the reconstruction utilizing (1) relative entropies between reduced density matrices, (2) wave function overlaps, and (3) correlation functions. In view of the discussion in section Sec. 2, we focus here on the regime 0 < α < 4; the regimes where the wave functions are captured by GHZ states are instead discussed in Appendix A.

Models for reconstruction
We consider two paradigmatic classes of operators as candidates for the parent Hamiltonian reconstruction. The first one are the k-local Hamiltonians constructed from the basis B N N (k) These Hamiltonians for k ≤ 4 are archetypal for the study of strongly correlated matter in 1D and 2D, and have been used for ab initio numerical studies of quantum spin liquid phases in different lattices [18-21, 51, 104]. We notice that these operators contains the XXZ and the J 1 − J 2 model as particular cases. The second class are long-range XXZ Hamiltonians constructed from the basis B LR in Eq. (36): The reason in the latter choice is twofold: on one hand J 1 = ∆ 1 is the Haldane-Shastry Hamiltonian, the exact parent Hamiltonian at α = 2. On the other hand, in Ref. [54] Shastry conjectured that α = 2 is the ground state of Eq. (43). We remark that the parent Hamiltonian is defined up to an overall multiplicative constant which sets the energy scales, and an additive zero energy value. Thus, without loss of generality, we factor out the J 1 term and we are interested in the values {w/J 1 }.

Numerical implementation
We search parent Hamiltonians of the above form through the BWA technique. The implementation is based on exact diagonalization (ED) routines in Fortran, using standard libraries and LAPACK [105]. We performed gradient descents with Figure 7: Relative entropy of the JG reduced density matrix and the BW converged one. We see that enlarging the domain of the operator involved, the quality of the results increases. The line α = 1 corresponds to a free fermions gas.
various threshold error th = 10 −3 − 10 −6 . In the considered region, we notice no qualitative change in the observable behavior, although a smaller threshold error requires more steps in the gradient descent convergence. For convenience, we present the results only for th = 10 −4 . At this value, the observables are determined with a precision of around 0.1%. The initial value of the couplings is drawn by a uniform random distribution on the interval [−2, 2]. Here the spreading plays a minimal role: since the optimal solution is unique (see Section 3), the only ambiguity is numerical and due to the truncation to th . The resulting uncertainty is in the last sensible digit of the relative entropy and of the other observables, which we lift through averaging over 50 initial configurations. As argued in Sec. 2, in the thermodynamic limit the system should exhibit a critical regime in the region α ∈ (0, 4.30). However, for the modest values considered L ∈ {4, 6, . . . , 20}, we chose to focus on the subregion α ∈ (0, 4), where finite-size effects are less severe.

Diagnostics for reconstruction
Let us introduce the observables we use to access the quality of the parent Hamiltonian reconstruction. Firstly, we evaluate the relative entropy S(ρ jas |ρ BW ) between the converged BW reduced density matrix ρ BW and the exact JG one ρ jas . Since this function is a "distance" in the density matrix space, it quantifies how much the BW density matrix approximates the input state.
We then introduce the module of the overlap | ψ jas |ψ rec | between the JG wave function |ψ jas and the ground state of the reconstructed Hamiltonian: We stress that this quantity is meaningful only for finite size systems, since it decays to zero in the thermodynamic limit, for any arbitrary small difference between two vector states (in analogy with orthogonality catastrophe [106]). Finally, we compute the following quantity, a cumulative estimate of how much the correlation functions over the reconstructed state differ from the exact ones: Here the first term is the correlation function respectively on the ground state of the reconstructed parent Hamiltonian and on the JG state eq (20). The 1/ (L) factor renders this object non-extensive, which is desirable when comparing different system sizes. For convenience, we call this operator the cumulative correlation difference. Equipped with these tools, in the following subsections we separately present the analysis for the previously introduced basis Eq. (42) and Eq. (43). On the former, we first discuss overlaps and relative entropies for different basis choices, and finally discuss correlation functions. On the latter, we focus the analysis only on the relative entropy.

Reconstruction with N N (k)
We begin by considering the models in Eq. (42) for k = 2, 3, 4. If a p-local Hamiltonian exists, we expect the terms k > p to be finite size terms and to decay to zero enlarging the system size. We anticipate that our result suggests that an exact local parent Hamiltonian exists only for α = 1 (see, e.g., the scaling of the overlap depicted in Fig. 8), which corresponds to free fermions 2-local Hamiltonian. At different values of α, the reconstruction is only approximate, although it improves considerably increasing the basis N N (k). We deduce that the exact parent Hamiltonian should involve long-range interactions.  Search for nearest-neighbor Hamiltonians. -Let us first restrict the easiest setting, that is choosing the N N (2) basis. In this case, the Hamiltonian Eq. (42) corresponds to the XXZ model. The value of interest is ∆ 1 /J 1 . When this is zero, the model reduces to the XX chain, which is a free fermion model up to a Jordan Wigner transformation. Moreover, it is interesting to compare our results with those of Ref. [56]. There, the authors considered the inverse variational problem, optimizing the parameter α with respect to the fixed ratio of ∆ 1 /J 1 . They argue that for α ∈ [0, 2] the wave functions are representatives of the critical phase ∆ 1 /J 1 ∈ [−1, 1] characterizing the spin-1/2 XXZ chain. Our results are compatible with their findings and the analytic results (Fig. 6).
For larger values of α, our results still indicate a very clear convergence to the thermodynamic limit. Moreover, the extrapolated values ( Table 1) always indicate that ∆ 1 > J 1 in this regime: this is compatible with an antiferromagnetic state with a very large correlation length. This finding is highly non-trivial, as there is no guarantee that our method shall return the correct parent Hamiltonian even in the presence of strong finite-volume effects, that have to be expected in this regime since, in the XXZ model, the transition to an antiferromagnetic phase belongs to the Berezinskii-Kosterlitz-Thouless universality class.
Search beyond nearest-neighbor Hamiltonians. -It is important to test the stability of these findings both with respect to enlarging the basis, considering N N (k > 2), and to system size. We thus considered the reconstruction also N N (3) and N N (4), and studied the behavior of the couplings {w α /J 1 }. As shown in Fig. 7 and Fig. 8, both the relative entropy and the overlap improve including higher-k terms. In addition, the magnitude of the couplings corresponding to the latter seems to increase with system size (see Fig. 9), suggesting that the exact Hamiltonians for the Jastrow-Gutzwiller states are long-ranged. An exception is the point α = 1, whose reconstructed Hamiltonian converges to the XX chain. As argued in Sec. 2, this is expected due to analytic arguments.
Ferromagnetic JG wave function. -Another particular point is α = 0. There, the corresponding JG wave function is the exact ground state of the ferromagnetic transition point XXZ. The BWA in principle should not work being this point described by a nonrelativistic field theory [80]. However, the converged coupling is flowing toward the correct ∆ 1 /J 1 = −1 enlarging the system size. Importantly, this result is strongly dependent on the basis chosen, and we see that it is unstable adding larger hopping terms (N N (3) and N N (4)).
Here the modulus of the couplings corresponding to (k > 2)-local terms increases, signal that a relativistic exact parent Hamiltonian for this point, if it exists, it is strongly long-range.
Correlation functions. -Finally, we present in Fig. 11 the results for the cumulative correlation difference V (rec|jas). At fixed system size L, it slightly increases when including higher k-terms. This is counterintuitive, since we observe that a larger basis N N (k) leads to states that are more similar to the JG wave functions (see Fig. 7 and Fig. 8). With the present analysis, we are not able to fully characterize if this trend is due to finite size effects or it has a more systematic nature. A possible explanation would be hidden in the BWA algorithm: since it optimizes over the short-k correlations (see Eq. (41)), the large distance correlators are less controlled and are subject to frustration effects. Within this interpretation, these discrepancies may suggest that longer range terms are required in the optimization to faithfully reconstruct an exact parent Hamiltonian.
Instead, at a fixed value of k, the cumulative correlation difference seems to saturate at some finite value. Being such an object deviation measure from a standard value (see Eq. (45)), it roughly gives how much on percentage the correlation functions change at a fixed site. In the worse scenario of our results, this has a value of around 10%. One may compare our findings with the exact results of the Haldane-Shastry model and the antiferromagnetic Heisenberg chain [30,31,107]: From the latter equations, we read the relative error of the nearest neighboring correlators and next-nearest neighboring ones, respectively of 2% and of 8%.
Combining the above reasonings, we state the reconstructed parent Hamiltonians are only approximate and the true parent Hamiltonians for the JG states require non-local terms. This further confirms our previous analysis. The exception is the point α = 1, where the cumulative correlation difference improves both with system size and by including larger N N (k). Figure 9: Ratios of the converged couplings ∆ 2 /J 1 and J 2 /J 1 versus inverse system size. We see that α = 1 is flowing toward the XX Hamiltonian (see also Fig. 6), while the other converged values are stationary in non-null values, suggesting long range 2-body physics for the JG states. Figure 10: Ratios of the converged couplings ratios versus inverse system size using the N N (4) basis. As for the previous cases, we see that α = 1 is flowing toward the XX Hamiltonian (see also Fig. 6,Fig. 9). Other values of α suggest a long range 2-body physics for the JG states. Relative error of the variational energy -As a last check we compute the variational energy of the parent Hamiltonian with respect to the Jastrow-Gutzwiller input state: and compare with the exact ground state energy E gs . The results are quantitatively compared via the relative error: We present the our results in figure Fig. 12. At fixed value of N N (k), our data suggest a mild linear growth of the relative error with system size. A linear extrapolation of the thermodynamic limit is given. All the considered cases lie within 1% of relative error in the energy landscape. Interestingly, at fixed L the relative error increases including larger N N (k), in a similar fashion to what we observe in the correlation functions. At present we cannot fully understand and characterize such counterintuitive behavior. As already mentioned in the previous paragraph, this may be due to the algorithm forcing the optimization on a finite size landscape and creating frustration effects. The latter likely explain the case α = 2, which should converge to the Haldane-Shastry pre-factors. Another possibility is that a new operator content is needed, and the chosen basis cannot grasp the thermodynamic properties of the systems. Further investigations on this problem are left for future studies.

Reconstruction with the long range model
We investigate the reconstruction when considering the model Hamiltonian Eq. (43), limiting our discussion to the relative entropy detector (see Sec. 4.2). The couplings are reported in Fig. 6, compared with the N N (k) cases. For the chain lengths considered, only at α = 2 the relative entropy shows a decreasing trend with system size (Fig. 13). This indeed corresponds to the exact Haldane-Shastry parent Hamiltonian. However, except at this fine-tuned point, At present we cannot infere a clear thermodynamic behavior for k → ∞, as the number of points is too modest to fit. Nonetheless, the error seems bounded within 1%-2%.
the relative entropy grows with system size, suggesting the parent Hamiltonian Eq. (43) is no the exact parent Hamiltonian for α = 0, and other more intricated terms must be added. Figure 13: Relative entropy between converged BW density matrix and the JG one for the long range model Eq. (43) for L = 8, 12, 16, 20. The results show a decreasing relative entropy for α = 2, which suggests the algorithm is approaching thermodynamic convergence. Instead, even points close to this Haldane-Shastry point exhibits increasing entropy, and certifying only an approximate reconstruction.

Conclusion and outlooks
In this work, we reconstructed approximate parent Hamiltonian for the one-dimensional Jastrow-Gutzwiller wave functions. We identified a region in parameter space where these wave functions display critical properties. Outside this interval, they are effectively described by Schrödinger cat states. Most likely, they are representatives of symmetry broken phases and their parent Hamiltonian is classical and constrained by the half-filling condition on the states. For the reconstruction technique, first we considered k-local Hamiltonians. We confirm the exact point α = 1 corresponding to free fermions, obtaining the XX Hamiltonian. At α = 0 the method fails to find local and relativistic parent Hamiltonians. This is due to a breakdown in the relativistic invariance in the wave function, whose exact parent Hamiltonian manifest gapless quadratic spectrum [56,80].
Our findings suggest the exact parent Hamiltonian for α = 1 should involve more complicated U(1)-invariant interactions, potentially with larger support. We checked the hypothesis of Shastry (Ref. [54]) of considering long-range XXZ chains with square secant couplings. Up to the considered system size there is a slow trend toward larger relative entropy, thus suggesting the ansatz is likely to be insufficient. Nevertheless, finite-size results are of value for Hamiltonian engineering and quantum simulations. Indeed, the BWA method provides inherently finite-size optimization and control on the basis chosen and on the quality of the outputs. In particular one can choose experimentally suitable operators in the basis, such as two-body operators. The fact that our technique is easily adaptable to include fully-longranged interactions may also be used in a different manner, that is, to certify and validate quantum simulators aimed at finding ground states of spin models including slowly-decaying power-law interactions, which are realized in both trapped ions [27] and Rydberg atom experiments [25,108].
It is of primary interest to apply similar techniques and considerations to two dimensional wave functions, such as the Laughlin wave functions. In fact, being the only computational demanding part of the algorithm the calculation of the ground state and the Bisognano-Wichmann expectation values, in principle one can tackle also higher dimensions by using Monte Carlo techniques. From the quantum engineering viewpoint, another intriguing perspective is to search for Liouvillians that have Jastrow-Gutzwiller wave functions as unique steady states [109,110]. In particular, dissipation may considerably soften the requirement for long-range couplings thanks to correlations induced by the bath.

A Correlation functions and parent Hamiltonian for the GHZ regimes
We argued that the JG states at α < 0 and α 1 corresponds to ferromagnetic and antiferromagnetic cat states. A first check is given by means of the participation spectrum and of the entanglement entropy (see Fig 2 in Section 2). Given the simple form of these GHZ states Eq. (9), we can compute their analytic correlation functions: In Fig. 14 we check the agreement between the above equations and the numeric correlation functions computed on the exact JG states. Our results suggest the state is in a symmetry broken phase [23]. Intuitively, we can guess classical parent Hamiltonians having these states as the ground state. For example, a ferro/antiferro-magnetic Ising model with the constraint of having zero magnetization. In practice, one can represent these states as MPS and use well-known results [11,23] to reconstruct local parent Hamiltonians. Figure 14: Difference between numerical correlation functions computed on the JG states and the analytic formulae Eq. (49). The different system sizes show a scaling to zero, confirming the correctness of the GHZ limit.