Momentum distributions in time-dependent density functional theory:
  Product phase approximation for non-sequential double ionization in strong
  laser fields
Momentum distributions in time-dependent density functional theory:
Product phase approximation for non-sequential double ionization
in strong laser fields
F. Wilken and D. Bauer
Max-Planck-Institut für Kernphysik, Postfach 103980, 69029 Heidelberg, Germany
(Dated: October 24, 2018)
We investigate the possibility to deduce momentum space properties from time-dependent density
functional calculations. Electron and ion momentum distributions after double ionization of a
model Helium atom in a strong few-cycle laser pulse are studied. We show that, in this case, the
choice of suitable functionals for the observables is considerably more important than the choice
of the correlation potential in the time-dependent Kohn-Sham equations. By comparison with
the solution of the time-dependent Schrödinger equation, the insufficiency of functionals neglecting
electron correlation is demonstrated. We construct a functional of the Kohn-Sham orbitals, which
in principle yields the exact momentum distributions of the electrons and the ion. The product-
phase approximation is introduced, which reduces the problem of approximating this functional
significantly.
PACS numbers: 31.15.Ew, 32.80.Rm
I. INTRODUCTION
Time-dependent density functional theory (TDDFT)
[1] is a remarkably successful approach to the study
of many-body systems in time-dependent external fields
[2]. The essential statement of TDDFT is the same as
that of the well-established ground state density func-
tional theory (DFT) [3]: all observables are, in princi-
ple, functionals of the particle density alone. Since the
latter is always a three-dimensional entity, independent
of the number of particles involved, the computational
cost of actual (TD)DFT calculations scales exponentially
more favorable than the solution of the many-body (time-
dependent) Schrödinger equation.
In practice, almost all (TD)DFT calculations are per-
formed using the (time-dependent) Kohn-Sham scheme
[(TD)KS] (see, e.g., [2]) where the density is calculated
with the help of auxiliary, non-interacting particles mov-
ing in an effective potential. The “art” of (TD)DFT is
two-fold, namely finding sufficiently accurate approxima-
tions to the density functionals of (i) the unknown effec-
tive potential and (ii) the observables of interest. For-
tunately, for many practical applications both items are
uncritical [2]. An example is the calculation of the opti-
cal response of bio-molecules where even the simple lo-
cal density approximation of the effective potential yields
reasonable results, and the observable can be calculated
from a known and explicit functional of the density (the
time-dependent dipole).
However, when it comes to the correlated motion of a
few particles in a strongly driven system, TDDFT faces
major challenges. In that respect, non-sequential dou-
ble ionization (NSDI) serves as the “worst case” scenario
for TDDFT. Theoretically, NSDI was addressed success-
fully using the strong-field approximation (see, e.g., [4]
and references therein) and classical methods [5, 6]. The
widely accepted mechanism behind NSDI relies on the
rescattering of the first electron with its parent ion, col-
lisionally ionizing (or exciting) the second electron.
In the recent publications Refs. [7, 8] significant
progress was made in the treatment of NSDI within
TDDFT as far as ionization yields are concerned. The
latter display as a manifestation of the electron-electron
correlation involved in NSDI the celebrated “knee” struc-
ture in the double ionization yield, which was, until re-
cently, not being reproduced within TDDFT. Reference
[7] addressed issue (i) above (the effective potential) while
Ref. [8] focused on item (ii), the functional for the observ-
able “double ionization”. It was shown that (i) taking
the derivative discontinuities at integer bound electron
numbers into account and (ii) using an adiabatic approx-
imation for the correlation function needed to calculate
the double ionization probability, the NSDI “knee” can
be reproduced.
In our current work we turn to the much harder prob-
lem of momentum distributions (or energy spectra [9]).
In the NSDI regime the ion momentum spectra, as mea-
sured in experiments employing “reaction microscopes”
(see, e.g., [4, 10]), show a characteristic “double-hump”
structure, i.e., maxima at non-vanishing ion momenta.
The maxima at non-zero ion momenta are easy to under-
stand within the rescattering scenario mentioned above:
the first electron preferentially returns to the ion, colli-
sionally ionizing the second electron, at times when the
vector potential of the laser field is non-zero. Since the
vector potential at the ionization time equals the final
drift momentum at the detector, non-vanishing electron
momenta (and, due to momentum conservation, non-
vanishing ion momenta) are likely. In a TDKS treat-
ment of NSDI in He starting from a spin-singlet state,
the rescattering scenario is “hidden” in a single, spa-
tial Kohn-Sham (KS) orbital. As we shall demonstrate,
taking the auxiliary KS particles for real electrons and
Fourier-transforming their position space product wave-
function to momentum space leads to ion momentum
http://arxiv.org/abs/0704.1423v1
spectra in very poor agreement with the exact ones. A
better approximation to calculate correlated electron mo-
mentum spectra in the NSDI regime is required. With
the present paper we aim at contributing to this goal by
showing that item (ii) above, namely the construction
of the functional for the observable, is the critical issue
while (i) known effective potentials are sufficient, at least
at the current level of accuracy.
In Sec. II the model Helium system used to study the
ionization process and the ensuing momentum distri-
butions is introduced. In Sec. III the method to cal-
culate electron and ion momentum distributions is ex-
plained. Results from the solution of the time-dependent
Schrödinger equation (TDSE) in Sec. IV serve as a refer-
ence for the results obtained using TDDFT in Sec. V:
The insufficiency of uncorrelated functionals to calcu-
late electron and ion momentum distributions (Sec. VA)
and the relative insignificance of the correlation poten-
tial (Sec. VB) lead us to the construction of correlated
functionals in Sec. VC. In Sec. VD we introduce the
product-phase approximation, which reduces the prob-
lem of approximating the correlated functionals to that
of approximating the exchange-correlation function.
For consistency, we restrict ourselves to the presenta-
tion of results for laser pulses with λ=780 nm and N=3
cycles. We stress, however, that the general conclusions
drawn hold also for λ= 614 nm, N = 3 and λ= 780 nm,
N=4 laser pulses, as we have checked explicitly.
II. MODEL SYSTEM
A Helium atom exposed to linearly polarized laser
pulses with N = 3 cycles and sin2-pulse envelopes is in-
vestigated. The length of the pulses with a frequency
of ω = 0.058 (corresponding to the experimentally used
λ = 780 nm) is T = 2Nπ/ω, and the vector potential
readsA(t) = Â sin2
sin(ω t) for 0 ≤ t ≤ T and zero
otherwise (atomic units are used unless otherwise indi-
cated). We use the dipole approximation, i.e., the spatial
dependence of the laser field is neglected. The linear po-
larization of the laser pulse thus allows us to describe the
system by a one-dimensional model Helium atom with
soft-core potentials for the Coulomb interactions. It is
known that the essential features of the ionization pro-
cess are described well by this model [7, 8, 11, 12, 13, 14].
Initially, the electrons are assumed to occupy the spin-
singlet groundstate of Helium, and due to the neglect
of magnetic effects in the dipole approximation the elec-
trons stay in the spin-singlet state during the interaction
with the laser pulse. Thus it is sufficient to study the
spatial wavefunction, which has to be symmetric under
exchange of the electrons.
The TDSE i∂t ψ(x1, x2, t) = Ĥ(x1, x2, t)ψ(x1, x2, t) is
solved for laser pulses with different effective peak inten-
sities I = I(Â). A trivial gauge-transformation cancels
the purely time-dependent A2-term and yields the Hamil-
tonian
i=1,2
∂2xi + V (xi, t)
+W (|x1 − x2|) , (1)
with Ĥ = Ĥ(x1, x2, t). The external potential is
V (x, t) = −iA(t) ∂x − 2/
x2 + ǫen, the electron-electron
interaction potential is given by W (x) = 1/
x2 + ǫee.
The soft-core parameters ǫen and ǫee are chosen to yield
the correct ionization potentials. Reproducing the ion-
ization potential of He
, I(2)p = 2.0 in a corresponding
model He
ion fixes ǫen = 0.5. The choice ǫee = 0.329
yields the ionization potential of Helium, I(1)p = 0.904.
All results presented in this work are qualitatively insen-
sitive to the precise values of the soft-core parameters.
As the two electrons constitute a spin-singlet state
for all times they are described by the same KS or-
bital. Therefore, in a TDDFT treatment, we have
only one time-dependent Kohn-Sham equation (TDKSE)
i ∂t φ(x, t) = Ĥ
KS(x, t)φ(x, t) with the Hamiltonian
ĤKS(x, t) = −1
∂2x + V (x, t) + vhx(x, t) + vc(x, t) . (2)
The Hartree-exchange potential vhx = vh + vx follows as
vhx(x, t) =
dx′ n(x′, t)/
(x− x′)2 + ǫKSee . We have
used the exact exchange term for Helium vx(x, t) =
−vh(x, t)/ 2, which is local as both electrons are de-
scribed by the same orbital.
Setting vc = 0 yields, in the special case of the
Helium atom or He-like ions, an identical description
as the time-dependent Hartree-Fock (TDHF) treatment
(due to the locality of vx). The LK05 potential v
[7] takes into account the discontinuous change in the
correlation potential when the number of bound elec-
trons N(t) =
dxn(x, t) passes integer numbers,
vLK05c (x, t) = [B(t)/ (1 + exp[C(B(t) − 2)])− 1] vhx(x, t),
where C is a sufficiently large constant (we set C = 50)
and B(t) = N0/N(t). In order to encompass all bound
states the parameter a is chosen as a = 6 a.u. throughout
this work, results being insensitive to the precise value of
a. We use ǫKSee = 0.343 in the Hartree-exchange potential
vhx to acquire I
p = 0.904 for the model Helium atom.
The TDSE and TDKSE are solved by a split-operator
time propagator on a numerical grid (see, e.g., [15] and
references therein).
Along the lines of Ref. [7] we construct from the
TDSE solution an exact KS orbital (EKSO). The
Schrödinger solution gives the exact density of our
model Helium atom n(x, t) = 2
dx2 |ψ(x, x2, t)|2 =
dx1 |ψ(x1, x, t)|2 and the exact probability current
j(x, t). From the equality of the exact and KS
currents in the case of a one-dimensional system,
the phase of the EKSO is determined as ϑ(x, t) =∫ x
dx′ j(x′, t)/n(x′, t) + α(t). The unknown purely
time-dependent phase factor α(t) does not affect the re-
sults presented in this work and is therefore set to zero.
The EKSO φ(x, t) =
n(x, t)/ 2 eiϑ(x,t) is thus identical
to the orbital a TDDFT calculation with the exact cor-
relation potential vc would yield via the TDKS scheme.
The EKSO allows us to separate the challenges facing
TDDFT calculations (cf. Sec. I): finding (i) a suitable
approximation of vc (where it serves as a reference for
the resulting orbital) and (ii) appropriate functionals for
observables (where it is the exact input).
III. MOMENTUM DENSITIES
We partition the two-electron space and associate
with single ionization the area A (He+) = {(x1, x2) |
|xi| > a, |xj 6=i| ≤ a ∀ i, j ∈ {1, 2}} and with dou-
ble ionization the area A (He2+) = {(x1, x2) | |x1| >
a, |x2| > a}. Integrating |ψ(x1, x2, t)|2 over these ar-
eas then yields the respective ionization probabilities,
with the double ionization probability given by P 2
(t) =∫∫
A (He2
dx1 dx2 |ψ(x1, x2, t)|2. This scheme to deter-
mine ionization probabilities from the two-electron wave-
function has been successfully used in numerous similar
calculations [11, 12, 14].
The wavefunction ψ(x1, x2, t) can be de-
scribed equivalently in momentum space by
its Fourier transform (2 π) ψ(k1, k2, t) =∫
dx2 ψ(x1, x2, t) e
−i (k1 x1+k2 x2). As the wave-
function in momentum space is normalized to one,
the pair density in momentum space is given by
ρ(k1, k2, t) = 2 |ψ(k1, k2, t)|2.
At times 0 < t < T during the laser pulse the velocity
of the electrons is actually given by ẋi(t) = ki(t) +A(t),
i.e., the sum of the canonical momentum ki and the value
of the vector potential at the respective time. In this
work we investigate properties of the system at t = T
after the laser pulse. As A(T )=0, canonical momenta k
and drift momenta are identical.
We are interested mainly in the double ionization pro-
cess and thus Fourier transform only the wavefunction
in the area A (He2+) associated with double ionization.
The resulting sharp step at the boundary of A (He2+) at
|xi| = a, |xj 6=i| ≥ a, with i, j ∈ {1, 2}, is a potential
source of artifacts when Fourier transformed. Hence, a
smoothing function f(x1, x2) =
i=1 1/
1 + e−c |xi−a|
is introduced. The factor c has to be of the order of
one, in this work we choose c = 1.25. The smoothing
function is constructed so that
dx1dx2 f
2(x1, x2) b =∫∫
A (He2
dx1 dx2 b for a constant b. This condi-
tion ensures that the wavefunction ψ(2
+)(x1, x2, t) =
f(x1, x2)ψ(x1, x2, t) gives to a good approximation the
same double ionization probability as the original wave-
function, i.e., that
dx2 f
2(x1, x2)|ψ(x1, x2, t)|2 ≃
. The correlated wavefunction of the electrons freed
in double ionization in momentum space is thus calcu-
lated as
(2 π) ψ(2
+)(k1, k2, t) =∫
dx2 ψ
(2+)(x1, x2, t) e
−i (k1 x1+k2 x2) . (3)
This approach is equivalent to projecting out the states
corresponding to single and no ionization and is known
to lead to accurate momentum distributions [13].
From the wavefunction we construct the momentum
pair density of the electrons freed in double ionization
+)(k1, k2, t) = 2 |ψ(2+)(k1, k2, t)|2 . (4)
The probability to find at time t an electron freed in dou-
ble ionization with momentum k1 in dk1 and an electron
with k2 in dk2 is then ρ
(2+)(k1, k2, t) dk1dk2.
In experiments, it is easier to measure the momentum
of the He2
ion kIon after double ionization instead of in-
dividual electron momenta. As the total photon momen-
tum involved is negligibly small, this provides informa-
tion about the sum of the electron momenta via momen-
tum conservation k1 + k2 = −kIon. The ion momentum
density then follows from the momentum pair density of
the electrons freed in double ionization (4) as
Ion (kIon, t) =
dk ρ(2
+)(−kIon− k, k, t)
dk ρ(2
+)(k,−kIon− k, t) , (5)
due to the symmetry of the electron momentum pair
density. The factor 1/2 ensures the correct normaliza-
tion since the system consists of only one ion but two
electrons. The ion momentum density n
Ion (kIon, t) dkIon
gives the probability to find at time t the He2
ion with
momentum kIon in dkIon.
IV. MOMENTUM DISTRIBUTIONS FROM
THE TDSE
From the numerical solution of the TDSE we obtain
ψ(x1, x2, T ) after the interaction with the laser pulse. In
the left hand side of Fig. 1 the momentum pair density
of the electrons freed in double ionization, as calculated
from Eq. (4), is shown.
For all but the highest intensity depicted, electrons
have the highest probability to move at different veloc-
ities |k1| 6= |k2| (ẋi(T ) = ki(T ) since A(T ) = 0, cf. dis-
cussion in Sec. III) but in the same direction (sgn(k1) =
sgn(k2)). Depending on the laser intensity the proba-
bility for the double ionization process is highest at dif-
ferent half-cycles of the laser pulse, i.e., different signs
of the vector potential. Therefore, the favored direction
in which the electrons leave the atom varies with inten-
sity. NSDI can be understood by a recollision mechanism
where one electron returns to the He
ion and frees the
second electron (see, e.g., [4]). The results of the TDSE
then imply that both electrons leave the atom in the same
direction but due to Coulomb repulsion their velocities
differ, in accordance with earlier results for a longer laser
pulse [13].
The “butterfly” shape of the momentum pair density of
the electrons freed in double ionization as shown in Fig. 1
is evidence that it is highly correlated, as it cannot be
reproduced by multiplying two orbitals for the respective
electrons.
For I = 6.96 × 1015W/cm2 both electrons have the
highest probability to leave the atom in the same direc-
tion with similar velocities k1 ≈ k2. This can only be the
case when the Coulomb repulsion between the electrons
is weak, i.e., when they are removed sequentially, result-
ing in a large spatial separation. The final non-vanishing
velocities are due to the high intensity of the laser pulse,
which ionizes the atom so rapidly that A(t) 6= 0 when the
first electron is freed. The grid-like structure typical for
a product wavefunction is seen, the electron correlation
being weak.
From the momentum pair density of the electrons freed
in double ionization ρ(2
+)(k1, k2, T ) (4) we calculate the
ion momentum density n
Ion (k1, k2, T ) (5). For different
effective peak intensities the density of the ion momen-
tum is depicted in Fig. 2. It exhibits peaks at non-zero
momenta. As explained in Sec. I, these are typical for rec-
ollision processes when the first freed electron recollides
close to the maximum of the vector potential, i.e., when
|A(t)| ≈ Â. Hence, the sum of the momenta of both elec-
trons is non-zero, and, by momentum-conservation, this
holds for the ion momentum as well [16].
For an infinitely long laser pulse of laser period T/N ,
Ĥ(t + T/N) = Ĥ(t) holds while this symmetry is bro-
ken in the case of few-cycle laser pulses. Hence, with
respect to the dislodged electrons there is no spatial in-
version symmetry, leading to asymmetric ion momentum
distributions [17, 18, 19]. This effect is clearly seen in
Fig. 2. For the three lowest intensities a process with
kIon ≥ 0 dominates while with increasing intensities pro-
cesses with kIon ≤ 0 become more likely. In addition, a
central peak gets more and more pronounced, showing
that the relative probability of sequential double ioniza-
tion increases. The fact that the peak is not centered
around kIon = 0 for I =6.96 × 1015W/cm2 is again due
to the high intensity and the short duration of the laser
pulse, as explained above.
V. MOMENTUM DISTRIBUTIONS FROM
TDDFT
DFT can be formulated in momentum space (see, e.g.,
[20]), and this seems to be the obvious path to follow
when one is interested in the calculation of momentum
spectra. However, momentum space DFT lacks the “uni-
versality” feature of the Hohenberg-Kohn theorem [3],
meaning that each system under study requires a dif-
ferent momentum space effective potential—an entirely
unattractive feature. We therefore prefer to make the
“detour” via standard, universal, position space TDDFT.
In the case of single ionization, a straightforward cal-
culation of the momentum or energy spectrum from
the Fourier-transformed valence KS orbital may be a
good approximation (see, e.g., the approach followed in
Ref. [21]). Instead, it is less obvious how to determine
correlated momentum spectra from position space TDKS
orbitals.
As explained in the Introduction, determining momen-
tum pair densities and ion momentum densities from a
TDDFT approach faces two challenges: The first is to
find an approximate correlation-potential vc in the TD-
KSE to reproduce the exact density n(x, t) with sufficient
accuracy. The second, more difficult one, amounts to as-
sign a suitable functional of the density to the respective
observable. As both the ion momentum density and the
momentum pair density (via their probability interpreta-
tions, cf. Sec. III) are observables, the Runge-Gross the-
orem assures that functionals of the density alone exist
A. Uncorrelated functionals
Treating the KS orbital as if it were a one-
electron wavefunction yields a product wavefunction
φ(x1, t)φ(x2, t). This is the same assumption frequently
made to derive uncorrelated ionization probability func-
tionals (see Ref. [8] and references therein).
The Fourier transformed KS orbital for |x| > a, i.e.,
with the bound states projected out (see Sec. III) is
2 π φ(
+)(k, t) =
dx f(x)φ(x, t) e−i k x , (6)
with f(x) = 1/
1 + e−c |x−a| the one-dimensional
smoothing function equivalent to the smoothing func-
tion used in Sec. III. Calculating the momentum pair
density (4) and the ion momentum density (5) from the
product wavefunction gives the uncorrelated functional
for the momentum pair density of the electrons freed in
double ionization
+)(k1, k2, t) = 2 |φ(+)(k1, t)φ(+)(k2, t)|2 (7)
and the uncorrelated functional for the ion momentum
density of He2
Ion (kIon, t) =
dk |φ(+)(−kIon− k, t)φ(+)(k, t)|2 . (8)
Equations (7) and (8) are not functionals of the den-
sity alone but due to the Fourier transformation they are
dependent on the density and on the phase of the KS
orbital.
The momentum pair density at t = T , as calculated
from the uncorrelated functional (7) using the EKSO,
is depicted in the right part of Fig. 1 for λ = 780 nm,
N =3-cycle laser pulses with different intensities. Com-
parison with the left hand side showing the momentum
pair density calculated from the correlated Schrödinger
wavefunction ψ(x1, x2, T ) confirms that only for the high-
est intensity a product wavefunction approach is reason-
able. For lower intensities the uncorrelated functional
FIG. 1: Contour plots of the momentum pair density
2+(k1, k2, T ) of the electrons freed in double ionization. Re-
sults calculated from the uncorrelated functional (7) using
the EKSO (right hand side) are compared to the TDSE (left
hand side) solution. Momentum pair densities for λ=780nm,
N =3-cycle laser pulses with different effective peak intensi-
ties are shown.
for the momentum pair density does not exhibit the
typical “butterfly”-shaped correlation structures of the
Schrödinger solution. Instead, the grid-like structure typ-
ical for a product wavefunction is clearly visible.
For the same system we calculate from Eq. (8) the ion
momentum density using the EKSO. In Fig. 2 the He2
ion momentum density is compared to the results from
the TDSE, which are scaled to enable the comparison of
qualitative features. The different values of the integrals
over the ion momentum densities are due to the differ-
ent double ionization probabilities, as can be seen from∫
dkIon n
Ion (kIon, t) ≃ P 2+, which follows from Eq. (3)
(see Ref. [8] and references therein for a discussion of this
particular problem). Apart from the highest intensity
the density is centered around a central peak at kIon ≈ 0.
This is evidence that correlations, which are not included
in the uncorrelated functionals for the observables, are
responsible for the distinct peaks of the ion momentum
density at non-zero momenta. This result is consistent
with the analysis of the results of the TDSE (Sec. IV),
which attributes the peaks at kIon 6= 0 to electron rescat-
tering, i.e., to an interaction between the electrons. For
the highest intensity shown in Fig. 2, sequential double
ionization becomes dominant (cf. Sec. IV), so that the de-
scription using the EKSO in the uncorrelated functional
reproduces the ion momentum density reasonably well.
FIG. 2: Ion momentum density of the model He2
ion after
interaction with λ=780 nm, N=3-cycle laser pulses with dif-
ferent effective peak intensities. The density calculated using
the EKSO in the uncorrelated functional (8) is compared to
results from the TDSE.
B. The role of the correlation potential
To underline the importance of the functional for the
ion momentum density we use the correlation potentials
vc = 0 (TDHF) and v
c (LK05) in the TDKSE for our
model He atom interacting with the λ=780 nm, N =3-
cycle laser pulses (cf. Sec. II).
In Fig. 3 the ion momentum densities obtained from us-
ing the respective orbitals in the uncorrelated functional
for the ion momentum density (8) are compared to the
results with the EKSO, i.e., the orbital which the exact vc
would yield. For the TDHF approach, results are similar
to the results using the LK05-potential. Both approx-
imations lead to uncorrelated ion momentum densities
which are close in qualitative terms to the EKSO results.
Only at the highest intensity I = 6.96×1015W/cm2 they
exhibit a single peak at kIon ≥ 0 and not, as the EKSO
solution, at kIon ≤ 0. In this intensity regime purely
sequential double ionization dominates, pointing to pos-
sible shortcomings in the description of this process with
both correlation potentials.
As the general deficiencies of the uncorrelated func-
tional described in the previous paragraph are entirely
due to the functional for the observable, these results
demonstrate the relative unimportance of the choice of
the correlation potential in the TDKSE for the observ-
ables of interest in this work.
FIG. 3: Ion momentum density of the model He2
ion after
interaction with λ=780 nm, N=3-cycle laser pulses with dif-
ferent effective peak intensities. The densities are calculated
from the uncorrelated functional (8) using the EKSO and the
orbitals obtained with vc = 0 (TDHF) and v
c (LK05).
C. Towards correlated functionals
In polar representation, the solution of the TDSE is
written as ψ(x1, x2, t) =
ρ(x1, x2, t)/ 2 e
iϕ(x1,x2,t) and
the KS orbital as φ(x, t) =
n(x, t)/ 2 eiϑ(x,t). We define
a time-dependent complex exchange-correlation function
κ(x1, x2, t) =
ψ(x1, x2, t)√
2 φ(x1, t)φ(x2, t)
gxc e
i [ϕ(x1,x2,t)−ϑ(x1,t)−ϑ(x2,t)] (9)
with the time-dependent exchange-correlation func-
tion gxc = gxc(x1, x2, t) given by gxc(x1, x2, t) =
ρ(x1, x2, t)/ n(x1, t)n(x2, t). Approximations to gxc =
|κ|2 have been used to construct correlated ionization
probability functionals [8, 22]. Note that while gxc is
an observable (and thus a functional of only the density
exists), the complex-valued κ is not an observable. Using
Eq. (9) to express the correlated wavefunction ψ(x1, x2, t)
in terms of the KS orbitals and the complex exchange-
correlation function, Eq. (4) gives the correlated func-
tional for the momentum pair density of the electrons
freed in double ionization
+)(k1, k2, t) = π
dx2 κ(x1, x2, t)
×φ(+)(x1, t)φ(+)(x2, t) e−i (k1 x1+k2 x2)
with φ(
+)(x, t) = f(x)φ(x, t). The correlated ion mo-
mentum density is calculated by using the correlated mo-
mentum pair density in Eq. (5). We thus have exact mo-
mentum distribution functionals, which depend only on
the complex exchange-correlation function κ and the KS
orbital φ.
The complex exchange-correlation function κ in turn
depends on the pair density and the phase of the
Schrödinger solution ψ(x1, x2, t). In order to derive mo-
mentum space properties for more complex atoms than
Helium from the KS orbitals directly through expressions
like Eq. (10), it is inevitable to approximate κ. However,
this is challenging since, due to the Fourier-integrals in
Eq. (10), the complex exchange-correlation function has
to be approximated in all A(He2+) (and not just for the
bound electrons, as in the calculation of ionization prob-
abilities [8, 22]).
D. Product phase approximation
The necessary approximation of the complex exchange-
correlation function κ (9) consists of approximating
gxc(x1, x2, t) and the phase-difference ϕ(x1, x2, t) −
ϑ(x1, t)− ϑ(x2, t).
Addressing the second part, the easiest approximation
follows from the assumption that the difference of the
sum of the phases of the KS orbitals and the phase of
the correlated wavefunction can be neglected when cal-
culating momentum distributions, i.e., we set
ϕ(x1, x2, t) = ϑ(x1, t) + ϑ(x2, t). (11)
Since ϑ(x, t) is the phase of the KS orbital we denote
this approach as the product phase (PP) approximation,
which yields
κPP(x1, x2, t) =
gxc(x1, x2, t) . (12)
It is noteworthy that knowledge of the exact κPP thus suf-
fices to calculate the exact double ionization probabilities
from the EKSO.
We calculate the ion momentum density using Eq. (12)
in Eq. (10) and in Eq. (5). Employing the EKSO, the
ion momentum densities shown in Fig. 4 for λ=780 nm,
N = 3-cycle laser pulses with different intensities are
obtained. The results from the TDSE are depicted as
well. For comparison of the qualitative features, they
are scaled, although the integrals over the ion momen-
tum densities are equal in both cases (note that the PP
approximation returns the exact double ionization proba-
bilities). A generally good qualitative agreement with the
Schrödinger solution is acquired. The asymmetric struc-
ture and distinct peaks are reproduced. For intensities
where NSDI is strongest, the quantitative agreement is
least convincing. Although the PP approximation does
not reproduce the exact kIon positions of the peaks, it
modifies the uncorrelated functionals in a way which al-
lows to deduce information about the underlying double
ionization processes at the different intensities. We can
therefore conclude that the difference between the phase
FIG. 4: Ion momentum density of the model He2
ion calcu-
lated from the correlated functionals in the PP approximation
using the EKSO. Results for λ = 780 nm, N = 3-cycle laser
pulses with different effective peak intensities are compared
to the ion momentum density obtained from the TDSE.
FIG. 5: Contour plots of the exchange-correlation function
gxc(x1, x2, t) for two effective peak intensities of λ=780 nm,
N =3 cycle laser pulses as acquired from the solution of the
TDSE. For clarity values larger than 10 are shown as 10.
of the correlated wavefunction and a product wavefunc-
tion is not as important for reproducing the structure
of the ion momentum density as is the correlation given
by gxc(x1, x2, t). This conclusion was verified by setting
gxc = 1 in Eq. (9) and using the exact phases in Eq. (10),
which did not yield the peaks present in the Schrödinger
solution. Using LK05 orbitals in the PP approximation
also reproduces distinct peaks while the general agree-
ment with the Schrödinger ion momentum density is not
as good as for the EKSOs.
The contour plots of the momentum pair density of the
electrons freed in double ionization ρ(2
+)(k1, k2, t) calcu-
lated from the correlated functional in the PP approx-
imation using the EKSO show a correlated structure,
while differences from the TDSE momentum pair den-
sities (Fig. 1) remain.
Using the PP approximation we obtain momentum
distributions which yield fundamental insight into the
double ionization processes. However, this still requires
knowledge of the exact gxc(x1, x2, t) at time t = T af-
ter the laser pulse, i.e., of the exact pair density in
real space. Approximating gxc(x1, x2, t) is a formidable
task itself. This can be seen from the highly correlated
structure in Fig. 5 where contour plots of the exchange-
correlation function gxc(x1, x2, T ) are shown for intensi-
ties where NSDI dominates. An adiabatic approximation
using the groundstate pair density [8] is not feasible as
the exchange-correlation function in the entire A(He2+)
is required in Eq. (10). An expansion for small inter-
electron distances [22, 23] will not include the correla-
tions for large |x1−x2|, which are clearly present in Fig. 5.
By multiplying the complex exchange-correlation func-
tion with a damping function F (|x1 − x2|) with F → 0
for large |x1 − x2|, we verified that short-range correla-
tions alone in the final wavefunction do not reveal the
characteristic peaks in the ion momentum density. It is
therefore of central importance to devise new strategies
of approximating gxc(x1, x2, t).
VI. SUMMARY
A model Helium atom in strong linearly polarized few-
cycle laser pulses was investigated. Solution of the time-
dependent Schrödinger equation yielded momentum pair
distributions of the electrons freed in double ionization
and corresponding ion momentum densities. They were
consistent with a recollision process and, at higher laser
intensities, with sequential double ionization. These re-
sults served as a reference for a time-dependent density-
functional treatment of the system. It was shown that
the choice of the correlation potential in the Kohn-
Sham equations is of minor importance compared to the
form of the functionals for calculating the momentum
distributions. An uncorrelated approach was found to
produce ion momentum densities differing significantly
from the Schrödinger solution in qualitative terms. We
constructed an exact correlated functional via the two-
electron wavefunction. The product-phase approxima-
tion reduces the problem of approximating this func-
tional.
This work was supported by the Deutsche Forschungs-
gemeinschaft.
[1] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997
(1984).
[2] M. A. L. Marques, C. A. Ullrich, F. Nogueira, A. Rubio,
K. Burke, and E. K. U. Gross, eds., Time-Dependent
Density Functional Theory (Springer, Berlin Heidelberg,
2006).
[3] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864
(1964).
[4] A. Becker, R. Dörner, and R. Moshammer, J. Phys. B:
At. Mol. Opt. Phys. 38, S753 (2005).
[5] L.-B. Fu, J. Liu, J. Chen, and S.-G. Chen, Phys. Rev. A
63, 043416 (2001).
[6] P. J. Ho, R. Panfili, S. L. Haan, and J. H. Eberly, Phys.
Rev. Lett. 94, 093002 (2005).
[7] M. Lein and S. Kümmel, Phys. Rev. Lett. 94, 143003
(2005).
[8] F. Wilken and D. Bauer, Phys. Rev. Lett. 97, 203001
(2006).
[9] V. Véniard, R. Täıeb, and A. Maquet, Laser Phys. 13,
465 (2003).
[10] J. Ullrich, R. Moshammer, A. Dorn, R. Dörner, L. P. H.
Schmidt, and H. Schmidt-Böcking, Rep. Prog. Phys. 66,
1463 (2003).
[11] D. Bauer, Phys. Rev. A 56, 3028 (1997).
[12] D. G. Lappas and R. van Leeuwen, J. Phys. B: At. Mol.
Opt. Phys. 31, L249 (1998).
[13] M. Lein, E. K. U. Gross, and V. Engel, Phys. Rev. Lett.
85, 4707 (2000).
[14] N. E. Dahlen and R. van Leeuwen, Phys. Rev. A 64,
023405 (2001).
[15] D. Bauer and P. Koval, Comput. Phys. Comm. 174, 396
(2006).
[16] A. Rudenko, K. Zrost, B. Feuerstein, V. L. B. Jesus,
C. D. Schröter, R. Moshammer, and J. Ullrich, Phys.
Rev. Lett. 93, 253001 (2004).
[17] X. Liu, H. Rottke, E. Eremina, W. Sandner, E. Gouliel-
makis, K. O. Keeffe, M. Lezius, F. Krausz, F. Lindner,
M. G. Schatzel, et al., Phys. Rev. Lett. 93, 263001 (2004).
[18] H. Rottke, X. Liu, E. Eremina, W. Sandner, E. Gouliel-
makis, K. O. Keeffe, M. Lezius, F. Krausz, F. Lindner,
M. G. Schätzel, et al., J. of Mod. Opt. 53, 149 (2006).
[19] C. Figueira de Morisson Faria, X. Liu, A. Sanpera, and
M. Lewenstein, Phys. Rev. A 70, 043406 (2004).
[20] R. M. Dreizler and E. K. U. Gross, Density Functional
Theory. An Approach to the Quantum Many-Body Prob-
lem (Springer, Berlin Heidelberg, 1999).
[21] A. Pohl, P.-G. Reinhard, and E. Suraud, Phys. Rev. Lett.
84, 5090 (2000).
[22] M. Petersilka and E. K. U. Gross, Laser Phys. 9, 105
(1999).
[23] A. D. Becke, J. Chem. Phys. 88, 1053 (1988).
