Study on deuteron formation mechanism in nucleon-induced reactions

NUCLEAR PHYSICS AND INTERDISCIPLINARY RESEARCH

Study on deuteron formation mechanism in nucleon-induced reactions

Ya-Jun He，

Chen-Chen Guo ，

Jun Su，

Long Zhu，

Zhen-Dong An

Nuclear Science and Techniques

Vol.31, No.8

Article number 84

Published in print 01 Aug 2020

Available online 08 Aug 2020

DOI：10.1007/s41365-020-00788-5

42202

The mechanism of deuteron formation in neutron-induced reactions is studied within the framework of the isospin-dependent quantum molecular dynamics (IQMD) model, using the GEMINI code. The influence of the n+p→d reaction channel is investigated by analyzing the deuteron production cross-sections in the neutron-induced reactions ¹²C(n,d), ¹⁶O(n,d), and ²⁸Si(n,d), with incident energies of 20-100 MeV. By including the n+p→d reaction channel when modeling the collision, the deuteron production cross-sections increase, optimizing the cross-section results and bringing them closer to the experimental data values. This indicates that the n+p→d reaction channel is an important mechanism for enhancing deuteron production.

Cluster mechanismDeuteron formation cross-sectionNucleon-induced reactions

1 Introduction

Spallation reactions, where high-energy, light particles collide with a heavy target, resulting in a lighter remnant nucleus and the ejection of numerous light particles, play an important role in a wide range of applications [1]. Ever since cosmic ray spallation reactions were studied in the 1930s [2], they have attracted much attention from researchers in applied and fundamental fields [3]. This interest only increased with the proposition of an optimum neutron resource, as well as the application of nuclear waste transmutation in an accelerator-driven system. Over the years, there have been extensive developments of applications entailing spallation reactions, including material physics [4], nuclear waste disposal [5-7], particle physics and nuclear physics [8], rare isotope production, accelerator radiation protection [9], cancer hadron therapy [11], and cosmic rays in the atmosphere [10].

A nucleon-induced spallation reaction is a two-stage process. First, the hard nucleon-nucleon (NN) collisions take place, dissipating the incident energy and resulting in a nuclei with high excitation energy. In the second stage, de-excitation takes place via evaporation and fission processes [12, 13]. In addition to fission and fusion, spallation is an effective process of releasing neutrons from nuclei; unlike the other two reactions, it is an endothermic process. To sustain a spallation reaction, a beam of high-energy particles, usually protons, must be supplied and focused on a heavy target. Spallation can be an important source of neutrons whose flux can easily be controlled by the drive beam [14]. However, it is hard to obtain the large quantity of reaction data necessary, in consideration of the low feasibility and high costs. Therefore, theoretical studies are imperative to improve our understanding of these reactions. There are many spallation models and codes currently used, including the particle and heavy-ion transport code system (PHITS) [15], the Boltzmann-Uehling-Uhlenbeck (BUU) transport model [16], the quantum molecular dynamics (QMD) model [17-19], the cascade model (JAM1.0) [20], the intranuclear cascade model (INCL)4 [21], ISABEL [22, 23], and the Monte Carlo code GEANT4 [24, 25]. These models and codes are used in conjunction with various statistical decay models, such as the generalized evaporation model (GEM) [26, 27], GEMINI [28], and the abrasion-ablation model (ABLA) [29]. More recently, a Bayesian neural network (BNN) method has been proposed for spallation reactions [30]. Nonetheless, the discrepancies between observables in theoretical calculations and experimental data are still relatively large, proving that the theoretical models are not well founded. One of the problems the theoretical models struggle with is being able to reproduce the production cross sections of the light cluster well. Attempts have been made to study the formation mechanism of clusters, such as deuteron, to improve upon the modeling. A surface coalescence mechanism has been used in INCL, Jaeri QMD (JQMD), and improved QMD (ImQMD) models to improve the description of the light complex particles emitted in various nucleon-induced reactions [21, 31, 32]. The pre-equilibrium light complex particle emission has been described reasonably well [18, 19]. An investigation into the formation mechanism of deuteron clusters is conducive to solving the current controversial clustering problem.

In this work, we study neutron-induced reactions by using the isospin-dependent QMD model (IQMD) model. The GEMINI model is applied to simulate the prefragment decays. The version of the IQMD model used is IQMD-BNU. The framework of the IQMD+GEMINI model has been applied to study nuclear reactions in an energy region from 10 MeV to GeV [33-42]. By using these models, we attempt to study the formation mechanism of deuteron clusters. The reaction channel of n+p→d is added into the model. In Sect. 2, a brief introduction of the IQMD model and cluster mechanism is made. Sect. 3 describes and analyzes the subsequent calculation results, with some brief discussion about the results. In Sect. 4, we summarize our findings, provide some concluding remarks, and discuss potential future directions that our work can take.

2 model

2.1 IQMD-BNU

The IQMD model is based on the same principles as the QMD model [43-45]. By considering the evolution of nucleons in mean field potentials and the binary NN collisions, the IQMD model has been widely and successfully used for the analysis of heavy-ion collisions (HICs) at incident energies below 2 GeV/nucleon.

The IQMD-BNU model was introduced and compared with other model versions in a transport-code-comparison project [46, 47]. In a QMD-type model [48], each nucleon is represented by a Gaussian wave packet [49]

ψ_{i} (r_{i}) = \frac{1}{{(2 π σ_{r}^{2})}^{3 / 4}} e^{- \frac{{(r_{i} - r_{i 0})}^{2}}{2 σ_{r}^{2}} + i (r_{i} - r_{i 0}) \cdot p_{i 0} / ℏ},

(1)

where σr and r_i0 are the width and centroid of the wavepacket, respectively. Usually, $σ_{r}^{2}$ is set at 1.0 fm². The total N-body wavefunction is considered to be a direct product of these coherent states,

ψ (r_{1}, \dots, r_{N}) = ϕ_{k_{1}} (r_{1}) ϕ_{k_{2}} (r_{2}) \dots ϕ_{k_{N}} (r_{N}),

(2)

ϕ_{k_{i}} (r_{i}) = \frac{1}{{(2 π σ_{r}^{2})}^{3 / 4}} \cdot \exp [- \frac{{(r_{i} - r_{i 0})}^{2}}{2 σ_{r}^{2}} + \frac{i p_{i 0} \cdot (r_{i} - r_{i 0})}{ℏ}],

(3)

where $ϕ_{k_{i}}$ is the wavefunction of the ith particle at a state ki (pi=ħki). $ϕ_{k_{i}}$ is chosen as the Gaussian wavepacket to avoid the negative values of the phase-space distribution (fN). By the Wigner transformation of the N-body wavefunction, the phase-space density is given by

\begin{array}{l} f_{N} (r_{1}, \dots, r_{N}; p_{1}, \dots, p_{N}) = \prod_{i = 1}^{N} f (r_{i}, p_{i}) \\ = & \prod_{i = 1}^{N} \frac{1}{{(π ℏ)}^{3}} \exp [- \frac{{(r_{i} - r_{i 0})}^{2}}{2 σ_{r}^{2}} - \frac{{(p_{i} - p_{i 0})}^{2}}{2 σ_{p}^{2}}], \end{array}

(4)

where, r_i0=〈r_i〉and p_i0=〈p_i〉 are the centroids of wavepackets in coordinate and momentum space, respectively.

The time evolution of the nucleons in the system is considered under self-consistent conditions. The mean field is governed by the Hamiltonian equations of motion,

\begin{array}{l} {\dot{r}}_{i} = \nabla_{p_{i}} H, \\ {\dot{p}}_{i} = - \nabla_{r_{i}} H, \end{array}

(5)

where the Hamiltonian, H, constitutes a kinetic energy component, T, a Coulombic potential energy term, and a nuclear potential energy term,

H = T + U_{Coul} + \int V_{n} (ρ, δ) d r .

(6)

The density can be calculated by taking the integral of the phase-space density [Eq. (4)] over momentum space. The nuclear potential V_n includes a two- and three-body Skyrme interaction term, as well as a symmetry potential term. The nuclear potential energy density of the asymmetric nuclear matter, with density ρ and asymmetry δ, is given by

\begin{array}{l} V_{n} (ρ, δ) = & \frac{α}{2} \frac{ρ^{2}}{ρ_{0}} + \frac{β}{γ + 1} \frac{ρ^{γ + 1}}{ρ_{0}^{γ}} + \frac{C_{sp}}{2} {(\frac{ρ}{ρ_{0}})}^{γ_{i}} ρ δ^{2}, \end{array}

(7)

where ρ₀ = 0.16 fm^-3 is the normal density. The constant parameters are taken to have values of α = -356 MeV, β = 303 MeV, γ = 7/6, C_sp = 38.06 MeV, and γi = 0.75. These values provide a compressibility of 200 MeV at saturation density for isospin symmetric nuclear matter.

The binary NN collisions are also included in our model. Specifically, the NN elastic scattering cross-section can be expressed by the following formula,

{(\frac{d σ}{d Ω})}_{i} = σ_{i}^{free} f_{i}^{angl} f_{i}^{med},

(8)

where σ^free is the elastic scattering cross-section of NN collisions in free space, f^angl is the angular distribution correction, and f^med is the in-medium correction. The isospin-dependent parametrization of σ^free and f^angl are taken from Ref. [50]. The in-medium factor for elastic scatterings is written as [51]

\begin{array}{l} f_{el}^{med} = σ_{0} / σ^{free} \tanh (σ^{free} / σ_{0}), \\ σ_{0} = 0.85 ρ^{- 2 / 3} . \end{array}

(9)

To compensate for the fermionic features, the phase-space density constraint (PSDC) method, taken from the constraint molecular dynamics (CoMD) model [52], is used. The phase-space occupation probability, ${\bar{f}}_{i}$ , is calculated by performing integration on an hypercube of volume h³ in the phase space centered around the ith nucleon at each time step. The effect of Pauli blocking is considered based on the phase-space occupation probability ${\bar{f}}_{i}$ ; only NN scattering producing final states with ${\bar{f}}_{i}$ < 1 is allowed.

2.2 GEMINI

The statistical code GEMINI, proposed by R. J. Charity in the 1980s [53], can be utilized to model the decay of a compound nucleus in fusion reactions and excited fragments in HICs. The GEMINI code is applied to the process of de-excitation of hot fragments [54, 55].

The IQMD code will stop running whenever the excitation energies of the two heaviest pre-fragments are less than a specified value, E_stop = 3 MeV/A. The charge number, mass number, and excitation energy of each pre-fragment are output by the IQMD code. If the excitation energy is greater than zero, the modeling of the sequential decay of the pre-fragment will be performed by the statistical code GEMINI. This statistical code not only allows light-particle evaporation and symmetric fission, but all possible binary decay modes. Subsequently, a Monte Carlo technique is used to track all the decay chains until the resulting products are unable to decay further. Part of the necessary decay width data comes from the Hauser-Feshbach formalism of the evaporation of light particles.

In summary, if the excited energy of the pre-fragments is lower than 3 MeV/A, the statistical decay model GEMINI is preferred. We start the simulation of GEMINI after a simulation time of 1000 fm/c in the IQMD model, also using the PSDC method. We assume that the evolution is stable enough until 1000 fm/c. For further details on GEMINI, refer to Ref. [53].

2.3 IQMD+GEMINI MODEL UPDATES

In the standard version of the IQMD model, only the elastic and inelastic NN collisions are considered. Taking the neutron-proton collision as an example, the scattering is performed if two requirements are satisfied. First, the nucleons must pass the point of closest approach in the next time step. Second, the collision parameter of the neutron-proton collision is less than the value of $\sqrt{f_{np}^{med} σ_{np}^{free} / π}$ , where $σ_{np}^{free}$ and $f_{np}^{med}$ are the cross-section in free space and the in-medium factor of the neutron-proton scattering, respectively.

In this work, the decay channel n+p→d is added, with the goal of studying the deuteron formation. The neutron-proton collision is performed according to the total cross-section,

σ_{t} = f_{np}^{med} σ_{np}^{free} + f_{d}^{med} σ_{d} .

(10)

where σ_d is the cross-section of the n+p→d reaction, and $f_{d}^{med}$ is the in-medium factor. The cross-section of the n+p→d reaction has been measured in various experiments, with the incident energy ranging from 0.01 to 100 MeV [56-65]. The experimental data is depicted as circles in Fig. 1. By fitting the data, an empirical formula of σ_d as a function of the incident energy is obtained,

Fig. 1.

Cross-section of the n+p→d reaction, in relation to the incident energy, E_lab. The circles show the experimental data taken from Ref. [56-65]. The solid curves are a result of the calculations from the empirical formula, Eq. (12).

σ_{d} (E) = {\begin{array}{l} - 0.14 \ln (E) - 0.25 \\ 0.035 \\ - 0.011 \ln (E) + 0.061 \end{array} \begin{array}{l} E < 0.14 M e V \\ 0.14 ⩽ E ⩽ 10 M e V \\ E > 10 M e V \end{array},

(11)

where the unit of the cross section is mb. The empirical value of the in-medium factor is,

f_{d}^{med} = {\begin{array}{l} 1 \\ 0 \end{array} \begin{array}{l} ρ < 0.2 ρ_{0} \\ ρ ⩾ 0.2 ρ_{0} \end{array},

(12)

where ρ is the density in the center of the colliding neutron and proton, and ρ₀ is the normal density. This in-medium factor value means that the deuteron can only be produced in the low-density region, i.e., the deuteron is generated when a neutron and a proton are close enough, and the medium density is less than 0.2ρ₀. When a deuteron evolves to a medium density greater than 0.2ρ₀, it will break up into a neutron and a proton.

It should be noted that the gamma rays were ignored in the n+p→d+γ decay channel in this work, since the beam energy is much larger than the energy. As a weak composite system, the propagation of the deuterons are approximately treated the same as the propagation of a neutron and proton under a nucleonic mean field potential. For the deuteron break-up mechanism, we adopt a phenomenological method that causes the deuteron to decay when its medium density is less than 0.2ρ₀. A long-term deuteron production mechanism in the QMD model is still an open question, which requires further study.

3 Results and discussion

The excited energy and charge distribution of pre-fragments for the n+²⁸Si reaction at 60 MeV are displayed in Fig. 2. The percentage of pre-fragments is represented by p. The excitation energy of the fragments can be calculated using,

Fig. 2.

The distribution of (a) the excited energy, E^*, and (b) the charge (atomic number, Z) of pre-fragments for n+²⁸Si at an incident energy of 60 MeV, calculated within the IQMD model. The percentage of pre-fragments in each distribution is represented by p.

E^{*} = \frac{\sum_{i} U_{i} + \sum_{i} \frac{{(p_{i} - p_{f})}^{2}}{2 m} - B (Z_{f}, A_{f})}{A_{f}} .

(13)

Here, Ui and p_i are the single-particle potential and momentum of the ith nucleon in the rest frame of the fragment, while p_f, Z_f, and A_f are the average momentum per nucleon, charge number, and mass number of the fragment, respectively. The variable B(Z_f, A_f) is the binding energy of a nucleus with charge number Z_f and mass number A_f. The summation is present to account for all the nucleons belonging to the same fragment. From Fig. 2 (a), one can see that the excited energy distribution of the pre-fragments has a center peak at around 1 MeV/A. We calculate that the excitation energy required to form ²⁹Si is 2.28 MeV; at this point the energy of the neutron is all transferred to ²⁹Si. While only 1 MeV/A is the approximate value of the center peak of the excited energy distribution; this is a result of peripheral collisions, where neutrons only transfer part of the energy to ²⁸Si. The interrupt time in our work is 1000 fm/c. Fig. 2 (b) shows the charge distribution of pre-fragments. More than 95% of pre-fragments have an atomic number of 12, 13, or 14. This means that the predominant reaction mechanisms of light nuclei with neutrons are scattering and nucleon trapping.

To study the formation of a deuteron cluster in detail, the formation time, collision energy, density at its position of formation, and excitation energy are shown in Fig. 3. The highest percentage of deuteron clusters are formed in the time frame 0-100 fm/c (Fig. 3(a)). The deuteron cluster is formed at this time because the incoming neutron captures a proton in the target, in the low-density region. From 100 to 1000 fm/c, the deuteron is formed from the excited nucleus in the low-density area during the reaction. It should be noted that the neutrons that form deuterons are not necessarily incident neutrons. Fig. 3(b) shows the relative kinetic energy of the neutron and proton in the excitation system when the deuteron cluster is formed. As the relative kinetic energy increases from 0 to 10 MeV, the probability of deuteron generation increases gradually. After 10 MeV, this probability gets smaller. At 60 MeV, the probability that a deuteron is formed is lower than one percent. This is because the cross-section of deuterons decreased gradually above 10 MeV, as shown in Fig. 1, in turn causing the probability of deuteron production to also decrease gradually. Fig. 3(c) shows the ratio between density and normal density, at the position where the deuteron cluster formed. The deuteron cluster prefers to form in a density ratio of 0.3–0.5. When the ratio is greater than 0.6, even though the deuteron cluster will form, the cluster will soon disappear owing to interactions with other nucleons. When the ratio is below 0.3, the density is too small for a deuteron cluster to form. Fig. 3(d) displays the distribution of the total energy of the deuteron cluster. The reason for negative total energies is a result of the choice to ignore the gamma photons.

Fig. 3.

The distribution of (a) formation time, (b) the relative kinetic energy, (c) density at the formation position, and (d) total energy, of a deuteron cluster for n+²⁸Si at 60 MeV. Here, p represents the percentage of the deuteron clusters in the distribution

Fig. 4 shows the deuteron cluster formation cross-section in ²⁸Si(n,d), ¹⁶O(n,d), and C12(n,d) reactions at beam energies of 20-100 MeV. The dotted line is the result of the deuteron cluster produced by fragmentation without the introduction of n+p→d reaction channel. The solid line is the result of the deuteron cluster produced by fragmentation and the neutron-proton collision, after the addition of this reaction channel. The hollow circles represent experimental data [66-68]. With the introduction of the n+p→d reaction channel into the modeling, the cross-section results are significantly optimized for each set of reactions, bringing the results much closer to the experimental data values. In Fig. 4(a), the prediction of the cross-section for ²⁸Si(n,d) reactions has the same trend as the experimental data at 20–50 MeV. By calculating the ratio of neutron-proton collisions forming a deuteron cluster, it is found that the contribution rate is greater than 52.6%. By the same method, the contribution rate of neutron-proton collisions for deuteron cluster formation is 36.5% for O16(n,d) reactions (Fig. 4(b)), and is 47.3% for C12(n,d) reactions (Fig. 4(c)). Hence, the addition of the n+p→d reaction channel is very significant when modeling the deuteron cluster cross-section. It is necessary to place greater focus on neutron-proton collisions when studying the effect on the formation and evolution of deuteron clusters.

Fig. 4.

(Color online) Measuring the cross-section of (a) ²⁸Si(n,d), (b) ¹⁶O(n,d), and (c) ¹²C(n,d) reactions at a beam energy range of 20-100 MeV. The lines are calculated by the IQMD+GEMINI model with and without the n+p→d reaction channel. The hollow circles are data measured from experiments [66-68].

4 Summary

A combination of the IQMD model and the GEMINI code was applied to study neutron-induced reactions. It was shown that the deuteron production cross-sections calculated by the standard version of the IQMD+GEMINI model are much smaller than the data. In an attempt to optimize the modeling of the deuteron production cross-sections, the n+p→d reaction channel was added to the IQMD model. The empirical formula of the cross-section and its in-medium factor of the n+p→d reaction were used. When the deuteron evolved to a medium density greater than 0.2ρ₀, it decayed into a neutron and a proton by the reaction d+N→n+p+N, where N is a nucleon. The contribution of neutron-proton collisions forming deuteron clusters was calculated. It was found that the contribution rate was greater than 52.6% for ²⁸Si(n,d) reactions, 36.5% for O16(n,d) reactions, and 47.3% for C12(n,d) reactions. The deuteron production cross-sections in the neutron-induced reactions ¹²C(n,d), ¹⁶O(n,d), and ²⁸Si(n,d) with incident energy 20–100 MeV were calculated and compared to the experimental data. The deuteron production cross-sections increased having incorporated the n+p→d reaction channel in the model, bring the values closer to the experimental data. The results indicated that the deuteron cluster is generated in a low-density region via the neutron-proton collision. In future work, the deuteron production mechanism in the QMD still requires further study, but our results hint that the n+p→d reaction channel should be considered in the transport model when studying deuteron clusters in HICs.

References

[1]

C.D. Bowman, G.J. Russell, J.W. Davidson et al.,

Nuclear energy generation and waste transmutation using an accelerator-driven intense thermal neutron source

. Nucl. Instrum. Meth. A 320, 336 (1992). https://doi.org/10.1016/0168-9002(92)90795-6