Mapping the power-law decay of high-harmonic spectra from laser-plasma interactions

Visible or near infra-red light can be manipulated to produce bursts of coherent extreme ultraviolet (XUV) or X-rays via the relativistic high-order harmonic generation process when a laser irradiates a solid plasma target. The intensity of the spectral components of the reflected signal decays with increase of harmonic order and the efficiency of this non-linear process largely hinges on how prompt this decay is. This is governed by the conditions of the laser-plasma interaction for which various models have been proposed. At relativistic intensities, a spectrum exhibiting a power-law decay with an exponent of 8/3 or 4/3 is often stated. Here, we analyse the dependence of this exponent on interaction parameters including the angle of incidence, the carrier envelope phase, intensity of the laser and the pre-plasma length, and discuss opportunities for optimization. Our simulations show that, rather than there being one universal exponent, the spectral decay is a continuous function of the laser-plasma interaction parameters.


I. INTRODUCTION
When a relativistically intense laser is focused on a plasma target, high-order harmonics of the incident laser light are generated around the specular direction. Highorder harmonic generation (HHG) by the laser-plasma interaction process can serve as a compact source of coherent XUV or soft X-ray pulses of attosecond duration. Unlike the HHG mechanisms in gaseous media [1], HHG mechanisms driven by laser-plasma interactions can withstand much higher field intensities and thereby offer an advantageous alternative [2]. The generated bright ultrashort pulses from HHG process have multiple applications ranging from studying warm-dense matter [3], non-destructive inspection techniques [4], studies of strong-field quantum electrodynamics [5] and capturing the ultra-fast physical, chemical or biological processes by serving as the fastest cameras [1].
Many experiments have demonstrated HHG from plasmas [2,[6][7][8][9][10][11] and distinct theoretical models explain the non-linear process responsible for the observed spectra [12][13][14][15][16][17][18][19]. The characteristics of the emission strongly depend on the initial parameters and configuration of laser-plasma interaction. Two forces, the laser driver force (∝ a 0 laser pulse amplitude in relativistic units) and plasma restoring force (∝ n e , electronic density normalised by the classical critical density n c ), compete with each other giving rise to a convenient similarity parameter S = n e /a 0 [20]. All mechanisms fundamentally depend on the coupling of laser's energy with that of the electrons at the plasma-vacuum interface and have distinct regions of dominance, efficiencies and signatures of the reflected spectrum. In the case of weakly relativistic laser intensities, a 0 1, coherent wake emission (CWE), where electrostatic plasma oscillations are excited by the wake of Brunel electrons [13] and the emission process is limited by the local plasma frequency, dominates. At higher intensities, the relativistic HHG mechanism is described by 'relativistic oscillating mirror' (ROM) model, where the harmonic emission is due to the Doppler upshifted laser frequency due to the relativistic motion of a reflecting oscillating point [14,18]. This model assumes the Leontovich boundary condition implying that the incoming and outgoing energy fluxes at some oscillating point are equal to each other at each instant of time, i.e. the energy is not accumulated by the plasma. At higher intensities and when moderate plasma densities are chosen (so as not to allow the onset of relativistic transparency), a significant part of incident laser energy (∼ 60%) can get deposited into plasma, due to the relativistic motion of electron nanobunches against relatively stationary ions. This process leads to release of short bursts of energy (with > 100 times the initial intensity) and can be modelled by 'relativistic electron spring' (RES) model [16,21]. The spectral trends of the bursts can be described by coherent synchrotron emission (CSE) [9,15] where the accelerated nanobunches generate coherent emission while following synchrotronlike trajectories. The high-order harmonics of the input laser in the reflected emission are phase locked and their spectral intensity decays with increasing frequency. Due to the non-linear nature of interaction between the electrons at the interaction surface of the solid plasma and intense laser driver, the generated harmonics show a long-tail frequency distribution. Theoretical modelling and PIC simulations suggest that this decay follows a power law [14,15,17]. The ROM model predicts that exponent of power-law decay of harmonic intensity is 8/3, based on the similarity theory highlighting that this law is universal and is independent of laser-plasma interaction conditions [14,17]. However, this is challenged in ultra-relativistic regime where the interaction is too complex and a range between 5/3 to 7/3 is predicted [22]. The RES model, on the other hand, predicts a slower exponential decay [16] which can be fitted to a power law with an exponent 1.31 [21]. The CSE model predicts a value in the range 4/3 to 6/5 [15]. A recent study on these exponents and cutoffs (that mark a faster decay) shows that this can even exceed 8/3 for a wide range of parameters [23]. The exponent is important because it determines how quick the decay is, which affects the IR to XUV/X-ray conversion efficiency and so the feasibility of HHG as a radiation source.
In this paper, we investigate the power-law exponent exhibited by the reflected spectrum for a laser-plasma system with a fixed similarity parameter, but different configuration parameters. These parameters include the laser's angle of incidence θ, the carrier envelope phase φ CEP and the length scales L p of an exponentially falling pre-plasma. These parameters determine the initial conditions of interaction and play a pivotal role in deciding which mechanism will dominate or interfere [7]. We show that the exponent of spectral power-law decay for S = 64.87 spans a range of values from 5.4 to 2.2, in contrast to specific values that are widely considered as universal (for ROM, 8/3, or CSE, 4/3). We identify an optimal region within this parameter space (θ, φ CEP and L p ), where the high-harmonic generation process is most efficient and has a slower decay, where the RES mechanism is in play. Subsequently, choosing the best value of φ CEP from this optimum region, we extend our investigation to higher laser intensities (lowering the S parameter down to 1.9), which further supports our findings.
The paper is organised in the following way. In Sec. II, we describe the details of a series of simulations performed with varying laser-plasma configuration parameters. Thereafter in Sec. III we discuss the characteristics of the reflected electric field as a function of the chosen parameters. The exponent of the spectral intensity decay and conversion efficiency are presented. Then we extend our parameter scan to higher intensities in Sec. IV. In Sec. V we summarize and discuss our findings. Appendix VI discusses the numerical reliability of our findings and presents analytical estimates of the optimal pre-plasma scale length.

II. SIMULATION DETAILS
A series of PIC simulations has been performed using ELMIS-1D [24] where a p-polarised laser of wavelength λ 0 =765 nm and intensity I 0 = 8.9 × 10 19 W/cm 2 (with vector potential a 0 = eE 0 /m e ω 0 c = 6.1659), irradiates a Si 4+ plasma at an angle of incidence θ. Here E 0 is the electric field strength of the laser, e, ω 0 ,m e and c are the elementary charge, the laser frequency, the electron mass and the velocity of the light in vacuum. This oblique-incidence problem is reduced to one dimension by the Bourdier technique [25], where a frame boosted by c sin θ along the plasma surface is chosen. The laser has the following profile: E y = E ′ 0 Θ(−ψ)Θ(−4π − ψ) sin(ψ/4)[0.5 cos(ψ/4) sin(ψ)+sin(ψ/4) cos(ψ+φ CEP )]. Here E ′ 0 = E 0 cos θ is the wave amplitude in the moving frame, Θ is the Heaviside step function, φ CEP is the carrier envelope phase, L m is the maximum length in laser's propagating direction x, ψ is the phase coordinate 2π(x − ct)/λ ′ , with λ ′ = λ 0 / cos θ being the wavelength of the laser in the moving frame. This configuration is close to a Gaussian profile of duration ∼3 fs at full-width-half-max, FWHM and has been chosen to ensure that the vector potential of a few-cycle laser has the same value as the coordinate approaches infinity. A strong dependence of the peak electric field with the φ CEP of the ultra short pulse laser can be observed in Fig. 1. The plasma is initialised with a temperature of 1 keV and the ions in the plasma are mobile with masses m i /m e = 28 × 1836, (m i being the mass of the ion). The plasma slab (with a step function profile) is 2 µm wide and has a peak density of 400 n c , where n c is the classical critical density for a laser pulse (n c = m e ω 2 0 /4πe 2 ) leading to S = 64.87. This is preceded by a preplasma that decays exponentially from the peak with a scale length L p . These simulations have been performed varying the angle of incidence θ ∈ [0 • , 75 • ] with a step size of 2.5 • , the carrier envelope phase φ CEP ∈ [0 • , 180 • ] with a step size of 6 • and the pre-plasma scale length is varied from L p ∈ [0.08, 0.73]µm with a step size of 0.05µm. This amounts a total of 14415 simulations. The length of the simulation box is L x = 16λ ′ extending from −8λ ′ to 8λ ′ . Fig. 2 shows the variation of the the spatial cell size (∆ x = L x /n x , with n x being the number of cells) with the angle of incidence of laser. It can be clearly seen that with increasing angle of laser incidence, the cell size increases very rapidly in the boosted frame. Thus, these automated simulations have been performed in two sets with different resolutions and are limited to a θ max = 75 • so as to trade off numerical accuracy against run-time. The number of cells in these simulations for angle of incidence from [0 • , 50 • ] is 64 × 1024 = 65536 and in those from [50 • , 75 • ] is 128 × 1024 = 131072. This has been done in order to ensure that the spatial resolution ∆ x remains small (between 0.24 nm for 0 • to 0.47 nm for 75 • ) even at larger angles to avoid numerical instabilities, as the size of the simulation box in the transformed frame is θ-dependent. The number of particles per cell is 85 and there are 256 timesteps per plasma period.

A. Peak electric field
The top row of Fig. 3 shows the the highest factor by which the incident laser field has been amplified in the reemitted spectrum, based on these simulations. Here, each subplot (panels [(a) − (f )]) shows the maximum strength of the reflected electric field, normalised by the incident laser field amplitude (η = |E ref /E 0 | max ), as a function of all angles of incidence (θ ∈ [0 • , 75 • ]) and carrier-envelope phase (φ CEP ∈ [0 • , 180 • ]) for certain pre-plasma scale lengths (labelled on the top of each panel). The preplasma scale length is increased from left to right with L p ∈ [0.08, 0.73]µm in the each row (corresponding to [0.105, 0.954]λ 0 ). Interestingly in each subplot, one can distinctly observe the formation of an optimum region in θ − φ CEP space, where the field amplification is significantly larger (η >> 1) than in the other regions. Beginning with the first subplot of shortest pre-plasma length (L p ∈ 0.08µm), where the plasma is closest to a perfectly sharp mirror-like boundary, this region is relatively small, 40 • θ 60 • and 30 • φ CEP 100 • . Moreover, in panel Fig. 3(a) η is not remarkably above unity in most of the region, and reaches 1.8 in the optimum region. The color gradient around this region indicates a smooth increase from the lower θ and φ CEP values followed by a similar subsequent decay. In low L p simulations, ROMlike conditions dominate the laser-plasma interaction and as the L p is gradually increased, the optimum region expands in θ − φ CEP space (notably wider in φ CEP ) yet approaches a much narrower region of maximum field strength gain (η ∼ 3 times). This augmentation, related to pre-plasma length, is expected as longer pre-plasma lengths would lower the effective S parameter, allowing an improved energy accumulation in plasma favouring the RES mechanism [16,26]. One can see in (panel (e)) that a very narrow peak of η within this optimum region at [θ, φ] = [55 • , 90 • ] which starts showing clearly from L p = 0.28µm and reaches its maximum at L p = 0.68µm, after which it disappears.
This can be better seen in Fig. 4, which shows the variation of η with different pre-plasma lengths for a range of angle of incidence θ ∈ [40 • , 50 • ] and a particular value of φ CEP = 102 • from the optimum region. Here we see a sharp increase in η as the pre-plasma length is increased up to ∼ 0.2µm beyond which η grows very slowly. A black star on each line of Fig. 4 marks a theoretical prediction of the optimal pre-plasma scale length, as given in Sec. VI B.
which ensure to S eff ∼ 1. The length scales below L opt indicate the interaction scenarios when S eff > 1. For an ultra-short laser pulse, like the one used in this work, φ CEP has a significant impact on the amplitude and position of the first peak of the electric-field profile, which drives the initial stages of the laser-electron interaction.
The effectiveness of the electron nanobunching subsequently affects the efficiency of the high-harmonic emission. Simulations indicate that φ CEP ∼ [90 • −100 • ] allow for a favourable interaction scenario for the case considered.
There is an appearance of a small well-defined island ( Fig. 3(f)) of potentially enhanced emission at lower an- where the field amplification is much enhanced (η = 2 for longer pre-plasma length L p = 0.73µm). Also, at very high θ > 70 • at L p = 0.73µm, there is yet another increase in the value of η ∼ 2. However, these could also be an outcome of numerical effects due to lower resolution at large θ. In order to convincingly identify this, we plot the peak field amplitude within the plasma region after the interaction as a measure of numerical effects in Fig. 9 (Appendix.VI). The impact of resolution at highangle of incidence leading to rising numerical instability is captured in Fig. 9, in the top row (panels [(a − f )]), where each subplot shows the maximum value of lasernormalised field in the spatial region of plasma (E err ), far away from the reflected signal whose peak is shown in Fig. 3. In Fig. 9, one can clearly see that the error increases with θ and is close to zero at lower values. This clarifies that the increased value of η in the top layer of Fig. 3 at θ > 70 • at L p = 0.73µm is potentially an outcome of numerical effects. Nevertheless, the distinct broad optimal region for enhanced emission as well as the small island (at about [θ, φ] = [5 • , 70 • ]) may be relied upon.

B. Decay exponent
The spectrum obtained after reflection of the laser in each configuration is unique as its intensity falls with increasing laser harmonic order. A slower decay of intensity or a smaller exponent is naturally favourable for effectively up-converting the laser frequency. In order to relate the optimum regions, in the first row of Fig. 3, where the reflected electric field boost is maximised, we plot the exponent of the decaying spectral power-law in the second row of Fig. 3 (panels(g) − (l)). This exponent is calculated by finding the slope of a linear fit to the log-scaled intensity ln I(ln ω) spectrum obtained from PIC simulations. It should be pointed out that the spectrum obtained from these simulations excludes the field contribution due to numerical effects shown in the top row of Fig. 9 shown in App.VI. The corresponding ad-justed R 2 value to the fit on each spectrum is plotted in the lower row of Fig. 9 indicating the goodness of fit captured in the exponent of the power law in Fig. 3. Clearly, the adjusted R 2 value is above 0.9 in most of the cases (meaning 90% of the spectral data can be explained by this fit) until about θ = 65 • , above which it varies from 0.8 to 0.7 with φ CEP . It should be noted that the range of harmonic spectrum chosen for data fitting with θ < 65 • is 2 < ω/ω 0 < 100 (corresponding to the energy being 2.62 eV to 162 eV, ω 0 = 1.62eV). For θ > 65 • a bit longer range is taken, 2 < ω/ω 0 < 250 as this leads to an overall improvement of the fit (R 2 min from 0.6 to 0.7) on the spectral data implying a subtle deviation from power-law behaviour. Having clarified the regions of a good fit on decaying spectral power-law, the exponent trends in the second row of Fig. 3 (panels (g) − (l)) can be analysed more meaningfully.
In panels (g) − (l) of Fig. 3, it can be seen that there is a smooth dependence of the exponent of the powerlaw on the experimental parameters. The clear global observation is the convergence of the exponents to a relatively lower value in the same θ − φ CEP space as the optimum region of high-η in top row of the same Fig. 3. This indicates the region of slower decay of harmonics with the potential to find rare events of high-frequencyhigh-intensity radiation. Interestingly, the range of exponents widely deviates from 8/3, the value expected by the universal law even in cases with smaller preplasma lengths, where ROM-like conditions dominate (see L p = 0.08, 0.13 µm). In the Fig. 3 (g) showing L p = 0.08 µm, the exponent in the optimum region converges to a moderate value of 3.5. Nevertheless, this optimum region is widely spread out in θ − φ CEP . In the non-optimal regions in Fig. 3 (g), the exponent value can also reach up to 5.8 (at around θ = 0, φ CEP = 156 • ) implying an undesired fast decay of harmonics and a significant deviation from the universal value even where the plasma configuration would allow ROM like conditions. As the pre-plasma lengths are increased towards the right from L p = 0.13 µm, the RES mechanism of HHG begins to dominate as the effective S parameter of interaction begins to get closer to unity. Here, the optimum region starts to converge to a more specific region of a lower absolute value of exponent. In Fig. 3-(k) (L p = 0.68 µm), the value of exponent reaches upto 3.1 for a very specific θ = 50 • and φ CEP = 80 • . This precisely coincides with the point of maximum η in the optimum region in top row of Fig. 3.
In order to crosscheck our simulation results, we also numerically solve the RES equations [21] for a subset of the parameter scans under consideration. The dependence of the decay exponent on φ cep , observed in the PIC simulations, is consistent with that predicted by the RES theory, which can be seen in Fig. 5. Here we restrict the fit to the frequency range corresponding to the first ten harmonics, which yields a decay exponent m 10 . This is because the frequency range in which the RES model is accurate is itself intensity-dependent [21].
In other regions of Fig. 3-(k), the furthest exponent value reaches up to 5.4. Moreover, at high θ−high φ CEP regions (particularly Fig. 3-(h, i)) the value of exponent for the spectral power law reaches around 2.2. However, the corresponding adjusted R 2 value ( Fig. 9-(h, i) in App.VI) reduces sharply to 0.68 in this region from 0.95 in the optimum region, which indicates a clear deviation from power-law behaviour. Interestingly, the small island of enhanced emission at lower angles (at about [θ, φ] = [5 • , 70 • ]) for L p > 0.4 also has a lower exponent (m = 3.1 for L p = 0.68µm, same as in the optimum region) implying the possibility of stable high harmonics at near-normal incidence. Overall for all simulations with a fixed S-parameter, the exponent of spectral power-law decay is found to vary in the range 5.8 > m > 2.2. The lower limit of the values, that mostly appears at low-θlow-φ and high-θ-high-φ regions, is very close to the ones found in a recent study that reports the exponent value to be in the range 6 > m > 4 for a single cycle laser with θ = 30 • , S ∼ 60, L p = 0, φ CEP = 0 [23]. The narrowing down of the optimum region of stable harmonics with increasing L p can be related to HHG generation by the RES mechanism.

C. Conversion efficiency
The third row of Fig. 3 [panels (m−r))] shows the fraction of laser energy (ǫ inc ) going into the entire reflected spectrum (ǫ ref ). It can clearly be seen that at lower preplasma lengths (esp. L p = 0.08µm, panel (m) ) almost all of the laser energy goes into the reflected signal. This is understandable as this configuration with a relatively sharp boundary allows for ROM-like conditions. As the pre-plasma lengths are increased, one can observe that in the optimum region of enhanced and stable re-emission, the laser energy in the re-emitted signal reduces in a small area which precisely coincides with the optimum regions in the corresponding panels above. This points to the region of RES dominance, where this reduction in laser energy going into the re-emitted spectrum can only be because of its deposition into plasma-electron nanobunches. Clearly, even at lower pre-plasma lengths L p = 0.08µm (Fig. 3(m)), there is a certain amount of laser energy that is being deposited into plasma (∼ 10 − 12%). Later, a clear trend of lowering of the laser energy imparted into the re-emitted spectrum in the optimum region can be seen in Fig. 3(o − r). In Fig. 3(q), at L p = 0.68µm, this total laser energy being imparted into the reflected signal is reduced to 0.62 at the same location (of θ−φ CEP ) where we observe the maximum field amplification and the lowest exponent of power law (θ = 50 • and φ CEP = 90 • ) in the corresponding panels above Fig. 3(a − l). This would imply that about 38% of the incident energy could be used up in electron heating at the solid boundary. Nevertheless, the fraction of incident energy (ǫ inc ) going into higher order harmonics is increased. This is captured in the lowest row Fig. 3[panels (s − x))] which clearly shows that the energy imparted into the reflected spectrum above the 10th harmonic (ǫ ref,≥10 ) also increases precisely in the same optimum region. Here as well, one can observe the expected trend of RES enhancement of energy going into higher harmonics at longer pre-plasma lengths within the same θ − φ CEP region. So, longer preplasma lengths that modify the effective S-parameter are favourable in compressing laser-plasma energy into higher harmonics as also seen in Ref. [26].

IV. HIGHER LASER INTENSITIES
Fixing the carrier-envelope phase φ CEP = 100 • from the optimum region, we now investigate the dependence of the exponent of spectral power-law decay for different intensities I 0 ∈ [8.9 × 10 21 , 1.0089 × 10 23 ] W/cm 2 . The other parameters are same as before, leading to a decreased bulk S parameter. This is shown in Fig. 6 which shows the peak electric field boost η (top row), powerlaw exponents (second row), and the fraction of laser energy (ǫ inc ) into the entire reflected spectrum (ǫ ref ) and spectra above the 10th harmonic (ǫ ref,≥10 ) (lower two rows respectively). This has been arranged in the same fashion as Fig. 3. One may notice small arrows systematically pointed at certain angles for all panels above 0.13µm. These arrows point to those precise angles for each pre-plasma length where the physical quantities in each row display an unusual trend (seen as a sharp change in colour) and shall be discussed later in this section.
First we see in the top row of Fig. 6 that, like the previous section, there is a clear optimum around 50 • where η is higher. In Fig. 6(a), we see a larger η ∼ 4 than before that is also θ−dependent. This is expected, as the relativistic motion of electrons in the plasma enhances RES. The electric-field increase η quickly increases with L p and peaks around L p = 0.13 − 0.23µm (Fig. 6(b − d)). Unlike in Fig. 3, this optimal region of high η decreases rapidly with increasing L p after reaching its peak (see Fig. 6(e − f )), as seen in Ref. [26]. This can largely be due to the fact that the laser intensity is held constant in the previous set of simulations, whereas here the laser intensities are varied over orders of magnitude. Around L p = 0.18µm (Fig. 6(c)), η can be as high as 8. This pre-plasma length of maximum gain is also closer to the optimal length predicted by RES theoretical estimates. This can be seen in Fig. 7, which shows a distinct peak. The black star in the figure is the theoretical estimate of L opt computed as before by Eq. (1). This length is also consistent with the widely reported finding that highharmonic generation is optimized for a pre-plasma with a similar length [6,7,11].
The second row shows that, even for higher intensities, the decay exponent does not appear to have a universal value. At the optimum angle of 50 • and L p = 0.13 − 0.18, the value of the exponent reaches as much as 1.16. Furthermore, similar to the row above, the exponent at this angle reduces with L p to around 2.5. This points to the efficient conversion region that can be seen in the two rows below. The conversion efficiency is reduced to around ∼ 0.8 in this optimal region (50 • and L p = 0.13 − 0.18). In the same region of lowest panel ( Fig. 6(t,u)), the energy in higher harmonics is higher, pointing to XUV conversion efficiencies of around 20%.
Beyond L p = 0.13 ( Fig. 6(o-r)), the loss in the reflected spectrum energy (up to 0.2) at lower angles is not directly related to but is an offshoot of RES mechanism. The arrows mark certain lower angles for each length scale, that break up the plasma surface. At these angles (with φ cep = 100 • ), a part of electron nano-bunches that drive RES are extracted from the surface and co-propagate with the re-emitted radiation. The radiation therefore loses some of its energy to electrons and the spectrum deviates from a power law. This can be seen consistently in all the panels. This lower angle θ r points to some L p -dependent resonance, and can be distinguished as the abrupt line of lower values of all parameters under discussion [η, m, ǫ ref and ǫ ref >10 ]. These have also been corroborated with higher resolution simulations and are also in agreement with RES theory which can be seen in Fig. 8. Here, the RES theory and PIC simulations predict a similar trend and a resonant peak in the exponent values is clear (This seems shifted by few degrees potentially due to sparse PIC data). Excluding these odd lines at θ r , one may clearly see the distinct stable radiation signatures at 50 • and L p = 0.13 − 0.18 with an exponent going beyond 4/3 [15] or 1.31 [16].

V. SUMMARY AND CONCLUSIONS
To summarize, we have carried out a systematic study of the decay exponent of the spectral intensity of high harmonics, for a fixed laser and plasma (same bulk S parameter), as well as higher intensities (lowering bulk S parameter), for different configuration parameters, includ- ing the angle of incidence and the carrier envelope phase of the laser and the pre-plasma length of the plasma. We have shown that the exponent of the power law can span a range of values for different parameters (5.8 to 1.1). Even for configurations where plasma mirror reflects the majority of the incident energy, the exponent non-universally deviates from 8/3. Optimal parameter spaces of enhanced emission and their dependence on the laser-plasma configuration parameters are identified and related to the mechanisms at play. We show that at lower pre-plasma lengths the region of moderate exponent is spread out widely in θ − φ CEP space while, at longer pre-plasma lengths, a lower exponent value can be reached in a more precise θ − φ CEP space. We have investigated these exponents in relation to other features of the re-emitted spectrum, such as reflected peak field amplitude and laser-energy conversion efficiency. We have highlighted regions in the parameter space that are not optimal for applications, despite the fact the spectral decay is relatively weak, because the conversion efficiency and peak-field strength are not large. Thorough understanding of the dependence of the exponent on laserplasma parameters, to which this study contributes, is paramount to the usefulness of high-harmonic generation as a potentially rich source of XUV or X-rays.

A. Numerical accuracy
Simulations are always subject to possible numerical inaccuracy. In order to gauge their reliability, we plot the peak field amplitude within the plasma region after the interaction as a function of the same parameter space in the top row of Fig. 9. This serves as a measure of numerical error and clearly, the results in the optimal region are devoid of the numerical instabilities in the plasma region. We also present the goodness of fit of the slopes (shown in the lower row of Fig. 3 that is found by a linear fit of the simulation spectra) in the lower row of Fig. 9. The adjusted R 2 value is close to 0.9 in the majority of the fits performed, while its reduction at high θ indicates a deviation from power-law behavior. It must be pointed out that the re-emitted spectrum on which the fitting is performed, precludes any involvement of the fields in the plasma region originating as numerical noise (shown in the top row of Fig. 9).

B. Preplasma dynamics
Here we repeat the analysis in Section 3 (Preplasma) in Blackburn et al. [26], but for an exponential, rather than linear, density ramp. The incident laser pulse penetrates the preplasma up to the point where the force of the charge separation field, established between the electron layer and the background ions, is balanced against the total magnetic force, established by the laser field and the uncompensated ion current. Assuming that all the electrons in the region z ′ < z ′ max are displaced, and that the ions are not displaced, we find that this balance occurs at The condition on L p in 3 ensures that the laser pulse does not break through the density ramp. This will be the case for all the results presented here. This result gives the maximum displacement of the electron layer, which occurs when the plasma interacts with the peak intensity of the laser pulse. We can derive the instantaneous displacement of the electrons ξ(t ′ ) by assuming that the temporal profile of the laser varies adiabatically and therefore that force balance applies at every instant of time between the charge separation field and the instantaneous magnetic field of the laser. We therefore replace a 0 in 3, which is implicit in the S factor, with a(φ) = a 0 exp[−φ 2 /(2τ 2 )], where the phase φ = ωt ′ cos θ and τ = ωT /(2 √ ln 2). The plasma density encountered by the laser at phasefront φ is n 0 exp[ξ(φ)/L], from the definition of the density ramp. The value of the normalized density at this point follows as s(φ) = n 0 exp[ξ(φ)/L]/[a(φ)n c ]. The dependence on phase cancels out, leaving s(φ) equal to a constant we label the "effective S" parameter: This is, rather interestingly, independent of the bulk value of S, due to the self-similar properties of the exponential function [11]. If we increase n 0 at fixed L, 3 tells us that the point at which the pulse is reflected, z ′ max , moves further out into the preplasma. This compensates for the increase in n 0 , because the local value of the density n e (z ′ max ) = n 0 exp[z ′ max /L p ]. . The top row shows the error defined as peak field value within the plasma region after its interaction with the laser (normalised by the peak field of the laser). The bottom row shows the adjusted R 2 value indicating the goodness of fit for the exponents of the spectral power-law scaling for high-order harmonics (shown in Fig. 3).
The relativistic electron spring model predicts that XUV generation is optimized at S ∼ 1, θ ∼ 60 • . We can identify that the scale length L necessary to reduce 4 to unity is eq.1. This is consistent with the widely-reported finding that high harmonic generation is optimized for a preplasma with L ≃ 0.25λ [6,7,11].