Position–momentum correlations in electron–ion scattering in strong laser fields

Electron–ion scattering in strong laser fields is of fundamental importance in the process of strong field ionization of atoms [1, 2]. The rescattering of the electron after ionization gives rise to scattering rings in the angular distribution [3] and to a plateau region in the energy distribution [4] in high-order above-threshold ionization (ATI). Laser-driven electron–ion scattering also plays an important role for inverse bremsstrahlung absorption of laser light in plasmas [5–8]. Numerical solutions to the time-dependent Schrödinger equation (TDSE) provide detailed information on the spatial probability distribution as well as the energy and momentum distributions of the scattered electrons [9]. However, while the energy distribution shows regular ATI peaks separated by one-photon energy, the spatial distributions apparently miss a simple physical interpretation. It is only in some cases that spatial wavepackets have been successfully identified with classical electron orbits [10].

In order to gain a better understanding of the position–momentum correlations of the scattered electrons, it is worth studying the evolution of the scattering problem in the joint phase space. Although phase-space trajectories are excluded in quantum mechanics by the Heisenberg uncertainty relation, a generalized phase-space representation of a quantum system can be given in terms of the Wigner function [11]. The Wigner function represents a quasi-probability distribution function, assuming positive and negative values, which approaches the usual phase-space probability distribution in the classical limit. In the context of strong-field laser–atom interactions, the Wigner function has apparently been studied only occasionally to visualize, e.g., stabilization [12], tunneling [13] and double ionization [14].

In this work, we calculate the Wigner function numerically from the underlying quantum Liouville equation to study its evolution during the scattering of an incident electron beam by an ion over several periods of the applied laser field. We focus our attention on the scattering by a 1d soft-core Coulomb potential in the presence of a linearly polarized laser field. Basic strong-field features of ionization [15, 16] and scattering [17, 18] have already been studied in this framework. Here we generalize previous work by considering the evolution in 2d phase space based on a fast Fourier transform (FFT) split operator method. It will be shown that the phase-space trajectories of the scattered electrons indicate both classical and nonclassical features and that the observed position–momentum correlations provide further insight into the formation of ATI spectra.

The evolution of the statistical operator ρ(t) of a quantum system with the Hamilton operator H is governed by the von Neumann equation

$$\begin{equation} \mathrm{i}\partial_t\rho=[H,\rho]. \end{equation}\noindent \tag{ 1 }$$

Here and in the following, atomic units are used, where the electron mass, the elementary charge and the Planck constant are equal to unity (m = e = ℏ = 1). We consider a 1d scattering model with the Hamiltonian

$$\begin{equation} H=-\textstyle\frac{1}{2}\Delta+U(x). \end{equation}\noindent \tag{ 2 }$$

The scattering potential is taken as a soft-core Coulomb potential

$$\begin{equation} U(x)=-\frac{1}{\sqrt{x^2+\epsilon^2}}. \end{equation}\noindent \tag{ 3 }$$

In ionization problems, the parameter is often taken in the range ⪆1. For  = 1.4 the ground state energy of the model potential assumes approximately the value −0.5 of hydrogen. However, for typical impact velocities of the order of unity the potential varies relatively slowly and its reflectivity therefore is quite small. In scattering problems, it is therefore useful to consider the parameter regime $\epsilon \lessapprox 1$. In the present calculations we have set  = 0.1, which accounts for about 10% reflection at impact velocities of the order of unity [17].

The interaction with the laser field is described in the Kramers–Henneberger (KH) framework [19]. In an oscillating laser field, ${\cal E}(t)={\cal E}_0\sin (\omega t)$, an electron is subject to the quiver motion

$$\begin{equation} \xi(t)=\xi_0\sin(\omega t),\quad \xi_0=\frac{v_0}{\omega},\quad v_0=\frac{{\cal E}_0}{\omega}. \end{equation}\noindent \tag{ 4 }$$

In the KH frame, the quiver motion is removed and, instead, the scattering potential in (2) is oscillating as U = U(x + ξ(t)). The KH frame is convenient, since the asymptotic velocities for |x| → ∞ correspond to the drift velocities of the particle measured outside the laser field.

The initial wavefunction for an incident electron is chosen as a Gaussian wavepacket,

$$\begin{equation} \psi(x)=\frac{1}{\sqrt{\sqrt{2\pi}\sigma}}~\mathrm{e}^{-\frac{(x-x_0)^2}{4\sigma^2}}\,\mathrm{e}^{\mathrm{i}p_0(x-x_0)}. \end{equation}\noindent \tag{ 5 }$$

The probability distribution of x is w(x) = |ψ(x)|². It is normalized to unity. The expectation values for position and momentum are x₀ and p₀, respectively, and the standard deviation of x from x₀ is σ.

The Wigner function is defined in terms of the density matrix 〈x|ρ|x'〉 by a coordinate transformation x = q + r/2, x' = q − r/2 and a Fourier transformation with respect to the difference coordinate r,

$$\begin{equation} f(p,q)=\int_{-\infty}^{+\infty} \mathrm{d}r~\left\langle q+\frac{r}{2}|\rho|q-\frac{r}{2}\right\rangle \,\mathrm{e}^{-\mathrm{i}pr}. \end{equation}\noindent \tag{ 6 }$$

By definition, the Wigner function is real; however, in contrast to a classical phase-space distribution function it can assume negative values. The Wigner function of (5) is a Gaussian probability distribution with respect to both p and q,

$$\begin{equation} f(p,q)=2\,\mathrm{e}^{-\frac{(q-x_0)^2}{2\sigma^2}}\mathrm{e}^{-2\sigma^2(p-p_0)^2}, \end{equation}\noindent \tag{ 7 }$$

being normalized according to

$$\begin{equation} \int_{-\infty}^{\infty}\frac{\mathrm{d}p}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}q ~f(q,p)=1. \end{equation}\noindent \tag{ 8 }$$

According to (1) and (6), the Wigner function satisfies the quantum Liouville equation

$$\begin{equation} \partial_t f(p,q,t) +{\cal L} f(p,q,t)=0, \end{equation}\noindent \tag{ 9 }$$

with the Liouville operator

$$\begin{equation} {\cal L}=p\partial_q+\mathrm{i}\left[U\left(q+\frac{\mathrm{i}}{2}\partial_p\right)-U\left(q-\frac{\mathrm{i}}{2}\partial_p\right)\right]. \end{equation}\noindent \tag{ 10 }$$

The Liouville equation of classical mechanics is obtained by expanding the potentials up to linear order around q.

The Liouville equation (9) has been solved numerically by an FFT-split-operator method similar to those used in [20–22]. For a sufficiently small time step, Δt = tⁿ⁺¹ − tⁿ, the Wigner function can be propagated by the time-evolution operator

$$\begin{equation} U^{n+1/2}=\mathrm{e}^{- {\cal L}^{n+1/2}\Delta t},\quad {\cal L}^{n+1/2}={\cal L}(p,q,t+\Delta t/2). \end{equation} \tag{ 11 }$$

The Liouville operator can be split into parts, ${\cal L}={\cal L}_q + {\cal L}_p$, which include only position or momentum derivatives,

$$\begin{equation} {\cal L}_q=p\partial_q,\quad {\cal L}_p=\mathrm{i}\left[U\left(q+\frac{\mathrm{i}}{2}\partial_p\right)-U\left(q-\frac{\mathrm{i}}{2}\partial_p\right)\right]. \end{equation} \tag{ 12 }$$

Up to second order in Δt, the time-evolution operator can then be factorized in the form

$$\begin{equation} U=U_q(\Delta t/2)U_p(\Delta t) U_q(\Delta t/2) \end{equation} \tag{ 13 }$$

with

$$\begin{eqnarray*} &&U_q(\alpha)= \mathrm{e}^{- {\cal L}_q\alpha},\quad U_p(\alpha)= \mathrm{e}^{- {\cal L}_p\alpha}. \end{eqnarray*}$$

Applying alternating Fourier transforms with respect to q and p, each factor of (13) can be reduced to a simple multiplication operator. For N time steps it is sufficient to perform first one half-step with U_q(Δt/2), then N full-steps with the product U_q(Δt)U_p(Δt) and finally one half-step backwards with U_q(−Δt/2). In our calculations the time step was chosen as Δt = 0.01 and phase-space grids with 2048 or 4096 points in each direction have been used. The separation of grid points was approximately 0.1 in the q-direction and 0.01 in the p-direction, which proved to be sufficient for resolving all interference fringes with appreciable visibility.

To gain a physical understanding of the Wigner function of the scattered wavepacket, it is helpful to first consider simple superpositions of quantum states and their interference patterns in phase space. Let |ψ〉 = c₁|k₁〉 + c₂|k₂〉 be a superposition of two momentum eigenstates. The wavefunction in the position representation is

$$\begin{equation} \psi(x)=c_1\,\mathrm{e}^{\mathrm{i}k_1x}+c_2\,\mathrm{e}^{\mathrm{i}\int k_2x} \end{equation} \tag{ 14 }$$

and the Wigner function of the corresponding density matrix 〈x|ρ|x'〉 = ψ(x)ψ*(x') is found to be

$$\begin{equation} \fl \frac{ f(p,q)}{2\pi}=|c_1|^2\delta(p-k_1)+|c_2|^2\delta(p-k_2)+2\Re\{c_1c_2\,\mathrm{e}^{\mathrm{i}(k_1-k_2)q}\}\delta\left(p-\frac{k_1+k_2}{2}\right). \end{equation}\noindent \tag{ 15 }$$

The first two terms correspond to classical streams with definite momenta k₁ and k₂. The third term is a nonclassical contribution that corresponds to the quantum interference in the spatial probability distribution. It varies along the q-direction with the difference wavenumber k₁ − k₂ and it propagates with the mean wavenumber (k₁ + k₂)/2. This example shows nicely how the Wigner function contains both classical particle properties and nonclassical wave properties being resolved in phase space by their different momenta.

In the semiclassical framework, the plane waves can be generalized with the Wentzel-Kramers-Brillouin (WKB) approximation to WKB waves with local wavenumbers k_1,2(q). A pair of branches can be continued up to a common turning point where they merge with a vertical slope. The generic behavior of the Wigner function of a superposition of momentum states is two branches originating from a turning point with vertical slope and enclosing on the concave side an interference pattern with the local wavenumber k₁(q) − k₂(q) centered at the mean momentum (k₁(q) + k₂(q))/2. For instance, figure 1 shows the Wigner function of the wavefunction

$$\begin{equation} \psi(x)=\mathrm{e}^{\mathrm{i}\int_0^x \mathrm{d}x' k_1(x')}-\mathrm{e}^{\mathrm{i}\int_0^x \mathrm{d}x' k_2(x')},\quad \end{equation} \tag{ 16 }$$

with local wavenumbers $k_{1,2}(x)=\pm \sqrt {-x}$ in the half-space x < 0. It represents an approximation to the stationary state with energy E = 0 in a linearly increasing potential V (x) = x/2. For simplicity the amplitudes are taken constant and a constant phase π/4 has been neglected. The well-known exact wavefunction can be expressed in terms of the Airy function Ai(x).

Zoom In Zoom Out Reset image size

Another important example is a superposition of two position eigenstates |ψ〉 = c₁|x₁〉 + c₂|x₂〉. A position eigenstate is a simple model of a particle localized at a definite position with a corresponding wide spread in momentum. In the momentum representation, the wavefunction becomes

$$\begin{equation} \psi(k)=c_1\mathrm{e}^{-\mathrm{i}kx_1}+c_2\mathrm{e}^{-\mathrm{i}kx_2}. \end{equation} \tag{ 17 }$$

The Wigner function can now be calculated from the density matrix in the momentum representation 〈k|ρ|k'〉 = ψ(k)ψ*(k') by setting k = p + s/2, k' = p − s/2 and taking the inverse Fourier transform

$$\begin{eqnarray} f(p,q)&=&\int_{-\infty}^{\infty} \frac{\mathrm{d}s}{2\pi}~\left\langle p+\frac{s}{2}|\rho|p-\frac{s}{2}\right\rangle \,\mathrm{e}^{\mathrm{i}sq}\nonumber\\ &=&|c_1|^2\delta(q-x_1)+|c_2|^2\delta(q-x_2) +2\Re\{c_1c_2\,\mathrm{e}^{\mathrm{i}p(x_1-x_2)}\}\delta\left(q-\frac{x_1+x_2}{2}\right).\nonumber\\ \end{eqnarray}\noindent \tag{ 18 }$$

It consists of the populations at the positions x_1,2 and an interference pattern at the mean position (x₁ + x₂)/2 varying along the momentum direction with the wavenumber x₁ − x₂. Generalizing again the plane waves to the WKB waves, two branches with local wavenumbers x_1,2(p) merging at a turning point with horizontal slope will enclose an interference pattern with local wavenumber x₁(p) − x₂(p) along the p coordinate. Schematically, the Wigner function has the same form as in figure 1 with q and p interchanged. The interference pattern along the p-direction is of particular interest, since it determines directly the momentum and energy spectra.

Basic properties of ATI spectra can simply be understood from this picture. First, consider jets of particles that originate from scattering at the ion in subsequent periods of the laser field. The particles with a given momentum p will be found at time t at the positions x_n(p) = p(t − nT), where n = 0,1,2,... denotes the laser cycle of the emission and T = 2π/ω the laser period. The interference pattern between neighboring branches x_n(p) and x_n+1(p) is directed along the p-direction and has the wavelength Δp = 2π/(x_n − x_n+1) = ω/p. In the energy spectrum, corresponding peaks occur at the distance

$$\begin{equation} \Delta E=p\Delta p=\omega. \end{equation}\noindent \tag{ 19 }$$

As a conclusion of this simple argument, one obtains the familiar ATI peaks, separated by the photon energy.

There is another interference pattern in ATI spectra that is associated with the particles scattered within one period. Neglecting for simplicity the drift motion, the backscattered particles are emitted at time t_s with the momentum p = −2v₀ cos(ωt_s) at the position x = ξ₀ sin(ωt_s) in the KH frame. The phase-space trajectory of the ensemble of emitted particles is given by

$$\begin{equation} \frac{[x-v(t-t_{\mathrm{s}})]^2}{\xi_0^2}+\frac{v^2}{(2v_0)^2}=1. \end{equation}\noindent \tag{ 20 }$$

On the x-axis, v = 0, the ellipse has the width 2ξ₀. From this width one can estimate the period of the momentum spectrum inside the ellipse as Δp = 2π/(2ξ₀) = π/ξ. Accordingly, the ATI peaks in the energy spectrum will be modulated with a period of approximately

$$\begin{equation} \Delta E=\pi\frac{ p}{\xi}=\pi\frac{\sqrt{2E}}{\xi}. \end{equation}\noindent \tag{ 21 }$$

In contrast to the constant separation of the ATI peaks (19), the modulational wavelength (21) depends on the quiver amplitude and increases with momentum.

Finally, we consider the Wigner function for a superposition |ψ〉 = c₁|x₀〉 + c₂|p₀〉 of position and momentum eigenstates. It describes the interference of free and bound parts of the wavefunction. The density matrix is given by

$$\begin{eqnarray} &&\fl \langle x|\rho|x'\rangle=|c_1|^2\delta(x-x_0)\delta(x'-x_0)+|c_2|^2\,\mathrm{e}^{\mathrm{i}p_0(x-x')}+c_1c_2^*\delta(x-x_0)\,\mathrm{e}^{-\mathrm{i}p_0x'} +c_1^*c_2\delta(x'-x_0)\,\mathrm{e}^{\mathrm{i}p_0x}.\nonumber\\ \end{eqnarray} \tag{ 22 }$$

Taking the Fourier transform with respect to the difference coordinate r = x − x' yields the Wigner function

$$\begin{equation} \fl f(p,q) = |c_1|^2\delta(q-x_0)+2\pi|c_2|^2\delta(p-p_0) +4\Re\{c_1c_2^*\,\mathrm{e}^{2\mathrm{i}(p-p_0)(q-x_0)-\mathrm{i}p_0x_0}\}. \end{equation} \tag{ 23 }$$

It shows hyperbolic interference fringes n = 0,1,2,... at the positions (p − p₀)(q − x₀) = nπ + const.

We first focus our attention on the scattering of a Gaussian wavepacket (x₀ = −50, p₀ = 1, $\sqrt {2}\sigma =20$) by a soft-core Coulomb potential in the absence of a laser field. Looking at this simple example, one can clearly recognize the classical and nonclassical features of the Wigner phase-space representation. In figure 2, the Wigner function is presented at different stages of the scattering process. Red color corresponds to positive and blue color to negative values of the quasi-probability. At time t = 10, the initial distribution has become slightly sheared due to the different distances moved by particles with different momenta, but it appears closely classical. At time t = 30, the center of the wavepacket has moved to x = −20 and the front has reached the scattering center at x = 0. Within the interaction region, the distribution becomes nonclassical. The population localized within the range of the potential interferes with the incoming beam, leading to hyperbolic fringes as in (23). At time t = 50, the center of the Gaussian has reached the scattering center. In the upper half-space (p > 0) one can recognize the quasi-classical behavior of particles crossing the potential with a given energy, thereby gaining momentum near the scattering center. In addition, nonclassical hyperbolic fringes arise from the incoming and transmitted beams interfering with the localized distribution around x = 0. In the lower half-space (p < 0) the distribution is dominant in the region q < 0 of the reflected wave. Along the negative q-axis one can recognize the interference pattern between the incoming and reflected waves in accordance with (15). At later times (t = 90), the incident packet has passed the scattering center and becomes spatially separated from the reflected one. Within the scattering region the distribution becomes rapidly oscillating in both the p- and q-directions. Due to the rapidly varying phase, the projected position and momentum probability distributions average to zero in this region, although the amplitudes of the Wigner function are still appreciable at the final time. This suggests that the spatial wavepackets decay after the interaction predominately by phase averaging of the Wigner function along the momentum direction.

Zoom In Zoom Out Reset image size

We now consider laser-driven scattering in a laser field with amplitude ${\cal E}_0=0.2$ and frequency ω = 0.2. The impact velocity (p₀ = 2) is chosen larger than the quiver velocity (v₀ = 1) to prevent the minimum collision velocity p₀ − v₀ from becoming small or negative. Furthermore, the laser period is larger than the typical duration of the collision (ω/v₀ = 0.02 ≪ 1). These parameters correspond to the regime of instantaneous collisions [17]. If the duration of the interaction is much smaller than the laser period, the laser field is approximately fixed during the collision. Therefore the basic effect of the laser field consists in the modulation of the impact velocities of the Coulomb collision. From instantaneous collision models one expects backscattered particles with momenta p = −p₀ − 2v₀ cos(ωt_s), depending on the time t_s of the collision.

The evolution of the Wigner function of a Gaussian wavepacket (x₀ = −75, p₀ = 2, $\sqrt {2}\sigma =30$) is shown in figure 3. The period of the laser field is T = 2π/ω ≈ 31. Within the first period the magnitude of the velocity of the scattered particles decreases from p₀ + 2v₀ = 4 to p₀ − 2v₀ = 0 (t = 10) and increases again to 4 (t = 20,30). The corresponding phase-space curve is indicated in the figure as the classical path of backward scattered particles. The two branches formed within the first and the second half-cycle interfere. One can already recognize transient interference fringes (denoted as momentum maxima) that evolve later to the ATI spectrum (19). Within the second period (t = 40, 50, 60) this process is repeated. The modulation of the phase-space curve by the quiver motion as indicated schematically in (21) now becomes clearly visible (denoted as modulation maxima). Within the area enclosed there appear the interference fringes that cause the modulation of the energy spectrum. After the third period (t = 90: left panel) one can observe three emitted classical jets with corresponding areas of nonclassical interference patterns. In the forward direction, the beam is periodically bunched by accelerated and decelerated particles (t = 90: right panel). The interference pattern (momentum maxima) produced between subsequent bunches corresponds to the ATI peaks in the forward peak of the energy spectrum (figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The energy spectrum has been calculated from the momentum presentation of the wavefunction $\hat {\psi }(k)$. The kinetic energies E(k) = k²/2 occur with the probabilities $F(E(k))=|\hat {\psi }(k)|^2/k$. Formally, the energy is taken positive for k > 0 and negative for k < 0 to distinguish between the spectra of forward and backward scattered particles, respectively. Energy spectra corresponding to the parameters of figure 3 are shown in figure 4. The transient spectrum after the first period (t = 30) shows peaks arising from a single jet of emitted particles interfering with trapped particles. The modulation is not yet present at this time and the spectrum decreases nearly exponentially. After another half-period (t = 40), a second jet of fast particles has formed and the spectrum now develops a plateau region. In the final stage (t = 100), the plateau clearly shows ATI peaks with a peak separation of ω = 0.2 and modulations by the quiver motion. The wavelength of the modulation according to (21) is $\sqrt {2E}\pi /\xi _0\approx 2.5$ around the maximum energy E ≈ 8. This is a rough estimate of the wavelength seen in figure 4; however, it is clearly observed that the wavelength decreases with decreasing energy.

In conclusion, the temporal evolution of the Wigner function shows nicely both the periodic emission of scattered particles and the associated formation of the momentum and energy distributions. In general, the spatial probability distribution at a fixed position contains contributions of electrons from different emission cycles. It is therefore generally not possible to identify spatial interference maxima with particles of a definite energy. In contrast, in phase space both the classical trajectories and the nonclassical interference structures can be well resolved. In more realistic two- and three-dimensional ionization and scattering problems, a full calculation of the Wigner function on a grid for every time step may be too demanding. In these cases, it will be more advantageous to propagate representative wavefunctions of an ensemble by the TDSE and to calculate the Wigner function of the ensemble only at a restricted sequence of times. Such procedures have recently been demonstrated for the density matrix of a quantum plasma [23].

The author thanks G M Fraiman for helpful discussions on the Wigner representation of electron–ion scattering. This work was funded by the Marie Curie Grant, project number 230777.