JENDL/DEU-2020: deuteron nuclear data library for design studies of accelerator-based neutron sources

ABSTRACT

Intensive fast neutron sources using deuteron accelerators have been proposed for various applications. To contribute to the design study of such neutron sources, a deuteron nuclear data library for ^6,7Li, ⁹Be, and ^12,13C up to 200 MeV, JENDL/DEU-2020 is developed. Evaluation of JENDL/DEU-2020 is performed by employing the code system DEURACS with particular attention to neutron production data. Toward the evaluation of JENDL/DEU-2020, some modifications are made to DEURACS. Validation of the library is performed through simulation with the Monte Carlo transport calculation codes. From the simulation, it is shown that calculation results based on JENDL/DEU-2020 reproduce measured neutron production data well in incident energies up to 200 MeV. The new library is expected to make a large contribution to diverse design studies of deuteron accelerator neutron sources.

KEYWORDS:

• Deuteron
• nuclear data
• double-differential cross section
• thick-target neutron yield
• accelerator-based neutron source

1. Introduction

Since deuteron is a weakly bound system consisting of a proton and a neutron, it easily breaks up and emits a neutron through interaction with a target nucleus. The neutron generated by this process has a broad energy peak around half the incident deuteron energy. Utilizing this property, highly efficient fast neutron sources using deuteron accelerators have been proposed. Typical examples in the fields of nuclear data and nuclear physics are Nuclear Data Production System (NDPS) at RAON [1] and Neutron For Science (NFS) at SPIRAL-2 [2], respectively. Moreover, application to the production of medical radioisotope [3] and transmutation of long-lived radioactive nuclear waste [4] has also been considered.

In the above-mentioned accelerator-based neutron sources, $(d,xn)$ reactions on Li, Be, or C are employed to generate neutron beams. As for deuteron beam energy, a value of less than 50 MeV is typically adopted but one up to 400 MeV is considered as a candidate in the transmutation application [4]. In the design of accelerator facilities using beams with energies more than several hundred MeV, an approach based on the intra-nuclear cascade (INC) model has often been adopted. Many attempts have so far been made to extend the INC model framework to lower incident energies. As a result, recent benchmark studies have concluded that an improved INC model works reasonably well for nuclear reactions induced by light particles up to ${ }^{12}$C when the incident energy is higher than 100 MeV/nucleon [5,6]. Considering this situation, deuteron nuclear data for stable isotopes of Li, Be, and C up to 200 MeV, especially reliable data for neutron production, are required for diverse design studies of neutron sources using deuteron accelerators.

The current status of deuteron nuclear data for Li, Be, and C isotopes is as follows: The evaluated data for ${ }^{6,7}$Li are given in ENDF/B-VII.1 [7] and ENDF/B-VIII.0 [8], but the upper limits of incident energies are 5 MeV for ${ }^6$Li and 20 MeV for ${ }^7$Li. Other than the ENDF/B series, TENDL-2017 [9], which is produced based on the calculation of the TALYS code [10], stores the deuteron nuclear data files for ${ }^{6,7}$Li, ${ }^9$Be, and ${ }^{12,13}$C up to 200 MeV.¹ However, the values of TENDL-2017 do not necessarily show good agreement with experimental data as shown later in this paper. Other than the above published data, the evaluated data for neutron and $\gamma$-ray production from ${ }^{6,7}\mathrm{L}\mathrm{i}+d$ reactions have been developed [11]. These data are bundled with McDeLicious [12], an extension of the MCNP Monte Carlo code [13]. McDeLicious has been adopted to the engineering design of the International Fusion Materials Irradiation Facility (IFMIF) [14], which is the intensive neutron sources using the interaction of a liquid lithium target and a 40-MeV deuteron beam [15]. However, the evaluated data are given only up to 50 MeV since they were developed focusing on the design of IFMIF.

Under these circumstances, the purpose of this work is to develop a deuteron nuclear data library up to 200 MeV for ${ }^{6,7}$Li, ${ }^9$Be, and ${ }^{12,13}$C with particular attention to neutron production data. The new library is named JENDL/DEU-2020 as one of the series of JENDL special-purpose files. Upon the evaluation of JENDL/DEU-2020, we employ DEURACS, which is the code system dedicated to deuteron-induced reactions we have developed. DEURACS was successfully applied to the analyses of the production of nucleons [16,17], composite particles up to $A=4$ [18,19], and residual nuclei [20]. From these results, it is expected that DEURACS describes the mechanism of the deuteron-induced reaction well and is suitable for completing deuteron nuclear data through interpolation and extrapolation of available experimental values.

The remainder of this paper is organized as follows: Section 2 describes the evaluation method of JENDL/DEU-2020, focusing especially on neutron production. Toward the evaluation of JENDL/DEU-2020, we make some modifications to DEURACS so as to address issues found in our previous study. These modifications are also explained here. Next, we develop two application libraries based on JENDL/DEU-2020 for use in the Monte Carlo transport calculation codes. Outlines of these libraries are presented in Section 3. Then, validations of JENDL/DEU-2020 are performed through comparison with experimental data in Section 4. Finally, a summary and conclusions are given in Section 5.

2. Evaluation method

2.1. Calculation method for neutron production

This section outlines the theoretical models and methods in DEURACS to calculate the double differential cross sections (DDXs) of $(d,xn)$ reactions that give fundamental information on neutron production. Calculation methods in DEURACS for other quantities are described in Refs. [18–20]. Note that the terminologies relating to breakup processes in the following description follow the latest ones defined in a recent review paper [21] and are slightly different from those used in our previous paper [17].

In DEURACS, the DDXs of $(d,xn)$ reactions are expressed by incoherent summation of the following components:

(1) $\frac{d^2{\sigma }_{(d,xn)}}{dEd\mathrm{\Omega }}=\frac{d^2{\sigma }_{\mathrm{E}\mathrm{B}}}{dEd\mathrm{\Omega }}+\frac{d^2{\sigma }_{\mathrm{N}\mathrm{E}\mathrm{B}}}{dEd\mathrm{\Omega }}+\frac{d^2{\sigma }_{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}}{dEd\mathrm{\Omega }},$(1)

where $d^2{\sigma }_{\mathrm{E}\mathrm{B}}/(dEd\mathrm{\Omega })$, $d^2{\sigma }_{\mathrm{N}\mathrm{E}\mathrm{B}}/(dEd\mathrm{\Omega })$, and $d^2{\sigma }_{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}/(dEd\mathrm{\Omega })$ correspond to the DDXs for elastic breakup, nonelastic breakup, and pre-equilibrium and compound nucleus processes, respectively.

First, as for the breakup processes, the elastic breakup component is directly calculated by the continuum-discretized coupled-channels (CDCC) method [22].

Next, the nonelastic breakup component is calculated by the Glauber model with $S$-matrices given by the optical model [23]. However, the Glauber model cannot properly treat the $(d,n)$ transfer reactions to the specific bound states in the residual nucleus, which is a part of the nonelastic breakup. To deal with this problem, we separately calculate the transfer reaction by a conventional zero-range-distorted wave Born approximation (DWBA) using the DWUCK4 code [24]. This approach results in overlapping of the components calculated by the Glauber model and DWBA in the emission energies where the DWBA components exist. Therefore, we cut off the component calculated by the Glauber model above the emission energy corresponding to the highest excitation energy considered in the DWBA calculations. The remaining Glauber model component is normalized so that the total cross section of the nonelastic breakup is conserved. This means that the sum of the cross sections of the normalized Glauber model and DWBA components is equal to that of the initial Glauber model component.

Finally, we mention the calculation of pre-equilibrium and compound nucleus processes. In deuteron-induced reactions, three types of composite nuclei can be formed by the absorption of either neutron or proton in the incident deuteron or the deuteron itself. In DEURACS, a calculation taking these effects into account is performed by combining several models. The DDXs for the pre-equilibrium and compound nucleus processes are calculated as follows:

(2) $\begin{array}{cc}\frac{d^2{\sigma }_{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}}{dEd\mathrm{\Omega }}=\int d{E}_n{R}_n(E_n)& \frac{d^2{\sigma }_{(n,xn)}^{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}}{dEd\mathrm{\Omega }}\\ +& \int d{E}_p{R}_p(E_p)\frac{d^2{\sigma }_{(p,xn)}^{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}}{dEd\mathrm{\Omega }}+R_d\frac{d^2{\sigma }_{(d,xn)}^{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}}{dEd\mathrm{\Omega }},\end{array}$(2)

where $E_n$ and $E_p$ are energies of the neutron and proton absorbed in the target, and $R_n$, $R_p$, and $R_d$ denote the formation fractions of three different highly excited compound nuclei, which are calculated with the Glauber model. $d^2{\sigma }_{(n,xn)}^{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}/(dEd\mathrm{\Omega })$, $d^2{\sigma }_{(p,xn)}^{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}/(dEd\mathrm{\Omega })$, and $d^2{\sigma }_{(d,xn)}^{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}/(dEd\mathrm{\Omega })$ are the pre-equilibrium and compound nucleus components for the DDXs of neutron production from absorption of neutron, proton, and deuteron, respectively. They are calculated by the exciton model and the Hauser–Feshbach statistical model. We employed the subroutines of the CCONE code [25,26] for these calculations. The formulation of the two-component exciton model by Kalbach [27,28] and the Hauser–Feshbach model with the width fluctuation correction [29] are implemented in CCONE. The angular distribution of the pre-equilibrium component is calculated with the systematics proposed by Kalbach [30].

2.2. Particle decay from discrete levels

In our previous study on the $(d,xn)$ reactions on ${ }^9$Be and ${ }^{12}$C [17], the DEURACS calculation systematically underestimated experimental data in a few MeV emission energy regions. This underestimation is attributed to the problem that particle decay from discrete levels in residual light nuclei [e.g. ${ }^9$Be($E_x$ = 2.43 MeV)] and unstable nuclei (e.g. ${ }^5$He) is absent in the framework of DEURACS. Toward the evaluation of JENDL/DEU-2020, we modify the subroutines of CCONE to take into account the contribution from the above-mentioned particle decay.

We make each unbound state decayed according to the decay schemes provided in RIPL-3 [31], which are based on the data of ENSDF [32]. The angular distribution of the emitted particle is assumed to be isotropic in a center-of-mass system (CMS) and then converted into a laboratory system (LS). To obtain the energy spectrum in LS, the velocity distribution of the compound system is necessary. In the present work, the distribution is assumed to be a Gaussian function centered at an average velocity, whose width is determined so that the CMS average energy of the recoiled nuclei is equal to that obtained by the Hauser-Feshbach formulation at each particle emission step. Since Gaussian functions are closed under convolution, this assumption would greatly simplify the calculation for the velocity distribution of the compound system created after the emission of multi-particles. The LS energy spectrum for an emitted particle is obtained by convolution of the Gaussian distribution and delta-function in a three-dimensional velocity space.

2.3. Determination of model parameters

Some input parameters are necessary for the models integrated in DEURACS. First, nucleon optical potentials (OPs) at half the incident deuteron energy are necessary for the CDCC method and the Glauber model. We chose the global nucleon OPs by Koning and Delaroche (KD) [33] for both neutrons and protons, even though their lower limit of target mass range is $A$ = 24. The validity of applying the KD OPs to 1$p$-shell nuclei in the two models was discussed in our previous study [17].

Next, the input parameters used in the DWBA calculation are summarized in Table 1. They are the same as those used in Ref [17]. Discrete states of the residual nuclei considered in the DWBA calculation are listed in Tables 2–4 for the $(d,n)$ reactions on ${ }^{6,7}$Li, ${ }^9$Be, and ${ }^{12,13}$C targets, respectively. In the tables, some quantities corresponding to the considered states, namely, the excitation energy $E_x$, the spin-parity $J^{\pi }$, and the angular momentum transfer $l$ are summarized. A conventional DWBA calculation for nucleon transfer reactions requires spectroscopic factors (SFs) for single-nucleon orbits. We extracted the SFs corresponding to the individual discrete states presented in the tables. The SFs were extracted by fitting calculated DWBA cross sections to experimental differential cross sections. In the case where experimental differential cross sections are scarce (e.g. ${ }^7\mathrm{L}\mathrm{i}(d,n{)}^8\mathrm{B}\mathrm{e}$ reactions), we integrated the sharp peak seen in experimental DDXs and obtained corresponding differential cross sections. In our previous study for the $(d,p)$ reactions on ${ }^{12}$C, ${ }^{27}$Al, ${ }^{40}$Ca, and ${ }^{58}$Ni up to 100 MeV [16], we reported that the SFs extracted from our DWBA analyses have a weak incident energy dependence and the trend is similar among the four target nuclei. The incident energy dependence of the SF value extracted from the analyses for the ${ }^{12}\mathrm{C}(d,p{)}^{13}{\mathrm{C}}_{\mathrm{g}.\mathrm{s}.}$ reaction is given as follows [16]:

(3) $\begin{array}{cc}{S}_{\mathrm{C}-12,0}(E_d)=-2.18& \times {10}^{-6}{E_d}^3+3.19\times {10}^{-4}{E_d}^2\\ \hspace{0.17em}& -1.56\times {10}^{-2}{E}_d+8.20\times {10}^{-1},\end{array}$(3)

Table 1. Input parameters used in the DWBA calculation

Neutron potential	Koning–Deraloche (KD) [33]
Deuteron potential	Adiabatic potential [34] from KD
Proton binding potential	Woods–Saxon form
	$r_0$ = 1.25 [fm], $a_0$ = 0.65 [fm]
Finite range	0.7457
correction factor [fm]
Zero range constant	1.5 $\times$ 10${ }^4$
$D_0$ [MeV${ }^2$fm${ }^3$]
Nonlocality parameters	Neutron:0.85, deuteron:0.54

Table 2. Discrete states considered in the DWBA calculation for the ${ }^{6,7}\mathrm{L}\mathrm{i}(d,n)$ reactions

	$i$	$E_x$ [MeV]	$J^{\pi }$	$l$	$F_{k,i}$
${ }^6\mathrm{L}\mathrm{i}(d,n{)}^7\mathrm{B}\mathrm{e}$	0	0(g.s.)	${\frac{3}{2}}^{-}$	1	1.35
	1	0.43	${\frac{1}{2}}^{-}$	1	1.84
${ }^7\mathrm{L}\mathrm{i}(d,n{)}^8\mathrm{B}\mathrm{e}$	0	0(g.s.)	$0^{+}$	1	2.40
	1	3.04	$2^{+}$	1	2.25
	2	16.62	$2^{+}$	1	1.16

Table 3. Discrete states considered in the DWBA calculation for the ${ }^9\mathrm{B}\mathrm{e}(d,n)$ reactions

	$i$	$E_x$ [MeV]	$J^{\pi }$	$l$	$F_{k,i}$
${ }^9\mathrm{B}\mathrm{e}(d,n{)}^{10}\mathrm{B}$	0	0(g.s.)	$3^{+}$	1	1.14
	1	0.72	$1^{+}$	1	2.12
	2	1.74	$0^{+}$	1	0.24
	3	2.15	$1^{+}$	1	0.22
	4	3.59	$2^{+}$	1	0.14
	5	4.77	$3^{+}$	3	0.05
	6	5.11	$2^{-}$	0	0.21
	7	5.16	$2^{+}$	1	0.63
	8	5.92	$2^{+}$	1	0.60
	9	6.13	$3^{-}$	2	0.55
	10	6.56	$4^{-}$	2	0.25

Table 4. Discrete states considered in the DWBA calculation for the ${ }^{12,13}\mathrm{C}(d,n)$ reactions

	$i$	$E_x$ [MeV]	$J^{\pi }$	$l$	$F_{k,i}$
${ }^{12}\mathrm{C}(d,n{)}^{13}\mathrm{N}$	0	0(g.s.)	${\frac{1}{2}}^{-}$	1	1.00
	1	2.36	${\frac{1}{2}}^{+}$	0	0.76
	2	3.50	${\frac{3}{2}}^{-}$	1	0.23
	3	3.55	${\frac{5}{2}}^{+}$	2	0.97
${ }^{13}\mathrm{C}(d,n{)}^{14}\mathrm{N}$	0	0(g.s.)	$1^{+}$	1	1.14
	1	2.31	$0^{+}$	1	1.47
	2	3.95	$1^{+}$	1	0.43

where $E_d$ is the incident deuteron energy. As with Ref. [17], we assume that the SFs for the $(d,n)$ reactions have the same energy dependence as $S_{\mathrm{C}-12,0}(E_d)$. Since $S_{\mathrm{C}-12,0}(E_d)$ is determined for the ${ }^{12}\mathrm{C}(d,p{)}^{13}{\mathrm{C}}_{\mathrm{g}.\mathrm{s}.}$ reaction, an energy-dependent SF, $S_{k,i}(E_d)$, for each target nucleus $k$ and residual state $i$ is obtained as follows:

(4) $\begin{array}{c}{S}_{k,i}(E_d)=F_{k,i}{S}_{\mathrm{C}-12,0}(E_d),\end{array}$(4)

where $F_{k,i}$ is the scaling factor depending on $k$ and $i$. We determined $F_{k,i}$ for each $(d,n)$ reaction in the same manner as in Ref. [17], and each $F_{k,i}$ is summarized in Table 2–4. We assumed no incident energy dependence in the range above 100 MeV since $S_{\mathrm{C}-12,0}(E_d)$ was derived from the analyses up to 100 MeV.

Finally, in the calculation of the pre-equilibrium and compound nucleus processes, we employed the OPs of KD [33] for both neutrons and protons, those of An and Cai [35] for deuterons, and those of Avrigeanu et al. [36] for $\alpha$-particles. For tritons and ${ }^3\mathrm{H}\mathrm{e}$, the simplified folding potential with the nucleon OPs of Kunieda et al. [37] has been proposed [38], and we adopted the folding potential. Pre-equilibrium model parameters were taken from the work of Koning and Duijvestijn [39] except for those related to the squared matrix element $M^2$. In Ref. [39], $M^2$ is given by the following equation:

(5) $M^2=\frac{1}{A^3}\left\{6.8+\frac{4.2\times {10}^5}{{(\frac{E}{n}+10.7)}^2}\right\},$(5)

where $A$ is the target mass number, $E$ is the excitation energy, and $n$ is the number of excitons. Referring to the experimental data for neutron spectra shown in Section 4, we optimized the parameters 6.8 and 10.7 in Equation (5)(5) $M^2=\frac{1}{A^3}\left\{6.8+\frac{4.2\times {10}^5}{{(\frac{E}{n}+10.7)}^2}\right\},$(5) to 3.4 and 21.4, respectively. For level densities of the highly excited nucleus, we adopted the composite formula of Gilbert and Cameron [40] together with the systematics for the Fermi gas model parameters by Mengoni and Nakajima [41].

3. Application libraries based on JENDL/DEU-2020

JENDL/DEU-2020 provides the evaluated deuteron nuclear data for ${ }^{6,7}$Li, ${ }^9$Be, and ${ }^{12,13}$C in the incident deuteron energies from 200 keV to 200 MeV. The evaluated data were compiled according to the ENDF-6 format [42]. In addition to the ENDF-6 format files, we developed two application libraries based on JENDL/DEU-2020 for use in transport calculation with the Monte Carlo codes such as MCNP [13] and PHITS [43]. One is an ACE (A Compact ENDF) format [44] file for MCNP and the other is a ‘Frag Data’ format [45] file for PHITS. The ACE files were generated by processing the ENDF-6 files of JENDL/DEU-2020 using the NJOY-99 code [46] with some modifications [47]. Unfortunately, the present version of the PHITS code cannot treat an ACE file for deuteron. On the other hand, the Frag-Data format defined uniquely in PHITS is available also for deuteron. We thus adopted the format and developed the Frag-Data files almost equivalent to the ACE files.

Figure 1 shows neutron yields from a thick ${ }^7$Li target calculated by MCNP-6.2 with the ACE file and by PHITS-3.20 with the Frag-Data file. Statistical errors are less than 20$\mathrm{\%}$ in all energy bins. As shown in the figure, the two calculations are in good agreement in the emission energies below 30 MeV. Although deviations of several percents are observed, we have confirmed that these are mainly attributed to the difference of the calculation methods for deuteron stopping power in the two codes.

Figure 1. Neutron yields from a thick ${ }^7$Li target calculated by MCNP-6.2 with the ACE file and by PHITS-3.20 with the Frag-Data file. The incident deuteron energy is 25 MeV, and the emission angle is 15${ }^{\circ }$. Statistical errors are less than 20$\mathrm{\%}$ in all energy bins

However, unphysical peaks are seen in the high-emission energy region of the PHITS calculation. This is due to the interpolation method called ‘unit-base interpolation’ provided in the ENDF-6 format [42] which is not available in the Frag-Data format. In the unit-base interpolation, the energy spectrum at each incident energy point is first normalized so that the maximum emission energy is unity, and then the energy spectrum at an arbitrary incident energy is obtained by interpolation. On the other hand, only the simple interpolation is available in the Frag-Data format. The energy-integrated values are conserved between both calculations presented in Figure 1, but it is desirable that PHITS is improved to handle an ACE file for deuteron. For this reason, in the rest of this paper, a combination of MCNP-6.2 and ACE files is used for the validation of JENDL/DEU-2020 by transport calculations. Users should note this point when they do the PHITS calculation with the Frag-Data files.

4. Results and discussions

4.1. Analysis of neutron production mechanism

To understand the relations among the reaction processes, we first perform a component-by-component analysis of neutron production. The results are presented in Figure 2. The DDXs of ${ }^7\mathrm{L}\mathrm{i}(d,xn)$ reactions calculated with DEURACS are decomposed into three components as expressed in Equation (1)(1) $\frac{d^2{\sigma }_{(d,xn)}}{dEd\mathrm{\Omega }}=\frac{d^2{\sigma }_{\mathrm{E}\mathrm{B}}}{dEd\mathrm{\Omega }}+\frac{d^2{\sigma }_{\mathrm{N}\mathrm{E}\mathrm{B}}}{dEd\mathrm{\Omega }}+\frac{d^2{\sigma }_{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}}{dEd\mathrm{\Omega }},$(1) . To make the analysis clearer, the target in the calculation is assumed to be 100% ${ }^7$Li. In the figure, each component and the sum of the three components are shown and compared with the experimental data [48]. Note that the experimental data are those for natural lithium (92.5% ${ }^7$Li and 7.5% ${ }^6$Li). As presented in the figure, the sums of the three components well reproduce both the shape and magnitude of the experimental data regardless of the emission angles.

Figure 2. Calculated and experimental DDXs for the Li$(d,xn)$ reactions at 40 MeV. The dotted, dash-dot-dotted, and dashed curves represent the component of elastic breakup (EB), nonelastic breakup (NEB), and pre-equilibrium and compound nucleus processes (PE+CN), respectively. The solid curves are sums of each component. The squares are the experimental data taken from Ref [48]. The number at the top of each plot denotes the emission angle. The target in the calculation is ${ }^7$Li but that in the experiment is natural lithium (92.5% ${ }^7$Li and 7.5% ${ }^6$Li)

The sharp peak observed around 50 MeV at 0${ }^{\circ }$ is attributed to the $(d,n)$ transfer reactions, which is part of the nonelastic breakup. We folded the calculated peaks of the transfer reactions with a Gaussian function corresponding to the experimental energy resolution. Slight differences are seen in the peak position between the experimental and calculated values. This is expected because the effect of deuteron energy loss in the target is included in the experimental data. In the experiment by Hagiwara et al. [48], the target with 0.86-mm thickness was used. According to the calculation with the SRIM-2013 code [49], the stopping power in a lithium target is estimated as 1.2 MeV/mm for 40-MeV deuteron. The experimental small peak seen around 40 MeV at 0${ }^{\circ }$ is a contribution from the ${ }^6\mathrm{L}\mathrm{i}(d,n{)}^7\mathrm{B}\mathrm{e}$ reactions, which is not considered in the calculation in Figure 2. This will be discussed later in Section 4.2.

As for the broad peaks seen around half the deuteron incident energy at forward angles, they are formed by the breakup processes, namely, elastic and nonelastic breakup. The nonelastic breakup component is dominant at 0${ }^{\circ }$ but it has a stronger angle dependence than the elastic breakup component, and consequently the former is smaller than the latter at 30${ }^{\circ }$. This result demonstrates that it is necessary to consider the two breakup components for the accurate prediction of the DDXs of the ($d,xn$) reaction at various angles.

Almost all of low-energy components below 10 MeV are due to the pre-equilibrium and compound nucleus processes. Especially at 90${ }^{\circ }$, the two breakup components become very small and almost all spectra are contributions from the pre-equilibrium and compound nucleus processes. This indicates that the calculation method taking into account the three compound systems shown in Equation (2)(2) $\begin{array}{cc}\frac{d^2{\sigma }_{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}}{dEd\mathrm{\Omega }}=\int d{E}_n{R}_n(E_n)& \frac{d^2{\sigma }_{(n,xn)}^{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}}{dEd\mathrm{\Omega }}\\ +& \int d{E}_p{R}_p(E_p)\frac{d^2{\sigma }_{(p,xn)}^{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}}{dEd\mathrm{\Omega }}+R_d\frac{d^2{\sigma }_{(d,xn)}^{\mathrm{P}\mathrm{E}+\mathrm{C}\mathrm{N}}}{dEd\mathrm{\Omega }},\end{array}$(2) and the parameter setting presented in Section 2.3 are working well.

In addition, we assess the effect of the particle decay from discrete levels described in Section 2.2. Figure 3 shows the components of the pre-equilibrium and compound nucleus processes given by the DEURACS calculations with and without particle decay from discrete levels. In the figure, the horizontal axes are presented up to 10 MeV and the components of the other reaction processes are omitted. As shown in the figure, the reproducibility for the low-energy part of the experimental data obviously improves by considering the particle emission from discrete levels. One may notice that the difference between the two calculations becomes smaller at 90${ }^{\circ }$. This is because in the present calculation, the particles are assumed to be isotropically emitted from the complex system moving forward.

Figure 3. Same as Figure 2, but the ranges of the vertical and horizontal axes are changed. The dashed and dotted curves represent the component of pre-equilibrium and compound nucleus processes with and without particle decay from discrete levels (PD-DL), respectively

4.2. Validation on the DDXs for the Li($d,xn$) reactions

Next, we perform a validation on the DDXs for the ($d,xn$) reactions on natural lithium in a wide incident-energy range. Although the number of experimental data on the DDXs for the $(d,xn)$ reactions is limited, those for natural lithium targets are available at several incident energies up to 200 MeV and we can perform a systematic validation utilizing these data. We used the experimental values stored in the database of Experimental Nuclear Reaction Data (EXFOR) [50] except for 200 MeV. The experimental data of 200 MeV were provided by the authors of Ref. [51].

Figure 4 illustrates the experimental and calculated DDXs for the $(d,xn)$ reactions on natural lithium at 40 MeV. Comparisons are made also at the angles not shown in Figure 2, and calculations are performed on a natural lithium target. As mentioned in Section 4.1, the effect of deuteron energy loss in the target is seen in the experimental data. To take this effect into account, we obtained the calculated DDXs from the MCNP calculation for a thin lithium target using the ACE file of JENDL/DEU-2020. The calculated results are presented by solid lines in the figure. The thickness of the target is set to 0.86 mm according to the experimental condition. For comparison, we also present the calculation results with the models implemented by the PHITS code. In the PHITS calculation, the approach combining the Intra-Nuclear Cascade of Liège (INCL) [52] and DWBA developed by Hashimoto et al. [53] is adopted and the MWO formula [54] is chosen to calculate deuteron total reaction cross sections. The target thickness is the same also in the PHITS calculation. We also plot the values of TENDL-2017. For comparison with the experimental data, we folded the calculated peaks in the high-energy region with a Gaussian function corresponding to the experimental energy resolution. This treatment was commonly performed in the three calculations shown in Figure 4.

Figure 4. Calculated and experimental DDXs for the $(d,xn)$ reactions on natural lithium at 40 MeV. The solid and dashed lines represent the calculated results with MCNP-6.2 using the ACE file of JENDL/DEU-2020 and with the models implemented in PHITS-3.20. The dash-dot-dotted lines are the values of TENDL-2017. The squares are the experimental data taken from Ref [48]. The number at the top of each plot denotes the emission angle

As shown in the figure, the calculation results based on JENDL/DEU-2020 reproduce experimental data better than those based on the models in PHITS and the values of TENDL-2017 in a wide range of emission energies and angles. As discussed in the previous section, the experimental smaller peaks in the second highest energy region are contributions from the ${ }^6\mathrm{L}\mathrm{i}(d,n{)}^7\mathrm{B}\mathrm{e}$ reactions. They are reproduced by performing the calculation on a natural lithium target. In the calculation with the models in PHITS, the broad peaks can be seen around half the deuteron incident energy at forward angles even though the magnitudes and positions of peaks are different from the experimental ones. This suggests the application limit of the INCL model in PHITS for a low incident energy. TENDL-2017 underestimates the experimental values considerably at forward angles. This indicates that the empirical model by Kalbach [55] describing the breakup reactions implemented in TALYS does not work well under the present condition. In fact, the values of TENDL-2017 are in reasonable agreement with the experimental data above 10 MeV at angles larger than 60${ }^{\circ }$ where the breakup reactions are negligibly small. On the other hand, TENDL-2017 underestimates the experiment in the energy below 5 MeV. This is due to the fact that the particle emission from discrete levels discussed in Section 2.2 is not considered in the TALYS calculation.

Next, we make a validation in the lower incident energies. Figures 5 and 6 show the comparisons at the incident energies of 16.4 and 25 MeV, respectively. As mentioned in Section 1, the evaluated data for ${ }^7$Li of ENDF/B-VII.1 and -VIII.0 are available at incident energies up to 20 MeV and thus they are also plotted in Figure 5 assuming that the target is 100$\mathrm{\%}$ ${ }^7$Li. As in the case of 40 MeV, the results of the MCNP calculation with the ACE file of JENDL/DEU-2020 show the best agreement with the experimental data at 16.4 and 25 MeV. Note that the calculated values were obtained from the calculation for the same thick targets as the experimental conditions. Regarding the ENDF/B series, the evaluated values change considerably between the two versions but both underestimate the experimental data. The underestimation is due to that the evaluations of them were performed with the $R$-matrix theory and only the compound nucleus process was taken into account. Although only the ${ }^7\mathrm{L}\mathrm{i}(d,n\alpha {)}^4\mathrm{H}\mathrm{e}$ reaction was considered as a neutron emission channel in the evaluation of ENDF/B-VII.1, the ${ }^7\mathrm{L}\mathrm{i}(d,2n{)}^7\mathrm{B}\mathrm{e}$ reaction was taken into account in addition to the ${ }^7\mathrm{L}\mathrm{i}(d,n\alpha {)}^4\mathrm{H}\mathrm{e}$ reaction in ENDF/B-VIII.0. Therefore, the DDXs at low energies of ENDF/B-VIII.0 are larger than those of ENDF/B-VII.1. As for the models in PHITS and the values of TENDL-2017, the tendencies are almost the same as 40 MeV.

Figure 5. Same as Figure 4 but for 16.4 MeV. The dash-dotted and dotted lines are the evaluated values on ${ }^7$Li of ENDF/B-VIII.0 and ENDF/B-VII.1. The experimental data were taken from Ref [56]

Figure 6. Same as Figure 4 but for 25 MeV. The experimental data were taken from Ref [57]

The comparisons at higher incident energies are presented in Figures 7 and 8. As with the results in lower incident energies, the calculated values based on JENDL/DEU-2020 are in overall good agreement with the experimental data up to 200 MeV. However, in the case of 200 MeV, underestimation in the neutron energies between 20 and 60 MeV is seen and the reproducibility of the experimental data worsens as the emission angle increases. For further validation, it is desirable that experiment for emission angles larger than 30${ }^{\circ }$ is performed at incident energies above 100 MeV. Regarding the results with the models in PHITS, the agreement with the experimental data becomes better as the incident energy increases. This is because the picture of the INC model becomes more appropriate at the high incident energies. Consequently, the prediction accuracy of PHITS is almost the same as that of JENDL/DEU-2020 at 200 MeV. As mentioned in Section 1, 200 MeV (=100 MeV/nucleon) may be expected to be valid as an incident energy point for switching from nuclear data to an approach based on the INC models.

Figure 7. Same as Figure 4 but for 102 MeV. The experimental data were taken from Ref [58]

Figure 8. Same as Figure 4 but for 200 MeV. The experimental data were taken from Ref [51]

4.3. Validation on neutron yield from thick Li, Be, and C targets

For incident energies below 50 MeV, experimental data of neutron yields from deuteron bombardment on thick Li, Be, and C targets exist more than those of the DDXs for the $(d,xn)$ reactions. In this section, we present some examples of the validation results of JENDL/DEU-2020 using thick-target neutron yields (TTNYs). In the succeeding discussion, the results of the two calculations will be shown. One is the MCNP calculation with the ACE files based on JENDL/DEU-2020, and the other is the PHITS calculation with the models implemented in it. On the other hand, we omit the comparison in terms of the TENDL-2017 and the ENDF/B series, considering the validation results on the DDXs discussed in the last section. Note that there is no significant difference among the DDX data of TENDL-2017 regarding the ($d,xn$) reactions on Li, Be, and C isotopes, especially for the components of the breakup reactions. All experimental values were obtained from the EXFOR database. Target density, deuteron incident energy $E_d$, and target thickness used in the TTNY calculation are summarized in Table 5. The target thickness is set to be the same as that used in the corresponding TTNY experiment presented in the table. In all calculations, we assumed that the target is a cylinder with a radius of 10 mm and the deuteron beam is a cylindrical one with a radius of 0.1 mm.

Table 5. Input parameters used in the TTNY calculation. The target thickness is set to be the same as that used in the TTNY experiment presented in the rightmost column

Target	Density [g/cm${ }^3$]	$E_d$ [MeV]	Thickness [mm]	Experiment
Lithium	0.534	9	1.3	[59]
		25	7.5	[60]
		40	21.4	[48]
Beryllium	1.848	25	3.0	[60]
		50	10.0	[61]
Carbon	2.253	18	1.0	[59]
		40	6.0	[62]

Figures 9–11 illustrate the experimental and calculated TTNYs on natural lithium at 9, 25, and 40 MeV, respectively. As shown in the figure, the calculation results based on JENDL/DEU-2020 reproduce the experimental data better than the PHITS calculation results at each incident energy. Some features seen in neutron spectra reflect those of the DDXs for the ($d,xn$) reactions. The broad bump structure around half the deuteron incident energy is attributed to the components of the breakup reactions. As demonstrated in Figure 2, the consideration of both the elastic and nonelastic breakup components is necessary to predict this bump structure accurately at various angles. Since the yields of these bumps are large, it is expected that modeling of the breakup reactions is essential to predict total neutron production yields. On the other hand, the steplike structure in the high-emission energy region is formed mainly by the $(d,n)$ transfer reactions. In the case of a natural lithium target, this part covers a wider energy region than that of beryllium or natural carbon target shown later. This is because the $Q$-value of the ${ }^7\mathrm{L}\mathrm{i}(d,n{)}^8\mathrm{B}\mathrm{e}$ reaction is 15.03 MeV and higher than that of the ${ }^9\mathrm{B}\mathrm{e}(d,n{)}^{10}\mathrm{B}$ reaction (4.36 MeV) and the ${ }^{12}\mathrm{C}(d,n{)}^{13}\mathrm{N}$ reaction ($-0.28$ MeV). In the neutron irradiation applications, it is also important to evaluate the radioactivity of residual nuclei produced from reaction channels opened by only high-energy neutrons. Therefore, accurate prediction of neutron yields in the steplike structure is also important even though they are much smaller than those in the bump structure.

Figure 9. Experimental and calculated neutron yields from a thick natural lithium target bombarded by a 9-MeV deuteron. The solid and dashed lines represent the calculated results with MCNP-6.2 using the ACE file of JENDL/DEU-2020 and with the models implemented in PHITS-3.20, respectively. The squares are the experimental data taken from Ref [59]. The number at the top of each plot denotes the emission angle

Figure 10. Same as Figure 9 but for a 25-MeV deuteron on natural lithium. The experimental data were taken from Ref [60]

Figure 11. Same as Figure 9 but for a 40-MeV deuteron on natural lithium. The experimental data were taken from Ref [48]

Next, the results for the beryllium target at 25 and 50 MeV are presented in Figures 12 and 13, respectively. As in the case of a natural lithium target, the calculated values based on JENDL/DEU-2020 are in good agreement with the experimental data in wide emission energies and angles. One may note that the systematic underestimation seen in our previous calculation [17] in a few MeV emission energy regions is considerably improved. This is because the particle decay from discrete levels is taken into account in the evaluation of JENDL/DEU-2020. As for the reason of underestimation in the high-emission energy region at larger angles seen in the data of 50 MeV, the angular distribution of the transfer reaction may be unfavorable in both calculations. Aside from this, there is a possibility that the uncertainties associated with digitalizing experimental values of literature may be large in the EXFOR data since the neutron yields are shown in the linear scale in a figure in Ref [61].

Figure 12. Same as Figure 9 but for a thick beryllium target bombarded by a 25-MeV deuteron. The experimental data were taken from Ref [60]

Figure 13. Same as Figure 12 but for a 50-MeV deuteron on beryllium. The experimental data were taken from Ref [61]

Finally, the results for a natural carbon target (98.9% ${ }^{12}$C and 1.1% ${ }^{13}$C) at 18 and 40 MeV are presented in Figures 14 and 15, respectively. It should be noted that the experimental data at 40 MeV by Hagiwara et al. [62] are uniformly multiplied by a factor of 2.0 according to the information provided in EXFOR. The calculations with JENDL/DEU-2020 reproduce the experimental data better than the PHITS calculations.

Figure 14. Same as Figure 9 but for a thick natural carbon target bombarded by a 18-MeV deuteron. The experimental data were taken from Ref [59]

Figure 15. Same as Figure 14 but for a 40-MeV deuteron on natural carbon. The experimental data were taken from Ref [62]

The validation results demonstrate that transport simulation based on JENDL/DEU-2020 well reproduces the experimental data on neutron production for various combinations of target nuclei and incident energies. From these results, it is expected that neutron production data of JENDL/DEU-2020 are reliable and the library makes a large contribution to design studies of neutron sources with deuteron accelerator.

5. Summary and conclusions

We developed JENDL/DEU-2020, a deuteron nuclear data library for ${ }^{6,7}$Li, ${ }^9$Be, and ${ }^{12,13}$C up to 200 MeV. The evaluation of JENDL/DEU-2020 was carried out by employing the code system DEURACS with particular attention to neutron production data. DEURACS was modified to take into account the contribution of particle decay from discrete levels. To be used in the Monte Carlo transport calculation codes, two application libraries based on JENDL/DEU-2020 were developed. One is the ACE files for MCNP and the other is the Frag-Data files for PHITS. Validation of JENDL/DEU-2020 was performed through transport simulation by MCNP-6.2 with the ACE files. From the simulation, it was demonstrated that the calculation results based on JENDL/DEU-2020 reproduced the measured neutron production data well in the incident energies up to 200 MeV. From these results, JENDL/DEU-2020 is expected to make a large contribution to diverse design studies of deuteron accelerator neutron sources. The ENDF-6 files of JENDL/DEU-2020 and the ACE files and the Frag-Data files based on JENDL/DEU-2020 will be available from the website of the Nuclear Data Center of Japan Atomic Energy Agency.

Acknowledgments

The authors wish to thank C. Konno for providing helpful information on the development of an ACE format file and on the calculation with the MCNP code. They are also grateful to S. Hashimoto for the fruitful discussion on the development of a Frag-Data format file.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was partially funded by JSPS KAKENHI Grant Number 19K15483 from Japan Society for the Promotion of Science.

Notes

1. In TENDL-2019, the latest version of TENDL series, the deuteron nuclear data for ${ }^{6,7}$Li are taken from ENDF/B-VII.1 and the data for ${ }^9$Be and ${ }^{12,13}$C are almost identical with those of TENDL-2017.