As discussed in Section 1, we performed two 3D PIC simulations with the PICCANTE code [22] to obtain reliable proton energy spectra from both a single-layer target and a DLT. In these simulations, we consider realistic laser parameters achievable with a 150 TW class laser. They are the same exploited in a recently published experimental work by [17]. The p-polarized pulse has a Gaussian envelope, 800 nm wavelength λ, a time duration of 30 fs, 2.8 μm FWHM focal spot and normalized laser intensity a ₀ = 16. The simulation box is x = 90λ along the laser propagation direction and the lateral dimensions are 60λ in the y and z directions. The spacial resolution is equal to 20 points per λ and the simulation duration is 100λ/c = 267 fs. We impose periodic boundary conditions in all directions.

The front surface of the targets is placed at x = 40λ. The bare target has a thickness of 1.875λ = 1.5 μm and 40n _c density ( n c = m e ω 0 2 ϵ 0 / e 2 is the critical density, m _e the electron mass, ω ₀ the laser frequency, ϵ ₀ the vacuum permittivity and e the electron charge). The density is sampled with 40 macro-electrons and 4 macro-ions with Z/A = 0.5 per cell. On the rear side of the target, there is a fully ionized hydrocarbon contaminant layer of thickness 0.1λ. Its density of 10n _c is partitioned in 5n _c for species with Z/A = 0.5 and 5n _c for species with Z/A = 1. The same number of macro-electrons, macro-ions with Z/A = 0.5 and Z/A = 1 per cell is set equal to 125.

For the simulation with the DLT, while keeping for the substrate the same parameters set for the single-layer target, we included a realistic three-dimensional nanostructured layer in front of the solid foil. It was obtained through the application of a Diffusion Limited Cluster-Cluster Aggregation (DLCCA) model [23]. The thickness of this layer is equal to 5λ = 4.0 μm and the resulting density is 3n _c , partitioned in 2.85n _c for a species with Z/A = 0.5 and 0.15n _c for a species with Z/A = 1. To simulate the carbon foam, we set 40 macro-electrons, 4 macro-ions with Z/A = 0.5 and 2 macro-ions with Z/A = 1 per cell. A snapshot of the simulation showing the interaction of the laser with the nanostructured DLT is reported in Figure 1A.

As shown in Figure 1B, with single-layer targets about 17% of the laser energy is absorbed by hot electrons at maximum (i.e., 25λ/c after the start of the simulation). Similarly, electrons from the DLT substrate absorb about 20% of the laser energy. On the other hand, the nanostructured layer absorbs 40% of the laser energy at 28λ/c. Overall, the presence of the low-density layer allows increasing by a factor of ∼ 3 the conversion efficiency of laser energy into hot electrons. Note that, after ∼ 50 λ / c the absorption reaches a plateau for both the single-layer target and DLT.

The proton spectra are retrieved at 70λ/c, and they are reported in Figure 1C with continuous lines. According to the PIC simulations, we can achieve an energy enhancement of about 2 in the presence of the nanostructure, being 11.6 MeV and 23.5 MeV the maximum energies achieved with the single-layer target and DLT, respectively. Moreover, considering protons having energies exceeding 2 MeV, their number is increased by a factor of 2.2 exploiting the DLT.

The proton energy spectra from PIC simulations are compared in Figure 1C with those reported in [17] and obtained with RCF stacks. The PIC results predict exactly the temperatures of the proton spectra for both kinds of targets. The experimental maximum energy of the accelerated protons with the single-layer target is equal to 10.6 MeV, thus resulting in a relative error of ∼ 8.6 % . In the DLT case, the maximum proton energy is slightly overestimated by PIC. Indeed, the experimental maximum energy of 18.8 MeV differs from the predicted value of ∼ 20 % . Overall, we can conclude that there is a quite fair agreement between our results and experimental data, proving the reliability of the 3D PIC simulations.

In order to quantify the number of ⁶⁴Cu isotopes that can be produced with monoenergetic and laser-driven proton sources, we performed Monte Carlo Geant4 [24] simulations. The goal is to compare the number of generated isotopes per unit incident protons (i.e., the yield) with laser-driven proton sources and conventional accelerators. For all irradiation conditions, we consider a pure Ni-64 sample of 5 cm thickness. We neglect the presence of other isotopes since enrichment levels higher than 95% can be achieved in samples for radioisotope production [25, 26]. It is placed 5 cm far from the source.

We performed two simulations with laser-driven protons providing the input proton spectra for single-layer target and DLT obtained from the PIC. Thus, the primary proton energies are sampled from the distributions reported in Figure 1C. The minimum energy is set to 2 MeV since, according to the ⁶⁴Ni(p, n)⁶⁴Cu reaction cross section reported in Figure 2A, no radioisotopes are produced below this threshold. Moreover, we carried out nine Monte Carlo simulations with monoenergetic protons having energies between 4 and 20 MeV. The energies provided by cyclotron machines for radioisotope productions usually lie in this range [1]. In the simulations involving laser-driven protons, we neglect the angular divergence of the beam since the irradiated material is large enough so that all particles reach the surface.

As far as the hadronic processes are concerned, we implemented the G4HadronElasticPhysicsHP and G4HadronPhysicsQGSP_BIC_AllHP physics lists to model elastic and inelastic interactions, respectively. In particular, the second one allowed us to exploit the TENDL-2019 cross section [27] for ⁶⁴Cu radioisotope production (i.e., the continuous filled line in Figure 2A). For each simulation, a total number of 4 × 10⁷ primary particles were simulated.

The number of produced ⁶⁴Cu isotopes per unit incident proton is reported in Figure 2. As expected, the yield achieved with monoenergetic protons monotonically increases with the primary particle energy. Around 20 MeV, it reaches a plateau since the proton energy is beyond the region where the cross section lies. Considering the laser-driven source with the single-layer target, the yield is Y = 3.3 × 10⁻⁵, thus comparable with that achievable with 4.5 MeV monoenergetic protons. Exploiting the DLT, the yield increases up to Y = 1.3 × 10⁻⁴, which is equivalent to the case of a 6.4 MeV monoenergetic proton source.

Exploiting the results obtained from Monte Carlo simulations, we can evaluate the activity achievable during irradiation with the laser-driven proton source. To this aim, we assume the same repetition rate of RR = 10 Hz that can be achieved with the laser system [28] considered for PIC simulations. For the single-layer target case, we assume N _p = 2.6 × 10¹⁰ accelerated protons per shot with energies higher than 2 MeV. This value has been obtained from the experimental scaling reported in [20]. The scaling is valid for 10–100 fs laser pulses and simple foils as targets. Then, we apply the gain factor of 2.2 for the proton number achievable with DLT. Therefore, the number of incident protons per shot with the near-critical nanostructured target is assumed to be 5.74 × 10¹⁰.

In order to retrieve the activity A during irradiation, we apply the following simple equation: A = R R ⋅ N p ⋅ Y ⋅ 1 − e − λ t (1) where λ = ln(2)/T _1/2 is the decay constant and T _1/2 = 12.7 h is the half life of ⁶⁴Cu. The results are reported in Figure 3.

First of all, we can notice that the activation rate and the saturation activity are one order of magnitude higher with the DLT compared to the single layer target, achieving approximately 4 MBq/h against 0.5 MBq/h and 73 MBq against 9 MBq, respectively. Overall, the highest activity we can achieve with the considered laser parameters would not allow performing clinical studies on human patients. Indeed, the values usually needed for treatments are of the order of 100 s MBq or higher [1, 29, 30]. They are achieved by exploiting cyclotrons providing 5–20 MeV monoenergetic protons and currents of the order of 10 s of μA. To achieve comparable performances, a more intense laser source and DLTs should be considered, as will be discussed in the next Section.

On the other hand, the activity values usually considered for pre-clinical studies exploiting mouse models lie in the range 5–100 MBq [31–36]. In particular, the indicative threshold of 5 MBq is achieved after 15.6 h of irradiation by exploiting single-layer targets. Vice versa, the same activity can be obtained in 78 min of elapsed time by adopting DLTs. These results show how near-critical nanostructured targets could play a fundamental role to achieve the activities required for pre-clinical studies by increasing the performances of the laser-driven proton source. In principle, a 150 TW class laser equipped with DLTs could provide sufficient ⁶⁴Cu isotopes for studies with mouse models. Notably, laser systems with parameters coherent with those considered for this study have become commercially available in recent years [37, 38]. On the other hand, important effort must be put in the future to develop target delivery systems compatible with DLTs and high repetition rate operation. This point is currently subject of research [39].

In the previous Sections, we have studied in detail the process with a fixed set of laser and target parameters by making use of 3D PIC simulations coupled with a Monte Carlo code. This approach is aimed at modeling the laser-driven radioisotope generation with the highest accuracy and demonstrating its potential in the view of practical applications. On the other hand, it would also be very interesting to extend the analysis to different laser and target configurations. In this Section, we aim at finding a suitable estimation of the total yield Y _T as a function of laser intensity a ₀ for optimized DLT conditions.

The total ⁶⁴Cu production yield Y _T in a semi-infinite, isotopically pure ⁶⁴Zn converter can be expressed by: Y T = N C u N p = N A v M ρ ∫ 0 ∞ f E p ∫ 0 E p σ E S p E d E d E p (2) where N _n is the total number of radioisotopes produced upon converter irradiation with N _p protons, N _Av is the Avogadro’s number, M is the converter atomic mass, ρ is the converter density, f(E _P ) is the proton spectrum normalized to the total number of protons N _p , S _p (E) is the proton linear stopping power of the converter material and σ(E) is the total cross section for the ⁶⁴Zn(p, n)⁶⁴Cu reaction, as a function of the proton energy E.

In order to develop an explicit analytical expression for Y _T we look for some suitable approximation of the integrals in Equation 2. Firstly, let’s define the differential yield g(E) and the monoenergetic yield h(E _p ) as: g E = N a t σ E S E (3) and h E p = ∫ 0 E p g E d E = ∫ 0 E p N a t σ E S E d E (4) in which N _at = ρ N _Av /M in the atomic density of the converter. We decided to approximate the g(E) with a Gaussian function Γ(E): g E → Γ E = Y m a x w π exp − E − E 0 2 w 2 (5)

The values of the parameters Y _max (2.023 × 10⁻³), E ₀ (10.25 MeV) and w (4.38 MeV) are obtained with a least square fitting of Γ(E) over the reference g(E). To this purpose, g(E) (Equation 3) is calculated by assuming N _at ≈ 8.379 cm⁻³ as the atomic density of a pure ⁶⁴Ni target, S(E) according to the Bethe’s formula for proton stopping power [40], and σ(E) from the TENDL nuclear data library [27].

According to the literature concerning TNSA [3], the normalized proton distribution is modeled as an exponential function characterized by a slope T _p and with a cut-off corresponding to the maximum proton energy ϵ _M : f ( E p ) = exp ( ϵ M / T p ) T p ( exp ( ϵ M / T p ) − 1 ) exp ( − E p T p )  if  E p ≤ ϵ M f ( E p ) = 0  if  E p > ϵ M

The integral ∫ 0 ϵ M f ( E p ) ∫ 0 E p Γ ( E ) d E d E p can be solved analytically, but the result can be simplified by assuming that ϵ _M ≫ T _p , E ₀, which is verified for DLT targets and a ₀ > 15. Under this assumption and taking only the leading term in T _p one gets the following approximation for the total yield: Y T ≈ Y m a x ⁡ exp w 2 − 4 E 0 T p T p 2 (6)

Since Y _max , w and E ₀ are known, the last step is to express T _p as a function of laser and target parameters. By considering the PIC spectra shown in Figure 1C one can get T _p = 1.86 MeV for the bare foil target and T _p = 3.24 MeV for the DLT. These values are close to the electron temperature in both cases (1.84 MeV and 3.36 MeV respectively), and therefore we assume that T _p ≃ T _e , where T _e can be estimated through the theoretical model presented in [14]. Although this approximation may seem rather crude it yields reasonable results nonetheless, as detailed in the following.

The comparison between the approximated ⁶⁴Cu yield achieved with bare targets and DLTs as a function of the normalized laser intensity is reported in Figure 4A. For a ₀ < 20, there is a remarkable difference between the yields obtained with the two kinds of targets. Indeed, the proton maximum energies are close to the energy range 5–15 MeV where the maximum of the cross section is located. Thus, a strong enhancement of the proton energies with DLTs results in a relevant increment of the yields in comparison with bare targets. On the other hand, the curves of the yields start having the same slope for a ₀ > 30, highlighting that the maximum proton energies have largely exceeded 15 MeV for both kinds of targets.

We notice a fair agreement between the yields from the model and those obtained from the Monte Carlo. They are practically superimposed for the bare-target case. The yield resulting from the Monte Carlo considering the DLT slightly underestimates the analytical estimation. The discrepancy can be ascribed to the approximations introduced in the model to treat the physical processes, as well as to the different shapes of the proton spectrum (resulting from the PIC in the Monte Carlo and purely exponential in the model). In addition, there is a small difference between the target parameters exploited in the PIC simulations and those obtained from the optimization with the model.

By exploiting the analytical yields, we can apply Equation 1 to obtain the activity as a function of the elapsed (i.e., irradiation) time for all values of a ₀ (corresponding to laser energies in the range 0.7–17.5 J). The results are reported in Figure 4B. Here, we keep constant RR = 10 Hz for all values of normalized laser intensities. While this high repetition rate has not yet been achieved for the higher intensities considered here, we can reasonably assume that laser technology will provide such systems in the future. Indeed, state-of-the-art multi-100 s TW and PW class lasers can work at ∼ 0.1 − 1 Hz repetition rate [37, 38, 41, 42]. For the number of accelerated protons per shot N _p , we have assumed the experimental scaling already exploited in Section 2.3 and presented in [20] for bare targets. Then, we multiply all values of N _p by the enhancement factor of 2.2 provided by the use of DLTs and obtained from the PIC simulations.

The isolines in Figure 4B lie in correspondence with constant values of activity. By way of comparison, we identified the 150 TW (i.e., a ₀ ∼ 16), DLT case study presented in the previous Sections as a dotted red line. It intercepts the 5 MBq activity isoline at 72 min, which is quite in agreement with the 78 min reported in Section 2.3. We noticed that, considering intensity values between 20 and 50, activities compatible with pre-clinical studies could be achieved in minutes rather than hours. To obtain A ≫ 100 MBq relevant for clinical studies, elapsed times of the order of 1 h or more should be considered. The same result can be deduced also from the activity per unit of elapsed time at the start of the irradiation. By looking at the activity rate as a function of a ₀ (Figure 4C) one should conclude that a ₀ > 30 will be sufficient to exceed 100 MBb/h.