Perturbation calculation of the uniform electron gas with a transcorrelated Hamiltonian

Electronic structure calculations usually suffer from numerical problems originating from the Coulomb potential $r12−1$⁠. On the one hand, the short-range singularity leads to non-smoothness in the many body wave function, characterized by a cusp at the electron coalescence point r₁₂ = 0,¹ which causes a slow convergence with respect to the basis-set in configuration expansions of the wave function. On the other hand, the slow decay of the long-range tail of the Coulomb potential also causes a major problem in calculations of metallic systems. Any straightforward perturbation treatment leads to divergent values in the thermodynamic limit, and in order to get a meaningful perturbation result, one has to sum over some sub-sequences of the perturbation expansion up to infinite order, a technique which is usually referred to as the random phase approximation (RPA).^2–6 

One way to deal with the short-range singularity of the wave function is to make use of the Jastrow ansatz,⁷ 

to split the cusp into a pair correlation factor

With a factor u fulfilling the asymptotic condition

the smoothness of the function Φ is one order higher than that of the function Ψ.⁸ This product form is widely used in quantum Monte Carlo methods,^9–11 where the expectation value of the Hamiltonian can be evaluated by random sampling. One can, thus, minimize this expectation value of energy by adjusting the correlation factor τ. Due to the involved high-dimensional integrals, it is difficult to use such kind of variational treatment of the Jastrow wave function in conventional quantum chemistry, where the function Φ needs to be further approximated by configuration expansions. Alternatively, the transcorrelated (TC) method of Boys and Handy^12–19 offers a fairly efficient way to deal with the Jastrow ansatz. With this method, the original Schrödinger equation is transformed into a non-Hermitian eigenvalue problem,

via an exact similarity transformation, which removes the involved exponential factor, and results in an exact effective Hamiltonian containing up to three body terms,

This method was initially studied by Boys and Handy for the single determinant ansatz and has only recently been combined with configuration expansion methods of quantum chemistry.^20–28 These works reveal a high level of efficiency of the transcorrelated Hamiltonian in handling short range correlations, that the resolving of the cusp leads to significant speed up of the basis convergence, and furthermore, more speed up can be achieved when more generalized factors are used (instead of the simple r₁₂-type or F(r₁₂)-type factors).

Regarding the study of long-range correlation with the transcorrelated method, there exist only a few works dealing with the single determinant ansatz for uniform electron gases (UEGs).^29–31 By choosing a proper correlation factor τ, the leading order singularity of the Coulomb potential can be removed from the effective Hamiltonian for both the short-range and the long-range singularities. We, therefore, expect that the divergence of the perturbation energy for metallic systems can be cured to some extent. It should at least be possible to get non-divergent results in the thermodynamic limit at low order perturbation theory without any RPA-type treatment. It would, then, be interesting to know whether this kind of non-RPA treatment leads to meaningful results. In this work, we will test the performance of the many body perturbation (MBPT) treatment of the transcorrelated Hamiltonian on UEG in the high density regime, which has a well-established RPA result,^5,32

In Sec. II, we describe the perturbation treatment of the TC Hamiltonian, and in Secs. III and IV, we present the details of the calculations of the first and the second order energies, where a Jastrow factor optimized for the Hartree–Fock reference function is used. In Sec. V, we also present some calculation results with the widely used Coulomb–Yukawa factor, which will be followed by some concluding remarks in Sec. VI.

Uniform electron gases are described by the Hamiltonian

The system is characterized by the charge density ρ (or, equivalently, the Wigner–Seitz radius r_s or the Fermi wavevector k_F). In order to get the explicit dependence of the physical quantities on r_s, we re-scale the spatial coordinates

so that the system keeps a fixed density,

but has a re-scaled Hamiltonian,

The formal volume Ω and the particle number N will finally be taken as ∞ in the thermodynamic limit. In this paper, we will only use this re-scaled form and, thus, will ignore the prime on the variables. We only need to remember that the final energies should be divided by $α2rs2$ according to Eq. (10). Then, we use Jastrow ansatz (1) for the ground state wave function of $Ĥ$⁠, and the two body correlation factor u(r_i, r_j) can be written as u(r_i − r_j) due to the translational symmetry. The transcorrelated Hamiltonian (5) can be derived straightforwardly in the second quantization formalism³¹ in terms of plane wave orbitals,

where $ũ(k)$ is the Fourier transformation of u(r). $F(▿u)2(k)$ is the Fourier transformation of (▿u(r))² and can be calculated via the convolution theorem,

Here, we have used the fact that in the thermodynamic limit, the function $ũ(k)$ is actually a 1D function of k due to the rotational symmetry.

It is of crucial importance to take a proper Jastrow factor u. It has to fulfill asymptotic conditions at both the limits,³¹ 

where $w(k)=4παrsk2$ is the Fourier transformation of the Coulomb potential and ρ is the constant density given in Eq. (9). These conditions remove the leading order singular term of the Coulomb potential and lead to finite second order perturbation energy. However, they do not guarantee the quality of the perturbation energy. Based on some recent numerical studies of the TC method on small molecular systems,^25,26 the correlation factor optimized for a single determinant reference Φ (e.g., the Hartree–Fock wave function) also performs well when Φ is further treated by configuration expansions, at least for weakly correlated systems.

Optimizations of the correlation factor for a single determinant Jastrow ansatz (e^τΦ₀) can be realized by solving either a variational equation or a transcorrelated equation for the factor.³³ Based on a minimization of the variational energy, one can derive a variational equation,

The explicit presence of the exponential factor e^τ in this equation usually leads to major numerical difficulties in solving the equation. For small- or medium-sized systems, it can be directly treated using the variational quantum Monte Carlo method. For large or extended systems, however, the exponential factor can only be treated approximately; for example, the linked cluster expansion^34–36 offers one of such kind approximate treatment. This is essentially a RPA-type treatment of the single determinant Jastrow ansatz. Alternatively, it is possible to avoid the treatment of the exponential factor by using the transcorrelated equation for the correlation factor,^12,14

The exponential factor is removed from the above equation by the similarity transformation, and thus, Eq. (19) can be solved more easily compared with Eq. (18), especially for large or extended systems. Despite this remarkable simplification, the solutions of Eq. (19) are found to be very close to those of Eq. (18) both for molecular systems³³ and for the UEG.³¹ 

For UEG in the high density limit, it is interesting to note that both the variational equation [Eq. (18)] and the transcorrelated equation [Eq. (19)] lead to the same final equation,^{29,31,34–36}

where T₂(k) is a volume factor defined by

with Θ being the Heaviside step function. In the thermodynamic limit, T₂(k) can be calculated as

The correlation factor optimized for the Hartree–Fock reference is, then, solved as

With this optimized factor $ũ$⁠, we can build the transcorrelated Hamiltonian (12), which is, then, split into two parts for a perturbation treatment,

where for simplicity, we take only the kinetic terms in $Ĥ0$⁠,

With this partition of the TC Hamiltonian, the perturbation energies can be written as

where E₀ is simply the total kinetic energy and E₁ contains the exchange energy $⟨Φ0|Ŵ|Φ0⟩$ and the first order correlation energy $⟨Φ0|K̂+L̂|Φ0⟩$⁠. The second order energy contributes only to the correlation energy.

The first order correlation energy is composed of the expectation values of the two body operator $K̂$ and the three body operator $L̂$⁠. For calculations of these values, we need to find all possible contractions of the creation and annihilation operators. From the two body operator, we get a direct contribution and an exchange contribution,

Here, only the $F(▿u)2$ term of $K̂$ makes a non-vanishing contribution, while the other terms coming from $[Ĥ,τ]$ have no contribution to the expectation value. The expectation value of the three body operator also has two possible contractions,

The volume factor O₃ is defined as

which can be calculated analytically in the thermodynamic limit.³¹ 

In the thermodynamic limit, the summation should be replaced by an integration

This replacement will be used in the following context without any notification.

The first term in Eq. (31) is partly canceled by the first term in Eq. (30), and we finally have three terms in the expression of the first order correlation energy,

For a more transparent treatment, we reformulate expression (23) of $ũ(k)$ as

The expression of the first order energy contains only quadratic terms of $ũ(k)$⁠, and the prefactor αr_s of $ũ$ in Eq. (38) will be canceled by the prefactor $1(αrs)2$ in the final energy expression according to Eq. (10). If r_s is set to be 0 in the remaining integrands in the first order energy, it is found that $E11c$ will diverge, while the expression of $E12c$ and $E13c$ still gives finite results. The singular term in $E11c$ originates from the integration over small k. In order to get an analytic result of the singular term, we rearrange the expression of $E11c$ as

In the first line of expression (40), T₂(k) is simplified with its leading term $34k$ in the small k region, and this simplified integrand is subtracted out from the original integrand in the second and third lines. The simplified integral in the first line can be calculated analytically as

In the second and third lines, the singularities of the two terms cancel each other, and hence, they contribute at most only to the constant term in the correlation energy. For the calculation of the contribution to the constant, we can simply take out the common $α2rs2$ factor and set r_s = 0 in the remaining integrand, which gives

In a similar way, the integral in the last line can be easily calculated as

$E12c$ and $E13c$ also contribute at most only to the constant term of the correlation energy and can also be simplified with the above method. After this simplification, it is still too difficult to be calculated analytically due to the complicated integrand. There is, however, no problem for a numerical calculation of such 1D integrals, which leads to

By adding all the results together, we get the final result for the first order correlation energy,

The first order correlation energy is essentially the correlation energy produced by the single determinant Jastrow ansatz. It is interesting to find that the leading logarithmic term $932π2ln(rs)$ is the same as that obtained by Talman based on the linked cluster expansion.³⁶ This logarithmic term makes up roughly 90% of the exact RPA logarithmic term $1−ln⁡2π2lnrs$⁠, which is also a typical portion of correlated energy one may expect from variational quantum Monte Carlo calculations based on the single determinant Jastrow ansatz.

The transcorrelated effective potential $Ŵ1$ depends on r_s nonlinearly. Besides the original Coulomb potential w(k), which is proportional to r_s, there are some terms that rely on $ũ(k)$ linearly and others quadratically. According to Eq. (38), in the high density regime where r_s ∼ 0, $ũ(k)$ scales roughly linearly with r_s. To get an estimation of the contribution to the second order energy E₂ from each quadratic term in the effective potential $Ŵ1$⁠, we can first set r_s = 0 in the denominator of expression (38) of $ũ(k)$⁠,

If this only leads to finite integrals in the calculation of E₂, we can ignore this quadratic term for the current calculation since it gives at most a $O(rs3)$ contribution to E₂ according to Eq. (29). Based on this estimation, we find there is only one quadratic term left, where approximation (47) leads to divergent integrals in the calculation of E₂. This term is given as

which is a contraction of the three body potential $L̂$ in Eq. (14) by taking k = k′. Consequently, for the calculation of the two leading terms in the second order energy E₂, we only need to take the following approximation of the effective potential:

By taking this two body potential approximation of $Ŵ1$⁠, the second order energy consists of two contributions,

obtained by a direct contraction and an exchange contraction, respectively.

The direct contraction term of E₂ can be written as

Q(k, t) and A(k, t) are given as integrals,

Detailed calculations dealing with Q(k, t) and A(k, t) are presented in  Appendixes A and  B. The results are written as

where X(k) is piecewisely defined for k < 2,

and for k ≥ 2,

By using expressions (56) and (57) in Eq. (53), we have

where we find that the last term equals $−E11c/N$⁠. Therefore, we have

where we have used a simple relation

which can be obtained from Eqs. (51) and (20). As far as the singularity at k ∼ 0 is concerned, the first integrand is similar to the integrand of $E11c$⁠,

We can reformulate the expression of $E21c$ as

where the remaining integrals are not any more singular at k ∼ 0 in the high density limit. Then, we can use the simple expression (47), and the integrals can be numerically evaluated as