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²⁺ 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)|². 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²⁺ are also given and compared with experimental data (http://physics.nist.gov/PhysRefData/contents-atomic.html). Without laser, the

Table 1. Stationary Ne configurations as obtained from our TDDFT solver. All values in Hartrees. For conversion to eV multiply with 27.21.

View Table

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²⁺, and Ne³⁺ 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²-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²⁺ 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.

 

Fig. 1. Ionization probability of the KS orbitals after a 20 cycle pulse vs. the peak intensity of the pulse. In (a) the result for the 2p0 orbital from the full TDKS run (red) is compared with the DFT-SAE result one would expect for sequential ionization (blue: first electron; green: second electron). In (b) the full TDKS results for the 2p1 and 2s orbitals are shown. Results obtained on a different numerical grid, appropriate for higher intensities, are drawn in orange.

Download Full Size | PDF

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₁ for the probability to find the 2_p1 orbital ionized, and $\overline{p_1}=1-p_1$)

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²-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).

 

Fig. 2. Ne ion yields after a 20 cycle laser pulse vs. its peak intensity. The red solid, dashed, and dashed-dotted lines are results for singly, doubly, and triply ionized Ne, respectively. The blue line is the DFT-SAE result for the first 2_p0 KS electron. The green line is the sequential double ionization result for the second 2_p0 KS electron, again from a DFT-SAE run. Results obtained on a different numerical grid are drawn in orange. The symbols are experimental yields from [6] (diamonds: Ne⁺, triangles: Ne²⁺).

Download Full Size | PDF

Again, the DFT-SAE results are plotted also (blue: Ne⁺, green: Ne²⁺). Although a NSI-knee is not present, the red dashed Ne²⁺-curve clearly predicts an enhancement of two-fold ionized Ne for intensities <10¹⁵ W/cm² 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}^{+}⁄2\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 ⁿ⁺ 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²⁺-yield (dashed, red) and the DFT-SAE Ne²⁺-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²⁺-curve at the right intensity is clearly present.

Another problem appears at low ionization probabilities: the TDDFT Ne²⁺ curve has a slope too steep when compared to experiment. It seems that at low intensities the experimental Ne²⁺ 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²⁺-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¹⁴W/cm². 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²⁺ yield. This was discussed in [10] but no effective potential was found which produced the NSI knee automatically.

 

Fig. 3. Focus averaged results for singly (solid) and doubly (dashed and dotted, red) ionized Ne, calculated as described in the text. The green, dashed curve is the DFT-SAE Ne²⁺-result. The symbols are from experiment, like in Fig. 2 (no relative shift in intensity was necessary).

Download Full Size | PDF

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²⁺ 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²⁺, 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.