Superintense Laser-driven Ion Beam Analysis

Ion beam analysis techniques are among the most powerful tools for advanced materials characterization. Despite their growing relevance in a widening number of fields, most ion beam analysis facilities still rely on the oldest accelerator technologies, with severe limitations in terms of portability and flexibility. In this work we thoroughly address the potential of superintense laser-driven proton sources for this application. We develop a complete analytical and numerical framework suitable to describe laser-driven ion beam analysis, exemplifying the approach for Proton Induced X-ray/Gamma-ray emission, a technique of widespread interest. This allows us to propose a realistic design for a compact, versatile ion beam analysis facility based on this novel concept. These results can pave the way for ground-breaking developments in the field of hadron-based advanced materials characterization.

www.nature.com/scientificreports www.nature.com/scientificreports/ PIXE analysis can be applied by ideally splitting the sample in a finite number J of fictitious layers having mass thicknesses ρ j R j and homogeneous compositions W i,j , as shown in Fig. 2(a). Assuming to perform K measurements at different proton energies E p k , the following system of equations can be written 22,23 :  refer to the k-th proton energy spectrum. Equation (4) implies that also in the case of broad spectrum sources it is necessary to use at least K different spectra to solve the system.
The formalism developed here for PIXE with arbitrary ion sources extends naturally also to PIGE 24 : the ionization cross section σ i (E) and the fluorescence yield ω i must be replaced with the proper nuclear reaction cross section and isotopic abundance, while γ-ray attenuation (i.e. the exponential terms in all the equations) can be neglected.
If a laser-driven proton source is used, the shaping function f p (E p ) can take different forms, depending on the specific ion acceleration process 10 . Target Normal Sheath Acceleration (TNSA) is arguably the most well understood laser-driven ion acceleration mechanism: it stands out for its robustness among other schemes 25 , it has been used extensively as diagnostic tool in laser plasma interaction experiments 26 and has been demonstrated in a repetitive regime 27 . TNSA can be outlined as follows: a micrometric solid foil is irradiated with an ultra intense laser pulse, which accelerates the electrons of the target towards the back side. The expansion of the hot electron cloud generates then an intense electric field, which drives ion acceleration. Accelerated ions (mostly protons from the hydrocarbon contaminants at the back side of the target) show an exponential energy distribution, with cut-off energies ranging from few MeVs for multi-TW class lasers 13,14 up to ~100 MeV 28,29 for high energy large scale laser systems. The energy distribution can also depend significantly on the angular selection as well as on the subsequent propagation toward the sample region. Other non-TNSA acceleration schemes, such as Collisionless Shock Acceleration (CSA) 30 or Radiation Pressure Acceleration (RPA) [31][32][33] have been reported in the literature. These acceleration schemes can provide a narrower energy spectrum with respect to TNSA, though still much broader than conventional electrostatic accelerators, and thus still incompatible with the standard algorithms for PIXE with monochromatic sources. Therefore, CSA or RPA cannot be considered more suitable than TNSA for PIXE. The approach described here could nonetheless be used in principle also for CSA or RPA.
Quantitative laser-driven PIXE analysis, both considering homogeneous samples and differential PIXE, requires to implement an iterative algorithm to solve the system of equations (3) or (4) to retrieve the sample composition, analogous to those used for PIXE with monoenergetic protons. We will test the capability to retrieve material properties in unknown samples by performing complete simulations of laser-driven PIXE measurements, as detailed in the following. Reliable and complete experimental datasets including both energy spectra and angular distribution for several laser intensities are not available in the literature. Therefore, we will exploit both a simplified analytical modeling and a realistic numerical description of the laser-driven proton source.

Simulation of a laser-driven PIXE experiment
The goal of the present work is to assess the potential of laser-driven PIXE in realistic scenarios, i.e. demonstrate the possibility to retrieve the composition of complex, realistic samples from X-ray emission induced by a laser-driven proton source. To this aim, we consider a selection of artifacts directly relevant for cultural heritage studies (e.g. paintings, jewelry…), a research field in which PIXE finds widespread use. The complexity of these artifacts, which naturally have a multi-elemental multi-layer composition, is an ideal testbed to demonstrate the general feasibility of laser-driven PIXE and differential PIXE for a wide range of applications. For the purpose outlined here, no useful experimental data are available in the literature, with only one, partial, exception. In 34 the possibility to obtain an X-ray spectrum irradiating ceramic and metallic samples with a laser-driven proton source was demonstrated, assessing also the absence of any damage to the irradiated samples. However, the full composition was not retrieved and the reported data are too simplified to be considered for the goals of this work.
Therefore, by means of Geant4 Monte Carlo simulations, we generate "synthetic" realistic X-ray spectra (i.e. the experimental X-ray yields ⁎ Y i or ⁎ Y i k ) in a variety of laser-driven PIXE scenarios (different artifacts and different laser-driven ion sources). For each scenario, the x-ray spectra are analyzed with an iterative code based on equations (3) and (4), as well as on an efficient minimization algorithm (see Methods section) designed to retrieve the sample compositions. These compositions are then compared with the sample composition initially set in the Monte Carlo simulations. In the following we will refer to the elemental concentrations set in the Monte Carlo code as the real concentrations and to the elemental concentrations retrieved by the iterative code as the retrieved concentrations. Specifically, we investigate a total of three scenarios (involving two different experimental set-ups).
For the first scenario we consider PIXE analysis of a homogeneous metallic sample in vacuum, performed with an idealized proton source. In this case, we choose a simple analytical TNSA-like exponential for the energy spectrum of the proton source in the Monte Carlo: p p p p , characterized by a maximum cut-off energy E p,max and a parameter T p (related to the mean proton energy). We choose E p,max = 5MeV and T p = 0.7MeV, which is realistic for a proton source driven by a 20 TW laser (see e.g. ref. 14 , where TNSA is performed with a pure hydrogen target). The simple model adopted for the ion source allows also to perform a sensitivity analysis to address the effects of typical fluctuations of both E p,max and T p on a PIXE measurement. Figure 3(a) shows these proton spectra. We consider a sample composed by iron, copper, zinc, tin and lead, which is representative of a Roman sword-scabbard 35 .
For the second scenario we consider differential PIXE analysis of a non-homogeneous metallic sample in vacuum performed with an idealized proton source. In this case we use 6 different analytical proton spectra with properties (E p,max and T p ) compatible with a laser system having a maximum power of 40 TW 14 (highest E p,max equal to 6 MeV). Figure 4(f) shows these spectra at the sample surface. In this scenario, the sample composition is inspired by the elemental concentration profiles of a medieval broach 36 : three layers (a 1.2 μm superficial layer, a 1.3 μm intermediate layer and a thick substrate), with different concentrations of copper, zinc, silver, gold and mercury in each layer.
For the third scenario we consider differential PIXE analysis of a non-homogeneous organic sample in air, performed with a realistic proton source obtained by means of Particle-In-Cell simulations. In particular, we model (2019) 9:9202 | https://doi.org/10.1038/s41598-019-45425-3 www.nature.com/scientificreports www.nature.com/scientificreports/ a realistic laser-driven ion source based on an advanced TNSA scheme, in which a double-layer target (a solid foil coated with a low density foam 37 ) is used to increase laser-target coupling and thus to enhance the maximum energy and the number of accelerated ions [16][17][18][19][20][21] . We exploit only the momentum distributions obtained from fully three-dimensional (3D) Particle-In-Cell (PIC) 38 simulations of such laser-driven ion source. The total number of simulated protons is selected in accordance with experimental data from the literature 14,27 , considering to perform multiple laser shots. In order to provide a realistic energy and angular distribution of the accelerated ions, we perform fully three-dimensional (3D) Particle-In-Cell (PIC) 38 simulations of such laser-driven ion source. We consider a 800 nm thick solid foil coupled with a near-critical 5 μm plasma layer, irradiated with a 30 femtoseconds, ultra-intense (up to ~5 × 10 19 W/cm 2 ) laser pulse. These parameters are realistic for a tightly focused (waist of 3 μm) 20 Terawatt Ti:Sapphire laser-system (e.g. they can be obtained scaling down the parameters reported in 39 for a 200 Terawatt laser, reducing the intensity while keeping the waist constant at 3 μm). Figure 2(b-d) show the main properties of this laser-driven proton source. Figure 2(d) shows a ring-like spatial distribution of the accelerated protons, which results into a hole in the energy spectrum at 0° (Fig. 2(b)). Similar distributions have been observed frequently in the literature (see 40 and references therein). In order to test our technique on a broader distribution we decide to tilt our source by ~3°. This means that the laser is still at normal incidence  www.nature.com/scientificreports www.nature.com/scientificreports/ with respect to the target, but the laser axis (which is also the target normal axis) is tilted by 3° with respect to the beam handler axis. We use 6 different proton spectra whose properties are changed varing the intensity of the laser beam in the PIC simulation. We choose a complex multi-layer sample, whose stratigraphic structure is representative of an oil painting 22,41 . The superficial layer has a thickness of 10 μm and it is composed by a protective organic varnish. The intermediate layer is 20 μm thick and contains a red pigment (HgS, the so-called Cinnabar), a lead-based white pigment (the so-called lead-white), some impurities (Ca and Fe) and the organic binder. The thick substrate (the so called imprimitura) is composed only by the lead-white and the organic binder.
The whole experimental apparatus assembled with Geant4 for the third scenario (PIXE in air) is illustrated in Fig. 1 (the set up for the first and second scenario neglects the presence of the air and exit window, being otherwise identical). Primary particles are generated at a punctual source placed in vacuum ( Fig. 1(a)). In principle, in case of a double layer target, we should consider protons, electrons and heavier ions. However we neglect heavier ions in the Monte Carlo simulations since in TNSA regime their number is usually orders of magnitude lower than for protons 42,43 . Moreover their cut-off energy per nucleon is usually lower, further reducing the X-ray yield due to heavier ions. Laser-driven proton sources generate also a copious amount of energetic electrons 10 , which would produce an unwanted X-ray signal from the sample. Thus the apparatus in Fig. 1 is designed to remove these energetic electrons. The particles are made to pass across an aperture slit ( Fig. 1(b1)) and a dipole magnet ( Fig. 1(b2)), which is a standard component to manipulate the laser-driven protons path 44,45 . The magnetic field steers electrons and protons trajectories. The electrons are completely deflected and they are stopped by a second slit ( Fig. 1(b3)). On the other hand, protons undergo a smaller deflection and they partially cross the second slit. Then, protons pass through an exit window ( Fig. 1(b4)) and they travel across the air before reaching the sample ( Fig. 1(b5)). The proton spectra at the sample surface are shown in Fig. 5(f). It is worth to point out that compact target chambers for laser-driven ion acceleration are already commercially available 46 and could be easily adapted for beam handlers. The X-ray detector has the features of a commercial Charged Coupled Device (CCD) (Fig. 1(b7)). The simulations account for the detection efficiency of the instrument. Under proper operating conditions (i.e. number of incident photons lower than the number of pixels), CCDs are able to perform single shot X-ray spectrum measurements 47 . This feature makes CCDs an interesting choice in order to record the X-ray spectra emitted by the sample during irradiation with laser-driven protons ( Fig. 1(b8)).

Results
We discuss the results of the three scenarios described in the previous section: PIXE analysis of a metallic homogeneous sample, differential PIXE analysis of a metallic non-homogeneous sample and differential PIXE analysis of a painting. We show that in all these cases the iterative algorithm is able to retrieve elemental concentration profiles in very good agreement with the real ones. All the details concerning the χ 2 minimization procedure are reported in the Methods section.
Typical PIXE measurements on delicate artifacts are performed with beam currents of the order of tens of pA for 100 s seconds 6 . This means ~10 10 protons on sample. Considering protons with energy greater than 0.2 MeV www.nature.com/scientificreports www.nature.com/scientificreports/ at the source, the efficiency of our beam-handling setup is ~10% in the most challenging case (i.e. the in-air laser-driven PIXE analysis and a 0 = 2). Therefore, we estimate that the laser-driven ion source should produce ~10 11 protons with energy greater than ~0.2 MeV at the source point. This can be achieved with ~10s of shots with a sub-100 TW class laser 27 (which can be operated at ~1 Hz). In this case the number of generated X-rays per shot would be low enough to avoid saturation of the CCD. It is also worth to point out that single-shot PIXE measurements can be performed with ~1 PW class laser, provided that a passive X-ray detector is used 34 .
In all the scenarios presented in this work we always perform Monte Carlo simulations with 10 11 protons with energy greater than ~0.2 MeV at the source. A typical sub-50 TW laser system accelerates ~10 9 protons per shot 14,27 . Therefore, our Monte Carlo simulations are representative of an experimental scenario in which the signal collected with several tens of laser shots is integrated. This applies also to the third case presented here, where PIC simulations are coupled to the Monte Carlo code to model a realistic source. Indeed, we use the momentum distributions from the PIC simulations to extract the initial properties of ~10 11 protons used in each Monte Carlo simulation.
As shown in Figs 4(g) and 5(g), in all the cases considered here the probed region on the samples is ~1 mm 2 (the probed region is not shown for the homogeneous sample analysis, however it is identical to that of Fig. 4(g)). This is perfectly adequate for most applications 48 . pIXe analysis of a metallic homogeneous sample. In the first scenario, an idealized proton source having a pure exponential energy distribution (see Fig. 3(a)) and a few degrees divergence is considered. Figure 3(b) shows that the agreement between the real concentrations and the retrieved ones is very good. For all the elements the discrepancy is less than 1%, which is comparable to the precision expected for traditional PIXE. However this excellent agreement has been obtained assuming a perfect knowledge of the proton energy spectrum (i.e. the black curve in Fig. 3(a)). On the other hand, it is well known that laser-driven proton sources are less stable than conventional particle accelerator sources, showing shot to shot fluctuations of the energy spectrum. These fluctuations might in principle hinder the applicability of laser-driven PIXE, thus it is imperative to asses their effect on the retrieved concentrations. This can be done straightforwardly for the simple scenario considered here. Keeping the X-ray yields obtained from the Monte Carlo simulations, we applied the iterative code considering a proton energy spectrum described by M C (i.e. α = 1 means that the iterative code uses a proton energy spectrum identical to that of the Monte Carlo simulation). α is varied between 0.8 and 1.2 in order to test the effect of fluctuations up to ±20%, as shown in Fig. 3(a). Figure 3(c) reports the ratio between the concentrations obtained with a given α and those obtained with α = 1. The graph shows a mild dependence of the relative error on the fluctuations of the proton energy spectrum: even for large (±20%) fluctuations, the relative error remains approximatively between ±10%.
Differential PIXE analysis of a metallic non-homogeneous sample and of a painting. Also in the second scenario, an idealized proton source is considered: six different exponential distributions with cut-off energies between 1 and 6 MeV (Fig. 4(f)). The concentration profiles ( Fig. 4(a,b)), reconstructed from the X-ray yields, show a good agreement with the original ones. In the third scenario (the painting) the output of six PIC simulations is used as a proton source in the Monte Carlo code (Fig. 5(f) shows the energy spectra at the sample surface). Concerning the analysis discussed here, the main difference with respect to the metallic case is that now the recorded X-ray peaks are representative of complex compounds rather than individual elements. The analysis clearly identifies the superficial protective varnish, the intermediate region containing the pigment, and the substrate. The concentration profiles of the elements are in satisfactory agreement with the real ones. Therefore, the system considered in this scenario (based on a 20 TW system and on the use of double layer targets) represents a realistic design for a laser-driven PIXE apparatus.

Conclusions
Laser-driven ion sources have been a thriving research topic in the last 2 decades, thanks to their many foreseen applications in a number of scientific, technological and societal fields. While for most of such applications major challenges still need to be tackled 49,50 , the results presented in this work for laser-driven PIXE clearly demonstrate the feasibility of a compact, versatile, cost-saving and portable laser-driven apparatus for IBA. This can be obtained combining a thorough theoretical description of the process with the most recent advanced laser-driven ion acceleration schemes and state-of-the-art laser systems. One can easily imagine further exciting developments, for example the possibility to exploit the same ion source for the production of neutrons with unique properties 51 , paving the way for a true revolution in the field of hadron-based advanced materials characterization techniques.

Methods
We developed an iterative code, based on χ 2 minimization algorithm, which solves the system of equations (3) or (4). We used this code to analyze the X-ray spectra obtained from Monte Carlo simulations of laser-driven PIXE experiments. For the third scenario, we performed Particle-In-Cell simulations to generate a realistic proton source.  55,56 . For X-ray mass www.nature.com/scientificreports www.nature.com/scientificreports/ absorption coefficients we used the online interface of the XCOM code 57 . For the proton stopping powers we used the Stopping and Range of Ions in Matter (SRIM) code 58 .
For the homogeneous sample, as initial guess, the elements are assumed to be equally concentrated. The differential PIXE analysis demands an initial guess also on the mass thickness of the fictitious layers. This choice should take into account the expected penetration depth of the protons and the expected scalelength of the inhomogeneities of the sample. We chose [0.65, 1.3, 2.6, 5.2,∞] mg/cm 2 both for the metal sample and for the painting. However, we introduced a free-parameter to readjust the mass thicknesses of the fictitious layers (e.g. [0.65ξ, 1.3ξ, 2.6ξ, 5.2ξ, ∞] mg/cm 2 with ξ constrained between 0.5 and 1.5). The fact that the analysis converges to good results for both the metallic sample and the painting, despite their very different structure, confirms the generality of this approach. As far as the initial guess on the concentrations is of concern, we followed a procedure analogous to that described in 22 .
For all the cases considered in this work, the χ 2 parameter has the general form: where the terms  Y i k and  ⁎ Y i k, are explained below. The mass concentrations W i (where i is the element index) are bounded between 0 and 1. In the case of the homogeneous sample analysis, only one measurement is necessary (i.e. the sum over k disappears).  Y i and  ⁎ Y i in equation (5) are, respectively, the normalized calculated yield and the normalized experimental yield (i.
This formulation allows to perform the χ 2 minimization without knowledge of the number of incident protons, which is particularly interesting from the experimental point of view.
As far as the differential PIXE analysis is concerned, χ 2 is a function of the free-minimizing parameter ξ and of the elemental concentrations W i,j of all the fictitious layers, which means 26 variables in the cases of interest here. Since the direct minimization of such function is particularly challenging we subdivided the procedure in multiple steps. We minimize χ 2 one element at a time (together with ξ) for all the elements and we repeat this procedure until the convergence is reached. Moreover, we did not used the normalized X-ray yields, but  Finally, in the case of the painting, the presence of a compound which does not provide detectable X-ray yields (i.e. the organic matrix) demands to impose the additional con- . Thus, we added a term λ∑ − ∑ ( ) to equation (5), where λ is a large number.
particle in cell (pIC) simulation. 3D Particle-In-Cell simulations were performed with the open-source code piccante 59 . We used a computational box of 120λ × 80λ × 80λ with a spatial resolution of 40 points per λ along x and 15 points per λ along ŷ and ẑ, where λ = 800nm is the laser wavelength. Time resolution is set at 98% of the Courant Limit. The P-polarized laser pulse had a cos 2 temporal profile (intensity FWHM of 30 fs) and a transverse Gaussian profile (i.e. E ~ E 0 exp(−r 2 /w 2 ) where E is the electric field, E 0 is the maximum electric field, r is the radius and w is the waist of 3 μm). The normalized peak laser intensity a 0 was set to 2, 2.5, 3, 3.5, 4, 4.5 (corresponding to a focused intensity between 8.7 ⋅ 10 18 W/cm 2 and 4.4 ⋅ 10 19 W/cm 2 ). The incidence angle was 0°. The target consisted in a 5 μm low-density layer with a density equal to 1 n c (sampled with 4 macro-electrons and 1 macro-ion with Z/A = 0.5 per cell), an overdense 800 nm thick foil with a density of 40 n c (sampled with 40 macro-electrons and 8 macro-ions with Z/A = 0.5 per cell) and a 80 nm thick contaminant layer with a density of 5 n c . In order to simulate a fully ionized hydrocarbon contaminant layer, the 5 n c ion charge was partitioned as follows: 4.285 n c for a species with Z/A = 0.5 and 0.715 for a species with Z/A = 1. 64 macro-electrons per cell and 125 macro-ions per species per cell where used for the contaminant layer. The electron population was initialized with a small (~eV) temperature. Periodic boundary conditions are used for both EM field and particles, however the box is large enough to avoid unphysical effects on the ion acceleration process. The duration of the simulations was 100 λ/c = 267 fs. A total of 6 PIC simulations was performed, one for each value for the normalized intensity a 0 . While in principle 3D PIC simulations can provide an absolute number for the total accelerated charge, this number depends on the choice of some parameters (namely the thickness and density of the contaminant layer) which might differ from those of a real experiment. For this reason in Fig. 2(c) we reported the energy spectra using arbitrary units. In order to use PIC results for the proton source in the Monte Carlo code, we calculated the distribution function d 3 N/dp x dp y dp z from the macro-proton phase space (p x , p y and p z are the Cartesian components of the momentum).
Geant4 monte carlo simulation toolkit. Geant4 (Geometry and Tracking) 60 is an abstract C++ base classes' toolkit dedicated to the Monte Carlo simulation of particles' transport through matter. It allows to implement the particle transport physics and the creation of secondary particles, covering a wide range of energies (from eV to TeV). Geant4 allows to simulate PIXE 61 reliably, as validated by several recent works [62][63][64] .
Our Monte Carlo simulations take into account protons, electrons, positrons, and photons. The secondary charged particles are tracked considering Bremsstrahlung, ionization and multiple scattering, while photoelectric effect, Compton scattering and pair production are activated for photons. The production cut is set equal to 0.5 μm for secondary particles in order to avoid infrared divergence. As far as these physical processes are concerned, we used the EmStandardPhysics_option3 module recommended by the Geant4 documentation for high accuracy with electrons and ions tracking. To properly model the proton ionization cross sections, the Energy-Loss Coulomb-Repulsion Perturbation-Stationary-State Relativistic Theory (ECPSSR) 53 www.nature.com/scientificreports www.nature.com/scientificreports/ In all the considered cases, the primary proton source is punctual and the primary particles can be electrons or protons.
Electrons are generated extracting their energy and angular divergence from an exponential and a uniform distributions respectively. Their energy spectrum was fitted from the PIC simulation with a 0 = 4.5. This is the worst case scenario since the highest energy electrons are generated.
As far as protons are concerned, for the metallic samples analysis a similar approach has been followed: the energy is extracted from an exponential distribution and a uniform angular distribution is assumed.
In the last case (i.e. the painting analysis), the proton source is modeled with the actual PIC simulations results. A PIC simulation provides a distribution of momenta. The primary proton momentum components are extracted with the Inverse Transform Sampling method from the PIC distribution and they are used to define the initial energy and propagation direction. To this aim we created a specific C++ class.
Several Geant4 simulations with different seeds of the Random Number Generator (RNG) were run in parallel and the final results were aggregated. A total of few 10 10 -10 11 primary protons were simulated for each case. We relied on the "HepJamesRandom" generator, which implements the Marsaglia-Zaman RANMAR algorithm. This RNG has the essential property of providing a large number (~8 ⋅ 10 8 ) of independent and very long sequences of pseudo-random numbers 65 .
It worth to point out that, since the Monte Carlo code allows to simulate one particle at a time, possible collective effects are neglected. This should not affect our results since collective effects have been reported to be negligible in this regime 66 .

Data Availability
The data supporting the findings of this work are available from the corresponding authors on request.