Effects of inner electrons on atomic strong-field ionization dynamics
aa r X i v : . [ phy s i c s . a t o m - ph ] M a r Effects of inner electrons on atomic strong-field ionization dynamics
J. Rapp and D. Bauer
Institut für Physik, Universität Rostock, 18051 Rostock, Germany (Dated: August 14, 2018)The influence of inner electrons on the ionization dynamics in strong laser fields is investigated in a wave-length regime where the inner electron dynamics is usually assumed to be negligible. The role of inner electronsis of particular interest for the application of frozen-core approximations and pseudopotentials in time-dependentdensity functional theory (TDDFT) and the single-active-electron (SAE) approximation in strong-field laserphysics. Results of TDDFT and SAE calculations are compared with exact ones obtained by the numerical abinitio solution of the three-electron time-dependent Schrödinger equation for a lithium model atom. It is foundthat dynamical anti-screening, i.e., a particular form of dynamical core polarization, may substantially alter theionization rate in the single-photon regime. Requirements for the validity of the approximations in the singleand multiphoton ionization domain are identified.
PACS numbers: 32.80.Fb , 32.80.Rm , 31.15.ee , 31.15.es
I. INTRODUCTION
Density functional theory (DFT) simulations (see, e.g., [1–3]) have become popular tools for electronic structure calcu-lations. Compared to the exact solution of the many-bodySchrödinger equation, discrepancies in the electron densityobtained from DFT-based Kohn-Sham (KS) [4] schemes are,by construction, caused by the unavoidable approximationto the generally unknown exact exchange-correlation (XC)potential. However, in practice it is common to apply ad-ditional approximations, most notably “pseudopotentials” orthe “frozen-core approximation” (see, e.g., [5, 6]) in order toavoid the numerical effort of treating tightly-bound inner elec-trons. The justification for this neglect is that core electronsdo not take directly part in, e.g., the formation of bonds.DFT has been extended to systems in time-dependent ex-ternal potentials. In principle, time-dependent density func-tional theory (TDDFT) (see, e.g., [7, 8]) allows to study many-electron systems such as atoms, molecules, or clusters instrong laser fields, even beyond linear response. It is known,however, that switching from DFT to TDDFT makes the un-known exact XC functional even more inaccessible becauseof memory effects and the initial-state dependence it shouldcontain [8–12] but all practicable approximations to it do not.In this paper, we study inner-electron dynamics induced bytime-dependent external fields. In contrast to typical applica-tions of DFT concerning the ground state electronic structureof the system at hand, even the lowest KS orbitals may un-dergo a significant modification if the system is subjected to astrong external laser field in TDDFT beyond linear response.Evidently, it is invalid to freeze those KS orbitals which di-rectly contribute to, e.g., the outgoing electron density of anatom being ionized. If, on the other hand, solely the KS va-lence orbital dominates ionization, it is an eligible questionif the essential dynamics can be reproduced by a frozen-core,pseudo, or single-active-electron (SAE) potential. In fact, theSAE approximation is ubiquitous in the strong-field ionizationcommunity (see, e.g., [13–16]). Only recently it has been rec-ognized that in multi-electron molecules it is often not permis-sible to consider only the highest occupied molecular orbitalin strong-field processes [17, 18]. The question we address in this paper is whether core elec-trons in atoms can be considered “frozen” or not in the inter-action with long wavelength radiation. Here, “long” meansthat the photon energy ~ ω should be small compared to theenergy by which the core electrons are bound. Given the en-ergy level spacings of the Li atom, we thus need to considerthe multiphoton regime, and the single-photon regime up tophoton energies well below values where the core electronsare accessed “directly” by the applied laser field.We employ a model Li atom, for which we are able tonumerically solve the time-dependent Schrödinger equation(TDSE) ab initio . Lithium is the simplest element with coreelectrons in the ground state configuration and thus serves as aperfect testing ground. However, strong-laser driven Li in fulldimensionality is well beyond nowadays computational capa-bilities. Even with the spatial degrees of freedom restrictedto one dimension (i.e., the laser polarization direction) perelectron, the numerical demand is enormous for strong laserfields. We have been able to speed up the calculations by em-ploying properties of the time-dependent, spatial three-bodywavefunction in the ionization regime considered and by op-timizing the TDSE solver for graphics processing units.The paper is organized as follows. In Sec. II the Li modelsystem is described. In Sec. III the methods and approxima-tions used in this work are introduced. Results are presentedin Sec. IV, a conclusion and outlook are given in Sec. V. Re-marks on numerical details are attached as an Appendix.Atomic units (a.u.) are used throughout. II. ONE-DIMENSIONAL LITHIUM MODEL
The Li atom is the simplest atom with “inner” and “outer”electrons in the ground state configuration. The reduction toone dimension per electron is required for the exact numeri-cal treatment, as the computational effort grows exponentiallywith both particle number and dimension. One-dimensionalatom models have been successfully used in the case of he-lium for more than 20 years [19] and more recently for Li aswell [20].Applying the dipole approximation, the Hamiltonian inlength gauge reads H ( t ) = X i T ( i ) + V ( i ) + H ( i )L ( t ) + 12 X j = i W ( ij ) (1)with indices i, j ∈ { , , } , the kinetic energy operator T ( i ) = 12 (cid:16) p ( i ) (cid:17) , (2)the core potential V ( i ) = − Z (cid:20)(cid:16) x ( i ) (cid:17) + ε (cid:21) − / , Z = 3 , (3)the coupling to the laser field H ( i )L ( t ) = E ( t ) x ( i ) , E ( t ) = − ˆ E sin ωt, (4)and the electron-electron interaction operator W ( ij ) W ( ij ) = (cid:20)(cid:16) x ( i ) − x ( j ) (cid:17) + ε (cid:21) − / . (5)The smoothing parameter ε = 0 . is tuned such that thetotal energy of the “real,” three-dimensional Li atom is repro-duced. III. METHODS AND APPROXIMATIONS
In this Section, we introduce the three methods used(TDSE, Floquet, and TDDFT) and the various approximations(frozen-core, pseudopotentials, and SAE). Particular empha-sis is put on the structure of the three-electron wavefunction,which can be decomposed into a sum of three terms, each fac-torizing in a spin and a spatial part.
A. Time-dependent Schrödinger equation
The TDSE i ∂ t | Ψ ( t ) i = H ( t ) | Ψ ( t ) i (6)is the fundamental equation describing the non-relativistictime-evolution of a many-particle quantum state | Ψ ( t ) i . Dueto the unavailability of an analytical solution for the Hamil-tonian (1) we solve the TDSE numerically. In that way exactbenchmark results are obtained to which results from approx-imate approaches will be compared.
1. Three-electron state | Ψ ( t ) i Let us expand the state | Ψ ( t ) i in orthonormal single-particle states {| n i} n ∈ N , | n i = | φ n i ⊗ | χ n i , (7) where | φ n i and | χ n i are spatial and spin components, respec-tively. Suppressing the time-argument, the expansion for threeparticles reads | Ψ i = X n " X k a kn | k i (1) ! ⊗ X l b ln | l i (2) ! ⊗ X m c mn | m i (3) ! (8) = | ABC i . (9)For brevity, the ⊗ -sign denoting the tensor product will beomitted from now on. The shorthand notation (9) allows toconcisely formulate the correct exchange antisymmetry in thecase of fermions | ABC i = −| ACB i = | CAB i = −| BAC i = | BCA i = −| CBA i . (10)The antisymmetry (10) can be enforced on a general three-particle state | A ′ B ′ C ′ i by the antisymmetrization operator A , | Ψ i = | ABC i = N ′ A| A ′ B ′ C ′ i , (11)where the normalization factor N ′ has to be chosen such that h Ψ | Ψ i = 1 , and A| A ′ B ′ C ′ i = 13! (cid:0) | A ′ B ′ C ′ i − | A ′ C ′ B ′ i + | C ′ A ′ B ′ i−| B ′ A ′ C ′ i + | B ′ C ′ A ′ i − | C ′ B ′ A ′ i (cid:1) . Introducing the abbreviation | xσ i = | x i| σ i| x i| σ i| x i| σ i , the expansion of | Ψ i in position-spin space reads | Ψ i = X σ σ σ Z Z Z d x d x d x | xσ ih xσ |N ′ A| A ′ B ′ C ′ i . (12)The configuration is chosen such that the total spin is S = 1 / and M S = +1 / at all times. This can be justified by the factthat interaction Hamiltonians that could induce spin-flips arenot considered in this paper. The primed state | A ′ B ′ C ′ i canbe chosen to have separable spin components, e.g., the corre-sponding expansion coefficients a ′ kn are only non-vanishingfor spin-down while the other two coefficients always resultin spin up. As a consequence, the function Ψ ( x , x , x ) = X σ σ σ | σ ih xσ |N ′ A| A ′ B ′ C ′ i (13)can be written as Ψ ( x , x , x ) = N ′ h |↓↑↑i (cid:0) φ ( x , x , x ) − φ ( x , x , x ) (cid:1) + |↑↓↑i (cid:0) φ ( x , x , x ) − φ ( x , x , x ) (cid:1) + |↑↑↓i (cid:0) φ ( x , x , x ) − φ ( x , x , x ) (cid:1)i , (14)with correlated spatial functions φ ( x , x , x ) = P klmn a ′ kn b ′ ln c ′ mn h x | φ k ih x | φ l ih x | φ m i . Defining φ ( x , x , x ) = N ′′ (cid:2) φ ( x , x , x ) − φ ( x , x , x ) (cid:3) , (15)which is antisymmetric with respect to the exchange of its sec-ond and third argument, one obtains the compact form Ψ ( x , x , x ) = N h |↓↑↑i φ ( x , x , x )+ |↑↓↑i φ ( x , x , x )+ |↑↑↓i φ ( x , x , x ) i , (16)where N = N ′ N ′′ . The full three-electron state | Ψ i is—at alltimes—completely determined by φ , | Ψ i = N h P (12) P (23) + P (23) P (12) i |↓↑↑i | φ i . (17)Here, P ( ij ) is the two-particle permutation operator which ex-changes the indices of particles i and j , and | φ i = Z Z Z d x d x d x | x i φ ( x , x , x ) . (18)We assume that | φ i is normalized to unity, h φ | φ i = 1 ,so that N = √ .
2. Spatial TDSE
Inserting a time-dependent state | Ψ ( t ) i of the form (17)into the TDSE (6) yields h P (12) P (23) + P (23) P (12) ih i ∂ t − H ( t ) i |↓↑↑i | φ ( t ) i = 0 (19)because both ∂ t and H ( t ) commute with any two-particle per-mutation operator P ( ij ) .Assuming a spin-diagonal Hamiltonian, multiplication of(19) from the left by, e.g., h↓↑↑| , yields a TDSE for the evolu-tion of | φ ( t ) i in time, i ∂ t | φ ( t ) i = h↓↑↑| H ( t ) |↓↑↑i | φ ( t ) i . (20)Although the Hamiltonian (1) does not act on spin compo-nents at all, the TDSE (20) still holds for spin-diagonal Hamil-tonians. This will be utilized in Sec. IV C. The TDSE (20) isthe one that is actually solved numerically in position space ona discretized x x x -grid. More details about the numericalsolution are described in the Appendix.
3. Observables for both spin projections
Although electrons are indistinguishable, the partial wave-function φ allows to extract information about inner andouter electrons separately. Given a one-particle operator a ( i ) acting on the spatial component only, one can construct a spin-spatial operator a σ of the form a σ = X i | σ i ( i ) h σ | ( i ) a ( i ) , (21)where | σ i (1) h σ | (1) = X σ σ | σ i | σ i | σ i h σ | h σ | h σ | , | σ i (2) h σ | (2) = X σ σ | σ i | σ i | σ i h σ | h σ | h σ | , | σ i (3) h σ | (3) = X σ σ | σ i | σ i | σ i h σ | h σ | h σ | . The choice of either | σ i = |↓i or | σ i = |↑i then yields an op-erator for calculating observables for the single spin-down in-ner electron on the one hand and for the two spin-up electronson the other hand, respectively. The latter will be referred toas inner-outer spin component [21].
4. Ionization x x x FIG. 1. (Color online) Schematic view of a cubic simulation boxaround the nucleus in the center. Different colors indicate those re-gions where zero (neutral Li), one (Li + ), two (Li ), or all threeelectrons (Li ) are located at positions far from the nucleus, re-spectively. The ionization probability is chosen as the primary observ-able for our investigations because it is well-defined and com-parable among all considered methods. In the TDSE simula-tion, position space ( x , x , x ) is divided into four regions,differing by the number of electrons which are located faraway from the nucleus (see Fig. 1). The respective ioniza-tion regions are (i) no ionization: single small cube aroundthe nucleus, (ii) single ionization: six channels pointing to thesurface centers of the simulation box, (iii) double ionization:twelve cuboids, four lying in each of the three central plains,(iv) triple ionization: eight cubes in the corners of the simula-tion box.The laser parameters considered throughout this paper aresuch that multiple ionization is negligible. Thus, the (single-)ionization probability p ( t ) reduces to p ( t ) = 1 − N ( t ) where N ( t ) is the norm inside the cube around the nucleus (repre-senting neutral Li). B. Floquet method
The whole purpose of applying the Floquet approach in thiswork is the determination of resonances, taking the AC Starkeffect into account.In general, the Floquet method (see, e.g., [22, 23]) allowsto study time-dependent problems without an explicit time-propagation. This is possible if the Hamiltonian H ( t ) is peri-odic in time, i.e. H ( t ) = H + H L ( t ) , H L ( t + T ) = H L ( t ) , (22)because the so-called Floquet theorem then allows to obtaina time-independent eigenvalue equation for the field-dressedstates and the corresponding quasienergies. Quasienergyspectra are useful for predicting resonance enhancements inthe ionization rate as a function of the laser frequency with theAC Stark effect automatically included. The Floquet methodwill thus be used to follow the quasienergies ǫ ( n ) m for vary-ing photon energies ω and a fixed electric field amplitude ˆ E .Here, the index m ∈ N refers to the (unperturbed) atomic en-ergy level, the index n ∈ Z to the “Floquet block” (note that ǫ ( n ) m + ω = ǫ ( n +1) m ). C. Time-dependent density functional theory
DFT [1–3] and its time-dependent extension TDDFT [7, 8]are approaches to overcome the exponential scaling of the nu-merical effort with the number of particles. They are basedon the fact that all information about the system is included inthe (time-dependent) single-particle density n ( x ) , as provenby the Hohenberg-Kohn theorem [24] and its time-dependentanalog, the Runge-Gross theorem [25]. (TD)DFT calcula-tions are performed in practice within a KS scheme, i.e., thesingle-particle density is reproduced with the help of a ficti-tious, much easier-to-solve noninteracting system evolving inan effective KS potential.In this work, the spin-polarized Li-system is studied. Wetherefore consider the spin-densities n σ ( x ) , σ ∈ {↓ , ↑} .The exchange-only local-spin-density approximation (LDA)is employed for the XC functional. The exchange functionalfor the three-dimensional electron gas is considered becausethe one-dimensional model introduced in Sec. II is meant tomimic a three-dimensional three-electron atom in a linearlypolarized laser field rather than a true one-dimensional sys-tem. The time-dependent KS equation reads (spatial argumentssuppressed) i ∂ t ϕ i ( t ) = H ( σ i )KS ( t ) ϕ i ( t ) (23)with the KS Hamiltonian H ( σ i )KS ( t ) = − ∂ ∂x + v ( t ) (24) + v (H) [ n ( t )] + v (XC) [ n σ i ( t )] and n σ ( t ) = P i n i ( t ) δ σ i σ , n i ( t ) = | ϕ i ( t ) | , n ( t ) = P σ n σ ( t ) . The same external potential v ( t ) as in the many-body TDSE (i.e., binding potential plus laser in our case) ap-pears here, v (H) is the Hartree potential, and v (XC) is the XCpotential (to be approximated).A known problem of LDA is the wrong asymptotic behav-ior of the KS potential v + v (H) + v (XC) . Each KS orbital in an,e.g., unperturbed, neutral atom should experience a potential − /r far away from the nucleus, representing one unscreenednuclear charge. This correct behavior can be enforced bythe so-called Perdew-Zunger (PZ) self-interaction correction(SIC) [26]. PZ SIC corrects the Hartree and XC potentials foreach orbital i by subtracting the potentials evaluated for the“own” single-particle density n i . In the case of PZ-correctedLDA, the final corrected XC potential is calculated as (timeand space arguments suppressed) v ( LDA+PZ ) i = v (LDA) [ n σ i ] − v (H) [ n i ] − v (LDA) [ n i ] . (25)It is easy to see that this leads to the correct SIC in the single-electron limit. However, due to the nonlinearity of the KSpotential, i.e., v (XC) [ n ] − v (XC) [ n i ] = v (XC) [ n − n i ] , the self-interaction is not completely removed by PZ-SIC in general.The PZ-corrected KS Hamiltonian is orbital-dependent, i.e., itis in general different for each orbital i (i.e., not only differentfor orbitals with different σ i , as in “conventional” spin-DFT).An unpleasant consequence of this orbital-dependence of theKS Hamiltonian is the non-orthogonality of the PZ-SIC KSorbitals (although in practice they are usually very close toorthogonal). Positive consequences of the PZ-SIC are that,besides the correct asymptotic behavior of potentials and den-sities, the values for the total energy, and ionization energies(or electron affinities of negative ions) typically improve. D. Frozen-core approximation
All-electron (TD)DFT calculations often become numeri-cally too demanding. Hence, it is common to apply additionalapproximations in order to reduce the numerical effort further.One may employ the fact that chemical bonds and reactionsare governed by valence electrons so that the relaxation ofelectronic core shells may be considered negligible. Conse-quently, tightly localized KS orbitals may be self-consistentlydetermined for the initial state configuration but “frozen” dur-ing the actual TDDFT time-propagation. This approach willbe referred to as “frozen-core approximation” (FCA). One ofthe issues in this work is the validity of the FCA for atoms in-teracting with a strong laser field, i.e., in TDDFT calculationsbeyond linear response.
E. Pseudopotentials
Consider a one-dimensional (e.g., radial) KS orbital, whichis orthogonal to n other mutually orthogonal orbitals. If theseother orbitals are constructed with the minimum number ofnodes the considered orbital must have at least n nodes. Anumerical grid thus requires a fine spatial resolution in orderto resolve all orbital nodes, including the behavior in-betweenwhere the second derivative can reach high absolute values.A popular tool which aims at circumventing the numeri-cal demand caused by orbitals with many, densely-distributednodes are “pseudopotentials” (see, e.g., [5, 6]). After applyingthe FCA, an artificial potential is constructed such that thoseorbitals which are not frozen yield the same single-particledensity outside a certain cutoff radius r c as in the full calcu-lation but have less nodes within [0 , r c ] . Moreover, the pseu-dopotential is “designed” to reproduce the KS energies of theunfrozen orbitals (and possibly also those of the unpopulated,excited states). As in the motivation of the FCA, the argumentfor using pseudopotentials in chemistry is that only valenceelectron densities are important for the questions of interestsuch as molecular binding properties and chemical reactions.In other words, only the electron density far away from thenuclei is important.In the hierarchy of approximations, pseudopotentials residebelow the FCA. Hence we do not test particular pseudopoten-tials in this work. If the FCA fails, pseudopotentials will failtoo (unless there is a lucky cancellation of errors caused bythe removal of the nodes for r < r c ). F. Single-active-electron approximation
In the single-active-electron (SAE) picture it is assumedthat the electron under investigation can be described as a sin-gle particle moving in an effective, external potential. It thusmay be also viewed as an FCA. DFT provides one option toapproximate this effective, external potential: one performsan all-electron DFT calculation for the desired initial elec-tron configuration (usually the ground state) and subsequently“freezes” all KS orbitals but the one for the SAE of interestfor the real-time propagation [27]. In that way, a dependenceon the XC potential is introduced only indirectly through theinitial state. In the strong-field community, often simple ana-lytical expressions for screened Coulomb potentials with ad-justable screening parameters are used [28, 29].The question arises which state the SAE should populatein the effective potential. Common choices are (i) the groundstate (corresponding to a pseudopotential approach with onevalence electron) or (ii) some excited state (corresponding topure FCA). The latter choice is more relevant for applicationsto intense laser-atom interaction, as the orbital symmetry ofthe initial state of the valence electron is important and can bemeasured in ionization experiments [30]. We will thereforeconcentrate on this case.
IV. RESULTS
The results in this work are organized as follows. InSec. IV A the lowest-lying stationary states of the model Liatom introduced in Sec. II are determined. Exact results forthe ionization rate are considered in Sec. IV B in order toidentify different mechanisms behind the ionization processfor different regimes of laser parameters. In Sec. IV C, a ge-dankenexperiment is performed, revealing the mechanisms bywhich seemingly passive “inner” electrons can influence theionization probability of the “outer” electron. In Sec. IV Dthe exact results are compared with TDDFT in various ap-proximations.
A. Stationary states
The highly optimized TDSE solver for propagating thefull three-particle wavefunction (details in the Appendix) isused in the imaginary-time mode for calculating the unper-turbed eigenstates of the model Li. The ground state energy is E = − . . As we are investigating single ionization inthe present work, it is useful to determine the single-ionizationcontinuum threshold. This can be done by comparing withthe ground state of Li + or by following the Rydberg series ofsingly-excited states E m of the neutral Li towards E ∞ , bothleading to an ionization potential E ip = E ∞ − E ≃ . .Table I lists the energies for the lowest eight excited states. TABLE I. Energies E m of the energetically lowest eigenstates m =0 , , , . . . of the Li model atom. m E m − . − . − . − . − . − . − . − . − . As an example, Fig. 10 in the Appendix shows a cut at x = 0 of | φ ( x , x , x ) | for the seventh excited state.The single spin-down component is oriented spatially along x in the partial wavefunction φ . As the spin-down com-ponent necessarily belongs to an inner electron, the extent ofthe probability density | φ | is small in x -direction. Hencethe cut at x = 0 . We further observe the antisymmetry plane x = x and the preference of spatial regions where no morethan one electron is located at a position comparatively farfrom the nucleus. B. Ionization in different laser regimes
The ionization rate Γ ω for a certain photon energy ω wasdetermined by fitting N ( t ) to p ω ( t ) ≃ − exp( − Γ ω t ) (seethe definition of p ( t ) in Sec. III A 4) during the flat-top partof a trapezoidal laser pulse. Given a maximum simulationtime τ ≥ t ≥ , the electric field amplitude ˆ E during theflat-top part of the pulse is chosen high enough to make therelevant ionization timescales for a “numerically measurable”ionization yield shorter than τ . Besides the inverse ioniza-tion rate Γ − ω also the inverse n -photon Rabi frequency Ω − ,n matters here. On the other hand, the ponderomotive energy U p = ˆ E ω should be smaller than both the ionization poten-tial E ip and the photon energy ω . Otherwise, (strong field)effects such as AC Stark shifts, above-threshold ionization,and stabilization [31] could influence the ionization dynamicsdramatically [32]. The electric field amplitude was thereforeset to ˆ E = 0 . .It is obvious that FCAs fail for photon energies ω ≫ E ip high enough to produce core holes. Too low laser frequenciesare numerically too demanding. Hence, the frequencies con-sidered in this work are restricted to ω ∈ [0 . , . . Within thisregime we encounter single-photon ionization for ω ≥ E ip and multiphoton ionization for ω < E ip .
1. Ionization rates
Ionization rates Γ ω obtained by the exact solution of theTDSE are shown in Fig. 2. As expected, one can qualitativelydistinguish between the single-photon and multiphoton ion-ization regime. Photon energy ω (a.u.) R a t e Γ ω ( a . u . ) − − E (2)ip E ip ω (2)4 ω (1)3 ω (1)7 FIG. 2. Logarithmic plot of the ionization rate Γ ω of the Li modelatom vs the laser frequency ω for ˆ E = 0 . , as obtained from the ab initio solution of the TDSE. Ionization thresholds for one andtwo photons are given by the dotted, vertical lines. n -photon reso-nances with the m th excited state are denoted by ω ( n ) m (dashed, ver-tical lines). The regime of multiphoton ionization ω < E ip is dom-inated by resonances at ω ( n ) m and E (2)ip whereas the photoionizationprobability decreases exponentially (i.e., linearly on the logarithmicscale) for ω > E ip , starting at its maximum for ω ≃ E ip . a. Single-photon ionization. If the photon energy in-creases beyond ω = E ip , the ionization rate drops expo-nentially. Note that this behavior can only partially be ex-plained by the decreasing number of photons per time andarea for the fixed laser intensity in the simulation. In a sim-ple picture, the ionization rate Γ ω should be the product ofthe photoionization cross section σ = σ ( ω ) and the photon impact rate per area Γ photon /A . Hence one would expect Γ ω = σ Γ photon /A = σI/ω , where I denotes the (in our caseconstant) laser intensity. Instead, the almost linearly decreas-ing slope in the logarithmic plot of Γ ω in Fig. 2 shows that Γ ω is not ∼ ω − . Hence there must be an exponential dependencein σ ( ω ) . b. Multiphoton ionization. Peaks in the ionization rate Γ ω for ω < E ip can be categorized into two groups. If the en-ergy of n photons is just sufficient to free the outer electron,the ionization probability is particularly high. Consequently,one finds a peak just above the two-photon ionization thresh-old E (2)ip ≃ E ip / . However, the AC Stark shift increases forsmaller laser frequencies ω . Thus, it becomes more importantto consider field-dressed states in order to predict n -photonionization thresholds for n > .Another mechanism that leads to peaks in the multiphotonregime is excited-state-assisted ionization where (a) the en-ergy of n photons matches the energy gap E m − E betweenthe ground state and the m th excited state, (b) the n -photontransition between the ground state and the m th excited stateis allowed, and (c) the binding energy E ∞ − E m of the m th excited state is smaller than the photon energy so that the ab-sorption of one additional photon leads to ionization. Thisscenario may be viewed as n -photon Rabi oscillations, ac-companied by ionization. Laser-dressed states have to be con-sidered in order to precisely predict the peak positions ω ( n ) m ,especially for small laser frequencies.
2. Time-dependent ionization probability
Ionization just above any n -photon ionization thresholdshould depend solely on a single timescale given by the rate Γ ω . Instead, in the case of excited-state-assisted ionizationthe ionization probability vs time should be modulated on thetimescale of (multiphoton) Rabi floppings. This is indeed thecase, as is shown in Fig. 3 where the inverse ionization prob-ability − p ω ( t ) is plotted vs t for four values of the photonenergy ω .
3. Position of resonance peaks
In order to quantitatively predict the position of the peaks inthe ionization rate Γ ω it is required to consider the AC Starkeffect. The laser parameters used in this work are such thatthe coupling of the ground state to states with excited innerelectrons is negligible. Hence, one can obtain Floquet spectraby considering states below the first ionization threshold only.The results of the Floquet solver are shown in Fig. 4. a. Avoided crossings of the shifted ground state. Avoided crossings of the field-dressed ground state are of par-ticular interest here because the system should be describedby this state after an appropriate ramping of the laser field inthe TDSE solution [33].Following the perturbed ground state m = 0 , n =0 from the high-frequency limit to lower frequencies, onefinds avoided crossings with the odd excited states m ∈ t (a.u.) N o r m N ω ( t ) = − p ω ( t ) . . . . . - ph o t o n r e s o n a n ce ω = . = ω ( ) - ph o t o n i o n i z a t i o n ω = . & E i p - ph o t o n r e s o n a n ce ω = . = ω ( ) ph o t o i o n i z a t i o n ω = . > E i p FIG. 3. (Color online) Comparison of time-dependent ionizationprobabilities p ω ( t ) for photon energies ω ∈ { . , . , . , . } and a fixed electric field amplitude ˆ E = 0 . , ramped over five cy-cles. In the case of resonances, i.e., excited-state-assisted ionization,the slope changes periodically with the Rabi frequency. In one- andtwo-photon ionization without excited-state assistance, the only rel-evant timescale is given by Γ − ω , leading to a straight slope in thelogarithmic plot (disregarding the transient behavior at small times t caused by laser ramping). { , , , , } of the next lower Floquet block, as expectedfrom the dipole selection rule. The minimum level distancesin these one-photon avoided crossings is given by the Rabi fre-quency Ω R , , which decreases for increasing m so that thosefor m = 9 and m = 7 are too close to be resolved. This isexpected because Ω R , is proportional to the dipole transitionamplitude, which decreases with increasing m .Following in Fig. 4 m = 0 , n = 0 below the two-photonionization threshold E (2)ip ≃ . , two-photon avoided cross-ings of the ground state with even states m = 8 (unresolved), m = 6 (unresolved), m = 4 (hardly resolved) and m = 2 (clearly resolved) show up. In the case of n -photon crossingsfor n ≥ , the identification of states becomes cumbersome,as Floquet blocks approach each other and the AC Stark shiftincreases. b. Prediction of peaks in the ionization rate. In Fig. 4,photon energies with a high ionization rate are markedby dashed (excited-state-assisted ionization) and dotted ( n -photon ionization thresholds) vertical lines. For each of thesephoton energies, the responsible mechanism can be identi-fied by inspecting the behavior of the state m = 0 , n = 0 .The other way around it is not that straightforward. There areavoided crossings at photon energies ω (1)1 and ω (1)5 which donot give a significant peak in the ionization rate Γ ω . However,for most of the avoided crossings one can quantitatively pre-dict a peak position in the ionization rate, which is supportingthe mechanisms introduced in Sec. IV B. C. Coupling of inner electrons
In this Section, a gedankenexperiment is performed.Halfway between freezing the inner electrons and taking theirdynamics fully into account lies a treatment where only theirinteraction with the laser is neglected while the electron-electron interaction W ij is included in the simulation. How-ever, switching-off the interaction with the laser for both in-ner electrons is not possible in the exact TDSE calculationbecause this would break the exchange symmetry discussedin Sec. III A. An interaction that does not break the exchangesymmetry is a “spin-selective laser” coupling P i H ( i )L ( t ) with[see Eq. (21)] H ( i )L ( t ) = |↑i ( i ) h↑| ( i ) E ( t ) x ( i ) . (26)When solving the TDSE (20) for φ , this hypothetical laseracts on all electrons except the single inner spin-down elec-tron. In that way we can investigate the role of the (seem-ingly passive) inner spin-down electron during the ionizationprocess by discriminating its reaction only to the spin-up elec-trons’ laser-induced dynamics from its full interaction with both the laser and the other electrons.
1. Ionization rates
Ionization rates for the spin-selective case are comparedwith the previous full-laser results in Fig. 5. In the multi-photon ionization regime ω < E ip both curves are in a goodagreement. Positions, heights, and widths of the peaks in bothcases match quantitatively. However, in the single-photon ion-ization regime ω ≥ E ip significant differences in the ioniza-tion rate are observed. Most notably, the ionization rate for thespin-selective laser is too high and shows a wrong asymptoticbehavior with increasing frequency. Hence, the interaction ofinner electrons with the laser field affects the ionization rateeven though all inner electrons stay bound.
2. Dipole expectation values in the laser-driven case
The excursions of loosely outer and tightly inner boundelectrons driven by an oscillating electric field E ( t ) = − ˆ E sin( ωt ) are expected to be opposite in phase if the laserperiod falls in between their respective time scales. In a har-monic binding potential, for example, the phase depends onthe sign of ω − ω with ω the eigenfrequency of the har-monic potential. In the case of a high-frequency driver ω ω < the electron is displaced opposite to the driving force −E ( t ) .On the other hand, a bound electron with a faster timescalethan the driver, ω ω > , is displaced in the direction of thedriving force. In terms of inner and outer electrons this meansthat the position expectation value of an outer electron is morelikely to oscillate in phase with the electric field, whereas in-ner electrons tend to have the opposite phase. Note that “inphase with the electric field” means “opposed to the force”due to the negative charge of the electron. Time-dependent - . - . - . - . - . - . - . Photon energy ω (a.u.) Q u a s i e n e r g i e s ǫ ( n ) m ( ω )( a . u . ) Excitedstateindex m , , , , |{z} F l o q u e t b l o c k i nd e x n = E (2)ip E ip ω (2)4 ω (1)3 ω (1)7 E E E ... FIG. 4. (Color online) Quasienergies ǫ ( n ) m ( ω ) of Floquet states vs photon energy ω for ˆ E = 0 . . Indices denote the state m and the Floquetblock n . Those Floquet states with a sizable projection on the atomic ground state m = 0 are marked by symbols (shaped and coloreddifferently for each Floquet block). Photon energies ω leading to a peak in the ionization rate Γ ω are indicated by gray vertical lines, as inFig. 2. The energies E m with m = 0 , , on the right hand side denote unperturbed atomic energy levels. Photon energy ω (a.u.) R a t e Γ ω ( a . u . ) − − exactspin-selective FIG. 5. Ionization rate Γ ω vs ω in the case of an artificial, spin-selective laser (dashed) compared to the previous results where allelectrons “see” the laser (solid). position expectation values for both spin components of theLi model interacting with a ramped sinusoidal laser field areshown in Fig. 6. Results were obtained for the case of an “or- dinary” laser on the one hand and for the case of the artificialspin-selective laser on the other hand.The amplitude of the single spin-down inner electron is oneorder of magnitude smaller than the amplitude of the inner-outer spin component. It is thus a good approximation toassign the position expectation value of the inner-outer spincomponent to the outer electron. As predicted by the har-monic oscillator the loosely bound outer electron oscillatesin phase with the electric field whereas the inner spin-downelectron oscillates with the opposite phase.The comparison of position expectation values for ordinaryand spin-selective laser shows that switching-off the laser forthe single spin-down inner electron does not affect the oscil-lation amplitude of the outer electron (see both bold curves ontop of each other in Fig. 6). In contrast, the amplitude for thesingle spin-down inner electron itself decreases by a factor oftwo if the ordinary laser (thin solid black) is replaced by thespin-selective (thin dashed red). Time t (a.u.) h x σ i ( t )( a . u . ) , E ( t )( a . u . ) - . . . laser type spin-selectiveusual ↓ (inner) ↑ (inner-outer) E ( t ) FIG. 6. (Color online) Time-dependent position expectation value ofthe spin-down component x ↓ ( t ) (inner electron; thin) compared withthe position expectation value of the other spin component x ↑ ( t ) (in-ner and outer electron; bold). The sinusoidal electric field E ( t ) withfrequency ω = 0 . (thin dashed-dotted, purple) is ramped over fiveperiods and interacts either with all electrons (usual case; solid) orwith all electrons except the single spin-down inner electron (artifi-cial spin-selective case; dashed). In the case of the spin-selective laser, the single spin-downinner electron oscillates with the laser frequency ω althoughit is not directly interacting with the laser. It only couples in-directly to the laser field via the electron-electron interaction W ( ij ) . Both spin-up electrons are directly coupled to the laserfield and repel the spin-down electron. The latter is there-fore slightly displaced in the direction of the electric field bythe other inner electron and in the opposite direction by theouter electron. The net result for the quantum mechanical ex-pectation value is an excursion in the opposite direction (thindashed red), but less so as if it were “seeing” the laser as well(thin solid black).As proven by the significantly differing single-photon ion-ization rates in the gedankenexperiment , the coupling betweeninner and outer electrons strongly affects the ionization pro-cess despite a seemingly harmless approximation: only thelaser interaction of one of the inner electrons is neglected.One can think of the decrease in the ionization rate in the fullsimulation as dynamical anti-screening of the nuclear charge.This is a particular form of dynamical core polarization, whichcould be modeled by adding a polarization potential to theSAE Hamiltonian (see, e.g., [34, 35]).One might wonder why the ionization rate is so differentfor the two laser types in Fig. 5 while the position expectationvalues for the inner-outer spin-alike electrons in Fig. 6 are vir-tually equal. The explanation is that the position expectationvalues are dominated by the bound part of the wavefunctionwhile the ionization rate is determined by the (small) escaping(and numerically absorbed) part.
3. Polarization by a constant electric field
In the limit ω ω ≫ a dipole expectation value x σ ( E ) = h Ψ ( E ) | x σ | Ψ ( E ) i with the sign opposite to E is expected forboth spin components if the electrons were non-interacting.However, the electron-electron interaction W ( ij ) modifies E (a.u.) h x σ i ( E )( a . u . ) − . . ↓ (spin down inner electron) ↑ (inner-outer spin component) FIG. 7. Dipole expectation values for both spin components of the Limodel in the presence of a constant electric field with field strength E (“seen” by all electrons). The inner-outer spin component is dis-placed in the direction opposite to the electric field, as expected fornegatively charged particles. The displacement of the opposite-spincomponent (corresponding to the spin-down inner electron here) isin the same direction as the electric field, as a response to the outerelectron. this.The numerical results in Fig. 7 support the picture of innerelectrons reacting to the displacement of the outer electron.The expectation value of the single spin-down inner electron x ↓ ( E ) has the same sign as the electric field E . D. Ionization rates obtained with different SAE and TDDFTapproximations with and without SIC
After having obtained insight into the role of inner elec-trons in ionization and polarization from ab initio solutions ofthe TDSE, results from full TDDFT and frozen-core calcula-tions are presented and interpreted in this Section. The SAEapproximation, as explained in Sec. III F, is also counted asan FCA in which only the KS valence orbital is propagated inthe frozen ground state KS potential (plus the potential due tothe laser). The SAE results do not suffer from self-interactionintroduced during the propagation of orbitals in time. How-ever, a self-interaction error may originate from the groundstate KS potential.For two approximations we considered, the results havebeen so unreasonable that they are not even shown here but arejust mentioned. First, the pure LDA TDDFT approach com-pletely fails in generating a reasonable behavior of ionizationprobability vs time in the multiphoton regime ω < E ip so thata rate p ω ( t ) ≃ − exp( − Γ ω t ) could not be extracted. Sec-ond, the TDDFT approach using PZ SIC leads to resonancesof inner electrons at high frequencies ω > . , leading to anon-monotonous behavior of the ionization rate not seen inthe exact result.
1. SAE results
Ionization rates from SAE calculations in which frozenKS potentials were used are shown in Fig. 8a. In the low-frequency regime ω < . both SAE ionization rates changerapidly with the frequency, as the exact result does due tothe many avoided crossings discussed in Sec. IV B 3. To the0 R a t e Γ ω ( a . u . ) − − − exact spin-selectivePZ SAE LDA SAE (a) SAE approaches Photon energy ω (a.u.) R a t e Γ ω ( a . u . ) − − − exact spin-selectivePZ TDDFTLDA TDDFT (b) full TDDFT calculations FIG. 8. (Color online) Ionization rate Γ ω obtained by SAE approaches (a) and by full TDDFT calculations (b). For reference, in both panelsthe result from the exact TDSE (labeled “exact”, drawn solid with bullets) and the “spin-selective” TDSE calculation (dotted with bullets) inthe single-photon regime ω > . where it is different from the full calculation (see Fig. 5) are included. The SAE results in (a) are labeled“LDA SAE” (i.e., SAE with frozen LDA groundstate KS potential, dashed) and “PZ SAE” (i.e., SAE with frozen LDA PZ-SIC groundstate KSpotential, solid). The full TDDFT results in (b) are labeled “LDA TDDFT” (dashed) and “PZ TDDFT” (i.e., LDA with PZ-SIC, solid). LDATDDFT results are omitted in the multiphoton regime ω < . , PZ TDDFT in the high-frequency regime ω > . , because a single-ionizationrate cannot be determined in these cases. right of the two-photon ionization peak a sharp minimum at ω min ∈ [0 . , . shows up in all results. In the remainingmultiphoton ionization part up to ω = E ip ≃ . one or twopeaks are visible, corresponding to one-photon excited-state-assisted ionization.The LDA SAE approach yields too few peaks and an in-correct curvature around ω ≃ . . Furthermore, the excited-state-assisted ionization peak around ω ≃ . is blue-shiftedwhile the ionization threshold is red-shifted. This can be par- tially explained by the wrong asymptotic behavior of the KSpotential originating from self-interaction in pure LDA. TheSAE approximation in the PZ-corrected case leads to the cor-rect number of peaks. If each of the n -photon-peaks is shiftedto the right by ∆ ω ≃ . n one finds a striking agreement withthe exact TDSE result with respect to peak positions, widths,and heights. This improvement over pure LDA is due to theasymptotically correct KS potentials in the case of PZ SIC. Infact, three-dimensional DFT calculations applying the PZ SIC1are often useful to quantitatively reproduce experimental val-ues for ionization thresholds and excitation energies. Hence,the required shift ∆ ω may be due to the one-dimensionalityof the Li model system considered in this work.In the single-photon ionization regime ω > E ip a mono-tonic decrease of the ionization rate, starting from its max-imum value for ω ≃ E ip , is observed. However, a convexcurvature as for the spin-selective laser gedankenexperiment arises. Hence, all approaches neglecting inner electron dy-namics completely (pure LDA SAE and PZ-corrected SAE)or partially (spin-selective laser gedankenexperiment ) yield aconvex curvature in the logarithmic plot, presumably becausethe anti-screening effect is neglected.
2. Full TDDFT results
Ionization rates from the full TDDFT calculations areshown in Fig. 8b. As mentioned above, plain LDA TDDFTdoes not allow to extract a rate in the multiphoton ionizationregime at all. Compared to the SAE results, the PZ TDDFTrate in the multiphoton ionization regime appears to be cal-culated with less spectral resolution. This is expected because“unfrozen” KS potentials do not support stationary energy lev-els that could aid ionization via resonant excitations.In the single-photon regime the full TDSE and the LDATDDFT calculations yield an exponential decrease of the ion-ization rate for increasing photon energies. The correspond-ing slopes approximately equal each other. However, the ratepredicted by LDA TDDFT is too high. A possible explana-tion for the overestimated rate is the following: the down shiftof the KS energy during ionization in LDA without SIC (seeSec. IV E below) leads to an increased ionization probabilityin the single-photon regime (because the ionization probabil-ity drops with increasing excess energy ω − E ip ). Hence, onemay view the LDA rate as rather being blue-shifted than up-shifted.Surprisingly, the PZ-corrected TDDFT calculations yieldthe worst rate, which is too small at the threshold, decreasestoo rapidly with increasing ω , and has a concave curvature.The reason for this wrong behavior is discussed in Sec. IV Fbelow.Concluding this subsection, we can state that none of theconsidered approximations is able to yield correct ionizationrates over the frequency interval [0 . , . , the best performingbeing the LDA PZ-SIC SAE approximation in the multipho-ton and the pure LDA TDDFT in the single-photon ionizationregime. E. Importance of SIC for resonances in the multiphotonregime
The PZ-SIC leads to “better” KS energies of the populatedlevels in the sense of being closer to the respective ioniza-tion energies. The energies of the unpopulated levels in thegroundstate KS potential also benefit from the SIC, leading to“better” excitation energies. Both are important ingredients for a correct ionization rate in the multiphoton regime. More-over, in a time-dependent calculation SIC also helps takinginto account the discontinuity of the KS potential at integer or-bital occupation numbers (this is the so-called “derivative dis-continuity” in the XC energy, cf. [36]). The ionization poten-tial must not depend on the ionization probability p ( t ) . Hence,the KS energy of the orbital from which predominantly ioniza-tion occurs should be independent of p ( t ) ∈ [0 , . Only when p ( t ) = 1 is reached, the KS potential, and thus the orbital en-ergy, should change discontinuously. The single-particle den-sity in the vicinity of the nucleus decreases as p ( t ) increases.As a result, the repulsive Hartree potential is weakened so thatall orbital energies would shift to lower values if the XC po-tential did not counteract. In fact, pure LDA does not coun-teract properly so that the KS level energies vary continuouslywith p ( t ) . PZ SIC implements the derivative discontinuity, atleast approximately. This is illustrated in Fig. 9 for a valenceKS orbital occupation n val ∈ [0 , for pure LDA and PZ SICapplied to LDA. By considering fractional occupations in thestationary groundstate calculations we are mimicking adiabat-ically evolving ionization, i.e., − n val = p with p the instan-taneous ionization probability. Fractional valence occupation number n val = 1 − p E n e r g y ( a . u . ) - . - . . positive ion ground state ←− ionization ←− PZPZ LDALDA ε val ∆ E tot FIG. 9. (Color online) Orbital energy of the “valence” KS orbital ε val (solid) and the total energy difference ∆ E tot (dashed) vs thefractional valence occupation number n val = 1 − p , where p is theionization probability, for pure LDA and for LDA with PZ SIC. For the neutral atom ground state n val = 1 the PZ SIClowers the valence orbital energy compared to the uncorrectedLDA case, − .
355 = ε (PZ)val (1) < ε (LDA)val (1) = − . . (27)In the Li + -limit n val = 0 both approaches almost agree, − .
671 = ε (PZ)val (0) ≃ ε (LDA)val (0) = − . . (28)In pure LDA the orbital energy shifts almost linearly duringionization down to small values of n val . With PZ SIC the or-bital energy is shifted much less during ionization as long as n val > . . In the region n val ∈ [0 , . the PZ-SIC orbitalenergy describes a “smoothed jump” down to the Li + -value.Hence, we find that the PZ SIC smoothes the step-function-like behavior the unknown, exact SIC would yield. The im-provement over pure LDA regarding the constancy of the or-bital energy of the “ionizing KS level” is essential.2 F. Problem of LDA PZ-SIC in single-photon ionization
With all the benefits from PZ SIC concerning KS level en-ergies and the asymptotic behavior of the KS potential, it is anobvious question why PZ-corrected LDA fails so badly in thesingle-photon ionization regime. The gedankenexperiment inSec. IV C indicates that dynamical coupling effects betweenelectrons such as anti-screening become increasingly impor-tant as the photon energy ω rises. In TDDFT, the coupling be-tween KS orbitals is mediated by the Hartree part and the XCpart of the KS potential. With PZ SIC the orbital-dependentKS potential reads v i ( x ) = v ( x ) + v (LDA) [ n σ i ]( x ) + v (H) [ n ]( x ) − v (LDA) [ n i ]( x ) − v (H) [ n i ]( x ) . (29)The Hartree potential is a linear functional of the total electrondensity so that v i ( x ) = v ( x ) + v (LDA) [ n σ i ]( x ) − v (LDA) [ n i ]( x ) + v (H) [ n − n i ]( x ) . (30)The effective Hartree term v (H) [ n − n i ]( x ) for the valence KSorbital after SIC is solely determined by the core KS orbitaldensity. As a result, anti-screening is stronger than withoutSIC of the Hartree potential. The SIC of the LDA-part acts inthe opposite direction. However, the SIC to the LDA XC po-tential is not exact, so that a net overestimated anti-screeningmay remain, leading to a lower ionization rate at higher pho-ton energies. Moreover, the dependence of the PZ TDDFTionization rate as a function of the laser frequency is wrongin Fig. 8b, pointing to a deficiency in the dynamics of the XCpotential (note that the PZ SAE rate bends in the opposite di-rection).The fact that pure LDA leads at least to the correct slopeof the ionization rate in the single-photon ionization regimeis thus likely due to a lucky cancellation of errors, i.e., thesuppressed anti-screening is compensated by the error in theLDA XC potential at all frequencies ω ∈ [0 . , . V. CONCLUSION AND OUTLOOK
In this work, we benchmarked various density functional-based approximate approaches to laser ionization with a nu-merically exactly solvable three-electron model atom.In the apparently simple photoeffect regime where only onephoton is required for ionization, a surprisingly pronounceddependence of the ionization dynamics on the correct treat-ment of the inner electrons is found. These inner electrons areusually assumed to be passive, justifying frozen-core, single-active-electron, and pseudopotential approaches. However,because of the opposite timescales of inner and outer elec-trons with respect to the laser period anti-screening of the nu-clear charge by the inner electrons occurs, which is ignoredin such approaches (but could be modeled by a dynamical po-larization potential). For instance, frozen-core orbitals lead toionization rates too high, with an erroneously curved slope of the ionization rate as a function of the laser frequency. Unfor-tunately, the more advanced Perdew-Zunger self-interaction-corrected local density approximation fails as well at high fre-quencies because of an overemphasized anti-screening.In order to correctly describe ionization in the multiphotondomain, energy levels as well as dipole transition probabilitiesmust be reproduced by the simulation method, which is verydemanding for pseudopotential and single-active-electron ap-proaches. Moreover, the energy levels should not change asexcited states get populated because this would move the sys-tem out of resonance. On the other hand, it is known thatonce the population is inverted, the density is the ground statedensity of a “new,” discontinuously changing Kohn-Sham po-tential [12, 37]. Only proper self-interaction free Kohn-Shampotentials may capture such multiphoton ionization effects in-volving resonances. We found that the Perdew-Zunger self-interaction-corrected local density approximation performswell for the lithium model atom in this respect, at least quali-tatively.In future work, it is worth to compare our exact numeri-cal model-Li results with Kohn-Sham results using more ad-vanced exchange-correlation functionals than we did in thecurrent work [38]. There might well be exchange-correlationpotentials “out there” that perform well in both the multipho-ton and single-photon ionization regime.
ACKNOWLEDGMENT
This work was supported by the SFB 652 of the GermanScience Foundation (DFG). Discussions with V. Kapoor aregratefully acknowledged.
Appendix A: Numerical solution of the three-particle TDSE
In this Appendix, numerical details concerning the solu-tion of the TDSE (20) are given. The full solution of thethree-particle TDSE for processes involving ionization is nu-merically demanding even if the spatial degrees of freedomare reduced to one per particle. For that reason a CartesianTDSE solver in three dimensions (i.e., one per particle) wasimplemented on a graphics processing unit (GPU) using theNVIDIA R (cid:13) CUDA TM [39] platform. The solver is highly opti-mized for the purpose of this paper. A speedup of three ordersof magnitude over a single-core CPU implementation cover-ing the full 3D Cartesian grid is achieved on a desktop com-puter featuring an NVIDIA R (cid:13) GeForce R (cid:13) GTX 580 GPU. Thespeedup originates equally from an efficient implementationutilizing the high single-precision performance of the GPUand the adjustment of the simulation grid. “Mixed precision”techniques allow to obtain results in double precision althoughmost of the numerical effort consists of single-precision float-ing point operations.The kinetic energy operator is discretized using the implicitNumerov expression, which is accurate up to fourth order inthe spatial gridspacing ∆ x . For the propagation in time weemploy the unitary Crank-Nicolson method consisting of an3 -100 -50 0 50 100 - - (cid:12)(cid:12) φ (7)23 ( x ; x , x ) | x =0 (cid:12)(cid:12) x (a.u.) x ( a . u . ) − − − − − − double ionization,core excitations:not treated same here s a m e h e r e a n t i s y m m e t r y p l a n e n o tt r e a t e ddu e t o a n t i s y m m e t r y FIG. 10. (Color online) Spin-resolved probability density | φ (7)23 ( x , x , x ) | of the seventh excited state cut at x = 0 . Theadapted grid omits those regions with a white background. (explicit) forward and an (implicit) backward step, accurateup to second order in the timestep ∆ t . The explicitly time-dependent Hamiltonian H ( t ) is evaluated at midpoints τ = t + ∆ t .The linear equations to propagate the discretized wavefunc-tion in time for one timestep can easily be arranged to requirethe solution of no more than a single implicit equation with a × × stencil. For a sufficiently small timestep ∆ t the corre- sponding coefficient matrix is diagonal-dominant. Hence, theJacobi method can be used to determine its solution. Fasterconvergence of long-wavelength errors is achieved by apply-ing a multigrid scheme. Finally, a high throughput of floatingpoint operations is obtained by using mixed precision tech-niques and different data caching stages combined with mas-sive parallelization.The numerical grid was optimized, exploiting what isknown about the probability density dynamics during long-wavelength, single ionization. Recall that the laser is tunedsuch that solely the outer electron is removed from the atom.As a consequence, the single spin-down inner electron istightly bound to the core at all times. This means that | φ ( t ; x , x , x ) | will only yield non-vanishing probabili-ties for small | x | . Hence, one may choose a small box size inthe x -direction. In contrast to that, the inner-outer spin com-ponent in the x and x -directions may be spatially extended.However, if both | x | and | x | are large, | φ ( t ; x , x , x ) | must be negligible for the allowed laser parameters, i.e., thosethat do not lead to double (or triple) ionization. Consequently,the grid regions where | x | and | x | are large are omitted. Fi-nally, by making use of the antisymmetry φ ( t ; x , x , x ) = − φ ( t ; x , x , x ) only the region x ≤ x needs to be con-sidered numerically.Denoting the width of an ionization channel in grid pointsby N small and its total length by N large , the total number ofgridpoints compared to an N small × N large × N large cuboidis reduced by a factor of N large /N small . The results pre-sented in this work have been checked to be converged for N small = 64 and N large = 1024 (for ∆ x = 0 . and ∆ t = ∆ x/ ), which corresponds to a speedup of causedby the grid adjustment alone. As an example, the probabilitydensity | φ (7)23 ( x = 0 , x , x ) | of the seventh excited state onthe optimized grid is shown in Fig. 10. [1] R. G. Parr and W. Yang, Density Functional Theory of Atomsand Molecules (Oxford Science Publications, 1989).[2] R. M. Dreizler and E. K. U. Gross,
Density Functional Theory:An Approach to the Quantum Many-Body Problem (SpringerBerlin Heidelberg, 1990).[3] E. Engel and R. M. Dreizler,
Density Functional Theory: AnAdvanced Course (Springer Berlin Heidelberg, 2011).[4] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[5] L. Szasz,
Pseudopotential Theory of Atoms and Molecules (Wi-ley New York, 1985).[6] F. Nogueira, A. Castro, and M. Marques, in
A Primer in DensityFunctional Theory , Lecture Notes in Physics, Vol. 620, editedby C. Fiolhais, F. Nogueira, and M. Marques (Springer BerlinHeidelberg, 2003) pp. 218–256.[7] M. A. L. Marques, N. T. Maitra, F. M. S. Nogueira, E. K. U.Gross, and A. Rubio, eds.,
Fundamentals of Time-DependentDensity Functional Theory , Lecture Notes in Physics, Vol. 837(Springer Berlin Heidelberg, 2012).[8] C. Ullrich,
Time-dependent density-functional theory: conceptsand applications (Oxford University Press, 2012).[9] C. A. Ullrich, The Journal of Chemical Physics , 234108(2006). [10] P. Hessler, N. T. Maitra, and K. Burke, The Journal of ChemicalPhysics , 72 (2002).[11] Y. Kurzweil and R. Baer, Phys. Rev. B , 075413 (2006).[12] N. T. Maitra, K. Burke, and C. Woodward,Phys. Rev. Lett. , 023002 (2002).[13] M. Protopapas, C. H. Keitel, and P. L. Knight,Reports on Progress in Physics , 389 (1997).[14] D. B. Miloševi´c, G. G. Paulus, D. Bauer, and W. Becker, Jour-nal of Physics B: Atomic, Molecular and Optical Physics ,R203 (2006).[15] K. C. Kulander, K. J. Schafer, and J. L. Krause,International Journal of Quantum Chemistry , 415 (1991).[16] M. Awasthi, Y. V. Vanne, A. Saenz, A. Castro, and P. Decleva,Phys. Rev. A , 063403 (2008).[17] O. Smirnova, Y. Mairesse, S. Patchkovskii, N. Du-dovich, D. Villeneuve, P. Corkum, and M. Y. Ivanov,Nature , 972 (2009).[18] A. E. Boguslavskiy, J. Mikosch, A. Gijsbertsen, M. Spanner,S. Patchkovskii, N. Gador, M. J. J. Vrakking, and A. Stolow,Science , 1336 (2012).[19] R. Grobe and J. H. Eberly, Phys. Rev. Lett. , 2905 (1992).[20] C. Ruiz, L. Plaja, and L. Roso, Phys. Rev. Lett. , 063002 (2005).[21] Had we chosen the Li configuration S = 1 / and M S = − / instead of M S = +1 / , the role of single-particle spin-up andspin-down components would be reversed, of course.[22] C. J. Joachain, N. J. Kylstra, and R. M. Potvliege, Atoms inIntense Laser Fields (Cambridge University Press, 2012).[23] S.-I. Chu and D. A. Telnov, Physics Reports , 1 (2004).[24] P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964).[25] E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997 (1984).[26] J. P. Perdew and A. Zunger, Phys. Rev. B , 5048 (1981).[27] Rigorously speaking, the fictitious KS particles should not beviewed as “real” electrons. However, it seems not unreasonableto take the valence KS orbital as the starting point for a SAEcalculation.[28] H. G. Muller, Phys. Rev. A , 1341 (1999).[29] C.-T. Chen and F. Robicheaux, Phys. Rev. A , 3261 (1996).[30] M. Meckel, D. Comtois, D. Zeidler, A. Staudte, D. Paviˇci´c,H. C. Bandulet, H. Pépin, J. C. Kieffer, R. Dörner, D. M. Vil- leneuve, and P. B. Corkum, Science , 1478 (2008).[31] P. Mulser and D. Bauer, in High Power Laser-Matter Interac-tion , Springer Tracts in Modern Physics, Vol. 238 (Springer-Verlag, Berlin Heidelberg, 2010) pp. 267–330.[32] Note, however, that all methods used in this work can be appliedto the nonperturbative regime as well. We restrict ourselves torather small field strengths here for purely methodical reasons.[33] V. Kapoor and D. Bauer, Phys. Rev. A , 023407 (2012).[34] Z. Zhao and J. Yuan, Phys. Rev. A , 023404 (2014).[35] B. Zhang, J. Yuan, and Z. Zhao,Phys. Rev. Lett. , 163001 (2013).[36] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Phys. Rev.Lett. , 1691 (1982).[37] M. Ruggenthaler and D. Bauer, Phys. Rev. Lett. , 233001(2009).[38] S. Kümmel and L. Kronik, Rev. Mod. Phys. , 3 (2008).[39] NVIDIA Corporation, “Parallel programmingand computing platform | cuda | nvidia,”