The ability to manipulate light at relativistic intensities1 (>1018 W cm⁻² for a laser wavelength of ∼1 μm) is an important aspect of the development of compact laser-driven particle and radiation sources.2,3 Such sources have wide-ranging potential applications including medicine,4–6 ultrafast imaging,7 industrial processing,8 and inertial confinement fusion.9,10 Tailoring the spatial and temporal profile of the laser pulse can be used to optimize focusing behavior in laser–plasma interactions,11 which can enhance laser-driven particle acceleration12,13 and has the potential to produce radiation sources with orbital angular momentum (OAM).14

The conventional approach to manipulating light is based on using solid-state optical media, but there is a limit to the energy density that these media can withstand before nonlinear effects distort the laser pulse and the optics themselves are damaged. This is typically circumvented by increasing the size of the optics as the laser energy is increased so that the overall energy density is below the critical value. However, large-area optics are very expensive, with long manufacture times and limited maneuverability because of their volume and weight. By contrast, plasma optical components are inherently compact, easily achieving energy densities of more than two orders of magnitude higher. In addition, the ultrafast evolution of the optical properties of laser-excited plasma enables other properties of the laser pulse to be tailored. By using plasma-based optics, it is possible to enhance the temporal intensity contrast of the laser pulse,15–17 refocus it to enhance the laser peak intensity significantly,18–20 compress its temporal profile,21 generate high harmonics,22–24 and amplify the pulse.25–27 Through the interaction of intense light with dense plasma oscillations, the production of transient plasma gratings,28–30 holograms,31 and coherent emission of light32 can also be achieved.

In recent studies involving intense laser-pulse interactions with ultrathin (tens of nanometers) foils, we introduced the concept of a self-generated relativistic plasma aperture that enables the laser light to pass through a region of an initially opaque target, resulting in diffraction and manipulation of the particle dynamics in the corresponding field structures.33,34 This plasma aperture forms because a region of the target becomes relativistically transparent35 to the laser light, which occurs when the plasma frequency ω_p [=nee2/(γϵ0me), where n_e is the electron density, e is the electron charge, γ is the electron Lorentz factor, ϵ₀ is the permittivity of free space, and m_e is the electron rest mass] is reduced to below the laser frequency ω_L. The electron density at which this occurs is known as the relativistically corrected critical density γn_crit, and at currently achievable laser intensities, the process is aided by target heating and expansion reducing the peak electron density in the focal-spot region. As relativistically induced transparency occurs above an intensity threshold, the spatial profile for a Gaussian focal spot is circular, forming a relativistic plasma aperture in the target.33 As the transmitted portion of the laser pulse interacts with the edge of this plasma aperture, intense light is generated at ω_L and 2ω_L with higher-order transverse spatial mode structures and different polarization states compared with the drive laser pulse.36 An example result from a 3D particle-in-cell (PIC) simulation of the interaction of an intense pulse with an initially solid-density 10-nm-thick aluminum target is shown in Fig. 1(a), in which can be seen the self-generated relativistic plasma aperture and 2ω_L light emitted in the TEM₀₁ mode. This was driven by a linearly polarized (in Y) 6-J, 40-fs-duration laser pulse focused to a peak intensity of 6 × 10²⁰ W cm⁻², in a fundamental TEM₀₀ mode; see Duff et al.36 for a full description of the high-order spatial mode generation mechanism. When the incoming laser light is circularly polarized, light with OAM can be produced at the higher harmonics of the laser frequency.37

The laser–plasma-aperture interaction effectively acts to transfer energy from the incoming laser pulse—which has a fundamental Gaussian spatial mode structure—to intense harmonic light in higher-order modes. This could potentially be developed as a plasma optical device. Although the use of ultrathin targets is challenging because of their fragility, which can limit the shot repetition rate, the scheme could potentially be simplified by using micrometer-thick targets incorporating a preformed circular aperture to produce the same effect.36 This is illustrated in Fig. 1(b), which is an example 3D simulation result showing the generation of relativistically intense 2ω_L light in the TEM₀₁ mode in a target with thickness l = 2 μm and aperture diameter d = 3 μm. These targets could then also be used as a platform to test more readily the underlying physics that occurs during the formation of a relativistic plasma aperture.

Herein, we investigate the potential use of preformed-aperture targets for generating frequency-doubled light with higher-order spatial mode structures. Using 2D and 3D PIC simulations, we determine the optimum target parameters to maximize the efficiency of energy transfer, and we show that this is related to the fraction of light transmitted through the aperture, the strength of the longitudinal electric field at the aperture edge, and the absorption of the laser pulse. Using 3D simulations, we compare the generation efficiency and mode structure formation when the polarization is varied from linear to circular. The results highlight the feasibility of using a preformed-aperture target—resulting in a plasma aperture—as a type of plasma optical device for mode conversion of intense light.

To investigate and optimize the generation of frequency-doubled light in higher-order spatial modes, 2D and 3D simulations were performed using the fully relativistic PIC code EPOCH.38 For the 2D case, the simulation XY grid was defined as 1000 × 720 mesh cells, corresponding to 20 ×   20 μm². For the 3D simulations, the XYZ grid was defined as 1000 × 720 × 720 cells, corresponding to 20 ×   20 ×  20 μm³. The targets comprised a slab of Al¹³⁺ ions of thickness l with a circular aperture of diameter d. This was neutralized with an electron population with a peak density of 30n_crit and an initial temperature of 10 keV. Each species was initialized with an average of 80 particles per cell in the 2D case and five particles per cell in the 3D case. All the boundaries of the simulation box were defined as free-space ones, and the ions were mobile in all cases.

In both the 2D and 3D simulations, the pulse of laser light with wavelength λ_L = 800 nm entering from the left (X) boundary and was focused at the plane corresponding to the front surface of the target (X = 0). Both the laser-pulse temporal and spatial profiles were defined in terms of intensity as Gaussian with pulse duration τ_L = 40 fs [full width at half maximum (FWHM)] and focal-spot diameter ϕ_L = 2.1 μm (FWHM). Time t = 0 was defined as when the peak of the pulse reached the front surface of the target, i.e., the front of the aperture. The 3D simulations were performed using linear polarization in the Y direction and circular polarization. In the 2D simulations, only linear polarization in the Y direction was used because circular polarization cannot interact with the full aperture edge. The peak laser-pulse intensity I_L was varied from 10²¹ to 10²³ W cm⁻² by increasing the energy in the pulse. In the 2D simulations, a simulation automation code called BISHOP39—developed by our group—was used to systematically scan a 2D parameter grid of target thickness and aperture diameter; l was varied from 0.5 to 10 μm and d was varied from 0.5 to 5 μm.

Figure 2(a) shows the electron density for a target with l = 0.5 μm at t = −4.8 fs from a 2D simulation, demonstrating the generation of electron bunches from the predefined (d = 3 μm) aperture edges as the laser interacts and propagates through; at each side, the electron bunches are emitted with frequency ω_L, π out of phase with the opposite side (because of the oscillation of the laser electric field vector), with an overall frequency of 2 ω_L when considering the bunches from both sides of the aperture. Figure 2(b) shows an example 2D map of the longitudinal electric field 6 fs thereafter; the pattern is produced by the interference from the two point sources at the opposite edges of the aperture. To investigate the temporal generation of 2 ω_L light, an arbitrary electron bunch—circled in Fig. 2(a) with its corresponding position shown 6 fs later in Fig. 2(b)—is tracked through the simulation. In Fig. 2(c), the total kinetic energy K_e in the bunch is plotted as a function of X (the axis of laser propagation), together with the generated 2 ω_L intensity I_2ωL and the longitudinal electric field experienced E_X (both averaged 400 nm around the bunch). The electron bunch is accelerated immediately to relativistic velocities from the aperture edge, before experiencing a decelerating longitudinal electric field as it propagates with the laser and into the sheath field at the target rear (formed because of the charge separation of earlier fast electrons and rear surface plasma surrounding the aperture40,41). This can be observed in Fig. 2(b) as the overall positive electric field (in the positive X direction) at the rear of the target, together with the longitudinal field associated with the laser pulse. Both the acceleration and deceleration of the electron bunch result in broadband emission that interferes constructively with the emission from neighboring bunches, resulting in the generation of 2 ω_L light. This is due to the periodic nature of the accelerated bunches produced at either side of the aperture, which have an overall frequency of 2 ω_L. As a result of this, the intensity of the 2 ω_L light along the longitudinal axis increases with distance from the aperture for ∼5 μm, after which it saturates and then reduces as the generated 2ω_L light expands away from the sampled region in the vicinity of the bunch.

Figures 2(d)–2(f) show the corresponding results for the case l = 10 μm (d = 3 μm), for which the aperture effectively becomes a channel. In this case, Fig. 2(e) shows the longitudinal electric field 24 fs after the electron density plot in Fig. 2(d), when the bunch has reached the rear of the target. As can be seen in Fig. 2(f), the initial acceleration of the electron bunch and generation of 2ω_L light are similar to those in case of l = 0.5 μm. However, as the electron bunch and laser co-propagate within the extended plasma aperture/channel, the electrons experience the positive longitudinal field associated with the focused Gaussian pulse (due to the divergence of the laser electric field profile, i.e., ∇⋅E⃗=0 in free space), resulting in a deceleration of the bunch. The bunch never comes to rest, and it continues to propagate in the positive X direction with a relativistic velocity at all times. As the bunch propagates, it begins to dephase with the laser pulse and experiences the negative longitudinal field associated with the subsequent half-wave period, accelerating the bunch again. In both the deceleration and acceleration phases, the electrons emit, and the constructive interference with the neighboring emission produces a strong 2ω_L signal. As the bunch leaves the aperture/channel at X = 10 μm, it experiences the sheath field at the target rear and again decelerates, resulting in further 2ω_L emission.

The aperture diameter and target thickness influence the energy conversion to 2ω_L. To determine the influence of these parameters, a detailed parameter scan was performed to observe the changes to the generation efficiency, based on the ratio of the total energy of the generated 2ω_L light detected 5 μm from the target rear to the energy of the input laser pulse. Figure 3 shows the conversion efficiency as a function of target thickness and aperture diameter for three different laser intensities. For I_L = 10²¹ W cm⁻², the optimal efficiency is found to be d ∼ 2 μm, shown in Fig. 3(a) as region I. At this aperture diameter, the efficiency saturates for l ≳ 6 μm. When the aperture diameter is increased, the edges of the aperture experience less interaction with the longitudinal field of the laser pulse and the efficiency of 2ω_L thus reduces, as observed in region II.

Region III indicates that for small d, the efficiency decreases with increasing l. This is due to poor matching conditions between the free-space Gaussian laser beam focus and the close-to-cut-off modes that can exist within the aperture when the thickness is increased beyond the wavelength; this acts to attenuate the beam within the aperture. As d is reduced from the optimal at 2 μm (for the 2.1-μm FWHM laser focal spot), the generation efficiency therefore becomes very dependent on l. For the thickest targets, the efficiency reduces to negligible levels because the beam is attenuated greatly. However, reducing l enables appreciable efficiency to be regained, as indicated by region IV. This is due to the aperture evolving with the laser pulse. While small apertures inhibit the propagation of laser light because of poor matching conditions, the transverse ponderomotive force of the focused light is intense enough to expel the plasma from the focal-spot region. If the target is thin enough, this expulsion is sufficient to increase the size of the aperture and continue to generate 2ω_L light, similar to that of the interaction with an ultrathin foil target as reported by Duff et al.36 This effect can be observed more readily as the peak intensity is increased. At I_L = 10²² W cm⁻², as shown in Fig. 3(b), the generation efficiency for the smallest d occurs over a much larger l range, and region IV widens in the 2D parameter map. As the intensity is increased further to I_L = 10²³ W cm⁻², as shown in Fig. 3(c), region IV becomes broad across the entire d and l ranges, resulting in no observable region III. This effect is likely to be exaggerated by the reduced peak density (in comparison to the solid density of ∼430n_crit) used in the simulations. In this scenario, I_L is high enough at the edges of the aperture that the heating of the plasma by the laser is sufficient for the interaction region to become relativistically transparent near the peak of the pulse, regardless of d. This acts to reduce the efficiency of the generation mechanism because the effective aperture diameter is now much larger than the starting value, in all cases. The optimal condition at d ∼ 2 μm still exists because the rising edge of the laser pulse (with less intensity) still interacts with the initial aperture diameter.

Next, we explore the observed dependence of the 2ω_L generation efficiency on d and l. It is shown above that during the interaction of the laser pulse with the aperture edge, the longitudinal electric field acts to accelerate and decelerate the bunched electrons. It is therefore assumed that the magnitude of the field close to the aperture edge plays a role in defining the efficiency of 2ω_L generation. The longitudinal electric field for a pulse with a focused Gaussian spatial profile can be determined from the spatial differential of the electric field along the polarization direction. The maximum absolute value of this field is at a transverse position of Y ∼ ϕ_L/2. As such, as d is increased, the magnitude of this longitudinal electric field at the aperture edge, E_X, is maximized when d ∼ ϕ_L; this is shown in Fig. 4(a). For a pulse that is linearly polarized in the Y dimension, E_X is defined asEX(Y)=ddYe−2ln(2)YϕL2.(1)

The fraction of the laser light that can enter the aperture is another factor that influences 2ω_L generation. The greater the energy that can propagate through the aperture, the greater the potential conversion into second-harmonic light. The energy in the pulse that can enter the aperture, K_in, is shown in Fig. 4(b) for the 2D simulations and is given byKin(d)=∫−d/2d/2e−2ln(2)YϕL2dY.(2)This varies with d because more of the laser pulse is blocked for smaller apertures. As I_L is increased, relativistic transparency at the aperture edges can occur, enabling more light to enter, but this is a second-order effect and is not included in this simple model.

A third parameter influencing conversion to 2ω_L light is the overall absorption of the energy of the laser pulse into the electron bunches, Abs_L, as it propagates through the aperture. We assume that this varies linearly with l with an absorption coefficient α because the number of electrons exposed to the propagating laser pulse increases with l. We also assume an exponential attenuation of the laser pulse with thickness with an attenuation constant β. Because the actual attenuation of the laser pulse is complicated by effects such as self-focusing and reflection, the attenuation constant is treated as a free parameter, approximated based on the level of transmitted light in the simulation results (β = 0.1). This is shown in Fig. 4(c), and the overall absorption is given asAbsL(l)=αle−βl.(3)

Assuming that these parameters are related linearly gives the following first-order approximate relationship for the 2ω_L conversion efficiency η2ωL:η2ωL∝|EX(d/2)|AbsL(l)Kin(d).(4)The variation of this relationship with d and l is plotted in Fig. 4(d), which matches well with the simulation results in Fig. 3. This simple model is useful only to demonstrate the behavioral trend of the conversion efficiency. Factors such as relativistic transparency and the reabsorption of 2ω_L generated light are not included: relativistic transparency will affect only the thinnest targets and the highest-intensity pulses, and the reabsorption of 2ω_L light will become significant only for very thick targets.

Because harmonic light from the plasma aperture can be produced in non-symmetric modes,36,37 it is necessary to extend the investigation to include 3D simulations. Because of the substantially increased computational requirements in doing so, a reduced set of simulations was performed to check whether the efficiency dependencies observed in 2D are the same in 3D. The results are presented in Fig. 5(a). The overall trends in the dependencies on l and d are similar to the 2D results, noting that the maximum conversion efficiency is reduced to ∼3.5% and the optimal d is found to be slightly larger at 3 ± 0.5 μm (compared to 2 μm for 2D). The reduction in conversion efficiency is likely due to the fact that in 2D, the intensity is simulated only in the plane of the laser polarization where the polarization vector is perfectly normal to the aperture edge. In 3D, there is no significant interaction when the polarization vector is parallel to the aperture edge, resulting in proportionally less energy available to produce 2ω_L light. The additional degree of freedom also means that the accelerated electrons are less confined to the laser field as they propagate. The total area of the laser focus, as well as the area of the aperture, must be considered in 3D, which modifies Eq. (3) as follows:Kin(d)=∫−ππ∫0d/2e−2ln(2)rϕL2rdrdϕ,(5)where the integral is given in cylindrical coordinates. Figure 5(b) shows the calculated η2ωL with Eq. (5) used to account for the changes in 3D. This results in an increase of the optimal d by ∼0.5 μm.

For a direct comparison of the dependencies of η2ωL on d and l in 2D and 3D, line-outs from both sets of simulation results are presented in Fig. 5(c) for fixed l = 10 μm and Fig. 5(d) for fixed d = 2 μm in 2D and d = 3 μm in 3D (i.e., the optimum diameter in each case). The model calculations based on the approximate conversion relationships are also shown for the 2D and 3D cases and are found to be in good overall agreement with the trends observed in the simulation results.

The 2D and 3D simulation results exhibit similar overall trends. The efficiency increases rapidly with d to an optimum value, and as the aperture size is increased beyond this, η2ωL decreases more slowly with d in the 3D case. The optimum d is also ∼0.5 μm higher in 3D, as predicted by the simple model. As l is varied (for the optimal d), the efficiency rises slightly faster in the 3D case compared to the corresponding results in 2D, before saturating at ∼3.5%. These discrepancies may be due to the optimal d being slightly larger in the 3D case. As the laser pulse expands after focus, the intensity decreases as r² in 3D compared to r in 2D, which will change the effective intensity. Second-order effects such as self-focusing42,43 and reflection of the laser pulse as it moves through the aperture may also induce discrepancies between 2D and 3D, because these process can change the intensity and longitudinal electric field.

To illustrate the self-focusing and scalloping behavior in the plasma aperture, Fig. 6 shows the laser intensity profile of the drive laser pulse in two example targets with differing l.

Figures 6(a) and 6(b) show the ω_L intensity profile in the XY plane along Z = 0 from the 3D simulations for l = 0.5 μm at example times t = 16 and 32 fs, respectively. The laser pulse is observed to experience self-focusing as it passes through the aperture to a position of X ≈ 2.5 μm before expanding out into free space. Note that the intensity increases from the nominal intensity of 1 × 10²¹ W cm⁻² (without plasma) to ∼4 × 10²¹ W cm⁻². This is twice as high as can be achieved with the same focus in the 2D simulations because of dimensionality effects (the pulse energy can be focused only in a plane in 2D but radially in 3D). When l is increased to 10 μm, as shown for the same example times in Figs. 6(c) and 6(d), the self-focusing behavior is very similar, but the expansion is inhibited by the extended aperture/channel in this case, resulting in scalloping behavior. At the later time step [Fig. 6(d)], the pulse is observed to have refocused toward the rear of the target.

In the 3D simulations, the spatial structure of the generated 2ω_L light takes the form of a TEM₀₁ mode. The exact spatial distribution of the mode depends on the aperture geometry. Figures 7(a)–7(d) show the change in the mode structure as d is increased for fixed l = 6 μm. This is temporally integrated over the full pulse duration. For the smallest aperture (d = 1 μm), low levels of 2ω_L light are generated (note the change in fluence scale) but the TEM₀₁ mode is still observed. As d is increased, the mode structure becomes very clear and the overall beam size is observed to reduce. As the aperture size is increased, there are also low amounts of 2ω_L light produced outside of this central mode, with even higher-order structures. This is caused by the interference of the generated light varying as the sources are moved farther apart.37 These can be seen more readily on a logarithmic scale, as shown in Fig. 8(a).

All the simulation results discussed thus far involve linear polarization. When the drive laser pulse is circularly polarized, high harmonic light is produced with a Laguerre–Gaussian (LG) mode structure.37 To explore the effect of polarization on the 2ω_L conversion efficiency, a further subset of simulations was also performed for circularly polarized drive laser light. In this case, the peak intensity was halved to maintain the same pulse energy as in the linear pulse case.

Figures 5(c) and 5(d) show the conversion efficiency as a function of d and l with circular polarization for the same fixed parameters as in the linear cases. The overall trends are similar to the linear polarization cases, but with a general increase in overall efficiency in conversion to 2ω_L light. This is due to the polarization vector always being normal to the aperture edge in the case of circular polarization, resulting in 2ω_L light being produced over the entire wave period.

Figures 7(e)–7(h) show the variation in spatial distribution of the temporally integrated 2ω_L light for circular polarization as d is varied. The generated mode takes the form of a ring distribution,37 as expected from an LG mode. Similar to the linear polarization case, the aperture diameter affects the size of the main central ring and low-level higher-order ring-like distributions are evident as separation of the electron bunches increases (i.e., with increasing d), resulting in a change to the interference pattern. This is seen more clearly on the logarithmic scale in Fig. 8(b).

In practice, depending on the temporal intensity contrast of the laser pulse, some degree of ionization and expansion of the target may be induced early in the interaction, prior to the arrival of the main pulse. Various techniques can be applied to improve the contrast of the pulse, including using plasma mirrors,15–17 but depending on the final intensity profile of the leading edge of the laser pulse, preplasma may still be generated, resulting in an expansion of plasma toward the center of the aperture. This will result in a reduction of the effective aperture diameter as the critical surface moves inward, causing the 2ω_L light generation to have the characteristics of that produced by a smaller-aperture target.

The rising edge of the laser may also produce a degree of underdense plasma within the aperture, which may affect how the laser pulse propagates. To test the impact of this, a series of 3D simulations was conducted in which the inside of a d = 3 μm, l = 6 μm aperture target was filled with plasma of varying density. Figure 9(a) shows the influence of this on the efficiency of conversion to 2ω_L light. For densities less than 0.1n_crit, there is very little (a few percent) reduction in conversion efficiency. The reduction is more significant at higher densities (∼60% at n_crit), but the targets are not expected to be used under a temporal intensity contrast so low as to induce such a closure of the aperture. A full investigation of this would require 3D magnetohydrodynamic simulations, which are beyond the present scope.

The results presented so far involve targets for which the entrance and exit apertures are the same size, creating a parallel channel in cases with a large thickness-to-diameter ratio. It is possible to optimize the 2ω_L light generation further by shaping the target geometry to control the self-focusing of the laser pulse. An additional series of 3D simulations was run in which the walls of the channel created within the aperture were angled (forming a conical target) by defining different entrance and exit diameters d₁ and d₂, respectively. Figure 9(b) shows the efficiency of conversion to 2ω_L light as a function of d₁ for fixed l = 6 μm and d₂ = 2 μm. The conversion efficiency is improved by increasing the size of the entrance diameter, reaching a maximum of ∼6.9% for an optimum d₁ = 8 μm. The optimum occurs when the convergence angle of the laser undergoing self-focusing matches the angle of the cone walls. This can be seen readily when comparing Figs. 9(c) and 9(d). For the case d₁ = 3 μm [Fig. 9(c)], self-focusing of the propagating light results in a focal position inside the target, similar to Fig. 6(c). For d₁ = 8 μm, the propagating laser light self-focuses close to the exit aperture. For higher d₁, the focal position moves increasingly beyond the exit aperture and the laser pulse interacts with a shorter length of plasma wall, because the laser pulse needs to propagate deep within the aperture before interacting with it. It may be possible to tailor the geometry of these targets further to improve the conversion efficiency, and the optimum may also depend on the temporal intensity contrast of the laser pulse.

Following our previous work on the self-generation of a relativistic plasma aperture during the interaction of intense laser light with an ultrathin foil33 and the resultant harmonic generation in higher-order spatial modes,36 we have performed a numerical investigation to explore the use of a preformed plasma aperture to achieve a similar output with a simpler target. The results showed that incoming laser light in a fundamental TEM₀₀ mode is converted to relativistically intense second-harmonic radiation emission in higher-order modes. The efficiency of conversion was explored as a function of aperture diameter and target thickness, and it was found that there is an optimal diameter at which the conversion efficiency saturates with increasing thickness. This was shown to be proportional to the intensity of light that can enter the aperture, the strength of the longitudinal electric field at the aperture edge, and the absorption of the laser energy. This study was performed for linear polarization in both 2D and 3D, which overall exhibit the same trends, allowing a multi-dimensional parameter space to be explored more readily using computationally efficient 2D simulations. Extension to circular polarization in 3D produced broadly similar trends in conversion efficiency.

For the parameter range explored, the optimum target conditions for generating 2ω_L light in the TEM₀₁ mode yielded a conversion efficiency of ∼6.9%, while the optimal efficiency for generating a 2ω_L LG mode was found to be ∼4%. There is potential for this to be increased further by tuning the profile of the tapered targets. As the intensity of the input pulse is increased, the efficiency is reduced slightly, but conversion efficiency of this order still produces 2ω_L light that is sufficiently intense to induce highly relativistic plasma oscillations. It has also been shown that there is scope for tuning the size of the mode structure and the mode order by varying the aperture diameter. This demonstrates that a preformed-aperture target irradiated by an intense laser pulse can function as a plasma optical device for generating 2ω_L light with a spatial mode structure that can be changed from non-OAM to OAM by varying the drive laser polarization.

Acknowledgment. This work was supported financially by EPSRC (Grant Nos. EP/R006202/1 and EP/V049232/1) and STFC (Grant No. ST/V001612/1). It involved the use of the ARCHIE-WeSt and ARCHER2 high-performance computers, with access to the latter provided via the Plasma Physics HEC Consortia (Grant No. EP/R029148/1), and the University of Cambridge Research Computing Service (funded by Grant No. EP/P020259/1). EPOCH was developed under EPSRC Grant No. EP/G054940/1. The research has also received funding from Laserlab-Europe (Grant Agreement No. 871124, European Union’s Horizon 2020 research and innovation program).