A two-dimensional pseudospectral Hartree-Fock method for low-Z atoms in intense magnetic fields

A two-dimensional pseudospectral Hartree-Fock method for low-Z atoms in intense magnetic fields

Anand Thirumalai anand.thirumalai@asu,edu. School of Earth and Space Exploration, Arizona State University, Tempe, Arizona, USA, 85287    Jeremy S. Heyl heyl@phas.ubc.ca. University of British Columbia, Vancouver, British Columbia, V6T 1Z1
July 6, 2019
Abstract

The energy levels of the first few low-lying states of helium and lithium atoms in intense magnetic fields up to  T are calculated in this study. A pseudospectral method is employed for the computational procedure. The methodology involves computing the eigenvalues and eigenvectors of the generalized two-dimensional Hartree-Fock partial differential equations for these two- and three-electron systems in a self-consistent manner. The method exploits the natural symmetries of the problem without assumptions of any basis functions for expressing the wave functions of the electrons or the commonly employed adiabatic approximation. It is seen that the results obtained here for a few of the most tightly bound states of each of the atoms, helium and lithium, are in good agreement with findings elsewhere. In this regard, we report new data for two new states of lithium that have not been studied thus far in the literature. It is also seen that the pseudospectral method employed here is considerably more economical, from a computational point of view, than previously employed methods such as a finite-element based approach. The key enabling advantage of the method described here is the short computational times which are on the order of seconds for obtaining accurate results for heliumlike systems.

pacs:

I Introduction

The motivation to study atoms in magnetic fields of strength beyond the perturbative regime was in a large part due to the discovery of such fields being present in white dwarf stars Kemp et al. (1970); Angel (1978); Angel et al. (1981) and neutron stars Trümper et al. (1977); Truemper et al. (1978). The most commonly observed neutron stars - pulsars, have been observed to have magnetic fields on the order of - T Ruder et al. (1994). Magnetars Duncan and Thompson (1992), which are strongly magnetized neutron stars, can have magnetic field strengths well in excess of T. White dwarf stars on the other hand have somewhat less extreme fields, albeit still high, - T Ruder et al. (1994). At such high field strengths, a Zeeman-type perturbative treatment of the field Landau and Lifshitz (2003) is not possible. The structure of atoms is considerably altered from the low field case.

The problem of atoms in magnetic fields has been tackled by various researchers since the 1970’s using a variety of different methods. In the literature, there exist numerous studies of hydrogen Canuto and Kelly (1972); Praddaude (1972); Simola and Virtamo (1978); Friedrich (1982); Wunner and Ruder (1980); Roesner et al. (1983, 1984); Ivanov (1988); Baye et al. (2008, 2002) and many recent studies of helium Proeschel et al. (1982); Thurner et al. (1993); Ivanov (1994); Jones et al. (1996a, 1997, 1999); Heyl and Hernquist (1998); Mori and Hailey (2002); Mori and Ho (2007); Al-Hujaj and Schmelcher (2003, 2003); Becken and Schmelcher (2002, 2001, 2000); Becken et al. (1999); Luhr et al. (2007); Ivanov and Schmelcher (1998, 2000); Wang and Qiao (2008) atoms in strong magnetic fields. There have also been studies conducted for molecules and chains of atoms for both hydrogen and helium atoms in strong to intense magnetic fields Lai et al. (1992, 1993); Lai and Salpeter (1996); Medin and Lai (2006a, b); Detmer et al. (1998); Schmelcher et al. (2000, 2001). Moreover, our recent investigation Thirumalai and Heyl (2009) using single-configuration Hartree-Fock (HF) theory Hartree (1957) was seen to yield accurate upper bounds for the binding energies of hydrogen and helium in strong magnetic fields. Our later study Heyl and Thirumalai (2010), obtained accurate binding energies for helium and lithium atoms in strong magnetic fields using a pseudospectral method. This approach was seen to be computationally far more economical than using our earlier finite-element based approach Thirumalai and Heyl (2009).

In sharp contrast to the somewhat simpler two-electron systems, there is very limited work available in the literature for atoms with more than two electrons. One of the first studies to investigate atoms in intense magnetic fields, in particular the iron atom, was by Flowers et al Flowers et al. (1977) in 1977. This variational study extended the work due to the authors in Ref. Glasser and Kaplan (1975) and obtained binding energies of iron atoms and condensed matter in magnetic fields relevant to neutron stars. Errors in this study were later corrected by Muller Mueller (1984). Other methods included density functional studies Jones (1985, 1986) and also employed the Thomas-Fermi-Dirac method Mueller et al. (1971); Skjervold and Ostgaard (1984) for estimating binding energies of atoms in intense magnetic fields. The first comprehensive HF studies of atoms with more than two electrons were carried out by Neuhauser et al Neuhauser et al. (1986, 1987) for magnetic fields greater than T, thus being directly relevant to neutron stars. Elsewhere, HF studies of atoms and molecules in intense magnetic fields were conducted by Demuer et al Demeur et al. (1994), with results consistent with previous findings. All of the above treatises, Refs. Flowers et al. (1977); Glasser and Kaplan (1975); Mueller (1984); Jones (1985, 1986); Mueller et al. (1971); Skjervold and Ostgaard (1984); Neuhauser et al. (1986, 1987); Demeur et al. (1994), concern themselves with magnetic fields in excess of T, well into the so-called intense magnetic field regime. At these field strengths, the interaction of the electron with the nucleus of the atom becomes progressively less dominant, in comparison to its interaction with the field itself. One of the first studies to carry out a rigorous HF treatment of atoms with more than two electrons in strong or intermediate field strengths was Ref. Jones et al. (1996a). Therein, they obtained estimates of the binding energies of a few low-lying states of lithium and carbon atoms, in low to strong magnetic fields. Elsewhere, Ivanov Ivanov (1998) and Ivanov & Schmelcher Ivanov and Schmelcher (1998); Schmelcher et al. (1999); Ivanov and Schmelcher (1999, 2000, 2001a, 2001b) have over recent years, carried out detailed HF and post-HF studies of atoms with more than two electrons using a numerical mesh-method for solving the unrestricted HF equations Ivanov and Schmelcher (1998). The special meshes were so constructed as to facilitate finite-difference calculations in a two-dimensional domain using carefully selected mesh node points Ivanov and Schmelcher (2001c). They were able to ascertain the binding energies of the first few low-lying states of low-Z atoms such as lithium, beryllium and mid-Z atoms such as boron and carbon etc., using this method. Moreover, using a gaussian basis of functions for expressing the wave functions of the electrons Al-Hujaj and Schmelcher (2003, 2003); Becken and Schmelcher (2002, 2001, 2000); Becken et al. (1999); Luhr et al. (2007), adopting a full configuration-interaction method, Al-Hujaj & Schmelcher Al-Hujaj and Schmelcher (2004a, b) have been able to estimate the binding energies of lithium and beryllium atoms in strong or intermediate magnetic fields, thereby improving upon previously obtained results. The sodium atom in a strong magnetic field has also been studied by Gonzalez-Ferez & Schmelcher González-Férez and Schmelcher (2003) obtaining estimates for the binding energies. Elsewhere, low lying states of the lithium atom have also been studied in strong magnetic fields using a configuration-interaction method, employing the so-called freezing full-core method both with Qiao and Li (2000) and without Guan and Li (2001) correlation between electrons. Recently, Medin & Lai Medin and Lai (2006a, b) have also studied atoms and molecules and infinite chains of condensed matter in magnetic fields greater than T, using density-functional-theory. Mori et al Mori and Hailey (2002); Mori and Ho (2007) have studied mid-Z atoms in strong to intense magnetic fields using perturbation theory as well, obtaining results consistent with previous findings. In recent years Engel and Wunner and co-workers Engel and Wunner (2008); Engel et al. (2009); Schimeczek et al. (2013, 2012); Klews et al. (2005) have computed accurate results for several atoms in magnetic fields relevant to neutron stars with a variety of techniques involving finite-element methods with B-splines both in the adiabatic approximation and beyond the adiabatic approximation with more than one Landau level. These highly accurate formulations employ a fast parallel Hartree-Fock-Roothan code, in which the electronic wave functions are solved for along the direction, with Landau orbitals (and combinations of more than one level in the latter studies) describing the remaining parts of the wave functions. Elsewhere, different ab initio Quantum Monte-Carlo approaches Bücheler et al. (2007, 2008) have also been successfully employed for determining the ground states of atoms and ions in strong magnetic fields including lithium. Recently excited states of helium have also been computed quite accurately in intense magnetic fields using a fixed-phase Quantum Monte-Carlo approach Meyer et al. (2013). The recent results of Schimeczek et al. (2012, 2013) for obtaining the ground state energies of atoms up to , were obtained within only a few seconds of computing time for helium and helium like atoms. Such speeds are essential for coupling atomic structure codes with atmosphere models and spectral analysis codes for magnetized white dwarfs and neutron stars. However, the work of Schimeczek et al. (2013) only concerns itself with the ground state configurations. The primary aim of the current study is to provide an accurate method for investigating the ground and excited states of atoms, which can be computed in a matter of seconds on run-of-the-mill computer architectures. The method described here provides a speed-up of a factor of when compared with our earlier investigation using a finite-element discretization Thirumalai and Heyl (2009). Computational times are comparable to the recently developed quantum monte-carlo methods due to Schimeczek et al. (2013). The method described in this study can be run on a distributed architecture for investigating multi-electron atoms.

It is also well known that post-HF methods such as configuration interaction (CI) or multiconfigration Hartree-Fock (MCHF) methods, yield considerable improvements with regard to the estimates of the upper bounds for the energies of various states. In the intermediate range of magnetic field strengths, where both the nucleus of the atom and the magnetic field have interactions with the electrons that are approximately equal in magnitude, the single configuration approximation then becomes increasingly ineffective with greater number of electrons. However, these methods are computationally more intensive than a single configuration calculation. Thus far, the most accurate CI methods involve decomposing the wave functions into a Gaussian basis set relying upon separation of variables in cylindrical coordinates (e.g. Al-Hujaj and Schmelcher, 2004a, b). These methods do however require a large set of basis functions. On the other hand, MCHF methods would require fewer basis functions, as the orbitals get optimized during the computation with the coefficients (e.g. Fischer, 1997). Separation of variables and/or basis decompositions speed up the computation in these post-HF methods considerably. However, there do not exist hitherto, any fully two-dimensional (2D) post-HF studies of multi-electron atoms in intense magnetic fields. This is partly due to the computational overhead associated with adopting a fully 2D picture. Central to the development of such a method would be the fast and accurate computation of the single-configuration problem in a full 2D framework without any basis expansions and separation of variables. Wave functions so determined could be used directly in 2D configuration-interaction calculations or the problem could be cast into a MCHF framework. Moreover, obtaining accurate estimates of the energy levels of atoms, in particular low-Z atoms, in strong and intense magnetic fields will ultimately facilitate a proper understanding of the spectra of neutron stars and white dwarf stars. In order to achieve this, a large amount of accurate data needs to be made available for not only binding energies of several different states, but also oscillator strengths, as well as estimates of the energies associated with both bound-free and free-free electron transitions in a variety of magnetic field strengths. There is also the added complication of strong electric fields affecting the spectra of atoms moving perpendicular to the magnetic field; however that being said, estimating the binding energies due to intense magnetic fields is part of that picture. Thus the central aim of the current work is to take a step in that direction and provide accurate and readily calculable estimates of the binding energies of the first few low-lying states of the simplest low-Z atoms; helium and lithium, in strong and intense magnetic fields using a two dimensional single-configuration pseudospectral method of solution.

The method outlined in the current study is an extension of the method developed in our two previous studies Refs. Thirumalai and Heyl (2009); Heyl and Thirumalai (2010).

Ii The HF Equations

We shall begin with the generalised single-configuration HF equations for an atom with -electrons and nuclear charge , in a magnetic field that is oriented along the -direction. A derivation of the single-configuration HF equations can be found in our earlier work Thirumalai and Heyl (2009), here we shall present only the salient points. The single configuration HF equation can be written in cylindrical coordinates as, where the length scale is in units of Bohr radii and the energy is scaled in units of Rydberg energy in the Coulomb potential of charge (see below for definitions).

(1)

where and Please note that the three-dimensional momentum operator has been split into two parts; which are the and parts of the Laplacian and which is the azimuthal part. The total Hartree-Fock energy of the state is given by

(2)

The direct () and exchange () interactions are determined according to the method outlined in Ref. Thirumalai and Heyl (2009), as the solutions of the elliptic partial differential equations for the potentials given by,

(3)

and

(4)

where and are the wave functions of the and electrons. The wave function of a given configuration of electronic orbitals is assumed to be given by a single Slater determinant as,

(5)

where is the anti-symmetrization operator. The individual electronic wave functions are given by,

(6)

where labels each of the electrons. The two-dimensional single particle wave functions are taken to be real functions.

Integration with respect to the azimuthal coordinate, , has been carried out, prior to writing the result in Eq. (1) above. The contribution due to electron spin has also been averaged out a priori. It is to be mentioned in this regard, that in the current study we shall only be concerned with fully spin-polarised states (FSP), in other words all the electrons of the atom are assumed to be anti-aligned with the magnetic field. Such states have an exchange interaction between the electrons providing an extra coupling term in the HF equations, . Additionally, FSP states are seen to be the most tightly bound states in the intense field regime. The extension to partially spin-polarised configurations is easily achieved by eliminating the exchange term in the HF equations. In the current study, we have chosen to work in units of Bohr radii along with the definitions given below.

The Bohr radius for an atom of nuclear charge is given by , where is the Bohr radius of the hydrogen atom. The magnetic field strength parameter , is given by the expression , where is the critical field strength at which point the transition to the intense magnetic field regime occurs Ruder et al. (1994). This is defined as T. Thus, beyond a value of , the interaction of the electron with the nucleus becomes progressively less dominant, as increases. Based upon the above definition of , it is convenient to classify the field strength Jones et al. (1996b) as low (), intermediate, also called strong () and intense or high (). These definitions of the different magnetic field strength regimes are useful to remember when discussing the results in the latter part of this paper and for distinguishing between “strong” and “intense” magnetic field strengths. The energy parameter of the electron is defined as , with , the Rydberg energy of the hydrogen atom. For brevity we shall refer to the units of energy as , which should be remembered as the Rydberg energy in the Coulomb potential of charge . The quantity is the fine structure constant. In the current study, all the physical constants were used in SI units. Additionally, the magnetic field , is taken to be in units of Tesla. Eq. (1) represents the -coupled Hartree-Fock equations in partial differential form for an -electron system with nuclear charge . The equations are coupled through the exchange interaction term between the electrons and as such the system of equations is solved iteratively. The discussion that follows is arranged in the following manner. In Section III, we shall describe the numerical methodology employed in the current study. For solving the system of partial differential equations we adopted a pseudospectral approach. The interested reader will find the entire method of discretizing and setting up the discrete problem using this approach in the appendix. Readers familiar with the method may wish to only read Sections III-V. In Section IV the results of the current study are provided alongside a discussion. Section V contains the the conclusions of the current study and avenues for further work are also discussed therein.

Iii Numerical Details

The numerical solution of the coupled eigenvalue problem in Eq. (1) proceeds via the so-called self consistent field (SCF) method due to Hartree Hartree (1957). First we find a solution to the hydrogenic problem, i.e., Eq. (1) without the direct and exchange interactions. This yields ionic single electron hydrogenic wave functions in the Coulomb potential of charge forming the initial estimates for the HF iterations. Second, using these estimates the elliptic partial differential equations for the direct and exchange interaction potentials in Eqs. (3) and (4) are solved. With these potentials now obtained, the coupled HF problem including the direct and exchange interactions in Eq. (1), is solved as an eigensystem. The exchange interactions which couple the equations are expressed using wave-functions from the previous iteration to solve the eigenvalue problem for each electron Slater (1951). The eigenvalues obtained are the individual particle energies and the normalized eigenvectors are the wave functions, . The SCF iterations then proceed with the updated electron wave functions and the steps from the second step described above, are repeated until convergence.

iii.1 Domain Discretization

The very first step in this direction is the discretization of the physical domain of the problem. By virtue of azimuthal symmetry and parity with respect to the plane, it is sufficient to restrict the domain of the problem to Ruder et al. (1994); Thirumalai and Heyl (2009). However, we solve the problem not in this semi-infinite domain, but rather in a finite but sufficiently large domain of size . Whenever the domain of a problem is restricted, this introduces an error since the computational problem then becomes an approximation of the actual problem. Therefore in order to eliminate this error associated with truncating the domain, we employ a sequence of domains of increasing sizes, obtaining a converged result in the limit of the computational domain approaching the size of the physical domain of the problem.

In our computations, the size of the computational domain and (in units of Bohr radii) are given by,

(7)

where is a scaling factor used for setting up computations in a sequence of increasing domain sizes. The effect of the logarithmic term , in the denominator is that it naturally makes the domain larger or smaller, depending on whether or , respectively. With the maximum domain size thus defined, we can then compactify the finite domain to with the transformation,

(8)

and

(9)

where and . Note that in our calculations we employed a square domain for achieving the best possible internally consistent convergence. Therefore in our work and therefore , but the possibility remains for using different sizes and scalings in the two orthogonal directions for optimizing computational effort particularly in the intense field regime.

With this compactification scheme, it is then possible to employ a Chebyshev-Lobatto spectral collocation method Trefethen (2000) with discrete points on the domain . The details regarding the pseudospectral approach developed for solving the HF equations in Eq. (1) can be found in the appendix. An atomic structure software package based on this method was written in the high level programming language MATLAB making particular use of its fast matrix manipulation algorithms. The eigenvalue problems described in Eqs. (36) and (44) are solved by discretizing the equations and solving the resultant algebraic eigenvalue problem. This ultimately produces a sparse matrix for the coupled eigenvalue problem (see the appendix and Fig. (8) therein). Thus we can take advantage of this fact and employ a sparse matrix generalized eigensystem solver; the widely used package ARPACK which utilizes the implicitly restarted Arnoldi method (IRAM) Arnoldi (1951); Sorensen (1992); Lehoucq et al. (1998). The key advantage of employing IRAM is that the memory storage requirements are considerably reduced since the Arnoldi factorization generates a small Krylov subspace using a few basis vectors and only a handful of eigenvalues are computed in a given portion of the spectrum Lehoucq et al. (1998); Saad (1984); Sorensen (1992). Since the Hamiltonian matrix that we are solving only has a few bound state solutions, employing IRAM for computing only a portion of the spectrum is therefore highly desirable and saves considerable computational effort, particularly for a two-dimensional problem where the discretisation can generate large matrices. This method was found to yield accurate results for the energy eigenvalues of the first few eigenstates of helium and lithium in intense magnetic fields. It was seen that generating a Krylov subspace with about to basis vectors was sufficient for determining around to eigenvalues in the vicinity of a given shift (), by employing the shift-invert algorithm Lehoucq et al. (1998). Runs were carried out for different values of the magnetic field strength parameter , in the range , for the cylindrical pseudospectral code. A typical tolerance of around was employed for the internal errors of ARPACK. It was observed during our runs that fast convergence was achieved; within about HF iterations. A convergence criterion for the HF iterations was employed wherein the difference between the HF energies for two consecutive iterations was tested. Typically, a tolerance on the order of was employed. Once the HF iterations attained convergence for a given level of mesh refinement, the total energy of the Hartree-Fock state under consideration is reported according to Eq. 2. Additionally for testing convergence of the pseudospectral method, we employed up to six different levels of mesh refinement for lithium and seven for helium, ranging from coarse to fine mesh i.e., mesh points in each direction, for a given domain size, . Once a converged solution was obtained for a given domain size, using a sequence of meshes, a different value for the scaling parameter was employed and the calculation was carried out again.

iii.2 Computational speeds

One of the key enabling advantages of the method adopted in this study is the reduced computational time. The gain in speed is largely due to two reasons. First, since the eigenvalue problem represents a smooth elliptic one whose solution does not have discontinuities, a pseudospectral approach yields accurate results with the use of a small number of discretisation points or grid elements. Second, our code is implemented on a parallel architecture, where the number of processors (or cores) employed is equal to the number of electrons in the atom. With this trivial parallelization, the eigenvalue computation then proceeds in parallel for all the electrons, and the computational time is then determined by the electron for which the eigenvalue computation takes the longest time. The calculations for lithium for example took a total of about seconds of computing time running on three Intel Xeon E5620  GHz processors, for obtaining accuracy . The fastest calculations we observed were for helium where the calculations converged on the order of about seconds, with computations carried out on two processors, for achieving the same level of accuracy. These speeds are comparable to the recent results of Schimeczek et al. (2013) who in their calculations achieve considerable speed up by first pre-calculating the inter-electron potentials and second by utilizing Landau states to describe the wave functions in the directions, thus only having to solve for the unknown -part of the wave functions. In our approach, we are not restricted to the adiabatic approximation and the unknown wave function that is solved for is a two-dimensional one, and the eigenvalue problem is therefore a fully coupled two-dimensional problem. Additionally, the inter-electron potentials in our approach are recalculated as the wave function gets optimized in each HF iteration, by solving elliptic partial differential equations, which adds to the computational times. This combination of speed and accuracy achieved in our computations is a result of spectral convergence. Figure 1 shows the dependence of computational time on both the level of domain discretization and the number of electrons in the problem, which is equal to the number of processors employed. The computational time increases roughly as , where is the number of grid points in each direction. The lines drawn through the data in Fig. 1 correspond to polynomials of degree .

Figure 1: Computational times shown as a function for the number of grid points, for typical calculations for both helium and lithium. The dotted vertical line corresponds to the number of grid points typically required to attain an accuracy of . The number of cores or processors employed by the parallel code is equal to the number of electrons in the atom. The dependence of the computational time on the number of grid points is roughly , while the dependence on the number of electrons is more or less linear, i.e., .

It can also be seen that the compute times increase with the number of electrons (). However, it can be seen that this dependence scales more or less linearly, i.e., . As a result, we expect that the method can be extended to mid-Z atoms with , without incurring a large computational overhead.

iii.3 Limitations of the approach

The atomic structure package developed in this study also has certain limitations. First, since the method employed presently requires a finite domain, computations are required to be carried out on a sequence of domains to eliminate errors due to domain truncation. This could be circumvented by employing an adaptive scheme, wherein the wave functions at the exterior domain points () are required to fall below a certain threshold or tolerance (e.g. Schimeczek et al., 2013) while employing an exponential grid. This would however complicate the pseudospectral approach, although a key feature of using Chebyshev-Lobatto points is that the density of points is greatest at the end points of the domain, which may aid in resolving the wave function at the outer limits of the domain, should such a scheme be implemented. Second, the current work does not include relativistic corrections to the energies. For the magnetic field strengths considered herein, the relativistic corrections to the energies were estimated using the scaling formula in Ref Poszwa and Rutkowski (2004). Their results for the hydrogen atom were used for this purpose and the corrections were estimated to be on the order of . This was seen to be smaller than the numerical errors arising from convergence of the entire numerical method including the extrapolation to the limit of a semi-infinite domain. Thus, while relativistic corrections are important, it was not possible to account for them accurately in the current study. Presently, the results are presented alongside a discussion in the following section.

Iv Results and Discussion

The atomic structure software package developed in the current study extends our earlier computations for atoms in strong magnetic fields Thirumalai and Heyl (2009); Heyl and Thirumalai (2010), towards the intense field regime, . The states that were considered in this study are labelled using both the field-free and strong-field notations for the sake of clarity; these can be found in Table 1 which lists these different states of helium and lithium. In the presence of a magnetic field states can be characterized using the notation , where is the total component of angular momentum. The summation is over all the electrons in the atom. This then forms a manifold within which different sub-spaces exist. The quantum number counts the excitation level within a given manifold and sub-space symmetry. The spin multiplicity is given in the usual way as . Finally, the parity of the state is indicated using , indicating positive or negative parity. We studied the six most tightly bound states of each of these two atoms in the intense magnetic field regime (). Within a given parity sub-space, typically there are crossovers that occur as the magnetic field is reduced and the reader is referred to Ref. Ivanov and Schmelcher (2000) for an excellent discussion regarding ground state crossovers. Within each given parity sub-space, we considered the most tightly bound state in the intense field regime.

Intense-field   Field-free
Helium
Lithium
Table 1: The different states of helium and lithium considered in this study, listed using both intense-field and field-free notation.

iv.1 The Helium Atom

For the states of helium listed in Table 1, eigenvalues were determined using the numerical method described in Section III. We began with the lowest value of the domain scaling parameter . This yielded a domain with dimensions given according to Eq. (7), and this domain size depends on . HF energies were then calculated using up to six different levels of mesh refinement in the domain. This enabled us to extrapolate the results to the limit of infinitely fine mesh, for a given domain size. Thereafter, the domain was rescaled to larger and larger values, corresponding to and the computations repeated once again with different levels of mesh refinement for each value of so chosen. Then, using the extrapolated values of the HF energy corresponding to infinitely fine mesh for each of the four domain sizes, a subsequent extrapolated value of the the HF energy () was obtained, in the limit of the domain size approaching infinity. These are then the “converged ” values reported in Figs. 2 and  3 as well as in Table 2. These figures show typical examples of spectral convergence that was achieved using our pseudospectral implementation, wherein the internal errors of the procedure diminish in an exponential fashion with mesh refinement.

Figure 2 shows HF energies obtained from calculations for the state of helium, in a magnetic field of strength . Meanwhile Fig. 3 shows convergence data for the state of helium in a magnetic field corresponding to . The spectral method converged to better than typically with mesh refinement within each domain size. At the higher end of the intense field regime, there was some loss of accuracy in the convergence. In this region, the spectral method only converged to with mesh refinement, particularly for the largest domain scaling parameter employed; . These HF energies for different mesh refinements were extrapolated to the limit of infinitely fine mesh using an exponential function which typically took the form . The errors associated with the extrapolation procedure were typically on the order of with a normalized squared value typically for the interpolating function employed. Again, at the upper end of the intense magnetic field regime, we noticed some loss of accuracy as the states become tightly bound, and for the extrapolation procedure had an error on the order of with a normalized squared of on average. For the extrapolation to infinitely fine mesh, the average area per unit grid size in the domain (), was taken as the independent variable and the energies extrapolated to the limit of , corresponding to infinitely fine mesh.

It can also be seen in Figs. 2 and 3 that by increasing the domain size (by increasing the domain scaling parameter ), the binding energies decrease. This clearly shows that a confinement energy is introduced when the domain of the problem is truncated which can lead to an overestimation of the binding energy. Additionally, in Figs. 2 and 3, it is readily seen that when the area of the domain is doubled, the difference in the binding energies are roughly halved, with each doubling. For example, , illustrating quadratic convergence of the method with increasing domain size. Graphically this can be seen in the approximate halving of the distance between the horizontal lines joining the data points, when the domain area is doubled. This convergence can be expressed as,

(10)

where . This can be readily seen in Figs. 2 and 3, where the converged value is approximately twice as displaced from the data corresponding to in comparison to the data corresponding to . In order to numerically obtain a converged HF energy in the limit of the domain size approaching infinity, we employed an extrapolating function of the form . The ordinates in this case were the four different converged HF energies in the limit of infinitely fine mesh in each of the four different domains, and the abscissae were the inverse domain areas, i.e. . Thus, extrapolating to zero inverse area corresponding to an infinite domain size yields the final converged HF energy indicated in Figs. 2 and 3, and it is these extrapolated values that are reported in Tables 2 and 3. The error in the extrapolation to the limit of an infinite domain size was on the order of with a normalized squared value of for the interpolating quadratic function. This extrapolation to the limit of infinite domain size eliminates the effect of introducing a confinement energy when the domain of the problem is truncated.

Inspection of Figs. 2 and 3 reveals that the single-configuration converged HF energy in the limit of an infinite domain size is in good agreement with the corresponding results due to Ref. Becken et al. (1999) and Ref. Becken and Schmelcher (2001). The results obtained in this study therefore represent an upper bound for uncorrelated HF energies. These estimates can be improved using correlated post-HF methods, such as CI or MCHF methods, and the D wave functions computed by our method can form the initial inputs. The difference between our results and those of Ref. Becken et al. (1999) or Ref. Becken and Schmelcher (2001) therefore convey a sense of how large the effect of electron correlation can be even for helium. This effect typically increases with the addition of more electrons, becoming increasingly more pronounced, and therefore it becomes important that for accurate binding energies, electron correlation be taken into account. However, this is outside the scope of the current single-configuration calculation and we leave such a full D post-HF computation for a future undertaking.

Moreover, by comparing Figs. 2 and 3 it can be seen that convergence with regards to mesh refinement is much cleaner when the magnetic field strength is lower. In Fig. 2 it can be seen that regardless of the domain scaling parameter , a converged solution (to well within ) is obtained with as few as about points whereas, for larger , the number of points required is typically greater, , for obtaining a similarly converged result. This is due to the fact that the geometry of the wave functions becomes increasingly extreme with increasing magnetic field strength and a greater number of discretization points are therefore required to accurately resolve the wave functions. Overall, we observed that on average about discretisation points were enough for achieving convergence to well within .

Figure 2: Convergence of the binding energy with mesh refinement for four different domain sizes for the state of Helium. The level of mesh refinement corresponds to points in each of the two orthogonal directions. The size of matrix of the coupled eigenvalue problem for a given level of mesh refinement is given by , where , is the number of electrons, which for helium is two. The levels of mesh refinement employed correspond to and points in each of the and . The exponential convergence of the spectral method can readily be seen.
Figure 3: Convergence of the binding energy with mesh refinement for four different domain sizes for the state of Helium. The level of mesh refinement corresponds to points in each of the two orthogonal directions. The size of matrix of the coupled eigenvalue problem for a given level of mesh refinement is given by , where , is the number of electrons, which for helium is two. The levels of mesh refinement employed correspond to and points in each of the and directions for helium. The exponential convergence of the spectral method can readily be seen for computations in the different domain sizes.

The converged HF energies for the positive parity states of helium are given in Table 2, while those for the negative parity states are shown in Table 3, alongside data from the correlated configuration interaction calculations of Schmelcher et al Becken et al. (1999); Becken and Schmelcher (2000, 2001) as well as Quantum Monte Carlo results due to Jones et al. (1999) and the recent results due to Schimeczek et al. (2013). The absolute values of the binding energies are provided therein. In this study we investigated the three most tightly bound states within each parity sub-space. The corresponding weak field orbitals of these states are those listed in Table 1. The values given in parentheses are from eigenvalue computations using spherical coordinates. These calculations were carried out using our pseudospectral atomic structure software developed in an earlier study Heyl and Thirumalai (2010). An improved and faster version of the code was employed, again using the same levels of mesh refinement for maintaining consistency. The cylindrical pseudospectral method begins to lose accuracy as the magnetic field decreases, in the weak field region (), while in contrast, the spherical pseudospectral method loses accuracy in the upper end of the strong field regime, i.e., . Therefore using a combination of the two types of codes, we can explore the entire range in .

Refs. Refs. Refs. Refs. Refs. Refs. Refs. Refs. Refs.
Here Ruder et al. (1994)/Thirumalai and Heyl (2009) Jones et al. (1999) Becken et al. (1999) Here Ruder et al. (1994)/Thirumalai and Heyl (2009) Jones et al. (1999) Becken and Schmelcher (2001)/Schimeczek et al. (2013) Here Ruder et al. (1994)/Thirumalai and Heyl (2009) Jones et al. (1999) Becken and Schmelcher (2000)/Ivanov and Schmelcher (2000)/Schimeczek et al. (2013)
0 (1.0899) 1.0871d 1.0876 (1.0307) 1.0179a 1.0278 (1.0686) 1.0668a 1.0657c 1.0666
0.01 (1.1256) 1.1213 1.1213 1.1220 (1.0872) 1.0852a 1.0830 1.0833 (1.1223) 1.1183a 1.1177 1.1193
0.05 (1.2103) 1.1911 1.2056 1.2064 (1.2212) 1.2175a 1.2160 1.2167 (1.2734) 1.2691a 1.2683 1.2704
0.1 (1.2914) 1.2133 1.2860 1.2868 (1.3495) 1.3510a 1.3436 1.3450 (1.4209) 1.4189a 1.4151 1.4178
0.125 (1.3301) 1.3253 (1.4063) 1.4016 (1.4859) 1.4828
0.2 (1.4389) 1.4330 1.4338 (1.5572) 1.5598a 1.5525 (1.6579) 1.6585a 1.6508 1.6544
0.25 (1.5112) 1.4999 (1.6459) 1.6411 (1.7583) 1.7545
0.5 1.7722 1.4670 1.7718 2.0035 2.0009a 1.9945 (2.1582) 2.1550a 2.1490
0.625 1.8848 1.8841 2.1426 2.1425 (2.3183) 2.3128
0.7 1.9457 1.6690 1.9447 2.2213 2.2246a 2.2171 2.2210c (2.4060) 2.4029a 2.3956
1 2.1607 1.9116 2.1601 2.4983 2.4981a 2.4933 2.4978c 2.7047 2.7026a 2.7000 2.7043c
1.25 2.3141 2.3137 2.6943 2.6940 2.9198 2.9197
2 2.6863 2.4793 2.6840 3.1691 3.1685a 3.1634 3.1691c 3.4387 3.4384a 3.4333 3.4386c
2.5 2.8868 2.8862 3.4275 3.4274 3.7203 3.7203
5 3.6311 3.4672 3.6272 4.3785 4.3740a 4.3693 4.3783c 4.7539 4.7502a 4.7441 4.7523c
6.25 3.9082 3.9076 4.7383 4.7380 5.1425 5.1421
7 4.0618 4.9310 4.9268a 4.9203 4.9310c 5.3514 5.3474a 5.3408 5.3507c
10 4.5742 4.4362 4.5693 5.5856 5.5851a 5.5770 5.5899c 6.0548 6.0543a 6.0506 6.0624c
12.5 4.9219 4.9215 6.0441 6.0443 6.5527 6.5524
20 5.7538 5.6367 5.7473 7.0968 6.9867 7.0938 7.1120c 7.6863 7.5750 7.6845 7.7018c
25 6.1896 7.6563 8.2896 8.2895b
50 7.7403 7.6473 7.7334 9.6683 9.5806 9.6643 9.6949c 10.4447 10.3567 10.4438 10.4726c
62.5 8.3074 10.4022 11.2348 11.2333b
70 8.6085 10.8018 10.7977 10.8306c 11.6577 11.6533 11.6879c
100 9.6191 9.5434 9.6143 12.1214 12.0449 12.1142 12.1584c 13.0635 12.9904 13.0632 13.1055c
125 10.3040 13.0077 14.0163 14.0161b
200 11.8785 11.8111 15.0689 15.0117 16.2018 16.1516
250 12.6878 16.1375 17.3497 17.3495b
500 15.5033 15.4537 19.8595 19.5816 21.2918 21.2521
700 17.0409 21.8937 23.4470
1000 18.7920 18.7519 24.2400 24.2004 25.9295 25.8917
  • Ref. Thirumalai and Heyl (2009)

  • Ref. Ivanov and Schmelcher (2000)

  • Ref. Schimeczek et al. (2013)

  • Ref. Jones et al. (1996a)

Table 2: Absolute value of the binding energies of the positive parity states of helium. Energies are in units of Rydberg energies in the Coulomb potential of nuclear charge for helium. Accurate data from other work is also provided for comparison. (). The values given in parentheses are obtained from a faster version of our spherical atomic structure code developed earlier Heyl and Thirumalai (2010).

It can be seen upon examining the data in Tables 2 and 3, that the fully converged results obtained in the current study are in good agreement with values obtained elsewhere, given that the current study is a single configuration calculation. However the recent results due to Schimeczek et al. (2013) are on average about more bound for . The positive parity states of helium are more bound than the negative parity states in intense magnetic field strengths. For the positive parity states in Table 2, over the entire range of magnetic field strengths investigated, our estimates of the binding energies agree with estimates elsewhere Jones et al. (1999); Becken et al. (1999); Becken and Schmelcher (2000, 2001); Schimeczek et al. (2013) to on average , and , for the states , and , respectively. We noticed loss of accuracy of the cylindrical pseudospectral method in the lower magnetic field regime ( 1) and therefore employed our spherical code for the lower magnetic field strengths while in the range we employed the cylindrical code. We observed that the cylindrical code (and the extrapolation method described above) maintained accuracy to within to in this latter range. In Table 3 the agreement of our results with those of Refs. Jones et al. (1999); Becken et al. (1999); Becken and Schmelcher (2000, 2001); Schimeczek et al. (2013) is , and , for the negative parity states , and , respectively.

One of the aims of the current study is to provide a fast method for the calculation of the energy landscape of atoms in intense magnetic fields; therefore, we have additionally calculated fits to the data provided in Tables 2 and 3. The model fits are rational functions whose analytic form is given by,

(11)

where and . The fitting was carried out using a non-linear least squares Levenberg-Marquardt algorithm with line searches Press et al. (1992). The coefficients and the maximal fitting errors over the entire range to are given in Table 4. These fitting functions could be employed directly in atmosphere models of neutron stars rather than incorporating a code that calculates the binding energy. Thus, atmosphere models which are computationally intensive to begin with, need not be further complicated with the addition of an atomic structure calculation module, even though the software developed in this study is compact and computationally efficient.

iv.2 The Lithium Atom

We investigated the six most tightly bound states of the lithium atom in intense magnetic fields. The binding energies obtained by solving the eigenvalue problem are shown in Tables 5 and 6. These tables show the results for the positive and negative parity states, respectively. As in the case of the helium atom, the HF binding energies are results that were obtained after extrapolating to the limit of infinite domain size.

In contrast to the helium atom, lithium has been investigated far less in the literature, and data is scarce for the binding energies of the different states, particularly in the intense field regime for the negative parity states. Table 5 shows the converged HF binding energies for the three most tightly bound states in the positive parity sub-space of the lithium atom. Over the entire range of magnetic field strengths investigated in the current single configuration study, the results obtained here agree with results obtained elsewhere Al-Hujaj and Schmelcher (2004a); Ivanov and Schmelcher (2000); Ivanov (1998); Schimeczek et al. (2013) to on average , and for the positive parity states , and , respectively.

Table 6 shows the converged HF binding energies of the negative parity states of lithium. Of the three most tightly bound states in this parity sub-space, only one has been investigated in the literature so far; the state . The average agreement between our results and those of Al-Hujaj and Schmelcher (2004a); Ivanov and Schmelcher (2000) is on average , over the entire range of magnetic field strengths considered in this study. We observed that the extrapolation to the limit of infinite domain size maintained accuracy to within to over the entire range of magnetic field strengths. Moreover, since our results are upper bounds for the uncorrelated single-configuration binding energies, the data for the states and given in Table 6 can be improved using post-HF methods. The D wave functions of these states calculated herein could be employed as initial estimates for improvements.

Moreover, with increasing magnetic field strength in the intense magnetic field regime, the most tightly bound negative parity state of lithium is seen to be the state , which is comprised of the field-free orbitals . This crossover occurs at around , below this field the is the most tightly bound of the three negative parity states of lithium shown in Table 6. To the best of our knowledge, this crossover has not been reported elsewhere in the literature. In addition, it can be seen by comparing the binding energies reported in Tables 5 and 6, that the state , is also among the most tightly bound states of lithium in intense magnetic fields. Moreover, the third state shown therein, the state, comprised of the orbitals , also has not been investigated in the literature. This latter state also becomes tightly bound with increasing magnetic field strength in the intense field regime.

Thus overall we see that for the six most tightly bound states of lithium in intense magnetic fields, two of the states have not been investigated earlier at all ( and ) and a third state () has not been investigated in the intense field regime (). Therefore the results presented here appear to be the first of such studies and represent upper bounds to the uncorrelated single-configuration binding energies. In addition we see that for the remaining three states that were investigated, the binding energies obtained in the current study are in relative good agreement with estimates obtained elsewhere at the level overall.

Furthermore, for the sake of facilitating atmosphere and crustal models of neutron stars, we have also carried out rational function fits to the data. Once again these analytic forms can be implemented directly in such codes, thereby circumventing the need for atomic structure calculations altogether. The rational functions have the same functional form as those described in Eq. (11) above. The coefficients of these rational functions are given in Table 7 alongside estimates of the fitting errors.

In the following section, we summarize the findings alongside a brief discussion of further avenues for investigation.

Refs. Refs. Refs. Refs. Refs. Refs. Refs. Refs. Refs.
Here Ruder et al. (1994)/Thirumalai and Heyl (2009) Jones et al. (1999) Becken et al. (1999)/Schimeczek et al. (2013) Here Ruder et al. (1994)/Thirumalai and Heyl (2009) Jones et al. (1999) Becken and Schmelcher (2000) Here Ruder et al. (1994)/Thirumalai and Heyl (2009) Jones et al. (1999) Becken and Schmelcher (2001)
0 (1.0686) 1.0641a 1.0657c 1.0665 (1.0307) 1.0278 (1.0189) 1.0156
0.01 (1.1056) 1.1029a 1.1016 1.1026 (1.0745) 1.0693 1.0702 1.0705 (1.0643) 1.0579 1.0601 1.0603
0.05 (1.2147) 1.2135a 1.2099 1.2112 (1.1808) 1.1565 1.1756 1.1761 (1.1684) 1.1364 1.1634 1.1636
0.1 (1.3229) 1.3231a 1.3174 1.3191 (1.2846) 1.2270 1.2764 1.2795 (1.2709) 1.1979 1.2650 1.2655
0.125 (1.3727) 1.3669 (1.3308) 1.3255 (1.3166) 1.3110
0.2 (1.4982) 1.4988a 1.4914 1.4936 [1.4938]c (1.4539) 1.3174 1.4469 1.4481 (1.4384) 1.2723 1.4319 1.4326
0.25 (1.5721) 1.5676 [1.5677]c (1.5262) 1.5202 (1.5044) 1.5042
0.5 1.8648 1.8647a 1.8593 1.8131 1.5131 1.8078 1.7918 1.4937 1.7912
0.625 1.9798 1.9796 1.9263 1.9259 1.9092 1.9090
0.7 2.0465 2.0461a 2.0404 1.9935 1.7202 1.9877 1.9717 1.7012 1.9709
1 2.2689 2.2685a 2.2641 2.2159 1.9678 2.2097 2.1939 1.9493 2.1931
1.25 2.4245 2.4243 2.3689 2.3687 2.3524 2.3522
2 2.8069 2.8060a 2.7989 2.7502 2.5433 2.7441 2.7295 2.5269 2.7287
2.5 3.0059 3.0057 [3.0058]c 2.9512 2.9511 2.9364 2.9361
5 3.7524 3.7522a 3.7466 3.7004 3.5369 3.6950 3.6836 3.5242 3.6827
6.25 4.0297 4.0297 3.9799 3.9795 3.9684 3.9680
7 4.1819 4.1817a 4.1762 4.1315 4.1264 4.1163 4.1155
10 4.6921 4.6920a 4.6862 4.6434 4.5067 4.6387 4.6304 4.4971 4.6275
12.5 5.0406 5.0400  [5.0401]c 4.9947 4.9947 4.9866 4.9860
20 5.8589 5.7495 5.8585 5.8195 5.7056 5.8158 5.8103 5.6989 5.8091
25 6.2942 6.2940b 6.2545 6.2462
50 7.8349 7.7491 7.8346 7.8016 7.7118 7.7958 7.7962 7.7080 7.7940
62.5 8.3991 8.3682 8.3632
70 8.7004 8.6998 8.6684 8.6651 8.6646 8.6620
100 9.7078 9.6370 9.7075 9.6799 9.6039 9.6748 9.6756 9.6015 9.6724
125 10.3885 10.3617 10.3586
200 11.9501 11.8970 11.9339 11.8675 11.9330 11.8661
250 12.7593 12.7427 12.7398
500 15.5573 15.5307 15.5449 15.5050 15.5290 15.5043
700 17.0580 17.0524 17.0427
1000 18.8466 18.8230 18.8275 18.7997 18.8198 18.7993
  • Ref. Thirumalai and Heyl (2009)

  • Ref. Ivanov and Schmelcher (2000)

  • Ref. Schimeczek et al. (2013)

Table 3: Absolute value of the binding energies of the negative parity states of helium. Energies are in units of Rydberg energies in the Coulomb potential of nuclear charge for helium. Accurate data from other work is also provided for comparison. ().
State Coefficients State Coefficients
Table 4: Coefficients of the different rational functions for fitting the six states of helium discussed. The absolute maximum fractional error of the eigenvalue relative to the fit from to is reported in the variable
Here Ref. Al-Hujaj and Schmelcher (2004a) Here Elsewhere Here Elsewhere
0 (1.1492) 1.1491 (1.1968) 1.1926a (1.1357) 1.1427a  [1.1299]d
0.00056 (1.1541) 1.1544 (1.2024) 1.1969a (1.1425) 1.1487a
0.0028 (1.1780) 1.1720 (1.2194) 1.2121c (1.1652) 1.1663a
0.0056 (1.1961) 1.1901 (1.2390) 1.2334a (1.1897) 1.1869a
0.0111 (1.2267) 1.2203 (1.2735) 1.2674a (1.2324) 1.2278a
0.0278 (1.2964) 1.2886 (1.3530) 1.3463a (1.3354) 1.3294a
0.0556 (1.3912) 1.3930 (1.4529) 1.4432a (1.4699) 1.4627a
0.5 2.1745 2.3607 2.5903
0.5556 2.4146 2.4145 2.4284 2.4280a 2.6579 2.6572a
1 2.7158 2.9551 3.2818
1.1111 2.8124 3.0619 3.0432b 3.4056 3.3695b
2 3.4404 3.7473 4.1986
2.3636 4.4581 4.4133b  [4.4203]e
2.7778 3.8590 4.2030 4.1781b 4.7247 4.6779b
5 4.7467 5.1669 5.8362
5.5556 4.9262 5.3616 5.3304c 6.0606 6.0043b
7 5.3444 5.8140 6.5822
10 6.0571 6.5846 7.4703
11.1111 6.2835 6.8301 6.7909c 7.7532 7.6856c
11.8178 7.9232 7.8544b  [7.8711]e
20 7.7046 8.3606 9.5172
23.6356 10.0807 9.9989b  [10.0238]e
27.7778 8.6194 9.3450 9.2936c 10.6480 10.5685b
50 10.5017 11.3643 12.9842
55.5556 10.8742 11.7638 11.7000c 13.4456 13.3464b
70 11.7297 12.6825 14.5023
100 13.1685 14.2197 16.2849
118.1780 17.1804 17.0622b  [17.1231]e
200 16.3929 17.6559 20.2679
277.7778 18.1335 19.5079 22.4181 22.2774b
500 21.6227 23.2111 26.7255
555.5556 22.2335 23.9240 27.5636 27.4029b
700 23.8351 25.5639 29.4767
1000 26.3906 28.2108 32.6287
  • Ref. Al-Hujaj and Schmelcher (2004a)

  • Ref. Ivanov and Schmelcher (2000)

  • Ref. Ivanov (1998)

  • Ref. Ralchenko et al. (2008), experimental result.

  • Ref. Schimeczek et al. (2013)

Table 5: Absolute value of the binding energies of the positive parity states of lithium. Energies are in units of Rydberg energies in the Coulomb potential of nuclear charge for lithium. Accurate data from other work is also provided for comparison. ().
Here Elsewhere Here Here
0 (1.1687) 1.1652c (1.1400) (1.1306)
0.00056 (1.1735) 1.1695c (1.1457) (1.1361)
0.0028 (1.1909) 1.1865c (1.1644) (1.1556)
0.0056 (1.2110) 1.2065c (1.1827) (1.1755)
0.0111 (1.2469) 1.2417c (1.2183) (1.2104)
0.0278 (1.3354) 1.3297c (1.3073) (1.2964)
0.0556 (1.4500) 1.4463c (1.4248) (1.4113)
0.5 2.4358 2.4048 2.3804
0.5556 2.4899 2.4898c 2.4853 2.4605
1 3.0313 3.0074 2.9807
1.1111 3.1376 3.1035a 3.1150 3.0881
2 3.8180 3.8029 3.7757
2.7778 4.2693 4.2319a 4.2585 4.2316
5 5.2238 5.2196 5.1945
5.5556 5.4168 5.3767b 5.4136 5.3889
7 5.8654 5.8641 5.8405
10 6.6301 6.6313 6.6095
11.1111 6.8738 6.8298b 6.8757 6.8544
20 8.3958 8.3998 8.3819
27.7778 9.3762 9.3242a 9.3808 9.3646
50 11.3892 11.3941 11.3809
55.5556 11.7877 11.7269a 11.7922 11.7799
70 12.7043 12.7090 12.6974
100 14.2387 14.2445 14.2331
200 17.6700 17.6795 17.6601
277.7778 19.5196 19.5324 19.5173
500 23.2201 23.2339 23.2195
555.5556 23.9385 23.9592 23.9414
700 25.5784 25.6165 25.5923
1000 28.2764 28.3062 28.2948
  • Ref. Ivanov (1998)

  • Ref. Ivanov and Schmelcher (2000)

  • Ref. Al-Hujaj and Schmelcher (2004a)

Table 6: Absolute value of the binding energies of the negative parity states of lithium. Energies are in units of Rydberg energies in the Coulomb potential of nuclear charge for lithium. Accurate data from other work is also provided for comparison. ().
State Coefficients State Coefficients
Table 7: Coefficients of the different rational functions for fitting the six states of lithium discussed. The absolute maximum fractional error of the eigenvalue relative to the fit from to is reported in the variable

V Conclusion

In the current study we have investigated low-Z atoms, helium and lithium in intense magnetic fields. A two-dimensional single-configuration Hartree-Fock method, as described in Ref. Thirumalai and Heyl (2009), was adopted. A key feature of the method is that the potentials for the inter-electronic interactions are obtained as solutions to the elliptic partial differential equations as given in Eqs. (3) and (4). The HF equations in Eq. (1) are solved using the self consistent field method. The system size grew as , with the number of electrons in the coupled problem and the number of grid points in each direction.

A pseudospectral approach was adopted for the numerical solution of the problem using cylindrical coordinates so as to facilitate calculations in the intense field regime. Domain discretization was achieved using the commonly employed Chebyshev-Lobatto spectral collocation method. The resulting discretized and coupled eigenvalue problem problem was solved using standard sparse matrix methods using the software package ARPACK. The key enabling advantage of the psuedospectral approach is the immensely reduced computational time particularly, since we have adopted here an unrestricted two-dimensional approach to the problem Thirumalai and Heyl (2009). The latter has the advantage that it does not require a basis of functions to describe the wave functions. Thus the wave functions obtained in the current unrestricted 2D approach can effectively be thought of as those arising from the superposition of a large number of basis functions.

We presented data for the six most tightly bound states of the helium atom in intense magnetic fields, in Tables 2 and 3. These were seen to be consistent with findings elsewhere. Similarly we investigated the six most tightly bound states of the lithium atom as well. However, we found that the data in the literature to be rather scarce for lithium. As a result we could only compare our results for four of the six states characterized in this study. We obtained, apparently for the first time, calculations for the binding energies for the states and of lithium. We find that the the state is also the most tightly bound negative parity state of lithium in the limit of intense magnetic fields.

The work described herein was motivated primarily by the need to have accurately determined upper bounds for the binding energies of atoms in intense magnetic fields employing a computationally straight-forward implementation. As the atomic structure software developed here is compact and computationally economical, it produces accurate results within a short amount of computing time. As a result, it can be incorporated directly into atmosphere and crustal models for neutron stars. However, while this may be desirable, it may present an additional layer of computational complexity. The user may wish to circumvent this by employing the rational function fits given in Tables 4 and 7. These analytic forms, model the data in the range and thus may simplify atmosphere and crustal models considerably. Estimates of binding energies and oscillator strengths are ultimately needed to correctly interpret the spectra of neutron stars and magnetized white dwarfs. However, there is an additional complication that in order to do so, it also becomes necessary to account for the effect of strong electric fields which are also present in the atmospheres of these objects. Therefore, the energies and wave functions obtained herein are only part of the solution, and there is yet work to be done before this goal can be achieved. Moreover, as the magnetic field strength increases in the intense magnetic field regime, effects due to finite nuclear mass become increasingly relevant. In the current study, the mass of the nucleus is assumed to be infinite, and as such we have not carried out a suitable correction. One way to account for finite nuclear mass effects is to employ a scaling relationship wherein the energies determined at a certain magnetic field strength for an infinite nuclear mass, would be related to the corresponding binding energies for a finite nuclear mass at a different value of the magnetic field strength Becken et al. (1999). Such a correction becomes increasingly important at that upper end of magnetic field strengths investigated in this study. However, while accounting for this correction is important, it was not possible to do so since the errors of the pseudospectral method at such field strengths was which is typically greater than the magnitude of the correction Becken et al. (1999). Additionally, the software developed here could be extended to tackle atoms with greater number of electrons, such as carbon or oxygen. It can also be extended towards a 2D post-HF framework such as CI of MCHF which would no doubt yield more accurate results. In either case, the method developed herein would be central to such enhancements and as such, the current study represents the very first implementation of a cylindrical pseudospectral method for atomic structure calculations in intense magnetic fields.

*

Appendix A The Pseudo-spectral Approach

For the numerical solution of the HF equations we employ a discretization based on pseudospectral methods extending our earlier investigation Heyl and Thirumalai (2010) to intense magnetic fields. However in contrast to our earlier work, in the current study we employ a cylindrical coordinate system. As a result the methodology for setting up the problem is considerably different from that described in Ref. Heyl and Thirumalai (2010). This section is arranged as follows. First we describe the methodology employed for solving the hydrogen atom using pseudospectral methods in cylindrical coordinates. Particular emphasis is placed on the implementation of boundary conditions. Thereafter, the treatment is extended to the particular case of the helium atom in a single configuration and a generalization of the scheme is then provided for multi-electron atoms. In the sections that follow, a considerable amount of detail is provided for the benefit of readers not familiar with pseudospectral methods, others may wish to only read Section A.2 of the appendix.

a.1 The Hydrogenic Problem

We begin with the Hamiltonian for the hydrogenic problem (single-electron) in a strong magnetic field,

(12)

The solution of the eigenvalue problem in Eq. (12) yields the individual electron energies and wave functions in a given configuration. The domains of both the radial and the axial coordinates are . The problem maintains azimuthal symmetry and thus, a solution of Eq. (12) in this domain, when reflected about the plane (respecting parity of course) and revolved about the z-axis through , gives the solution in three-dimensional cylindrical coordinates.

With the transformations given in Section III in Eqs. (8) and (9), we can re-write Eq. (12) as

(13)

where, we have dropped the subscripts on the coordinate labels. For brevity we have introduced the constant . The discretisation points are thereafter taken to be the commonly used Chebyshev-Lobatto points Heyl and Thirumalai (2010); Trefethen (2000); Boyd (2001) given by

(14)

where . As is customary, we employ monic polynomials of degree as the cardinal functions to interpolate between these points and are given by Heyl and Thirumalai (2010); Trefethen (2000),

(15)

with

(16)

Derivatives of these interpolating polynomials at the discretisation points then yield the so-called Chebyshev differntiation matrix, whose elements are given by,

(17)

and

(18)

In writing Eq. (13), we have removed the co-ordinate singularity at , by replacing it with , where in units of Bohr radii. This approximation produced acceptable results within error tolerances. The outer boundary conditions of the domain, at (corresponding to ), are taken care of by imposing Dirichlet boundary conditions since the wave function must vanish at infinity (see below). The remaining inner boundaries of the compactified domain at can have either Dirichlet or Neumann boundary conditions, depending upon the wave function in question. The following discussion delineates the methodology for the 2D problem.

a.1.1 An Explicit Example - Domain Discretization

We consider here an explicit example to illustrate the use of pseudospectral methods for solving an eigenvalue problem. The method developed here is a non-trivial extension of the one developed by the authors in Ref. Talbot and Crampton (2005).

Let us begin with a domain which is discretised using points in each of the two Cartesian directions; and , with in this explicit example. This is illustrated in Fig. (4)

Figure 4: Pictorial representation of the domain discretised using points in each direction, with . The outer boundaries have Dirichlet conditions imposed.

The partial differential equation in Eq. (13) is two-dimensional therefore, following the procedure outlined in Refs. Heyl and Thirumalai (2010); Talbot and Crampton (2005), we can construct two-dimensional operators by employing Kronecker products of matrices, for example we can define,

(19)

and

(20)

where and are Chebyshev differentiation matrices in the and directions respectively, defined similar to Eqs. 17 and 18. The matrices and are identity matrices of dimension and respectively. The matrices and are discrete representations of the and operators, respectively. In these matrices the entries corresponding to the outer boundaries at have already been excised due to Dirichlet boundary conditions (Heyl and Thirumalai, 2010; Trefethen, 2000, see Figure 4). As a result, the indices of these matrices are limited to an upper value of rather than . In our example, the number of node points in the and -directions are equal, therefore . Thus, using these Chebyshev differentiation matrices we can write down Eq. (13) in matrix form as,

(21)

where