Variational wavefunctions for Sachdev-Ye-Kitaev models

Given a class of $q$-local Hamiltonians, is it possible to find a simple variational state whose energy is a finite fraction of the ground state energy in the thermodynamic limit? Whereas product states often provide an affirmative answer in the case of bosonic (or qubit) models, we show that Gaussian states fail dramatically in the fermionic case, like for the Sachdev-Ye-Kitaev (SYK) models. This prompts us to propose a new class of wavefunctions for SYK models inspired by the variational coupled cluster algorithm. We introduce a static ("0+0D") large-$N$ field theory to study the energy, two-point correlators, and entanglement properties of these states. Most importantly, we demonstrate a finite disorder-averaged approximation ratio of $r \approx 0.62$ between the variational and ground state energy of SYK for $q=4$. Moreover, the variational states provide an exact description of spontaneous symmetry breaking in a related two-flavor SYK model.

Introduction.-Variational wavefunctions are at the heart of our understanding of a variety of condensed matter systems like quantum Hall systems [1], superconductors [2], and correlated metals [3]. These wavefunctions provide an intuitive description of these phases, and are often useful for numerics. Working with pure states also makes it possible to study entanglement, a property which has been crucial to characterize exotic phases of matter [4]. Further, with the advent of quantum simulators [5,6], and in particular of hybrid quantumclassical variational algorithms [7,8], it is desirable to find preparable states that can reach low energy regimes of strongly correlated Hamiltonians.
A related topic of recent interest is Hamiltonian complexity [9], which studies the computational complexity of approximating the ground state of certain classes of Hamiltonians. These problems belong to the quantum Merlin-Arthur (QMA) class since a verifier can check a solution (i.e. a quantum state) efficiently on a quantum computer by measuring its energy [10][11][12]. Whereas approximating the ground state energy within a small additive error was shown to be QMA-complete for a wide range of Hamiltonians, the complexity of approximating the ground state energy density within finite relative error is still undecided, and is closely related to the quantum PCP [13-17] and NLTS conjectures [18]. Proving these conjectures would, roughly speaking, require finding classes of Hamiltonians for which not only the ground state but all states below a finite energy density are impossible to reach with a simple ansatz.
Given a class of traceless Hamiltonians and a class of ansatz wavefunctions, one can define a figure of merit called approximation ratio, given by r ψ ≡ E ψ /E GS , where E GS is the energy of the ground state, and E ψ = min ψ ψ| H |ψ , where ψ belongs to the class of ansatz wavefunctions. For non-trivial Hamiltonians, simple wavefunctions (e.g. product states) are of course not expected to reach an approximation ratio very close to 1. The question we aim to answer instead is whether they * arijit.haldar@utoronto.ca can at least achieve r ψ > 0 in the thermodynamic limit.
Our work is motivated by the following question: can similar results be obtained for q-local fermionic Hamiltonians [21]? For fermionic systems, a natural analog of product states are Gaussian states, which include the Slater determinants calculated with Hartree-Fock. However, for q > 2, we will provide strong evidence that the approximation ratio of Gaussian states goes to 0 in the thermodynamic limit: r Gauss → 0 for N → ∞. This highlights a fundamental difference between the bosonic and fermionic case. It also motivates the following question: if Gaussian states are not up to the task, is there any other class of tractable wavefunctions that could provide a finite approximation ratio?
Model and ansatz.-The q-SYK model is defined as: with i, j ∈ [1, N ] and with g = (q/2)!/ (q/2)( N 2 ) q 2 − 1 2 . The symbolsĉ † i ,ĉ i denote fermionic creation, annihilation operators. The couplings J i1···i (q/2) ;j1···j (q/2) are Gaussian random numbers which satisfy appropriate symmetrization conditions [51]. The variance is represented as J 2 , and will be set to one except when written explicitly. This Hamiltonian has an extensive energy bandwidth which is symmetric around zero due to particle-hole symmetry [52].
The simplest variational wavefunctions for a fermionic model are Gaussian states (which include Slater determinants), and the corresponding optimization procedure is the celebrated Hartee-Fock [53]. In quantum chemistry, this technique typically recovers 99% of the electronic energy, and is the basis for a variety of more sophisticated approaches. By contrast, for SYK models with q > 2, an elementary calculation (see [54], SM. S1)) shows that the energy bandwidth of Gaussian states, which is also symmetric around zero, scales subextensively with N . This means that, in the large-N limit, Gaussian states only reach a vanishing fraction of the full many-body bandwidth, and are effectively at infinite temperature. This is a strong indication that the approximation ratio of Gaussian states for q-local fermionic Hamiltonians with q > 2 goes to 0 in the large-N limit, in contradiction to the conjecture found in Ref. [21]. Intuitively, this happens since minimizing the energy requires optimizing over the value of q-point correlators, but these correlators are over constrained for a Gaussian state: due to Wick's theorem, all higher-order correlators are simple functions of two-point correlators.
Since Hartee-Fock does not produce any useful result, we take a different approach: let us look for a subset of terms in H which commute with each other, and for which the energy can be minimized easily. The selected subset of terms should be extensive in order for the state to have a finite approximation ratio, i.e. it should contain a number of terms which scales as N q . We propose to construct such a set by partitioning the system into two subsystems (see fig. 1), with N L sites on the left and N R = N − N L sites on the right, and by keeping only terms with creation operators on the left side, and annihilation operators on the right side: The N orbitals are partitioned into left, right subsystems. The operatorT † moves (q/2)-fermions (with q = 4 in the figure) at a time from the right side to the left side, with the same amplitude Ji 1 ···i (q/2) ;j 1 ···j (q/2) as the corresponding term in the Hamiltonian. Starting from a state |0 in which the right-side is filled with fermions and the left is empty, the variational wavefunction is constructed by repeated applications ofT † .
where L = [1, . . . , N L ] and R = [N L + 1, . . . , N ]. The parameter p = N R /N L can be tuned at will, but we will focus on p = 1 for now. It will be useful to define the partitioned-SYK Hamiltonian, which contains an extensive subset of the terms of H SY K , and which is an example of the systems studied in Ref. 35. Using this notation, the ansatz wavefunction is defined as where |0 is the state for which all states on the right (resp. left) are full (resp. empty), a is a real variational parameter, and where the normalization is given by N (a) = 0 | exp(−aT ) exp(−aT † )|0 . The intuition behind this state is as follows: starting from a state that is empty on the left and fully occupied on the right, we create a population of particles on the left and holes on the right by applying the corresponding terms from the Hamiltonian. Interestingly, this wavefunction belongs to the class of variational coupled cluster (VCC) states developed for quantum chemistry [55][56][57][58]. This algorithm has the advantage of being variational (as opposed to regular coupled cluster [59,60]), but is usually limited to a very small number of orbitals due to the factorial complexity of the method. By contrast, we were able to perform VCC directly in the large-N limit for a class of SYK models.
The disorder-averaged energy density for the state is given by E(a) = 1 N ψ(a)|H SY K |ψ(a) and can be calculated using where we used the fact that the expectation value of the terms which are present in H SY K but not in H pSY K vanishes after disorder averaging.
Large-N theory.-To enable the computation of log(N ), we introduce a field-theoretic approach similar to the fermionic path integral(see SM. S2 for details). First, we perform a particle-hole transformation on the right side, wherebyĉ i∈R =ĥ † i∈R andĉ † i∈R =ĥ i∈R . We then define the fermionic-coherent states |c i∈L , |h i∈R for left and right, characterized by the Grassmann numbers c i , The disorder averaging is implemented using the replica-trick log(N ) = lim R→0 [N R − 1]/R. This results in a "static" action involving Grassmann numbers c i , h i with no imaginary time dynamics. Introducing the static Green's functions along with the self-energies Σ c , Σ h , into the action, allows us to integrate the fermions c i , h i . The particle densities in the left and right subsystems are simply given by For p = 1, particle conservation implies ρ L +ρ R = 1, and thus G c = G h ≡ G and Σ c = Σ h = Σ. At the saddle point, one finds which are polynomial equations for G(a) and Σ(a) that can easily be solved numerically. These relations derive from the generating function log(N ), which takes the form − log(N ) = −N log(1 + Σ) + ΣG + a 2 J 2 q G q , (9) at the large-N saddle-point. Interestingly, this generating functional can be interpreted as a static limit of the free-energy for the SYK model [25, 28], given by: where τ ∈ [0, β] denotes the imaginary-time variable and β = T −1 is the inverse temperature. Indeed, if the imaginary time dynamics is eliminated by substituting ∂ τ → 1, G(τ ) → G and Σ(τ ) → Σ, an expression similar to −log(N ) in eqn. 9 is recovered. The energy density E(a) is calculated using eqn. 6 to give where G is obtained by solving the saddle point equations (see fig. 2). The most important point is that E(a) does not decay with N , which means the variational states have an extensive bandwidth, and thus a finite approximation ratio in the large-N limit. The variational energy has a single minimum as a function of a (see fig. 2 bottom), with the following properties: We can now compare E min with the energy density of the ground-state of the SYK model (E SY K ). The latter can be obtained by taking the zero temperature limit of F SY K (see eqn. 10). We give a comparison as a function of q in Fig. 3(a). For example, for q = 4, we find E min = 8/25 √ 5 ≈ −0.143 and E SY K ≈ −0.2295. Since we expect both E min and E SY K to be self-averaging, we define the disorder-averaged approximation ratio as r ψ = E min /E SY K . We thus find r ψ ≈ 0.62 for q = 4. To put things into perspective, we have calculated that E min has the same energy density as the thermal ensemble of the SYK model at temperature T /J ≈ 0.455.
A peculiarity of |ψ(a) is that it does not in general have a uniform density of particles in the left and right subsystems, in contrast to the original SYK model for which all orbitals are at half-filling. As seen in fig. 2, the densities ρ L and ρ R evolve monotonically with a, and are only equal at a = a s = 2 (q/2)−1/2 . This discrepancy arises from the fact that our construction aims at minimizing H pSY K , which contains only a subset of the terms in H SY K , and creates an artificial distinction between the two subsystems. The Hamiltonian H pSY K is actually interesting in its own right as it can be understood as an example of two-flavor SYK models, in which two SYK quantum dots are coupled by q-body interactions, as studied in Ref. 35. For q ≥ 4, H pSY K was shown to have a low temperature phase which exhibits phase separation: one subsystem (say the one on the left) has density 1/(q + 1), and the other one has density q/(q + 1). This phase has a gap to single-particle excitations, and spontaneously breaks particle-hole symmetry and left-right interchange symmetry, but conserves their product.
Interestingly, |ψ(a min ) reproduces this density imbalance perfectly: we find ρ L,min = 1 − ρ R,min = 1/(q + 1). Further, we find E min to be equal to the ground state energy of H pSY K (which can be obtained in a similar fashion at E SY K , see SM. S3)) within numerical accuracy for q ≥ 4. We checked that this agreement even extends to the asymmetric case of p = N R /N L = 1. In the context of two-flavor SYK models, this ratio gives the relative size of the two dots [35]. The comparison for ρ L , ρ R , E min with the exact values is shown in fig. 3(b), for q = 4. The only discrepancy appears as p → 0, which is expected since H pSY K undergoes an additional phasetransition to a gapless phase at p c 0.072 [35]. Another discrepancy appears for q = 2, in which case the variational wavefunction fails to describe the Fermi liquid phase of H pSY K which survives down to T = 0 (see SM. S4 for more details).
Entanglement.-The entanglement properties of |ψ(a) can also be calculated using a recently developed formalism [36]. Some of the earlier studies on entanglement in the SYK model can be found in Ref. [61][62][63][64]. We focus on the second Rényi entropy, , for a biparition of the system into regions A and B, and whereρ A = Tr B |ψ(a) ψ(a)| is the reduced density matrix. The partition is parametrized by x ∈ [0, 1], which gives the proportion of orbitals in A. For x ≤ 0.5 we take region A to be entirely comprised of the left-side fermions, while x > 0.5 also includes a portion (x − 0.5) of the right-side fermions. The large-N limit for S (2) is obtained using an approach similar to calculating log(N ) (see SM. S5).
The results for S (2) are shown in Fig. 4 for q = 4. The x dependence of S (2) resembles the one obtained for KM states in SYK [50], with a small-x linear behavior indicative of a volume law of entanglement, and a maximum at x = 0.5. Starting from 0 at a = 0, the entanglement grows until the left-right symmetric point a = a s is reached, after which it decays monotonically (see Fig. 4 inset). Remarkably, we find S (2) (x) = min(x, 1−x) log(2) at a = a s , which means |ψ(a s ) is maximally entangled between the left and right subsystems.
Discussion.-In this work, we have highlighted a fundamental difference between bosonic (or qubit) and fermionic q-local Hamiltonians, as regards to the complexity of finding wavefunctions with a finite approximation ratio (E ψ /E GS > 0) in the thermodynamic limit. We showed that, for a prototypical fermionic model, the SYK model, the bandwidth of Gaussian states scales subextensively with system size, leaving a parametrically large gap between the ground state and Gaussian states. This raises the question of whether other classes of tractable wavefunctions could (partially) fill this gap. We took a step in that direction by proposing a wavefunction inspired by the variational coupled cluster algorithm with a disorder-averaged approximation ratio of r ≈ 0.62.
From a physical perspective, this wavefunction is easily tractable, since it is described by a static large-N field theory for which saddle point equations are simply given by polynomial equations. It remains however unpractical from a computational point of view since a "bruteforce" calculation of its properties would have factorial complexity on a classical computer. Further, to the best of our knowledge, there is no efficient algorithm to prepare a VCC state on a quantum computer. It is therefore desirable to find other classes of wavefunctions with r > 0 which could efficiently be studied with a classical or quantum computer. Unitary coupled cluster states are particularly promising regarding the latter possibility [7].
Moreover, our approach of focusing on a subset of terms in the SYK Hamiltonian could be transposed to other versions of SYK models with a reduced number of terms, like low-rank SYK [65] and sparse SYK [66]. More generally, we surmise that large-N techniques and SYK models could prove a useful tool in the search for new variational wavefunctions.
Acknowledgments In this section we show that the energy bandwidth of Gaussian states is subextensive for the SYK model with q > 2. We start with a derivation specific to q = 4, and then treat the more general case by mapping it to a classical spin glass model.
Case of q = 4.-We use the Majorana version of SYK for convenience, written as where γ i represent the Majorana fermions with {γ i , γ j } = δ ij . Using Wick's theorem and the permutation properties of J ijkl , the expectation value of the Hamiltonian for an arbitrary Gaussian state can be written as (S1.2) Interpreting J i<j,k<l as a real symmetric matrix and L i<j ≡ i γ i γ j as a vector, we go to the eigenbasis of J, leading to where λ µ are the eigenvalues of J i<j,k<l , and L µ are its eigenvectors. Minimizing H now amounts to minimizing this quadratic form, but with an extensive number of constraints on the values of L µ in order for them to be consistent with a Gaussian state. In order to obtain a non-trivial bound on H , it is sufficient to take into account the simplest of such constraints, which sets the norm of the vector L: Minimizing the quadratic form under this single constraint is straightforward, and leads to the following bound: where λ min is the smallest eigenvalue of J.
We now need to find the scaling of λ min . In the large N limit, we expect the matrix J i<j,k<l to behave as a random matrix of dimension O(N 2 ) × O(N 2 ), and thus to have a semi-circle distribution of eigenvalues with radius O(N ) (this was verified numerically for N up to 200). We therefore expect λ min to be a negative number of order N . From Eq. S1.5, this means that the bandwidth of Gaussian states scales at most like √ N (whereas the full bandwidth scales like N since it is extensive).
General case.-We can show that the above sub-extensive scaling also holds when q > 4 by mapping the problem to the p-spin spherical spin glass model [67]. To do this, we start from the following Hamiltonian: (S1.6) We compute the expectation value for a Gaussian state in a similar way as above, leading to: (S1.7) Denoting a 1 = (i 1 < i 2 ) and similarly for the other indices, we rewrite the expectation value as H = 1 N (q−1)/2 q! (q/2)! a1,...,a q/2 J a1...a q/2 L a1 . . . L a q/2 , (S1.8) where L a is again understood as a N (N − 1)/2-dimensional vector. Even though there exists a large number of constraints on the vector L, we find again that it is sufficient to impose the simplest one ( a L 2 a = N/8) to obtain a non-trivial bound. This will provide the spherical constraint for the mapping to the spherical p-spin model.
The p-spin spherical model is defined as [67] H p-spin = 1 M (p−1)/2 We can now make the following identifications: in order to relate the two models. This finally leads to where E p-spin is the ground state energy of an instance of the spherical p-spin model for which the couplings J a1,...,ap are given by the corresponding J (i1<i2),...,(iq−1,iq) of the SYK Hamiltonian. We now make the assumption that these instances of the p-spin spherical model are typical, or in other words that the correlations present in J (i1<i2),...,(iq−1,iq) due to permutation symmetries can be neglected. If that is the case, we can use the fact that the spherical p-spin model is extensive to deduce that E p-spin scales like M ∼ O(N 2 ). By using this relation, the right-hand side of Eq. S1.13 can be shown to scale like N 3 2 − q 4 . The bandwidth of Gaussian states therefore scales at most like N 3 2 − q 4 , which is subextensive for q > 2. Setting q = 2, we find a Gaussian state bandwidth which is extensive, as expected since in that case the ground state is a Gaussian state. For q = 4, we find √ N as previously shown. For larger q, the Gaussian states' bandwidth gets narrower and narrower.

S2: Large-N analysis of the variational wavefunction
In this section, we discuss the details pertaining to the computation of log(N ) (see eqn. 9) in the large-N limit. As stated in the main text, the said quantity works as a generating functional for computing observables and correlation functions for the variational wavefunction. Since calculating the disorder average of the log-term directly is hard, we use the replica trick to represent the term as where R denotes the number of replicas. The normalization N (a) = 0 | exp(−aT ) exp(−aT † )|0 (see eqn. 3 for the definition ofT ) can be written as an integral over the fermionic-coherent states |c i∈L , |h j∈R , representing the particles and holes, such that where the Grassmann-numbersc ir , c ir ,h ir , h ir are indexed by the replica index r and the site-index i. Contrary to usual thermal-field theory, the Grassmann-numbers do not require an imaginary-time τ index since the terms in the cluster-operatorT commute. Disorder averaging eqn. S2.2 over all possible realizations of J i1···i (q/2) ;j1···j (q/2) , gives us (S2.3) To obtain the large-N limit of the above integral, we introduce the static Green's function G c , G h , and demand that they must satisfy at the large-N saddle point. The above constraints can be incorporated into eqn. S2.3 using the static self-energies Σ c , Σ h such that where Σ c,h act as Lagrange multipliers. We integrate out the fermions from the above to get where we have introduced the effective action S[G, Σ], and 1 represents the identity matrix in the replica-space. We evaluate the integral in eqn. S2.6 at the saddle-point for the action S. Furthermore, we shall consider a replicadiagonal ansatz for G c,h , Σ c,h , i.e. G c,h (r 1 , r 2 ) = δ r1,r2 G c,h and Σ c,h (r 1 , r 2 ) = δ r1,r2 Σ c,h . This results in the following simplified form for the effective action where we have used the site-ratio p = N R /N L and the total number of sites N = N R + N L . Minimizing the above replica-diagonal action with respect to G c,h , Σ c,h we get the saddle-point conditions The value for log(N ) at the saddle-point is given by while the energy-density is obtained from The expression for log(N ) given in eqn. 9 of the main-text is then obtained by setting p = 1 and Similarly, the saddle-point conditions in eqn. S2.8 take the form as reported in eqn. 8 of the main-text.
Density of particles and holes The density of particles, say for the left-side fermions, is obtained by calculating expectation value (S2.13) in the large-N limit. Instead of evaluating the above expression directly, we use a chemical-potential-like source term µ, such that (S2.14) The advantage of using a source-term is that we can repeat the same analysis used for computing the energy earlier (eqn. S2.10) in this case as well. At the end of which we get the following saddle point equations that give back the saddle-point conditions of eqn. S2.8 in the µ → 0 limit. The corresponding replica-diagonal action for the log(· · · ) term in eqn. S2.14 is found to be (S2.16) Using the fact log 0 | exp(−aT ) exp(µĉ † iĉ i ) exp(−aT † )|0 = S ρ L (µ)/R, we compute the derivative of S ρ L (µ)/R w.r.t µ as shown below which according to eqn. S2.14 gives us the density of particles on the left side Similarly, the density of holes on the right-side, i.e. ĥ † iĥ i , can be calculated by using the source-term exp(µĥ † iĥ i ) in place of exp(µĉ † i c i ) in eqn. S2.14 to get ĥ † iĥ i = 1 + G h , from which the density of right-side fermions (particles) can be determined to be

S3: Thermal field theory for the partitioned-SYK model
We now discuss the the thermal field theory for the partitioned-SYK model. This will allow us to compute the exact properties for the ground-state when the temperature T is extrapolated to zero. We reiterate the Hamiltonian for the partitioned-SYK model for ease of access (S3.2) using the fermionic coherent states |c i=1···N [35] described by anti-periodic Grassmann fieldsc i (τ ), c i (τ ) living on the imaginary-time interval τ ∈ [0, β]. In the above equation, we have defined the action S[c, c] whose saddle-point would give us access to the large-N limit. Since, we are interested in the disorder averaged free-energy we use the replica trick, log(Z) = lim R→0 (Z R − 1)/R, yet again, to perform the averaging over J i1···i (q/2) ;j1···j (q/2) . The replica-partition function Z R is found to be where S R denotes a new action over replicas and the Grassmann fieldsc i,r1 (τ ), c i,r2 (τ ) have picked up the additional replica indices r 1 , r 2 . Introducing the large-N Green's functions for the left-side, right-side fermions along with their associated self-energies Σ (r1,r2) L (τ 1 , τ 2 ), Σ (r1,r2) R (τ 1 , τ 2 ), and subsequently integrating out the fermionic-fieldsc i , c i , etc., we arrive at the following action . Assuming time-translational invariance and a replica-diagonal ansatz for the saddle-point, i.e. G (r1,r2) L,R (τ 1 , τ 2 ) ∝ δ r1,r2 G L,R (τ 1 − τ 2 ) and same for Σ (r1,r2) L,R (τ 1 , τ 2 ), we obtain a simplified form for the replica-action where we have used the site-ratio p = N R /N L and defined the action-per-replica S. The saddle-point conditions are found to be by minimizing S w.r.t. G, Σ. The above equations were solved iteratively [35], for a given value of T and p, after discretizing the imaginary-time interval [0, β]. The free-energy can then be calculated by plugging the solutions of eqn. S3.8 in eqn. S3.7 and using The ground-state energy density E pSY K is calculated using the thermodynamic relation which follows from the usual definition of the two-point Green's functions. We access the energy-density and particledensity for the ground-state by numerically extrapolating the values for small but finite T to T → 0.
S4: The non-interacting (q = 2) partitioned-SYK model In this section, we study the non-interacting limit of the partitioned-SYK model on a system of 2N sites. The Hamiltonian for the model is obtained by setting q = 2 in eqn. 3, which is 2) The single-particle spectrum of the above model was checked to be gapless via exact-diagonalization of the Hamiltonian. The energy for our variational ansatz was also calculated using the large-N approach discussed in the main-text (see SM. S2) and the minimized energy was found to be E min = −0.3849.... Interestingly, due to the non-interacting nature of the Hamiltonian, we can calculate this value for energy analytically. We now discuss the analytical approach. For simplicity, let us take J ij to be a real N by N symmetric matrix, with Gaussian matrix elements having variance 1. We can then diagonalizeT † , leading toT with µ ∼ O(1) distributed according to semi-circle law. The operators C † µ , C µ and D † µ , D µ represent the singleparticle eigen-states (orbitals) obtained after diagonalization and obey fermionic anti-commutation relations. We express the variational ansatz in the following way where µ represents the direct product operation and |01 µ denotes the state where the C µ -orbital is occupied and the D µ -orbital is empty, while the reverse is true for the state |10 µ . The energy for the above ansatz is then obtained as where we have defined the matrix J = 1 √ N J ij . We can now take the disorder average using random matrix theory where C n are the Catalan numbers. Using which we find FIG. S1. Comparison of the exact ground-state energy-density EpSY K (circles), obtained using thermal-field theory by extrapolating T → 0, with the prediction from the variational ansatz (eqn. 5) Emin (squares). The match is excellent for q ≥ 4 since the partitioned-SYK model breaks PH-symmetry and develops a gap in the single-particle spectrum. However, when q = 2, i.e. the non-interacting limit of the partitioned-SYK model, the single-particle spectrum is gapless and the prediction from the ansatz deviates from the exact value. The ground-state energy ESY K (denoted with crosses) for the full q-body SYK model shown for comparison.
where F (x) is the ordinary generating function of the Catalan numbers, given by This leads to an expression for the disorder-averaged energy The minimum of E/2N occurs at a = √ 3/2, and the minimum value is E min = E/2N = −2/(3 √ 3) −0.38, which is equal to the value reported in the beginning of this section.
Comparison with exact ground state (GS) The exact ground state (GS) at half-filling is given by where the states |01 µ , |10 µ have the same meaning as described below eqn. S4.4. The ground-state energy-density is given by Since µ are distributed according to the semicircle law, taking the disorder average leads to Comparing the energy-density of the wavefunction (E min ≈ −0.38) with E pSY K , we find E min > E pSY K (only slightly). More importantly, we see that, unlike q ≥ 4 case, the wavefunction does not predict the energy exactly when the ground-state is gapless, see fig. S1. In order to estimate entanglement within the variational wavefunction, we divide the system into parts A (subsystem) and B (rest) and compute the reduced density matrixρ A from the full density matrix We denote the fraction of sites in A as x. We demonstrate the computation of Rényi-entropy for x = 0.5, i.e. A is comprised of the left-side fermions, and report the result for arbitrary x at the end. Additionally, we also set the site-ratio p = 1.0. The reduced-density matrixρ A , when x = 0.5, is found to bê where we have used the fermionic-coherent state |h i∈R for the holes and their corresponding Grassmann numbersh i , h i . The symbol |0 c denotes the vacuum for theĉ i∈L fermions. Since, the 2-nd Rényi-entropy is related to the second moment of the reduced-density matrix, i.e.