Identification of coupling mechanisms between ultraintense laser light and dense plasmas

The interaction of intense laser beams with plasmas created on solid targets involves a rich non-linear physics. Because such dense plasmas are reflective for laser light, the coupling with the incident beam occurs within a thin layer at the interface between plasma and vacuum. One of the main paradigms used to understand this coupling, known as Brunel mechanism, is expected to be valid only for very steep plasma surfaces. Despite innumerable studies, its validity range remains uncertain, and the physics involved for smoother plasma-vacuum interfaces is unclear, especially for ultrahigh laser intensities. We report the first comprehensive experimental and numerical study of the laser-plasma coupling mechanisms as a function of the plasma interface steepness, in the relativistic interaction regime. Our results reveal a clear transition from the temporally-periodic Brunel mechanism to a chaotic dynamic associated to stochastic heating. By revealing the key signatures of these two distinct regimes on experimental observables, we provide an important landmark for the interpretation of future experiments.

High-density plasmas can be created by focusing intense laser pulses on initially-solid targets. The interaction of such plasmas with laser light, investigated for several decades, involves a rich non-linear physics which is not only of fundamental interest, but also of high relevance for a wide range of applications over a large interval of laser intensities, and spanning thermonuclear fusion [1], laboratory astrophysics [2], or laser-driven particle acceleration [3,4]. For most if not all applications, depositing laser energy into the plasma is essential. Due to their high density, largely in excess of the so-called critical density where the local electron plasma frequency equals the laser frequency, these plasmas however tend to reflect a large fraction of the laser light. The actual coupling with the incident light field can only occur either in the undercritical part of the density gradient at the plasma-vacuum interface, where the laser wave propagates, or within the skin depth of the overcritical plasma, where the laser wave is evanescent.
At such plasma densities, physicists initially anticipated the main mechanism of energy deposition to be collisional absorption [5]: electron-ion collisions disrupt the regular quivering motion of the plasma electrons in the light field, statistically leading to a net kinetic energy gain from the laser field. Soon after the invention of lasers, experimental studies on the feasibility of laser-driven thermonuclear fusion however revealed the importance of non-collisional light absorption mechanisms, coming into play for moderate laser intensities (Iλ 2 10 13 W µm 2 /cm 2 ) [6]. For interactions at the surface of dense plasmas, these processes are expected to be most relevant when the laser beam impinges the * henri.vincenti@cea.fr † fabien.quere@cea.fr target at oblique incidence (angle of incidence θ i = 0) and in p-polarization, such that an electric field component efficiently drives electron motion along the normal to the target surface. Among these so-called 'anomalous absorption' mechanisms, the resonant excitation and subsequent damping of collective electronic plasma waves at the critical plasma density has been the first key process to be identified both theoretically and experimentally, and is commonly known as resonant absorption [7][8][9][10].
Since the coupling with laser light occurs at the interface of the plasma with vacuum, the characteristic spatial length L of the plasma density gradient across this interface is a crucial parameter. This density gradient is generally not step-like, due to the unavoidable expansion of the plasma into vacuum, either during the main laser pulse driving the interaction, or even before this pulse when laser pre-pulses are present (either accidentally, or voluntarily introduced). In an influential paper from 1987, Brunel predicted a transition from resonant absorption to a new coupling mechanism, that he ironically called 'not-so-resonant, resonant absorption' [11], depending on the value of L. He anticipated this mechanism to come into play when the laser intensity becomes so strong that the quivering amplitude of the plasma electrons in the field along the surface normal gets larger than L, such that the plasma-vacuum interface can be modeled as step-like.
This simple and intuitive mechanism, now known as Brunel absorption or vacuum heating, is qualitatively analogous to the intensively studied three-step model of atomic and molecular strong-field physics, where an intense laser field (Iλ 2 10 13 W µm 2 /cm 2 ) drives the recollision of ionized electrons with their parent ions [12]. Here, in each optical cycle, electrons at the target surface are dragged out of the plasma into vacuum when the component of the laser E-field normal to the surface points to the target (see simulation results in Fig.1). 11 12 13 t/T Temporal dynamics of a dense plasma exposed to an ultraintense laser field in the Brunel regime. This graph displays results from a Particle-In-Cell simulation performed for a0 = 2 (I = 8.5 × 10 18 W/cm 2 for λ = 800 nm), θi = 55 o , and a density gradient scale length L = λ/10 (with λ the laser wavelength). The gray-scale color map shows the temporal evolution, during three laser optical periods T , of the plasma electron density around the target surface, while the purple color scale shows the attosecond light pulses emitted by this plasma (harmonics 8 to 15). Two representative trajectories for the particles forming the expelled electron beam (blue color) and the 'recolliding' electron flux (red color) are also displayed. The right panel shows a typical distribution of electrons in the x − px phase space (x spatial coordinate along the target normal) at time t/T = 12.7.
Later in the cycle, when this driving field changes sign, some of these expelled electrons are pushed back toward the 'parent plasma': as they penetrate into this dense plasma, they escape the influence of the laser field due to plasma screening, and propagate ballistically into the target (red trajectory in Fig.1). The initial model by Brunel for a step-like surface of perfectly-conducting plasma focused on this returning electron population. However, numerical simulations show that for smoother interfaces and/or higher laser intensities, another fraction of the electrons escapes into vacuum (blue trajectory in Fig.1), typically in the form of periodic attosecond bunches. In both cases, the fast electrons resulting from this suboptical-cycle dynamics carry away energy acquired from the laser: for convenience, the terms 'Brunel electrons' or 'Brunel absorption' used in this paper will encompass these two populations.
A few years later, with the development of high-power femtosecond lasers, Brunel mechanism appeared as an ideal 'toy-model' to understand the interaction of dense plasmas with these ultrashort pulses. First, their ultrahigh intensities result in large electron quivering amplitudes. Second, these pulses are so short that plasma expansion during the interaction is very limited, potentially leading to sharp density gradients at the plasma surface, i.e. small values of L. However, from an experimental point of view, reaching this regime turned out to be much more challenging than expected. The key difficulty arose from the unavoidable light pedestal ahead of ultrashort laser pulses [13]: this pedestal, if too intense, leads to the premature creation and expansion of the plasma, and thus to long and largely uncontrolled density gradients at the plasma-vacuum interface when the main laser pulse hits the target, making the Brunel regime inaccessible. More generally, this major issue has considerably complicated the interpretation of most early experiments on the interaction of intense ultrashort lasers with dense plasmas.
It took more than an additional decade to find methods to efficiently reduce the light pedestal accompanying ultrashort laser pulses [14][15][16][17], and thus obtain temporal contrasts that at last made extremely sharp plasma surfaces accessible and compatible with ultrahigh laser intensities [18,19]. Nowadays, Brunel mechanism is most likely at play in experiments performed on solid targets with ultraintense laser pulses of suitably high-contrast. Yet, direct experimental evidence is still elusive, and its range of validity is not precisely known so far. Furthermore, following the historical development of this topic, the 'common wisdom' tends to be that when L is increased, a transition from Brunel to resonant absorption should at some point occur [20][21][22], but no clear experimental evidence of this transition has been reported yet.
A broad range of topical experiments are now performed worldwide on the interaction of ultraintense laser pulses (Iλ 2 > 10 18 W µm 2 /cm 2 ) with dense plasmas, driven by applications such as laser-driven ion [3,4,23,24] and electron acceleration [25][26][27][28][29][30][31][32], or the generation of intense harmonics and/or attosecond light pulses [22,33]. Clearly identifying the laser-plasma coupling mechanisms at play in this interaction regime, and determining the range of physical parameters where they are relevant, is essential for the proper understanding of such experiments. This is what we achieve in this article, by focusing ultraintense femtosecond laser pulses on a dense plasma with a sharp, controlled and measured density gradient scale length L, which we systematically vary from L λ to L ≈ λ (with λ the laser wavelength). We show how performing and correlating measurements of the high-order harmonics and relativistic electrons emerging from the target provide clear signatures of these couplings mechanisms, and relate these observations to the underlying physics through an advanced analysis of 2D and 3D Particle-In-Cell (PIC) simulations, solving the coupled Vlasov-Maxwell equation system. This comprehensive study shows that Brunel mechanism is indeed the relevant physical process for sharp enough plasma-vacuum interfaces. As expected, a transition occurs to a different mechanism when the density gradient scale length L is increased. Measurements of this transition as a function of the laser incidence angle provide confirmation of Brunel's transition criteria based on the comparison of the electron quivering amplitude with the typical spatial extent of the interface. However, we establish that in the regime of ultrahigh laser intensities considered here, resonant absorption plays no significant role in the regime of large L (L ≈ λ). The coupling is rather dominated by another kinetic mechanism, so far known as stochastic heating, in which collective plasma a L ~ λ / 1.5 L ~ λ / 15 The target is impinged by a controlled prepulse, followed by the main pulse after an adjustable delay τ . Two of the main diagnostics are displayed, the lanex screen for the measurement of the spatial profile of the high-energy electron beam, and the angularly-resolved xuv spectrometer. These two diagnostics can be used either separately, or simultaneously when small holes are made in the electron detection assembly, as shown in the figure. They can also be replaced by an angularly-resolved electron spectrometer. The main experimental findings for a p-polarized laser field are summarized in the two lower panels b and c: left images, angular emission pattern of relativistic electrons; bottom right images, angularly-resolved energy spectrum of electrons in the incidence plane (θy = 0); top right images, angularly-resolved harmonic spectrum. These illustrate the major changes occurring on these three observables as the density gradient scale length is increased from L λ to L ∼ λ (θi = 55 o , a0 = 3.5, τ = 1 ps for panel b, leading to L = λ/15, τ = 10 ps for panel c, leading to L = λ/1.5). These very contrasted features constitute signatures of the different underlying laser-plasma coupling mechanisms. effects play little role: as initially suggested in Ref. [34], electrons in the underdense part of the density gradient gain energy in the interference pattern resulting from the superposition of the incident laser field with the field reflected by the overdense part of the plasma. It has been established theoretically that at the laser intensities considered here, electron dynamics in such an interference pattern is not integrable, gets chaotic, and can lead to high energy transfer from the laser wave to the electron population [35].
In this paper, the amplitude of the incident laser field, which determines the intensity on target, is characterized by the dimensionless potential vector at the peak of the pulse, a 0 = eE 0 /mcω = λ[µm](I[W/cm 2 ]/1.37 × 10 18 ) 1/2 , with c the speed of light, e the elementary charge, m the electron mass, ω the laser frequency, and E 0 the amplitude of the laser electric field. All experiments and simulations presented here have been performed with a 0 > 1, which corresponds to the interaction regime where relativistic effects play an important role on electron motion. We define n c as the critical plasma density associated to the laser frequency ω (n c = m 0 ω 2 /e 2 = 1.74 × 10 21 cm −3 for 800 nm laser light).

I. DESCRIPTION OF THE EXPERIMENT
A sketch of the experiment is presented in the upper part of Fig.2. A high-power femtosecond laser beam is focused on a silica target, which it fully ionizes on a thin surface layer, thus producing a dense plasma (maximum plasma density n 0 6.10 23 cm −3 , i.e 400n c for 800 nm light).
We use the UHI100 laser at CEA Saclay, a commercial system delivering 20 to 25 fs pulses with a peak power of 100 TW. After correction of its wavefront by an adaptive optic system, the beam is focused by an off-axis parabolic mirror with a f -number of f /6, leading to a focal spot of 5 µm diameter (FWHM in intensity) and to an estimated peak intensity of 2.10 19 W/cm 2 (a 0 ≈ 3.5) on target. By default, the laser beam is p-polarized on target, but the polarization can be switched to s by inserting a thin zeroorder mica half wave-plate in the beam.
The first key aspect of the experiment is that it was carried out with high degree of control and an accurate knowledge of the plasma density gradient scale length L at the target surface, which is a prerequisite for the study of the laser-plasma coupling mechanisms. This implies the use of laser pulses of very high temporal contrast, so that premature creation of the plasma on target is avoided: this is achieved thanks to a double plasma mirror system placed before the main experimental chamber [18], which increases the contrast by about 4 decades, from 10 9 to 10 13 on a 100 ps time scale. This ultrahigh temporal contrast is of paramount importance for all experiments presented here. Starting from the very steep density gradient allowed by this ultrahigh temporal contrast, L is then varied in a controlled way, thanks to the introduction of a 'weak' prepulse (fluence ≈ 1 kJ/cm 2 ) at an adjustable delay τ before the main pulse (0 ≤ τ ≤ 15 ps), produced from an edge of the main beam using the optical layout described in Ref. [36]. This prepulse is strong enough to ionize the target and initiate plasma expansion, at a typical velocity in the 40 to 60 nm/ps range. For all experimental conditions considered herein, and in particular for all incidence angles, we systematically measured L(τ ) using the recently-introduced technique of Spatial Domain Interferometry [37] (SDI, see supplementary material).
The second key aspect of the experiment is the combination of diagnostics that were implemented to study the interaction. We concentrated on two types of observables: the relativistic electron beam emitted by the target towards vacuum, and the beam of high-order harmonics generated around the specular reflection direction.
Two diagnostics were used for the electron beam. First, a lanex screen, placed behind a 13 µm thick aluminum foil (to eliminate laser light and its harmonics) and a 2 mm thick glass plate (to filter out low-energy electrons), and which fluorescence was imaged on a CCD camera, provided the spatial profile of the emission of electrons with energies higher than 1 MeV, at a distance of ≈ 10 cm from the target (left images in Fig.2b-c). Second, we designed a new type of magnetic spectrometer for relativistic electrons (see supplementary material), which provided, for every laser shot, the angularly-resolved kinetic energy spectrum of electrons, in the incidence plane (i.e. for θ y = 0). This (θ x , E) distribution, with θ x the angle in the plane of incidence and E the electron kinetic energy, was measured with a very large angular acceptance of ∆θ x = 500 mrad around the specular reflection direction (θ x = 0) (bottom right images in Fig.2b-c).
The harmonic beam was characterized using an angularly-resoled XUV spectrometer [36], with an angular acceptance of 200 mrad around the specular direction (top right images in Fig.2b-c). The harmonic spectrum and the electron beam spatial profile were initially measured on the same laser shots, thanks to small holes in the aluminum foil, glass plate and lanex screen that let the harmonic beam go through (Fig.2). However, we observed an excellent shot-to-shot reproducibility of the experimental results, so that these multiple diagnostics were finally implemented on different laser shots performed under identical interaction conditions. A simple additional diagnostic, implemented on separate laser shots, consisted in measuring the spatial profile of the laser beam reflected by the target, by inserting a frosted glass plate in this beam 20 to 30 cm after the target, and measuring the image of the laser light scattered by this plate on a camera placed behind a bandpass filter centered on the fundamental laser frequency. This can be exploited to determine the plasma reflectivity for the incident laser beam.

II. EXPERIMENTAL RESULTS
The lower part of Fig.2 summarizes the main finding of the experiment for a p polarization of the incident laser, by presenting the electron beam angular profiles and (θ x , E) distributions, as well as the angularly-resolved harmonic spectra, measured for two different density gradient scale lengths L, L 1 λ and L 2 ≈ λ: as L is increased from L 1 to L 2 , the measured signals for all these observables radically change. Three main differences are observed: (i) When L λ, the electron emission is predominantly peaked at θ x ≈ 100 mrad, i.e. close to the direction of laser specular reflection (θ x = 0 mrad), with a slight shift towards the target normal. As L is increased, it then switches to θ x ≈ −200 mrad, a direction between specular direction and the tangent to the target surface (θ x = −600 mrad), and simultaneously slightly broadens angularly. (ii) Electrons reach energies about twice higher for large L (spectral peak around 10 MeV), with a (θ x , E) distribution that significantly changes. In the short gradient regime, a clear correlation is observed between emission angle θ x and electron energy E, especially in the most intense part of the distribution (0 ≤ θ x ≤ 200 mrad): the electron energy increases as one gets closer to the specular direction. In contrast, in the long gradient regime, the electron spectrum hardly varies angularly, i.e. no significant correlation is observed on this (θ x , E) distribution. (iii) Harmonic emission is clearly observed for small L, but it drops below the experimental detection threshold for large L.
The details of the transition between these two regimes are presented in Fig.3, which displays the evolution with L of the electron beam angular profile in the incidence plane, and of the harmonic spectrum. The most important point is the quantitative correlation, observed at short gradients, between the emission of relativistic electrons and the harmonic signal (see curves in Fig.3c). As L is gradually increased, the electron signal around θ x = 100 mrad and the harmonic signal reach a common optimum around L = λ/15, and then both quickly decrease. The electron signal on the other side of the specular direction then grows, but is not associated with any harmonic signal. The transition between these two regimes occurs around L ≈ λ/5. Another important difference between these two interaction regimes is the dependence of the observables on laser polarization direction, illustrated in Fig.4. In the short gradient regime, the electron and harmonic signals are totally suppressed when the polarization is switched from p to s. By contrast, for longer gradients, the electron signal is still observed for s-polarization, although it gets about five times weaker.
These observations on the electron and harmonic beams clearly point to a complete change in the coupling mechanism between the laser field and the plasma, which we will analyze in the rest of this article. In the following, we will refer to these two interaction regimes as the short-gradient and long-gradient regimes for convenience.
This transition also has consequences on even simpler observables: Fig.5 (left panels) displays the spatial intensity profiles of the reflected laser beam, measured on a scattering screen in these two distinct coupling regimes. In the short gradient regime, a smooth beam is observed, which is almost unaltered compared to the incident laser beam: this is the so-called plasma mirror regime [19], where the plasma acts as a usual high-quality mirror, specularly reflecting the fundamental frequency, despite the ultrahigh intensity on target. By contrast, in the long gradient regime, the beam profile is strongly perturbed and starts exhibiting spatial structures that were not present on the incident beam. The term plasma mirror is thus no longer appropriate to this regime, although the laser field still interacts with a dense -and hence reflective-plasma. Experimentally, the spatial profile of the reflected laser beam might thus also be used as an alternative and very simple signature of the transition in the laser-plasma interaction. By spatially integrating these images, the variation of the plasma reflectivity at the fundamental laser frequency as a function of L can be determined, and is displayed in the right panel of Fig.5, for both p and s polarizations of the incident laser field. All these measurements have been repeated over four different experimental campaigns on the UHI100 experimental facility over the last four years, and all effects described above have been observed to be very reproducible and robust.

III. 3D PARTICLE-IN-CELL SIMULATIONS
To interpret these experimental observations, we will now turn to Particle-In-Cell (PIC) simulations of the laser-plasma interaction. This requires ensuring that (i) these PIC simulations are performed in the actual physical conditions of the experiment, (ii) they are reliable and properly reproduce the key experimental findings. To check these two critical points, we first carry out full 3D simulations of the interaction, so that we can directly confront the numerical and experimental results -especially the full angular pattern of the electron emission, which is only accessible by 3D simulations.
All simulations reported herein have been performed using the recently-developed WARP+PXR code [38][39][40][41][42]. The specificity of this code is the use of a massivelyparallel high-order spectral solver for Maxwell's equations [43], which greatly reduces numerical dispersion of electromagnetic waves as well as numerical noise. This ensures convergence of the simulations for larger spatial and temporal mesh steps than in most other PIC codes, and thus makes physically-realistic and reliable 3D simulations of the interaction with dense plasmas compu-tationally tractable [43,44]. Each 3D simulations reported here required 6.3 millions computation hours on a massively-parallel machine [45,46]. More detailed information on the numerical parameters of these simulations is provided in the supplementary material.
We performed two 3D simulations, for the same physical conditions as in the experiments of Fig.2b and c, corresponding to fixed laser parameters but different density gradients scale lengths, L 1 λ and L 2 ≈ λ. From these simulations, we extracted the exact same observables as those measured in the experiment: the calculated angular profiles and angle-energy distributions of the emitted electron beam, as well as the angularly-resolved harmonic spectra, are displayed in Fig.6. Comparison with the lower panels of Fig.2 shows that these simulations very well reproduce the two distinct interaction regimes observed experimentally for all these observables.
These 3D benchmark simulations clearly demonstrate both the reliability of the PIC simulations, as well as our excellent control of the interaction conditions in the experiment, which ensures a relevant choice of the physical parameters used in the theoretical study. In the next section, we will further exploit simulations performed with this WARP+PXR code to get detailed insight into the physical processes underlying these distinct interaction regimes. To this end, we will rely on more tractable 2D simulations [47]. Detailed information on the numerical parameters of these simulations is also provided in the supplementary material.

IV. PHYSICAL ANALYSIS AND DISCUSSION
The starting point for our analysis of the experimental results is the joint measurement of the relativistic electron and harmonic signals, which provides information on the temporal structure of the electron emission by the plasma (section IV A). We then discuss the spatial structure of the observed electron beams (section IV B). In order to understand the electron heating and ejection mechanism in the long gradient regime, we turn to a detailed analysis of PIC simulations, presented in section IV C. With the physical insight provided by this analysis, we finally discuss the influence of the laser polarization (Fig.4) and the evolution with L of the plasma reflectivity at the fundamental laser frequency (Fig.5) in section IV D.
A. Temporal structure of the electron emission The clear correlation observed between the high-energy electron signal and the harmonic signal for short gradients (Fig.3c) suggests that in this regime, the relativistic electrons are involved in the harmonic emission. Then, the fact that a highly-contrasted harmonic comb is produced would be an indication that this electron emission is periodic in time, being locked to the driving laser field. To support this tentative interpretation, we consider the PIC simulations of Fig.1 and Fig.7b: they show that temporal periodicity is indeed a characteristic of the Brunel mechanism. Electron emission occurs in the form of bunches that are initially extremely short (in the attosecond range), emitted once every optical period. For a 0 > 1, these electrons reach relativistic velocities when they escape the plasma, and thus induce a Doppler effect on the reflected laser field: this results in the generation of a train of attosecond light pulses, spaced by one laser period, which are clearly observed on Fig.1 and Fig.7. This well-identified process is known as the Relativistic Oscillating Mirror (ROM) effect [19,[48][49][50][51][52][53][54][55]. In simulations, the resulting periodic light emission has a spectrum consisting of a comb of high-order harmonics of the laser frequency: this is the origin of the harmonic signal observed in our experiment. The electron-harmonic correlations observed experimentally is therefore a signature of the periodicity of the dynamics of electrons.
As described in section II, the harmonic signal is observed to collapse for longer density gradients L. One possible interpretation for this collapse can be that the electron emission ceases to be periodic in time. And indeed, PIC simulations for longer density gradients (Fig.7d) strikingly show that, in contrast to the Brunel mechanism, electron emission is no longer periodic in this regime. The absence of harmonic signal in conjunction with the relativistic electron emission can thus be considered as a signature of the transition to a new coupling mechanism, associated to the very different plasma temporal dynamics observed in Fig.7. This mechanism will be described in section IV C.

B. Spatial structure of the electron emission
We now discuss the spatial properties of the outgoing electron beams, with the support of the simulation results of Fig.7. We show that in the short gradient regime, this structure is mostly determined by the interaction of expelled electrons with the reflected laser field in vacuum (section IV B 1), while in the long gradient regime, it is rather imposed by large quasi-static surface fields that develop in the vicinity of the plasma surface during the interaction (section IV B 2).

Short-gradient regime
In the case of short density gradients, the peculiar angular structure of the electron beam has recently been analyzed experimentally and theoretically in Ref. [32]. In this Brunel regime, electrons are expelled from the plasma as a very laminar beam, with relativistic velocities initially quasi-parallel to the direction of specular reflection (Fig.7a). These relativistic electrons thus co-propagate with the intense reflected laser field, with which they interact in vacuum over a distance of the order of the Rayleigh length. This interaction always results in the ejection of electrons out of the laser beam, and therefore digs a hole in the electron beam, centered on the specular direction, as observed in our experiment.
There are two typical scenarios for this ejection, depending on the exact initial conditions with which electrons are expelled from the plasma into vacuum. Some electrons explore multiple optical cycles of the laser field, and thus oscillate in the field and get expelled from the laser focal volume by the so-called ponderomotive effect, isotropically and with a limited energy gain [56,57]. They form the ring-shaped halo observed on the electron beam. But most electrons actually remain around a given phase of the reflected field and rather 'surf' a single wavefront of the reflected field, thus escaping the laser beam laterally along the laser polarization direction, and forming the bright peak observed next to the specular direction. The side on which this peak forms is determined by the laser phase at which electrons are expelled from the plasma into vacuum: the observation of a peak on one side only of the 'ponderomotive hole' (between the specular direction and the target normal) is a clear indication that electrons are ejected periodically once every laser period, when the laser field drags them out of the plasma, in the form of a sub-optical cycle bunches.
This second set of electrons experiences a quasiconstant E-field from the laser in vacuum until they escape the focal volume, leading to a greater energy gain than in ponderomotive scattering: this 'vacuum laser acceleration' (VLA) process accounts for the observed asymmetry of the (θ x , E) distribution (see Fig.2b), where higher energies are observed on one side of the hole (VLA electrons) than on the other ('ponderomotive' electrons), as well as for the angle-energy correlations on this distribution [58]. An important consequence is that the Doppler upshift factor induced by outgoing electrons on the reflected field, which leads to the generation of high-order harmonic (ROM mechanism), cannot be directly deduced from the electron spectra measured experimentally, since electrons keep gaining energy after they escaped the target and emitted high-order harmonics, and before being detected. For instance, simulations show that in the present experiment, the electron Lorentz factor typically varies from γ ≈ 2 − 3 as they are ejected from the plasma and emit harmonics, to γ 15 after their interaction with the reflected laser field [32], when they are detected.

Long-gradient regime
The spatial properties of the electron beam observed in the 'long-gradient' regime described here have not been explained in detail yet, to the best of our knowledge. The electron trajectories displayed in Fig.7 show that the conditions of electron ejection from the plasma are already very different in the short and long gradient regimes. In the second case, the expelled electron beam is no longer laminar, and rather has a complex velocity distribution. As will be shown in the next section, this feature can be attributed to the chaotic character of the electron heating mechanism leading to ejection from the plasma.
Furthermore, these 2D PIC simulations show that a quasi-static magnetic field develops at the plasma surface (see map of the magnetic field in Fig.7c, where a red area indicating a high magnetic field is present close to the plasma surface). This field grows during the laser-plasma interaction, reaches an amplitude typically of the same order of magnitude as the laser magnetic field, and then persists even after the laser pulse has been reflected by the plasma. A detailed analysis of electron trajectories in these simulations shows that this surface field, which is much larger than the one occurring in the short gradient regime (compare maps of the magnetic field in Fig.7a and  c), strongly deflects the escaping electrons toward the target surface. This deflection accounts for the fact that the electron angular distributions observed in this regime are essentially centered between the specular direction and target surface. Test simulations have been performed to check that the reflected laser field plays no role in the deflection of electrons after their ejection from the plasma, in strong contrast to what is observed in Brunel regime.
Such surface quasi-static fields have already been reported in multiple studies of the interaction of intense lasers with dense plasmas (see e.g. [59][60][61][62]), and can be induced by a variety of physical processes (see e.g. [63][64][65][66]). In the present case, our simulations indicate that their development can be attributed to the 'fountain effect' described in [66] where they originate from the plasma cold return current that compensates for the lateral charge ejection from the laser focal volume. This is supported by the fact that these fields are not observed at all in plane wave simulations, where the plasma surface is homogeneously illuminated by the laser field.
In contrast to the short-gradient regime, here the spatial properties of the electron beam do not provide much insight into the involved electron heating mechanism. To identify this mechanism, we now turn to a more detailed numerical investigation based on PIC simulations.
C. Electron heating mechanism in the long gradient regime

Importance of the reflected laser field
Our analysis of the electron heating mechanism in the long gradient regime is based on a set of 2D plane wave PIC simulations [47] for three different physical configurations. Their key results are summarized in Fig.8, and shed light on the underlying physical mechanism. The upper panels (case A) correspond to the interaction of an ultraintense laser pulse with a dense plasma in the long gradient regime, i.e. the same physical configuration as in Fig.7 and as in our experiment: they display the temporal evolution of the plasma electron density (panel a), and the x−p x phase space distribution of electrons (panel d) at the time when electron ejection from the plasma is observed to start (blue dashed line in panel a).
The middle panels (case B) display the same quantities, but now in a situation where the plasma profile has been truncated for densities n > 0.4n c cos 2 θ, i.e. keeping only the underdense part of the plasma, such that there is hardly any reflection of the incident laser by the plasma. This underdense plasma layer, surrounded by vacuum on both sides, is irradiated by the same laser beam as before, but also by a second beam of slightly lower intensity (80%), symmetrically arriving from the other side of the plasma. The role of this second laser is to emulate the beam reflected by the dense part of the plasma in case A.
The key point here is that both the temporal dynamics of the plasma density profile and the electron phase-space distributions look very similar for cases A and B.
By contrast, if the truncated plasma layer is irradiated by one laser beam only (case C, lower panels), the plasma dynamics becomes totally different. More specifically, while similar upward electron ejections are observed in cases A (corresponding to the electron signal observed in our experiment) and B, this electron emission is strongly reduced in case C: less electrons are emitted, and they have much weaker velocities. This toy-model study leads to two important conclusions. (i) The comparison of cases A and B indicates that in the long gradient regime, the coupling mechanism leading to electron ejection mostly occurs in the underdense part of the density gradient. (ii) The comparison of cases B and C indicates that the overdense part of the plasma nonetheless plays a key role, by producing a reflected beam that, when crossing and interfering with the incident beam, strongly modifies the dynamics of electrons in the underdense plasma layer. As we explain in the next section, the electron heating mechanism coming into play in such a physical situation is already well-identified in the existing literature.

Introduction to stochastic heating
Let us start by considering a free electron (i.e. without any collective plasma effects involved) exposed to two non-collinear ultraintense laser beams. Using a quantum description of the field as an ensemble of photons provides a simple way to understand that this electron can gain more energy than when exposed to a single laser beam. When a single laser beam (assumed to be a plane wave) is present, it is well-known that photon absorption processes are hindered because they do not enable to conserve both energy and momentum of the total system. By contrast, when two non-collinear beams are present, the combined absorption of multiple photons simultaneously from both beams is allowed, because the availability of photons with different k vectors makes it possible to conserve both energy and momentum of the total system. In other words, the presence of a second beam allows for energy absorption by the electron from the laser field, through a process that can be defined as stimulated multiphoton Compton scattering [67].
For large field amplitudes, the laser field can be treated classically, and many previous studies in the literature have shown that electron dynamics in these combined . case A (upper line) corresponds to an overdense plasma with a density gradient of scale length L = λ, irradiated by a single laser beam with an incidence angle θi = 55 o and a0 = 2.5. In case B, this plasma has been truncated for n > 0.4nc cos 2 θ, keeping only the underdense part of the plasma, which is now irradiated by two almost identical laser beams (same parameters as in case A) arriving symmetrically from both sides of the plasma layer. In case C, this same plasma layer is irradiated only by the upper laser beam. In these three cases, plots of the temporal evolution of the plasma electron density (panels a-c) and snapshots of the x − px phase space distribution of electrons (panels d-f, taken at the time indicated by the blue line in panels a-c) are displayed. The multilayered phase space distributions of panels d and e strongly contrast with the smooth regular distribution of panel f, and are typical of chaotic dynamics resulting from a repetitive stretching and folding effect in phase space.
noncollinear fields is not integrable and gets chaotic for high enough laser amplitudes (typically a 0 0.15 for at least one of the two beams). This results in large energy gains, and this efficient regime of energy absorption by electrons is known as stochastic heating [35,[68][69][70][71] -although the name chaotic heating would probably be more appropriate here, since the system is perfectly de-terministic and involves no stochastic process.
This effect is obviously not restricted to isolated free electrons: it can equally occur for electrons in an underdense plasmas, leading to an energy absorption process where neither collisions nor collective plasma effects play any major role. Such a coupling of the plasma with multiple laser beams has been studied experimentally in [72] by exposing an underdense plasma to two laser beams (like in case B of Fig.8). It is also known to play a role in electron injection in laser-driven plasma wakefield accelerators by the colliding pulse scheme [73,74].
To the best of our knowledge, Ref. [34] was the first to point out that this mechanism should also come into play when an small underdense plasma layer is present at the surface of a dense plasma exposed to a single laser beam. In this case, the required second non-collinear laser beam results from the reflection of the single input beam by the dense plasma. Electrons in the underdense plasma are then exposed to the standing wave formed in front of the dense plasma by the superposition of the incident and reflected beams, and can gain energy by stochastic heating: this is precisely how we interpret our present experimental results in the 'long' gradient regime.

Numerical evidence for stochastic heating in the long gradient regime
With the help of simulations, we now support this interpretation by providing evidence that stochastic heating indeed occurs in the long gradient regime. To this end, we analyze the electron phase space distributions in cases A, B and C (panels d-f in Fig.8) and their temporal evolution during the laser-plasma interaction (see movie M1). When a single laser beam is present (case C), electrons are observed to simply oscillate non-linearly in the laser field, leading to a smooth and regular phase space distribution (panel f and movie M1). In striking contrast, in cases A and B, electron dynamics in the standing wave resulting from the superposition of two non-collinear laser beams is complex: the key point is that we observe a very strong local 'stretching and folding' effect on the phase-space distribution, around each node of the standing wave (see movie M1).
Such a stretching and folding effect results in very different trajectories for particles that are initially very close in phase space: this is known to be one of the most typical routes to chaotic dynamic [75], exemplified in the wellknown horseshoe map models. These repetitive stretching and folding eventually result in a highly-structured, multilayered phase space distribution (panels d-e), where electrons at a given spatial position have a complex momentum distribution, typical of chaotic dynamics. The striking contrast between these highly-structured phase space distributions, and the smooth distribution observed in case C (Fig.8f)  of these phase space distributions with the type of distribution observed for the Brunel mechanism (Fig.1) is illustrative of the major difference in the dynamics of the system between the short and long gradient regimes.
The chaotic character of the electron dynamics can be further supported by the calculation of the Lyapunov exponents of plasma electrons (see supplementary material), which should be positive in the case of chaotic dynamics. The Lyapunov exponent σ px for the p x variable, obtained from 2D plane wave PIC simulations of case A of Fig.8, are displayed in Fig.9, as a function of the incident laser amplitude. This exponent is negative at low intensity, and gets positive when a 0 0.15, thus pointing to chaotic dynamics. This threshold in laser intensity is consistent with early theoretical investigations of stochastic heating [35,68].

D. Effect of the polarization, plasma reflectivity
With the support of the previous physical analysis, we finally discuss the experimental observations on the influence of the laser polarization direction (Fig.4), and the evolution of the plasma reflectivity with L, for both s and p polarizations (Fig.5).
In the short gradient regime, switching the polarization from p to s is expected to suppress both electron emission and harmonic generation, as observed experimentally ( Fig.4a-b), because the laser E-field component normal to the plasma surface is the main driving force of Brunel absorption at the laser intensities considered here. In contrast, when stochastic heating is involved, the polarization direction of the incident laser beam is not expected to have a strong influence, since the plasma surface only comes into play by producing a reflected wave: this explains why, in the long gradient regime, the electron signal is experimentally observed to still persist in s-polarization (Fig.4c-d). This experimental observation alone makes a strong case against an interpretation of the long gradient regime in terms of resonant absorption, which should rather be very sensitive to the laser polarization.
In the long gradient regime, a spatial degradation of the reflected laser beam wavefronts is observed right after the target in these simulations (Fig.7c), while the beam wavefront is preserved in the short gradient regime (Fig.7a). This is qualitatively consistent with experimental observations (Fig.5), where the laser beam intensity profile far from the target is observed to become degraded for long gradients. The chaotic character of the electron dynamics, identified in the previous section, affects the laser beam propagation in the underdense plasma layer, and thus provides a possible interpretation for this degradation of the reflected laser wavefronts.
From a more quantitative point of view, 2D PIC simulations also properly reproduce the evolution with L of the plasma reflectivity at the fundamental laser frequency, for both s and p polarizations (Fig.5). For short density gradients, the reflectivity is lower in p (≈ 50%) than in s (close to 100%) polarization, while similar values (≈ 70%) are observed for both polarizations in the long gradient regime. This is consistent with the results of Fig. 4 on relativistic electron emission, and with our previous interpretation of the inter-action: in s-polarization, for short L, Brunel absorption is suppressed, while for long L stochastic heating is still present.
These 2D PIC simulations can be exploited to get more insight into the redistribution of the initial laser energy after the laser-plasma interaction. The different distributions obtained for both short and long gradient regimes, and for both p and s laser polarizations, are displayed as pie charts in Fig. 10. For short gradients and ppolarization, around 25 % of the laser is converted into harmonics of the laser frequency (mostly low-orders), and about one third is deposited as kinetic energy of the plasma particles (among which the relativistic electrons observed in our experiment). When the polarization is switched to s, these two contributions get considerably reduced, down to around 5% each, leading to a much higher reflectivity of the fundamental laser frequency. For long gradients, by contrast, very little energy is converted into harmonics, regardless of polarization. As already emphasized, the energy stored into quasi-static fields around the plasma surface significantly increases compared to the short gradient regime. In p-polarization, the fraction of energy going into particle kinetic energy is only slightly weaker than in the short gradient case, and gets reduced by about 50% in s-polarization. The reflectivity for the fundamental frequency remains similar for both polarizations, while more energy goes into scattered light in s-polarization.

E. Conclusion of the physical analysis
The previous combination of multiple experimental observables and PIC simulations has provided strong evidence for the transition from Brunel absorption to stochastic heating as the density gradient scale length L is increased, while no evidence of resonance absorption has been observed in this ultrahigh laser intensity regime.
The intuitive physical insight underlying this transition is that the Brunel mechanism requires a sharp interface, such that the amplitude of the quivering motion of electrons in the laser field exceeds the length scale of this interface. On the opposite, stochastic heating can only occur if electrons are present within the standing wave interference pattern formed by the incident and reflected fields: it will therefore be favored by longer density gradients, for which this interference zone contains a larger number of particles provided by the underdense part of the plasma located just in front of the laser reflection point.
It however remains difficult to predict analytically the value L t of the density gradient scale length for which this transition occurs. In the next section, we determine this transition length experimentally as a function of the laser incidence angle on target θ i -an essential physical parameter, since it affects the point of the density gradient where laser reflection occurs. We also discuss the effect of this angle on the properties of the emitted electron beams.

A. Evolution of the transition gradient scale length
We have repeated the previous measurements, carried out for θ i = 55 o , for five other incidence angles θ i ranging from 40 o to 65 o . For each angle, the density gradient scale length was systematically varied by changing the prepulse delay, and was measured using the SDI technique. The main outcomes of this experiment are summarized in Fig.11. The left panels show how the electron beam angular profile in the incidence plane evolves as a function of L, for three different incidence angles. In all cases, the same type of transition as reported in Fig.2 and  3 is observed, and it occurs for shorter values of L as the incidence angle is increased. As before, we observe that this transition is correlated with major changes in the electron angle-energy distribution and in the harmonic emission. The experimental transition length L e t deduced from all these data is plotted in Fig.11b (blue dots), and clearly decreases with θ i .
To get some qualitative insight into this angular dependance, we consider the starting point of Brunel's model [11], which is that the physics of the laser-plasma coupling changes when the quivering amplitude ∆z e of plasma electrons starts exceeding the typical spatial extent d of the plasma-vacuum interface. We then definesomewhat arbitrarily-this spatial extent as d = |z s − z c |, the distance between the effective reflective surface of the plasma, located at z s (θ i ) such that n(z s ) = n c cos 2 θ i , and the fixed point corresponding to the location of the critical plasma density, z c such that n(z c ) = n c (see sketch in Fig.11b). d can be calculated by using the standard assumption of an exponential density profile at the target surface [76]: n(z) = n 0 exp(−z/L), with n 0 the maximum plasma density in the target, reached in z = 0. This leads to d = −2L ln(cos θ i ).
In the simplest possible approach, the electron quivering amplitude can be considered to be the one of free electrons in the laser field, which is ∆z e ≈ λ/2π×a 2 0 /(a 2 0 +1). At sufficiently high intensity, well into the relativistic regime (a 0 1), ∆z e ≈ λ/2π is quasi-independent of the laser field amplitude. The condition ∆z e = d(θ i ) then leads to a simple theoretical expression for the transition gradient scale length (see movie M1): This simple analysis indeed predicts an angular dependence of L t on θ i : physically, this is because for a given plasma density profile at the surface (i.e. a fixed L), the distance d increases with θ i , due to the wellknown angular dependence of the effective critical den-  11. Effect of the laser incidence angle on the laser-plasma coupling. Panel a shows the evolutions with L of the high-energy electron beam angular profile in the incidence plane, measured for three different incidence angles, for a laser field amplitude a0 ≈ 3.5. For each angle, L is measured experimentally using the SDI technique. The value L e t where a transition occurs is clearly observed to decrease with θi. The measured evolution of L e t with θi is plotted in panel b, where it is compared to the prediction of the simple model described in the text. The predictions of this model are shown for a range of laser amplitudes 2 ≤ a0 ≤ 5, to account for the experimental uncertainty on this amplitude. These predictions are very close to the high-amplitude limit given by Eq.(1). sity, n s = n c cos 2 θ i . Panel b in Fig.11 displays a quantitative comparison of the experimental values L e t with the prediction L t of this model, showing a good agreement. Despite its extreme simplicity, this model can thus account for the angular dependence observed experimentally. We note however that it becomes inappropriate for small angles -a regime that could not be investigated in our experiment for practical reasons-since it predicts that L t −→ ∞ as θ i −→ 0.

B. Evolution of the electron number and energy
Several other effects are observed when the incidence angle is changed, depending on the coupling mechanism. This is summarized in Fig.12, where the evolutions of the electron signal and average energy are displayed as a function of the gradient scale length L and incidence angle θ i , over a range that covers the short and longgradient interaction regimes described before. In these plots, the areas associated to these two regimes are separated by a white zone, which corresponds to the transition curve predicted by the simple model of the previous section, and displayed in Fig.11b. In the following, we describe these evolutions and suggest tentative interpretations, which will need further investigations to be validated in detail.
Both in the Brunel absorption and stochastic heating regimes, the number of emitted electrons grows with θ i (see Fig.12a) in the angular range investigated here. This might be attributed to the fact that when θ i is increased, the target area covered by the laser focal spots increases, so that more electrons get involved in the interaction.
As far as energy is concerned, it is hardly affected by the gradient scale length L in each interaction regime, while the effect of incidence angle is different for the these two regimes. In the Brunel regime, the electron average energy only weakly changes with incidence angle, slightly decreasing for larger angles. Physically, in this regime, most of the electron energy gain is due to the VLA process, which is not expected to depend on the incidence angle since it occurs in vacuum. Changing θ i might however modify the conditions of injection of electrons in the reflected field (e.g. their initial energy), and this might explain the observed slight angular dependence of the energy of VLA electrons.
By contrast, in the stochastic heating regime, the electron spectrum clearly shifts to higher energies as θ i is increased (see Fig.12b), up to about 20 M eV . This might be due to the fact that, all laser parameters being kept fixed, the 'life time' of the interference pattern formed by the incident and reflected fields increases with θ i . A simple analytical calculation indeed gives the following equation for the duration τ i of this standing wave: where τ L is the laser pulse duration, and w the laser focal spot size. As a result, for larger θ i , the stochastic heating process driven by the interference pattern can last longer, leading to a larger final energy gain for electrons, which could account for the experimental observation of Here again, the electron spectra used for this calculation were selected on the right of the ponderomotive hole for short gradients, and on its left for long gradients.

VI. CONCLUSION
In conclusion, we have combined state-of-the-art experiments and PIC simulations to provide the first unambiguous experimental evidence of the transition from Brunel absorption to a different laser-plasma coupling mechanism, as the steepness of the plasma surface is varied. This mechanism has been identified to be stochastic heating of electrons in the underdense plasma layer at the target surface, driven by the standing wave formed by the incident and reflected laser waves. This work has enabled the identification of clear signatures of these two coupling mechanisms, carried by the relativistic electron emission towards vacuum, the generated harmonic signal, and even the spatial profile of the reflected laser beam. At the laser intensities considered here, no evidence of the process known as resonant absorption has been found.
These signatures should prove extremely useful for the interpretation of a broad range of topical experiments performed with high-power ultrashort lasers, related e.g. to ion and electron acceleration, or to the generation of short-wavelength light and attosecond pulses. In the later case, ultrahigh contrast pulses are required, generally obtained with plasma mirrors or by frequencydoubling after the final compression stage, and our results confirm that such experiments involve the Brunel mechanism. But in many other experiments -e.g. for ion acceleration from dense plasmas-preserving very steep interfaces is not strictly necessary. Such experiments are therefore often performed without contrast improvement after temporal compression: the laser pulses then typically have a very high temporal contrast up a few ps before the main pulse, but their short-time contrast is often not as good, due e.g. to remaining high-order terms in the pulse spectral phase. With typical plasma expansion velocities in the 50 nm/ps range, this so-called coherent pedestal on the ps time scale will lead to density gradient scale lengths of the order of λ at the arrival of the main pulse. In these conditions, stochastic heating should be the dominant coupling mechanism, and can now be readily be identified through the multiple signatures described in this work.