Deep learning approaches for modeling laser-driven proton beams via phase-stable acceleration

Laser-driven ion acceleration is a topic of significant importance in the field of laser–plasma interaction. The generation of high-energy ion beams holds great potential for applications in cancer therapy as well as for fundamental science studies.^1–3 Various acceleration schemes have been proposed. One widely studied scheme is target-normal sheath acceleration (TNSA), which involves irradiating a solid target with an intense laser pulse. The laser–plasma interaction creates a plasma sheath at the target rear, causing ions to be accelerated away from the target.^4–7 Another notable scheme is radiation pressure acceleration (RPA), where the radiation pressure of the laser pulse drives the acceleration of ions. The interaction between the laser and the plasma generates a strong electrostatic field, which aids in the ion acceleration process.^8–12 Breakout afterburner (BOA) is another acceleration scheme that involves a two-stage process using intense laser pulses and skin-depth scale targets.^13,14 During the initial phase of laser absorption at the front surface of the target, rapid evolution occurs, leading to the volumetric heating of electrons throughout the target volume beneath the laser focal point. As a consequence, the nanometer-scale target quickly transitions to a state of relativistic transparency, allowing the laser to penetrate through the dense target plasma. During the subsequent phase, the co-motion of the accelerating electric field and ions throughout the acceleration period is instrumental. This mechanism results in a substantial enhancement in beam energy and conversion efficiency when compared to traditional TNSA. Magnetized-vortex-acceleration (MVA) is a scheme that exploits the interaction between a magnetized plasma and a laser pulse. The presence of a magnetic field enhances the ion acceleration by enabling self-focusing of the laser pulse and reducing the electron energy loss.^15,16 Collisionless shock acceleration (CSA) is a scheme where a shock wave is generated in a plasma through the interaction of an intense laser pulse. This shock wave can efficiently accelerate ions to high energies without significant heating.^17,18

The existing acceleration schemes are primarily based on one-dimensional (1D) theory and demonstrate good agreement with particle-in-cell simulations (PICs), but multi-dimensional PIC simulations reveal that the complex physical processes involved¹⁹ and transverse effects introduced instabilities significantly reduce the efficiency of ion acceleration.^20–23 Moreover, numerous simulations indicate that the ion energy and spectrum are highly sensitive to the laser and target parameters, such as beam size, energy, focal position, wavelength, target thickness, density, and structure. Despite achieving the highest recorded proton energy of 94 MeV in 2018,²⁴ which nearly meets the requirements for radio-biological in vivo studies, further enhancements in proton energy and numbers are necessary for general tumor treatment. In 2022, an experiment of tumor irradiation in mice with a laser-accelerated proton beam was conducted,²⁵ and thus, clinical experiments by laser-driven proton beam will be expected in the coming future.²⁶ 

The progress of laser-driven ion acceleration has been limited and Badziak discussed the challenges and prospects,²⁷ e.g., understanding and controlling ion beam properties, mitigating instabilities, optimizing laser and target parameters, and developing diagnostic techniques to characterize the ion beams accurately. Therefore, finding an optimized parameter set and understanding the underlying physics are critical. However, the consumed-cost associated with conducting experiments or performing multi-dimensional PIC simulations presents a challenge in parameter optimization.

Recently, deep learning (DL) has emerged as a powerful tool for scientists to address the challenges of the optimization problem because²⁸ DL possesses the potential to explore the vast parameter space associated with ion acceleration. By leveraging large datasets generated from simulations and experiments, DL algorithms can capture intricate patterns, correlations, and nonlinear relationships that may not be easily discernible through traditional analytical approaches. This data-driven approach allows scientists to extract valuable insight and develop predictive models for ion acceleration.²⁹ For example, in recent work by Djordjevic et al., in 2021, an ensemble of over 1000 1D-particle-in-cell (PIC) simulations was used to develop a neural network (NN) surrogate model to study the dependence of ion energy on the pre-plasma gradient length scale in target-normal sheath acceleration (TNSA) and conclude that the pre-plasma is critical issue in the interaction between the laser and the target.³⁰ This paper suggests carefully measuring the laser pre-pulse and amplified spontaneous emission during experiments for modeling the pre-plasma. Hence, varying the pre-plasma profile would be their future study. Therefore, the approach allowed for a comprehensive exploration of the parameter space and provided valuable insight into the underlying physics of ion acceleration.

In addition to the NN surrogate model, Dolier et al. demonstrated the optimization of laser-driven proton acceleration PIC simulations controlled by a Bayesian algorithm.³¹ Moreover, Goodman et al. also indicated the utility of applying a machine learning-based optimization approach providing new insight into the physics of radiation generation in laser–foil interactions, which will inform the design of experiments in the quantum electrodynamics (QED)-plasma regime.³² An essential development in the integration of physics and machine learning is physics-informed neural networks (PINNs).³³ A PINN is a specialized class of neural networks that leverage the principles of physics, frequently expressed in the form of ordinary differential equations or partial differential equations, to act as a regularization mechanism. This ensures that the model's solutions align harmoniously with the fundamental laws of physics. Furthermore, Scheinker and Pokharel recently presented a physics-constrained neural network (PCNN) approach to solving Maxwell's equations for the electromagnetic fields of intense relativistic charged particle beams.³⁴ In 2023, Döpp et al. published a review article that presents an overview of relevant machine learning methods with a focus on applicability to laser-plasma physics and its important sub-fields of laser-plasma acceleration and inertial confinement fusion.³⁵ Hence, data-driven science and machine-learning methods as a new approach to research have been a trend.

Our NN surrogate model can predict the maximum proton energy by the given parameter set $( a L , D / λ L , n 0 / n c )$ to reveal the validity of the PSA theory. To verify the predictions, we performed 1D-PIC to compare with the prediction to demonstrate the successful establishment of a PSA scenario, where the radiation pressure and the electrostatic field achieve a delicate equilibrium, resulting in highly efficient and controlled ion acceleration. Furthermore, our NN surrogate model also indicates a discrepancy between the PSA theory and PIC simulations at larger target thickness region (⁠ $D / λ L > 1$⁠). By exploring scenarios with larger target thicknesses, we observed the disruption of the electron layer due to an imbalance between the radiation pressure and the electrostatic force, ultimately leading to the termination of the PSA process. We performed 1D-PIC to investigate the physical reason and found that the collapse of the electron layer occurs when the plasma skin depth surpasses the thickness of the pulse-compressed electron layer during the interaction process, causing the termination of the PSA process.

In our study, we constructed a dataset consisting of scalar outputs extracted from 355 one-dimensional (1D) particle-in-cell (PIC) simulations performed using the EPOCH PIC code.³⁷ The simulation setup is illustrated in Fig. 1. We considered a circularly polarized laser pulse with a wavelength of $λ L = 800 nm$⁠. The laser pulse was modeled as a square pulse with a duration of $100 T L$ and a Gaussian profile with a width of $0.3 T L$ on both the rising and falling edges, where T_L represents the laser period. The target was a pure hydrogen plasma with a step density profile, and the front edge was located at x = 0. The target thickness (D) and density n₀ were normalized to the laser wavelength λ_L and the critical density n_c, respectively. We employed a grid number with N_X = 50 000 cells and a grid size of $Δ X = 2.4 nm$⁠. Each cell contained 500 macro-particles. Here, we aimed to use DL to verify the PSA formula so the pulse duration $100 T L$ is purposely fixed, which is the same as Yan's PSA paper.³⁶ The simulation time for each run was $500 fs$ because it is long enough to distinguish the successes and the failures of the PSA scheme.

To construct the dataset, we varied the laser amplitudes a₀ as 5, 10, 15, and 20. The target thickness $D / λ L$ ranged from 0.01 to 1.5, while the target density $n 0 / n c$ ranged from 0.5 to 15. At $t = 500 fs$⁠, we recorded the maximum kinetic energy E_k of protons in each simulation run. The 1D-PIC simulation results for the dataset of deep learning and the phase-stable acceleration (PSA) theory (illustrated by black-dashed lines) are shown in Fig. 2. The higher proton energies were obtained near the region of the black-dashed line showing that the proton was efficiently accelerated when the PSA criterion was satisfied.³⁶ Moreover, the maximum kinetic energy of protons E_k increases with laser strength a_L. In Subsection II B, we aim to develop a neural network (NN) surrogate model to verify the PSA theory and predict the maximum proton energy E_K in laser-driven ion acceleration instead of running the consumed-cost PIC.

To construct the NN surrogate model, we employ the PyTorch framework, a powerful tool for deep learning that has gained widespread recognition in the scientific community.³⁸ The input dataset, as shown in Fig. 1, used for training and validation consists of input parameters (a_L, $D / λ L , n 0 / n c$⁠) paired with their corresponding maximum proton energy E_k. By leveraging this dataset, we can guide the NN surrogate model to learn the intricate relationships between the input parameters and the target variable. Our NN architecture follows a four-layer structure, as depicted in Fig. 3. PyTorch's flexible and intuitive interface enables us to easily define the layers and activation functions within the network, allowing for seamless customization and fine-tuning.

The initial layer serves as the input layer, where the dataset's three features are fed into the network. Leveraging the SELU activation function (Scaled Exponential Linear Unit), this layer outputs nine features. Subsequently, the second layer applies the LeakyReLU activation function (leaky rectified linear unit) to the nine input features, resulting in the transformation into 27 output features. In the third layer, the 27 input features are further expanded to 81 output features using the ELU activation function (exponential rectified linear unit). Finally, the output layer takes the 81 input features and generates a single target value E_k employing the ReLU activation function (rectified linear unit). ReLU, widely adopted in regression tasks, effectively handles non-linearity by passing positive values while disregarding negative ones which guarantees none of the negative proton energy output.

To further validate the performance of our NN surrogate model, we conducted a comparative analysis by contrasting its predicted values with those derived from an ensemble of 355 1D-PIC simulations, depicted as the blue data points in Fig. 5. This comparative examination serves as a pivotal evaluation of the predictive capabilities inherent in our NN surrogate model. In Fig. 5, the reference line with a slope of 1, representing the ideal predictions by the NN surrogate model, indicates a scenario where predicted values align perfectly with the simulation data. The blue data points exhibit exceptional concordance between the PIC-derived values and the NN-predicted values, resulting in a corresponding R-squared value of $R 2 = 0.997$⁠.

To further scrutinize the reliability of our NN model, we conducted additional PIC simulations using datasets both within and beyond the confines of the training domain, defined as $5 < a L < 20 , 0.5 < n 0 / n c < 15$⁠, and $0.01 < D / λ L < 1.5$⁠. These are represented by the red (within) and green (beyond) data points depicted in Fig. 5. The computed R² values for the red and green datasets are 0.83 and –1.56, respectively. In the case of the red data points, the majority of predictions exhibit a favorable alignment with the PIC simulations, albeit with some instances of underestimation or overestimation. The scattering of red points corresponds to the predictive accuracy of our model. This accuracy, defined as the percentage of sources with errors less than 10%, reaches 0.6 in the final epoch, as shown in Fig. 4(b). Out of a total of 270 red points, 162 exhibit errors below 10%. Consequently, the calculated predictive accuracy for the red points is 162/270 = 0.6, aligning with the overall model accuracy.

Conversely, for the green data points, most predictions tend to be overestimated, and the negative R² value signifies a challenging extrapolation task for the NN surrogate model, where the data points fall outside the training domain, defined as $20 < a L < 40 , 15 < n 0 / n c < 30$⁠, and $1.5 < D / λ L < 3$⁠, which may cause the green points are always beyond the ideal model prediction (black line). Notwithstanding the rudimentary NN surrogate model presenting the inherent limitations of extrapolation predictions in this paper, numerous scientists are still dedicated to the refinement of extrapolation techniques and the reduction of extrapolation error.^39–42 

The well-established NN surrogate model enables the construction of a continuous data map based on a discrete ensemble of datasets derived from sparse PIC simulations. Following the training of our NN surrogate model using the 355 input PIC datasets, we leverage it to conduct an extensive parameter sweep, accompanied by the corresponding predicted proton energy E_k, as depicted in Fig. 6. In the figure, a discernible layer becomes evident, signifying the presence of an optimal parameter set. This narrow layer exhibits two distinctive characteristics:¹ an increase in the maximum proton energy E_k with rising a_L;² a symmetric hyperbolic layer presents due to the structure of Eq. (1), but the predicted energy distribution is asymmetric on the hyperbolic layer. For example, in the hyperbolic curve at a_L = 20, it is clear that the higher energy (red) and the lower energy (yellow-green) are obtained by the thin-high density target and the thick-low density target, respectively.

Considering $( n 0 / n c ) ( D / λ L )$ as a variable, the PSA formula can be represented as a linear equation with a slope of π, as depicted in Fig. 7. It becomes evident that the PSA theory aligns well with the predictions generated by our NN surrogate model, indicating efficient proton acceleration when Eq. (1) is satisfied. Consider a_L = 20 and $( n 0 / n c ) ( D / λ L ) = 7$ to illustrate. In this case, increasing target density while reducing thickness is equivalent to decreasing target density while increasing thickness. However, when comparing Fig. 7 with Fig. 6, it becomes apparent that the symmetry property of the parameter set exists only within a limited region of the parameter space. This implies that some information remains obscured in the representation of Fig. 7.

To compare the predictions generated by the NN surrogate model with the input PIC data, as illustrated in Fig. 2, we have performed parameter slice analyses on the parameter cube $( a L , D / λ L , n 0 / n c )$ for values of $a L =$ 5, 10, 15, and 20, as demonstrated in Fig. 8. Some filament patterns in the region below the black lines are observed because of the sparse input PIC data.

For a specific slice of a_L = 5, the NN surrogate model predicts a maximum proton energy of approximately E_k = 1200 MeV at the coordinates $( D / λ L , n 0 / n c ) = ( 0.40 , 4.02 )$⁠. To validate this prediction, an additional PIC simulation was conducted using these parameter settings, resulting in an energy measurement of $E k = 664.5$ MeV. This observation highlights the potential for the NN surrogate model to occasionally overestimate or underestimate proton energy, aligning with the presentation in Fig. 5. In addition, Fig. 8 showcases a promising level of agreement between the predicted values and the PSA theory, particularly in cases, where $D / λ L < 1$⁠. Nevertheless, certain discrepancies between the predictions and the theoretical expectations are apparent especially in the region of $D / λ L > 1$⁠, warranting further investigation.

In the case of $n 0 = 1.24 n PSA$⁠, the density distribution exhibits a double-layer structure, where the electron layer is propelled forward by radiation pressure, while the proton layer is consistently drawn by the established electrostatic field, as illustrated in Fig. 9(a). We introduce a negative sign to the electron density to differentiate between the proton and electron layers. The accumulated densities of the proton and electron layers (⁠ $30 n c$⁠) are much higher than the initial target density. Moreover, the mono-energetic proton energy spectrum with a central energy of 1300 MeV is observed which is evidence of the efficient acceleration of the proton by the PSA mechanism. Conversely, for $n 0 = 1.23 n PSA$⁠, the density distribution demonstrates layer destruction without a mono-energetic proton energy spectrum, leading to the failure of the PSA mechanism, as depicted in Fig. 9(b).

To elucidate the underlying physical mechanisms leading to the success or failure of the PSA mechanism, we examine the density distribution snapshots for both scenarios, as depicted in Fig. 10.

Figure 11 showcases the density distributions at t = 105 fs for both $n 0 = 1.24 n PSA$ and $n 0 = 1.23 n PSA$⁠. In this depiction, we define D as the initial target thickness and $D ′$ as the pulse-compressed electron layer's thickness. It is important to acknowledge that precisely defining $D ′$ proves challenging due to its dynamic nature throughout the pulse-pushing process. In Fig. 11, the designated region represents the skin depth (depicted in gray), computed using the maximum electron density. A comparison is made between the calculated skin depth l_s and the compressed electron layer thickness $D ′$⁠. For $n 0 = 1.24 n PSA$⁠, the compressed electron layer thickness $D ′$ surpasses the calculated skin depth l_s. As a result, the laser field (depicted by the black dashed line) experiences a steep drop to near-zero values after traversing the electron layer (⁠ $x > 4.8 λ L$⁠), as exemplified in Fig. 11(a). This occurrence signifies complete reflection of the laser field, effectively channeling energy into the electron layer and indirectly facilitating proton bunch acceleration via the electrostatic field between the double-layer structure. Conversely, for $n 0 = 1.23 n PSA$⁠, the dimensions of the sharp electron density structure prove smaller than the skin depth l_s, thereby enabling the transmission of the laser field (indicated by the black dashed line) past the electron layer (⁠ $x > 4.5 λ L$⁠), as depicted in Fig. 11(b). Consequently, the radiation pressure resulting from partial reflection fails to counterbalance the electrostatic force originating from the sharp electron density structure.

By employing the PSA mechanism as the benchmark, we have already ensured that the DL can extract the physical feature instead of performing numerous consumed-cost PIC if the NN model is trained appropriately. However, the situation is predictably different from 1D to 2D/3D PIC simulations because a critical parameter, focal spot size, is included. Even though we obtain the optimal conditions based on the 1D PSA scheme and apply the parameter set to 2D/3D PIC simulations, the 2D/3D simulation result would not present the same performance (e.g., energy, spectra, or density distributions) as the 1D PIC simulation since the transverse effects or instabilities cause the damage of the accumulated electron layer, such as Rayleigh-Taylor instability,^20,43 Weibel instability,^44,45 half-wavelength instability,⁴⁶ and transverse instability.^21–23 Although the 1D scheme cannot be directly compared with 2D/3D PIC, 1D theory still plays an important role in establishing the physical insight for the complex process of laser-driven proton acceleration. Hence, we expect that the NN surrogate model using sparse training data based on this paper will be a powerful tool to allow us to figure out the dependencies between the accelerated proton energy and each parameter.

After studying Fuchs' scaling by DL, we will switch the 2D scenarios for the hybrid acceleration scheme proposed by Isayama.⁴⁷ This hybrid acceleration scheme is a combination acceleration scheme of RPA,^8–11 laser wakefield acceleration (LWFA),^48,49 and TNSA^4–6 by the interaction of the tabletop dual pulses with a tandem solid density and a near-critical density target. By choosing an appropriate parameter set, the achievement of 100 MeV (or even higher) protons with a tabletop laser system is possible, e.g., NCU 100-TW laser system.⁵⁰ In this hybrid acceleration scheme,⁴⁷ multiple parameters are necessary to consider for the optimization of the maximum proton energy, such as the individual laser parameters (energy, duration, spot size), the pulse time delay between two lasers, the density and the size of the solid target and the near-critical density target, and the space between the two targets.

We developed a neural network (NN) surrogate model, constructed from an ensemble of 355 one-dimensional (1D) particle-in-cell (PIC) simulations, to investigate the phenomenon of phase-stable acceleration (PSA) in laser-driven ion acceleration. This work highlights that even with a sparse and discrete dataset, a well-trained NN model can effectively generate a continuous data map, enabling accurate predictions. Significantly, the widely accepted PSA theory serves as a benchmark for evaluating the reliability and robustness of the predictions made by the trained NN model. Moreover, our NN surrogate model not only validates the well-established PSA theory but also uncovers a noteworthy disparity between theory and predictions when considering targets with larger thicknesses. Through comprehensive 1D PIC simulations, we identify that the collapse of the electron layer occurs when the plasma skin depth surpasses the thickness of the pulse-compressed electron layer during the interaction process, causing the termination of the PSA process. Although our study predominantly focuses on the 1D effects, it is essential to underscore the significance of the skin depth in multi-dimensional simulations. The skin depth plays a pivotal role in determining the dynamics of laser–plasma interactions and the mechanisms governing ion acceleration in both one and higher-dimensional scenarios.

The integration of deep learning (DL) methodologies into multi-dimensional PIC simulations or experiments holds great promise for advancing the field of laser-driven ion acceleration. Future investigations will leverage DL techniques to explore intricate phenomena, such as Rayleigh–Taylor instability,^20,43 Weibel instability,^44,45 half-wavelength instability⁴⁶ and transverse instability,^21–23,51 measurement technique,^52,53 target design,^54,55 and novel acceleration schemes,^56,57 especially the hybrid acceleration scheme proposed by Isayama as mentioned in future work.⁴⁷ These breakthroughs will have far-reaching implications in areas, such as cancer therapy, fundamental science research, and other practical applications, relying on high-energy ion beams.

See the supplementary material for two GIF video files for n 0 = 1.24 n PSA and n 0 = 1.23 n PSA in Fig. 10. There are three subfigures in each GIF file: the first is the electrostatic field (Ex), second is the laser field (Ey), and third is the electron and the proton densities ( n / n c).

This work was supported by the National Science and Technology Council (Ministry of Science and Technology), Taiwan, under Grant Nos. MOST 111-2112-M-006-011-MY3, 111-2111-M-006-005-MY2, and 112-2811-M-006-027. This work was also supported by JSPS KAKENHI under Grant Nos. 19H00668, 20KK0064, and 22H01195 and Core-to-Core Program (No. JPJSCCA2019002). We thank the National Center for High-performance Computing (NCHC) for providing computational and storage resources.

The authors have no conflicts to disclose.

Yao-Li Liu: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (equal); Project administration (lead); Resources (lead); Software (equal); Supervision (equal); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Yen-Chen Chen: Data curation (equal); Formal analysis (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Chun-Sung Jao: Formal analysis (supporting); Methodology (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Mao-Syun Wong: Data curation (supporting); Methodology (supporting); Software (supporting); Visualization (supporting). Chun-Han Huang: Data curation (supporting); Formal analysis (supporting); Software (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Han-Wei Chen: Data curation (supporting); Formal analysis (supporting); Methodology (supporting); Software (supporting); Visualization (supporting). Shogo Isayama: Writing – original draft (supporting); Writing – review & editing (supporting). Yasuhiro Kuramitsu: Conceptualization (supporting); Project administration (equal); Supervision (equal); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.