Integrating sheath and radiation-based a

Due to the ultra-thin nature of the graphene target, simulating the interaction of a prepulse with the target is of interest. However, due to limitations in conventional hydrodynamic and particle-in-cell (PIC) methods, the prepulse simulations are carried out in the current study. The results from these simulations are used as the initial conditions for the subsequent main pulse simulation using the PIC method. The schematic of this method is shown in Fig. 1.

Three dimensional molecular dynamics simulation

A graphene target is modeled and simulated within the 3D molecular dynamics (MD) environment called LAMMPS¹⁵. A laser prepulse profile similar to that from the experiments^10,16 is applied to the MD graphene target. An elementary description of this MD methodology for simulating the prepulse-target interaction is provided as supplementary information 1.

The MD simulations were carried out for both the 4-layer and 8-layer graphene targets. The prepulse leads to the compression of the graphene target at the spot area, resulting in a slight increase in target density. The average density of carbon at the spot area, after prepulse interaction, is $n_\mathrm{C^{6+}}$=1.2 $\times$ 10$^{29}$$\hbox {m}^{-3}=120n_c$, which is nearly identical to the target’s initial density. Here, $n_c$ is the critical density obtained for $\lambda =1054$ nm. A similar compression behavior was also observed for the 4-layer graphene target. The averaged carbon density from the spot area in MD is considered as the initial density of the fully ionized carbon species ($n_\mathrm{C^{6+}}$) in the 2D PIC simulation. The corresponding atomic kinetic energy is used as the initial temperature for the electrons.

Two dimensional particle-in-cell simulation

Since impurities such as protons ($\hbox {H}^+$) are also expected to be included in the PIC simulation, a target was created in the EPOCH PIC environment¹⁷ with an electron density of $n_e=1.65\times 6n_\mathrm{C^{6+}}$ and a proton density $n_\mathrm{H^{+}}=0.05n_e$. The factor 1.65 arises from the experimental measurement of the contaminant contribution in suspended graphene¹⁰. The effect of target deformation from MD simulations is neglected, as separate simulations confirmed that the deformation profile had a negligible contribution to the final ion energies. The target is made charge neutral through the introduction of proton contaminants, considered as 1 nm thick layers sandwiching the carbon core. The electrons have a temperature of 10 keV, equivalent to the average kinetic energy of atoms from the MD simulations, and the ion species is assumed to be cold. Currently, this setup causes the entire target to undergo pre-expansion uniformly before the main pulse arrives. While this uniform pre-expansion can lead to more realistic target parameters within the laser spot area during the main pulse interaction, it’s important to note that electron pre-expansion should occur only within the spot area. The widespread pre-expansion currently is an artifact resulting from transitioning from MD simulations to PIC simulations. In future work, we will incorporate the kinetic energy distribution profiles from MD simulations to ensure that pre-expansion in the PIC simulations is realistically localized to the laser spot area.

The main pulse is a LP pulse with $\lambda =1054$ nm, pulse duration $\tau =770$ fs, and spot diameter (FWHM) $d=10$$\mu$m with a peak intensity of $I=3.2\times 10^{20}$ W/$\hbox {cm}^2$. These pulse parameters are similar to the VULCAN laser pulse parameters used in experiments¹⁴. The pulse is Gaussian both spatially and temporally. The temporal Gaussian profile for the laser with a time shift of w, where w is the laser FWHM, is used. The laser enters the domain from the left boundary, while the other boundaries are configured as simple outflow boundaries. A 2D simulation box of dimension 130 $\mu$$m \times 30$$\mu$m comprising of 390,$000 \times 3$,000 cells is defined, and it is ensured that there are at least 100 (700) ions (electrons) per cell. The proton contaminants were modelled as 1 nm layers on either side of the $\hbox {C}^{6+}$ core. The PIC simulation was carried out with the above laser and simulation setup for 4-layer (4 nm), 8-layer (6 nm), 12-layer (8 nm), 16-layer (10 nm), and 24-layer (14 nm) graphene targets. The thickness of the carbon core in these targets are defined based on the MD simulation results. The laser contacts the target at $t \approx 133\,\textrm{fs}$.

TNSA-RPA hybrid acceleration

The maximum $\hbox {C}^{6+}$ and $\hbox {H}^+$ ion energy with respect to time for all the cases is shown in Figs. 2(a) and (b), respectively. The time evolution of the energy is similar for the 8, 12, and 16-layer cases, with the maximum $\hbox {C}^{6+}$ ($\hbox {H}^{+}$) energies reaching over > 1 GeV (> 100 MeV) around 1 ps. However, the energy evolution profile for the 4-layer and the 24-layer cases differs significantly from the other cases. Additionally, one-dimensional (1D) simulations under similar conditions were performed to confirm the energy evolution trends for the various target cases. The energy evolution plot from the 1D simulations is provided in the supplementary, and the results generally align with the results from the 2D simulations. The 1D cases were also extended to a longer duration, where a change in the slope of energy evolution indicating ion energy saturation could also be observed.

4-layer graphene

The initial electron temperature causes electron pre-expansion on either side of the target. As the pulse interacts with the target, the electrons in the spot area expand faster than in the surrounding regions, which remain electron-rich by comparison. This pre-expansion and subsequent laser-driven expansion are highlighted in Fig. 3(d), and neither the pre-expanded electron population nor the electron-rich regions outside the spot area contribute significantly to the final ion energies. For the 4-layer case, the target thickness is below the optimal thickness to achieve RPA for the given laser intensity. The condition $a \le \zeta$, required for obtaining maximum radiation pressure, is violated. Here, a is the normalized laser amplitude, and $\zeta =\pi (n_e/n_c)(l/\lambda )$ characterizes the optical property of the target^9,18, $n_e$ is the electron density, $n_c$ is the critical density, l is the target thickness, and $\lambda$ is the laser wavelength. Consequently, the electrons are blown out from the target, interrupting the acceleration process. For the current target and laser parameters, $a \approx 15$ and $\zeta \approx 13.4$. The electron blowout could be verified from the simulations by examining the time evolution of the electron density profiles shown in Figs. 3(a), (b), and (c). The RPA electric field peak R (black solid line in Figs. 3(d), (e), and (f)) causes an electron blowout, and there is no hot electron cloud formation behind the target as seen in Fig. 3(b).

An electron depletion region is therefore formed in the target rear as seen in Fig. 3(c). This reduces the intensity of the TNSA sheath acceleration electric field peak T that is setup behind the target as seen in Figs. 3(d), (e), and (f). The peaks R and T collectively form a double peaked electric field (DPEF) structure, initially proposed for the linearly polarized pulse by Zhuo etal.¹². However, in this context, the radiation pressure dominates over the electrostatic pressure resulting in sub-optimal acceleration of the $\hbox {C}^{6+}$ ions in the DPEF. This contrasts with the 24-layer case, where the radiation pressure is unable to penetrate through the target before the electrostatic pressure becomes dominant. This leads to hotter electrons, and a more pronounced TNSA-driven acceleration.

As for the protons, it is not accelerated by the DPEF field but rather by the shielded Coulomb repulsion force, as discussed by Liu etal.¹⁹ for a multi-species foil. This is evident from Figs. 3(d), (e), and (f), where the proton (green) is accelerated as monoenergetic layer over an extended duration, without overlapping the carbon (red) layer. The electric field peak P further signifies the Coulomb repulsive force exerted on the proton layer by the carbon ions.

8-layer graphene

Similar to the 4-layer case, the electron pre-expansion and the laser-driven expansion regions are highlighted in Fig. 4(d). The acceleration mechanisms in the 8, 12, and 16-layer graphene cases share many similarities. Therefore, this discussion focuses solely on the 8-layer case. This is confirmed by the electron density distribution shown in Figs. 4(a), (b), and (c). In the 8-layer target, the condition $a \le \zeta$ is met, resulting in a ponderomotive force-driven electron compression layer, few nanometer thick, at the target’s rear as evident in Figs. 4(a), (b), and (c). The compression from the ponderomotive force in combination with the laser heating blows away the electrons from the target front leading to a high density electron region at the target rear as seen in Fig. 4. A brief discussion regarding the evolution of the particle densities and the electric field under these conditions is provided in the supplementary.

The proton layer continues to be accelerated by the shielded Coulomb field, but the acceleration is more effective in this case. This enhancement can be attributed to the increase in the carbon ion concentration behind the proton layer, as the Coulomb repulsive force increases in proportion to the total charge of the carbon layer¹⁹. This leads to an intense peak P at earlier times, as seen in Fig. 4(d), which eventually subsides as shown in Figs. 4(e) and (f), due to the expanding population of hot electrons that effectively shields the carbon ions. The $\hbox {C}^{6+}$ ions also experience more efficient acceleration by the DPEF structure, particularly between 150 fs to 250 fs. Beyond this period, the DPEF structure diminishes over time, as the intensity of the peak T reduces due to target expansion. The diminishing of the clear distinction between the peaks T and P beyond 200 fs can be observed from the electric field evolution shown in the supplementary. The emergence of this DPEF structure, particularly at 150 fs, aligns well with the observed sudden increase in particle energies, as seen in Figs. 2(a) and (b).

The target remains opaque during the nearly 100 fs long DPEF acceleration phase. Throughout this phase, the condition $n_e > \gamma n_c$²⁰ is consistently observed as detailed in the supplementary. Here, $\gamma n_c$ represents the relativistically corrected density. This condition precludes the possibility of enhanced acceleration due to the relativistic induced transparency (RIT) mechanism, as discussed in a past experimental study by Higginson etal.¹⁴, where DPEF coexisted with RIT, leading to 100 MeV protons. Past studies^21,22,23 have also found further evidence that RIT could be avoided for ultra-thin targets by coating a lighter species target with a heavier species. This allows for a stable LS-RPA acceleration of the lighter species by the excess electrons from the heavier species, similar to the proton acceleration observed in this study. The LS-RPA ion energy scaling thus obtained, accounting for the laser pulse and target properties, is independent of laser polarization since the cycle-averaged laser pressure is the same for both CP and LP pulses^23,24,25. The appropriate target thickness required to avert RIT and achieve bulk RPA throughout the pulse duration could also be inferred from these scaling relationships.

Additionally, there was no onset of the Bunemann instability, that couples the electron momentum to the ion momentum. This absence further negates the the possibility of breakout afterburner (BOA) acceleration^{26,27,28,29,30}, a mechanism that can accelerate the carbon ions to higher energies. A past experimental study³¹ has attributed the occurrence of BOA to the generation of 1 GeV carbon ions.

According to the hybrid acceleration model proposed by Qiao et al.¹³ RPA dominancy for graphene is achieved for target thickness within the range 4.54 nm < l < 54 nm, for the laser and target parameters used in this study. The 8, 12, and 16-layer graphene falls within this range, resulting in the generation of extremely energetic ions. Conversely, the 4-layer graphene, being below the lower limit, generates ions with less energy. While the simulation results generally align with this model, the 24-layer graphene, despite being within the model’s proposed range, shows different energy scaling and potentially lower energies.

Hybrid acceleration framework

To better understand the acceleration dynamics of energetic carbon ions in the DPEF, the hybrid acceleration framework has been established. This framework primarily considers only the TNSA and RPA mechanisms. The shielded Coulomb repulsion acceleration can be neglected, as it primarily affects protons. Previous research by Liu et al.¹⁹ demonstrated that, due to the complete separation of carbon and proton layers and the protons being a small minority (both points applicable to the current study), they have almost no effect on the motion of carbon ions. The current study focuses on the 8, 12, and 16 layer cases, as these targets achieve higher carbon ion energies likely through the HA process. The maximum velocity achieved by a particle undergoing HA at a given time $\tau$ is defined by Qiao et al.³² as the Lorentz transformed velocity summation of TNSA and RPA. This velocity is termed as the hybrid velocity $v_h$. In our study, we adopt a similar definition for $v_h$, but with the introduction of the scaling coefficients $\alpha$ and $\beta$ as follows.

$$\begin{aligned} v_{h}=\frac{\alpha v_t+\beta v_r}{1+\alpha \beta v_tv_r/c^2}, \end{aligned}$$

(1)

where $v_t$ is the ion TNSA velocity, $v_r$ is the ion RPA velocity, and $\alpha$ and $\beta$ are the TNSA and RPA specific scaling constants, respectively. The particle energy is then obtained from $v_h$ as $E_h=(\gamma _h-1 )m_ic^2$, where $\gamma _h=1/\sqrt{1-v_h^2/c^2}$ is the Lorentz factor, $m_i$ is the ion mass, and c is the speed of light. $v_t$ can be estimated based on the model proposed by Schreiber et al.³³ given by

$$\begin{aligned} \frac{\tau }{\tau _0}=X\bigg (1+\frac{1}{2(1-X^2)}\bigg )+\frac{1}{4}\textrm{ln}\frac{1+X}{1-X}, \end{aligned}$$

(2)

where $\tau$ is the laser pulse duration, $\tau _0=(r_L+l\,\mathrm{tan\theta })/v_{\infty }$, and $X=(E_t/E_{t,\infty })^{1/2}$. $E_t$ is the ion TNSA energy, $E_{t, \infty }=Zk_BT_e(r_L+l\,\mathrm{tan\theta })/\lambda _D$ is the theoretical maximum energy of an ion with charge Z, and therefore the maximum achievable ion velocity is $v_{\infty }=(2E_{t,\infty }/m_{i})^{1/2}$. Here, $r_L$ is the laser spot radius, $\theta$ is the half-angle of the electrons traveling through the target, $k_B$ is the Boltzmann constant and $\lambda _D$ is the Debye length. The hot electron temperature $T_e$ is given by $T_e=(\gamma n_c/n_e)^{1/2}(\gamma -1)m_ec^2$,³⁴ where $\gamma =(1+a^2)^{1/2}$, and $v_t$ can be obtained from $E_t$ through $v_{t}=(2E_{t}/m_{i})^{1/2}$.

The RPA velocity can be obtained from the equation describing the LS motion of the target^7,9:

$$\begin{aligned} \frac{du}{dt}=\frac{a^2}{n_0MD}\frac{\sqrt{1+u^2}-u}{\sqrt{1+u^2}+u}, \end{aligned}$$

(3)

where u=P/mc is the normalized momentum with $P=mv_r$, t is the time normalized by the laser period, $n_0=n_e/n_c$ is the normalized electron density, $M=m/m_e$, $m_e$ is the electron mass, and $D=l_\textrm{eff}/\lambda$ is the normalized effective target thickness. The radiation pressure accelerates only a thin layer of thickness $l_\textrm{eff}=l-d$ on the target rear, where $d=a\pi \lambda {n_c}/{n_e}$ is the electron depletion depth of the target due to the ponderomotive force⁹.

Based on the curve fit of Eq. (1) to the temporal energy evolution of 8L, 12L, and 16L in Fig. 2(a), for unique values of $\alpha$ and $\beta$, we propose the definition of these scaling coefficients in terms of the optical parameter ($\xi =l\omega _{\mathrm pe}^2/2c\omega$)^9,18, where $\omega _{\mathrm pe}$ is the plasma frequency and $\omega$ is the laser frequency.

$$\begin{aligned} \begin{aligned} \alpha&= 3.8\frac{c}{aC_h}\frac{n_c}{n_e}\xi ,\\ \beta&= \frac{\sqrt{3}a}{\xi }\bigg (\frac{l}{d}-1\bigg ), \end{aligned} \end{aligned}$$

(4)

where $C_h=(Zk_BT_e/m_i)^{1/2}$ is the hot-electron sound speed. For a given set of laser and target parameters, by combining Eqs. (1) – (4), specific values of $v_h$ at each time could be obtained. $\xi$ determines the opacity of the target to the laser and, for $\alpha$, it is modified by $C_h$ and $n_e$ indicating the role of laser heating on target transparency. Meanwhile in $\beta$, $\xi$ is modified by $l_\textrm{eff}$ indicating the radiation pressure’s role on target transparency. In combination, Eq. (4) provides the fractional roles of TNSA and RPA on ion acceleration through the coefficients $\alpha$ and $\beta$. Note that while Eq. (2) predicts ion energy saturation, it does not account for the effects of plasma cooling due to expansion. In contrast, Eq. (3) does not incorporate energy saturation. Since Eq. (1) combines these equations, it similarly lacks a saturation limit. The coefficients $\alpha$ and $\beta$ solely modify the slope of the temporal energy evolution, making the energy curves in Figs. 2(a) and (b) sufficient to determine if a target is undergoing hybrid acceleration (HA) before reaching saturation. Nonetheless, the ion energies in the 1D simulations provided in the supplementary materials are nearing saturation, and applying this theory demonstrates good agreement with the 2D results.

When $\alpha =\beta =1$, Eq. (1) presents the TNSA and RPA velocities under their respective ideal scenarios. The Schreiber definition of $v_t$ and the equation of motion of LS for $v_r$ treats TNSA and RPA as mutually independent acceleration mechanisms. However, in reality, these mechanisms can potentially influence each other. This paper aims to explore this coupling behavior through the introduction of non-dimensional scaling coefficients $\alpha$ and $\beta$. The primary assumption is that TNSA and RPA are considered to be the predominant acceleration mechanisms, and due to total energy conservation, an increase in particle energy due to RPA implies a decrease in contribution from TNSA, and vice versa.

The scaling coefficients are useful in tuning this TNSA-RPA energy balance based on the target/laser conditions. The definitions of $\alpha$ and $\beta$ contain parameters that influence the energy gained by particles through the respective acceleration mechanisms. The coupling between the acceleration mechanisms (and between $\alpha$ and $\beta$) is achieved through the introduction of $\chi$, which controls the target’s transparency. In the case of $\alpha =1$, Eq. (4) can be rewritten as $C_h \approx 2l\omega /a$, which represents efficient energy transfer from the laser to the target electrons. The condition $\alpha =1$ therefore indicates perfect heating of the target by the laser. Similarly, for $\beta =1$, Eq. (4) can be rewritten as $\xi =\sqrt{3}a l_\textrm{eff}/d$, which balances laser penetration against the target’s resistance to the penetration. Thus, the condition $\beta =1$ indicates perfect reflectance of the laser, efficiently transferring momentum to the target and accelerating it effectively as a light-sail.

The coefficients $\alpha$ and $\beta$ are related through $\xi$, and we find that $\alpha \propto \beta ^{-1}$. Since $\beta \ge 0$, the target fraction condition $l/d \ge 1$ should be satisfied, and by setting $\beta = 0$, we obtain $l_{r,min}=d\sim 4.5$ nm as the lower limit of target thickness required to prevent electron blowout, in agreement with the previous studies^9,13. Here, $l_{r,min}$ is the minimum target thickness where RPA could be achieved. For $\alpha =0$, $l_{t,min}=0$, and an upper limit on the values of $\alpha$ and $\beta$ is imposed by setting $\alpha =\beta =1$, which gives a maximum target thickness of $l_{r,max}=10.75$ nm for $\beta =1$ and $l_{t,max}=59$ nm for $\alpha =1$. Here, $l_{t,min}$ is the minimum target thickness where TNSA is achieved. The RPA dominant region for graphene thus obtained from the $\beta$ limits is 4.5 nm $<l<$ 10.75 nm, which is much narrower than the range proposed by the previous model¹³.

Note that the numerical values of 3.8 and $\sqrt{3}$ in Eq. (4) are derived from modifications to the scaling equations based on curve fitting of particle energy scaling observed in the simulations. It is important to note that two-dimensional (2D) simulations can overestimate the maximum ion energies almost by a factor of 2 compared to more realistic three-dimensional (3D) simulations. This overestimation has been noted in several previous studies^35,36,37 and may explain the adjustments made to the scaling equations. While the maximum ion energies can potentially become lower in the actual 3D phenomenon, the slope of the energy evolution over time could remain similar to that observed in the 2D case. Future work could benefit from incorporating 3D simulations to provide a more accurate representation of the energy dynamics and to further validate Eq. (4).

The RPA dominant region, within the hybrid framework, for the graphene target under different laser intensities for the current simulation target parameters is shown in Fig. 5. The RPA dominant region ($0\le \beta \le 1$) co-exists with two other regions, where for $\beta <0$ there is electron blowout that disrupts the RPA and for the $\beta >1$ region, the TNSA starts dominating over RPA. In this region, the ions undergo an acceleration mechanism similar to the relativistic induced transparency,¹⁴ where volumetric heating causes the ion energies to increase (energies less than achievable by RPA dominancy) by enhancing TNSA. For $\alpha >1$, the acceleration mechanism falls outside the hybrid framework described by the scaling coefficients in Eq. (1). The 24-layer case falls in this region. The region $\beta < 0$ will be termed as the blowout region, while the region $0\le \beta \le 1$ will be termed as the RPA region, and the region $\beta >1$ and $\alpha <1$ will be termed as the transparency region termed after the relativistic induced transparency mechanism. All the three regions experience hybrid acceleration to different extents, and our definition of $\beta$ is applicable only to the region where RPA is dominant.

Table 1 Coefficients $\alpha$ and $\beta$ values obtained from theory (Eq. (4)) are compared with the simulation values obtained by fitting Eq. (1) to the $\hbox {C}^{6+}$ energy evolution shown in Fig. 2(a). The error is much less and the $\hbox {C}^{6+}$ energy scales to GeV energies for the targets (8, 12, and 16 layers) undergoing RPA dominant HA.

Table 1 shows the coefficients predicted by Eq. (4) along with $\alpha$ and $\beta$ values obtained by fitting Eq. (1) to the $\hbox {C}^{6+}$ energy evolution from the simulation. From the maximum percent error that arises by comparing the scaling coefficients from theory and simulations, it is clear that the theoretical $\hbox {C}^{6+}$ hybrid energy evolution is in good agreement with the simulations for the targets undergoing RPA dominant acceleration in the HA framework. The error is large for targets outside the RPA region but can provide qualitative information about the type of acceleration. The energy of $\hbox {C}^{6+}$ ions for targets which are exclusive to the RPA dominant region achieves GeV energies (around twice the other cases). To verify the upper limit of RPA dominancy (10.75 nm), we choose a 24-layer graphene target whose thickness of 14 nm ($\alpha =0.583,\, \beta >1$ from Eq. (4)) lies in the transparency region. The corresponding $\hbox {C}^{6+}$ energy evolution profile is shown in Fig. 2(a), and it could be seen that the maximum energy obtained is much less than the cases which lie in the RPA region. The hot electrons generated by the oscillating laser field penetrates through the target before the hole boring from RPA can, and establishes a strong sheath field at the target rear, thus making the acceleration TNSA dominant.

Table 2 The theoretical $\alpha$ and $\beta$ obtained from Eq. (4) for the past simulation and experimental conditions. Cases are in blowout region when $\beta <0$, in RPA region when $0\le \beta \le 1$, and in transparency region when $\alpha \le 1$ and $\beta >1$. The particle energies are maximum when they are in the RPA region.

For further validation of the coefficients, a comprehensive literature review has been conducted, and the $\alpha$ and $\beta$ values for the past experiments and simulations relating to HA have been obtained as shown in Table 2. It is observed that, the most energetic particle generation cases lie in (800 nm¹², 150 nm¹⁴) or close to (80 nm¹³, 75 nm, 90 nm, 300 nm, 500 nm¹⁴) the RPA region. The energy of the ions from these simulations and experiments increases when they are closer to the RPA region, and decreases for an increase in distance (200 nm¹³, 1000 nm¹⁴) from the RPA region in either direction. Finally, the ion energy is the lowest for cases above the transparency region (1500 nm¹⁴).

Finally, the RPA region for different materials is compared by maintaining a constant laser intensity while varying the target’s mass density. The coefficient $\beta$ is inversely related to $\chi$, which includes the electron number density $n_e$ via the plasma frequency $\omega _e$. The values of $n_e$ for different materials can then be derived from their mass densities. The acceleration regions obtained for graphene (C), aluminum (Al), prepulse compressed contaminated graphene ($\hbox {C}^{*}$) (used in the PIC simulations in this study) , copper (Cu), and gold (Au) are shown in Fig. 6. It is found that for an increase in the mass density, the RPA region becomes narrower, thus necessitating the use of thinner targets. However, producing such thin targets is often experimentally unfeasible. This observation underscores the advantage of using low mass density targets for RPA applications, as they not only facilitate the establishment of RPA dominance but are also relatively easier to fabricate.

Integrating sheath and radiation-based acceleration using scaling coefficients for tailoring radiation dominant hybrid accelerat

Three dimensional molecular dynamics simulation

Two dimensional particle-in-cell simulation

TNSA-RPA hybrid acceleration

4-layer graphene

8-layer graphene

Hybrid acceleration framework