Three-electron correlations in strong laser field ionization: Spin induced effects

Strong field processes in the non-relativistic regime are insensitive to the electron spin, i.e. the observables appear to be independent of this electron property. This does not have to be the case for several active electrons where Pauli principle may affect the their dynamics. We exemplify this statement studying model atoms with three active electrons interacting with strong pulsed radiation, using an ab-initio time-dependent Schr\"odinger equation on a grid. In our restricted dimensionality model we are able, for the first time, to analyse momenta correlations of the three outgoing electrons using Dalitz plots. We show that significant differences are obtained between model Neon and Nitrogen atoms. These differences are traced back to the different symmetries of the electronic wavefunctions, and directly related to the different initial state spin components.

Introduction: Interaction of strong laser light with matter still brings new surprises for more than half century as reviewed by [1]. It was realized quite early that the impact of short but strong laser radiation may lead to multiple ionization of atoms [2][3][4][5][6][7] and subsequently discussion on collective, non-sequential or sequential character of electrons emission has begun [3,[8][9][10][11][12][13][14][15][16]. It soon became apparent that the physics involved seems to be more interesting at intermediate intensities when the existence of the characteristic knee in the ionization yield triggered a lively discussion about the importance of the so-called re-scattering mechanism on the non-sequential double ionization [17][18][19]. The theoretical picture came with the simple-man three-step model [20], that provided a simple and intuitively explanation of above-threshold ionization (with electrons leaving the atom with energy far exceeding the threshold [21]) as well as high harmonic generation (HHG) [22] and re-scattering mechanism leading to enhanced double ionization [17][18][19]. The coldtarget recoil-ion-momentum spectroscopy (COLTRIMS) [23][24][25][26][27] allowed one to detect the momenta of an ion and ejected electrons during the ionization event opening the path for the detailed study of both the electron and ion laser-induced dynamics. Further in-depth studies of electron-electron correlation during strong-field evolution led to deeper understanding of ongoing processes and the development of more sophisticated models [28][29][30]. For an additional resolution of ionization delays, the coincidence scheme was assisted by a streaking technique [27,31,32]. An alternative path for accessing realtime dynamics implies a combination of coincidence spec-troscopy with an attosecond interferometry [33,34].
Importantly, recent coincidence experimental schemes [33][34][35][36][37][38] are kinematically complete, that is, momenta, position and detection time are measured for each charged particle explicitly. With increased energy resolution [33,34] one is able to catch fine interference structures that form the quantum fingerprint of the laser-driven electron dynamics. Those recent advances removed the limitations on the number of particles detected: there is no need need to use ionic momenta for recalculating electronic ones.
Although challenging, the experimental arena is thus ready for three-electron coincidence experiments. Specially considering that, in ion-impact induced ionization, the pioneering experiments were already done long time ago [39]. Such an experiment for strong field ionization would definitely open a new page in the understanding of electronic correlations. Its power is of particular importance for studies of recollisional dynamics induced by IR femtosecond laser pulses. In particular, different double ionization channels of multi-electron atoms (affecting different electrons) cannot be resolved in a two-electron coincidence experiment [40]; the three-electron coincidence scheme, on the other hand, potentially could allow one to track the relative contributions of different channels [41].
The theoretical support needed for these prospective experiments are, however, not still available partly because numerical simulations for many electron ionization are extremely challenging. A full quantum-mechanical simulations are, at best, possible for two electrons thanks to the efforts of Ken Taylor's group in Belfast [42,43]. Furthermore, models for three electrons are at their infancy being limited to either purely classical Monte Carlo methods [44][45][46][47][48][49] or to time-dependent orbital approaches which, while successful for description of HHG [50], can hardly provide detailed information on the ionization dynamics and momenta distributions.
Recently, we presented introductory studies of strongfield triple ionization, however, due to the intrinsic complexity of the problem, inevitably we restricted ourselves to only investigate the dependence of ionization yields on the laser pulse amplitude [40,41,51,52]. In the present work, we raise the analysis to a significantly higher level, addressing three electrons momenta distributions. Our method is based on the Eckhardt-Sacha model of reduced dimensionality, derived from the analysis of saddles in the effective adiabatic potential for electrons in an instantaneous electric field [55]. Let us emphasize, that this is a distinctly different approach from the standard Rochester reduction, which restrain the motion of electrons to the laser polarization axis [53,54]. The model: The Hamiltonian of the model studied reads [51]: where r i and p i correspond to the i-th electron's position and momentum, respectively and is a parameter softening the Coulomb singularity. The driving time dependent field F (t) = −∂A/∂t is given via the vector potential with the pulse length T p = 2πn c /ω 0 and carrier-envelope phase (CEP) φ which we put to zero. We limit here to n c = 3 optical cycles case and assume ω 0 = 0.06 a.u., that corresponds to a laser wavelength of about 760 nm. The smoothing parameter = √ 0.83 a.u. and effective electric charge q ee = 1 assures that the ground state energy is equal to the triple ionization potential of Neon, I p = 4.63 a.u., i.e. 126 eV. (for appropriate symmetry of the wavefunction, see below). Starting from the ground state, found with imaginary time propagation method, we directly solve the time-dependent Schrödinger equation (TDSE) on a large 3-dimensional grid using a standard propagation scheme [51]. The algorithm used for simulating momentum distributions is a direct extension of the approach introduced for the two-electron case [56,57]. For the sake of brevity, will describe it very shortly: for each electronic coordinate we distinguish the "bounded motion" region where the numerical solution of the TDSE is used and the "outer" region where each electron is assumed to move freely and the interactions with other electrons and the nucleus are neglected. A distinction between "bounded motion" and "outer" regions is done by setting the threshold distance, defined by a classical distance of no return, r t = F 0 /ω 2 0 from the center of the coordinate system -when a given electron coordinate value exceeds r t it is evolved in the "outer" region; the electron never comes back to the "bounded motion" region. The solutions from the different regions are added coherently (including the possible probability fluxes between the regions). The conceptual idea was created in [56] for two-electron problems (see also [57] for a detailed description) and is generalized here to the three-electron case.
Additionally, and novel for three electrons, a twist comes from the fact that for three electrons one has to be careful about the spatial symmetries of the wavefunction. As is known they are directly related to electronic spin and Pauli exclusion principle. This has been already noticed in the pioneering studies of Li atom ionization [58,59]. In its 1s 2 2s 1 initial state the corresponding Slater wavefunction is a sum of 3 terms (each corresponding to electron permutations) with two electrons with spin up (U) and one down (D). Since the Hamiltonian is spin-independent (neglecting spin-orbit coupling) for non-relativistic strong laser-atom interactions all three components evolve identically and a single combination may be time evolved only. The same would occur for species with electronic ns 2 np 1 configuration, e.g. for Boron or Aluminum, however, in these cases a real multiphoton regime falls within very low frequencies. As a model system with electronic ns 2 np 1 configuration we consider an artificial three-electron atom with the first three ionization thresholds corresponding to Neon [40]. Thus, we will refer to that model as Neon. The ground state wavefunction in the model atom is spatially partially antisymmetric (assuming the spin configuration is (UUD) meaning Up-Up-Down for X, Y and Z electrons): [51]. For comparison, we will consider a model atom that has p 3 configuration with all spins oriented in the same direction (say UUU) and, thus, a totally antisymmetric spatial wavefunction : . The corresponding ionization potentials are set to match the first three ionization potentials of Nitrogen, i.e. 0.52 a.u., 1.61 a.u. and 3.92 a.u. with = √ 1.02 and q ee = √ 0.5. We will refer, therefore, to that model as Nitrogen. The Data Representation: The time evolution of TDSE [60] leads to events in which the three electrons are ionized. Can we extract a useful information on the dynamics from momenta distribution in such a case? We are interested in the triple ionization events, i.e. those events where, once the laser pulse have ceased, all three electrons are in the "outer" region. The 3D wavefunction obtained after the integration of the TDSE allows us to learn about the momentum of outgoing electrons employing the technique of Dalitz plots. These representation is instrumental in particle physics and have been proved to be very useful to disentangle electron-electron and electron-ion correlations in ion atom collisions [39,61]. In its essence, the 3D wavefunction is projected onto the surrounding sphere by integrating over the radial coordinate, and then each octant of this sphere is projected onto a plane, forming an equilateral triangle. Each projection plane is chosen perpendicular to the diagonal line belonging to the particular octant; the type of projection is gnomonic with a tangent equal to the radius of the sphere. The borders between octants are formed by the XY , XZ and Y Z plains, thus the vertices of the Dalitz plots correspond to exclusive motion of the X, Y or Z electron with the other two electrons at rest. We refer to these vertices as X, Y and Z, respectively. For each Dalitz plot, an internal position is defined by a set of distances {π X , π Y , π Z } to the sides of the triangle opposite to the X, Y, Z vertices, correspondingly (see Fig. 6 for illustration): with p X , p Y and p Z being the point coordinates in the momentum space. The particular octant, in turn, is defined by a set of quantities (ξ X , ξ Y , ξ Z ): (Color online) Experimental-type (spin direction averaged) Dalitz plots for triply-ionized state of Neon depict the relative momentum distribution of electrons released in different spatial directions with respect to the external field vector ("+" for positive, "−" for negative). In the notation ( * * * ) the first, second and third signs correspond to X, Y and Z electrons. The interpretation of the position in the plots is presented in Fig. 6.
for the purpose of presentation we substitute −1 by − and +1 by +. One needs to notice that by their very definition, Dalitz plots map ratios of the ejected electrons' momenta to the triple-ion momentum value, while information about momenta absolute values is lost. Still, however, one is able to determine the correlated escape of electrons, e.g. let us consider the plot representing the (+ + +) octant, then two electrons escaping with similar momenta, in the same direction along polarization axis, should contribute to the maxima that are equally distant to two chosen sides of the triangle, i.e. the maxima should appear along each of altitudes forming a structure with trifold symmetry due to electron indistinguishability.
Let us now see the connection between the simulated data and the data of a prospective experiment. While an all-angles experiment allowing the detection of all ejected electrons independently of their direction in 3D space is a standard tool these days, the restricted-dimensionality code we use is not capable of resolving such an information. The numerical data do allow one to determine the positive or negative direction of electronic motion with respect to the laser polarization axis only. On the other hand, in our calculation, we are able to extract spin-resolved information -a knowledge inaccessible in real experiments because the strong field phenomena in the non-relativistic regime are insensitive to the electron spin. Thus, in order to have a compatible set of data, we provide a Dalitz plots transformation to wash out any spin-resolved fingerprint. Such spin-averaged plots would correspond to angular-integrated (with respect to forth-"+" and back-"−" propagation) momentum-resolved 3-electron coincidence experiment.
The corresponding transformation that averages over spin directions is simple: for each experimentallyresolved combination of electrons [(+ + +), (+ + −), (− − +) and (− − −)] one collects all three possible orientation of electrons (vertices XY Z, Y ZX, ZXY ) and then sums them up. As a result, one is left with four plots per one simulation, compare Fig. 2. Such a procedure is necessary for a system consisting of different spin directions. On the contrary, in the p 3 case no averaging is necessary [60]. Results: In order to show the capabilities of an anticipated coincidence experiment on triple strong-field ionization, we simulated the dynamics within three-activeelectron models of Nitrogen and Neon atoms -the examples of systems with spatially fully antisymmetric wavefunction and with partially antisymmetric wavefunction, correspondingly. These two models were previously studied in the context of the ionization yield dependence on the laser field amplitude [40,41,51,52]. In essence the results are as follows. At intermediate laser intensities, when the characteristic knee in the ionization yield is observed, there exist several channels leading to a triple ion, i.e. a sequential ionization, when each electron is liberated independently of others; a direct ionization, when all three electrons are set free at once; and a mixed ionization, when two electrons escape from the atom simultaneously and the third one is ejected separately (either single ionization follows double or vice versa). For atoms with fully antisymmetric wavefunction in their ground state the channels involving any kind of simultaneous escape (either the double ionization in the mixed triple ionization or the direct triple ionization) are suppressed in comparison to the sequential ionization channel. For atoms with a partially antisymmetric wavefunction, however, the mixed ionization channel is not suppressed [41]. Basing on the fact that in the case of two-electron events the non-sequential ionization manifest its correlated character on the electron's momentum distributions in a form of the famous finger-like structures [62,63], we expect to recognize indications of correlated motion in Dalitz plots as well. Furthermore, as the significance of diverse ionization channels in the case of Nitrogen and Neon is different, the correlated motion traces in Dalitz plots for Nitrogen and Neon, respectively, are expected to differ. Fig. 2 shows a set of four Daliz plots obtained for Neon exposed to a pulse with amplitude F = 0.12 a.u. One recognizes in panels (a) and (d), which collect momenta of electrons escaping in the same direction of the polarization axis, a structure centered at the triangle centroid that exhibits a tri-fold symmetry. The structure from panel (a) has three equivalent maxima corresponding to the motion of two electrons with approximately equal momenta and the third one with significantly lower momentum -these maxima point to a mixed ionization: a direct double ionization accompanied by a single ioniza- tion event. There is also a maximum right at the triangle centroid that corresponds to the direct triple escape. One could expect a similar structure in panel (d), as this panel collects momenta of electrons escaping in opposite directions. It is not the case, because we limit ourselves to few-cycle optical cycles. On the contrary, for longer pulses, this should be observed. Panels (b) and (c) collect events in which two electrons are moving in the same direction and the third one in the opposite one -again the correlated motion of the former two electrons manifests itself in a structure that is aligned with and symmetric about the bisector falling from the vertex corresponding to the third electron moving in the opposite direction.
Let us now compare the Dalitz plots for Neon and Nitrogen. We limit our discussion to the octants representing electrons moving in the same direction with respect to the polarization axis. Such plots are shown in Fig. 3. The difference in the observed structures is evident. Although in both panels the observed structures have tri-fold symmetry, in the case of Neon one observes clear maxima on three bisectors and at the triangle centroid (see panel (a)), in striking contrast to the zeros in the distributions along each of the bisectors that are visible for Nitrogen (see panel (b)). The spin symmetry clearly manifests itself in Dalitz plots. For Nitrogen, a spatial totally antisymmetric wavefunction has Ψ = 0 for r i = r j , i, j = {1, 2, 3}; this property is conserved during the transformation to the momentum space, leaving the diagonals empty. Such diagonals map to bisectors during the Dalitz plot creation. In the case of Neon, however, there is a single "zero" bisector in the Dalitz plots of the 8-plot set, thus it does not survive rotations and additions provided to obtain the "experimental" spinaveraged plot. The observed difference in Dalitz plots for Neon and Nitrogen is in line with recalled suppression of correlated escapes for the case of Nitrogen [41]. At the same time both the N and Ne plots share a common feature: a qualitative change of Dalitz plot structures around the field amplitude F = 0.1 a.u. (see Fig. 4)distributions become narrower, centered at, aligned with and symmetric about the bisector falling from vertex corresponding to the electron moving in the opposite direction to the motion of the other two electrons. Comparison with ionization yields curves [41,51], directly implies a clear correlation between such a change and entering into the "knee" regime, where nonsequential ionization processes are instrumental.
There is one more remarkable difference between the Dalitz plots for N and Ne. In the intensity region below the knee, that is, about F < 0.1 in our case, a noticable difference is seen in the general shapes of the structures in (− − +) octants (see Fig. 5). The maxima for N are all concentrated below the triangle centroid O -within the triangle XOY ; the maxima for Ne, on the contrary, are located in the upper part of the plot, within triangles ZOX and ZOY . The effect is robust and repeats for a number of intensities. Importantly, the absence of similar differences in (+ + −) octants (see Fig. 4 panels (a) and (c)) suggests a significant dependence of the signal on the carrier envelope phase (CEP), as should be expected for a few cycle pulse -study of such dependencies are in progress. Indeed, for the low field amplitude triple ionization would almost necessarily imply at least two cycles of ionization with subsequent recollision of an electron or two electrons. Efficiency of these sub-processes is inherently CEP-dependent. From the present data, the two same-direction electrons appear to be fastest in the case of N, while the opposite directions electrons are most probable to have largest velocity for Ne. Summary: The three-electron coincidence scheme appears to be an important tool for studying multiple electron ionization dynamics. We have proposed a possible experimental scheme for detecting electronic momenta after triple ionization of atoms by strong femtosecond laser pulses and have shown how the obtained data can be processed, visualized and analyzed. We have demonstrated how ionization of atoms with different spin symmetry of electronic valence shell can be visualized experimentally. The work in progress will analyze in detail the dependence of the Dalitz plots on the carrier-envelope phase which would provide a more detailed understanding of the process. Our simulations indicate a need for a theoretical model allowing one for a detailed interpretation of the Dalitz plots for three electron strong-field ionization.
We provide here more technical details on extracting momenta distributions from the wavefunctions obtained from the TDSE solution on a spacial grid. An extended description of the algorithm and its implementation will be published elsewhere. We show also exemplary Dalitz plots obtained directly from the wavefunction analyzed, before the spin averaging procedure described in the main text.

Extracting momenta distributions
As mentioned in the main text the ground states within an appropriate symmetry classes (different for Neon and Nitrogen) are obtained by an imaginary time propagation technique. Ground states were computed on a 512 points grid for up to 100 a.u. of length in each direction. While this region seems quite large, it is necessary to correctly represent the tails of each ground state. At each time step, during the imaginary time propagation, the wavefunction was symmetrized back following the rules for a given symmetry class. In this way, the numerical inaccuracies would not drive us, e.g.. to lower lying symmetric (nonphysical) states. For real time propagation the grid has been expanded to 1024 points in each direction.
Typically, in a real time propagation one introduces some complex absorbing potential close to the edges of the grid to avoid spurious reflection. The loss of the norm may then serve as an indicator of the ionization yield. Since we aim at momenta distributions, such a complex potential could artificially affect the electron momenta so we need to use a different approach. One possible way of approaching the problem would be to extend the grid sufficiently far so that the wavefunction never reaches the boundaries during the real time evolution. This is not feasible while keeping a reasonable grid step size. An alternative strategy has been devised in a seminal paper by Lein, Gross and Engel [56] for the two-electron case. We have also followed this scheme in our earlier studies [57]. The present algorithm is a simple (albeit technically hard) extension of this approach to three dimensions (electrons). The idea is as follows, all the space is divided into regions corresponding nominally to the Neutral atom (N), the Single (S), the Double (D) and the Triple (T) ions. During the time evolution one keeps track of the probability fluxes between different regions and that allows one to perform ionization dynamics simulations without absorbing boundaries at all.
The basic physical idea employed in the method of Ref. [56] is that an electron located far away from the parent ion, far beyond the classical point of no return r t = F/ω 2 , would never again propagates close to the nucleus, thus will never again experience an influence of the ionic core as well as the companion electrons. Because of that its evolution is governed, to a good precision, exclusively by the kinetic energy operator p 2 i /2, where: where p i is kinetical momentum,p i is canonical momentum and A(t) is vector potential: F (t) = −∂A(t)/∂t. This way we neglect the influence of the nucleus as well as the brother electrons. Solution of the time-dependent Schrödinger equation in the dimension corresponding to this electron is analytic in the momentum representation; the wavefunction can be stored in momentum representation forever.
Following this idea we define a neutral atom region (N) as ∀i : |r i | < r t ; single ion region (S k ) with k-th electron away from nucleus: ∀i = k : |r i | < r t , |r k | ≥ r t ; double ion region (D kl ) with k-th and l-th electrons away from nucleus: ∀i = k, l : |r i | < r t , |r k | ≥ r t , |r l | ≥ r t ; triple ion region (T ) with all electrons away from nucleus: ∀i : |r i | ≥ r r . Graphical representation of the space division is presented in Fig. 6.
Since we partially neglect interactions, a different effective Hamiltonian is solved in each region. Hamiltonians are obtained by dropping the potential terms. They appear to be small enough due to the large values of the (Color online) Spin sensitive plots (distinguishing the direction to which each spin is ejected -(a,b,c)) are obtained by straightforward projection of the wavefunction; X, Y, Z electrons are attributed with a particular spin. They are further concatenated to (d) after a preliminary rotation applied to each of them: we check that "+" becomes in position Z, so (a) is rotated counter clockwise, (b) -clockwise and (c) does not change its orientation.
distances, taken as a parameter: Note that the interaction with the laser field is hidden in the p i 's. Therefore, in the implementation of the model, we use the velocity gauge Hamiltonian. The Hamiltonians contain the kinetic and (restricted in some regions) potential terms. For calculating the wavefunction evolution during an incremental time step ∆t, the split operator technique is used as in [51], implying that the total solution of the TDSE is expressed as consequent solutions of the TDSE for the kinetic and potential operators. As solutions for these operators are analytical in momentum and position representations, respectively, the Fourier transform is applied to switch between them. In respective representations the actions of the kinetic and potential terms appear as a simple phase factor. If for a given region the i-coordinate the potential operators are absent, the i-dimension of the wavefunction is stored in the momentum representation exclusively. Thus, in the T region the wavefunction is ψ(p X , p Y , p Z ) and directly yields the dominant ingredient of the threedimensional momentum distribution as |ψ(p X , p Y , p Z )| 2 .
The continuity between different regions and the natural flow of the wavefunction through their boundaries is provided by a set of wavefunction-transfer operators, that are applied at each time step. In essence, during the action of such an operator, in each region, starting from (N) to (D)-s and excluding (T), the wavefunction is profiled. In this way, its parts close to the boundaries with the neighboring regions are separated and, after the corresponding partial Fourier transform, are coherently added to the wavefunctions of the neighboring regions.
After ending the time evolution of the wavefunction, i.e. when the external laser field ceased, a mask is applied to all the regions incorporated in the calculations. The mask spatially extracts the ionized part of the wavefunction and transfers it to the neighboring regions: from N to S-regions, from S-regions to D-regions and from Dregions to the T one. Technically, masking is identical to the wave-transfer described previously with the only difference being the shape of the profile: the mask utilizes the boundary between areas corresponding to predominantly bound and predominantly free states. After application of the mask, one ends up with one N region, three S regions, three D regions and eight T regions corresponding to physically neutral, single-ionized, doubleionized and triple-ionized states. Thus, each region can be used for plotting momentum distributions independently: distributions for single ions from S regions, for double ions from D regions and for triple ions from T region with the corresponding one-, two-and three-electron distributions.

Dalitz plots
The construction of the Dalitz plots is described in the main text. As pointed out there, the experiments cannot yield spin-sensitive data, thus the data from simulations have to be spin averaged as shown and described in the main text. Here we present an example of spin-sensitive numerical simulation -compare Fig. 7.
The resulting plots imply that the unique nonsymmetric pattern of Dalitz plots for Neon (partially antisymmetric spin configuration) cannot be accessed with the experimental data that provide a two-fold symmetry for (+−−) and (++−) data for both Neon and Nitrogen -see example in Fig. 7.