Laser–target interactions have been demonstrated to provide a highly versatile source of secondary radiation, of interest for many applications^[1^,2^] as well as the study of fundamental science^[3–5^]. Specifically, laser-driven ion accelerators^[6^,7^] have desirable characteristics pertaining to applications in medicine^[8^], material science^[9^], nuclear fusion^[10^] and imaging^[11–13^]. The need for stable, reproducible beams that can be tuned presents a necessary, yet challenging, goal towards the realization of many of these applications^[6^].

The most extensively researched mechanism driving laser-driven proton acceleration is sheath acceleration (often termed target normal sheath acceleration, or TNSA)^[14^]. Here, a laser pulse is focused to a relativistic intensity on the surface of a solid foil, ionizing the material and heating electrons to MeV temperatures. As the accelerated electrons escape the rear target surface, they build an electrostatic sheath field that can reach approximately more than 1 TV/m and accelerate ions to much greater than 1 MeV energies over just a few micrometres^[15^]. Although target materials vary, the accelerated ions are typically dominated by protons from hydrocarbon surface contaminants, which are preferentially accelerated due to their high charge–mass ratio^[16^].

Due to their dependence on the rear-surface electrostatic sheath field, the specific characteristics of these MeV proton beams are strongly influenced by the laser–electron energy coupling, electron transport through the bulk of the target and disruption to the target rear surface. Experiments have demonstrated the dependence of the electron and proton beam on various experimental control parameters (e.g., laser intensity and contrast or target thickness). Of particular importance is the plasma scale length at the front surface, which affects the laser–electron coupling mechanisms^[17–19^]. In pre-plasmas with long scale lengths (>100 μm) the laser beam has been observed to filament, subsequently reducing the coupling efficiency, while for optimal scale lengths the beam can undergo relativistic self-focusing effects that enhance laser energy coupling^[18–20^]. The amplified spontaneous emission (ASE) pedestal and prepulses, common to short pulse lasers, can pre-heat the target and lead to significant plasma expansion before the main pulse arrives. For a given ASE pedestal duration, an optimal thickness exists at which this enhanced coupling and electron recirculation^[21^] will be advantageous for the acceleration process, while for thinner targets the inward travelling shock-wave launched by the rapid surface pre-heating can disrupt the accelerating sheath field at the rear surface^[22–24^]. Radiative heating from X-rays, generated in the focus of a prepulse incident on the target front-side, can similarly induce rear-surface expansion of thin targets, impacting TNSA for interaction parameters in which the ASE-induced shock may not reach the rear surface during the acceleration window^[22^].

While broad trends within TNSA proton beams have been established, comparison of experimental results from different experiments highlights variation in measured beam parameters and is indicative of the nuanced relationship between laser parameters and proton beam characteristics. For example, while maximum proton energy has been observed to increase with laser intensity, the scaling follows an approximately ${I}^{1/2}$ dependence for long ( ${>}300$ fs) duration laser pulses and a linear dependence of approximately $I$ for ultra-short (40–150 fs) pulses^[7^,25^,26^]. A number of numerical and experimental studies have explored the impact of laser pulse duration on maximum proton energy due to TNSA in interactions with moderate laser contrast. These indicate an optimal pulse duration for proton acceleration associated with a fixed laser intensity or laser energy^[27–30^] (e.g., optimal duration between 100 and 300 fs for laser pulse energies of 1 J). The dependence is attributed to differences in laser energy to electron coupling efficiency and acceleration time relative to the rear-surface expansion timescale. For higher laser contrast, a similar trend is observed for targets with thicknesses of tens of micrometres^[31^]. More advanced temporal pulse shaping (e.g., shaping of the rising and falling edge of the pulse) has been explored recently for the interaction of high-contrast lasers with ultra-thin targets, which indicates significant enhancement of proton maximum energies over those observed for best laser compression^[32^,33^].

With the proliferation of multi-Hz high-power laser pulses^[34^] and the development of high repetition-rate (HRR)-compatible solid-density targetry^[35–41^], it is now possible to quickly obtain large datasets from laser-driven ion acceleration experiments. This opens up the possibility to perform extensive multi-dimensional parameter scans to elucidate the interdependence of different experimental control parameters, as well as to apply machine learning techniques to optimize ion beam properties – within complex multi-dimensional parameter spaces – in automated experiments^[42–44^] and simulations^[45^,46^].

Here, we describe the first experimental demonstration of real-time Bayesian optimization (BO) of a laser-driven ion source, using a closed-loop algorithm. The fully automated control system operated the laser, analysed the diagnostic results and made changes to the experiment control parameters. This enabled rapid and efficient optimization of the accelerator performance through simultaneous tuning of up to six different input parameters, producing proton beams with equivalent peak energy using 57% of the laser energy of the manually optimized interaction.

The experiment (see Figure 1 for the setup) was performed at the Gemini TA2 facility, using a Ti:sapphire laser that forms part of the Central Laser Facility at the Rutherford Appleton Laboratory. The laser pulses contained up to 500 mJ in a transform-limited pulse duration full width at half maximum (FWHM) of approximately 40 fs, with a central wavelength of 800 nm and an FWHM bandwidth of 30 nm. The laser was focused to a high intensity ( ${I}_{\mathrm{L}}>{10}^{19}$ W/cm²) using an $f/2.5$ off-axis parabolic mirror and interacted with the target at an angle of incidence of 30° with p-polarization. The target was Kapton tape of 12 μm, spooled continuously during shots using a motorized tape drive^[35^].

The interaction was diagnosed using a suite of particle diagnostics, including a scintillator (EJ-440), positioned along the rear-surface target normal, to measure the proton spatial profile, two point measurements of the proton energy spectrum using a time-of-flight (TOF) diamond detector^[47^] and fibre-coupled Thomson parabola spectrometer (at 3° to the target normal and along the target normal axis respectively) and a 0.15 T permanent magnet electron spectrometer in the laser-forward direction. The near-field of the specularly reflected laser light was also measured at the first and second harmonics of the drive laser.

The laser spectral phase was controlled by an acousto-optic programmable dispersive filter (DAZZLER) and measured using a small central sample of the compressed pulse and a SPIDER diagnostic. The laser parameters (six wavefront aberrations generated through Zernike polynomials, second-, third- and fourth-order temporal phases, energy and polarization) and target position relative to the laser focus were controlled using a fully automated control and acquisition system. This enabled data scans consisted of bursts of shots (up to 20) at fixed input values in parameter space. Following a burst of shots, the control code performed analysis of the measured data from the online diagnostics and adjusted the laser or target parameters for the next burst. Although the controls enabled adjustment of the laser temporal pulse shape, the data presented here corresponds to pulse shapes close to best compression ( $\sim 40$ fs) with variations due to the day-to-day variation in laser tuning.

For high-energy, multi-Hz laser facilities, prolonged HRR operation can affect the laser pulse parameters, leading to a degradation in peak intensity^[48^]. To ensure that our setup was not subject to these effects, measurements of laser parameters were made over periods of extended 1 Hz operation. These measurements concluded that the effect of prolonged HRR operation on the quality of the temporal pulse shape was negligible, with the standard deviation of the random fluctuation in pulse FWHM measured as approximately 3 fs over 1400 shots.

With the automated setup, parameter scans can be readily obtained by following a pre-programmed procedure. In doing this, the control algorithm moved through a series of equally spaced locations, taking a number of repeat shots at each configuration to quantify shot-to-shot fluctuations. Figure 2 shows proton and electron spectra for a 1D scan of the target position through the laser focus with a 12 μm Kapton tape and a pulse length of ${\tau}_{\mathrm{FWHM}}=49\pm 3$ fs. The burst-averaged 95th percentile proton energy (hereafter referred to as the maximum energy) and average electron energies are overlaid in magenta with the standard deviation for each burst indicated by error bars. The proton and electron spectra are seen to extend to higher energies as the target position approaches the laser best focus, as would be expected due to the increasing laser intensity at the target surface. While the electron spectra peak around the best focus, where the laser intensity is highest, a characteristic dip in the maximum proton energy and flux is observed. Around the best focus ( $\left|{z}_{\mathrm{T}}\right|<25\;\unicode{x3bc} \mathrm{m}$ ), a comparatively small number of protons are still observed at high ( $\approx 3.5$ MeV) energy, with the spectrum dominated by lower (sub-MeV) energies, as seen in Figure 2(c). A second signal is also seen in the TOF spectrum, appearing as a band peaking at approximately equal to 0.5 MeV per nucleon in Figure 2(a). This is most likely due to heavy ions that were accelerated to lower velocities due to their lower charge-to-mass ratio. Together with the spectrally peaked proton spectrum, this appears similar to observations of ‘buffered’ proton acceleration for higher intensities and thinner targets^[49^], suggesting that the sheath field is dynamic during the acceleration process.

An increase in the laser focal spot size has previously been observed to increase the number of accelerated particles, although with a reduced maximum energy^[50^]. While this is consistent with our measurements, the strong suppression of proton flux at the highest intensity may indicate that the acceleration process is further compromised at the highest laser intensities by the contrast levels of our laser, with the prepulses and amplified spontaneous emission (ASE) causing adverse pre-heating of the target. Similar disruption has previously been attributed to rear-surface deformation by ASE-driven shock break-out, which can effectively steer a high-energy component of the proton beam emission towards the laser axis^[23^,24^], modifying the spectrum measured at a single angular position, or the presence of a long scale-length plasma on the rear surface, which has been shown to suppress the production of ions through TNSA in experiments and simulations^[51–53^].

A laser pulse contrast measurement (using an Amplitude Sequoia) showed an ASE intensity contrast of better than ${10}^{-9}$ up to $t=-20$ ps, after which the laser intensity in the coherent pedestal^[54^] increased exponentially. Individual prepulses with a relative intensity of ${10}^{-6}$ were also observed between $t=-50$ ps and $t=-65$ ps. The measured contrast, starting at $t=-150$ ps, was used to perform 2D cylindrical hydrodynamic modelling of the target evolution ahead of the arrival of the peak intensity. The modelling was performed using the FLASH code (v4.6.2). The ASE and coherent pedestal from $t=-150$ ps to $t=-1$ ps were coupled to the target electrons using ray-tracing with inverse bremsstrahlung heating, and Lee–More conductivity and heat exchange models were used. This indicated the formation of an approximately 2 μm exponential scale-length pre-plasma. For the target thicknesses used in the experiment ( $\gg$ 1 μm), the measured prepulse was not large enough for the generation of a shock moving quickly enough to perturb the density step of the target rear surface. This matches previous results^[22^,24^] in similar interaction conditions, which indicates that the ablation-launched density shock does not have time to affect the rear surface during the acceleration process. This dip in signal at the highest intensities is a surprising result for targets with tens of micrometre thickness. It is more commonly observed for experiments using ultra-thin targets, where it is attributed to ASE shock break-out^[55^]. For the interaction presented here, the rear surface may be affected by poor long-timescale contrast (before the start of the measurement window at $t=-150$ ps), or through X-ray heating of the target bulk^[22^]. Determination of the specific processes driving the disruption of proton acceleration for this interaction requires additional measurements of long-timescale contrast and pre-plasma scale length, and is beyond the scope of this optimization demonstration.

The laser intensity can be varied by adding wavefront aberrations to the laser, which changes the focal spot shape as well as the peak intensity. Figures 3(a)–3(d) show the burst-averaged electron and proton flux as well as the relative specular reflectivity at the fundamental and second harmonic wavelengths for a varying target plane, ${z}_{\mathrm{T}}$ (Figures 3(a) and 3(c)) and 45-degree astigmatism, ${Z}_2^{-2}$ (Figures 3(b) and 3(d)). Again, a characteristic drop was observed in the proton flux and maximum energy for the best focus. Moderately decreasing the peak intensity by either defocusing or adding astigmatism decreased the electron flux and average energy, but maximized the proton acceleration. At the low intensity limit the particle acceleration drops to zero, as expected.

Evidence of disruption to the front surface during the interaction can be seen from the sharp drop in the fundamental and second harmonic laser reflectivity at high intensities. This matches previous observations of target reflectivity and harmonic generation being adversely affected for high-intensity low-contrast laser–plasma interactions as a result of the formation of a large scale-length pre-plasma^[17^,56^,57^].

To explore the interplay between astigmatism and defocus, an automated 2D grid scan was performed. Burst-averaged measurements of the electron and proton flux, mean electron energy and maximum proton energy are displayed in Figures 3(e)–3(h). The electron generation was maximized for both zero defocus and zero astigmatism, monotonically decreasing as ${z}_{\mathrm{T}}$ and ${Z}_2^{-2}$ were increased. The proton flux and energy are approximately maximized for a ring around the origin, indicating a threshold intensity for disrupting the proton acceleration that can be achieved either by defocusing or increasing the focal spot size through optical aberrations. There also appears to be enhanced proton flux for zero defocus, ${z}_{\mathrm{T}}=0$ μm, but with the application of significant astigmatism ${Z}_2^{-2}=\pm 1.2\;\unicode{x3bc} \mathrm{m}$ when compared with the increased flux achieved just through defocusing with no astigmatism. This indicates that the proton acceleration process is not just intensity dependent, but rather is also sensitive to the spatial intensity profile on-target^[58^].

For all results in Figure 3, the laser pulse had a shorter pulse duration of ${\tau}_{\mathrm{FWHM}}=46\pm 4$ fs and was skewed with a slower rising edge than for the results in Figure 2, as can be seen in Figure 4. The pulse shape appears to affect the range of ${z}_{\mathrm{T}}$ over which the proton acceleration was suppressed. For the case of Figure 2, proton acceleration was maximized at ${z}_{\mathrm{T}}=\pm 30\;\unicode{x3bc} \mathrm{m}$ , while with the slower rising edge used for Figure 3(a), the maximum occurs at ${z}_{\mathrm{T}}=\pm 67\;\unicode{x3bc} \mathrm{m}$ , with very low flux obtained at ${z}_{\mathrm{T}}=\pm 33\;\unicode{x3bc} \mathrm{m}$ . For a slower rising edge, there is more time for any disruption of the target to occur prior to the arrival of the peak of the pulse, meaning that a larger defocus is required to maintain proton acceleration. A lower maximum proton energy of $0.8\pm 0.1$ MeV is observed in Figure 3(a) compared to $2.8\pm 0.4$ MeV for Figure 2. This shows the benefit of using temporal pulse shaping to minimize target pre-heating, as it allows for efficient proton acceleration at higher intensity interaction, leading to higher energy protons. Despite the very different interaction conditions (here – moderate laser contrast and micrometre thick target), this indicates a similar trend to that observed in recent experiments using high-contrast ultra-thin targets that demonstrated the enhancement of particle energies and numbers by modification of the pulse shape, away from nominal pulse compression, with a steepened rising edge in comparison with the falling edge of the pulse^[32^,33^].

The interaction parameter space is multi-dimensional. While 1D and 2D slices of this parameter space provide valuable insights, true mapping of the parameter space is required for deeper understanding and control of sheath-accelerated proton beams, and the location of global optima requires the inclusion of more dimensions. This is particularly important when individual parameters are coupled in complex relationships, as is evident in the case of proton acceleration in Figure 3. While grid scanning is feasible for mapping up to 2D slices of parameter space with multi-Hz laser systems, for higher numbers of dimensions it becomes prohibitively time consuming. In addition, a grid scan covers regions of parameters with a high signal and no signal with the same resolution. Given that measurements in large regions of the accessible parameter space will return no appreciable proton acceleration, this is uneconomical in terms of target usage, debris production and laser operation. It is therefore desirable to use more intelligent algorithms for probing and optimization of the beam in high-dimensional parameter spaces.

In BO, experimental data is used to update a prior model to more accurately fit observations, thereby obtaining a posterior model. The model is then used to make predictions over the experimental parameter space and select the parameter set for the next measurement, with the goal of efficiently finding the optimum within the parameter space. A commonly used type of model is Gaussian process regression (GPR)^[59^], which is a well-suited technique for modelling multi-dimensional experimental data. A key advantage for experimental science is that GPR can naturally include uncertainty quantification on the data, and also can be used to estimate the uncertainty when making predictions. BO is widely used for the optimization of noisy processes that are expensive to evaluate and for which there is no adequate analytical description, as is typically the case in complex non-linear systems.

In the field of laser–plasma acceleration, BO has been used in laser-wakefield acceleration to optimize the generated electron and X-ray beam properties^[43^,44^] and in simulated laser-driven ion beams for maximizing proton energy^[45^]. To demonstrate its applicability in a laser-driven ion acceleration experiment, we have adapted the algorithm from Shalloo et al.^[43^]. Using the automated control of the experiment and online analysis of the experimental diagnostics, we were able to perform multi-dimensional optimization of any fitness function that outputted a scalar property of interest, such as the maximum proton energy.

The input values passed to the model for each burst were taken from the parameters set by the control algorithm, with the exception of the plane of the interaction, ${z}_{\mathrm{T}}$ . It was found that errors in positioning the tape target, while significantly smaller than the Rayleigh range of the focusing laser (15 μm), were large compared to the sensitivity of the plasma accelerator, and so the target position input values were taken from a spatially resolved measurement of the self-emission region at the target surface, collected at 60° to the front surface normal of the tape. This was found to greatly improve the confidence of the model and its ability to find the optimum.

To demonstrate the BO algorithm, a 6D optimization was performed to maximize the maximum proton energy (as measured by the TOF). Five Zernike mode coefficients – namely ${Z}_2^{-2}$ , ${Z}_2^2$ (oblique and vertical astigmatism), ${Z}_3^{-1}$ , ${Z}_3^1$ (vertical and horizontal coma) and ${Z}_4^{-2}$ (oblique second astigmatism) – were used to change the spatial phase of the laser pulse, affecting the focal spot shape and peak intensity. In addition, the tape surface position relative to the focal plane of the laser, ${z}_{\mathrm{T}}$ , was varied using a linearly motorized stage and an online self-emission diagnostic to measure its position, as mentioned. This optimization started from an initially flat wavefront (optimized using a feedback loop with a HASO wavefront sensor) and the tape target was initially positioned at the estimated focal plane of the laser. This represents the typical starting point for conventional optimization of manual experiments, with the highest intensity being assumed to be optimal.

Figure 5 shows the measured maximum proton energy along with the variation of each input parameter as a function of the burst number. By chance, one of the initial randomly selected points produced a large enhancement in maximum proton energy, with every parameter apart from ${Z}_2^2$ close to its eventual optimum value. With each additional measurement, the model gained more knowledge of the parameter space and adjusted its prediction of the global optimum of proton energy (red line) and its location in parameter space. After 61 bursts, the model optimum was $2.30\pm 0.10$ MeV, compared to a starting point of $1.22\pm 0.04$ MeV. This optimized maximum energy is close to the value shown in Figure 2 and significantly greater than that seen in Figure 3, despite being limited to only 260 mJ of on-target laser energy, compared to more than 430 mJ for the parameter scans.

The optimum found involved a shift of 70 μm from the initial position (best focus), as well as the addition of significant wavefront aberrations compared to the initially flat wavefront. The focal spot fluence distribution at the target plane was calculated from on-shot measurement of the near-field phase and fluence profiles using a HASO wavefront sensor. The resulting intensity maps are shown in Figure 6 for ${Z}_2^{-2}=-\mathrm{1,0,1}$ at zero defocus and for the optimized wavefront found through the optimization. The peak intensity was ${I}_0=5\times {10}^{19}$ W cm⁻² for the case of a flat wavefront at the focal plane. Each of the modified pulses had a similar peak intensity of ${I}_0\approx 3\times {10}^{19}$ W cm⁻², with the spot shape appearing close to an ellipse for the optimized pulse case. From analysis of the previous parameter scans it is inferred that this intensity distribution was optimal, maintaining a maximum possible intensity, while limiting disruption to the rear target surface.

The reason for the focal spot shape found in this optimization would require expensive 2D–3D numerical simulations investigating how subtle changes in the wavefront affect the various energy transfer processes and plasma dynamics in sheath acceleration. This would be valuable for understanding how to further optimize laser-driven proton acceleration, but is beyond the scope of this paper, which is focused on the utility of online optimization for proton beam parameter optimization as a particle source for applications as well as its use to identify interesting regions for deeper study within a complex non-linear system.

In conclusion, we have demonstrated the automation and optimization of laser-driven proton acceleration from a solid tape drive. The ability to take high-fidelity parameter scans in one or two dimensions will be of great benefit to the field in understanding what is a highly complex and dynamic interaction. BO of the generated proton beam was demonstrated, by using real-time analysis of experimental diagnostics to create a closed-loop system with limited human intervention. This can find optima that would elude manual optimization or single-parameter scans, due to the complex interplay among the large number of control parameters. In the case of this low-contrast interaction, the temporal pulse shape was shown to play an important role in determining the intensity threshold at which the proton acceleration process was disrupted. Including temporal and spatial pulse shaping simultaneously in the optimization process may lead to further improvement.

Automated BO can quickly find regimes of stable and optimized operation without requiring the constant attention of laser–plasma experts. It is anticipated that this development will be essential for efficient utilization of laser-driven ion acceleration for its many applications in future user facilities^[60^]. BO could also be extremely valuable for optimizing radiation pressure acceleration for which laser–plasma instabilities^[61^,62^] typically limit the acceleration process, as well as for optimizing enhanced acceleration through relativistic transparency, which has already demonstrated highly desirable near-100 MeV proton energies^[63^] and for which target evolution plays a central role. Fine tuning of the laser parameters may be able to mitigate these instabilities or further tailor the target evolution respectively, significantly enhancing the accelerated proton beam.