Numerical studies of a layered lithium-boron target for laser-driven aneutronic fusion reactions

This paper explores a novel target design for laser-driven, aneutronic, proton-boron and proton-lithium fusion reactions consisting of a stack of boron and lithium foils. In contrast to a homogeneous target, this multi-layer setup provides additional fusion channels in the diﬀerent materials. The composition of the layers is chosen in descending order of the fusion reactions’ thresholds, facilitating the fusion of protons that penetrate further into the material despite their energy losses due to electronic and nuclear stopping power. We employ a combination of Fluka simulations and additional numerical computations to evaluate thousands of target con-ﬁgurations. Four diﬀerent beam energy distributions are considered: two Gaussian distributions with 6 MeV and 10 MeV mean energies, respectively, a Maxwell-Boltzmann distribution and a power law distribution. We explore the production of energy in a range of layer thicknesses motivated by the proton ranges based on ionization losses. The conﬁguration which maximizes the produced energy for each beam type is reported. The production of fusion energy ranges from hundreds to thousands of millijoules for proton bunches of 10 15 having mean energies between 2-10 MeV.


Introduction
The advent of high-power lasers has led to the production of particle beams through laser-matter interactions.Proton acceleration with lasers has been well studied and different regimes have been identified for a range of laser parameters (e.g.Target Normal Sheath Acceleration [1, 2, 3], Collisionless shock [4,5], Coulomb explosion [6,7,8], in increasing order of the laser intensity).Furthermore, simple relations have been derived for Target Normal Sheath Acceleration which relates the laser pulse properties to the maximum energy of the protons [9].Such relation describes very well existing data which cover a range of proton energies from a few MeV up to ∼ 50 MeV.For a review of the different acceleration regimes and experimental evidence see e.g.Ref. [10].The energies achieved have motivated further applications such as ion therapy with lasers and fusion reactions with targets.Early observations by Belyaev et al. [11] showed that aneutronic fusion reactions can be initiated by high-power lasers.In this reference a 10-12 J laser was used to produce 1.5 ps pulses with a maximum intensity of 2 × 10 18 W/cm 2 , which was directed onto a 11 B + CH 2 target to yield about 1000 alpha particles per shot (a followup paper found this yield to have been underestimated by at least 100 times [12]).Subsequent works have studied different aspects of laser-acceleration-driven fusion both theoretically [13,14,15] and experimentally [16,17].Of particular interest are aneutronic fusion reactions (e.g.proton-boron reactions), where an unexpectedly high yield of secondary alphas has been reported and explained with fusion chain reactions [15,18,19].It also has been argued [20] that chain reactions between fusion products (i.e.alpha-particles) and 11 B nuclei in densely packed solid fusion targets could be the cause.Another explanation proposed involves collisions between the secondary alphas and beam protons, as pointed out in a very recent paper by Giuffrida et al. [21].The authors reported a record yield of alpha particles from fusion reactions of a laser-accelerated proton beam and boron nitride.The laser power in this experiment was rather low at about 10 16 W/cm 2 .However, the peak proton current was up to 2 A during a pulse duration of 10 ns (i.e. 10 11 alpha particles per pulse).In recent years a new type of small-scale fusion reactor was proposed [22] along with a possible new target design [23].The target suggested consists of alternating layers of lithium and boron with the aim of maximizing the number of fusion reactions.Hence, the following aneutronic fusion reactions were considered: (3) Furthermore, the optimal geometric form of the proton source that minimizes the spread of the proton beam was studied.It was found that the optimal shape for a solid-state proton source is a slab geometry.This paper expands on the ideas from [22,23], focusing on the optimization of the layered target.A target with three layers and four proton beams with different energy distributions are considered.The optimal energy production was sought by exploring a range of layer thicknesses motivated by the theoretical range of protons in the materials affected by ionization losses only.The method combined Fluka [24,25] simulations to account only for ionization losses and particle transport, and additional numerical computations to estimate the number of protons undergoing the aneutronic reactions (1), ( 2) and ( 3).The four different beam distributions are a Gaussian energy distribution with mean 6 MeV and with mean 10 MeV, a power law distribution and a thermal distribution.These are distributions expected from laser acceleration and for quasi mono-energetic beams, the spreads can be ∼ 10 % [26].Even narrower energy spreads of ∼ 1% have been demonstrated with proton energies reaching ∼ 20 MeV [27].It should be noted that our results are also valid when other means for proton acceleration are used, such as conventional accelerators.However, laser-driven acceleration seems a more promising method of acceleration due to its compactness and its increasing efficiency in bringing protons up to tens of MeV energies.This paper is organized as follows: Section II details the target composition and material properties employed.Section III discusses the methodology employed to combine Fluka simulations and numerical integration to obtain the reaction number and energy produced.This section also presents the simulation parameters and the cross-section data employed.Section IV presents the results with the combination of thicknesses of the target layers which yield the most energy for each of the proton beams.

II. Target design
The Li-B target under consideration consists of layers of 6 Li , 7 Li , and 11 B .Each of the layers has a constant particle density and we assume that only the proton beam is hitting the target (i.e. the laser beam is decoupled from the protons before they reach the target).
Protons moving through the target will progressively lose energy through electron stopping in addition to undergoing the desired fusion reactions.This allows controlling the mean energy of the protons at each depth in the target by choosing the thickness of each layer appropriately.For example, employing Bethe's formula for the energy losses due to ionization, the range of 10 MeV proton in lithium is around 3 mm, while for 1 MeV it is 85 µm.
Thus, a suitable choice of layer thickness would allow limiting the average energy loss of the proton to any desired value.In our case, where different materials are involved, we can estimate the required thickness by using that the stopping power scales with ρZ A where ρ is the density and Z, A the atomic and mass number of the material respectively.Since the range is the integral of the inverse of the stopping power, its scaling would then be A ρZ .Hence, the range in boron for a 10 MeV proton can be estimated from the range in lithium by multiplying it by which yields 759 µm employing the values in Table 1.Thus, to fully stop 10 MeV protons it would be necessary a layer of Li of at least 3 mm of thickness, or for example 1 mm lithium layer and the equivalent of 2 mm of lithium in boron, that is 0.4 mm.The order of the layers is presented in Table 1 in addition to the densities and the proton energies where the cross-section is maximal in each case.The cross sections employed and sources are suffer enough losses to fall in the energy range covering the maximum cross-section of the corresponding reactions.This target design is expected to increase the number of fusion reactions in contrast with a target with only one layer where the fusion number decreases as soon as the proton energies are lower than the value where the crosssection peaks.

III. Methodology
The energy distribution of protons N (x k , E) is affected by a differential thickness x k of material k as expressed by the following differential equation: where Σ k (E) is the macroscopic cross section for reactions and µ k (E) fraction of protons of energy E that losing an amount of energy of at least dE.
Hence, the term Σ k (E)N (x k , E) represents the number of protons that undergo nuclear reactions, and the term µ k (E)N (x k , E) represents the number of protons losing energy and unable to react at energy E. These quantities are dependent on both the thickness x k and the proton energy E. The macroscopic cross section for each material is com- with the densities ρ k and atomic masses m k given in Table 1.The cross sections σ k (E) were obtained and discretized for folding with the proton distribution as explained in Appendix VI..
The form of µ k (E) is more complex and requires the inclusion of different effects.For example, protons undergo nuclear elastic collisions, which causes non-straight trajectories.This has the effect of spreading the initial energy distribution of protons as they penetrate the target since different paths lead to different losses.Such an effect depends non-linearly on the layer thickness and it is important in predicting the entry energy for subsequent layers.
To resolve these difficulties, we have employed the code Fluka [24,25], which models the diffusion of protons through the target including the elastic collisions and electronic and nuclear stopping powers very well.However, to better control the fusion reactions for each material, the computation is done outside of Fluka using the cross sections for each of the reactions considered that are available in the literature (see Appendix VI.).
The number of interactions in a material Γ k requires integrating within the range of energies considered (E min , E max ) and throughout the layer's thickness d k .Having calculated the number of interactions per material, the total energy produced is given by where the energies produced in each material Q k are shown in Table 2.

Numerical integration using Fluka
In order to integrate Eq. 4 we employed a numerical procedure.The method consists of recursively setting thicknesses for the target layers, simulating the proton diffusion with Fluka, which yields the proton energy distributions after each layer.The number of reactions is thus computed by estimating the proton energy distribution inside each layer and convolving with the energy-dependent cross sections.
To numerically integrate Eq. 4 is convenient to express it in finite differences N (x k , E) can be integrated choosing sufficiently small and equal steps of ∆x k = ∆x 0 , leading to a monotonically increasing sequence of thicknesses After a layer, the redistribution of protons due to energy losses are modeled with Fluka and reflected on the exiting distribution {M l k }.However, the fusion reactions considered here are not included in the Fluka simulations, and therefore need to be accounted for in the first term of Eq. 7. The total energy produced was obtained for each combination of layer thicknesses, represented in Figs.7-9.

Parameters of the Fluka simulation
The target geometry employed is shown in Fig. 1.The simulation volume is represented by a cylinder of 10 mm radius and variable length determined by the sum of the thickness of all layers.The proton emitting surface is a circle with a 50 µm radius positioned at 10 mm from the face of the target, having an impact diameter of ∼ 1 mm on the surface of the target.The beam directions follow a flat distribution centered along the symmetry axis of the target and a spread of 100 mrad (from -50 mrad to 50 mrad).
The densities shown in Table 1 are given as parameters to Fluka.In Fluka the materials are specified by providing the composition and the density.The material of each layer was set to a pure composition of the corresponding isotope and density as tabulated.
In addition to varying the target thicknesses, four different beam types were considered (see Fig. 2): a) a Gaussian-distributed beam with mean energy 6 MeV and a full width at half maximum of 0.5 MeV All additional products from proton interactions can be neglected in this study, and we have suppressed their production by using appropriate settings in Fluka.In some cases, where a high energy tail of the incoming beam is present, additional reactions may occur that lead to products like alpha particles with sufficient energy to produce additional fusion reactions.Such cases will also not be considered in this work, and consequently, they have not been scored in the simulations.

IV. Results
Tab. 3 presents for each of the beam scenarios, the values of the parameters computed for the combination of layers corresponding to the maximum Q tot .Figs. 3, 4, 5, 6 present the proton energy distributions crossing each layer boundary for the optimal target configuration.It should be noted that the proton number is given in MeV −1 and hence, there are fewer protons per MeV when the energy of the incoming beam is higher.The cases that produce the most energy are the Gaussian beams both around 3 J which is remarkable given that in the 6 MeV case the proton energies are 60% of the 10 MeV case.This energy differ- Gauss  Table 3: Values of the parameters for the target thicknesses that maximize the total energy, separated by beam case and by layer.The columns from the left, are the beam case, the layer order, the layer thickness, the number of fusion reactions per layer, the energy produced per layer and the total energy produced.ence is reflected in the much thicker layers for the 10 MeV case and the higher mean energy of protons entering the second layer.In both cases, most of the energy is generated in the first layer, but the energy losses differ as Figs 3 and 4 show.The protons entering the second layer (orange line) have a mean energy of 2.8 MeV in the 6 MeV case and 5 MeV in the 10 MeV case, which indicates a ∼ 50% mean relative energy loss in both cases, but the absolute is almost double in the latter.This indicates that for the 6 MeV case, a comparable amount of fusion energy can be generated with much fewer losses due to ionization.
The thermal beam distribution has a mean energy comparable to that of the Gaussian 6 MeV case, however, the maximal energy produced is almost three times smaller.This is in part caused by the broader distribution of the initial beam, with about half the number of protons above the mean value.
The lower energy protons are largely stopped in the first layer without interacting given that the fusion cross section for the first layer is much smaller below 1 MeV (see Figs. 11 and 14).This also produces the effect of increasing the mean proton energy as they penetrate further into the target.Indeed, Fig. 5 shows how the initial distribution's 6.1 MeV mean energy, changes to ∼7 MeV entering the second layer, and to ∼ 8 MeV as they enter the last layer.This effect conspires against the desired lowering of the average energy to match the descending fusion optimal energies of the target design.Additionally, the thermal distribution extends to higher energies than the Gaussian 10 MeV which explains the much larger thicknesses found, 2-3 times compared to the 10 MeV Gaussian case.This is also contributing to the increase of protons' mean energy with depth.An inversion of the ordering of the layers could be a solution for the thermal case, but this was not explored in this work.
In contrast, the power-law case has a mean proton energy of 2.3 MeV and a maximum energy of 10 MeV and the total thicknesses of the first two layers are similar to the 10 MeV case (within 6 %), but the energy produced is almost six times smaller.
The presence of 10 MeV protons in the power-law (around 50 times that of the Gaussian 10 MeV case) determines the need for equivalent thickness to reduce their energy, but this also causes most of the protons to be stopped in the first layer.Fig. 6 shows how only a small fraction of protons enters the last layer, leading to a negligible fusion number.In this case, the increase of the mean energy as the protons reach further into the target occurs as well, given the large fraction of protons at lower energies, and it is possible that a different layer ordering could produce more energy than the optimal configuration found here.The slices are all intercepting the point of the phase space with the maximal energy produced, which is represented with a white star.This means that in each section, the value of the missing parameter is that of the maximal energy configuration.In all phase space distribution sections L1-L2 show a relation between the thicknesses, with smaller values of one thickness correlated with higher values of the other to produce a banana-shaped region of higher energy values.In Fig. 7 this is more clearly shown, where is visible in all panels.In Fig. 9 this relation is not appreciable due to the strong dependence of the produced energy on the thickness of layers 1 and 2, causing the vertical lines in sections L1-L3 and L2-L3.Such a relation is to be expected since the protons have similar energy losses in both lithium layers.In other words, making the first layer shorter and increasing the second one by a similar value would lead to comparable losses with penetration in the target.In addition, the smaller mean fusion cross section in layer 1 compared to that in layer 2 is compensated by the larger value in produced energy per reaction.In Fig. 7 the strong dependence on the thickness of layer 1 leads to sections L2-L3 practically independent of the thickness of the last layer.It is notable that in the Gaussian 6 MeV case, other regions of significant energy production are visible (see Fig. 7).In section L1-L3 (layer 2 thickness is 350 µm) very large values of energy can be achieved for layer 1 in the range 270-360 µm and layer 3 above 60 µm.Similarly, in sections L2-L3 (layer 1 thickness is 750 µm) even larger values of energy are reached for layer 2 of ∼200 µm and layer 3 also above 60 µm.However, the optimal combination favors values of L1 and L2 such that the energy produced in layer 3 is negligible.

V. Conclusions
In this paper, we showed simulations for highenergy proton beams traversing a hybrid fusion target that is comprised of lithium and boron isotopes.The computer calculations were done for proton beams with different types of energy distribution functions.The functions, which were used for the numerical simulations were: a distribution following a power law, a thermal beam and a beam with a Gaussian energy distribution function.The layer thicknesses that maximize the energy production in each case were reported.Energy productions of  0.4-3 J for each proton beam were found (10 15 protons per bunch, 2.3-10 MeV mean energies).The Gaussian beam cases (quasi-monochromatic) produced the most energy, but the Gaussian 6 MeV case would be preferable to reach the same energy production with fewer losses and lower proton energies.The thermal and power-law beams although yielding lower energies show indications  that larger production could be achieved with a different ordering of the layers than for the Gaussian beams.The layer thicknesses explored range from millimeters to tens of microns, meaning the targets are very compact with less than 6 mm total thickness.In all cases, the production of energy is mostly occurring in the first layer, with a strong reduction in the second layer (at least 24 times from the first) and a negligible contribution in the third layer.This seems to indicate that the Li-B targets don't seem to benefit as much from the multiple layers, and hence, don't perform significantly better than a onelayer target.Nevertheless, the phase space of thicknesses explored here indicates that other regions might be at least as well-performing, and could lead to better production.Furthermore, the behavior of the thermal and power-law case, with a large number of protons in lower energies, prompts a different ordering of the layers and would need examination in future works.This work, on the other hand, demonstrates that each beam shape has unique features as its development within the target has a non-trivial relation with its energy distribution, layers' thicknesses and even ordering.Another point, which has to be considered, is the conversion efficiency from the laser beam to the protons.Recent experiments by Brenner et al. [29] have yielded conversion efficiencies of about 15 percent.However, these numbers were obtained for proton bunches of only 10 13 protons, with average kinetic energies of 5 to 30 MeV.Hence, the assumption used in this paper seems to be realistic and achievable currently, and given the current rate of progress, better efficiencies are to be expected in the near future.

VI. Appendix
The cross sections for the reactions used in this paper were obtained by interpolating available experimental data.Figs.11-13 present the experimental values employed and the resulting interpolations.The interpolation values were subsequently discretized into the energy bins employed in the simulations.This facilitated computing the number of reactions as the sum of the bin-by-bin products between the discretized cross sections and the proton energy distribution.Fig. 14 shows the discretization of the continuous interpolations.

Figure 1 :
Figure 1: Schematic representation of the target and proton beam.The geometrical dimensions are indicated with arrows and different colors.The protons start at the red circle and impinge on the target within the gray circle.The angular spread of the protons is shown by the conical opening on the leftmost side of the image, and the conical edges represent the maximal angular deviation of the protons.

17.346 93 & 65 2 6→ 4 3 Table 2 :
Li+p He+ He Fusion reactions triggered by the protons in different materials.The kinetic energy of the products is given in brackets in MeV units, and based on the masses form Ref.[28].The corresponding energy output Q k , cross-section peaks σ max k and the proton energy where the cross-section peaks E kin k are also given.
x l k = l∆x 0 with l ∈ [0, N ].The evaluation of the first term in the right-hand side of Eq. 7 is independent of the thickness except on the factor N (x k , E), which needs to be evaluated for each thickness x l k , ı.e. a sequence {N l k} : N l k (E) = N (x l k , E) is needed.The second term in the right-hand side of Eq. 7 will be evaluated as a whole leading to a se-quence of {M l k } : M l k = µ k (E)N (x l k , E)∆x k .The determination of the mentioned sequences {N l k } and {M l k } proceeded by using Fluka.The target was implemented with a fixed surface and varying the thicknesses of each layer.The proton distribution after each layer was recorded and used to obtain the sequences {N l k } and {M l k }.The values {N l k } are directly obtained as the proton distribution entering a layer.The values {M l k } are directly obtained as the proton distribution exiting a layer.

Figure 2 :
Figure 2: Comparison between the different beam profiles employed in terms of the proton energy distributions.The Gaussian cases appear as only 1-2 bins due to the smaller energy spread compared to the bin sizes employed.Compared to the gaussian cases, the power-law and thermal beams have fewer protons of energies above 6 MeV, with the former having substantially less than the latter.

Figure 3 :
Figure 3: The number of protons crossing each layer as a function of energy for a monochromatic beam of 6 MeV.The different colors represent the protons entering each layer, which are indicated in the legend (e.g. in Li7 indicates the protons entering the 7 Li layer.)

Figure 4 :
Figure 4: Analogous to Fig. 3 but with a beam energy of 10 MeV.

Figure 5 :
Figure 5: Analogous to Fig. 3 but in the case of a thermal distribution of the proton beam.

Fig . 7 , 8 , 9 ,
Fig .7,8, 9, 10 represent the total energy produced for different configurations of the target.Each panel is a 2D slice of the 3D phase space constituted by the thicknesses of the three layers.The color scale represents the total energy produced by fusions through all the targets.The slices are all intercepting the point of the phase space with the maximal energy produced, which is represented with a white star.This means that in each section, the value of the missing parameter is that of the maximal energy configuration.In all phase space distribution sections L1-L2 show a relation between the thicknesses, with smaller values of one thickness correlated with higher values of the other to produce a banana-shaped region of higher energy values.In Fig.7this is more clearly shown, where is visible in all panels.In Fig.9this relation is not

Figure 6 :
Figure 6: Analogous to Fig. 3 but in the case of a powerlaw distributed beam.

Figure 7 :
Figure 7: Dependency of Qtot versus the widths of different layers for the Gaussian 6 MeV case.Each panel shows Qtot values averaged over the third layer's width.The star represents the layer thicknesses where the maximum Qtot is achieved.

Figure 8 :
Figure 8: Dependency of Qtot versus the widths of different layers for the Gaussian 10 MeV case.Each panel shows Qtot values averaged over the third layer's thicknesses.The star represents the layer width where the maximum Qtot is achieved.

Figure 9 :
Figure 9: Dependency of Qtot versus the widths of different layers for the thermal beam case.Each panel shows Qtot values averaged over the third layer's thicknesses.The star represents the layer width where the maximum Qtot is achieved.

Figure 10 :
Figure 10: Dependency of Qtot versus the widths of different layers for the power law case.Each panel shows Qtot values averaged over the third layer's thicknesses.The star represents the layer thicknesses where the maximum Qtot is achieved.

Figure 11 :
Figure 11: Experimental values and interpolation of cross section for the reaction in the 7 Li layer.Data reported in references[30,31,32].

Figure 12 :
Figure 12: Experimental values and interpolation of cross section for the reaction in the 6 Li layer.Data reported in references[33,34].

Figure 13 :
Figure 13: Experimental values and interpolation of cross section for the reaction in the 11 B layer.Data reported in reference [35].

Figure 14 :
Figure 14: Comparison of interpolated cross section and the discretization values in the chosen logarithmic grid.

Table 1 :
The composition, order and densities of the layers in the target.The order corresponds to the sequence in which the protons cross each layer.discussed in the Appendix VI..A bunch of laseraccelerated protons impinging onto the Li-B target would first encounter the7Li , where a fraction of the protons undergo fusion reactions and the rest lose energy continuously until reaching the following layer.The surviving protons that reach the subsequent layers with materials6Li and 11 B should