Plasma density limits for hole boring by intense laser pulses

High-power lasers in the relativistic intensity regime with multi-picosecond pulse durations are available in many laboratories around the world. Laser pulses at these intensities reach giga-bar level radiation pressures, which can push the plasma critical surface where laser light is reflected. This process is referred to as the laser hole boring (HB), which is critical for plasma heating, hence essential for laser-based applications. Here we derive the limit density for HB, which is the maximum plasma density the laser can reach, as a function of laser intensity. The time scale for when the laser pulse reaches the limit density is also derived. These theories are confirmed by a series of particle-in-cell simulations. After reaching the limit density, the plasma starts to blowout back toward the laser, and is accompanied by copious superthermal electrons; therefore, the electron energy can be determined by varying the laser pulse length.

In the interaction of intense laser fields with overdense targets, the laser radiation pressure pushes electrons at the plasma surface, which sets up an electrostatic field originating from the charge separation that ultimately accelerates ions in the forward direction. The laser light proceeds to push the plasma surface, which has the relativistic critical density γn c into the target. Here, γ is the relativistic factor of electrons, n c = m e ω 2 L /(4πe 2 ) is the non-relativistic critical density, ω L the laser frequency, m e the rest mass of electron, and e the fundamental charge. This process is referred to as the laser hole boring (HB) [26][27][28][29][30][31] . The HB or the surface steepening has been considered to be important for applications, e.g., laser channeling [32][33][34] , high harmonic generation 35,36 , and plasma mirror 37,38 . For these applications, how long the steepened clean interface is sustained is an essential question.
Laser absorption and hot electron generation also depend on the steepening of the plasma surface in the HB process. In the conventional model of the HB 26,[28][29][30] , which is based on the momentum transfer equation from laser to ions, laser lights bore holes as long as the laser pulse continues, which implies there is no density limit for the HB.
Recently, for the laser-solid interactions in the multi-ps regime, it is reported that the plasma on the laser-irradiated front surface expands significantly in the order of 10 μm during the laser irradiation, and superthermal electrons beyond the conventional ponderomotive scaling 26 are generated [39][40][41] . A recent LFEX laser experiment demonstrated that the hot electron temperature increases drastically when the laser pulse duration is extended from 1 to 4 ps while keeping the peak intensity same 42 . Ion acceleration in such a multi-ps laser interaction cannot be described by the conventional target normal sheath acceleration model that assumes plasma expansion with isothermal electron temperature [43][44][45] . The nonisothermal plasma expansion theory is proposed to explain the multi-ps laser-driven ion acceleration experiments where the time evolution of hot electron temperature is important 46,47 . In order to control such a hot and/or superthermal electron generation in the multi-ps regime, the cause of the plasma blowout to the front side is essential to be figured out.
In this study, we find that the HB is stopped in the ps time regime even while the laser pulse is still on, and the surface plasma eventually starts to blowout to the front side. We here develop a theory that explains the transition from the HB to the plasma blowout regime. Based on a pressure balance relation between laser radiation pressure and electron thermal pressure, we derive the limit density for the HB as 8Ra 2 0 n c $ 5:8RI 18 λ 2 μm n c , above which the laser field cannot push beyond. Here, a 0 = eE 0 /(m e cω L ) is the normalized laser field amplitude, E 0 the amplitude of the laser electric field, R the reflectivity, I 18 the laser intensity normalized by 10 18 W cm −2 , λ μm the laser wavelength normalized by 1 μm. We show the validity of the derived model by using particle-in-cell (PIC) simulations. In addition, we obtain the time scale for the transition from the HB to the blowout regime based on the momentum transfer equation in a preformed plasma 48 . The transition time scale is found in the ps regime, and therefore, the prediction by this theory will appear in multi-ps laser experiments.

Results
Process of hole boring. As an introduction of the HB, we show an interaction of plasma with linearly polarized high-intensity laser field by using Fig. 1. Here, the laser intensity is in the relativistic regime, I =5×10 18 W cm −2 so thatâ 0 ¼ 2 where the hat indicates the peak value. We consider a thick plasma with initial ion and electron densities n i0 and n e0 = Zn i0 , respectively, where Z is the ion charge state. The target electron density n e0 is an order of magnitude higher than the relativistic critical density γn c . We assume a pre-plasma of scale length L at the front side where the electron density n e increase linearly from zero to n e0 . Figure 1a-f are the results we obtained in a PIC simulation in one-dimensional (1D) geometry. The laser field transmits through the underdense region n e < n c . After reaching the critical density n c , the laser field can further penetrate into the plasma n e ≤ γn c due to relativistic transparency. Since γ fluctuates in the interaction, the plasma is not transparent completely in this regime, and thus, the laser field starts to push electrons at the pulse front. In this phase, the pile up of electrons swept up by the laser is fast enough to create a strong electric field that can accelerate ions at a speed exceeding the sound velocity. Consequently, a collisionless shock is formed at the front as seen in the ion phase plot for the longitudinal direction in Fig. 1a. Above the relativistic critical density, n e > γn c , the plasma is opaque and the speed of the interaction front decreases significantly. The laser is incapable of driving shocks by pushing the overdense electrons. We refer to this stage to as the 'HB' phase. In the HB phase, electrons in the front surface are also pushed by the laser light making a charge separation. If the laser pulse can sustain the charge separation, ions start to move forward. Note that for the HB, the pulse has to be longer than the ion response time scale, 2πω À1 pi ; which is in the order of 100 fs at the critical density. Here, ω pi is the ion plasma frequency. As seen in Fig. 1b, the front of the shock generated in the relativistic transparency phase proceeds forward in a faster speed than the laser interaction front, that is the HB front.
As the HB proceeds, we find that the average ion momentum at the HB front p xf changes from positive as in (b) to zero as in (c). This indicates that the HB front cannot go further when it reaches the condition that the average longitudinal density flux of ions is zero. The asymmetric relation between positive and negative ion density flow can be seen in (d) and (e) where the distributions of ion longitudinal momentum in the region of the pulse front, sampled within Δx = 0.2 μm, are shown for the corresponding times for (a) and (b), respectively. One can see that the distribution that has a notable amount of positive p xf component in (e) turns to a symmetrical distribution with respect to the p xf axis in (f). Therefore, at t = 0.6 ps ((c) and (f)), a stationary state of the plasma front is established, which corresponds to the stopping condition for the HB. At this time, the laser radiation pressure balances with the plasma pressure. The establishment of the stationary state is owing to the electron heating during the HB, by which a hot dense electron cloud is formed behind the HB surface, and thus the negative ion density flow is driven. Since the plasma is heated by the laser continuously, the plasma will start to expand when the electron pressure exceeds the laser radiation pressure. This corresponds to the transition from the HB to the plasma blowout. Note that the stationary state never appears in the conventional aspect of the HB where the electron heating is not taken into account.
In Fig. 2a, we show the time evolution of the positions of the critical density n c and relativistic critical density γn c for the same simulation of Fig. 1. The relation γ 1 þ ð1 þ RÞϵ 2â2 0 =2 À Á 1=2 is used where ϵ = 1 for linear polarization and ϵ ¼ ffiffi ffi 2 p for circular polarization. We here use the reflectivity R at the transition time from the HB to the blowout. Figure 2b shows the reflectivity R, which keeps an almost constant value around the transition time from the HB phase to the blowout phase. In early time t < 0.4 ps, positions of n c and γn c are pushed toward positive x direction, however, during the time t = 0.5-0.7 ps that is shaded in Fig. 2a, b, the plasma front hardly moves from x~9.45 μm. After this time, the plasma starts to blowout to the negative x direction.
In Fig. 2d, e, we show the electron distributions in the phase space for the longitudinal x direction at t = 0.4 ps, that is the HB phase, and t = 1.0 ps, that is the blowout phase, respectively. In the HB phase (d), electrons are experiencing J × B acceleration at the peripheral of the HB surface, x~9.4 μm. In the low-density region in front of the HB surface, mainly the electrons in the region where n e < 0.1n c interact directly with the incident laser field and gain momenta p xe /m e c~5 (2.1 MeV). On the other hand, in the blowout phase (e), the region where 0.1n c < n e < n c expands about 1 μm, which causes the enhanced acceleration of electrons 40 . The spectra of electron energy ε e at t = 0.4 and 1.0 ps are shown in Fig. 2c. The number of hot electrons in the MeV range are increased by an order of magnitude, and the high-energy slope is enhanced, resulting in temperatures higher than the ponderomotive temperature T p and also than the Beg scaling T B 49 .
Hole boring limit density. We now consider the pressure balance at the laser-plasma (LP) interaction front in the stationary state, where the HB cannot proceed further. We assume that a plane laser field with the normalized amplitude a 0 is irradiating an overdense plasma distributed in the region x ≥ 0. At the stationary state, the electrons are pushed slightly by the laser radiation pressure with the distance ' s =2 where ' s ¼ m e c 2 = 4πn e e 2 ð Þ ð Þ 1=2 is the skin depth for the incident laser field. Note here that since the laser amplitude decays in ' s , the intensity decreases in the scale length of ' s =2. Then, only ions remain in the front surface 0 x<' s =2 with the density n i , while the inside x ! ' s =2 is filled by both ions with density n i and electrons with density n e = Zn i . The positive electrostatic field generated by the charge separation at the front surface, E s ¼ 2πen e ' s , drags ions to the positive x direction. However, in the stationary state of the LP interaction front, the ion density is maintained by the negative density flow of ions as seen in Fig. 1c, f. The pressure balance relation at the interaction surface between laser field and electrons is derived by integrating the electron fluid equation of motion in the stationary state as described in the Methods as The left-hand side represents the laser radiation pressure with intensity I, the first term on the right-hand side (RHS) corresponds to the electron pressure with temperature T e , and the second term on the RHS denotes the sheath electrostatic potential energy density, which corresponds to the surface tension in the width of the skin depth. Note that some of the previous papers have also included the electron pressure term 28,30 . However, the stationary state sustained by the surface tension Þis irradiated from the left boundary with a 0.15 ps Gaussian pulse leading edge. A fully ionized deuteron plasma of an overdense density n e0 ¼ 40n c and a 1 μm pre-plasma with linear density gradient is initially distributed. a-c are the phase plots of ions in the longitudinal direction at times t = 0.18, 0.40, and 0.60 ps, respectively, where the color represents the number of PIC particle. In d-f, the longitudinal momentum distributions at the laser-plasma (LP) interaction front at the corresponding times are presented. Here, the momentum p xf is normalized by M i c where M i is the ion rest mass of the field has never been considered, which is the critical difference of Eq. (1) from the conventional descriptions of the HB. We assume that the plasma is composed of hot and bulk electron components as n e = n h + n b . We here assume n e T e~nh T h where T h is the hot electron temperature, by neglecting the bulk electron pressure.
The hot electron density and velocity can be determined by the conservation of energy density flux, Here, β e = v e /c is the ratio between electron drift velocity v e and the speed of light c, and α ≡ ir/2 is the geometrical factor where r = 1 for the non-relativistic Maxwell momentum disrtibution, r = 2 for the relativistic Maxwell (Maxwell-Jüttner) momentum distribution, and i = 1, 2, or 3 represents the dimension of momentum distribution. Hereafter, we consider r = 2. For quasi-1D relativistic interactions, we can assume i = 1, and thus, α = 1. Equation (2) indicates that the absorbed laser energy flux is carried by electrons. We simplify Eq. (2) by replacing n e T e β e on the RHS by n h T h . This approximation is valid for the relativistic condition, a 0 > 1 and thus β h $ 1. When the laser field amplitude becomes super relativistic (a 0 ≳ 300), other energy loss mechanisms such as the radiation damping will come into play, and hence, the pressure balance at the LP interaction surface will change 50 . From Eqs. (1) and (2), we can derive the front density n e for the stationary state. We define this stationary density as the limit density of the HB, n s , which is obtained as Here, we used the relation I=c ¼ n c m e c 2 ϵ 2 a 2 0 =2. Equation (3) represents the limit density for the HB, above which the laser light cannot proceed forward. For this limit density, the laser light with intensity I is incapable of sustaining a charge separation that is sufficient for driving positive mean density flux of ions. When one assumes a 1D relativistic Maxwell distribution α = 1, the relativistic limit for the hot electron velocity β h = 1, and linear polarization ϵ = 1, Eq. (3) reduces to Equation (4) presents the dependence of the HB limit density on normalized laser amplitude a 0 and reflectivity R. Note that for the ideal condition for the HB, we generally consider the density regime n e~γ n c~a0 n c . On the other hand, the HB limit density n s given by Eq. (4) corresponds to a higher density than γn c for a 0 > 1 when the reflectivity is not so low satisfying R>a À2 0 =8. In other words, the HB proceeds beyond the density regime n eã 0 n c , and stops at the density n e ¼ 8Ra 2 0 n c . The maximum HB limit density is obtained for R = 1, that is n smax =n c ¼ 8a 2 0 . In the case of multi-dimensional interactions (α > 1), due to the additional freedom of lateral energy diffusion, the plasma pressure becomes lower even with the same R as seen from Eq.
(2). Consequently, laser lights can proceed higher density in multi-dimensional geometry than in the 1D case. However, the maximum limit density 8a 2 0 n c is not changed by the dimension effect as can be confirmed by Eq. (3).
In Fig. 3, we plot Eq. (4) for various reflectivities R by dotted lines. The bold solid line for R = 1 corresponds to the maximum HB density. Hence, the shaded area above the solid line is the prohibited area where laser light cannot reach that density. We also present the maximum density to which the laser light can reach in 1D and two-dimensional (2D) (quasi-1D) PIC simulations. We see that all the simulation points are below the maximum HB limit density. The maximum density obtained in each simulation agrees with the theoretically obtained limit density n s for R with a difference limited in the range of ±0.1, except for the result of 2D circular polarization (C-pol), which will be addressed below. Note that in the 2D geometry, the reflectivity tends to be lower than the 1D geometry, so that the maximum density in the 2D simulation is typically far below the theoretical limit, that is n smax for R = 1. Although the reflectivity R changes in time, the variation of R in the shaded interval in Fig. 2b, that is ±0.1 ps, around the HB stationary state is small about 0.01 as shown in Fig. 3 for each simulation. Hereafter, we estimate the time scale t s to reach the HB limit density. As seen in Fig. 2, superthermal electrons appear after the transition to the blowout phase. Hence, t s is important to control the electron heating, which is essential for various applications of intense LP interactions.
During the HB stage, the momentum transfer from laser light to ions, in the frame moving with the ion front velocity v f , is expressed by where M i is the ion mass, and the electron pressure term is neglected for simplicity. Here we assume an exponential density profile with the scale length L as n e = Zn i = γn c exp(x/L). The front position of the HB (x f ) is obtained by integrating v f in time as as in ref. 48 . Then, the electron density at the HB front n f is given as The transition time scale t s is obtained by solving Eqs. (6) and (7) for t with substituting n s to n f in Eq. (7) as where A ≡ (Zm e /(2M i )) 1/2 /2, and n s /n c on the RHS is given by Eq.
(3). Here, we define the time when the laser front reaches to the position of n e = γn c as t = 0. In the above derivation, we assume the laser field amplitude a 0 is constant in time. Note that this assumption is valid for the present simulations where the pulse profile is constant in most of the interaction time. In the case of laser pulses with Gaussian temporal profiles, the time scale for reaching the HB limit density n s is estimated approximately to be 2t s . Equation (8) implies that the LP interaction front proceeds by boring the exponentially distributed plasma with the initial scale length L, and then subsequently reaches the HB limit density n s in the time scale of t s . In Fig. 4a, the transition times t s given by Eq. (8) are plotted as a function of ϵa 0 by black lines for various reflectivities R. Here, the relativistic 1D momentum condition, that is α = 1 and β h = 1, and pre-plasma scale length L = 2 μm are used. The derived time scales are in the order of ps. The times when the HB stops in the simulations agree well with the theoretical lines. In the case of circular polarization, which is the combination of P-and Spolarizations, the P-polarization has higher absorption than the S-polarization, resulting that the polarization of the laser field during the interaction is not circular anymore. Namely, the interaction becomes a two-beam-like interaction, so that the result departs from the single beam scaling, e.g., n s and t s .
As an example of the HB process, Fig. 4b shows the electron density profile in position-time space obtained by the 2D simulation shown in Fig. 4a at ϵa 0 = 4. The laser front gets to the position of n e = γn c at t = 0 by the process of relativistic transparency. Up to time t = t s , the HB surface moves to the higher density direction. After passing t = t s , the HB surface is pushed backward (lower density side), which corresponds to the plasma blowout.
To check the multi-dimensional effect, we demonstrated a 2D simulation for a tightly focused laser pulse. The results are plotted in Figs. 3 and 4 by white circles. Both results are in a good agreement with values derived by Eqs. (3) and (8), assuming the 2D condition α = 2 and reflectivity R obtained in each simulation. The obtained scalings are therefore confirmed to be applicable in multi-dimensional situations.
In conclusion, we derived the limit density for the HB n s based on a balance relation between laser radiation pressure and plasma pressure. The HB limit density is found to be n s ¼ 8Ra 2 0 n c $ 5:8RI 18 λ 2 μm n c . After reaching the HB limit density n s , the laser light is incapable of sustaining the charge separation that is sufficient for driving the forward mean density flux of ions in the HB surface, and thus, the laser light can no longer proceed into the higher density region. By using the PIC simulation, we demonstrated that the HB reaches the stationary state when the laser pulse front proceeds to the density equal to  4), that is for the 1D relativistic condition α = 1 and β h = 1, are presented by black solid and dotted lines for various reflectivities R. The shaded area corresponds to that prohibited by the theory for any reflectivity R. Blue triangles and blue square show the PIC simulation results for one-dimensional geometry with linear and circular polarizations, respectively. Red circles and red square are the 2D PIC simulation results for the quasi-1D geometry with widely focused (spot size of 60 μm) linear-P and circular polarizations, respectively. The white circle corresponds to the PIC simulation result for a tightly focused (spot size of 1.5 μm) P-polarized laser pulse in two-dimensional geometry, for which the corresponding theoretical value of n s assuming the 2D relativistic condition α = 2 and β h = 1 in Eq. (3) is presented. Reflectivity R obtained in each simulation is also shown the limit density n s . We derived the transition time t s , at which the plasma turns from the HB to the blowout, as Eq. (8).
The transition time t s is critical for applications such as laser channeling, high harmonic generation, and plasma mirror, which require steepened clean interface for efficient operation. The blowout plasma interacts with the laser field resulting copious superthermal electrons, which are important to increase the efficiency of the applications, e.g., ion acceleration 46,47 and pairplasma creation 51,52 with multi-ps kiloJoule laser lights.

Methods
Numerical simulations. Simulation results in 1D geometry shown in Figs. 1,2-4 are obtained using a fully relativistic collisional PIC code EPIC3D 53 . The calculations are executed in 2D spatial dimension where the size of simulation box is L x = 40.96 μm in the laser propagation direction, while we use four meshes in the transverse y direction with the mesh size of 10 nm. A laser field with the normalized amplitude of a 0 ðtÞ ¼â 0 f a ðtÞ is excited by an antenna at the left boundary. Here,â 0 is the peak value and f a is the pulse shape factor whose time dependence at the antenna position is given by f a ¼ exp t À t 0 ð Þ 2 =τ 2 L Â Ã for t < t 0 and f a = 1 for t ≥ t 0 , where t 0 = 0.1 ps and τ L = 0.15 ps. The laser wavelength is λ L = 1.05 μm. In the calculations for Figs. 1 and 2, the peak normalized amplitude isâ 0 ¼ 2; which corresponds to the intensity of 5 × 10 18 W cm −2 . A uniform fully ionized neutral deuteron plasma is distributed initially from x = 10-20 μm with a linear pre-plasma of length L = 1 μm whose electron density increases from zero at x = 9 μm to the uniform plasma density n e0 >8â 2 0 at x = 10 μm. In the calculations for Figs. 1 and 2, n e0 = 40n c is assumed. The initial plasma distribution and the incident laser field are uniform in the transverse y direction.
Simulation results in 2D geometry shown in Figs. 3 and 4 are obtained using a fully relativistic PIC code PICLS 54 . The size of simulation box is L x = 40 μm in the laser propagation direction and L y = 120 μm in the transverse direction with the mesh size of 20 nm. The laser pulse function and the wavelength are same as those used in the EPIC simulations. A uniform fully ionized neutral deuteron plasma is distributed initially from x = 20-40 μm with an exponential pre-plasma of scale length L = 2 μm whose electron density increases from 0.6n c at x = 10 μm to the uniform plasma density 100n c at x = 20 μm. The initial plasma distribution is uniform in the transverse y direction. The laser field is focused at x = 20 μm to the spot sizes of 60 μm in the wide focus case, and 1.5 μm in the tight focus case, which is close to the diffraction limit to reduce the self-focusing effect. Since the spot size for the wide focus cases is much larger than the plasma scale length L, the interaction around the beam center can be regarded as quasi-1D, that is α~1, where the effect of multi-dimensional energy diffusion is small.
The reflectivities R for the simulations referred in Figs. 3 and 4 are observed around the transition time t s ± Δt, where Δt = 0.1 ps, which corresponds to the scale of the ion response time 2πω À1 pi .
Derivation of the pressure balance equation. The pressure balance relation for the stationary state of the HB Eq. (1) is derived from the electron fluid equation of motion in the stationary state given by where ϕ is the electrostatic potential, and all quantities are averaged over the laser period. Here, we consider a 1D geometry where an overdense plasma is distributed in the region x ≥ 0. Electrons are pushed by the laser radiation pressure with the distance ' s =2; where ' s is the skin depth, while ions remain in the front surface 0 x<' s =2 with the density n i ( = n e /Z). We assume the scales of spatial variation ∇n e $ 2n e =' s and ∇a 0 $ Àa 0 =' s in the interface, and obtain ∇ n e a 2 0 =γ À Á ¼ Àn e ∇a 2 0 À Á =ð2γÞ. Here, γ 1 þ ð1 þ RÞϵ 2â2 0 =2 À Á 1=2 is used, and we approximated ð1 þ RÞϵ 2â2 0 =2 $ γ 2 in the calculation of ∇γ. Using the above assumptions for spatial variation, integrate Eq. (9) from x = −∞ to x = +∞, then which is the pressure balance relation Eq. (1). Here, the relation I=c ¼ n c m e c 2 ϵ 2 a 2 0 =2 is used, and the electron density at the laser front is given by γn c . E 2 s =8π is the sheath electric field energy density generated by the charge separation in the interface 0 x<' s =2. In obtaining this sheath field term, we used the Poisson equation and neglected the ion pressure.
Note that in the HB stage where the ion front has not reached the stationary state, the charge separation field eventually moves the ion front in the forward (positive x) direction. Therefore, for describing the HB stage, the last term on the RHS of Eq. (10) is replaced by the kinetic energy density of ions moving with the HB velocity v f , and then the momentum transfer equation in the frame moving with the velocity v f is obtained as Eq. (5).
Data availability. The data that support the findings of this study are available from the corresponding author upon reasonable request. for linear and circular polarizations, respectively. Here, the relativistic 1D condition, that is α = 1 and β h = 1, and pre-plasma scale length L = 2 μm are used. Theoretical lines given by Eq. (8) are presented by black lines. Red circles and red square are the 2D PIC simulation results for the quasi-1D geometry with widely focused (spot size of 60 μm) linear-P and circular polarizations, respectively. The blue triangle (square) is the PIC simulation result for one-dimensional geometry with linear (circular) polarization. The white circle corresponds to the PIC simulation result for a tightly focused (spot size of 1.5 μm) P-polarized laser pulse in twodimensional geometry, for which the corresponding theoretical value of t s assuming the 2D relativistic condition α = 2 and β h = 1 in Eq. (8) is presented. Reflectivity R obtained in each simulation is also shown. b Electron density profile in position-time space obtained by the 2D simulation shown in a at ϵa 0 = 4 with the wide laser focal spot. The presented density is the averaged value over the central region y = ±1 μm with respect to the laser axis. Black solid and dashed curves are contour lines for critical density n c and relativistic critical density γn c , respectively