High-Harmonic Generation in the Water Window from mid-IR Laser Sources
Keegan Finger, David Atri-Schuller, Nicolas Douguet, Klaus Bartschat, Kathryn R. Hamilton
HHigh-Harmonic Generation in the Water Window from mid-IR Laser Sources
Keegan Finger, David Atri-Schuller, Nicolas Douguet, Klaus Bartschat, and Kathryn R. Hamilton ∗ Department of Physics and Astronomy, Drake University, Des Moines, Iowa 50311, USA Department of Physics, Kennesaw State University, Kennesaw, Georgia 30144, USA (Dated: January 28, 2021)We investigate the harmonic response of neon atoms to mid-IR laser fields (2000 − ab initio all-electron R -Matrix with Time-dependence (RMT) method. The laser peak intensity and wavelength are varied to find optimalparameters for high-harmonic imaging in the water window. Comparison of the SAE and RMTresults shows qualitative agreement between them as well as parameters such as the cutoff frequencypredicted by the classical three-step model. However, there are significant differences in the details,thus indicating the importance of multi-electron effects. I. INTRODUCTION
The “Water Window” is a region of the electro-magnetic (EM) spectrum spanning from the K -absorption edge of carbon at 282 eV to the K -edge ofoxygen at 533 eV. Importantly, water is nearly transpar-ent to soft X-rays in this region, while carbon, and thusorganic molecules, are absorbing. This enables the imag-ing of biologically important molecules in their typicalaqueous environments [1].Light sources in this region of the EM spectrum aregenerally produced by synchrotron radiation from an ac-celerator [2]. However, a major problem with this ap-proach is the inability to conduct time-dependent imag-ing of live cells without freezing due to the typically lowinstantaneous power of synchrotron radiation [3]. Analternative approach is high-order harmonic generation(HHG), which employs a few-cycle laser system to gener-ate coherent soft X-rays by producing odd harmonics ofa mid-range infrared (IR) fundamental frequency [4–6].Such HHG sources could lead to the generation of single-shot absorption spectra, and live-cell imaging with a fem-tosecond time resolution [7].HHG can be understood in terms of the semi-classical“three-step model” suggested by Corkum [8]. In thismodel, an electron is 1) tunnel-ionized and 2) acceleratedby a strong electric field. As the electric field changessign, the electron is accelerated back towards its parention, where it 3) recollides and releases some of its energyin the form of a high-energy photon. The recollision canalso excite another electron with a higher energy, and theprocess can be repeated over several cycles of the electricfield.From the above model, one can obtain a formula forthe maximum photon energy, or cutoff energy, producibleby HHG from a single laser source. This cutoff energy isgiven by E c = I p + 3 . U p , (1) ∗ [email protected] where I p is the atomic ionization potential and U p is theponderomotive potential. The ponderomotive potentialcan be approximated in terms of the laser peak intensityand wavelength as U p = 2 e c(cid:15) m · I · λ π . (2)Here e is the electron charge, (cid:15) the vacuum permittivity, c the speed of light, m the electron mass, I the laser peakintensity, and λ the (central) laser wavelength. The pon-deromotive potential scales linearly with intensity andquadratically with wavelength.At first sight, one might just increase the intensity orthe wavelength to generate harmonics in the water win-dow. However, there are a number of problems with thisassumption. Experimentally, increasing the intensity re-sults in more ionization, which leads to depletion of thetarget. On the other hand, if the wavelength is increasedinstead, the conversion efficiency lowers, following ap-proximately a λ − or even λ − scaling [9]. Furthermore,long-wavelength lasers are relatively rare compared totheir mid-IR counterparts.The computational side of HHG is not without chal-lenges either. The largest of these is an extensive angular-momentum expansion needed to obtain partial-wave con-vergence in a fully quantum-mechanical treatment of theprocess. The computational load for the calculation re-quires large computers and significant run times. Fur-ther, longer wavelengths mean larger excursion distancesfor the electron, thus requiring a large configuration spaceto correctly describe the electron’s motion far away fromthe nucleus.These aforementioned computational challenges aretypically overcome by applying some approximations toreduce the overall complexity of the problem. They in-clude the semi-classical strong-field approximation (SFA)[10] or the single-active electron (SAE) models [11], bothof which neglect multi-electron effects. However, sucheffects have been shown to significantly change the har-monic cutoff and/or the intensity of the generated har-monics through phenomena such as the giant resonancein xenon [12].The R -Matrix with Time-dependence (RMT) method a r X i v : . [ phy s i c s . a t o m - ph ] J a n accounts for the multi-electron effects in HHG by employ-ing a highly efficient code specifically designed for high-performance computing [13]. RMT is indeed suitable forserious computational challenges such as HHG from mid-IR lasers, and has previously been used to describe HHGfrom 1800 nm sources [14]. In the present project, wesignificantly extend previous RMT calculations by em-ploying driving laser wavelengths of more than 2000 nm,resulting in harmonic cutoffs within the water window.Harmonic spectra from RMT calculations are also com-pared with those obtained by an SAE approach, provid-ing us with means to determine the influence of multi-electron effects on the spectra in this energy region.Section II gives a brief overview of the numerical meth-ods and the parameter space investigated. This is fol-lowed by the presentation and discussion of the predictedharmonic spectra for a neon atom in laser fields of vari-ous intensities and wavelengths in Sect. III. Conclusionsand an outlook are given in Sect. IV. Unless specifiedotherwise, atomic units (a.u.) are used below. II. NUMERICAL METHODS
Below we describe the methods used in the presentstudy. We start with a brief summary of our own SAEapproach in Section II A, which is followed by a more de-tailed description of the RMT method in Section II B.Section II C summarizes the parameter space investi-gated.
A. The Single-Active-Electron Approach
We use the same basic method and the associated com-puter code as described by Pauly et al. [15]. Briefly, wesolve the time-dependent Schr¨odinger equation (TDSE)for the active electron, which is initially in the 2 p orbitalof the neon ground state. We use the potential suggestedby Tong and Lin [16], in which we calculate the 2 p andother bound orbitals (only used here to check the qualityof the structure description), as well as the continuumorbitals to represent the ionized electrons. The abovepotential produces the first ionization potential of neon,i.e., the binding energy of the 2 p electron, very well. Weemploy a finite-difference method with a variable radialgrid, with the smallest stepsize of 0.01 near the nucleusand a largest stepsize of 0.05 for large distances. Thetimestep is held constant at 0.005. The calculations arecarried out in the length gauge of the electric dipole op-erator. Partial waves up to angular momenta (cid:96) = 300were included in order to ensure converged results. Testscarried out by varying the time step, the radial grid, andthe number of partial waves give us confidence that anyinaccuracies in the predictions are due to the inherentlimitations of the SAE model rather than numerical is-sues. We then calculate z ( t ) ≡ (cid:104) Ψ( r , t ) | z | Ψ( r , t ) (cid:105) , i.e., effec-tively the dipole moment for linearly polarized light withthe electric field vector along the z direction, at each timestep and then numerically differentiate it once with re-spect to time to obtain the dipole velocity v ( t ) and againto obtain the dipole acceleration a ( t ). Following Joachain et al. [17] and Telnov et al. [18], the spectral density S ( ω )is obtained as S ( ω ) = 23 πc | ˜ a ( ω ) | = 23 πc ω | ˜ v ( ω ) | , (3)where ˜ a ( ω ) = (cid:90) + ∞−∞ a ( t ) e iωt dt (4)with a similar definition for ˜ v ( ω ). Without any furthermanipulation, we cannot use ω | ˜ z ( ω ) | , since the distri-bution of the ejected electrons may lead to a remainingnet dipole moment at the end of the pulse, which causesproblems when calculating the Fourier transform. Sincethis final dipole moment is constant, its time derivativevanishes, and hence neither the velocity nor the acceler-ation form are affected. We verified that, indeed, eitherof those two forms yield the same result for S ( ω ).Finally, we need to account for the fact that the ac-tual ground state of Ne is (1 s s p ) S . Hence, if weonly look at the 2 p electron, we need to perform calcula-tions for the possible magnetic sublevels, i.e., m = 0 and m = ±
1, respectively. Due to the symmetry propertiesof the process, the results for m = − m = +1 fora linearly polarized driving field are identical. An un-polarized initial state is then simulated by weighting theresults for those sublevels by 1/3 for m = 0 and 2/3 for m = 1, respectively. We will see below, however, that forthe process of interest in this manuscript, the contribu-tion from m = 0 dominates that for m = ±
1. Hence, inpractice it is sufficient to only perform the calculation for m = 0 and divide the result by 3 to obtain the appropri-ate absolute scale. Even though such an absolute scale,representing the single-atom response, is not relevant forexperimental applications, we believe it is important fortheory to put the predictions on an absolute scale usingthe procedure outlined above. B. The R-Matrix with Time-Dependence Method
RMT is an ab-initio technique that solves thetime-dependent Schr¨odinger equation by employing the R -matrix paradigm, which divides the configurationspace into two separate regions over the radial coordi-nate of a single ejected electron [19].In the inner region, close to the nucleus, the time-dependent N -electron wave function is represented byan R -matrix (close-coupling) basis with time-dependentcoefficients. In the outer region, the wave function is ex-pressed in terms of residual-ion states coupled with theradial wave function of the ejected electron on a finite-difference grid. The solutions in the two separate re-gions are then matched directly at their mutual bound-ary. The wave function is propagated in the length gaugeof the electric dipole operator, which, for the typicalatomic structure descriptions used in RMT, convergesmore quickly than the velocity gauge [20].In the RMT calculations presented here, the neon tar-get is described within an R -matrix inner region of ra-dius 20 a , and an outer region of up to 65,000 a . Thefinite-difference grid spacing in the outer region is 0.08 a and the time step for the wave-function propagation is0.01 a.u.. The description is similar to that employed byBurke and Taylor [21], and includes all available 2 s p (cid:15)(cid:96) and 2 s p (cid:15)(cid:96) channels up to a maximum total orbital an-gular momentum of L max = 250. The continuum func-tions are constructed from a set of 50 B -splines of or-der 13 for each angular momentum of the outgoing elec-tron. For the calculations presented in this manuscript,the RMATRX-II suite of codes was used to generate theatomic structure description [22].The harmonic spectra from RMT calculations are ob-tained by Fourier transforming and squaring the time-dependent expectation value of the dipole operator [23],and scaling by the appropriate power of ω . Both thelength and velocity gauges can be applied to obtain the fi-nal spectrum, as the length and velocity matrix elementsgenerated from the time-independent structure calcula-tion can be utilized by the RMT code. Thus we candirectly determine the time-varying expectation valuesof both the dipole operator and the dipole velocity. Un-like in the SAE code, a Gaussian mask is applied to thewave function when determining the expectation value ofthe dipole length, so that the calculation is limited to aregion of almost no ionization. This results in a minimalnet dipole moment and improved agreement between thespectra calculated from either form. C. Calculation Parameters
The primary laser pulse is a 6-cycle, 2000-3000 nmpulse with peak intensities ranging from 1 . × W/cm to 2 . × W/cm . A 3-cycle, sine-squared ramp on/off profile approximates the more re-alistic Gaussian profile sufficiently well, while providingsome numerical advantages over the Gaussian. The shortduration of the pulse was chosen to isolate the generationof the cutoff harmonics to a single trajectory of a singleIR cycle [24], and to mimic a recent experiment thatused HHG from mid-IR sources to generate soft-X-raypulses [25].Calculations were performed on the Stampede2 [26]and Frontera [27] supercomputers, both based at theTexas Advanced Computing Center at the University ofTexas at Austin, and Bridges-2 [28] at the PittsburghSupercomputing Center. Each RMT calculation requiredabout 500 cores for 8-10 hours, with variations depend- ing on intensity and wavelength. Calculations performedwith the SAE code took at most a few hours on a singlenode using all available cores of that node via OpenMPparallelization. III. RESULTS AND DISCUSSION
The first step in determining the efficacy of RMT inproducing water-window harmonics is to ensure that thelength and velocity forms of the produced spectrum agreereasonably well. -20 -19 -18 -17 -16 -15 -14
0 50 100 150 200 250 300 350 400 S ( ω ) [ a . u . ] ω [eV]LengthVelocity2400 nm; 1.0 x 10 W/cm FIG. 1. (Harmonic spectrum of neon in a 2400 nm, 2 . × W/cm laser pulse, as obtained in the length (red) andvelocity (blue) form of the electric dipole operator. The har-monic plateau spans from 35 eV to a cutoff of 372 eV, puttingthe upper limit well into the water window. Figure 1 depicts typical spectra obtained in bothgauges. While the results for individual energies mayvary by 50% or even more, smoothing them (as wouldbe done in practice) suggests excellent agreement for thehigh intensities of the plateau region.Due to the large angular-momentum expansion neededfor the calculation of multi-electron effects, it is impor-tant to verify that the chosen value of the largest term( L max ) does not affect the convergence of the results.This test is performed by running several calculationswhile increasing L max . A fully converged calculation isshown in Fig. 2, where the results obtained with differentvalues of L max are very similar.Having verified that the RMT calculations are bothgauge-invariant and partial-wave converged, we can nowinvestigate the results obtained using the SAE and RMTmethods, as well as the trends when varying the phys-ical parameters the of the calculations. Figure 3 showsa comparison of SAE and RMT results for a wavelengthof 2000 nm and a peak intensity of 2 . × W/cm .We notice good qualitative agreement between the twosets of predictions, although there are substantial differ-ences in the details, especially in light of the logarithmicscale. As mentioned above, the SAE results need to be -23 -22 -21 -20 -19 -18 -17 -16 -15
0 20 40 60 80 100 120 140 160 S ( ω ) [ a . u . ] ω [eV]L max = 200L max = 250L max = 3202000 nm; 1.0 x 10 W/cm FIG. 2. Harmonic spectra for neon in a 2000 nm, 1 . × W/cm laser pulse obtained with different values of L max . The harmonic plateau extends from 15 eV to a cut-off of 145 eV. -21 -20 -19 -18 -17 -16 -15 -14 -13 -12
0 50 100 150 200 250 300 S ( ω ) [ a . u . ] ω [eV] RMT1 x SAE (m=0) / 32 x SAE (m=1) / 32000 nm; 2.0 x 10 W/cm FIG. 3. Harmonic spectra for neon in a 2000 nm laser pulsewith a peak intensity of 2 . × W/cm , as obtained in theRMT and SAE models. For the latter, the results for m = 0and m = 1 are shown with the appropriates weighting factors.See text for details. weighted by the initial magnetic quantum number of the2 p orbital. However, the figure shows that the contri-butions from m = ± m = 0. We also see that both sets of results (RMTand SAE for m = 0) show the plateau reaching up toabout 260 eV, in excellent agreement with the predictionof 258 eV from the three-step model.Since our principal goal is to test the RMT model,we now move to a discussion of the RMT results as afunction of various experimentally controllable parame-ter. Figure 4 shows that the harmonic cutoff changesapproximately linearly with intensity, as expected fromEq. (2). Furthermore, the calculated harmonic cutoffsagree well with those expected from Eq. (1). The largestdisagreement is about 5 eV for the calculation with an -20 -19 -18 -17 -16 -15 -14 -13 -12
0 50 100 150 200 250 S ( ω ) [ a . u . ] ω [eV]1.6 x 10 W/cm W/cm W/cm FIG. 4. Harmonic spectra for neon in a 2000 nm laserpulse with intensities between 1 . × W/cm and 2 . × W/cm . The harmonic cutoffs range from 215 eV to264 eV, following an approximately linear increase with in-tensity. -21 -20 -19 -18 -17 -16 -15 -14 -13
0 50 100 150 200 250 300 350 400 S ( ω ) [ a . u . ] ω [eV]1.6 x 10 W/cm W/cm W/cm FIG. 5. Harmonic spectra for neon in a 2400 nm laserpulse with intensities between 1 . × W/cm and 2 . × W/cm . The harmonic cutoffs range from 295 eV to364 eV, following an approximately linear increase with in-tensity. intensity of 2 . × W/cm .Though the calculations shown in Fig. 4 exhibit thecharacteristics expected from Eqns. (1) and (2), the cut-offs are not yet within the water window. Increasingthe wavelength to 2400 nm, as shown in Fig. 5, pushesthe cutoffs into the water window while maintaining thetrends exhibited in lower wavelength calculations.Another expected trend based on Eq. (2) is an approx-imately quadratic increase in the harmonic cutoff withincreasing wavelength. The cutoffs in Fig. 6, in general,agree with this expectation. However, the disagreementsbetween expected and calculated harmonic cutoffs arelarger than in the calculations in Fig. 4, with disagree-ments up to 6 eV. This is not unexpected, as Eq. (2) does -23 -22 -21 -20 -19 -18 -17
0 50 100 150 200 250 300 S ( ω ) [ a . u . ] ω [eV]2600 nm2800 nm3000 nm1.0 x 10 W/cm FIG. 6. Harmonic spectra for neon in a 1 . × W/cm laserpulse with wavelengths ranging from 2600 nm to 3000 nm. Asexpected, the harmonic cutoffs vary from 222 eV to 289 eV. -22 -21 -20 -19 -18 -17 -16 -15 -14
0 50 100 150 200 250 300 350 400 450 S ( ω ) [ a . u . ] ω [eV]2600 nm2800 nm3000 nm1.5 x 10 W/cm FIG. 7. Harmonic spectra for neon in a 1 . × W/cm laserpulse with wavelengths ranging from 2600 nm to 3000 nm. Asexpected, the harmonic cutoffs vary from 332 eV to 424 eV. not account for the pulse envelope.For practical applications, it is also very important toinvestigate how the harmonic yield depends on the inten-sity and wavelength of the laser. We expect the yield todecrease as the wavelength increases for fixed intensity.This expectation is confirmed in the lowering of the har-monic plateau in Figs. 6 and 7. On the other hand, theHHG intensity increases with increasing laser intensity,as seen in Figs. 4 and 5.A prominent feature in some of the RMT spectra arethe large, localized increases in harmonic yield occur-ring at energies above 250 eV, and most easily visible inFigs. 5 and 7. These features are pseudoresonances, aris-ing from the inclusion of certain ( N +1)-electron boundconfigurations that are included in the close-couplingexpansion to compensate for orthogonality constraintsimposed on the continuum orbitals. These pseudo- resonances could be reduced by a variety of methods:either by improving the underlying atomic structure de-scription in RMATRX-II (although this may result inprohibitively large calculations), or by generating thenecessary atomic structure in RMATRX-I [29], and thenusing a procedure such as those described in [30, 31]to remove the pseudoresonances. A further option is togenerate the required atomic structure using the atomic B -Spline R -Matrix code (BSR) [32], which is currentlybeing interfaced with the RMT code. BSR has theadvantage of using non-orthogonal continuum orbitals,allowing us to avoid the introduction of these ( N +1)-electron terms completely, thereby generally eliminat-ing the appearance of spurious pseudoresonances. Whilethese pseudoresonances are certainly undesirable, theirnet effect is comparatively small, since in practice onewould smooth the predicted spectra. IV. CONCLUSIONS AND OUTLOOK
We have investigated the efficacy of the R -Matrixwith Time-dependence (RMT) method to generate high-order harmonics in the water window using mid-IR laserlight. Previous approaches relied on significant simplifi-cations that neglect both quantum-mechanical and multi-electron effects on the predicted spectra. Since it is notclear how accurate the results of these simplified mod-els are, the present project was designed to carry out athorough test.We successfully demonstrated that RMT is capable ofgenerating harmonics in the water window. We foundqualitative agreement with the trends predicted in clas-sical calculations regarding the cutoff frequency, with thelargest disagreement being a few eV.Comparison with results from SAE calculations showsgood agreement between the two methods in the plateauregion and cutoff regions of the spectrum. At lower en-ergies, the spectra exhibit some significant differences,which could be due to the influence of multi-electron ef-fects. We intend to pursue the source of this discrep-ancy in future calculations. Work is also underway toreduce the effect of pseudoresonances, which manifestthemselves in the spectral region above 250 eV in theRMT calculations. Finally, as an ultimate test of theRMT method, we hope to compare our predictions withexperimental data. ACKNOWLEDGEMENTS
This work was supported by the United States Na-tional Science Foundation under grants No. OAC-1834740 and PHY-1803844. The calculations were car-ried out on Stampede2 and Frontera at the Texas Ad-vanced Computing Center in Austin, with access pro-vided by the XSEDE supercomputer allocation No. PHY-090031 (Stampede2) and Frontera Pathways allocationPHY20028, and on Bridges-2 at the Pittsburgh Super-computing Center through the Early User Program of the allocation No. PHY-180023.The RMT code is part of the UK-AMOR suite and canbe accessed free of charge at [33]. [1] R. Neutze, R. Wouts, D. van der Spoel, E. Weckert, andJ. Hajdu, Nature , 752 (2000).[2] W. Ackermann et al. , Nature Photonics , 336 (2007).[3] Y. Fu, K. Nishimura, R. Shao, A. Suda, K. Midorikawa,P. Lan, and E. J. Takahashi, Communications Physics , 92 (2020).[4] J. Rothhardt, S. H¨adrich, A. Klenke, S. Demmler,A. Hoffmann, T. Gotschall, T. Eidam, M. Krebs,J. Limpert, and A. T¨unnermann, Optics Letters , 5224(2014).[5] G. J. Stein, P. D. Keathley, P. Krogen, H. Liang, J. P.Siqueira, C.-L. Chang, C.-J. Lai, K.-H. Hong, G. M. Lau-rent, and F. X. K¨artner, Journal of Physics B AtomicMolecular Physics , 155601 (2016).[6] J. Li, X. Ren, Y. Yin, K. Zhao, A. Chew, Y. Cheng,E. Cunningham, Y. Wang, S. Hu, Y. Wu, M. Chini, andZ. Chang, Nature Communications , 186 (2017).[7] Y. Fu, K. Nishimura, R. Shao, A. Suda, K. Midorikawa,P. Lan, and E. J. Takahashi, Communications Physics , 92 (2020).[8] P. B. Corkum, Phys. Rev. Lett. , 1994 (1993).[9] A. D. Shiner, C. Trallero-Herrero, N. Kajumba, H.-C.Bandulet, D. Comtois, F. L´egar´e, M. Gigu`ere, J.-C. Ki-effer, P. B. Corkum, and D. M. Villeneuve, Phys. Rev.Lett. , 073902 (2009).[10] M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier,and P. B. Corkum, Phys. Rev. A. , 2117 (1994).[11] I. A. Ivanov and A. S. Kheifets, Phys. Rev. A , 053827(2009).[12] T. Mazza, A. Karamatskou, M. Ilchen, S. Bakhtiarzadeh,A. J. Rafipoor, P. O’Keeffe, T. J. Kelly, N. Walsh, J. T.Costello, M. Meyer, and R. Santra, Nat. Commun. ,6799 (2015).[13] A. C. Brown, G. S. Armstrong, J. Benda, D. D. Clarke,J. Wragg, K. R. Hamilton, Z. Maˇs´ın, J. D. Gorfinkiel,and H. W. van der Hart, Comp. Phys. Commun. ,107062 (2020).[14] O. Hassouneh, A. C. Brown, and H. W. van der Hart, Phys. Rev. A , 043418 (2014).[15] T. Pauly, A. Bondy, K. R. Hamilton, N. Douguet, X.-M.Tong, D. Chetty, and K. Bartschat, Phys. Rev. A ,013116 (2020).[16] X. M. Tong and C. D. Lin, Journal of Physics B: Atomic,Molecular and Optical Physics , 2593 (2005).[17] C. J. Joachain, N. Kylstra, and R. M. Potvliege, Atomsin Intense Laser Fields (Cambridge University Press,2012).[18] D. A. Telnov, K. E. Sosnova, E. Rozenbaum, and S.-I.Chu, Phys. Rev. A , 053406 (2013).[19] L. R. Moore, M. A. Lysaght, L. A. Nikolopoulos, J. S.Parker, H. W. van der Hart, and K. T. Taylor, J. Mod.Optics , 1132 (2011).[20] S. Hutchinson, M. A. Lysaght, and H. W. van der Hart,J. Phys. B. At. Mol. Opt. Phys. , 095603 (2010).[21] P. G. Burke and K. T. Taylor, J. Phys. B. At. Mol. Opt.Phys. , 2620 (1975).[22] https://gitlab.com/Uk-amor/RMT/rmatrixII .[23] A. C. Brown, D. J. Robinson, and H. W. van der Hart,Phys. Rev. A , 053420 (2012).[24] K. R. Hamilton, H. W. van der Hart, and A. C. Brown,Phys. Rev. A , 013408 (2017).[25] T. Gaumnitz, A. Jain, Y. Pertot, M. Huppert, I. Jordan,F. Ardana-Lamas, and H. J. W¨orner, Opt. Express ,27506 (2017).[26] .[27] .[28] .[29] K. A. Berrington, W. B. Eissner, and P. H. Norrington,Computer Physics Communications , 290 (1995).[30] T. W. Gorczyca, F. Robicheaux, M. S. Pindzola, D. C.Griffin, and N. R. Badnell, Phys. Rev. A , 3877 (1995).[31] K. Bartschat, Computer Physics Communications ,168 (1998).[32] O. Zatsarinny, Computer Physics Communications174