Using QFDFT to obtain ground states of systems

This article has been published in Journal of Physics A, Volume 40 pages 5981-5993, 2007.

**1. Introduction**

In recent years there have been many efforts to calculate the static as well as dynamic (time- dependent) electron density of many-electron systems through the use of a single equation within the framework of density functional theory (DFT) [1,2,3]. Amongst the earliest attempts is the Thomas-Fermi method [4,5] in which the system is considered as a gas of almost-free electrons and the static electron densities are calculated making use of a single equation. However the time-dependent analogue of Thomas-Fermi method cannot be constructed with a single equation. Such a single equation formulation has been made possible by combining the principles of Quantum Fluid Dynamics (QFD) and DFT [6]. It has been shown that the two basic QFD equations, viz., the equation of continuity and Euler-type equation of motion, can be merged to yield a single time-dependent generalized nonlinear Schrödinger equation (GNLSE). There have been other time-dependent formulations like the self-interaction free time-dependent DFT [7,8] and time-dependent-current DFT [9], which, however, deal with individual occupied spin-orbitals which automatically render themselves very-demanding computationally, especially when dealing with large systems. The QFD-DFT approach has been subsequently used to obtain the ground-state densities in spherical coordinates [10,11,12], as well as to study dynamic properties of two-electron systems [13,14]. Of late there has also been a spurt of interest in studying the static and dynamic properties of noble gas atoms and clusters, especially the phenomenon of Above Threshold Ionization (ATI) and High Harmonic Generation (HHG) in very intense laser fields [16,17,18]. The QFD-DFT approach has been tested on such systems to study their ground state properties by Roy and Chu [12] using a time-dependent generalized pseudospectral technique. They have subsequently applied their technique to study multiphoton ionization (MPI) and High Harmonic Generation (HHG) of He and Ne atoms in intense laser fields [15]. However their calculations having been done in one dimension, much of the actual picture of the dynamics seems to have been obscured.

In this paper we have recast the GNLSE as a diffusion equation in scaled cylindrical coordinates and evolved it in imaginary time to obtain on convergence the ground state electron densities, energies and other properties of the full spectrum of noble gas atoms from Helium to Xenon. We have compared the results with those obtained from other methodologies and which were performed in one dimension with a view to test the efficacy of our formulation in studying dynamical processes involving the same systems in real time and in three dimensions.

**2. Methodology**

The QFD equations involving the two dynamical variables in 3-dimensional coordinate space, viz., the electron density
(**r**, *t*) and the current density
**j**(**r**, *t*) are expressed as

where
**j**(**r**, *t*) = (**r**, *t*) being the velocity potential. The term
*v*_{(}**r**, *t*) comprises the electron-nuclear attraction potential plus the potential due to any external
field. The fourth term is the classical inter-electronic Coulomb repulsion potential *v*_{el-el}. *G*[] is a universal density functional which
comprises the kinetic, exchange and correlation energy functionals. It can be written thus:a

The first term on the right-hand-side of equation (3) is the Weizsäcker kinetic energy functional. The second term is a correction to the first term to properly account for not only the total kinetic energy of an n-electron (n 2) system but also its global and local behaviour, including atomic shell structure [19]. The remaining two terms are the exchange and correlation energy functionals respectively.

If we define a complex hydrodynamical wave function as

and combine equations (1) and 2 we would arrive at a single time dependent generalized nonlinear Schrödinger equation (GNLSE) [6]:

where
*v*_{eff}(;**r**, *t*) contains the sum of the contributions of all the potential terms.

Equations (1) and 2 can be written in imaginary time instead of real time *t* by substituting
= - *it*. Proceeding as above we find that the GNLSE takes the form of a diffusion type equation:

in real time. *R* in equation (6) is a real diffusion function different from the complex function .

Solving the diffusion equation (6) for sufficiently long time and forcing normalization after each iteration, we obtain the minimum energy ground state of the system characterized by the effective potential
*v*_{eff}(;**r**, *t*).

*v*_{eff}(;**r**, *t*) can be written as the sum of different potential energy terms as follows:

*v*_{el-el} = *E*_{el-el}/ is the inter-electronic coulomb repulsion potential, *E*_{el-el} being the corresponding energy.

*v*_{nu-el} = *E*_{nu-el}/ is the nuclear-electron attraction potential, *E*_{nu-el} being the corresponding energy.

*v*_{x} = *E*_{x}/ and
*v*_{c} = *E*_{c}/ are the electron-electron exchange and correlation potentials respectively,

*E*_{x} and *E*_{c} being the corresponding energies.

*v*_{Tcorrec} = *T*_{correc}/ is the potential corresponding to the kinetic energy
correction term
*T*_{correc}.

Finally *v*_{ext} is the potential corresponding to the externally applied field.

The functional form of the different components are as follows:

Z being the atomic number.

For He the exact form of the exchange energy and potential have been employed, viz.,

For the other atoms in question we have used parametrized local exchange energy functionals [10,12,20]. A local parametrized Wigner-type functional [10,12,21] has been used for the correlation energy for all the atoms.

The Kinetic energy correction term is taken as a modified Thomas-Fermi energy term [10,12,19]. The corresponding potential is found therefrom.

It may be noted that evaluation of the inter-electronic coulomb repulsion potential in cylindrical coordinates
and *z* presents its own challenges, especially for atoms with higher Z values. The Green's function in equation (9) has been expressed in terms of the half-integer degree Legendre function of the second kind
(*Q*_{m-}) [22]:

where

Equation (14) can be equivalently expressed as

where is the Neumann factor [23].

Substituting equation (16) in equation (9) and noting that m=0 if we consider an axis symmetric system we obtain

*Q*_{-1/2}() is expressed in terms of the complete elliptic integral of the first kind *K*())[24]

where

We have used the above procedure to calculate the value of the classical coulomb potential (*v*_{el-el}) at the grid boundaries. For the interior points of the cylindrical mesh *v*_{el-el} values are obtained from the solution of the Poisson equation

with the values at the boundaries serving as Dirichlet boundary conditions.

**3. Numerical solution of the diffusion equation in scaled cylindrical coordinates**

In order to solve equation (6) in cylindrical coordinates and *z* we use scaled cylindrical coordinates [25,26] and where and .

It may be noted that the equal spacings in the scaled coordinates ensure more number of and *z* points near the origin. This helps in tackling steep coulomb potentials near the origin, especially for systems with high atomic number. The diffusion function *R* must satisfy the normalization condition

where *N* is the total number of electrons. In order that the above normalization condition is satisfied in the scaled coordinate system we employ a transformed diffusion function
*R*^{}
where
. We also note that the transformed function
*R*^{} is zero at the origin thereby ensuring we do not encounter numerical difficulties at this Coulomb singular point.

We denote the non-linear operator in equation (6) as so that equation (6) looks as follows:

So . The Taylor expansion of
*R*^{}(**r**, *t* + *t*) about
*R*^{}(**r**, *t*) can then be expressed as

The evolution operator is real and non-unitary.

is splitted along the lines of Peaceman and Rachford [25] as follows:

where

Writing
*R*^{}(**r**, *t* + *t*) as
*R*^{n+1} and
*R*^{}(**r**, *t*) as
*R*^{n} we obtain a symmetric form of the above equation in terms of *Â* and :

Now using Peaceman Rachford (PR) Alternate Direction Implicit (ADI) scheme we recast the above equation into two equations involving the value of
*R*^{} at the fictitious intermediate time step
*t*/2 as follows:

and

We approximate the first and second order partial derivatives in equations (30) and (31) with two-point and three-point central difference formulae respectively and therefrom write each of the above two equations as a set of tridiagonal system of equations which are then solved with the knowledge of the boundary conditions [27] to yield
*R*^{n+1} from
*R*^{n}.
After each iteration normalization is forced on
*R*^{} and a fresh iteration is started.

The values of the parameters in the scaled coordinates have been taken as follows:

= 1.5 0.02 and = 14.

The number of points *N*1 and *N*2 in the and directions respectively have been chosen 401 and 401 for He, 801 and 801 for Ne and 1001 and
1001 for Ar, Kr and Xe. The grid boundaries were taken as 1.0e-06 and 0.6 and the grid boundaries were taken as -16.0 and +16.0 for all
the atoms. The temporal spacing *t* is chosen to have a value of 0.0005 for He, 0.00005 for Ne and 0.00001 for Ar, Kr and Xe. It may be mentioned
that the choice of a larger number of points and finer grids, especially for species with larger atomic numbers, yields better accuracy but at the cost of
computation speed.

There are various ways of solving the Poisson equation (20) subject to Dirichlet boundary conditions. We have however taken recourse to a fast and efficient program named SEPELI which is included in the FISHPACK [28] package of partial-differential-equation-solvers freely available at the `*www.netlib.org*' repository.

The solution has been launched by taking Slater type functions as tabulated in reference [10]. This ensures a faster rate of convergence compared to that using any other trial function.

**Results and discussion**

The most remarkable aspect of the results of our calculation is that all the quantities studied maintain perfect spherical symmetry, not only in their converged values at very long times, but also throughout the evolution process, in spite of the fact that the evolution is carried out in cylindrical coordinates and not in spherical coordinates. All the figures presented in this paper bear testimony to this aspect.

Figure 1.
shows the spatial profile of the radial densities of He, Ne, Ar, Kr and Xe as a function of the square root of the radial distance *r* calculated
at all the points of the scaled cylindrical grid after a sufficient number of iterations. The radial densities closely resemble
not only the Hartree Fock (HF) radial density profiles but also those of Deb and Ghosh[19] as also of Roy and Chu [12]. All of them
exhibit their inherent shell structures, just as expected. Figure 1. also shows how we may obtain completely spherically-symmetric densities
after evolving the GNLSE for a sufficiently long time, in spite of our calculation being done with the help of finite-differencing in (scaled)
cylindrical coordinates. To get a better insight into this aspect we present in Figure 2.
the radial density plots for *Xe* at intermediate steps of
the numerical evolution. Here again we observe how the spherical symmetry of the electron density is maintained throughout the evolution process.
It is important to mention that this aspect of maintainance of spherical symmetry may be thought of as a check of the correctness of the algorithm as also
the correct choice of the grid-parameters , the grid boundaries and the spatial and temporal grid spacings, , and *t* respectively. It has also been observed that the correctness of the algorithm is extremely sensitive to the correlation
between and with *t*. Although, no analytical criterion for stability has been made use of [25], we have
found from a lot of experimentation that *t* when chosen to be less than , yields sufficiently accurate results.

We have also studied the variation of the different potential terms constituting *v*_{eff} viz., *v*_{nu-el}, *v*_{el-el}, *v*_{x}, *v*_{c} and
*v*_{Tcorrec} with the radial distance from the nuclear site. The nuclear-electron attraction potential *v*_{nu-el} has been found to shoot down
remarkably fast as one approaches the Coulomb-singularity at the origin. The use of appropriately scaled coordinates does ensure that this potential term
is properly calculated at points close to the origin. Same holds true for the electron-electron repulsion potential *v*_{el-el} as well as the exchange and
correlation potentials *v*_{x} and *v*_{c} respectively. We have presented in Figure 3.
the plot of *v*_{el-el} with *r* for Xe. It is very gratifying to
note that the plot too shows spherical symmetry, just like those for the radial electron density, in spite of having been calculated using finite
differences along scaled cylindrical coordinates. The other potential terms also show spherical symmetry, which is however to be expected as they are
calculated as functions of the radial distance *r*, and hence their plots are not presented. And finally
*v*_{Tcorrec} reveals a marked shell structure
in its variation with *r* [19], just as expected. It may be recalled that this potential term is responsible for the shell structure in the
radial density and proper choice of the corresponding kinetic energy correction term (
*T*_{correc}) is very much essential in order to obtain the
correct electron density [19].

In Table 1. we have presented the nonrelativistic ground state energies and different energy components as obtained after evolving the GNLSE (6) for sufficiently large number of iterations in scaled cylindrical coordinates. The total energy in all cases does go beyond the corresponding Hartree-Fock values and are also very close to those found elsewhere [10,12]. It is important to keep in mind that our endeavour is not to find very accurate energy and ground state density of spherically-symmetric systems through our cylindrical-coordinates-based algorithm. Such accurate calculations can and have been performed by solving equations in only one (radial) coordinate [10,12]. Instead our aim has been to put to test our algorithm on the spherically-symmetric systems in question and see if they match the results obtained through the radial-coordinate-based algorithms [10,12]. Even if the results obtained are close, we would have gained confidence in using our algorithm on such systems where spherical symmetry breaks down and cylindrical symmetry holds, as for example, noble gas atoms subjected to intense laser fields along the axial directions. As is evident from the the figures and the table, we have been more or less successful to this end.

The values of the individual energy terms, though close to those from literature [10] for He and Ar appear to diverge somewhat as we move towards species with higher atomic number. We have however observed that the discrepancies can be reduced if we take larger and larger number of spatial grid points near the origin and also use smaller spatial grid spacings which in turn forces us to take yet smaller temporal spacing leading to remarkable decrease in the speed of our numerical calculation. What we have made use of is a trade-off between accuracy and the speed. We have settled for optimally dense grids which yield exactly spherically-symmetric densities with the expected shell structure and closely matching, if not very accurate energy values. Another point to be noted is that the rate of convergence becomes slower and slower as we approach the minimum energy ground state configuration. This constraint, though not an issue with He and Ne, does pose a problem with higher atomic-number species like Kr and Xe where we need to use much denser grids. Also to be noted is that the value of the Weizsäcker kinetic energy is higher in all the species in question. Consequently the value of the virial constant becomes less than the ideal value of 2.0. This is also somewhat to be expected as this term is calculated using finite differences in two-dimensions. It value can however be improved by following the prescription mentioned above, viz., by taking more points near the origin and smaller spatial grid spacings.

And finally a note on the correction term to the kinetic energy (
*T*_{correc}). The choice of the correction term is very crucial in obtaining the correct
density. Although the somewhat arbitrary choice of
*T*_{correc} in our work has yielded pretty good results, there is scope for finding improved and more
general forms of
*T*_{correc} so that we may make use of our algorithm to deal with larger systems of interest.

**Conclusions**

An imaginary-time evolution method has been made use of to obtain non-relativistic ground state electron-densities and energies for all the closed shell noble gas atoms through the solution of a single diffusion-type equation in appropriately scaled cylindrical coordinates. The radial plots of the radial- densities and associated energy terms are absolutely spherically-symmetric and reveal all the expected features. The calculated values of the non-relativistic energy terms are in close, if not exact agreement with those from the literature. The algorithm can be easily modified to study the dynamics in real-time of similar systems subjected to external polarized fields in terms of the electron density in 3-dimensional space. This would entail solving a single equation under the influence of a single effective potential as opposed to solving a number of equations as required in other methods. This becomes all the more important for atoms with higher atomic number and a greater number of electrons like, say Xenon, for which the number of equations needed to be solved would have been prohibitively large otherwise. The calculations being performed in (scaled) cylindrical coordinates, the 3-dimensional picture of the dynamics can also be visualized in real-time throughout the evolution process, thereby facilitating a mo

Figure 1. Radial density (
4*r*^{2}) versus where
*r* = for each point
of the scaled cylindrical coordinate mesh (clockwise from bottom) for He, Ne, Ar, Kr and Xe respectively at convergence.

Figure 2. Radial density (
4*r*^{2}) versus for Xe where
*r* = for each point of the scaled cylindrical coordinate mesh (clockwise from bottom-left) at the start, after 300 time steps, after 600 time steps, after 1100 time steps, after 7000 time
steps and at convergence respectively.

Table 1. Values of ground state properties (a.u) for different noble gas atoms

Property | Ref. | He | Ne | Ar | Kr | Xe |

- E |
PW | 2.9068 | 128.7422 | 528.3923 | 2754.1551 | 7234.8235 |

HF | 2.8617 | 128.5470 | 526.8174 | 2752.0546 | 7232.1302 | |

Ref [10] | 2.8973 | 128.9065 | 527.5486 | 2753.8809 | 7234.9742 | |

- E_{nu-el} |
PW | 6.7969 | 307.6406 | 1269.4524 | 6586.6530 | 17199.8669 |

HF | 6.7492 | 311.1333 | 1255.0504 | 6582.5412 | 17164.9821 | |

Ref [10] | 6.7850 | 311.0597 | 1245.5699 | 6533.8352 | 17038.2385 | |

E_{el-el} |
PW | 2.0675 | 61.7228 | 234.3324 | 1100.7495 | 2616.4803 |

HF | 2.0516 | 66.1476 | 231.6093 | 1172.3372 | 2880.0352 | |

Ref [10] | 2.0651 | 65.7129 | 220.6552 | 1119.3762 | 2744.7642 | |

- E_{x} |
PW | 1.0334 | 11.7012 | 30.5016 | 91.4425 | 172.1502 |

HF | 1.026 | 12.11 | 30.19 | 93.89 | 179.2 | |

Ref Ref [10] | 1.0325 | 12.1111 | 29.4850 | 91.5847 | 173.9435 | |

- E_{c} |
PW | 0.0423 | 0.3391 | 0.7371 | 1.7259 | 2.7123 |

HF | 0.042 | 0.390 | 0.787 | 1.835 | 2.921 | |

Ref [10] | 0.0423 | 0.3561 | 0.7011 | 1.7529 | 2.8407 | |

T_{W} |
PW | - | 95.1316 | 328.9982 | 1433.7140 | 3445.8170 |

HF | - | 90.1640 | 308.4206 | 1276.7349 | 2932.0548 | |

Ref [10] | - | 94.2068 | 322.0345 | 1377.5940 | 3226.9174 | |

T_{correc} |
PW | - | 34.0841 | 208.9682 | 1391.2027 | 4077.6085 |

HF | - | 37.3886 | 214.4033 | 1465.2484 | 4298.9068 | |

Ref [10] | - | 37.7006 | 205.5177 | 1376.3217 | 4008.3670 | |

Virial |
PW | 2.002 | 1.996 | 1.9822 | 1.975 | 1.962 |

HF | 2.000 | 2.000 | 2.000 | 2.000 | 2.000 | |

Ref [10] | 1.999 | 1.999 | 1.999 | 1.999 | 1.999 | |

*PW* stands for present work.
*HF* stands for values obtained using Hartree-Fock density [10].

2007-09-27