Neutron beams are commonly employed in scientific research, medicine, and industry for a wide range of applications.1 In practice, they are generated from nuclear reactions initiated by accelerator proton beams. Conventional neutron source facilities range from compact tubes to large-scale linacs like the Spallation Neutron Source (Oak Ridge, USA)2 and the European Spallation Source (Lund, Sweden)3 currently under construction, where 1–2 GeV protons hit a heavy-metal target to produce neutrons through spallation reactions. These consist of a cascade of binary collisions between the incident projectile and the nucleons inside the target nuclei, followed by de-excitation (or evaporation) of the excited nuclei, leading to the emission of neutrons, and also, to a smaller extent, protons, alpha particles, lighter heavy ions, gamma rays, etc.4

The production of bright neutron beams using high-power short-pulse lasers was demonstrated in the early 2000s (see Ref. 5 and references therein for an overview) and has since been actively investigated. Laser-generated neutrons have already been utilized for a variety of purposes, such as materials testing for fusion experiments,6 nondestructive imaging,7 and studies of equations of state via neutron resonance spectroscopy.8,9 Such neutron sources exploit laser-driven energetic protons10 (with current record high energies of ∼100 MeV),11,12 electrons,13 or gamma-ray photons14 as the primary drivers, with typical cross sections in the barn range.5,15 Most previous studies on this topic have been focused on improving the yield10,14 or the energy spectrum16 of the emitted neutrons.

Progress in high-brightness neutron sources is a necessary step toward the laboratory production of neutron-rich isotopes via rapid neutron captures (r-process), which would allow nuclear physics models to be tested17 and improve our understanding of the formation of heavy nuclei in the Universe.18,19 Half of the elements heavier than iron (Z = 26), and all those beyond bismuth (Z = 83), are indeed believed to originate from the r-process during cataclysmic astrophysical events such as supernova explosions or neutron star mergers.20 To compensate for the short lifespan (in the millisecond range) of the intermediate isotopes, a minimum neutron flux >10²⁰ n cm⁻² s⁻¹ is estimated to be necessary for the r-process to operate.21 This value is several orders of magnitude above the capability of conventional accelerator-based facilities (∼10¹⁶ n cm⁻² s⁻¹),22 but also significantly larger than the current record high flux (∼10¹⁸ n cm⁻² s⁻¹) obtained with intense short-pulse lasers.10,14 Neutron fluxes as high as ∼10²⁴ n cm⁻² s⁻¹ can be attained at large-scale laser fusion facilities,23 but with limited user access and very few shots per experiment. Systematic laboratory investigations of r-process nucleosynthesis therefore require further development of laser-based neutron sources.

Laser acceleration of proton beams should greatly benefit from next-generation petawatt or multi-petawatt facilities, delivering pulse intensities in excess of 10²¹ W cm⁻².24–31 At such intensities, the dominant ion acceleration mechanism is expected to transition from target normal sheath acceleration (TNSA)32 to radiative pressure acceleration (RPA) or light-sail acceleration (LSA).33 The accompanying increase in proton energy above the 100 MeV level should trigger spallation reactions in a secondary neutron-producing target, entailing the emission of multiple neutrons per incident proton (see Fig. 1).

In Fig. 1, the cross sections of various neutron-producing reactions (as given by the ENDF/B-VIII database34) are plotted as a function of the projectile (proton or photon) energy (solid lines). Already at relatively low proton energies of tens of MeV, easily produced by current high-intensity laser facilities, neutrons are efficiently produced. It is therefore clear that to access reactions with increased neutron multiplicity, higher projectile energies are needed.

In the present study, we will focus on exploiting protons as a driver to induce the desired neutron beam. This is justified, because the cross sections of photonuclear reactions in Pb (black dashed curve in Fig. 1), although comparable to those of proton-induced reactions around 10 MeV, are 100 times lower at higher energy. Thus, integrating the expected production of neutrons, induced by either photons or protons, in the cases that will be addressed at the LMJ-PETAL35 and Apollon36 laser facilities shows that the neutron yield induced by protons is at least ten times higher than that due to photons. The characteristics of the proton beams considered here will be detailed below, while to estimate the photons that can be generated by PETAL and by Apollon, we refer respectively to Ref. 37 and to Refs. 38 and 39. Hence, photonuclear reactions will be neglected in a first-order approach in the physical situations treated in the following. Yet, these could prove to be influential at more extreme laser intensities (≥10²³ W cm⁻², not presently achievable) for which massive high-energy synchrotron radiation is expected to arise.40

This paper reports on a numerical study aiming to characterize the neutron yield and flux achievable at the petawatt-class LMJ-PETAL and Apollon laser facilities. Our results are also relevant to similar laser systems, such as ELI Beamlines41 and ELI-NP.42 Overall, this is done by combining particle-in-cell (PIC) simulations of the laser-driven particle acceleration in a thin-foil target with Monte Carlo calculations of the proton-induced nuclear reactions in a secondary target. In Sec. II, we first investigate numerically how the production of spallation neutrons varies with the target material and incident proton energy. In Sec. III, we present the results of PIC simulations of laser-based proton acceleration from solid foils under conditions accessible at LMJ-PETAL and Apollon. The generation of spallation neutrons from a Pb convertor by the PIC-predicted proton beams is studied in Sec. IV via Monte Carlo simulations, which predict that neutron yields of 10⁸–10¹¹ n sr⁻¹ and neutron fluxes of 10²³–10²⁴ n cm⁻² s⁻¹ are achievable. Finally, our results are summarized and discussed in Sec. V.

We start by briefly examining the neutron yield from spallation as a function of the converter target material and the projectile proton energy in a range within the reach of present or near-future laser systems.11,12,33 For this purpose, we used the FLUKA 3D Monte Carlo code43,44 to simulate the nuclear reactions induced during irradiation of a convertor target by a monokinetic and monodirectional proton beam. The proton beam energy was varied from ϵ_p = 25 MeV to 1 GeV. Four target materials were considered: aluminum (Z = 13), copper (Z = 29), silver (Z = 47), and lead (Z = 82), a standard material for spallation purposes.

The simulated target was a 50-cm-radius cylinder of variable length L. Introducing λ, the energy-dependent projected range of a proton due to ionization and excitation,45 five L/λ values (0.2, 0.4, 0.6, 0.8, 1) were considered for each material and input proton energy. The values of λ corresponding to our parameter range are given in Table I.

Figure 2 shows the L/λ dependence of the neutron multiplicity M_n, i.e., the number of neutrons produced per incident proton as a function of its input energy ϵ_p, for the different materials under consideration. Collisional stopping and scattering of the projectiles and product neutrons are taken into account in these simulations. The main result is that, whatever the material, M_n rises sharply (i.e., approximately quadratically) with ϵ_p when L/λ is kept constant. In Ag and Pb, M_n approaches unity (a usual criterion for the onset of spallation) for ϵ_p ≃ 250 MeV and L/λ ≲ 1. At higher ϵ_p, M_n ≥ 1 can be achieved in lower L/λ targets. The maximum neutron multiplicity (M_n ∼ 10) is obtained in Pb with ϵ_p = 1 GeV and L/λ = 0.4–0.6. It should be noted that the increasing trend of M_n with ϵ_p ceases beyond ϵ_p ≃ 0.5–1 GeV when L/λ ≳ 0.5. This is a known behavior in spallation studies, ascribed to the increasingly significant contribution of pion production to proton energy losses.4

At fixed proton energy, M_n is predicted to rise by a relatively modest (∼3×) factor when the normalized target thickness is increased from L/λ = 0.2 to 1. Finally, at fixed ϵ_p ≤ 0.5 GeV and L/λ, M_n shows a moderate increase with atomic number, namely, a ∼3× enhancement between Al and Pb. At ϵ_p = 1 GeV, however, a ∼10× enhancement is obtained.

We now numerically characterize the neutron beams that could be produced through spallation reactions in the near future using the LMJ-PETAL and Apollon laser systems. Our methodology comprises two steps. First, we have employed the CALDER code46 to perform multidimensional PIC simulations of proton acceleration from laser-irradiated foil targets, under conditions relevant to the LMJ-PETAL and Apollon lasers, as detailed in Table II. Second, the proton distributions recorded in the PIC simulations have been used as input in 3D FLUKA Monte Carlo simulations,44 describing the proton transport and associated nuclear reactions through a secondary Pb convertor target. This procedure is sketched in Fig. 3. In this section, we present the results of the proton acceleration simulations for the LMJ-PETAL and Apollon cases.

Proton acceleration using LMJ-PETAL was investigated in a quasi-3D geometry with the CALDER-CIRC PIC code.47 In this code, the particles are advanced in 3D Cartesian space, but the fields and particle densities are computed in cylindrical coordinates (x, r, θ), using a reduced number of Fourier angular modes (only two angular modes were used here). This technique allows one to describe both the axisymmetric self-generated plasma fields (such as the TNSA field) and the non-axisymmetric (linearly polarized) laser fields at a much lower computational cost than in a fully 3D simulation. Based on experimental measurements,48 the PETAL laser beam was modeled as the sum of two superimposed, time-synchronized Gaussian laser waves, each of 610 fs FWHM duration, polarized along the y axis and propagating along the x axis. This was done in order to best describe the measured complex focal spot of PETAL, which is composed of a central, intense, spot, surrounded by low-intensity extended wings (see Fig. 1 in Ref. 48). The central component of the laser spot was modeled by a Gaussian having a FWHM width of 20 μm and a dimensionless field strength a_L = e E_L/m_e cω_L = 1.1. The wings were modeled by a Gaussian with a 130 μm FWHM width and a dimensionless field strength a_L = 1.3. The cumulative maximum laser intensity was of 8 × 10¹⁸ W cm⁻².48

The target consisted of a semitransparent Al preplasma and an overcritical plastic (CH) plasma (5 μm thick). Henceforth, nc=meϵ0ωL2/e2 (where ω_L is the laser frequency, m_e the electron mass, e the elementary mass, and ϵ₀ the vacuum permittivity) will denote the critical density beyond which the laser can no longer propagate in the plasma. For the λ_L = 1 μm wavelength of the PETAL laser pulse, one has n_c = 1.1 × 10²¹ cm⁻³. Based on hydrodynamic–radiative simulations,49 the electron density profile of the preplasma was taken to evolve as n_e(x) = 5n_c exp[5.45(x/150)^0.29], where the longitudinal position x is here expressed in micrometers. The minimum and maximum density values were set to 0.02n_c and 5n_c. The latter maximum density of 5n_c, which also characterized the uniform CH layer, was chosen such that it was high enough to accurately represent the absorption of the laser light and low enough to relax the constraints on the numerical discretization. The Al, C, and H ions were assumed to be fully ionized, and were initialized at a 10 eV temperature. The ionic species Al¹³⁺, C⁶⁺, H⁺, and the electrons were represented with 4, 4, 8, and 32 macroparticles per cell, respectively. Coulomb collisions between the plasma particles were neglected. The laser wave impinged normally onto the plasma. Absorbing boundary conditions were used for both particles and fields. The simulation domain, of dimensions L_x × L_r = 318 × 163 μm², was discretized with a Δx = Δr = 32 nm mesh size. The simulation was run during a ∼3.4 ps integration time.

With the above parameters, proton acceleration proceeds according to the standard TNSA mechanism.32,33 During the interaction, the laser wave propagates through the extended undercritical preplasma while driving the electrons to relativistic energies through various processes. These have been examined in detail by some of the present authors in a recent study,48 which revealed the importance of stochastic electron heating as a result of stimulated laser backscatter and laser filamentation in the preplasma.

After traversing the target, the laser-generated hot electrons form a negatively charged cloud at the backside. The associated electrostatic field then accelerates the plasma ions in the +x direction. Figure 4 displays the longitudinal phase space of the protons as measured at t = 3.4 ps. The proton distribution exhibits a linear shape typical of TNSA, extending over ∼100 μm longitudinally and (from analysis of the x − r proton density map) ∼300 μm transversely. The corresponding energy spectrum of the protons is plotted in Fig. 5. It has a decreasing exponential shape characteristic of TNSA, with a cutoff energy of ∼37 MeV.

To confirm the validity of our approach, we compare in Fig. 5 the proton distribution predicted by CALDER-CIRC with an experimental one, which was recorded by us at LMJ-PETAL. It was measured with the CRACC diagnostic, which uses a radiochromic film (RCF) stack as detector; the spectrum is shown here after angular integration over the whole surface of the films.50 In this shot, the PETAL laser irradiated a 7 μm-thick titanium foil with a pulse of 960 fs duration, 354 J energy, and ∼5.3 × 10¹⁸ W cm⁻² intensity, i.e., parameters relatively close to those of the simulation. Although the experimental proton spectrum is characterized by a harder slope and a larger cutoff energy (∼42 MeV) than the simulated spectrum, the two curves agree fairly well with each other. This result is particularly satisfactory, given the above simplifications made to the simulation in order to handle the spatiotemporal scales of the experiment. It should be noted that the experimental cutoff energy represents the last recorded point on the RCF stack, although there were more RCFs in the stack; such sharp cutoffs are common features of TNSA-accelerated proton beams.

Proton acceleration using Apollon was simulated in a 2D Cartesian geometry with the CALDER PIC code. The Apollon pulse was modeled as a 0.8 μm wavelength Gaussian electromagnetic wave linearly polarized along the (in-plane) y axis, with 20 fs FWHM duration and 5 μm FWHM transverse size. We considered two operating regimes of the Apollon laser, characterized by a peak intensity of 2 × 10²¹ W cm⁻² (0.6 PW regime) and 2 × 10²² W cm⁻² (6 PW regime), respectively; see Table II. Note that the temporal profile of the laser pulse is truncated at twice its FHWM duration ahead of the peak, i.e., the laser is switched on at that time. This starting point is chosen since it corresponds to an intensity that is 10⁻⁵ times the peak value. This ratio is close to what can be achieved using plasma mirrors as ultrafast switches, allowing one to virtually eliminate any ionizing light in the laser pulse prior to that level.51 In practice, however, modulations in the pulse spectrum36 could lead to the laser temporal profile deviating from the perfect Gaussian shape that we assume here. This may induce premature expansion of the target surface, thereby affecting the irradiation conditions at the pulse maximum and thus the overall performance of the target. This will have to be taken in consideration, using experimental data for the laser temporal profile, when carefully planning experimental campaigns at maximum power.

The target was a thin, solid-density CH foil with sharp gradients. It was assumed to be fully ionized (yielding an electron density n_e = 200n_c) and initialized at a temperature of 100 eV. Its thickness, l, was chosen based on the parametric simulation study of Ref. 29. The optimum foil thickness for RPA by femtosecond laser pulses was found to be l_opt ≃ 0.5a_L(n_c/n_e)λ_L. In the 0.6 PW (respectively, 6 PW) regime, corresponding to a_L = 30 (respectively, a_L = 96), we chose l = 64 nm (respectively, l = 192 nm), close to l_opt = 60 nm (respectively, l = 192 nm). Each of the plasma constituents (C⁶⁺, H⁺, electrons) was modeled by 100 macroparticles per cell. The laser pulse interacted at normal incidence with the target. Absorbing boundary conditions were adopted for the fields and particles. The simulation domain was set to L_x × L_y = 40 × 50 μm² at 0.6 PW (56 × 96 μm² at 6 PW), with a spatial resolution Δx = Δy = 3.2 nm.

The proton acceleration dynamics is illustrated by the longitudinal proton phase spaces shown at two successive times in Figs. 6 and 7 for the 0.6 and 6 PW irradiation cases, respectively. In each panel is also plotted (in blue) a lineout (along the laser axis) of the accelerating longitudinal electric field 〈E_x〉 (where 〈·〉 indicates averaging over the laser cycle). This field is scaled (in units of E₀ = 3.2 × 10¹² V m⁻¹) to fit within the p_x axis of the phase space.

At time t = −20 fs (here the time origin t = 0 is when the laser maximum reaches the target) in Figs. 6(a) and 7(a), proton acceleration originates from both RPA and TNSA,12,33 as demonstrated by the two 〈E_x〉 peaks at the front and rear sides of the target. At the front side, the electrons are pushed and compressed by the laser’s ponderomotive force. The ensuing charge separation generates an electrostatic field 〈E_x〉 ≃ 2E₀ at 0.6 PW (≃10E₀ at 6 PW), which in turn accelerates the front-side protons in the forward direction. These RPA protons have then reached a longitudinal momentum p_x/m_ic ≃ 0.015 (≃0.06 at 6 PW). Simultaneously, the backside protons have started expanding toward vacuum owing to TNSA triggered by the fast electrons. The associated (on-axis) electric field 〈E_x〉 ≃ 1.5E₀ (≃3E₀ at 6 PW), however, turns out to be weaker than the one induced by the radiation pressure (especially at 6 PW). At this early stage of the interaction, the maximum momentum of the TNSA protons is of p_x/m_ic ≃ 0.0011 (≃0.02 at 6 PW).

At t = +4 fs [Figs. 6(b) and 7(b)], the RPA protons have caught up with the TNSA protons, and the two previously observed field structures have merged into a single accelerating structure. Importantly, the expanding, lower-density plasma has then become transparent to the central part of the laser pulse. This causes the electrons to be further heated and the accelerator field strength to be boosted to 〈E_x〉 ≃ 7.5E₀ at 0.6 PW and 〈E_x〉 ≃ 25E₀ at 6 PW, which, in both cases, represents about a quarter of the laser field amplitude.

As a result of this sequence of processes, the protons eventually attain high cutoff energies (≃115 MeV at 0.6 PW and ≃660 MeV at 6 PW) as demonstrated by the energy spectra displayed in Fig. 8. Note that these spectra are recorded when the protons reach the right-hand side of the box: their spatial distribution then has a ∼20–30 μm transverse size, comparable to the distance traveled. The electrostatic field seen by the fastest protons should therefore be relatively well captured by our 2D simulation. We acknowledge, however, that the combination of a 2D geometry and laser polarized along the (in-plane) y-axis, likely leads to a significant overestimation of the electron heating, and hence of the accelerator field compared with a real-world 3D configuration.52 This, on top of the idealized temporal laser profile and target conditions (i.e., the neglect of preplasma formation) that we considered, is bound to degrade the performance in the experiments compared with our simulation results.

In the CALDER simulations, the properties (statistical weight, position, momentum, and time of arrival) of the macro-protons crossing a virtual detector plane near the right-hand side of the simulation box were all recorded. The resulting output files contained between about 10⁵ (PETAL) and 10⁷ (Apollon) macro-protons. The proton distributions obtained from the 2D Apollon simulations had to be post-processed in order to be used as input in the FLUKA 3D Monte Carlo code. For this purpose, they were converted into cylindrically symmetric distributions. Specifically, the position and momentum of each macro-proton were rotated around the x axis by a random azimuthal angle. Moreover, the transverse density profile of the proton distribution was interpreted as a radial density profile: the statistical weight of each macro-proton (a linear density in a 2D simulation) was therefore multiplied by its transverse radius to obtain a dimensionless quantity, corresponding to the number of physical protons represented by the macro-proton. These numbers are reflected in the spectra shown in Fig. 8 and are further detailed below.

The convertor target was taken to be a lead cylinder of fixed radius 5 cm and varying length (10 μm ≤ l ≤ 10 cm), located 0.5 cm behind the proton-generating target. To get good statistics on the simulated events, we carried out 1000 independent Monte Carlo simulations for each set of initial conditions. Only those neutrons crossing the target rear side were characterized.

The energy-angle spectra of the outgoing neutrons from a 0.3-mm-thick converter Pb target as predicted by FLUKA are shown for the LMJ-PETAL case in Fig. 9(a) and the 0.6–6 PW Apollon cases in Figs. 9(b) and 9(c). Overall, the neutron energy spectra show an exponentially decreasing shape up to a maximum energy close to that of the incident protons. This is more clearly seen in Fig. 12(b), obtained from angular integration of the spectra of Fig. 9. It is worth noting that the lower-energy part of the neutron distribution is essentially isotropic, while its higher-energy part is preferentially emitted in the initial direction of the proton beam.

Figure 10 plots, as a function of the Pb target length l, the fraction of kinetic energy lost by the primary protons through nuclear reactions (blue curves) and that transmitted through the target (green curves). In the PETAL case, the proton beam energy is wholly dissipated for l ≥ 0.1 cm. In the 0.6 PW Apollon case, this takes place for l ≥ 1 cm, while in the 6 PW case, the absorption is limited to ∼80% at l = 10 cm.

Figure 11(a) shows how the total number (per unit solid angle) of outgoing neutrons varies with l. Under the LMJ-PETAL conditions, the neutron number is seen to rise with l up to l ≃ 0.1 cm and to saturate at a ∼10⁸ n level in the range 0.1 cm ≤ l ≤ 3 cm, before dropping at larger l as a result of reabsorption. In the 0.6 and 6 PW Apollon cases, owing to higher proton energies, saturation occurs in thicker targets, namely, at ∼5 × 10⁸ n for 1 ≤ l ≤ 3 cm and at ∼10¹¹ n for 3 ≤ l ≤ 10 cm, respectively (the plateau observed at 6 PW may actually extend beyond the range of thicknesses considered here). These trends are consistent with the evolution of the dissipated proton energy as discussed above. Interestingly, the neutron yield is predicted to be quite similar for the LMJ-PETAL and 0.6 PW Apollon lasers. At first glance, this result may seem surprising, given that LMJ-PETAL generates about 50 times more fast protons than 0.6 PW Apollon (∼5 × 10¹² in a 2–40 MeV energy range vs ∼10¹¹ in a 10–120 MeV range), owing to its relatively large spot size and long pulse duration. Yet the lower yield of protons achieved at 0.6 PW Apollon is compensated for by their higher cutoff energy, which increases their neutron generation efficiency. Such an enhancement is even more dramatic using the 6 PW Apollon parameters, in which case a sharp rise in the proton energies is observed. The ∼5 × 10¹¹ protons then produced above 20 MeV are predicted to translate into a two orders of magnitude higher neutron yield.

Figure 11(b) plots the corresponding variations in the maximum neutron flux at the backside of the Pb target. Note that FLUKA takes account of the time of injection of each proton into the convertor, so that the temporal profile of the total neutron flux across a given surface can be computed. The maximum flux appears to be reached in l ≲ 100 μm targets and to drop in increasingly thick targets. Values in excess of ∼1 × 10²³ n cm⁻² s⁻¹ are expected at LMJ-PETAL and 0.6 PW Apollon, while a maximum flux as high as ∼4 × 10²⁴ n cm⁻² s⁻¹ is found with the 6 PW Apollon parameters. This trend results from the increase in duration and transverse size of the neutron distribution when the converter target is made thicker. These variations originate from the energy dispersion of the incident proton beam (which leads to an elongation of the proton beam, and therefore of the generated neutron beam) as well as from elastic scatterings of both the protons and neutrons throughout the target (which mainly account for the transverse size of the neutron source). The temporal dependence of the neutron flux is illustrated in Fig. 12(a) for LMJ-PETAL. It is observed that upon thickening the target from l = 10 μm to 1 cm, the neutron pulse is lengthened from ∼3 ps to ∼6 ns. The correlation between the neutron source duration and transverse size is clearly shown in Fig. 13. The three laser configurations give rise to similar behavior: very compact neutron sources, of a few picoseconds duration and ∼50–100 μm width only, are expected from l ≤ 100 μm Pb targets, which evolve into sources a few nanoseconds long and a centimeter wide when centimeter-thick Pb targets are employed.

We have assessed here the possibility of exploiting spallation reactions to generate high-flux neutron sources using petawatt-class lasers as the primary drivers. We have tested this scenario by combining two simulation codes: the first to simulate the proton generation by the laser and the second to simulate the neutron generation in a lead converter. Most notably, our results highlight the interest in using ultraintense femtosecond laser pulses to push the maximum proton energy. In particular, irradiation at the 6 PW level using Apollon, although yet untested experimentally, should produce protons beyond the 100 MeV threshold, i.e., well above the current performance of higher-energy picosecond lasers such as the 0.5 PW LMJ-PETAL system. Such high proton energies translate into much larger neutron multiplicity from the converter targets, and therefore allow multi-petawatt short-pulse lasers to make up for their lower proton output.

Regarding the quantitative accuracy of our study, satisfactory agreement was demonstrated in the LMJ-PETAL case between an experimental proton spectrum and that obtained from a quasi-3D PIC simulation. Although, as of now, such a comparison cannot be made in the Apollon setting with ultrathin targets, since the facility is still undergoing commissioning,53 we acknowledge the limitations of our 2D PIC simulations, and the fact that they may appreciably overestimate the proton cutoff energy (particularly in the 6 PW regime) as claimed by previous works.52 This leaves room for a refined (but much more computationally demanding) simulation study based on 3D simulations, and using more realistic (i.e., non-Gaussian) temporal laser profiles,36 to be carried out in the future.

To conclude, we note that the >10²³ n cm⁻² s⁻¹ peak neutron fluxes predicted by our numerical study could be appropriate to laboratory studies of r-process nucleosynthesis. In particular, the short duration (in the picosecond-to-nanosecond range) of the neutron source [shown in Figs. 12(a) and 13] is adequate to perform nucleosynthesis experiments, since the β-decay of the created isotopes resulting from multiple neutron absorption occurs over much longer timescales (longer than milliseconds).54 Beyond the practical achievement of such high instantaneous flux, which we will soon be able to verify using the Apollon laser facility53 and other multi-petawatt facilities like ELI, an evaluation of the overall amount of isotopes that could be produced per laser shot or in a cumulative mode (over several shots) needs to be conducted using presently available, i.e., theoretically estimated, cross sections; this is an ongoing task that will be the focus of a separate publication.

Acknowledgment. This work was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 787539). It was also supported by Grant No. ANR-17-CE30-0026-Pinnacle from the Agence Nationale de la Recherche. We acknowledge GENCI, France, for granting us access to HPC resources at TGCC/CCRT (Allocation No. A0010506129). S.N.C. acknowledges support from the Extreme Light Infrastructure Nuclear Physics (ELI-NP) Phase II, a project co-financed by the Romanian Government and the European Union through the European Regional Development Fund–the Competitiveness Operational Programme (1/07 July 2016, COP, ID 1334), and by the project ELI-RO-2020-23 funded by IFA (Romania). The PETAL laser was designed and constructed by CEA under the financial auspices of the Conseil Régional d’Aquitaine, the French Ministry of Research, and the European Union. The CRACC diagnostic was designed and commissioned on the LMJ-PETAL facility as a result of the PETAL+ project coordinated by University of Bordeaux and funded by the French Agence Nationale de la Recherche under Grant No. ANR-10-EQPX-42-01. The LMJ-PETAL experiment presented in this article was supported by the Association Lasers et Plasmas and by CEA. The diagnostics used in the experiment have been realized in the framework of the EquipEx PETAL+ via Contract No. ANR-10-EQPX-0048.