# Strong-Field Double Ionization in a Three-Electron System:
Momentum Distribution Analysis

Dmitry K. Efimov<EMAIL_ADDRESS>Institute of Theoretical Physics,
Faculty of Fundamental Problems of Technology, Wrocław University of Science
and Technology, 50-370 Wrocław, Poland  Artur Maksymov  Institute of
Theoretical Physics, Faculty of Physics, Astronomy and Applied Computer
Science, Jagiellonian University, Stanisława Łojasiewicza 11, 30-348 Kraków,
Poland  Jakub Zakrzewski  Instytut Fizyki Teoretycznej, Wydział Fizyki,
Astronomii i Informatyki Stosowanej, Uniwersytet Jagielloński, ulica
Stanisława Łojasiewicza 11, 30-348 Kraków, Poland  Mark Kac Center for Complex
Systems Research, Jagiellonian University, Stanisława Łojasiewicza 11, 30-348
Kraków, Poland  Jakub S. Prauzner-Bechcicki  jakub.prauzner-
<EMAIL_ADDRESS>Jagiellonian University in Kraków, Faculty of Physics,
Astronomy and Applied Computer Science, Marian Smoluchowski Institute of
Physics, Łojasiewicza 11, 30-348 Krakow, Poland

###### Abstract

We study strong-field double ionization in a three-electron system by applying
a simplified, reduced-dimensionality model with three active electrons. The
influence of the spin-induced symmetry of the spatial part of the wavefunction
on the final two-photoectron momentum distribution is discussed. We identify
partial momentum distributions originating from different sets of spins of
outgoing electrons providing in this way a quantum support connection between
V-structure and direct ionization typically explained classically. Changes in
the momentum distribution with increasing field amplitude obtained in our
simplified model are shown to be well-correlated with experimental data known
from the literature. The possible relation between the observed dependencies
and different ionization mechanisms is discussed.

††preprint: APS/123-QED

##  I Introduction

Multiple ionization of atoms and molecules is an important branch of strong-
field and attosecond physics primarily due to the fact that it serves as
scientific playground to study electron-electron correlations. The best
example of such a correlation is the non-sequential multiple ionization [1, 2]
with its characteristic “knee” in the yield versus laser intensity plot [3, 4,
5] and the V-structure in the photoelectron momentum distribution (PMD) for
double ionization [6, 7]. When two electrons are involved in a process, the
well-known three-step model [8] provides the widely accepted mechanism: (i) a
tunneling ionization of a single electron, (ii) its rescattering with a parent
ion, and (iii) a double ionization. When more than two electrons are involved,
the first two steps are similar but, as a result, in the third step one
observes double, triple, and higher ionization events, all of which exhibit
correlations.

The measurement of ionized particles’ momenta is a challenging task,
especially when triple and higher ionization is considered [9, 10, 11, 12, 13,
14, 15]. Yet, both numerical simulations and analytical (including semi-
analytical) descriptions of the process remain a difficult task too,
particularly when a quantum-mechanical explanation is requested. Even for two-
electron systems, the results of full-dimensional quantum simulations are
scarce [16, 17, 18, 19, 20]. The methods of choice in the multiple ionization
case are then reduced-dimensionality quantum models [21, 22, 23, 24, 25, 26,
27] or classical/semi-classical calculations [28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41]. In the following, we shall use the reduced-
dimensionality model for three active electrons introduced earlier [42].

When one constructs theoretical models for double ionization it is usually
assumed that two electrons are active only, even for targets with higher
number of electrons possibly involved (eg. Ne, Ar, Kr). Consequently, it is
often taken for granted that these two active electrons have opposite spins.
The metastable ${}^{3}S$ state of He was analysed in that contexts to expose a
clear suppression of correlated escape of electrons as a main result [43, 44,
45]. When a larger number of electrons is considered, the spin degrees of
freedom become indispensable even though the strong-field Hamiltonian does not
couple to the spin part; the spin part of the wavefunction determines symmetry
of the spatial part of the function. For the sake of simplicity we restrict
the discussion to three-electron systems. Let us consider atoms that have
$s^{2}p^{1}$ electronic configuration - the spatial part of the wavefunction
will be partially symmetric with respect to exchange of electrons, because the
spin part is partially antisymmetric too. On the contrary, elements with
$p^{3}$ electronic configurations have symmetric spin part of the ground state
wavefunction and, thus, the spatial part is completely antisymmetric. We have
addressed the role that is played by symmetry of the spatial part of the
wavefunction in the strong-field ionization of three-electron atoms in our
previous studies focusing on the triple ionization yields and three-electron
momenta distributions [42, 46, 47] and the double ionization yields only [48,
49]. The present study fills up a missing part, i.e. we present momenta
distribution for double ionization events from the three-electron atom.

The description of double ionization of a three-electron atom is a seemingly
simpler task than dealing with triple ionization, especially when it comes to
presenting electron momenta distributions [47, 40]. However, the task is not
that easy when the spin degrees of freedom are properly included. We have
shown that in the general case, double ionization of a three-electron atom
cannot be reproduced by a straightforward combination of judiciously chosen
two-active electron models [48]. Thus, the two-electron (2E) PMDs of the
three-electron model cannot be reproduced by combinations of PMDs from the
two-electron models. Let us stress that the latter conclusion has been made
solely by examination of ionization potentials diagram together with the
results of two-electron ionization yields simulations.

In the following, we simulate interaction of a three-electron atom with a
strong-laser field while focusing on analysis of double ionization events in
terms of 2E–PMDs. Since our simulation can reveal the momentum distribution
for different pairs of electrons, we ask ourselves whether they can be used
for uncovering some properties of a strong-field dynamics. As we show later,
there is a kind of structure in PMD, namely the so-called V-structure, that
can be put into direct correspondence with the particular channel for
ionization event. Another finding that concerns the asymmetry property of the
PMDs is discussed extensively as well.

The work is organized as follows. First, we describe the physical and
numerical models used in the simulations, then results are presented and
discussed. The article ends with conclusions and acknowledgments. Atomic units
are used throughout the paper unless otherwise noted.

##  II Description of the model

Numerical treatment of tree-electron atoms in a full configuration space is a
formidable task. Therefore, we resort to a restricted-space model capable of
capturing the essence of strong field electron dynamics [29, 42]. In the
model, each of the electrons move along one-dimensional (1D) track inclined at
a constant angle $\alpha$ ($\tan\alpha=\sqrt{2/3}$) with respect to the
polarization axis of the laser pulse and forming a constant angle $\pi/6$
between each pair of tracks. The chosen geometry results from the local
stability analysis applied to study the adiabatic full-dimensional potential
[29, 50, 51] and has been applied previously to study three-electron dynamics
yielding plausible results [42, 48, 49, 46, 47].

The Hamiltonian reads:

$H=\sum_{i=1}^{3}\frac{p^{2}_{i}}{2}+V+V_{int},$ (1)

where $V$ and $V_{int}$ are the atomic and interaction potentials,
respectively. The atomic potential in the restricted-space has a form:

$V=-\sum_{i=1}^{3}\frac{3}{\sqrt{r^{2}_{i}+\epsilon^{2}}}+\sum_{i,j=1;i<j}^{3}\frac{q^{2}_{ee}}{\sqrt{(r_{i}-r_{j})^{2}+r_{i}r_{j}+\epsilon^{2}}}$
(2)

where a smoothing factor $\epsilon=\sqrt{0.83}$ and an effective charge
$q_{ee}=1$ allow for reproduction of the triple ionization potential of Neon
atom ($I_{p}=4.63$ a.u.). The interaction term is described as follows:

$V_{int}=\sqrt{\frac{2}{3}}A(t)\sum_{i=1}^{3}p_{i},$ (3)

where the vector potential is given as

$A(t)=\frac{F_{0}}{\omega}\sin^{2}\left(\frac{\pi
t}{T_{p}}\right)\sin\left(\omega t+\phi\right)$

for $t\in[0,T_{p}]$. The laser pulse parameters are: field amplitude, $F_{0}$;
the carrier frequency, $\omega=0.06$ a.u.; the pulse length $T_{p}=2\pi
n_{c}/\omega$; the number of optical cycles, $n_{c}=3$; and the carrier-
envelope phase, $\phi$, which is set to zero.

As we are interested in studying the interaction of three-electron atoms with
a laser pulse we need a wavefunction that describes the motion of three
electrons (the nucleus is assumed to be infinitely heavy). The wavefunction is
composed of spatial and spin parts and, as a whole, has to be antisymmetric
with respect to exchange of electrons. Depending on the spin configuration,
the spatial part will have different symmetry properties with respect to the
exchange of electrons (see detailed discussion in Appendix A). For the system
described with simplified Hamiltonian, eq. (1), and the atom with $s^{2}p^{1}$
electron configuration the spin part is such that two electrons have the same
spin, i.e. the spin part needs to be symmetric with respect to exchange of
electrons in that pair. Consequently, the spatial part needs to be
antisymmetric with respect to exchange of electrons in that pair.

Without loss of generality we may write the wavefunction for such atom in the
restricted-space model as (see eq. (13)):

$\Psi\propto\Psi_{12}(r_{1},r_{2},r_{3})|{\rm UUD}\rangle\\\
+\Psi_{23}(r_{1},r_{2},r_{3})|{\rm DUU}\rangle\\\
+\Psi_{13}(r_{1},r_{2},r_{3})|{\rm UDU}\rangle.$ (4)

Here $U$ and $D$ stand for electron with spin-up and spin-down, respectively.
The subscripts denote the pair of electrons with respect to exchange of which
the wavefunction is antisymmetric.

The Hamiltonian, eq. (1), does not influence the spin part during evolution,
therefore for the sake of simplicity, in numerical implementation it is enough
to evolve one of the three terms on the right side of the eq. (4) only.
However, as explained in detail in the Appendix, special caution must be taken
with the symmetry properties of the wave function. In the following, we shall
evolve the wavefunction $\Psi_{12}(r_{1},r_{2},r_{3})$, i.e. a function that
is antisymmetric with respect to exchange of electrons $1$ and $2$, but with
respect to exchange of electrons $1\leftrightarrow 3$ and $2\leftrightarrow 3$
it is neither symmetric nor antisymmetric. The initial wavefunction, mimicking
the ground state, is found using imaginary time propagation of the
Hamiltonian, eq. (1), without the interaction with the field and enforcing the
proper symmetry [42].

The numerical solution of the time-dependent Schrödinger equation (TDSE) with
the Hamiltonian, eq. (1), is essentially an extension of the method used in
the two-dimensional case described in extenso elsewhere [24, 47]. Here, we
recall this method only very briefly. The TDSE is solved on a large
3-dimensional grid with the use of a split-operator technique and the fast
Fourier transform algorithm. The goal is to obtain momenta distributions of
outgoing electrons, therefore the wavefunction cannot be absorbed at edges of
the grid as it is typically implemented. In our approach, we divide the whole
evolution space into regions, i.e., the “bounded motion” and “outer” regions.
In the “bounded motion” region the evolution is done without further
simplifications, whereas in the “outer” regions the Hamiltonian is simplified,
i.e. the interaction of electrons with nucleus and other electrons is
successively neglected resulting in regions describing a single ion, double
ion and triple ion. Thanks to these simplifications, in the “outer” regions
the wavefunction can be represented and evolved in the momentum representation
- the evolution simplifies to the multiplication by an appropriate phase
factor. The transfer between the “bounded motion” and “outer” regions is done
by a smooth cutting and coherent adding of the wavefunction following the
procedure introduced in [21].

At the end of the simulation, the wavefunction from “outer” regions can be
integrated over the “bounded” part, thus leaving the momentum distributions
corresponding to single ionization, double ionization and triple ionization.
More details of the algorithm for the simulation of momentum distributions of
the three-electron atom can be found in the supplementary material of the work
[47].

Strong-field double ionization is a quite convoluted process. For a purpose of
the present work we shall divide the processes into two categories, i.e.
recollision-impact ionization (RII) and time-delayed ionization (TDI). The RII
is a process in which the rescattering electron ionizes the bound electron
upon collision and is usually connected with the V-structure (also named a
finger-like structure) in PMD [6]. The TDI includes sequential ionization
(independent electron releases) [52], recollision-excitation with subsequent
ionization (RESI) of all kinds [53, 54] and other possible mechanisms (eg. the
slingshot non-sequential ionization [38]).

The ionization yields ascribed to different ionization channels are calculated
in our model by an integration of probability fluxes over the border of the
appropriate regions. The RII yield is calculated by integration of the direct
probability flux over the border between the regions describing a neutral atom
and that for the double ion, whereas TDI yield is calculated by an integration
of the flux over the border between the regions describing a single ion and
double ion. Detailed procedure is described elsewhere [24, 27, 46].

The grid has $n=1024$ nodes in each direction, the step in the grid is ${\rm
d}r=0.195$. The time-step is ${\rm d}t=0.05$ and total number of steps amounts
to $6500$.

##  III Results and discussion

In the following we shall analyze the final wavefunction of the system. Once
the laser pulse is ended after long enough time one may plausibly assume that
the wavefunction can be written as:

$\Psi=\Psi_{N}+\Psi_{SI}+\Psi_{DI}+\Psi_{TI},$ (5)

where $\Psi_{N}$, $\Psi_{SI}$, $\Psi_{DI}$ and $\Psi_{TI}$ are wavefunctions
of an neutral atom - three bound electrons, a single ion - a single freed and
two bounded electrons, a double ion - two freed electrons and a single bounded
one, and a triple ion - three freed electrons, respectively. We refer to freed
electrons as we are focused on electron dynamics. Thanks to the division into
the “bounded motion” and “outer” regions in our numerical approach we have
indeed access to these four wavefunctions. Before we continue our analysis,
let us note that it is justified to expect that the wavefunctions, eq. (5),
which the final wavefunction is decomposed to, are at least approximately
orthogonal to each other after enough time passed since the action of the
laser field, i.e. one has a set of properties:
$\langle\Psi_{N}|\Psi_{SI}\rangle\simeq
0,\langle\Psi_{N}|\Psi_{DI}\rangle\simeq 0,$ etc. Therefore, the final
observable, a PMD may be written:

$|\Psi|^{2}\simeq|\Psi_{N}|^{2}+|\Psi_{SI}|^{2}+|\Psi_{DI}|^{2}+|\Psi_{TI}|^{2}.$
(6)

The last term, $|\Psi_{TI}|^{2}$, is a three-electron (3E) PMD and have been
already studied elsewhere [47]. Here we concentrate on a 2E–PMD, i.e.
$|\Psi_{DI}|^{2}$. Let us stress that our model allows us to study a three-
electron atom and extract 2E–PMDs, so the influence of the third electron is
treated explicitly in the evolution and may be traced in the distributions.
Typically, when double ionization of multi-electron atoms is considered one
refers to some kind of two-active electron model, where the influence of the
third electron is neglected.

The 2E–PMD is an experimental observable. In coincidence measurements,
momentum of a double ion is measured together with momentum of one of the
electrons (see for example [55]). The momentum of the remaining electron is
retrieved from the momentum conservation principle. Our approach allows to
shed some more light on components contributing to the observed PMDs, and
originating in different processes.

The final wavefunction, eq. (5), is a three-electron wavefunction. Its part
that describes double ions, that is, $\Psi_{DI}$, has, in principle, a form of
eq. (4). The first term, i.e. $\Psi_{12}(r_{1},r_{2},r_{3})$, is antisymmetric
with respect to the exchange of electrons 1 and 2, while it does not possess
symmetries with respect to the exchange of other pairs of electrons. The other
terms posses similar properties as applied to different pairs of electrons. To
obtain 2E–PMD we need to write the wavefunction in momentum representation,
therefore the Fourier transform is applied. The Fourier transform does not
change the symmetry properties of the wave function. Next step is to take
modulus squared of the wavefunction, and here the fact that spin parts are
orthogonal allows us to write it in a formal form:

$|\Psi_{DI}(p_{1},p_{2},p_{3}))|^{2}\propto\\\
|\Phi^{DI}_{12}(p_{1},p_{2},p_{3})|^{2}+|\phi^{DI}_{23}(p_{1},p_{2},p_{3})|^{2}\\\
+|\phi^{DI}_{13}(p_{1},p_{2},p_{3})|^{2}.$ (7)

This operation changes symmetry, so for example the first term of eq. (7)
becomes symmetric with respect to exchange of electrons 1 and 2. The last step
in pursue to 2E–PMD is to integrate eq. (7) with respect to the variable
describing the bound electron. Let us assume the electron 3 is still bound,
the integration over $p_{3}$ variable leaves us with:

$PMD(p_{1},p_{2})\propto\\\
|\Phi^{DI}_{12}(p_{1},p_{2})|^{2}+|\phi^{DI}_{23}(p_{1},p_{2})|^{2}\\\
+|\phi^{DI}_{13}(p_{1},p_{2})|^{2}.$ (8)

The first term is symmetric with respect to exchange of electrons 1 and 2 but
the other two terms are not necessarily. Overall, 2E–PMD has to be symmetric
with respect to the exchange of electrons because the electrons are
indistinguishable. Therefore, the sum of the second and third terms in eq. (8)
has to be symmetric. It can be expected that when the first term of eq. (8) is
plotted alone it should have zeros on the $p_{1}=p_{2}$ axis and be symmetric
with respect to that axis, while the other two terms if plotted separately
could be asymmetric with respect to $p_{1}=p_{2}$, but the axis itself does
not need to be filled with zeros.

Figure 1:  Two-photoelectron momentum distributions for field amplitude F=0.12
a.u. (a) partial distribution - symmetric part, corresponding to the first
term on the right side of eq. (8), (b) partial distribution - asymmetric part,
corresponding to the third term on the right side of eq. (8), (c) the sum of
the second and third terms on the right side of eq. (8), (d) the full 2E-PMD,
eq. (8).

These symmetry expectations are confirmed in Fig. 1 where we present an
exemplary 2E–PMD. In panel (a) the partial 2E–PMD is clearly symmetric with
respect to the $p_{1}=p_{2}$ and the axis has pronounced zeros. The zeros on
the axis are remnants of antisymmetry of that part of the wavefunction and
have their origin in the Pauli exclusion principle (a similar effect was
discussed in Helium in [45]). In Fig. 1(b), where the third term of the full
2E–PMD, eq. (8), is depicted alone, the distribution is clearly asymmetric
with respect to $p_{1}=p_{2}$ axis and has non-zero values on the axis;
however, at the same time it contains elements resembling the famous V
structure. The sum of the second and third terms is shown in Fig. 1(c). It
exhibits the V-structure undoubtedly - it is more prominent in the first
quadrant because the simulations are done for a very short laser pulse (only 3
cycles long) and for a single carrier-envelope phase. The V structure is still
visible in the entire 2E–PMD; see Fig. 1(d). From the inspection of the full
distribution, one sees that the V-structure is dominant feature, and the
signal originating in the $|\Phi^{DI}_{12}(p_{1},p_{2})|^{2}$ part is barely
visible. The partial PMDs cannot be measured experimentally, as the electrons
are indistinguishable in the detector. The above discussion is valid
exclusively within the framework of the theoretical model. However, the
observed differences in the partial PMDs are the consequence of the initial-
state symmetry, i.e. for atoms with totally antisymmetric spatial part (the
spin part is symmetric as for $p^{3}$ valence shell systems, [49]) all the
partial PMDs look the same.

Furthermore, the 2E–PMD shown in Fig. 1(a) can be interpreted according to the
analysis presented by Staudte et al. in [6]. The suppression of correlated
escape of electrons connected with the Pauli exclusion principle is reserved
to binary collisions, whereas the observed features can be connected with
recoil collisions exclusively.

The lack of symmetry in relation to the electron exchange visible in Fig.
1(b), which we will call an asymmetry, follows directly from the symmetry
properties of the initial state and the decomposition of the final
wavefunction, see eq. (8). Although the existence of such an asymmetry is
clear, the degree of it does not follow directly from the definition of the
model and should rather depend on the particular properties of ionization
dynamics. Thus, it is an interesting task to reveal the connection between
ionization and partial momentum distribution asymmetry.

The partial PMD shown in Fig. 1(b) possesses a slight asymmetry, but gives no
clue whether the asymmetry depends on the physical parameters such as the
field amplitude. One could suspect that interaction with pulses having larger
amplitudes may introduce more asymmetric partial distributions, mainly because
the time-delayed double ionization, such as sequential ionization or RESI, is
more pronounced for these conditions. In fact, such an increase in asymmetry
was observed in our calculations.

Let us introduce the asymmetry parameter of the partial 2E–PMD obtained by
providing a comparison between the above-diagonal and below-diagonal regions
of the distribution. In particular, we collect information about the
difference of points located symmetrically in relation to the main diagonal.
The asymmetry parameter $\beta$ is defined as follows:

$\beta=\cfrac{\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\left||\Psi(p_{1},p_{2})|^{2}-|\Psi(p_{2},p_{1})|^{2}\right|\,dp_{1}\,dp_{2}}{\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}|\Psi(p_{1},p_{2})|^{2}\,dp_{1}\,dp_{2}}.$
(9)

The $\beta$ parameter calculated for different field intensities for the part
corresponding to the first term in eq. (8) was unsurprisingly equal to 0. The
dependence of $\beta$ on the field value for the third term in eq. (8) is
shown in Fig. 2. Within the standard laser field amplitude interval, in which
one observes the famous “knee” structure in double ionization yield, the
$\beta$ parameter increases considerably, therefore, let us contrast it with
ionization yields for different channels as a function of field amplitude, see
Fig. 3.

Figure 2: Asymmetry parameter computed for the partial 2E–PMD (the third term
of eq. (8)) for different field values. Figure 3: Yields for different
channels of double ionization. The TDI stands for Time-dependent ionization,
that includes sequential ionization and recollision-excitation with subsequent
ionization channel (RESI); the RII, or Recollisionally Induced Ionization
stands for a set of channels that assume direct ionization after recollision.
The order in which electrons ionize with respect to the spin degree of freedom
is denoted by capital letters from left to right, ex. (0-U-D) stands for TDI
in which the electron with spin up (U) is ionized first; if a direct
ionization happened, the letters denoting spins are written close to each
other, ex. (0-DU).

As discussed in the previous section, the Hamiltonian, (eq. (1)), does not
influence the spin part during the evolution, and therefore in our numerical
implementation we evolve only the first term of the three terms on the right
side of eq. (4). Consequently, the spin degree of freedom is “encoded” in
electron’s index and it is possible to trace the yield of spin-resolved
ionization channels [42]. In Fig. 3 we show the dependence of ionization yield
on the field amplitude split into TDI and RII channels, each of which is
further assigned with respect to the order of the escaping electrons.
Eventually, we calculate five partial yields, i.e. RII (0-DU) – the direct
double ionization of electrons with opposite spins, RII (0-UU) – the direct
double ionization of electrons with the same spin, TDI (0-U-D) – the time-
delayed double ionization in which the first freed electron has spin up and
the second has spin down, TDI (0-D-U) – the time-delayed process with reversed
sequence of escaping electrons, and finally TDI (0-U-U) – the time-delayed
ionization of electrons with the same spin. Within the discussed model, these
partial yields can be easily connected with the partial 2E–PMDs depicted in
Fig. 1 (a)–(c). Indeed, the fluxes that allow calculation of yields RII (0-UU)
and TDI (0-U-U) are also responsible for filling the region of the integration
grid from which the partial 2E–PMD with symmetric distribution is retrieved
(see Fig. 1(a)). Similarly, the fluxes allowing calculation of RII (0-DU), TDI
(0-D-U) and TDI (0-U-D) yields are connected with partial PMDs depicted in
Fig. 1(b) and (c).

From comparison between Fig. 2 and Fig. 3 it is hard to discriminate which of
the identified ionization channels is responsible for changes in $\beta$. All
channels that contribute to partial 2E–PMD for which the asymmetry parameter
is calculated increase monotonically with increasing field amplitude. The only
important fact in the ionization yield plot is the change of order, i.e. for
low amplitudes the RII (0-DU) channel dominate over the two TDI channels,
whereas for higher amplitudes, when sequential double ionization start to
dominate the TDI (0-U-D) channel prevails.

It is a known fact that for amplitudes below the knee structure in the
ionization yield plot the non-sequential ionization is the main channel for
double ionization with its acknowledged trace in a form of the V-structure
that, with increasing field intensity, changes into cross-shaped feature [56].
This fact is both reflected in Fig. 2 and in the 2E–PMDs in Fig. 4. The
increase of $\beta$ with increasing field amplitude may be interpreted as a
signature of unequal energy sharing upon recollision.

Figure 4:  Two-photoelectron momentum distributions for field amplitudes (from
left to right) F=0.08, 0.10.0.12, 0.15, 0.18, 0.22, 0.27, 0.32, 0.40, 0.50
a.u. The upper, middle and lower raws correspond to (i) partial distribution -
symmetric part, corresponding to the first term on the right side of eq. (8),
(ii) partial distribution - asymmetric part, corresponding to the third term
on the right side of eq. (8), (iii) the sum of the second and third terms on
the right side of eq. (8). The range of each numerical box is (-3,3) a.u.

The connection between the probability fluxes and partial PMDs allows us to
conclude that the V-structure observed in double ionization of three-electron
atoms corresponds to the ionization via RII (0-UD) channel that has a
distinctively higher yield than RII (0-UU) yield. It is of particular interest
to recall that the V-structure in two-electron momentum distribution is known
for about two decades [6] and, indeed, is clearly associated to the direct
impact channels of double ionization. However, such correspondence has been
set up after simulations employing a purely classical picture of ionization-
recollision processes. The observation we made nicely compliments the
classical argument by adding a quantum evidence of the V-structure connection
to the RII channels even in ionization of three-electron atoms.

Let us now point out to one more ingredient that has its role in the formation
of final PMDs, namely, RESI process. Unfortunately, within the model that we
use we cannot separate RESI channel from the other possible ionisation paths.
Note, we use integration of probability fluxes over borders between regions
defined as atom, single, double and triple ion to calculate ionization yields
and therefore RII yield can be separated from all TDI events, that is
sequential ionisation and all flavors of RESI. Typically, when discussing RESI
one assumes that the ionisation involves excitation of a bound electron to a
singly-excited state (SES) of an ion. However, Liao et al. [54] have shown in
their study that there is a group of channels that involves population of a
doubly excited state (DES) of an atom. Indeed, when the energy of the
returning electron is not sufficient to support direct ionisation that
electron can be trapped by the parent ion and share part of its energy with
the bound electron to form DES. The lifetime of such states relative to
autoionization is typically of order of several femtoseconds thus they may
play an important role in strong-field dynamics. What is more, the relative
participation of DES in ionization is noticeable for low field amplitudes and
decreases monotonically with its increase, whereas the role of SES in RESI
increases with increasing field amplitude [54]. Although the symmetry of the
excited state plays the key role in the final shape of PMD [57, 58] it is
impossible, at the present stage of the analysis performed, undoubtedly to
indicate its role in the change of the $\beta$ parameter. Further research in
that direction is anticipated.

##  IV Conclusions

We have studied strong-field double ionization in a three-electron system by
applying a simplified, reduced-dimensionality model with three active
electrons. Special attention has been paid to the spin-induced symmetry of the
spatial part of the wavefunction and its influence on the final two-
photoelectron momentum distribution. Partial momentum distributions
originating from different sets of spins of outgoing electrons were
identified. Thereafter, the distribution exhibiting the famous V-structure was
undoubtedly linked with the escape of electrons with opposite spins; we have
shown that it corresponds to the RII channel thus providing a quantum support
for the classically grounded connection between the V-structure and direct
ionization. Furthermore, the latter distribution was contrasted with the one
resulting from the escape of electrons with same spins and displaying
suppressed correlation. Changes in the photoelectron momentum distribution
with increasing field amplitude was discussed and found coherent with data
known from literature. Finally, possible relationship between the observed
changes and various RESI mechanisms was discussed and the direction of future
research was indicated.

###### Acknowledgements.

Support by PL-Grid Infrastructure is acknowledged. This work was realized
under National Science Centre (Poland) project Symfonia No.
2016/20/W/ST4/00314 (DKE, AM, JSPB) as well as the OPUS call within the WEAVE
programme 2021/43/I/ST3/01142 (JZ). The research has been also supported by a
grant from the Priority Research Area (DigiWorld) under the Strategic
Programme Excellence Initiative at Jagiellonian University. Views and opinions
expressed in this work are, however, those of the authors only and do not
necessarily reflect those of the European Union, nor any other granting
authority. Neither the European Union nor any granting authority can be held
responsible for them. No part of this work was written with the help of
artificial intelligence software.

##  Appendix A Symmetry properties of wavefunction for three electrons

Here we shall discuss the symmetry properties of the spatial part of the
wavefunction for three electrons. Let us denote all the quantum numbers needed
to describe the state as $Q_{k}=(q_{k},m_{s,k})$, where $m_{s}$ is a spin
projection and $q$ other necessary quantum numbers; by $R_{k}=(r_{k},k)$ let
us denote all the variables needed to describe the electron, that is, spatial
coordinates and the spin variable $k$. Because electrons are fermions, the
wavefunction has to be antisymmetric with respect to exchange of any pair of
electrons, and therefore it can be written as a Slater determinant. For three
electrons, it reads:

$\Psi\propto\left|\begin{array}[]{ccc}\Psi_{Q_{1}}(R_{1})&\Psi_{Q_{1}}(R_{2})&\Psi_{Q_{1}}(R_{3})\\\
\Psi_{Q_{2}}(R_{1})&\Psi_{Q_{2}}(R_{2})&\Psi_{Q_{2}}(R_{3})\\\
\Psi_{Q_{3}}(R_{1})&\Psi_{Q_{3}}(R_{2})&\Psi_{Q_{3}}(R_{3})\end{array}\right|$
(10)

Each term of the determinant is composed of three single-electron
wavefunctions, each possessing a spatial and spin part, for example:

$\Psi_{Q_{1}}(R_{1})\Psi_{Q_{2}}(R_{2})\Psi_{Q_{3}}(R_{3})=\\\
\psi_{q_{1}}(r_{1})\phi_{m_{s,1}}(1)\psi_{q_{2}}(r_{2})\phi_{m_{s,2}}(2)\psi_{q_{3}}(r_{3})\phi_{m_{s,3}}(3)=\\\
\Psi(r_{1},r_{2},r_{3})|m_{s,1}m_{s,2}m_{s,3}\rangle$ (11)

In the last line, we collected all the spatial parts and spin parts into two
separate functions keeping the order. Other quantum numbers than $m_{s}$ are
omitted in subscript for clarity. When calculating the evolution with our
numerical models we need to take only the spatial part with proper symmetry,
i.e. when the spin part is symmetric with respect to the exchange of electrons
(for example, when it is $|{\rm DDD}\rangle$), the spatial part has to be
totally antisymmetric.

In the paper, we discuss the configuration where $m_{s,1}={\rm U}$,
$m_{s,2}={\rm U}$ and $m_{s,3}={\rm D}$. So the whole wavefunction has a
following form:

$\Psi\propto\left(\Psi(r_{1},r_{2},r_{3})-\Psi(r_{2},r_{1},r_{3})\right)|{\rm
UUD}\rangle\\\
+\left(\Psi(r_{2},r_{3},r_{1})-\Psi(r_{3},r_{2},r_{1})\right)|{\rm
DUU}\rangle\\\
+\left(\Psi(r_{3},r_{1},r_{2})-\Psi(r_{1},r_{3},r_{2})\right)|{\rm
UDU}\rangle.$ (12)

The structure of the wavefunction is simple, it is a sum of three products of
spatial and spin functions.

As the Hamiltonian describing the evolution of a system does not have a spin
operator part, the spin function remains unchanged, and it is enough to use
the spatial function in studying the dynamics. Furthermore, it is not
necessary to evolve three separate functions, since proper results are
obtained by simple permutation of variables. However, caution should be
exercised when discussing the symmetry properties of the partial
wavefunctions.

The entire wavefunction is antisymmetric with respect to the exchange of each
pair of electrons, but each term on the right of eq. (12) separately behaves
differently. The first term in eq. (12) has spatial part that is antisymmetric
with respect to exchange of electrons labeled with 1 and 2 subscripts, and is
accompanied with a symmetric spin part with respect to the same exchange.
However, if one exchanges electrons 2 and 3 instead, the spatial part is
neither symmetric nor antisymmetric, similarly the spin part has no symmetry
with respect to that exchange. Likewise the second term in eq. (12) has the
property with respect to electrons 2 and 3, and the third term with respect to
electrons 1 and 3, respectively.

Now, let us consider exchange of electrons 1 and 3: in the first term, the
spin part changes into the one present in the second term and the spatial part
changes accordingly into the one present in the second term but with a minus
sign. Here again, if the first term is analyzed only, the spatial part is
neither symmetric nor anti-symmetric with respect to that exchange. But the
whole wavefunction remains antisymmetric because the discussed exchange of
electrons is done in all terms at once, i.e. terms the first and second change
into each other with minus sign, the third term just changes sign.

Therefore, for the sake of brevity, let us write the wavefunction in the
following form:

$\Psi\propto\Psi_{12}(r_{1},r_{2},r_{3})|{\rm UUD}\rangle\\\
+\Psi_{23}(r_{1},r_{2},r_{3})|{\rm DUU}\rangle\\\
+\Psi_{13}(r_{1},r_{2},r_{3})|{\rm UDU}\rangle.$ (13)

Here, the subscript of each term denotes the pair of coordinates associated
with antisymmetry of the spatial part. In our simulations we used only the
spatial part of the first term. Let us stress, starting with the spatial
function antisymmetric with respect to the change of the electrons 1 and 2
protects the symmetry of the whole wavefunction.

## References

  * [1] L’Huillier A, Lompre L A, Mainfray G and Manus C 1982 Phys. Rev. Lett. 48(26) 1814–1817
  * [2] L’Huillier A, Lompre L, Mainfray G and Manus C 1983 Phys. Rev. A 27 2503
  * [3] Fittinghoff D N, Bolton P R, Chang B and Kulander K C 1992 Phys. Rev. Lett. 69(18) 2642–2645
  * [4] Kondo K, Sagisaka A, Tamida T, Nabekawa Y and Watanabe S 1993 Phys. Rev. A 48(4) R2531–R2533
  * [5] Walker B, Sheehy B, DiMauro L F, Agostini P, Schafer K J and Kulander K C 1994 Phys. Rev. Lett. 73(9) 1227–1230
  * [6] Staudte A, Ruiz C, Schöffler M, Schössler S, Zeidler D, Weber T, Meckel M, Villeneuve D, Corkum P, Becker A et al. 2007 Phys. Rev. Lett. 99 263002
  * [7] Rudenko A, de Jesus V L B, Ergler T, Zrost K, Feuerstein B, Schröter C D, Moshammer R and Ullrich J 2007 Phys. Rev. Lett. 99(26) 263003
  * [8] Corkum P B 1993 Phys. Rev. Lett. 71(13) 1994–1997
  * [9] Basnayake G, Fernando S, Lee S K, Debrah D A, Stewart G A and Li W 2022 Molecular Physics 120 e1931722
  * [10] Grundmann S, Serov V V, Trinter F, Fehre K, Strenger N, Pier A, Kircher M, Trabert D, Weller M, Rist J, Kaiser L, Bray A W, Schmidt L P H, Williams J B, Jahnke T, Dörner R, Schöffler M S and Kheifets A S 2020 Phys. Rev. Research 2 033080
  * [11] Larimian S, Erattupuzha S, Baltuška A, Kitzler-Zeiler M and Xie X 2020 Phys. Rev. Research 2 013021
  * [12] Henrichs K, Eckart S, Hartung A, Trabert D, Fehre K, Rist J, Sann H, Pitzer M, Richter M, Kang H, Schöffler M S, Kunitski M, Jahnke T and Dörner R 2018 Phys. Rev. A 98 043405
  * [13] Bergues B, Kübel M, Johnson N G, Fischer B, Camus N, Betsch K J, Herrwerth O, Senftleben A, Sayler A M, Rathje T, Pfeifer T, Ben-Itzhak I, Jones R R, Paulus G G, Krausz F, Moshammer R, Ullrich J and Kling M F 2012 Nat. Commun. 3 813
  * [14] Mikaelsson S, Vogelsang J, Guo C, Sytcevich I, Viotti A L, Langer F, Cheng Y C, Nandi S, Jin W, Olofsson A, Weissenbilder R, Mauritsson J, L’Huillier A, Gisselbrecht M and Arnold C L 2020 Nanophotonics 10 117–128
  * [15] Zhong S, Vinbladh J, Busto D, Squibb R J, Isinger M, Neoričić L, Laurell H, Weissenbilder R, Arnold C L, Feifel R, Dahlström J M, Wendin G, Gisselbrecht M, Lindroth E and L’Huillier A 2020 Nat. Commun. 11 5042 ISSN 2041-1723 (Preprint eprint 2005.12008)
  * [16] Parker J S, Smyth E S and Taylor K T 1998 J. Phys. B-At. Mol. Opt. 31 L571
  * [17] Parker J S, Glass D, Moore L R, Smyth E S, Taylor K and Burke P 2000 J. Phys. B-At. Mol. Opt. 33 L239
  * [18] Parker J S, Doherty B J S, Taylor K T, Schultz K D, Blaga C I and DiMauro L F 2006 Phys. Rev. Lett. 96(13) 133001
  * [19] Feist J, Nagele S, Pazourek R, Persson E, Schneider B, Collins L and Burgdörfer J 2008 Phys. Rev. A 77 043420
  * [20] Hao X, Chen J, Li W, Wang B, Wang X and Becker W 2014 Phys. Rev. Lett. 112(7) 073002
  * [21] Lein M, Gross E K U and Engel V 2000 Phys. Rev. Lett. 85(22) 4707–4710
  * [22] Ruiz C, Plaja L, Roso L and Becker A 2006 Phys. Rev. Lett. 96 053001
  * [23] Prauzner-Bechcicki J S, Sacha K, Eckhardt B and Zakrzewski J 2007 Phys. Rev. Lett. 98(20) 203002
  * [24] Prauzner-Bechcicki J S, Sacha K, Eckhardt B and Zakrzewski J 2008 Phys. Rev. A 78 013419
  * [25] Chen S, Ruiz C and Becker A 2010 Phys. Rev. A 82 033426
  * [26] Eckhardt B, Prauzner-Bechcicki J S, Sacha K and Zakrzewski J 2010 Chem. Phys. 370 168 – 174
  * [27] Efimov D K, Maksymov A, Prauzner-Bechcicki J S, Thiede J H, Eckhardt B, Chacón A, Lewenstein M and Zakrzewski J 2018 Phys. Rev. A 98(1) 013405
  * [28] Grobe R, Rzazewski K and Eberly J H 1994 J. Phys. B-At. Mol. Opt. 27 L503–L508
  * [29] Sacha K and Eckhardt B 2001 Phys. Rev. A 64(5) 053401
  * [30] Ho P J and Eberly J 2005 Phys. Rev. Lett. 95 193002
  * [31] Ho P J, Panfili R, Haan S and Eberly J 2005 Phys. Rev. Lett. 94 093002
  * [32] Ho P J and Eberly J 2006 Phys. Rev. Lett. 97 083001
  * [33] Emmanouilidou A and Rost J M 2006 J.Phys. B At. Mol. Opt. 39 4037–4048
  * [34] Ho P J and Eberly J H 2007 Opt. Express 15 1845–1850
  * [35] Zhou Y, Liao Q and Lu P 2010 Optics Express 18(15) 16025
  * [36] Emmanouilidou A 2011 Phys. Rev. A 83(2) 023403
  * [37] Tang Q, Huang C, Zhou Y and Lu P 2013 Optics Express 21(18) 21433–43
  * [38] Katsoulis G P, Hadjipittas A, Bergues B, Kling M F and Emmanouilidou A 2018 Phys. Rev. Lett. 121(26) 263203
  * [39] Peters M B, Katsoulis G P and Emmanouilidou A 2022 Phys. Rev. A 105(4) 043102
  * [40] Jiang H, Efimov D, He F and Prauzner-Bechcicki J S 2022 Phys. Rev. A 105(5) 053119
  * [41] Emmanouilidou A, Peters M B and Katsoulis G P 2023 Phys. Rev. A 107(4) L041101
  * [42] Thiede J H, Eckhardt B, Efimov D K, Prauzner-Bechcicki J S and Zakrzewski J 2018 Phys. Rev. A 98(3) 031401
  * [43] Ruiz C, Plaja L, Vázquez de Aldana J R and Roso L 2003 Phys. Rev. A 68(2) 023409
  * [44] Ruiz C, Plaja L, de Aldana J V and Roso L 2004 Applied Physics B: Lasers and Optics 78 829–833
  * [45] Eckhardt B, Prauzner-Bechcicki J S, Sacha K and Zakrzewski J 2008 Phys. Rev. A 77(1) 015402
  * [46] Prauzner-Bechcicki J S, Efimov D K, Mandrysz M and Zakrzewski J 2021 Journal of Physics B: Atomic, Molecular and Optical Physics 54 114001
  * [47] Efimov D K, Maksymov A, Ciappina M, Prauzner-Bechcicki J S, Lewenstein M and Zakrzewski J 2021 Opt. Express 29 26526–26537
  * [48] Efimov D K, Prauzner-Bechcicki J S, Thiede J H, Eckhardt B and Zakrzewski J 2019 Phys. Rev. A 100 063408
  * [49] Efimov D K, Prauzner-Bechcicki J S and Zakrzewski J 2020 Phys. Rev. A 101 063402
  * [50] Eckhardt B and Sacha K 2006 J. Phys. B: At. Mol. Opt. Phys. 39 3865
  * [51] Arnold V I 2013 Mathematical methods of classical mechanics vol 60 (Springer Science & Business Media)
  * [52] Geltman S 1985 Phys. Rev. Lett. 54(17) 1909–1912
  * [53] Rudenko A, Zrost K, Feuerstein B, de Jesus V L B, Schröter C D, Moshammer R and Ullrich J 2004 Phys. Rev. Lett. 93(25) 253001
  * [54] Liao Q, Winney A H, Lee S K, Lin Y F, Adhikari P and Li W 2017 Phys. Rev. A 96 023401
  * [55] Jiang H, Mandrysz M, Sanchez A, Dura J, Steinle T, Prauzner-Bechcicki J S, Zakrzewski J, Lewenstein M, Fend H, Biegert J and Ciappina M F 2023 New Journal of Physics 25(3) 033002
  * [56] Kübel M, Burger C, Kling N G, Pischke T, Beaufore L, Ben-Itzhak I, Paulus G G, Ullrich J, Pfeifer T, Moshammer R, Kling M F and Bergues B 2016 Phys. Rev. A 93(5) 053422
  * [57] Maxwell A S and Faria C F d M 2016 Phys. Rev. Lett. 116(14) 143001
  * [58] Chen A, Kübel M, Bergues B, Kling M F and Emmanouilidou A 2017 Sci. Rep.-UK 7 7488