Electron acceleration from transparent targets irradiated by ultra-intense helical laser beams

The concept of electron acceleration by a laser beam in vacuum is attractive due to its seeming simplicity, but its implementation has been elusive, as it requires efficient electron injection into the beam and a mechanism for counteracting transverse expulsion. Electron injection during laser reflection off a plasma mirror is a promising mechanism, but it is sensitive to the plasma density gradient that is hard to control. We get around this sensitivity by utilizing volumetric injection that takes place when a helical laser beam traverses a low-density target. The electron retention is achieved by choosing the helicity, such that the transverse field profiles are hollow while the longitudinal fields are peaked on central axis. We demonstrate using three-dimensional simulations that a 3 PW helical laser can generate a 50 pC low-divergence electron beam with a maximum energy of 1.5 GeV. The unique features of the beam are short acceleration distance (∼100 μm), compact transverse size, high areal density, and electron bunching (∼100 as bunch duration). Plasma mirrors have become the preferred method for electron injection to laser-based accelerators, but the optimal configuration is difficult to achieve. Here, an alternative injection method employing a low-density foam target and a helical laser pulse is investigated numerically.

T he interest in helical light structures was triggered by a groundbreaking study showing that they carry orbital angular momentum 1 . There has been an ongoing effort to utilize this property for various applications 2 . In recent years, the attention has shifted towards higher intensity laser pulses with helical wave fronts. The shift in interest can be credited to the realization that the helical wave front structures can be achieved by reflecting a conventional high-intensity laser beam off a helical step-like mirror [3][4][5][6] . Given the high conversion efficiency afforded by these mirrors, there is potential for these beams to be created on the latest generation of extreme-intensity (~10 23 W cm −2 ) laser facilities 7,8 .
High-intensity helical beams represent drivers with features that are qualitatively different from those of conventional laser beams. Manifestations of these features include the generation of large magnetic fields inside a plasma [9][10][11] , production of helical electron beams in laser-wakefield accelerators 12,13 , and improved laser-driven ion acceleration 14 . The helicity or the wave front twist introduces an extra 'control knob' for adjusting the field topology. For example, the radial profile of the transverse electric and magnetic field can be made hollow by twisting their wave fronts. For a properly chosen twist, the helical laser beams can have strong longitudinal electric and magnetic fields that peak on the axis where there are essentially no strong transverse fields. Such a unique field structure offers an exciting possibility of generating collimated beams of ultra-relativistic electrons via direct acceleration of electrons by the longitudinal laser electric field.
The idea of using intense laser pulses to directly accelerate electrons in a vacuum to ultra-relativistic energies 15,16 is as old as the capability to generate such pulses. However, successful experimental implementation of the concept has been elusive, as it requires a method for efficient electron injection into the laser beam and a mechanism for counteracting transverse electron expulsion that is typical for conventional laser beams. It has recently been shown that no expulsion occurs if, rather than using a conventional laser beam, one uses a helical laser beam with dominant longitudinal electric and magnetic fields in the nearaxis region 17 . In this case, the electrons obtain their energy from the longitudinal electric field. The longitudinal magnetic field is essential for electron retention. Radially polarized laser beams can also have a dominant on-axis longitudinal electric field 18,19 , but they have no on-axis magnetic field. As a result, electron acceleration by these beams is rather sensitive to perturbative effects 19 .
In search of a solution to the electron injection problem, multiple approaches have been examined and laser reflection off a plasma mirror has emerged as the preferred method [17][18][19] . However, this method is sensitive to the gradient of the plasma density. Achieving a necessary steep density gradient at high laser intensities is extremely difficult. It is also challenging to reliably measure the electron density profile all the way to the reflection point since the reflection at relativistic intensities occurs inside a classically opaque region. Another difficulty arises due to the constraint that experiments at ultra-high intensities must be performed using oblique incidence. It has been shown experimentally for radially polarized beams that the laser reflection in the case of oblique incidence negatively impacts electron injection and their subsequent acceleration 19 .
In this paper, we consider a qualitatively different approach, offered by the recent progress in target fabrication, where the injection takes place during volumetric interaction between the laser and a solid target. It is now possible to produce inorganic foam targets (e.g., SiO 2 foam targets 20,21 with the pore size of just a few nm) whose electron density, after the target is fully ionized and homogenized, is close to the classical cutoff density for an optical laser. For a laser with an 800 nm wavelength, the cutoff density (also often referred to as the critical density) is n c ≈ 1.7 × 10 21 cm −3 . This value is much higher than the electron density in an expanding gas jet and much lower than the density in a solid metal target. The advantage of near-critical targets, as those made of inorganic foams 22 , is that they become transparent when irradiated by a high-intensity laser, so they can transmit the laser pulse and enable a volumetric interaction of the pulse with a dense electron population. In the context of the electron acceleration problem, the benefit of the volumetric interaction compared to the reflection is that the electron injection takes place at an electron density that is known relatively well. Additionally, the density is set by the target fabrication process rather than by the laser prepulse, so it is possible to choose the desired density 21 .
Using three-dimensional particle-in-cell (3D PIC) simulations we show that both the electron injection and the electron retention problems can be successfully solved by employing a~10-μmthick near-critical target irradiated by a helical laser beam. We consider a 3 PW, 29 fs laser beam focused to a 3 μm spot, but our results are rather generic and not specific to these exact parameters. The electrons are "injected" into the laser beam while it is being transmitted by the near-critical target. We find that a proper choice of the target electron density n e is important for preserving the topology of the laser beam. At maxðn e Þ ¼ 0:5n c , the considered laser beam retains its structure without generating high-amplitude higher-order radial modes whose presence is detrimental to the generation of an ultra-relativistic collimated electron beam. The simulation for the optimal regime shows that, upon exiting the target, the laser beam carries dense bunches of electrons, which solves the injection problem. As shown in Fig. 1a, b, the bunches remain close to the axis of the laser beam due to the hollow profiles of the transverse fields, which solves the retention problem. The electrons are accelerated by the longitudinal electric field that has a peak on the axis. The electron energies seen in Fig. 1c reach 1.5 GeV, which is similar to the maximum energy achieved using a plasma mirror and a 6 PW laser pulse 17 . The divergence of the beam, whose charge is about 50 pC, is less than 3 ∘ . Such a low divergence is the advantage of the acceleration by the longitudinal, rather than the transverse, electric field. The presented electron injection and acceleration mechanisms are qualitatively different from those employed by more established schemes that can achieve similar electron energy and charge, e.g., laser wakefield acceleration 23,24 . The differences lead to several unique features that include short acceleration distance (~100 μm); high areal density (~10 16 cm −2 ); the compact transverse size of the beam (≲1 μm); and bunching of energetic electrons, with the bunches being~100 as in duration.

Results
We consider a setup where a thin (~10 μm) near-critical target is irradiated by an ultra-high-intensity laser beam and becomes very transparent due to the effect of the relativistically induced transparency. The purpose of the target is to "inject" electrons into the transmitted laser beam for their subsequent acceleration in a vacuum. In our setup, the majority of the energy gain occurs during the vacuum propagation of the transmitted laser. In order to make this process effective, one must use a beam structure that can accelerate electrons without the transverse expulsion that is typical for conventional laser beams whose transverse fields peak on the central axis. Here we focus specifically on helical beams and, in the next two subsections, we identify the most suitable helical mode and make analytical predictions for the electron energy gain. The desired mode can be created (even at a high intensity) from a conventional laser beam using recently developed techniques involving reflective optics. However, the considered target must preserve the helicity of the beam during its transmission without generating high-amplitude higher-order radial modes. Our 3D PIC simulations for different target densities explore the changes induced during the beam transmission. We use these results to identify the optimal regime. This section concludes by considering a case of oblique incidence for the optimal density aimed to show that the considered setup can be used in those laser systems that do not permit experiments at normal incidence.
Vacuum field structure of a helical laser beam. As shown in Fig. 1a, we consider a laser beam propagating in the positive direction along the x-axis. In the case of a beam with helical wave fronts, it is convenient to represent the spatial structure of each field component as a superposition of Laguerre-Gaussian modes. The modes describing vacuum propagation are well-known, so we provide the corresponding expressions without derivation. The purpose is to identify the key features and set the stage for the discussion of electron acceleration.
The polarization of a laser beam determines the orientation of its fields, whereas the helicity or twist index l determines the topology of the wave fronts. We illustrate the difference by considering a linearly-polarized beam with wavelength λ 0 and beam width w 0 that has E y and E x electric field components and B z and B x magnetic field components. It is convenient to use cylindrical coordinates (x, r, ϕ) to describe the spatial structure of each component, where r is the distance and ϕ is the polar angle in the (y, z) plane. The transverse electric field E y of a mode with a twist index l and radial index p is given by E y ðe x;e r; ϕ; tÞ ¼ E 0 ψ p;l gðξÞ expðiξÞ; ð1Þ where E 0 is a constant that sets the field amplitude, C p,l is a normalization coefficient, g(ξ) is the temporal envelope, L jlj p is the generalized Laguerre polynomial. Additionally, we introduced the beam divergence angle θ d = w 0 /x R , the Rayleigh length x R ¼ πw 2 0 =λ 0 , and the laser frequency ω = 2πc/λ 0 . The tilde indicates the following normalization: e r r=w 0 : ð6Þ We set C l;p ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2p!=πðp þ jljÞ! p , so that R 2π 0 dϕ R 1 0 ψ p;l ψ Ã p;l e rde r ¼ 1. The transverse magnetic field is simply B z = E y . It follows from Eqs. (1) and (2) that, even though the polarization of the modes with different l is the same, the surface of the constant phase changes its shape with the change of l. Specifically, it becomes helical for l ≠ 0, whereas the l = 0 case corresponds to a conventional laser beam whose wave fronts have no twist.
The helicity profoundly alters the radial dependence of the laser fields. For simplicity, we examine the field structure in the vicinity of the focal plane (e x ¼ 0). The transverse fields of a conventional beam (l = 0) peak on the central axis. In contrast to that, E y and B z of a helical beam vanish on the central axis, with E y ¼ B z / e r jlj at e r ( 1. It follows from Eqs. (1) and (4) that for l ≠ 0 the surface of constant phase for the transverse electric and magnetic fields is a helix. The sign of l determines the direction of the twist and the value of |l| determines the number of full rotations/twists per wavelength: Δx/λ 0 = − lΔϕ/2π.
The changes in the topology of the transverse fields lead to similarly profound changes in the topology of the longitudinal electric and magnetic fields. The changes result from the fact that the electric and magnetic fields must be divergence-free, with (∇ ⋅ E) = 0 and (∇ ⋅ B) = 0. We are particularly interested in the field structure close to the axis where E y and B z vanish for helical beams. It follows from (∇ ⋅ E) = 0 that for the considered linearly-polarized wave we have For a helical beam (i.e., l ≠ 0), we have ∂E y =∂e r % jljE y =e r close to the central axis. We take this into account to find from Eq. (7) that at e r ( 1 the longitudinal electric field is approximately given by where the superscript on the left-hand side represents the sign of l. The longitudinal magnetic field B x has the same radial dependence. The modes with |l| = 1 are evidently unique because their longitudinal fields peak on the central axis. The longitudinal fields of these modes dominate in the near-axis region, so these modes are well-suited for longitudinal electron acceleration. These results can be generalized to the case of a circularlypolarized laser beam. We construct such a beam out of two linearly-polarized beams with a phase offset. Its transverse electric field is set by the relation E z = iσE y , where E y is given by Eq. (1). Here σ = 1 corresponds to a right-circularly-polarized wave and σ = −1 corresponds to a left-circularly-polarized wave. Only for |l| = 1 the longitudinal rather than transverse fields peak on axis, but the two circular polarizations (right and left) are not equivalent. In the case of σ = −l, the longitudinal fields reach their highest amplitude at e r ! 0. On the other hand, in the case of σ = l, the longitudinal fields vanishes on the central axis. Figure 1 compares the topology of the transverse and longitudinal electric fields in the circularly-polarized beam with σ = −l and l = −1. Figure 1d, e show that the transverse field for this beam indeed vanishes on the central axis. In the near-axis region, the longitudinal electric field dominates, as seen in Fig. 1f. The structure of the longitudinal field is shown in Fig. 1g, h.
Close to the axis, E x and B x are axis-symmetric for both linearly and circularly-polarized pulses. The linearly-polarized beam loses this symmetry away from the axis. In contrast to that, the fields of the circularly-polarized beam with σ = −l retain their symmetry [e.g., see Fig. 1h]. Motivated by this observation, we focus in the remainder of this paper on the more symmetric circularly-polarized beam with σ = −l= 1. Its longitudinal electric field is iθ d E y e iϕ jlj=e r À e rf Â Ã for p ¼ 0: The longitudinal magnetic field has the same radial dependence, but it has an additional phase shift: B x = iE x . We have intentionally retained higher mode numbers corresponding to p ≥ 1 in our analysis. Even if the irradiating laser beam contains only a single radial mode (e.g., p = 0), the transmitted beam is likely to consist of various modes with different amplitudes (this assertion is later confirmed via mode decomposition). The mode coupling during laser propagation through the thin targets occurs due to laser focusing by the plasma. By definition, the radial mode number p changes the radial field structure of E x away from the axis. However, the onaxis field changes with p as well. In order to clearly identify the corresponding changes, we take the limit of e r ! 0 in Eq. (9). Note that the expressions in square brackets reduce to jlj=e r for p ≥ 1 and for p = 0. We then have is the peak value of the longitudinal electric field. The explicit dependence of the phase of the on-axis field on the mode number via the (1 + p) multiplier deserves particular attention. We define the phase velocity as v ph = dx/dt for the constant phase given by the expression in the square brackets in Eq. (10). We take the derivative of this expression with respect to t and assume that θ d is small to find the following approximate expression for the relative superluminosity: The key feature here is that the relative degree of superluminosity is given by Eq. (12) increases with the radial mode number. Even though the phase velocity is very close to the speed of light, the explicit dependence on p can have a significant impact on the dynamics of ultra-relativistic electrons with Estimates for electron acceleration by the longitudinal electric field of a helical beam. In order to obtain analytical scalings for the energy and momentum gain, we consider an electron that is moving along the central axis of a circularly-polarized helical beam with σ = −l = 1. The transverse fields of this beam vanish on the axis, so such a simple model is a reasonable approximation for the dynamics of electrons that start their motion close to the axis. Even though most of the energy gained by electrons in our setup takes place in a vacuum, we anticipate that the electrons undergo some initial acceleration while moving with the laser through the target. We mimic this by assuming that electron starts its vacuum acceleration with an initial relativistic long- Under our assumptions, the change in electron momentum results from the electron interaction with E x ðe r ¼ 0Þ whose explicit expression is given by Eq. (10): where t 0 is the time when the electrons starts its vacuum acceleration with the initial momentum p 0 . In order to compute this integral, one needs to know how the phase of E x evolves along the electron trajectory. Even though a general dependence is complicated, a simplified expression can be derived for a regime for a forward-moving electron with γ ¼ 1= . The irradiating laser beam in the simulation whose results are shown in Fig. 1 has w 0 = 3 μm, λ 0 = 0.8 μm, and p = 0. This means that the requirement c − v x ≪ v ph − c is satisfied after the electron γfactor (at e x ≲ 1) exceeds γ Ã % 12. This is not a particularly constraining condition, considering that the peak value of γ exceeds 3000. The condition eventually breaks down at large e x as v ph decreases and approaches c, with the corresponding value of e x increasing with γ. In the considered example, this occurs at a very large distance of e x % 17 even for a relatively modest γ ≈ 100. At such a large distance, the laser amplitude is very low due to the beam divergence, so no energy gain takes place and the breakdown of the considered condition is inconsequential.
In the regime where c − v x ≪ v ph − c, the change in phase is predominantly determined by the superluminosity of the wave fronts. We can thus set x = x 0 + c(t − t 0 ) and eliminate the explicit time dependence of the phase. For simplicity, we consider an electron that undergoes its acceleration in a region close to the peak of the field envelope, so that g(ξ) ≈ 1. The resulting expression for the longitudinal on-axis electric field is where and e and m e are the electron charge and mass. We have introduced Φ 0 , which is the phase of E x at e x ¼ e x 0 at t = t 0 , to explicitly take into account the initial conditions. In our analysis, Φ 0 is an input parameter that can be interpreted as the 'injection phase' at the start of the vacuum acceleration. In order to calculate the integral in Eq. (14) we change the integration variable from t 0 to e x 0 by taking into account that de x 0 =dt 0 % c=x R . The resulting expression for the change of momentum is Based on the obtained expression, there are two distinct regimes of electron acceleration: acceleration by modes with an even radial index p and acceleration by modes with an odd radial index p. In order to make the difference more evident, we consider the terminal change in momentum. The corresponding expression is obtained by taking the limit of e x ! 1. We further simplify the expression by assuming that e x 0 ( 1, so that tan À1 ðe x 0 Þ % e x 0 . In the case of an even p value, we obtain In the case of an odd p value, we obtain The difference in the prefactor that is independent of the injections phase Φ 0 means that the momentum gain for odd values of p can be significantly suppressed if ðp þ 1Þe x 0 ( 1. The p = 1 mode is the one that is most impacted for e x 0 ( 1. The qualitative difference in acceleration between different p values can be understood by examining the terminal phase slip. According to (17), the phase slip is ΔΦ ¼ 2ðp þ 1Þ½tan À1 ðe x 0 Þ À tan À1 ðe xÞ. For an electron that starts its acceleration at e x 0 % 0, the phase slip at e x ! 1 is ΔΦ = − π(p + 1). The dependence on p here is a result of the dependence of v ph on p [see Eq. (12)]. The key point is that the electron is slipping with respect to E x faster at higher p. Since ΔΦ = − π for p = 0, then a properly injected electron can remain in the accelerating phase of the beam almost the entire time. In contrast to that, we have ΔΦ = −2π for p = 1. The electron then experiences not only acceleration but also deceleration by E x . The deceleration leads to the reduction of ðΔp x Þ term given by Eq. (19). At p = 2, ΔΦ = −3π and the electron can experience the accelerating phase twice. As a result, the momentum gain increases compared to the case with p = 1. However, it is reduced compared to the p = 0 case because the acceleration alternates with deceleration.
The key conclusion from our results is that, for given a 0 and θ d , the biggest momentum increase and thus the biggest energy gain should be expected for a mode with p = 0. Excitation of the modes with higher p values during the beam propagation through the target is likely to be detrimental. The plasma can also focus the beam, reducing w 0 . The size of the focal spot w 0 however has no impact on the terminal momentum gain at fixed power 17 . At fixed power, a 0 scales as 1/w 0 . On the other hand, the divergence angle θ d also scales as 1/w 0 . As a consequence, the ratio of a 0 /θ d that enters the expressions for ðΔp x Þ term remains unaffected by changes in w 0 , provided that the beam power remains the same.
3D PIC simulation results. In order to perform a self-consistent analysis of the electron acceleration in the considered setup, we have performed a 3D PIC simulation with a moving window where a uniform plasma slab is irradiated by a normally incident 3 PW, 29.4 fs helical laser beam that is focused to a 3 μm spot to achieve a 0 ≈ 89. Detailed simulation and target parameters are provided in the Methods section. The incident beam is chosen to be circularly-polarized with l = −1, σ = 1, and p = 0 because this mode has been identified as the optimal mode for the vacuum acceleration of electrons. The laser focal plane (in the absence of the target) is located at x = 0 μm. The target is 12 μm-thick with a central plane located at x = 0 μm. The target is uniform in the central region that is 8 μm-thick. There are 2 μm density ramps on both sides to mimic target pre-expansion that likely takes place in experiments with ultra-high intensity due to the presence of a prepulse. The electron density in the target is subcritical, n max ¼ 0:5n c , so the target is transparent to the laser pulse. Figure 1a shows the topology of the longitudinal laser electric field E x and electron density in the beam transmitted by the target at t = 88.3 fs. The time t = 0 fs is defined as the time when the peak of the laser envelope passes through x = 0 in the absence of the target. In the incident beam, E x peaks on the axis of the beam, since l = −1 and σ = 1. The target evidently preserves this key feature, with E x remaining peaked on the axis of the transmitted beam. The surface plot of the electron density past the initial surface of the target (x > 6 μm) confirms that the laser propagation through the transparent target serves as an effective mechanism for electron injection. The peak electron density in the transmitted beam exceeds 10% of the initial electron density in the target. There are two key features that are evident from the plot of n e : the electrons are confined in the transverse direction and they are bunched longitudinally with a period roughly equal to λ 0 . The confinement results from a combination of vanishing transverse electric and magnetic fields near the central axis and a peaked longitudinal magnetic field on the central axis. The electron bunching is a manifestation of the electron acceleration by the longitudinal laser electric field. Similar bunching can be observed for test electrons that are initially uniformly distributed along one wavelength by integrating electron equations of motion with E x given by Eq. (15). Figure 1b, c show several representative electron trajectories and the time evolution of the electron spectrum for all injected electrons. Rather than moving only along the x-axis, the electrons in Fig. 1b move at an angle of θ ≈ 3.3 × 10 −3 to the central axis during the acceleration process. In this case, the longitudinal velocity is determined mainly by the angle θ and not by the difference between v and c, so that ðc À v x Þ=c % 1 À cos θ % θ 2 =2. On the other hand, the laser divergence angle is at least θ d = 8.5 × 10 −2 (the transmitted beam is narrower, which means its divergence angle is greater than this value). It follows from Eq. (12) that, for e x ≲ 1 and p = 0, we have ðv ph À cÞ=c ≲ θ 2 d =2. We then conclude that c − v x ≪ v ph − c for the considered particle trajectories. Therefore, the change in phase during the acceleration is predominantly determined by the superluminosity of the wave fronts and the expressions derived for the electron momentum gain are applicable to these trajectories in Fig. 1b.
In order to estimate the maximum possible momentum gain, we assume that the transmitted beam is only a p = 0 mode. We have already shown that the beam narrowing or widening is inconsequential in terms of the momentum gain, so we take the original values of θ d = 8.5 × 10 −2 and a 0 = 88.6. It follows from Eq. (18) that This corresponds to a maximum energy gain of 1.2 GeV during the acceleration in vacuum by the longitudinal electric field of the helical beam with p = 0. The electrons that start at the rear side of the target, e.g. x 0 = 5 μm, can achieve an energy gain very close to this value because the target is relatively thin compared to the Rayleigh length x R . Indeed, we have cos e x 0 À Á cos Φ 0 þ e x 0 À Á % 0:99 for Φ 0 ≈ −0.047π and e x 0 ¼ x 0 =x R % 0:14, where x 0 = 5 μm and x R ≈ 35 μm.
According to the time evolution of the electron spectrum shown in Fig. 1c, the maximum terminal electron energy in the 3D PIC simulation is however about 25% or 300 MeV higher than what is predicted by Eq. (20). We have performed particle tracking that shows that the electrons leaving the target have appreciable energy, with a peak kinetic energy of~300 MeV at x ≈ 5 μm. Using this value as an initial condition, we can then reconcile the result of the PIC simulation with the prediction given by Eq. (20). Specifically, ε term where p 0 is the initial momentum corresponding to ε k = 300 MeV and ðΔp x Þ term is the contribution from the vacuum acceleration for Φ 0 ≈ −0.047π and x 0 = 5 μm. The conclusion then is that the acceleration within the target, which is likely to differ from that in a vacuum due to the presence of the plasma electric and magnetic fields, provides an additional mechanism for the energy gain in the considered setup.
Target density scan. We have performed two additional 3D PIC simulations to examine the impact of the target density on the discussed electron acceleration process. In both simulations, the electron density in the target is higher than the electron density in the original simulation (n max ¼ 0:5n c ). These two densities are n max ¼ 2:5n c and n max ¼ n c . The results of the two additional runs are shown in Fig. 2, where we also show the results of the original simulation to facilitate the comparison between different target densities. Figure 2a-f show the electron density and the topology of the longitudinal electric field shortly after the laser pulse transmission in all three simulations (t = 48.3 fs). We find from Fig. 2a-c that the divergence of the electron "beam" generated by the laser pulse outside of the target (x > 6 μm) reduces as we reduce the target density. The fact that the electrons in Fig. 2c have traveled further along the x-axis than in Fig. 2a is a result of a higher group velocity in the lower density target. We also observe noticeable changes to the structure of the longitudinal electric field E x in the transmitted laser beam. As we reduce the target density, we find that the radial dependence in Fig. 2d-f changes, which suggests that higher radial modes (the original beam is a p = 0 mode) are generated during the laser beam transmission by the target. Figure 2g-i show the electron density profiles at t = 248.3 fs for all three targets. By this point in time, the electrons shown in the plots have experienced significant acceleration. We find that the already discussed increased divergence at higher target density at t = 48.3 fs translates into the higher divergence of the accelerated Fig. 2 Target density scan for a circularly-polarized helical beam. a-c Electron density n e normalized to the critical density n c for three different target densities (n max ¼ 2:5n c , 1.0n c , and 0.5n c ), shortly after laser pulse transmission (t = 48.3 fs). d-f Normalized amplitude of the longitudinal electric field, |e||E x |/ m e cω, for the same three target densities also at t = 48.3 fs. g-i Density of accelerated electrons for the same three target densities at t = 248.3 fs. j-l Time evolution of the energy distribution for forward-moving electrons selected using the same criterion as that used for Fig. 1c, with ε k = (γ − 1)m e c 2 . In panels (a-i), the quantities are shown in the (x, y)-plane at z = 0. Panels (a), (d), (g), and (j) correspond to a target with n max ¼ 2:5n c ; panels (b), (e), (h), and (k) correspond to a target with n max ¼ n c ; panels (c), (f), (i), and (l) correspond to a target with n max ¼ 0:5n c . electrons at t = 248.3 fs. In the case of the target with an initial density of n e = 0.5n c , the accelerated electrons remain within roughly 1 μm distance of the central axis, shown with the vertical dashed lines in Fig. 2i. In the case of the target with n max ¼ 2:5n c , there are electrons that are at least 6 μm away from the central axis, which means that the divergence angle for the electrons at this density is roughly six times higher than that in the simulation with n max ¼ 0:5n c .
The target density increase also impacts the energy gain by the electrons during their vacuum acceleration by the transmitted laser beam. Figure 2j -l show the time evolution of the energy spectra for the electrons that remain close to the central axis (|y| < 1 μm and |z| < 1 μm). As we increase the target density from n max ¼ 0:5n c to n max ¼ 1:0n c , the energy gain remains relatively unchanged. However, a further increase from n max ¼ 1:0n c to n max ¼ 2:5n c causes a visible reduction. We conjecture that the difference in the energy gain is, in part, caused by the excitation of higher-order radial modes, as these modes are less effective in terms of electron acceleration.
In order to find the changes in the structure of the transmitted beam in each simulation, we perform a mode decomposition detailed in the Methods. We perform our analysis for a snapshot of the magnetic field B z at t = 88.3 fs. In general, it is not possible to recover the mode structure from a single snapshot due to insufficient information regarding the propagation of the field. However, the inherent ambiguity can be resolved in our case by taking into account that there are no reflected modes behind the target. This information is not automatically encoded in the real field values outputted by the PIC code. We use the Hilbert transform to convert the real field B z provided by the PIC code into a complex field B H z , where the imaginary part is such that the field represents a beam propagating in the positive direction along the x-axis.
The field is given by the following mode decomposition: jejB H z ðe x;e r; ϕ; tÞ m e cω ¼ ∑ The transverse coordinate is normalized to w 0 0 rather than to w 0 that is associated with the irradiating beam. We determine the value of w 0 0 by minimizing the number of modes with considerable amplitude in the decomposition, as detailed in the Methods. To understand the need for minimization, consider the original beam. It is represented by a single-mode with l = −1 and p = 0 if w 0 0 ¼ w 0 . On the other hand, it is represented by multiple modes if w 0 0 ≠w 0 . The amplitudes of the additional modes increase with the increase of the mismatch in the transverse scale, jw 0 0 À w 0 j. The two representations are equivalent, but the decomposition that contains only a single mode is much easier to use when considering electron acceleration. Figure 3 shows jB H z j of the transmitted beam in the simulations with n max ¼ 2:5n c and n max ¼ 0:5n c at t = 88.3 fs [see Fig. 3a, d]. The discussed minimization procedure yields w 0 0 =w 0 % 0:8 for n max ¼ 2:5n c and w 0 0 =w 0 % 0:9 for n max ¼ 0:5n c . Additional details are provided in Methods. In both simulations, the modes with l = −1 and p = 0, which are the l and p values of the original beam, dominate the decomposition. The corresponding field profiles for these modes in the (x, y) plane are shown in Fig. 3b, e, where we use the following notation: jejB p;l =m e cω ¼ b p;l ðξÞψ p;l ðe x;e r; ϕÞ expðiξÞ. We find that max jb 0;À1 j % 92:43 for n max ¼ 2:5, whereas max jb 0;À1 j % 94:37 for n max ¼ 0:5.
For comparison, max jb 0;À1 j % 88:6 in the original beam. The amplitude increase is a result of the beam narrowing that occurs during the propagation inside the target. In both cases, the amplitudes of the modes with l ≠ −1 are negligible.
The narrowing of the beam inside the plasma generates modes with higher radial numbers. In both cases, the generated mode with the biggest amplitude is the p = 2 mode (see Methods for a plot of mode amplitudes). However, the mode amplitude increases with the increase of the target density. We find that max jb 2;À1 j % 30:68 for n max ¼ 2:5, whereas max jb 2;À1 j % 16:51 for n max ¼ 0:5. The structure of the modes with p = 2 and l = −1 is shown in Fig. 3c, f. The mode decomposition indicates that there is a phase shift with respect to the dominant mode, with Δξ ≈ −2.73 for n max ¼ 2:5 and Δξ ≈ −2.40 for n max ¼ 0:5 close to the peak of the envelope. The phase shift varies only gradually along each transmitted beam compared to the oscillations associated with λ 0 , so we can treat it as being constant for our estimates.
The described changes in the structure of the transmitted beams impact the acceleration of electrons in vacuum. Our primary focus is on the relative differences between the simulations with n max ¼ 2:5 and n max ¼ 0:5. We use Eq. (18) to find that for the acceleration by the main mode, i.e., the mode with p = 0 and l = −1, where we took into account that the amplitude of the transverse electric field is equal to the amplitude of the magnetic field and that θ d / 1=w 0 0 . The estimate assumes that the acceleration starts close to the focal plane and that the electrons are injected into the most favorable phase, which is Φ 0 = 0. The conclusion from Eq. (23) is that there is less power associated with the dominant mode at higher density, which contributes to the reduction in energy gain.
The impact of the p = 2 mode can also be assessed using Eq. (18). The analysis that we used to derive Eq. (18) neglects changes to the electron velocity during the acceleration. Therefore, the contributions from the modes with p = 2 and p = 0 can be added up while taking into account their phase shift. For the sake of simplicity, we again assume that the acceleration starts close to the focal plane and that the electrons are injected into the p = 0 mode with Φ 0 = 0. Using the parameters provided earlier for the p = 2 mode in both simulations, we find that The difference between the two cases has further increased with the inclusion of the p = 2 mode that is generated while the beam passes through the target. At n max ¼ 2:5, the amplitude of this mode is almost two times higher than at n max ¼ 0:5. This mode is decelerating the electron due to its phase shift, which explains why the difference between the two cases has increased.
Our analysis neglects higher-order harmonics generated during the considered laser-plasma interaction. The amplitude of the second harmonic in our simulations is roughly two orders of magnitude smaller than the amplitude of the main mode. We have performed two additional simulations with a higher spatial resolution (see the Methods section for details): one for the target with n e = 2.5n c and one for the target with n e = 0.5n c . In these simulations, the amplitude of the second harmonic is three orders of magnitude smaller than the amplitude of the main mode, which means that the amplitude of the second harmonic is lower than in the original simulations. We can therefore conclude that the harmonics generated in our setup are too weak to significantly influence electron acceleration. The mode decomposition shows higher radial modes with amplitudes that arẽ 30% of the zeroth mode, far larger than the contributions of the second harmonic. This is the reason why our focus is on the impact of the higher-order radial modes.
Our estimates explain the trend observed in the simulations and, more importantly, clarify the negative impact of higher radial modes on electron acceleration. Another feature of the modes with the high p values is that their transverse fields peak much closer to the axis than for the mode with p = 0. For example, the field in Fig. 3c peaks close to the vertical dashed line, whereas the field in Fig. 3b peaks at a distance from the central axis that is about three times greater. The same is the case for the transverse electric field. The transverse fields of the higher p modes close to the central axis deteriorate electron confinement and increase the rate of the phase slip by adding transverse momentum. These aspects are likely to further reduce the energy gain in the higher density case where the modes with p ≥ 2 have a much higher amplitude (see Methods).
Optimal regime at normal and oblique incidence. We now consider just the lowest density case, n max ¼ 0:5n c , that delivers the best performance with the goal of providing additional details regarding the electron acceleration. Figure 4 shows the case of normal incidence, where panel a shows the electron density during the laser beam transmission and panels b through e show various characteristics of the accelerated electrons much later in time. During the laser beam transmission, it produces a ring-like channel inside the target, with a narrow mostly straight compressed region at the center, as seen in Fig. 4a. The electrons injected into the laser pulse are accelerated in the positive direction along the x-axis. Figure 4e shows the number density n e in the (x, y)-plane at t = 238.3 fs and demonstrates that the electrons remain tightly grouped in bunches close to the laser axis. Figure 4b shows the corresponding areal density obtained by integrating n e of the entire electron beam along x. The areal density is concentrated in the region where the longitudinal laser electric field dominates. The average divergence angle for the electrons with r < 1 μm is 47 mrad or 2.7 ∘ . The electron energy distribution as a function of the longitudinal coordinate x and divergence angle θ ⊥ is shown in Fig. 4d, c, respectively. The profile of the maximum kinetic energy in Fig. 4d is similar to the envelope of the laser pulse, though the higher energy electrons (ε k ≳ 1 GeV) can be seen to be grouped in bunches with a length of ≲0.2 μm, which is equivalent to a temporal bunch length of ≲0.6 fs. The peaks correspond to the oscillations of the longitudinal laser electric field.
In order to emphasize the advantages of using the helical beam to inject and accelerate electrons, we have repeated the simulation discussed here using a plane-polarized Gaussian beam that has the same energy, peak power, and temporal envelope. Figure 4f shows a snapshot of the electron density while the Gaussian beam is traversing the target. In contrast to Fig. 4a for the helical beam, no dense population of electrons is being injected in the near-axis region. Instead, the electrons appear to be strongly divergent even at this early stage. This aspect, together with fundamental differences in the field topology, leads to a strongly divergent and very diffuse electron population at later times. The areal density of laser-accelerated electrons in Fig. 4g shows that the electron profile is more reminiscent of a result associated with a point source rather than with a beam.
The electron energy spectra are compared between the two runs by using all of the electrons generated by the Gaussian beam. In the case of the helical beam, we select only the electrons in the near-axis region. This is done by limiting the transverse dimensions to |y, z| < 5 μm at t ≈ 398 fs. The Gaussian laser beam produces so few electrons in the same region, marked with a dashed square in Fig. 4g, that a direct comparison would discard most of the laser-accelerated electrons. This is the reason why the plots in Fig. 4h, i that corresponds to the Gaussian beam include all of the electrons in the simulation box. Our box is a square that is 28 μm wide, so the corresponding selection criterion is |y, z| < 14 μm. It is evident from the electron spectra in Fig. 4i that the cut-off energy for the Gaussian beam is roughly 1.7 lower than that for the helical beam. Moreover, the energetic electrons generated by the Gaussian beam have a much wider angular distribution, as seen in Fig. 4h. Another way to interpret the result is by pointing out that the energetic collimated part (ε k > 800 MeV; θ ⊥ < 80 mrad) from Fig. 4c is missing in Fig. 4h.
It is also worth noting that even though both beams have the same peak power (3 PW), the helical beam has a lower peak electric field. Indeed, due to both the circular polarization and the intrinsically hollow beam profile, the peak electric field in the helical beam, jejE max =m e cω ' 42, is roughly two times lower than the peak electric field in the plane-polarized Gaussian beam, jejE max =m e cω ' 88. Despite this, the helical beam produces electrons with peak energy close to twice that of the Gaussian beam. This underlines the fact that our findings are a manifestation of qualitative differences in physics between the two types of laser beams. Fig. 4 Optimal regime of electron acceleration at normal incidence and comparison with a plane-polarized Gaussian beam. a-e Electron data from a simulation where a target with n max ¼ 0:5n c is irradiated by the helical beam. f-h Electron data from a simulation where the same target is irradiated by a plane-polarized Gaussian beam that has the same peak power of 3 PW. a, f Electron density n e normalized to the critical density n c during laser pulse transmission. b, g Areal density at t = 238.3 fs calculated by integrating n e along x. e Density of electrons accelerated by the helical beam at t = 238.3 fs. The vertical dashed lines mark the region shown in (b). c, h Electron energy distribution at t = 398.3 fs for different divergence angles The distribution in (c) is calculated for |y, z|≤5 μm, whereas the distribution in (h) is calculated for |y, z|≤14 μm. d Electron energy distribution along the x-axis at t = 398.3 fs for the case with the helical beam, where ε k is the kinetic energy. Only the electrons with |y, z|≤5 μm are shown. i Electron spectra (t = 398.3 fs), where the red curve shows electrons with |y, z|≤5 μm in the simulation with the helical beam and the black curve shows electrons with |y, z|≤14 μm in the simulation with the Gaussian beam.
Firing a laser pulse at normal incidence to a target surface, even one that has low reflectivity is not a procedure that is regularly performed in high-intensity high-power laser systems and so it is necessary to examine the electron acceleration mechanism for a case of oblique incidence. An incidence angle of 15 ∘ (to the normal direction to the target surface) is often sufficient to eliminate the danger of reflection. At the same time, this angle is sufficiently small not to fundamentally alter the interaction discussed in the manuscript. In order to implement this regime, we rotate the target rather than the irradiating laser beam, so that the laser axis is at 15 ∘ to target normally. The simulation parameters are provided in the Methods section. Figure 5 shows the results for the oblique incidence case. The electron density during the laser pulse transmission by the target is shown in Fig. 5a. The shape of the channel is asymmetric, which translates into changes in the radial position of the accelerated electrons and their density compared to the case of normal incidence. Specifically, the density of the accelerated electron population is lower by a factor of~7. However, the average divergence angle of the electrons close to the axis remains relatively unchanged. It is 50 mrad or 2.9 ∘ for the electrons within 1 μm of the central axis. The electron energy distribution as a function of the divergence angle θ ⊥ is shown in Fig. 5b. Figure 5c compares the electron energy spectra for the normal and oblique cases at the end of the simulation. Here we select only those electrons that are closer than 1 μm to the laser propagation axis in order to examine the spectrum of the collimated part of the electron population. Below 800 MeV, the two spectra are similar. The spectrum for the oblique incidence case has a slightly steeper slope and, as a result, a somewhat lower number of energetic electrons. Above 800 MeV, the spectrum for the oblique case appears to be under-resolved, which precludes us from using it to make a reliable comparison. The likely cause is the reduction in electron density. One would need to increase the number of macroparticles per cell by an order of magnitude in order to better resolve the energetic tail above 1 GeV. We can conclude based on the provided results that the investigated mechanism is sufficiently robust and that it can be used at oblique incidence without a significant increase in the divergence of the electron beam.

Discussion
We have shown that the injection and retention problems of the electron vacuum acceleration can be successfully solved by employing a thin near-critical target and a high-intensity helical laser beam with a helical index l = −1 and a radial index p = 0. The injection takes place while the laser pulse goes through the target. In contrast to a conventional beam, the retention is provided by the field of the beam itself, because the transverse fields vanish on the central axis while the longitudinal fields, including the magnetic field, peak on axis. The electrons are first accelerated inside the target and then behind the target in the vacuum by the longitudinal electric field of the transmitted pulse. Our simulations predict an electron beam with a peak energy of 1.5 GeV for a 3 PW laser. The divergence of the beam, whose charge is about 50 pC, is less than 3 ∘ . The following unique features distinguish our mechanism from other more mature approaches (e.g., laser wakefield acceleration): short acceleration distance (~100 μm); high areal density (~10 16 cm −2 ); compact transverse size (~1 μm) of the beam; and bunching of energetic electrons, with the bunches being 100 as in duration.
From an experimental point of view, the extremely short acceleration distance is a significant advantage. It opens up the possibility of performing experiments with energetic electron beams using very compact setups. Laser wakefield acceleration experiments with high-power lasers require parabolic mirrors of long focal length, typically 10s of meters, and gas targets of 10s of centimeters, whereas the setup proposed in this manuscript can be achieved with parabolic mirrors of short focal length, typically below 2 m, and very small targets. This is an important aspect because high-power laser facilities would need to have large space for very long focusing optics and such space may not be always available. Therefore, having the possibility of performing particle acceleration over short acceleration distances can be beneficial for a wide number of facilities. The resulting electron beam can be used to create a compact collider for quantum electrodynamicsbased experiments 25 and applied research (e.g., generation of highly energetic photons via inverse scattering and electronpositron beam generation).
There are also advantages to our setup when compared to a more commonly used setup for vacuum electron acceleration where the laser is reflected off a plasma mirror. In our setup, the interaction of the laser with the target is volumetric, which eliminates the sensitivity to the surface density gradient. The terminal energy gain is more efficient. Our simulations predict peak energy of 1.5 GeV for a 3 PW beam, which is comparable to the prediction for a 6 PW laser that is reflected off a plasma mirror 17 . However, our simulations also show that the laser-target interaction can become detrimental if it leads to the excitation of higher-order radial modes. This consideration imposes an upper limit on the target density. For the considered laser parameters, we determined that the electron density of the~10 μm target should remain at or below the critical density.  5 Optimal regime of electron acceleration at oblique incidence. a Electron density n e normalized to the critical density n c at t = 26.6 fs for a target with an initial maximum electron density n max ¼ 0:5n c . b Electron energy distribution at t = 386.6 fs for different divergence angles Only those electrons that have |y|≤5 μm and |z|≤5 μm are shown in this plot. c Electron energy spectra for the normal and oblique cases at the end of each simulation (t = t end ). The spectra are for the electrons that are closer than 1 μm to the axis of the laser beam.
Here t end = 398.3 fs for the case of normal incidence and t end = 386.6 fs for the case of oblique incidence.
The parameters selected for this study are such that both the laser and target specifications are within reach of the existing laser systems and target fabrication facilities. The production of helical laser beams has already been demonstrated 3 , whereas the laser pulse intensity and duration are similar to those of the Extreme Light Infrastructure-Nuclear Physics (ELI-NP) laser facility 7,8 . The capability to generate circularly-polarized beams is also already available at ELI-NP and was successfully demonstrated in recent experimental campaigns. This method is based on commercially available large aperture mica plates. The production of foams at densities as low as 2 mg cm −3 is almost routine at this point 20,21 . The considered thickness can potentially be achieved by employing a concave structured foam within a small opening (of the order of an mm) on a thin sheet. One advantage of inorganic foams such as SiO 2 is that they combine particularly fine pore structures (typically a few nm) with extremely low densities (as low as 2 mg cm −3 ) 20,21,26 . This allows for some reduction in the complexity when considering high-power highintensity laser pulses. Specifically, the pore collapse time 27,28 is likely to occur during laser prepulse due to the fine pore structure.

Methods
Particle-in-cell simulations. Table 1 provides detailed parameters for the simulations with normal incidence presented in the manuscript. Simulations were carried out using the fully relativistic PIC code Smilei 29,30 . All our simulations are 3D-3V.
The x-axis is aligned with the axis of the irradiating laser beam that propagates in the positive direction. The beam is launched at the boundary (x = −8 μm) of the simulation box by specifying the transverse magnetic field. The beam is focused (in the absence of the target) at x = 0 μm. The wavelength λ 0 and the divergence angle θ d are given in Table 1. The corresponding Rayleigh length is x R ¼ λ 0 =πθ 2 d % 35:33 μm and the corresponding beam width is w 0 = λ 0 /πθ d = 3 μm. The envelope of the pulse is Gaussian, with a pulse duration of 29.4 fs (full width at half maximum for the intensity). We limit the total duration of the envelope to 150 fs (the envelope is truncated and set to zero 75 fs before and after the peak of the envelope). The laser beam is circularly-polarized with l = −1, p = 0, and σ = 1. It is created by two linearly-polarized beams (l = −1, p = 0) with a phase offset. We use a moving window in our simulation, which reduces the size of the computational domain needed to simulate the electron acceleration for~400 fs. The window moves with the speed of light along the x-axis. The window starts moving after the beam has been fully in injected into the box. The start time is t = 28.3 fs, with t = 0 fs being the time when the peak of the envelope reaches x = 0 μm.
The target is initialized as a fully ionized silicon dioxide plasma slab. The electron density profile along the x-axis is n e ¼ n max exp½ðjxj À Δx 1 Þ=Δx 2 for 4 μm<jxj<6 μm; where Δx 1 = 4 μm and Δx 2 = 0.5 μm, so that the central plane is located at x = 0 μm. The electron density is independent of y and z. The three values for n max used in our simulations are listed in Table 1. In each case, the densities of oxygen and silicon ions are chosen as specified in Table 1 to ensure that the target is quasineutral.
In order to examine the generation of high harmonics, we have performed two simulations with a higher spatial resolution than that listed in Table 1. The new resolution is increased to 40 cells/micron in the x-direction and to 27.4 cells/ micron in both y and z directions. These simulations have a larger box size of 36.8 μm × 28 μm × 28 μm. They are run without using the moving window.
One additional simulation is performed for a case of oblique incidence. In this simulation, the target is tilted such that the front surface normal is at 15 ∘ to the xaxis. One of the initial domain boundaries is moved by 3 μm from x = −8 μm to x = −11.5 μm and the width of the simulation box along x is increased to 31.5 μm to keep the original spatial resolution. All other parameters remain the same.
We compare our results for the optimal regime at normal incidence with those for a linearly-polarized Gaussian pulse. The comparison is performed using the simulation setup for a target with n e = 0.5n c where the original beam is replaced by a linearly-polarized Gaussian beam that has the same energy, peak power, and temporal envelope. The transverse electric field of this beam is directed along the yaxis, whereas its magnetic field is directed along the z-axis. The normalized peak amplitude is a 0 = 88.6.
Mode decomposition. A snapshot of an electric or magnetic field generated by the 3D PIC code can be presented as a superposition of the ψ p;l ðe x;e r; ϕÞ modes with different twist indices l and radial indices p. In the transmitted beam, all of the modes propagate in the positive direction along the x-axis, which means that their complex field has the expðiξÞ ¼ expðikx À iωtÞ multiplier, where k = 2π/λ 0 is the amplitude of the wave vector. However, this information is not automatically encoded in the real field values provided by the code. To illustrate the ambiguity, consider a real electric field with a longitudinal profile given by E 0 cosðkxÞ at t = 0. This field can be the real part of E 0 expðikx À ωtÞ or the real part of E 0 expðÀikx À iωtÞ, which corresponds to a mode propagating in the negative direction along the x-axis. This field can even be the field of two counterpropagating modes, with the complex field given by 0:5E 0 expðikx À iωtÞ þ 0:5E 0 expðÀikx À iωtÞ.
The ambiguity can be removed by recovering the missing information using the Hilbert transform along x. The Hilbert transform converts the real field provided by the PIC code into a complex field, where the imaginary part is such that the field represents a beam propagating in the positive direction along the x-axis. We denote the field obtained via the Hilbert transform using the superscript H. For example, the Hilbert transform converts the real longitudinal field (at time t) E x ðe xÞ ¼ À g½ξðtÞE max k 1 þ e x 2 sin ξðtÞ À 2ð1 þ pÞtan À1 ðe xÞ 1 þ e x 2 cos ξðtÞ À 2ð1 þ pÞtan À1 ðe xÞ The resulting complex field E H x is a snapshot of the field of a forward propagating mode given by Eq. (10).
The fields obtained via the Hilbert transform can be presented as the following superposition: jejB H z ðe x;e r; ϕ; tÞ m e cω ¼ ∑ where, without any loss of generality, we used the z-component of the magnetic field to provide an explicit expression. Here b p,l (ξ) is a dimensionless complex amplitude that has both real and imaginary parts. If B H z ðe x;e r; ϕÞ is known at time t, then one can recover all complex amplitudes by performing a double integral, b p;l ðξÞ ¼ for all possible combinations of l and p. The meaning of b p,l can be readily understood when considering a simple signal B H z ðe x;e r; ϕÞ ¼ B 0 ψ p 0 ;l 0 gðξÞ expðiξÞ. This signal consists of an amplitude of B 0 , a single Laguerre-Gaussian mode ψ p 0 ;l 0 , an envelope g(ξ), and an oscillating term expðiξÞ. In this simple case the resulting decomposition, when trying all values of p and l, from Eq. (29) will give a single nonzero amplitude of b p 0 ;l 0 ðξÞ B 0 gðξÞ, where g(ξ) is the temporal envelope of the original signal. The right-hand side of Eq. (29) is calculated for ξ at time t, so the result is a function of e x. In order to find the corresponding dependence on ξ, we take into account that ξ ¼ 2e x=θ 2 d À ωt at a given time t. We have benchmarked this method using a 3D PIC simulation (without a target) for a helical circularlypolarized beam with p = 0, l = −1, and σ = −1. Following our procedure, we have recovered b 0,−1 (ξ) from a single snapshot of the real field B z ðe x;e r; ϕÞ outputted by the PIC code.
The representation given by Eq. (28) is unique for a given w 0 . The value of w 0 sets a characteristic transverse scale. A single helical mode with a beam width w b is represented by a single term if w 0 = w b . On the other hand, if w 0 ≠ w b , then, in general, it is represented by multiple terms. The number of terms with an appreciable amplitude can be adjusted by adjusting w 0 . We use this approach when post-processing the transmitted beam whose w b is not known to minimize the number of terms in the mode decomposition. Specifically, we minimize the sum M ∑ by varying the width that we denote as w 0 0 to distinguish it from the width of the original beam w 0 . This approach assumes that the mode with l = −1 and p = 0 is the dominant mode in the transmitted beam, which is indeed the case for all our simulations. It is worth noting that M = 0 at w 0 0 ¼ w 0 for the irradiating beam because it can be represented by a single mode with l = −1 and p = 0. Figure 6 shows the mode decomposition for two simulations discussed in the main text: a simulation with n max ¼ 2:5n c and a simulation with n max ¼ 0:5n c . The values of the sum defined by Eq. (30) for different w 0 0 are shown in Fig. 6. We find that M reaches its local minimum at w 0 0 =w 0 ¼ 0:8 for n max ¼ 2:5n c and at w 0 0 =w 0 ¼ 0:9 for n max ¼ 0:5n c . In both cases, the modes with l = −1 dominate. The plots of 〈b p,−1 〉 for these optimal values of w 0 0 are shown in Fig. 6a, b, d.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.