Finite elements and the discrete variable representation in nonequilibrium Green’s function calculations. Atomic and molecular models
In this contribution, we discuss the finite-element discrete variable representation (FE-DVR) of the nonequilibrium Green’s function and its implications on the description of strongly inhomogeneous quantum systems. In detail, we show that the complementary features of FEs and the DVR allows for a notably more efficient solution of the two-time Schwinger/Keldysh/Kadanoff-Baym equations compared to a general basis approach. Particularly, the use of the FE-DVR leads to an essential speedup in computing the self-energies.
As atomic and molecular examples we consider the He atom and the linear version of H in one spatial dimension. For these closed-shell models we, in Hartree-Fock and second Born approximation, compute the ground-state properties and compare with the exact findings obtained from the solution of the few-particle time-dependent Schrödinger equation.
In the last decade, the application of the nonequilibrium Green’s function (NEGF) to describe strongly inhomogeneous quantum systems has started to become an actively considered subject. Thereby, various finite and localized systems have challenged attention and different many-body approximations have been applied. Recent state-of-the-art approaches discuss small atoms and molecules [1, 2], few-electron quantum dots [3, 4] and quantum dots coupled to leads , molecular junctions , and Hubbard nanoclusters . In part, these works also include the monitoring of the system’s temporal evolution which, accounting for correlation and memory effects, requires to solve the two-time Schwinger/Keldysh/Kadanoff-Baym equations [8, 9, 10] (SKKBE).
To solve the SKKBE for homogeneous quantum systems [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22] has become routine. However, this still does not hold for finite, localized and inhomogeneous systems. The reason is, that in contrast to homogeneous systems where one spatial coordinate or momentum drops out of the NEGF, the inhomogeneity claims adequate resolution in both coordinates or momenta. To meet these requirements all above mentioned works rely on an expansion of the NEGF in one-particle orbitals: for atomic and molecular systems, linear combinations of Slater-type or Gauss-type orbitals are being used (also in connection with tight-binding models), whereas for other classes of problems e.g. potential-eigenstate basis functions are being utilized. Beyond a basis ansatz, alternatives such as direct grid (finite-difference) methods do not exist as they are computationally very expensive. Nevertheless, also basis approaches are so far limited to relatively small basis sets. Generally, the numerical complexity involved in the description of the binary interactions does not permit to extend nonequilibrium calculations to much larger basis dimensions than a guide number of orbitals. This explains the low resolution in the description of photoionization processes of model atoms [23, 24] using nonequilibrium Green’s functions.
In this contribution, we—for the first time—develop a grid-based approach in the frame of the finite-element discrete variable representation (FE-DVR), e.g. Refs. [25, 26, 27, 28] and references therein. This method leads to specially-designed, flexible basis sets which are capable to combine the advantages of pure grid and standard basis approaches, see Sec. 2 and also Ref. . In particular, it allows for a very efficient treatment of the binary interactions and, in turn, a drastic simplification of higher-order self-energy expressions, which, generally, require the main effort in all NEGF calculations. As a result, only semi-analytical matrix elements of the interaction energy operator are required to compose the second Born self-energy. This has to be compared to a general basis representation of the NEGF: There, matrix elements are involved, cf. . As a consequence, the use of the FE-DVR enables more efficient calculations at less storage memory and computing time and provides the basis to consider also spatially extended hamiltonians, where particles may occupy large domains in coordinate space.
After the outline of the method, we, in Sec. 3, apply the FE-DVR representation of the NEGF to compute the ground-state properties of atomic and molecular models. First, in Sec. 3.1, we discuss the one-dimensional (1D) helium atom [30, 31, 32, 33, 34, 35, 36, 37], where we focus on technical details such as grid size, DVR basis size, and convergence in the case of the Hartree-Fock and second Born approximation. The benchmarking results shown are also of relevance for time-dependent calculations, as they define and border the requirements to resolve explicit correlation effects (within a NEGF approach) e.g. the two-electron resonances  which are embedded in the one-electron continuum of the atom and are going to be occupied during laser-atom interactions—see also  in the present volume. Finally, we consider the linear molecular ion H in the symmetric 1D singlet configuration [40, 41] as an example of a two-electron molecule and vary the interatomic distance to record the binding-energy curves in Hartree-Fock and second Born approximation, see Sec. 3.2. From this we can extract the minimum ground-state energies and the respective bond-lengths, which are compared to the findings from the few-particle time-dependent Schrödinger equation.
2 The nonequilibrium Green’s function in FE-DVR representation
For the description of the NEGF using finite elements together with the discrete variable representation  (DVR), we consider the general -electron Hamiltonian
with the kinetic energy , the time-dependent potential energy , and the pair interaction energy formulated in atomic units. The one-particle nonequilibrium Green’s function with space-time arguments and reads
with addition of its adjoint equation. Further, is the one-electron (kinetic plus potential) energy, denotes the self-energy, and equilibrium initial correlations are treat in the mixed Green’s function approach [43, 44, 45], cf. Sec. 2.3.
The favorable aspects of the FE-DVR representation have been successfully used for the time-dependent Schrödinger equation (TDSE) by Rescigno  and others, e.g. . There, the accuracy of the DVR  on the one hand and the sparse character of FEs on the other led to an efficient TDSE code which is well parallelizable. In our case, this hybrid approach allows us to rewrite the SKKBE (3) in a highly effective matrix notation using optimal combinations of a grid and a local basis.
2.1 Grid-based ansatz
On a predefined spatial 1D interval , we expand the nonequilibrium Green’s function of hamiltonian (1) as
with time-dependent complex coefficients and real basis functions . Outside the interval , we assume the NEGF to vanish. All indices (superscripts), ranging in Eq. (4), are linked to a grid of length composed of finite elements, see Fig. 1. All indices (subscripts), ranging , are connected to locally defined basis functions which are constructed as follows: First, we divide the interval into finite elements with boundaries . To each FE , we then attach a local DVR basis using the generalized Gauss-Lobatto (GGL) points  and weights :
with the standard Gauss-Lobatto points (and weights ). For the special case of Legendre interpolating functions, the points are defined as roots of the first derivative of Legendre polynomials according to
and the weights are given by
where denotes the total number of basis functions per element. In our approach, we keep the number of basis functions constant in each FE . However, in order to combine the locally defined DVR basis functions to a continuous basis set, the FE-DVR space is being spanned by two classes of functions—’bridge’ and ’element’ functions, see also Ref. . The bridge function (the case in Eq. (9)) extends over two adjacent FEs and ensures communication between the grid domains and . In particular, it guarantees the continuity of the NEGF. The ’element’ functions are zero at and outside the respective element boundaries. Using definitions (5) and (6), the basis functions have the explicit form
and obey . As in the last finite element (i.e. in the element ) no bridge function is needed due to the boundary condition of vanishing outside the interval , the total FE-DVR set consists of
2.2 Matrix elements of the kinetic, potential and interaction energy operators
With representation (4) of the NEGF, the SKKBE will transform into an equation of motion for the matrix with elements , cf. Sec. 2.3. Obviously, this equation involves also the kinetic, potential and interaction energy operator of Eq. (1) in matrix form, which we specify in the following. Thereby, integrations over coordinate space are performed by using the generalized Gauss-Lobatto quadrature rule, and case differentiations arise from the basis functions being split into element and bridge functions.
First, let us consider the potential and the kinetic energy in FE-DVR representation: The potential-energy matrix is given by
This implies that the potential energy is diagonal with respect to elements and local DVR basis indices , and that, consequently, it can be represented by a vector of dimension . Moreover, Eq. (2.2) holds true also for any other local operator. As the operator of the kinetic energy is non-local in coordinate space, the matrix elements are not diagonal. We follow the derivation of Ref.  and obtain the kinetic-energy matrix as
Here, the case differentiations lead to the matrix having a block-diagonal form , and the quantity is given by
which involves the first derivatives of the Lobatto shape functions at the GGL points, see also Ref. .
Next, let us focus on the matrix elements of the binary-interaction operator (the two-electron integrals) which are carrying a set of four index-pairs and are defined by
In a general basis approach, the two-electron integrals [Eq. (16) with all index-pairs replaced by single indices] often require a careful analysis, as they are not analytically accessible and have to be numerically precomputed for all combination of indices, e.g. . Although symmetry relations  help to restrict oneself to a subset of indices, the effort scales with and thus can be huge for larger basis sets. On the contrary, using the FE-DVR, the evaluation of the two-electron integrals turns out to be much simpler. In particular, the integrals can be performed in a semi-analytical way such that Eq. (16) reduces to
where the kernel matrix is symmetric and follows as
To obtain Eq. (18), we have used the separable form of the discretized interaction potential , and, correspondingly, the quantities denote the eigenvalues of the matrix
and are related to the eigenvectors :
In comparison with any single-electron matrix element (such as the kinetic or potential energy), in FE-DVR, the calculation of the binary-interaction matrix elements involves just an additional but numerically elementary matrix diagonalization. Furthermore, besides the fact that with Eq. (17) the two-electron integrals attain a very simple form independent of the specific pair-interaction potential, only a single matrix of dimension needs to be stored in the code. This memory-friendly property is based on the high degree of diagonality determined by the product of Kronecker deltas and represents a main attractive feature of the FE-DVR representation. Particularly, this aspect opens the way towards efficient NEGF calculations, since Eq. (17) has direct consequences for the structure of the self-energies, see the following Section.
2.3 Equations of motion
Once all relevant matrix elements are known with respect to the chosen FE-DVR basis (9), we can start to solve the Schwinger/Keldysh/Kadanoff-Baym equations (3) for the one-particle Green’s function expanded in the form of Eq. (4). This implies, though, the SKKBE in the matrix form of the finite-element discrete variable representation:
where we have denoted , , and Eq. (21) has to be supplied with its adjoint equation. Further, Eq. (22) separates the self-energy matrix into a time-local Hartree-Fock part () and a contribution that accounts for electron-electron (-) correlation and memory effects. However, as an exact treatment of - correlations is impractical, we have to apply many-body approximations for which the second Born diagrams provides one of the most basic models. Hence, besides the general form of the HF self-energy
where denotes we, in second Born approximation, have
Eqs. (23) and (24) involve the spin degeneracy factor , the matrix elements of Eq. (17) and show a very simple form compared to the situation when a general basis is applied, e.g. . The reason for this is the subtle structure of the FE-DVR. In detail, the Hartree term is completely diagonal [including a single sum over elements] and the exchange term involves only a product of two matrix elements. For the second-order Born terms, the degree of simplification is even more drastic: In a general basis representation, two sums are required for each full vertex point in the second-order self-energy diagrams and, additionally, a single sum is needed for the start- and the end-point. This leads to a scaling with . However, in our case, due to the diagonality of the two-electron integrals, cf. Eq. (17), the evaluation reduces remarkably to a scaling with per matrix element.
For the atomic and molecular model calculations to be outlined as first benchmarks of the method in Sec. 3, we, in this contribution, restrict ourselves to the (equilibrium) ground-state properties. For completeness, we give a short summary of the main computational steps. First, in the FE-DVR picture, the Hartree-Fock equilibrium Green’s function, denoted , follows from
where is the HF approximation of with and is the inverse temperature, and is defined via Eq. (23), but with being replaced by matrix . Thereby, Eq. (25) has to be solved to self-consistency by iteration and as result we obtain
Here, the vector contains the energy eigenvalues of , the matrix summarizes the corresponding eigenvectors, and the chemical potential is determined by normalization of the Fermi distribution: . Eq. (26) solves a simple differential equation, see e.g. , and corrections due to - correlations in second Born approximation are obtained by insertion into the Dyson equation  (the SKKBE in the limit ) for the full Matsubara Green’s function :
where with instead of in Eqs. (23) and (24). The convolution integrals in Eq, (27) are performed in sequence by direct integration. This has been found to be more stable and controllable than the method applied before in Refs. [2, 4]. Once the self-consistent is computed, we have direct access to many observables, e.g. the one-electron density is obtained as . The total energy is computed similar as in Refs. [43, 4]. Overall, in order to ensure the atoms and molecules being in the ground state, we set the inverse temperature .
3 Performance for atomic and molecular models
As first benchmarks and preparatory work for the investigation of the temporal evolution of small atoms and/or molecules following an external (e.g. laser-induced) perturbation, we here consider their equilibrium-(initial-)state preparation within the FE-DVR context. As examples, we focus on the He atom and the linear molecular ion H, modeled in one spatial dimension. For both two-electron systems, the Coulomb potential is considered in the regularized form which, from the physical point of view, allows for a transverse extension of the few-particle wave function. Furthermore, - correlations are treated in second Born approximation.
3.1 The 1D helium atom
The helium atom is the most elementary closed-shell two-electron system. In one spatial dimension, it is well modeled by the nucleus potential , where the atomic number is , , and is a regularization parameter. In this setup, the 1D helium atom serves as a fundamental ’testing ground’ for multi-electron calculations [36, 37, 39] and provides many features of the single- and double-ionization dynamics including the so-called knee structure [34, 35]. Considering the singlet state, we refine the model and set in Eqs. (23) and (24), and used .
Further, for the equilibrium calculations, we have used a FE-DVR basis that covers a total interval of a.u. length. This corresponds to a domain that is about times larger than the characteristic extension of the ground-state wave function or density, cf. Ref. . Such a grid extension is more than adequate to resolve the ground-state features of the model (to be discussed here) but will become crucial, when the helium atom is perturbed by external fields and, in turn, electrons start to occupy highly excited or continuum states. In this sense, our results are benchmarks also with relevance for the computation of the system’s temporal evolution. In particular, a grid with an extension of about a.u. should be well capable to resolve the two-electron resonance states embedded within the one-electron continuum of dipole spectra. This follows from the performance of the few-particle time-dependent Schrödinger equation (TDSE) using absorbing potentials, e.g.  and references therein. We note, that also in the FE-DVR approach, such imaginary one-electron potentials that damp reflections at the interval boundaries are easily implemented, just allowing the matrix elements in Eq. (2.2) to be complex.
The explicit partitioning of the interval into finite elements has been organized as follows: (case I) the interval is divided into equidistant segments, (case II) the central FE has a width of one atomic unit and the width of the surrounding elements is linearly increasing towards the interval boundaries. The effect of these segmentations on the Hartree-Fock ground-state energy convergence of the helium atom is displayed in Fig. 2 a) as function of the total basis size. In case I, for equidistant FEs, the convergence is relatively slow with and strongly depends on the number of elements as well as on the number of local DVR basis functions used (see the different symbols and lines). Consequently, more than functions are needed for the HF energy to deviate by less than Ha from the converged HF result Ha. In case II, the situation is completely different as essentially more basis functions are available in the center region of the interval. This enables a superior representation of the Matsubara Green’s function in coordinate space, and leads to adequate convergence at to FE-DVR basis functions with an error reduced by several orders of magnitude compared to case I. Furthermore, the ground-state energy depends less on the number of elements. For completeness, we note that, besides the total energy, an additional indicator for the basis quality, is to look at how well the potential can be expanded into the chosen FE-DVR basis.
For case II with and , we have computed the ground-state energy in second Born approximation. Thereby, all quantities are found to be well converged with respect to the basis size, compare with the HF case in Fig. 2 a). The time-argument in the Matsubara Green’s function () has been discretized using a uniform power mesh  with parameters and —for definition see e.g. Refs. [4, 43]. Fig. 2 b) indicates the convergence with respect to these parameters, where the total number of -grid points is given by . At , a mesh parameter ensures the particle number being sufficiently stable during iteration of the Dyson equation (27). Particularly, with more than grid points it is possible to compute the ground-state energy to relatively high precision, Ha. As result, the electron-electron correlations lower the total energy accounting for % of the correlation energy and, hence (improving the HF result), approaches the exact ground-state energy which is Ha. For the discussion of other observables such as the one-electron density in HF and second Born approximation, the reader is referred to Ref. .
3.2 The linear molecule H
With more than two nuclei, the molecular ion H can, in 1D, only be realized in its linear version , where the H-H bonds are oriented parallel to each other. Hence, the electrons move along the molecular axis and the one-electron potential is modeled as
where denotes the interatomic distance, and the last term (offset) collects all internuclei interaction energy contributions. For the ground-state NEGF calculations, the coordinate space has been constrained to an interval of a.u. with a grid of FE being linearly increasing, starting from a one atomic unit wide central element. In total, FE-DVR basis functions have been used.
The total binding energy of the singlet state is shown in Fig. 3 a) against distance for the exact solution of the TDSE (dotted [and triangles]), the Hartree-Fock (dashed) and the second Born approximation (solid). Over a broad range of internuclear distances, the second Born approximation, thereby, accounts for about -% of the correlation energy, and—leading overall to a larger bond—essentially improves the HF result.
Furthermore, with a value of Ha, H has the same dissociation threshold as the 1D hydrogen molecule , but, in addition, leaves behind a positively charged hydrogen ion. However, we note, that the Hartree-Fock and the second Born approximation, cannot resolve this threshold. The reason is that the H molecule dissociates into open-shell fragments—two hydrogen atoms and a single hydrogen ion. These cannot be represented within a spin-restricted calculation (with ), and, hence, lead to a strong deviation of in the limit of large . Regardless of this failure of the ansatz, the NEGF calculations are well capable to describe the behavior of around the minima in the binding energy curves. Also, the equilibrium positions are consistently reproduced, and the correct trend is observed when - correlations are being included: the bond-length shifts to larger nuclear separations, for the specific values obtained see caption of Fig. 3.
The electron ground-state density of the molecular ion is displayed in Fig. 3 b) together with schematic curves for the spatial potentials , where is twice the equilibrium internuclear distance, cf. Eq. (28). In HF approximation (dashed curve), the density shows a pronounced maximum in the central region of the molecule, whereas the exact density (dotted curve) is essentially less peaked. However, also the exact result does not indicate onset of electron localization, i.e. does not show separated maxima in . In second Born approximation, corresponding to a lower total energy (cf. Fig. 3 a)), we resolve the correct trend of this density reduction. Moreover, replacing the self-consistent bond-lengths by the corresponding length obtained from the TDSE (dash-dotted curves) only slightly improves the results for the Hartree-Fock and the second Born approximation. This explains that the substantial differences in the density profiles are unambiguously due to - correlation effects.
The FE-DVR ansatz (4) provides an elegant and very efficient way to treat binary interactions in NEGF calculations for inhomogeneous quantum systems. To this end, the method uses a flexible combination of grid (FE) and basis (DVR) strategies, which allows for simple, (semi)-analytical formulas for the matrix elements of the kinetic-, potential- and, especially, the interaction-energy operator, cf. Sec. 2.2. Further, for the most basic model that accounts for particle-particle correlations—the second Born approximation—the use of a FE-DVR basis enables remarkable scaling properties: Only summations are required for the computation of a single matrix element of the second-order self-energy instead of summations that are needed in a general basis approach.
Also, we emphasize that the FE-DVR space can, e.g. via finite-element variations, be well adjusted to the problem considered and, thus, allows for an efficient expansion of the NEGF [also for spatially extended hamiltonians] and quick convergence of the observables of interest. This has been exemplified for the He atom and the triatomic molecule H in Sec. 3.1 and 3.2. In addition, the method is found to be stable also for large grids [large basis sets with ] and also for larger particle numbers , considering interacting fermions in a harmonic trap potential.
Finally, we believe that the FE-DVR method is attractive also for other classes of many-body approximations, such as or -matrix calculations, as it will likewise simplify the computation of self-energy contributions of higher than second order. Moreover, though we, in this contribution, focused on the solution of the Dyson equation, the formalism presented is, in particular, well applicable in nonequilibrium situations solving the full two-time SKKBEs which will be demonstrated in a forthcoming publication.
-  A. Stan etal, Europhys. Lett. 76, 298 (2006).
-  N.E. Dahlen, and R. van Leeuwen, Phys. Rev. Lett. 98, 153004 (2007).
-  K. Balzer, and M. Bonitz, J. Phys. A: Math. Theor. 42, 214020 (2009).
-  K. Balzer, M. Bonitz, R. van Leeuwen, A. Stan, and N.E. Dahlen, Phys. Rev. B 79, 245306 (2009).
-  P. Myöhänen, A. Stan, G. Stefanucci, and R. van Leeuwen, Europhys. Lett. 84, 67001 (2008).
-  K.S. Thygesen, Phys. Rev. Lett. 100, 166804 (2008).
-  M. Puig von Friesen, C. Verdozzi, and C.-O. Almbladh, Phys. Rev. Lett. 103, 245306 (2009).
-  P.C. Martin, and J. Schwinger, Phys. Rev. 115, 1342 (1959).
-  L.V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1515 (1964) [Sov. Phys. JETP 20, 235 (1965)].
-  L.P. Kadanoff, and G. Baym, Quantum Statistical Mechanics (Benjamin, Inc., New York, 1962).
-  P. Danielewicz, Ann. Phys. (N.Y.) 152, 305 (1984).
-  H.S. Köhler, Phys. Rev. C 51, 3232 (1995).
-  P. Bożek, Phys. Rev. C 56, 1452 (1997).
-  N.H. Kwong, and M. Bonitz, Phys. Rev. Lett. 84, 1768 (2000).
-  M. Bonitz, D. Kremp, D.C. Scott, R. Binder, W.D. Kraeft, and H.S. Köhler, J. Phys.: Cond. Mat. 8, 6057 (1996).
-  D. Semkat, D. Kremp, and M. Bonitz, Phys. Rev. E 59 1557 (1999).
-  D. Semkat, D. Kremp, and M. Bonitz, J. Math, Phys. 41, 7458 (2000).
-  R. Binder, H.S. Köhler, and M. Bonitz, Phys. Rev. B 55, 5110 (1997).
-  N.H. Kwong, M. Bonitz, R. Binder, and S. Köhler, phys. stat. sol. (b) 206, 197 (1998).
-  W. Schäfer, J. Opt. Soc. Am. B 13, 1291 (1996).
-  L. Banyai, H. Haug, and P. Gartner, Eur. Phys. J. B 1, 209 (1998).
-  P. Gartner, J. Seebeck, and F. Jahnke, Phys. Rev. B 73, 115307 (2006).
-  D. Hochstuhl, K. Balzer, S. Bauch, and M. Bonitz, arXiv:0902.0768, Physica E, in press (2009).
-  M. Bonitz, D. Hochstuhl, S. Bauch, and K. Balzer, arXiv:0909.1964, submitted to Contributions to Plasma Physics (2009).
-  T.N. Rescigno, and C.W. McCurdy, Phys. Rev. A 62, 032706 (2000).
-  L.A. Collins, S. Mazevet, J.D. Kress, B.I. Schneider, and D.L. Feder, Physica Scripta T110, 408 (2004).
-  B.I. Schneider, L.A. Collins, and S.X. Hu, Phys. Rev. E 73, 036708 (2006).
-  J. Feist, R. Pazourek, S. Nagele, E. Persson, B.I. Schneider, L.A. Collins, and J. Burgdörfer, J. Phys. B: At. Mol. Opt. Phys. 42, 134014 (2009).
-  K. Balzer, S. Bauch, and M. Bonitz, arXiv:0910.5458, submitted to Phys. Rev. A (2009).
-  M.S. Pindzola, D.C. Griffin, and C. Bottcher, Phys. Rev. Lett. 66, 2305 (1991).
-  R. Grobe and J.H. Eberly, Phys. Rev. A 48, 4664 (1993); S.L. Haan, R. Grobe, and J.H. Eberly, Phys. Rev. A 50, 378 (1994).
-  D. Bauer, Phys. Rev. A 56, 3028 (1997).
-  W.-C. Liu, J.H. Eberly, S.L. Haan, and R. Grobe, Phys. Rev. Lett. 83, 520 (1999).
-  N.E. Dahlen, and R. van Leeuwen, Phys. Rev. A 64, 023405 (2001).
-  M. Lein, E.K.U. Gross, and V. Engel, Phys. Rev. Lett. 85, 4707 (2000).
-  J. Zanghellini, M. Kitzler, T. Brabec, and A. Scrinzi, J. Phys. B: At. Mol. Opt. Phys. 37, 763 (2004).
-  M. Ruggenthaler, D. Bauer, Phys. Rev. Lett. 102, 233001 (2009).
-  G. Tanner, K. Richter, and J.-M. Rost, Rev. Mod. Phys. 72, 497 (2000).
-  D. Hochstuhl, and M. Bonitz, proceedings of conferencne ”PNGF IV”, this issue (2009).
-  I. Kawata, H. Kono, A.D. Bandrauk, Phys. Rev. A 64 043411 (2001).
-  N. Suzuki, I. Kawata, and K. Yamashita, Chemical Physics 338, 348-353 (2007).
-  J.C. Light, I.P.Hamilton, and J.V. Lill, J. Chem. Phys. 82, 1400 (1985).
-  N.E. Dahlen, and R. van Leeuwen, J. Chem. Phys. 122, 164102 (2005).
-  N.E. Dahlen, A. Stan, and R. van Leeuwen, J. Phys: Conf. Ser. 35, 324 (2006); N.E. Dahlen, R. van Leeuwen, and A. Stan, J. Phys: Conf. Ser. 35, 340 (2006).
-  A. Stan, N.E. Dahlen, and R. van Leeuwen, J. Chem. Phys. 130, 224101 (2009).
-  J.C. Light, and T. Carrington Jr., Advances in Chemical Physics 114, 263 (2007)
-  The generalization to different numbers of basis functions per element is straightforward, and only slightly alters the matrix elements involved.
-  D.E. Manolopoulos and R.E. Wyatt, Chem. Phys. Lett. 152, 23 (1988).
-  The two-electron integrals are symmetric with respect to interchange of , and tupels .
-  S. Bauch, and M. Bonitz, Phys. Rev. A 78, 043403 (2008).
-  W. Ku and A.G. Eguiluz, Phys. Rev. Lett. 89, 126401 (2002).
-  J.D. Alexander, C.R. Calvert, R.B. King, O. Kelly, L. Graham, W.A. Bryan, G.R.A.J. Nemeth, W.R. Newell, C.A. Froud, I.C.E. Turcu, E. Springate, I.D. Williams, and J.B. Greenwood, J. Phys. B: At. Mol. Opt. Phys. 42, 141004 (2009).