1 Introduction

Since being first recognised in the classical context of tides [1,2,3] vortex phenomena have held an iconic status across a diverse range of disciplines [4,5,6,7]. Of particular interest are the quantum mechanical versions of vortices, which can be found in a wide range of systems [8,9,10,11,12]. Current intense interest in phase vortices follows the first experimental observations of this phenomenon for photons [8] and unbound electrons [10] (for reviews on vortices in electrons see [13,14,15]). The unique topological properties of vortex states [15] render them fundamental in the study of structured wave fields, and have led to the inception of whole research areas, such as singular optics [16,17,18,19]. For instance, a vortex state cannot transform to another by simple deformation such as stretching and compressing, or the addition of noise. Furthermore, at its centre along the propagation axis, a vortex beam has a zero amplitude and an ill-defined phase. Currents around singularities imply that such states carry intrinsic orbital angular momenta (OAM), and thus can be used to influence the dynamical properties of physical systems [20, 21].

The study of vortices in attosecond physics has great potential in controlling light and matter. In high-harmonic generation (HHG) it has been shown that optical vortices in the IR driving fields lead to optical vortices in the resulting UV light produced [22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]. This has been exploited to allow a high degree of control over the light. In one example, UV light has been produced exhibiting torus knot topology [34, 35] by sophisticated trefoil IR pulses, for which the orientation of the trefoil varies with the azimuthal angle. In another study  [36], extreme UV (EUV) light was imparted with time-varying OAM, leading to self-torque by employing two time-delayed IR pulses with different OAM. This was the first demonstration of self-torque in light [37] and could aid in probing systems with naturally time-varying OAM.

Work understanding the OAM of photoelectrons emitted in attosecond processes is still in its infancy. Initial studies include the exploration of high OAM values for quasi-relativistic field intensities [38] and terahertz fields [39] as well as using the OAM in rescattering to probe bound state structures [40]. The OAM has also been studied indirectly via the coherent combination of pairs of vortex states that lead to the interference vortices [41,42,43,44,45,46,47,48,49,50,51], which occur for two counter rotating circularly polarised fields separated by a time delay. However, ideas as advanced as knots or self-torque in electron vortices have not been explored. This can partly be attributed to the difficulty in experimental implementation; for example, no measurement scheme has been devised to detect the OAM of photoelectrons emitted in strong-field experiments. Initial ideas on how to achieve this have been suggested in our recent publication [51], using both interferometric schemes and adapting existing methods [52] used in electron beams. Such experimental development needs sufficient theoretical backing to guide the implementation and justify the cost, which at present is lacking. In this study we introduce a theoretical framework to derive analytical conservation laws and understand the basic dynamics of the OAM during strong field ionisation. Note in all cases we consider the OAM of the electron but not that of the laser fields. The OAM we refer to in this work is associated with the OAM carried by free particles, described by Bessel functions and is the eigenvalue of \(\hat{L}_z\), that is, the quantum magnetic number. The outgoing photo electron wavepacket in strong-field ionisation does not naturally form a beam, so the z-direction can be arbitrarily chosen. Throughout this work we use the convention chosen in Fig. 1.

We utilise the strong field approximation (SFA) [53]—often considered the workhorse of strong field physics—to describe the basic ionisation dynamics. The SFA in the form used here [51, 54, 55] is an approximation, which primarily neglects the effects of the Coulomb potential in the continuum. The use of the SFA, over a more accurate model, is justified by the possibility of analytical solutions and the ease of interpreting the results. Thus, the SFA allows the basic laws of the OAM to be derived for strong field ionisation. Furthermore, we focus mostly on circularly polarised light, where the Coulomb effects are less significant, than those observed in linear polarisation [56]. In a recent publication [51], we numerically computed OAM distributions using the SFA, QProp [57, 58] and the R-matrix with time dependence method [59,60,61] for circularly polarised fields and confirmed that the SFA was able to qualitatively reproduce the key features found in these two more accurate numerical methods. This work differed from the present in that it focused on the interference effects due to employing two time-delayed counter-rotating circularly polarised IR fields. Furthermore, the OAM computations were performed numerically without full derivation of the analytical expressions and the SFA calculations employed a monochromatic field, whereas in this work we extend the model to include a \(\sin ^2\) envelope.

The paper is structured as follows. In Sect. 2 key theoretical results of the SFA (Sect. 2.1) and vortex states (Sect. 2.2) are given. Next, in Sect. 3 we incorporate the OAM into the SFA, starting with expressions for a general field (Sect. 3.2). Following this, we derive analytical conditions for the cases of a linear field (Sect. 3.3) and a monochromatic circular field; in the latter we do this both without (Sect. 3.4) and with (Sect. 3.5) the use of the saddle point approximation. In Sect. 4 numerical results are presented for a circular \(\sin ^2\) laser field. Therein we derive an effective frequency to extend the analytic condition for a monochromatic field to work for the \(\sin ^2\) case. Finally, in Sect. 5 we present our conclusions.

2 Key results from the strong field approximation and vortex states

In this section we provide some key results from the SFA and for electron vortex states, necessary for understanding the OAM derivation. Throughout the article we use atomic units unless otherwise stated. Both cylindrical \((\hat{e}_{\parallel },\hat{e}_{\perp },\hat{e}_{\phi } )\) and Cartesian \((\hat{e}_{x},\hat{e}_{y},\hat{e}_{z} )\) coordinates will be used in this article, they will always be aligned such that \(\hat{e}_{\parallel }=e_{z}\) and any radial quantity will be given by e.g. \(r_{\perp }=\sqrt{r_x^2+r_y^2}\). We will consider specific cases where the laser field polarisation is parallel to \(\hat{e}_{\parallel }\) [linearly polarised] and in the xy-plane [circularly polarised].

2.1 Key results from SFA

Within the SFA S-matrix formalism [54, 62] the transition amplitude for direct strong field ionisation from the bound state \(|{\psi _0(t')}\rangle \) to a final continuum state \(|{\psi _{f}(t)}\rangle \) is given by the following expression [54, 62]

$$\begin{aligned} {} M_{f}=-i\lim \limits _{t\rightarrow \infty }\int _{-\infty }^{t} \mathrm{d}t'\langle {\psi _{f}(t)}|{U_{v}(t,t')V}|{\psi _{0}(t')}\rangle , \end{aligned}$$
(1)

where the time evolution is approximated by the Volkov operator

$$\begin{aligned} U_{v}(t,t')=\int d^{3}\mathbf {p}e^{\frac{-i}{2}\int _{t'}^{t}(\mathbf {p}+\mathbf {A}(\tau )^{2})\mathrm{d}\tau }|{\mathbf {p}+\mathbf {A}(t)}\rangle \langle {\mathbf {p}+\mathbf {A}(t')}|. \end{aligned}$$
(2)

Combining Eqs. (2) and (1) and taking \(\mathbf {A}(t)=0\) we get the following for the transition amplitude

$$\begin{aligned} M_{f}&=\lim \limits _{t\rightarrow \infty } \int d^{3}\mathbf {p}'\; e^{-iS(\mathbf {p},t)}\langle {\psi _{f}|\mathbf {p}'}\rangle M(\mathbf {p}) \end{aligned}$$
(3)

with

$$\begin{aligned} M(\mathbf {p})&= -i\int _{-\infty }^{\infty }\mathrm{d}t'e^{iS(\mathbf {p},t')}d(\mathbf {p},t'), \end{aligned}$$
(4)

where

$$\begin{aligned} d(\mathbf {p},t')&=\langle {\mathbf {p}+\mathbf {A}(t')|V|\psi _{0}}\rangle \end{aligned}$$
(5)

and \(S(\mathbf {p},t)\) and \(S(\mathbf {p},t')\) are the upper and lower limit of the semi-classical action, respectively, both given by

$$\begin{aligned} S(\mathbf {p},t)=\mathrm {I}_\mathrm {p}t +\frac{1}{2}\int ^{t}_{-\infty } \mathrm{d}\tau (\mathbf {p}+\mathbf {A}(\tau ))^2. \end{aligned}$$
(6)

A plane wave momentum state is commonly used for the final continuum state \(|{\psi _f}\rangle \rightarrow |{\psi _{\mathbf {p}}}\rangle \), which leads to \(\langle {\psi _f|\mathbf {p}'}\rangle \rightarrow \langle {\mathbf {p}|\mathbf {p}'}\rangle =\delta (\mathbf {p}'-\mathbf {p})\) and therefore \(M_f\rightarrow M(\mathbf {p})\), where Eq. (4) gives the definition of \(M(\mathbf {p})\). However, in this work, a Bessel beam vortex state will be used as final continuum state instead of a plane wave. This will enable the development of an analytical model for the OAM of the outgoing photoelectron. In order to proceed we will introduce some properties of the Bessel beam.

2.2 Key results from electron vortex state

Vortex states are topologically distinct both from plane waves and vortex states with a different orbital angular momentum l [15]; this means two vortex states cannot be transformed into one another via continuous deformations. At its centre along the propagation axis, the vortex has a zero amplitude and undefined phase. In order to enforce the single-valued wave function, a quantised phase factor \(e^{il\phi }\) is needed, where \(\phi \) is the azimuthal angle and l is a topological charge with integer value known as the orbital angular momentum (OAM). The general form of the Bessel beam electron-vortex is [14]

$$\begin{aligned} \langle {\mathbf {r}|\psi _{l}(t)}\rangle =N_{l}J_{l}(p_{\perp }r_{\perp })e^{il\phi }e^{i p_{\parallel }r_{\parallel }}e^{-i\omega t}. \end{aligned}$$
(7)

Here, \(N_{l}\) is a normalisation factor and \(J_l(p_{\perp }r_{\perp })\) is the Bessel function of the first kind. Ignoring the time dependence, the Fourier transform is

$$\begin{aligned} \langle {\mathbf {p}'|\psi _{l}}\rangle =\frac{i^{-l}e^{i l\phi '}}{2\pi p_{\perp }}\delta (p'_{\parallel }-p_{\parallel }) \delta (p'_{\perp }-p_{\perp }), \end{aligned}$$
(8)

where \((p_{\parallel },p_{\perp },\phi )\) are the cylindrical coordinates of \(\mathbf {p}\). Inserting this into Eq. (3) gives

$$\begin{aligned}&M_{l}(p_{\parallel },p_{\perp },t) \nonumber \\&\quad =\frac{i^{l}}{2\pi }\int _{-\pi }^{\pi } \mathrm{d}\phi 'e^{-i S(p_{\parallel },p_{\perp },\phi ',t)} e^{-il\phi '} M(p_{\parallel },p_{\perp },\phi '), \end{aligned}$$
(9)

where \(M(p_{\parallel },p_{\perp },\phi ')\) is the transition amplitude Eq. (5) written in cylindrical coordinates, which have been used as they are natural for vortex states. Now the vortex state is incorporated into the SFA framework. Next, we will compute more explicit expressions for the transition amplitude.

3 Matrix element calculations

In this section, the transition matrix element will be derived for a general laser field. We will examine the specific cases of a linear field and a circular monochromatic field, as well as the saddle point approximation.

3.1 Coordinate systems

In Fig. 1 the two laser fields considered in this work are depicted. The circular field propagates along the z or \(\hat{e}_{||}\) direction, while the linear field propagates along the x direction. The angle relating to the final OAM of the photoelectron is chosen to be fixed for both fields around the z-axis, while its conjugate variable, the angle \(\phi \), is marked on the figure. This convention will be used throughout the article. Note that, for any vector \(\mathbf {v}\), the notation \(v_{\perp }\) is used for the radial component in a cylindrical coordinate system, which is the projection of the vector onto the xy-plane given by \(v_{\perp }=\sqrt{v_x^2+v_y^2}\).

Fig. 1
figure 1

Depiction of the two laser fields and the coordinate system used throughout this article. The circularly polarised field propagates along the z-axis and its electric field vector rotates in the xy-plane (marked by a square grid). In contrast, the linearly polarised field propagates along the x-axis and its electric field vector oscillates along the z/\(\hat{e}_{||}\) direction. The angle \(\phi \) is defined to be in the xy-plane (i.e. the angle about the z-axis) and starts from the x-axis. The orbital angular momentum is defined to be around the z-axis as marked on the figure by \(\hat{L}_z\)

Full size image

3.2 Matrix element calculations for a general field

In order to proceed we will collect \(\phi \) dependent and independent parts to analytically perform the integral in Eq. (9). The upper limit of the action contains \(\phi \) independent terms that contribute only a phase and thus may be neglected as well as the following \(\phi \) dependent term

$$\begin{aligned} p_{\perp } (\cos (\phi ) \alpha _x(t) + \sin (\phi ) \alpha _y(t)), \end{aligned}$$
(10)

where \(\varvec{\alpha }(t)=\int _{-\infty }^{t} \mathbf {A}(\tau ) \mathrm{d}\tau \). As \(t\rightarrow \infty \) for a ‘well-behaved’ pulse or monochromatic field this term will vanish. For the case of \(\sin ^2\) pulses, which will be considered later, this holds for an N-cycle pulse, where N is an integer greater than one. To simplify matters we will only consider such ‘well-behaved’ laser fields. Thus, the upper limit of the action can be entirely neglected here.

The lower limit of the action, which governs the dynamics, can be split into two parts:

$$\begin{aligned}&S(p_{\parallel },p_{\perp },\phi ,t')= S_A(p_{\parallel },p_{\perp },t')+S_B(p_{\perp },\phi ,t'), \end{aligned}$$

where

$$\begin{aligned}&S_A(p_{\parallel },p_{\perp },t')\nonumber \\&\quad = \left( \mathrm {I}_\mathrm {p}+\frac{1}{2}p_{\parallel }^2+\frac{1}{2}p_{\perp }^2\right) t' +p_{\parallel }\alpha _{\parallel }(t') +\frac{1}{2}\int \mathrm{d} t'A^2(t') \end{aligned}$$
(11)

is independent of \(\phi \), and

$$\begin{aligned} S_B(p_{\perp },\phi , t')&= p_{\perp } \left( \alpha _x(t')\cos (\phi )+\alpha _y(t')\sin (\phi )\right) \nonumber \\&= p_{\perp }\alpha _{\perp }(t')\sin (\phi -\nu (t')), \end{aligned}$$
(12)

encodes the dependence on the momentum azimuth angle \(\phi \). This dependence is in terms of the angle

$$\begin{aligned} \nu (t')&= \mu (t') - \pi /2, \end{aligned}$$
(13)

where

$$\begin{aligned} \mu (t')&=\arctan (\alpha _y(t')/\alpha _x(t')), \end{aligned}$$
(14)

which is related to the rotation of the laser field and magnitude

$$\begin{aligned} \alpha _{\perp }(t')&=\sqrt{\alpha _x(t')^2+\alpha _y(t')^2} \end{aligned}$$
(15)

of the field integral \(\varvec{\alpha }(t')\).

Finally, in order to remove the \(\phi \) dependence from the bound-state matrix element, \(d(\mathbf {p},t')=\langle {\mathbf {p}+\mathbf {A}(t')|V|\psi _{0}}\rangle \), we decompose it into its Fourier series

$$\begin{aligned} d(\mathbf {p},t')=\sum _{m}e^{i m \phi } V_{m}(p_{\parallel },p_{\perp },t'). \end{aligned}$$
(16)

If the bound state is that of an atomic target, then the matrix element will have only a limited number of nonzero relevant Fourier terms, and these will be easy to compute numerically applying the fast Fourier transform (FFT). The transition amplitude from Eq. (9) may be written as

$$\begin{aligned}&M_{l}(p_{\parallel },p_{\perp },t)\nonumber \\&\quad = -i\sum _{m}\int _{-\infty }^{t}\mathrm{d}t' e^{i S_A(p_{\parallel },p_{\perp },,t')} V_{m}(p_{\parallel },p_{\perp },t') I^{\phi }_l(p_{\perp },t'), \end{aligned}$$
(17)

where \(I^{\phi }_l(p_{\perp },t')\) contains the \(\phi \) integral and corresponding terms and is given by

$$\begin{aligned} I^{\phi }_l(p_{\perp },t')&=\frac{i^l}{2\pi }\int _{-\pi }^{\pi }\mathrm{d}\phi ' e^{-i(l-m)\phi '} e^{i p_{\perp }\alpha _{\perp }(t')\sin \left[ \phi '-\nu (t')\right] }\nonumber \\&=i^{l}e^{-i (l-m) \nu (t')} J_{l-m}(p_{\perp }\alpha _{\perp }(t')). \end{aligned}$$
(18)

Now the full transition amplitude can be written as

$$\begin{aligned}&M_{l}(p_{\parallel },p_{\perp })= i^{l-1}\sum _m \int _{-\infty }^{\infty }\!\!\!\mathrm{d}t' e^{i S_{l-m}(p_{\parallel },p_{\perp },t')}\nonumber \\&\quad V_{m}(p_{\parallel },p_{\perp },t') J_{l-m}(p_{\perp }\alpha _{\perp }(t')), \end{aligned}$$
(19)

where

$$\begin{aligned} S_{l}(p_{\parallel },p_{\perp },t')&=\left( \mathrm {I}_\mathrm {p}+\frac{1}{2}p^2_{\parallel }+\frac{1}{2}p^2_{\perp }\right) t' -l \nu (t')\nonumber \\&\quad +\,p_{\parallel }\alpha _{\parallel }(t') +\frac{1}{2}\int \mathrm{d} t'A^2(t'). \end{aligned}$$
(20)

Note that the OAM l only appears in the action coupled with \(\nu (t')\), which mediates interaction of the laser field with the OAM of the electron. The time-varying ‘angle’ \(\nu (t')\) can be interpreted as the dynamical rotational action of the field on the OAM of the photoelectron. For a monochromatic field, we will see that \(\nu (t')=\omega t\). We can rewrite \(\nu (t')\) in terms of its derivative to gain more insight

$$\begin{aligned} \nu '(t') = \frac{\left( \varvec{\alpha }(t')\times \mathbf {A}(t')\right) \cdot \hat{e}_{\parallel }}{\alpha _{\perp }(t')^2}. \end{aligned}$$
(21)

This makes the rotational nature and link to the spin of the field [63] more apparent. We will now examine a few special cases that will significantly simplify the transition amplitude.

3.3 Matrix element calculations for a linear field

In this section we will consider a laser field, which only has components in the \(\hat{e_{\parallel }}\) direction in the cylindrical coordinate system (i.e. perpendicular to the xy-plane see, Fig. 1 for more details). Thus, the field has no \(\phi \) dependence and only \(S_A(p,\theta _p,t,t')\) contributes to the action. Furthermore, given that \(\alpha _{\perp }(t')\rightarrow 0\), the Bessel function will become a Kronecker delta

$$\begin{aligned} J_{l-m}(p_{\perp }\alpha _{\perp }(t'))\rightarrow \delta _{lm}, \end{aligned}$$

and this can be demonstrated by re-evaluating Eq. (18)

$$\begin{aligned} I^{\mathrm {linear}}_{\phi }=\frac{1}{2\pi }\int _{-\pi }^{\pi } e^{-i(m-l)\phi '}\mathrm{d}\phi ' =\delta _{l,m}. \end{aligned}$$
(22)

This is a selection rule, enforcing \(l=m\), given by the symmetry of the problem. Substituting this back into the transition amplitude leaves

$$\begin{aligned} M_{l}(p_{\parallel },p_{\perp },)=i^{l-1}\int _{-\infty }^{\infty }\mathrm{d}t' e^{i S_A(p_{\parallel },p_{\perp }, t')} V_{l}(p_{\parallel },p_{\perp },t'). \end{aligned}$$
(23)

The selection rule comes from the exponential factor of the bound state and OAM of the free electron, which indicates that there will be a one-to-one correspondence between the two. Thus, if there is only one Fourier term \(m=m_0\) in Eq. (16), the OAM of the photoelectron will be \(l=m_0\). Aside from this conservation of angular momentum, the form of the OAM distribution is unchanged from the original SFA formalism for the plane wave momentum state.

3.4 Matrix element calculations for a circular field

In this section we turn to the OAM transition amplitude for a monochromatic circular field, where there is now a \(\phi \) dependence in the semi-classical action. The form of vector potential for a circular monochromatic field is given by

$$\begin{aligned} \mathbf {A}(t)=-\sqrt{2 \mathrm {U}_\mathrm {p}} \left( \cos (\omega t)\varvec{e}_{x} +\sin (\omega t)\varvec{e}_{y}\right) , \end{aligned}$$
(24)

for more details on the orientation of the field see Fig. 1. For this field in Eq. (19) \(\alpha _{\perp }(t') = \sqrt{2\mathrm {U}_\mathrm {p}}/\omega \) and \(\nu (t')=\omega t'\). Thus, the action from Eq. (20) becomes linear in \(t'\) such that

$$\begin{aligned} S_{l-m}(p_{\parallel },p_{\perp },t')&=\chi _{l-m}(p_{\parallel },p_{\perp }) t', \end{aligned}$$
(25)

where

$$\begin{aligned} \chi _{l-m}(p_{\parallel },p_{\perp })&=\mathrm {I}_\mathrm {p}+\mathrm {U}_\mathrm {p}+\frac{1}{2}(p^2_{\parallel }+p^2_{\perp }) -(l-m)\omega . \end{aligned}$$
(26)

The transition amplitude is then

$$\begin{aligned} M_{l}(p_{\parallel },p_{\perp })&=i^{l-1}\sum _m J_{l-m}\left( p_{\perp }\sqrt{2\mathrm {U}_\mathrm {p}}/\omega \right) \nonumber \\&\quad \times \int _{-\infty }^{\infty }\mathrm{d}t' e^{i \chi _{l-m}(p_{\parallel },p_{\perp }) t'} V_{m}(p_{\parallel },p_{\perp },t'). \end{aligned}$$
(27)

The integral acts to Fourier transform the prefactor term \(\tilde{V}_{\mathbf {p}0}\) to give

$$\begin{aligned} M_{l}(p_{\parallel },p_{\perp })&=i^{l-1}\sum _m J_{l-m}\left( p_{\perp }\sqrt{2\mathrm {U}_\mathrm {p}}/\omega \right) \nonumber \\&\quad \times \hat{V}_{m}(p_{\parallel },p_{\perp },\chi _{l-m}) \end{aligned}$$
(28)

where \(\hat{V}_{m}(p_{\parallel },p_{\perp },\chi _{l-m})\) is the Fourier transform of \(V_{m}(p_{\parallel },p_{\perp },t')\) with the frequency \(\chi _{l-m}(p_{\parallel },p_{\perp })\).

To further simplify the problem we will now consider the case with a very simple bound state, where the matrix elements have only one Fourier series term for \(m=0\), as would be the case for a simple s-state or the bound state for a zero range potential. In the example of a zero range potential [64] the matrix element of the bound state is given simply by a constant dependent on the ionisation potential \(V_0 = \sqrt{2\pi \sqrt{2 \mathrm {I}_\mathrm {p}}}\). Now Eq. (27) becomes

$$\begin{aligned} M_{l}(p_{\parallel },p_{\perp })&=2\pi V_0 i^{l-1} J_{l-m}\left( p_{\perp }\sqrt{2\mathrm {U}_\mathrm {p}}/\omega \right) \nonumber \\&\quad \times \delta \left( \chi _{l-m}(p_{\parallel },p_{\perp })\right) . \end{aligned}$$
(29)

Thus, for the case of monochromatic field with a simple bound state we arrive at another conservation equation

$$\begin{aligned} \mathrm {I}_\mathrm {p}+\mathrm {U}_\mathrm {p}+\frac{1}{2}(p^2_{\parallel }+p^2_{\perp }) -(l-m)\omega =0. \end{aligned}$$
(30)

This is directly related to the semi-classical condition for ATI peaks [65], which can be interpreted forming due to additional photons absorbed by the photoelectron in the continuum, beyond that required for ionisation. In this case, however, each ATI peak will correspond to a different value of OAM, which will be shifted by the quantum magnetic number m. This has interesting implications as it suggests that different values of the OAM will be localised to specific energy regions in photoelectron emission spectra for ionisation via circularly polarised light. It is this idea that was exploited to produce interference vortices in our recent work [51].

3.5 Exploiting the saddle point approximation

In this subsection we discuss an alternative way to compute the OAM-SFA transition amplitude. In order to compute the OAM transition amplitude we employ the saddle point approximation for Eq. (5) and then analytically perform the integral over \(\phi \). Applying the saddle point approximation and ignoring the bound state prefactor Eq. (5) becomes

$$\begin{aligned} M(p_{\parallel },p_{\perp },\phi )&= \sum _{s=0}^{N-1}\sqrt{\frac{2\pi i}{\partial ^2 S/\partial t^2|_{t=t_s}}} \exp \left[ i S(p_{\parallel },p_{\perp },\phi ,t_s) \right] , \end{aligned}$$
(31)

where N is the number of laser cycle considered and \(t_s\) is given by

$$\begin{aligned} \left( \mathbf {p}+\mathbf {A}(t_s) \right) ^2=-2\mathrm {I}_\mathrm {p}. \end{aligned}$$
(32)

For a general field \(t_s\) will depend on \(\phi \) in a nontrivial way so the integral over \(\phi \) cannot be performed analytically. However, for a monochromatic circular field the dependence on \(\phi \) is linear and we can write \(t_s\) in terms of \(t'_s\), which is independent of \(\phi \)

$$\begin{aligned} \omega t_s(p_{\parallel },p_{\perp },\phi )&= \phi +\omega t'_s(p_{\parallel },p_{\perp })&\quad (\mathrm {mod}\quad 2\pi ). \end{aligned}$$
(33)

Substituting this into the action leads to

$$\begin{aligned} S(p_{\parallel }, p_{\perp }, \phi ,t'_s)&= \left( \mathrm {I}_\mathrm {p}+\mathrm {U}_\mathrm {p}+\frac{1}{2}(p_{\parallel }^2 + p_{\perp }^2) \right) \left( \frac{\phi }{\omega }+t'_s \right) \nonumber \\&\quad -\,\frac{\sqrt{2\mathrm {U}_\mathrm {p}}}{\omega }p_{\perp } \sin \left( \omega t'_s\right) \end{aligned}$$
(34)
$$\begin{aligned}&=S(p_{\parallel },p_{\perp },0,t'_s)\nonumber \\&\quad +\,\left( \mathrm {I}_\mathrm {p}+\mathrm {U}_\mathrm {p}+\frac{1}{2}(p_{\parallel }^2 + p_{\perp }^2) \right) \frac{\phi }{\omega }. \end{aligned}$$
(35)

Thus, substituting the action and transition amplitude into Eq. (19) leads to the expression

$$\begin{aligned} M_{l}(p_{\parallel },p_{\perp })&=\sum _{s=0}^{N-1} \sqrt{\frac{2\pi i}{\partial ^2 S/\partial {t'}_{s}^2}} \exp (i S(p_{\parallel },p_{\perp },0,t'_s))\nonumber \\&\quad \times \frac{1}{2\pi }\int _{-\pi }^{\pi }d\phi ' \exp \Bigg [ i (\mathrm {I}_\mathrm {p}+\mathrm {U}_\mathrm {p}+1/2 p^2)\nonumber \\&\quad \times (\phi '/\omega )-i\phi ' l \Bigg ] \end{aligned}$$
(36)
$$\begin{aligned}&=\sum _{s=0}^{N-1} \sqrt{\frac{2\pi i}{\partial ^2 S/\partial {t'}_{s}^2}} \exp (i S(p_{\parallel },p_{\perp },0,t'_s))\nonumber \\&\quad \times \mathrm {sinc}\left[ (\mathrm {I}_\mathrm {p}+\mathrm {U}_\mathrm {p}+1/2 p^2 -\omega l)(\pi /\omega ) \right] . \end{aligned}$$
(37)

Note, the prefactor in the square root can be determined to be independent of \(\phi '\), thus it is outside of the \(\phi '\) integral. The sum leads to the contribution of one identical ionisation event per laser cycle, due to the field being monochromatic. Thus, the sum can be evaluated analytically to give the following

$$\begin{aligned} M_{l}(p_{\parallel },p_{\perp })&= \mathrm {\Omega }_{N}(p_{\parallel },p_{\perp }) \sqrt{\frac{2\pi i}{\partial ^2 S/\partial {t'}_{s}^2}} \exp (i S(p_{\parallel },p_{\perp },0,t'_0))\nonumber \\&\quad \times \mathrm {sinc}\left[ \frac{\chi _l\pi }{\omega }\right] , \end{aligned}$$
(38)

where \(t'_0\) denotes the solution in the first laser cycle and

$$\begin{aligned} \mathrm {\Omega }_{N}(p_{\parallel },p_{\perp })&= \sum _{n=0}^{N-1} \exp \left( 2\pi i n \chi _0/\omega \right) \nonumber \\&=\frac{\exp \left[ 2\pi i N\chi _0/\omega \right] -1}{\exp \left[ 2\pi i \chi _0/\omega \right] -1} \end{aligned}$$
(39)

as previously demonstrated in [66, 67]. In the limit \(N\rightarrow \infty \) \(\mathrm {\Omega }_{N}(p_{\parallel },p_{\perp })\) becomes a Dirac comb

$$\begin{aligned} \lim _{N\rightarrow \infty } \mathrm {\Omega }_{N}(p_{\parallel },p_{\perp }) =\sum _{n=0}^{\infty } \delta \left( \chi _0- n\omega \right) . \end{aligned}$$
(40)

The Dirac delta functions lead to the following replacement in the argument of the \(\mathrm {sinc}\) function

$$\begin{aligned} \mathrm {sinc}\left[ \frac{\chi _l\pi }{\omega }\right]&\rightarrow \mathrm {sinc}\left[ \pi (n-l)\right] =\delta _{nl}. \end{aligned}$$
(41)

Thus, this leads to

$$\begin{aligned} M_{l}(p_{\parallel },p_{\perp })&= \sqrt{\frac{2\pi i}{\partial ^2 S/\partial {t'}_{s}^2}} \exp (i S(p_{\parallel },p_{\perp },0,t'_0))\nonumber \\&\quad \times \delta \left( \chi _l(p_{\parallel },p_{\perp })\right) , \end{aligned}$$
(42)

which gives the same condition as Eq. (30) but has different prefactors due to the saddle point approximation. Taking N as a finite fixed value we are able to plot the result, which can be considered a rough approximation to employing a laser pulse of N cycles, which has been done in Fig. 2. The plot is over energy to make clear that each peak is separated with equal energy spacing as dictated by Eq. (30). As the number of laser cycles goes from 2 to 8 the OAM peaks get sharper, approaching the delta functions predicted by Eq. (29). The bound state prefactor is neglected in Fig. 2, but it would be expected that including this should shift the peak by the value of the magnetic quantum number. This behaviour was in fact observed in [51]. The analytical formalism we have derived is able to give the core properties of the OAM for strong field ionisation; however, in an actual experiment the laser will have a pulse envelope and thus a distribution of frequencies. In the next section, using numerical computations, we explore the effect this has on the OAM distribution.

Fig. 2
figure 2

The result of Eq. (38) is plotted for 2 (a), 4 (b) and 8 (c) ionisation events (given by N). The laser field intensity employed is \(I=2\times 10^{14}\) W/cm\(^2\) and a wavelength \(\lambda =800\) nm, corresponding to \(\mathrm {U}_\mathrm {p}=0.44\) a.u. and \(\omega =0.057\) a.u. for the ponderomotive energy and angular frequency, respectively. The distributions are plotted versus energy in order to easily see condition (30), which is marked by the vertical dashed lines on each figure and has constant spacing in energy. The OAM l values 20–37 are included as shown in the figure legend. The bound state prefactor has been neglected in the calculation

Full size image

4 Numerical computations

In this section we will examine the effect of a pulse envelope on the OAM distributions. We will employ 2-cycle, 4-cycle and 8-cycle laser pulses as in the previous section. Employing a \(\sin ^2\) envelope we set the vector potential to be

$$\begin{aligned} \mathbf {A}(t)=-\sqrt{2\mathrm {U}_\mathrm {p}}\sin \left( \frac{\omega t}{2 N}\right) ^2 \left( \cos (\omega t)\varvec{e}_{x} +\sin (\omega t)\varvec{e}_{y}\right) . \end{aligned}$$
(43)

The numerical computation of OAM distributions closely follows the methodology of the previous section, using the saddle point approximation to compute the plane-wave transition amplitude [as in Eq. (31)] and then numerically computing the \(\phi \) integral from Eq. (19). The \(\phi \) integral can be computed very efficiently using the FFT algorithm. The saddle point solutions are computed using Eq. (32). The solutions can also be found very efficiently by transforming the equation to a polynomial of order \(2N+2\) and finding the roots, as outlined in [68]. This means there will be \(N+1\) valid solutions, see Fig. 3.

Fig. 3
figure 3

The laser pulse (first row), imaginary parts of ionisation times (second row) and momentum distributions (final row). In the first row the 2-cycle, 4-cycle and 8-cycle laser fields are plotted in a–c, respectively. The x [y] component of the laser field is plotted by the blue [orange] line. The real part of the times of ionisation are given by the vertical dashed lines and coloured spots are used to mark where this intersects with the x component of the laser field. These times are found by solving Eq. (32). The imaginary parts of the solutions are plotted in the middle row, d–f versus the perpendicular momentum coordinate \(p_\perp \), the colours of each line corresponds to the colour of the spots in the panels above, the real parts are also explicitly given (in units of laser cycles) on the right-hand side of each panel. The minima of the imaginary part of the times of ionisation are marked by dotted lines. The photoelectron momentum distributions are plotted along the bottom row (h–i), for the momentum components in xy-plane and \(p_z=0.1\) a.u., computed using Eq. (31) with saddle point from Eq. (32) for each laser pulse considered

Full size image

In Fig. 3 the 2-cycle, 4-cycle and 8-cycle results are shown. In the first column the laser field is plotted with the real parts of the times of ionisation indicated by the vertical dashed lines. In the second column the imaginary parts of the ionisation times are plotted versus the radial momentum coordinate. The minima of these times (see the dotted line in Fig. 4) is a predictor for where the momentum distribution will be maximum. The momentum distributions are plotted in the final row for each laser pulse. As suggested by the imaginary part of the times, the peak of these distributions is away from the centre forming doughnut shapes. Interference between different paths can be seen as faint circular fringes. The following references [68,69,70] can give some further insight into these solutions and the method. Such distributions bear some similarity with the attoclock [71,72,73,74,75,76,77,78], where an elliptical nearly circular field is used to relate the electron emission angle to the tunnelling time. In this work we are relating the photoelectron OAM to the emission energy, which are the conjugate variables of angle and time, respectively.

Fig. 4
figure 4

The OAM distributions are shown for a 2-, 4- and 8- \(\sin ^2\) laser pulse (a–c), respectively. The peak intensity, wavelength and ionisation potential is the same as that used in Fig. 2. The bound state prefactors are neglected. The dot-dashed lines correspond to the condition Eq. (30) with the following effective frequencies used \(3/4 \omega \), \(15/16 \omega \) and \(63/54 \omega \), for (a–c), respectively. The condition has been slightly shifted in each case so that the central peak lines up with the corresponding l value. The OAM values utilised are from \(l=20\) to \(l=40\)

Full size image

In Fig. 4 the OAM distributions for the three laser pulses used in Fig. 3 are shown. The distribution of frequencies in the \(\sin ^2\) pulse means that we should not expect Eq. (30) to hold. In panel (a), the OAM distribution is plotted for a 2-cycle pulse. There are still specific OAM peaks in each energy region; however now they significantly overlap. Furthermore, the central peak at around 0.6 a.u. corresponds to \(l=31\) instead of \(l=27\) as in the monochromatic case. This is because the spacing between the peaks is reduced. The dot-dashed lines in Fig. 4a use a spacing of \(\frac{3}{4}\omega \). This closely matches most of the peaks, but for increasing l above \(l=31\), the spacing drifts to higher values, while for decreasing l below that of \(l=31\) the spacing drifts to lower values. For the longer pulses, 4 and 8 cycles, the distributions move closer to the monochromatic case, with the \(l=27\) peak moving towards its previous position of 0.6 a.u. and the spacing between peaks getting closer to \(\omega \). However, the spacing is still below this, with spacing of \(\frac{15}{16}\omega \) and \(\frac{63}{64}\omega \) being used for the dot-dashed lines for 4 and 8 cycles, respectively. Note that in all cases as well as altering the frequency in condition Eq. (30) a small shift was required to align the dot dashed lines to the correct central peak (reducing for longer pulses), see the caption of Fig. 4 for more details.

A clear pattern is visible in the spacing of the OAM distribution of 2-, 4- and 8-cycle pulses and it is possible to analytically derive this spacing dependence. In the monochromatic case the spacing between OAM peaks is given by \(\nu (t')=\omega t'\), which leads to an \(\omega \) spacing in energy in the condition Eq. (30). There is not such a simple relation for \(\nu (t')\) when employing a \(\sin ^2\) pulse, but through a Taylor expansion we can derive such an expression. A first-order Taylor expansion about the peak of the pulse (\(t=\pi N/\omega )\) gives \(\nu (t')=\omega ^*_{N} t'\) in terms of an effective frequency \(\omega ^*_{N}\), dependent on the number of laser cycles. Note a constant term has been discarded as this only contributes an overall phase and does not affect the OAM peak separation. The effective frequency can be calculated to be

$$\begin{aligned} \omega ^*_{N} = \nu '(\pi N/\omega )= {\left\{ \begin{array}{ll} \left( 1-\frac{1}{N^2} \right) \omega , &{} \quad \text {if}\,N\,\text {is even}\\ \omega &{} \quad \text {otherwise} \end{array}\right. }, \end{aligned}$$
(44)

which matches the peak spacing seen in 2-cycle, 4-cycle and 8-cycle pulse in Fig. 4.

Fig. 5
figure 5

The same as Fig. 4 except 3-, 5- and 9- cycle \(\sin ^2\) laser pulses have been employed. The dot-dashed lines use the carrier frequency for the spacing and are again slightly shifted so that Eq. (30) to match the central peak

Full size image

It is interesting to note that, if the number of cycles N is odd, the peaks remain separated at the carrier frequency. This is exemplified in Fig. 5, which shows the OAM distributions for 3-cycle, 5-cycle and 9-cycle laser pulses. The dot-dashed lines in this figure are all separated by the carrier frequency and the corresponding peaks in the distributions line up with this well. It is possible to see some drift away from the central peak in (a), as the spacing is slightly smaller than the carrier frequency. These shifts are simply because the condition, although accurate, it is no longer exact and thus there is a increasing drift for higher l values.

Beyond the spacing between peaks there are further differences between Figs. 4 and 5. Namely, in Fig. 4 considerable secondary peaks can be observed, which for in cases where the OAM peaks are low [see \(l=23\) in (c)] the secondary peaks can be nearly as high as the primary peak [or even higher in extreme cases e.g. \(l=22\) in (c)]. However, the secondary peaks always occur in the regions where there is a more dominant OAM peak corresponding to another l value. In contrast, for an odd number of cycles, Fig. 5, the secondary peaks are considerably lower, not playing much of a role at all. These features suggest that, if a clean well separated (in energy) OAM distribution is required, it would be much better to utilise an odd number of laser cycles in the laser pulse envelope. Furthermore, in the even-cycle case the peaks are slightly asymmetric, with a longer tail on one side than the other, while for an odd number of cycles the peaks are symmetric. These differences suggest one may be able use the OAM of the photoelectrons to detect if the laser field had a closer to odd or even number of cycles. However, additional effects such as the intensity variation over the laser focal volume (see Ref. [79] for an example) would need to be considered first, and as previously stated, the measurement of the OAM of photoelectrons in strong field experiment is an open problem [51, 52].

5 Conclusions

In this study we have developed a new version of the strong field approximation (SFA) to explore the orbital angular momentum of photoelectrons undergoing strong field ionisation. Employing this model we are able to derive analytic conditions relating to selection rules/conservation laws about the system. In the case of a linear field we demonstrate that the photoelectron OAM is determined by the magnetic quantum number of the initial bound state, while for a circular field we show that a range of OAMs are possible, which will occur in well-defined energy regions. We derive an analytic condition that relates their position to the ATI peaks. Computations using a \(\mathrm{sin}^2\) laser pulse demonstrate that this condition continues to be accurate even for very short pulses but the separation of the peaks is now described by an analytically derived effective frequency as opposed to the carrier frequency of the laser field. We find that employing an odd or even number of laser cycles has a marked effect on the OAM distributions, leading to two classes of effective frequency for even and odd pulses.

The present work is an interesting example of how the ellipticity and time profile of the field add dynamic shifts \(\nu (t')\) to angular variables intrinsic to the target, such as the OAM of the electron. This effect bears some similarity with previous work on the angular properties of the photoelectrons. In one example, the tunnelling angle, derived in [80], was used to show the preferential tunnel ionisation of ‘counter-rotating’ electrons with circularly polarised fields. An electron’s angle of return also leaves imprints in HHG spectra [81, 82]. These angular shifts were incorporated in a purely structural two-centre interference condition for HHG in diatomic molecules [82], and could be made visible by exploiting macroscopic propagation [83]. A key difference is that these articles computed angles related to the velocity of the electron, whereas in the present work \(\nu (t')\) is related only to the rotation of the field and represents the interaction of the laser field with the OAM of the electron. It could be said to be the ‘stirring’ action of the field upon the electron at the moment of ionisation.

The method proposed here is general and could be extended to initial states with arbitrary angular momentum and other types of tailored fields. However, an open question is the role of the residual binding potential. In [51] good agreement was found between the SFA and the TDSE solvers Qprop and RMT but the longer wavelengths used in this study may lead to more significant Coulomb effects. This question has been addressed in the attoclock setup, which (as previously stated) deals with the conjugate variables of the emission angle and corresponding ionisation time. In this setting the Coulomb potential has been shown to shift the photoelectron emission angle (see e.g. [76]), which may have a profound effect on the interpretation of the tunnelling time (see [78, 84] for reviews on the attoclock). Thus, it should be expected that the Coulomb potential will also shift the OAM of the photoelectron. This may be addressed by incorporating the OAM into SFA-like models which account for the binding potential, such as the Coulomb quantum-orbit strong field approximation [56, 66, 85] as well as performing comparison with TDSE solvers.

The OAM decomposition of the final photoelectron wavepacket presented here reveals the angular momentum content along a chosen direction, associated with the OAM that may be carried by freely propagating particles. The electron wavepacket does not form a beam. Therefore, the wavepacket cannot be envisaged as single vortex beam propagating along the OAM axis but rather a superposition of many different Bessel beam vortex states with different OAMs. This can be exploited either in strong field measurement [51, 52], where such distributions over OAM may reveal unique properties about the target system, or alternatively could be used to produce a custom electron vortex beam [15]. Such as scheme could exploit the relationship between the ATI peaks and the OAM to filter photoelectrons from a specific ATI peak in order to select a single OAM.

Another very important question is: How about other types of fields, for more complex vortices? The response of matter to structured laser fields is a central subject of intense laser-matter physics, and is inevitably related to the OAM of matter. Recently, a novel pump-probe scheme incorporating OAM was theoretically demonstrated using a vortex IR beam and XUV pulse allowing for time-resolved photoionisation [86], while DeNinno [87] demonstrated experimentally that a free electron matter wave, produced by an XUV pulse, is sensitive to a vortex IR beam. Similarly, angle-resolved attosecond streaking of twisted attosecond pulses has been recently proposed [88]. In a condensed matter context, it was shown recently that THz laser pulses with circular polarisation induce transient Chern insulator in graphene [89]. Similarly, linearly polarised pulses with OAM induce non-uniform Chern insulators [90]. All these examples clearly illustrates that we are entering an era of OAM and more complex textures in laser-matter physics.