Single-active electron calculations of high-order harmonic generation from valence shells in atoms for quantitative comparison with TDDFT calculations

We present a reproducible ab-initio method to produce benchmark tests between time-dependent Schrödinger equation (TDSE) in the single-active-electron approximation (SAE) and time-dependent density functional theory (TDDFT) in the highly nonlinear multiphoton and tunneling regime of strong-field physics. To this end we compare results for high-order harmonic generation from valence shells in atoms using the SAE-TDSE approach and TDDFT calculations. As key to the benchmark comparison we obtain an analytic form of SAE potentials based on density functional theory, which we applied for different atoms and ions. The ionization energies of atomic ground and excited states, as well as the energies of inner shells, for the SAE potentials agree well with experimental data. Using these potentials we find remarkable agreement between the results of the two independent numerical approaches (TDDFT and SAE-TDSE) for the high-order harmonic yields in helium, demonstrating the accuracy of the SAE potentials as well as the predictive power of SAE-TDSE and TDDFT calculations for the nonperturbative and highly nonlinear strong-field process of high harmonic generation in the ultraviolet and visible wavelength regime. Finally, as another application of the SAE potentials, high harmonic spectra from outer and inner valence shells are calculated and it is shown that unphysical artifacts in the SAE-spectra from the individual shells are removed once all the amplitudes are considered.

To study nonlinear effects in laser driven larger atoms, numerical solutions of the time-dependent Schrödinger equation (TDSE) within the single-active electron (SAE) approximation are a widely used method.
In the SAE approximation it is assumed that the laser field only interacts with one-typically the most weaklybound-electron in the atom or molecule, considering all other electrons as frozen spectators to the evolving physics. SAE-TDSE calculations have the advantage of being a straightforward extension from common hydrogen TDSE codes. These calculations are computationally efficient and scale well to parameter regimes (e.g., long wavelength and/or high intensity) that are often infeasible to explore with multielectron codes. In cases where multielectron calculations are possible, single-electron approximations still have value, since SAE-TDSE simulations can serve in initial exploratory parameter sweeps before running more expensive calculations to identify multielectron effects.
Furthermore, ab-initio SAE calculations can provide important benchmark results for more sophisticated multielectron methods in parameter regimes where electron correlation effects are small. In the perturbative interaction regime, such a comparison between results of TDDFT simulations and ab-initio calculations using lowest-order perturbation theory has been done recently and a good agreement has been found [28]. This established the applicability of TDDFT calculations for perturbative multiphoton processes. We are not aware of a similar study at higher intensities for highly nonlinear and nonperturbative processes. A few rare studies reporting a comparison between TDDFT and SAE-TDSE results for nonperturbative strong-field processes typically attribute differences to the impact of multielectron effects [27,29] without establishing a benchmark test between SAE-TDSE and TDDFT calculations. The absence of such benchmark tests partially contributes to the misconception in parts of the strong-field community that TDDFT studies cannot provide quantitative predictions for highly nonlinear processes involving the continuum. In this paper we show by carefully constructing a SAE potential that indeed agreement between SAE-TDSE benchmark results and TDDFT results can be found. Hence, we establish a method to check the accuracy of TDDFT calculations in view of multielectron effects and disprove the misconception about the quantitative predictive power of TDDFT calculations.
In section 2 we present the method to obtain analytical forms of single-active-electron potentials based on DFT calculations (which is essential for the benchmark tests). We then briefly discuss the approaches to perform the time-dependent calculations for high-order harmonic generation. In section 3 we proceed by establishing the accuracy of the SAE potentials by presenting predictions for energies of the ground state and singly excited states of various atoms within the single-active-electron approximation (section 3.1). On this basis, we then compare results for high harmonic generation from helium within the single-active-electron approximation with those from TDDFT calculations based on the same DFT approach for a few highly nonlinear cases in 3.2. Finally, we show as another application of the SAE potentials results for high-order harmonic generation from outer and inner valence shells (section 3.3). We summarize our results in section 4. Hartree atomic units (e=ÿ=m=1) are used throughout unless stated otherwise.

Theory
In this section, we discuss the construction of SAE potentials for atomic systems, based on fitting to the effective Kohn-Sham potential. First, we discuss the concept in constructing the SAE potential for the benchmark tests before we briefly reiterate the density functional theory formalism used to obtain the Kohn-Sham potential. Then we present the approach used to approximate the Kohn-Sham potential via analytical functions and discuss strategies to determine the coefficients in the model potential. Finally, we present the methods used to perform the ab-initio SAE-TDSE and TDDFT calculations.

Methodology
The construction of a potential to model single-active-electron behavior has been pursued through several different strategies (e.g., [30][31][32][33][34][35][36][37]). SAE models are frequently developed using analytic functions with free parameters that are fit to accurately simulate experimental behavior in the scenario of interest; often by matching the ionization potential for the single active valence electron to the experimental value since many strong field processes-such as tunneling ionization-strongly depend on it. Besides the long-range behavior of the single-active-electron potential-which is essential for the determination of the ionization behaviorelectron dynamics driven by a strong laser pulse also probe the short-range properties of the potential (e.g., excited states). Constructing a SAE potential, that reproduces well ground, inner and excited state energies, is therefore useful to study laser driven processes from the valence as well as inner shells. Furthermore, it provides a reliable ground for benchmarking TDDFT calculations against ab-inito SAE-TDSE results, as we will demonstrate below. Conceptually, the SAE potential should be therefore based on the same DFT approach that is used in the TDDFT calculations without any free parameter adjustment.
To this end, we have obtained analytical forms of single-active-electron potentials, using density functional theory (DFT) and physical reasoning regarding the electronic structure of an atomic system. Fits to the effective Kohn-Sham potential are made without using ad hoc free parameters. We establish the accuracy of the potentials by comparing ionization and excitation energies. Basing the approaches on the same functional we are set up to make an accurate quantitative comparison between SAE-TDSE and TDDFT calculations for high harmonic generation (HHG) from helium at intensities in the nonperturbative regime. We further note that by providing analytical forms of the SAE potentials the SAE-TDSE versus TDDFT comparisons can be reproduced or extended in the future.

Analytical single-active-electron potential fits: Modeling atomic species
In the Kohn-Sham density functional theory (DFT) formulation of multielectron systems [38,39], each orbital ( ) y s r i KS , (i=1, K N) involved in the ground state of the spin-polarized N-electron system is determined by solving a set of Schrödinger-like eigenvalue equations given by is the effective Kohn-Sham potential and σ denotes the spin of the electron. The terms in equation (2) represent the external Coulomb potential . An exact expression for V xc is not known and several popular approximations to this quantity exist. We make use of the strategy proposed by Krieger, Li, and Iafrate [40] and extended in [26]. We choose a sum of local (spin-)density approximation (LDA) exchange and correlation functionals: , , as defined in appendix C of [41], and apply a self-interaction correction [41]: This results in an orbital-dependent potential, but the Optimized Effective Potential procedure from [40] resolves this via where s X i, is the solution to the set of linear equations: where m is the index of the highest occupied orbital (or indices for a degenerate set), so those terms are ignored. The potential ( ) has been found self-consistently by guessing a set of orbitals, i.e. r s i, , then using the orbitals to calculate an exchange correlation potential (equation (5)), which is then used to update the orbital set as the lowest energy eigenstates of the effective Kohn-Sham potential (equation (2)). The DFT calculations have been completed on a 1D spherical grid with spacing dr=0.0025 a.u. and length L=80 a.u., requiring that the valence orbital energy s E m KS , converges to a difference below 10 −10 a.u. between iterations. One can proceed by using the DFT potential numerically; it is however instructive and useful to provide an analytical form of the single-active-electron potential as: where V long is the long-range Coulomb term: In test calculations we have left C 0 and Z c as free parameters in the fitting procedure. The results are close to the expected values of ( , where Z is the charge of the nucleus and N e the number of electrons in the atom or ion, and = -Z Z C c 0 for the remainder of the charge. For this reason we did not vary these parameters in the remainder of our studies.
V shell captures the effective SAE potential at intermediate distances with one exponential term for each of the n orbital shells: This is motivated by the comparison of electron density and exchange correlation potential based on DFT calculations (see figures 1 and 2) where we identify one smooth region per shell. We have applied the above potential form to both atoms and ions up to argon. For a given species the same potential is found for both valence and inner orbitals. We note that for larger atoms the present approach based on an orbital shell picture may need to be modified, since in such atoms subshell densities are more strongly overlapping. The effective single-active-electron potential V SAE depends on n linear parameters (a i ) and n+1 nonlinear parameters using the Levenberg-Marquardt (LM) least-squares curve fitting algorithm [42,43] implemented in the LMFIT library [44]. The LM algorithm involves methods to avoid local minima, resulting in nonlocal fitting based on initial parameter guesses. As an extra check, the parameters b i in the exponentials have been scanned across a range of initial guesses. In table 1 we show the results for the  parameters of our fits for a number of atomic and ionic species. The relative error between the fit and the actual DFT results for  V eff , is presented in figure 3 for helium and argon.

High-order harmonic calculations
To obtain HHG spectra in the SAE-TDSE calculations, we have propagated the wavefunction using the static potential ( ) V r SAE and the time-dependent laser interaction potential described in length gauge as: where E is the electric field of the laser. In the present calculations, the initial state was set as the ground state of the helium SAE potential. We have discretized the wavefunction and (radial) potential on a cylindrical 2D grid with dρ=dz=0.03 a.u, and grid sizes up to r = 50 max a.u. and =  z 100 max a.u.. Using the Caylee form of the propagation exponential, we have applied the 2nd order Crank-Nicolson method and the split-operator approximation to propagate the wavefunction starting from the initial state. As absorbing boundary, we have used exterior complex scaling (ECS), where the edges of the grid are rotated into complex space by an angle η=π/4.
To study the accuracy of the SAE potentials we compare results of SAE-TDSE calculations with those of TDDFT calculations. For the latter we have propagated the wavefunction using the Kohn-Sham orbitals obtained by solving the time-dependent Kohn-Sham equations    where 1 7 ext I and the Hartree and exchange correlation terms depend now on the time-dependent density ρ(r, t). The TDDFT calculations have been performed using the OCTOPUS code [45][46][47] on a Cartesian 3D grid with spacing dx=0.3 a.u. in a cylindrical box of length at least 70 a.u. in the z-direction, and diameter at least 40 a.u. in the xy-plane. An additional 20 a.u. in every direction outside the box was added for the complex absorbing potential. Note that use of 4th-order finite difference allows for coarser dx here.
We have taken the laser to be linearly polarized in the z-direction, with a sin 2 pulse envelope, variable central frequency ω 0 , and full-width duration T=8π/ω 0 , corresponding to four optical cycles. That is, for 0tT. We have evaluated the dipole acceleration a(t) using the Ehrenfest theorem where V=V SAE (for SAE-TDSE) and = -V Z r (for TDDFT). In the case of TDDFT, the dipole accelerations from all orbitals are added together, and other interaction terms cancel. The harmonic spectrum P(ω) is then obtained by taking the windowed Fourier transform of the dipole acceleration, 0.08 cos 4 is the Blackman window.

Results and discussion
In this section we first test the single-active-electron potential fits via comparison of ionization and excitation energies for different atomic species obtained within the approach as well as comparison of SAE-TDSE and TDDFT high-order harmonic generation calculations for helium atom. We then proceed by studying the impact of the different terms in the SAE potentials on the HHG spectra and finally consider HHG spectra from outer and inner valence shells.

Ionization and excitation energies
We have obtained energies for the ground state and (singly) excited states of the atoms and ions for the analytical SAE potential fits using the Implicitly Restarted Arnoldi Method [48]. For the analytical SAE potentials we have performed calculations on the same 1D spherical grid as for the DFT calculations with dr=0.0025 a.u. and L=80 a.u.. In table 2 we compare the ground state energies for the atoms and ions studied, obtained in the DFT calculations and with the analytical SAE potential fits, with the experimental data from the NIST Atomic Spectra Database [49]. The SAE potential results reproduce the DFT results remarkably well, showing the success of the fitting procedure. The same conclusion holds for the ionization energies of inner orbitals, as the example for neon atom in table 3 shows. Further improvement concerning the agreement between the theoretical results and experimental data requires more advanced DFT calculations. In tables 4 and 5 we present a comparison between experimental data [49] and the predictions obtained using the SAE potentials for some singly excited state energies of the helium and argon atoms, which shows a reasonable agreement.

Comparison between SAE-TDSE and TDDFT calculations for HHG spectra
The good agreement between analytical fit and DFT results provides the grounds for the key goal of our study, namely to proceed with the quantitative comparison between TDDFT and SAE-TDSE calculations for HHG. In figure 4, we compare results of SAE-TDSE and TDDFT calculations for the HHG spectrum driven by a 4-cycle laser pulse at central wavelengths of (a) l = 267 nm (ω 267 =4.65 eV), (b) λ=400 nm (ω 400 =3.10 eV), (c) λ=600 nm (ω 600 =2.07 eV) and peak intensity of = I 10 peak 15 W cm −2 . For the comparison we have matched the two results at one point, namely the energy of the driving laser pulse. As a target we used the helium atom, this being the simplest atomic system with electron correlation effects. For laser wavelengths in the ultraviolet and optical regime we expect that the correlation effects on the nonlinear and nonperturbative dynamics may be small. Thus, the single-active-electron approximation should be applicable and, hence, constitute a test bed for the comparison.
At the shortest wavelength (panel (a)), it requires six photons to exceed the field-free ionization potential I p =25.7 eV (c.f., table 2), while the energy of seven photons is needed to exceed the field-shifted ionization threshold of where =´-U 9.3 10 is the ponderomotive energy of a free electron in a laser field at the peak intensity. Furthermore, the HHG cut-off energy of + = U I 3. 17 46.7 eV p p w »10 267 . The results of the TDDFT and SAE-TDSE calculations are in excellent agreement over the whole spectrum, up to the expected cut-off and even beyond. The position of the maxima and minima, the form of the odd high-order harmonic lines as well as the cut-off in the spectrum are nearly identical in the results of both independent calculations. We note small deviations between the results in the region of the 3rd and 5th harmonics, which may be either related to differences in the excited state spectra for the two calculations or indicate minor electron correlation effects. However, overall the comparison shows that electron correlation effects appear to be negligible for this nonlinear process in helium at a UV driver laser wavelength. Similarly good agreement is found at the longer central wavelengths of 400 nm (b) and 600 nm (c), The other laser parameters are held the same as for the calculations at 267 nm. We note that while the 267 nm calculation was in the so-called multiphoton regime, the Keldysh parameter γ=0.62 [50] indicates that the process at 600 nm occurs in the tunneling interaction regime. Due to the longer wavelength the ponderomotive energies are significantly larger. This results in larger field-shifted ionization potentials and corresponding HHG cutoffs of 73.1 eV (b) and 132.2 eV (c), corresponding to the 24th and 64th harmonic of the driver frequency, respectively. While the differences between the results of the two independent calculations are slighter larger than at the UV driver wavelength, they remain overall clearly less than a factor of two. The largest deviations are again found in the energy region near the ionization threshold. We note that at higher energies even subtle features are in excellent agreement for the two calculations.
The agreement between the results demonstrates the accuracy and applicability of the single-atom potentials and shows the capability of TDDFT to provide reliable results for nonperturbative strong-field processes in this parameter regime. The deviations in the ionization threshold region may be attributed to a minor role of electron correlation and/or excited states during the process. We note that studies reporting a comparison between TDDFT and SAE-TDSE results for nonperturbative strong-field processes are rare [27,29] and we are not aware that a similar benchmark test between the two types of numerical calculations, which are based on independent computer codes, in the strong-field and tunneling regime has been performed and an agreement in parameter regimes of small correlation effects has been achieved. Typically, differences between results have been interpreted as the effect of correlation effects [27], even in parameter regimes in which we find excellent agreement between SAE-TDSE and TDDFT results. This demonstrates the importance of the method in constructing SAE potentials used in the present study and the general benchmark tests performed. It also clarifies the sometimes expressed misconception that TDDFT calculation cannot provide accurate results in the highly nonlinear multiphoton and tunneling regime.

HHG from outer and inner valence shells
Beyond the application for benchmark-testing TDDFT calculations, the present SAE potential can of course be applied for other studies as well. Recently, the influence of inner valence shells on high harmonic generation has been studied [51][52][53][54]. The effects in the spectra are often interpreted by multielectron or intershell correlation, so accurate single-active-electron calculations can serve as a useful tool for comparison and, hence, identification of the correlation effects.
Using the potentials for neon we have obtained high harmonic spectra from the outermost 2p shell (m=0) and and the inner 2s valence shell. In the calculations we have applied a laser pulse of 4 optical cycles at 800 nm and a peak intensity of 10 14 W cm −2 . Both spectra (c.f., figure 5(a) and (b)) show the characteristic features of a HHG spectrum with plateau and cut-off. However, in both spectra there occurs a prominent peak between the 15th and 17th harmonic, at the energy corresponding to the energy difference between the 2p and 2s states in our SAE potentials for neon. The peak is an artefact of the SAE-TDSE calculations which allow for transitions between the sub-shells, although they are filled in the real multielectron atom. The unphysical feature is however completely removed once the amplitudes for the processes driven from the two sub-shells are added and the full Figure 5. Results of SAE calculations for high-order harmonic spectra generated from (a) 2p-shell and (b) 2s-shell in neon by a 4-cycle laser pulse at peak intensity of 10 14 W cm −2 and wavelength of 800 nm. The vertical dashed line indicates thes p 2 2 transition energy at which an artificial peak occurs in the spectra from the two sub-shells, which is removed once the the amplitudes for both processes are added (panel (c)).

Summary
To summarize, we have presented an analytical form for atomic SAE potentials, which are based on results of DFT calculations. The potential fits are shown to provide a good agreement with the DFT results for the energy of the ground state as well as with experimental data for the energies of ground, excited states, and inner shells, of different atoms and ions, without any ad hoc adjustment of parameters in the analytical forms Using the potentials, a benchmark test between between SAE-TDSE and TDDFT calculations has been performed and a remarkable agreement throughout the whole HHG spectra at ultraviolet and visible wavelengths has been achieved. Since the calculations are based on two different approaches and different numerical codes, the agreement indicates the accuracy of the SAE potentials as well as the capability of TDDFT to provide reliable results for a nonperturbative, highly nonlinear process such as high harmonic generation in this parameter regime. Finally, results of SAE calculations for HHG spectra from different shells in neon reveal an artificial large peak at the transition energy between the 2s and 2p-states, which is removed once the amplitudes of the processes from the two sub-shells are added up.