Progress in the development of high-power laser technology has enabled the construction of laser facilities with laser beam power exceeding the petawatt (PW) level.1 After tight focusing, laser pulses generated by such systems may achieve extreme peak intensities of 10²² W cm⁻² and higher. Although the number of PW class facilities is gradually increasing, achieving ultrahigh peak intensities remains a difficult technical challenge. In 2004, a peak intensity of 10²² W cm⁻² was demonstrated on the 0.3 PW HERCULES facility (USA),2 and, after an upgrade, this record was surpassed by an intensity of 2 × 10²² W cm⁻².3 It was a long time before similar results could be achieved elsewhere. Only toward the end of the last decade were intensities exceeding 10²² W cm⁻² obtained on the 0.3 PW J-KAREN-P (Japan),4 Texas PW (USA),5 4.0 PW CoReLS (South Korea),6 and 5.4 PW SULF (China)7 facilities. More recently, in 2021, with careful wavefront control, a higher peak intensity exceeding 10²³ W cm⁻² was reported on CoReLS.8 With these increasing intensities, their measurement is becoming more and more important and challenging, owing to the many concurrent processes that occur in the highly nonlinear regime of extreme intensity.

A rough method of measurement is to estimate the peak laser intensity based on the peak laser power P₀, which can be calculated knowing the total energy contained in the laser pulse and its temporal profile. The peak intensity I₀ may then be obtained for a given spatial distribution of the laser beam in the focal plane I(r) = I₀f(r) as $I_0=P_0/{\int }_0^{\infty }f(r)2\pi rdr$. For example, for a laser beam with Gaussian spatial profile given by $I(r)=I_0\exp [-4r^2/(D_F^2\ln 2)]$, where D_F is the focal-spot diameter at half-maximum, the expression for the peak intensity is $I_0=4P_0/(\pi D_F^2\ln 2)$. Thus, estimation of the peak intensity with this approach requires information on the beam size and its spatial profile in the focal plane. At lower laser energies, these characteristics may be measured with conventional techniques by imaging the focal plane on diagnostic cameras and beam profilers. After that, the result needs be extrapolated to the full laser energy, but it may differ substantially from the real experimental values owing to nonlinear effects. The record-making results mentioned above have been obtained by a far-field reconstruction at the beam focus, which involves the use of expectations based on mathematical prediction.5 Thus, new approaches based on direct intensity measurement, which could probe what actually happens at the focus of an ultrahigh-intensity laser beam, are of great importance.

Several methods have been proposed in this context, including the use of multiple tunneling ionization of high-Z atoms with high ionization potentials,9–11 nonlinear Compton/Thomson scattering,12–19 electron–positron pair production,20 and ponderomotive scattering of electron bunches.21 Some of these can be used to diagnose other pulse parameters (such as pulse duration and focal spot size) along with the peak intensity estimation. To achieve the desired accuracy with the proposed methods proposed, certain difficulties need to be overcome. While the major problem of the multiple tunneling ionization method at intensities above 10²² W cm⁻² is the importance of ponderomotive acceleration of high-Z ions in the laser field, the other approaches mentioned require perfect synchronization and spatial overlap of both particle and laser beams, as well as full characterization of the electron beam. At the same time, achieving the highest laser intensity requires tight focusing of the laser pulse, which results in strong gradients of the ponderomotive force that can accelerate even some heavy particles without synchronized pre-acceleration. The resulting angular energy distributions contain information about the laser pulse parameters, which in principle may be extracted as a solution of an inverse problem. There have been some investigations considering electrons as detecting particles,22–28 but, for a tightly focused laser pulse with high peak intensity, electrons are expelled from the focal volume before interaction with the main laser peak, which is one of the main disadvantages of this approach for full laser pulse characterization. In a recent work,29 it has been proposed to use protons for diagnostic purposes in the case of laser pulses with intensity from 10²¹ to 10²⁴ W cm⁻². A further development of this method is to use both types of particles (i.e., protons and electrons),30 which makes it possible to simultaneously estimate the peak intensity of a high-power laser pulse, the size of its focal spot, and its duration. In this work, we address the principal approaches to the inverse problems of focal parameter reconstruction for a tightly focused extremely intensive laser beam, based on ponderomotively accelerated protons.

Although, as mentioned above, some methods may be applied to estimate several laser parameters, they need information concerning the spatial–temporal distribution of the laser pulse, because of the dependence of the results on the focal laser distribution.28,31 It has been shown28 that the form of the electron energy spectrum depends on the focal distribution of the laser beam. Similarly, some relationship is expected between the distributions of accelerated protons and laser profiles that are studied in this paper. The inverse problem of reconstruction of laser focal profiles based on particle angular energy distributions, in turn, may be targeted with use of an artificial neural network (ANN). Here, we use convolutional neural network (CNN) architectures, which have gained wide recognition as useful tools for computer vision and image analysis tasks.32,33 Moreover, with the development of deep CNNs, such architectures have become the dominant approach for almost all recognition and detection tasks and are approaching human performance.34 In addition to general technical tasks (e.g., face and text recognition), they have also been widely used to solve various classification and regression problems in physical sciences.35–42

This article focuses on the interplay between the physics of proton dynamics in the extremely highly focused electromagnetic field and the key featured imprints on the obtained spectra, with the aim of reconstructing both the laser intensity distribution and the peak value in the focal region. In our previous work,29 we have already shown that in the regime of nonrelativistic proton acceleration, the angular characteristics of proton distributions depend on the focal spot size, while their energy characteristics are also determined by the pulse peak intensity and its duration. It has been pointed out that a nonisotropic distribution of the laser intensity in the focal plane leads to nonisotropic proton spectra, i.e., the dependence of the laser spatial distribution on the polar angle results in a corresponding dependence of the particle distribution. At the same time, the important question of how the dependence of the laser spatial distribution on the radius in the focal plane impacts on the proton distributions has not been considered, while non-Gaussian focal distributions with different energies in diffraction half-rings are sometimes observed in experiments.6,7 To analyze this effect, here, for the purpose of demonstration, we consider axially symmetric laser beams. As in our previous work,29 the proton dynamics are calculated using the nonrelativistic ponderomotive force approximation,43 and the laser pulses are simulated by Stratton–Chu integrals,44 which model laser beams with different spatial profiles focused by an off-axis parabolic mirror down to the diffraction limit. As shown in Ref. 29, in a certain range of laser intensities and focal spot sizes, the proton angular energy distribution can be scaled based on a simple expression. In this range, this allows us to prepare an augmented dataset using data calculated for only one intensity value for each focal spot size and spatial distribution. Using such a dataset, we present a neural network approach that enables reconstruction of the laser pulse focal distribution (including its quantitative characteristics such as focal spot size) as well as the laser peak intensity using the angular energy distributions of the particles. We use the calculated angular energy distributions and focal spatial profiles as samples for ANN training and testing. Although our approach is tested here for axially symmetric beams, it can be applied in more general cases with changes in the shape of the input data covering additional variables upon which the focal distribution may depend, such as the polar angle ϕ.

The method presented in this paper for reconstructing the intensity profile of a laser pulse is based on analysis of the angular energy distributions of protons after their acceleration by the diagnosed tightly focused laser beam. Figure 1 shows the main setup of this laser diagnostic method. A laser pulse focused by an off-axis parabolic mirror interacts with protons from a rarefied gas, after which those particles expelled from the focal volume are detected at certain angles counted from the direction of laser propagation, and their angular energy distributions are used to predict the focal profile of the laser beam and its peak intensity. Such a diagnostic approach is possible if the distribution of particles is determined by the parameters of the laser pulse. This limits the parameters of diagnostic gases by neglecting the effects of both the interaction of protons with the volume of the residual charge (after expelling electrons) and multiple proton scattering due to collisions of protons with hydrogen nuclei (ion–ion collisions). Both factors have been discussed in detail in Ref. 29.

We recall the main expressions for estimating the experiment parameters to satisfy the conditions for applicability of the method. The first effect is related to the insignificance of the Coulomb interaction of an individual proton with the volume charge compared with the ponderomotive force of the laser pulse. In this case, the proton concentration is limited by $n_p\ll 10^{23}{a}_{0p}^2{D}_{F\lambda }^{-8/3}{(\lambda /0.8\mu \mathrm{m})}^{-2}$ cm⁻³, where a_0p = qE/mωc is the normalized amplitude of the laser field in terms of the proton mass m, λ is the laser wavelength, and D_Fλ = D_F/λ. For the laser parameters considered in this paper (a minimum peak intensity I = 10²¹ W cm⁻² and a maximum focal spot size D_F = 4λ), the proton concentration should be less than 10¹⁷ cm⁻³. Proton–proton scattering leads to energy losses and distortion of emission angles. Taking the latter into consideration, the proton concentration and thickness l of the gas jet should be limited by $n_{p}l<1.5\times 10^{19}{(\mathrm{\Delta }{\vartheta }_{\text{er}})}^2{(2{\epsilon }_p/m_e{c}^2)}^2$ cm⁻², where Δϑ_er is the acceptable angular error in degrees, which should be as low as the accuracy of the angular measurement, and m_e is the electron mass. For example, for Δϑ_er = 1° and ɛ_p = 2 keV, the gas parameters are limited by n_pl < 10¹⁵ cm⁻². This condition can be satisfied in experiments, for instance, by using a gas jet of low concentration ($<10^{16}$ cm⁻³) and 1 mm thickness. For low values of Δϑ_er, the energy losses of protons become negligible compared with the final particle energy. This condition may require a specifically designed nozzle backed by a pressure of a few kilopascals. With optimal parameters of the nozzle and sufficiently low pressure, densities as low as 10¹⁵–10¹⁶ cm⁻³ can be obtained in a jet with a diameter of a few millimeters.45

Following the general approach presented in Ref. 29, we use the approximation of a nonrelativistic ponderomotive force to describe the interaction of a proton with a laser pulse whose intensity is in the extreme range from 10²¹ to 10²⁴ W cm⁻²:$$\frac{d{\vec{p}}_{\text{drift}}}{dt}={\vec{F}}_p=-\frac{q^2}{4m{\omega }^2}\vec{\nabla }|\vec{E}{|}^2,$$(1)where q and m are the charge and mass of the particle, and ω and $\vec{E}$ are the laser frequency and the electric component of its field, respectively. The force describes the dynamics of the average proton momentum related to the slow proton drift due to high laser intensity gradients in the case of tight laser focusing. To predict the spatial distribution of real laser pulses, the ANN must be trained on data similar to those obtained in experiments. For generating such data, we normally use a model based on Stratton–Chu integrals,44 which can simulate laser pulse focusing by an off-axis parabolic perfect mirror,46 as used in real-life experiments. Occasionally and if necessary, an arbitrary phase shift can be introduced at any point on the mirror surface. This model accurately describes particle dynamics in the intensity range from 10²¹ to 10²⁴ W cm⁻², although the upper limit of the peak intensity may be easily extended. In this work, it is set by the nonrelativistic ponderomotive force approximation, which provides accurate results up to 10²⁴ W cm⁻². However, Eq. (1) may be modified to take relativistic effects into account. The lower limit is set by factors that may affect particle dynamics, such as proton–proton scattering. The conditions on the gas jet described at the beginning of this section were obtained for 10²¹ W cm⁻². For reliable applicability at lower intensities, more challenging conditions may be required. In principle, the constraints on the gas density and the lower intensity limit could be lifted by switching to a different model in which particle dynamics were calculated via kinetic simulations accounting for collective effects.

We let the incident laser pulse propagate along the z axis and take the origin of the coordinate system at the position of the parabola focus, and thus the electric field of the laser pulse focused by the mirror (with radius ρ, off-axis angle ψ_off, and parent parabola focal length F) is represented as$$\begin{array}{cc}\hfill \vec{E}(x_1,y_1,z_1)& =\frac{1}{2F\lambda }\iint \left(i{\mathbf{A}}_e+\frac{{\mathbf{a}}_e}{k{r}_{QS}}\right)\cdot \frac{{\vec{E}}_i{e}^{\mathit{ikl}}}{r_{QS}^2}dxdy,\hfill \\ \hfill l& =r_{QS}-F(s-1)\hfill \\ \hfill & =\sqrt{r_1^2+r^2-2(x{x}_1+y{y}_1+z{z}_1)}-F(s-1),\hfill \end{array}$$(2)where ${\vec{E}}_i(x,y,z)=[E_{0x}(x,y){\vec{e}}_x+E_{0y}(x,y){\vec{e}}_y]\exp [-i(kz+\omega t+{\varphi }_0)]$ is the electric field of the laser pulse incident on the mirror, k and λ are the wavenumber and the wavelength, respectively, ${\vec{r}}_{QS}=\{\mathrm{\Delta }{x}_{QS},\mathrm{\Delta }{y}_{QS},\mathrm{\Delta }{z}_{QS}\}=\{x-x_1,y-y_1,z-z_1\}$ is the vector from the observation point Q to the point S on the mirror surface [see Fig. 2(a)], r and r₁ are the respective distances between the laser focus and the points S and Q, and the tensors$$\begin{array}{c}{\mathbf{A}}_e=\left(\begin{array}{ccc}\hfill 2F{r}_{SQ}-x\mathrm{\Delta }{x}_{QS}\hfill & \hfill -y\mathrm{\Delta }{x}_{QS}\hfill & \hfill 0\hfill \\ \hfill -x\mathrm{\Delta }{y}_{QS}\hfill & \hfill 2F{r}_{SQ}-y\mathrm{\Delta }{y}_{QS}\hfill & \hfill 0\hfill \\ \hfill x(r_{SQ}-\mathrm{\Delta }{z}_{QS})\hfill & \hfill y(r_{SQ}-\mathrm{\Delta }{z}_{QS})\hfill & \hfill 0\hfill \end{array}\right),\\ \hfill {\mathbf{a}}_e={\vec{r}}_{SQ}\cdot \left(\begin{array}{ccc}\hfill x\hfill & \hfill y\hfill & \hfill 0\hfill \end{array}\right).\hfill \end{array}$$(3)The integration area is the projection of the paraboloid surface onto the xy plane, i.e., (x − h)² + y² ≤ ρ², where h is the distance between the centers of the mirror and the corresponding parent parabola on the xy plane. Numerical integration is carried out by the Simpson method,47 and the discretization scheme resolves the variations in the shape of the incident pulse.48 Owing to the asymmetry of the integration surface with respect to the center of the paraboloid, the laser beam reflected by the mirror propagates at an angle ψ_off relative to the z axis. Therefore, it is convenient to switch to a primed coordinate system obtained by applying a rotation matrix:$$\begin{array}{c}{x}_1=x^{\prime }\cos {\psi }_{\text{off}}-z^{\prime }\sin {\psi }_{\text{off}},\\ z_1=x^{\prime }\sin {\psi }_{\text{off}}+z^{\prime }\cos {\psi }_{\text{off}},\\ y_1=y^{\prime }.\end{array}$$(4)The components of the laser field in the primed system can then be derived from Eq. (2) by the inverse transformation: {E_x′, E_y′, E_z′} = {E_x cos φ + E_z sin φ, E_y, −E_x sin φ + E_z cos φ}; the laser pulse propagates along the z′ axis. The latter representation of the laser components is used in this work to simulate laser pulses. Changing the ratio f_# = F/ρ makes it possible to model laser beams with different focal spot sizes. The larger the ratio, the wider is the focal spot.

We neglect border effects from the mirror edge, as well as the influence of the temporal profile on the laser spatial distribution. The former is possible for mirrors with a large focal length F of the parent mirror compared with the laser wavelength, i.e., kF ≫ 1,49 and although spatial–temporal couplings must be taken into account to simulate single- or few-cycle tightly focused laser pulses,50–52 these remain beyond the ranges discussed in this article.

In this paper, we demonstrate the robustness of the proposed method while applying it to laser pulses linearly polarized along the x axis (hereinafter, we omit the primes on the coordinates: x′, y′, z′) with a Gaussian temporal profile and a FWHM duration of 30 fs and one of the following spatial profiles: flat-top, Gaussian distributions of the second, fourth, and sixth orders, and the Laguerre–Gaussian mode $(L_1^0)$. The pulses are focused to a spot with a FWHM size in the range of ∼(1–4)λ. This is achieved using a parabolic mirror with radius ρ = 15 cm and off-axis angle ψ_off = 60°, with its focal length being varied to obtain the desired beam size. The laser peak intensity is set in the range of 10²¹–10²³ W cm⁻².

Using the above laser configuration, we calculate the final energies and the emission angles of protons, which were initially at rest, by independently integrating Eq. (1) by the Adams method.53 With axially symmetric laser beams, the final distribution of accelerated particles also does not depend on the polar angle in the xy plane ϕ, measured from the direction of laser polarization.29 Note that the use of axially symmetric beams in this example facilitates the analysis, although it is not demanded by the method itself, which, with certain modifications of the training data set, namely, considering more general asymmetric cases, can be used to assess more complex beams with no axial symmetry. With polar-angle symmetry, we consider θE distributions of particles, where E is the final proton energy and θ is the azimuth angle between the particle emission direction and the z axis.

The laser intensity distribution in the focal spot, in both the transverse and longitudinal directions, has a direct impact on the angular and energy distributions of protons accelerated from the focal spot region. To demonstrate this, we consider, for example, two beams with Gaussian and Laguerre–Gaussian initial transverse intensity profiles, focused into a spot with the same FWHM size D_F = 1.0λ and peak intensity 10²² W cm⁻². Note that for the same peak intensity, the two considered profiles have different total energies: for the Laguerre–Gaussian beam, the total energy is ∼4 times greater than that of the Gaussian beam owing to differences in their intensity distributions. The transverse intensity profiles of the two beams in the focal plane are shown in Figs. 3(a1) and 3(a2). In this plane, the main and most noticeable difference between the profiles is the presence of a secondary maximum for the Laguerre–Gaussian profile in the form of a radial ring that extends from ∼1λ to ∼2λ. Within this ring, the differences between the two profiles are less pronounced. Moreover, their longitudinal distributions also make it possible to distinguish the Gaussian and Laguerre–Gaussian distributions from each other. These longitudinal cross-sections are presented in Figs. 3(b1) and 3(b2), without any time dependence being imposed on them. As can be seen, the Laguerre–Gaussian beam is much more prolonged along the propagation direction z in comparison with the Gaussian beam of the same diameter D_F = λ. The Rayleigh length of the former is about four times greater than that of the latter. This difference is reflected in the angular energy spectra of protons. For axial symmetric laser pulses, proton spectra do not depend29 on the polar angle ϕ and can be integrated over it. In a more general case, the data for each ϕ have to be preserved. In this work, we focus on assessment of the radial intensity profile from the angular energy distribution of particles, and we integrate over ϕ. The resulting spectra for the considered Gaussian and Laguerre–Gaussian beams are shown in Figs. 3(c1) and 3(c2). Qualitatively, both particle distributions are similar in the sense that the protons of maximum energy are accelerated in the direction perpendicular to the laser axis θ = 90°; the more they depart from this direction, the less energetic the particles become. The cutoff energy, i.e., the maximum energy observed for particles at θ = 90°, is about 18 keV in both cases. However, the characteristic width of the angular energy distribution at half the maximum energy differs qualitatively between the two cases. For the Gaussian beam, the characteristic width is ∼4 times greater. These results can be interpreted using the simplest analytical expressions. The maximum energy of particles, E_max, is observed along the direction of the highest intensity gradient, i.e., along the beam radius r, and, for an elementary estimate, the value of this cutoff can be considered29 to be proportional to the square of the intensity gradient in this direction:$$E_{\mathit{\text{max}}}\propto {\left({\left.\frac{\partial I(r,z)}{\partial r}\right|}_{r_0,z_0}\right)}^2,$$(5)where r₀ and z₀ are the coordinates of the position at which the maximum gradient is observed with the given beam spatial profile. Since the intensity variation is steepest in the inner area of the Laguerre–Gaussian beam, where it closely matches the Gaussian beam, and both beams have the same diameter, protons with the highest energy accelerate from the same region under the action of intensity gradients of approximately the same strength. Thus, their maximum energies almost coincide, which can be observed in Figs. 3(c1) and 3(c2). When studying the same distributions for beams with no axial symmetry, the value of the cutoff energy in each direction ϕ will be different, and so if the input is treated as a three-dimensional array with ϕ as the third variable, it is possible to estimate the maximum gradient and consequently the slope of the beam profile in this particular direction. The angular width of the spectrum, Δθ, in turn, depends on the gradient along the propagation direction z, which should be greater for a Gaussian beam owing to the smaller Rayleigh length in comparison with that of the Laguerre–Gaussian beam. According to an equation from Ref. 29, the angular width of the proton spectrum relates to the parameters of the beam as follows:$$\tan (\mathrm{\Delta }\theta /2)\propto \frac{D_F}{z_R}.$$(6)Substituting the corresponding beam parameters, namely, D_F = λ and z_R = 4.3λ for the Gaussian beam and D_F = λ and z_R = 17.8λ for the Laguerre–Gaussian beam, into this formula yields a difference by a factor of ∼4 in the angular width of their proton spectra, which is consistent with the results presented in Figs. 3(c1) and 3(c2).

The presence of a ring in the Laguerre–Gaussian beam also affects the proton spectrum, which is more evident if one considers ϕ-integrated spectra with a logarithmic scale on the energy axis, as shown in Figs. 4(b1) and 4(b2), in particular their low-energy parts, outlined with red rectangles. These distributions are shown together with plots of the squared intensity gradients in the xz plane in Figs. 4(a1) and 4(a2), from which it can be seen that as well as a more prolonged profile along the beam axis z, the Laguerre–Gaussian distribution has two additional areas where the gradient may reach ∼10⁻¹ of its maximum value. These correspond to the intensity slopes on the inner and outer edges of the ring. Protons accelerated from these regions form the low-energy part of the angular energy distribution; see Figs. S1 and S2 in the supplementary material for better insight into this. For the Laguerre–Gaussian beam, owing to the spatial arrangement of gradients in this region, the low-energy part of the proton spectrum, i.e., the part with energies 10⁻⁴ ≲ E/E_max ≲ 10⁻², has a peculiar shape with several pronounced branching wings [see Fig. 4(b2)]. By contrast, for the Gaussian profile, this part of the angular energy distribution contains significantly fewer particles and has a much simpler structure, owing to the less complex structure of the beam itself. These differences can also be noticed in the total θ, ϕ-integrated energy spectra of protons. These distributions are presented in Fig. 5(a). Here, for the Laguerre–Gaussian profile, a substantially greater number of protons accumulate in the low-energy region, owing to the large volume of the secondary ring region, while high-energy parts of the spectra look almost identical, since they are formed by protons accelerated from the central region with steep intensity changes, where the two profiles closely resemble one another.

As can be seen in Fig. 5(a), the distinct cutoff energies are almost the same for these two beam profiles. The general dependence of the cutoff energy on the maximum intensity gradient is a useful characteristic that, on the one hand, allows the validity of Eq. (5) to be checked and, on the other hand, allows the analysis to be simplified by making it possible to retrieve the maximum gradient and peak intensity for a given profile from the measured maximum energy of protons and vice versa. The curves characterizing this dependence were obtained in a set of simulations with three different beam profiles (flat-top, Gaussian, and Laguerre–Gaussian), beam sizes in the range ∼(1–4)λ, and two different maximum laser intensities I₀ = 10²² W cm⁻² and I₀ = 10²³ W cm⁻². These curves are shown for the two intensities in Figs. 5(b1) and 5(b2), respectively. At lower intensities, they appear to closely follow the dependence predicted by Eq. (5), while at higher intensities, nonlinear effects start to manifest themselves. The greatest deviation from Eq. (5) is observed in the case of the most tightly focused beam with a diameter of ∼λ and a peak intensity 10²³ W cm⁻² and constitutes about 16% [see Fig. 5(b2)]. The lower values for the real dependence in comparison with those predicted by Eq. (5) can be explained by the fact that in a sufficiently strong intensity gradient, a proton can attain energies high enough to cause significant shifts in its position during the time of its interaction with the laser pulse, which is normally much shorter than the time necessary for the proton to move from its initial position by a notable distance. As a result, the proton interacts with maximum gradient for only a fraction of the whole time of the interaction and attains a lower final energy. Thus, this analysis indicates that the simplified scaling given by Eq. (5) is applicable up to 10²³ W cm⁻². In this range, if the shape of the intensity profile remains the same, i.e., $I_2(\vec{r})=k{I}_1(\vec{r})$, then the final kinetic energy of the proton can be immediately found as E_i2 = k²E_i1, eliminating the need for additional simulations. Conversely, given two values of the cutoff energy E_max1 and E_max2, it is possible to retrieve the unknown scale of the second intensity profile with the same relative slope as $I_2(\vec{r})=\sqrt{E_{\text{max}2}/E_{\text{max}1}}{I}_1(\vec{r})$, which is extremely useful for the reconstruction task considered here.

For the development of a robust technique for diagnosing tightly focused laser pulses with extreme intensity based on the properties of directly accelerated proton spectra, it is necessary to find an effective solution of a multiparameter inverse problem. As the relations between intensity and particle distributions are nontrivial, assessment of experimental data requires an approach that takes this complex and nonlinear behavior into account. Here, to reconstruct the spatial profile of the laser pulse and estimate its peak intensity from the angular energy spectrum of accelerated protons, a method based on an ANN is proposed and analyzed. ANNs constitute a vast subset of machine learning techniques and are based on an algorithm that is “trained” in advance on labeled data sets. For a detailed explanation of the basic principles of ANNs, the reader can consult Refs. 54 and 55, and here we just provide a brief description of the dense feedforward neural network relevant in the context of this work. The training data set contains both the input and output variables—explicitly in the case of the former. The ANN is trained to model output variables for a given input vector. This vector is passed through one or more intermediate layers of the ANN, called hidden layers. Each of these generally contains several nodes, called hidden units, or “neurons.” A neuron takes an input vector, applies to it adjustable weight and bias vectors, and passes it through a nonlinear activation function. The output vector is compiled from linear combinations of all the resultant outputs of each hidden node and passed on in the forward direction. In the final hidden layer, the hidden units are connected elementwise to the nodes of the output vector. The procedure of training involves optimizing the weights and biases of each hidden node to provide the best fit to the output data available for the training set. The “best fit” metric is chosen for a particular task by the user. One of the common metrics for regression tasks, where the input variable is mapped to a continuous output variable (as in our case), is the mean squared residual.

The input of an ANN can have an arbitrary shape, in particular it can be a two-dimensional map, or, colloquially speaking, an image. In this case, using a regular dense feedforward neural network may not always be optimal, since groups of pixels may form distinctive features, making their relative positions important. If the image is converted into a one-dimensional vector, these features are ignored, since the information about mutual positions of pixels is lost. This problem can be avoided by using a CNN architecture, which introduces one or more convolutional layers before the first hidden layer. In each of them, an image is convolved with several optimizable kernels extracting different characteristic features. Convolutional layers may also be followed by so-called “pooling” layers, where the output of the convolutional layer is downsampled on a rougher grid by, for example, taking local maxima in nodes of the old grid as the values in each node of the new grid. This procedure reduces the size of the data with which the ANN works and thus reduces the number of parameters for optimization, improving the ANN training. It is interesting to note that the input image of a CNN may also contain multiple color channels, making it a three-dimensional array rather than a two-dimensional image. This can be used to factor in some extra information or include an additional variable, which in the context of this work can be employed for assessment of axially asymmetric beams that have an additional dependence on the angle ϕ. This work aims to show that a CNN combined with one fully connected layer is capable of retrieving radial focal distributions of the laser pulse from the angular energy distribution of protons, which represents a relevant task for the development of intense laser pulses diagnostics.

To reconstruct the axially symmetric focal distributions of laser pulses from the angular energy distributions of protons, we chose a CNN architecture with two convolutional layers, two pooling layers, and one fully connected (dense) layer [Fig. 6(a)]. Grayscale 8-bit images of the angular energy spectra with dimensions 100 × 100 serve as our input. These are similar to the distributions presented in Figs. 3(a1) and 3(a2), but the energy axis is normalized and converted into a logarithmic scale, since the latter appears to be informative in the context of the mandated task. The normalization is done such that the top row of pixels of the image corresponds to the cutoff energy E_max, while the bottom row corresponds to an energy that is four orders of magnitude lower than E_max. The initial absolute value of the cutoff energy is saved and treated separately to obtain the scale on the intensity range in accordance with the scaling law described above. The θ axis uses a linear scale such that the left column of the image corresponds to θ = 0° and the right column to θ = 180°. The intensity of each pixel is set on the basis of proton density in a particular bin in a logarithmic scale such that the maximum number of protons always corresponds to 255. For the intensity focal distributions, 200 × 200 8-bit images serve as the output. Note that in this work, for a proof-of-principle demonstration, we consider only axially symmetric beams, which means that the output can be reduced from a two-dimensional image to a one-dimensional vector. However, aiming at subsequent realization of a more general situation in which the output profile may be axially asymmetric, the two-dimensional shape of the retrieved intensity distributions is retained. The x and y axes have linear scales and cover the range from −7.5λ to 7.5λ. The color scale for intensity images, unlike that for particle distributions, is linear, since this has proved to be more robust for the reconstruction task. It is normalized such that the maximum intensity always corresponds to 255. The absolute value of the peak intensity is saved and treated separately in conjunction with E_max and the scaling discussed above. First, input images are passed through two convolutional layers [“Conv2D” in Fig. 6(a)] with the rectified linear unit (ReLU) f(x) = max(0, x). As a result of convolution with optimizable kernels, 32 and 64 feature maps are obtained after the first and second convolutional layers, respectively. In addition, images after each convolutional layer are downsampled using pooling layers [“MaxPool2D” in Fig. 6(a)]; in this case, max pooling is used, i.e., an image is divided into small regions, and local maxima in each region are taken for the new pixel value. This allows small local groups of similar features to be merged into one and significantly reduces the number of optimized parameters; also, it enhances CNN training. The output of the final pooling layer with 64 different feature maps is converted into a one-dimensional vector and submitted to the hidden layer, which is a fully connected (dense) layer that connects each of the 1600 nodes of this input vector to each of the 40 000 output vector nodes. Here, the sigmoid function $\sigma (x)={(1+e^{-x})}^{-1}$ is used as a nonlinear activation function. The adjustable kernels in the first and second convolutional layers have sizes of (5, 5) and (3, 3), respectively. These match the sizes of the first and second pooling kernels, which are also chosen as (5, 5) for the first pooling layer and (3, 3) for the second pooling layer. The dense layer consists of 100 neurons. The hyperparameters of the neural network, such as the number of layers of each type and the number of nodes in each layer, are adjusted on the basis of a grid-search algorithm.56 For future studies, the proposed architecture can be easily modified for assessment of axially asymmetric beams that depend on the additional variable ϕ by changing the dimensionality of the input from two to three, i.e., by moving from grayscale single-channel images to colored multichannel ones. In addition, for axially symmetric beams, the longitudinal profile may be retrieved by introducing an additional dimension in the plane containing the beam axis.

To train the CNN, a data set consisting of 1030 pairs of intensity distributions and their corresponding angular energy distributions of protons was created. The data were roughly evenly distributed among five different possible initial intensity profiles: flat-top, Gaussian of the second, fourth, and sixth orders, and Laguerre–Gaussian. For each profile, ∼20 different f-number values f_# = F/ρ were considered, yielding focal spots with sizes in the range ∼(1–4)λ. To make the CNN more robust to the quality of input data, the number of particles for each parameter combination alternated between 10³ and 10⁴ for all profiles except Laguerre–Gaussian. For the latter profile, the number of particles was varied between 2.1 × 10⁴ and 7 × 10⁴. The higher overall number of particles for the Laguerre–Gaussian beam was used to ensure that the density of the proton cloud was approximately the same in all cases. Also, since the Laguerre–Gaussian beam is significantly wider, taking into account the ring and large Rayleigh length, the total number of particles in this case had to be increased.

The mean squared error was chosen as the loss function for optimization, while the optimization itself was performed with a first-order gradient-based method using the Adam algorithm.57 To avoid overfitting, which is a common problem resulting from the high flexibility of ANNs,58 the training continued until the loss function for the validation subset of the data, which is not explicitly used for weight adjustment, did not decrease for more than 20 epochs. At this point, the training was stopped, and the best weights corresponding to the lowest loss for the validation subset were saved for subsequent analysis of the intensity distributions. The data were divided into training and validation subsets in the proportion 9:1. To obtain statistically significant estimates, a total of ten training runs were performed. This was done in a k-fold cross-validation setup55 with k = 10, which implies that for each training run a different 10% portion of the data was taken for validation. An example of the learning curves obtained in one of these training runs —one showing the decrease in the root-mean-square error with the number of iterations for the training and validation subsets—is presented in Fig. 6(b). As can be seen, after a rapid descent in the initial training stage, the curves at some point start to level out, especially for the validation subset, implying that subsequent training does not result in a significant decrease in the prediction error. It can also be seen that the final error for the training subset is somewhat lower than the error for the validation data. This can be attributed to the fact that the validation subset consists of examples that were not used for the ANN optimization. Thus, they are “new” to the ANN and are used to test its ability to generalize and make adequate predictions on the previously unseen data. Therefore, slightly higher errors are to be expected for the validation subset.

The trained CNN performed reasonably well on the validation data set, with an overall profile similar to its a priori known correct counterpart and with prediction errors for the characteristic parameters such as the beam diameter not exceeding a few percent. In the final test, the validity of its predictions was verified on a completely new small data set consisting of nine images. The images were obtained for the following profiles:•flat-top profiles with {f_# = 1.3, I₀ = 2 × 10²² W cm⁻²}, {f_# = 2.8, I₀ = 3 × 10²¹ W cm⁻²}, and {f_# = 3.3, I₀ = 6 × 10²¹ W cm⁻²};•Gaussian profiles with {f_# = 1.4, I₀ = 8 × 10²¹ W cm⁻²}, {f_# = 2.2, I₀ = 1.3 × 10²² W cm⁻²}, and {f_# = 2.6, I₀ = 4 × 10²¹ W cm⁻²};•Laguerre–Gaussian profiles with {f_# = 2.1, I₀ = 3 × 10²¹ W cm⁻²}, {f_# = 2.9, I₀ = 1.6 × 10²² W cm⁻²}, and {f_# = 3.5, I₀ = 3 × 10²² W cm⁻²}.

These profiles, both correct and predicted by the CNN, are shown in Fig. 7. For the plots presented here, the maximum number of particles in the proton angular energy spectra was used, i.e., N_tot = 10⁴ for the flat-top and Gaussian profiles and N_tot = 7 × 10⁴ for the Laguerre–Gaussian profile. The results are summarized in Table I, including the correct values of the beam diameter D_F and the maximum intensity I₀, the values predicted for the same parameters by the CNN, their relative errors, and the root-mean-square errors calculated for the whole two-dimensional profile. It appears that the CNN has achieved the ability to almost perfectly reproduce the correct intensity profile from the angular energy distribution of protons. In some cases, the predictions made by the CNN are nearly indistinguishable from the correct distributions, while the relative errors for both the beam diameter D_F and the maximum intensity I₀ do not exceed ∼6% for all the presented cases. One-dimensional cuts are shown for ease of comparison. As an example, two two-dimensional profiles, the “true” one for laser beam parameters in Fig. 7(b2) and the corresponding one retrieved by the CNN, are presented in Figs. 7(d1) and 7(d2). The retrieved profile is also axially symmetric and closely matches the correct one. It has been verified that the same applies to the rest of the retrieved two-dimensional distributions. Thus, we can conclude that the CNN-based approach developed here is viable for reconstructing axially symmetric intensity distributions in the focal spot region from the angular energy distributions of accelerated protons.

In real experiments, the number of particles in each bin of the angular energy spectra may be different and depends on the density of the plasma cloud, the sensitivity of detectors, and other factors. Thus, it is important to ensure that the developed approach is reasonably insensitive to these variations. For this purpose, the CNN was a priori trained on angular energy distributions with different numbers of particles in them, ranging between 10³ and 10⁴ for all profiles except Laguerre–Gaussian, for which the number of particles in the training data was varied between 2.1 × 10⁴ and 7 × 10⁴ for the reasons outlined above. The predictions shown in Fig. 7 were made using particle distributions with the maximum number of particles in them. To understand how a lower number of particles affects accuracy and what statistics are sufficient for the CNN to make adequate predictions, the number of particles was gradually decreased, and the obtained errors were compared with those presented in the first nine rows of Table I. The lower limit at which the accuracy of the CNN predictions appears to be independent of the particle statistics is about three times lower than the maximum value, i.e., 3 × 10³ for the flat-top and Gaussian profiles and 2.1 × 10⁴ for the Laguerre–Gaussian profile. Examples of angular energy profiles formed by 100% and 30% of the maximum number of particles in the simulation are shown for the flat-top intensity profile in Figs. 8(a1) and 8(a2). As can be seen, if only one-third of the total number of particles participate in forming the angular energy spectrum, it becomes quite discrete, clearly lacking statistics in certain regions. Nevertheless, as the analysis conducted here shows, it is quite sufficient to make accurate predictions. Predicted radial intensity profiles for the considered angular-spectral distributions are shown in Figs. 8(b1) and 8(b2). The CNN seems to make correct predictions, which has also been verified for other test profiles, and the results have been included in Table I: see rows 10–18. In certain cases, it was found that the number of particles can even be decreased by a factor of ∼10; however, for some profiles, the accuracy decreases, and the errors for the beam diameter and the maximum intensity may increase by more than 10%–20%.

Considering the limitation on the maximum density of protons introduced at the beginning of Sec. II, it is worth estimating how this fits with the obtained requirements on the number of particles for sufficient statistics. As was shown in the previous paragraph, ≳3 × 10³ protons are required for correct retrieval of the intensity profile if it does not have wide secondary rings (e.g., a flat-top or Gaussian profile). In simulations, we average the results over the polar angle ϕ. In real-life experiments, it will most likely be impossible to position detectors all around the target in a spherical formation. The ring-like arrangement, on the other hand, is more plausible. In this case, if the characteristic size of the detector is Δϕ ∼ 1°, then the total number of particles has to be higher by a factor of 360°/Δϕ = 360, i.e., there must be about 10⁶ particles, to obtain the same statistics as in simulations. In addition, detector efficiency has to be taken into account. For multichannel plate detectors, which are among the possible types of detectors that can be used for registering relatively low-energy protons, the absolute detection efficiency exceeds 10⁻³ for protons with energies greater than 10 eV.59 Thus, the number of particles in a focal volume should be greater by an additional factor of ∼10³. The resulting number of 10⁹ particles are distributed in a volume of ∼(100)³ µm³, yielding a concentration of ∼10¹⁵ cm⁻³. This value is consistent with the limitations discussed at the beginning of Sec. II (10¹⁶ cm⁻³ for a 1 mm thick jet). Furthermore, it is worth mentioning that some high-power laser facilities offer high repetition rates, enabling tens of shots to be performed per hour; see, for example, Ref. 60. In this case, if needed, the statistics of registered particles may be bolstered by accumulating data from multiple shots, provided that the shot-to-shot stability of the laser pulse energy is sufficiently high. In addition, they may also be improved by positioning more detectors around the target, provided that the experimental setup allows this.

In addition to particle statistics, it is important to verify the effect of detector misalignment on the predicted value. For example, in an experiment, the positions of the detectors corresponding to θ = 90° and other directions may be set incorrectly, leading to some distortion of the angular energy proton spectrum in comparison with the model distribution obtained in simulations. To find the range where the predictions are insensitive to this misalignment, particle distributions from the test set were artificially modified by rotation through a small arbitrary angle. The resultant intensity profiles reconstructed by the CNN from the modified particle distributions were compared with the correct ones. It was found that the optimal range of the misalignment angle is about ±1°. Beyond this range, the particle distributions become significantly distorted and are misinterpreted by the CNN, with the errors for D_F and I₀ increasing above ∼10%–20%.

Thus, we can conclude that the total number of particles and the misalignment angle should be viewed as important parameters that have to be considered at the stage of preparing the experiment to ensure that intensity profiles are retrieved with good accuracy.

As with other neural network-based algorithms, our approach is data-driven, and its applicability for reconstruction of a particular focal distribution depends on the extensiveness and diversity of the training data set. If the spatial profile of the laser pulse is qualitatively similar in appearance to those on which the neural network was trained, the reconstruction is expected to work correctly. However, this approach may fail in a completely new case that was not represented in the training data.

Consider an example of an elliptic beam, as shown in Fig. 9(a). In simulations, this was obtained by tight focusing of a Gaussian beam with an off-axis mirror, half of which was covered by an aperture in the area of x < 0 (before reflection). This leads to an elongation of the intensity profile along the x axis (ϕ = 0°). The ϕ-integrated angular energy distribution of protons accelerated by this laser beam is shown in Fig. 9(b). Since the beam has different widths in different directions ϕ, the angular energy profiles for these directions are also different, as can be seen in Figs. 9(c) and 9(d), where spectral distributions of protons along the major and minor ellipse axes are presented. One of the key distinctions between these is the value of the cutoff energy. Since the beam has a greater width and a lower gradient along ϕ = 0°, protons accelerated along this axis gain less energy than those accelerated along ϕ = 90°, and the cutoff for the former is a few times lower. Passing the ϕ-integrated spectral distribution through the CNN, which was trained purely on axially symmetric data, yields an axially symmetric reconstructed profile that has the same widths along the x- and y-axis directions (ϕ = 0° and ϕ = 90°, respectively), in contrast to the spatial distribution of the original laser beam with different widths [see Figs. 10(a) and 10(b)]. Interestingly, the width of the reconstructed axially symmetric profile is closer to the width of the original profile for ϕ = 90°, i.e., in the direction of the ellipse’s minor axis, which is least affected by the introduction of the aperture.

It can also be noted that the predicted profile has a pronounced peripheral ring that was not observed in the original laser beam. The dependence of the beam width on the angle ϕ leads to a redistribution of the ϕ-integrated spectrum of accelerated protons compared with the results obtained for the axially symmetrical Gaussian beam. The neural network was trained on data in which changes in the form of the angular energy spectra, namely, changes in the proportions of high-energy and low-energy particles, were the result of the appearance of peripheral rings, such as the Laguerre–Gauss mode ring. Thus, in the case of an elliptical laser beam, the neural network works in exactly the same way, predicting peripheral rings to describe the cause of the proton redistribution. Although separate treatment of spectral distributions for each angle ϕ for retrieval of the intensity profile in this direction may seem an apparent solution, this is also incapable of correct reconstruction, as can be seen from Figs. 10(c) and 10(d), where predictions made on the basis of ϕ = (0.0 ± 1.0)° and ϕ = (90.0 ± 1.0)° segments of the angular energy proton spectrum are presented. In this case, for the minor axis [ϕ = (90.0 ± 1.0)°], the beam width and intensity are reproduced with $\approx 24$% error, while the difference between the retrieved width and intensity constitutes a factor of ∼2 for the major axis [ϕ = (0.0 ± 1.0)°]. However, the characteristic gradient value is reproduced correctly in both cases, which is evident from a comparison of the slopes of the original and reconstructed profiles. The great inconsistency of the results obtained for ϕ = (0.0 ± 1.0)° can be explained by the impact of the longitudinal distribution of the laser beam on the proton angular spectra. Previously, the beam diameter and its Rayleigh length changed reciprocally for a given profile, and their relation determined the angular width of the spectral distribution of protons according to Eq. (6). For the beam partially covered by the aperture, the Rayleigh length stays almost the same, while the size of the beam along one of the dimensions increases. In accordance with Eq. (6), this leads to a significant broadening of the spectrum, which can be clearly seen in Fig. 9(c). In the training data set, such a great angular width was only obtained for a very tightly focused beam with a diameter of one wavelength; see, for example, the angular energy spectrum in Fig. 3(c1). Thus, the implemented neural network is biased toward retrieving smaller focal spot sizes given a spectral distribution with large angular width, which leads to the observed underestimation of the beam size in the x-axis direction.

Based on the analysis conducted here, we come to the familiar conclusion that the neural network may perform inaccurately for an input on which it was not trained. Thus, to assess more complex beam profiles, a sufficient amount of relevant data has to be included in the training data set. For example, if we expect elliptic beams in a real experiment owing to the use of one-half of an off-axis mirror, then the training set has to be chosen accordingly, so that the neural network can model this type of beam. Note that in this case, a transition to three-dimensional input data will be required, since the cutoff energies and the general shape of the two-dimensional angular energy spectrum depend on the variable ϕ. With necessary modifications and extension of the training data set, the approach can be applied to treat asymmetric and noncircular beam profiles. In some cases, some lower symmetry of the beam may probably be used for reduction of the necessary number of measurements at different angles ϕ. For example, square beam profiles have two lines of symmetry along the sides of the square and two diagonal lines of symmetry, and so full measurement of their angular energy distribution in the range Δϕ = 2π is not absolutely necessary for their characterization. Some information about possible symmetries may be obtained on a first step of rough estimation of focal properties with a reduced intensity.

It is also worth noting that the focal intensity distribution may be affected by wavefront phase distortions. If these distortions are significant, they may impede the correct retrieval of the intensity profile. To avoid this problem, a more diverse and broad training set with the inclusion of complex intensity profiles resulting from an imperfect laser beam wavefront would be required. In this regard, preliminary wavefront estimation with other techniques such as Shack–Hartmann sensors61 at reduced laser beam energy may provide useful information about possible wavefront distortions and their effect on the intensity profile.

The development of the high-energy laser facilities with laser beam power exceeding the PW level requires diagnostic approaches that can characterize the intensity distribution in the focal spot region and the peak intensity without the need for a reduction in laser power. In this work, a promising diagnostic technique for tightly focused laser pulses with extreme intensities has been discussed. It is based on measuring the angular energy spectrum of protons accelerated from the focal spot region. As has been shown here, the parameters of the particle distribution depend on the focal distribution of intensity, although in a nontrivial manner, which consequently complicates the analysis of the obtained data. Although it has been demonstrated that some parameters, such as the cutoff energy and the angular width of the particle energy angular spectrum, are defined by the maximum intensity and diameter of the focused beam, as well as its Rayleigh length, the whole particle distribution has to be considered to obtain a complete reconstruction of the intensity profile. So, for instance, its low-energy part may contain information about the presence of secondary rings with relatively low-intensity slope, making this part of the particle spectrum essential for distinguishing, for example, between a tightly focused beam with a simple Gaussian profile and a beam with sufficiently high-intensity slopes on the periphery, such as for the Laguerre–Gaussian profile considered here.

In this study, we have shown that by using a CNN, it is possible to reconstruct the radial distribution of intensity of a tightly focused high-energy laser beam in the focal plane with acceptable accuracy and with measurement errors of the main beam parameters that do not exceed a few percent, assuming that necessary conditions on the number of accelerated particles and the misalignment angle of the detectors are satisfied. Implementation of the proposed diagnostic involves the following:(i)positioning an array of proton detectors able to perform energy-resolved measurements around the focal spot region at equal distances from it, so that each of the detectors corresponds to a particular angle θ and the whole array of them encompasses a sufficiently large area around θ = 90°;(ii)obtaining a set of measurements of proton spectra in each direction θ;(iii)preprocessing the obtained data and integrating the one-dimensional energy spectrum from each detector into one two-dimensional angular energy distribution;(iv)passing the experimentally resulting angular energy spectrum through the trained CNN to obtain the intensity distribution in the focal spot region.

Although the current scheme has assumed axially symmetric beams, it is not limited to the study of such beams. In the absence of axial symmetry, additional arrays of detectors have to be placed around the target, so that each subset of detectors corresponds to a different ϕ. The data for each ϕ can then be treated separately or by a CNN with three-dimensional input to come up with the radial intensity profile for a particular direction ϕ. Moreover, the comparison of the proton spectra along different angles ϕ allows a test of the hypothesis of laser beam axial symmetry at the beginning of laser profile reconstruction.

The described method allows direct implementation of the focal distribution diagnostics in state-of-the-art and forthcoming facilities such as ELI, APPOLON, and XCELS. Preliminary measurements and expectations of the focused laser beam profile may facilitate the application of the method, but although the uncertainties of the former do not generally cause problems, more diverse training sets will be required to account for the more complex focal distributions resulting from imperfect laser beam wavefronts and other factors expected in real experiments.

For the further development of the proposed approach, it can also be adapted to angular energy distributions of particles with charge-to-mass ratio different from that of protons, i.e., heavier ions or electrons. The former would allow the approach to be extended without significant modifications to a wider intensity range, since heavy ions are expected to reach lower energies than protons, which justifies the use of the same nonrelativistic ponderomotive force approximation given by Eq. (1) and a simple linear scaling between the cutoff energy and maximum squared intensity gradient for intensities exceeding 10²³–10²⁴ W cm⁻². The use of electrons, in turn, may supply additional information about tightly focused laser pulses, since despite the great difference in their charge-to-mass ratio in comparison with protons and their much more efficient acceleration in the laser field, their angular energy distributions also seem to exhibit some relation to the laser spatial profile, as has already been shown in previous work.28 Finally, the approach presented here can be extended to deal simultaneously with several kinds of particles with different dynamics (protons, electrons, and possibly heavy ions), which can be used to estimate a wider range of laser pulse parameters, including duration.30

Acknowledgment. This work was supported by RFBR, ROSATOM (Grant No. 20-21-00023), and the Ministry of Science and Higher Education of the Russian Federation (Agreement No. 075-15-2021-1361). O.E.V. acknowledges the Foundation for the Advancement of Theoretical Physics and Mathematics (“BASIS”) for financial support (Grant No. 22-1-3-28-1). Numerical simulations were partially supported by resources of the NRNU MEPhI High-Performance Computing Center.