1. science
2. physics

Laser Physics, Vol. 12, No. 2, 2002, pp. 355–361.
STRONG FIELD
PHENOMENA
Original Text Copyright © 2002 by Astro, Ltd.
Copyright © 2002 by MAIK “Nauka /Interperiodica” (Russia).
Photo-Ionization of Two-Electron Atoms:
How to Solve a 5-Dimensional Schrödinger Equation
on a Desktop PC
H. G. Muller
FOM-Institute for Atomic and Molecular Physics, Kruislaan 407, SJ Amsterdam, 1098 The Netherlands
e-mail: muller@amolf.nl
Received October 1, 2001
Abstract—We discuss some numerical techniques that help to make numerical simulation of ionization of a
three dimensional, two-electron system more efficient. A method to obtain highly accurate ground state, and a
novel absorbing boundary condition are discussed in detail.
1. INTRODUCTION
Recently, efficient algorithms have been developed
for numerically solving the time-dependent Schrödinger
equation for an atom in a laser field [1]. These calculations make use of the single-active-electron approximation (SAE) [2], which means that the atom is modeled
as a single electron moving in a potential well describing the nuclear attraction and an average repulsion by
all other electrons. This well is assumed to be static,
i.e., the other electrons are supposed to not react to the
laser field.
For efficient computation of such problems, the
choice of coordinates plays a crucial role. Solutions to
the equations often are rapidly oscillating in one direction, with very little structure in others. In quantum
mechanics the rapid dependence coincides with the
direction of motion of the particle: the wave vector represents the momentum. In photo-ionization the velocity
is directed radially far away from the atom, and can be
quite large. A description in spherical coordinates thus
concentrates all rapid variation along the radial coordinate, requiring only a sparse grid to represent the angular dimensions. Note that, for instance, in Cartesian
coordinates the large radial wave vector would have
non-zero components along each of the coordinates,
and even adaptive step-size techniques would require
fastly more points than in the spherical-coordinate
case, especially in three dimensions.
In this paper we discuss how the techniques that
have proved so successful in the SAE case [3–6] can be
applied to calculations in two-electron systems. In particular we will focus on the question of non-sequential
double ionization [7]: the creation of doubly charged
ions at intensities where the second ionization step on
an isolated singly charged ion is not possible.
The most important issue in making a highly dimensional calculation manageable is the choice of grid.
Any savings here relax the memory requirement
(increasing the number of accessible PCs that could be
put on the job) and reduce the required CPU time
(which is proportional to the number of grid points).
A two-electron problem has six spatial dimensions,
three for each electron [8]. In the dipole approximation
the laser field has cylinder symmetry, and for linear
polarization this symmetry persists in time. This
enables one to separate the overall rotation around the
polarization axis from the other coordinates. The associated angle Φ thus is a “must” as one of the dimensions
of our coordinate system, and the dependence of the
wave function on this coordinate can be solved analytically as exp(iMΦ). The quantum number M is a constant of the motion, and in the helium ground state M =
0. In other words, the separability ensures that this
dimension collapses to a single point.
With the success of the SAE methods in mind, the
five remaining dimensions are described as polar coordinates (ri , ϑi ) for the individual electrons i = 1, 2,
within their (half-)planes through the polarization axis,
and the angle ϕ between those planes. For helium the
Hamiltonian is given by
H = H 1 + H 2 + 1/r 12
(1)
atom
,
(2)
1 2
= --- p i – 2/r i .
2
(3)
with
Hi = Hi
laser
+ Hi
and
atom
Hi
In the double ionization of a helium atom in the ground
state, the wave function is expected to have very weak
dependence on the ϕ-dimension. The initial state is
almost completely independent of it: when expanded as
a Fourier series in this dimension, only 0.5% of the
norm of the wave function resides in terms with a
ϕ dependence exp(imϕ) with m ≠ 0. (Note that ϕ is the
355
356
MULLER
difference of the azimuthal angles ϕi of the individual
electrons, so that the exponent factorizes as the product
of two one-electron wave functions. For each term in
the m expansion the ϕ motion of the electrons is thus
still independent, and we can say that each of the electrons has the quantum number m or – m.) The bulk of
the wave function (99.5%) has m = 0. In the low-frequency regime, ionization of the first electron is dominated by tunneling. This helps in purifying the m = 0
component even more, since all other components have
to leave one electron in an excited state of the ion, causing a vastly bigger barrier [9]. Once the two electrons
get separated far along the polarization axis, the potential energy is hardly dependent on ϕ, and correlation in
this direction can only establish itself after tunneling
when the electron recollides with its parent ion.
The other two angular dimensions (ϑi ) are treated as
an expansion over angular momenta li . This reduces the
angular part of the one-electron Hamiltonians Hi to a
simple multiplicative potential (the centrifugal barrier,
l(l + 1)/2r2), rather than a derivative operator with a singular prefactor. Otherwise, treating a function as a linear combination of L Legendre polynomials or as an
interpolation polynomial through L support points is
completely equivalent. The wave function is thus written as
Ψ =
∑ψ
m
l 1 l 2 ( r 1,
m
–m
r 2 ) ( Y l1 ( ϑ 1, ϕ 1 )Y l2 ( ϑ 2, ϕ 2 )
l 1, l 2, m
–m
(4)
m
+ Y l1 ( ϑ 1, ϕ 1 )Y l2 ( ϑ 2, ϕ 2 ) ),
where the radial dimensions are discretized on an equidistant grid. Only the operator 1/r12 couples the different m subspaces.
The time integration is done as a multiply splitoperator scheme. Each partial propagator of a piece HP
of the Hamiltonian over a time-interval 2τ is approximated as (1 + iHPτ)–1(1 – iHPτ). The propagators for the
one-electron operators Hi commute. Each of them is
treated like in the one-electron problem [1]: split into
an atomic and a laser part, half-step propagators for the
laser
laser interaction H i surrounding a full-step propagaatom
tor of H i . Their product is then surrounded by halfstep propagators for the electron repulsion. The latter is
itself split according to a multipole expansion.
To elucidate the mechanism for non-sequential double ionization [7], simulations have to cover only a limited region of the five-dimensional space. As soon as
both electrons have left the vicinity of the nucleus, say
beyond a distance R, it seems natural that the presence
of a strong laser field should in the end set those electrons free. This distance R can be pretty small: at the
laser fields the double ionization takes place (E ≈ 0.1),
even the Coulomb attraction of a doubly charged
nucleus can only compete with it to a distance of about
5 Bohr. It thus seems safe to set R to a value around
8 Bohr. To see how the wave function moves to the
region with both r1 and r2 > R, we only have to follow
it through a region where at least one of the ri is smaller
than R. We will take this “inner” electron to be electron 2.
Electron 1, on the other hand, will then have to be followed to a distance far enough until it can have no possible effect on the ejection of electron 2. In practice, this
means that we have to follow this electron to a distance
so large that even the quiver motion can not bring it
back close to the nucleus again, and this can mean several hundred Bohr.
As is the single-electron problem, the angular
dependence is very sensitive to the gauge in which the
problem is formulated. In the length gauge the laser
interaction is described by the scalar potential E · r,
which can grow very large far away from the atom,
where it dominates tile Hamiltonian. In this case it
induces a factor U = exp(iA · r), with ∂t A = E, into the
wave function. For large r this factor has both a very
rapid angular and temporal dependence. It is therefore
advantageous to eliminate this factor from the numerical calculation, by assuming its presence implicitly and
solve for the quantity χ = exp(–iA · r)Ψ. This could be
viewed as a kind of rotating-wave anzats, which happens to coincide with a gauge transform: The function χ
obeys the velocity-gauge form of the Schrödinger equation, where the laser field is represented by the operator
A · p.
It should be noted that the velocity gauge only simplifies the description for weakly-bound electrons,
where the motion is dominated by the laser field. For a
tightly bound ground state, the laser field only acts as a
small perturbation, and the factor U is not naturally
present in the wave function at all. Using velocity
gauge for the description of such states thus only complicates the description of such states. In the helium
problem, where ionization occurs around 1015 W/cm2,
the field in atomic units is 0.168. At 800 nm this makes
A = 3. The wave function is largely static up to the saddle point in the combined Coulomb laser potential,
which occurs at r = 2.5. Expanding the factor U in
spherical waves up to this distance would require about
20 angular momenta. For this reason we only incorporate the effect of E · r1 in the wave function, by making
the ansatz χ = exp(–iA · r1)Ψ, after which χ obeys the
equation
1 2 1 2
i∂ t χ = ⎛ --- p 1 + --- p 2 – 2/r 1 – 2/r 2 + 1/r 12
⎝2
2
(5)
--- + A ( t ) ⋅ p 1 + E ( t ) ⋅ r 2 ⎞ χ.
⎠
This ansatz is thus equivalent to treating the two
electrons in a different gauge. The fact that r2 is rather
small in the entire region where we do the simulation
ensures that this electron can simply never built up a
high velocity within that range through interaction with
LASER PHYSICS
Vol. 12
No. 2
2002
PHOTO-IONIZATION OF TWO-ELECTRON ATOMS
the laser. Collisions with the other electron that kick it
out are supposed to be only marginally above the
threshold and thus also do not leave much kinetic
energy in either of the electrons. Low velocity at close
range means also little angular momentum, and in practice the calculations converge to the 1% level with l2
running up to 4 for angle-integrated quantities, and l2
up to 6 for angular distributions. Electron 1, on-the
other hand, will be partly in the ground state, and partly
in the continuum. The quiver motion brings it quickly
far from the nucleus, without the opportunity to spread
much in the angular dimension. The narrow angular
confinement that results makes the use of a large number of angular momenta obligatory anyway, and with this
number of l1 the ground state can be accurately represented even in the velocity gauge. At 800 nm 40–50 angular momenta are usually sufficient to obtain convergence.
A classical estimate of the azimuthal (internal)
angular momentum can be easily made. The maximum
return energy of the first electron is 3.17 times [10] the
ponderomotive energy UP = I/4ω2. At 800 nm and 5 ×
1014 W/cm2 this return energy is 3.48 Hartree. The core
electron is on the average located 0.5 Bohr from the
nucleus, (moving radially with a speed 2), where the
potential energy is –4 Hartree. By the time the returning
electron gets there, its kinetic energy is thus 7.48 Hartree,
corresponding to a velocity 3.87. Kinematics of the
electron–electron collision tell us that only half this
energy can be converted into out-of-plane motion,
divided over two electrons. The maximum velocity perpendicular to the collision plane is thus 2.18 a.u., and at
0.5 Bohr from the nucleus this corresponds to an azimuthal angular momentum 1.09. If the collision would
take place at r = 1 Bohr (the outer turning point of the
1s orbital) the angular momentum can get up to 1.66.
Classically, there is thus no way for the azimuthal angular momentum to get up to 2, and it might thus be
expected that components of the wave function with
m ≥ 2 will hardly contribute.
In agreement with these estimates, practice shows
that only two points in the ϕ dimension are required for
convergence to the 0.1% level. One could view these as
the points ϕ = 0, π or as m = 0, 1, depending on if we
chose an angle or an angular momentum representation. We opted for the latter, although a transform
between the two of them with only two points is so trivial that this hardly matters. More interesting is that with
two ϕ points the problem becomes topologically equivalent to a two-electron atom in a two-dimensional
world. In such a four-dimensional problem the symmetry with respect to the polarization axis leads to an overall parity P, which is a constant of the motion, and an
internal parity p, that takes over the role of m and can
have only two values. After separating off these coordinates, the remaining space consists of the product of
two half-planes, just as in the 3d case after accounting
of the azimuthal variables.
LASER PHYSICS
Vol. 12
No. 2
2002
357
1. ABSORBING BOUNDARY
As the size of the grid has such a large impact on the
resource requirements, a lot can be gained from
exploiting localization of the wave function, not wasting any grid points on regions that the wave function
never enters. An even more aggressive implementation
of this strategy eliminates also regions where the
charge does get, but from which it is not expected to
return to the regions we are interested in. For instance,
an electron can only reach a point much more than two
quiver amplitudes from the atom if it has an outward
cycle-averaged (“drift”) velocity. For a free electron the
drift velocity becomes a constant of the motion, and
such an electron will continue to move outward and
never approach the nucleus again.
Eliminating an area where the electron gets from the
grid introduces an unnatural grid boundary that the
electron will touch. To achieve the intended effect,
namely, the disappearance of the electron into the nonsimulated region, it is essential that no reflections occur
from this artificial boundary. We thus need an absorbing boundary condition.
For the purpose of calculating the ion yields for double ionization, we do not have to follow the second
ejected electron once it becomes certain that it is
removed from the atom. Although in principle one can
never be sure if an electron in the presence of a longrange Coulomb potential will end up into some very
high Rydberg state, the probability for this quickly
decreases with increasing distance. The pragmatic
approach, adopted by many people, is therefore to simply assume that ionization has taken place after the
electron exceeds a certain distance. For instance, a
sphere of 8 Bohr completely contains the (n = 2)- and
(n = 3)-Rydberg states of He2+, so as soon as the second
electron crosses this distance we might consider it ionized.1
To simulate the inner electron only upto r = 8
reduces the grid size enormously: with the fourth-order
radial derivative operator [1] only about 32 points are
needed in this dimension to get an acceptable accuracy
(we aim at about 1%). But at the end of this range we
then need a boundary condition that absorbs outgoing
waves. The commonly practiced way of absorbing the
outgoing flux, by continuing propagation past the point
of interest in the presence of a smoothly growing complex (absorbing) potential until the wave amplitude
becomes negligible, spoils this advantage. To keep
reflections at a bearable level, this method requires and
absorber thickness of dozens of points, driving up the
number of points needed in this dimension by an unacceptably large factor.
1 Even
if the electron would get trapped in higher Rydberg states,
from a philosophical point this would not matter: what we are trying to fathom is how it is the mechanism through which the inner
electron can be kicked out of the ground state at intensities where
the laser alone can not do it, and where the electron ends up is of
minor importance.
358
MULLER
To achieve the desired result, we developed a
method of absorbing outgoing waves that requires only
a few points. The most simple one even reduces the
absorber size to zero points, in the sense that all points
of the grid do faithfully represent the wave function as
it would have been in absence of the absorbing boundary (i.e., if the grid was continued beyond R).
These absorbers are based on the idea that the radial
spacing of the grid is determined by accuracy considerations. As a result the grid is able to support very high
momenta (e.g., a grid spacing of 0.25 Bohr supports
electrons with kinetic energy upto 80 Hartree). Most of
these high momenta propagate very inaccurately, but
this does not matter, since they do not occur in the problem. The low energies that do occur propagate accurately, and as long as the integration scheme does not
generate the high momenta out of nothing (a numerical-stability issue), we don’t care how unrealistically
they propagate. In the same spirit, it is also not important if these high, non-occurring momenta are absorbed
by the boundary, as long as the boundary does not inject
them.
For convenience the following discussion ignores
any potential energies, and assumes the electron
behaves as if completely free in the asymptotic region
where it is absorbed. The zero-point absorber can be
constructed by simply require that the outgoing wave
exp(ik0r), for some selected momentum k0, is an eigenfunction of the Hamiltonian with boundary condition.
The discretized kinetic-energy Hamiltonian is tridiagonal and thus involves only relations between a point and
its neighbors. This makes most points of the grid oblivious of the boundary, and the exponential wave is automatically an eigenfunction, just as in the boundary-less
case. Only in the last row of the Hamiltonian matrix the
boundary is felt, but it is always possible to set the diagonal element on this row to a complex number such that
the outgoing exponential wave satisfies the eigenvalue
equation also in the last point. To absorb outgoing
waves the imaginary part of this lower-right corner then
automatically has the proper sign for making the norm
of the wave function decrease.
The zero-point method gives perfect absorption for
outgoing momentum k0, but the reflection coefficient
rapidly increases for momenta that deviate from this.
For high energies the reflection rises to nearly 1, but we
don’t care about this since such energies will not occur.
For vanishing energy, the reflection will also approach 1.
This is a much bigger nuisance, but unfortunately
unavoidable. The very definition of reflection coefficient (as the ratio of the waves of momentum k and –k)
implies that it will be 1 at k = 0. In a finite system the
reflection coefficient will be a continuous function of
the momentum, so some momenta close to 0 will be
poorly absorbed too. A complex-potential absorber will
also have very poor absorption as soon as the wavelength of the electron becomes larger than the thickness
of the absorber.
The energy window of good absorption of the zeropoint absorber is thus very narrow, and this makes it
only useful in problems where we can make a pretty
good guess about the energies that will hit the boundary. In a tunneling situation this might very well occur
somewhat beyond the exit of the tunnel. We could tune
the absorber to perfectly absorb the tunneling current at
the peak electric field. At lower field the tunneling current would have lower kinetic energy at the same point,
and would thus be absorbed less well. But the tunneling
current decreases so rapidly with field that it would
have effectively turned off before this starts to be a
problem anyway. Nevertheless, the applicability of the
zero-point absorber remains limited: in the nonsequential double-ionization problem the mechanism
that ejects the inner electron is not expected to be tunneling, and this makes it hard to predict the electron
energy at the boundary.
To improve on the zero-point absorber, we simply
extend the philosophy that underlies its design stepwise. For the first upgrade, we thus require exact
absorption of two independently chosen outgoing
momenta k1 and k2. Again such exponential waves satisfy the eigenvalue equation for the tridiagonal operator
automatically everywhere except in the boundary point.
Unfortunately it is not possible to make two different
exponentials eigenfunctions by just adapting the lowest
row of the Hamiltonian. With the same amplitude, outgoing waves of different momenta represent different
currents to be absorbed, while absorption of norm in
the final point only depends on the amplitude in that
point.
The desired result can thus only be achieved by giving the different outgoing waves eigenfunctions that
have an amplitude in the point(s) where the absorption
takes place proportional to the square-root of their
momentum, i.e., the final point can no longer in all
cases represent the wave function. Making the value of
the eigenfunctions in this point an arbitrary dynamical
variable introduces two complex degrees of freedom
(one for each momentum). The point, however, is referenced in the two bottom rows of the Hamiltonian, and
we thus have to satisfy four complex conditions (two
for each momentum). This imposes two complex (or
four real) conditions on the coefficients appearing in
the Hamiltonian on these rows. Maintaining the tridiagonal form, five coefficients are available there, so we
can even impose extra conditions. We elected to keep
the Hamiltonian matrix symmetric, leaving three coefficients to play with. To provide four real degrees of
freedom with three coefficients, only one of them has to
be complex, and we took this to be the lower-right corner. The absorption is then limited to this single point.
The appendix shows how to derive the required coefficients. Usually we pick k1 = 0.6 and k2 = 1.5 causes better than 99% absorption in the energy range 3 to 40 eV.
(Note that at R = 8 the depth of the potential well of a
doubly charged ion is 0.25 Hartree, or 6.8 eV, so that it
LASER PHYSICS
Vol. 12
No. 2
2002
PHOTO-IONIZATION OF TWO-ELECTRON ATOMS
is not evident that in the full problem electrons with an
energy below 3 eV would not reflect back to the
nucleus.)
With the described choices, the only deviation from
hermiticity is the imaginary part of the lower-right corner of the Hamiltonian. Using a split-operator scheme
to implement the time propagation, where we split the
Hamiltonian into the hermitian and anti-hermitian part,
reduces the partial propagator of the latter to a simple
mask multiplication, where the mask differs form 1 in
only a single point! To have as only non-unitary part of
the propagator a multiplication in a single point with a
factor smaller than 1 guarantees that the absorption process can only decrease the norm of the wave function.
The reflection coefficient is thus smaller than 1 for any
outgoing momentum. It is furthermore very easy to
keep track of how much probability is removed that
way from the grid, and the equation of continuity
(which the propagator for the hermitian part obeys)
tells us that this removed probability is a measure for
the current leaving the grid edge.
The absorber discussed above assumed a simple
three-point finite-difference representation Δ2 of the
radial Laplacian. In our calculations, however, we use
2
–1
h
the fourth-order correction ⎛ 1 + ------ Δ 2 ⎞ Δ 2 . The cor⎝
12 ⎠
rection brought about by the first factor is supposed to
be small, so to lowest order the boundary condition in
Δ2 would have to be the one discussed. We simply use
the hermitian part of this operator as Δ2 in the fourthorder operator, which guarantees that the corrected
operator remains hermitian. The anti-hermitian part is
not corrected at all. The correction of the hermitian part
turns out to hardly affect the absorber properties.
2. INITIAL STATE
The initial state for the calculations has to be chosen
with great care. Photo-ionization probabilities can be
extremely small, since they are exponentially sensitive
to the magnitude of the laser field. Double-ionization
probabilities are even smaller. An initial state that is not
the exact ground-state eigenfunction of the propagator
will contain components along other eigenstates of it. If
the propagator is any good as an approximation to reality, such eigenstates are either continuum states
(amongst with the double-ionization continuum), leaving the atom from the very beginning, or excited bound
states that will be completely stripped of the atom at
very low laser fields. With an inaccurate initial state,
such initial bursts of ionization can easily swamp the
true signal if the latter is weak.
In the SAE case, using a half-implicit scheme (1 –
iHt)/(1 + iHt) for the atomic propagator, this problem
is solved by using the ground state of the discretized
atomic Hamiltonian as initial state. Since the propagator is a function of the Hamiltonian, they have the same
LASER PHYSICS
Vol. 12
No. 2
2002
359
eigenfunctions. Note that in another popular way to
solve the SAE problem, where the propagator is split in
a potential- and kinetic-energy part, this method can
not be used, since even in the absence of the field the
propagator is not a function of the atomic Hamiltonian.
In the two-electron case we face a similar problem. The
atomic propagator is split in several parts: one part for
each electron, and one for their Coulomb repulsion.
Due to non-commutation of the two one-electron operators with the repulsion, this is not the same as the
exponent of the sum of these operators. One really has
to determine the eigenfunctions of the propagator.
A popular way to determine ground-state eigenfunctions is imaginary-time propagation. In this
method one integrates the Schrödinger equation in
imaginary time, i.e., one solves ∂tΨ = –HΨ. The solutions of this equation blow up exponentially, the ground
state (having the most negative eigenvalue) the fastest
of all. So after some time all other contributions are
negligible in comparison.
Unfortunately, this method does not solve the problem. When the atomic propagator is split, the splitting
error depends on the time step τ, usually as τ2. With
non-commuting operators this error will not only show
up in the propagator eigenvalues, but also in its eigenfunctions. Making τ imaginary flips the sign of the error
and alters the eigenfunctions. In fact the eigenfunctions
differ twice as much from the desired one as the eigenfunctions of the Hamiltonian itself (which coincide
with those of the propagator for infinitely small τ, since
in this limit the splitting error vanishes).
A simple and reasonably efficient solution to the
problem exists. The numerical solutions generated by
the (split) propagator (in real time) are superpositions
of the evolution of the contributing eigenstates. So in
principle they can be extracted from this evolution by
Fourier analysis. Since it is impractical to generate a
solution for an infinite amount of time, the various frequency components can not be resolved with infinite
precision. If we generate the solution Ψ(t) over a finite
time interval [0, 2T], we can approximately extract the
2T
Fourier component of frequency E as 1/2T 0 Ψ (t)eiEtdt.
∫
Since Ψ(t) = e–iHt Ψ(0), this extraction procedure is formally equivalent to applying the operator sin((E –
H)T)/(E – H)T to the initial state Ψ(0). From this we can
see that any contribution of an eigenstate of energy Ek
is suppressed by a factor sin(Ek – E)T/(Ek – E)T.
If the sought state is well separated (by ΔE) from all
other states in the spectrum, we can choose T large
enough to be sure that there are no other states within
the central lobe of sinc function (ΔET > π), and outside
this lobe the maximum of the sinc function is 0.21. This
means that all contaminations of the state are suppressed by a factor 5. Simply repeating this procedure
a number of times reduces all contaminations exponentially to zero.
360
MULLER
This method is trivial to program, since it can make
use of the ordinary propagator, and the only added part
has to calculate the time-average of the exponentially
weighted wave function. The integration time 2T has to
be on the order of 2π/ΔE, where ΔE is the lowest excitation energy. In fact convergence can be sped up by not
using the same value of T on every iteration, but varying it a little: this spreads the zeros of the sinc function
to reduce the contribution of any contaminating state
even faster.
The eigenvalue equation is
⎛
.
⎜
⎜ exp ( – 2ik i h )
( Δ 2 – Λ i ) ⎜ exp ( – ik i h )
⎜
⎜
1
⎜
⎝
γi
⎞
⎟
⎟
⎟ = 0,
⎟
⎟
⎟
⎠
(7)
where γi is the arbitrary value of the wave function in
2
3. RESULTS AND CONCLUSIONS
In conclusion, we can say that it is possible to solve
important aspects of the two-electron problem even on
a single PC. Using the highly efficient techniques
developed for integrating the single-electron problem
reduces the integration times to about 5 μs per grid
point per time step on a 333 MHz Celeron. Aggressive
trimming of the five-dimensional grid can reduce it to
500 × 40 × 32 × 5 × 2 = 6.4 M points (r1 × ϑ1 × r2 × ϑ2 ×
m), for a modest memory requirement of 52 Mbyte and
an execution speed of 120 steps per hour. With 2000
steps per cycle a short pulse can be simulated in four
days.
the last (auxiliary) point of the grid, and Λi ≈ – k i the
eigenvalue corresponding to momentum ki (i = 1, 2).
We assume everywhere that ki h 1, so that we can
approximate exp(–ikh) by 1 – ikh. As mentioned in the
text only the equations resulting from the last two rows
of the Laplacian matrix are relevant; the others are satisfied automatically by the choice of Λi .
To find the conditions for the Hamiltonian coefficients A, B, and C, we first eliminate γi from the equations using the last row:
2
γ i = B/ ( h Λ i – C ),
(8)
substitution of which in the equations of the forelast
row produces
2
2
2
– ik i h + A + B / ( h Λ i – C ) – h Λ i = 0.
ACKNOWLEDGMENTS
This work is part of the research program of FOM
(Fundamental Research on Matter), which is subsidized by NWO (Netherlands Organization for the
Advancement of Research).
(9)
We then proceed to eliminate C (which, like the γi , is an
arbitrary complex number):
2
2
2
C – h Λ i = B / ( A – ik i h – h Λ i ),
(10)
leaving the single (complex) equation
2
APPENDIX
Obtaining the coefficients for the absorbing termination on the Laplacian is not entirely trivial, due to the
fact that we are dealing with mixed real and complex
conditions and a set of non-linear equations. As a consequence, a brute-force elimination attempt quickly
leads us to unwieldly equations. With a little subtlety
solutions can be easily found analytically, as shown
below.
The simplest finite-difference approximation of the
Laplacian, written in matrix representation on a grid
with spacing h is
⎛ . .
⎞
⎜
⎟
⎜ . –2 1
⎟
–2 ⎜
⎟.
Δ2 = h
1
⎜ 1 –2
⎟
⎜
1 –1 + A B ⎟
⎜
⎟
B C⎠
⎝
h ( Λ2 – Λ1 )
2
1
1
= B ⎛ ------------------------------------– -------------------------------------⎞ .
⎝ A – ik h – h 2 Λ A – ik h – h 2 Λ ⎠
1
1
2
(11)
2
Both A and B follow uniquely from this single complex
equation because we required them to be real; to solve
for them we write separate equations for the real and
imaginary part.
2
2
2
2
( k 1 k 2 h – ( A – h Λ1 ) ( A – h Λ2 ) ) = B ,
2
2
(12)
2
h ( Λ2 – Λ1 ) ( – k 2 ( A – h Λ1 ) – k 1 ( A – h Λ2 ) )
2
= B ( k 1 – k 2 ).
(13)
Eliminating B2, we are left with a quadratic equation
for A, which is easily solved for any momenta k1 and k2
that we care to pick:
(6)
2
2
2 2
A = k 1 k 2 h – k 1 k 2 h ( ( k 1 + k 2 ) h + 1 ).
(14)
From this B2 and C then follow easily by backsubstitution in Eqs. (12) and (10), respectively. (The choice of
LASER PHYSICS
Vol. 12
No. 2
2002
PHOTO-IONIZATION OF TWO-ELECTRON ATOMS
sign of the square root in Eq. (14) guarantees that
B2 > 0.)
REFERENCES
1. Muller, H.G., 1999, Laser Phys., 9, 138.
2. Kulander, K.C., Schafer, K.J., and Krause, J.L., 1992,
Atoms in Intense Laser Fields, Gavrila, M., Ed. (Boston:
Academic), p. 247; Schafer, K.J. and Kulander, K.C.,
1990, Phys. Rev. A, 42, 5798.
3. Muller, H.G. and Kooiman, F.C., 1998, Phys. Rev. Lett.,
81, 1207.
4. Muller, H.G., 1999, Phys. Rev. A, 60, 1341.
5. Nandor, M.J., Walker, M.A., Van Woerkom, L.D., and
Muller, H.G., 1999, Phys. Rev. A, 60, R1771.
6. Muller, H.G., 1999, Phys. Rev. Lett., 83, 3158.
LASER PHYSICS
Vol. 12
No. 2
2002
361
7. Muller, H.G., 2001, Opt. Express, 8, 44.
8. Muller, H.G., 2001, Opt. Express, 8, 86.
9. Yang, B., Schafer, K.J., Walker, B., Kulander, K.C.,
Agostini, P., and DiMauro, L.F., 1993, Phys. Rev. Lett.,
71, 3770.
10. Smyth, E.S., Parker, J.S., and Taylor, K.T., 1998, Comput. Phys. Commun., 114, 4.
11. Muller, H.G., 2000, CP525, Multiphoton Processes:
ICOMP VIII, 8th International Conference, DiMauro, L.F.,
Freeman, R.R., and Kulander, K.C., Eds. (AIP), p. 257.
12. Corkum, P.B., 1993, Phys. Rev. Lett., 71, 1994.
13. Muller, H.G., 2001, Opt. Express, 8, 417.
14. Muller, H.G., 2001, Super-Intense Laser–Atom Physics,
Piraux, B. and Rzazewski, K., Eds. (Netherlands: Kluwer), p. 95.