The Hypernetted Chain Equations for Periodic Systems

Starting from the general inhomogeneous Fermi hypernetted chain equations, the equations for periodic systems are derived by simple Fourier transform. It is shown how the symmetry reduces the size of the involved quantities. First results for a one-dimensional (1D) model system are presented. The results allow a reliable estimation of the numerical demand even for realistic 3D systems, such as solids. It is shown that treatment of this systems is feasible with moderate computational resources.


Introduction
Jastrow correlated trial wave functions, i.e., a wave function which additionally to a Slater determinant includes two-particle correlations, are widely used in quantum Monte Carlo (MC) techniques [1,2]. Less well known are the analytic methods to calculate expectation values with this wave function [3][4][5][6]. The Fermi hypernetted chain method (FHNC) is an effective scheme to calculate a large class of cluster diagrams. The optimal correlation function is obtained by minimization of the energy expectation value. A drawback of the method, compared to MC techniques, is that certain diagrams, i.e., the elementary diagrams, are not covered by the scheme. Therefore, the obtained quantities and matrix elements are only approximations. However, the advantage of the method is that it is generally numerically less demanding than MC B Martin Panholzer martin.panholzer@polytechnique.edu 1 Laboratoire des Solides Irradies, Ecole Polytechnique, CNRS-CEA, Universite Paris-Saclay, 91128 Palaiseau cedex, France methods. The analytic representation of the wave functions further allows to deal with excited states [7,8].
An additional advantage of the used formulation of the FHNC method is that the functional form of the method allows a parameter-free optimization of the two-particle correlation function.
The FHNC method has been first developed for homogeneous systems, but has been generalized to inhomogeneous systems by Krotscheck et al. [9][10][11][12]. In order to keep the numerical demand low, one has to utilize the symmetries of the system. This has been done for slab [13] and for spherical geometries.
In this paper a possible discretization for periodic systems is presented which reduces the numerical demand considerably and leads to a simple form of the Euler equation. A different discretization and coordinates are proposed by Krotscheck [14]. Although it allows a further reduction in numerical demand by reducing the resolution of the center of mass coordinate, it sacrifices the simplicity of the relevant equations and reduces the resolution and consequently the accuracy of the result.
This work is related to an implementation of Fantoni and Schmidt [15], which formulates and solves the FHNC equations in a periodic box, but with a uniform density. The presented method extends this approach by explicitly including a nonconstant density.

The Inhomogeneous HNC Equations
Here the generalization to the general inhomogeneous case of the FHNC equations in the formulation of Kallio and Piilo [16] is given. The reason for this is that results are shown for a 1D model system and the simplified FHNC of Krotscheck doesn't work in 1D.
The inhomogeneous theory yields coupled equations for the one-particle and twoparticle correlations. The one-particle equation is a generalized Hartree-Fock equation (gHF) as given in [12], and the two-particle equation is actually a set of equations given below. These two equations are coupled and need to be solved self-consistently. However, it seems to be a reasonable approximation for certain systems to omit the coupling [13]. This is the justification for using a model for the single-particle states in the next section and not solving the gHF.
The inhomogeneous Euler-Lagrange equations can be written as with the definition of the convolution and where S is the, to be determined, static structure function. V ph is the particle hole irreducible interaction explicitly given below. The one-body operator H 1 is defined as where ρ 1 is the one-body density obtained from the gHF. The induced interaction is where the V ph from the previous iteration is used. (In order to start the iteration process, a reasonable guess for V ph is needed.) With these ingredients and by using the pair distribution function g(r 1 , r 2 ), the new particle-hole irreducible interaction is computed where ∇ 1 is short for the gradient with respect to r 1 and v is the interaction potential of the particles in the system. Note that here no convolutions are involved, i.e., the functions are simply multiplied. The Fermi potential v F is defined as with The free static structure function S F is given by where l is the non-interacting density matrix as result of the gHF and ν is the spin degeneracy and g F is the free pair distribution function analog of S F . In order to obtain the boson version of the equations v F has to be set to zero. (1) is solved by considering the eigenvalue problem

Solution of the Euler equation: Equation
The ψ l are orthonormalized according to The static structure factor is obtained by and its inverse

Definition of the Fourier Transform
In order to describe a periodic system a unit cell, spanned by the lattice vectors a 1 · · · a d , where d is the dimensionality of the system, is defined. One-body quantities like the density are periodic: ρ(r + T n ) = ρ(r). With the translation vector T n = n 1 a 1 · · · + n d a d . A similar relation also holds for two-body quantities, e.g., the pair distribution function: g(r 1 + T n , r 2 + T n ) = g(r 1 , r 2 ). This symmetry is used to constrain the first coordinate to the basic unit cell. This is denoted by a bar: g(r 1 , r 2 ). Further, we may decompose the second coordinate: g(r 1 ,r 2 + T n ). In order to describe pair properties properly we have to choose a cutoff in the translations. Therefore, we restrict the integers n i to the interval 0 · · · N − 1, 1 which we call now crystal on which we impose periodic boundary conditions. (Although this is non-physical, it becomes a good approximation if N is large enough.) This crystal consists of N d unit cells. The unit cell is discretized with n G points in each direction, so the resolution is r i = a i n G . A periodic function can be expressed as a sum of plane waves (13) where G = m 1 b 1 · · · + m d b d are the reciprocal lattice vectors, with By utilizing the periodicity of the system the Fourier transform is defined as Crystal and its inverse f (r 1 , r 2 ) = q,G,G e −i(q+G)r 1 e i(q+G )r 2 f (q, G, G ) (15) With this definition the convolution as defined above becomes in momentum space r 1 , r 3 )] (q, G 1 , G 3 With the above definitions the Euler equation is written in reciprocal space as The induced interaction is also calculated most effectively in reciprocal space, by Fourier transform or Eq. (4). However, the particle hole irreducible interaction is calculated in r-space, like in the homogeneous case. The most time-consuming part is the solution of the Euler equation. The Fourier transform reduces the eigenvalue problem Eq. (9) considerably since q appears as a parameter: so the computation timescales are linear with the number of q-vectors.

Results for a 1D Model System
The described method is applied to a 1D model system similar to that described by Asgari [17]. The interaction potential is derived from the electron gas in a homogeneous trap v(r ) = e 2 √ 2 2b exp r 2 4b 2 erfc |r| 2b where b characterizes the thickness of the wire. Since the single-particle equation is not solved, a simple model for the single-particle basis is used These single-particle states result in a sinusoidal modulated density with the amplitude determined by λ. G = 2π a is the reciprocal lattice vector and a the length of the unit cell.
Results of this model are plotted in Fig. 1 with the parameters a = 6r S a B , b = 0.1r S a B , λ = 0.02 and the Wigner-Seitz radius r S = 1 and a B = 0.529Å the Bohr radius. In order to converge the result a real space resolution of r = 0.1r S a B was sufficient for the homogeneous and the inhomogeneous implementation. This results in n G = 60 reciprocal lattice points, and the number of q-points needed is n q = 20 which determine the size of the large cell. It has been verified that the implementation reproduces the FHNC/0 result of [17] in the homogeneous limit. In Fig. 1 the deviation of the inhomogeneous pair distribution function from the homogeneous on is clearly seen. Further, the minimum in g(r, r ) at r = r follows approximately a local density approximation. The number of G vectors n 3 G for a realistic 3D system is now estimated. Sodium is used as an example. The lattice constant of the primitive unit cell is a = 1.75r s a B , where r s = 4 for sodium. From the treatment of the homogeneous electron gas and the results for the 1D model system it is seen that n G = 20 points are sufficient for that length. It is even possible to reduce this number further. This results in a 8000 × 8000 eigenvalue problem for each q-vector, which is numerically tractable. Only q-vectors in the irreducible Brillouin zone need to be calculated, which further reduces the numerical load.

Conclusions
The FHNC equations for periodic systems have been derived by simple Fourier transform of the general inhomogeneous equations. This considerably reduces the numerical demand so that realistic systems become numerically feasible. Further reduction in the numerical demand is possible, but only at the cost of simplicity of the theory.
First results for a 1D system have been shown. There is increasing interest in the theoretical description of 1D systems, e.g., [18][19][20]. The presented method could extend existing methods by applying it to more realistic systems. Also more fundamental questions could be addressed, like the extension of the local density approximation to pair quantities could be investigated. That is how to combine results of the homo-geneous system for different densities to obtain an estimation of the inhomogeneous result. Work in that direction is in progress.