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

Massimo Ostilli^1,2 and Carlo Presilla^1,2,3

¹Dipartimento di Fisica, Università di Roma `La Sapienza', Piazzale A. Moro 2, Roma 00185, Italy
²Center for Statistical Mechanics and Complexity, Istituto Nazionale per la Fisica della Materia, Unità di Roma 1, Roma 00185, Italy
³Istituto Nazionale di Fisica Nucleare, Sezione di Roma 1, Roma 00185, Italy

Email: carlo.presilla@roma1.infn.it

Received 16 April 2004
Published 12 August 2004

Contents

• 1. Introduction
• 2. Exact probabilistic representation of lattice dynamics
• 3. Canonical decomposition of expectation
• 4. Canonical averages via a central limit theorem
  • 4.1. Hard-core bosons
  • 4.2. Fermions
• 5. Ground-state correlation functions
• 6. Numerical results
• 7. Conclusions
• Acknowledgments
• References

1. 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 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 with 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]-[7] in the limit when the latter becomes exact [8].

In the present 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 and 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 [8]-[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 for 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.

2. Exact probabilistic representation of lattice dynamics

We illustrate our approach for 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_iσ²  =  0. We are interested in evaluating the matrix elements n|e ^- Ht|n₀, 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 neighbours and η_ij  =  0 otherwise. We will also assume γ_i  =  γ and δ_i  =  0. We shall call such a model the first neighbour 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_ijσ^t}, (i, j) Γ, be a family of 2|Γ| independent Poisson processes with jump rate ρ. At each jump of the process N_ijσ^t, if λ_ijσ ≠ 0 a particle moves from site i to site j or vice versa, whereas the lattice configuration n remains unchanged if λ_ijσ  =  0. The total number of jumps at time t is N_t  =  ∑_{(i,j)Γ,σ = ↑↓}N_ijσ^t. By ordering the jumps according to the times s_k, k  =  1,  ..., N_t, which take place in the interval [0,t), we define a trajectory as the Markov chain generated from the initial configuration n₀. Let us call and 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 is defined by

if N_t  >  0 and if N_t  =  0. Here V₀  =  V(n₀) and s₀  =  0.

Several quantities can be obtained from the matrix elements (4). The ground-state energy is given by

3. Canonical decomposition of expectation

To evaluate (6), we decompose the expectation as a series of conditional expectations with a fixed number of jumps (canonical averages)

where Ω_N  =  Ω_N(n₀) is the set of trajectories with N jumps branching from the initial configuration n₀ and

The weights (9) are obtained on multiplying (5) by the infinitesimal probability e ^- 2|Γ|ρtρ^Nds₁ds₂ ...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 n₀, n₁,  ..., n_N - 1. If A_k  =  ∑_{(i,j)Γ,σ = ↑↓}|λ_ijσ(n_k)| is the number of active links in the configuration n_k, the probability associated with the trajectory r is p_N^(r)  =  ∏_k = 0^N - 11/A_k^(r) and we have

Note that 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. ∑_AN_A  =  N. For the FNU model, the possible values of A depend on the number of particles and we have the bound A≤ min(2d(N_↑  +  N_↓), 2|Γ|).

For a generic trajectory, weights (9) satisfy the recursive differential equation

where . In terms of the Laplace transform we get

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. ∑_VN_V  =  N  +  1. For model (1), the possible values assumed by V are V  =  0, γ, 2γ, ..., N_pγ, where N_p  =  min(N_↑, N_↓).

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 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₀ is the solution of the equation

Note that in the case γ  =  0, equation (13) reduces to Stirling's approximation of the exact value .

4. Canonical averages via a central limit theorem

As evident from the explicit expression given in the case γ  =  0, the weights 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 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 behaviour of the stochastic variables N_V and N_A. In the following, we will indicate by m_V and m_A the number of different values assumed by the variables V and A, respectively. In this limit, we will not distinguish the different normalizations, N  +  1 and N, of N_V and N_A, respectively. For clarity, we consider separately the hard-core boson and fermion cases.

4.1. Hard-core bosons

In this case, we have 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

equation (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₀  +  V)/ε] ...; ...log A  ...), respectively. Note that u depends on ν through x₀  =  x₀(ν).

The probability density is given by the fraction of trajectories branching from the initial configuration n₀ and having after N jumps multiplicities N_V  =  ν_VN and N_A  =  ν_AN. For N large, it can be approximated in the following way. We rewrite the multiplicities as N_V  =  ∑_k = 0^Nχ_V(n_k) and N_A  =  ∑_k = 1^N χ_A(n_k - 1), where χ_V(n)  =  1 if V(n)  =  V and χ_V(n)  =  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 or  [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 - 1V-like components and m_A - 1A-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 is the normal density defined in terms of the vector having the m  -  2 chosen components of ν

Here and 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 , whereas rows and columns of the VV, VA, 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

To evaluate the expectation value , we need to sum the series (7). This can be done with a further saddle-point integration. According to equation (20), the terms of this series are exponentially peaked at N  =  N^sp, where N^sp satisfies

with v^T  =  ( ...(x₀  +  V) ^- 1 ...; ...0 ...). We observe that, according to equations (14) and (21), the term vanishes for ν  =  ν^sp so that the above condition reduces to

Equation (23) is a time-independent equation, which determines as a function of and Σ. According to equation (14), this means that for ν  =  ν^sp, the quantity N^sp increases linearly with time so that 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 equation (27) is a nontrivial formula for the ground state of a system of bosons interacting via a hard-core potential. With the above expression equation (26) can be written more compactly as

By using the bounds E_0B⁽⁰⁾  <  E_0B  <  0, the scalar equation (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 to take into account all the possible values of the generalized potential V corresponding to the operator ∑_iΛγ_ic_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 equation (10) consists in replacing ∏_k = 0^N - 1A_k with ∏_k = 0^N - 1(A_kη_k/ε). Explicitly, now the vectors and u are

and the generalized ground-state energy E_0B⁽⁰⁾ corresponding to γ  =  δ  =  0 reads

Similar considerations hold for Hamiltonians with arbitrary potential operators.

4.2. Fermions

In this case, we have 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 by . The approach developed for hard-core bosons can be extended also to fermions by including this further multiplicity N _- . We will report on this procedure elsewhere. In the present paper, 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. 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 ∑_nn| e ^- Ht|n₀ in powers of t and comparing with the expansion (7), for γ  =  0 and N large, we obtain for hard-core bosons and fermions, respectively,

where c_0B(n₀) and c_0F(n₀) are coefficients related to the initial configuration n₀ and E_0B⁽⁰⁾ and E_0F⁽⁰⁾ are the γ  =  0 ground-state energies. The average weighted sign after N jumps for γ  =  0 is then given by

For fermions the noninteracting energy, E_0F⁽⁰⁾, is known exactly, whereas for hard-core bosons E_0B⁽⁰⁾ can be computed with Monte Carlo simulations or analytically as shown above. In general, E_0F⁽⁰⁾/E_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 . Therefore, this can be evaluated as for hard-core 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 in the fermion case. Finally, as in the hard-core boson case, one finds that

and equation (35) becomes equal to equation (28) with E_0B and E_0B⁽⁰⁾ substituted by E_0F and E_0F⁽⁰⁾, respectively. However, in this case we do not have the analogue of equation (27) and E_0F⁽⁰⁾ must be provided by an independent calculation, i.e. by exact diagonalization of the separable many-particle Hilbert space.

5. 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₀(ξ) is the ground-state energy of H(ξ). Hereafter, we assume a normalized ground state, E₀(ξ)|E₀(ξ)  =  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. O  =  c_iσ^†c_jσ, with i ≠ j, O  =  c_i↑^†c_i↑c_i↓^†c_i↓ or O  =  c_iσ^†c_iσ if H is the generalized Hamiltonian (1). In this case, by using equation (37) for hard-core bosons, we have

The expressions in the second lines of (38)-(40) have been obtained by using the derivatives of equation (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 for fermions by using the derivatives of equation (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₀(ξ). 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 ith lattice point. The quantum expectation of the operators S_iS_j in the ground state of the hard-core boson FNU model, E_0B(ε, γ)|S_iS_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.

6. Numerical results

We now apply the approach developed in the previous sections to some example cases. In particular, we compare the ground-state energy obtained by equations (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 behaviour 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 with the measurement of and Σ are negligible on the scales considered.

In figure 1, we show the behaviour of E_0B for the hard-core boson FNU model as a function of the interaction strength γ for several lattice systems. The solution of equation (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 and , whereas the function g_N given by equation (17) depends on ν_VN and ν_AN. This implies that the tails of the probability density give a finite contribution to the integral (20). Furthermore, from the structure of g_N, 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 1. Ground-state energy per particle for the hard-core boson FNU model versus the interaction strength γ. The results from equation (28) are compared with those from exact diagonalizations ( + ) and QMC simulations (×) for two-dimensional systems with periodic boundary conditions and different Lx×Ly sites and N↑  =  N↓  =  Np particles. The statistical errors in the QMC data are negligible in this scale.

In figure 2, we show the behaviour of E_0F evaluated according to equation (35) as a function of the interaction strength γ for the same cases considered in figure 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 for hard-core bosons.

We stress that, once and Σ are known, our approach provides any other ground-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 figure 3, we report the spin-spin structure factor S(q_x, q_y) evaluated by using equations (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 . Note that in the 4×4 system with periodic boundary conditions considered in figure 3, we have only five different values for E₀|S_iS_j|E₀ corresponding to the five possible distances d_ij  =  |x_i  -  x_j|  +  |y_i  -  y_j|  =  0, 1, 2, 3, 4. In figure 3, we also show the behaviour of E₀|S_iS_j|E₀ 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.

7. 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 and . 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 ground-state correlation functions are particularly simple; and (ii) the corresponding results compare rather well with the exact ones in regions of physical interest.

In the present 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.

Acknowledgments

We thank F Cesi for stimulating discussions and G Jona-Lasinio for insightful comments and a critical reading of the paper. This work was supported in part by Cofinanziamento MIUR protocollo 2002027798_001.

References

[1]

Simon B 1979 Functional Integration and Quantum Physics  (New York:  Academic) 

[2]

De Angelis G F, Jona-Lasinio G and Sirugue M 1983 J. Phys. A: Math. Gen. 16 2433 
IOPscience

[3]

De Angelis G F, Jona-Lasinio G and Sidoravicius V 1998 J. Phys. A: Math. Gen. 31 289 
IOPscience

[4]

Beccaria M, Presilla C, De Angelis G F and Jona-Lasinio G 1999 Europhys. Lett. 48 243 
IOPscience

[5]

Trivedi N and Ceperley D M 1990 Phys. Rev. B 41 4552 
CrossRef

[6]

von der Linden W 1992 Phys. Rep. 220 53 
CrossRef

[7]

Calandra Buonaura M and Sorella S 1998 Phys. Rev. B 57 11446 
CrossRef

[8]

Sorella S and Capriotti L 2000 Phys. Rev. B 61 2599 
CrossRef

[9]

Zhang S, Carlson J and Gubernatis J E 1997 Phys. Rev. B 55 7464 
CrossRef

[10]

Carlson J, Gubernatis J E, Ortiz G and Zhang S 1999 Phys. Rev. B 59 12788 
CrossRef

[11]

van Bemmel H J M et al 1994 Phys. Rev. Lett. 72 2442 
CrossRefPubMed

[12]

Chandrasekaran S and Wiese U J 1999 Phys. Lett. 83 3116 
CrossRef

[13]

Beccaria M 2001 Phys. Rev. B 64 73107 
CrossRef

[14]

Billingsley P 1986 Probability and Measure  (New York: Wiley) 

[15]

Fano G, Ortolani F and Parola A 1990 Phys. Rev. B 42 6877 
CrossRef

[16]

Parola A, Sorella S, Baroni S, Car R, Parrinello M and Tosatti E 1989 Physica C 162-164 771 
CrossRef

[17]

Husslein T, Fettes W and Morgenstern I 1997 Int. J. Mod. Phys. C 8 397 
CrossRef