A theoretical model of laser-driven ion acceleration from near-critical double-layer targets

Laser-driven ion sources are interesting for many potential applications, from nuclear medicine to material science. A promising strategy to enhance both ion energy and number is given by Double-Layer Targets (DLTs), i.e. micrometric foils coated by a near-critical density layer. Optimization of DLT parameters for a given laser setup requires a deep and thorough understanding of the physics at play. In this work, we investigate the acceleration process with DLTs by combining analytical modeling of pulse propagation and hot electron generation together with Particle-In-Cell (PIC) simulations in two and three dimensions. Model results and predictions are confirmed by PIC simulations—which also provide numerical values to the free model parameters—and compared to experimental findings from the literature. Finally, we analytically find the optimal values for near-critical layer thickness and density as a function of laser parameters; this result should provide useful insights for the design of experiments involving DLTs. Laser-driven ion acceleration is an important topic for next-generation compact accelerators and material characterisation. The authors present a theoretical study on ion acceleration with near-critical double-layer targets that supports the experimental realisation of these targets and the interpretation of experiments of laser-ion acceleration.

L aser-driven ion acceleration is a well-established research topic [1][2][3] . Many distinct acceleration mechanisms have been explored so far (such as radiation pressure acceleration 4,5 , breakout-after-burner 6 , relativistic induced transparency 7,8 , collision-less shock acceleration 9 , magnetic vortex acceleration 10 ). Yet, Target Normal Sheath Acceleration (TNSA) 11 is arguably the most established and robust ion acceleration scheme. The unique properties of TNSA ions (e.g. broad exponential spectrum with tens of MeV cut-off energies, high bunch density, ultrafast duration) could be exploited in the near future for a number of interesting applications for materials characterization (e.g. particle induced x-ray emission 12,13 ), in nuclear science (e.g. bright neutron sources 14,15 , radioisotopes production 16 ) and in the study of harsh radiation environment effects in materials science 17,18 . Nonetheless, TNSA is affected by a low conversion efficiency of laser energy into energetic ions, which limits the use of compact laser systems for such applications. A viable route to overcome this limit may be to use double-layer targets (DLTs) consisting of a thin solid foil coated with a near-critical density layer, where the critical density n c ¼ m e ω 2 =4πe 2 marks the transparency threshold for the propagation of an electromagnetic wave with frequency ω (m e is the electron mass and e is the elementary charge). Both numerical simulations [19][20][21][22] and experiments [23][24][25][26][27][28][29][30][31] have demonstrated that a laser pulse strongly interacts with the near-critical layer, generating a larger number of energetic electrons and increasing the ions energy and number with respect to an uncoated target. This higher acceleration efficiency has been attributed to several phenomena: the laser selffocusing (SF) induced by the radial dependence of the refractive index within the channel generated via the ponderomotive force 32,33 , the generation of a strongly magnetized channel carrying high currents 34 and the Direct Laser Acceleration of electrons (DLA) through the betatron resonance [35][36][37] . Since in real experiments near-critical layers are usually nanostructured materials, the effect of realistic nanostructures was studied, highlighting differences from the ideal uniform case 38 . In addition, it has been demonstrated that the laser interaction with near-critical plasmas follows an ultra-relativistic scaling with respect to the transparency factor n ¼ n e =γ 0 n c (or the normalized plasma density), where γ 0 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ a 2 0 =} p is the mean Lorentz factor of the electron motion due to a laser with normalized amplitude a 0 and linear (} ¼ 2) or circular (} ¼ 1) polarization 39,40 .
DLTs have been extensively studied. However, the literature does not provide a comprehensive theoretical view able to account for the various effects at play and to provide an estimation for the optimal DLT parameters. In this work, we present a theoretical description of laser-DLT interaction and the consequential ion acceleration process. We develop a model which is comprehensive of all the main relevant phenomena (laser selffocusing, electron heating and ion acceleration) including the dependence of the quantities of interest on a large set of parameters altogether (e.g. near-critical layer density and thickness, laser waist and duration). Modelling the essential physical aspects of the interaction, we are able to unfold a relativistic scaling for the accelerated ions. We support our arguments with an extensive multi-dimensional (2D/3D) particle-in-cell (PIC) simulations 41 campaign. By exploiting suitable approximations, we identify a set of optimal DLT parameters (i.e. density and thickness of the near-critical layer) which maximize the ion energy enhancement with respect to the uncoated target case. To this purpose we develop simple estimations that can be easily carried out without having to perform many time-consuming PIC simulations. Our results provide a convenient guide both for the design of engineered DLTs and for the interpretation of laserdriven ion acceleration.

Results
Laser-DLT interaction. As mentioned in the introduction, the interaction between a super-intense laser pulse and a near-critical plasma leads to the onset of complex phenomena, characterized by a strong coupling between the electromagnetic field and the plasma. When the transparency factor is lower than one, n ¼ n e =γ 0 n c <1, the laser can propagate within the plasma and dig a channel in it, which in turn induces relativistic-ponderomotive SF (see Fig. 1a). In this process, in the whole interaction volume, the laser generates hot electrons (see Fig. 1b), which are characterized by a large particles density and super-ponderomotive mean energy. Furthermore, the strong currents induced in the channel generate a dipole magnetic vortex with fields easily exceeding 10kT (see Fig. 1c) which constrain the electrons to a directional motion. Finally, part of the pulse reaches the substrate, heating additional electrons from the surface and starting the TNSA process. In this complex framework, we account only for the essential mechanisms for a sufficiently accurate description of the interaction and we restrict our analysis to a simple scheme, depicted in Fig. 1d: normal incidence, linear polarization and homogeneous plasma. We also assume that the laser pulse duration is short enough that the plasma ions can be considered almost motionless during the interaction.
In this section, we focus first on the evolution of the pulse waist due to SF in the near-critical plasma for which we propose a simple law. Then, we describe the pulse energy loss and the amplitude amplification with equations valid in D-dimensional geometry (D ¼ 1; 2; 3) which depend on few free parameters. This formulation allows us to validate the model through both 2D and 3D PIC simulations (in particular, in this work, with a large amount of 2D ones and a limited number of more expensive 3D simulations). Moreover, we calculate the hot electrons mean energy and total particles number for both the near-critical layer and the substrate of the DLT.
Firstly, it should be mentioned that SF is one of the most important features of the laser-DLT interaction. The phenomenon occurs when ultra-high intensity lasers propagate within near-critical plasmas. In the literature 21,32 a simple expression for the minimum achievable laser waist w m has been proposed: where λ is the laser wavelength. In order to carry out quantitative calculations we propose the following relation for the laser waist evolution w(x) against the pulse propagation length x inside the plasma, which is based on a thin-lens approximation: ; ð2Þ where w 0 is the initial laser pulse waist, x R the Rayleigh length, and l f the SF focal length, namely the position at which the laser is most strongly focused in the plasma, w(l f ) = w m . Equation (2) describes the plasma focusing phenomenon as a thin lens. While a plasma should act more like a gradient-index lens 42 , a simple thin lens model is able to predict the waist measured in 2D/3D PIC simulations (see Fig. 2a, b). It should be pointed out that the equation should be reasonably valid for n≳0:05 (under this threshold we don't expect a focusing-defocusing behaviour, but a stable propagation 43 ) and within a pulse propagation length of the order of the SF focal length (x≲2l f ). At longer distances other phenomena are induced, as laser-plasma non-linear instabilities 42 and filamentation 44 , which distort the ideal Gaussian shape of the   w (a, b), the pulse energy ε p (c, d) and the pulse amplitude a (e, f) evolution inside a uniform plasma at different propagation lengths x, normalized to the transparency factor n and the laser wavelength λ, in 2D (a, c, e) and in 3D geometry (b, d, f). The full lines and points refer to particle-in-cell (PIC) simulations with intensity ranging in a 0 = 2-8 and plasma density n e /n c = 0.5-2, while the dashed lines refer to our model. The pulse energy in (c, d) is calculated from PIC simulations by integrating the total electromagnetic energy, including the reflected part of the pulse; these values are compared to the results of the model by adding the reflectance R D to the calculated energy. COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-020-00400-7 ARTICLE COMMUNICATIONS PHYSICS | (2020) 3:133 | https://doi.org/10.1038/s42005-020-00400-7 | www.nature.com/commsphys laser. Although we expect the model to fail below a sufficiently small waist, we have extensively validated it for w = 4 λ, which is a very common value for tight focused ultra-intense lasers. Indeed, focusing an ultra-intense beam closer to the diffraction limit poses considerable technological challenges.
It is worth mentioning that l f % w 0 = ffiffiffi n p when w 0 ) w m (which is a reasonable approximation). Under this condition, Eq. (2) can be approximated by Here we defined the relativistically normalized space variable x ¼ ffiffiffi n p x=λ (defined in ref. 39 within the framework of the ultrarelativistic similarity theory). This equation emphasizes that the evolution of the laser waist inside the near-critical plasma depends on x only through the x variable and leads to self-similar curves for equal initial waist w 0 /λ and for constant transparency factor n. We expect the pulse propagating inside the channel to lose part of its energy heating electrons and to increase its intensity due to SF. To calculate the laser energy loss in the propagation we can assume that all the electrons inside the channel are heated with the well-known ponderomotive scaling, with arguments similar to those presented in ref. 45 : where R ch is the plasma channel radius, V DÀ1 R ch ðxÞ DÀ1 is the plasma channel section, Þ-dimensional ipersphere with unitary radius (V 0 = 1, q is the electrons mean Lorentz factor in linear polarization at a given pulse propagation length, while C nc a constant allowing to estimate the mean electron energy as C nc γ x ð Þ À 1 ð Þ m e c 2 . Thus C nc absorbs the details of the electron heating process 46 , which may hold super-ponderomotive features as seen in other works 29,37,45 . Considering an ideal Gaussian pulse, both in time and space, we can express its initial energy in D dimensions as ε p0 ¼ π D=2 2 ÀD=2À1 m e c 2 n c a 2 0 w DÀ1 0 τc 47 , where τ is the fields temporal duration (1/e) and c the speed of light in vacuum (the full-width-half-maximum, FWHM, temporal duration and the focal spot over the intensity, frequently used in the experimental field, can be retrieved as τ I FWHM ¼ ffiffiffiffiffiffiffiffiffiffi ffi 2log2 p τ and w I FWHM ¼ ffiffiffiffiffiffiffiffiffiffi ffi 2log2 p w 0 ). Then, the normalized energy loss assumes the following form: where we introduced the ratio of the plasma channel radius to the waist r c ¼ R ch ðxÞ=w x ð Þ, assuming it constant. We will adopt r c as a free parameter to describe the channel radius as a function of the pulse waist. It is reasonable to think that the channel will be larger than the waist (hence r c > 1), but of the same order of magnitude.
To solve Eq. (7) we need γ(x) along the propagation length, which can be calculated from the pulse amplitude a(x), through the following equation: Here we have neglected the pulse temporal shaping effects because the intensity amplification during the propagation is mostly due to SF, rather than to the temporal compression. Indeed, the temporal compression takes place only in special conditions, and τ is at most halved 21 , while the SF waist can reach the diffraction limit, implying a waist reduction which can easly exceed a factor of 10. Equation (7) depends on the laser waist equation (Eq. (2)) and can be numerically solved, coupled with Eq. (8), with a finite difference method. In order to do so, the initial condition ε p ð0Þ=ε p0 must be imposed. While one could approximate this initial condition to 1, we also take into account that part of the laser pulse can be reflected by the plasma. We propose 2D/3D relations for the reflectance R D (Eqs. (29) and (30) in "Methods"), which are used to set the initial condition ε p ð0Þ=ε p0 ¼ 1 À R D . Since Eqs. (7) and (8) depend on the waist evolution and thus on the self-similar variable x, both the pulse energy and the laser amplitude evolutions along the path length can be scaled to self-similar curves, even for very different initial conditions (i.e. plasma density, initial waist, laser intensity, pulse temporal duration). Equation (7) can be solved for different values of the free parameters C nc and r c to find the pulse energy loss and the laser amplitude during the propagation. The free parameters will assume different constant values for different problem dimensionality and are fixed by fitting the numerical results with the theoretical model. Figure 2a-f shows the model results (dashed lines) compared with PIC simulations results (solid lines with points) in normalized units. A good agreement is observed in all cases for both dimensionalities. We note that the fitted values for C nc and r c (see Table 1) are consistent with their physical interpretation. Indeed we obtain C nc ≳1 and r c ≈ 2, which are compatible with that the hot electron heating may be slightly super-ponderomotive in near-critical plasmas and the channel radius is expected to be larger but not too large than the waist. Moreover, we observe in 2D (Fig. 1b) that the channel radius at x = 0 is about 6λ, corresponding to r c,2D~1 .5, to be compared with the fitted value of 2.0. The C nc value is found to be a little higher than one both in 2D and 3D, as expected. The fact that C nc,3D is lower than C nc,2D by a factor of about 1.5 could be explained by considering that, at equal peak amplitude, the mean laser Ratio of the channel radius (formed by ponderomotive force) to the pulse waist n 1.2 × 10 −3 n c 5 × 10 −2 n c Hot electron density normalization constant which scales the estimated proton energy in the quasistationary model The free parameters of our model for the two-dimensional and three-dimensional cases are reported. Cnc and Cs are proportionality constants which correct the ponderomotive scaling for the near-critical layer and substrate electrons, respectively (see Eqs. (6) and (13)). rc is a proportionality constant between the near-critical channel radius and the laser waist (see Eq. (7)).ñ is a normalization constant of the hot electron distribution function (see Eq. (26)).
amplitude is lower in 3D than in 2D. This dimensionality effect has been reported also in ref. 46 . Now that we have a model for the pulse propagation in a nearcritical plasma, we characterize the hot electron population that is generated in the same process. Solving Eq. (7) allows retrieving the fundamental properties of the hot electrons heated by the laser. Firstly, assuming that the electrons are the principal absorbers of the laser energy, we can define the absorption efficiency η nc , i.e. the fraction of the laser energy that is converted into hot electrons kinetic energy: Within this approximation we neglect that the pulse energy can be directly absorbed by plasma ions or emitted as secondary radiation, such as the synchrotron-like emission. It is worth noting that these hypotheses are reasonably valid when the pulse duration is short enough (tens of fs) and the intensity sufficienlty low (a 0 < 50 38 ). We confirm that these approximations hold by comparing the calculated η nc (x) with the 2D PIC simulations data (see Fig. 3a).
In this framework, it is also possible to calculate the total number of hot electrons N nc at a given x position, by integrating the electron density inside the plasma channel: Exploting Eqs. (9) and (10), the mean hot electron energy E nc is retrieved: The good agreement between the calculated E nc with the one obtained from the PIC simulations can be observed in Fig. 3b.
Up to now, we have described laser interaction with a semiinfinite near-critical plasma. Nonetheless, if we want to describe laser-DLT interaction, we have to take into account the effect of a thin solid substrate, with a given thickness d s , coupled with a near-critical plasma with length d nc . To do so, we have to consider that the laser pulse can reach the substrate with some residual energy and produce hot electrons at the substrate interface, with an absorption efficiency η s defined as where N s is the total number of hot electrons generated at the surface of the substrate, and E s is their mean energy which can be expressed by the ponderomotive scaling, as done in Eq. (6): where C s is a constant which includes the physical details of the surface interaction, which may not exactly follow the ponderomotive scaling, and γ is the mean Lorentz factor at a given nearcritical plasma thickness x = d nc . To calculate the total number of substrate electrons, the efficiency η s must be determined. The dependence of the absorption efficiency into hot electrons on the laser intensity is a topic adressed in several works [48][49][50][51] . Here, in order to assure the consistency with the 2D PIC simulation results, we fit the absorption efficiency in the range a 0 = 1-16 from bare solid target simulations. We obtain the linear relation η s = 0.00388 a 0 + 0.04257 (see Fig. 4a). Following the same method adopted for C nc and r c , we fix C s through the fitting of the mean substrate electrons energies calculated from Eq. (16) and the ones retrieved from PIC simulations with a bare target in the intensity range a 0 = 1-16. The retrieved value of 0.52 in 2D, less than one, is consistent with the analysis of ref. 46 . Again, the C s,2D value is 1.5 times the C s,3D due to the higher mean laser envelope amplitude. Finally, we note that the electron population in the near-critical layer is generated directionally in the magnetized plasma and it tends to mix toghether with that of the substrate. Thus, we define a new single population with a mean energy E DLT (d nc ) obtained as a weighthed average of the two: where N tot d nc ð Þ ¼ N s d nc ð ÞþN nc d nc ð Þ is the total number of hot electrons. Since x is the actual independent variable, we normalize the near-critical layer thickness in the same fashion, by introducing The comparison between the results of the model and the PIC simulations, in Fig. 4b, shows also in this case a good agreement for d nc ≲2 ffiffiffi n p l f =λ % 2w 0 =λ, which is another indication that the C nc,2D and r c,2D values can be considered acceptable. The figure indicates that the mean energy of DLT electrons grows when increasing the near-critical layer thickness until a maximum value is reached, at the SF focal length. At this length the energies of both the near-critical layer Fig. 3 Evolution of electron heating in a near-critical layer. a shows the absorption efficiency η nc of the laser energy into the hot electrons of the nearcritical layer for different laser propagation lengths x, normalized to the transparency factor n and the laser wavelength λ. The full lines and points refer to 2D particle-in-cell (PIC) simulations with intensity ranging in a 0 = 2-8 and plasma density in n e /n c = 0.5-2, while the dashed lines refer to our model. b shows the mean energy E nc of the hot electrons of the near-critical layer, normalized to the ponderomotive scaling γ 0 À 1 À Á m e c 2 , for different laser propagation lengths. The full lines and points refer to 2D PIC simulations with intensity ranging in a 0 = 2-8 and plasma density in n e /n c = 0.5-2, while the dashed lines refer to our model. and substrate hot electron populations are at their maxima due to the SF intensity amplification. Furthermore, it should be observed that E DLT tends, for increasing d nc values, to the near-critical layer electrons mean energy E nc . This is due to the high absorption efficiency of the near-critical layer exceeding the one of the substrate electrons, as can be seen by the comparison between Figs. 3a and 4a.
Ion acceleration with the near-critical DLT. In this section, we estimate the maximum ion energy ϵ max using a DLT target. To do so, we exploit a quasi-stationary TNSA model (see "Methods" section), combined with the DLT hot electrons mean energy that we deduced in the previous Section and compare the results with the 2D/3D PIC simulations. Moreover, we discuss the different features in the DLT proton acceleration in 2D and 3D.
To estimate the accelerated ions energy we use the approximate relation ϵ max (26)) for protons, and ϵ max ¼ Zϵ max p =A for ions with Z charge and A mass. Therefore, we need to express the hot electron temperature T h and density n h0 according to our model. T h is related to the electron energy E DLT with a functional dependence determined by the shape of the hot electron distribution function. While different kinds of distribution functions can be plugged into the quasi-static TNSA model 52 , here we consider a perfectly exponential spectrum, by which we have T h = E DLT (d nc ). To calculate n h0 we assume that the electrons are spread uniformly in a 'cylinder' with volume of Lastly, to carry out the proton energy estimation we have to fix theñ free parameter. The parameterñ comes from the adopted model for ion acceleration. In this modelñ represents the density of hot electrons far away from the target, where the electrostatic field driving the acceleration vanishes. Since this quantity does not represent a straightforward physical observable, we leave it as a free parameter to be fitted from simulation results with the bare solid target, d nc = 0. We obtain the value of 1.2 × 10 −3 n c and 5 × 10 −2 n c in 2D and 3D, respectively. Consistently with the physical interpretation ofñ, these values are well below n c .
The resulting maximum proton energy is compared to the 2D PIC data in Fig. 5a, where a remarkable agreement is obtained. The data are represented as a function of the normalized abscissa d nc ¼ ffiffiffi n p d nc =λ, as done in Fig. 4b. The fact that, at a given a 0 and for different values of n, the proton energies lie on almost the same curve is an indication that the proton energy is roughly proportional to the mixed population temperature and it follows the same relativistic normalization. Another indication of this point is that, referring to Fig. 5b, all the points tend to collapse to a self-similar curve when the maximum energy is normalized to the ponderomotive scaling γ 0 −1, similarly to what was done for the near-critical layer electrons (see Fig. 3b).
It should be noticed that the maximum of the self-similar curve is situated at about the initial waist value w 0 /λ. This behaviour is explained by the following considerations: since the proton energy linearly depends on the electrons temperature, it reaches its maximum value when the mean energy of DLT electrons reaches its maximum as well, at the SF focal length. For this reason it is quite straightforward to estimate the optimal normalized thickness for the near-critical layer as d opt nc % ffiffiffi n p l f =λ % w 0 =λ, which is similar to the results obtained in refs. 21,22,26,31,32 .
We also numerically solved the 3D equation set with a finite difference method for a 0 = 4 and the resulting proton energies are compared to 3D PIC simulations in Fig. 6a, b, against the normalized thickness and density, respectively. The model is remarkably accurate in this case as well, even if the n>0:5 cases suffer from a limited error, probably due to the overestimation of the reflectance at relatively high n (see "Methods"). We note that also in the 3D case the largest ϵ max p is obtained at about the SF focal length.
It should be noticed that in 3D the proton energies do not lie on a self-similar curve for different n values, as seen in the 2D case instead. This is due to the form of the equations which depend on n in a different fashion: in particular, Eq. (8) predicts a higher amplitude amplification with respect to 2D, indicating a stronger SF for n approaching 1. This has an effect also on the energy loss equation and the total number of near-critical layer electrons (Eqs. (7) and (10)), where the waist appears raised to the Here the maximum value of the hot electrons mean energy is observed at the Self-Focusing focal length l f , equal to 4.
second power in 3D. In addition, due to Eq. (30), the reflectance has a steeper trend on the normalized density (see "Methods"), suggesting that lower n values allow exploiting more efficiently the pulse energy for hot electrons heating. We therefore expect that an optimal density value exists, where the SF is sufficiently strong to produce high mean energy electrons, yet the reflected part of the pulse is at the same time reduced.
Finally, we compared the 3D model predictions, which should give more realistic results, to available experimental data. Among all the experimental data obtained so far with DLTs, only a limited number of cases satisfy all the hypotheses introduced in our model, namely normal incidence, short pulse duration and uniform near-critical layer (see Table 2). In particular, we take into account the experimental cut-off energies of protons and carbon ions from refs. 26,29 , obtained, respectively, with linear and circular polarization. To include circular polarization in the model we only adjusted the Lorentz factor to γ x ð Þ ¼ in all the calculations. The free parameters were set to the 3D values of Table 1, except forñ ¼ 3 10 À2 n c , fixed by fitting the maximum energy of the protons obtained with bare targets. Since we don't have a realistic 3D fit for the substrate electron absorption η s , we exploited for η s the scaling law presented in   Lower and upper bounds for the model parameters, which respect the model hypotheses. a 0 is the normalized laser amplitude, τ the laser duration, n ¼ n e =γ 0 n c the transparency factor (with n e and n c the electron and critical density, respectively), d nc the near-critical layer thickness, and d s the substrate thickness. l f is the Self-Focusing focal length (see Eq. (4)). λ is the laser wavelength. The slash indicates that the parameter bound is not limited by the model.
ref. 48 , based on experimental data. We also included a 10% error in the reported values of a 0 , τ, w 0 , n e /n c to explicitly take into account experimental errors and obtain a confidence area. Figure 7 represents these results, showing a good agreement with the protons and carbon ions experimental data. This confirms the validity of the 3D model in predicting the ions maximum energies for different species and also in different polarization conditions.
Determination of optimal near-critical layer parameters. In the previous section, we noted that in the more realistic 3D case the cut-off proton energy is sensitive not only to d nc but also to the n parameter. To support this point we calculated ϵ max p with the 3D model, for w 0 = 5 λ and a 0 = 32, as a function of n e /n c and d nc /λ and represented the results in Fig. 8a. The highest predicted proton energies lie on an island with optimal thickness well approximated by the SF length, corresponding to d opt nc ¼ ffiffiffi n p l f =λ % w 0 =λ, and for a limited range of densities. For this reason not only the thickness, but also the density of the near-critical layer must be carefully chosen to optimize the ion acceleration process. In order to find an explicit relation for the optimal near-critical layer parameters, we solve the 3D geometry model Eqs. (1)- (14), (26) and (30) in an approximate analytical way.
First, we observe that Eq. (26) predicts a linear dependence of the proton energy on the hot electron temperature, and weaker logarithmic dependence on the hot electrons density. Thus, the ion energy enhancement factor (defined as the ratio of the cut-off proton energy obtained with the DLT to the one obtained with the standard target) can be roughly estimated with the hot electrons enhancement factor, defined as E DLT d nc ; n e ð Þ=E 0 s , where E 0 s ¼ C s;3D γ 0 À 1 À Á m e c 2 is the hot electron mean energy obtained in the standard bare target case. Thus the DLT temperature equation (Eq. (14)) should be analytically solved. Exploting Eqs. (10) and (12) to re-write the denominator, we find the relation: To solve Eq. (15), an explicit expression for ε p x ð Þ=ε p0 should be written. To do so, we restrict the analysis to the ultra-relativistic limit (a 0 ≫ 1) where the normalized amplitude is proportional to the Lorentz factor (aðxÞ % ffiffi ffi 2 p γðxÞ). Under this approximation, Eq. (7) reduces to This equation can be solved by the variable separation method to obtain an analytical solution: Comparison between cut-off ions energy ϵ max from experimental data 26,29 (points) and our model (continuous line), as a function of the near-critical layer thickness d nc , normalized to the transparency factor n and the laser wavelength λ. The filled area represents the model predictions considering a 10% error in the laser intensity, waist and temporal duration, and in the near-critical layer density. As previously mentioned, the highest temperature of the DLT hot electron population is found at the SF length: By substituting this value into Eq. (17), the pulse residual energy at the SF length, which depends only on the normalized density n, is retrieved: where the R x 0 w x ð Þ=w 0 dx integral was solved explicitly from Eq. (1). The term in the square brackets tends to (l f /x R ) 2 when l f / x R increases (when l f /x R > 2, which is equivalent to n ≥ λ 2 =2w 2 0 , the relative error is under 50%), thus the energy loss can be expressed as Now Eq. (20) can be used to write the enhancement factor as a function of n only; in addition, owing to the fact that η s is often quite low compared to η nc at the SF length (see Figs. 3a and 4a), we can neglect its contribution, which is equivalent to state that E DLT tends to E nc at the SF focal length, as observed in Fig. 4b: Furthermore, we exploit that the reflectance R 3D approaches zero when the transparency factor is sufficiently low, approximately n<1=4 (this assumption is verified a posteriori in "Methods"). To find the normalized density which optimizes the enhancement factor, we calculate the derivative of Eq. (21) and impose it to zero. The derivative numerator is proportional to n 3=2 þ 3ϖ n 1=2 À 4ϖ=ρ, where ϖ ¼ 3λ 2 =π 2 w 2 0 and ρ ¼ C nc;3D r 2 c;3D w 0 = ffiffi ffi 2 p πτc are constants. Since ϖ approaches zero and n is assumed low, the term 3ϖ n 1=2 is an infinitesimal of higher order than 4ϖ=ρ and can be neglected (also this approximation is verified in "Methods") in order to easily find the zero of the derivative, which is as follows: Equations (18) where the numerical values are calculated from the free parameters of Table 1 and the intial laser waist and the temporal duration are expressed as the FWHM over the intensity. As a final step, the optimal near-critical parameters can be used to determine the value of the optimized enhancement factor. Substituting n opt in Eq. (15) The comparison between the optimal values analytically estimated by Eqs. (23)- (25) with the numerical solution of the 3D model are satisfactory, as shown in Fig. 8b.

Discussion
We derived explicit relations for the optimal near-critical layer parameters which depend on the pulse waist, temporal duration and intensity. In particular, we obtained from Eqs. (23) and (24) that a larger waist requires a thicker and less dense near-critical layer in order to efficiently focus the laser, keeping its energy loss limited, and heat the electrons to higher energies. An opposite behaviour is observed for the pulse temporal duration: the pulse energy increases linearly with τ, at fixed intensity, which means that the laser can be more strongly focused using higher densities and lower thicknesses, without an excessive energy loss. Equation (23) predicts also that n e should be increased as the laser intensity increases (since γ 0 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 þ a 2 0 =2 p ), because the relativistic effects make the plasma more transparent.
The proton energy enhancement factor (Eq. (25)) is quite straightforward to interpret: when the square brackets term approaches the unity (for higher temporal durations), E opt DLT tends to a constant, equal to about 3C nc;3D =2C s;3D % 4:6. This can be interpreted as the ratio of the super-ponderomotive hot electron temperature of the near-critical plasma to the ponderomotive energy for a bare solid foil; for this reason, within the validity ranges of the proposed model, the maximum enhancement value actually remains invariant with respect to the laser parameters. The obtained enhancement factor appears quite reasonable in light of the published experimental results, since enhancements in the range 1.5-3 have been reported in the literature [24][25][26][27][28][29][30][31] . Equations (23)- (25) can be regarded as a useful tool to carry out experiments which aim at optimizing the DLT performances and to scale the results to other laser sources.
The validity of our theoretical calculations is of course limited by the adopted initial assumptions and approximations: making reference to Table 2, the model is able to make accurate predictions when the pulse amplitude is relativistic, yet, in order to neglect the plasma ions motion and the synchrotron-like radion, the pulse duration is short enough (<100 fs) and the intensity sufficiently low (a 0 < 50). Moreover, as previously explained, our waist evolution equation, describing the self-focusing, is reasonably valid for near-critical density (0:05< n<1) and for nearcritical layer thickness lower than about two-times the selffocusing length. Since we are describing TNSA acceleration, the substrate should not be destroyed by the laser. A lower limit for its thickness can be given by the optimal thickness for RPA light sail, d s ¼ a 0 λ=π n c =n e ð Þ s 4 . Approaching this value we expect radiation pressure to distort the ions spectrum and ultimately, under this threshold, to disrupt the target and suppress ion acceleration.
In addition, we initially assumed linear polarization, homogeneous plasma and normal incidence. Nevertheless, our model lets us gain insights on non-ideal configurations as well. Firstly, we point out that the laser interaction with the near-critical plasma and, ultimately, the enhancement factor should be weakly dependent on the pulse polarization. It was demonstrated in a previous simulation work 38 that P and C polarized pulses produce, when the first layer is sufficiently transparent ( n≲0:3), similar electron temperatures and thus we expect comparable C nc,3D values in this density range. This is confirmed by the good agreement between experimental data and model results observed in Fig. 7. The independence of the DLT proton energy on the polarization was also observed experimentally in refs. 27,28 . Nonetheless, it is worth noting that in C polarization, the γ 0 factor differs from the one in linear polarization of a factor about ffiffi ffi 2 p . Secondly, our analysis allows making some considerations about the near-critical plasma homogeneity effects: PIC simulations works 38,40 reported that a nanostructured plasma, with an inhomogeneity scale greater than the laser wavelength, can suppress the DLA resonant mechanism, with a reduction in the electron temperature. Moreover, it has been shown that, when n≲0:3, the nanostructure is capable of increasing the number of mildly energetic electrons, and keeping the total pulse energy absorption similar to the homogeneous case. To take into account these effects a simple corrective factor α ns could be introduced (with 0 < α ns < 1) which adapts the near-critical layer hot electron temperature (Eq. (11)) and number (Eq. (10)) to the nanostructured case: E ns ðxÞ ¼ α ns E nc ðxÞ and N ns ðxÞ ¼ N nc ðxÞ=α ns . We emphasize that this point is beyond the scope of this work and it should be the aim of a deeper analysis.
Thirdly, we believe that the variation of the incidence angle from the normal could be the most crucial issue, with respect to the other two. Indeed, if the pulse interacts with a tilted nearcritical plasma, the self-focusing axial simmetry is destroyed, inducing other effects: such as the pulse refraction and a mismatch in the angular distributions of the hot electrons populations (the near-critical ones should be accelerated in the laser direction while the substrate ones along the normal of the target, eventually separating at the rear of the substrate).
In conclusion we have quantitatively described the essential aspects of ultra-intense laser interaction with near-critical DLTs, characterizing the pulse attenuation, the SF intensity amplification and the hot electrons populations mean energy and total number. The free parameters adopted in this theoretical description were fitted from 2D and 3D PIC simulation results, finding a reasonable agreement in both the trends and the absolute values of all the observed quantities. We could have let the free parameters vary depending on the specific configuration, however, we found a good agreement even by fixing them to the reported values once and for all. We coupled this model with a well-established quasi-stationary TNSA model in order to estimate the maximum energy of the accelerated ions for different near-critical layer densities and thicknesses. We observed both in the 2D/3D model and in 2D/3D PIC simulations a self-similar behaviour in the proton energy with respect to the normalized thickness d nc ¼ ffiffiffi n p d nc =λ, with a maximum at the self-focusing length d opt nc % ffiffiffi n p l f =λ % w 0 =λ. We used the 2D version of our model to validate the hypotheses of the model itself. We did this by comparing the model results with 2D PIC simulations results for a large number of target densities and thicknesses. On the other hand, the 3D version of the model is intended as a convenient tool for the interpretation of 3D simulations and experimental results and to guide their design. Finally, the explicit model solution, valid in the ultra-relativistic limit, can be exploited to explicitly calculate the optimal nearcritical layer density and thickness to maximize the proton energy. We showed that the retrieved values depend directly on the specific laser source used for the acceleration experiment, in particular on its intensity, its focal spot and its temporal duration. Moreover, we derived a theoretical maximum enhacement of the ion energy that can be obtained using a DLT with respect to a standard foil. We found that this only depends on the ratio of the near-critical super-ponderomotive electron temperature to the ponderomotive energy of the electrons in the substrate. We also discussed the validity ranges of our model and we suggested possible ways to widen them.
Our results provide an effective tool to design near-critical DLTs that are optimized for the purpose of laser energy conversion into hot electrons, hence for laser-driven ion acceleration. At this regard, we obtained a simple recipe for the optimal DLT properties. These properties (i.e. density and thickness of the near-critical layer) can be relatively easy controlled and manipulated during the DLT production phase, which usually relies on advanced synthesis techniques 28,30,[53][54][55][56] . Certainly, optimizing a DLT-based laser-driven ion source is also of great interest for the potential applications of TNSA, since it would allow to obtain higher ion energies without having to improve the laser system. Lastly, our theoretical approach could be used in other contexts than TNSA, for instance the DLT parameters could be suitably tuned to optimize other acceleration mechanisms (e.g. magnetic vortex acceleration with free-standing nearcritical plasmas or radiation pressure acceleration with ultrathin substrate DLTs) or even other physical processes (e.g. photon sources by synchrotron-like emission).

Methods
Particle-in-cell simulations. A total of 58 simulations were performed with the open-source, massively parallel code piccante 57 . The laser pulse had an idealized cos 2 temporal profile of the fields (to approximate an ultra high contrast laser) and a Gaussian transverse profile, and it was linearly polarization with the electric field lying in the simulation plane (P-polarization). The temporal duration was 15 λ/c (FWHM of the fields). The intensity was varied between a 0 = 2 and a 0 = 8 at fixed normal incidence. These parameters, if scaled to Ti:Sapphire lasers (λ = 800 nm), correspond to a 28.5 fs FWHM pulse, attaining a peak intensity in the range 8:7 10 18 W=cm 2 <I<1:4 10 20 W=cm 2 which is found in small-medium scale super-intense laser facilities 58 .
For the 2D simulations with the near-critical uniform plasma only (used to study the SF, the pulse energy loss and amplification), a 4 λ waist (corresponding to a spotsize FWHM about 3.7 μm) and a 100 λ × 60 λ box were used, with a resolution of 20 points per wavelength. The plasma filled a 50 λ × 60 λ box starting from x = 0; 9 macro-electrons per cell were used. We simulated the following densities: 0:5 n c ; 1 n c ; 2 n c . The laser was focused at the vacuum-plasma boundary.
For the 2D simulations with the DLT (used to study the DLT hot electrons mean energy and the proton cut-off energy), a 4 λ waist and a 200 λ × 120 λ box were used, with a resolution of 64 points per wavelength. The near-critical layer, starting at x = 0, had 10 macro-electrons per cell with densities of 0:5 n c ; 1 n c ; 2 n c and thickness varied in the range 0.5-32λ; 64 macro-electrons were used for the solid density layer, with density fixed to 64 n c and thickness 0.5 λ. The laser was focused at the near-critical layer-substrate boundary. The laser was focused at the near-critical layer-substrate boundary. It is worth considering that a laser waist of 4 λ implies a Rayleigh length of~50 λ, which is larger than the longest near-critical layer considered in our simulation campaign. For this reason, within the parameter range that we have explored, it is reasonable to expect a small dependence of the simulation results on the position of the focal spot of the beam.
For the 3D simulations with the DLT, a 5 λ waist and a 100 λ × 60 λ × 60 λ box were used, with a resolution of 20 points per wavelength. The near-critical layer, starting at x = 0, had 10 macro-electrons per cell with densities of 0.5, 1, 1.5, 1.8, 2.3, 3 n c and thickness varied in the range 4-24 λ; 40 macro-electrons were used for the solid density layer, with density fixed to 40 n c and thickness 0.5 λ. The laser was focused at the vacuum-plasma boundary.
We performed additional convergence tests to ensure that both the spatial resolution and the number of macro-electrons per cell were adequate for our physical scenario. The case a0 = 4, P-polarization, homogeneous plasma was selected and simulations with resolution of 20 and 40 points per wavelength were carried out. Negligible differences were observed for the electron energy spectra, the absorption efficiency and the proton energy in these cases.
In all the simulated cases the plasma was fully pre-ionized and the charge/mass ratio of the ions was 0.5 (e.g. C 6+ ). The electron population was initialized with a small temperature (few eVs) to avoid numerical artefacts. In order to verify that this choice did not affect our results, we performed an additional simulation with a higher initial temperature (~1 keV) in a reduced box, where no differences in the simulation were observed.
Several analyses were performed on the PIC simulations: the normalized amplitude of the pulse was calculated as ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where the propagation length x is calculated as the position of the amplitude half maximum. In a similar way the pulse waist was calculated as the radial width at 1/e threshold of the fields. The total energy of the laser was approximated by the integral over the whole box of the electromagnetic energy density relative to the B z and the E y component; the reflectance was calculated as the ratio of the reflected energy to the total one. The electron mean energies were calculated on the whole box excluding the non-relativistic electrons, namely the ones with energy lower than m e c 2 . Since in 2D simulations the proton energy does not saturate we set the time at which the maximal proton energy is calculated by imposing the energy time derivative to a constant 59 .
Quasi-stationary TNSA model. Several theoretical models have been proposed to estimate some of the most important features of laser accelerated ions. Three main branches of TNSA models exist, defined by a different treatment reserved to the ion dynamic: quasi-stationary 60 , dynamical 61 , and hybrid 62 . A model that provides a good agreement with the experimental results in a wide range of conditions is the quasi-stationary description 60,[63][64][65][66] . This model gives an estimation of the ions cutoff energy in TNSA, which reads as follows: where T h is the hot electron distribution temperature, φ * is the normalized p dp. The normalized potential is retrieved by solving the implicit equationñ I φ * ð Þ ζκ 1 ζ ð Þ e φ * ¼ n h0 , with n h0 the hot electron density, andñ a normalization constant of the hot electron distribution function that is used as a free parameter 66 . Note that the last approximated part of Eq. (27) is verified under the φ * ≫ 1 condition, which has the physical meaning of imposing that the hot electrons distribution cut-off energy is much higher than its temperature, which is a very common and reasonable condition.
Reflectance calculation in two and three dimensions. It is well known that an electromagnetic wave is reflected by an overcritical plasma while it is transmitted in an undercritical medium. In our case we can have a mixed behaviour since a plasma can be at the same time relativistically undercritical (n e /γ 0 n c < 1), but classically overcritical (n e /n c > 1). Indeed, if we consider a Gaussian pulse amplitude envelope in 2D, a t; y ð Þ ¼ a 0 e Àt 2 =τ 2 e Ày 2 =w 2 0 , near the laser peak the electrons move at relativistic speed allowing the pulse to propagate; while in correspondance of the envelope tails, the electrons move with a not-relativistic quiver motion eventually resulting in an overcritical reflecting plasma.
Making reference to Fig. 9a, to calculate the reflectance in 2D, we firstly find the threshold of this process given by the condition n e =γ t; y ð Þn c ¼ 1. Roughly approximating γ t; y ð Þ % a t; y ð Þ= ffiffi ffi 2 p we can rewrite it as n e γ 0 n c With a change of variables (t=τ ¼ ξ and y=w 0 ¼ χ) and taking the natural logarithm of the latter equation we obtain To calculate the fraction of energy which is not reflected, we have to integrate the electromagnetic energy density and we use the more convenient polar coordinates r 2 ¼ ξ 2 þ χ 2 , since Eq. (28) represents a circumference: This relation is in agreement with the trend given by 2D PIC simulation as shown in Fig. 9b even if it underestimates the absolute values at high n, when a 0 is low (since we have approximated the Lorentz factor in the ultra-relativistic limit).
In a similar way we can also evaluate the transmittance in 3D with the following integral: Approximations validity. We discuss the range of validity of the approximations underlying Eqs. (22)- (25). In order of appearance we assumed the following. n>λ 2 =2w 2 0 : substituting this inequality into Eq. (22), we get the condition τ>π 2 r 2 c;3D C nc;3D λ=12 ffiffiffiffiffi 2π p c % 1:6λ=c. This corresponds to a FWHM temporal duration longer than 5 fs which is generally verified for nearly all high intensity laser systems.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Code availability
The codes generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. Fig. 9 Calculation of the reflectance R. a shows the level plot of the amplitude gaussian in normalized units ξ ¼ t=τ and χ ¼ y=w 0 . The dashed black line marks a general threshold À log n obtained in Eq. (28); the dashed part of the plot represents the tails of the pulse which are reflected by the overcritical plasma. b represents the reflectance as obtained from Eq. (29) (2D, full line) and Eq. (30) (3D, dashed line); to be compared to 2D/3D particle-in-cell (PIC) simulation data (points).