Estimation of uncertainty in lead spallation particle multiplicity and its propagation to a neutron energy spectrum

ABSTRACT

This paper presents an approach to uncertainty estimation of spallation particle multiplicity of lead (^natPb), primarily focusing on proton-induced spallation neutron multiplicity (x_pn) and its propagation to a neutron energy spectrum. The x_pn uncertainty is estimated from experimental proton-induced neutron-production double-differential cross-sections (DDXs) and model calculations with the Particle and Heavy Ion Transport code System (PHITS). Uncertainties in multiplicities for (n, xn), (p, xp), and (n, xp) reactions are then inferred from the estimated x_pn uncertainty and the PHITS calculation. Using these uncertainties, uncertainty in a neutron energy spectrum produced from a thick ^natPb target bombarded with 500 MeV proton beams, measured in a previous experiment, is quantified by a random sampling technique, and propagation to the neutron energy spectrum is examined. Relatively large uncertainty intervals (UIs) are observed outside the lower limit of the measurement range, which is prominent in the backward directions. Our findings suggest that a reliable assessment of spallation neutron energy spectra requires systematic DDX experiments for detector angles and incident energies below 100 MeV as well as neutron energy spectrum measurements at lower energies below ∼1.4 MeV with an accuracy below the quantified UIs.

KEYWORDS:

• Uncertainty
• spallation particle multiplicity
• spallation neutron
• double-differential cross-section
• random sampling
• PHITS
• lead

1. Introduction

Monte Carlo particle transport codes are widely used to simulate the transport of incident and secondary particles for neutronic and shielding analyses in various applications (e.g. fission and fusion reactors, accelerators, medical radiation, and space exploration). Among secondary particles produced from nuclear reactions, energy spectrum intensity of spallation neutrons has a particularly significant impact on the neutronic performance, radiation shielding, and radioactive waste management of the accelerator-driven nuclear transmutation systems (ADSs) [1–3] and the spallation neutron sources [4–9]. Reliable assessment of these neutronic characteristics requires detailed information about their uncertainties.

As a method for quantifying the uncertainties of neutronic characteristics, the so-called random sampling technique [10,11] has been developed. In order to apply this technique to the Monte Carlo particle transport calculation, uncertainties in nuclear data or nuclear reaction model parameters must be evaluated beforehand; because the uncertainties are expressed in the form of variance–covariance matrices, they are often referred to as covariance data [12]. To meet this requirement, a great deal of effort has been devoted to assessing the covariance data of the nuclear data libraries (e.g. ENDF [13], JENDL [14], and JEFF [15]) for neutrons below 20 MeV, and various uncertainties in reactor physics parameters have been quantified for the neutronic design of innovative nuclear systems such as ADSs (e.g. effective neutron multiplication factor $k_{\mathrm{e}\mathrm{f}\mathrm{f}}$ [16–18], effective delayed neutron fraction ${\beta }_{\mathrm{e}\mathrm{f}\mathrm{f}}$ [19,20], burnup reactivity [21], decay heat [22], and ${ }^{210}$Po production [23]). In the current body of research, few attempts have been made to estimate the uncertainty in nuclear data for protons and neutrons above 20 MeV. Stankovskiy et al. [24] endeavored the first such attempt for lead (${ }^{204,206,207,208}$Pb) and bismuth (²⁰⁹Bi) isotopes; however, a satisfactory estimation methodology has yet to be established.

This study focuses on the uncertainty in the spallation particle multiplicity of lead (^natPb) and its propagation to a neutron energy spectrum produced and emitted from a thick ^natPb target, even though the various uncertainties (e.g. nonelastic and elastic scattering cross-sections, energy–angle distributions of produced neutrons, energy loss of charged particles in materials) each affect it. Furthermore, although various particles, such as neutrons, protons, deuterons, pions, and photons, are emitted from the spallation reactions, this study exclusively focuses on protons and neutrons, which are major contributors to the formation of the spallation neutron energy spectrum. While Stankovskiy et al. [24] estimated the uncertainties in the proton-induced spallation neutron multiplicity ($x_{pn}$) and neutron-induced spallation neutron multiplicity ($x_{nn}$) using various samples of nuclear reaction models [25–27] and high-energy nuclear data libraries [28–30], in our study, we estimate the $x_{pn}$ uncertainty on the basis of experimental proton-induced neutron-production $(p,xn)$ double-differential cross-sections (DDXs) and model calculations using the Particle and Heavy Ion Transport code System (PHITS) version 3.08 [31]. Because experimental $(n,xn)$, $(p,xp)$, and $(n,xp)$ DDX data are scarce, uncertainties in the multiplicities for these reactions (i.e. $x_{nn}$, $x_{pp}$, and $x_{np}$) were inferred from the estimated $x_{pn}$ uncertainty and the PHITS nuclear reaction calculations, for which the default nuclear reaction model, a combination of the Intra-Nuclear Cascade model of Liège version 4.6 (INCL4.6) [27] and Generalized Evaporation Model (GEM) [32], was employed. For neutrons below 20 MeV, JENDL version 4.0 (JENDL-4.0) [14] was used. Then, the uncertainty in the spallation neutron energy spectrum was quantified using the estimated $x_{pn}$, $x_{nn}$, $x_{pp}$, and $x_{np}$ uncertainties via the random sampling technique; hereafter, their collective quantity, or the spallation particle multiplicity, will be referred to as $x$.

In Section 2, we will describe our estimation method of the $x$ uncertainty. Section 3 describes our uncertainty quantification method for a neutron energy spectrum, elucidates our calculation conditions, and offers a detailed discussion of our results. A summary and suggestions for future work are given in Section 4.

2. Estimation of uncertainty in spallation particle multiplicity

Using $x_{pn}$, the proton-induced neutron-production cross-section (${\sigma }_{pn}$) and its DDX can be expressed as

(1) ${\sigma }_{pn}(E_p)={\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}(E_p)\cdot x_{pn}(E_p)$(1)

and

(2) $\mathrm{D}\mathrm{D}\mathrm{X}(E_p;E^{\prime },\theta )={\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}(E_p)\cdot x_{pn}(E_p)\cdot f(E_p;E^{\prime },\theta ),$(2)

respectively, where $E_p$ is the incident proton energy and ${\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}$ is the proton-induced nonelastic scattering cross-section. Meanwhile, $f$ is the energy–angle distribution of the outgoing neutrons for $E_p$ as a function of polar angle $\theta$ in the laboratory frame and outgoing energy ($E^{\prime }$), which obeys the following normalization condition:

(3) $\int \int f(E_p;E^{\prime },\theta )sin\theta \mathrm{d}\theta \mathrm{d}{E}^{\prime }=1.$(3)

For simplicity, in this study, we assume that the model calculation provides correct $f$ and ${\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}$ values; that is, the relative uncertainty in $x_{pn}$ is equivalent to that in the magnitude of ${\sigma }_{pn}$ and DDX. While the calculated ${\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}$ values should also involve uncertainty, our previous study confirmed that Pearlstein–Niita systematics [33] implemented in PHITS reproduces experimental ${\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}$ data for lead fairly well (see Iwamoto et al. [34]); in other words, this uncertainty is negligible.

In order to estimate the uncertainty in $x_{pn}$, experimental spallation neutron multiplicities and their uncertainties were assessed in advance using experimental DDX data taken from the EXchange FORmat (EXFOR) experimental nuclear reaction database [35] as well as experimental data recently measured by Satoh et al. [36]. As an example, Figure 1 compares the PHITS nuclear reaction calculation with the experimental DDX data for the $p$(800 MeV) + ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}}$Pb reaction (Amian et al. [37], Nakamoto et al. [38], Ishibashi et al. [39], Ledoux et al. [40], Satoh et al. [41], and Trebukhovsky et al. [42]). Although PHITS broadly reproduces the experimental data, discrepancies are observed not only between the calculated and experimental values, but also among the experiments that should have been measured under the same projectile, target, and detector angle conditions; this discrepancy is likely linked to unknown experimental uncertainties (e.g. detection efficiency, energy calibration, and geometric measurement conditions). Thus, we first calculated the relative differences between the experimental DDX values and the PHITS calculations for each experiment ($d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{\mathrm{D}\mathrm{D}\mathrm{X}}$), and counted them as samples, which are expressed as

(4) ${\left(d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{\mathrm{D}\mathrm{D}\mathrm{X}}\right)}_{ij}=\frac{\mathrm{D}\mathrm{D}{\mathrm{X}}_{\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}({E^{\prime }}_i,{\theta }_j)-\mathrm{D}\mathrm{D}{\mathrm{X}}_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}({E^{\prime }}_i,{\theta }_j)}{\mathrm{D}\mathrm{D}{\mathrm{X}}_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}({E^{\prime }}_i,{\theta }_j)},$(4)

Figure 1. Comparison of DDX for $p$(800 MeV) + ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}}$Pb reaction between the PHITS calculations and experimental data

where $E_i^{\prime }$ and ${\theta }_j$ stand for $i$th measurement energy and $j$th measurement angle of the outgoing neutrons, respectively; $\mathrm{D}\mathrm{D}{\mathrm{X}}_{\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}$ denotes the measured DDX value; and $\mathrm{D}\mathrm{D}{\mathrm{X}}_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}$ indicates the calculated DDX value. Based on the above assumption, a sample of relative $x_{pn}$ difference ($d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}$) can be written as

(5) ${\left(d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}\right)}_{ij}={\left(d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{\mathrm{D}\mathrm{D}\mathrm{X}}\right)}_{ij}.$(5)

Figure 2 shows histograms of the frequency distribution for the

Figure 2. Histograms of relative difference from the PHITS calculation of DDXs for 20, 34, 48, 62.9, 63, 78, 113, 256, 597, and 800 MeV, as well as 1.0, 1.2, 1.5, 1.6, and 3.0 GeV $p$ + ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}/208}$Pb reactions, where the values in each bin are normalized such that the integral value is unity. In each panel, the datapoint with the error bar indicates $q_{\mathrm{M}}$ and uncertainty width ${ }_{-v_x}^{+u_x}$ of the distribution

$d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}\left(=d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{\mathrm{D}\mathrm{D}\mathrm{X}}\right)$ samples in each experiment (Satoh et al. [36,41], Guertin et al. [43], Meier et al. [44,45], Amian et al. [37,46], Nakamoto et al. [38], Ishibashi et al. [39], Ledoux et al. [40], and Trebukhovsky et al. [42]). Because the sample sizes of $i$ and $j$ were small, the $d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{\mathrm{D}\mathrm{D}\mathrm{X}}$ values were in most cases asymmetrically and multimodally distributed. Therefore, as reference values of the statistical distributions, we used the median ($q_{\mathrm{M}}$) and interquartile range (IQR), which are more robust indices than the mean and standard deviation. Table 1 summarizes key information about the experimental data and the $q_{\mathrm{M}}$ values with an uncertainty width of ${ }_{-v}^{+u}$ for the above distributions, in which we define the superscript $u$ and subscript $v$ as

(6) $u=2\left(q_{\mathrm{U}}-q_{\mathrm{M}}\right)/1.349,$(6)

(7) $v=2\left(q_{\mathrm{M}}-q_{\mathrm{L}}\right)/1.349,$(7)

Table 1. Summary of DDX experiments for ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}/208}$Pb$(p,xn)$ reactions and evaluation results for each incident proton energy

$E_p$ (MeV)	Target	Reference no.	Detector angle (°)	Sample size	$d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{\mathrm{D}\mathrm{D}\mathrm{X}}$ dist.	$x_{pn,\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}(n/p)$
20	^natPb	[36]	0	13	$-{0.684}_{-0.279}^{+0.414}$	${0.263}_{-0.263}^{+0.344}$
34	^natPb	[36]	0	17	$-{0.761}_{-0.689}^{+0.142}$	${0.528}_{-1.520}^{+0.314}$
48	^natPb	[36]	0	35	$-{0.498}_{-0.471}^{+0.057}$	${1.40}_{-1.31}^{+0.16}$
62.9	^208Pb	[43]	24, 35, 55, 80, 120	272	${0.170}_{-0.506}^{+0.656}$	${3.92}_{-1.70}^{+2.20}$
63	^natPb	[36]	0	19	$-{0.320}_{-0.261}^{+0.079}$	${2.28}_{-0.88}^{+0.26}$
78	^natPb	[36]	0	29	$-{0.431}_{-0.191}^{+0.300}$	${2.14}_{-0.72}^{+1.13}$
113	^natPb	[44]	7.5, 30, 60, 150	292	${0.471}_{-0.276}^{+0.321}$	${6.72}_{-1.26}^{+1.47}$
256	^natPb	[45]	7.5, 30, 60, 150	339	${0.211}_{-0.231}^{+0.172}$	${8.25}_{-1.57}^{+1.17}$
597	^natPb	[46]	30, 60, 120, 150	395	${0.071}_{-0.152}^{+0.147}$	${11.1}_{-1.6}^{+1.5}$
800	^natPb	[37–42]	0, 10, 15, 25, 30, 40, 55, 60, 85,
			90, 120, 130, 145, 150, 160	1110	${0.027}_{-0.250}^{+0.266}$	${12.2}_{-3.0}^{+3.2}$
1000	^natPb	[42]	15, 30, 60, 90, 120, 150	231	${0.002}_{-0.196}^{+0.184}$	${13.1}_{-2.4}^{+2.6}$
1200	^natPb	[40]	0, 10, 25, 40, 55, 70, 85, 100,
			115, 130, 145, 160	563	$-{0.050}_{-0.332}^{+0.285}$	${13.3}_{-4.7}^{+4.0}$
1500	^natPb	[39]	30, 60, 90, 120, 150	124	${0.275}_{-0.347}^{+0.634}$	${19.2}_{-5.2}^{+9.5}$
1600	^natPb	[40,42]	0, 10, 15, 25, 30, 40, 55, 60,
			85, 90, 100, 115, 120, 130,
			145, 150, 160	806	$-{0.040}_{-0.468}^{+0.402}$	${14.7}_{-7.2}^{+6.2}$
3000	^natPb	[39]	15, 30, 60, 90, 120, 150	121	$-{0.264}_{-0.185}^{+0.420}$	${12.3}_{-3.1}^{+7.0}$

where

(8) $\frac{u+v}{2}=\mathrm{I}\mathrm{Q}\mathrm{R}/\mathrm{1.349.}$(8)

Here, $q_{\mathrm{U}}$ and $q_{\mathrm{L}}$ indicate the upper and lower quartiles, respectively. The IQR/1.349 is an index known as the normalized IQR (NIQR), which equals standard deviation (1$\sigma$) when the population is normally distributed. Figure 3 illustrates a comparison of the proton-induced

Figure 3. Proton-induced ^natPb spallation neutron multiplicity (top panel: linear scale, bottom panel: logarithmic scale). The dot-dashed lines are the calculations with PHITS. The solid line is the JENDL-4.0/HE evaluation; the dotted line is the TENDL-2017 evaluation; the datapoints with error bars indicate $x_{\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}$

${ }^{\mathrm{n}\mathrm{a}\mathrm{t}}$Pb spallation neutron multiplicity between the PHITS calculations and the estimated statistics (${q_{\mathrm{M}}^{x_{pn}}}_{-v^{x_{pn}}}^{+u^{x_{pn}}}$) of the experimental values ($x_{pn,\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}$), in which the evaluations of the high-energy nuclear data library of JENDL-4.0 (JENDL-4.0/HE) [28] and the ninth version of TALYS evaluated nuclear data library (TENDL-2017) [29] are also drawn for later discussion. In this figure, a sample of $x_{pn,\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}$ was calculated as

(9) ${\left(x_{pn,\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}\right)}_{ij}=x_{pn,\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}\left\{1+{\left(d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}\right)}_{ij}\right\}.$(9)

In the above expression, $x_{pn,\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}$ denotes the spallation neutron multiplicity calculated using PHITS. The estimated statistics (${q_{\mathrm{M}}^{x_{pn}}}_{-v^{x_{pn}}}^{+u^{x_{pn}}}$) are also presented in Table 1. Although the PHITS results are generally consistent with the estimated experimental values for incident proton energies above 100 MeV, significant overestimation is observed below 100 MeV. In addition, INCL4.6 corroborates the experiments of Guertin et al. [43] according to Boudard et al. [27], but it fails to reproduce the neutron energy distribution for incident protons below 100 MeV at 0^° (see Satoh et al. [36]). As suggested by Cugnon and Henrotte [47], these trends may be attributed to the limited application of the INC model, which treats the nuclear reactions as classical kinematics using the classical collision term. Although the threshold of validity of the INC model is uncertain, for the sake of convenience, the uncertainty in $x_{pn}$ was estimated separately for below and above 100 MeV.

Figure 4 shows a histogram of

Figure 4. Histogram of relative difference from the PHITS calculations for incident proton energies above 100 MeV. The dashed and solid curves indicate a fit to the histogram with Gaussian distribution and Student’s t-distribution, respectively (top panel: linear scale, bottom panel: logarithmic scale)

$d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}(=d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{\mathrm{D}\mathrm{D}\mathrm{X}})$ for incident protons above 100 MeV, together with two curves fitted by the maximum likelihood estimation method with Gaussian distribution and Student’s t-distribution:

(10) $\mathcal{N}(d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}|\mu ,{\sigma }^2)=\frac{1}{{(2\pi {\sigma }^2)}^{1/2}}exp\left(-\frac{{(d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}-\mu )}^2}{2{\sigma }^2}\right)$(10)

and

(11) $\mathrm{S}\mathrm{t}(d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}|\mu ,\lambda ,\nu )=\frac{\mathrm{\Gamma }(\nu /2+1/2)}{\mathrm{\Gamma }(\nu /2)}{\left(\frac{\lambda }{\pi \nu }\right)}^{1/2}{\left(1+\frac{\lambda {(d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}-\mu )}^2}{\nu }\right)}^{-\nu /2-1/2},$(11)

respectively, where $\mu$ is the mean, $\sigma$ is the standard deviation, $\nu$ is the degree of freedom, $\lambda$ is the precision of the t-distribution, and $\mathrm{\Gamma }$ indicates the gamma function. Student’s t-distribution is an alternative probability distribution to Gaussian distribution, which is often used for the robust regression analysis (e.g. see Bishop [48]). To find the maximum likelihood solution constrained by the condition that $\mu =q_{\mathrm{M}}$, this study employed the expectation-maximization algorithm [49]. It can be observed from this figure that the width of the Gaussian fit is excessively large (${\sigma }_{\mathrm{f}\mathrm{i}\mathrm{t}}=0.999$). This is likely caused by the presence of outliers originating from large differences at high-energy edges of the neutron energy spectra. In contrast, the robustness of the t-distribution to outliers provides a good fit to the histogram. Based on the statistical analysis and the fitting with Student’s t-distribution, $\mu =0.049$, $\lambda =5.72$, and $\nu =4.01$ were obtained.

Because PHITS significantly overestimates the $x_{pn,\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}$ below 100 MeV, we first considered a simple fitting function for the estimated $x_{pn,\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}$ values. Then, the most probable $x_{pn}$ values with a probability distribution were estimated from the function, which is expressed as

(12) $x_{pn,\mathrm{f}\mathrm{i}\mathrm{t}}(E_p)=a(E_p-c{)}^b(n/p),$(12)

where $E_p$ is the incident proton energy in MeV; $a$ and $b$ are fitting parameters; and $c$ is a free parameter in MeV, which is chosen so that the most probable $x_{pn}$ line below 100 MeV connects to that above 100 MeV as smoothly as possible. Owing to the much smaller sample size of Satoh et al.’s [36] experimental data compared to that of Guertin et al. [43], both the data were each arranged in 1,000 samples via the bootstrap resampling technique [50]. By fitting these data using the weighted least-squares method, we chose $c=15$ and optimal values of $a=0.0936$ and $b=0.741$ were obtained. Figure 5 depicts a histogram of relative difference from Equation (12)(12) $x_{pn,\mathrm{f}\mathrm{i}\mathrm{t}}(E_p)=a(E_p-c{)}^b(n/p),$(12) with the optimal values

Figure 5. Histogram of relative difference from Equation.(12)(12) $x_{pn,\mathrm{f}\mathrm{i}\mathrm{t}}(E_p)=a(E_p-c{)}^b(n/p),$(12) for incident energies below 100 MeV (top panel: linear scale, bottom panel: logarithmic scale)

(${\delta }_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}$) for incident protons below 100 MeV, where ${\delta }_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}$ was calculated as

(13) ${\left({\delta }_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}\right)}_{ij}=\frac{{\left(x_{pn,\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}\right)}_{ij}-x_{pn,\mathrm{f}\mathrm{i}\mathrm{t}}}{x_{pn,\mathrm{f}\mathrm{i}\mathrm{t}}}.$(13)

From this statistical analysis, $q_{\mathrm{M}}$ and NIQR of the ${\delta }_{\mathrm{r}\mathrm{e}\mathrm{l}}^{x_{pn}}$ distribution were measured as 0.485 and 1.109, respectively. Unfortunately, however, no proper parametric function was obtained owing to the presence of many outliers in the prominently asymmetric and multimodal distribution.

The above analytical results are summarized in Table 2. Meanwhile, Figure 6 displays the same information as Figure 3, but it incorporates the estimated uncertainty interval (hereafter, referred to as UI), in which 50% and 90% UIs are drawn with dark and light bands, respectively. Here the upper and lower boundaries of the 50% UI were calculated as 75th and 25th percentiles, respectively, whereas those of the 90% UI were calculated as 95th and 5th percentiles, respectively. It can be observed that the UIs below 100 MeV are much larger than those above 100 MeV. This large uncertainty is due to the lack of systematic experimental data. Among the available references, Satoh et al. [36] measured DDXs for five incident energies (i.e. 20, 34, 48, 63, and 78 MeV) all dedicated to 0^°, and Guertin et al. [43] measured DDXs for five detector angles (i.e. 24^°, 35^°, 55^°, 80^°, and 120^°) but focused on an incident energy of 62.9 MeV. We note that, to reduce this uncertainty, further systematic DDX experiments for detector angles and incident energies below 100 MeV are needed.

Figure 6. Same as Figure 3, but including the uncertainty estimated from $x_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}$, $x_{\mathrm{f}\mathrm{i}\mathrm{t}}$, and $x_{\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}$ (top panel: linear scale, bottom panel: logarithmic scale). The dotted lines and the dark and light bands are $q_{\mathrm{M}}$ and 50% and 90% UIs, respectively

Table 2. Summary of statistical analyses for samples of $d_{\mathrm{r}\mathrm{e}\mathrm{l}}^x$ ($E_p\ge 100$ MeV) and ${\delta }_{\mathrm{r}\mathrm{e}\mathrm{l}}^x$ ($E_p<100$ MeV)

Incident energy	Statistics	Probability density function
$E_p\ge 100$ MeV	$q_{\mathrm{L}}=-0.170$, $q_{\mathrm{M}}=0.049$, $q_{\mathrm{U}}=0.270$,	Student’s t ($\mu =0.049$, $\lambda =5.72$, $\nu =4.01$)
$E_p<100$ MeV	$q_{\mathrm{L}}=-0.060$, $q_{\mathrm{M}}=0.485$, $q_{\mathrm{U}}=1.436$,	N/A

In contrast to the $(p,xn)$ reaction, few experimental DDX data are available for the $(n,xn)$, $(p,xp)$ and $(n,xp)$ reactions. However, because the nucleon emission mechanism at higher incident energies is described well by Serber’s classical picture [51], it is expected that the reproducibility of the $(p,xn)$, $(n,xn)$, $(p,xp)$, and $(n,xp)$ spectra via the PHITS nuclear reaction model based on Serber’s picture should be similar. Therefore, we assume that uncertainties for the $(p,xn)$, $(n,xn)$, $(p,xp)$, and $(n,xp)$ reactions share the same probability density function. Figure 7 demonstrates the calculated

Figure 7. Neutron-induced spallation proton multiplicity ($x_{nn}$) (left panels), proton-induced spallation proton multiplicity ($x_{pp}$) (middle panels), and neutron-induced spallation proton multiplicity ($x_{np}$) (right panels) inferred from the estimated $x_{pn}$ uncertainty and the PHITS nuclear reaction calculations (top panels: linear scale, bottom panels: logarithmic scale). The dark and light bands are 50% and 90% UIs, respectively

$x_{nn}$, $x_{pp}$, and $x_{np}$ together with estimated uncertainties. No proper experimental $(n,xn)$ and $(n,xp)$ data are found in EXFOR, which is because of the technical difficulty in producing monoenergetic neutrons. On the other hand, two pieces of experimental data for $(p,xp)$ are available, and the evaluation results for these data are listed in Table 3. It can be observed in Figure 7 that the two pieces of experimental $(p,xp)$ data of Iwamoto et al. [52] and Beck et al. [53] are within the 90% and 50% UIs, respectively. Therefore, the above assumption appears to be reasonable, although an uncertainty remains for the estimated UIs. Moreover, Figures 6 and 7 indicate that neutrons are produced more easily from the spallation reaction than protons. Accordingly, the $x_{pn}$ and $x_{nn}$ uncertainties are larger than the $x_{pp}$ and $x_{np}$ uncertainties. This phenomenon is due to neutron evaporation in the deexcitation process, which implies that neutrons produced from the spallation reaction, especially evaporation neutrons, largely affect the intensity of neutron fluxes or neutron energy spectra.

Table 3. Summary of DDX experiments for ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}/208}$Pb$(p,xp)$ reactions and evaluation results for each incident energy

$E_p$ (MeV)	Target	Reference no.	Detector angle (°)	Sample size	$d_{\mathrm{r}\mathrm{e}\mathrm{l}}^{\mathrm{D}\mathrm{D}\mathrm{X}}$ dist.	$x_{\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{t}}(p/p)$
392	${ }^{208}$Pb	[52]	20, 25, 30, 40, 45, 50, 60, 75, 90, 105	287	$-{0.255}_{-0.284}^{+0.295}$	${1.76}_{-0.50}^{+0.52}$
558	${ }^{\mathrm{n}\mathrm{a}\mathrm{t}}$Pb	[53]	20, 30, 40, 50, 60	498	${0.111}_{-0.184}^{+0.170}$	${3.12}_{-0.57}^{+0.53}$

The estimated statistics of $x_{pn}$, $x_{nn}$, $x_{pp}$ and $x_{np}$ (i.e. $q_{\mathrm{M}}$ and 50% UI) for representative incident proton energies are listed in Table 4. Here, for incident energies below 100 MeV, we did not estimate in this study the uncertainties for $x_{nn}$, $x_{pp}$, and $x_{np}$ owing to the lack of systematic experimental data and theoretical rationale. The nuclear reaction mechanism and its uncertainty in this energy range could be accounted for by some quantum mechanical approach using other nuclear reaction models such as CCONE [54] and TALYS [29]. However, Figure 3 indicates that the CCONE and TALYS, which were utilized in evaluating the JENDL-4.0/HE and TENDL-2017, respectively, presumably have room for improvement.

Table 4. Estimates of $x_{pn}$, $x_{nn}$, $x_{pp}$, and $x_{np}$ for ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}}$Pb at representative incident energies

	$x_{pn}$ $(n/p)$		$x_{nn}$ $(n/n)$		$x_{pp}$ $(p/p)$		$x_{np}$ $(p/n)$
Incident energy (MeV)	$q_{\mathrm{M}}$	50% UI	$q_{\mathrm{M}}$	50% UI	$q_{\mathrm{M}}$	50% UI	$q_{\mathrm{M}}$	50% UI
35	1.3	(0.8, 2.1)	–	–	–	–	–	–
40	1.5	(1.0, 2.5)	–	–	–	–	–	–
50	1.9	(1.2, 3.2)	–	–	–	–	–	–
80	3.1	(1.9, 5.0)	–	–	–	–	–	–
150	5.7	(4.5, 6.7)	7.1	(5.6, 8.6)	2.1	(1.7, 2.5)	0.6	(0.5, 0.7)
200	6.5	(5.0, 7.7)	7.7	(6.1, 9.3)	2.0	(1.6, 2.4)	0.8	(0.6, 0.9)
400	9.0	(6.9, 10.6)	10.3	(8.1, 12.4)	2.5	(2.0, 3.0)	1.2	(1.0, 1.5)
600	11.1	(8.6, 13.2)	12.7	(10.1, 15.4)	3.0	(2.4, 3.7)	1.8	(1.4, 2.2)
800	12.8	(9.9, 15.1)	14.4	(11.4, 17.4)	3.5	(2.8, 4.3)	2.2	(1.8, 2.7)
1000	14.0	(10.9, 16.6)	15.7	(12.4, 19.0)	3.9	(3.1, 4.7)	2.6	(2.1, 3.2)
1200	15.1	(11.6, 17.8)	16.6	(13.1, 20.1)	4.2	(3.3, 5.1)	2.9	(2.3, 3.6)
1500	16.2	(12.5, 19.1)	17.6	(14.0, 21.4)	4.6	(3.6, 5.6)	3.3	(2.6, 4.0)
2000	17.2	(13.3, 20.4)	18.6	(14.7, 22.5)	5.0	(4.0, 6.1)	3.7	(3.0, 4.5)
3000	18.0	(13.9, 21.3)	19.2	(15.2, 23.3)	5.4	(4.3, 6.5)	4.1	(3.3, 5.0)

3. Uncertainty quantification of a neutron energy spectrum

3.1. Methodology

The variables $x_{pn}$, $x_{nn}$, $x_{pp}$, and $x_{np}$ would typically correlate with one another (each correlation coefficient [$\rho$] must be in the range between $-1$ and 1); however, because the nuclear reaction model treats the particle emission on an event-by-event basis, estimating the correlations is difficult in practice. For simplicity, therefore, in this study, we assume that each set of variables has a total positive correlation; namely, $\rho =1$.

Let $x_{pn}$ be written as

(14) $x_{pn}=x_{pn,\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}\left(1+x_{\mathrm{r}\mathrm{e}\mathrm{l}}\right),$(14)

where $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ indicates a random variable obeying a probability distribution. Then, from Equation (1)(1) ${\sigma }_{pn}(E_p)={\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}(E_p)\cdot x_{pn}(E_p)$(1) , ${\sigma }_{pn}$ can be expressed as

(15) ${\sigma }_{pn}={\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}\cdot x_{pn,\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}\left(1+x_{\mathrm{r}\mathrm{e}\mathrm{l}}\right).$(15)

However, for the aforementioned reason, it is practically impossible to directly treat $x_{pn}$ as a simple random variable. In order to circumvent this problem, we assume that the nuclear reaction is stochastically triggered on the basis of ${\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}(1+x_{\mathrm{r}\mathrm{e}\mathrm{l}})$ and conduct the usual nuclear reaction simulation, which is expressed as follows:

(16) ${\sigma }_{pn}={\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}\left(1+x_{\mathrm{r}\mathrm{e}\mathrm{l}}\right)\cdot x_{pn,\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}},$(16)

where ${\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}$ is calculated using the Pearlstein–Niita systematics [33]. Although this assumption does not satisfy the event-by-event nuclear reaction kinematics, Equation (16)(16) ${\sigma }_{pn}={\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}\left(1+x_{\mathrm{r}\mathrm{e}\mathrm{l}}\right)\cdot x_{pn,\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}},$(16) is formally the same as Equation (15)(15) ${\sigma }_{pn}={\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}\cdot x_{pn,\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}\left(1+x_{\mathrm{r}\mathrm{e}\mathrm{l}}\right).$(15) . Thus, we quantified the uncertainty in the neutron energy spectrum according to the following procedure:

1. Sample $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ according to a probability distribution.

2. Multiply ${\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}$ and ${\sigma }_{n,\mathrm{n}\mathrm{o}\mathrm{n}}$ described in a PHITS source file by $(1+x_{\mathrm{r}\mathrm{e}\mathrm{l}})$ and save the file.

3. Compile a set of PHITS source files and create a PHITS executable file.

4. Calculate the neutron energy spectrum using the executable file.

5. Repeat the above steps $N$ times, where $N$ is the sample size.

6. Calculate the uncertainty from the $N$ output files via statistical processing.

As for the probability distribution of the uncertainty $p(x_{\mathrm{r}\mathrm{e}\mathrm{l}})$ for incident energies above 100 MeV, we employed Student’s t-distribution assessed in Figure 3 as follows:

(17) $p(x_{\mathrm{r}\mathrm{e}\mathrm{l}})=\mathrm{S}\mathrm{t}\left(x_{\mathrm{r}\mathrm{e}\mathrm{l}}|\mu ,\lambda ,\nu \right),$(17)

where $\mu =0.049$, $\lambda =5.72$, and $\nu =4.01$. To avoid unphysical occurrences, that is, a negative ${\sigma }_{\mathrm{n}\mathrm{o}\mathrm{n}}$ value, samples with $x_{\mathrm{r}\mathrm{e}\mathrm{l}}<-1$ (3.3% of total samples) were omitted. While this procedure slightly shifts the median and mean values from 0.049 to 0.067 and 0.099, respectively, its impact on the distribution shape was negligible. Because no reasonable probability function for incident energies below 100 MeV could be found, we first focused on the propagation of the uncertainty in $x$ above 100 MeV in this section; then, the impact of the $x$ below 100 MeV on the neutron energy spectrum is discussed in the subsequent section.

As an example, responses of $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ random samples to ${\sigma }_{pn}$ and ${\sigma }_{pp}$ are presented in Figure 8. Here,

Figure 8. Responses of $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ random samples to ${\sigma }_{pn}$ and ${\sigma }_{pp}$ for the 800 MeV proton + ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}}$Pb reaction (top panel: histogram of the $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ random samples, right panel: histogram of the ${\sigma }_{pn}$ and ${\sigma }_{pp}$ output values, middle panel: relationship between $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ and ${\sigma }_{pn}$ [${\sigma }_{pp}$]). The datapoints with error bars are the calculation values with 1$\sigma$ statistical uncertainties; the solid line is the least-squares linear fit for a set of the points with error bars. The sample size is 300

${\sigma }_{pn}$ and ${\sigma }_{pp}$ are obtained by the PHITS simulations, assuming that a thin (1 mm thick) ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}}$Pb target is bombarded with 800 MeV protons. We can see that the PHITS calculation exhibits a strong linear relationship between ${\sigma }_{pn}$ and $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$, as well as between ${\sigma }_{pp}$ and $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$, as anticipated in Equation (16)(16) ${\sigma }_{pn}={\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}\left(1+x_{\mathrm{r}\mathrm{e}\mathrm{l}}\right)\cdot x_{pn,\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}},$(16) .

3.2. Calculation condition

This study is based on an experiment of ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}}$Pb thick-target neutron yields (TTNYs) conducted by Meigo et al. [55] at KEK. In this experiment, neutron energy spectra, which were produced from a thick ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}}$Pb target (H15 cm $\times$ W15 cm $\times$ D20 cm) bombarded with 500 MeV and 1.5 GeV proton beams, were measured using the time-of-flight method with organic liquid-scintillator detectors for six angles (i.e. 15${ }^{\circ }$, 30${ }^{\circ }$, 60${ }^{\circ }$, 90${ }^{\circ }$, 120${ }^{\circ }$, and 150${ }^{\circ }$). This study focuses on the 500 MeV proton-induced neutron energy spectrum. The range for the 500 MeV protons in ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}}$Pb is 19.9 cm, for which all incident protons stop in the material. A schematic of the experimental condition is illustrated in Figure 9. In the PHITS simulation, the experimental geometry and the proton beam profile are precisely modeled.

Figure 9. A schematic of the experimental condition [55]

Using the PHITS experimental model, the neutron energy spectrum ($\psi$) at each detector angle is calculated with $4.5\times {10}^5$ histories ($4.5\times {10}^4$ source histories per one batch times 10 batches). For details, see the PHITS user’s manual [56]. The energy range of the tally is largely divided into four groups, as listed in Table 5. Hereafter, we express $\psi$ with the tally group number $\mathrm{\ell }$ as ${\psi }_{\mathrm{\ell }}$.

Table 5. Energy group structure for tally

	Energy range (MeV)
Group number $\mathrm{\ell }$	lower boundary	upper boundary
1	$100$	$500$
2	$10$	$100$
3	$1$	$10$
4	$0.1$	$1$

3.3. Sample size of the random sampling

To determine a proper sample size for quantifying the $\psi$ uncertainty, we first examine its convergence for the sample size. Figure 10 displays the convergence of

Figure 10. Convergence of ${\psi }_{\mathrm{\ell }}$ (left panels) and their 1$\sigma$-equivalent uncertainties (right panels) as a function of the sample size ($N$) for the detector angle of 30${ }^{\circ }$. For the left panels, the dots with error bars indicate the calculated ${\psi }_{\mathrm{\ell }}$ samples, the solid line is median of the ${\psi }_{\mathrm{\ell }}$ sample distributions, and the dashed line is a reference calculated without random sampling. For the right panels, the dashed line indicates $\stackrel{‾}{{\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}}}$, which is the 1$\sigma$ statistical uncertainty of ${\psi }_{\mathrm{\ell }}$ averaged over $N$ samples; the dot-dashed line is ${\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{o}\mathrm{b}\mathrm{s}}$, which is the observed statistical uncertainty obtained from $N$ PHITS runs; the solid line is ${\left(\delta {\psi }_{\mathrm{\ell }}\right)}_x$, which is the uncertainty due to $x$

${\psi }_{\mathrm{\ell }}$ and its uncertainty for each energy group $\mathrm{\ell }$ of the 30${ }^{\circ }$ detector angle, as a function of sample size $N$. Here, we define the uncertainty due to $x$, ${\left(\delta {\psi }_{\mathrm{\ell }}\right)}_x$, whose interval approximately corresponds to 1$\sigma$, as follows:

(18) ${\left\{{\left(\delta {\psi }_{\mathrm{\ell }}\right)}_x\right\}}^2={\left\{{\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{o}\mathrm{b}\mathrm{s}}\right\}}^2-{\left\{\stackrel{‾}{{\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}}}\right\}}^2,$(18)

where $\stackrel{‾}{{\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}}}$ is 1$\sigma$ statistical uncertainty of ${\psi }_{\mathrm{\ell }}$ averaged over $N$ samples and ${\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{o}\mathrm{b}\mathrm{s}}$ is the observed statistical uncertainty obtained from $N$ PHITS runs. Here, ${\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{o}\mathrm{b}\mathrm{s}}$ is defined as the NIQR of the calculated ${\psi }_{\mathrm{\ell }}$ distribution:

(19) ${\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{o}\mathrm{b}\mathrm{s}}=\left(q_{\mathrm{U}}^{{\psi }_{\mathrm{\ell }}}-q_{\mathrm{L}}^{{\psi }_{\mathrm{\ell }}}\right)/1.349,$(19)

where $q_{\mathrm{U}}^{{\psi }_{\mathrm{\ell }}}$ and $q_{\mathrm{L}}^{{\psi }_{\mathrm{\ell }}}$ are the upper and lower quartiles, respectively. For $\mathrm{\ell }=2,3,$ and $4$, it is seen in this figure that the ${\left(\delta {\psi }_{\mathrm{\ell }}\right)}_x$ values sufficiently converge at $N=800$, and $\stackrel{‾}{{\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}}}$ are almost negligible; that is,

(20) ${\left(\delta {\psi }_{\mathrm{\ell }}\right)}_x\simeq {\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{o}\mathrm{b}\mathrm{s}}.$(20)

For $\mathrm{\ell }=1$, although $\stackrel{‾}{{\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}}}$ was so large that Equation (20)(20) ${\left(\delta {\psi }_{\mathrm{\ell }}\right)}_x\simeq {\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{o}\mathrm{b}\mathrm{s}}.$(20) was not satisfied, ${\left(\delta {\psi }_{\mathrm{\ell }}\right)}_x$ appears to sufficiently converge at $N=800$. Therefore, considering the calculation cost and time spent, the random PHITS runs were ended at $N=800$, and then the ${\psi }_{\mathrm{\ell }}$ UIs were quantified via the statistical processing.

3.4. Propagation to a neutron energy spectrum

The estimated $\psi$ with $50\mathrm{\%}$ and $90\mathrm{\%}$ UIs propagated from the estimated $x$ uncertainty for the six angles are pictured in Figures 11 and 12, along with the experimental values, which are drawn using logarithmic and linear scales, respectively. Here, the UIs were obtained by the statistical processing of Equation (18)(18) ${\left\{{\left(\delta {\psi }_{\mathrm{\ell }}\right)}_x\right\}}^2={\left\{{\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{o}\mathrm{b}\mathrm{s}}\right\}}^2-{\left\{\stackrel{‾}{{\left(\delta {\psi }_{\mathrm{\ell }}\right)}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}}}\right\}}^2,$(18) from the ${\left(\delta {\psi }_{\mathrm{\ell }}\right)}_x$ distributions for the $800$ samples. It can be observed in these figures that most measurement values are within the $50\mathrm{\%}$ and $90\mathrm{\%}$ UIs for the measurement scope above 1.4 MeV. However, relatively large UIs are seen around a broad peak of several hundred kiloelectron volts, which is out of the measurement scope. This result indicates that true $\psi$ values likely exist within these UIs. It is important to note that the UIs are derived from the estimated $x$ uncertainty at energies from 100 to 500 MeV; therefore, other uncertainty sources could broaden the UIs. We also note that neutron yields around the peak energy are important for the neutronic performance and radiation shielding of high-energy accelerator facilities, as suggested by Iwamoto et al. [34,57]

Figure 11. Estimated $\psi$ along with experimental data, which is drawn using the logarithmic scale. The dark and light bands indicate 90% and 50% UIs, respectively, propagated from the estimated $x$ uncertainty above 100 MeV

Figure 12. Same as Figure 11, but drawn with the linear scale

To understand how the $x$ uncertainty propagates to the ${\psi }_{\mathrm{\ell }}$ uncertainties, we examine the relationships between $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ random samples and ${\psi }_{\mathrm{\ell }}$ for all the energy groups and detector angles, the results of which are presented in Figure 13, where, for each panel, the solid line indicates the weighted least-squares linear fit to the samples. Furthermore, as an example of the ${\psi }_{\mathrm{\ell }}$ distributions, a response of $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ random samples to ${\psi }_4$ at 30${ }^{\circ }$ is shown in Figure 14. In contrast to Figure 8, in which the cross-section values increase in proportion to $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$, the $\psi$ values exhibit a nonlinear relationship with $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$, and thus, the resultant $\psi$ distributions are more asymmetric than those of $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$, which is most prominent at $\mathrm{\ell }=1$. As described in Section 3.5, this asymmetric distribution is likely linked to the balance of neutron production and consumption in the nuclear reactions.

Figure 13. Relationships between $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ random samples and ${\psi }_{\mathrm{\ell }}$ ($\mathrm{\ell }=1,..,4$) for the six detector angles in the sampling energy range from 100 to 500 MeV. For each panel, the datapoints with error bars indicate the calculation values with 1$\sigma$ statistical uncertainties; the solid line is the least-squares linear fit for a set of points with the error bars; the dashed line with the dark and light bands indicates $q_{\mathrm{M}}$ with observed 50% and 90% UIs of $\psi$, respectively. The sample size is 200

Figure 14. Response of $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ random samples to ${\psi }_4$ (top panel: histogram of the $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ random samples, right panel: histogram of the ${\psi }_4$ output values, middle panel: relationship between $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ and ${\psi }_4$). The datapoints with error bars are the calculation values with 1$\sigma$ statistical uncertainties; the solid line is the least-squares linear fit for a set of points with error bars. The sample size is 800

Moreover, it can be observed from Figure 13 that as the detector angle shifts forward, the median $\psi$ values at energies below 10 MeV decrease. This phenomenon is likely due to the neutron shielding of the target material. Concurrently, we can also see that the median $\psi$ values at energies above 10 MeV increase as the detector angle shifts forward, which is due to the forward directionality of neutron emission in the intra-nuclear cascade process [25].

3.5. Impact of spallation particle multiplicity on a neutron energy spectrum

In order to examine the impact of $x$ on ${\psi }_{\mathrm{\ell }}$ in detail, the sampling energy range for incident protons from 100 to 500 MeV is divided into four groups, and a group for 20–100 MeV is added. The energy group structure for the random sampling is summarized in Table 6. Then, responses of $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ random samples in each sampling group to ${\psi }_{\mathrm{\ell }}$ are calculated, where the random seed used for $m=1,2,3$, and $4$ is kept uniform. Hereafter, we express ${\psi }_{\mathrm{\ell }}$ with the sampling energy group $m$ as ${\psi }_{\mathrm{\ell }}^{(m)}$. For $m=5$, we use a continuous uniform distribution function as a noninformative probability distribution function:

(21) $p(x_{\mathrm{r}\mathrm{e}\mathrm{l}})=\mathcal{U}(x_{\mathrm{r}\mathrm{e}\mathrm{l}}|\alpha ,\beta )$(21)

(22) $=\left\{\begin{array}{cc}\frac{1}{\beta -\alpha }& (\alpha \le x_{\mathrm{r}\mathrm{e}\mathrm{l}}\le \beta )\\ 0& (\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e})\end{array}\right.$(22)

Table 6. Energy group structure for random sampling

	Energy range (MeV)
Group number $m$	lower boundary	upper boundary
1	$400$	$500$
2	$300$	$400$
3	$200$	$300$
4	$100$	$200$
5	$20$	$100$

where $\alpha =-1$ and $\beta =1$. Figure 15 demonstrates the energy breakdown of the response of $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ random samples to ${\psi }_{\mathrm{\ell }}^{(m)}$ for the detector angle of 30${ }^{\circ }$, where, for each panel, the solid line represents the weighted least-squares linear fit to the samples. Most slopes of the linear fits are positive, but negative slopes are observed in the panel of ${\psi }_1^{(4)}$ and ${\psi }_2^{(5)}$. These trends are likely attributed to the balance of neutron production and consumption originating from the assumption of Equation (16)(16) ${\sigma }_{pn}={\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}\left(1+x_{\mathrm{r}\mathrm{e}\mathrm{l}}\right)\cdot x_{pn,\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}},$(16) , which is interpreted as follows. As $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ increases, more neutrons are produced from the nuclear reactions; therefore, ${\psi }_{\mathrm{\ell }}^{(m)}$ is generally expected to increase, and the slopes are expected to be positive. However, if most of the produced neutron energies are lower than the lower energy boundary of $\mathrm{\ell }$, the slopes for ${\psi }_{\mathrm{\ell }}^{(m)}$ will possibly be negative. Let us take a case of $m=4$ as an example. As $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ for $m=4$ increases, more neutrons, which are mainly composed of evaporation neutrons with energies below 20 MeV, are produced from the nuclear reaction, whereas more neutrons from 100 to 200 MeV are consumed as projectiles. This production and consumption balance yields the positive slopes for ${\psi }_2^{(4)}$, ${\psi }_3^{(4)}$, and ${\psi }_4^{(4)}$ and a negative slope for ${\psi }_1^{(4)}$ and eventually involves a strong distortion of the $\psi$ distribution.

Figure 15. Energy breakdown of the response of $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ random samples to ${\psi }_{\mathrm{\ell }}^{(m)}$ for the 30${ }^{\circ }$ detector angle

Figure 16 depicts a relative sensitivity matrix of

Figure 16. Relative sensitivity matrix of ${\psi }_{\mathrm{\ell }}^{(m)}$ to $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ for each detector angle

${\psi }_{\mathrm{\ell }}^{(m)}$ to $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$ for each detector angle, where the relative sensitivity matrix $\mathit{S}(=S_{lm})$ is defined as

(23) $\mathit{S}\equiv {\left.\frac{\mathrm{\partial }\mathit{\psi }/\mathit{\psi }}{\mathrm{\partial }{x}_{\mathrm{r}\mathrm{e}\mathrm{l}}/x_{\mathrm{r}\mathrm{e}\mathrm{l}}}\right|}_{x_{\mathrm{r}\mathrm{e}\mathrm{l}}=\mu }$(23)

(24) $\simeq \frac{\mathit{a}\mu }{\mathit{a}\mu +\mathit{b}},$(24)

where $\mu =0.049$; $\mathit{\psi }(={\psi }_{\mathrm{\ell }}^{\mathit{(}m\mathit{)}}\mathit{)}$ represents the neutron flux matrix, and $\mathit{a}(=a_{lm})$ and $b(={\mathit{b}}_{lm})$ are the slope and $y$-intercept matrices of the linear fits, respectively. It can be observed that, as the detector angle shifts backward, the ${\psi }_{\mathrm{\ell }}^{(1)}$ values are most sensitive to $x_{\mathrm{r}\mathrm{e}\mathrm{l}}$, but this trend is mitigated as the detector angle shifts forward. From these phenomena, we can interpret that ${\psi }_{\mathrm{\ell }}$ at backward directions are formed primarily from neutrons produced from a small number of nuclear reactions, whereas the forward ${\psi }_{\mathrm{\ell }}$ are formed mainly from neutrons produced from many nuclear reactions when passing through the target materials. These suppositions imply that the $x$ uncertainty below 100 MeV should be considered in quantifying the total ${\psi }_{\mathrm{\ell }}$ uncertainty, particularly at the forward directions, which also suggests that the $x$ uncertainty below 100 MeV could impact the neutron energy spectra for large-scale spallation target systems involving a large number of nuclear reactions.

4. Conclusions

In this study, we have estimated the uncertainty in $x_{pn}$ from experimental DDXs and the PHITS calculations and inferred the uncertainties in $x_{nn}$, $x_{pp}$, and $x_{np}$ on the basis of the estimated $x_{pn}$ uncertainty. Using these estimated uncertainties at incident energies from 100 to 500 MeV, we quantified the uncertainty in $\psi$ produced from the thick ${ }^{\mathrm{n}\mathrm{a}\mathrm{t}}$Pb target bombarded with 500 MeV proton beams using the random sampling technique, and propagation to $\psi$ has been examined. Consequently, we obtained the following findings:

• The $x_{pn}$ and $x_{nn}$ uncertainties at incident energies below 100 MeV should be considered in quantifying the total $\psi$ uncertainty. However, no proper probability density function for these uncertainties was found in this energy range, which we attribute to the lack of systematic experimental data.

• For neutron energies below 100 MeV, the 1$\sigma$-equivalent $\psi$ UIs ranged from 15% to 25% in relative value, depending on the neutron energies and detector angles. For neutron energies above 100 MeV, the $\psi$ distribution was distorted, especially at forward directions, which is likely due to the neutron production and consumption balance in the nuclear reactions originating from the assumption of Equation (16)(16) ${\sigma }_{pn}={\sigma }_{p,\mathrm{n}\mathrm{o}\mathrm{n}}\left(1+x_{\mathrm{r}\mathrm{e}\mathrm{l}}\right)\cdot x_{pn,\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}},$(16) .

• Relatively large $\psi$ UIs were observed out of the lower limit of the measurement range, that is, at energies below 1.4 MeV, which is prominent in the backward directions. Although true $\psi$ values most likely exist within the estimated UIs, uncertainty remains in the estimated values.

Although this study has taken the TTNY experiment of Meigo et al. as an example of the propagation of the $x$ uncertainty to $\psi$, our approach and findings are expected to be applicable to other TTNY experiments using high-energy proton beams and applications such as ADSs and spallation neutron source facilities. Based on the above findings and expectations, we have drawn the following challenges to be addressed in our future work:

• Reliable assessment of $\psi$ requires systematic $(p,xn)$ DDX measurements at incident energies below 100 MeV; experiments for the $(n,xn)$ reactions are also important.

• Validation of $\psi$ assessment requires experimental $\psi$ values at energies below $\sim$1.4 MeV, with an accuracy below the quantified $\psi$ UIs; that is, development of neutron detector systems that enable the measurement of $\psi$ in this energy range is required.

• Our random sampling approach will be applied to other TTNY experiments, ADSs, and/or spallation neutron source facilities, and subsequent issues on the impact of the $x$ uncertainty on the neutronic and shielding design will be investigated.

Acknowledgments

We would like to thank Dr. F. Maekawa from the Japan Atomic Energy Agency for his review of the manuscript and his comments.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Japan Society for the Promotion of Science (JSPS) Grants-in-Aid for Young Scientists B (Grant no. 17K14916).