Prerequisite Numeric TDSE solver for Helium
, Coulomb wave function
, Identical particles
Basic data of helium atom
- parahelium: singlet spin, symmetric wave function, there is a 1s ground state (which we use). orthohelium: triplet spin, anti-symmetric wave function, there is no 1s ground state, but the energy of other bound states is relatively low.
- ground state energy (numerical calculation): -2.9037244au (-79.014366 eV)
- 1s2s energy: -2.14597, 1s3s energy: -2.06127, 1s4s energy: -2.03359 (refer to Feist graduation thesis)
- Helium ion Ground state energy: -2au (-54.422772eV) First excited state -0.5au (-13.605693eV) Second excited state -2/9=-0.22...au (-6.0469747eV)
- direct (n=1) Single ionization threshold: 0.9037244au (24.591594eV)
- shakeup (n=2) Single ionization threshold: 2.4037244au (65.408673eV)
- shakeup (n=3) Single ionization threshold: 2.68150218au (72.967391eV)
Fig. 2：Energy Level (Scri98)
90eV XUV ionization
- 90eV's xuv photon energy is 3.3074390au
- 90eV xuv direct (n=1) after ionization, the electron’s final kinetic energy is 2.4037146au (65.408405eV), and the final momentum is 2.1925850au
- 90eV xuv's shakeup (n=2) after ionization, the electron's end kinetic energy is 0.9037146au (24.591326eV), and the end momentum is 1.3444066au
- 90eV xuv's shakeup (n=3) after ionization, the electron's end kinetic energy is 0.62593682au (17.032608eV), and the end momentum is 1.1188716au
Energy
After the ground state is obtained by the virtual time method, we can use the average energy to calculate the ground state energy
\begin{equation}
E = \left\langle \Psi \middle| H_0 \middle| \Psi \right\rangle
\end{equation}
Momentum spectrum (6 dimensions)
The double ionization is considered here, because we are going to project both electrons to the continuous state. If any electron is in the ground state, then the projection will be zero. The scattering state of the hydrogen atom is the Coulomb wave function , while the (two-electron) scattering state of the helium atom has no analytical solution. Therefore, it is assumed that when the time is large enough, there is no interaction between the two electrons, you can directly combine the two The Coulomb wave function can be tensor product and then symmetric/antisymmetric. In this way, the scattering state and the bound state are obviously not orthogonal, so all the bound states need to be removed in the analysis.1. Aihua's method is to directly subtract the projection of the ground state, but in this way the component of the excited state still exists. I think that if the ionized wave packet and the unionized wave packet have been separated, the wave function within a certain radius can be directly mined like a hydrogen atom.
The Coulomb wave function (see eq. 11 , we only need the plane wave to exit) is ($\eta = Z/k$, $Z = -2$)
\begin{equation}
\psi_{ \boldsymbol{\mathbf{k}} }( \boldsymbol{\mathbf{r}} ) = \sqrt{\frac{2}{\pi}}\sum_{l, m} \frac{ \mathrm{i} ^l}{r k} \mathrm{e} ^{- \mathrm{i} \sigma_l(\eta)} F_l (\eta, kr) Y_{l, m}^* ( \hat{\boldsymbol{\mathbf{k}}} ) Y_{l, m}( \hat{\boldsymbol{\mathbf{r}}} )
\end{equation}
Let the symmetrization operator be
\begin{equation}
\hat{S} f( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) = \frac{1}{\sqrt{2}} \left[f( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) + f( \boldsymbol{\mathbf{r}} _2, \boldsymbol{\mathbf{r}} _1) \right]
\end{equation}
Then the symmetrized scattering state is
\begin{equation} \begin{aligned}
& \qquad \psi_{ \boldsymbol{\mathbf{k}} _1, \boldsymbol{\mathbf{k}} _2}( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) = \hat{S} \left[\psi_{ \boldsymbol{\mathbf{k}} _1}( \boldsymbol{\mathbf{r}} _1)\psi_{ \boldsymbol{\mathbf{k}} _2}( \boldsymbol{\mathbf{r}} _2) \right] \\
&= \hat{S} \frac{2}{\pi} \sum_{l_1, m_1} \sum_{l_2, m_2} \frac{ \mathrm{i} ^{l_1 + l_2}}{r_1 r_2 k_1 k_2} \times\\
&\qquad [( \mathrm{e} ^{- \mathrm{i} [\sigma_{l_1}(\eta_1) + \sigma_{l_2} (\eta_2)]} F_{l_1} (\eta_1, k_1 r_1) F_{l_2}(\eta_2, k_2 r_2) Y_{l_1, m_1}^* ( \hat{\boldsymbol{\mathbf{k}}} _1) Y_{l_2, m_2}^* ( \hat{\boldsymbol{\mathbf{k}}} _2))]\times\\
& \qquad Y_{l_1, m_1} ( \hat{\boldsymbol{\mathbf{r}}} _1) Y_{l_2, m_2} ( \hat{\boldsymbol{\mathbf{r}}} _2)
\end{aligned} \end{equation}
The wave function (assuming that the exchange symmetry has been satisfied) is
\begin{equation} \begin{aligned}
&\Psi( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) = \sum_{L, M, l_1, l_2} \frac{1}{r_1 r_2} \psi_{l_1, l_2}^{L, M}(r_1, r_2)\mathcal{Y}_{l_1, l_2}^{L, M}( \hat{\boldsymbol{\mathbf{r}}} _1, \hat{\boldsymbol{\mathbf{r}}} _2)\\
& \quad = \sum_{L, M, l_1, l_2} \frac{1}{r_1 r_2} \psi_{l_1, l_2}^{L, M}(r_1, r_2) \sum_{m_1, m_2} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} Y_{l_1, m_1} ( \hat{\boldsymbol{\mathbf{r}}} _1) Y_{l_2, m_2} ( \hat{\boldsymbol{\mathbf{r}}} _2)
\end{aligned} \end{equation}
Projected (the process is omitted, the symmetrization operator becomes $\sqrt{2}$, two $Y_{l_1, m_1} ( \hat{\boldsymbol{\mathbf{r}}} _1) Y_{l_2, m_2} ( \hat{\boldsymbol{\mathbf{r}}} _2)$ integrals will produce two $\delta$ functions)
\begin{equation}
\left\langle \psi_{ \boldsymbol{\mathbf{k}} _1, \boldsymbol{\mathbf{k}} _2} \middle| \Psi \right\rangle = \frac{2\sqrt 2}{\pi k_1 k_2} \sum_{L, M, l_1, l_2} (- \mathrm{i} )^{l_1 + l_2} \mathrm{e} ^{ \mathrm{i} [\sigma_{l_1}(\eta_1) + \sigma_{l_2} (\eta_2)]} I_{l_1, l_2}^{L, M}(k_1, k_2) \mathcal{Y}_{l_1, l_2}^{L, M}( \hat{\boldsymbol{\mathbf{k}}} _1, \hat{\boldsymbol{\mathbf{k}}} _2)
\end{equation}
among them
\begin{equation}
I_{l_1, l_2}^{L, M}(k_1, k_2) = \int_0^\infty \int_0^\infty F_{l_1} (\eta_1, k_1 r_1) F_{l_2}(\eta_2, k_2 r_2) \psi_{l_1, l_2}^{L, M}(r_1, r_2) \,\mathrm{d}{r_1} \,\mathrm{d}{r_2}
\end{equation}
The square of the projection is the momentum distribution
\begin{equation}
P( \boldsymbol{\mathbf{k}} _1, \boldsymbol{\mathbf{k}} _2) = \left\lvert \left\langle \psi_{ \boldsymbol{\mathbf{k}} _1, \boldsymbol{\mathbf{k}} _2} \middle| \Psi \right\rangle \right\rvert ^2
\end{equation}
If you fix the directions of the two vectors, you can draw a two-dimensional momentum/energy distribution diagram (note that you need to multiply $k_1^2k_2^2$).
Radial probability distribution
\begin{equation}
P(r_1, r_2) = 2\iint \left\lvert \Psi( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) \right\rvert ^2 r_1^2 r_2^2 \,\mathrm{d}{\Omega_1} \,\mathrm{d}{\Omega_2}
\end{equation}
Analogous to the case of hydrogen atoms (
eq. 15 ), let $\lambda$ represent $l_1, l_2, L, M$
\begin{equation}
\begin{aligned}
& \quad \iint P(r_1, r_2) \,\mathrm{d}{r_1} \,\mathrm{d}{r_2} \\
&= \iint 2\iint \left\lvert \frac{1}{r_1 r_2} \sum_\lambda \psi_\lambda(r_1, r_2) \mathcal{Y}_\lambda ( \hat{\boldsymbol{\mathbf{r}}} _1, \hat{\boldsymbol{\mathbf{r}}} _2) \right\rvert ^2 \,\mathrm{d}{\Omega_1} \,\mathrm{d}{\Omega_2} r_1^2 r_2^2 \,\mathrm{d}{r_1} \,\mathrm{d}{r_2} \\
&= \iint 2\sum_\lambda \left\lvert \psi_\lambda(r_1, r_2) \right\rvert ^2 \,\mathrm{d}{r_1} \,\mathrm{d}{r_2} = 2
\end{aligned} \end{equation}
and so
\begin{equation}
P(r_1, r_2) = 2\sum_\lambda \left\lvert \psi_\lambda(r_1, r_2) \right\rvert ^2
\end{equation}
This is the total modulus length in the $(r_1, r_2)$ subspace, which is equal to the sum of the modulus squares of the coefficients of each basis in the subspace.
Energy distribution
Continuing to consider double ionization first, the probability in eq. 8 is normalized to (because the integral is repeated for the whole space)
\begin{equation}
\int P( \boldsymbol{\mathbf{k}} _1, \boldsymbol{\mathbf{k}} _2) \cdot k_1^2 \,\mathrm{d}{\Omega_{k1}} \,\mathrm{d}{k_1} \cdot k_2^2 \,\mathrm{d}{\Omega_{k2}} \,\mathrm{d}{k_2} = 2
\end{equation}
So the distribution function of the absolute value of momentum is
\begin{equation}
P(k_1, k_2) = k_1^2 k_2^2 \int P( \boldsymbol{\mathbf{k}} _1, \boldsymbol{\mathbf{k}} _2) \,\mathrm{d}{\Omega_{k1}} \,\mathrm{d}{\Omega_{k2}}
\end{equation}
The energy distribution function is
\begin{equation}
P(E_1, E_2) = \frac{1}{k_1 k_2} P(k_1, k_2)
\end{equation}
So just substitute
eq. 8 into
eq. 13 to get
\begin{equation}
P(k_1, k_2) = \frac{8}{\pi^2} \sum_{L, M, l_1, l_2} \left\lvert I_{l_1, l_2}^{L, M} (k_1, k_2) \right\rvert ^2
\end{equation}
In fact, the derivation of this formula can directly use spherical waves instead of projecting to plane waves and then do angular integration. The process is analogous to hydrogen atoms
.
Double ionization probability and cross section
total probability of double ionization is
\begin{equation}
P_{DI} = \frac{1}{2}\int P(k_1, k_2) \,\mathrm{d}{k_1} \,\mathrm{d}{k_2} = \frac{1}{2}\int P(k_1, k_2) \,\mathrm{d}{k_1} \,\mathrm{d}{k_2}
\end{equation}
There is an assumption here that the tensor product space of two Coulomb functions is the eigenstate space of double ionization, that is, the space where the two electron Coulomb functions are located. This should be true from physical intuition, but I don't know how to prove it.
Let the average number of incident photons per unit area in the entire wave packet be $n$, and the total cross-section of double ionization $\sigma$ is defined as ($\sigma n$ equals the total probability of double ionization)
\begin{equation}
\sigma = \frac{P_{DI}}{n}
\end{equation}
JAD(joint angular distribution)
The JAD here is not the total directional distribution, but the directional distribution of a certain energy division, namely $k_2 - k_1 = \Delta k$, so there must be a delta function in the integral.
\begin{equation}
\begin{aligned}
P_{JAD}( \hat{\boldsymbol{\mathbf{k}}} _1, \hat{\boldsymbol{\mathbf{k}}} _2) &= \int_0^\infty \int_0^\infty P( \boldsymbol{\mathbf{k}} _1, \boldsymbol{\mathbf{k}} _2) k_1^2 k_2^2 \delta(k_1 - k_2 + \Delta k) \,\mathrm{d}{k_1} \,\mathrm{d}{k_2} \\
&= \frac{8}{\pi^2}\int_0^\infty \int_0^\infty \left\lvert f(k_1, k_2, \hat{\boldsymbol{\mathbf{k}}} _1, \hat{\boldsymbol{\mathbf{k}}} _2) \right\rvert ^2 \delta(k_1 - k_2 + \Delta k) \,\mathrm{d}{k_1} \,\mathrm{d}{k_2} \\
&= \frac{8}{\pi^2} \int_0^\infty \left\lvert f(k_1, k_1 + \Delta k, \hat{\boldsymbol{\mathbf{k}}} _1, \hat{\boldsymbol{\mathbf{k}}} _2) \right\rvert ^2 \,\mathrm{d}{k_1}
\end{aligned}
\end{equation}
Where $f$ is the sum in
eq. 6
\begin{equation}
f(k_1, k_2, \hat{\boldsymbol{\mathbf{k}}} _1, \hat{\boldsymbol{\mathbf{k}}} _2) = \sum_{L, M, l_1, l_2} (- \mathrm{i} )^{l_1 + l_2} \mathrm{e} ^{ \mathrm{i} [\sigma_{l_1}(\eta_1) + \sigma_{l_2} (\eta_2)]} I_{l_1, l_2}^{L, M}(k_1, k_2) \mathcal{Y}_{l_1, l_2}^{L, M}( \hat{\boldsymbol{\mathbf{k}}} _1, \hat{\boldsymbol{\mathbf{k}}} _2)
\end{equation}
Generally, $ \hat{\boldsymbol{\mathbf{k}}} _1$ and $ \hat{\boldsymbol{\mathbf{k}}} _2$ are taken on the same plane ($\phi_1 = \phi_2$ is a constant, such as 0), so that a two-dimensional distribution map $P_{JAD}(\theta_1, \theta_2)$ with the horizontal and vertical coordinates of $\theta_1$ and $\theta_2$ can be drawn.
In the paper, a parameter $\varepsilon$ is commonly used to measure the energy division, which is defined as the ratio of the smaller photoelectron end kinetic energy to the total end kinetic energy
\begin{equation}
\varepsilon = \frac{\min \left\{E_1, E_2 \right\} }{E_1 + E_2} \qquad (0 \leqslant \varepsilon \leqslant 0.5)
\end{equation}
Understanding JAD from the perspective of photons: The total energy is $E = k_1^2/2 + k_2^2/2$, plus $k_2 - k_1 = \Delta k$ to determine the energy of the two photons, and the remaining degrees of freedom are only angles. But in fact, it can be seen from the figure of $P(k_1, k_2)$ that the arc represented by $E = k_1^2/2 + k_2^2/2$ has a width, so it is equivalent to averaging this width.
TDCS(triple differential cross section)
The definition of the analog differential cross section, the double differential cross section can be defined as
\begin{equation}
\frac{\partial^2\sigma}{\partial\Omega_1\partial\Omega_2} = \frac{1}{n} \int P( \boldsymbol{\mathbf{k}} _1, \boldsymbol{\mathbf{k}} _2)k_1^2k_2^2 \,\mathrm{d}{k_1} \,\mathrm{d}{k_2}
\end{equation}
The triple differential section can be defined as
\begin{equation}
\frac{\partial^3\sigma}{\partial E_1\partial\Omega_1\partial\Omega_2} = \frac{k_1}{n} \int P( \boldsymbol{\mathbf{k}} _1, \boldsymbol{\mathbf{k}} _2)k_2^2 \,\mathrm{d}{k_2}
\end{equation}
In Aihua's paper, fix $ \boldsymbol{\mathbf{k}} _1$ and make $ \boldsymbol{\mathbf{k}} _2$ rotate on a certain plane containing $ \boldsymbol{\mathbf{k}} _1$ to obtain the experimental diagram of the angular distribution.
Single photon absorption
Analogous hydrogen atom absorbs a photon (first-order perturbation) after it has $\Delta l = \pm 1$ (it can be imagined that the photon has an angular momentum quantum number 1), and the total angular momentum will increase after helium atom absorbs a photon from the ground state. So to consider the single electron absorption, we can only observe all the partial waves of $L = 1$, and even use the so-called "weak field TDSE"2
\begin{equation}
\mathrm{i} \frac{\partial \Psi^{^1P}}{\partial t} = H_0 \Psi^{^1P} + H_E \Psi^{^1S} \mathrm{e} ^{- \mathrm{i} E_0 t}
\end{equation}
Among them, $H_0$ is the Hamiltonian without field, and $H_E = E(t) (z_1 + z_2)$ is the Hamiltonian for electric field.
\begin{equation}
\Psi^{^1P}( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2, t) = \sum_{l_1, l_2} \frac{1}{r_1 r_2} \psi_{l_1, l_2}^{1, 0} (r_1, r_2, t) \mathcal Y_{l_1, l_2}^{1, 0} ( \hat{\boldsymbol{\mathbf{r}}} _1, \hat{\boldsymbol{\mathbf{r}}} _2)
\end{equation}
The probability of absorbing a photon is to sum all the partial waves of $L = 1$. However, if you do Streaking, it usually absorbs multiple infrared photons, so many subwaves are still needed.
Angular basis transformation
In addition to using the generalized spherical harmonic function as the angular basis, another natural choice is the product of two spherical harmonic functions. If the total wave function is
\begin{equation}
\Psi( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) = \frac{1}{r_1 r_2} \sum_{l_1, m_1, l_2, m_2} \psi_{l_1, l_2}^{m_1, m_2}(r_1, r_2) Y_{l_1, m_1}( \hat{\boldsymbol{\mathbf{r}}} _1) Y_{l_2, m_2} ( \hat{\boldsymbol{\mathbf{r}}} _2)
\end{equation}
Then from
eq. 3 and
eq. 1
\begin{equation}
\psi_{l_1, l_2}^{m_1, m_2}(r_1, r_2) = \sum_{L,M} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \psi_{l_1, l_2}^{L,M}(r_1, r_2)
\end{equation}
If we have a list of $l_1, l_2, L, M$ and a list of $l_1, l_2, m_1, m_2$, the transformation can be represented by a matrix.
Single electron probability distribution
From ex. 1 , the probability density of finding any particle at $ \boldsymbol{\mathbf{r}} _1$ is
\begin{equation}
P( \boldsymbol{\mathbf{r}} _1) = 2\int \left\lvert \Psi( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) \right\rvert ^2 \,\mathrm{d}^{3}{r_2}
\end{equation}
If you use the wave function
3, the form of this integral becomes quite simple, shilling $ \boldsymbol{\mathbf{r}} _1$ is a constant, calculate the radial wave function about $r_2$
\begin{equation}
\begin{aligned}
\psi_{l_2, m_2}( \boldsymbol{\mathbf{r}} _1, r_2) &= \frac{1}{r_1} \sum_{l_1, m_1} \psi_{l_1, l_2}^{m_1, m_2}(r_1, r_2) Y_{l_1, m_1}( \hat{\boldsymbol{\mathbf{r}}} _1)\\
&= \frac{1}{r_1} \sum_{L,M,l_1} \begin{bmatrix}l_1 & l_2 & L\\ M - m_2 & m_2 & M\end{bmatrix} \psi_{l_1, l_2}^{L,M}(r_1, r_2) Y_{l_1, m_1}( \hat{\boldsymbol{\mathbf{r}}} _1)
\end{aligned} \end{equation}
Then
eq. 25 becomes
\begin{equation}
\Psi( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) = \frac{1}{r_2} \sum_{l_2, m_2} \psi_{l_2, m_2}( \boldsymbol{\mathbf{r}} _1, r_2) Y_{l_2, m_2} ( \hat{\boldsymbol{\mathbf{r}}} _2)
\end{equation}
The single electron probability is
\begin{equation}
\begin{aligned}
P( \boldsymbol{\mathbf{r}} _1) &= 2\sum_{l_2, m_2} \int \left\lvert \psi_{l_2, m_2}( \boldsymbol{\mathbf{r}} _1, r_2) \right\rvert ^2 \,\mathrm{d}{r_2} \\
& = \frac{2}{r_1^2} \sum_{l_2, m_2} \int \left\lvert \sum_{l_1, m_1} \psi_{l_1, l_2}^{m_1, m_2}(r_1, r_2) Y_{l_1,m_1}( \hat{\boldsymbol{\mathbf{r}}} _1) \right\rvert ^2 \,\mathrm{d}{r_2}
\end{aligned}
\end{equation}
Single ionization momentum spectrum
In the $P(r_1, r_2)$ distribution, single ionization can generally be distinguished from double ionization. You only need to take one box: $r_1 < r_{1,max}$, $r_2 > r_{2,min}$
\begin{equation}
\left\langle \delta_{ \boldsymbol{\mathbf{r}} _1} \right\rvert \left\langle C_{l_2,m_2}(k_2) \middle| \Psi \right\rangle = \frac{1}{r_1} \sum_{l_1, m_1} \psi_{l_1,m_1}(r_1) \cdot Y_{l_1, m_1}( \hat{\boldsymbol{\mathbf{r}}} _1)
\end{equation}
\begin{equation}
\psi_{l_1,m_1}(r_1) = \int_{r_{2,min}}^{+\infty} \sqrt{\frac{2}{\pi}}F_{l_2}(\eta_2, k_2 r_2)\psi_{l_1, l_2}^{m_1, m_2}(r_1, r_2) \,\mathrm{d}{r_2}
\end{equation}
However, what streaking trace needs is not the probability of the above formula but the probability of a certain direction. We can use
eq. 13 for linear combination
\begin{equation}
\left\langle \delta_{ \boldsymbol{\mathbf{r}} _1} \right\rvert \langle{\psi_{ \boldsymbol{\mathbf{k}} _2}^{(-)}}|{\Psi}\rangle = \sum_{l_2,m_2} \frac{ \mathrm{i} ^{-l_2}}{k_2} \mathrm{e} ^{ \mathrm{i} \sigma_{l_2}}Y_{l_2,m_2}( \hat{\boldsymbol{\mathbf{k}}} _2) \left\langle \delta_{ \boldsymbol{\mathbf{r}} _1} \right\rvert \left\langle C_{l_2,m_2}(k_2) \middle| \Psi \right\rangle
\end{equation}
So (
eq. 31 is substituted into
eq. 33 , and it is sorted into the form of $1/r_1\sum \dots Y_{l_1,m_1}( \boldsymbol{\mathbf{r}} _1)$)
\begin{equation}
\begin{aligned}
P( \boldsymbol{\mathbf{k}} _2) &= \int_0^{r_{1,max}} \left\lvert \left\langle \delta_{ \boldsymbol{\mathbf{r}} _1} \right\rvert \langle{\psi_{ \boldsymbol{\mathbf{k}} _2}^{(-)}}|{\Psi}\rangle \right\rvert ^2 \,\mathrm{d}^{3}{r_1} \\
&= \sum_{l_1,m_1} \int_0^{r_{1,max}} \left\lvert \sum_{l_2,m_2} \frac{ \mathrm{i} ^{-l_2}}{k_2} \mathrm{e} ^{- \mathrm{i} \sigma_{l_2}} Y_{l_2,m_2}( \hat{\boldsymbol{\mathbf{k}}} _2) \cdot \psi_{l_1,m_1}(r_1) \right\rvert ^2 \,\mathrm{d}{r_1}
\end{aligned}
\end{equation}
Shakeup ionization
Let the normalized ground state be
\begin{equation}
\phi_{n_1, l_1, m_1}( \boldsymbol{\mathbf{r}} _1) = R_{n_1,l_1}(r_1)Y_{l_1,m_1}( \hat{\boldsymbol{\mathbf{r}}} _1)
\end{equation}
The normalized scattering state is
\begin{equation}
C_{l_2, m_2}(k_2, \boldsymbol{\mathbf{r}} _2) = \frac{1}{r_2}\sqrt{\frac{2}{\pi}} F_{l_2}(\eta_2, k_2 r_2) Y_{l_2, m_2} ( \hat{\boldsymbol{\mathbf{r}}} _2)
\end{equation}
The following calculation is very similar to the single-electron probability distribution. We also don't need to care too much about the symmetry of the particle exchange. Just multiply the inner product by $\sqrt{2}$ and the probability by two. Similarly, first convert the total wave function to the form of
eq. 25 , and project to get
\begin{equation}
\left\langle \phi_{n_1, l_1, m_1} \right\rvert \left\langle C_{l_2, m_2}(k_2) \middle| \Psi \right\rangle = \sqrt{\frac{2}{\pi}} \iint r_1 R_{n_1,l_1}(r_1) F_{l_2}(\eta_2, k_2 r_2) \psi_{l_1, l_2}^{m_1, m_2}(r_1, r_2) \,\mathrm{d}{r_1} \,\mathrm{d}{r_2}
\end{equation}
However, what streaking trace needs is not the probability of the above formula but the probability of a certain direction. We only need to use
eq. 13 for linear combination
\begin{equation}
\left\langle \phi_{n_1, l_1, m_1} \right\rvert \langle{\psi_{ \boldsymbol{\mathbf{k}} _2}^{(-)}}|{\Psi}\rangle = \sum_{l_2,m_2} \frac{ \mathrm{i} ^{-l_2}}{k_2} \mathrm{e} ^{ \mathrm{i} \sigma_{l_2}} Y_{l_2,m_2}( \hat{\boldsymbol{\mathbf{k}}} _2) \left\langle \phi_{n_1, l_1, m_1} \right\rvert \left\langle C_{l_2, m_2}(k_2) \middle| \Psi \right\rangle
\end{equation}
\begin{equation}
P(n_1, l_1, m_1, \boldsymbol{\mathbf{k}} _2) = \left\lvert \left\langle \phi_{n_1, l_1, m_1} \right\rvert \langle{\psi_{ \boldsymbol{\mathbf{k}} _2}^{(-)}}|{\Psi}\rangle \right\rvert ^2
\end{equation}
When it is $M \equiv 0$, the sum of $m_2$ is actually only one item, $m_2 = - m_1$, which is also a constant. In theory, if eq. 39 is added to all $n_1, l_1, m_1$, it should be equal to eq. 34 .
Consider the solid angle of the detector
In $P(n_1, l_1, m_1, \boldsymbol{\mathbf{k}} _2)$, the detector actually captures the photoelectrons in a cone with a solid angle of $\Omega$. If we suppose $\Omega \to 0$, the total probability is
\begin{equation}
\int P(n_1, l_1, m_1, \boldsymbol{\mathbf{k}} _2) k_2^2 \Omega \,\mathrm{d}{k_2}
\end{equation}
So the actual measured photoelectron momentum spectrum should be
\begin{equation}
\Omega k_2^2 P(n_1, l_1, m_1, k_2)
\end{equation}
Relative delay analysis of single ionization
Osiander/Pazourek's interpretation of the delay analysis is: at the moment of single ionization, helium ions exist in electric dipole (almost not change with time), which exerts an attractive effect on the ionized electrons. But it does not seem to explain how to draw this dipole from the wave function.
Obviously, this cannot be drawn with the aforementioned single-electron probability distribution, because the ionization rate in the current calculation is only on the order of 1e-6 to 1e-5. This shows that the total wave function is mainly still in the bound state. I think the correct way is to calculate the location distribution $P( \boldsymbol{\mathbf{r}} _1) = P( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2)$, where $ \boldsymbol{\mathbf{r}} _2$ can take the average position of the ionization wave packet or integrate all the positions of the ionization wave packet
\begin{equation}
P( \boldsymbol{\mathbf{r}} _1) = \int P( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) \,\mathrm{d}^{3}{r_2}
\end{equation}
Alternatively, you can also calculate $P( \boldsymbol{\mathbf{r}} _1) = P( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{k}} _2)$, where $ \boldsymbol{\mathbf{k}} _2$ takes the central momentum of the ionization wave packet or the integral of all the momentum of the wave packet (make sure that only the ionization wave packet is included)
\begin{equation}
P( \boldsymbol{\mathbf{r}} _1) = \int P( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{k}} _2) \,\mathrm{d}^{3}{k_2}
\end{equation}
Or integrate the shakeup or direct interval of a vertical line on the streaking trace.
\begin{equation}
P( \boldsymbol{\mathbf{r}} _1) = \int P( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{k}} _2) \,\mathrm{d}{k_2}
\end{equation}
I think the last approach is the most reliable, after all, streaking trace is a momentum spectrum.
How to calculate $P( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{k}} _2)$ is discussed below. The wave function can be written as
\begin{equation}
\Psi( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) = \frac{1}{r_1}\sum_{l_1,m_1} \left[\frac{1}{r_2}\sum_{l_2,m_2} \psi_{l_1,l_2}^{m_1,m_2}(r_1, r_2) Y_{l_2,m_2}( \hat{\boldsymbol{\mathbf{r}}} _2) \right] Y_{l_1,m_1}( \hat{\boldsymbol{\mathbf{r}}} _1)
\end{equation}
The content in the brackets is projected to the Coulomb plane wave $ \left\lvert \boldsymbol{\mathbf{k}} _2 \right\rangle $ item by item to get $\psi_{l_1,m_1}(r_1)$, just follow the method of eq. 14 . then
\begin{equation}
P( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{k}} _2) = 2 \left\lvert \frac{1}{r_1}\sum_{l_1,m_1} \psi_{l_1,m_1}(r_1) Y_{l_1,m_1}( \hat{\boldsymbol{\mathbf{r}}} _1) \right\rvert ^2
\end{equation}
Multiplying by 2 considers identical particles.
Calculation result: The dipole is not seen, probably because the electron cloud near the nucleus contributes a large part when calculating the projection $\psi_{l_1,m_1}(r_1)$. It may be necessary to remove the nucleus and then project it.
In fact, it is better to directly calculate the distribution of $P( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2)$ if it is removed, and keep $ \boldsymbol{\mathbf{r}} _2$ unchanged (the exact center of the photoelectron wave packet can be taken). Just use the 2 times modulo square of eq. 45 directly (fixed at $ \boldsymbol{\mathbf{r}} _2$).
Fig. 3：Classical mechanics calculation delay
Use classical electromagnetics and mechanics to try to calculate the delay: $P( \boldsymbol{\mathbf{r}} _2)$ position distribution is two circular dipoles, the outer ring represents direct ionization, and the inner ring represents shake-up. Then we can take $ \boldsymbol{\mathbf{r}} 2$ as the positions of the center of the outer ring and the center of the inner ring in the directions of $z$, calculate the $P( \boldsymbol{\mathbf{r}} _1)$ distributions of the two positions at each time $t$, and calculate the gravitational force of the charge distribution at $ \boldsymbol{\mathbf{r}} $. Then compare the time delay of the gravitation to the two classical mechanics trajectories. Force calculation:
\begin{equation}
\begin{aligned}
&F_z( \boldsymbol{\mathbf{r}} _2) = \iiint \frac{P( \boldsymbol{\mathbf{r}} _1) \cos\alpha}{r_{12}^2} r_1^2\sin\theta_1 \,\mathrm{d}{r_1} \,\mathrm{d}{\theta_1} \,\mathrm{d}{\phi_1} \\
\qquad
&r_{12}^2 = r_1^2 + r_2^2 - 2r_1r_2 \cos\theta_1\\
&\cos\alpha = \frac{r_2^2 + r_{12}^2 - r_1^2}{2r_2 r_{12}}
\end{aligned}
\end{equation}
Fig. 4：estimated $\Delta k$
\begin{equation}
\Delta k_{max} = k_0 \sin\left(2\pi \tau/T_{IR}\right) \approx 2 \pi k_0 \tau/T_{IR}
\end{equation}
In test800, $k_0 \approx 0.042$, $\tau = 16$, $T_{IR} = 110.245$. Get $\Delta k_{max} \approx 0.038$.
1. ^ Theoretically, it doesn’t matter if the bound state is not removed, because the Coulomb wave function of the single electron and the bound state The states are inherently orthogonal. This may be due to numerical errors?
2. ^ refer to Colgan 2001, J Phys. B: 34 (2001) L457-L466 during propagation.
3. ^ in the form of eq. 25 , you don’t need to do this first in the program. You only need to substitute eq. 26 into the following equations and put the sum of $L, M, l_1, l_2$ to the outside.
