Manipulating twisted electrons in strong-field ionization

Abstract

We investigate the discrete orbital angular momentum (OAM) of photoelectrons freed in strong-field ionization. We use these ‘twisted’ electrons to provide an alternative interpretation on existing experimental work of vortex interferences caused by strong field ionization mediated by two counter-rotating circularly polarized pulses separated by a delay. Using the strong field approximation, we derive an interference condition for the vortices. In computations for a neon target we find very good agreement of the vortex condition with photoelectron momentum distributions computed with the strong field approximation, as well as with the time-dependent methods Qprop and R-Matrix. For each of these approaches we examine the OAM of the photoelectrons, finding a small number of vortex states localized in separate energy regions. We demonstrate that the vortices arise from the interference of pairs of twisted electron states. The OAM of each twisted electron state can be directly related to the number of arms of the spiral in that region. We gain further understanding by recreating the vortices with pairs of twisted electrons and use this to determine a semiclassical relation for the OAM. A discussion is included on measuring the OAM in strong field ionization directly or by employing specific laser pulse schemes as well as utilizing the OAM in time-resolved imaging of photo-induced dynamics.

---

I. Introduction

Phase vortices in light and matter are the product of the orbital angular momentum (OAM) carried by free particles resulting in a rotating wavefront. The OAM of a free particle is a quantum number that can be measured and manipulated, with many of the same properties as spin, yet it has only been observed in experiments fairly recently in photons (3 decades ago)¹ and even later for electrons (1 decade ago)² (see ref. 3 and 4 for reviews). It has been found that these twisted photons or electrons have huge potential in imaging, and recent experiments with electron vortex beams have revealed their role in chiral energy-loss spectroscopy and magnetic dichroism.³

The study of twisted light and electrons in strong field processes is a very new but rapidly developing topic. In high-harmonic generation (HHG) it has been shown that a twisted driving field can be used to produce twisted UV light.^5–9 Later this was made time-varying¹⁰ and further extended to the generation of UV light with torus knot topology.¹¹ The effect of twisted electrons has been investigated in HHG.¹² However, the effect of OAM on photoelectrons in strong-field ionization has not been studied in detail.

Recently, the OAM content of photoelectrons was computed using the strong field approximation (SFA). This work focused on reaching exceedingly high OAM values for quasi-relativistic field intensities¹³ and terahertz fields.¹⁴ Twisted electrons have also been studied in the context of rescattering with a well-defined OAM.¹⁵ However, currently there is very little work focused on the generation, measurement and manipulation of photoelectrons carrying OAM for strong field ionization. That said, there is a still a lot to learn from related strong field work. Through conservation of momentum, photoelectrons ionized by a linearly polarized laser field will adopt the same OAM as the magnetic quantum number of the bound state of origin.¹⁵ The linear field does not alter the magnetic quantum number or the photoelectron’s OAM, and there is a direct mapping between the two. Furthermore, it has been demonstrated that circularly polarized laser fields preferentially tunnel ionized p electrons, whose magnetic quantum number m dictates that they are ‘counter-rotating’ with respect to the laser field.^16,17 Depending on the helicity of the laser field, either m = 1 or m = −1 electrons dominate the signal, while ionization from orbitals with m = 0 is strongly suppressed. Thus, given the interaction with the magnetic quantum number, tunnel ionization via circularly polarized light may be a route to generate, control, and measure photoelectrons with non-trivial OAM.

A more explicit example of twisted electrons in multiphoton processes is the formation of interference vortices.^18–27 Computations of the momentum distribution of photoelectrons ionized via two time-delayed, counter-rotating circularly-polarized laser fields reveal a Fermat spiral interference pattern. Initially treated theoretically for single-photon ionization¹⁸ and later two-photon ionization¹⁹ of helium, such spirals have since been observed in both single and double ionization of molecules.^20,21 Experimental verification came in ref. 22, where the spiral interference pattern was demonstrated for three-photon ionization of potassium. Calculations using the multielectron ab initio R-Matrix with time dependence (RMT) method fully corroborated this effect,²⁷ and revealed the three-dimensional characteristics of the interference vortices. Experimental progress in this domain continues to gather pace, and vortex measurements at mid-infrared wavelengths are being carried out.²⁸ Interference vortices are indicative of photoelectrons carrying OAM, as the interference of two states with differing OAM produces such spirals. The above-mentioned multiphoton studies have provided explanations in terms of the magnetic quantum number in the bound states and scattering states to which the photoelectron is promoted, but do not go as far as discussing the OAM of the free electron or its measurement and control. In this work we examine the interference vortices in the non-perturbative regime, where a description of promotion to scattering states will not hold. However, analysis of co- and counter-rotating electrons^16,17 is valid and provides additional insight.

Twisted electrons and the OAM distributions are the focus of this article, and with these tools we present an alternative understanding of the vortices. Using the SFA,²⁹ as well as the time-dependent Schrödinger equation (TDSE)-based solvers, Qprop^30,31 and RMT,^32–34 we compute both the photoelectron’s momentum distribution and its constituent OAM components. The resulting momentum–space vortices may be understood by deriving an interference condition using the SFA, and by recreating the same pattern using the interference between states of different values of OAM. We also identify a semi-classical relationship that demonstrates a link between the OAM and above threshold ionization peaks of the photoelectron. Understanding the OAM in photoelectrons, which undergo strong-field ionization via circularly-polarized light, imbues the seemingly plain momentum distributions with previously unseen structure revealed only by such interferometric schemes. This opens the question: is there a generalised way to measure, directly or indirectly, the OAM of photoelectrons in strong-field experiments? Given the high correlation between the photoelectron OAM and the quantum numbers of the initial, bound state, consideration of the OAM could be extremely useful in relating physical observables to the bound state and disentangling interferences in time-resolved measurements. Furthermore, the potential for the outgoing twisted photoelectron states to be chiral could aid the detection and spectrographic measurement of chiral molecular targets.³⁵

The article is structured as follows. In Section II we outline the theoretical background for the three methods employed in this article, the SFA (Section II A), the single active electron (SAE) 3D-TDSE solver Qprop (Section II B), and the multielectron TDSE solver RMT (Section II C). In Section II D we present the methodology for computing the OAM distributions for each model. Next we focus on the interference vortices themselves. In Section III, we present the derivation of the vortex interference condition using the SFA (Section III A) and compare the results of all methods and that condition (Section III B). While in Section IV we turn our attention to the orbital angular momentum (OAM); in Section IV A we compare the photoelectron OAM distribution from all methods and in Section IV B we use this to construct the interference vortex entirely from twisted electron states. In Section V we make our concluding statements and discuss potential direct and indirect measurement schemes. We employ atomic units throughout (denoted a.u.), where the elementary charge, electron mass and the Planck constant are set to one, e = m = ℏ = 1.

II. Background and methods

1. Laser field

We consider a bicircular field comprised of two time-delayed, 400 nm pulses, described by the vector potential

A(t) = A₊(t) + A₋(t − δ), (1)

where ± denotes the rotation of the circular field. Each field is given bywhere is the peak vector potential strength, which is related to the peak electric field amplitude by A₀ = E₀/ω. The laser frequency is ω = 0.114 a.u. (corresponding to a 400 nm wavelength), N_c is the total number of laser cycles, and δ is the time delay between the pulses, which we will typically take to be so that the two field envelopes do not overlap. Each pulse reaches an associated laser peak intensity of 2.5 × 10¹³ W cm⁻² and the field retains this profile for time t ∈ [0, 2πN_c/ω], and is zero otherwise.

A. Strong field approximation. Our starting point in the SFA²⁹ is the transition amplitude for direct ATI from the initial bound state |ψ₀〉 to a final Volkov state with drift momentum p given by^36–38where d(p,t′) = 〈p + A(t′)|V|Ψ₀〉, |Ψ₀〉 is the ground state of the target and the action is given byHere, I_p is the ionization potential of our target. We employ the saddle point approximation, seeking stationary action for the integration variable t′,Now the probability distribution can be computed from eqn (3) aswhere the prefactor c(p,t_s,t), derived from application of the saddle point approximation and also includes the t′ independent phase from eqn (3), is given by

In order to capture the essential physics of the interference vortices, we make some additional assumptions. The minimal requirement is to have interference between two photoelectron orbits deriving from each pulse. Thus, in the SFA model we take only two ionization events from near the peak of each laser pulse, A₊(t) and A₋(t − δ), but otherwise neglect the pulse envelopes to allow for an analytic description. Thus, in the SFA model, the vector potentials are defined by

where the ± denotes the direction of rotation of the field. In this way we are able to separate the two contributions of ionization at the times t₊ and t₋ + δ, given byIn eqn (9) we have written the solutions in spherical momentum coordinates p, θ and ϕ as they are the most natural for this system. The transition amplitude can then be written as a sum of the two separate solutions

M(p) = c₊(p,t)d₊(p)exp(iS₊(p)) + c₋(p,t)d₋(p)exp(iS₋(p)). (10)

The actions S₊(p), S₋(p) and prefactors c₊(p,t), c₋(p,t), d₊(p) and d₋(p) have the times t₊ or t₋ + δ substituted into them such thatandFollowing the prescription of ref. 16, we split the prefactor d(p,t) into two parts

d(p,t) = e^imϕ_k(p,t)f_nlm(k(p,t),θ_k(p,t)), (13)

where k = p + A(t) and k(p,t), θ_k(p,t) and ϕ_k(p,t) are the spherical coordinates of k. The function f_nlm(k,θ) is dependent on the specific bound state, which in our case is the ground state of neon, while ϕ_k is the so-called tunneling angle,¹⁶ which leads to the enhanced ionization of ‘counter-rotating’ electrons. It will continue to play an important role in the description of the interference vortices.

B. Qprop. The 3D-TDSE calculations were performed using the latest version of Qprop,³¹ where a very accurate method for the calculation of photoelectron spectra (PES), dubbed i-SURFV, is implemented. In short, Qprop is a velocity gauge 3D-TDSE solver that allows studies within both the single active electron (SAE) approximation and many-electron systems via the solution of the time-dependent Kohn–Sham equations. From a computational viewpoint, the velocity gauge formulation of the strong laser–matter problem appears to be more efficient, providing results of accuracy equivalent to the length gauge with much less computational effort. For our SAE model of the Ne atom, we have employed the model potential of ref. 39 and started the time-dependent simulations from the 2p₀, 2p₋₁ and 2p₊₁ bound states. These states were computed via imaginary time propagation. The pulse parameters are those given in Section II 1 and, specifically, the effective potential has the formwhere Z = 1. For neon, the coefficients a_i are a₁ = 8.069, a₂ = 2.148, a₃ = −3.570, a₄ = 1.986, a₅ = 0.931, a₆ = 0.602,³⁹ which gives the correct ionization potential of I_p = 0.79 a.u. In this computation we considered angular momenta up to L = 14 and all magnetic sublevels.

C. R-Matrix with time dependence. The R-Matrix with time dependence method (RMT) is an ab initio approach that solves the time-dependent Schrödinger equation for single ionization of multielectron atoms, ions, and molecules in arbitrarily-polarized, strong laser fields.^32–34 The method divides position space into two distinct regions. The inner region, which usually extends to around 20 Bohr radii (20 a.u.), confines the target nucleus. In this region, a many-body wave function is treated, accounting for electron exchange and correlation. The outer region, which extends to larger distances, contains a single ejected electron that interacts with both the laser field and the singly-ionized residual target.

Our description of the neon target is provided in a previous work.⁴⁰ In short, we couple a single electron to the Ne⁺ residual ion, expanding the wave function in LM_LSπ symmetries up to L = 99, retaining both the 1s²2s²2p⁵ and 1s²2s2p⁶ Ne⁺ residual ion states. The photoelectron wave function is calculated on a finite-difference radial grid in the outer region, which extends to 3400 a.u. The laser pulse configuration is that given in Section II 1. The wave function is propagated in time for a total of 35 fs, which gives sufficient time for the ionizing wavepacket to extend to large distances within the outer region. The photoelectron momentum distribution is obtained by performing a Fourier transform of the ejected-electron radial wavefunction.

D. Computing OAM distributions. In order to compute the orbital angular momentum distributions we will take two approaches. Firstly, given a transition amplitude such as that in the SFA we can transform it into the ‘OAM basis’ by computing its Fourier series coefficients¹³The coordinates p⊥ and p_‖ are cylindrical momentum coordinates, where p_‖ is the axial coordinate parallel to the laser field propagation and p is the radial distance given by . However, given that this is a transformation into cylindrical co-ordinates, and many TDSE solvers employ a basis of spherical harmonics, some computational power can be saved by converting a distribution expanded in such a basis,so that the transformation of eqn (15) reduces towhere tan θ = p_⊥/p_‖, and L is the maximum value of the azimuthal quantum number l required for adequate convergence.

III. Vortex interference

A. Vortex interference condition from the standard SFA implementation

In this section we derive a condition for the interference vortices from the SFA description. Such conditions can help to give some insight into the physical dynamics of the process. The starting point to derive the condition is eqn (10). We can demonstrate that the prefactors c₊(p,t) and c₋(p,t) are the same and given byThus, c_±(p,t) can be factored out and will not play a role in the interference. For the prefactor d_±(p), co-ordinates k and θ_k of f_nlm(k,θ_k) can be shown to be the same when evaluated at either t₊ or t₋ + δ. From the saddle point eqn (5) we can determineThus, we can deduce

f_nlm(k(p,t₊),θ_k(p,t₊)) = f_nlm(k(p,t₋ + δ),θ_k(p,t₋ + δ)) =: f_nlm(k,θ_k). (20)

Hence, f_nlm(k,θ_k) can be factored out and we are left with the e^imϕ_k term, which (as presented in ref. 16) will act to weight the contribution from each laser pulse depending on the value of m. For non-zero m this will lead to a loss of contrast of the interference fringes as electrons ‘co-rotating’ with the field will be preferentially ionized. Now we may write the transition amplitude as

M(p) = Ce^imϕ_k+(p)e^iS₊(p)(1 + e^imΔϕ_k(p)e^iΔS(p)), (21)

where C contains all the prefactors. The term e^imΔϕ_k(p) acts to switch off or blur the interference in the way described. The interference condition can be derived by maximizing the term e^iΔS(p), whereSetting ΔS(p) = 2πn, where n is an integer we can find its maxima and thus a condition for the interference vortices

(2I_p + p² + 2U_p)(2ϕ + δω) = 4πnω. (23)

This can be re-arranged to give the equation of a spiralFrom eqn (24) for any particular angle ϕ there is a minimum value of n, below which there will be no real solutions for the radius p. This is reminiscent of similar interference conditions that may be derived in the SFA e.g. to describe above-threshold ionization (ATI) peaks.⁴¹ The condition given by eqn (24) is capable of describing a varying number of spiral arms. In Fig. 1 we plot the spiral for different values of n using both the positive and negative branches. In the three panels from left to right, we use the first 7, next 8 and the following 9 integer values of n, respectively, and both branches of eqn (24). Using more and higher values of n leads to the spiral arms filling a larger region, with some arms becoming crowded in the center, more evenly spaced for higher p and finally developing gaps for the highest p. Thus, this suggests that in specific regions there will be a particular number of evenly spaced spiral arms (similar to Fermat spirals), which increases with p. In short, we expect the interference vortices to have an increasing number of spiral arms with p.

Fig. 1 Plotting the interference vortex condition over the polarization p_xp_y-plane of the laser field. The three spirals are plotted using 3 different sets of n from eqn (24); the first 7 valid n, the next 8 and then the next 9 are used, from left to right. The minimum valid value of n is 37 for our field, target parameters and ϕ = π. Using a delay δ = 4 × 2π/ω, an angular frequency of ω = 0.114 a.u., U_p = 0.027 a.u., corresponding to a peak laser intensity of I₀ = 2.5 × 10¹³ W cm⁻². The ionization potential is I_p = 0.79 a.u. for a neon target.

B. Computation of interference vortices

In this section we present numerical calculations of strong field ionization of neon via the two counter-rotating fields described above. In Fig. 2 we plot the result of performing this calculation with the SFA [left column], Qprop [central column] and RMT [right column], for initial 2p orbitals of ground state neon with magnetic quantum number m = 0 [top row] and m = 1 [bottom row]. The vortices are clearly visible in all plots. The interference vortex condition is superimposed on all plots and matches very well, clearly capturing the core behaviour. It is notable that such a simplified SFA model is able to capture the features of the vortex. There are however differences, such as the hole in the center of the SFA results, which is less visible in the RMT and Qprop outcomes. This is likely due to the Coulomb force, which will decelerate the outgoing electron wavepacket and thus reduce its final momentum closing the hole. This can be roughly accounted for in the SFA if the radial final momentum is taken to be lower than the initial, which will reduce the size of the hole. One clear aberration is the discontinuity along ϕ = 0 in the SFA results. In the SFA, ϕ is fixed between ±π yielding a discontinuity. In the current construction the interference vortices can only be reproduced properly if ϕ can vary over a range larger than 2π in some regions, as it does in the interference condition.

Fig. 2 Momentum dependent photoelectron yield of 2p₀ [top row] and 2p₁ [bottom row] electrons ionized by two counter-rotating, 400 nm, 2.5 × 10¹³ W cm⁻², 4-cycle sin² pulses separated by a delay of 4 laser cycles, computed using the SFA [left], Qprop [centre] and the RMT method [right]. All figures show the vortex interferences, upon which we have superimposed the spiral condition derived above with dashed lines. The colours of the lines correspond to the value of n used in eqn (23). The Qprop and RMT momentum distributions sample states with l > 3 in the spherical harmonic expansion to resolve the processes involving more than 3 photons. All panels are plotted on a logarithmic scale. A different scale is used for each panel to highlight the suppression of the photoelectron yield from 2p₀ in comparison to 2p₁. The units of the scale are arbitrary but comparisons can be made between computations using the same model. The distributions are calculated in the p_∥ = 0.1 a.u. plane. This was chosen as at p_∥ = 0 there is a node for m = 0 leading to nearly zero signal.

In the RMT and Qprop momentum distribution calculation the spherical harmonic expansion was limited to angular momenta l > 3. The rationale for this is that the spiral interference arises from an ionization process involving at least 7 photons, which typically leads to population of states close to l = 7. Population in states of low angular momenta are likely due to processes involving intermediate excited states, which are not accounted for by the SFA treatment. At the specific field parameters used here (peak intensity of 2.5 × 10¹³ W cm⁻² and 4-cycle sin² pulses), some contributions from such intermediate excited states can be identified for neon. This leads to additional interference vortices [not shown] superimposed on those presented in Fig. 2.

We find that the intermediate state effects are sensitive to the laser pulse parameters, and their contribution may be reduced if alternative parameters are used. We have verified that for higher intensities (10¹⁴ W cm⁻²) or long pulses (10-cycle sin² pulses), the low values of angular momenta do not play a significant role, and the interference vortices arising from at least 7-photon ionization from the ground state are the dominant effect. In fact very similar spirals can be produced with QProp and RMT using all angular momenta for an intensity of 10¹⁴ W cm⁻². This is aided by the fact the spiral fringes do not vary considerably with intensity in this regime. For longer pulses, the interference vortices are visible along with the characteristic ATI rings/peaks, with the number of arms of each spiral increasing with the order of the ATI processes. A similar effect was observed in lower-order ATI processes in recent experiments.²²

Although it is difficult to resolve by eye, the number of spiral arms does vary, in Fig. 2, increasing with momentum. We have verified this by performing Fourier analysis on the momentum distributions across rings of fixed cylindrical radial momentum p_⊥ so that we can count the number of spiral arms intersecting each ring. This is also corroborated by the very good match with the interference condition. Thus, very good agreement can be achieved with a relatively simple SFA model, where many features can be identified. However, if we now examine the OAM of the outgoing photoelectron we can achieve additional insight and understanding.

IV. Interpretation of vortex in terms of OAM

A. Numerical OAM results

Using eqn (15) and (17) we compute the ‘OAM distributions’ for the three models presented. The vortex state basis enables a description of the interference vortices in terms of the OAM of the freed electron. In Fig. 3 we show the distribution of OAM over perpendicular momentum for a fixed p_∥ = 0.1 a.u., chosen as there is a node at p_‖ = 0 leading to a very low ionization signal for the case m = 0. The left column of Fig. 3 displays results for neon initially in the m = 0, p state. All distributions show a single peak corresponding to each of three momentum regions, with a small, fourth peak at the highest momentum shown. For each peak there are actually two overlapping sets of lines corresponding to OAM values with opposite signs (negative OAM is shown by a solid line while positive OAM is dotted), and it is the interference between these that leads to the interference vortices. The RMT and Qprop OAM distributions are very similar but they differ somewhat from the SFA. In particular, there is a slight shift in the peak positions and relative heights. The highest peak in the SFA is for l_v = ±9, while it is l_v = ±8 for RMT and Qprop. This is related to the larger ‘hole’ in the SFA results (Fig. 2) and most likely can be traced to Coulomb distortions. There is also a shift of around half a photon’s energy between the SFA peaks and those from RMT and Qprop. On each panel in Fig. 3 the ATI peaks are marked by vertical dashed lines, this can be interpreted as the number of photons that must be absorbed to reach a particular momentum. Despite the shift of the peaks between the models, in all cases the order of the ATI peak corresponds approximately to the OAM gained. Thus we are extending the interpretation of the multiphoton picture, where each additional photon gained can add ±1 to the OAM.

Fig. 3 Photoelectron yield |M_lv(p_⊥,p_‖)|² for differing orbital angular momenta l_v. Computed using the SFA [top], Qprop [middle] and the RMT method [bottom]. The same field and target parameters have been used as in Fig. 2. The vertical dashed lines denote OAM regions from solving eqn (29). The left and right columns shows computations for a 2p₀ and 2p₁ initial state, respectively. For positive values of l_v the distribution is marked by data points, while for negative values a line plot is used. This is so that positive and negative OAM can be distinguished as they completely overlap in the m = 0 case. In the case of m = 1, a marked asymmetry is apparent, with states of negative OAM dominating. This reflects the strong-field preference of the m = 1 electron to be ionized by a pulse of opposite (negative) helicity, which drives transitions to negative OAM values. As in Fig. 2 the OAM distributions were computed for p_‖ = 0.1 a.u.

On the right-hand side of Fig. 3 the results for the m = 1, p-state of neon are shown. The main effect (relative to the m = 0 case) is the suppression of the positive OAM values, which is what causes the blurring for the equivalent results in Fig. 2. This asymmetry reflects the preferential ionization of m = 1 electrons by a field of negative helicity. Such a field will tend to decrease the OAM of the ejected electron, and populate states of negative OAM. The asymmetry is a strong-field effect, quite different to the symmetric yields observed in previous studies of few-photon ionization of s electrons. We see in Fig. 3 that in this case states with OAM values symmetric about l_v = 1 contribute over a common range of momenta. In particular, both Qprop and RMT calculations find that the peaks for l_v = −7 and l_v = 9 are now aligned at p = 0.5 a.u. (and similarly l_v = −8 and l_v = 10 contribute at around p = 0.65 a.u.). Despite their asymmetric yields, their interference is still sufficient to give rise to the spiral interference pattern observed in Fig. 2. Thus, it is clear that at a given momentum, the photoelectron yield is contained in states whose OAM value is shifted by +1 relative to the m = 0 case, suggesting that the magnetic quantum number is simply added to the final OAM. For m = −1 electrons (not shown), the opposite tendency appears, and states of positive OAM dominate.

B. Building interference vortices from twisted electrons

Using the distributions presented in Fig. 3 the interference vortices can be reconstructed. This provides particular insight into the meaning of the OAM in a strong field context. To recreate this interference, for the case of m = 0, we will consider a distribution of outgoing photoelectrons with two opposite values of OAM ±l_v. For an initial state of general m one should consider two states with OAM m ± l_v. This can be written asHere |ψ_m±lv,p(t)〉 is the electron vortex state with OAM m ± l_v and momentum p = (p_‖,p_⊥) and w(p) provides a weighting over the cylindrical momentum coordinates. Now, if we project this state onto a 3D plane wave momentum state we can examine the interference in the p_xp_y-plane,The m cancels in the interference fringes, as they depend on the difference in the values of OAM. We have used the following momentum representation of Bessel electron vortex states⁴Thus, we identify the interference between the two vortex states as leading to the condition

4ϕl_v + δp² = 2π(2n + l_v). (28)

This result generalizes the expressions derived in ref. 18 for one-photon and two-photon ionization to the strong field regime, and describes a Fermat spiral with 2l_v arms. In Fig. 4[(a, b) and (d, e)], we plot examples l_v = 1 and l_v = 2, for different values of δ, which inversely sets the growth of the radius with respect to the angle ϕ. These spirals are reminiscent of those presented in ref. 18, and observe a similar delay dependence. In panels (c) and (f) we plot a reconstruction of the full interference vortices by combining pairs of OAM in different momentum regions as indicated by the distributions given in Fig. 3. The SFA and RMT distributions for m = 0 are recreated in panels (c) and (f), respectively. This interference construction will also hold for m = 1, the only difference being positive OAM are suppressed. This explains why the original spiral condition works so well for both cases. For each pair of OAM values the weighting function w(p) is set as a Gaussian over p_⊥ to reflect the position, heights and widths of the respective OAM distributions.

Fig. 4 Interference vortices from twisted electrons. Panels (a) and (b) show interference between photoelectrons with an OAM of ±1, panels (d) and (e) show interference between photoelectrons with an OAM of ±2, the left column uses a delay of 4 cycles, the middle column uses a delay of 1 cycle. The distributions are chosen to be Gaussian in momentum coordinates p_x and p_y. The right column shows the interference of the OAM distributions of the SFA [top] and RMT [bottom]. Gaussian distributions of twisted electron states are taken to reflect the OAM distributions show in Fig. 3. Each panel is normalized with respect to its peak value and plotted on a logarithmic scale.

This clearly demonstrates the link between interference vortices and the OAM of the photoelectrons. However, we have two separate spiral conditions. From the interference of photoelectrons of opposite OAM, eqn (28) leads to a condition for Fermat spirals with a fixed number of arms determined by the OAM. But due to the localization of OAM of the photoelectron, this condition is only valid in specific momentum regions and the number of arms of the interference vortex does not remain fixed. Alternatively, from the SFA, eqn (23), we derived a condition capable of describing a varying number of arms. It is, however, possible to reconcile these two, seemingly disparate descriptions. By employing the semiclassical relation

2ωl_v(p) := p² + 2I_p + 2U_p, (29)

in eqn (23) we recover the same form as eqn (28),

4ϕl_v(p) + δp² = 4π(n − N_cl_v(0)), (30)

where δ = 2πN_c/ω is the number of laser cycles in the delay. Thus, there is a p-dependent semiclassical expression for the OAM, l_v(p), which maps directly onto the ATI equation and can be interpreted as each additional photon contributing ±1 to the OAM. The integer n is shifted by N_c × l(0), where l(0) corresponds to lowest possible OAM given l(0) photons are required to overcome the barrier I_p + U_p. The expression given by eqn (30) maps the SFA condition onto vortex Fermat spirals. It demonstrates, via this semiclassical relation, in specific regions where l(p) is an integer, that eqn (30) will behave like the condition given by eqn (28).

V. Discussion and conclusions

In the previous sections the main dynamics of the interference vortices in the strong field regime were captured and new light was shed on their formation by analysing the orbital angular momenta (OAM). In this section we will place the interference vortices and photoelectron OAM in a wider context, making parallels with other systems as well as discussing the potential to measure the OAM and exploit it in time-resolved imaging.

Vortices and spirals in a physical system are indicative of rotational dynamics and symmetries. When we talk about vortices in quantum matter the physics of liquids is the first example that comes to mind.⁴² In the present context, however, the analogy to Bose Einstein condensates (BEC) of dilute atomic gases is more appropriate.⁴³ Vortices in BEC can be created by rotations (like in a bucket of liquid helium),⁴⁴ that is “stirring with an optical spoon”⁴⁵ or a “phase imprinting”.⁴⁶ While the vortices of charges higher than ±1 are not stable in these systems, rotating BECs form beautiful Abrikosov lattices of charge one vortices. Still, interference of matter waves was, from the very beginning, proposed as a way to observe topological defects,^47–49 stimulated by non-linear optics.^50–53 Interference of a vortex state with a plain BEC leads to the famous fork patterns.^42,43 The interference of two BECs with vortices leads to spiral patterns similar to those discussed here,^42,43 except that (i) the patterns in BEC are formed in the density profile in real space and (ii) they are formed by a macroscopic quantum state, and are thus robust with respect to decoherence.

The question of how to measure the OAM of photoelectrons in strong field ionization may also be answered via interference. The main difference here is that interference structures are formed in the electron momentum space. The interference vortices are in principle measurable in experiments with angular resolution of the electrons created in ionization. Indeed, in the multi-photon regime the experimental measurements clearly showed the Fermat spirals.^22,23 Still, the problem is how sensitive the interference patterns may be to the parameters of the system. Will incoherent effects in the strong field regime blur the interference? Preliminary simulations of focal averaging in our regime of strong and non-perturbative laser fields indicate that the feasibility of observing the interference pattern is robust. The interference fringes do not vary considerably with intensity until U_p > I_p, which for neon at λ = 400 nm occurs at intensities greater than 7 × 10¹⁴ W cm⁻². Going to higher intensities or longer wavelength would perhaps require the use of pre-designed laser pulses with a flat top, for instance,^54,55 using more sophisticated methods. It is also important that the strong field ionization leading to the interference vortices is the dominant process. For lower intensities and multielectron targets, coupling to excited states can occur. Among other effects this could cause superposed interference vortices from both the ground and excited states and the combination can be difficult to disentangle. For larger intensities and longer pulses these effects can be minimized. For more complex targets such as molecules, contributions from multiple states or molecular centres may also present difficulties as again the superposition of many interference vortices, associated with different ionization potentials and magnetic quantum numbers, could make them difficult to resolve.

For the case of circular fields we have presented a clear path to determine the OAM of the photoelectrons, via a second field of opposite helicity. What about elliptical, linear or more complex fields? In the case of linear fields, as they do not interact with the OAM, the same scheme as discussed in this work would not be possible. It would, however, for any field be possible to produce a well-characterized reference for measurement of OAM via interference using a circular pulse. Interference of pairs of twisted electron states will always yield a spiral in the p_xp_y-plane with the number of arms being equal to the difference between the two OAM values. Thus, a careful choice of circular pulses following any target pulse (after a delay) could reveal the OAM distribution. Exploiting the localisation of OAM in momentum space for a circular field would lead to different interference vortices in different regions, in effect ‘scanning’ the OAM of the photoelectrons ionized from the target pulse. However, this method would start to break down if multiple OAMs populate the same region of momentum space for a particular target pulse. Multiple vortices may be difficult to resolve and lead to interferences too complex to analyse.

Another approach is to measure the OAM of the outgoing photoelectrons directly. This would require significant alteration to standard strong field experimental setups, such as velocity map imaging (VMI)^56,57 or a reaction microscope (ReMi) detection system,^58–60 so needs some justification. However, it would not only be a much more generalized and robust way to determine the OAM. Direct measurement could also allow for measurement of incompatible observables, which is not obviously possible using interferometric schemes via tailored laser fields. Incompatible observables are exploited in various aspects of quantum information and metrology including quantum key distribution, Bell inequalities and quantum sensing.^61,62 Additionally, the measurement of OAM could enhance existing strong field procedures. For example, it would allow direct measurement of the quantum magnetic number m of the initial state. Furthermore, it has already been shown that recolliding twisted photoelectrons are sensitive to chiral molecular targets,¹⁵ which may be exploited for imaging. This may also have implications for photoelectron holography,⁶³ where recollision/reinteraction with the target plays an important role and chiral phases could be revealed. The inclusion of fields with high ellipticity and recollision probability, such as bi-circular fields, would enable manipulation of the OAM while retaining strong interaction with the target. Thus, OAM measurements could open up a range of possibilities for strong field research.

So what is the state of the art in OAM measurement and how can this be applied to strong field systems? In ref. 3 and 4 a review of the state of measurement in OAM is given. A range of methodologies are available, an example of a typical approach is using diffraction via holograms or particularly shaped apertures, converting the phase information into spatial information. The first/most well-known method uses a fork phase mask^64,65 but there are difficulties such as the limited transparency of the mask. In ^{ref. 66} the OAM sorter developed transforms the OAM into spatial position, reminiscent of a Stern–Gerlach style measurement of spin. Such a scheme could be envisaged in a strong field set-up. This would forgo the usual momentum information from the p_xp_y-plane in favour of a p_⊥ and l measurement. In theory this could be achieved in existing detection systems with the addition of the hologram/phase mask and corrector as described in the OAM sorter.⁶⁶ This would map one of the dimensions in the p_xp_y-plane to the OAM, leading to different spots along the detector for different values of OAM. However, there are a great many practical considerations including the kinetic energy of the photoelectrons, propagation distance, detector and mask efficiencies to name a small few. In trying to combine these systems together it may be that they are incompatible without major alterations. However, the additional control and information that may be gleaned from strong field systems as a result certainly makes it worth consideration.

Because the OAM gives a direct measurement of the quantum magnetic number as well as providing sensitivity to chiral molecules,^15,35 it is a useful tool for photo-induced time-resolved measurements. The first thought may be in exploiting OAM for such measurements through the incident light such as in twisted attosecond pulses. This is a non-trivial task but progress has been made in high-harmonic generation table top sources^67,68 as well as in free-electron lasers.⁶⁹ However, the focus of this work is on the OAM of the outgoing photoelectron, which may be exploited via detection. As such there would be many possibilities for pump–probe schemes. For example the two counter-rotating circular fields employed in this work could each act as a probe, while two pumps separated by the same delay could be (i) XUV attosecond pulses or (ii) few-cycle IR fields. In the case of (i) we have a situation reminiscent of the RABBITT (reconstruction of attosecond harmonic beating by interference of two-photon transitions) and attosecond streaking techniques.^70,71 In both cases, (i) and (ii), the idea would be to use VMI to exploit the interference, counting the number of spirals to determine the OAM. In order to apply such methods the interference vortices should be further studied in more complex targets to address the above-mentioned reservations for molecules. Direct measurement of the OAM could allow for any pump–probe configuration as well as potentially removing the difficulties of using more complex targets. The ability to employ a variety pump–probe scheme with any ellipticity could allow for non-trivial coupling of the laser pulses to the OAM and novel time-resolved measurements of photo-induced dynamics.

In conclusion, we have presented a description of interference vortices in the strong field regime. Using the strong field approximation we have captured the main physical mechanism, deriving an interference condition, which closely matches calculations performed with the time-dependent Schrödinger equation-based solvers Qprop and R-Matrix with time dependence. We can explain the blurring of interference in the m = 1 case by the asymmetric yield from circular fields of opposite helicity, which is an exclusively strong field effect. By examining the system using orbital angular momentum (OAM) of twisted photoelectron states, we find a new interpretation of the electron vortices in terms of the interference of pairs of OAM states. Not only do we find good agreement between all models here but we also uncover a semiclassical relationship for the OAM by considering the interference between pairs of vortex states and linking it to the condition derived using the SFA. This has consequences, opening up the possibility to measure the OAM directly, or indirectly using interferometric schemes. The OAM of photoelectrons can open the possibility to a range of novel measurements and control in strong field systems.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

This work was in part funded by the UK Engineering and Physical Sciences Research Council (EPSRC). ASM acknowledges grant EP/P510270/1, which is within the remit of the InQuBATE Skills Hub for Quantum Systems Engineering. CFMF would like to acknowledge EPSRC grant EP/T019530/1. GSJA and ACB acknowledge EPSRC grants EP/P022146/1, EP/P013953/1 and EP/R029342/1. The RMT code is part of the UK-AMOR suite, and can be obtained for free at ^{ref. 72}. This work benefited from computational support by CoSeC, the Computational Science Centre for Research Communities, through CCPQ. This work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk), for which access was obtained via the UK-AMOR consortium funded by EPSRC. ICFO group acknowledges support from ERC AdG NOQIA, Spanish Ministry of Economy and Competitiveness (“Severo Ochoa” program for Centres of Excellence in R&D (CEX2019-000910-S), Plan National FISICATEAMO and FIDEUA PID2019-106901GB-I00/10.13039/501100011033, FPI), Fundació Privada Cellex, Fundació Mir-Puig, and from Generalitat de Catalunya (AGAUR Grant No. 2017 SGR 1341, CERCA program, QuantumCAT _U16-011424, co-funded by ERDF Operational Program of Catalonia 2014–2020), MINECO-EU QUANTERA MAQS (funded by State Research Agency (AEI) PCI2019-111828-2/10.13039/501100011033), EU Horizon 2020 FET-OPEN OPTOLogic (Grant No 899794), and the National Science Centre, Poland-Symfonia Grant No. 2016/20/W/ST4/00314.

References

1. L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw and J. P. Woerdman, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 45, 8185.
2. K. Y. Bliokh, Y. P. Bliokh, S. Savel’ev and F. Nori, Phys. Rev. Lett., 2007, 99, 190404.
3. K. Y. Bliokh, I. P. Ivanov, G. Guzzinati, L. Clark, R. Van Boxem, A. Béché, R. Juchtmans, M. A. Alonso, P. Schattschneider, F. Nori and J. Verbeeck, Phys. Rep., 2017, 690, 1.
4. S. M. Lloyd, M. Babiker, G. Thirunavukkarasu and J. Yuan, Rev. Mod. Phys., 2017, 89, 035004.
5. M. Zürch, C. Kern, P. Hansinger, A. Dreischuh and C. Spielmann, Nat. Phys., 2012, 8, 743.
6. G. Gariepy, J. Leach, K. T. Kim, T. J. Hammond, E. Frumker, R. W. Boyd and P. B. Corkum, Phys. Rev. Lett., 2014, 113(15), 153901.
7. R. Géneaux, A. Camper, T. Auguste, O. Gobert, J. Caillat, R. Taïeb and T. Ruchon, Nat. Commun., 2016, 7, 12583.
8. F. Kong, C. Zhang, F. Bouchard, Z. Li, G. G. Brown, D. H. Ko, T. J. Hammond, L. Arissian, R. W. Boyd, E. Karimi and P. B. Corkum, Nat. Commun., 2017, 8, 14970.
9. D. Gauthier, P. R. Ribič, G. Adhikary, A. Camper, C. Chappuis, R. Cucini, L. F. DiMauro, G. Dovillaire, F. Frassetto, R. Géneaux, P. Miotti, L. Poletto, B. Ressel, C. Spezzani, M. Stupar, T. Ruchon and G. De Ninno, Nat. Commun., 2017, 8, 14970.
10. L. Rego, K. M. Dorney, N. J. Brooks, Q. L. Nguyen, C. T. Liao, J. S. Román, D. E. Couch, A. Liu, E. Pisanty, M. Lewenstein, L. Plaja, H. C. Kapteyn, M. M. Murnane and C. Hernández-García, Science, 2019, 364, 1253.
11. E. Pisanty, G. J. Machado, V. Vicuña-Hernández, A. Picón, A. Celi, J. P. Torres and M. Lewenstein, Nat. Photonics, 2019, 13, 569.
12. S. Gemsheim and J.-M. Rost, Phys. Rev. A, 2019, 100, 43408 , arXiv:1909.00728.
13. F. C. Vélez, K. Krajewska and J. Z. Kamiński, Phys. Rev. A, 2018, 97, 043421.
14. F. Cajiao Vélez, J. Z. Kamiński and K. Krajewska, Phys. Rev. A, 2020, 101, 053430.
15. O. I. Tolstikhin and T. Morishita, Phys. Rev. A, 2019, 99, 063415.
16. I. Barth and O. Smirnova, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 84, 063415.
17. I. Barth and O. Smirnova, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 87, 013433.
18. J. M. Ngoko Djiokap, S. X. Hu, L. B. Madsen, N. L. Manakov, A. V. Meremianin and A. F. Starace, Phys. Rev. Lett., 2015, 115, 113004.
19. J. M. Ngoko Djiokap, A. V. Meremianin, N. L. Manakov, S. X. Hu, L. B. Madsen and A. F. Starace, Phys. Rev. A, 2016, 94, 013408.
20. K. J. Yuan, S. Chelkowski and A. D. Bandrauk, Phys. Rev. A, 2016, 93, 053425.
21. J. M. Ngoko Djiokap, A. V. Meremianin, N. L. Manakov, L. B. Madsen, S. X. Hu and A. F. Starace, Phys. Rev. A, 2018, 98, 063407.
22. D. Pengel, S. Kerbstadt, D. Johannmeyer, L. Englert, T. Bayer and M. Wollenhaupt, Phys. Rev. Lett., 2017, 118, 053003.
23. D. Pengel, S. Kerbstadt, L. Englert, T. Bayer and M. Wollenhaupt, Phys. Rev. A, 2017, 96, 043426.
24. M. Li, G. Zhang, X. Kong, T. Wang, X. Ding and J. Yao, Opt. Express, 2018, 26, 878.
25. S. Kerbstadt, K. Eickhoff, T. Bayer and M. Wollenhaupt, Adv. Phys.: X, 2019, 4(1) DOI:10.1080/23746149.2019.1672583.
26. S. Kerbstadt, K. Eickhoff, T. Bayer and M. Wollenhaupt, Nat. Commun., 2019, 10, 1.
27. G. S. J. Armstrong, D. D. A. Clarke, J. Benda, J. Wragg, A. C. Brown and H. W. van der Hart, Phys. Rev. A, 2019, 100, 63416.
28. T. Bayer, C. Philipp, K. Eickhoff and M. Wollenhaupt, Phys. Rev. A, 2020, 102, 013104.
29. K. Amini, J. Biegert, F. Calegari, A. Chacón, M. F. Ciappina, A. Dauphin, D. K. Efimov, C. F. de Morisson Faria, K. Giergiel, P. Gniewek, A. S. Landsman, M. Lesiuk, M. Mandrysz, A. S. Maxwell, R. Moszyński, L. Ortmann, J. A. Pérez-Hernández, A. Picón, E. Pisanty, J. Prauzner-Bechcicki, K. Sacha, N. Suárez, A. Zaïr, J. Zakrzewski and M. Lewenstein, Rep. Prog. Phys., 2019, 82, 116001.
30. V. Mosert and D. Bauer, Comput. Phys. Commun., 2016, 207, 452.
31. V. Tulsky and D. Bauer, Comput. Phys. Commun., 2020, 207, 107098.
32. L. R. Moore, M. A. Lysaght, L. A. A. Nikolopoulos, J. S. Parker, H. W. van der Hart and K. T. Taylor, J. Mod. Opt., 2011, 58, 1132.
33. D. D. A. Clarke, G. S. J. Armstrong, A. C. Brown and H. W. van der Hart, Phys. Rev. A, 2018, 98, 053442.
34. A. C. Brown, G. S. J. Armstrong, J. Benda, D. D. A. Clarke, J. Wragg, K. R. Hamilton, Z. Mašín, J. D. Gorfinkiel and H. W. van der Hart, Comput. Phys. Commun., 2020, 250, 107062.
35. A. Asenjo-Garcia and F. J. García de Abajo, Phys. Rev. Lett., 2014, 113, 066102.
36. W. Becker, F. Grasbon, R. Kopold, D. B. Milošević, G. G. Paulus and H. Walther, Adv. At., Mol., Opt. Phys., 2002, 48, 35.
37. C. Figueira de Morisson Faria, H. Schomerus and W. Becker, Phys. Rev. A: At., Mol., Opt. Phys., 2002, 66, 043413.
38. L. V. Keldysh, Sov. Phys. JETP, 1965, 20, 1307; F. H. Faisal, J. Phys. B: At. Mol. Phys., 1973, 6, L89; H. R. Reiss, Phys. Rev. A: At., Mol., Opt. Phys., 1980, 22, 1786.
39. X. M. Tong and C. D. Lin, J. Phys. B: At., Mol. Opt. Phys., 2005, 38, 2593.
40. P. G. Burke and K. T. Taylor, J. Phys. B: At., Mol. Opt. Phys., 1975, 8, 2620.
41. M. Lewenstein, K. C. Kulander, K. J. Schafer and P. H. Bucksbaum, Phys. Rev. A: At., Mol., Opt. Phys., 1995, 51, 1495.
42. A. Fetter, Rev. Mod. Phys., 2009, 81, 647.
43. L. Pitaevskii and S. Stringari, Bose-Einstein Condensation and Superfluidity, Oxford Univ. Press, Oxford, 2nd edn, 2016.
44. K. W. Madison, F. Chevy, W. Wohlleben and J. Dalibard, Phys. Rev. Lett., 2000, 84, 806.
45. R. Srinivasan, Pramana, 2006, 66, 3.
46. Ł. Dobrek, M. Gajda, M. Lewenstein, K. Sengstock, G. Birkl and W. Ertmer, Phys. Rev. A: At., Mol., Opt. Phys., 1999, 60, R3381.
47. E. L. Bolda and D. F. Walls, Phys. Rev. Lett., 1998, 81, 5477.
48. J. Tempere and J. T. Devreese, Solid State Commun., 1998, 108, 993.
49. F. Dalfovo, S. Giorgini, L. P. Pitaevskii and S. Stringari, Rev. Mod. Phys., 1999, 71, 463.
50. K. Staliunas, C. O. Weiss, and G. Slekys, Horizons of World Physics, Nova Science Publishers, Commack, NY, 1999, vol. 228.
51. I. Basistiy, V. Bazhenov, M. Soskin and M. Vasnetsov, Opt. Commun., 1993, 103, 422.
52. L. V. Kreminskaya, M. S. Soskin and A. I. Khizhnyak, Opt. Commun., 1998, 145, 377.
53. J. F. Nye and M. V. Berry, Proc. R. Soc. London, Ser. A, 1974, 336, 165.
54. W. Boutu, T. Auguste, O. Boyko, I. Sola, P. Balcou, L. Binazon, O. Gobert, H. Merdji, C. Valentin, E. Constant, E. Mével and B. Carré, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 84, 063406.
55. E. S. Toma, P. Antoine, A. de Bohan and H. G. Muller, J. Phys. B: At., Mol. Opt. Phys., 1999, 32, 5843.
56. A. T. J. B. Eppink and D. H. Parker, Rev. Sci. Instrum., 1997, 68, 3477,  DOI:10.1063/1.1148310.
57. M. Takahashi, J. P. Cave and J. H. D. Eland, Rev. Sci. Instrum., 2000, 71, 1337,  DOI:10.1063/1.1150460.
58. R. Moshammer, M. Unverzagt, W. Schmitt, J. Ullrich and H. Schmidt-Böcking, Nucl. Instrum. Methods Phys. Res., Sect. B, 1996, 108, 425.
59. R. Dörner, V. Mergel, O. Jagutzki, L. Spielberger, J. Ullrich, R. Moshammer and H. Schmidt-Böcking, Phys. Rep., 2000, 330, 95.
60. J. Ullrich, R. Moshammer, A. Dorn, R. Dörner, L. P. H. Schmidt and H. Schmidt-Böcking, Rep. Prog. Phys., 2003, 66, 1463.
61. T. Heinosaari, T. Miyadera and M. Ziman, J. Phys. A: Math. Theor., 2016, 49, 123001.
62. K. Bongs, M. Holynski, J. Vovrosh, P. Bouyer, G. Condon, E. Rasel, C. Schubert, W. P. Schleich and A. Roura, Nat. Rev. Phys., 2019, 1, 731.
63. C. F. d. M. Faria and A. S. Maxwell, Rep. Prog. Phys., 2020, 83, 034401.
64. G. Guzzinati, L. Clark, A. Béché and J. Verbeeck, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 89, 025803.
65. K. Saitoh, Y. Hasegawa, K. Hirakawa, N. Tanaka and M. Uchida, Phys. Rev. Lett., 2013, 111, 074801.
66. V. Grillo, A. H. Tavabi, F. Venturi, H. Larocque, R. Balboni, G. C. Gazzadi, S. Frabboni, P. H. Lu, E. Mafakheri, F. Bouchard, R. E. Dunin-Borkowski, R. W. Boyd, M. P. J. Lavery, M. J. Padgett and E. Karimi, Nat. Commun., 2017, 8, 15536.
67. K. M. Dorney, L. Rego, N. J. Brooks, J. San Román, C.-T. Liao, J. L. Ellis, D. Zusin, C. Gentry, Q. L. Nguyen, J. M. Shaw, A. Picón, L. Plaja, H. C. Kapteyn, M. M. Murnane and C. Hernández-García, Nat. Photonics, 2019, 13, 123.
68. J. W. Wang, M. Zepf and S. G. Rykovanov, Nat. Commun., 2019, 10, 5554.
69. P. Rebernik Ribič, B. Rösner, D. Gauthier, E. Allaria, F. Döring, L. Foglia, L. Giannessi, N. Mahne, M. Manfredda, C. Masciovecchio, R. Mincigrucci, N. Mirian, E. Principi, E. Roussel, A. Simoncig, S. Spampinati, C. David and G. De Ninno, Phys. Rev. X, 2017, 7, 031036.
70. M. J. J. Vrakking, Phys. Chem. Chem. Phys., 2014, 16, 2775.
71. L. Cattaneo, J. Vos, M. Lucchini, L. Gallmann, C. Cirelli and U. Keller, Opt. Express, 2016, 24, 29060.
72. The RMT repository, https://gitlab.com/Uk-amor/RMT/rmt, 2020.

---