The X-Ray Emission Effectiveness of Plasma Mirrors: Reexamining Power-Law Scaling for Relativistic High-Order Harmonic Generation

Ultrashort pulsed lasers provide uniquely detailed access to the ultrafast dynamics of physical, chemical, and biological systems, but only a handful of wavelengths are directly produced by solid-state lasers, necessitating efficient high-power frequency conversion. Relativistic plasma mirrors generate broadband power-law spectra, that may span the gap between petawatt-class infrared laser facilities and x-ray free-electron lasers; despite substantial theoretical work the ultimate efficiency of this relativistic high-order-harmonic generation remains unclear. We show that the coherent radiation emitted by plasma mirrors follows a power-law distribution of energy over frequency with an exponent that, even in the ultrarelativistic limit, strongly depends on the ratio of laser intensity to plasma density and exceeds the frequently quoted value of −8/3 over a wide range of parameters. The coherent synchrotron emission model, when adequately corrected for the finite width of emitting electron bunches, is not just valid for p-polarized light and thin foil targets, but generally describes relativistic harmonic generation, including at normal incidence and with finite-gradient plasmas. Our numerical results support the ω−4/3 scaling of the synchrotron emission model as a limiting efficiency of the process under most conditions. The highest frequencies that can be generated with this scaling are usually restricted by the width of the emitting electron bunch rather than the Lorentz factor of the fastest electrons. The theoretical scaling relations developed here suggest, for example, that with a 20-PW 800-nm driving laser, 1 TW/harmonic can be produced for 1-keV photons.

The first high-order harmonics were generated from solid-density plasmas in the sub-relativistic regime [32][33][34][35][36][37][38] , where the highest frequencies produced are limited by the plasma frequency [39][40][41] . For relativistic intensities the nonlinear frequency conversion mechanism is fundamentally different, allowing much higher efficiencies and generating harmonics well above the plasma frequency 5,6 ; experimentally demonstrated photon energies have reached 3.8 keV ( ≈3200 harmonics) 23 . Advances in laser technology are now making attosecond sources based on RHHG feasible 42 . The relativistic harmonic generation mechanism is often described as a relativistically oscillating mirror (ROM), where the incident laser field is Doppler upshifted by the longitudinal motion of the reflective plasma surface 6,10 , though for strong laser-plasma interactions sub-cycle energy storage in the plasma becomes important and a more detailed treatment of the electron motion is required, e.g. the relativistic electron Figure 1. Schematic of relativistic harmonic generation with a p-polarized ultrafast laser incident on a solid target. The spectrum was experimentally measured using the Princeton 20 TW Ti:sapphire laser system (Amplitude Technologies), with 70 mJ incident on a BK7 target in 25 fs to produce a peak intensity of 8 × 10 19 W/cm 2 . The central fundamental wavelength is 800 nm. The reflected single-shot spectrum was passed through a 150 nm aluminum filter to remove the fundamental and imaged with a flat-field diffraction-grating spectrometer. The plasma frequency limits sub-relativistic harmonics to 40 nm. Note that the harmonic structure appears because the driving laser pulse is multi-cycle; a single-cycle driver would produce an isolated attosecond pulse and a continuum spectrum. the full spectrum of reflected radiation can be higher than that of the driving laser, making RHHG a possible route to intensities beyond the limits of CPA 12,50,51 . The power of an attosecond pulse which could be achieved from a 10 PW driver based on simulation results is also marked 15 .
We show that the power-law exponent p of the reflected spectrum is a continuous function of the experimental parameters, even in the relativistic limit, and that p can be less than 8/3 across a wide range of conditions. Drawing on particle-in-cell simulations, we delineate scaling relationships for the power-law exponent and cutoff and relate them to observed properties of surface electrons, including the Lorentz factor, in an effort to understand both the range of validity of proposed models and the reachable efficiencies of high-order-harmonic generation.

The Power-Law Spectrum of RHHG
RHHG is primarily characterized by the spectrum of the reflected light. Most theoretical models and simulation results are consistent with a power-law spectral shape, although there is some disagreement regarding the appropriate exponent 5,8,10,45,49 . A power-law fit seems to reasonably approximate the computed spectra we observe. For each spectrum in Fig. 2, the initial power-law [ ω ω ω ∝ − I( ) ( / ) L p , where ω L is the laser frequency] eventually gives way to a steeper falloff. We will define the cuffoff frequency ω c as the frequency at which the spectral intensity has dropped by a factor of e below that predicted by continuing the power-law scaling, a definition which is reasonably robust to the variation in spectral shape that can occur near the cutoff; ω c is marked for the two spectra in Fig. 2. The value of p is important because it can be used to calculate an approximate value for the energy in each harmonic with ω < ω c . In this section, we consider how p varies due to the laser intensity and plasma density and relate this variation to the microscopic motion of plasma electrons.
Without atomic effects, the interaction physics may be scaled by the laser frequency, removing explicit dependence on ω L . The strength of the incident laser is described by its normalized field amplitude: a 0 = eE/m e ω L c, where E is the maximum value of the electric field envelope, m e and e are the mass and charge of an electron, and c is the speed of light. The plasma density (n e ) is normalized by the critical density ( ω = n m c e L 2 /4πe 2 ) to give N = n e /n c . Since the two forces of interest in the problem come from the laser and the plasma, the dynamics of relativistic RHHG are governed by the balance between the laser driving force (∝a 0 ) and the plasma restoring force (∝N) 15 . The similarity parameter (S = N/a 0 ) 52-54 arises from nondimensionalization of the Vlasov-Maxwell equations in the relativistic limit (a 0 ≫ 1) and generally plays an important role in relativistic laser-plasma interactions, although its exact relevance for RHHG has been questioned 49,[55][56][57] . The generation of high frequency light via relativistic RHHG is a subcycle process, so to isolate the detailed mechanism we will primarily show results for single-cycle (Gaussian) driving pulses. Figure 4a shows the variation in p for 1 ≤ a 0 ≤ 1000 and 0.01 ≤ a 0 /N ≤ 5 at angle of incidence θ L = 0. The value of p is found from a linear fit to the logarithm of the spectrum vs the logarithm of the normalized frequency in the range 1 < ω/ω L < 100; the fitting procedure is described in more detail in the Methods section. There are two key points of interest in Fig. 4a. First, for a 0 → ∞ at fixed a 0 /N, p approaches a constant value that can be approximated as proportional to the logarithm of a 0 /N, i.e. ∝ →∞ p a N a N lim , for a wide range of values. Second, for a 0 /N > 0.15, p exceeds the predicted maximum scaling at normal incidence (p = 8/3) 10 . This corresponds to where the maximum reflected intensity exceeds the incident intensity. Note in particular that as a 0 /N is varied, a very large range of values for p can be reached even in the limit a 0 → ∞, from p larger than 6 to less than 2. The example spectra at right show how the spectrum of reflected radiation varies at fixed a 0 /N and www.nature.com/scientificreports www.nature.com/scientificreports/ varied a 0 for selected values of a 0 /N. As a 0 increases, the spectra follow the power law over a larger frequency range; the cutoff increases with a 0 .
At moderate values of a 0 the fitted value of p can be far from the asymptotic limit. Neglecting deformation of the plasma surface, the electric field at the vacuum-plasma interface for normal incidence and a 0 ≪ N can be approximated as: To derive this, note that at the plasma-vacuum boundary, the complex valued incident E I,0 , reflected E R,0 and transmitted (E T,0 ) electric fields must satisfy E I,0 + E R,0 = E T,0 and the magnetic fields satisfy B I,0 + B R,0 = B T,0 . In vacuum, E 0 = B 0 , and the plasma dispersion relation gives ,0 . The surface field E S,0 is simply E T,0 at the interface. Since the transverse Lorentz factor of electrons at the surface will be proportional to the surface field, the condition for relativistic motion is a 0 /  N 1, not simply a 0 ≫ 1. When we filter the plotted points by this condition, and fit p in the range 0 < ω < 20 to avoid the spectral cutoffs, we find that p is a function of only a 0 /N (Fig. 4a inset); there is no variation with a 0 for lower-order harmonics. This suggests the relevance of the similarity model to RHHG and means we can attribute the deviations in p observed for smaller a 0 , particularly at smaller values of a 0 /N, to fact that the motion is not fully relativistic.
The same basic dependence of the spectral power law on a 0 and N persists at non-zero angle of incidence for P-polarized light, shown, for example, at θ L = 30° in Fig. 5b. As for θ L = 0°, p approaches a constant value for a 0 → ∞ at fixed a 0 /N and decreases for a 0 /N < 0.5. Although for both θ L = 0° and θ L = 30° the fitted value of p does not quickly increase for a 0 /N > 1, Fig. 5a shows that the attosecond pulse generation efficiency (I a /I L ) reaches a maximum near a 0 /N ≈ 0.5 and falls for larger a 0 /N. This decrease occurs for all spectral filters, e.g. 0 < ω/ω L (red), 4 < ω/ω L < 100 (black, solid), 50 < ω/ω L (black, dashed), as well as for the total reflected energy (solid blue) and reflected fundamental energy (dashed blue). This corresponds to the relativistic transparency regime, where a 0 ≫ N and the high-intensity laser pulse will propagate into plasma which is overdense for less-intense light. In this regime, the laser force overpowers the plasma restoring forces, leading to inefficient RHHG. The relatively shallow values of p for a 0 /N > 1 do not capture the total decrease in the spectrum as a result of the transparency of the plasma in this regime. The most efficient attosecond pulse generation is achieved when the laser and plasma forces are close to balanced.
As the angle of incidence (θ L ) changes, the dependence of p on a 0 /N as a 0 → ∞ remains relatively consistent. Some variation in the exact values of p that are found at a given a 0 /N can be seen in Fig. 6, where the a 0 → ∞ limit is approximated by a 0 = 1000. For a 0 /N < 0.5, p can be very roughly estimated, to within ±1, by ≈ + p N a C 2ln ( / ) 0 , shown by the diagonal dashed line. More importantly, at all angles considered (θ L = 0°, 30°, 45°, and 60°) the shallowest values of p lie in the interval 0.1 < a 0 /N < 1, with a limit near p = 4/3. Apart from a few points for θ L = 45°, this −4/3 line appears to be a limit for the power-law exponent across all of the angle of incidence, intensities, and densities considered here.

Electron Trajectories
The −4/3 spectral scaling arises from the CSE model of high-order harmonic generation 45 by assuming that the radiation of interest is emitted by a single bunch of electrons over a time when its position in the direction of specular reflection can be approximated for small t as: x L L 1 3 www.nature.com/scientificreports www.nature.com/scientificreports/ and its transverse current varies linearly as: where t = 0 is the time of emission, β x = v x /c is the normalized velocity component in the x direction, and α 0 and α 1 are constants. Using the stationary phase approximation, this leads to 45 : where Ai′ is the derivative of the Airy function of the first kind, ω α γ = is the Fourier transform of the shape function f describing the finite width of the electron bunch. For ω ≪ ω γ and an electron bunch much narrower than 1/ω, we recover the power law spectrum: 4/3 For this to be valid, ω/ω L must also be sufficiently large that t/T L is small in the expansions of x(t) and j y (t) and higher order terms can be safely neglected.
The mean and individual longitudinal and transverse velocities of an XUV-emitting electron bunch are shown in Fig. 7, along with the individual Lorentz factors in the direction of emission and the advanced time [t a = t + x(t)/c] for each electron. To reach Eq. (3) we must fit the longitudinal velocity with a parabola and the transverse velocity with a line, as marked by the dashed curves. The width of the region where these approximate the underlying trajectories sets the lowest harmonic frequency which will follow Eq. (3); 4/3-scaling comes  www.nature.com/scientificreports www.nature.com/scientificreports/ from the third-order approximation of x and it is in principle possible to reach shallower power laws for a single electron if higher-order terms cannot be neglected. However, under most observed conditions, the third order term in x dominates even for relatively large t/T L (low-order harmonics) and the spectra produced by individual electrons tends to follow Eq. (3) with f = 1.
If the spectra associated with single electrons are well-predicted by Eq. (3) with , deviations of the total reflected spectrum from this shape can be attributed to either the finite distribution in time and space of the electrons as they emit or the range in maximum Lorentz factors, and therefore ω γ , of the radiating electrons.
Appropriate treatment of f allows the finite distribution of electrons to be included. The curves in Fig. 7 represent a selection of particles that reach the highest velocities, rather than all macroparticles in the simulation. These are the electrons that will dominate the emitted radiation, particularly at high frequencies, but the contribution of slower electrons to the total spectrum can result in relatively higher emission in low-order harmonics and a steeper power-law slope.
The electron-bunch kinematics shown in Fig. 7 produce characteristic synchrotron-like trajectories in the x-y plane 14,46 . The key attribute of these trajectories is that velocity and acceleration are perpendicular when the electrons are traveling in the specular-reflection direction, a characteristic of synchrotron-like motion that can be found under many variations of the laser and plasma conditions. Figure 8 shows synchrotron-like trajectories from PIC simulations of a thin foil at normal incidence in one (a) and two (b) dimensions, (c) of a semi-infinite target at oblique incidence, and (d) of a finite plasma density gradient at oblique incidence. In the oblique incidence examples, the direction of emission is 30° above the y = 0 axis. In each of the subplots, the specular direction Lorentz factor peaks while the transverse acceleration is maximized, producing instantaneously circular trajectories. These characteristic synchrotron-like trajectories appear under a broad range of interaction conditions, so understanding their dynamics is critical to quantifying relativistic RHHG.
Similarity theory suggests that in the limit a 0 → ∞ at fixed a 0 /N the electron trajectories will approach an asymptote. Figure 9 shows the (a) electron bunch velocity components, (b) trajectories, (c) Lorentz factors for varied a 0 at fixed a 0 /N, and (d) the maximum Lorentz factors over a range of a 0 and N. From these plots it is apparent that the electron kinematics and trajectories for fixed a 0 /N approach a limit; the velocities and displacements for a 0 = 80 and a 0 = 100 are almost identical. This is in agreement with the collapse in Figs. 4 and 5b of p, which depends on lower-order harmonics and thus the shape of the electron trajectories over a relatively long time, to an asymptotic value in the same limit. For example, since the shape of the electron bunch velocity curves are close to identical for a 0 = 80 and a 0 = 100 over the interval 0.05 < t/t L < 0.15, the tenth harmonic should be the same proportion of the fundamental for both.
However, the Lorentz factor of the electron bunch (γ = 1/ − v c 1 / 2 2 or, when calculated only for the velocity component in the direction of reflection, γ r = 1/ − v c 1 / r 2 2 ) does not collapse and instead scales linearly with a 0 (Fig. 9c,d), which is consistent with the prediction from similarity theory that γ ∝ α S a 0 S 53 . Note that because emission occurs when the velocity points in the specular reflection direction, γ = γ r at the time of emission. To be consistent with the literature, where not otherwise noted the variable γ will indicate the value of the Lorentz factor at the emission time (which is generally the largest value of γ r ). To be valid for moderate a 0 , we write this as: The cutoff frequency is dependent on γ so the efficiency of the highest frequencies may depend on a 0 independent of a 0 /N, Both the ROM 10 and CSE 45,46 models predict scaling of the cutoff with γ 3 . The electron trajectories are therefore in agreement with the observation in Fig. 4 that p approaches an asymptotic limit for large a 0 , but the point at which the spectra deviate from p continues to vary with a 0 . www.nature.com/scientificreports www.nature.com/scientificreports/ Since we have shown that γ − 1 ∝ a 0 at fixed a 0 /N, the final question is how (γ − 1)/a 0 varies with a 0 /N. In Fig. 10 the normalized emitting-direction maximum Lorentz factor for electrons is plotted against a 0 /N for varied a 0 and N at θ L = 30°. The quantity (γ − 1)/a 0 collapses to a single line with a double power-law shape, with (γ − 1)/a 0 approximately proportional to a N ( / ) 0 3 in the region a 0 /N < 0.3 and to a 0 /N in the region a 0 /N > 0.3. Note that due the transition of α s , other normalizations of γ will not produce a single collapsed curve. If the values of γ are compared to I a /I L in Fig. 5, we can see that γ continuously increases with a 0 /N over the entirety of the optimal harmonic efficiency region, so the peak in attosecond pulse generation efficiency cannot be tied to a peak in γ. This suggests that an analysis that includes only γ in harmonic generation is missing key physics; in the next section we will address this gap by considering the spectral cutoff.

The Harmonic Cutoff
To determine the efficiency of harmonic conversion, especially for x rays, we must know the highest frequency for which power-law scaling with p is valid. The spectral cutoff sets the limit of p-scaling and is thus critical for determining the x-ray generation efficiencies that can be reached. We have chosen here to define the observed cutoff frequency ω c as the point where the spectrum has dropped below ω −p by a factor of 1/e, in part because ≈ ′ ′ [Ai (0)/Ai (1)] 1 2 /e, so that when ω = ω γ , the value of I(ω) from Eq. (3) is 1/e times that from the −4/3 scaling law. Therefore, if ω γ governs the cutoff, ω γ will approximately equal ω c .
We have previously established that γ ∝ a 0 at fixed a 0 /N, and ω α γ = γ 8 1 3 , so for a 0 ≫ 1: The similarity-scaling of electron trajectories suggests that because α 1 is derived from the shape of the electron velocity, it is a function of a 0 /N and does not strongly depend on a 0 independently. Consider Fig. 11, where reflected spectra calculated for a 0 /N = 0.33, θ = 30°, and selected a 0 between 5 and 100 are compared to Eq. (3) (gray line). We calculate the predicted cutoff ω γ from the maximum γ of the emitting electron bunch observed in the PIC simulations and fitted α 1 to the longitudinal velocity traces, which the insets show for randomly-selected particles as a function of the advanced time t a . The predicted cutoff ω γ scales favorably with a 0 : at a 0 = 100, we find ω γ = 22300, implying that for p = 4/3 the process would produce 20 keV photons with 10 −6 efficiency in each harmonic, or greater than 10 −3 energy conversion efficiency to the spectral band www.nature.com/scientificreports www.nature.com/scientificreports/ 20-22 keV. However, the a 0 = 100 spectrum actually begins to deviate from the −4/3 law around ω/ω L = 200, and the 10 keV efficiency is 10 −11 per harmonic, five orders of magnitude below the optimistic power-law-ω γ prediction. Though this still makes RHHG an extremely bright x-ray source, the substantial discrepancy requires explanation. The spectral cutoffs ω c are much smaller than ω γ in each of these plots, and, like the cutoffs in Fig. 4, appear to scale close to linearly with a 0 , not as a 0 3 . Since ω γ is often assumed to predict the maximum frequencies that can be efficiently generated, it is crucial to understand why RHHG spectra deviate from −4/3 scaling at much lower frequencies.
The trajectory analysis in Fig. 7 substantiates Eqs. (1) and (2) for the average bunch trajectory, and it has been established elsewhere that Eq. (3) accurately predicts the spectra of individual high-γ electrons 58 , so it is initially surprising that ω γ fails to match ω c so dramatically. To address the discrepancy, we consider the width of the emitting electron bunch by accurately including the shape factor term ω  f ( ) 2 that appears in the original formulation of Eq. (3). There are two important differences with how we treat the emitting electrons in comparison to the previous literature. First, we must consider the emission width as a function of the advanced time, not just approximate the width of the electron nanobunch in space at a fixed laboratory time. Since the electrons move relativistically and emission time in the laboratory frame varies with position in the bunch, the dependence of advanced time (t a ) on time (t) is complicated. Second, the distribution is not Gaussian; the electron emission distribution has a sharp leading edge, so a better approximation is the Heaviside step function, resulting in a power-law decay rather than an exponential cutoff at ω b . This sharp leading edge can be seen in the insets of Fig. 11, which show the individual particle Lorentz factors as a function of advanced time. The distribution of emission events -which is related to, but not exactly the same as, any instantaneous electron density profilemore directly affects the total emitted radiation than the electron number density. The number density of electrons is proportional to the density of drawn lines in these plots. In each plot, emission from the fastest electrons arrives at the observation point first. After the peak, the rate of emission events remains relatively constant as the maximum Lorentz factor of the involved electrons slowly decays, which is due to electrons further from the surface feeling weaker fields. The step function is an appropriate approximation because the rate of emission events does not slowly ramp up before the peak arrives; in each of these plots no emission peak is produced before t A = 0. The difference between the orange line (Gaussian) and red line (step function) spectra for a 0 = 40 in Fig. 11 illustrates the enormous impact the actual distribution of charge in the electron nanobunch has on the high-frequency components of the spectrum.
The distributions of peak γ-factors, which indicate when the highest frequency radiation is emitted by each electron, can have widths that are a substantial fraction of the laser period. For a given harmonic (n) with frequency ω n = nω L , only electrons within about half a harmonic period (T n = T L /n) of the leading edge will contribute. The advanced time electron bunch width τ a therefore sets a frequency cutoff ω b approximated as: At a 0 = 10, for example, the emission lasts for around 2% of the laser period, which implies that all of the electrons in the bunch only contribute constructively to the first 25 harmonics. From a physical perspective, if we assume that the electron emission times are uniformly distributed after the leading edge, the condition that only www.nature.com/scientificreports www.nature.com/scientificreports/ those within T n /2 constructively add means the fraction of the electron bunch that contributes to the emission of frequency ω n scales as 1/n. For coherent radiation, the radiated power will be proportional to 1/n 2 , or, equivalently, ω − n 2 . Mathematically, the square of the Fourier transform of the Heaviside step function will lead to an additional ω −2 factor for ω > ω b , i.e.: 1 . We apply this correction to the spectra in Fig. 11 to plot the red curves, using detailed particle trajectories to calculate τ a and ω b . The corrected predictions much more closely approximate the actual spectra, and the correct treatment of the bunch width becomes critically important for keV photon energies. Given that much of the detailed structure of the electron bunches has been approximated away without much justification, the analytic spectra match the PIC results remarkably well. Note that although the analytic fit has a sharp angle between the two scaling regimes, we expect both experimental and simulation spectra to exhibit a smoother transition because the onset of destructive interference is gradual.
The importance of ω b extends further than the single value of a 0 /N presented in Fig. 11 and is starkly apparent in Fig. 12, where spectra and electron γ-factors are shown for selected a 0 at a 0 /N = 1.1 and θ L = 45°. Higher values of a 0 /N produce higher γ factors (see Fig. 10), and a 0 /N = 1.1 gives γ ≈ 0.3a 0 , producing ω γ at MeV levels Figure 11. HHG spectra with CSE model fits and comparison to bunch width (ω b ) and Lorentz factor (ω γ ) cutoffs at a 0 /N = 0.33 and θ = 30°. Spectra (black) for a 0 varied between 5 and 100 driven by single-cycle laser pulses. The CSE Airy-function spectrum neglecting bunch shape is shown in gray. In red, the analytic spectrum is corrected for the observed spread of emission times by assuming that only electrons within a half wavelength of the bunch leading edge emit constructively and that electron emission events are evenly distributed in time with a sharp initial edge. This last assumption is confirmed by the inset plots of electron γ-factor against the advanced time of the individual electrons, radiation from the fastest electrons arrives first, with minimal interference, after which the bunch slowly trails off. The inset plots contain a random selection of electron trajectories. The incoherent combination of higher-order frequencies results in an additional ω −2 factor in the spectral scaling above ω b = (T L /2τ a )ω 0 . ω γ is the cutoff calculated from the Lorentz factor: ω α γ = γ 8 1

3
. The bunch width is the time between the peak γ-factor of the first and last electrons to emit. For a 0 = 40, an additional line (orange) shows the spectrum for a Gaussian emitting bunch shape, which leads to a factor ω ωω for a 0 = 100 and at 5 GeV for the, admittedly unrealistic, a 0 = 4000. However, ω b is even smaller for this larger value of a 0 /N, leading to bunch-width-dependent cutoffs at much lower frequencies. Note that for this a 0 /N the fundamental and lower-order harmonics are not strongly generated, so the first few harmonics appear to scale above the −4/3 power law. The advanced time γ-factor shapes for individual electrons are much narrower than the total time over which the emission arrives, which can be viewed as the physical interpretation of ω γ ≫ ω b .
To present a more general picture of how the different measures of cutoff vary, the variation in ω b , ω γ and ω c with a 0 is presented in Fig. 13 Fig. 11, the PIC-derived spectra are compared to analytic fits based on the CSE model with and without bunch-width correction, with the cutoffs associated with both the bunch width and the electron Lorentz factor marked. Inset: the electron Lorentz factors in the emission direction as a function of the advanced time. Note that for these conditions the width in γ r for an individual electron is even narrower in comparison to the total width of the bunch than for the trajectories in Fig. 11. for C γ found from Fig. 10 and α 1 found from the average velocity (dashed lines), both ω b and ω c increase much less rapidly, with the solid lines linearly proportional to γ. The linear γ scaling of ω b is simply an approximate empirical fit to this data and a derived relationship between the two may point to a slightly different functional relationship. The linear fit may also not be general. We can, however, say with some confidence that neither ω b nor ω c is proportional to γ 3 . In plots where ω γ lies below ω/ω L = 1000 the exponential cutoff can be identified directly from PIC spectra in addition to the ω b cutoff; for each of these simulations, ω γ calculated from the electron γ-factors closely agrees with ω γ PIC found by direct fitting to the observed spectrum. Note that the Doppler mirror cutoff 8 at 4γ 2 does not appear as a feature in any of the spectra which contribute to this plot.
We must also consider whether these relationships hold under less-efficient conditions, where p is substantially larger than 4/3. In Fig. 14, where θ = 30° and a 0 /N = 0.1, p lies around 2.5 for both a 0 = 40 and a 0 = 100. For these parameters, a significant fraction of the electrons which contribute to the reflected radiation do not reach the peak Lorentz factor associated with the overall interaction. Since the γ achieved by each electron determines the highest frequencies it will emit, the number of electrons which contribute to a particular harmonic decreases for higher harmonic. This causes a difference between the 4/3 scaling of the CSE model and the steeper scaling of the observed spectrum for a less ideal ratio of a 0 to N. In these interactions ω b and ω γ are comparable, so although the spectral fit for the bunch-width correction appears to better match the data the difference is relatively small.
The interpretation of the steep power-law-like behavior of spectra in the regime ω > ω b as due to destructive interference of electrons emitting more than λ n /2c after the sharp leading edge can be checked by filtering the reflected radiation and examining it in the time domain. Figure 15 shows the intensity of reflected radiation as a function of time for three different filters. Even for the lowest frequency filter (10 < ω/ω L < 50), we have created an attosecond pulse: the approximate duration produced by an 800 nm driver is 100 as. For 600 < ω/ω L < 3000, the process produces an attosecond pulse with 4 as duration. For filters with ω LF > ω b , the attosecond pulse is much shorter than the duration of the emitting electron bunch in the advanced time (Fig. 15d); the circles in (b) and (c) indicate the times at which electrons with individual ω γ greater than ω LF reach their maximum γ r . Without destructive interference between electrons after the leading edge of the bunch, this would be the width of the attosecond pulse for this filter. The inset in Fig. 15d shows the full-width-half-maximum duration (T a ) of the attosecond pulse as a function of the filter frequency ω LF , where ω UF = 5ω LF , which follows the scaling T a /T L = ω L /ω LF closely. ω LF and ω UF define the filter passband; the filter is transparent for ω LF < ω < ω UF . Therefore, although the efficiency of harmonic generation drops for ω > ω b , the duration of the attosecond pulses that are produced continues to shorten.

Beyond Single-Cycle Semi-Infinite Interactions
Although the use of single-cycle interactions on plasma targets with steep density transitions between vacuum and solid density allows study of the underlying RHHG dynamics, real interactions often involve multi-cycle laser pulse on targets which have developed finite exponential density gradients. We consider here how both gradients and multiple-cycle pulses affect the scaling of the reflected spectrum.
Finite plasma density gradients. Over a picosecond to nanosecond timescale after ionization of a solid density target by a laser pulse, the plasma formed at the surface will expand outwards self-similarly, forming an exponential plasma density gradient 59,60 . For relativistic surface interactions, even prepulses orders-of-magnitude weaker than the main pulse arriving tens-to-hundreds of picoseconds earlier can perturb the target enough to create a substantial plasma density gradient. The effect of pre-plasma gradients on RHHG has been studied 28,29 , and though there is some disagreement over the optimal scale length, most experimental and theoretical results point to enhanced harmonic generation for L/λ up to a few tenths of the laser wavelength. We must therefore establish how the analysis developed for plasma targets with sharp edges extends to the plasma density gradient case: here we show that (1) small-to-moderate gradients increase efficiency for a 0 /  N 1 max , (2) Figure 14. HHG spectra with CSE model fits and comparison to bunch width and Lorentz factor cutoffs at a 0 /N = 0.1 and θ = 30°. In this lower efficiency regime, the slope p is not close to 4/3. The difference between the 4/3 power law scaling (red) and the actual scaling (dashed blue) results from the substantial contribution in these spectra of electrons which have low maximum velocities compared to the peak of the electron bunch.
Let us first examine the effect of a gradient on attosecond pulse generation efficiency. In Fig. 16a, the ratio of attosecond intensity to incident intensity is plotted for a range of a 0 values and selected L/λ between 0 and 0.5. The target maximum density is = N 200 max and θ L = 30°. For L/λ = 0, this means a 0 /N = a 0 /200 and the peak efficiency of I a /I L ≈ 30 occurs around a 0 = 70. The attosecond pulses are filtered from the reflected field in the interval 20 < ω/ω L < 100. As the gradient is increased from 0 to L/λ = 0.1, the efficiency increases for a 0 / < . N 0 2 max , which can be attributed to the interaction happening at a much lower density (N ≈ a 0 ) than N = 200. Under these conditions, for L/λ > 0.1 the efficiency does not continue to increase; for gradients larger than 0.5λ the attosecond pulse become weaker. In the regime a 0 > 100 the gradient has little effect, since the interaction is in the transparency regime for the maximum plasma density. These simulations are configured so that the gradient decreases exponentially to a minimum value of 0.01n c , where it drops to zero. Weighted particles are used to provide adequate resolution across the entire range of densities with reasonable computational expense.
The key point of Fig. 16a is that although a small (L/λ ≈ 0.1) gradient can improve RHHG efficiency by orders of magnitude for a 0 /N < 0.2 -which is an important regime experimentally -the peak efficiencies reached with a gradient are similar to or less than the maximum efficiencies reached for the optimal values of a 0 /N without a gradient. In this sense, gradients are simply a way of correcting target densities that are too high for currently reachable laser intensities. This, along with the observation that electron trajectories are similar for both efficient flat targets and gradients (Fig. 8), suggests that the corrected synchrotron model is also an adequate method to treat targets with moderate gradients.
If we compare the reflected spectra produced by targets with gradient scale lengths in the interval 0 ≤ L/λ ≤ 0.2 (Fig. 16b), we find that where the zero-gradient case is close to optimum (a 0 /N = 0.4), the gradient does not substantially change RHHG efficiency. For larger gradients, the reflectivity in the interval 2 < ω/ω L < 30 drops, which flattens the apparent initial scaling law, but does not increase η ω 1 2 is the fraction of energy in the between harmonics n 1 and n 2 , inclusive. In all finite-gradient cases significant modulations appear in the spectrum for higher frequencies, indicating complex plasma dynamics. These modulations produce significant apparent fluctuations in efficiency for particular intervals of the harmonic spectrum, but the large-scale behavior of each of the spectra shown here is a decrease at a similar rate.
To consider a broader range of parameters, we plot η ω [30,50] for both varied gradient scale length with = N 1000 max (Fig. 16c) and varied a 0 /N without a gradient (Fig. 16d). The color scheme for both plots is chosen  Fig. 12c. For frequencies higher than 1000ω L numerical dispersion of the propagating radiation has a measurable widening effect on the duration of the attosecond pulse. Note that for (a-c) we consider light passing the observation point, so the advanced time is equivalent to the lab frame.
www.nature.com/scientificreports www.nature.com/scientificreports/ so that blue regions represent efficiencies greater than 50% of those predicted by the −4/3 power law. Although the zero-gradient case reaches its optimum around a 0 /N ≈ 0.5 and a 0 > 100, the finite gradient case is restricted to lower efficiencies, even for a 0 up to 200. This suggests that although a gradient is useful for increasing efficiency for moderate intensity lasers, the ideal condition for RHHG is a semi-infinite target with density appropriately balanced to the laser field strength.
The trajectories of emitting electron bunches in finite gradients are similar to those for a high-efficiency flat target: as shown in Fig. 17, the width of the emitting electron bunch in the advanced time frame is much longer than the attosecond pulses that are produced for high frequencies. The destructive interference in the trailing edge of the electron bunch means that attosecond pulses are still produced, but the harmonic conversion efficiency drops below the ideal power-law scaling. Provided the gradient is reasonably short (L/λ < 0.2), the corrected synchrotron emission model is still applicable. The most efficient gradients are those with L/λ < 0.2, so the breakdown of the synchrotron model for very long gradients is less important for high-order harmonic generation.
Multi-cycle driving pulses. The emphasis on single-cycle driving pulses in the previous sections has been justified by the fact that the frequency conversion process is intrinsically sub-cycle; the emission occurs within a fraction of the driving laser period. Here we briefly suggest how the single-cycle analysis extends to multi-cycle pulses, arguing that a long driving laser can be represented as a sequence of single cycle interactions provided that the initial conditions of each interaction can be treated appropriately.
The carrier envelope phase (CEP), represented here by φ, describes the alignment of the envelope and underlying phase. We take the convention φ = 0 for a sine-like pulse. Since the generation of attosecond pulses is associated with transitions of the laser electric field during particular optical half-cycles 15 , we can assign each attosecond pulse a phase φ a based on the CEP of the pulse that produced it and the optical cycle from which it was generated. For example, an attosecond pulse generated during the central optical cycle for a φ = 0 incident beam has phase φ a = 0 and additional pulses will be generated during different cycles of the same interaction with φ a = nπ for normal incidence. We can in general write φ a = πn + φ, noting that for oblique incidence, attosecond pulses at odd values of n will be suppressed, a similar role to that fulfilled by the differently-defined ψ g used by Ma et al. 61 .
In the limit of short driving pulses (full-width-half-maximum duration τ), where the parasitic effects of preceding cycles can be neglected, the strength an emitted attosecond pulse will depend on the effective normalized potential ã during the optical half-cycle which leads to its emission. A Gaussian incident laser www.nature.com/scientificreports www.nature.com/scientificreports/ actions with steep surfaces, we may often write -as can be seen in Fig. 5a-I 2 , so that for a Gaussian driving laser the attosecond pulse intensity scales with a 0 , φ a and τ as 2 . This equation predicts the attosecond pulse train shapes in Fig. 18a,b, where the relative intensities of attosecond pulses are presented for −0.8 < t/τ < 0.8. In Fig. 18a,b, the attosecond pulse trains produced by τ = 3.5 fs and τ = 21.2 fs incident pulses agree with the predicted attosecond pulse envelope (red and blue dashed lines), based on q = − 6.6 calculated from a 0 /N scaling for single-cycle pulses under these conditions. Figure 18d summarizes attosecond pulse intensities calculated for a range of incident τ and φ at a 0 = 100, N = 1000, and θ L = 15°, showing that for these parameters the attosecond pulse intensity depends only on position under the envelope φ a /ω L τ, and that the relative strength of pulses in the train matches the expectation from the single-cycle scaling. In Fig.18e-g the attosecond pulse intensity (I a /I L ) is plotted against the envelope position φ a /ω L τ for varied θ L and L/λ, showing a reasonable collapse of η a to a single function of φ a /ω L τ for varied τ across a range of conditions. The second line for L/λ = 0.1 and θ L = 45° are from attosecond pulses emitted on alternative half-cycles from the main pulse.
Generation efficiency becomes more complicated for multi-cylce pulses under the most efficient conditions, where, for example at a 0 /N = 0.25, the maximum attosecond pulse intensity decreases for longer pulse duration (Fig. 18c). Since intensities for different values of τ coincide at early times on our normalized scale, the decrease in maximum intensity for longer τ can be attributed to the parasitic effect of strong early pulses on the achievable intensity in later cycles, through disruption of the plasma surface. This effect is not observed for smaller values of a 0 /N because the plasma response time is much faster than the laser period, so the disruption due to the generation of an attosecond pulse is damped before the next pulse arrives. In contrast, at higher a 0 /N, the relativistic plasma frequency is closer to the laser frequency and plasma disruption may persist for a full optical cycle.
The total spectrum will be the sum of the spectrum contributions the individual pulses. In all multi-cycle pulse cases, a range of intensities under the envelope drive different attosecond pulses, so the individual spectral contributions vary under the envelope. The total spectrum can therefore appear more complex, even if the individual attosecond pulse generation effects are well described by an analytic CSE model.

Conclusion
In this article we have presented results from a large number of particle-in-cell simulations, focusing on how the reflected high-order harmonic spectra varies with plasma density, laser intensity, and angle of incidence, and we have explained key features in terms of a bunch-width corrected coherent synchrotron emission model. Specifically, we have shown that under most conditions the energy spectrum decays as a power of frequency with exponent −p up to a rolloff set by the width of the emitting electron bunch. The spectrum then approximately decreases as ω −p−2 , often with substantial modulation due to the interference between different parts of the electron bunch, until an exponential falloff at a frequency proportional to the cube of the electron Lorentz factor. The initial power-law exponent p appears limited to around 4/3, the prediction of the CSE model, so the maximum harmonic peak powers reachable in the 0.1-10 keV range are likely ω ω .
0 in each harmonic, where P 0 is the initial power and the 0.33 factor arises from normalizing the total reflected power to the incident power. for www.nature.com/scientificreports www.nature.com/scientificreports/ a 10-PW driver at 800 nm, this leads to 13 TW/harmonic at 100 eV, 0.5 TW/harmonic at 1 keV, and 25 GW/harmonic at 10 keV.
Outside of the most efficient conditions, we have found that p is a continuous function of a 0 and N which approaches an asymptotic limit for a 0 → ∞ at fixed a 0 /N and decreases monotonically for a 0 /N < 0.5. The bunch-width cutoff ω b normally lies at smaller frequencies than the Lorentz-factor dependent cutoff ω γ , leading to a steep power-law decay for ω b < ω < ω γ before the exponential falloff above ω γ . Since ω b approximately scales linearly with a 0 rather than with the a 0 3 dependence of ω γ , the increase in the range of high-order harmonics that can be efficiently generated is not as rapid as the naive CSE model predicts. Nonetheless, for experimentally reachable conditions, efficient keV photon production is possible.
Although previous work has limited CSE models to thin foils in reflection 46 or transmission 30,47,62 , or highly-oblique p-polarized incidence 45 , we have found that when the bunch width is treated appropriately, the CSE model is valid across a wide range of conditions and that both ω b and ω γ are important for many reachable parameters. In short, we have generalized the CSE model to explain the evolution of the reflected harmonic spectrum observed in a large number of particle-in-cell simulations spanning a comprehensive set of interaction parameters, developing scaling relationships between the interaction parameters and the efficiency of harmonic and attosecond pulse generation. The scaling laws developed here point to relativistic high-order harmonic generation driven by petawatt-class lasers as the brightest extreme-ultraviolet and soft x-ray source, and competitive with free-electron lasers for peak power at few-keV photon energies.

Methods
Using the PIC codes EPOCH 63 and BOPS 64 , we examined the interaction of single-cycle pulses with fully-ionized overdense (N > 1) semi-infinite plasma at relativistic intensities (a 0 > 1) to find p and ω c . To determine the underlying scaling behavior, target densities and laser intensities were allowed to vary somewhat outside of realistic parameters. Here we neglect quantum electrodynamic (QED) effects, even for light fields with very large nominal a 0 . Although the neglect of QED for large a 0 is not physical, it allows us to determine the "relativistic limit" of a particular set of conditions; if we fix a 0 /N, electron trajectories approach an asymptote as a 0 → ∞. This limit provides an idealized case which more rigorously validates the underlying mechanism. At more realistic conditions, a 0 ≈ 1 effects can obscure the relativistic mechanism and simulations results may not unambiguously distinguish models, even when those models make different predictions. With a mechanism validated by the relativistic limit, www.nature.com/scientificreports www.nature.com/scientificreports/ we can treat more realistic values of a 0 -where there are no QED effects -as small deviations from the ideal relativistic case. The key physics of relativistic RHHG can be captured with a one-dimensional treatment; for non-zero laser angle of incidence (θ L ) a relativistic transform to a boosted reference frame (Bourdier method 65 ) removed dependence on the second spatial dimension. Although the spatial axis is one-dimensional, all three components of the velocity are calculated, and two-dimensional trajectories are found by integration of the transverse velocity components. The simulations presented here used spatial resolutions between 3,000 and 50,000 cells per laser wavelength, with at least 10 cells/wavelength for the shortest harmonic wavelength from which any conclusions were drawn. The number of particles/cell was varied based on noise requirements and ensuring that the critical density was resolved; up to 5000 particles/cell were used. The reflected spectra were found by taking the Fourier transform of the reflected fields. Filtered attosecond pulses were found in a spectral range ω LF < ω < ω UF from the inverse Fourier transform. In all cases, the pulse envelopes were taken as Gaussian.
Calculation of power law exponent. The exponent p of the spectral power law [I(ω) ∝ ω −p ] is found by a linear fit of the PIC-calculated spectrum. The slope of the logarithm of the intensity vs the logarithm of the frequency ( ω = − + I p C log log 10 10 ) is the exponent associated with the power-law. The coefficient of proportionality becomes the additive constant C. The fit is based on a particular interval of the spectrum: e.g. in Fig. 4 p is calculated based on the range 1 < ω/ω L < 100. Since the entire spectrum -not just the harmonic peaks -is used, spectra with distinct harmonic structure will be fit by a line that passes below their harmonic maxima. Note that the range of the fit also means that the fundamental frequency does not dominate or constrain the fit, and in some cases there may be a noticeable discrepancy between the power-law fits and the calculated spectra at ω/ω L ≈ 1.