## Abstract

Time-dependent density functional theory (TDDFT) is employed to study the interaction of a Ne atom with short and strong 800 nm laser pulses. In the intensity regime covered (10^{14}–10^{16} W/cm^{2}) up to triply ionized Ne is observed. Good quantitative agreement with the experimental Ne^{+} ion-yield (and the Ne^{2+}-yield near saturation) is obtained. Nonsequential ionization (NSI) leads to a strong increase of the probability for double and triple ionization when compared to a single active electron (SAE)-approach. A NSI-“knee” is observed but the agreement with its experimental counterpart is not satisfactory.

© Optical Society of America

## 1 Introduction

With the discovery of the so-called nonsequential ionization (NSI) “knee” [1] it became clear that electron-correlation can lead to measurable effects in intense laser atom interaction not accessible by means of single active electron (SAE)-approaches. Since full ab initio simulations, i.e., solving the many-body time-dependent Schrödinger equation for more than two active electrons is beyond the capabilities of even the best super-computers, alternative approaches are indispensable. Density functional theory (DFT), an in principle exact mean field approach to the many-body problem [2], is extremely successful in atomic physics, solid state physics, and chemistry. Therefore it seems reasonable to suppose that time-dependent density functional theory (TDDFT) [3] may give useful insight also in the interaction of strong laser light with many electronatoms. Time-dependent Thomas-Fermi theory can be considered the simplest version of TDDFT and has been applied to multiple ionization of noble gases in laser fields [4].

In this paper we perform a TDDFT study of Ne in an 800 nm laser pulse, and we try to reproduce measured Ne^{+} and Ne^{2+} ion yields.

## 2 Method

The time-dependent Kohn-Sham (TDKS) equation for the orbital Ψ* _{iσ}*(

**r**,

*t*) reads (in atomic units)

Here, *σ*=↑, ↓ indicates the spin polarization, *V*(*r*)=-*Z*/*r* is the ionic potential, *V _{I}*(

*t*) is the laser in dipole approximation, and ${V}_{{\mathrm{ee}}_{\sigma}}[{n}_{\uparrow},{n}_{\downarrow}]$ is the effective electron-electron interaction potential which is a functional of the electron spin densities

*n*(

_{σ}**r**,

*t*)=∑

^{N}^{σ}_{i=1}|Ψ

_{iσ}(

**r**,

*t*)|

^{2}.

*N*is the number of orbitals occupied by Kohn-Sham (KS) particles with spin σ. The total electron density is

_{σ}*n*(

**r**,

*t*)=∑

_{σ}

*n*(

_{σ}**r**,

*t*). The electron-electron potential is splitted, ${V}_{{\mathrm{ee}}_{\sigma}}[{n}_{\uparrow},{n}_{\downarrow}]=U\left[n\right]+{V}_{\mathrm{xc}\sigma}[{n}_{\uparrow},{n}_{\downarrow}]$, where

*U*[

*n*] is the Hartree potential

and ${V}_{{\mathrm{xc}}_{\sigma}}[{n}_{\uparrow},{n}_{\downarrow}]$ is the exchange correlation-part, preferably with a self-interaction correction (SIC), for which an approximation has to be used in practice. We chose the Slater expression

where ${u}_{{\mathrm{xc}}_{i\sigma}}={V}_{{\mathrm{xc}}_{\sigma}}^{\mathrm{XLDA}}[{n}_{\uparrow},{n}_{\downarrow}]-U\left[{n}_{i\sigma}\right]-{V}_{{\mathrm{xc}}_{\sigma}}^{\mathrm{XLDA}}[{n}_{i\sigma},0]$, i.e., the self-interaction is correctly removed, and the exchange-only local density approximation (XLDA) was employed.

In our numerical code, the KS orbitals are expanded in spherical harmonics. The radial wave functions are discretized in position space. Each time step, the effective potential has to be calculated self-consistently, which makes a TDKS solver more time-consuming than an ordinary time-dependent Schrödinger equation (TDSE) code running a corresponding SAE problem. The effective potentials were expanded up to the dipole. We are confident that this approximation does not affect the validity of our conclusions since we could reproduce all results presented in this paper even stopping after the monopole. We will come back to this issue later on. The actual propagation is performed in velocity gauge using an algorithm similar to the one for the TDSE described by Muller in Ref. [5]. As usual, probability density reaching the boundary of the numerical grid at a radius *r*≈100 a.u. was removed by an imaginary potential. The eventually decreasing norm then can be interpreted as one minus the ionization probability of that orbital.

## 3 Results

#### 3.1 Stationary configuration

In order to obtain quantitatively acceptable results for ionization yields and rates it is essential to start from the proper ground state configuration. The stationary orbital configurations were obtained by imaginary time propagation (with the laser switched off) in combination with Gram-Schmidt orthonormalization each time step. For real time propagation, the 1s electrons were kept frozen since they cannot be removed by laser intensities we study in this work. This relaxed the requirements concerning the grid spacing, and a Δ*r*=0.05 was found to be sufficient to obtain the orbital energies and the total energies E displayed in Table 1. Ionization potentials *I _{p}* for neutral Ne, Ne

^{+}, and Ne

^{2+}are also given and compared with experimental data (http://physics.nist.gov/PhysRefData/contents-atomic.html). Without laser, the

two 2_{p0} and four 2_{p1} electrons are degenerate. Therefore, only three KS orbitals are actually needed to obtain the results of Table 1. Spin-polarization effects were neglected in the calculations of the Ne^{+}, Ne^{2+}, and Ne^{3+} open shell-results so that three KS orbitals were still sufficient. As expected, the 2p electrons are most loosely bound. With the laser switched on, *four* KS orbitals need to be propagated since the 2_{p0} and 2_{p1} KS electrons behave differently in the laser field.

## 3.2 Orbital ionization probabilities

In Fig. 1 we present results for the ionization probability of the individual KS orbitals after a 20 cycle 800 nm pulse vs. the peak intensity of the pulse. The electric field had a sin^{2}-shape. The orbital ionization probability was calculated as one minus the norm of the orbital after the pulse. Note that this cannot be directly interpreted as the probability to find Ne, Ne^{+}, Ne^{2+} etc. because each orbital represents two KS electrons (the 2_{p1} even four) so that some statistics has to be performed in order to construct functionals of the orbitals which give the correct single, double, triple etc. ionization probability. How to do this rigorously is, to our knowledge, an unsolved problem. However, from Fig. 1a useful information can be inferred: when compared to the DFT-SAE run (with only one 2_{p0} KS electron active and all others frozen, drawn blue) the 2_{p0} KS electron in the full run (drawn red) ionizes less likely. This is expected because it can share energy with the other KS electrons through their common mean field. At higher intensities it even ionizes less likely than the second 2p0 electron in the DFT-SAE run for Ne^{+} (green curve) indicating even stronger collective effects.

In Fig. 1b the results for the two other active orbitals are presented. Comparing the ionization probability in (b) with the one in (a) one sees that the 2_{p1} and 2s KS electrons will play a minor role in single and double ionization. It is interesting that the 2s result shows a knee-structure. However, since the 2p orbitals are almost exclusively responsible for up to triply ionized Ne this knee has nothing to do with the one observed in experiment [6]. The orange curves in Fig. 1 stem from runs we performed with a smaller time step but with all effective electron-electron potentials expanded up to the monopole only. For higher orbital ionization probabilities there is no difference at all, and the knee structure in the 2s result is still clearly visible. This is in contrast to calculations [7] where the dipole term of the Hartree potential was found to be essential to produce a knee within the so-called CRAPOLA model.

## 3.3 Yields

As mentioned before it is not clear how to calculate from the KS orbital data ionization yields which can be compared with experiment. One reason is that the KS particles are auxiliary entities which must not be confused with the “real” electrons. However, the easiest way to construct functionals for multiple ionization probabilities directly from the KS orbital ionization data is to adopt an independent particle-like statistical viewpoint. The same method was applied in [8] where Ne in a 248 nm laser field was studied. In this approximation, we have (using the notation, e.g., p_{1} for the probability to find the 2_{p1} orbital ionized, and $\overline{{p}_{1}}=1-{p}_{1}$)

$$+6{\overline{s}}^{2}{\overline{{p}_{0}}}^{2}{\overline{{p}_{1}}}^{2}{{p}_{1}}^{2}+8\overline{s}\phantom{\rule{.2em}{0ex}}\overline{{p}_{0}}\phantom{\rule{.2em}{0ex}}{\overline{{p}_{1}}}^{3}\left(s\overline{{p}_{0}}{p}_{1}+\overline{s}{p}_{0}{p}_{1}\right),$$

as the probabilities to find neutral, singly and doubly ionized Ne, respectively. Similar, although more lengthy expressions for higher ionization states can be easily written down.

In Fig. 2 results for Ne ion yields after 20 cycle 800 nm laser pulses (again sin^{2}-shaped), calculated according Eqs. (5) and (6), are presented. The red and orange curves are single, double, and triple ionization yields (solid, dashed, and dashed-dotted, respectively). The symbols are experimental data for a longer (200 fs) pulse from [6] (experimental ion yields typically show a continuing increase after saturation due to focal expansion; we account for this in the following Subsection).

Again, the DFT-SAE results are plotted also (blue: Ne^{+}, green: Ne^{2+}). Although a NSI-knee is not present, the red dashed Ne^{2+}-curve clearly predicts an enhancement of two-fold ionized Ne for intensities <10^{15} W/cm^{2} by several orders of magnitude compared to the DFT-SAE result in green. In TDDFT, NSI, when calculated from (6), has its origin rather in statistics than in a particular scenario such as rescattering. This can be seen by comparing the leading terms in (5) and (6): for low intensities ${P}^{2+}\approx {p}_{0}{P}^{+}\u20442\overline{{p}_{0}}\approx {p}_{0}{P}_{\mathrm{SAE}}^{+}$ holds.

If a NSI-knee is present it is due to the transition from NSI to the sequential behavior because of the depletion of neutral Ne. Therefore we can draw the preliminary conclusion that either the pulse was too short to yield depletion of neutral Ne, or there are intrinsic problems, e.g., with functional (6).

## 3.4 Rates

From our TDDFT studies, by plotting 1-*P*
^{n+} vs. time, we determined ionization rates as a function of the laser intensity. For sequential single and double ionization we determined the rates using our DFT-SAE results. The nonsequential double ionization rate was determined from the full TDDFT result. Then, to simulate the experimental conditions of [6], those rates were plugged into the rate equations for the singly and doubly ionized Ne ion yields which we solved for a set of 200 fs laser pulses with Gaussian temporal shape and various peak intensities. Subsequently, a focal averaging was performed, leading to the results presented in Fig. 3.

The Ne^{+}-yield (solid, red) agrees very well with the experimental data, as the Ne^{2+}-yield (dashed, red) and the DFT-SAE Ne^{2+}-yield (dashed, green) do as far as the saturation intensity is concerned. However, the NSI-knee is not well reproduced although a small dip in the red, dashed Ne^{2+}-curve at the right intensity is clearly present.

Another problem appears at low ionization probabilities: the TDDFT Ne^{2+} curve has a slope too steep when compared to experiment. It seems that at low intensities the experimental Ne^{2+} yield is almost “parallel” to the singly ionized Ne. The same appears to be true for He [9]. The only theory which reproduces also this behavior correctly is, to our knowledge, the one in [11]. Instead, the empirical NSI rate proposed in [12] leads to results similar to ours in that respect.

The red, dotted curve is the Ne^{2+}-result when the NSI rate was set to zero for intensities greater than where the DFT-SAE single ionization yield crosses the full TDDFT result in Fig. 2, i.e., around 6.8×10^{14}W/cm^{2}. Of course, this artificial cutting of the rate is not satisfactory, and a better approach is required, i.e., improved multiple ionization functionals are needed. It is an interesting question whether, with a more accurate effective potential, the functionals (4)–(6) would yield better results concerning the Ne^{2+} yield. This was discussed in [10] but no effective potential was found which produced the NSI knee automatically.

## 4 Conclusion

We presented TDDFT studies of Ne in short 800 nm laser pulses. The agreement with experiment concerning the ion yields for Ne^{+} (and Ne^{2+} near saturation) is very good when the DFT ground state is used, and only the 2_{p0} KS electron is propagated (DFT-SAE). The TDDFT runs with all 2s and 2p KS electrons active show NSI, i.e., an increased probability for Ne^{2+}, when calculated using the simple functionals (4)–(6). Although a NSI-knee is clearly visible the agreement with experiment is rather poor.

Finally, we want to mention that, concerning NSI of He [9], we obtained with our TDDFT code results similar to those presented in this paper. However, we plan to perform more studies, especially with respect to the frequency dependence of NSI in TDDFT.

## Acknowledgements

This work was supported in part by the Deutsche Forschungsgemeinschaft within the SPP “Wechselwirkung intensiver Laserfelder mit Materie.”

## References and links

**1. **D. N. Fittinghoff, P. R. Bolton, B. Chang, and K. C. Kulander, “Observation of nonsequential double ionization of helium with optical tunneling,” Phys. Rev. Lett. **69**, 2642–2645 (1992). [CrossRef]

**2. **R. M. Dreizler and E. K. U. Gross, “Density Functional Theory: An Approach to the Quantum Many-Body Problem,” (Springer, Berlin, 1990).

**3. **Erich Runge and E. K. U. Gross, “Density-Functional Theory for Time-Dependent Systems,” Phys. Rev. Lett. **52**, 997–1000 (1984). [CrossRef]

**4. **Miroslaw Brewczyk, Kazimierz Rzazewski, and Charles W. Clark, “Appearance intensities for multiply charged ions in a strong laser field,” Phys. Rev. A **52**, 1468–1473 (1995). [CrossRef]

**5. **H. G. Muller, “An Efficient Propagation Scheme for the Time-Dependent Schrödinger Equation in the Velocity Gauge,” Laser Physics **9**, 138–148 (1999).

**6. **S. Larochelle, A. Talebpour, and S. L. Chin, “Non-sequential multiple ionization of rare gas atoms in a Ti:Sapphire laser field,” J. Phys. B: At. Mol. Opt. Phys. **31**, 1201–1214 (1998). [CrossRef]

**7. **J. B. Watson, A. Sanpera, D. G. Lappas, P. L. Knight, and K. Burnett, “Nonsequential Double Ionization of Helium,” Phys. Rev. Lett. **78**, 1884–1887 (1997). [CrossRef]

**8. **C. A. Ullrich and E. K. U. Gross, “Many-electron atoms in strong femto-second laser pulses: A density-functional study,” Comm. At. Mol. Phys. **33**, 211 (1997).

**9. **B. Walker, B. Sheehy, L. F. DiMauro, P. Agostini, K. J. Schafer, and K. C. Kulander, “Precision Measurement of Strong Field Double Ionization of Helium,” Phys. Rev. Lett. **73**, 1227–1230 (1994). [CrossRef]

**10. **M. Petersilka and E. K. U. Gross, “Strong-Field Double Ionization of Helium: A Density-Functional Perspective,” Laser Physics **9**, 105–114 (1999).

**11. **A. Becker and F. H. M. Faisal, “*S*-matrix analysis of ionization yields of noble gas atoms at the focus of Ti:sapphire laser pulses,” J. Phys. B: At. Mol. Opt. Phys. **32**, L335–L343 (1999). [CrossRef]

**12. **U. Eichmann, M. Dörr, H. Maeda, W. Becker, and W. Sandner, “Collective Multielectron Tunneling Ionization in Strong Fields,” Phys. Rev. Lett. **84**, 3550–3553 (2000). [CrossRef]