Introduction

One of the most notorious departures from classical physics is quantum entanglement, a subtle combination of classical correlation and quantum superposition, which both caused a seismic shift in how the world is viewed, and also provides a resource for quantum computation and metrology1,2,3, including the possibility of entanglement enhanced imaging4,5, with much interest in atomic and molecular physics, e.g., properties of entangled photons from such systems6,7,8. However, the potential for entanglement to optimize attosecond (10−18s) imaging processes is unexplored, and therefore the role and physical insight afforded by entanglement for such processes remains unclear.

Attosecond and strong-field physics deal with processes in matter on the scale of attoseconds9,10. The promise of resolving atomic and molecular electron dynamics on its natural timescale has led to the development of a host of imaging procedures boasting attosecond time resolution; including strong-field initiated methods like high-order high harmonic spectroscopy11, laser-induced electron diffraction12,13, photoelectron holography14,15 and attosecond pump-probe techniques like attosecond streaking16,17 and reconstruction of attosecond harmonic beating by interference of two-photon transitions18,19. Although most of these imaging protocols depend on quantum processes, none explicitly exploit entanglement. Recently, the quantum nature of attosecond processes was exhibited through the generation of non-classical states of light in the laser driving field20,21.

From as early as 199422, there has been a growing interest in the role of entanglement in attosecond processes. Many studies have focused on entanglement between photoelectrons and ions23,24,25,26,27,28,29,30,31, an essential part of understanding decoherence32, while other studies have focused on electron–electron entanglement33,34,35. However, the majority of studies involve the calculation of a continuous variable density matrix (exceptions include studies focused on entanglement involving discrete vibrational states, see, e.g., Refs. 28,30,31) and entanglement measures such as the purity26,29,30,34,35 or von Neumann entropy26,27,29,35. These quantities, when derived from continuous variables, have some drawbacks: (i) They are challenging to compute, and often approximations must be imposed, or in limited cases analytical approximations can be found36. Furthermore, these methods are restricted to pure states, while in strong-field experiments we must always consider mixed states from incoherence averaging. (ii) The physical interpretation of these quantities is difficult, and may not improve understanding of the process. (iii) Direct experimental evidence of the entanglement is often practically impossible, requiring the measurement of incompatible observables, such as momentum and position.

One solution to these difficulties is to use quantized observables. All free particles have such a quantum observable in the form of orbital angular momentum (OAM)37,38. Photons carrying OAM have received significant attention in producing extreme ultraviolet (XUV) high-order harmonics with OAM, see e.g., Refs. 39,40,41,42. Strong-field studies on OAM in photoelectrons include exploiting OAM in rescattering electrons to probe bound state structures43, and recent work44,45,46,47, providing insight into the role of OAM for circularly polarized fields and conservation laws for OAM in strong-field ionization48. Measurement of OAM is a rapidly expanding field, with a host of techniques becoming available for electron beams49,50,51,52,53,54,55,56. Conservation between the initial quantum magnetic number and final OAM, which occurs for systems with rotational symmetry around the quantization axis, may be exploited in strongly-correlated two-electron processes, where entanglement could allow for enhanced photoelectron imaging.

Non-sequential double ionization (NSDI) is a highly correlated two-electron ionization process, the details of which are depicted in Fig. 1. Despite strong electron–electron correlation and rescattering being confirmed in NSDI as early as 200057,58, there has been little focus on the quantum entanglement between the two electrons. This is primarily for the following reasons, (i) classical models have been very successful in modelling NSDI59, (ii) early work suggested entanglement would not play a decisive role. These studies focused on the momentum coordinate parallel to the laser field and found a small degree of quantum correlation33, later it was shown classical correlation was sufficient for field intensities greater than 1014 W/cm260. (iii) Furthermore, computation of NSDI is a very arduous task and computation of continuous variable entanglement measures is even more difficult.

Fig. 1: Diagram of NSDI and the resulting entanglement of OAM.
figure 1

a Orientation of the linearly polarized laser, with polarization e along the z axis taken for the OAM denoted by \({\hat{L}}_{\parallel }\) and the target atom. b NSDI process depicted for the EI and RESI mechanisms. Interaction via the field (e.g. tunnel ionization) is depicted by a dashed line, excitation in the singly charged ion is denoted by dotted lines, while the recollision and OAM sharing is denoted by the yellow spark. The two electron states are defined in the text, following the convention of ref. 66. c The excitation pathways in RESI, which lead to different final OAM states and an entangled superposition. The final states are given by OAM states, \(\left|{l}_{{{{{{{{\rm{e}}}}}}}}},-{l}_{{{{{{{{\rm{e}}}}}}}}}\right\rangle\) as used in (8). The dashed-dotted lines are to denote, which two-electron state the first electron ends up in. The numbers −1, 0, 1 refer to the values of the quantum magnetic number mη of the intermediate excited state, \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\, \eta \right\rangle\).

In this work we address (ii) and (iii), by exploiting the quantized nature of OAM to clearly demonstrate entanglement in NSDI, which we show may occur most easily through the RESI pathway. The use of a quantized degree of freedom enables a simple analysis through the logarithmic negativity61 and entanglement witnesses62,63,64, which enables the inclusion of incoherent effects, as well as a search over a wide range of parameters. The interplay of channels of excitation allows photoelectrons to approach maximally entangled states for some final momenta, which could be investigated as a source of OAM entangled electrons. We show that the entanglement exhibited is robust to incoherent averaging of laser intensities over the focal volume. By decomposing an entanglement witness, we strongly reduce the difficulty of detecting entanglement, by avoiding full tomographic measurements. Furthermore, the OAM entanglement could be used to perform correlated OAM and momentum measurement on the two electrons. Thus, paving the way for a unique kind of entanglement enhanced attosecond imaging technique.

Results

Orbital angular momentum in non-sequential double ionization

A pictorial depiction of non-sequential double ionization (NSDI) is given is Fig. 1. The process follows the three-step mechanism [panel (b)]65: (i) The first electron is removed from the two-electron ground state \(\left|0\right\rangle\) by the laser via tunnel ionization into the state \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,0\right\rangle\) (one electron in the continuum and the other in its ground state). (ii) The continuum electron subsequently undergoes a laser driven recollision with its parent ion. (iii) The energy imparted by the collision allows for two pathways, a second electron is directly ionized in the electron-impact ionization (EI) mechanism, or the second electron is excited, resulting in the state \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,\eta \right\rangle\)—here η is used to label the excited state—and the second electron subsequently ionizes due to the laser field in the recollision with subsequent ionization (RESI) mechanism. In both cases, the final state of the photoelectrons is \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}\right\rangle\), the two electron continuum state. For more information on this notation, see ref. 66, while for reviews of NSDI and these mechanisms, see refs. 67,68. The laser field is polarized in the z direction, in the same direction as the total OAM operator \({\hat{L}}_{|\vert }\), so the laser field cannot change the total OAM since \([{\hat{L}}_{|\vert },\,\hat{H}(t)]=0\), where \(\hat{H}(t)\) is the total Hamiltonian of the system44,48.

We expand the NSDI two-electron continuum wave function in a basis of electron vortex states, the one-electron vortex state is denoted \(\left|{{{{{{{\bf{p}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}\right\rangle\) and a plane wave by \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}}\right\rangle\). Note, for the two-dimensional vectors we use p = (p, p) and for three-dimensional vectors we include a tilde, \(\tilde{{{{{{{{\bf{p}}}}}}}}}=({p}_{|\vert },\,{p}_{\perp },\,\phi )\), written in cylindrical coordinates. Here, \({p}_{\perp }=\sqrt{{p}_{x}^{2}+{p}_{y}^{2}}\) is the radial coordinate, p = pz is the momentum coordinate along the cylindrical axis (parallel to the laser field polarization) and ϕ the azimuthal angle \(\phi=\arctan ({p}_{y}/{p}_{x})\), while le is the topological charge or azimuthal OAM. We will employ atomic units throughout, unless otherwise stated. The spatial representation of the vortex state is given by52

$$\left\langle \tilde{{{{{{{{\bf{r}}}}}}}}}|{{{{{{{\bf{p}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}\right\rangle=\frac{1}{{(2\pi )}^{3/2}}{J}_{{l}_{{{{{{{{\rm{e}}}}}}}}}}({p}_{\perp }{r}_{\perp }){e}^{i\phi {l}_{{{{{{{{\rm{e}}}}}}}}}}{e}^{i{p}_{|\vert }{r}_{|\vert }},$$
(1)

while the momentum representation is

$$\langle \tilde{{{{{{{{\bf{k}}}}}}}}}|{{{{{{{\bf{p}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}\rangle=\frac{{i}^{-{l}_{{{{{{{{\rm{e}}}}}}}}}}{e}^{i\phi ^{\prime} {l}_{{{{{{{{\rm{e}}}}}}}}}}}{2\pi {p}_{\perp }}\delta ({k}_{\parallel }-{p}_{\parallel })\delta ({k}_{\perp }-{p}_{\perp }).$$
(2)

We assume in the asymptotic limit of large distances from the atomic nucleus that the final state of the two electrons in the continuum can be written as a product of vortex states

$$\left|{{{{{{{\bf{p}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime},\,{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}\right\rangle=\left|{{{{{{{\bf{p}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}\right\rangle \otimes \left|{{{{{{{\bf{p}}}}}}}}^{\prime},\,{l}_{e}^{\prime}\right\rangle .$$
(3)

Typically for NSDI we consider momentum measurement and the associated transition amplitude

$$M({\tilde{{{{{\bf{p}}}}}}},{\tilde{{{{{\bf{p}}}}}}}^{\prime})=\mathop{{{{{\mathrm{lim}}}}}}\limits_{t\to \infty}\left\langle {\tilde{{{{{\bf{p}}}}}}},{\tilde{{{{{\bf{p}}}}}}}^{\prime}|\psi (t)\right\rangle,$$
(4)

where \(\left|\psi (t)\right\rangle\) denotes the continuum two-electron wavepacket after the interaction with the external field, and \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}\right\rangle\) denotes the scattering state with asymptotic momenta \(\tilde{{{{{{{{\bf{p}}}}}}}}}\) and \({\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}\). Here, however, we consider the transition amplitude corresponding to OAM measurement, including only the doubly ionized portion of the system, which may be expressed as

$${M}_{{l}_{{{{{\rm{e}}}}}},{l}_{{{{{\rm{e}}}}}}^{\prime}}({{{{{\bf{p}}}}}},\ {{{{{\bf{p}}}}}}^{\prime}) =\mathop{{{{{\mathrm{lim}}}}}}\limits_{t\to \infty}\left\langle {{{{{\bf{p}}}}}},\ {l}_{{{{{{{{\rm{e}}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime},\,{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}|\psi (t)\right\rangle \\ =\frac{{i}^{{l}_{{{{{{{{\rm{e}}}}}}}}}+{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}}{{(2\pi )}^{2}}\iint d\phi d\phi ^{\prime} {e}^{-i\phi {l}_{{{{{{{{\rm{e}}}}}}}}}-i\phi ^{\prime} {l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}M(\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}),$$
(5)

where \({M}_{{l}_{{{{{{{{\rm{e}}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )\) is also the two-dimensional Fourier series coefficient of \(M(\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime})\). Note, we will always construct \({M}_{{l}_{{{{{{{{\rm{e}}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )\) such that electron indistinguishability is accounted for.

The total azimuthal OAM of the two electrons will be conserved at all times

$${l}_{{{{{{{{\rm{e}}}}}}}}}+{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}=m+m^{\prime},$$
(6)

where le and \({l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}\) are the azimuthal OAM of the final vortex state and m and \(m^{\prime}\) are the initial quantum magnetic numbers of the two-electron ground state \(\left|0\right\rangle\). For the recolliding electron in NSDI m ≠ 0 is strongly suppressed, so we consider m = 0. Similarly, for the second electron contributions from \(m^{\prime}\) with opposite signs will destructively interfere, thus we take \(m^{\prime}=0\). The OAM conservation is now trivially \({l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}=-{l}_{{{{{{{{\rm{e}}}}}}}}}\), see Fig. 1(b) and (c). For the RESI mechanism, the second electron leaves from an excited state \(\left|\eta \right\rangle\), and thus we have the additional conservation \({l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}={m}_{\eta }\), thus le = − mη, which again can be seen in Fig. 1(b) and (c). It is also important to consider the role of the ion, if different ionic states were associated with different final OAM states, this would lead to decoherence that could reduce entanglement. However, given the selection rules lead to \(m=m^{\prime}=0\), the residual ion will only feasibly have one final OAM state. Thus, the ion may be traced out without affecting the electron–electron entanglement.

Due to the recollision there is OAM sharing in NSDI, while, in the RESI mechanism, the OAM is tunable via the excited state. Furthermore, the excited electron may occupy a superposition of states69,70,71,72, which means entanglement can emerge in the OAM degree of freedom. From this point onwards, we only consider the excited state populates mη = 0, ± 1. This reduction on the OAM space captures all necessary physics, while limiting the complexity of measurement/ implementation. OAM measurement across considerably larger ranges has been achieved in electron beams, see, e.g., refs. 49,53,55,56. The final two-electron continuum state corresponding to this scenario can be described by,

$$\left|\psi \right\rangle= \iint {d}^{2}{{{{{{{\bf{p}}}}}}}}{d}^{2}{{{{{{{\bf{p}}}}}}}}^{\prime} \left({M}_{1-1}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )\left|{{{{{{{\bf{p}}}}}}}},\,1,\,{{{{{{{\bf{p}}}}}}}}^{\prime},\, {-}1\right\rangle \right.\\ +\left.{M}_{00}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )\left|{{{{{{{\bf{p}}}}}}}},\,0,\,{{{{{{{\bf{p}}}}}}}}^{\prime},\,0\right\rangle+{M}_{-11}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )\left|{{{{{{{\bf{p}}}}}}}},\, {-}1,\,{{{{{{{\bf{p}}}}}}}}^{\prime},\,1\right\rangle \right),$$
(7)

where we consider a time after the end of the pulse and have suppressed t in our notation for convenience, and we will follow this convention throughout the remainder of this work. Here, \(\left|\psi \right\rangle\) is a maximally entangled qutrit if \({M}_{1-1}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )={M}_{00}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )={M}_{-11}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )\). In fact, even for arbitrary M’s, it is possible to show, see Methods, that after the momentum coordinates are traced out the resulting mixed state is always entangled via the positive partial transpose (PPT) criterion73,74. With the condition \(\iint {d}^{2}{{{{{{{\bf{p}}}}}}}}{d}^{2}{{{{{{{\bf{p}}}}}}}}^{\prime} {M}_{k-k}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} ){M}_{k^{\prime} -k^{\prime} }^{*}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )\,\ne\, 0\)\(k,\,k^{\prime} \in [-1,\,1]\). This means OAM entanglement can survive integration over all the momentum coordinates, as long as the probability of excitation to states with two or more values of mη is non-zero. Entanglement in the EI mechanism is also possible, however, in general the final OAM \({l}_{{{{{{{{\rm{e}}}}}}}}}={l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}=0\) dominates, keeping the entanglement associated with EI relatively low, see Supplementary Information for details.

Entanglement measure and witness

In order to quantify and measure entanglement, we consider the density matrix, \(\rho=\left|\psi \right\rangle \left\langle \psi \right|\). To greatly simplify matters we will assume that the experimentalists are only interested in measuring the OAM, thus we will compute the reduced density matrix, tracing over all the continuous momentum components, see (36), leaving an entangled mixed state,

$${\rho }_{{{{{{\rm{OAM}}}}}}}=\mathop{\sum}\limits_{{l}_{{{{{{\rm{e}}}}}}},\,{l}_{{{{{{\rm{e}}}}}}}^{\prime}}{\alpha }_{{l}_{{{{{{\rm{e}}}}}}}{l}_{{{{{{\rm{e}}}}}}}^{\prime}}\left|{l}_{{{{{{\rm{e}}}}}}},\,-{l}_{{{{{{\rm{e}}}}}}}\right\rangle \left\langle {l}_{{{{{{\rm{e}}}}}}}^{\prime},\,-{l}_{{{{{{\rm{e}}}}}}}^{\prime}\right|,$$
(8)

with

$${\alpha }_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}=\int {d}^{2}{{{{{{{\bf{p}}}}}}}}\int {d}^{2}{{{{{{{\bf{p}}}}}}}}^{\prime} {M}_{{l}_{{{{{{{{\rm{e}}}}}}}}}-{l}_{{{{{{{{\rm{e}}}}}}}}}}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} ){M}_{{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}-{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}^{*}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} ).$$
(9)

Note, for all computations the density matrix will be normalized by its trace.

The logarithmic negativity (\({E}_{{{{{{{{\mathcal{N}}}}}}}}}\))61 is a measure of entanglement, that is valid for mixed state systems and exploits the PPT separability criterion73,74 and is easy to compute, which makes it a good choice for our purposes. \({E}_{{{{{{{{\mathcal{N}}}}}}}}}\) is given by

$${E}_{{{{{{{{\mathcal{N}}}}}}}}}={\log }_{2}\left[|\vert {\rho }_{{{{{{{{\rm{OAM}}}}}}}}}^{{T}_{{{{{{{{\rm{A}}}}}}}}}}|{|}_{1}\right],$$
(10)

where the trace norm reads

$$|\vert \rho|{|}_{1}={{{{{{{\rm{Tr}}}}}}}}\left[\sqrt{{\rho }^{{{{\dagger}}} }\rho }\right]$$
(11)

and \({\rho }^{{T}_{{{{{{{{\rm{A}}}}}}}}}}\) is the partial transpose, i.e., the transpose with respect to one subsystem. \({E}_{{{{{{{{\mathcal{N}}}}}}}}}\) may vary between 0 and \({\log }_{2}(3) \ \approx \ 1.58\) for the qutrit system in (7).

This measure exploits the partial transpose (\({\rho }^{{T}_{{{{{{{{\rm{A}}}}}}}}}}\)), which can be related to time-reversal of one of the subsystems. As such, \({\rho }^{{T}_{{{{{{{{\rm{A}}}}}}}}}}\) may not always correspond to a physical system, which manifests by negative eigenvalues. The partial transpose preserves positive eigenvalues for a separable system, but not for entangled systems. There exist entangled systems, where the partial transpose does not yield negative eigenvalues. These will yield a logarithmic negativity of zero. Correspondingly, \({E}_{{{{{{{{\mathcal{N}}}}}}}}}={\log }_{2}(2{{{{{{{\mathcal{N}}}}}}}}-1)\), where \({{{{{{{\mathcal{N}}}}}}}}\) (the negativity) is given by the absolute value of the sum of the negative eigenvalues of \({\rho }^{{T}_{{{{{{{{\rm{A}}}}}}}}}}\). \({E}_{{{{{{{{\mathcal{N}}}}}}}}}\) gives an upper bound to the distillable entanglement, i.e., quantifying the number of copies of ρ required to transform it into a maximally entangled state. Other commonly used entanglement measures, such as entropy of entanglement and the purity, exploit that the reduced density matrix (tracing over one of the subsystems) cannot easily be used if the input state is mixed, for more information, see the Supplementary Information.

To measure if there is entanglement, one approach is to use an entanglement witness, which can distinguish a subset of entangled states as non-separable. They can commonly be associated with an observable for which the expectation value is negative for some entangled states, i.e., for a witness \({{{{{{{\mathcal{W}}}}}}}}\) and state ρ, the condition \({{{{{{{\rm{Tr}}}}}}}}\left[{{{{{{{\mathcal{W}}}}}}}}\rho \right] \, < \, 0\) implies ρ is entangled. It can be shown that

$${{{{{{{\mathcal{W}}}}}}}}(\theta )=\frac{1}{d}{\mathbb{1}}-|\nu (\theta )\rangle \langle \nu (\theta )|$$
(12)

with

$$\vert \nu (\theta )\rangle=\frac{1}{\sqrt{d}}\mathop{\sum }\limits_{{l}_{{{{{{{{\rm{e}}}}}}}}}=-1}^{1}{e}^{i\theta {l}_{{{{{{{{\rm{e}}}}}}}}}}|{l}_{{{{{{{{\rm{e}}}}}}}}},\,-{l}_{{{{{{{{\rm{e}}}}}}}}}\rangle,$$
(13)

is a valid entanglement witness. Here, the dimension d = 3. This witness is useful as the state \(\left|\nu (\theta )\right\rangle\) mirrors the anti-correlated (\({l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}=-{l}_{{{{{{{{\rm{e}}}}}}}}}\)) entangled superposition of the OAM and also contains a system-dependent tuning parameter θ, which can be set depending on the target and field parameters to enhance the detectability. The parameter θ can be tuned depending on the phase information in the state, in particular depending on \({\theta }_{\alpha } {:} {=} \arg ({\alpha }_{01})\), where α01 represents the coherence between \(\left|0,\,0\right\rangle\) and \(\left|\pm 1,\,\mp 1\right\rangle\) in the density matrix ρOAM given by (8). We will set θ to two extremes θ = 0 or π and discuss the correspondence with θα.

Optimizing entanglement

In Fig. 2 we show the logarithmic negativity and the expectation value of the entanglement witness for many targets and laser parameters. We do this using the strong-field approximation (SFA), the details of which are described in the Methods section and refs. 70,71. The SFA is an approximate method, but it captures the basic OAM correlation and is rapid to compute, so is well-suited to a broad parameter search. The parameters are chosen such that the return energy of the first electron, ≈3.17Up, is not above the ionization potential of the second, to ensure that the RESI mechanism is dominant. Here, Up is the ponderomotive or quiver energy of the electron in the laser field. Also, we stay approximately in the tunneling regime, with the Kelydsh75 parameter \(\gamma=\sqrt{{I}_{{{{{{{{\rm{p}}}}}}}}}/(2{U}_{{{{{{{{\rm{p}}}}}}}}})}\) spanning 0.87–1.20. Note that the entanglement is not dependent on tunnel ionization, i.e., there will still be entanglement if the second electron is ionized through multiphoton absorption.

Fig. 2: Extensive search over targets and parameters.
figure 2

a Quantifies the entanglement using the logarithmic negativity (\({E}_{{{{{{{{\mathcal{N}}}}}}}}}\) of (10)). b and c Expectation value of the entanglement witnesses (\(\langle {{{{{{{\mathcal{W}}}}}}}}(0)\rangle\) and \(\langle {{{{{{{\mathcal{W}}}}}}}}(\pi )\rangle\), see (12)), respectively. Here, \(\langle {{{{{{{\mathcal{W}}}}}}}}\rangle : ={{{{{{{\rm{Tr}}}}}}}}[{{{{{{{\mathcal{W}}}}}}}}{\rho }_{{{{{{{{\rm{OAM}}}}}}}}}]\). Blue text indicates λ = 400 nm and red λ = 800 nm. d, e, f Estimated contribution of intermediate excited states for three extremal cases Mg, Be and Ne, respectively. The field values are given on the panels. The contribution is computed by taking the peak value of the momentum dependent probability of each channel individually. The values are normalized by the largest channel and in arbitrary units. We show \({E}_{{{{{{{{\mathcal{N}}}}}}}}}\), \(\langle {{{{{{{\mathcal{W}}}}}}}}(0)\rangle\), \(\langle {{{{{{{\mathcal{W}}}}}}}}(\pi )\rangle\) and the phase \({\theta }_{\alpha }=\arg ({\alpha }_{10})\), see (8). The excited state labels, at the top and bottom of the panels, are in the format \(|\eta \rangle=\vert n{\ell }_{{m}_{\eta }}\rangle\), where n is the principle quantum number, is the angular quantum number and mη is the magnetic quantum number.

Immediately from Fig. 2a, it is clear that magnesium and beryllium have the highest logarithmic negativity, while the nobel gases argon and neon have the lowest. We may write the logarithmic negativity in terms of the coefficients \({\alpha }_{{l}_{{{{{{{{\rm{e}}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}\), \({E}_{{{{{{{{\mathcal{N}}}}}}}}}={\log }_{2}({\alpha }_{00}+4|{\alpha }_{10} \vert+4{\alpha }_{11})\), due to indistinguishability of the electrons, which means αij = αij = αij = αij. The coefficients α00 and α11 are determined by the population of intermediate states with mη = 0 and mη = ± 1, respectively, while the term α10 is a measure of the coherence between these states. In the bottom row, panels (d)–(f), the contribution of channels via different excited states is estimated for the two highest and lowest cases. In the case of the nobel gases, states with mη = 0 dominate, thus reducing the entanglement as a single OAM prevails, see panel (f). This selectivity provides information on recollision dynamics: For nobel gases, the second electron initial p-state (with \(m^{\prime}=0\)) is aligned along the OAM axis and transitions to excited states which keep this alignment are more probable. For magnesium and beryllium in panels (d) and (e), there is a more balanced superposition across mη = − 1, 0, 1, the initial s-state of these targets is spherically symmetric, so there will be no directional preference for the excited state, leading to higher logarithmic negativity and a higher degree of entanglement. Thus, the logarithmic negativity is a direct probe of excited state population and geometry.

In the panels (b) and (c), the expectation values of \({{{{{{{\mathcal{W}}}}}}}}(0)\) and \({{{{{{{\mathcal{W}}}}}}}}(\pi )\) are shown. Values below and above zero correspond to entangled states that may and may not be distinguished from separable states, respectively. The witness \({{{{{{{\mathcal{W}}}}}}}}(\pi )\) outperforms \({{{{{{{\mathcal{W}}}}}}}}(0)\), which can be traced back to the phase θα, which varies with the target and field parameters. The expectation value of the witnesses can be expressed in terms of the coefficients αij, \(\langle {{{{{{{\mathcal{W}}}}}}}}(0/\pi )\rangle=-\frac{2}{3}({\alpha }_{11}\pm 2|{\alpha }_{10}|\cos ({\theta }_{\alpha }))\). Thus, the witness probes the coherence term, α10 and its phase θα. The witness \({{{{{{{\mathcal{W}}}}}}}}(\pi )\) improves the detectability of entanglement for λ = 800 nm, while less difference is seen for λ = 400 nm. E.g. beryllium in panel (e) θα ≈ π so the witness \({{{{{{{\mathcal{W}}}}}}}}(\pi )\) is much lower than \({{{{{{{\mathcal{W}}}}}}}}(0)\). Furthermore, for magnesium in panel (d) at λ = 400 nm θα ≈ 0.5π, which is why \({{{{{{{\mathcal{W}}}}}}}}(0)\) and \({{{{{{{\mathcal{W}}}}}}}}(\pi )\) give nearly identical results. Thus, the witness also provides key information on coherence between \(\left|0,\,0\right\rangle\) and \(\left|\pm 1,\,\mp 1\right\rangle\) in ρOAM via α10, which is not usually accessible given that α10 is not an observable. Thus, tuning the witness via the parameter θ, enables efficient entanglement detection for specific targets and fields, as well allows the phases between channels of excitation to be determined.

In Fig. 3, we plot the correlated momentum distribution \({{\Omega }}({p}_{|\vert },\,{p}_{|\vert }^{\prime})\) over p and \({p}_{|\vert }^{\prime}\), where we have integrated over the components perpendicular to the laser field polarization.

$${{\Omega }}({p}_{|\vert },\,{p}_{|\vert }^{\prime})\propto \iint d{p}_{\perp }d{p}_{\perp }^{\prime}{p}_{\perp }{p}_{\perp }^{\prime}|M(\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}){|}^{2}.$$
(14)

Interferences can be seen, which is a hallmark of the superposition of excited states69,70,71,72,76. In the second column, we plot momentum-dependent logarithmic negativity \({{{{{{{{\rm{E}}}}}}}}}_{{{{{{{{\mathcal{N}}}}}}}}}({p}_{|\vert },\,{p}_{|\vert }^{\prime})\). Here, the logarithmic negativity is computed at specific momentum values, see caption, quantifying the entanglement between the photoelectrons given a specific final momentum.

Fig. 3: Momentum dependent probability, logarithmic negativity, and density matrices.
figure 3

a and d The correlated momentum distribution \({{\Omega }}({p}_{|\vert },\,{p}_{|\vert }^{\prime})\) for magnesium and beryllium, respectively, with respect to the parallel momentum coordinates p and \({p}_{|\vert }^{\prime}\). We have integrated over the components perpendicular to the laser field polarization, given by (14). A logarithmic scale is given, normalized to the peak value in arbitrary units. The target is given in the title, while the field parameters are listed in the top left. b and e Momentum-dependent logarithmic negativity \({{{{{{{{\rm{E}}}}}}}}}_{{{{{{{{\mathcal{N}}}}}}}}}({p}_{\vert \vert},\,{p}_{\vert|}^{\prime})\) for the same respective targets. The perpendicular components are taken to be (nearly) zero. c and f Display the density matrices, ρOAM (8), for the indicated targets. The non-zero complex element of ρOAM, are colored blue and represented by the phase \(\theta=\arg ({\alpha }_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}})\) (recall that \({\theta }_{\alpha }=\arg ({\alpha }_{01})\)) and the modulus \(r=\vert {\alpha }_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}} \vert\). The bra and ket labels, c and f, correspond to the states \(\left|{l}_{{{{{{{{\rm{e}}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}\right\rangle\), where the le and \({l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}\) is the OAM of the two photoelectrons.

We may write, \({{{{{{{{\rm{E}}}}}}}}}_{{{{{{{{\mathcal{N}}}}}}}}}({p}_{|\vert },\,{p}_{|\vert }^{\prime})={\log }_{2}(|{M}_{-11}+{M}_{00}+{M}_{1-1}{|}^{2})\), where we have used \({M}_{ij}\equiv {M}_{ij}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} ){|}_{{p}_{\perp }={p}_{\perp }^{\prime}=\delta p}\), see (5). This coherent mix of OAM channels reveals the importance of their relative phases. Such phases are determined by an interplay between recollision dynamics and ionization from the excited state. As such, \({{{{{{{{\rm{E}}}}}}}}}_{{{{{{{{\mathcal{N}}}}}}}}}({p}_{|\vert },\,{p}_{|\vert }^{\prime})\) varies considerably, forming entanglement fringes, which reveal the phase between the OAM channels that follow the fringes in the momentum distribution. For \({{{{{{{{\rm{E}}}}}}}}}_{{{{{{{{\mathcal{N}}}}}}}}}({p}_{|\vert },\,{p}_{|\vert }^{\prime})\) of magnesium, panel (b), there are entanglement maxima either side of the p and \({p}_{|\vert }^{\prime}\) axes, with diagonal fringes and a central ring of high logarithmic negativity. This behavior can be explained through the excitation channels 4s0, 4f0 and 4f±1, which contribute the most, see Fig. 2. In Fig. 3 there are peaks in the entanglement when there is a balanced superposition of channels with mη = −1, 0 and 1. This occurs where 4f0 and 4s0 are maximum, which is adjacent to the axis and in the center, respectively. For \({{{{{{{{\rm{E}}}}}}}}}_{{{{{{{{\mathcal{N}}}}}}}}}({p}_{|\vert },\,{p}_{|\vert }^{\prime})\) of beryllium, panel (e), has the same entanglement maxima either side of the p and \({p}_{|\vert }^{\prime}\) axes as magnesium. This arises from the states 4d0, 4d±1, 4f0 and 4f±1, which leads to the same effect due to on-axis nodes of 4d0 and 4f0. In this case, the mixing of dominant channels with different l leads to the large peaks on the diagonal (around \({p}_{|\vert },\,{p}_{|\vert }^{\prime} \approx \pm 0.8\) a.u.), where the photoelectrons have large correlated momentum.

Through the interplay of excitation channels, we have identified three cases where the two photoelectrons approach the maximum logarithmic negativity. (i) Combination of three channels with equal and mη = ±1, 0, which results in off-axis maxima, i.e., fast and slow entangled photoelectrons. (ii) The combination of an s state with higher channel with mη = ±1, leads to central maxima or two slow entangled photoelectrons. (iii) The combination of two sets of differing channels each with mη = ±1, 0, this leads to two fast pairs of entangled photoelectrons, with a correlated direction. The different types of pairs of highly entangled photoelectrons could be useful as imaging probes accessing different momentum regions, or alternatively as a source of OAM entangled electrons, which can be optimized by tuning the interplay between different channels of excitation, for instance with tailored fields.

In Fig. 3, we display ρOAM for beryllium and magnesium targets. The phases of the complex entries of ρOAM are written on each element, it is clearly only θα, the phase between \(\left|{\pm} 1,\,\mp 1\right\rangle\) and \(\left|0,\,0\right\rangle\), that plays a role. The closer the tuning parameter θ is to θα the more effective the entanglement witnesses. The measure \({E}_{{{{{{{{\mathcal{N}}}}}}}}}\), on the other hand, is independent of this phase and is determined by the relative magnitude of the non-zero elements of ρOAM. We can see in the figure, it is the elements related to \(\left|{\pm} 1,\,\mp 1\right\rangle\) and \(\left|0,\,0\right\rangle\) that are most reduced. The coherence between these states is reduced after tracing over the momentum coordinates. The coherence between \(\left|{\pm} 1,\,\mp 1\right\rangle\) and \(\left|\mp 1,\,\pm 1\right\rangle\) is robust, however, as the channels that leads to these states such as 4d±1 are degenerate, so will have the same behavior over momentum. Thus, the final mixed state keeps a reasonably high logarithmic negativity.

Measurement considerations

There are additional incoherent averaging effects leading to mixed states, which will primarily be in the form of fluctuations and averaging over the carrier envelope phase and the laser intensity. For simplicity, we are assuming long enough pulses for the former to not be an issue. The latter takes place through intensity variation over focal volume (focal averaging) as well as intensity variation from shot-to-shot. In the Supplementary Information, we compute the focal averaged momentum distributions and density matrices, following the procedure set out in refs. 71,77. We find that this does not have a significant effect on the overall entanglement, as the coherence between the states relating to mη = ±1 is robust to intensity variation.

Aside from incoherent averaging, ion–photoelectron entanglement would lead to decoherence. However, the strict selection rules, prevent decoherence by ensuring there is only one final state of the ion. One could speculate whether a multielectron treatment is required, beyond the 2e + ion model employed here? The remaining core electron in the cases of beryllium and magnesium are tightly bound and unlikely to play a role. In argon these effects may have to be considered. It is possible this would provide additional decoherence channels to final OAM states, which would reduce the entanglement. Full inclusion of the Coulomb potential is also an important consideration. Recent work demonstrated the inclusion of the Coulomb potential introduced recolliding trajectories for the second electron78. However, this treatment uses the same S-matrix treatment for the recollision-excitation step, allowing for entanglement via the same mechanism. The ion could couple to the photoelectrons via its angular momenta. Given the large mass of the ion, this coupling is expected to have a vanishing effect.

In its current form, (12), the expectation value of the witness cannot be easily measured as it would require a combined measurement on both particles. However, any witness of this type may be decomposed79 into a series of local measurements performed on each particle separately62,63. Practically, this will mean doing multiple experiments to compute expectation values and then combining these results with suitable weights determined by the decomposition. For the specific values of θ = 0 and θ = π, the witness decomposition is

$${{{{{{{\mathcal{W}}}}}}}}(0/\pi )= \bigg[\frac{8}{3}{{\mathbb{1}}}^{\otimes 2}\,-\,{\lambda }_{3}^{\otimes 2}\,-\,2({\lambda }_{4}^{\otimes 2}\,+\,{\lambda }_{5}^{\otimes 2})\,+\,{\lambda }_{8}^{\otimes 2}\,+\,\sqrt{3}({\lambda }_{3}\,\otimes \,{\lambda }_{8}\,+\,{\lambda }_{8}\,\otimes \,{\lambda }_{3})\\ \pm 2({\lambda }_{1}\,\otimes \,{\lambda }_{6}\,+\,{\lambda }_{6}\,\otimes \,{\lambda }_{1}\,+\,{\lambda }_{2}\,\otimes \,{\lambda }_{7}\,+\,{\lambda }_{7}\,\otimes \,{\lambda }_{2})\bigg]/12,$$
(15)

where − is taken for θ = 0 and + is taken for θ = π. Here, \({\mathbb{1}}\) corresponds to the identity, i.e., performing no measurement, while λi are matrices that label 8 different single-particle measurements. In the case of a two-level system (qubit/ SU(2)) the local measurement may correspond to the Pauli matrices plus the identity, which have a simple geometrical interpretation. For our three-level system (qutrit/ SU(3)) we have used the eight Gell-Mann matrices. As an example, \({\lambda }_{3}={L}_{|\vert }^{(1)}\) corresponds to a direct OAM measurement, while λ4 corresponds to measuring in the basis (\(\left|{-} 1\right\rangle \pm \left|1\right\rangle\)) and \(\left|0\right\rangle\). For a full definition of the Gell–Mann matrices, see the Methods section or ref. 80.

We can use the decomposition (15) to determine how it could reduce the total number of measurements required. There are 11 combinations of λxλy. However, we can discount \({{\mathbb{1}}}^{\otimes 2}\) as this corresponds to doing nothing. Furthermore, [λ3, λ8] = 0 (i.e., they share eigenstates) so all parts containing only λ3 and λ8 can be determined with the same measurement. Additionally, measurements of the form λxλy and λyλx may be collected simultaneously, as we do not distinguish the two electrons. This leaves five measurement settings, a large reduction to the 81 combinations of measurement for a general qutrit system.

Entanglement-enhanced attosecond imaging

In the non-perturbative regime, we have shown entanglement can be generated from a recollision process, which protects the electrons from decoherence with the ion. Now we will examine the possibility of exploiting this entanglement for enhanced attosecond imaging, to access information beyond that available in a classically correlated system. Quantum enhancement in attoscience is relatively unexplored, but coherence has been shown to improve measurement precision81. To exploit the entanglement, we may measure each electron in a different basis, e.g., OAM (vortex state) and momentum (plane wave) measurements. The probability of this mixed measurement is \({{{{{{{\rm{Prob}}}}}}}}({{{{{{{\bf{p}}}}}}}},\,l,\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime})=\vert {M}_{l-l}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} ){|}^{2}\), where due to the OAM correlation, we may disentangle channels of excitation. This is in contrast to the purely momentum-based measurement, which leads to \({{{{{{{\rm{Prob}}}}}}}}(\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime})=\vert \mathop{\sum }\nolimits_{l=-1}^{1}{M}_{l-l}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} ){|}^{2}\); note for a classically correlated system, this would be an incoherent sum. We may combine such measurements, e.g.,

$$\chi (\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime})\equiv {{{{{{{\rm{Prob}}}}}}}}(\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime})-\mathop{\sum}\limits_{l}{{{{{{{\rm{Prob}}}}}}}}({{{{{{{\bf{p}}}}}}}},\,l,\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}).$$
(16)

In a separable classically correlated system \(\chi (\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime})=0\), while in our entangled system we have \(\chi (\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime})={\sum }_{l\ne l^{\prime} }{M}_{l-l}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} ){M}_{l^{\prime} -l^{\prime} }^{*}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )\), giving access to the correlation terms. Mixed measurements in other bases, such as those defined by the Gell–Mann matrices, would provide yet further correlated measurements. Switching between a coherent and incoherent sum bears similarity to entangled double slit interference observed in doubly ionized H282. By measuring the OAM, we are able to (i) disentangle the channels of excitation with differing quantum magnetic numbers, not possible in an uncorrelated system, and (ii) access coherence terms, yielding phase information not accessible in a classically correlated system. Results, which include the Coulomb potential in the electron propagation for RESI78, show that recollision of the second electron plays an important role. Thus, trajectories will interfere with differing amounts of interaction with the Coulomb potential, providing holographic interferences15. The combination of correlated mixed measurement, exploiting OAM entanglement and photoelectron holography provides two ways to extract phase information, making a powerful attosecond imaging tool.

Discussion

In this work we use the correlated process of non-sequential double ionization (NSDI), as a backbone, to generate entanglement in the orbital angular momentum (OAM) of two photoelectrons. The photoelectron OAM in NSDI has some general rules, which show that entanglement occurs in the recollision with subsequent ionization (RESI) mechanism of NSDI, prevalent at lower intensities60. The OAM in NSDI has not been explored before, and previous entanglement studies have focused on the momentum coordinates parallel to the laser field33. We find that OAM entanglement can be maximal in specific momentum regions, and that this is controlled by the interplay of channels of excitation, which may be investigated as a source of OAM-entangled electrons. Furthermore, the entanglement is robust as it survives incoherent averaging over the focal volume and it is not generated by the symmetrization resulting from the indistinguishability of the photoelectrons83.

The use of the OAM has many benefits, firstly, the entanglement is simply understood as a consequence of angular momentum sharing during recollision, coupled with a superposition over OAM states due to contribution from excitation-channels with a differing quantum magnetic number [Fig. 1c]. Secondly, the quantization enables a clear and simple analysis using the logarithmic negativity, which enables an extensive search over targets and parameters to maximize the entanglement, where it is clear ideal targets are those with two s-state valence electrons, given that this will promote a balanced superposition across OAM states [Fig. 2]. Thirdly, the reduction in computational difficulty allows the density matrix and momentum-dependent logarithmic negativity to be computed [Fig. 3] and enables incoherent averaging. Finally, we can construct an entanglement witness, which may be decomposed into local measurements, avoiding full state tomography or the measurements of incompatible continuous observables, like position and momentum, reducing the difficulty of experimental implementation for the detection of entanglement.

A key question is how to perform such an experiment? The measurement of OAM for electron vortex beams, has received a great deal of attention and a variety of methods have been demonstrated49,50,51,52,53,54,55,56. This includes diffractive methods such as a fork hologram49,50 or Dammann vortex grating55 or conformal mapping techniques using phase plates53 or electrostatic fields56. This arsenal of techniques is capable of detecting an OAM range far in excess of what is required here. The latter electrostatic OAM sorter56 boasts high efficiency and was retrofitted as an OAM analyzing element to a transmission electron microscope. Thus, such an element could be conceivably added to a typical reaction microscope (ReMi)84,85 employed for NSDI—a ReMi measures correlation of the ion and the two electrons, thus, it is already well-suited for studying entanglement86. Beyond the measurement of OAM, there will be some difficulty in using alkaline earth metals in experiment over the typical nobel gas targets. Works87,88 have employed a range of alkaline earth metals in studies on ion rates, even at a longer wavelength of 2000 nm88, however, a ReMi-style measurement of the electrons’ momentum would be challenging. An alternative route could be to investigate larger parameter ranges and tailored fields in order to boost the entanglement in the noble gases or other simpler targets. This search would be aided if the theoretical model included the electron impact mechanism to accurately address the effect on overall entanglement.

The entangled OAM of the photoelectrons in NSDI may be exploited in various ways, e.g., through interferometric schemes exploiting tailored fields, as a source of OAM-entangled electron pairs, or correlated measurements of OAM and momenta. An interesting consideration is the spin of the electrons. Previous work observed spin polarization from the ground states of heavier targets, see e.g.,89,90. However, for the targets we have employed this effect will be small. The spin–orbit coupling during strong-field dynamics, e.g., during recollision, may play some role but has been generally neglected in strong-field studies. OAM in attosecond processes provides a rich burgeoning research area, to help achieve the aim of imaging and controlling matter on ultrafast times scales and photoelectron OAM-entanglement plays an important role in achieving this goal. Beyond this, like ref. 20, it further demonstrates the fundamental non-classicality of such processes.

Methods

Strong-field approximation

The wave function for NSDI, which describes both electron-impact ionization (EI) and recollision with subsequent ionization (RESI) can be written using an SFA flavor ansatz66 in the following way

$$\left|\psi (t)\right\rangle= a(t)\left|0\right\rangle+\int {d}^{3}\tilde{{{{{{{{\bf{p}}}}}}}}}\,b(\tilde{{{{{{{{\bf{p}}}}}}}}},\,t)\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,0\right\rangle+\mathop{\sum}\limits_{\eta }\,\,\int {d}^{3}\tilde{{{{{{{{\bf{p}}}}}}}}}\,c(\tilde{{{{{{{{\bf{p}}}}}}}}},\,\eta,\,t)\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,\eta \right\rangle \\ +\iint {d}^{3}\tilde{{{{{{{{\bf{p}}}}}}}}}{d}^{3}{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}\,d(\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime},\,t)\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}\right\rangle,$$
(17)

where \(\left|0\right\rangle\) is the two electron ground state, \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,0\right\rangle\) corresponds to one electron in the continuum and the other in its ground state, \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,\eta \right\rangle\) is the same as the latter with the second electron excited and \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}\right\rangle\) is the two electron continuum state.

The final symmetrized transition amplitude is related to (17) via

$${M}^{{{{{\rm{RESI}}}}}}({\tilde{{{{{\bf{p}}}}}}},\,{\tilde{{{{{\bf{p}}}}}}}^{\prime}) =\mathop{{{{{\mathrm{lim}}}}}}\limits_{t\to \infty}d({\tilde{{{{{\bf{p}}}}}}},\,{\tilde{{{{{\bf{p}}}}}}}^{\prime},\,t)\\ =\frac{1}{\sqrt{2}}\left({M}_{{{{{{{{\rm{unsym}}}}}}}}}^{{{{{{{{\rm{RESI}}}}}}}}}({\tilde{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime})+{M}_{{{{{{{{\rm{unsym}}}}}}}}}^{{{{{{{{\rm{RESI}}}}}}}}}({\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime},\,{\tilde{{{{{{{\bf{p}}}}}}}}})\right),$$
(18)

where \({M}_{{{{{{{{\rm{unsym}}}}}}}}}(\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime})\) is the unsymmetrized transition amplitude. This is valid for an initial singlet state, as considered here. In the SFA, in atomic units, for the RESI mechanism and using the assumptions listed in refs. 70,71, the unsymmetrized transition amplitude can be written as

$${M}_{{{{{{{{\rm{unsym}}}}}}}}}^{{{{{{{{\rm{RESI}}}}}}}}}({\tilde{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime})=\mathop{\sum}\limits_{\eta }\int {d}^{3}t\int {d}^{3}{\tilde{{{{{{{\bf{k}}}}}}}}}{V}_{\tilde{{{{{{{{\bf{p}}}}}}}}}{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime},{\tilde{{{{{{{\bf{p}}}}}}}}}\eta }{V}_{\tilde{{{{{{{{\bf{p}}}}}}}}}\eta,{\tilde{{{{{{{\bf{k}}}}}}}}}0}{V}_{\tilde{{{{{{{{\bf{k}}}}}}}}}0,0}\exp [iS({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime},\,{{{{{{{\bf{k}}}}}}}},\,t,\,t^{\prime},\,t^{\prime\prime} )],$$
(19)

where

$$\int {d}^{3}t\equiv \int\nolimits_{-\infty }^{\infty }dt\int\nolimits_{-\infty }^{t}dt^{\prime} \int\nolimits_{\infty }^{t^{\prime} }dt^{\prime\prime}$$
(20)

and

$$ S({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime},\,{{{{{{{\bf{k}}}}}}}},\,t,\,t^{\prime},\,t^{\prime\prime} ) \\ \quad={I}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{10}}}}}}}}}t^{\prime\prime}+{I}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{20}}}}}}}}}t^{\prime}+{I}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{\eta }}}}}}}}}(t-t^{\prime} )-\int\nolimits_{t^{\prime\prime} }^{t^{\prime} }\frac{{[{{{{{{{\bf{k}}}}}}}}+{{{{{{{\bf{A}}}}}}}}(\tau )]}^{2}}{2}d\tau \\ \qquad -\int\nolimits_{t^{\prime} }^{\infty }\frac{{[{{{{{{{\bf{p}}}}}}}}+{{{{{{{\bf{A}}}}}}}}(\tau )]}^{2}}{2}d\tau -\int\nolimits_{t}^{\infty }\frac{{[{{{{{{{\bf{p}}}}}}}}^{\prime}+{{{{{{{\bf{A}}}}}}}}(\tau )]}^{2}}{2}d\tau$$
(21)

denotes the semiclassical action and \({I}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{10}}}}}}}}}\), \({I}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{20}}}}}}}}}\) and \({I}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{\eta }}}}}}}}}\) are the one-electron ionization potentials corresponding to removing a bound electron from \(\left|0\right\rangle\), \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,0\right\rangle\) and \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,\eta \right\rangle\), respectively. Note that for three-dimensional vectors we include a tilde \(\tilde{{{{{{{{\bf{p}}}}}}}}}=({p}_{|\vert },\,{p}_{\perp },\,\phi )\), while for two-dimensional vectors we do not p = (p, p), where p = pz, \({p}_{\perp }=\sqrt{{p}_{x}^{2}+{p}_{y}^{2}}\) and ϕ is the azimuthal angle. The prefactors are given by

$${V}_{\tilde{{{{{{{{\bf{k}}}}}}}}}0,\,0} = \langle \tilde{{{{{{{{\bf{k}}}}}}}}}(t^{\prime\prime} ),\,0|V|0\rangle \\ =\frac{1}{{(2\pi )}^{3/2}}\int {d}^{3}\tilde{{{{{{{{\bf{r}}}}}}}}}\,V(\tilde{{{{{{{{\bf{r}}}}}}}}}){e}^{-i\tilde{{{{{{{{\bf{k}}}}}}}}}(t^{\prime\prime} )\cdot \tilde{{{{{{{{\bf{r}}}}}}}}}}{\psi }_{10}(\tilde{{{{{{{{\bf{r}}}}}}}}}),$$
(22)
$${V}_{\tilde{{{{{{{{\bf{p}}}}}}}}}\eta,\,\tilde{{{{{{{{\bf{k}}}}}}}}}0}= \langle \tilde{{{{{{{{\bf{p}}}}}}}}}(t^{\prime} ),\,\eta|{V}_{12}|\tilde{{{{{{{{\bf{k}}}}}}}}}(t^{\prime} ),\,0\rangle \\= \frac{1}{{(2\pi )}^{3}}\iint {d}^{3}{\tilde{{{{{{{{\bf{r}}}}}}}}}}^{\prime}{d}^{3}\tilde{{{{{{{{\bf{r}}}}}}}}}\exp [-i(\tilde{{{{{{{{\bf{p}}}}}}}}}-\tilde{{{{{{{{\bf{k}}}}}}}}})\cdot \tilde{{{{{{{{\bf{r}}}}}}}}}]\\ \times {V}_{12}(\tilde{{{{{{{{\bf{r}}}}}}}}},\,{\tilde{{{{{{{{\bf{r}}}}}}}}}}^{\prime}){[{\psi }_{\eta }({\tilde{{{{{{{{\bf{r}}}}}}}}}}^{\prime})]}^{*}{\psi }_{20}({\tilde{{{{{{{{\bf{r}}}}}}}}}}^{\prime})$$
(23)

and

$${V}_{\tilde{{{{{{{{\bf{p}}}}}}}}}{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime},\,\tilde{{{{{{{{\bf{p}}}}}}}}}\eta }= \left\langle \tilde{{{{{{{{\bf{p}}}}}}}}}(t),\,\tilde{{{{{{{{\bf{p}}}}}}}}}^{\prime} (t)|{V}_{{{{{{{{\rm{ion}}}}}}}}}|\tilde{{{{{{{{\bf{p}}}}}}}}}(t),\,\eta \right\rangle \\= \frac{1}{{(2\pi )}^{3/2}}\int {d}^{3}{\tilde{{{{{{{{\bf{r}}}}}}}}}}^{\prime}{V}_{{{{{{{{\rm{ion}}}}}}}}}({\tilde{{{{{{{{\bf{r}}}}}}}}}}^{\prime}){e}^{-i\tilde{{{{{{{{\bf{p}}}}}}}}}^{\prime} (t)\cdot {\tilde{{{{{{{{\bf{r}}}}}}}}}}^{\prime}}{\psi }_{\eta }({\tilde{{{{{{{{\bf{r}}}}}}}}}}^{\prime}),$$
(24)

where \(\tilde{{{{{{{{\bf{p}}}}}}}}}(t)\), \({\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}(t)\) and \(\tilde{{{{{{{{\bf{k}}}}}}}}}(t)\) are defined according to \(\tilde{{{{{{{{\bf{k}}}}}}}}}(t)=\tilde{{{{{{{{\bf{k}}}}}}}}}+\tilde{{{{{{{{\bf{A}}}}}}}}}(t)\) or \(\tilde{{{{{{{{\bf{k}}}}}}}}}(t)=\tilde{{{{{{{{\bf{k}}}}}}}}}\) in the length or velocity gauge, respectively. In this work, for simplicity, we employ the velocity gauge. This formalism describes the RESI process, in which an electron is ionized by the laser field from the ground state \(\left|0\right\rangle\) at time t into \(|\tilde{{{{{{{{\bf{k}}}}}}}}},\,0 \rangle\) with intermediate momentum \(\tilde{{{{{{{{\bf{k}}}}}}}}}\), it recollides at \(t^{\prime}\) and excites a second electron into the state \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,\eta \right\rangle\) with a final momentum p for the initial electron. The second electron is ionized via the laser field at time t into the state \(\left|\tilde{{{{{{{{\bf{p}}}}}}}}},\,{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}\right\rangle\) with a final momentum \({\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}\). The prefactors give the information about all the bound states91 and interactions, for which we employ V the singly charged binding potential for the first electron, Vion the doubly charged binding potential for the second electron and V12 the electron–electron interaction. In this approximation, electron–electron correlation is described by the prefactor \({V}_{\tilde{{{{{{{{\bf{p}}}}}}}}}\eta,\,\tilde{{{{{{{{\bf{k}}}}}}}}}0}\). The transition amplitude (19) is computed using the steepest descent method. In which, we look for values of the variables t, \(t^{\prime}\)t and \(\tilde{{{{{{{{\bf{k}}}}}}}}}\) such that the action is stationary. This leads to the following saddle-point equations

$${\left[\tilde{{{{{{{{\bf{k}}}}}}}}}+\tilde{{{{{{{{\bf{A}}}}}}}}}(t^{\prime\prime} )\right]}^{2}=-2{I}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{10}}}}}}}}},$$
(25)
$$\tilde{{{{{{{{\bf{k}}}}}}}}}=-\frac{1}{t^{\prime} -t^{\prime\prime} }\int\nolimits_{t^{\prime\prime} }^{t^{\prime} }d\tau \tilde{{{{{{{{\bf{A}}}}}}}}}(\tau ),$$
(26)
$${\left[\tilde{{{{{{{{\bf{p}}}}}}}}}+\tilde{{{{{{{{\bf{A}}}}}}}}}({t}^{\prime})\right]}^{2}={\left[\tilde{{{{{{{{\bf{k}}}}}}}}}+\tilde{{{{{{{{\bf{A}}}}}}}}}(t^{\prime} )\right]}^{2}-2({I}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{20}}}}}}}}}-{I}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{\eta }}}}}}}}}),$$
(27)

and

$${\left[{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime}+\tilde{{{{{{{{\bf{A}}}}}}}}}(t)\right]}^{2}={{{{{{{\boldsymbol{-}}}}}}}}2{I}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{\eta }}}}}}}}}.$$
(28)

Eqs. (25) and (28) give the energy conservation of the first and second electron at the instant of tunnel ionization, (26) enforced the first electron will return and (27) describe energy sharing between the electrons.

Conservation laws

The conservation laws exploited in the main text are general and can be arrived at with few assumptions. The Hamiltonian of the two-electron system may be written as

$$\hat{H}(t)=\frac{{\hat{p}}^{2}+{\hat{p}}^{{\prime} 2}}{2}+({r}_{|\vert }+{r}_{|\vert }^{{\prime} }){E}_{|\vert }(t)+{V}_{{{{{{\rm{ion}}}}}}}(\tilde{{{{{{\bf{r}}}}}}})+{V}_{{{{{{\rm{ion}}}}}}}(\tilde{{{{{{\bf{r}}}}}}}^{\prime} )+{\hat{V}}_{12}.$$
(29)

All terms are independent of the coordinates ϕ and \(\phi ^{\prime}\) except for \({\hat{V}}_{12}\), which depends on the relative distance between the electrons \(|\tilde{{{{{{{{\bf{r}}}}}}}}}-{\tilde{{{{{{{{\bf{r}}}}}}}}}}^{\prime} \vert=\sqrt{{({r}_{|\vert }-r^{\prime}_{|\vert } )}^{2}+{r}_{\perp }^{2}{r^{\prime} }_{\perp }^{2}\cos (\phi -\phi ^{\prime} )}\), thus this is still invariant to a rotation to both particles. Hence, \([{\hat{L}}_{|\vert },\,\hat{H}]=0\), where \({\hat{L}}_{|\vert }=-i{\partial }_{{{{{{{{\rm{\phi }}}}}}}}}-i{\partial }_{{{{{{{{\rm{\phi }}}}}}}}^{\prime} }\), and total OAM is conserved, \(m+m^{\prime}={l}_{{{{{{{{\rm{e}}}}}}}}}+{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}\). If, during ionization by the laser field, the single-active-electron approximation is employed (as is widespread) orbital angular momentum will be conserved during ionization, given that \([{\hat{L}}_{|\vert }^{(1)},\,{\hat{H}}^{(1)}]=[{\hat{L}}_{|\vert }^{(2)},\,{\hat{H}}^{(2)}]=0\), where \({\hat{L}}^{(1)}\), \({\hat{L}}^{(2)}\), \({\hat{H}}^{(1)}\) and \({\hat{H}}^{(2)}\) are the OAM operators and Hamiltonians for each individual electron, including only one-particle terms. Thus, this lead to \({l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}={m}_{\eta }\).

The same conservation laws plus additional constraints are encoded in the SFA via the prefactors, which may be written in terms of their dependence on the azimuthal angles,

$${V}_{\tilde{{{{{{{{\bf{k}}}}}}}}}0,\,0}={e}^{im{\phi }_{{{{{{{{\bf{k}}}}}}}}}}{\tilde{V}}_{{{{{{{{\bf{k}}}}}}}}0,\,0}$$
(30)
$${V}_{\tilde{{{{{{{{\bf{p}}}}}}}}}\eta,\,\tilde{{{{{{{{\bf{k}}}}}}}}}0}={e}^{i(m^{\prime} -{m}_{\eta }){\phi }_{{{{{{{{\bf{p}}}}}}}}}}{\tilde{V}}_{{{{{{{{\bf{p}}}}}}}}\eta,\,{{{{{{{\bf{p}}}}}}}}^{\prime\prime} 0}$$
(31)
$${V}_{\tilde{{{{{{{{\bf{p}}}}}}}}}{\tilde{{{{{{{{\bf{p}}}}}}}}}}^{\prime},\,\tilde{{{{{{{{\bf{p}}}}}}}}}\eta }={e}^{i{m}_{\eta }{\phi }_{{{{{{{{\bf{p}}}}}}}}^{\prime} }}{\tilde{V}}_{{{{{{{{\bf{p}}}}}}}}{{{{{{{\bf{p}}}}}}}}^{\prime},\,\tilde{{{{{{{{\bf{p}}}}}}}}}\eta },$$
(32)

where the tilde indicates quantities independent of the azimuthal angle of all coordinates. Now using (5) and the above equations we may write the SFA OAM transition amplitude

$${M}_{{l}_{{{{{{{{\rm{e}}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}^{{{{{{{{\rm{RESI}}}}}}}}}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )={i}^{-({l}_{{{{{{{{\rm{e}}}}}}}}}+{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime})}{\delta }_{m,\,0}{\delta }_{m^{\prime} -{m}_{\eta },\,{l}_{{{{{{{{\rm{e}}}}}}}}}}{\delta }_{{m}_{\eta },\,{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}\tilde{M}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )$$
(33)

with

$$\tilde{M}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )\,=\,\int \,{d}^{3}t\,\int \,{d}^{2}{{{{{{{\bf{k}}}}}}}}{\tilde{V}}_{{{{{{{{\bf{p}}}}}}}}{{{{{{{\bf{p}}}}}}}}^{\prime},\,{{{{{{{\bf{p}}}}}}}}\eta }{\tilde{V}}_{{{{{{{{\bf{p}}}}}}}}\eta,\,{{{{{{{\bf{k}}}}}}}}0}{\tilde{V}}_{{{{{{{{\bf{k}}}}}}}}0,\,0}{e}^{iS({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime},\,{{{{{{{\bf{k}}}}}}}},\,t,\,t^{\prime},\,t^{\prime\prime} )}.$$
(34)

Thus, we recover the above stated conservation equations along with the condition m = 0. The second electron in its ground states has quantum magnetic number \(m^{\prime}\), and from the behavior of the spherical harmonic \({Y}_{{\ell }_{\eta }}^{-m^{\prime} }({{\Omega }})={(-1)}^{m^{\prime} }\overline{{Y}_{{\ell }_{\eta }}^{m^{\prime} }}({{\Omega }})\) we can deduce that, for odd values of \(m^{\prime}\) we will get opposite signs in the final transition amplitude and thus odd pairs will cancel, as the initial states will be degenerate. Thus, for the s and p initial states employed here, we can assume \(m^{\prime}=0\).

Density matrix

The full density matrix \(\rho=\left|\psi \right\rangle \left\langle \psi \right|\), where \(\left|\psi \right\rangle\) is given by (7), is

$$\rho=\mathop{\sum}\limits_{{l}_{{{{{{\rm{e}}}}}}},\,\,{l}_{{{{{\rm{e}}}}}}^{{\prime} }}\iint {d}^{2}{{{{{\bf{p}}}}}}{d}^{2}{{{{{\bf{p}}}}}}^{\prime} {d}^{2}{{{{{\bf{p}}}}}}^{\prime\prime} {d}^{2}{{{{{\bf{p}}}}}}^{\prime \prime \prime} {M}_{{l}_{{{{{\rm{e}}}}}},\,-{l}_{{{{{\rm{e}}}}}}}({{{{{\bf{p}}}}}},\,\,{{{{{\bf{p}}}}}}^{\prime} ){M}_{{l}_{{{{{\rm{e}}}}}}^{{\prime} }-{l}_{{{{{\rm{e}}}}}}^{{\prime} }}^{\ast }({{{{{\bf{p}}}}}}^{\prime\prime},\,{{{{{\bf{p}}}}}}^{\prime \prime \prime} )|{{{{{\bf{p}}}}}},\,\,{l}_{{{{{\rm{e}}}}}},\,{{{{{\bf{p}}}}}}^{\prime},\,-{l}_{{{{{\rm{e}}}}}}\rangle \langle {{{{{\bf{p}}}}}}^{\prime \prime},\,{l}_{{{{{\rm{e}}}}}}^{{\prime} },\,\,{{{{{\bf{p}}}}}}^{\prime \prime \prime},\,-{l}_{{{{{\rm{e}}}}}}^{{\prime} }|.$$
(35)

We do not compute this explicitly due to the continuous momentum coordinates, instead most commonly, we will trace out the momentum coordinates

$${\rho }_{{{{{{{{\rm{OAM}}}}}}}}}=\iint {d}^{2}{{{{{{{\bf{k}}}}}}}}{d}^{2}{{{{{{{\bf{k}}}}}}}}^{\prime} \left\langle {{{{{{{\bf{k}}}}}}}},\,{{{{{{{\bf{k}}}}}}}}^{\prime} \right|\rho \left|{{{{{{{\bf{k}}}}}}}},\,{{{{{{{\bf{k}}}}}}}}^{\prime} \right\rangle,$$
(36)

where we assume \(\left\langle {{{{{{{\bf{k}}}}}}}}|{{{{{{{\bf{p}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}\right\rangle=\delta ({{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{p}}}}}}}})\left|{l}_{{{{{{{{\rm{e}}}}}}}}}\right\rangle\). Here \(\left|{l}_{{{{{{{{\rm{e}}}}}}}}}\right\rangle\) is an OAM state with the property \(\left\langle r|{l}_{{{{{{{{\rm{e}}}}}}}}}\right\rangle \propto {e}^{i{l}_{{{{{{{{\rm{e}}}}}}}}}\phi }\). Applying this rule results in (8). In order to compute the momentum dependent logarithmic negativity \({{{{{{{{\rm{E}}}}}}}}}_{{{{{{{{\mathcal{N}}}}}}}}}({p}_{|\vert },\,{p}_{|\vert }^{\prime})\) we need to compute the density matrix conditioned on some specific final momentum

$$\rho ({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )=\mathop{\sum}\limits_{{l}_{{{{{{{{\rm{e}}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}\,{M}_{{l}_{{{{{{{{\rm{e}}}}}}}}},\,-{l}_{{{{{{{{\rm{e}}}}}}}}}}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} ){M}_{{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime},\,-{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}^{*}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} )\left|{{{{{{{\bf{p}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime},\,-{l}_{{{{{{{{\rm{e}}}}}}}}}\right\rangle \left\langle {{{{{{{\bf{p}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime},\,{{{{{{{\bf{p}}}}}}}}^{\prime},\,-{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}\right|.$$
(37)

Positive partial transpose

The positive partial transpose (PPT) or Peres-Horodecki criterion73,74 is a necessary condition on density matrices to determine if a system is separable. It is valid for both pure and mixed states. The approach is to take the partial transpose—i.e., transpose one subsystem—and compute the eigenvalues, if any are negative the state is non-separable and thus, entangled. For 2  2 and 2  3 the condition is also sufficient, so no negative eigenvalues imply separability, for higher dimensional systems (such as our 3  3 system), this is not the case, however, it is still a very powerful method and witnesses can be constructed to detect any PPT entangled state92, which we can exploit to our advantage.

For our NSDI qutrit mixed state

$${\rho }_{{{{{{\rm{OAM}}}}}}}=\mathop{\sum }\limits_{{l}_{{{{{{\rm{e}}}}}}},\,{l}_{{{{{{\rm{e}}}}}}}^{\prime}=-1}^{1}{\alpha }_{{l}_{{{{{{\rm{e}}}}}}}{l}_{{{{{{\rm{e}}}}}}}^{\prime}}\left|{l}_{{{{{{\rm{e}}}}}}},\,-{l}_{{{{{{\rm{e}}}}}}}\right\rangle \left\langle {l}_{{{{{{\rm{e}}}}}}}^{\prime},\,-{l}_{{{{{{\rm{e}}}}}}}^{\prime}\right|,$$
(38)

where α is defined as in (8). In taking the partial transpose we swap the indices for one of the subsystems (in this case the second electron),

$${\rho }_{{{{{{{{\rm{OAM}}}}}}}}}=\mathop{\sum }\limits_{{l}_{{{{{{{{\rm{e}}}}}}}}},\,{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}=-1}^{1}{\alpha }_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}\vert {l}_{{{{{{{{\rm{e}}}}}}}}},\,-{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}\rangle \langle {l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime},\,-{l}_{{{{{{{{\rm{e}}}}}}}}}|.$$
(39)

The eigenvalues can be computed analytically as

$${\alpha }_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}}=\int {d}^{2}{{{{{{{\bf{p}}}}}}}}\int {d}^{2}{{{{{{{\bf{p}}}}}}}}^{\prime}|{M}_{{l}_{{{{{{{{\rm{e}}}}}}}}}-{l}_{{{{{{{{\rm{e}}}}}}}}}}({{{{{{{\bf{p}}}}}}}},\,{{{{{{{\bf{p}}}}}}}}^{\prime} ){|}^{2}\,\,{{{{{{{\rm{for}}}}}}}}\,\,{l}_{{{{{{{{\rm{e}}}}}}}}}\in [-1,\,1],$$
(40)

which are positive and \(\pm|{\alpha }_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}|\) for \({l}_{{{{{{{{\rm{e}}}}}}}}}\,\ne\, {l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}\), which provides three negative eigenvalues as long as \(|{\alpha }_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}|\,\ne\, 0\) for \({l}_{{{{{{{{\rm{e}}}}}}}}}\,\ne\, {l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}\). Thus, the electrons from the RESI mechanism of NSDI are always PPT entangled as long as there is non-zero population across different excited state with differing mη. This also means there will always be an entanglement witness, which may be used to experimentally verify this entanglement.

Logarithmic negativity

The logarithmic negativity is well-suited for quantifying entanglement in PPT entangled states, as it is equal to the entanglement cost to create this entanglement via PPT operations. In the main article we construct the logarithmic negativity from the reduced density matrix, traced over momentum coordinates, however, we may instead define a momentum dependent logarithmic negativity

$${E}_{{{{{{{{\mathcal{N}}}}}}}}}({p}_{|\vert },\,p^{\prime} )={\log }_{2}\left[{\left|\left|{\rho }^{{T}_{{{{{{{{\rm{A}}}}}}}}}}({p}_{|\vert },\,\delta p,\,p^{\prime},\,\delta p)\right|\right|}_{1}\right],$$
(41)

where the perpendicular momentum coordinates are set to nearly to zero, δp = 0.05 a.u., in order allow for 2D visualization. The value is chosen to be where the distribution will have high probability but no nodes due to the geometry of the excited states.

The logarithmic negativity can be related to the sum of the negative eigenvalues (λi) of the density matrix

$${E}_{{{{{{{{\mathcal{N}}}}}}}}}= {\log }_{2}\left[|\vert {\rho }_{{{{{{{{\rm{OAM}}}}}}}}}^{{T}_{{{{{{{{\rm{A}}}}}}}}}}|{|}_{1}\right]\\= {\log }_{2}\left(1+2\left|\mathop{\sum}\limits_{{\lambda }_{i} < 0}{\lambda }_{i}\right|\right)$$
(42)

which can be written as minus the sum of the negativive eigenvalues

$$={\log }_{2}\left(1+2\left(|{\alpha }_{-10} \vert+\vert {\alpha }_{-11} \vert+\vert {\alpha }_{01}|\right)\right)$$
(43)

Using the Cauchy–Schwartz inequality we can show \({\alpha }_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}\le \sqrt{{\alpha }_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}}{\alpha }_{{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}}\), leading to

$$ |{\alpha }_{-10} \vert+\vert {\alpha }_{-11} \vert+\vert {\alpha }_{01}|\\ \le \sqrt{{\alpha }_{-1-1}{\alpha }_{00}}+\sqrt{{\alpha }_{-1-1}{\alpha }_{11}}+\sqrt{{\alpha }_{00}{\alpha }_{11}}\\ \le \sqrt{{\alpha }_{-1-1}+{\alpha }_{11}+{\alpha }_{00}}\sqrt{{\alpha }_{00}+{\alpha }_{-1-1}+{\alpha }_{11}}\\ =1,$$
(44)

where other forms of the Cauchy–Schwartz inequality were used in the additional inequalities. With this we can place a bound upon the logarithmic negativity

$${E}_{{{{{{{{\mathcal{N}}}}}}}}}\le {\log }_{2}(1+2)\,\approx\, 1.58\,.$$
(45)

Entanglement witnesses

Here, we will show that the witness used in this work (see (12)) is a valid entanglement witness. We do this by showing that (i) it is positive for all separable states and (ii) the trace is negative for at least one entangled state. For a general separable pure state

$$\left|\psi \right\rangle=\mathop{\sum}\limits_{i}{a}_{i}\left|i\right\rangle \otimes \mathop{\sum}\limits_{j}{b}_{i}\left|j\right\rangle$$
(46)

the density matrix may be written as

$${\rho }_{{{{{{{{\rm{s}}}}}}}}}=\mathop{\sum}\limits_{i,\,j,\,m,\,n}{a}_{i}{a}_{m}^{*}{b}_{j}{b}_{n}^{*}\left|i,\,j\right\rangle \left\langle m,\,n\right|$$
(47)

and we can compute the trace with the witness of (12)

$${{{{{{{\rm{Tr}}}}}}}}[{\rho }_{{{{{{{{\rm{s}}}}}}}}}{{{{{{{\mathcal{W}}}}}}}}(\theta )]=\frac{1}{d}\mathop{\sum}\limits_{i,\,j}|{a}_{i}{|}^{2}|{b}_{j}{|}^{2}-\left\langle \nu (\theta )\right|{\rho }_{{{{{{{{\rm{s}}}}}}}}}\left|\nu (\theta )\right\rangle$$
(48)
$$\quad\quad=\frac{1}{d}-{\left|\mathop{\sum}\limits_{ij}{a}_{i}{b}_{j}\left\langle i,\,j|\nu (\theta )\right\rangle \right|}^{2}$$
(49)
$$\quad=\frac{1}{d}-\frac{1}{d}{\left|\mathop{\sum}\limits_{{l}_{{{{{{{{\rm{e}}}}}}}}}}{a}_{{l}_{{{{{{{{\rm{e}}}}}}}}}}{e}^{i\theta l}{b}_{-{l}_{{{{{{{{\rm{e}}}}}}}}}}\right|}^{2}$$
(50)
$$\quad\quad\quad\quad\ge \frac{1}{d}-\frac{1}{d}\mathop{\sum}\limits_{{l}_{{{{{{{{\rm{e}}}}}}}}}}|{a}_{{l}_{{{{{{{{\rm{e}}}}}}}}}}{e}^{i\theta {l}_{{{{{{{{\rm{e}}}}}}}}}}{|}^{2}\mathop{\sum}\limits_{{l}_{{{{{{{{\rm{e}}}}}}}}}}|{b}_{{l}_{{{{{{{{\rm{e}}}}}}}}}}{|}^{2}=0.$$
(51)

Thus, we have demonstrated (i) and the entangled state we use for (ii) is \(\left|\nu (\theta )\right\rangle\)

$${{{{{{{\rm{Tr}}}}}}}}[\left|\nu (\theta )\right\rangle \left\langle \nu (\theta )\right|{{{{{{{\mathcal{W}}}}}}}}(\theta )]=\frac{1}{d}-|\left\langle \nu (\theta )|\nu (\theta )\right\rangle {|}^{2}$$
(52)
$$\hskip65pt=-\frac{d-1}{d} \, < \, 0.$$
(53)

Hence, given d = 3, we have shown that we are using a valid witness.

Witness decomposition

The decomposition of entanglement witness into a series of local measurement is described in refs. 62,63,79. As described in the main text, for a qutrit, a convenient decomposition is in terms of the Gell-Mann matrices. These are defined by the following construction80,

$$\begin{array}{lll}{\chi }_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}^{\pm }=\vert {x}_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}^{\pm }\rangle \langle {x}_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}^{\pm }|&{{{{{{{\rm{and}}}}}}}}&|{x}_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}^{\pm }\rangle=\frac{1}{\sqrt{2}}(\left|{l}_{{{{{{{{\rm{e}}}}}}}}}\right\rangle \pm \left|{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}\right\rangle ),\\ {{{\Upsilon }}}_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}^{\pm }=\left|{y}_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}^{\pm }\rangle \langle {y}_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}^{\pm }\right|&{{{{{{{\rm{and}}}}}}}}&|{y}_{{l}_{{{{{{{{\rm{e}}}}}}}}}{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}}^{\pm }\rangle=\frac{1}{\sqrt{2}}(\left|{l}_{{{{{{{{\rm{e}}}}}}}}}\right\rangle \pm i\left|{l}_{{{{{{{{\rm{e}}}}}}}}}^{\prime}\right\rangle ),\end{array}$$
(54)

the matrices are then given by

$$\begin{array}{ll}{\lambda }_{1}={\chi }_{01}^{+}-{\chi }_{01}^{-}&{\lambda }_{2}={{{\Upsilon }}}_{01}^{+}-{{{\Upsilon }}}_{01}^{-}\\ {\lambda }_{4}={\chi }_{0-1}^{+}-{\chi }_{0-1}^{-}&{\lambda }_{5}={{{\Upsilon }}}_{0-1}^{+}-{{{\Upsilon }}}_{0-1}^{-}\\ {\lambda }_{6}={\chi }_{1-1}^{+}-{\chi }_{1-1}^{-}&{\lambda }_{7}={{{\Upsilon }}}_{1-1}^{+}-{{{\Upsilon }}}_{1-1}^{-}\\ {\lambda }_{3}=\left|0\right\rangle \left\langle 0\right|-\left|1\right\rangle \left\langle 1\right|&\\ {\lambda }_{8}=\frac{1}{\sqrt{3}}(\left|0\right\rangle \left\langle 0\right \vert+\left|1\right\rangle \left\langle 1\right|\! - \!2\left|-1\right\rangle \left\langle -1\right|).&\end{array}$$
(55)

These matrices can be related to the well-known Pauli matrices. For example, the measurement \({\lambda }_{3}={\hat{L}}_{|\vert }\), corresponds directly to an OAM measurement, similar to the σz Pauli matrix, while λ1, λ4, λ6 can be related to σx as they are all pairwise superposition of two OAM states with the phase±. The measurements λ2, λ5, λ7 relate to σy as they are pairwise superposition with the phase ± i. Measurements λ1λ7 have eigenvalues of −1, 0, 1, as with OAM, while λ8 has eigenvalues \(-2/\sqrt{3}\) and \(1/\sqrt{3}\).

The full entanglement witness decomposition given in terms of the Gell-Mann matrices is

$${{{{{{{\mathcal{W}}}}}}}}(\theta )= \frac{1}{12}\bigg[\frac{8}{3}{{\mathbb{1}}}^{\otimes 2}-{\lambda }_{3}^{\otimes 2}-{\kappa }^{2}({\kappa }^{-4}\,+\,1)({\lambda }_{4}^{\otimes 2}+{\lambda }_{5}^{\otimes 2})+{\lambda }_{8}^{\otimes 2}\\ +\sqrt{3}({\lambda }_{3}\,\otimes \,{\lambda }_{8}\,+\,{\lambda }_{8}\,\otimes \,{\lambda }_{3})+i{\kappa }^{2}({\kappa }^{-4}\,\,-\,1)({\lambda }_{4}\,\otimes \,{\lambda }_{5}\,-\,{\lambda }_{5}\,\otimes \,{\lambda }_{4})\\ -\kappa ({\kappa }^{-2}\,+\,1)({\lambda }_{1}\,\otimes \,{\lambda }_{6}\,+\,{\lambda }_{6}\,\otimes \,{\lambda }_{1}\,+\,{\lambda }_{2}\,\otimes \,{\lambda }_{7}\,+\,{\lambda }_{7}\,\otimes \,{\lambda }_{2})\\ -i\kappa ({\kappa }^{-2}\,-\,1)({\lambda }_{1}\,\otimes \,{\lambda }_{7}\,-\,{\lambda }_{7}\,\otimes \,{\lambda }_{1}\,-\,{\lambda }_{2}\,\otimes \,{\lambda }_{6}\,+\,{\lambda }_{6}\,\otimes \,{\lambda }_{2})\bigg],$$
(56)

where κ = eiθ. This decomposition is specific to the Gell-Mann matrices, however, a general procedure may be outlined for any complete set of measurement bases. Thus, if it is simpler to measure in another way, use another entanglement witness or perform this for a higher-dimensional system, the following recipe can be used to achieve the decomposition. For a complete set of one-particle measurements defined by the operators λi for i [1, d2] and a known entanglement witness \({{{{{{{\mathcal{W}}}}}}}}\), the decomposition in terms of the local one-particle observables may be written as

$${{{{{{{\mathcal{W}}}}}}}}=\mathop{\sum}\limits_{ij}{c}_{ij}{\lambda }_{i}\otimes {\lambda }_{j},$$
(57)

where the coefficients cij can be used to combined with experimental local expectation values to determine \({{{{{{{\mathcal{W}}}}}}}}\). This sum can be inverted to find cij via vectorization, i.e., flattening one dimension. We may define

$${{{{{{{{\mathcal{W}}}}}}}}}_{{{{{{{{\rm{v}}}}}}}}} ={{{{{{{\rm{vec}}}}}}}}({{{{{{{\mathcal{W}}}}}}}})\\ ={({{{{{{{{\mathcal{W}}}}}}}}}_{0,\,0},\,{{{{{{{{\mathcal{W}}}}}}}}}_{0,\,1},\,...,\,{{{{{{{{\mathcal{W}}}}}}}}}_{0,\,{d}^{2}-1},\,{{{{{{{{\mathcal{W}}}}}}}}}_{1,\,0},\,...,\,{{{{{{{{\mathcal{W}}}}}}}}}_{{d}^{2}-1,\,{d}^{2}-1})}^{T}$$
(58)

and \({c}_{{{{{{{{\rm{v}}}}}}}}}={({c}_{0,\,0},\,{c}_{0,\,1},\,...,\,{c}_{0,\,{d}^{2}-1},\,{c}_{1,\,0},\,...,\,{c}_{{d}^{2}-1,\,{d}^{2}-1})}^{T}\). This reduces the d2 × d2 matrix operators, \({{{{{{{\mathcal{W}}}}}}}}\) and c, to d4 dimensional vectors. Now we defined a d4 × d4 matrix

$${M}_{{{{{{{{\rm{v}}}}}}}}}=\left({{{{{{{\rm{vec}}}}}}}}({\lambda }_{0}\otimes {\lambda }_{0}),\,..,\,{{{{{{{\rm{vec}}}}}}}}({\lambda }_{d-1}\otimes {\lambda }_{d-1})\right),$$
(59)

where each row is the vectorized matrix λiλj for specific values of i and j. Now using these definitions we may rewrite (57) as a linear matrix equation \({{{{{{{{\mathcal{W}}}}}}}}}_{{{{{{{{\rm{v}}}}}}}}}=M{c}_{{{{{{{{\rm{v}}}}}}}}}\), thus we may obtain the coefficients by inverting, such that \({c}_{{{{{{{{\rm{v}}}}}}}}}={M}^{-1}{{{{{{{{\mathcal{W}}}}}}}}}_{{{{{{{{\rm{v}}}}}}}}}\). For our finite-dimensional system with d = 3, this can be done quickly and straightforwardly.