A twodimensional pseudospectral HartreeFock method for lowZ atoms in intense magnetic fields
Abstract
The energy levels of the first few lowlying 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 twodimensional HartreeFock partial differential equations for these two and threeelectron systems in a selfconsistent 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 finiteelement 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 Zeemantype 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); AlHujaj 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 singleconfiguration HartreeFock (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 finiteelement based approach Thirumalai and Heyl (2009).
In sharp contrast to the somewhat simpler twoelectron 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 ThomasFermiDirac 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 socalled 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 lowlying 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 postHF studies of atoms with more than two electrons using a numerical meshmethod for solving the unrestricted HF equations Ivanov and Schmelcher (1998). The special meshes were so constructed as to facilitate finitedifference calculations in a twodimensional domain using carefully selected mesh node points Ivanov and Schmelcher (2001c). They were able to ascertain the binding energies of the first few lowlying states of lowZ atoms such as lithium, beryllium and midZ 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 AlHujaj and Schmelcher (2003, 2003); Becken and Schmelcher (2002, 2001, 2000); Becken et al. (1999); Luhr et al. (2007), adopting a full configurationinteraction method, AlHujaj & Schmelcher AlHujaj 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 GonzalezFerez & Schmelcher GonzálezFé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 configurationinteraction method, employing the socalled freezing fullcore 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 densityfunctionaltheory. Mori et al Mori and Hailey (2002); Mori and Ho (2007) have studied midZ 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 coworkers 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 finiteelement methods with Bsplines both in the adiabatic approximation and beyond the adiabatic approximation with more than one Landau level. These highly accurate formulations employ a fast parallel HartreeFockRoothan 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 MonteCarlo 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 fixedphase Quantum MonteCarlo 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 runofthemill computer architectures. The method described here provides a speedup of a factor of when compared with our earlier investigation using a finiteelement discretization Thirumalai and Heyl (2009). Computational times are comparable to the recently developed quantum montecarlo methods due to Schimeczek et al. (2013). The method described in this study can be run on a distributed architecture for investigating multielectron atoms.
It is also well known that postHF methods such as configuration interaction (CI) or multiconfigration HartreeFock (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. AlHujaj 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 postHF methods considerably. However, there do not exist hitherto, any fully twodimensional (2D) postHF studies of multielectron 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 singleconfiguration problem in a full 2D framework without any basis expansions and separation of variables. Wave functions so determined could be used directly in 2D configurationinteraction calculations or the problem could be cast into a MCHF framework. Moreover, obtaining accurate estimates of the energy levels of atoms, in particular lowZ 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 boundfree and freefree 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 lowlying states of the simplest lowZ atoms; helium and lithium, in strong and intense magnetic fields using a two dimensional singleconfiguration pseudospectral method of solution.
Ii The HF Equations
We shall begin with the generalised singleconfiguration HF equations for an atom with electrons and nuclear charge , in a magnetic field that is oriented along the direction. A derivation of the singleconfiguration 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 threedimensional momentum operator has been split into two parts; which are the and parts of the Laplacian and which is the azimuthal part. The total HartreeFock 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 antisymmetrization operator. The individual electronic wave functions are given by,
(6) 
where labels each of the electrons. The twodimensional 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 spinpolarised states (FSP), in other words all the electrons of the atom are assumed to be antialigned 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 spinpolarised 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 HartreeFock 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 IIIV. 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 socalled 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 wavefunctions 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 semiinfinite 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 ChebyshevLobatto 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 twodimensional 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 shiftinvert 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 HartreeFock 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 precalculating the interelectron 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 twodimensional one, and the eigenvalue problem is therefore a fully coupled twodimensional problem. Additionally, the interelectron 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 .
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 midZ 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 ChebyshevLobatto 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 semiinfinite 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 fieldfree and strongfield 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 subspaces exist. The quantum number counts the excitation level within a given manifold and subspace 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 subspace, 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 subspace, we considered the most tightly bound state in the intense field regime.
Intensefield  Fieldfree  

Helium  
Lithium  
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 singleconfiguration 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 postHF 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 singleconfiguration calculation and we leave such a full D postHF 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 .
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 subspace. 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.0871^{d}  1.0876  (1.0307)  1.0179^{a}  1.0278  (1.0686)  1.0668^{a}  1.0657^{c}  1.0666  
0.01  (1.1256)  1.1213  1.1213  1.1220  (1.0872)  1.0852^{a}  1.0830  1.0833  (1.1223)  1.1183^{a}  1.1177  1.1193 
0.05  (1.2103)  1.1911  1.2056  1.2064  (1.2212)  1.2175^{a}  1.2160  1.2167  (1.2734)  1.2691^{a}  1.2683  1.2704 
0.1  (1.2914)  1.2133  1.2860  1.2868  (1.3495)  1.3510^{a}  1.3436  1.3450  (1.4209)  1.4189^{a}  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.5598^{a}  1.5525  (1.6579)  1.6585^{a}  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.0009^{a}  1.9945  (2.1582)  2.1550^{a}  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.2246^{a}  2.2171  2.2210^{c}  (2.4060)  2.4029^{a}  2.3956  
1  2.1607  1.9116  2.1601  2.4983  2.4981^{a}  2.4933  2.4978^{c}  2.7047  2.7026^{a}  2.7000  2.7043^{c}  
1.25  2.3141  2.3137  2.6943  2.6940  2.9198  2.9197  
2  2.6863  2.4793  2.6840  3.1691  3.1685^{a}  3.1634  3.1691^{c}  3.4387  3.4384^{a}  3.4333  3.4386^{c}  
2.5  2.8868  2.8862  3.4275  3.4274  3.7203  3.7203  
5  3.6311  3.4672  3.6272  4.3785  4.3740^{a}  4.3693  4.3783^{c}  4.7539  4.7502^{a}  4.7441  4.7523^{c}  
6.25  3.9082  3.9076  4.7383  4.7380  5.1425  5.1421  
7  4.0618  4.9310  4.9268^{a}  4.9203  4.9310^{c}  5.3514  5.3474^{a}  5.3408  5.3507^{c}  
10  4.5742  4.4362  4.5693  5.5856  5.5851^{a}  5.5770  5.5899^{c}  6.0548  6.0543^{a}  6.0506  6.0624^{c}  
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.1120^{c}  7.6863  7.5750  7.6845  7.7018^{c}  
25  6.1896  7.6563  8.2896  8.2895^{b}  
50  7.7403  7.6473  7.7334  9.6683  9.5806  9.6643  9.6949^{c}  10.4447  10.3567  10.4438  10.4726^{c}  
62.5  8.3074  10.4022  11.2348  11.2333^{b}  
70  8.6085  10.8018  10.7977  10.8306^{c}  11.6577  11.6533  11.6879^{c}  
100  9.6191  9.5434  9.6143  12.1214  12.0449  12.1142  12.1584^{c}  13.0635  12.9904  13.0632  13.1055^{c}  
125  10.3040  13.0077  14.0163  14.0161^{b}  
200  11.8785  11.8111  15.0689  15.0117  16.2018  16.1516  
250  12.6878  16.1375  17.3497  17.3495^{b}  
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)
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 nonlinear least squares LevenbergMarquardt 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 subspace 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 AlHujaj 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 subspace, only one has been investigated in the literature so far; the state . The average agreement between our results and those of AlHujaj 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 singleconfiguration binding energies, the data for the states and given in Table 6 can be improved using postHF 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 fieldfree 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 singleconfiguration 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.0641^{a}  1.0657^{c}  1.0665  (1.0307)  1.0278  (1.0189)  1.0156  
0.01  (1.1056)  1.1029^{a}  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.2135^{a}  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.3231^{a}  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.4988^{a}  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.8647^{a}  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.0461^{a}  2.0404  1.9935  1.7202  1.9877  1.9717  1.7012  1.9709  
1  2.2689  2.2685^{a}  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.8060^{a}  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.7522^{a}  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.1817^{a}  4.1762  4.1315  4.1264  4.1163  4.1155  
10  4.6921  4.6920^{a}  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.2940^{b}  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 
State  Coefficients  State  Coefficients 

Here  Ref. AlHujaj and Schmelcher (2004a)  Here  Elsewhere  Here  Elsewhere  

0  (1.1492)  1.1491  (1.1968)  1.1926^{a}  (1.1357)  1.1427^{a} [1.1299]^{d} 
0.00056  (1.1541)  1.1544  (1.2024)  1.1969^{a}  (1.1425)  1.1487^{a} 
0.0028  (1.1780)  1.1720  (1.2194)  1.2121^{c}  (1.1652)  1.1663^{a} 
0.0056  (1.1961)  1.1901  (1.2390)  1.2334^{a}  (1.1897)  1.1869^{a} 
0.0111  (1.2267)  1.2203  (1.2735)  1.2674^{a}  (1.2324)  1.2278^{a} 
0.0278  (1.2964)  1.2886  (1.3530)  1.3463^{a}  (1.3354)  1.3294^{a} 
0.0556  (1.3912)  1.3930  (1.4529)  1.4432^{a}  (1.4699)  1.4627^{a} 
0.5  2.1745  2.3607  2.5903  
0.5556  2.4146  2.4145  2.4284  2.4280^{a}  2.6579  2.6572^{a} 
1  2.7158  2.9551  3.2818  
1.1111  2.8124  3.0619  3.0432^{b}  3.4056  3.3695^{b}  
2  3.4404  3.7473  4.1986  
2.3636  4.4581  4.4133^{b} [4.4203]^{e}  
2.7778  3.8590  4.2030  4.1781^{b}  4.7247  4.6779^{b}  
5  4.7467  5.1669  5.8362  
5.5556  4.9262  5.3616  5.3304^{c}  6.0606  6.0043^{b}  
7  5.3444  5.8140  6.5822  
10  6.0571  6.5846  7.4703  
11.1111  6.2835  6.8301  6.7909^{c}  7.7532  7.6856^{c}  
11.8178  7.9232  7.8544^{b} [7.8711]^{e}  
20  7.7046  8.3606  9.5172  
23.6356  10.0807  9.9989^{b} [10.0238]^{e}  
27.7778  8.6194  9.3450  9.2936^{c}  10.6480  10.5685^{b}  
50  10.5017  11.3643  12.9842  
55.5556  10.8742  11.7638  11.7000^{c}  13.4456  13.3464^{b}  
70  11.7297  12.6825  14.5023  
100  13.1685  14.2197  16.2849  
118.1780  17.1804  17.0622^{b} [17.1231]^{e}  
200  16.3929  17.6559  20.2679  
277.7778  18.1335  19.5079  22.4181  22.2774^{b}  
500  21.6227  23.2111  26.7255  
555.5556  22.2335  23.9240  27.5636  27.4029^{b}  
700  23.8351  25.5639  29.4767  
1000  26.3906  28.2108  32.6287 

Ref. AlHujaj and Schmelcher (2004a)

Ref. Ivanov and Schmelcher (2000)

Ref. Ivanov (1998)

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

Ref. Schimeczek et al. (2013)
Here  Elsewhere  Here  Here  

0  (1.1687)  1.1652^{c}  (1.1400)  (1.1306) 
0.00056  (1.1735)  1.1695^{c}  (1.1457)  (1.1361) 
0.0028  (1.1909)  1.1865^{c}  (1.1644)  (1.1556) 
0.0056  (1.2110)  1.2065^{c}  (1.1827)  (1.1755) 
0.0111  (1.2469)  1.2417^{c}  (1.2183)  (1.2104) 
0.0278  (1.3354)  1.3297^{c}  (1.3073)  (1.2964) 
0.0556  (1.4500)  1.4463^{c}  (1.4248)  (1.4113) 
0.5  2.4358  2.4048  2.3804  
0.5556  2.4899  2.4898^{c}  2.4853  2.4605 
1  3.0313  3.0074  2.9807  
1.1111  3.1376  3.1035^{a}  3.1150  3.0881 
2  3.8180  3.8029  3.7757  
2.7778  4.2693  4.2319^{a}  4.2585  4.2316 
5  5.2238  5.2196  5.1945  
5.5556  5.4168  5.3767^{b}  5.4136  5.3889 
7  5.8654  5.8641  5.8405  
10  6.6301  6.6313  6.6095  
11.1111  6.8738  6.8298^{b}  6.8757  6.8544 
20  8.3958  8.3998  8.3819  
27.7778  9.3762  9.3242^{a}  9.3808  9.3646 
50  11.3892  11.3941  11.3809  
55.5556  11.7877  11.7269^{a}  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 
State  Coefficients  State  Coefficients 

V Conclusion
In the current study we have investigated lowZ atoms, helium and lithium in intense magnetic fields. A twodimensional singleconfiguration HartreeFock method, as described in Ref. Thirumalai and Heyl (2009), was adopted. A key feature of the method is that the potentials for the interelectronic 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 ChebyshevLobatto 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 twodimensional 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 straightforward 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 postHF 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 Pseudospectral 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 multielectron 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 (singleelectron) 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 zaxis through , gives the solution in threedimensional cylindrical coordinates.
With the transformations given in Section III in Eqs. (8) and (9), we can rewrite 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 ChebyshevLobatto 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 socalled Chebyshev differntiation matrix, whose elements are given by,
(17) 
and
(18) 
In writing Eq. (13), we have removed the coordinate 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 nontrivial 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)
The partial differential equation in Eq. (13) is twodimensional therefore, following the procedure outlined in Refs. Heyl and Thirumalai (2010); Talbot and Crampton (2005), we can construct twodimensional 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