Calculation of the Coulomb barrier and the level density parameter at the fission saddle point in spallation nuclear reactions

The knowledge of level densities is crucial for testing the validity of models of nuclear structure, and a lot of theoretical and experimental work has been made [1, 2], leading to several models to fit experimental data. These models take into account the shell and pairing effects, as well as the rotational–vibrational character of nuclei, which, in one way or another, has been implemented in several computational codes [3], especially those dedicated to the study of spallation reactions.

Based on the close relationship between the nuclear structure and the level density, the latter must determine, in some way, the behavior of most of the quantities, such as the Coulomb barrier (CB) for charged particles' emission and the level density parameter at the fission saddle point. In this respect, CB measures the amount of energy ΔΛ needed to put or remove one unit of electric charge from the nucleus, so that for the ground state the electrostatic formula corrected for barrier penetration should be enough to calculate ΔΛ, while for excited states (when the nuclear structure has changed), a dependence on the excitation energy must appear.

A change of the excitation energy means a variation of nuclear temperature, which is related to the level density through its thermodynamic definition. Then it is possible to think about the existence of some relationship reflecting the fact that a change in CB is due to a change in the nuclear structure. Since the level density parameter characterizes the behavior of the level density, it is almost intuitive to try to establish a relationship between the level density parameter at the saddle point and the level density model used in the calculations.

Nevertheless, both quantities commonly receive the treatment of adjustable parameters in the codes devoted to studying spallation nuclear reactions [4, 5] and their relationship with the level density remains unexplored. Thus, CB has the form V_C = V₀f(U), with V₀ = κV_E, where V_E comes from electrostatics and κ is a correction that takes into account the barrier penetration (tunneling effect). The factor f(U) is another correction introduced to reflect the expected reduction of the barrier due to the expansion and surface oscillations of an excited nucleus, which is treated in several ways [6, 7] that, in general, are not related with the level density. Something similar happens with the level density parameter at the fission saddle point a_f.

That is why the objective of this work is to study the relationship between the level density and both the Coulomb barrier for charged particles' emission and the level density parameter at the saddle point. In this regard, expressions corresponding with the partial widths for light particles' emission and occurrence of fission will be used together with the level density model of Egidy and Bucuresco [8], based on the back shifted Fermi gas (BSFG), to derive two differential equations that allow the calculation of these magnitudes.

The structure of the present work is as follows: sections 2 and 3 develop new methods for calculation of CB and the level density parameter at the fission saddle point, respectively. Section 4 describes the features of the statistical module of the CRISP code, as well as the changes made to incorporate the ideas of sections 2 and 3. Section 5 shows the results of the simulations using the new version of the code. Finally, section 6 highlights the results of the work, together with an assessment of them, including a comparison with other de-excitation models and establishing some directions for future work and development of the code.

Historically, spallation nuclear reactions have been treated by a mechanism which includes two stages [9]: a sequence of binary collisions initiated by the incident nucleon, known as the intranuclear cascade (INC), followed by the decay, through several channels, of the excited compound nucleus. The second step of the model takes into account the evaporation-fission competition, and the partial width associated with the emission of light particles is [10, 11]

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{l}={\hslash }{\int }_{{V}_{l}}^{U-{B}_{l}}\pi {r}_{0}^{2}{A}^{2/3}(\epsilon -{V}_{l})\displaystyle \frac{{g}_{l}{m}_{l}{\rho }_{F}}{{\pi }^{2}{{\hslash }}^{3}{\rho }_{P}}\,{\rm{d}}\epsilon ,\end{eqnarray} \tag{ 1 }$$

where B_l is the light particle's separation energy, ρ_P and ρ_F are the level densities of the nucleus before and after the emission, respectively, r₀ is the Bohr radius, A the mass number before the emission, m_l the mass of the emitted particle and g_l the occupation number, this expression strongly depends on both the level density ρ and the Coulomb barrier height V_l.

After some operations, it is possible to write (see the appendix)

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{l}=\displaystyle \frac{\gamma }{{a}_{l}}{\int }_{0}^{{X}_{l}}({X}_{l}-x){\rho }_{{Fx}}(x){\rm{d}}{x},\end{eqnarray} \tag{ 2 }$$

where $x={a}_{l}(U-{E}_{{Fl}}-{B}_{l}-\epsilon )$, ${X}_{l}={a}_{l}(U-{E}_{{Fl}}-{B}_{l}-{V}_{l})$ and

$$\begin{eqnarray}&&{\rho }_{{Fx}}(x)=\displaystyle \frac{\exp (2\sqrt{x}\ \ )}{{x}^{3/4}{[1+\sqrt{1+4x}]}^{1/2}},\end{eqnarray} \tag{ 3 }$$

being U the excitation energy, while a_l and E_Fl are the level density parameter and the energy correction after emission of the particle l (using the subscripts n, p and a for neutrons, protons and alpha particles, respectively).

The light charged particles' to neutrons' emission ratio is

$$\begin{eqnarray}&&{\rm{\Theta }}=\displaystyle \frac{{{\rm{\Gamma }}}_{l}}{{{\rm{\Gamma }}}_{n}}=\displaystyle \frac{{a}_{n}{g}_{l}{m}_{l}{A}_{n}^{5/6}}{{a}_{l}{g}_{n}{m}_{n}{A}_{l}^{5/6}}\displaystyle \frac{{I}_{l}}{{I}_{n}},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&{I}_{l}={\int }_{0}^{{X}_{l}}({X}_{l}-x){\rho }_{{Fx}}(x){\rm{d}}{x},\end{eqnarray} \tag{ 5 }$$

while I_n is the same as equation (5), changing X_l for ${X}_{n}={a}_{n}(U-{E}_{{Fn}}-{B}_{n})$. Now, considering Θ = Θ(U, V_l) leads to (see the appendix)

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Theta }}/\partial U}{\partial {\rm{\Theta }}/\partial {V}_{l}}=-1+\displaystyle \frac{{a}_{n}({X}_{l}-\overline{{X}_{l}})}{{a}_{l}({X}_{n}-\overline{{X}_{n}})},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&\overline{{X}_{l}}=\displaystyle \frac{{\int }_{0}^{{X}_{l}}x{\rho }_{{Fx}}(x){\rm{d}}{x}}{{\int }_{0}^{{X}_{l}}{\rho }_{{Fx}}(x){\rm{d}}{x}}.\end{eqnarray} \tag{ 7 }$$

For X_l values such that ${X}_{l}\gg 1$ arises the following expression,

$$\begin{eqnarray}&&\overline{{X}_{l}}\approx {X}_{l}\exp \left(-\displaystyle \frac{1}{\sqrt{{X}_{l}}}\right),\end{eqnarray} \tag{ 8 }$$

as long as the derivative of V_l respect to U is (see the appendix)

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{V}}_{l}}{{\rm{d}}{U}}=\displaystyle \frac{\partial {\rm{\Theta }}/\partial U}{\partial {\rm{\Theta }}/\partial {V}_{l}}+{\rm{\Xi }}(U,{V}_{l}).\end{eqnarray} \tag{ 9 }$$

Equation (9) has two unknowns, namely, dV_l/dU and Ξ(U, V_l), so more information is needed to solve the problem. This information comes from the assumptions made on the asymptotic behavior of V_l. In this way, assuming that V_l is a monotonically decreasing function which tends to a stationary value at high energies results in

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{V}}_{l}}{{\rm{d}}{U}}\lt 0,\quad \mathop{\mathrm{lim}}\limits_{U\to \infty }\displaystyle \frac{{{\rm{d}}{V}}_{l}}{{\rm{d}}{U}}=0,\end{eqnarray} \tag{ 10 }$$

generalizing the accepted idea that CB decreases with increasing excitation energy [6, 7]. If one accepts the conditions given in (10), then it is possible to estimate Ξ(U, V_l) as described below.

At high energies equation (9) takes the form

$$\begin{eqnarray}&&\mathop{{\rm{lim}}}\limits_{U\to \infty }\displaystyle \frac{{{\rm{d}}{V}}_{l}}{{\rm{d}}{U}}=0\to \mathop{{\rm{lim}}}\limits_{U\to \infty }{\rm{\Xi }}(U,{V}_{l})=1-{\left[\displaystyle \frac{{a}_{n}}{{a}_{l}}\right]}^{1/2}.\end{eqnarray} \tag{ 11 }$$

Combining equation (11) and the exponential introduced in equation (8) leads to the following proposal:

$$\begin{eqnarray}&&{\rm{\Xi }}(U,{V}_{l})=\left[1-{\left(\displaystyle \frac{{a}_{n}}{{a}_{l}}\right)}^{1/2}\right]\exp \left(-\displaystyle \frac{1}{\sqrt{{X}_{l}}}\right),\end{eqnarray} \tag{ 12 }$$

finally getting that

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{V}}_{l}}{{\rm{d}}{U}}=-1+\displaystyle \frac{{a}_{n}({X}_{l}-\overline{{X}_{l}})}{{a}_{l}({X}_{n}-\overline{{X}_{n}})}+\left[1-{\left(\displaystyle \frac{{a}_{n}}{{a}_{l}}\right)}^{\tfrac{1}{2}}\right]\exp \left(-\displaystyle \frac{1}{\sqrt{{X}_{l}}}\right),\end{eqnarray} \tag{ 13 }$$

which allows the calculation of the Coulomb barrier height as a function of the excitation energy, starting from the initial condition V_l(U₀) = V₀ = κV_E, with ${U}_{0}={V}_{0}+{B}^{* }$, being ${B}^{* }$ the highest value between B_l and B_n, whilst κ is 0.70 for protons and 0.83 for alpha particles.

Figure 1(a) shows CB for proton emission as a function of the excitation energy for several nuclei, obtained from equation (13), while figure 1(b) shows the same, but using a different approach given by equation (28) (see section 4). Comparing both figures it is concluded that for a given nucleus at a certain value of the excitation energy, equation (13) produces lower values of CB, what increases X_l in equation (2) and consequently Θ(U, V_l). So that, if implemented within a computational code, equation (13) should lead to higher values of the probability of emission of light charged particles (see section 5).

Zoom In Zoom Out Reset image size

Of course, the choice of Ξ(U, V_l) remains, to some extent, arbitrary, because the total derivative of Θ(U, V_l) is unknown. Nevertheless, regardless of the selection, Ξ(U, V_l) will always become important at high energies, to compensate for the behavior of the first term of the right side of equation (9), so that it is zero at low energies. Then it is possible to say that the curves Ξ_m(U, V_l) are, at least, close in the sense of zero-order proximity, and the absolute value of the difference between any two functions of the set is small [12].

A natural choice would be Ξ(U, V_l) = 0, but this option, together with the conditions of equation (10) would force the level density parameter to be a function of the excitation energy, in such a way that

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{U\to \infty }{a}_{n}=\mathop{\mathrm{lim}}\limits_{U\to \infty }{a}_{l},\end{eqnarray} \tag{ 14 }$$

which is not the case with this work, where it is independent of the excitation energy. The same need for an energy-dependent level density parameter has been previously pointed out as a necessary condition to account for the damping of the shell effects [13].

Now it is possible to understand why the function Ξ(U, V_l) remains with a certain degree of arbitrariness. It has to guarantee compliance with the conditions given in (10), regardless of the choice of the level density parameter.

The partial width associated with the occurrence of fission has the form [14]

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{f}=\displaystyle \frac{1}{2\pi {\rho }_{P}({U}_{P})}{\int }_{0}^{{U}_{P}-{B}_{f}}{\rho }_{P}^{* }({a}_{f},{U}_{P}-{B}_{f}-\epsilon )\,{\rm{d}}\epsilon ,\end{eqnarray} \tag{ 15 }$$

where ${\rho }_{P}^{* }$ is the nuclear level density at the saddle and B_f is the fission barrier height, while the ratio between the widths of fission and neutron emission is

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{f}=\displaystyle \frac{{{\rm{\Gamma }}}_{f}}{{{\rm{\Gamma }}}_{n}}=\displaystyle \frac{{a}_{n}{(A-1)}^{5/6}{{\hslash }}^{2}}{4{m}_{n}{r}_{0}^{2}{A}^{3/2}}\displaystyle \frac{{I}_{f}}{{I}_{n}},\end{eqnarray} \tag{ 16 }$$

with

$$\begin{eqnarray}&&{I}_{f}={\int }_{0}^{{X}_{f}}{\rho }_{{Fx}}(x){\rm{d}}{x},\end{eqnarray} \tag{ 17 }$$

being X_f = a_f(U_P − B_f). Considering Θ_f = Θ_f(U, a_f) and following a similar path to the one taken in section 2 it turns out that

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Theta }}}_{f}/\partial U}{\partial {{\rm{\Theta }}}_{f}/\partial {a}_{f}}=\displaystyle \frac{{a}_{f}}{{U}_{P}-{B}_{f}}\left[1-\displaystyle \frac{{I}_{f}}{{I}_{n}}\displaystyle \frac{\tfrac{\partial {I}_{n}}{\partial U}}{\tfrac{\partial {I}_{f}}{\partial U}}\right],\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {I}_{f}}{\partial U}={a}_{f}{\rho }_{{Fx}}({X}_{f}),\end{eqnarray} \tag{ 19 }$$

whilst

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{a}}_{f}}{{\rm{d}}{U}}=\displaystyle \frac{\partial {{\rm{\Theta }}}_{f}/\partial U}{\partial {{\rm{\Theta }}}_{f}/\partial {a}_{f}}+{{\rm{\Xi }}}_{f}(U,{a}_{f}).\end{eqnarray} \tag{ 20 }$$

This time, the condition imposed to a_f is that it must be a stationary function at high energies; that means

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{U\to \infty }\displaystyle \frac{{{\rm{d}}{a}}_{f}}{{\rm{d}}{U}}=0,\end{eqnarray} \tag{ 21 }$$

which is fulfilled by equation (18) if the additional condition

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{U\to \infty }{a}_{f}={a}_{n}\end{eqnarray} \tag{ 22 }$$

is also taken into account. Then it is possible to establish that

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{a}}_{f}}{{\rm{d}}{U}}=\displaystyle \frac{\partial {{\rm{\Theta }}}_{f}/\partial U}{\partial {{\rm{\Theta }}}_{f}/\partial {a}_{f}}=\displaystyle \frac{{a}_{f}}{{U}_{P}-{B}_{f}}\left[1-\displaystyle \frac{{a}_{n}}{{a}_{f}}\displaystyle \frac{{I}_{f}}{({X}_{n}-\overline{{X}_{n}}){\rho }_{{Fx}}({X}_{f})}\right].\end{eqnarray} \tag{ 23 }$$

Figures 2(a) and (b) show the level density parameter at the fission saddle point and the Γ_f/Γ_n ratio for several nuclei. In figure 2(a) the behavior of a_f in the low energy limit is determined by the relationship between the neutron separation energy B_n and the fission barrier B_f, obtaining that (see the appendix)

$$\begin{eqnarray}&&\mathop{{\rm{lim}}}\limits_{U\to {B}_{n}^{+}}\displaystyle \frac{{{\rm{d}}{a}}_{f}}{{\rm{d}}{U}}=-\mathop{{\rm{lim}}}\limits_{U\to {B}_{n}^{+}}\displaystyle \frac{2}{(U-{B}_{f})}\displaystyle \frac{{\int }_{0}^{{X}_{f}}{\rho }_{{Fx}}(x){\rm{d}}{x}}{(U-{B}_{n}){\rho }_{{Fx}}({X}_{f})}\lt 0,\quad ({B}_{n}\gt {B}_{f}),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\mathop{{\rm{lim}}}\limits_{U\to {B}_{f}^{+}}\displaystyle \frac{{{\rm{d}}{a}}_{f}}{{\rm{d}}{U}}=\mathop{{\rm{lim}}}\limits_{U\to {B}_{f}^{+}}\displaystyle \frac{{a}_{f}}{(U-{B}_{f})}\gt 0\qquad ({B}_{f}\gt {B}_{n}),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\mathop{{\rm{lim}}}\limits_{U\to {B}_{f/n}^{+}}\displaystyle \frac{{{\rm{d}}{a}}_{f}}{{\rm{d}}{U}}=-\mathop{{\rm{lim}}}\limits_{U\to {B}_{f/n}^{+}}\displaystyle \frac{{a}_{f}}{(U-{B}_{f})}\lt 0\qquad ({B}_{f}\approx {B}_{n}),\end{eqnarray} \tag{ 26 }$$

where B_f/n is the highest value between B_n and B_f.

Zoom In Zoom Out Reset image size

Equations (24)–(26) were obtained considering that for values of X_n and X_f less than one the level density becomes constant [8]. These equations, together with the information presented in table 1 allow to say that for those nuclei with B_f > B_n the values of a_f decrease when $U\to {B}_{f}^{+}$ (da_f/dU > 0). For ²³⁸U, although B_f > B_n, it also happens that B_f ≈ B_n; then a mixed behavior occurs, initially da_f/dU < 0 and a_f grows while $U\to {B}_{f}^{+}$ reaching a maximum, after which an extremely abrupt fall begins (not shown in the figure). Finally, for ²³²Pu the value of B_n is appreciably greater than B_f, then da_f/dU < 0 and a_f grows while $U\to {B}_{n}^{+}$, being the growth much faster than in the case of ²³⁸U.

Shell effects are also observable in figure 2(a), which shows that at high energies a_f seems to be directly proportional to the mass of the nucleus, while below 80 MeV the proximity to the closed shell reduces these values, in such a way that for a given energy, the value of a_f for ²⁰⁸Pb is the lowest when compared with ¹⁹⁷Au and ¹⁸⁴W. In the same way, the ratio Γ_f/Γ_n for ²⁰⁸Pb is always the lowest in figure 2(b).

Other calculations of a_f and Γ_f/Γ_n can be found in [13, 15, 16] and references therein. In [15] the general behavior of a_f for ²⁴⁰Pu is similar to the curve for ²³²Pu presented here in figure 2(a). The same happens with the ratio Γ_f/Γ_n for ²⁰²Pb in [15] and ²⁰⁸Pb in figure 2(b), with the particularity, that the values corresponding to ²⁰⁸Pb are much lower, which could be related to its doubly magical character. References [13, 16] offer data on the level density parameters as well as on the models to calculate them. The level densities at the fission saddle points are calculated according to the microscopic HFBCS model [17], (which includes shell, pairing and deformation effects related to the fission barrier structures in a consistent manner) and can be quantified in terms of the ratio a_f/a between the level density parameter at the fission saddle point and the level density parameter of the equilibrium ground-state configuration.

The use of a microscopic model of the level densities should lead to a more consistent description of the saddle points than that offered by phenomenological models such as the BSFG [13]. Then it is not surprising that when making a comparison between the results of both models some discrepancies are observed, mainly at low energies. Such is the case of figure 3, which shows the ratio a_f/a for ²³⁸U reported in [13] together with the one obtained in the present work. At high energies, the results of both models should be the same. It is important to note that the model used in [13] does not consider the possibility of the damping of the collective effects at increasing excitation energies, which could lead to a greater than expected difference between the models.

Zoom In Zoom Out Reset image size

CRISP (Colaboração Río-Ilhéus-São Paulo) [18] is a computational code to simulate nuclear reactions at intermediate energies (40 MeV–4 GeV). It is composed by several modules dedicated to the study of the different stages in which a nuclear reaction is divided (INC [19], evaporation–fission competition [20] and fission fragments distribution [21, 22]). There is also possible to study multi-fragmentation reactions. All modules are under continuous study and development [23–29].

The statistical module of CRISP works with an excited nucleus coming from the INC, calculating several magnitudes, such as the separation energies of the light particles to be emitted (neutrons, protons and alpha particles), the Coulomb and the fission barrier height, the level density parameters corresponding with the initial nucleus, the fission saddle point and the residual nucleus after emission. Then the probabilities of fission and light particles' emission are calculated, and one of the outputs is randomly selected. In the case of fission, the process stops, saving the data of the fissioning nucleus. For light particles' emission, the program recalculates all the magnitudes. This process continues until the excitation energy of the residual nucleus is below the neutron binding energy, when events such as gamma emission become competitive.

Although it is possible to find a lot of information about the code in the references, it is useful to highlight some of its characteristics. In this way, for nuclei such that $78\leqslant Z\leqslant 120$ the fission barrier height is calculated using data reported at the Reference Input Parameter Library-2 (RIPL2) [16], whereas that the following parameterization is valid when $68\leqslant Z\lt 78$

$$\begin{eqnarray}&&{B}_{f}={f}_{0}{E}_{S}+C(A){{\rm{e}}}^{-0.40{(Z-{Z}_{0})}^{2}}.\end{eqnarray} \tag{ 27 }$$

Where ${Z}_{0}=A/(1.98+0.015{A}^{2/3})$. This parameterization is based on the work of Möller and Nix [30, 31]. Here, f₀ is a function of the critical fissility and also other parameters that take into account the competition between nuclear and Coulomb forces, E_S defines the surface energy of the nucleus, and C(A) is introduced to include the shell effects. In the case Z < 68 fission is not considered.

The code includes two different forms of the level density: one very similar to the one used by Dostrovsky [5–7] and another based on the BSFG model, presented by Egidy and Bucuresco [8], in which the nuclear masses allows the calculation of the level density parameter in two ways: independent and dependent on the excitation energy. In this work, for the case of BSFG level density, was adopted the variant in which the level density parameter is independent of the excitation energy.

For Dostrovsky parameterization (DP), the level density parameter at the fission saddle point is calculated according to a_f = r_f a_n where r_f is a dimensionless parameter depending on Z²/A. On the other hand, a_f for BSFG parameterization (BsP) is determined using equation (23). Coulomb barrier height for DP comes from

$$\begin{eqnarray}&&{V}_{D}={V}_{0}\left(1-\displaystyle \frac{U}{{B}_{{\rm{lig}}}}\right),\end{eqnarray} \tag{ 28 }$$

where B_lig is the nuclear binding energy. Whereas that in BsP the Coulomb barrier height is calculated using equation (13), being also possible to link equation (28) with BsP.

The excitation energy after $n+1$ steps of the evaporation chain is

$$\begin{eqnarray}&&{E}_{n+1}={E}_{n}-{B}_{l}-\epsilon ,\end{eqnarray} \tag{ 29 }$$

where is the kinetic energy of the particle, which is obtained attending to the level density parameterization. In the case of DP, is determined considering a Maxwell distribution [10]. On the other hand, to be consistent with the changes made, a different way of calculation of the kinetic energy was adopted for BsP. First, the program generates a random number δ. Then, the integral given by equation (5) is calculated, keeping the same upper limit and changing the lower limit to ${X}_{{\rm{low}}}={a}_{l}(U-{E}_{{Fl}}-{B}_{l}-{E}_{m})$, where E_m is the maximum value for the kinetic energy of the emitted particle, given by the condition [10] ${E}_{m}\ll U-{E}_{{Fl}}-{B}_{l}$. In the present work it was considered that ${E}_{m}={V}_{l}+0.03(U-{E}_{{Fl}}-{B}_{l}-{V}_{l})$. At this moment an iterative process starts, in which the integral of equation (5) is evaluated, increasing the lower limit x_lw, until the moment when ${I}_{l}({x}_{{lw}},{X}_{l})/{I}_{l}({X}_{{\rm{low}}},{X}_{l})=\delta$. Finally, the kinetic energy is determined from the condition

$$\begin{eqnarray}&&\epsilon =U-{E}_{{Fl}}-{B}_{l}-\displaystyle \frac{{x}_{{lw}}}{{a}_{l}}.\end{eqnarray} \tag{ 30 }$$

In order to know the influence of the new parameterization, equations (13), (23) and (30) were implemented into the CRISP code, and used to simulate the reactions $p+{}^{208}{\rm{Pb}}$ and $p+{}^{238}{\rm{U}}$ at 1000 MeV. Firstly, the INC, whose resultant nuclei were used as entry of the evaporation phase. Figure 4 shows some general features of the nuclei at the beginning of evaporation for these reactions.

Zoom In Zoom Out Reset image size

Figures 5(a) and (b) show the average values of the number of emitted neutrons, protons and alpha particles as a function of the excitation energy of the initial residual nucleus, which correspond to three parameterizations of CB: V₀, with no energy correction, V_D, given by equation (28) [6, 7], and V_N, obtained from equation (13). In both cases, the number of neutrons remains more or less the same, no matter the parameterization, partly because the barrier height for neutrons is always zero. V₀ and V_D parameterizations lead to a slightly higher number of neutrons until energies around 500 MeV when V_N starts to give higher values. For protons, V_N offers the highest multiplicity, while V₀ and V_D give similar results. In the case of the alpha particles, V_N also leads to the highest values. So the use of the V_N parameterization produces the highest values of the average number of emitted protons and alpha particles, what is in agreement with the predictions made in section 2. Table 2 confirms these ideas, showing the multiplicity of the particles emitted during the evaporation.

Zoom In Zoom Out Reset image size

The content of table 3 is the same as in table 2, but considering fission during the simulations. It is possible to observe how fission reduces the number of emitted particles (indeed, fission represents a sudden cut in the evaporation chain), mainly in the case of ²³⁸U.

Table 4 and figure 6 show data corresponding with the kinetic energy distribution of the emitted particles. It is possible to observe that, for neutrons, the distributions are very similar; nevertheless, the values given by V_N are slightly lower. In the case of protons, the mean values of V₀ and V_D are very close and higher than V_N. Finally, the mean value of V₀ for alpha particles is the highest, V_N the lowest and V_D is between them.

Zoom In Zoom Out Reset image size

This information, together with the one given in table 2 and figures 5(a) and (b) allows establishing that, the higher the multiplicity of one kind of particle, the lower the mean value of the kinetic energy associated to it when emitted. That's why V_N, who gives the highest values of the multiplicity for protons and alphas, leads to the lowest ones for the mean kinetic energy. Then, the almost constant value of the neutron multiplicity occurs because any change in the charged particles' multiplicity, is directly related with a variation on the kinetic energy of this kind of particles, without affecting the neutrons' behavior.

Figures 7 and 8 show the isotopic distribution of the spallation products, obtained using both the previous (V_D-dostrovsky) and the new (V_N-BSFG) parameterization included in the CRISP code, which are compared with experimental data [32, 33]. Here the junction V_N-BSFG determines CB from equation (13), a_f from equation (23) and the level density parameter of the nucleus before and after the emission is calculated according to [8]; while the coupling V_D-dostrovsky obtains V_C from equation (28), a_f = r_fa_n and the level density parameter of the nucleus before and after the emission comes from [6, 7]. The combination V_N-BSFG offers the best results, mainly for the reaction $p+{}^{238}{\rm{U}}$ at 1000 MeV. The pair V_D-dostrovsky does not correctly reproduce the experimental data because it overestimates the occurrence of fission, which immediately leads to an underestimation of the isotopic distribution of the spallation products.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The combination of the BSFG model for the level density with the partial widths for fission and particles' emission allowed the creation of two new parameterizations, one for the Coulomb barrier and another for the level density parameter at the fission saddle point. These parameterizations are explicitly dependent on the excitation energy and the level density parameter of the nucleus after the emission of light particles. As a result, the number of light charged particles emitted is higher, and the values of a_f reflect the influence of the shell effects at low energies.

The work performed in section 2 lies in the framework of the Weisskopf–Ewing formalism (WE) [11] for the emission of light particles. Then a question arises on whether to use or not this formalism, knowing that the Hauser–Feshbach (HF) formalism [34] allows a more rigorous treatment of angular momenta and spins of nuclei and outgoing particles involved in a nuclear reaction.

It happens that, despite offering a more rigorous treatment of angular momentum, HF is not the most common choice between the codes devoted to studying spallation [3] (ABLA [35–38], CRISP [5, 18–28], SMM [39], EVAP [7, 40–42], GEM [43, 44]), where WE is used instead. A remarkable exception to this is the GEMINI code [35, 45].

The origin of this behavior lies in the high computational cost of the implementation of the HF formalism (simulations are pretty much longer than the WE equivalents) as well as in the fact that the intrinsic angular momentum of the compound nucleus formed in spallation reactions is typically low in comparison with fusion reactions. Indeed, the original version of GEMINI studied fusion, that is why HF remained in the new version for studying spallation (GEMINI++) [35].

There are several studies in which it is possible to compare different de-excitation models, such is the case of [46], where GEMINI++ (uses HF for evaporation) and SMM (works with WE) used data coming from the intranuclear-cascade model INCL4.5 [47] for simulating several reactions, including spallation. Both combinations showed strengths and weaknesses, although the INCL4.5/GEMINI++ model gave better general results when compared with experimental data.

Nevertheless, it was not possible to identify the formalism for evaporation as the only source of the discrepancies. Other possible causes are

• the parameterization of some ingredients of the models, such as nuclear level densities, Coulomb barriers, and collective enhancement;
• the proper adjustment of the free parameters included in the codes;
• the inclusion or not of a pre-equilibrium stage between cascade and de-excitation;

whose influence will depend on the observable under study. Then the origin of the discrepancies lies in all the set of models that make up a computer code, so the individual influence of each model is overlapped, being quite unlikely that just one of them is responsible for the differences observed when comparing with experimental data, which is common in the codes studying nuclear reactions.

Some of these ideas are confirmed by figures 9 and 10, which compare the new parameterization of the statistical module of CRISP with other de-excitation models that were coupled to the INCL4.5 code to simulate the reactions. The data corresponding to the different couplings come from the IAEA Benchmark of Spallation Models [48], which offers extensive theoretical and experimental information on spallation reactions. The junction ${incl}45-{abla}07$ gives the best results. However, there are regions where all models reasonably describe the experimental data and others where they diverge, showing that all codes can be improved [48, 49]. Thereby the description of the isotopic distribution of the spallation products offered by CRISP is in correspondence with that given by most of the de-excitation models, confirming the reliability of the new parameterizations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At this point, it is important to stress again that the objective of this work is to study the relationship between the level density and magnitudes such as CB and a_f. So that to implement the new models only the sections of code directly related to them were modified, without changing any other parameters, so no attempt was made to achieve better agreement with the experimental data, which could be the goal of future works.

On the other hand, it is also possible to get information from the general expression corresponding to the probability per unit time W()d that a nucleus P with excitation energy U, emits a light particle l with kinetic energy between and  + d, resulting in a nucleus F with excitation energy (U − B_l − ) [10]. Integrating this expression for all the possible values of results in the partial width corresponding to the emission of the particle l:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{l}={\int }_{{V}_{l}}^{U-{B}_{l}}\sigma (U,{V}_{l},\epsilon ...)\displaystyle \frac{{g}_{l}{m}_{l}\epsilon {\rho }_{F}}{{\pi }^{2}{{\hslash }}^{2}{\rho }_{P}}\,{\rm{d}}\epsilon ,\end{eqnarray} \tag{ 31 }$$

where $\sigma (U,{V}_{l},\epsilon ...)$ is the cross section for the reverse process.

For incorporating equation (31) into the codes is necessary to determine σ(U, V_l, ...), whose form will depend on the formalism for studying evaporation (WE or HF) and will always include U, V_l and , so that the partial width for light particle's emission will have the form Γ_l = Γ_l(U, V_l, ...) while the ratio Θ = Γ_l/Γ_n will be of the form Θ = Θ(U, V_l ...), leading to the results of section 2.

The final equations would vary, but the main conclusion remains the same: there is a relationship between CB and the level density. The same line of thinking applies to the level density parameter at the fission saddle point.

The inclusion of the HF formalism would help to determine which observables are sensitive to a rigorous treatment of the conservation of both the angular momentum and the parity. Currently, this is not possible in CRISP, so one of the future objectives would be its implementation, which is, in fact, a demanding work.

The new methods for calculation of V_C and a_f were implemented into the CRISP code and used in the simulation of spallation nuclear reactions, leading to an acceptable description of the isotopic distribution of the spallation products of the reactions p + ²⁰⁸Pb and p + ²³⁸U at 1000 MeV. Nevertheless, some discrepancies appeared. In the case of p + ²⁰⁸Pb, there is a displacement between the theoretical and the experimental spallation parabolas; as well as in the case of p + ²³⁸U, which offers poorer results for nuclei with an atomic number near to 92. The source of the differences lies in the causes previously exposed, especially the treatment of collective effects [13].

The BSFG was the base to develop the models, although it is possible to adapt them to any other level density parameterization, which is currently being doing with DP as well as with the variant of BSFG with the level density parameter dependent on the excitation energy [8]. Then it is possible to say that the main result of the present work is to show that there is a relationship between the level density and magnitudes such as CB and the level density parameter at the fission saddle point, which was not studied before.

O Tumbarell would like to thank CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior) for the financial support. F Guzmán would like to thank the SBF (Sociedade Brasileira de Física) for the financial support of PLAF/SBF/CNPq in the frame of the Cátedra Itinerante program. And finally, the authors thank the financial support of Conselho Nacional de Desenvolvimento Científico e Tecnológico(CNPq).

A.1. Deduction of equation (2)

The level density model used in this work is the BSFG, defined as [50]

$$\begin{eqnarray}&&{\rho }_{P}(U)=\displaystyle \frac{\exp (2\sqrt{a(U-{E}_{P})})}{12\sqrt{2}\sigma {a}^{1/4}{(U-{E}_{P})}^{5/4}},\end{eqnarray} \tag{ 32 }$$

where U is the excitation energy, E_P the energy correction, a the level density parameter and σ the spin-cutoff parameter, given by [8]

$$\begin{eqnarray}&&\sigma =0.121{A}^{5/6}\displaystyle \frac{{[1+\sqrt{(1+4a(U-{E}_{P})}]}^{1/2}}{{[2a(U-{E}_{P})]}^{1/2}},\end{eqnarray} \tag{ 33 }$$

then, equation (32) becomes

$$\begin{eqnarray}&&{\rho }_{P}(U)=\displaystyle \frac{0.69a\exp (2\sqrt{{{aU}}_{P}})}{{A}^{5/6}{({{aU}}_{P})}^{3/4}{[1+\sqrt{1+4{{aU}}_{P}}]}^{1/2}},\end{eqnarray} \tag{ 34 }$$

with U_P = U − E_P. Using (34), equation (1) takes the form

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{l}=\gamma {\int }_{{V}_{l}}^{{U}_{F}-{B}_{l}}\displaystyle \frac{(\epsilon -{V}_{l}){a}_{l}\exp (2\sqrt{{a}_{l}{U}_{{Fl}}})\,{\rm{d}}\epsilon }{{({a}_{l}{U}_{{Fl}})}^{3/4}{[1+\sqrt{1+4{a}_{l}{U}_{{Fl}}}]}^{1/2}},\end{eqnarray} \tag{ 35 }$$

being

$$\begin{eqnarray}&&\gamma =\displaystyle \frac{0.69{g}_{l}{m}_{l}{r}_{0}^{2}{A}^{2/3}}{\pi {{\hslash }}^{2}{A}_{l}^{5/6}{\rho }_{P}({U}_{P})},\end{eqnarray} \tag{ 36 }$$

A_l the residual nucleus' mass, ${U}_{F}=U-{E}_{{Fl}}$ and ${U}_{{Fl}}={U}_{F}-{B}_{l}-\epsilon$, while a_l and E_Fl are the level density parameter and the energy correction after emission of the particle l (using the subscripts n, p and a for neutrons, protons and alpha particles, respectively). The substitutions $x={a}_{l}(U-{E}_{{Fl}}-{B}_{l}-\epsilon )$ and ${X}_{l}={a}_{l}(U-{E}_{{Fl}}-{B}_{l}-{V}_{l})$ transforms (35) into equation (2).

A.2. Deduction of equation (6)

Considering Θ = Θ(U, V_l) its partial derivatives are

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Theta }}}{\partial U}={\rm{\Theta }}\left[\displaystyle \frac{1}{{I}_{l}}\displaystyle \frac{\partial {I}_{l}}{\partial U}-\displaystyle \frac{1}{{I}_{n}}\displaystyle \frac{\partial {I}_{n}}{\partial U}\right],\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Theta }}}{\partial {V}_{l}}={\rm{\Theta }}\left[\displaystyle \frac{1}{{I}_{l}}\displaystyle \frac{\partial {I}_{l}}{\partial {V}_{l}}\right],\end{eqnarray} \tag{ 38 }$$

which leads to

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Theta }}/\partial U}{\partial {\rm{\Theta }}/\partial {V}_{l}}=\displaystyle \frac{\left[\tfrac{\partial {I}_{l}}{\partial U}-\tfrac{{I}_{l}}{{I}_{n}}\tfrac{\partial {I}_{n}}{\partial U}\right]}{\tfrac{\partial {I}_{l}}{\partial {V}_{l}}}.\end{eqnarray} \tag{ 39 }$$

Besides,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {I}_{l}}{\partial U}={a}_{l}{\int }_{0}^{{X}_{l}}{\rho }_{{Fx}}(x){\rm{d}}{x},\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {I}_{l}}{\partial {V}_{l}}=-\displaystyle \frac{\partial {I}_{l}}{\partial U},\end{eqnarray} \tag{ 41 }$$

and $\tfrac{\partial {I}_{n}}{\partial U}$ is the same as equation (40), changing l for n. Therefore,  (39) turns into equation (6).

A.3. Deduction of equation (9)

Let Θ be a function such that Θ = Θ(U, V_l), then, considering V_l = V_l(U) the total derivative of Θ is

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\rm{\Theta }}}{{\rm{d}}{U}}=\displaystyle \frac{\partial {\rm{\Theta }}}{\partial U}+\displaystyle \frac{\partial {\rm{\Theta }}}{\partial {V}_{l}}\displaystyle \frac{{{dV}}_{l}}{{\rm{d}}{U}},\end{eqnarray} \tag{ 42 }$$

leading to

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{V}}_{l}}{{\rm{d}}{U}}=\displaystyle \frac{\tfrac{{\rm{d}}{\rm{\Theta }}}{{\rm{d}}{U}}}{\partial {\rm{\Theta }}/\partial {V}_{l}}-\displaystyle \frac{\partial {\rm{\Theta }}/\partial U}{\partial {\rm{\Theta }}/\partial {V}_{l}},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{V}}_{l}}{{\rm{d}}{U}}=\displaystyle \frac{\partial {\rm{\Theta }}/\partial U}{\partial {\rm{\Theta }}/\partial {V}_{l}}+\displaystyle \frac{\tfrac{{\rm{d}}{\rm{\Theta }}}{{\rm{d}}{U}}}{\partial {\rm{\Theta }}/\partial {V}_{l}}-2\displaystyle \frac{\partial {\rm{\Theta }}/\partial U}{\partial {\rm{\Theta }}/\partial {V}_{l}},\end{eqnarray} \tag{ 44 }$$

where the last two terms of the right side of equation (44) are functions of U and V_l. Then it is possible to establish that

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{V}}_{l}}{{\rm{d}}{U}}=\displaystyle \frac{\partial {\rm{\Theta }}/\partial U}{\partial {\rm{\Theta }}/\partial {V}_{l}}+{\rm{\Xi }}(U,{V}_{l}).\end{eqnarray} \tag{ 45 }$$

Once the function ${\rm{\Xi }}(U,{V}_{l})$ has been estimated, V_l is calculated by numerical integration of the equation (45). Following a similar path to the one exposed before it is possible to obtain the equations corresponding to the level density parameter at the saddle point.

A.4. Deduction of equations (24)–(26)

Explicitly, equation (23) has the form

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{a}}_{f}}{{\rm{d}}{U}}=\displaystyle \frac{{a}_{f}}{U-{B}_{f}}\left[1-\displaystyle \frac{{a}_{n}}{{a}_{f}}\displaystyle \frac{{\int }_{0}^{{X}_{f}}{\rho }_{{Fx}}(x){\rm{d}}{x}}{\left({X}_{n}-\tfrac{{\int }_{0}^{{X}_{n}}x{\rho }_{{Fx}}(x){\rm{d}}{x}}{{\int }_{0}^{{X}_{n}}{\rho }_{{Fx}}(x){\rm{d}}{x}}\right){\rho }_{{Fx}}({X}_{f})}\right],\end{eqnarray} \tag{ 46 }$$

with ρ_Fx(x) given by equation (3). For the sake of simplicity, the value of the energy correction was included inside B_f/n.

When $U\to {B}_{f/n}^{+}$, ${X}_{f/n}={a}_{f/n}(U-{B}_{f/n})\to 0\lt 1$. For ${X}_{f/n}\lt 1$ it is possible to assume that the level density takes a constant value [8], then the function ρ_Fx takes the constant value ρ_Fx(X_f/n) and it can be extracted from the corresponding integral in (46). Here B_f/n is the highest value between B_f and B_n, while ${X}_{f/n}={X}_{f},{X}_{n}$ and a_f/n = a_f, a_n; depending on B_f/n.

If B_n > B_f, ${X}_{n}\to {0}^{+}$ when $U\to {B}_{n}^{+}$ , while X_f > 0. Taking into account the previous paragraph equation (46) becomes

$$\begin{eqnarray}&&\mathop{{\rm{lim}}}\limits_{U\to {B}_{n}^{+}}\displaystyle \frac{{{\rm{d}}{a}}_{f}}{{\rm{d}}{U}}=\mathop{{\rm{lim}}}\limits_{U\to {B}_{n}^{+}}\displaystyle \frac{{a}_{f}}{(U-{B}_{f})}\left[1-2\displaystyle \frac{{\int }_{0}^{{X}_{f}}{\rho }_{{Fx}}(x){\rm{d}}{x}}{{a}_{f}(U-{B}_{n}){\rho }_{{Fx}}({X}_{f})}\right],\end{eqnarray} \tag{ 47 }$$

where the absolute value of the second term of the right side of equation (47) is much higher than the first, leading to

$$\begin{eqnarray}&&\mathop{{\rm{lim}}}\limits_{U\to {B}_{n}^{+}}\displaystyle \frac{{{\rm{d}}{a}}_{f}}{{\rm{d}}{U}}=-\mathop{{\rm{lim}}}\limits_{U\to {B}_{n}^{+}}\displaystyle \frac{2}{(U-{B}_{f})}\displaystyle \frac{{\int }_{0}^{{X}_{f}}{\rho }_{{Fx}}(x){\rm{d}}{x}}{(U-{B}_{n}){\rho }_{{Fx}}({X}_{f})}\lt 0,\quad ({B}_{n}\gt {B}_{f})\end{eqnarray} \tag{ 48 }$$

which is the equation (24). In a similar way it is possible to obtain the expressions given in equations (25) and (26).