Ground state of many-body lattice systems: an analytical probabilistic approach

On the grounds of a Feynman–Kac-type formula for Hamiltonian lattice systems, we derive analytical expressions for the matrix elements of the evolution operator. These expressions are valid at long times when a central limit theorem applies. As a remarkable result, we find that the ground-state energy as well as all the correlation functions in the ground state are determined semi-analytically by solving a simple scalar equation. Furthermore, explicit solutions of this equation are obtained in the noninteracting case.


Introduction
The Feynman-Kac formula [1] provides a powerful connection between the imaginary time evolution of quantum systems in the continuum and probabilistic expectations over Wiener classical trajectories.Its extension to the case of multicomponent wavefunctions requires the introduction of a further expectation over stochastic Poisson processes [2].The role of Poisson processes is central in the probabilistic representation of Berezin integrals over anticommuting variables and, in general, of the time evolution of discrete systems [3].In the simplified formulation [4] it is shown that the real time or the imaginary time dynamics of systems described by a finite Hamiltonian matrix, representing bosonic or fermionic degrees of freedom, is expressed in terms of the evolution of a proper collection of independent Poisson processes.For a lattice system the Poisson processes are associated to the links of the lattice and their jump rates can be arbitrary.In [4] it is demonstrated that when the rates of the Poisson processes are chosen equal to the hopping coefficients of the system the probabilistic representation leads to an optimal algorithm which coincides with the Green Function Quantum Monte Carlo (QMC) method [5,6,7] in the limit when the latter becomes exact [8].
In this paper we exploit the probabilistic representation [4] to derive analytical expressions for the matrix elements of the evolution operator in the long time limit.Our approach is based on a series expansion of the probabilistic expectation in terms of conditional expectations with a fixed number N of jumps of the Poisson processes.This resembles the expansion of the grand canonical partition function in statistical mechanics in terms of canonical averages with a fixed number of particles.By integrating out the N stochastic jump times, we show that the conditional expectations become averages of functions which depend only on the multiplicities N V and N A of the values assumed by the potential and hopping energies, respectively, of the configurations visited by the system.According to a central limit theorem, at large values of N, i.e. long times, the rescaled multiplicities N V / √ N and N A / √ N become Gaussian distributed and the corresponding averages can be evaluated analytically.The parameters of the Gaussian probability density do not depend on the Hamiltonian parameters and are easily determined statistically.Finally, the series of the conditional expectations can be re-summed with a saddle-point method.As a remarkable result, we find that the ground state energy is semi-analytically determined by solving a simple scalar equation.Once this equation is solved the quantum expectation in the ground state for other operators can be determined analytically by using the Hellman-Feynman theorem.The result is valid for boson or fermion systems.As regards fermions the well known sign problem [9,10,8,11,12,13] is avoided by introducing an approximation related to the exact counting of the positive and negative contributions in the noninteracting case.
The paper is organized as follows.In Section 2 we review the probabilistic representation of quantum dynamics in the case of a generalized Hubbard Hamiltonian.In Section 3 we decompose the expectation in canonical averages of weights which are calculated analytically.The canonical averages are evaluated at long times via a central limit theorem in Section 4. The equation for the ground-state energy is discussed in Section 4.1 for hard-core bosons and in Section 4.2 for fermions.In Section 5 we show how to calculate ground-state correlation functions within our approach.Finally in Section 6 we show some example cases and compare with the results of exact numerical calculations.General features of our approach are summarized and discussed in Section 7.

Exact probabilistic representation of lattice dynamics
We illustrate our approach in the case of imaginary time dynamics for a system of hard-core bosons or fermions described by a generalized Hubbard Hamiltonian where Λ ⊂ Z d is a finite d-dimensional lattice with |Λ| sites and c iσ the commuting or anticommuting destruction operators at site i and spin index σ with the property c 2 iσ = 0. We are interested in evaluating the matrix elements n|e −Ht |n 0 where n = (n 1↑ , n 1↓ , . . ., n |Λ|↑ , n |Λ|↓ ) are the lattice occupation numbers taking the values 0 or 1.The total number of particles is N σ = i∈Λ n iσ for σ =↑↓.In the following we shall use the mod 2 addition n ⊕ n ′ = (n + n ′ ) mod 2. Let Γ be the set of system links, i.e. the unordered pairs (i, j) with i, j ∈ Λ such that η ij = 0.For simplicity, we will start by assuming η ij = ǫ if i and j are first neighbors and η ij = 0 otherwise.We will also assume γ i = γ and δ i = 0. We shall call such a model the first neighbor uniform (FNU) model.For a d-dimensional lattice the number of links per spin component is |Γ| = d|Λ|.Let us introduce where 1 iσ = (0, . . ., 0, 1 iσ , 0, . . ., 0), and let {N t ijσ }, (i, j) ∈ Γ, be a family of 2|Γ| independent Poisson processes with jump rate ρ.At each jump of the process N t ijσ , if λ ijσ = 0 a particle moves from site i to site j or vice versa, while the lattice configuration n remains unchanged if λ ijσ = 0.The total number of jumps at time t is N t = (i,j)∈Γ,σ=↑↓ N t ijσ .By ordering the jumps according to the times s k , k = 1, . . ., N t , at which they take place in the interval [0, t), we define a trajectory as the Markov chain n 1 , n 2 , . . ., n Nt generated from the initial configuration n 0 .Let us call λ 1 , λ 2 , . . ., λ Nt and V 1 , V 2 . . ., V Nt the values of the matrix elements (2) and (3) occurring along the trajectory.As proved in [4], the following representation holds where the stochastic functional M t is defined by if N t > 0 and M t = e 2|Γ|ρt e −V 0 t if N t = 0.Here V 0 = V (n 0 ) and s 0 = 0. Several quantities can be obtained from the matrix elements (4).The ground-state energy is given by

Canonical decomposition of expectation
To evaluate (6) we decompose the expectation E(M t ) as a series of conditional expectations with a fixed number of jumps (canonical averages) where Ω N = Ω N (n 0 ) is the set of trajectories with N jumps branching from the initial configuration n 0 and The weights (9) are obtained on multiplying (5) by the infinitesimal probability e −2|Γ|ρt ρ N ds 1 ds 2 . . .ds N to have N jumps and integrating over the jump times.
A link ij with spin σ is called active if λ ijσ = 0. From (8) it is clear that only trajectories formed by a sequence of active links contribute to (7).Hereafter we restrict Ω N to be the set of these effective trajectories with N jumps.The sum over the set Ω N in (7) can be rewritten as an average, • , over the trajectories with N jumps generated by extracting with uniform probability one of the active links available at the configurations is the number of active links in the configuration n k , the probability associated to the trajectory r is p (r) and we have Note that N −1 k=0 A k = A A N A depends only on the multiplicities N A of the values A assumed by the number of active links; these multiplicities are normalized to N, i.e.
A N A = N.For the FNU model the possible values of A depend on the number of particles and we have the bound For a generic trajectory the weights (9) satisfy the recursive differential equation where W −1 (t) = 0.In terms of the Laplace transform In (12) we see that the weights depend only on the multiplicities N V of the values V assumed by the potential; these multiplicities are normalized to N + 1, i.e.
For an assigned set of multiplicities of the potential, the antitransform of ( 12) can be evaluated with the residue method.In this way one obtains an exact recursive expression of W N (t) which, however, for large N must be evaluated numerically with multi-precision algebra.On the other hand, for N large a complex saddle-point method can be used which provides the following asymptotically exact explicit expression where x 0 is the solution of the equation Note that in the case γ = 0, Eq. ( 13) reduces to Stirling's approximation of the exact value W N (t) = ǫ N t N /N!.

Canonical averages via a central limit theorem
As is evident from the explicit expression given in the case γ = 0, the weights W N (t) have a maximum at some N which increases by increasing t.This remains true also in the general case, as shown in the following.Therefore in the long time limit the most important contributions to the expansion (7) of the expectation E(M t ) come from larger and larger values of N. In this section we will evaluate the canonical averages (10) analytically for N large by using the asymptotic behavior of the stochastic variables N V and N A .In this limit we will not distinguish the different normalizations, N + 1 and N, of N V and N A , respectively.In the following, we will indicate with m V and m A the number of different values assumed by the variables V and A, respectively.For clarity, we consider separately the hard-core boson and fermion cases.

Hard-core bosons
In this case we have S N = 1 and the canonical averages ( 10) are averages of a function which depends only on the multiplicities, N V and N A (besides a parametric dependence on time).In terms of the corresponding frequencies, ν V = N V /N and ν A = N A /N, which for N large become continuously distributed in the range [0, 1] with the constraints Eq. ( 10) can be rewritten as where ν and u are vectors with m = m V + m A components defined as ν T = (. . .ν V . . .; . . .ν A . ..) and u T = (. ..−log[(x 0 +V )/ǫ] . . .; . . .log A . ..), respectively.Note that u depends on ν through x 0 = x 0 (ν).
The probability density P N (ν) is given by the fraction of trajectories branching from the initial configuration n 0 and having after N jumps multiplicities For N large, it can be approximated in the following way.We rewrite the multiplicities as 0 otherwise, and similarly for χ A .Since the configurations n k form a Markov chain with finite state space, a central limit theorem applies to each rescaled sum N V / √ N or N A / √ N [14].However, due to the constraints (15), the joint probability for these m rescaled sums is not Gaussian.Given an arbitrary set of m V − 1 V -like components and m A − 1 A-like components, the joint probability density is the product of a Gaussian density for this set of variables and two delta functions which take into account the constraints (15).For the frequencies ν, therefore, we have where F N ( ν) is the normal density defined in terms of the vector ν having the m − 2 chosen components of ν Here ν and ΣN −1 are the (m − 2)-component subvector and submatrix, respectively, of the mean value, ν, and the covariance matrix, ΣN −1 , of ν.As discussed in Section 6, the quantities ν and Σ are easily measured by sampling over trajectories with a large number of jumps.By using (18), the m-dimensional integral which appears in ( 16) can be performed by the saddle-point method.Note that this integration method is asymptotically exact for N large.Due to the constraints (15), ν satisfies the property V ν V = A ν A = 1 while rows and columns of the V V , V A, AV , and AA blocks of Σ are normalized to zero.By using these properties, in terms of ν and Σ we get where ν sp is the saddle-point frequency defined by the equation In order to evaluate the expectation value E(M t ) we need to sum the series (7).This can be done with a further saddle-point integration.According to Eq. ( 20) the terms of this series are exponentially peaked at N = N sp , where with v T = (. . .(x 0 + V ) −1 . . .; . . .0 . ..).We observe that, according to Eqs. ( 14) and ( 21), the term {t − N [(ν, v) + (Σu, v)]} vanishes for ν = ν sp so that the above condition reduces to Equation ( 23) is a time independent equation which determines x 0 | ν=ν sp ,N =N sp as a function of ν and Σ.According to Eq. ( 14), this means that for ν = ν sp , the quantity N sp increases linearly with time so that x 0 | ν=ν sp ,N =N sp becomes independent of time.
In conclusion, a saddle-point integration with respect to N of the series (7) provides By taking the time derivative of this expectation and using (6), we obtain that the ground-state energy of the hard-core boson system is Equation ( 23) is, therefore, the equation for the ground-state energy.It defines E 0B in terms of ν and Σ and explicitly reads In the case γ = 0 the ground-state energy E 0B can be solved explicitly and one has Note that Eq. ( 27) is a nontrivial formula for the ground state of a system of bosons interacting via a hard-core potential.With the above expression Eq. ( 26) can be written more compactly as By using the bounds E (0) 0B < E 0B < 0, the scalar Eq. ( 28) can be easily solved with the bisection method.
The generalization of the above results to Hubbard Hamiltonians (1) with arbitrary parameters η, γ, δ is straightforward.Equations (20-25) remain formally unchanged, however the vectors ν, u, v and the covariance matrix Σ are modified in order to take into account all the possible values of the generalized potential V corresponding to the operator i∈Λ γ i c † i↑ c i↑ c † i↓ c i↓ + i∈Λ σ=↑↓ δ iσ c † iσ c iσ and all the possible values of the generalized kinetic quantities T = Aη/ǫ, where now ǫ is a unity of energy and η/ǫ is the dimensionless hopping value corresponding to the current jumping link.In fact the generalization of Eq. ( 10) consists in replacing N −1 k=0 A k with N −1 k=0 (A k η k /ǫ).Explicitly, now the vectors ν and u are ν T = (. . .ν V . . .; . . .ν T . ..) and the generalized ground-state energy Similar considerations hold in the case of Hamiltonians with arbitrary potential operators.

Fermions
In this case we have S N = ±1 and the canonical averages (10) are averages of a function which depends not only on the multiplicities, N V and N A but also on N − , the sign multiplicity related to S N by S N = (−1) N − .The approach developed in the case of hard-core bosons can be extended also to fermions by including this further multiplicity N − .We will report on this procedure elsewhere.Here, we introduce an approximation which allows to reduce the calculation of the fermion ground-state energy to that of an effectively modified hard-core boson system.This approximation is motivated by the observation that the correlations between N − and N V are smaller than those between N − and N A .
Let us consider again the FNU model.In order to evaluate (10) for a fermion system we introduce the average weighted sign s N after N jumps The quantity s N is a function of the interaction strength.For γ = 0 it can be evaluated in the following way.Expanding n n|e −Ht |n 0 in powers of t and comparing with the expansion (7), for γ = 0 and N large we get in the case of hard-core bosons and fermions, respectively where c 0B (n 0 ) and c 0F (n 0 ) are coefficients related to the initial configuration n 0 and E (0) 0B and E (0) 0F are the γ = 0 ground-state energies.The average weighted sign after N jumps for γ = 0 is then given by In the case of fermions the noninteracting energy, E 0F , is known exactly while for hardcore bosons E (0) 0B can be computed with Monte Carlo simulations or analytically as shown above.In general 0B < 1 so that s N vanishes exponentially for N large.Approximating s N with its value (34) at γ = 0 removes effectively the negative signs in the expectation E(M t ).Therefore this can be evaluated as in the case of hardcore bosons with the same Gaussian probability density.In particular the saddle-point condition for N = N sp now becomes which is a time independent equation determining x 0 | ν=ν sp ,N =N sp in the fermion case.Finally, as in the hard-core boson case one finds that and Eq. ( 35) becomes equal to Eq. ( 28) with E 0B and E (0) 0B substituted by E 0F and E (0) 0F , respectively.However, in this case we do not have the analogue of Eq. ( 27) and E (0) 0F must be provided by an independent calculation, i.e. by exact diagonalization of the separable many-particle Hilbert space.

Ground-state correlation functions
The equation obtained in the previous Section for the determination of the ground-state energy depends on the parameters of the Hamiltonian only explicitly through the values of the generalized potentials V and of the kinetic quantities T .In fact, the statistical moments ν and Σ are determined by the structure of the Hamiltonian not by the values of the Hamiltonian parameters.Therefore we are able to evaluate the derivatives of the ground-state energy with respect to any parameter ξ of the Hamiltonian H(ξ).This allows the determination of arbitrary ground-state correlation functions via the Hellman-Feynman theorem where E 0 (ξ) is the ground-state energy of H(ξ).Hereafter, we assume a normalized ground state, E 0 (ξ)|E 0 (ξ) = 1.Suppose that we want to evaluate the quantum expectation of an operator O in the ground state of the Hamiltonian H.We have two possibilities.
i) The operator O is a term of the Hamiltonian itself, e.g.
. In this case, by using Eq.(37) for hard-core bosons we have The expressions in the second lines of (38-40) have been obtained by using the derivatives of Eq. ( 23), which for a generic parameter ξ read where ν sp is given by ( 21) and is determined once E 0B (η, γ, δ) has been solved.Similar expressions hold in the case of fermions by using the derivatives of Eq. ( 35) ii) If the operator O is not a term of the Hamiltonian H we consider a new Hamiltonian H(ξ) = H +ξO and calculate the corresponding ground-state energy E 0 (ξ).Note that, since the used probabilistic representation holds for any system described by a finite Hamiltonian matrix [4], the nature of the operator O is arbitrary.As an example we study the spin-spin structure factor where S i = c † i↑ c i↑ − c † i↓ c i↓ and x i and y i are the coordinates of the i-th lattice point.The quantum expectation of the operators S i S j in the ground state of the hard-core boson FNU model, E 0B (ǫ, γ)|S i S j |E 0B (ǫ, γ) , can be obtained by considering the Hamiltonians where H represents the FNU model.For these Hamiltonians, the possible values of the potential are V = mγ + kξ ij , with m = 0, 1, . . ., N p and k = −1, 0, 1, and we have where E 0B (ǫ, γ, ξ ij ) is the ground-state energy of the Hamiltonian (44) and is calculated as explained in Section 4.1.

Numerical results
Now we apply the approach developed in the previous Sections to some example cases.
In particular, we compare the ground-state energy obtained by Eqs. ( 28) and (35) with that from exact numerical calculations.
In our approach the starting point is the evaluation of the statistical moments ν and Σ.These are obtained by generating trajectories in the lattice configuration space and counting the multiplicities N V and N A .The length of the trajectories is chosen to be sufficiently large for the asymptotic behavior to be established.The determination of ν and Σ with good statistical precision requires a number of trajectories which increases no more than linearly with the number of lattice sites |Λ|.Therefore the evaluation of these moments is feasible even for large systems.In the following applications the statistical errors associated to the measurement of ν and Σ are negligible on the scales considered.
In Fig. 1 we show the behavior of E 0B for the hard-core boson FNU model as a function of the interaction strength γ for several lattice systems.The solution of Eq. (28) compares rather well with the results of exact diagonalizations or QMC simulations.There is a small systematic error which grows with increasing γ and/or the system size and which is maximum at half density.This error is related to the Gaussian shape (19) assumed for the asymptotic probability density.In fact, the mentioned central limit theorem for Markov chains applies to the variables ν V √ N and ν A √ N while the function g N given by Eq. ( 17) depends on ν V N and ν A N. This implies that the tails of the probability density give a finite contribution to the integral (20).From the structure of g N , furthermore, it is evident that this error becomes large when the components of the vector u assume large values, i.e. for large values of γ and/or large system sizes.Figure 2. As in Fig. 1 for the fermion FNU model.In this case no QMC simulations are available.Exact diagonalization data for the 4 × 4 systems are taken from [15] (N p = 8) and [16,17] (N p = 5).
In Fig. 2 we show the behavior of E 0F evaluated according to Eq. (35) as a function of the interaction strength γ for the same cases considered in Fig. 1.Compared with the hard-core boson case we observe a further systematic error due to the approximation s N (γ) ≃ s N (0).Depending on the particle density of the system, this error adds to or subtracts from the systematic error due to the Gaussian tails of the probability density discussed in the case of hard-core bosons.
We stress that, once ν and Σ are known, our approach provides any other ground- frag replacements Figure 3. Spin-spin structure factor S(q x , q y ) for the hard-core boson FNU model at γ = 0 in a lattice with 4 × 4 sites, N ↑ = N ↓ = 8 particles, and periodic boundary conditions.For the same system, in the right plot we show the values of the groundstate expectation of the operators S i S j for d ij = 0, 1, 2 as a function of the interaction strength γ (solid lines) compared with Monte Carlo results (dots with error bars).
state quantity analytically as a function of the Hamiltonian parameters.On the other hand, QMC methods require, due to the unavoidable branching or reconfiguration techniques [7], different simulations for different values of the parameters.
As an example of calculation of correlation functions in the ground state, in Fig. 3 we report the spin-spin structure factor S(q x , q y ) evaluated by using Eqs.(43-45) for the hard-core boson FNU model.The value S(π, π) represents a maximum for S(q x , q y ) and for large values of |Λ| is related to the staggered magnetization m s through m s = S(π, π)/|Λ|.Note that in the 4 × 4 system with periodic boundary conditions considered in Fig. 3 we have only five different values for E 0 |S i S j |E 0 corresponding to the five possible distances d ij = |x i − x j | + |y i − y j | = 0, 1, 2, 3, 4. In Fig. 3 we also show the behavior of E 0 |S i S j |E 0 as a function of γ for d ij = 0, 1, 2. These terms provide the most important contributions to S(π, π).For small values of γ, the results compare rather well with Monte Carlo data.

Conclusions
By using saddle-point techniques and a central limit theorem, we have exploited an exact probabilistic representation of the quantum dynamics in a lattice to derive analytical approximations for the matrix elements of the evolution operator in the limit of long times.For both hard-core boson and fermion systems, this development yields to a simple scalar equation for the ground-state energy.This equation depends on the values of the generalized potentials V and of the kinetic quantities T , and on the statistical moments ν and Σ of their asymptotic multiplicities N V and N T .In turn, these moments depend only on the structure of the system Hamiltonian, not on the values of the Hamiltonian parameters.This implies that the statistical moments must be measured una tantum for a given Hamiltonian structure and, once ν and Σ are known, our approach provides the ground-state energy analytically as a function of the Hamiltonian parameters.
In the long time limit, the saddle-point integrations used in our approach are asymptotically exact and the central limit theorem evoked applies rigorously to the rescaled multiplicities N V / √ N and N T / √ N .However, since functions depending on N V and N T are involved, we have a small systematic error related to the finite contributions from the tails of the probability density.This systematic error could be reduced by a large deviation analysis.In fact, equations ( 27) and (28) suggest that the Gaussian approximation for the probability density corresponds to a second order truncation of a cumulant expansion.Anyway, the present Gaussian approach has the following relevant features: i) the equations derived for the ground-state energy and the groundstate correlation functions are particularly simple; ii) the corresponding results compare rather well to the exact ones in regions of physical interest.
In this paper we have considered two-dimensional lattice models at imaginary times.Our approach, however, is valid in any dimension.Similar analytical expressions can be obtained also for the real time evolution.