Superintense Laser-driven Ion Beam Analysis

TL;DR

This paper explores the use of superintense laser-driven proton sources as a compact, versatile alternative to traditional ion beam analysis tools for advanced material characterization.

Contribution

It develops an analytical and numerical framework for laser-driven ion beam analysis and proposes a realistic design for a portable, flexible facility.

Findings

Framework successfully models laser-driven ion analysis

Design demonstrates potential for portable analysis facilities

Enhances flexibility and accessibility of ion beam techniques

Abstract

Ion beam analysis techniques are among the most powerful tools for advanced material 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 material characterization.

Equations

$$Y_i=N_p\frac{\Delta \Omega }{4\pi }{\epsilon }_i\frac{N_{av}}{M_i}{W}_i{\int }_{E_0}^{E_f}{\sigma }_i(E){\omega }_i{e}^{-{\mu }_i{\int }_{E_0}^{E^{\prime }}\frac{d{E}^{\prime }}{S(E^{\prime })}\frac{\cos \theta }{\cos \varphi }}\frac{dE}{S(E)}\equiv N_p{\Gamma }_i(E_p)$$

$$Y_i^k=N_p^k\frac{\Delta \Omega }{4\pi }{\epsilon }_i\frac{N_{av}}{M_i}\sum _{j=1}^{J-1}{W}_{i,j}{e}^{-{\sum }_{l=1}^{j-1}(\frac{\mu }{\rho }{)}_{i,l}\frac{{\rho }_l{R}_l-{\rho }_l{R}_{l-1}}{\cos \varphi }}{\int }_{E_{p,j+1}^k}^{E_{p,j}^k}{\sigma }_i(E){\omega }_i{e}^{-{\mu }_{i,j}{\int }_{E_p}^{E^{\prime }}\frac{d{E}^{\prime }}{S_j(E^{\prime })}\frac{\cos \theta }{\cos \varphi }}\frac{dE}{S_j(E)}$$

$$Y_i=\frac{\Delta \Omega }{4\pi }{\epsilon }_i\frac{N_{av}}{M_i}{W}_i{\int }_{E_{p,min}}^{E_{p,max}}{f}_p(E_p){\int }_{E_p}^0{\sigma }_i(E){\omega }_i{e}^{-{\mu }_i{\int }_{E_p}^{E^{\prime }}\frac{d{E}^{\prime }}{S(E^{\prime })}\frac{\cos \theta }{\cos \varphi }}\frac{dE}{S(E)}d{E}_p$$

$$Y_i^k=\frac{\Delta \Omega }{4\pi }{\epsilon }_i\frac{N_{av}}{M_i}\sum _{j=1}^{J-1}{W}_{i,j}{e}^{-{\sum }_{l=1}^{j-1}(\frac{\mu }{\rho }{)}_{i,l}\frac{{\rho }_l{R}_l-{\rho }_l{R}_{l-1}}{\cos \varphi }}{\int }_{E_{p,min}^k}^{E_{p,max}^k}{f}_p^k(E_p){\int }_{E_{p,j+1}^k}^{E_{p,j}^k}{\sigma }_i(E){\omega }_i{e}^{-{\mu }_{i,j}{\int }_{E_p}^{E^{\prime }}\frac{d{E}^{\prime }}{S_j(E^{\prime })}\frac{\cos \theta }{\cos \varphi }}\frac{dE}{S_j(E)}d{E}_p$$

$$\chi^{2}=\sum_{i,k}\Big{(}\frac{\tilde{Y}_{i}^{k}-\tilde{Y}_{i}^{k,*}}{\sqrt{\tilde{Y}_{i}^{k,*}}}\Big{)}^{2}$$

Full text

Introduction

Materials characterization is of crucial importance for basic research as well as for a wide number of applications, ranging from virtually every field of advanced technology to cultural heritage purposes. Ion Beam Analysis (IBA)[1, 2] is an important family of modern analytical techniques to probe the composition and the surface structure of solid samples with MeV ion beams. IBA has found use in an impressive number of applications [3], including biomedical elemental analysis [4], semiconductor industry [5], cultural heritage studies[6], forensic analysis[7], nuclear fusion research[8]. The appeal of these techniques is due to their non-destructive nature and their unparalleled detection capabilities.

The layout of an IBA apparatus involves an ion source, a target chamber (where the sample is placed) and suitable diagnostics. The interaction of the ions with the sample generates a secondary emission, whose properties allow to retrieve the structure and the elemental composition of the sample.

Despite the growing interest for IBA and its widespread adoption by industries and research institutes, in most cases the oldest accelerator technologies, mostly Van de Graaf and Tandem systems, are adopted. As a consequence, al the available IBA techniques have been developed considering only monochromatic, conventional accelerator-based ion beams. Conventional ion sources are invariably affected by a number of limitations: non-tunable energy, radioprotection issues [9], high costs, non-portable size. Ion sources driven by superintense lasers represent a very interesting alternative[10]. The investigation of laser-driven particle sources is a thriving research area, as demonstrated by the growing worldwide investments in large-scale facilities[11], such as the pan-European Extreme Light Infrastructure[12].

The most relevant properties of laser-driven ion sources (maximum ion energy and ion number per laser shot) already exceed, in principle, the requirements for several IBA techniques. In addition, recent progress in exploiting non-conventional targets allows to achieve such beam parameters with reduced laser requirements, which may now realistically fall into the class of compact, table-top Terawatt systems [13, 14]. Moreover, laser-driven ion sources offer the possibility to easily and rapidly change the energy of the accelerated ions. This can be done simply by tuning the laser pulse energy, with concrete, major advantages in terms of flexibility. As far as radioprotection is concerned, laser-driven ion sources should be more manageable since the radiation fields originate from a small region surrounding the laser-target interaction point. Thus, laser-driven ion accelerators for IBA offer a concrete route for increased flexibility and more compact, cheaper systems with reduced radioprotection concerns, which could revolutionize the entire field. On the other hand, the peculiar nature of these sources, invariably characterized by a broad energy distribution[10], requires to address novel, so-far unsolved, challenges to make laser-driven IBA a reality.

In this work we demonstrate the full potential of laser-driven ion sources for IBA, by means of a detailed theoretical-numerical investigation of a particularly interesting technique: Proton Induced X-ray/Gamma-ray Emission, PIXE/PIGE. Specifically, we address and solve the following points: i) extension of the theoretical models to interpret PIXE/PIGE data, in order to make them suitable for laser-driven ions ii) development of realistic numerical simulations of a laser-driven PIXE experiment to show the capability of retrieving material properties with the developed numerical methods iii) design of a proof-of-concept, transportable, laser-driven source for IBA (see figure 1), obtained combining a compact, TW-class laser[15, 16] with a non-conventional targetry solution based on nanostructured double-layer targets [17, 18, 19, 20, 21, 22].

PIXE/PIGE theory with laser-driven ions

Laser-driven ions are invariably characterized by a broad-spectrum energy distribution[10], whose shape depends on the specific acceleration mechanism. Since PIXE/PIGE formalisms have been developed only for monochromatic ion sources[2], extending the theory for arbitrary energy distributions is a crucial step. Consider first conventional PIXE analysis, performed on a thick homogeneous sample composed by $I$ elements of atomic weights $M_{i}$ and mass concentrations $W_{i}$, irradiated with a monoenergetic proton beam of initial energy $E_{0}$ and total number of protons $N_{p}$. The number of X-rays $Y_{i}$, associated to one characteristic line for the i-th element, is [2, 24]:

[TABLE]

where $\Delta\Omega$ and $\varepsilon_{i}$ are the subtended solid angle and the efficiency of the detector, $N_{av}$ is the Avogadro’s number, $E_{f}$ is the final proton energy after slowdown in matter, $\sigma_{i}(E)$ is the ionization cross section, $\omega_{i}$ is the fluorescence yield, $S(E)$ is the proton stopping power in the sample, $\mu_{i}$ is the X-ray attenuation coefficient, $\theta$ is the proton impact angle and $\phi$ is the X-ray emission angle. The integral takes into account the production of X-rays along the proton range, as well as their attenuation inside the material, and it is therefore computed over the projectile energy.

By neglecting the secondary fluorescence contribution, quantitative PIXE analysis can be performed by solving a system of $I$ equations like equation (1), where the X-ray yields $Y_{i}^{*}$ are measured experimentally and the unknowns are the mass concentrations $W_{i}$ of the elements. The system is non-linear because both the X-ray attenuation coefficients $\mu_{i}$ and the proton stopping power $S(E)$ depend upon the sample composition. Accordingly, the solution can be found by minimizing a parameter like $\chi^{2}=\sum_{i}[(Y_{i}-Y_{i}^{*})/\sqrt{Y_{i}^{*}}]^{2}$, making use of an iterative procedure.

If the sample is non-homogeneous, the goal of the analysis is to retrieve the concentration profiles of the element as a function of the thickness. In general, X-rays are emitted along the whole path of the protons inside the sample and the emission cannot be directly related to the deposited energy. In other words, it is not possible to simply assume that X-ray emission is representative of the sample composition at the Bragg peak. Therefore, a dedicated procedure must be followed as far as a non-homogeneous sample is of concern. In this case, Differential PIXE analysis can be applied by ideally splitting the sample in a finite number $J$ of fictitious layers having mass thickness $\rho_{j}R_{j}$ and homogeneous composition $W_{i,j}$, as shown in figure 2(a). Assuming to perform $K$ measurements at different proton energies $E_{p}^{k}$, the following system of equations can be written [25, 24]:

[TABLE]

where $N_{p}^{k}$ is the number of protons with energy $E_{p}^{k}$, $(\frac{\mu}{\rho})_{i,j}$ and $S_{j}(E)$ are the X-ray attenuation coefficient and proton stopping power associated to the $j$-th layer, $E_{p,j}^{k}$ and $E_{p,j+1}^{k}$ are the proton energies at the boundaries of the $j$-th layer (see figure 2(a)). The system of $I\times K$ equations can be solved if $K\geq J$ (i.e. if a number of measurements greater than or equal to the number of fictitious layers is performed) by $\chi^{2}$ minimization. The solution is represented by the set of discrete values of $W_{i,j}$ and mass thicknesses $\rho_{j}R_{j}$ which best fit the experimental X-ray yields $Y_{i}^{k*}$.

We extended this formalism to account for a non-monoenergetic proton energy spectrum, represented by the function $f_{p}(E_{p})=dN_{p}(E_{p})/dE_{p}$, which we assume to be different from zero in the energy interval from $E_{p,min}$ to $E_{p,max}$. The number of X-rays $dY_{i}$ generated by the incident protons with energy between $E_{p}$ and $E_{p}+dE_{p}$ is $dY_{i}=dN_{p}\Gamma_{i}=f_{p}\Gamma_{i}dE_{p}$. Therefore, integrating over the incident proton energy, equation (1) can be generalized as follows:

[TABLE]

With the choice $f_{p}=N_{p}\delta(E_{p}-E_{0})$ (i.e. a monochromatic distribution) equation (3) reduces to equation (1).

Also the case of Differential PIXE can be generalized for $K$ arbitrary energy distributions. Equation (2) becomes:

[TABLE]

Here, the quantities $f_{p}^{k}(E_{p})$, $E_{p,max}^{k}$ and $E_{p,min}^{k}$ 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[26]: the ionization cross section $\sigma_{i}(E)$ and the fluorescence yield $\omega_{i}$ must be replaced with the proper nuclear reaction cross section and isotopic abundance, while $\gamma$-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 [27], it has been used extensively as diagnostic tool in laser plasma interaction experiments [28] and has been demonstrated in a repetitive regime [29]. 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 $\sim 100$ MeV [30, 31] 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) [32] or Radiation Pressure Acceleration (RPA) [33, 34, 35] 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, preventing the use of 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 measurement, 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 [36] 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 investigated 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: $f_{p}(E_{p})\sim\exp\left({-E_{p}/T_{p}}\right)$, 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}=5\textrm{ MeV}$ and $T_{p}=0.7\textrm{ MeV}$, 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. We consider a sample composed by iron, copper, zinc, tin and lead, which is representative of a Roman sword-scabbard [37].

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). In this scenario, the sample composition is inspired by the elemental concentration profiles of a medieval broach [38]: three layers (a 1.2 $\mu$m superficial layer, a 1.3 $\mu$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 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[39]) is used to increase laser-target coupling and thus to enhance the maximum energy and the number of accelerated ions[17, 18, 19, 20, 21]. In order to provide a realistic energy and angular distribution of the accelerated ions, we perform fully three-dimensional (3D) Particle-In-Cell (PIC) [40] simulations of such laser-driven ion source. We consider a 800 nm thick solid foil coupled with a near-critical 5 $\mu$m plasma layer, irradiated with a 30 femtoseconds, ultra-intense (up to $\sim 5\times 10^{19}\textrm{W/cm}^{2}$) laser pulse. These parameters are realistic for a tightly focused (waist of $3\mu$m) 20 Terawatt Ti:Sapphire laser-system (e.g. they can be obtained scaling down the parameters reported in [41] for a 200 Terawatt laser, reducing the intensity while keeping the waist constant at 3 $\mu$m). Figures 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^{\circ}$ (fig. 2(b)). Similar distributions have been observed frequently in the literature (see [42] and references therein). In order to test our technique on a broader distribution we decided to tilt our source by $\sim 0^{\circ}$. This means that the laser is still at normal incidence with respect to the target, but the laser axis (which is also the target normal axis) is tilted by $3^{\circ}$ 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[43, 24]. The superficial layer has a thickness of 10 $\mu$m and it is composed by a protective organic varnish. The intermediate layer is 20 $\mu$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 figures 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 (figure 1 (a)). In principle, in case of a double layer target, we should consider protons, electrons and heavier ions. However we neglected heavier ions in the Monte Carlo simulations since in TNSA regime their number is usually orders of magnitude lower than for protons [44, 45]. 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 figure 1 is designed to remove these energetic electrons. The particles are made to pass across an aperture slit (figure 1 (b1)) and a dipole magnet (figure 1 (b2)), which is a standard component to manipulate the laser-driven protons path [46, 47]. The magnetic field steers electrons and protons trajectories. The electrons are completely deflected and they are stopped by a second slit (figure 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 (figure 1 (b4)) and they travel across the air before reaching the sample (figure 1 (b5)). The proton spectra at the sample surface are shown in figure 5 (f). It is worth to point out that compact target chambers for laser-driven ion acceleration are already commercially available [48] and could be easily adapted to beam handlers. The X-ray detector has the features of a commercial Charged Coupled Device (CCD) (figure 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[49]. 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 (figure 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 $\chi^{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 100s seconds [6]. This means $\sim 10^{10}$ protons on sample. Considering protons with energy greater than 0.2 MeV at the source, the efficiency of our beam-handling setup is $\sim 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 $\sim 10^{11}$ protons with energy greater than $\sim 0.2$ MeV at the source point. This can be achieved with $\sim 10$s of shots with a sub-100 TW class laser [29] (which can be operated at $\sim 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 $\sim 1$ PW class laser, provided that a passive X-ray detector is used [36].

As shown in figures 4(g) and 5(g), in all the cases considered here the probed region on the samples is $\sim 1\textrm{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 [50].

PIXE analysis of a metallic homogeneous sample

In the first scenario, an idealized proton source having a pure exponential energy distribution (see figure 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 $f_{p}(E_{p},\alpha)\sim\exp\left({-E_{p}/(\alpha T_{MC})}\right)\Theta(E_{p}-\alpha E_{MC})$ (i.e. $\alpha=1$ means that the iterative code uses a proton energy spectrum identical to that of the Monte Carlo simulation). $\alpha$ was varied between 0.8 and 1.2 in order to test the effect of fluctuations up to $\pm 20\%$, as shown in figure 3(a). Figure 3(c) reports the ratio between the concentrations obtained with a given $\alpha$ and those obtained with $\alpha=1$. The graph shows a mild dependence of the relative error on the fluctuations of the proton energy spectrum: even for large ($\pm 20\%$) fluctuations, the relative error remains approximatively between $\pm 10\%$.

Differential PIXE analysis of a metallic non-homogeneous sample and of a painting

Also in the second scenario, an idealized proton source was considered: six different exponential distributions with cut-off energies between 1 and 6 MeV (figure 4(f)). The concentration profiles (figures 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 (figure 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 18 years, 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 [51, 52], 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 [53], paving the way for a true revolution in the field of hadron-based advanced material characterization techniques.

Methods

We developed an iterative code, based on $\chi^{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.

$\chi^{2}$ minimization algorithm

The minimization algorithm used to retrieve the elemental composition of the samples relies on the open source dlib C++ library [54], in particular on the BOBYQA algorithm. All the physical parameters used in the algorithm are taken from standard libraries. For the ionization cross sections we used the Energy-Loss Coulomb-Repulsion Perturbed-Stationary-State Relativistic (ECPSSR) theory [55, 56]. For fluorescence yields we used values tabulated in the Evaluated Atomic Data Library (EADL) [57, 58]. For X-ray mass absorption coefficients we used the online interface of the XCOM code [59]. For the proton stopping powers we used the Stopping and Range of Ions in Matter (SRIM) code [60].

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,\infty]$ $\textrm{mg}/\textrm{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\xi,1.3\xi,2.6\xi,5.2\xi,\infty]$ $\textrm{mg}/\textrm{cm}^{2}$ with $\xi$ 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 [24].

For all the cases considered in this work, the $\chi^{2}$ parameter has the general form:

[TABLE]

where the terms $\tilde{Y}_{i}^{k}$ and $\tilde{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). $\tilde{Y}_{i}$ and $\tilde{Y}_{i}^{*}$ in equation (5) are, respectively, the normalized calculated yield and the normalized experimental yield (i.e. $\tilde{Y}_{i}=Y_{i}/\sum_{i}Y_{i}$ and $\tilde{Y}_{i}^{*}=Y_{i}^{*}/\sum_{i}Y_{i}^{*}$). This formulation allows to perform the $\chi^{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, $\chi^{2}$ is a function of the free-minimizing parameter $\xi$ 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 $\chi^{2}$ one element at a time (together with $\xi$) for all the elements and we repeat this procedure until the convergence is reached. Moreover, we do not used the normalized X-ray yields, but $\tilde{Y}_{i}^{k}$ and $\tilde{Y}_{i}^{k,*}$ coincide with calculated $Y_{i}^{k}$ and experimental $Y_{i}^{k,*}$ ones. 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 constrain $\sum_{i}W_{i,j}=1$. Thus, we added a term $\lambda\sum_{i}{(1-\sum_{i}W_{i,j})}^{2}$ to equation (5), where $\lambda$ is a large number.

Particle In Cell (PIC) simulation

3D Particle-In-Cell simulations were performed with the open-source code piccante[61]. We used a computational box of $120\lambda\times 80\lambda\times 80\lambda$ with a spatial resolution of 40 points per $\lambda$ along $\hat{x}$ and 15 points per $\lambda$ along $\hat{y}$ and $\hat{z}$, where $\lambda=800$nm 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\sim E_{0}\exp\left({-r^{2}/w^{2}}\right)$ where $E$ is the electric field, $E_{0}$ is the maximum electric field, $r$ is the radius and $w$ is the waist of 3 $\mu$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\cdot 10^{18}\leavevmode\nobreak\ \textrm{W/cm}^{2}$ and $4.4\cdot 10^{19}\leavevmode\nobreak\ \textrm{W/cm}^{2}$). The incidence angle was $0^{\circ}$. The target consisted in a 5 $\mu$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 ($\sim$ 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 $\lambda/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 figure 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) [62] 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[63] reliably, as validated by several recent works[64, 65, 66].

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 $\mu$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) [67] is selected.

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 ($\sim 8\cdot 10^{8}$) of independent and very long sequences of pseudo-random numbers [68].

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 [69].

Data availability

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

Author contribution statement

M.P. conceived the project and supervised all the activities. F.M. wrote the iterative code with the assistance of L.F, developed the theoretical model and performed the Monte Carlo simulations. L.F. performed the Particle-In-Cell simulations. All the authors contributed equally to the preparation of the manuscript.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (ENSURE grant agreement No 647554). We also acknowledge LISA and Iscra access schemes to MARCONI HPC machine at CINECA(Italy) via the projects LAST, LIRF and EneDaG.

Competing interests

The authors declare no competing interests.