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

Abstract

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.

Introduction

Laser-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⁶, relativistic induced transparency^7,8, collision-less shock acceleration⁹, magnetic vortex acceleration¹⁰). Yet, Target Normal Sheath Acceleration (TNSA)¹¹ 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¹⁶) 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\omega ^2{\mathrm{/}}4\pi 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 self-focusing (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³⁴ 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³⁸. 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 $\bar n = n_e{\mathrm{/}}\gamma _0n_c$ (or the normalized plasma density), where $\gamma _0 = \sqrt {1 + a_0^2{\mathrm{/}}\wp }$ is the mean Lorentz factor of the electron motion due to a laser with normalized amplitude a₀ and linear ($\wp = 2$) or circular ($\wp = 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 self-focusing, 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⁴¹ 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 laser-driven 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, $\bar n = n_e{\mathrm{/}}\gamma _0n_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.

Fig. 1: Laser interaction with a near-critical plasma.

a shows the transverse component of the magnetic field (the polarization plane is xy) of a super-intense laser with a₀ = 8, propagating inside a uniform near-critical plasma with n_e/n_c = 1 ($\bar n = 0.17$) at the time 20 λ/c after the start of the interaction. The self-focusing of the pulse is evident from the figure. b shows the hot electrons density (with energy higher then m_ec²) at the time 20 λ/c. c shows the transverse magnetic field at the time 36 λ/c after the start of the interaction. The magnetization of the channel should be noticed. Figure (d) shows the schematization of our model.

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 ${\cal{D}}$-dimensional geometry (${\cal{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:

$$\begin{array}{*{20}{c}} {w_m = \frac{\lambda }{\pi }\frac{1}{{\sqrt {\bar n} }}}, \end{array}$$

(1)

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:

$$\begin{array}{*{20}{c}} {w\left( x \right) = w_m\sqrt {1 + \left( {\frac{{x - l_f}}{{x_R}}} \right)^2} }, \end{array}$$

(2)

$$\begin{array}{*{20}{c}} {x_R = \frac{\pi }{\lambda }w_m^2 = \frac{\lambda }{{\pi \bar n}}}, \end{array}$$

(3)

$$\begin{array}{*{20}{c}} {l_f = \sqrt {x_R\left( {\frac{\pi }{\lambda }w_0^2 - x_R} \right)} = \sqrt {\frac{{w_0^2 - w_m^2}}{{\bar n}}} }, \end{array}$$

(4)

where w₀ 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⁴², 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 $\bar n \gtrsim 0.05$ (under this threshold we don’t expect a focusing-defocusing behaviour, but a stable propagation⁴³) and within a pulse propagation length of the order of the SF focal length ($x \lesssim 2l_f$). At longer distances other phenomena are induced, as laser-plasma non-linear instabilities⁴² and filamentation⁴⁴, which distort the ideal Gaussian shape of the 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.

Fig. 2: Evolution of laser beam parameters.

The figure shows, respectively, the laser waist w (a, b), the pulse energy ${\it{\varepsilon _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 $\bar 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₀ = 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 ${\cal{R}}_{\cal{D}}$ to the calculated energy.

It is worth mentioning that $l_f \approx w_0/\sqrt {\bar n}$ when $w_0 \gg w_m$ (which is a reasonable approximation). Under this condition, Eq. (2) can be approximated by

$$\begin{array}{*{20}{c}} {\frac{{w\left( {\bar x} \right)}}{\lambda } \approx \sqrt {\frac{1}{{\pi ^2\bar n}} + \left( {\bar x - \frac{{w_0}}{\lambda }} \right)^2} }. \end{array}$$

(5)

Here we defined the relativistically normalized space variable $\bar x = \sqrt {\bar n} x/\lambda$ (defined in ref. ³⁹ 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 $\bar x$ variable and leads to self-similar curves for equal initial waist w₀/λ and for constant transparency factor $\bar 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. ⁴⁵:

$$\begin{array}{*{20}{c}} {{\rm{d}}\varepsilon _p = - V_{{\cal{D}} - 1}R_{ch}(x)^{{\cal{D}} - 1}n_eC_{nc}\left( {\gamma \left( x \right) - 1} \right)m_ec^2{\rm{d}}x}, \end{array}$$

(6)

where R_ch is the plasma channel radius, $V_{{\cal{D}} - 1}R_{ch}(x)^{{\cal{D}} - 1}$ is the plasma channel section, $V_{{\cal{D}} - 1} = \pi ^{\frac{{{\cal{D}} - 1}}{2}}{\mathrm{/}}\Gamma \left( {\frac{{{\cal{D}} + 1}}{2}} \right)$ is the volume of a $\left( {{\cal{D}} - 1} \right)$-dimensional ipersphere with unitary radius (V₀ = 1, V₁ = 2, V₂ = π), and Γ the Euler gamma function. $\gamma \left( x \right) = \sqrt {1 + a\left( x \right)^2{\mathrm{/}}2}$ 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}\left( {\gamma \left( x \right) - 1} \right)m_ec^2$. Thus C_nc absorbs the details of the electron heating process⁴⁶, 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 ${\cal{D}}$ dimensions as $\varepsilon _{p0} = \pi ^{{\cal{D}}/2}2^{ - {\cal{D}}/2 - 1}m_ec^2n_ca_0^2w_0^{{\cal{D}} - 1}\tau c$⁴⁷, 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 $\tau _{{{FWHM}}}^I = \sqrt {2{\mathrm{log}}2} \tau$ and $w_{{{FWHM}}}^I = \sqrt {2{\mathrm{log}}2} w_0$). Then, the normalized energy loss assumes the following form:

$$\begin{array}{*{20}{c}} {\frac{1}{{\varepsilon _{p0}}}\frac{{{\rm{d}}\varepsilon _p(x)}}{{{\rm{d}}x}} = - 2\left( {\frac{2}{\pi }} \right)^{\frac{{\cal{D}}}{2}}V_{{\cal{D}} - 1}C_{nc}\frac{1}{{\tau c}}\frac{{n_e}}{{a_0n_c}}\frac{{\gamma \left( x \right) - 1}}{{a_0}}\left( {\frac{{r_cw\left( x \right)}}{{w_0}}} \right)^{{\cal{D}} - 1}}, \end{array}$$

(7)

where we introduced the ratio of the plasma channel radius to the waist $r_c = R_{ch}(x)/w\left( x \right)$, 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:

$$\begin{array}{*{20}{c}} {\frac{{a\left( x \right)}}{{a_0}} = \sqrt {\frac{{\varepsilon _p(x){\mathrm{/}}\varepsilon _{p0}}}{{\left( {w\left( x \right){\mathrm{/}}w_0} \right)^{{\cal{D}} - 1}}}} }. \end{array}$$

(8)

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²¹, 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 $\varepsilon _p(0){\mathrm{/}}\varepsilon _{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 ${\cal{R}}_{\cal{D}}$ (Eqs. (29) and (30) in “Methods”), which are used to set the initial condition $\varepsilon _p(0){\mathrm{/}}\varepsilon _{p0} = 1 - {\cal{R}}_{\cal{D}}$. Since Eqs. (7) and (8) depend on the waist evolution and thus on the self-similar variable $\bar 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} \gtrsim 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 amplitude is lower in 3D than in 2D. This dimensionality effect has been reported also in ref. ⁴⁶.

Table 1 Model free parameters.

Now that we have a model for the pulse propagation in a near-critical 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:

$$\begin{array}{*{20}{c}} {\eta _{nc}\left( x \right) = 1 - \frac{{\varepsilon _p\left( x \right)}}{{\varepsilon _{p0}}} - {\cal{R}}_{\cal{D}}}. \end{array}$$

(9)

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₀ < 50³⁸). We confirm that these approximations hold by comparing the calculated η_nc(x) with the 2D PIC simulations data (see Fig. 3a).

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 near-critical layer for different laser propagation lengths x, normalized to the transparency factor $\bar n$ and the laser wavelength λ. The full lines and points refer to 2D particle-in-cell (PIC) simulations with intensity ranging in a₀ = 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 $\left( {\gamma _0 - 1} \right)m_ec^2$, for different laser propagation lengths. The full lines and points refer to 2D PIC simulations with intensity ranging in a₀ = 2–8 and plasma density in n_e/n_c = 0.5–2, while the dashed lines refer to our model.

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:

$$\begin{array}{*{20}{c}} {\frac{{{\rm{d}}N_{nc}\left( x \right)}}{{{\rm{d}}x}} = V_{{\cal{D}} - 1}R_{ch}\left( x \right)^{{\cal{D}} - 1}n_e = V_{{\cal{D}} - 1}n_e\left( {r_cw\left( x \right)} \right)^{{\cal{D}} - 1}}. \end{array}$$

(10)

Exploting Eqs. (9) and (10), the mean hot electron energy E_nc is retrieved:

$$\begin{array}{*{20}{c}} {E_{nc}\left( x \right) = \frac{{\eta _{nc}\left( x \right)\varepsilon _{p0}}}{{N_{nc}\left( x \right)}}} \end{array}$$

(11)

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 semi-infinite 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

$$\begin{array}{*{20}{c}} {\eta _s = \frac{{N_sE_s}}{{\varepsilon _p(d_{nc})}}}, \end{array}$$

(12)

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):

$$\begin{array}{*{20}{c}} {E_s\left( {d_{nc}} \right) = C_s\left( {\gamma \left( {d_{nc}} \right) - 1} \right)m_ec^2}, \end{array}$$

(13)

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 near-critical 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₀ = 1–16 from bare solid target simulations. We obtain the linear relation η_s = 0.00388 a₀+ 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₀ = 1–16. The retrieved value of 0.52 in 2D, less than one, is consistent with the analysis of ref. ⁴⁶. Again, the C_s,2D value is 1.5 times the C_s,3D due to the higher mean laser envelope amplitude.

Fig. 4: Electron heating in the Double-Layer-Target (DLT).

a shows the absorption efficiency η_s of the laser energy into hot electrons, in the bare target case, retrieved from 2D particle-in-cell (PIC) simulations, as a function of the laser intensity a₀. The dashed line refers to the fit η_s = 0.00388 a₀+ 0.04257. b shows the DLT hot electron mean energy E_DLT, with a near-critical layer of thickness d_nc normalized to the transparency factor $\bar n$ and the laser wavelength λ. The full points refer to 2D PIC simulations with intensity a₀ = 8 and plasma density in n_e/n_c = 0.5–2, while the lines refer to the hot electron energies calculated by our model relative to the substrate (dashed lines), the near-critical layer (dotted lines) and the weighted average (continuous lines). Here the maximum value of the hot electrons mean energy is observed at the Self-Focusing focal length l_f, equal to 4.

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:

$$\begin{array}{*{20}{c}} {E_{DLT}\left( {d_{nc}} \right) = \frac{{\eta _s\varepsilon _p(d_{nc}) + \eta _{nc}\left( {d_{nc}} \right)\varepsilon _{p0}}}{{N_{tot}\left( {d_{nc}} \right)}}}, \end{array}$$

(14)

where $N_{tot}\left( {d_{nc}} \right) = N_s\left( {d_{nc}} \right) + N_{nc}\left( {d_{nc}} \right)$ is the total number of hot electrons. Since $\bar x$ is the actual independent variable, we normalize the near-critical layer thickness in the same fashion, by introducing $\bar d_{nc} = \sqrt {\bar n} d_{nc}{\mathrm{/}}\lambda$. The comparison between the results of the model and the PIC simulations, in Fig. 4b, shows also in this case a good agreement for $\bar d_{nc} \lesssim 2\sqrt {\bar n} l_f{\mathrm{/}}\lambda \approx 2w_0{\mathrm{/}}\lambda$, 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 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 $\bar 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 $\it{\epsilon ^{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 ${\it{\epsilon }}_p^{max} = T_h\left[ {{\mathrm{log}}\left( {n_{h0}{\mathrm{/}}\tilde n} \right) - 1} \right]$ (Eq. (26)) for protons, and ${\it{\epsilon }}^{max} = Z{\it{\epsilon }}_p^{max}{\mathrm{/}}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⁵², 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 $V_{{\cal{D}} - 1}w_0^{{\cal{D}} - 1}d_s$, where d_s is the substrate thickness: $n_{h0} = N_{tot}\left( {d_{nc}} \right){\mathrm{/}}V_{{\cal{D}} - 1}w_0^{{\cal{D}} - 1}d_s$. Lastly, to carry out the proton energy estimation we have to fix the $\tilde n$ free parameter. The parameter $\tilde n$ comes from the adopted model for ion acceleration. In this model $\tilde n$ 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⁻³ n_c and 5 × 10⁻² n_c in 2D and 3D, respectively. Consistently with the physical interpretation of $\tilde n$, 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 $\bar d_{nc} = \sqrt {\bar n} d_{nc}{\mathrm{/}}\lambda$, as done in Fig. 4b. The fact that, at a given a₀ and for different values of $\bar 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 γ₀−1, similarly to what was done for the near-critical layer electrons (see Fig. 3b).

Fig. 5: Proton cut-off energy in two dimensions: comparison between model and simulations.

Cut-off proton energy ${\it{\epsilon }}_p^{max}$ obtained in 2D particle-in-cell (PIC) simulations with fixed substrate thickness and spot size and variable near-critical layer thickness d_nc (normalized to the transparency factor $\bar n$ and the laser wavelength λ) and density n_e, with intensity ranging in a₀= 2–8. In (a) the simulation results (points) are compared with the cut-off proton energy predicted by our model (dashed lines). In (b) the PIC proton energies are normalized to the ponderomotive scaling $(\gamma _0 - 1)m_ec^2$. Here the maximum value of the proton energy is observed at the Self-Focusing focal length l_f, equal to 4.

It should be noticed that the maximum of the self-similar curve is situated at about the initial waist value w₀/λ. 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 $\bar d_{nc}^{opt} \approx \sqrt {\bar n} l_f{\mathrm{/}}\lambda \approx w_0{\mathrm{/}}\lambda$, 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₀ = 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 $\bar n > 0.5$ cases suffer from a limited error, probably due to the overestimation of the reflectance at relatively high $\bar n$ (see “Methods”). We note that also in the 3D case the largest ${\it{\epsilon }}_p^{max}$ is obtained at about the SF focal length.

Fig. 6: Proton cut-off energy in three dimensions: comparison between model and simulations.

The cut-off proton energy ${\it{\epsilon }}_p^{max}$ obtained by 3D particle-in-cell (PIC) simulations are compared to the predictions of our model. The substrate thickness, spot size and intensity are fixed (a₀ = 4), while the near-critical layer thickness d_nc, normalized to the transparency factor $\bar n$ and the laser wavelength λ (a) and density n_e, normalized with the transparency factor $\bar n$ (b) are varied.

It should be noticed that in 3D the proton energies do not lie on a self-similar curve for different $\bar n$ values, as seen in the 2D case instead. This is due to the form of the equations which depend on $\bar n$ in a different fashion: in particular, Eq. (8) predicts a higher amplitude amplification with respect to 2D, indicating a stronger SF for $\bar 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 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 $\bar 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 $\gamma \left( x \right) = \sqrt {1 + a\left( x \right)^2}$ in all the calculations. The free parameters were set to the 3D values of Table 1, except for $\tilde n = 3 \times 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 ref. ⁴⁸, based on experimental data. We also included a 10% error in the reported values of a₀, τ, w₀, 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.

Table 2 Model bounds.

Fig. 7: Ions cut-off energy: comparison between model and experiments.

Comparison between cut-off ions energy ${\it{\epsilon }}^{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 $\bar 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.

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 $\bar n$ parameter. To support this point we calculated $\epsilon _p^{max }$ with the 3D model, for w₀ = 5 λ and a₀ = 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 $\bar d_{nc}^{opt} = \sqrt {\bar n} l_f{\mathrm{/}}\lambda \approx w_0{\mathrm{/}}\lambda$, 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.

Fig. 8: Optimization of Double-Layer-Target parameters.

a shows the cut-off proton energy ${\it{\epsilon }}_{p,DLT}^{max}$, calculated solving numerically the 3D model from Eqs. (1)–(14), 26 and (30) (with a₀ = 32, w₀ = 5 λ) as a function of the density n_e/n_c and thickness d_nc/λ of the near-critical layer; on top of the image the Self-Focusing (SF) focal length is represented (Eq. (18)). b shows that the optimal thickness $d_{nc}^{opt}{\mathrm{/}}\lambda$ (solid blue curve) as a function of the density, calculated numerically from the data of Figure (a), superimposes to the SF focal length l_f (dashed blue curve), calculated by Eq. (18); the orange solid curve represents the enhancement factor calculated by Eq. (21) to be compared to the enhancement factor ${\it{\epsilon }}_{p,DLT}^{max}{\mathrm{/}}{\it{\epsilon }}_{p,bare}^{max}$ (orange dashed curve), calculated numerically from Figure (a); the orange diamond represents the optimal density value and the relative enhancement factor obtained by Eqs. (23) and (25).

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}\left( {d_{nc},n_e} \right){\mathrm{/}}E_s^0$, where $E_s^0 = C_{s,3D}\left( {\gamma _0 - 1} \right)m_ec^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:

$${\frac{{E_{DLT}\left( x \right)}}{{E_s^0}} = \frac{{1 + \left( {\eta _s - 1} \right)\frac{{\varepsilon _p\left( x \right)}}{{\varepsilon _{p0}}} - {\cal{R}}_{3D}}}{{\eta _s\sqrt {\frac{{\varepsilon _p\left( x \right)}}{{\varepsilon _{p0}}}} \frac{{w\left( x \right)}}{{w_0}} + \frac{{2\sqrt 2 C_{s,3D}r_{c,3D}^2}}{{\sqrt \pi }}\frac{{\bar n}}{{\tau c}}{\int}_0^x {\left( {\frac{{w\left( {x\prime } \right)}}{{w_0}}} \right)^2} {\rm{d}}x\prime }}}. $$

(15)

To solve Eq. (15), an explicit expression for $\varepsilon _p\left( x \right){\mathrm{/}}\varepsilon _{p0}$ should be written. To do so, we restrict the analysis to the ultra-relativistic limit (a₀ ≫ 1) where the normalized amplitude is proportional to the Lorentz factor ($a(x) \approx \sqrt 2 \gamma (x)$). Under this approximation, Eq. (7) reduces to

$$\begin{array}{*{20}{c}} {\frac{1}{{\varepsilon _{p0}}}\frac{{{\rm{d}}\varepsilon _p\left( x \right)}}{{{\rm{d}}x}} \approx - \frac{{2\sqrt 2 C_{nc,3D}r_{c,3D}^2}}{{\sqrt \pi }}\frac{{\bar n}}{{\tau c}}\frac{{w\left( x \right)}}{{w_0}}\sqrt {\frac{{\varepsilon _p\left( x \right)}}{{\varepsilon _{p0}}}} }. \end{array}$$

(16)

This equation can be solved by the variable separation method to obtain an analytical solution:

$$\begin{array}{*{20}{c}} {\frac{{\varepsilon _p\left( x \right)}}{{\varepsilon _{p0}}} \approx \left( {\sqrt {1 - {\cal{R}}_{3D}} - \frac{{\sqrt 2 C_{nc,3D}r_{c,3D}^2}}{{\sqrt \pi }}\frac{{\bar n}}{{\tau c}}\mathop {\int}\limits_0^x {\frac{{w\left( {x\prime } \right)}}{{w_0}}} {\rm{d}}x\prime } \right)^2}. \end{array}$$

(17)

As previously mentioned, the highest temperature of the DLT hot electron population is found at the SF length:

$$\begin{array}{*{20}{c}} {\bar d_{nc}^{opt} = \sqrt {\bar n} \frac{{l_f}}{\lambda } \approx \frac{{w_0}}{\lambda }}. \end{array}$$

(18)

By substituting this value into Eq. (17), the pulse residual energy at the SF length, which depends only on the normalized density $\bar n$, is retrieved:

$$\begin{array}{*{20}{c}} \begin{array}{c}\frac{{\varepsilon _p\left( {\bar n,\bar d_{nc}^{opt}} \right)}}{{\varepsilon _{p0}}} \approx \left( {\sqrt {1 - {\cal{R}}_{3D}} - \frac{{\sqrt 2 C_{nc,3D}r_{c,3D}^2}}{{\sqrt \pi }}\frac{{\bar n}}{{\tau c}}\frac{{w_mx_R}}{{2w_0}}\left[ {\frac{{l_f}}{{x_R}}\sqrt {1 + \left( {\frac{{l_f}}{{x_R}}} \right)^2} } \right.} \right.\\ \hskip -82pt\left. {\left. { + \sinh ^{ - 1}\left( {\frac{{l_f}}{{x_R}}} \right)} \right]} \right)^2\end{array} \end{array}$$

(19)

where the ${\int}_0^x w \left( x \right){\mathrm{/}}w_0dx$ integral was solved explicitly from Eq. (1). The term in the square brackets tends to (l_f/x_R)² when l_f/x_R increases (when l_f/x_R > 2, which is equivalent to $\bar n \ge {\uplambda}^2{\mathrm{/}}2w_0^2$, the relative error is under 50%), thus the energy loss can be expressed as

$$\begin{array}{*{20}{c}} {\frac{{\varepsilon _p\left( {\bar n,\bar d_{nc}^{opt}} \right)}}{{\varepsilon _{p0}}} \approx \left( {\sqrt {1 - {\cal{R}}_{3D}} - \frac{{C_{nc,3D}r_{c,3D}^2}}{{\sqrt {2\pi } }}\frac{{\sqrt {\bar n} w_0}}{{\tau c}}} \right)^2}. \end{array}$$

(20)

Now Eq. (20) can be used to write the enhancement factor as a function of $\bar 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:

$$\frac{{E_{DLT}\left( {\bar n,\bar d_{nc}^{opt}} \right)}}{{E_s^0}} \approx \frac{{E_{nc}\left( {\bar n,\bar d_{nc}^{opt}} \right)}}{{E_s^0}}\\ = \frac{{1 - \left( {\sqrt {1 - {\cal{R}}_{3D}} - \frac{{C_{nc,3D}r_{c,3D}^2}}{{\sqrt {2\pi } }}\frac{{\sqrt {\bar n} w_0}}{{\tau c}}} \right)^2 - {\cal{R}}_{3D}}}{{\frac{{2\sqrt 2 C_{s,3D}r_{c,3D}^2}}{{3\sqrt \pi }}\frac{{\sqrt {\bar n} w_0}}{{\tau c}}\left( {1 + \frac{{3\lambda ^2}}{{\pi ^2w_0^2\bar n}}} \right)}}.$$

(21)

Furthermore, we exploit that the reflectance ${\cal{R}}_{3D}$ approaches zero when the transparency factor is sufficiently low, approximately $\bar 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 $\bar n^{3/2} + 3\varpi \bar n^{1/2} - 4\varpi {\mathrm{/}}{\it{\upvarrho }}$, where $\varpi = 3\lambda ^2{\mathrm{/}}\pi ^2w_0^2$ and ${\it{\upvarrho }} = C_{nc,3D}r_{c,3D}^2w_0{\mathrm{/}}\sqrt 2 \pi \tau c$ are constants. Since ϖ approaches zero and $\bar n$ is assumed low, the term $3\varpi \bar n^{1/2}$ is an infinitesimal of higher order than $4\varpi {\mathrm{/}}{\it{\upvarrho }}$ 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:

$$\begin{array}{*{20}{c}} {\bar n^{opt} \approx \left( {\frac{{4\varpi }}{{\it{\upvarrho }}}} \right)^{2/3} = \left( {\frac{{12\sqrt {2\pi } }}{{\pi ^2C_{nc,3D}r_{c,3D}^2}}\frac{{\lambda ^2\tau c}}{{w_0^3}}} \right)^{2/3}}. \end{array}$$

(22)

Equations (18) and (22) can be reformulated in dimensional units to obtain

$$n_{e}^{opt} \approx 0.91\gamma_0 n_c\frac{{\lambda ^2}}{{w_{FWHM}^I{\,}^2}}\left( {\tau _{FWHM}^I\frac{c}{\lambda }} \right)^{2/3},$$

(23)

$$\begin{array}{*{20}{c}} {d_{nc}^{opt} \approx 0.84w_{FWHM}^I\sqrt {\frac{{\gamma _0n_c}}{{n_e^{opt}}}} \approx 0.88\frac{{w_{FWHM}^I{\,}^2{\mathrm{/}}\lambda }}{{\root {3} \of {{\tau _{FWHM}^Ic{\mathrm{/}}\lambda }}}}}, \end{array}$$

(24)

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 $\bar n^{opt}$ in Eq. (15) and neglecting again the factor ϖ, we obtain

$$\frac{{E_{DLT}\left( {n_e^{opt},d_{nc}^{opt}} \right)}}{{E_s^0}} \approx \frac{3}{2}\frac{{C_{nc,3D}}}{{C_{s,3D}}}\left[ {1 - \frac{1}{\pi }\left( {\frac{{\sqrt 3 C_{nc,3D}r_{c,3D}^2}}{2}\frac{\lambda }{{\tau c}}} \right)^{2/3}} \right]\\ \approx 4.58\left[ {1 - 0.92\left( {\frac{\lambda }{{\tau _{FWHM}^Ic}}} \right)^{2/3}} \right].$$

(25)

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 $\gamma _0 = \sqrt {1 + a_0^2{\mathrm{/}}2}$), 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), $\bar E_{DLT}^{opt}$ tends to a constant, equal to about $3C_{nc,3D}{\mathrm{/}}2C_{s,3D} \approx 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₀ < 50). Moreover, as previously explained, our waist evolution equation, describing the self-focusing, is reasonably valid for near-critical density ($0.05 < \bar n < 1$) and for near-critical layer thickness lower than about two-times the self-focusing 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\lambda {\mathrm{/}}\pi \left( {n_c{\mathrm{/}}n_e} \right)_s$⁴. 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³⁸ that P and C polarized pulses produce, when the first layer is sufficiently transparent ($\bar n \lesssim 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 γ₀ factor differs from the one in linear polarization of a factor about $\sqrt 2$.

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 $\bar n \lesssim 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) = \alpha _{ns}E_{nc}(x)$ and $N_{ns}(x) = N_{nc}(x){\mathrm{/}}\alpha _{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 near-critical 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 $\bar d_{nc} = \sqrt {\bar n} d_{nc}{\mathrm{/}}\lambda$, with a maximum at the self-focusing length $\bar d_{nc}^{opt} \approx \sqrt {\bar n} l_f{\mathrm{/}}\lambda \approx w_0{\mathrm{/}}\lambda$. 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 near-critical 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 near-critical 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⁵⁷. The laser pulse had an idealized cos² 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₀ = 2 and a₀ = 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 \times 10^{18}{{\rm{W}}/{\rm{cm}}^2} < I < 1.4 \times 10^{20}{{\rm{W}}/{\rm{cm}}^2}$ which is found in small-medium scale super-intense laser facilities⁵⁸.

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⁶⁺). 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 $\sqrt {{\mathbf{E}}({\boldsymbol{x}})^2 + {\mathbf{B}}({\boldsymbol{x}})^2} {\mathrm{/}}\sqrt {{\mathbf{E}}_0^2 + {\mathbf{B}}_0^2}$ 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_ec². 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⁵⁹.

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⁶⁰, dynamical⁶¹, and hybrid⁶². 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:

$$\begin{array}{*{20}{c}} {{\it{\epsilon }}_p^{max} = T_h\left[ {\varphi ^ \ast - 1 + \frac{{\beta \left( {\varphi ^ \ast } \right)}}{{I\left( {\varphi ^ \ast } \right)e^{\zeta + \varphi ^ \ast }}}} \right] \approx T_h\left[ {{\mathrm{log}}\left( {\frac{{n_{h0}}}{{\tilde n}}} \right) - 1} \right]}, \end{array}$$

(26)

where T_h is the hot electron distribution temperature, φ^* is the normalized potential inside the substrate $\varphi ^ \ast = \phi {\mathrm{/}}T_h$, $\zeta = m_ec^2{\mathrm{/}}T_h$, $\beta \left( {\varphi ^ \ast } \right) = \sqrt {\left( {\varphi ^ \ast + \zeta } \right)^2 - \zeta ^2}$ and $I\left( {\varphi ^ \ast } \right) = {\int}_0^{\beta \left( {\varphi ^ \ast } \right)} {{\rm{e}}^{ - \sqrt {\zeta ^2 + p^2} }} {\rm{d}}p$. The normalized potential is retrieved by solving the implicit equation $\tilde n\frac{{I\left( {\varphi ^ \ast } \right)}}{{\zeta \kappa _1\left( \zeta \right)}}e^{\varphi ^ \ast } = n_{h0}$, with n_h0 the hot electron density, and $\tilde n$ a normalization constant of the hot electron distribution function that is used as a free parameter⁶⁶. 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/γ₀n_c < 1), but classically overcritical (n_e/n_c > 1). Indeed, if we consider a Gaussian pulse amplitude envelope in 2D, $a\left( {t,y} \right) = a_0e^{ - t^2/\tau ^2}e^{ - y^2/w_0^2}$, 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/\gamma \left( {t,y} \right)n_c = 1$. Roughly approximating $\gamma \left( {t,y} \right) \approx a\left( {t,y} \right){\mathrm{/}}\sqrt 2$ we can rewrite it as

$$\begin{array}{*{20}{c}} {\frac{{n_e}}{{\gamma _0n_c}} = \frac{{\gamma \left( {t,y} \right)}}{{\gamma _0}} = {\rm{e}}^{ - \frac{{t^2}}{{\tau ^2}}}{\rm{e}}^{ - \frac{{y^2}}{{w_0^2}}}}. \end{array}$$

(27)

With a change of variables ($t{\mathrm{/}}\tau = \xi$ and $y{\mathrm{/}}w_0 = \chi$) and taking the natural logarithm of the latter equation we obtain

$$\begin{array}{*{20}{c}} {\xi ^2 + \chi ^2 = - \log \bar n}. \end{array}$$

(28)

Fig. 9: Calculation of the reflectance ${\mathbf{{\cal{R}}}}$.

a shows the level plot of the amplitude gaussian in normalized units $\xi = t{\mathrm{/}}\tau$ and $\chi = y{\mathrm{/}}w_0$. The dashed black line marks a general threshold $- \log \bar 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).

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 = \xi ^2 + \chi ^2$, since Eq. (28) represents a circumference:

$$\begin{array}{*{20}{c}} {1 - {\mathbf{{\cal{R}}}}_{2D} = \frac{{{\int}_0^{\sqrt { - \log \left( {\bar n} \right)} } {r{\rm{e}}^{ - 2r^2}} {\rm{d}}r}}{{{\int}_0^{ + \infty } {r{\rm{e}}^{ - 2r^2}} {\rm{d}}r}} = 1 - \bar n^2}. \end{array}$$

(29)

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 $\bar n$, when a₀ 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:

$$\begin{array}{*{20}{c}} {1 - {\mathbf{{\cal{R}}}}_{3D} = \frac{{{\int}_0^{\sqrt { - \log \left( {\bar n} \right)} } {r^2{\rm{e}}^{ - 2r^2}} {\rm{d}}r}}{{{\int}_0^{ + \infty } {r^2{\rm{e}}^{ - 2r^2}} {\rm{d}}r}} = {\mathrm{erf}}\left( {\sqrt { - 2\log \bar n} } \right) - \frac{4}{{\sqrt {2\pi } }}\bar n^2\sqrt { - \log \bar n} }. \end{array}$$

(30)

Approximations validity

We discuss the range of validity of the approximations underlying Eqs. (22)–(25). In order of appearance we assumed the following.

$\bar n > {\uplambda}^2{\mathrm{/}}2w_0^2$: substituting this inequality into Eq. (22), we get the condition $\tau > \pi ^2r_{c,3D}^2C_{nc,3D}\lambda {\mathrm{/}}12\sqrt {2\pi } c \approx 1.6\lambda {\mathrm{/}}c$. This corresponds to a FWHM temporal duration longer than 5 fs which is generally verified for nearly all high intensity laser systems.

$\bar n < 1{\mathrm{/}}4$: from this condition we find $w_0{\mathrm{/}}\lambda > \root {3} \of {{24.4\,\tau c{\mathrm{/}}C_{nc,3D}r_c^2\lambda }}$. Fixing the FWHM temporal duration to 28.5 fs, the latter inequality reads as $w_0 > 3.5\lambda = 2.8\,{\rm{\mu m}}$, which is often the case for high intensity laser systems. Note that the inequality weakly depends on the temporal duration of the pulse, because of the cubic root.

$3\varpi \bar n^{1/2} \ll 4\varpi {\mathrm{/}}{\it{\upvarrho }}$: by substituting Eq. (22) into the inequality, we obtain that $\tau \gg 3^{3/2}C_{nc,3D}r_c^2\lambda {\mathrm{/}}2^{5/2}\pi ^{3/2}c$; we find $\tau \gg 0.8\,\lambda {\mathrm{/}}c$ which is implied by the first condition.