2 Experimental setup

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^[Reference Xu, Streeter and Ettlinger^35].

Figure 1 Illustration of the experimental setup, showing the orientation of the laser–plasma interaction and the main diagnostics. The laser was focused, with an f/2.5 90° diamond-turned off-axis parabolic mirror (OAP), to a 1.6 μm radius focal spot containing a median of $35\%\pm 3\%$ of pulse energy. The plane of the laser–plasma interaction was monitored by imaging self-emission at 800 nm at ${60}^{\circ }$ to the laser propagation axis.

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^[ Reference Margarone, Krsa, Giuffrida, Picciotto, Torrisi, Nowak, Musumeci, Velyhan, Prokpek, Lska, Mocek, Ullschmied and Rus^{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^[ Reference Li, Tsubakimoto, Yoshida, Nakata and Miyanaga^{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.

3 Automated grid scans

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^[ Reference Dover, Palmer, Streeter, Ahmed, Albertazzi, Borghesi, Carroll, Fuchs, Heathcote, Hilz, Kakolee, Kar, Kodama, Kon, MacLellan, McKenna, Nagel, Neely, Notley, Nakatsutsumi, Prasad, Scott, Tampo, Zepf, Schreiber and Najmudin^{49 ]}, suggesting that the sheath field is dynamic during the acceleration process.

Figure 2 (a) Proton and (b) electron energy spectra from the rear side of the target during an automated target position scan ( ${z}_{\mathrm{T}}$ ) with a 12 μm Kapton tape and an on-target laser energy of $438\pm 32$ mJ. (c) Average proton spectra (and standard deviation) for different ${z}_{\mathrm{T}}$ positions as indicated in the legend. The proton spectra are recorded by the time-of-flight diamond detector. Each column of the waterfall plots is the average of the 10 shots from each burst. The scan comprises 31 bursts at different target positions spaced at 7.3 μm intervals along the laser propagation axis. Negative values of ${z}_{\mathrm{T}}$ are when the target plane is closer to the incoming laser pulse and ${z}_{\mathrm{T}}=0$ is the target at the best focus of the laser pulse. The magenta data points, connected with a guide line, indicate the burst-averaged 95th percentile energy as well as the standard deviation of this value across the burst.

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^[ Reference Green, Carroll, Brenner, Dromey, Foster, Kar, Li, Markey, Mckenna, Neely, Robinson, Streeter, Tolley, Wahlström, Xu and Zepf^{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^[ Reference Lindau, Lundh, Persson, Mckenna, Osvay, Batani and Wahlström^{23 ,} Reference Lundh, Lindau, Persson, Wahlström, Mckenna and Batani^{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^[ Reference Hooker, Tang, Chekhlov, Collier, Divall, Ertel, Hawkes, Parry and Rajeev^{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^[ Reference Kaluza, Schreiber, Santala, Tsakiris, Eidmann, Meyer-Ter-Vehn and Witte^{22 ,} Reference Lundh, Lindau, Persson, Wahlström, Mckenna and Batani^{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^[ Reference Morrison, Feister, Frische, Austin, Ngirmang, Murphy, Orban, Chowdhury and Roquemore^{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^[ Reference Kaluza, Schreiber, Santala, Tsakiris, Eidmann, Meyer-Ter-Vehn and Witte^{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.

Figure 3 One-dimensional scans of (a) and (c) target z-position ${z}_{\mathrm{T}}$ and (b) and (d) astigmatism ${Z}_2^{-2}$ for 12 μm thickness Kapton tape and a pre-plasma laser energy of $453\pm 40$ mJ. The electron and proton flux are plotted in (a) and (b), and the specularly reflected fundamental and second harmonic laser signals are plotted in (c) and (d). All fluxes are normalized to their observed maxima over the 2D parameter scans. Two-dimensional scans of electron and proton flux are shown in (e) and (f), with the average detected electron energy and the maximum (95th percentile) proton energies shown in (g) and (h), respectively. The 2D scan is a result of 143 bursts of 15 shots and the datapoints are the mean of each individual burst.

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^[ Reference Chopineau, Leblanc, Blaclard, Denoeud, Thévenet, Vay, Bonnaud, Martin, Vincenti and Quéré^{17 ,} Reference Pirozhkov, Choi, Sung, Lee, Yu, Jeong, Kim, Hafz, Kim, Pae, Noh, Ko, Lee, Robinson, Foster, Hawkes, Streeter, Spindloe, McKenna, Carroll, Wahlström, Zepf, Adams, Dromey, Markey, Kar, Li, Xu, Nagatomo, Mori, Yogo, Kiriyama, Ogura, Sagisaka, Orimo, Nishiuchi, Sugiyama, Esirkepov, Okada, Kondo, Kanazawa, Nakai, Akutsu, Motomura, Tanoue, Shimomura, Ikegami, Daito, Kando, Kameshima, Bolton, Bulanov, Daido and Neely^{56 ,} Reference Streeter, Foster, Cameron, Borghesi, Brenner, Carroll, Divall, Dover, Dromey, Gallegos, Green, Hawkes, Hooker, Kar, McKenna, Nagel, Najmudin, Palmer, Prasad, Quinn, Rajeev, Robinson, Romagnani, Schreiber, Spindloe, Ter-Avetisyan, Tresca, Zepf and Neely^{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^[ Reference Wilson, King, Butler, Carroll, Frazer, Duff, Higginson, Dance, Jarrett, Davidson, Armstrong, Liu, Hawkes, Clarke, Neely, Gray and McKenna^{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^[ Reference Tayyab, Bagchi, Chakera, Khan and Naik^{32 ,} Reference Ziegler, Albach, Bernert, Bock, Brack, Cowan, Dover, Garten, Gaus, Gebhardt, Goethel, Helbig, Irman, Kiriyama, Kluge, Kon, Kraft, Kroll, Loeser, Metzkes-Ng, Nishiuchi, Obst-Huebl, Püschel, Rehwald, Schlenvoigt, Schramm and Zeil^33].

Figure 4 Laser pulse temporal profiles as measured by the on-shot SPIDER diagnostic for the results of the 1D scan (Figure 2), 2D scan (Figure 3) and optimization (Figure 5). The integrals of the signals are set by independent measurements of the on-target laser energy, which were $438\pm 32$ mJ (1D scan), $453\pm 40$ mJ (2D scan) and $258\pm 22$ mJ (optimization). The corresponding measured FWHM pulse widths were $49\pm 3$ fs, $45\pm 4$ fs and $39\pm 1$ fs.

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.

4 Bayesian optimization

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)^[ Reference Rasmussen and Williams^{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^[ Reference Shalloo, Dann, Gruse, Underwood, Antoine, Arran, Backhouse, Baird, Balcazar, Bourgeois, Cardarelli, Hatfield, Kang, Krushelnick, Mangles, Murphy, Lu, Osterhoff, Põder, Rajeev, Ridgers, Rozario, Selwood, Shahani, Symes, Thomas, Thornton, Najmudin and Streeter^{43 ,} Reference Jalas, Kirchen, Messner, Winkler, Hübner, Dirkwinkel, Schnepp, Lehe and Maier^{44 ]} and in simulated laser-driven ion beams for maximizing proton energy^[ Reference Dolier, King, Wilson, Gray and McKenna^{45 ]}. To demonstrate its applicability in a laser-driven ion acceleration experiment, we have adapted the algorithm from Shalloo et al. ^[ Reference Shalloo, Dann, Gruse, Underwood, Antoine, Arran, Backhouse, Baird, Balcazar, Bourgeois, Cardarelli, Hatfield, Kang, Krushelnick, Mangles, Murphy, Lu, Osterhoff, Põder, Rajeev, Ridgers, Rozario, Selwood, Shahani, Symes, Thomas, Thornton, Najmudin and Streeter^{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.

Figure 5 Optimization of the 95th percentile proton energy determined by the rear-surface time-of-flight diagnostic through the adjustment of the laser wavefront and position of target along the laser propagation direction ( ${z}_{\mathrm{T}}$ ). The top panel shows the measured values of the proton energy (median and median absolute difference of each burst) as a function of the burst number (black points and error bars, respectively), together with the model predicted optimum after each burst (red line and shaded region) as well as the final optimal value from the model (blue horizontal line). The variation of each control parameter (given in micrometres) is shown in the lower plots (black points) along with the final optimized values (blue horizontal line), also as functions of the burst number. The best individual burst is indicated by the vertical magenta line in each plot and it can be seen that, for all parameters, the experimental parameters fall very close to the optimum value predicted by the model (e.g., they are close to the horizontal blue line). For this data series, each burst contained 20 shots, the target was 12 μm Kapton tape and the laser energy was $258\pm 22$ mJ.

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.

Figure 6 Reconstructed laser intensity profiles at ${z}_{\mathrm{T}}=0\;\unicode{x3bc} \mathrm{m}$ for (a) ${Z}_2^{-2}=-1.2\;\unicode{x3bc} \mathrm{m}$ , (b) ${Z}_2^{-2}=0\;\unicode{x3bc} \mathrm{m}$ , (c) ${Z}_2^{-2}=1.2$ μm and (d) for the optimal pulse (burst 53) from the optimization shown in Figure 5. The peak intensity of each focus was 2.7 × 10¹⁹, 5.1 × 10¹⁹, 2.9 × 10¹⁹ and 3.2 × 10¹⁹ W cm⁻², respectively.

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.