text:
# 1 Introduction

## 1 Introduction

Neutrinos play an important role during Big Bang Nucleosynthesis (BBN) due to their contribution to the total energy density. Furthermore, the electron neutrinos directly influence the neutron-proton reactions. In this context, the Standard Model with its three active massless neutrinos agrees very well with the observations of the primordial light element abundances: a recent analysis  of observational data claims that any particles beyond the Standard Model must not contribute more than 20% to the total energy density of a single neutrino flavor.

The recent discovery of $`\nu _\mu `$-$`\nu _\tau `$ oscillations at the Super-Kamiokande detector  has proven the existence of neutrino masses. Reasonable neutrino masses of order $`m_\nu 1\mathrm{eV}`$ will not spoil the standard BBN picture, because BBN takes place at temperatures $`T0.1\mathrm{MeV}m_\nu `$. Also, mixing between the active neutrino flavors plays no role since they are equally populated during the relevant epoch.

But the situation dramatically changes if one supposes the existence of one or more sterile neutrinos: even a partial excitation of more than 20% of a sterile neutrino state, e.g. via oscillations, will contradict the observational data of light element abundances. Actually, introducing sterile neutrinos is quite natural, since right-handed neutrinos would necessarily be sterile in the Standard Model sense!

Another important place in which sterile neutrinos can play a role is the creation of lepton asymmetry in the active neutrino sector. For example, if a significant asymmetry between $`\nu _\mathrm{e}`$ and $`\overline{\nu }_\mathrm{e}`$ is established, this could change the neutron-proton ratio and lead to a contradiction with the observational data.

Early considerations of the simplified kinetic equations found no significant increase of the asymmetry in the case where the sterile neutrino has a mass $`m_{\nu _\mathrm{s}}1\mathrm{MeV}`$ and any vacuum mixing angle and concluded that this asymmetry was always small, $`\eta 1`$ . However, later it was found that in the region of small mixing angles the asymmetry can increase significantly and reach values $`𝒪(0.1)`$ . This statement was confirmed by numerical calculations  where a large final value of the asymmetry was found for the first time. Recently, more accurate numerical calculations also confirmed the increase of asymmetry to large values $`𝒪(0.1)`$ . The reason for such a big increase of the asymmetry is that it grows according to a power law close to $`T^4`$.

Contrary to this statement is the result of an analytical approach by Dolgov et al. , who found an additional counterterm in the effective equation for the asymmetry, such that the power law would change from $`T^4`$ to $`T^1`$. Thus the final asymmetry would be $`𝒪(10^4)`$.

In this paper, we analyze the method of  and show that it is valid only when the evolution of the kinetic equations dictated by the collision terms, i.e. when the temperature of the primordial plasma is large enough, $`T1\mathrm{MeV}`$. For smaller temperatures, $`T`$ 2-3 MeV, collisions become less important, and the evolution of the neutrinos can be described by the well-known Mikheyev-Smirnov-Wolfenstein (MSW) effect . In this small temperature region the approach of  is not valid. This effect has been studied in .

The counterterm found in , which leads to the power law $`T^1`$, only becomes compatible to the main term when the applied perturbative expansion breaks down. This fact suggests that the power law $`T^1`$ could be wrong, but strictly speaking it does not prove this. Therefore, we develop a new independent analytical method which can be applied in the regions where the method of  does not work. Instead of expanding in the coefficients of the kinetic equations, we solve these equations pertubatively for the case when the creation rate of sterile neutrinos is not significantly large. Using this new method we show that the $`T^1`$ power law is not a self-consistent solution for the lepton asymmetry. Instead, the power law $`T^{11/3}`$ is approximately valid, and therefore the asymmetry can reach large values of order $`0.1`$, although our new method in its present form does not allow us to calculate the final value of the asymmetry.

In section 2 we define the kinetic equations. Then in section 3 we shortly discuss the evolution of the lepton asymmetry before it reaches the power-law regime and explain the power-law solution known in literature. In section 4 we discuss some recent papers on this topic, especially the one claiming the $`T^1`$ power law . In section 5 we present an improved approach similar to  and discuss its validity. In section 6 we present a new analytic solution of the kinetic equations which we use to check the validity of the $`T^1`$ power law. In section 7 we conclude and summarize our results. In the Appendix we comment on some numerical fits.

## 2 Kinetic equations

In order to describe the evolution of the system of oscillating active $`\nu _\mathrm{a}`$, $`\mathrm{a}=\mathrm{e},\mu ,\tau `$, and sterile neutrinos $`\nu _\mathrm{s}`$ with a mass difference $`\delta m^2`$ and a vacuum mixing angle $`\theta `$ in the Early Universe we need to consider the density matrix formalism  and take into account second order effects in the Fermi coupling constant . Then the following set of equations for the components of the density matrix can be found:

$`i(_tHp_p)\rho _{aa}`$ $`=`$ $`F_0(\rho _{sa}\rho _{as})/2i\mathrm{\Gamma }_0(\rho _{aa}f_{\mathrm{eq}}),`$
$`i(_tHp_p)\rho _{ss}`$ $`=`$ $`F_0(\rho _{sa}\rho _{as})/2,`$
$`i(_tHp_p)\rho _{as}`$ $`=`$ $`W_0\rho _{as}+F_0(\rho _{ss}\rho _{aa})/2i\mathrm{\Gamma }_1\rho _{as},`$
$`i(_tHp_p)\rho _{sa}`$ $`=`$ $`W_0\rho _{sa}F_0(\rho _{ss}\rho _{aa})/2i\mathrm{\Gamma }_1\rho _{sa},`$ (1)
where $`a`$ and $`s`$ mean “active” and “sterile” respectively, $`F_0=\delta m^2\mathrm{sin}2\theta /2E`$, $`W_0=\delta m^2\mathrm{cos}2\theta /2E+V_{\mathrm{eff}}^a`$, $`H=\sqrt{8\pi \rho _{\mathrm{tot}}/3m_{\mathrm{pl}}^2}`$ is the Hubble parameter, $`m_{\mathrm{pl}}=1.22\times 10^{22}`$ MeV is the Planck mass. Furthermore, since the mass of the neutrinos is negligible in comparison with their momentum, the neutrino momentum and energy equal, $`p=E`$. Finally, $`f_{\mathrm{eq}}`$ is the equilibrium Fermi distribution function, given by

$$f_{\mathrm{eq}}=\frac{1}{e^{(E\mu )/T}+1},$$
(2)
where $`T`$ is the photon temperature and $`\mu `$ is the chemical potential for the active neutrinos.

The effective potential for (anti)-neutrinos is

$$V_{\mathrm{eff}}^a=C_1\eta G_FT^3+C_2^a\frac{G_F^2T^4E}{\alpha },$$
(3)
where $`E`$ is the neutrino energy, $`T`$ is the temperature of the plasma, $`G_F=1.166\times 10^5`$ GeV-2 is the Fermi coupling constant, $`\alpha =1/137`$ is the fine structure constant, $`C_1=0.345`$, and $`C_2^\mathrm{e}=0.61`$ for $`\nu _\mathrm{e}`$-$`\nu _\mathrm{s}`$ mixing and $`C_2^{\mu ,\tau }=0.17`$ for $`\nu _{\mu ,\tau }`$-$`\nu _\mathrm{s}`$ mixing. The signs “$``$” refer to neutrinos and anti-neutrinos, respectively. The individual contributions to the effective asymmetry $`\eta `$ from different particle species are

$`\eta `$ $`=`$ $`2\eta _{\nu _\mathrm{e}}+\eta _{\nu _\mu }+\eta _{\nu _\tau }+\eta _\mathrm{e}\eta _\mathrm{n}/2\text{for}\nu _\mathrm{e},`$
$`\eta `$ $`=`$ $`2\eta _{\nu _\mu }+\eta _{\nu _\mathrm{e}}+\eta _{\nu _\tau }\eta _\mathrm{n}/2\text{for}\nu _\mu ,`$ (4)
and $`\eta `$ for $`\nu _\tau `$ is obtained from Eq. (4) by the interchange $`\mu \tau `$.

Equations (1) account exactly for the first-order terms described by the refractive index, while the second order terms describing the coherence breaking are approximately modelled by the damping coefficients $`\mathrm{\Gamma }_j`$. The latter are equal to :

$$\mathrm{\Gamma }_0=2\mathrm{\Gamma }_1=g_a(p)\frac{180\zeta (3)}{7\pi ^4}G_F^2T^4p\left[1z_a\eta _{\nu _\mathrm{e}}+O(\eta _{\nu _\mathrm{e}}^2)\right],$$
(5)
where the expression in brackets stems from the chemical potential, and again “$``$” denotes $`\nu `$ and $`\overline{\nu }`$, respectively. The coefficients of the expansion in powers of the asymmetry are small, $`z_{\nu _\mathrm{e}}=0.1`$ and $`z_{\nu _{\mu ,\tau }}=0.04`$ , and in the following we will neglect these terms. The coefficients $`g_a(p)`$ are in general momentum-dependent, but in order to simplify Eqs. (1) we can take their momentum averaged values, which are $`g_{\nu _\mathrm{e}}=3.56`$ and $`g_{\nu _{\mu ,\tau }}=2.5`$ .

The anti-neutrino density matrix satisfies a set of equations similar to Eqs. (1), but with opposite sign in the antisymmetric term of $`V_{\mathrm{eff}}^a`$ and with a slight difference in the damping factors that is proportional to the lepton asymmetry.

For convenience, we will introduce new variables, following the notation of Dolgov et al.  to a certain degree, and we will restrict ourselves to treating $`\nu _\mathrm{e}`$ only. First we note that for $`1\mathrm{MeV}<T<100\mathrm{MeV}`$ the universe is radiation dominated so that the scale factor evolves as $`a(t)T^1`$. This implies $`H=\dot{T}/T`$. Therefore, if we use the comoving frame given by the variables $`x=1\mathrm{MeV}/T`$ and $`y=p/T`$, the LHS of Eqs. (1) simplify to $`ixH_x\rho _{ij}`$. Using $`\rho _{\mathrm{tot}}=\frac{\pi ^2}{30}g_{}T^4`$ with the constant $`g_{}=10.75`$, we get $`H5.5T^2/m_{\mathrm{pl}}`$.

Furthermore, we substitute

$`\rho _{aa}`$ $`=`$ $`f_{\mathrm{eq}}(y)[1+a(x,y)],`$
$`\rho _{ss}`$ $`=`$ $`f_{\mathrm{eq}}(y)[1+s(x,y)],`$
$`\rho _{as}=\rho _{sa}^{}`$ $`=`$ $`f_{\mathrm{eq}}(y)[h(x,y)+il(x,y)].`$ (6)
Next, we introduce a new time variable $`q`$ which is defined so that for $`\eta =0`$ the resonance condition, given by $`W_0=0`$, is fulfilled when $`q=y`$. We find

$$q=\sqrt{\frac{|\delta m^2|\mathrm{cos}2\theta }{2C_2^\mathrm{e}G_F^2\alpha ^1\mathrm{MeV}^6}}x^36.6\times 10^3\sqrt{|\delta m_{\mathrm{eV}}^2|\mathrm{cos}2\theta }x^3,$$
(7)
where $`\delta m_{\mathrm{eV}}^2=\delta m^2/\mathrm{eV}^2`$. For the evolution of a single momentum mode, we will use the variable $`\tau =q/y`$.

With all these modifications, the Eqs. (1) have the form

$`s^{}`$ $`=`$ $`Fl,`$
$`a^{}`$ $`=`$ $`Fl2\gamma a,`$
$`h^{}`$ $`=`$ $`Wl\gamma h,`$
$`l^{}`$ $`=`$ $`{\displaystyle \frac{F}{2}}(as)Wh\gamma l,`$ (8)
where the prime means differentiation with respect to $`\tau `$. The coefficients are derived using $`W=W_0\left(Hx\frac{\tau }{x}\right)^1`$, $`\gamma =\mathrm{\Gamma }_1\left(Hx\frac{\tau }{x}\right)^1`$, and so forth. We also split $`W`$ into its symmetric and anti-symmetric parts, so $`W=UV\eta `$. Then

$$F=Q\mathrm{tan}2\theta Q\mathrm{sin}2\theta ,\gamma =\delta /\tau ^2,U=Q(1/\tau ^21),V=V_0y^{1/3}\tau ^{4/3}.$$
(9)
The dimensionless constants are given by

$`Q`$ $`=`$ $`{\displaystyle \frac{\delta m^2\mathrm{cos}2\theta }{2yT}}\left(Hx{\displaystyle \frac{\tau }{x}}\right)^15.6\times 10^4\sqrt{|\delta m^2|\mathrm{cos}2\theta },`$ (10)
$`\delta `$ $`=`$ $`Q{\displaystyle \frac{g_{\nu _\mathrm{e}}90\zeta (3)}{7\pi ^4C_2^\mathrm{e}/\alpha }}Q/148,`$ (11)
$`V_0`$ $`=`$ $`{\displaystyle \frac{m_{\mathrm{pl}}C_1G_F\mathrm{MeV}}{5.5\times 3}}\left({\displaystyle \frac{|\delta m^2|\mathrm{cos}2\theta }{2C_2^\mathrm{e}G_F^2\alpha ^1\mathrm{MeV}^6}}\right)^{1/6}5.6\times 10^{10}(|\delta m_{\mathrm{eV}}^2|\mathrm{cos}2\theta )^{1/6}.`$ (12)
Finally, the effective asymmetry has the form:

$$\eta =\frac{1}{2\zeta (3)}𝑑yy^2f_{\mathrm{eq}}(a\overline{a})\eta _0,$$
(13)
where $`\eta _0=\eta _\mathrm{n}/2\eta _\mathrm{e}\eta _{\nu _\mu }\eta _{\nu _\tau }`$.

Note that merely $`V`$ depends directly on the momentum $`y`$. Also note that the system of Eqs. (8) is completely equivalent to Eqs. (1), i.e. we have not used any approximations yet.

## 3 Evolution of the lepton asymmetry

We now briefly discuss the different regimes of the evolution of the lepton asymmetry. Let us first sketch the reason why the lepton asymmetry can change at all. Due to the mixing of $`\nu _\mathrm{e}`$ with $`\nu _\mathrm{s}`$, the number densities of the $`\nu _\mathrm{e}`$ and $`\overline{\nu }_\mathrm{e}`$ can alter. The sterile neutrinos that are thus created have no effect on $`\eta `$. Now if $`\eta 0`$ initially, the $`\nu _\mathrm{e}`$ and $`\overline{\nu }_\mathrm{e}`$ number densities will change differently, thus changing $`\eta `$. Thus the neutrino mixing works as a backreaction on $`\eta `$.

In Fig. 2 we have plotted an example of the lepton asymmetry evolution as a function of the inverse background medium temperature for the case $`(\delta m^2,\mathrm{sin}2\theta )=`$ $`(1\mathrm{eV}^2,5\times 10^5)`$ and an initial effective asymmetry $`\eta _0=10^{10}`$. Figure 2 shows the evolution of the resonance momentum $`y_{\mathrm{res}}`$, which is given by

$$y_{\mathrm{res}}=\sqrt{C(q)^2+q^2}\pm C(q),$$
(14)
where $`C(q)=V_0\eta (q)q^{2/3}/(2Q)`$, and $`\pm `$ is valid for neutrinos and anti-neutrinos, respectively. For definiteness, we have taken $`\eta `$ to be positive.

¿From Figs. 2 and 2 we can see that there are three different stages of lepton asymmetry evolution. During the first stage the asymmetry is so small that the variable $`C(q)`$ can be neglected in Eq. (14). The neutrinos and the anti-neutrinos simultaneously pass through the resonance. Furthermore, any initial asymmetry decreases due to the negative back-reaction.

The second stage is the stage of exponential growth of the lepton asymmetry, which occurs when the back-reaction becomes positive. From the point of view of the resonance condition, the coefficient $`C(q)`$ becomes important so that the two resonance conditions for $`\nu `$ and $`\overline{\nu }`$ are driven in opposite directions (see Fig. 2). One of them (in our example $`\nu `$) quickly passes through all the momentum modes to very high momenta where the Fermi distribution is negligible, the other one ($`\overline{\nu }`$ in our example) decreases to small momenta $`y_{\mathrm{res}}<1`$. Though the duration of this stage is short in time, the asymmetry grows many orders of magnitude and reaches values of order of $`10^6`$. (For a recent analytical investigation on these two regimes, see .)

However, we should mention that the behavior of the asymmetry in the exponential regime is still unsettled. Several publications  have claimed that the asymmetry starts oscillating in this regime, thereby changing its sign. We believe that this effect appears when the momentum distribution is neglected, since then all neutrinos pass the resonance simultaneously. Still, there might be a region in the parameter space $`(\delta m^2,\mathrm{sin}2\theta )`$ where these oscillations also occur when the momentum distribution is taken into account . In the present paper, we will not further investigate this question and will assume that no such asymmetry oscillations occur.

In the third stage the lepton asymmetry grows according to a power law, $`\eta T^\alpha `$. Unfortunately, in the literature there exist different results on the exact power $`\alpha `$. The two main results are $`\alpha =4`$, see e.g. , and $`\alpha =1`$, derived in . In the case $`\alpha =4`$, the resonance momentum is constant (see Fig. 2). Since in this case the $`\overline{\nu }_\mathrm{s}`$ at the resonance momentum will soon be fully excited, the resonance momentum of the $`\overline{\nu }`$ will in reality slowly increase, thus passing through the whole spectrum. Therefore, $`\alpha \mathrm{\Gamma }>4`$. In the end, the $`\overline{\nu }_\mathrm{s}`$ sector will essentially be excited, so that $`\eta =𝒪(0.1)`$. If $`\alpha =1`$, the resonance condition passes much faster through all the momentum modes so that the final asymmetry value remains small, $`\eta 1`$. In other words, the $`\overline{\nu }_\mathrm{s}`$ are only excited to a fractional amount.

The $`\eta T^1`$ law arrives due to the counter-term in Eq. (51) of . If this counter-term is neglected, we find a different power law for the asymmetry, $`\alpha =11/3`$, a result which  also found when neglecting their counter-term. We also show this power law in Fig. 2. This power law is quite similar to the one with $`\alpha =4`$.

Evidently the final amplitude of the asymmetry strongly depends on the exact power-law dependence. It is therefore very important to find the correct answer.

## 4 Discussion of the paper by Dolgov et al.

We now discuss the treatment of the system of differential equations as done by Dolgov et al. in . For convenience, we will denote references to equations in their paper in the form (D.51). Furthermore, to avoid confusion, we will strictly use their notation in this section.

First we would like to pay attention to some minor misprints, which however do not change the results. First, the sign in Eq. (D.20) is wrong. As a consequence, the total sign of Eq. (D.51) changes. However, since the initial value of $`b_0(0)=+1`$ instead of $`1`$, the physical content of Eq. (D.51) remains unchanged.

Next, in Eqs. (D.47) and (D.48), the signs of the first terms, i.e. with the factor $`F`$, are wrong.

Also, a correct derivation of Eq. (D.50) yields111In terms of our notation, $`F`$, $`\gamma `$ and $`U`$ must be divided by $`Q`$. Furthermore, $`VZ=D`$ corresponds to our $`V\eta /Q`$. Finally $`\sigma ^2=\gamma ^2+U^2`$ and $`\stackrel{~}{\sigma }^2=\gamma ^2(1D^2/\sigma ^2)^2+U^2(1+D^2/\sigma ^2)^2`$.

$$L=10^{10}\frac{FVZ\gamma U}{\sigma ^2\stackrel{~}{\sigma }^2}b_0\left[1+\frac{F^2(\sigma ^2+D^2)}{4\sigma ^2\stackrel{~}{\sigma }^2}\left(14\gamma ^2\frac{\sigma ^2+D^2}{\sigma ^2\stackrel{~}{\sigma }^2}\right)\right],$$
(15)
where the total sign is the same as in  when taking into account the changed sign of $`Z`$. Furthermore, the second term has changed sign with respect to , and there has appeared a new term. Still, the $`T^1`$ power law holds, since the third term now acts as the counter-term. To see this, we take the square brackets in (15) at resonance, i.e. when $`D^2=U^2`$, and use the fact that $`\gamma ^2D^2`$, which is certainly true after the exponential regime. We get $`[\mathrm{}]=[1+\frac{F^2}{8\gamma ^2}(12)]`$, which is exactly what we get from an equivalent treatment of Eq. (D.50). Thus, the counter-term dominates at resonance when $`\gamma |F|/\sqrt{8}`$.

As a consequence of Eq. (20), Eq. (D.51) becomes

$`{\displaystyle \frac{1}{Z}}{\displaystyle \frac{dZ}{dq}}`$ $`=`$ $`\delta Bq^{5/3}{\displaystyle _0^{\mathrm{}}}𝑑t{\displaystyle \frac{t^4(t^21)f_{\mathrm{eq}}(tq)b_0(1/t)}{\sigma ^2\stackrel{~}{\sigma }^2}}`$ (16)
$`\times \left(1+B_1{\displaystyle \frac{\sigma ^2+D^2}{\sigma ^2\stackrel{~}{\sigma }^2}}\left[14\gamma ^2{\displaystyle \frac{\sigma ^2+D^2}{\sigma ^2\stackrel{~}{\sigma }^2}}\right]\right).`$
We also checked Eq. (D.55), and found it to be

$$\zeta ^{}=\underset{j=1,2}{}\frac{\pi B}{2}b_0\left(1/t_j\right)\frac{q^{8/3}(t_j^21)f_{\mathrm{eq}}(qt_j)}{2\zeta \sqrt{4q^{2/3}+\zeta ^2}}\left[1\frac{B_1}{4\delta ^2t_j^4}\right].$$
(17)
Still, the general physical conclusions do not change, since the missing factor in the second term $`(t_j^21)^21`$ anyway. We will show later on how to derive this type of equation, see Eq. (56ff).

We would also like to comment on the criticism of  given in . We agree that neglecting the term $`\frac{F}{2}(AS)`$ in Eq. (D.28) is not allowed when the asymmetry $`Z`$ is essentially damped to zero. However, the effect of this term is merely to slow down the decrease of $`Z`$. In fact, it is the asymmetry in chemical potential that prevents $`Z`$ from becoming arbitrarily small, see section 5.2. In any case, during the exponential regime, any information on the amplitude of $`Z`$ before the exponential regime is lost. Thus, one can safely neglect the term in the formalism of  when one is interested in the final value of the asymmetry.

As a second point,  criticized that Dolgov et al. neglected the terms $`H^{}/Q`$ and $`L^{}/Q`$ in Eqs. (D.27,D.28). To check the validity of this approximation, we calculated these two terms numerically in the formalism used by , and compared the results with the other terms in Eqs. (D.27,D.28). We found that $`H^{}/Q`$ and $`L^{}/Q`$ were indeed negligible at all times. The differing results in  appear due to the neglect of the momentum distribution. We are of the opinion that one cannot neglect the momentum distribution, since this would imply that all neutrinos were at resonance at the same time, which is simply not the case.

Third, the author of  suggested that the $`T^1`$ power law obtained in  stemmed from neglecting the influence of the lepton asymmetry on the fermion distribution function $`f_{\mathrm{eq}}`$. By expanding $`f_{\mathrm{eq}}`$ in $`Z\zeta (3)/\pi ^2\times 10^{10}\eta _{\nu _\mathrm{e}}`$, we get

$$f_{\mathrm{eq}}(Z)f_{\mathrm{eq}}(0)\left(1\pm \frac{Z}{\frac{\zeta (3)}{\pi ^2}\times 10^{10}}\frac{12\zeta (3)}{\pi ^2}\frac{1}{1+e^y}\right),$$
(18)
where $`\pm `$ denotes neutrinos and anti-neutrinos, respectively. The correction solely leads to an additional term in $`A^{}/Q`$ which is given by

$$2\gamma Z\frac{12}{1+e^y}.$$
(19)
We compare this with the first term in $`A^{}/Q`$, i.e. $`FL`$, using the leading order in Eq. (15):

$$2\gamma Z\frac{10^{10}F^2VU}{\sigma ^2\stackrel{~}{\sigma }^2}.$$
(20)
One can clearly see that the additional term can safely be neglected near the resonance in the power-law regime.

Finally, we discuss the problems we found in the analysis of . The first occurred when the authors diagonalized their set of equations for the symmetric functions, Eq. (D.33). When diagonalizing and expanding it in $`F`$, they get as the second eigenvalue

$$\mu _22\gamma +\frac{F^2\gamma }{2[\sigma ^2(VZ)^2]},$$
(21)
where $`\sigma ^2=\gamma ^2+U^2`$. Here, the denominator of the second term diverges close to the resonance in the case that $`VZ\gamma `$. We found that this condition is already fulfilled before the end of the exponential regime. Furthermore, we found that this divergence continued to cause problems in the subsequent analysis.

The second problem occurs when $`F=𝒪(\gamma )`$. Then the expansion in $`F`$ is no longer valid. Unfortunately, the counter-term in Eq. (D.50), i.e. our Eq. (15), becomes of order 1 only when $`F\gamma `$, so that one cannot say whether the counter-term is physical or simply an artifact of the expansion.

Finally, the validity of Eq. (D.49) breaks down when $`Q\gamma 1`$. Thus, assuming that the counter-term is of order 1, (D.49) is invalid for $`QF1`$, which corresponds to $`\mathrm{sin}2\theta \sqrt{|\delta m^2|\mathrm{cos}2\theta }1.8\times 10^5`$.

## 5 Improved expansion

In this section we develop a similar approach to that of , using our own notation unless quoted explicitly. However, we want to circumvent the divergent behavior that occurs in their analysis. To this end, we treat the neutrinos and anti-neutrinos separately.

The system of equations (8) can be rewritten in matrix form as $`𝒱^{}=𝒱`$:

$`\left(\begin{array}{c}s^{}\\ a^{}\\ h^{}\\ l^{}\end{array}\right)=\left(\begin{array}{cccc}0& 0& 0& F\\ 0& 2\gamma & 0& F\\ 0& 0& \gamma & W\\ F/2& F/2& W& \gamma \end{array}\right)\left(\begin{array}{c}s\\ a\\ h\\ l\end{array}\right).`$ (34)
We can apply a non-degenerate transformation

$$𝒱=𝒳.$$
(35)
For the new vector $``$ we get the equation

$$^{}=\left(𝒳^1𝒳+𝒳^1𝒳^{}\right).$$
(36)
Now, if $`𝒳`$ is a matrix, constructed from the eigenvectors of $``$, then the first term on the RHS of Eq. (36) will be diagonal. In the case where $``$ is a constant matrix, the eigenvectors will also be constants and the second term on the RHS of Eq. (36) will be zero. Then the equations for each component of $``$ will be independent of each other and could easily be solved. Then the general solution of the matrix equation will be a linear combination of the eigenvectors of $``$ with coefficients $``$.

In the general case, when $``$ is not constant, there is no straightforward way to solve the matrix equation. But in our case in some regimes we can hope that the second term on the RHS of Eq. (36) will be a small correction in comparison to the main term, and then we can solve Eq. (36) perturbatively.

To begin, we write down the exact eigenvalues of $``$:

$`\mu _1`$ $`=`$ $`\gamma {\displaystyle \frac{1}{\sqrt{2}}}\sqrt{F^2W^2+\gamma ^2+\sqrt{4F^2\gamma ^2+(F^2+W^2+\gamma ^2)^2}},`$
$`\mu _2`$ $`=`$ $`\gamma +{\displaystyle \frac{1}{\sqrt{2}}}\sqrt{F^2W^2+\gamma ^2+\sqrt{4F^2\gamma ^2+(F^2+W^2+\gamma ^2)^2}},`$
$`\mu _3`$ $`=`$ $`\gamma {\displaystyle \frac{1}{\sqrt{2}}}\sqrt{F^2W^2+\gamma ^2\sqrt{4F^2\gamma ^2+(F^2+W^2+\gamma ^2)^2}},`$
$`\mu _4`$ $`=`$ $`\gamma +{\displaystyle \frac{1}{\sqrt{2}}}\sqrt{F^2W^2+\gamma ^2\sqrt{4F^2\gamma ^2+(F^2+W^2+\gamma ^2)^2}}.`$ (37)
Note that at resonance $`W=0`$, while out of resonance $`WF,\gamma `$. We will consider three different regimes, $`F\gamma `$, $`F\gamma `$ and $`F\gamma `$. We will concentrate on the first regime and will afterwards stress the problems occurring in the two other regimes.

### 5.1 Expansion for $`𝑭\mathbf{}𝜸`$

In this case we expand Eqs. (37) in $`F`$. We have done so to second order:

$`\mu _1`$ $`=`$ $`2\gamma +{\displaystyle \frac{F^2\gamma }{2\rho ^2}},`$
$`\mu _2`$ $`=`$ $`{\displaystyle \frac{F^2\gamma }{2\rho ^2}},`$
$`\mu _{3,4}`$ $`=`$ $`\gamma iW\left(1+{\displaystyle \frac{F^2}{2\rho ^2}}\right),`$ (38)
where we have used $`\rho ^2=\gamma ^2+W^2`$.

As a next step, we can construct the eigenvectors, and thus the matrix $`𝒳`$, from them. Then we re-express the original wave function vector $`(s,a,h,l)`$ through the new functions $`b_j(\tau )`$ according to (35), again keeping only terms up to order $`F^2`$:

$`s`$ $`=`$ $`\left(1{\displaystyle \frac{F^2(3W^2+\gamma ^2)}{2\rho ^4}}\right)b_1+b_2{\displaystyle \frac{F^2}{4\rho ^2}}`$ (39)
$`{\displaystyle \frac{F}{\rho ^2}}\left[\gamma \left(b_3\mathrm{cos}\mathrm{\Omega }b_4\mathrm{sin}\mathrm{\Omega }\right)W\left(b_3\mathrm{sin}\mathrm{\Omega }+b_4\mathrm{cos}\mathrm{\Omega }\right)\right],`$
$`a`$ $`=`$ $`b_1{\displaystyle \frac{F^2}{4\rho ^2}}+\left(1{\displaystyle \frac{F^2(3W^2+\gamma ^2)}{2\rho ^4}}\right)b_2`$ (40)
$`{\displaystyle \frac{F}{\rho ^2}}\left[\gamma \left(b_3\mathrm{cos}\mathrm{\Omega }b_4\mathrm{sin}\mathrm{\Omega }\right)+W\left(b_3\mathrm{sin}\mathrm{\Omega }+b_4\mathrm{cos}\mathrm{\Omega }\right)\right],`$
$`h`$ $`=`$ $`b_1{\displaystyle \frac{FW}{2\rho ^2}}+b_2{\displaystyle \frac{FW}{2\rho ^2}}+\left(b_3\mathrm{sin}\mathrm{\Omega }+b_4\mathrm{cos}\mathrm{\Omega }\right),`$ (41)
$`l`$ $`=`$ $`b_1{\displaystyle \frac{F\gamma }{2\rho ^2}}b_2{\displaystyle \frac{F\gamma }{2\rho ^2}}+\left(b_3\mathrm{cos}\mathrm{\Omega }b_4\mathrm{sin}\mathrm{\Omega }\right).`$ (42)
Here, we have at the same time introduced the fast-oscillating variable $`\mathrm{\Omega }`$, with $`\mathrm{\Omega }^{}=W`$. Also, while the functions $`b_1`$ and $`b_2`$ are the eigenfunctions of the eigenvalues $`\mu _1`$ and $`\mu _2`$, respectively, the functions $`b_3`$ and $`b_4`$ are superpositions of the eigenfunctions to the remaining eigenvalues: the actual eigenfunctions are $`b_2(\mathrm{cos}\mathrm{\Omega }+i\mathrm{sin}\mathrm{\Omega })+`$ $`b_3(\mathrm{sin}\mathrm{\Omega }+i\mathrm{cos}\mathrm{\Omega })`$ and $`b_2(\mathrm{cos}\mathrm{\Omega }i\mathrm{sin}\mathrm{\Omega })+b_3(\mathrm{sin}\mathrm{\Omega }i\mathrm{cos}\mathrm{\Omega })`$ for $`\mu _3`$ and $`\mu _4`$, respectively. We have chosen these definitions to circumvent imaginary values.

Inserting these definitions into the system of differential equations (34) we derive

$`b_0^{}`$ $`=`$ $`{\displaystyle \frac{F^2\gamma }{2\rho ^2}}b_0,`$ (43)
$`b_2^{}`$ $`=`$ $`\gamma b_2{\displaystyle \frac{1}{2}}b_0\left[\left({\displaystyle \frac{FW}{\rho ^2}}\right)^{}\mathrm{sin}\mathrm{\Omega }+\left({\displaystyle \frac{F\gamma }{\rho ^2}}\right)^{}\mathrm{cos}\mathrm{\Omega }\right],`$ (44)
$`b_3^{}`$ $`=`$ $`\gamma b_3{\displaystyle \frac{1}{2}}b_0\left[\left({\displaystyle \frac{FW}{\rho ^2}}\right)^{}\mathrm{cos}\mathrm{\Omega }\left({\displaystyle \frac{F\gamma }{\rho ^2}}\right)^{}\mathrm{sin}\mathrm{\Omega }\right],`$ (45)
where we have neglected all terms of insignificant order, and $`b_1`$ can be neglected totally. Furthermore, the initial values are given by $`b_0(0)=+1`$ and $`b_{1,2,3}(0)=0`$. Following the corrected Eqs. (D.45–D.51), we get

$$l(\tau )=\frac{F\gamma }{2\rho ^2}b_0\left[1+\frac{F^2}{2\rho ^2}\left(1\frac{2\gamma ^2}{\rho ^2}\right)\right].$$
(46)
Also in this case we find a counter-term which is similar to that derived by .

To obtain the evolution of the lepton asymmetry, we also have to include the anti-neutrinos, i.e.

$$\overline{l}(\tau )=\frac{F\overline{\gamma }}{2\overline{\rho }^2}\overline{b}_0\left[1+\frac{F^2}{2\overline{\rho }^2}\left(1\frac{2\overline{\gamma }^2}{\overline{\rho }^2}\right)\right].$$
(47)
Then, by using the differential equation for $`a`$ and Eq. (13), we get

$$\frac{d\eta }{dq}=\frac{|F|}{2\zeta (3)}_0^{\mathrm{}}𝑑yy\left[f_{\mathrm{eq}}(y,\mu )l\left(\frac{q}{y}\right)f_{\mathrm{eq}}(y,\overline{\mu })\overline{l}\left(\frac{q}{y}\right)\right].$$
(48)
Note that one power of $`y`$ cancels when we change from $`d/d\tau `$ to $`d/dq`$. Furthermore, the $`\gamma a`$ term in the differential equation for $`a`$ disappears due to the conservation of leptonic charge

$$𝑑yy\left[f_{\mathrm{eq}}(y,\mu )\gamma af_{\mathrm{eq}}(y,\overline{\mu })\overline{\gamma }\overline{a}\right]=0.$$
(49)
Equation (48) is valid in the case $`F<\gamma `$. For $`F\gamma `$, the perturbative expansion breaks down. We should also check the validity of Eq. (D.49), i.e. the approximation of the integral. This approximation is valid when $`\gamma 1`$, which is a weaker condition than $`F<\gamma `$ if $`F>1`$, i.e. $`\mathrm{sin}2\theta \sqrt{|\delta m_{\mathrm{eV}}^2|\mathrm{cos}2\theta }>1.8\times 10^5`$. Thus, in the case of small mixing angles and small $`\delta m^2`$, the condition that Eq. (48) is valid is given by $`\gamma >1`$. From solving Eq. (48) numerically, we found that these conditions hold until the power-law regime is established for $`\mathrm{sin}2\theta <5\times 10^5`$ when $`\delta m^2=1\mathrm{eV}^2`$. If we neglect the counter-term in Eq. (48), we get the same solution as , Eqs. (79-80), up to some minor corrections.

To compare our result Eqs. (4648) with , we neglect the chemical potential $`\mu `$, use $`b_0\overline{b}_0`$ and expand $`(l\overline{l})`$ in $`V\eta `$. We get

$`(l\overline{l})`$ $``$ $`{\displaystyle \frac{2F\gamma UV\eta }{\sigma ^4}}b_0(1+{\displaystyle \frac{F^2(U^22\gamma ^2)}{\sigma ^4}}`$
$`+(V\eta )^2\{{\displaystyle \frac{2(U^2\gamma ^2)}{\sigma ^4}}+{\displaystyle \frac{F^2(5U^426U^2\gamma ^2+9\gamma ^4)}{\sigma ^8}}\}),`$
where $`\sigma ^2=U^2+\gamma ^2`$. We compare it with the expansion of (15)

$`2\times 10^{10}L`$ $``$ $`{\displaystyle \frac{2F\gamma UVZ}{\sigma ^4}}b_0(1+{\displaystyle \frac{F^2(U^23\gamma ^2)}{4\sigma ^4}}`$
$`+(VZ)^2\{{\displaystyle \frac{2(U^2\gamma ^2)}{\sigma ^4}}+{\displaystyle \frac{F^2(5U^430U^2\gamma ^2+13\gamma ^4)}{4\sigma ^8}}\}).`$
This equation is given in terms of the variables as defined in . Still we can compare the two equations, since the differences to our definitions cancel. We see that, up to some minor corrections, the main difference is in a factor of $`4`$ in the $`F^2`$ terms. This factor has no fundamental consequences.

### 5.2 Discussion of $`𝑭\mathbf{}𝜸`$

Now we discuss the solutions of Eq. (48) in the region of its validity in more detail. First, let us consider the region where the asymmetry decreases (see Fig. 2). In this region, the corrections to the main terms in Eqs. (46) and (47) are negligible. We expand the integrand of Eq. (48) in terms of $`\eta `$ and $`\mu /T`$ and get

$`f_{\mathrm{eq}}(y,\mu )lf_{\mathrm{eq}}(y,\mu )\overline{l}f_{\mathrm{eq}}(y){\displaystyle \frac{F\gamma _0b_+}{2}}{\displaystyle \frac{4UV\eta }{\rho ^2\overline{\rho }^2}}f_{\mathrm{eq}}(y){\displaystyle \frac{F\gamma _0b_{}}{2}}\left({\displaystyle \frac{1}{\rho ^2}}+{\displaystyle \frac{1}{\overline{\rho }^2}}\right)`$ (52)
$`+`$ $`f_{\mathrm{eq}}(y){\displaystyle \frac{F\gamma _0b_+\eta _{\nu _\mathrm{e}}}{2}}\left({\displaystyle \frac{1}{\rho ^2}}+{\displaystyle \frac{1}{\overline{\rho }^2}}\right)\left\{3[1f_{\mathrm{eq}}(y)]0.2\right\}+O(\eta _{\nu _\mathrm{e}}^2,\eta ^2,\eta _{\nu _\mathrm{e}}\eta ),`$
where we used the notations

$$b_+=\frac{b_0+\overline{b}_0}{2}\text{and}b_{}=b_0\overline{b}_0.$$
(53)
On the RHS of Eq. (52), the first term corresponds to the main term in Eq. (D.50), the second term takes into account the $`AS`$ contribution neglected in Eq. (D.50), and the last term takes into account the non-zero chemical potential both for the $`\gamma `$-terms, see Eq. (5), and for $`f_{\mathrm{eq}}`$, which we rewrite through $`\eta _{\nu _\mathrm{e}}1.5\frac{\mu }{T}`$. Note that we have used $`\overline{\mu }=\mu `$, which is valid for large temperatures. Before substituting Eq. (52) into Eq. (48), let us write the differential equations for $`b_\pm `$, in which we can neglect the small difference between $`\gamma `$ and $`\overline{\gamma }`$:

$`{\displaystyle \frac{db_+}{d\tau }}`$ $`=`$ $`{\displaystyle \frac{F^2\gamma }{2}}\left({\displaystyle \frac{\gamma ^2+U^2+(V\eta )^2}{\rho ^2\overline{\rho }^2}}b_++{\displaystyle \frac{UV\eta }{\rho ^2\overline{\rho }^2}}b_{}\right),`$
$`{\displaystyle \frac{db_{}}{d\tau }}`$ $`=`$ $`{\displaystyle \frac{F^2\gamma }{2}}\left({\displaystyle \frac{\gamma ^2+U^2+(V\eta )^2}{\rho ^2\overline{\rho }^2}}b_{}+{\displaystyle \frac{4UV\eta }{\rho ^2\overline{\rho }^2}}b_+\right).`$ (54)
In the case of small asymmetry $`\eta `$ we approximately get from these equations $`b_+1`$ and $`b_{}\mathrm{sin}^22\theta \eta `$. Using this fact, we can substitute Eq. (52) into Eq. (48) and get the following equation:

$$\frac{d\eta }{dq}=C\eta \left[I_1(q)\mathrm{sin}^22\theta I_2(q)\right]+\eta _{\nu _\mathrm{e}}I_3(q),$$
(55)
where the functions $`I_i(q)`$ are integrals over the momentum which do not depend on the asymmetries. The term $`I_1(q)`$ corresponds to the main term in (D.51), and without $`I_2`$ and $`I_3`$ the asymmetry decreases to infinitely small values.

The second term has opposite sign and is proportional to an additional power of $`\mathrm{sin}^22\theta `$. For $`\mathrm{sin}^22\theta 10^4`$, this term becomes of the order of the first one, but it can not cancel it completely. We would like to mention that, although this term seems to be related to the expansion in $`F`$, it is actually of a different nature, and should therefore not be neglected here. This term corresponds to the $`AS`$ contribution in terms of  and it is proportional to $`b_{}`$ in our notations.

But the actual reason that the lepton asymmetry does not decrease to very small values is the non-zero chemical potential, which corresponds to $`I_3`$. When the asymmetry, which evolves according to Eq. (55), becomes small enough, $`\eta =2\eta _{\nu _\mathrm{e}}\eta _0\eta _{\nu _\mathrm{e}}\eta _0/2`$, the last term in Eq. (55) starts to compensate the first two terms, and prevents the asymmetry to get tiny values.

We will not consider the regime of exponential growth of the asymmetry (see Fig. 2) separately, let us just note again that during this stage the resonance condition for anti-particles (particles) go to very large momenta, while the resonance condition for particles (anti-particles) go to small momenta, $`y<1`$, for $`\eta >0`$ ($`\eta <0`$), see Fig. 2. Therefore, in the power-law regime, one may neglect either neutrinos or anti-neutrinos in (48), since the contribution of one of them is negligible due to the exponential damping of the Fermi distribution function.

Finally, we discuss the power-law regime. In this regime, we can neglect either neutrinos or anti-neutrinos. For definiteness, we consider the case $`\eta <0`$, where the resonance momentum of the anti-neutrinos is very large. Thus we can neglect $`\overline{l}`$ in Eq. (48). Also, we will assume that the asymmetry is still $`\eta 1`$, so that we can also neglect the chemical potential. This is in contrast to the case described by Eq. (55), where the main terms of the $`\nu `$ and $`\overline{\nu }`$ almost canceled each other, so that the even the small chemical potential would become important. Here, however, the $`\overline{\nu }`$ have no influence at all, so that no cancellation of the main terms appears. Furthermore, we assume $`b_01`$, which is valid for small $`F/\gamma `$. To be able to integrate analytically, we expand $`W`$ around $`y_{\mathrm{res}}`$, given in Eq. (14):

$$W(y)c_W(yy_{\mathrm{res}})+\mathrm{},\text{with}c_W=\frac{2y_{\mathrm{res}}Q}{q^2}+V\eta .$$
(56)
Then the leading order in (48) gives

$`{\displaystyle \frac{d\eta }{dq}}`$ $`=`$ $`{\displaystyle \frac{F^2\delta }{4\zeta (3)q^2}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑yf_{\mathrm{eq}}(y){\displaystyle \frac{y^3}{\gamma ^2+W^2}}`$ (57)
$``$ $`{\displaystyle \frac{F^2\delta }{4\zeta (3)q^2}}f_{\mathrm{eq}}(y_{\mathrm{res}}){\displaystyle \frac{y_{\mathrm{res}}^3}{\gamma ^2(y_{\mathrm{res}})}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑y{\displaystyle \frac{1}{1+\frac{c_W^2}{\gamma _{\mathrm{res}}^2}(yy_{\mathrm{res}})^2}}`$
$`=`$ $`{\displaystyle \frac{F^2\delta }{4\zeta (3)q^2}}f_{\mathrm{eq}}(y_{\mathrm{res}}){\displaystyle \frac{y_{\mathrm{res}}^3}{\gamma ^2(y_{\mathrm{res}})}}{\displaystyle \frac{\gamma (y_{\mathrm{res}})}{c_W}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑z{\displaystyle \frac{1}{1+z^2}},`$
where we have substituted $`z=c_W(yy_{\mathrm{res}})/\gamma _{\mathrm{res}}`$. The final integral gives $`\pi `$. Analogously, we find the next order corrections, so that we get

$$\frac{d\eta }{dq}=\frac{F^2\pi }{4\zeta (3)}f_{\mathrm{eq}}(y_{\mathrm{res}})\frac{y_{\mathrm{res}}}{c_W}\left(1\frac{F^2}{8\gamma ^2(y_{\mathrm{res}})}\right),$$
(58)
which is our differential equation for the asymmetry $`\eta `$ in the power-law regime. We can compare this to Eq. (17), which is the corresponding equation for the approach by Dolgov et al.. The interesting point is that the counter-term, compared with the leading term, is the same in both approaches up to a factor of 2 (the counter-term in Eq. (17) corresponds to $`F^2/(16\gamma ^2(y_{\mathrm{res}}))`$ in our notation).

¿From Eq. (58) it is easy to derive the power law for $`\eta `$. Let us first neglect the counter-term. If we further assume that $`C(q)q`$ in (14), we get $`y_{\mathrm{res}}Q/(V_0\eta q^{4/3})`$ and can neglect the first term in $`c_W`$. Then

$$\frac{d\eta }{dq}\frac{F^2\pi }{4\zeta (3)}f_{\mathrm{eq}}(y_{\mathrm{res}})\frac{Q}{V_0^2}\eta ^2q^{8/3}C_\eta \eta ^2q^{8/3}.$$
(59)
If we multiply by $`dq`$ and integrate, we get

$$\eta =\left(\frac{9}{11}C_\eta \right)^{1/3}q^{11/9}.$$
(60)
Here, we have neglected the boundary conditions, which become negligible for large $`q`$ in any case. We see that this corresponds to the $`\eta T^{11/3}`$ power law found in  when neglecting the counter-term. We should however mention that even for all the simplifications made above, this power law is only approximate, since $`C_\eta `$ contains the Fermi distribution $`f_{\mathrm{eq}}(y_{\mathrm{res}})`$. Since $`y_{\mathrm{res}}q^{1/12}`$ for the derived power law, $`f_{\mathrm{eq}}(y_{\mathrm{res}})`$ will slowly decrease, thus changing Eq. (60).

If we take the counter-term into account, the condition for the power law is that the brackets in Eq. (58) equal 0. Using again the approximation $`y_{\mathrm{res}}Q/(V_0\eta q^{4/3})`$, we can derive

$$\eta =\frac{Q}{V_0}\left(\frac{8\delta ^2}{F^2}\right)^{1/4}q^{1/4},$$
(61)
which is exactly the $`\eta T^1`$ power law found in .

Let us still have a look on the evolution of $`l`$ for the $`\eta T^1`$ power law. In Fig. 3 we have plotted $`l`$ as given in Eq. (46), once with and once without the counter-term, for the case where $`F=\sqrt{2}\gamma `$ at resonance. Then the counter-term merely compensates the leading order term near the resonance, but the integration over the momentum still yields a similarly large value when taking the counter-term into account. Thus, it is safe to say that the counter-term is of minor importance for $`F\sqrt{2}\gamma `$. This corresponds to the result in Eq. (58).

To summarize, this ansatz gives similar results as the approach of , so that the divergence in their approach seems not to be critical. However, our approach has shown two further advantages: First, we need not make the adiabatic approach which they used for $`L`$ and $`H`$. Second, our result yields the correct behavior for very small $`\eta `$ before the exponential regime, since we included the term $`(F/2)(as)`$ and the chemical potential.

### 5.3 $`𝑭\mathbf{}𝜸`$

This region can occur in two cases: either when $`F`$, and thus the mixing angle $`\mathrm{sin}2\theta `$, is large, or when $`\gamma `$ is very small, which is the case for $`\tau 1`$. Since we only discuss small mixing angles in our analysis, the first case is not important. However, the second case appears, since for large asymmetries, the (anti-)neutrinos at resonance have $`\tau =q/y`$, where $`q>1`$ is large and $`y0.1`$ is small.

In this limit the eigenvalues Eqs. (37) should be expanded in $`\gamma `$. However, in this regime, the approximation in Eq. (D.49) does not hold, since $`\gamma `$ is no longer very large. Therefore, we did not investigate further in this case.

### 5.4  $`𝑭\mathbf{}𝜸`$

In this case let us rewrite $`\gamma =F(1+\chi )`$ and suppose that $`\chi 1`$. Then we can expand the eigenvalues (37) in $`\chi `$. To the main order, which does not contain $`\chi `$, we get:

$$\mu _i=F\pm \frac{1}{\sqrt{2}}\sqrt{W^2\pm W\sqrt{W^2+4F^2}}+O(\chi ).$$
(62)
The important feature of Eq. (62) is the fact that at resonance, $`W=0`$, all 4 eigenvalues coincide. This means that the resonance matrix $``$ changes its rank, and thus diagonalization becomes impossible.

We conclude that the approach of Dolgov et al. is useful and correct as long as $`F<\gamma `$, but that it seems very difficult to develop this kind of approach for $`F\mathrm{\Gamma }>\gamma `$.

## 6 Solution for $`\mathbf{|}𝒂\mathbf{}𝒔\mathbf{}\mathrm{𝟏}\mathbf{|}\mathbf{}\mathrm{𝟏}`$

### 6.1 Solution for one momentum mode

In this section, we present a different analytic solution of the system of differential equations (8) in the power-law regime for a given momentum $`y`$. We introduce two simplifications. First, to circumvent the non-linearity, we give the lepton asymmetry an explicit form,

$$\eta (T)=\eta _iT^\alpha ,$$
(63)
where $`\alpha `$ represents the power-law behavior. Second, we treat the change of the variables $`a`$ and $`s`$, expressed by the function $`f_{(as)}(as1)`$, as a perturbation in the differential equation of $`l`$. In terms of the physical content, $`f_{(as)}1`$ means that there is only very little conversion from active to sterile neutrinos. The advantage of this approximation is simple: provided that the $`T^1`$ power law holds, the condition $`f_{(as)}1`$ holds for all times222This condition corresponds to $`P_z11`$ in the formalism used by  and others..

To start with, we write down the differential equations for $`h`$ and $`l`$:

$`l^{}`$ $`=`$ $`{\displaystyle \frac{F}{2}}(1+f_{(as)})Wh\gamma l,`$ (64)
$`h^{}`$ $`=`$ $`Wl\gamma h.`$ (65)
Subsequently, we will assume that $`f_{(as)}`$ is an explicit function. Then the above equations form a complete set of differential equations for $`h`$ and $`l`$.

By solving (64) for $`h(\tau )`$, taking its derivative $`h^{}(\tau )`$, and inserting both into (65), we find the differential equation of second order for $`l(\tau )`$:

$`Wl^{\prime \prime }+\left(2W\gamma W^{}\right)l^{}+\left(W^3+\gamma ^2W+W\gamma ^{}\gamma W^{}\right)l`$
$`+{\displaystyle \frac{F}{2}}\left[(1+f_{(as)})\left(W^{}\gamma W\right)Wf_{(as)}^{}\right]`$ $`=`$ $`0.`$ (66)
Analogously, we find

$`Wh^{\prime \prime }+\left(2W\gamma W^{}\right)h^{}+\left(W^3+\gamma ^2W+W\gamma ^{}\gamma W^{}\right)h`$
$`{\displaystyle \frac{F}{2}}(1+f_{(as)})W^2`$ $`=`$ $`0,`$ (67)
where the only difference to (66) is given in the inhomogeneous $`F`$-term.

We can now solve the equation for $`l`$ analytically. To this end, we first write down the solution for the homogeneous differential equation, i.e. $`F=0`$:

$`l_{\mathrm{hom}}(\tau )`$ $`=`$ $`\left(k_c^l\mathrm{cos}({\displaystyle \underset{\tau _i}{\overset{\tau }{}}}W_2𝑑\tau _2)+k_s^l\mathrm{sin}({\displaystyle \underset{\tau _i}{\overset{\tau }{}}}W_2𝑑\tau _2)\right)`$ (68)
$`\times \mathrm{exp}\left({\displaystyle \underset{\tau _i}{\overset{\tau }{}}}\gamma _2𝑑\tau _2\right).`$
Here, the coefficients $`k_c^l`$ and $`k_s^l`$ are given by the boundary conditions at time $`\tau _i`$, and the subscript 2 means that the parameters are evaluated at $`\tau _2`$. The solution for the inhomogeneous differential equation can then easily be found, e.g. by using ,

$`l_{\mathrm{inhom}}(\tau )`$ $`=`$ $`{\displaystyle \frac{F}{2}}{\displaystyle \underset{\tau _i}{\overset{\tau }{}}}{\displaystyle \frac{(1+f_{(as)})_1\left(W_1^{}\gamma _1W_1\right)W_1(f_{(as)}^{})_1}{W_1^2}}`$ (69)
$`\times \mathrm{sin}\left({\displaystyle \underset{\tau _1}{\overset{\tau }{}}}W_2𝑑\tau _2\right)\mathrm{exp}\left({\displaystyle \underset{\tau _1}{\overset{\tau }{}}}\gamma _2𝑑\tau _2\right)d\tau _1.`$
The complete solution is then $`l_{\mathrm{analytic}}=l_{\mathrm{hom}}+l_{\mathrm{inhom}}`$.

Analogously, we get $`h=h_{\mathrm{hom}}+h_{\mathrm{inhom}}`$, where $`h_{\mathrm{hom}}`$ is identical to the homogeneous solution for $`l`$, (68), now with coefficients $`k_c^h`$ and $`k_s^h`$. The inhomogeneous solution is

$$h_{\mathrm{inhom}}=\frac{F}{2}\underset{\tau _i}{\overset{\tau }{}}(1+f_{(as)})_1\mathrm{sin}\left(\underset{\tau _1}{\overset{\tau }{}}W_2𝑑\tau _2\right)\mathrm{exp}\left(\underset{\tau _1}{\overset{\tau }{}}\gamma _2𝑑\tau _2\right)𝑑\tau _1.$$
(70)
Now we turn to the variables $`a`$ and $`s`$. It is straightforward to solve the differential equations describing them, see (8), provided that $`l`$ is an explicit function:

$`a(\tau )`$ $`=`$ $`a_i\mathrm{exp}\left({\displaystyle \underset{\tau _i}{\overset{\tau }{}}}\gamma _2𝑑\tau _2\right)F{\displaystyle \underset{\tau _i}{\overset{\tau }{}}}l_{\mathrm{expl}}(\tau _1)\mathrm{exp}\left({\displaystyle \underset{\tau _1}{\overset{\tau }{}}}\gamma _2𝑑\tau _2\right)𝑑\tau _1`$ (71)
$`s(\tau )`$ $`=`$ $`s_i+F{\displaystyle \underset{\tau _i}{\overset{\tau }{}}}l_{\mathrm{expl}}(\tau _1)𝑑\tau _1.`$ (72)
We still need to settle the boundary conditions: From  we know that for $`\tau _i=0`$, $`h=l=a=0`$ and $`s=1`$. It follows that $`k_c^l=k_s^l=k_c^h=k_s^h=a_i=0`$ and $`s_i=1`$ as initial conditions.

The next step is to connect the two systems of differential equations, $`(h,l)`$ and $`(a,s)`$, using perturbation theory. The leading order, labeled with the subscript $`0`$, is defined by $`f_{(as)}l_{\mathrm{expl}}0`$. We find $`a_0(\tau )s_0(\tau )0`$, while $`l_0`$ and $`h_0`$ are given by (69) and (70), respectively, using $`f_{(as)}f_{(as)}^{}0`$.

The first-order corrections, $`a_1`$ and $`s_1`$, are found by inserting $`l_{\mathrm{expl}}=l_0`$ into the analytic solutions, (71) and (72), respectively, and for $`h`$ and $`l`$ by inserting $`f_{(as)}=a_1s_11`$ into (70) and (69), respectively. For higher order corrections, this procedure is repeated recursively. It is clear that this expansion holds only provided that $`f_{(as)}1`$ for all earlier times.

### 6.2 How to check the consistency of a power law

In the previous subsection, we derived an analytic solution for the density matrix of a single momentum $`y`$. The analytical equation for $`\eta `$ in the power-law region is found by using the differential equation for $`a`$ and Eq. (13),

$$\left(\frac{d\eta }{dq}\right)_{\mathrm{analytic}}=\underset{0}{\overset{\mathrm{}}{}}\frac{dy}{2\zeta (3)}f_{\mathrm{eq}}(y)yFl[q,W(q,y,\eta )],$$
(73)
where we have neglected the anti-neutrinos and have used the conservation of leptonic charge, Eq. (49). We can insert the solution of $`l`$ as given in the previous subsection by using the relation $`\tau =q/y`$. However, note that the resulting equation for $`\eta `$ is valid only as long as $`f_{(as)}1`$.

Furthermore, we need to define $`W`$, which is the only parameter in the analytic solution of $`l`$ that is dependent on $`\eta `$ and $`y`$. If we define $`\eta (q)=\eta _0q^{\alpha /3}`$, where $`\alpha `$ is the temperature behavior of $`\eta `$, $`\eta T^\alpha `$, we can write

$$W(\tau ,y)=U(\tau )V_0\eta _0\tau ^{(4+\alpha )/3}y^{(1+\alpha )/3}.$$
(74)
Having given $`W`$ in this form, we now have two unknown parameters: the power-law behavior $`\alpha `$ and the asymmetry factor $`\eta _0`$.

Checking the consistency of a power law is now simple in principle: given the correct values for $`\alpha `$ and $`\eta _0`$, the integration of the rhs of (73) for any $`q`$ in the power-law regime should yield

$$\left(\frac{d\eta }{dq}\right)_{\mathrm{analytic}}=\frac{\alpha }{3}\eta _0q^{1\alpha /3}.$$
(75)
### 6.3 Results

We now consider the case derived by , i.e. $`\alpha =1`$. In this case, (74) does not depend on $`y`$ directly, while the $`y`$-dependence persists indirectly through $`\tau =q/y`$. The remarkable consequence is that $`\tau _{\mathrm{res}}`$ is independent on $`y`$, and thus is the ratio of $`F`$ to $`\gamma `$ at resonance, $`RF/\gamma _{\mathrm{res}}`$.

We found it convenient to substitute the integration over $`y`$ in Eq. (73) with an integration over $`\tau `$:

$$\frac{d\eta }{dq}=\frac{Fq^2}{2\zeta (3)}\underset{0}{\overset{\mathrm{}}{}}\frac{d\tau }{\tau ^3}f_{\mathrm{eq}}(q/\tau )l(\tau ).$$
(76)
Also, we use the parameter $`R`$ instead of $`\eta _0`$. The advantage of this is obvious: the $`T^1`$ power law in  originates in the counter-term compensating the leading order term, which occurred when $`R=\mathrm{const}=𝒪(1)`$. This yields $`\tau _{\mathrm{res}}=\sqrt{F/(R\delta )}`$. Furthermore, by solving $`W(\tau _{\mathrm{res}})=0`$, we get

$$|\eta _0(R)|=QV_0^1\tau _{\mathrm{res}}(R)\left[\tau _{\mathrm{res}}^2(R)1\right].$$
(77)
The following calculations were done using the fits presented in Appendix. To begin, we have calculated the evolution of $`l_{\mathrm{analytic}}`$ using $`f_{(as)}=0`$ for $`0.1<R<1000`$. We have given an example for $`R=100`$ in Fig. 5. Within this range, $`l`$ showed a similar behavior, the only differences occurring in the scaling of $`\tau `$ and the amplitude of $`l`$. Then we calculated the first-order contribution from $`f_{(as)}`$. In terms of the formalism used in , this corresponds to including contributions from $`b_0^{}`$, i.e. the origin of the counter-term they found. We found that the correction was negligible for $`R<10`$, and that it became of order of the main term only when $`R1000`$, and also then stayed smaller than the LO term.

Thus, we conclude that the counter-term found in  must be an artifact of the expansion in $`F`$, which is only valid for $`R1`$.

Looking at Fig. 5, one might argue that even without the NLO term, the oscillations after the resonance might be the origin of a large counter-term. However, we found this not to be the case. We calculated Eq. (76) for $`0.2<R<10`$. According to (75), the result should be equal to $`\eta _0q^{2/3}/3`$ for all $`q`$, with $`\eta _0`$ as given in (77). We found that this was not the case, even when including the first-order contribution from $`f_{(as)}`$: as an example, we show $`\eta _0q^{2/3}/3`$ and (76) in Fig. 6, with $`R=1`$.

By the way, we found that $`\left(\frac{d\eta }{dq}\right)_{\mathrm{analytic}}`$ could be fitted very well by the simple function

$$\left(\frac{d\eta }{dq}\right)_{\mathrm{analytic}}=5.4\times 10^6\frac{q^2}{R}f_{\mathrm{eq}}\left(\frac{q}{\tau _{\mathrm{res}}}\right),$$
(78)
in the region $`0.2<R<10`$, $`0.1<q<60`$. This corresponds to our understanding that only the momenta close to resonance are significant. From this fit, we see that $`\left(\frac{d\eta }{dq}\right)_{\mathrm{analytic}}`$ has a totally different $`q`$ behavior than the expected $`q^{2/3}`$.

The reason that the oscillations do not bring a counter-term can be seen in the evolution of $`s`$, shown in Fig. 5. The better part of the asymmetry is created in the region around the resonance, where $`l`$ varies slowly. The oscillations after this region very soon become fast, and each half period approximately compensates the previous one.

For the case $`\alpha 1`$, the situation is more complicated, since then the evolution of $`l`$ depends on $`y`$. We intend to eventually present the results for this case in a later publication.

However, we decided to consider here the special case $`\alpha =11/3`$, i.e. the power law derived in Eq. (60), numerically. To this end, we assumed this power law, with $`\eta _0=(9/11C_\eta )^{1/3}`$, where $`C_\eta `$ was defined in Eq. (59). Then we calculated the evolution of $`l`$ for a grid of momenta $`y`$ numerically, using the exact differential equations (8), and integrated the results for several $`q`$ according to Eq. (73). Then we compared the result with the expected value, in this case given by

$$\left(\frac{d\eta }{dq}\right)_{\mathrm{expected}}=\frac{11}{9}\eta _0q^{2/9}.$$
(79)
We have plotted the ratio in Fig. 7. If our result was perfect, the ratio would have been $`1`$ for all $`q`$. The discrepancy has several origins: the result depends on $`\alpha `$ and $`\eta _0`$. As an example, choosing a different $`\eta _0`$ will shift the curve upward or downward. Hereby, already a change of a few percent would suffice to make the curve in Fig. 7 cross unity. Furthermore, a change in $`\alpha `$ alters the curvature of the plot. In particular, the downward slope at small $`q`$ is due to the fact that we have used the approximation $`C(q)q`$ when deriving the power law, see Eqs. (5860). In this context, the result presented in Fig. 7 is fairly good. We conclude that it is safe to say that the real power law is around $`\alpha 4`$.

## 7 Conclusions

In this work we have analytically investigated the system of kinetic equations which governs the creation of lepton asymmetry in the Early Universe due to active-sterile neutrino oscillations with small vacuum mixing angle $`\mathrm{sin}2\theta <10^2`$ and negative $`\delta m^2`$. We have improved the analytical approach of , taking into account the non-zero chemical potential of neutrinos and relaxing the adiabatic approximations made in .

In analogy to , we have derived an equation for the asymmetry evolution \[Eq. (48)\], which we found to be valid only if the collision terms dominate the kinetic equations (8), or in our notation if $`F<\gamma `$. Except for very small values of the asymmetry, our approach gave the same results as found in . Also, the main term in Eq. (48) coincides with the so-called static limit derived in . However, Eq. (48) contains an additional term which compensates the main term when $`F\gamma `$.

It was found by the authors of  that this additional term resulted in a growth of the asymmetry according to a $`T^1`$ power law. However, we found that the additional term in Eq. (48) becomes important only when the expansion breaks down, or $`F/\gamma =𝒪(1)`$. But strictly speaking, this fact does not prove that $`\eta T^1`$ is wrong, because the additional term already exists when the perturbative expansion is valid.

In order to check the validity of the $`T^1`$ power law, we developed a new analytical approach which is independent from perturbations in $`F/\gamma `$. We expand the kinetic equations not in the coefficients (as done in ), but by expanding in the combination of functions $`(as1)`$, which remains small as long as the number density of the sterile neutrinos remains small for all momenta.

This condition strongly depends on the power-law index. For the $`T^1`$ power law, the number density of the sterile neutrinos remains small for all momenta. This fact allowed us to check the consistency of $`\eta T^1`$, and we found that this power law is not a solution of the kinetic equations. Thus, we proved analytically that $`\eta T^1`$ is erroneous.

For the $`T^4`$ power law, which corresponds to a full transition of active to sterile neutrinos, the condition $`(as)1`$ breaks down already in the region $`F/\gamma <1`$. For the $`T^{11/3}`$ power law, this condition holds until $`F/\gamma 1`$. Thus, our new method is not useful beyond this point.

We conclude that our Eq. (48) provides a correct analytical solution, which is in qualitative agreement with the numerical results of  in the region prior to the formation of a power law. For $`\mathrm{sin}2\theta <5\times 10^4`$, our solution shows the formation of the power law $`T^{11/3}`$. On the other hand, we can not apply our Eq. (48) in the region where the MSW effect dominates and can not decide if this power law changes to $`T^4`$ or not. Either way, the final value of asymmetry should be large.

When this work was nearly complete we received a manucript from P.Di Bari, R. Foot, R.R. Volkas and Y.Y.Y. Wong , who also conclude that the power law $`\eta T^1`$ is erroneous by showing that the approach of  does not work in the region where the MSW effect dominates over the second-order effects.

## Acknowledgment

We are grateful to S. Hansen and G. Raffelt for helpful comments and reading the manuscript, and would like to thank M. Gorbahn, F. Gubarev, A. Neronov and B. Pötter for fruitful discussions.

We also like to thank P. Di Bari, R. Foot, R.R. Volkas and Y.Y.Y. Wong for making their manuscript  available before publication. We thank R. Foot for comments on our manuscript.

This work was partly supported by the Deutsche Forschungsgemeinschaft under grant No. SFB 375 and in part by INTAS grant 1A-1065.

## Appendix: Fits for l

Since the exact analytic function for $`l`$ is an integral, and we are forced to integrate over $`l`$, we have fitted $`l`$ to speed up the numerics. Here we present the fits used in our calculations in subsection 6.3.

At very small $`\tau `$, the adiabatic solution, given by $`h^{}=l^{}=0`$, is a good approximation. However, for small $`\gamma `$, the condition $`h_{\mathrm{adiab}}^{}Wl_{\mathrm{adiab}},\gamma h_{\mathrm{adiab}}`$ breaks down long before the resonance. We therefore relax the condition $`h^{}=0`$, while we maintain $`l^{}=0`$ (which is valid much longer due to the term $`F/2`$). Thus, by setting $`l^{\prime \prime }=l^{}=0`$ in (66), we find

$$l_{\mathrm{semiad}}=\frac{F}{2}\frac{\gamma WW^{}}{W^3+W\gamma ^2+W\gamma ^{}\gamma W^{}}.$$
(80)
The error of this solution is of order $`l_{\mathrm{semiad}}^{}/(F/2)`$.

For the region before and during the resonance, we decided to calculate the analytic solution of $`l`$ for a sample of points and fit a function to it. We found that before the resonance, $`l`$ showed an exponential-like behavior, which is why in this region we used a fit of the form

$$l_{\mathrm{expfit}}=e^{\alpha (\tau )},\text{where}\alpha (\tau )=\alpha _0+\alpha _1\tau +\mathrm{}$$
(81)
Using $`n=30`$ points and a polynomial of order $`(n1)`$, we found it necessary to do two such fits to achieve an accuracy of better than $`1\%`$.

Around the resonance and shortly after, we applied two linear fits, also with $`n=30`$.

Finally, we fitted the oscillating behavior after the resonance. We found that it is a good approximation to add the homogeneous solution (68) to the semi-adiabatic fit (80),

$$l_{\mathrm{oscfit}}(\tau )=l_{\mathrm{hom}}(\tau )+l_{\mathrm{semiad}}(\tau ).$$
(82)
Then the semi-adiabatic fit accounts for the inhomogeneous term

When using the analytic solution for $`\tau >\tau _0`$, we found that the integrand diverges at the resonance, $`\tau _1=\tau _0`$. However, this problem can easily be circumvented: between $`\tau _0`$ and some point $`\tau _ϵ\mathrm{\Gamma }>\tau _0`$ we used a quadratic fit for $`l`$ with boundary conditions given by the analytic solution at $`\tau _i=\tau _0`$, while for $`\tau >\tau _ϵ`$ we used the analytic solution with boundary conditions given by the quadratic fit at $`\tau _i=\tau _ϵ`$. We took $`\tau _ϵ=\tau [W(\tau )=F]`$.

We should also mention that for $`\tau \tau _0`$ the integrand of the analytic solution is oscillating fast, so that it is difficult to calculate the integral starting from $`\tau =0`$. Therefore, we have used the semi-adiabatic solution up to some $`\tau _i\tau _0`$, where the semi-adiabatic solution still is a very good fit, e.g. has an error of less than $`10^3`$.


format: text

text:
# Abstract

## Abstract

113 residual product nuclide yields in a 1.0 GeV proton-irradiated thin monoisotopic 208Pb sample and 107 residual product nuclide yields in a 2.6 GeV proton-irradiated natW sample have been measured and simulated by 8 different codes. The irradiations were made using proton beams extracted from the ITEP synchrotron. The nuclide yields were $`\gamma `$-spectrometered directly using a high-resolution Ge-detector. The $`\gamma `$-spectra were processed by the GENIE-2000 code. The ITEP-developed SIGMA code was used together with the PCNUDAT nuclear decay database to identify the $`\gamma `$-lines and to determine the cross-sections. The 27Al(p,x)22Na reaction was used to monitor the proton flux. The measured yields are compared with calculation the LAHET (with ISABEL and Bertini options), CEM95, CEM2k, CASCADE, CASCADE/INPE, INUCL, and YIELDX codes. Estimates of the mean deviation factor are used to demonstrate the predictive power of the codes. The results obtained may be of interest in studying the parameters of the Pb and W target modules of the hybrid Accelerator-Driven System (ADS) facilities.

## Introduction

At present, the Pb-Bi eutectic and W are regarded as the most promising target materials for ADS facilities -. As a result, the high-energy irradiation mode of using the materials necessitates additional studies of the nuclear-physics characteristics of Pb, Bi, and W, particularly the yields of residual product nuclei under proton irradiation in a broad range of energies from a few MeV to 2-3 GeV. Results of such studies are extremely important when designing even demonstration versions of the ADS facilities.

Undoubtedly, computational methods will play an important role when forming a set of nuclear constants for ADS facilities. Therefore, verification of the most extensively used simulation codes has proved to be of a high priority.

## Basic definitions and computational relations

The formalism of representing the reaction product yields (cross sections) in high-energy proton-irradiated thin targets is described in sufficient detail in . In terms of the formalism, the variations in the concentration of any two chain nuclides produced in an irradiated target ($`N_1\stackrel{\lambda _1}{}N_2\stackrel{\lambda _2}{}`$) may be presented to be a set of differential equations that describe the production and decays of the nuclides. By introducing a formal representation of the time functions, $`F_i`$, of the form $`F_i=\left(1e^{\lambda _i\tau }\right)\frac{1e^{\lambda _iKT}}{1e^{\lambda _iT}}\text{ ,}`$ (i=1, 2, and Na or another monitor product; $`\tau `$ is the duration of a single proton pulse; $`T`$ is the pulse repetition period; $`K`$ is the number of pulses within the irradiation period), which characterize the nuclide decays within the irradiation time, and by expressing (similar to the relative measurements) the proton fluence size via the monitor reaction cross section, $`\sigma _{st}`$, we can present the cumulative and independent yields as:

$$\sigma _1^{cum}=\frac{A_0}{\eta _1\epsilon _1F_1N_{Na}}\frac{N_{Al}}{N_T}\frac{F_{Na}}{\lambda _{Na}}\sigma _{st}$$
(1)
$$\sigma _1^{cum}=\frac{A_1}{\nu _1\eta _2\epsilon _2F_1N_{Na}}\frac{N_{Al}}{N_T}\frac{\lambda _2\lambda _1}{\lambda _2}\frac{F_{Na}}{\lambda _{Na}}\sigma _{st}$$
(2)
$$\sigma _2^{ind}=\left(\frac{A_2}{F_2}+\frac{A_1}{F_1}\frac{\lambda _1}{\lambda _2}\right)\frac{1}{\eta _2\epsilon _2N_{Na}}\frac{N_{Al}}{N_T}\frac{F_{Na}}{\lambda _{Na}}\sigma _{st}$$
(3)
$$\sigma _2^{cum}=\sigma _2^{ind}+\nu _1\sigma _1^{cum}=\left(\frac{A_1}{F_1}+\frac{A_2}{F_2}\right)\frac{1}{\eta _2\epsilon _2N_{Na}}\frac{N_{Al}}{N_T}\frac{F_{Na}}{\lambda _{Na}}\sigma _{st}\text{ ,}$$
(4)
where $`\sigma _1^{cum}`$ is the cumulative cross section of the first nuclide; $`\sigma _2^{ind}`$ and $`\sigma _2^{cum}`$ are the independent and cumulative cross sections of the second nuclide; NAl and NT are the numbers of nuclei in the monitor (standard) and in experimental sample, respectively; $`\eta _1`$ and $`\eta _2`$ are the $`\gamma `$-line yields; $`\epsilon _1`$ and $`\epsilon _2`$ are the spectrometer efficiencies at energies E$`_{\gamma _1}`$ and E$`_{\gamma _2}`$; $`\nu _1`$ is the branching ratio of the first nuclide; $`\lambda _1`$, $`\lambda _2`$, and $`\lambda _{Na}`$ are, respectively, the decay constants of the first and second nuclides and of the monitor product (22Na and/or 24Na).

The factors A0, A1, and A2 are calculated through fitting the measured counting rates in the total absorption peaks, which correspond to energies E$`_{\gamma _1}`$ (the first nuclide) and E$`_{\gamma _2}`$ (the second nuclide), by exponential functions. It should be noted that formulas (1)–(4) were derived on assumption that the $`\gamma `$-intensities of the two nuclides produced under irradiation are recorded up to the desired accuracy within an interval of time from the irradiation end to the moment of the ultimate detectable intensity. If, for some reasons, the factor A1 cannot be found, then the factor A2 will be used together with expression (14) from  to determine the quantity $`\sigma _2^{cum^{}}`$, which we called the supra cumulative yield:

$`\sigma _2^{cum^{}}`$ $`=`$ $`\sigma _2+{\displaystyle \frac{\lambda _2}{\lambda _1\lambda _2}}\nu _1\sigma _1^{cum}=`$ (5)
$`=`$ $`{\displaystyle \frac{A_2}{\eta _2\epsilon _2F_2N_{Na}}}{\displaystyle \frac{N_{Al}}{N_T}}{\displaystyle \frac{F_{Na}}{\lambda _{Na}}}\sigma _{st}`$
The resultant value of $`\sigma _2^{cum^{}}`$ may prove to be very different from $`\sigma _2^{cum}`$. Nevertheless, the supra cumulative yield can be either used directly to verify the codes, or determined further up to $`\sigma _2^{cum}`$ if the needed data are obtained elsewhere (for example, from inverse-kinematics experiments).

## Experimental techniques

A 10.5 mm diameter, 139.4 mg/cm2 monoisotopic 208Pb metal foil sample (97.2% 208Pb, 1.93% 207Pb, 0.87% 206Pb, $`<`$0.01% 204Pb, $`<`$ 0.00105% of chemical impurities) and a 38.1 mg/cm2 natW metal foil sample (99.95% W, $`<`$0.05% of chemical impurities), both of 10.5-cm diameter, were proton-irradiated. 139.6 mg/cm2 and 139.1 mg/cm2 Al foils of the same diameter were used as monitors. Chemical impurities of the monitor did not exceed 0.001%.

The samples were irradiated by the external proton beam from the ITEP U-10 synchrotron . The average flux densities during irradiation of Pb and W samples were of 1.4$`\times `$1010 p/cm2 and 2.8$`\times `$1010 p/cm2, respectively.

Our measurements were supported by extra researches aimed at reducing the systematic errors in the experimental results. Theses researches included:

* experiments to specify the neutron component in the extracted proton beams,
* experiments to specify the 27Al(p,x)24Na monitor reaction cross section,
* studies to specify the dependence of the $`\gamma `$-spectrometer detection efficiency on the position geometry of irradiated sample,
* studies to optimize the $`\gamma `$-spectrum simulation codes.

Figs. 2 and 2 show the results of measuring the neutron component in the extracted proton beams, i.e., the neutron-to-proton flux density ratio, $`\mathrm{\Phi }_n/\mathrm{\Phi }_p`$. Fig. 4 presents the monitor reaction cross sections measured here and in other works111 MI85 – Nucl. Phys. A 441 (1985) 617; MI86 – NIM B 16 (1986) 61; MI89 – Analyst 114 (1989) 287; Mi90 – NEANDC(E)-312-U(1990) 46; MI93 – INDC(GER)-037/LN(1993) 49; MI95 – NIM B 103 (1995) 183; MI96 – NIM B 114 (1996) 91; MI97 – NIM B 129 (1997) 153.. The height-energy dependence of the detection efficiency is displayed in Fig. 4.

The discrepancies between the two sets of high-energy data in Fig. 4 have yet to be studied. The analytical expression of the spectrometer detection efficiency as a function of energy and sample position height is

$`\epsilon (E,H)=\epsilon _{base}\left(E\right)\left[{\displaystyle \frac{\left(q_1+q_2\mathrm{ln}E+H_{base}\right)}{\left(q_1+q_2\mathrm{ln}E+H\right)}}\right]^2,`$ (6)
where $`q_1`$ and $`q_2`$ are parameters defined by fitting the experimental results. An analysis of the $`\gamma `$-spectrum processing codes has shown that the GENIE-2000 code is superior to the others because of its interactive mode of fitting the peaks, which permits correction of the automated computer-aided processing; therefore, we chose it for our work.

## Experimental results and measurement errors

Tables 1 and 2 present the results of measuring the reaction product yields in the 1 GeV proton-irradiated 208Pb and 2.6 GeV proton-irradiated natW samples. 113 yields from 208Pb have been obtained, of which, 6 independent yields (i), 17 independent yields of metastable states (m), 15 independent yields of metastable and ground states ($`\mathrm{\Sigma }m_j`$+g), 64 cumulative yields (c), and 11 supra cumulative yields, when the addend may exceed the determination error (c). 107 yields from natW are presented, of which, 6 independent yields (i), 9 metastable state yields (m), 5 yields of metastable and ground states ($`\mathrm{\Sigma }m_j`$+g), 86 cumulative yields (c), and 1 supra cumulative yield (c).

From Tables 1 and 2 one can see that the experimental errors are ranging within $``$(6-35)%. The main contribution to the total error is from uncertainties in the nuclear data, namely, in the absolute quantum yields and cross sections of the monitor reactions.

## Comparison with experimental data obtained elsewhere

Table 3 and Fig. 9 compare some of the present results with experimental data of other laboratories published in .

| Continuation of Table 1 | | | |
| --- | --- | --- | --- |
| Product | T1/2 | Type | Yield (mb) |
| 201Tl | 72.912h | c | 43.7$`\pm `$ 2.9 |
| 200Tl | 26.1h | c | 40.6$`\pm `$ 2.6 |
| 200Tl | 26.1h | i(m+g) | 22.7$`\pm `$ 1.5 |
| 199Tl | 7.42h | c | 38.5$`\pm `$ 5.2 |
| 198m1Tl | 1.87h | i(m1+m2) | 17.6$`\pm `$ 3.6 |
| 198Tl | 5.30h | c | 35.9$`\pm `$ 5.0 |
| 196mTl | 84.6m | i(m) | 34.8$`\pm `$ 4.4 |
| 203Hg | 46.612d | c | 4.03$`\pm `$ 0.27 |
| 197mHg | 23.8h | i(m) | 10.7$`\pm `$ 0.7 |
| 195mHg | 41.6h | i(m) | 13.6$`\pm `$ 2.0 |
| 193mHg | 11.8h | i(m) | 18.9$`\pm `$ 2.5 |
| 192Hg | 4.85h | c | 35.2$`\pm `$ 2.8 |
| 198mAu | 54.48h | i(m) | 1.01$`\pm `$ 0.14 |
| 198Au | 64.684h | i(m+g) | 2.11$`\pm `$ 0.22 |
| 198Au | 64.684h | i | 1.09$`\pm `$ 0.30 |
| 196Au | 6.183d | i(m1+m2+g) | 4.13$`\pm `$ 0.35 |
| 195Au | 186.098d | c | 48.7$`\pm `$ 5.5 |
| 194Au | 38.020h | i(m1+m2+g) | 7.06$`\pm `$ 0.75 |
| 192Au | 4.94h | c | 46.9$`\pm `$ 6.6 |
| 192Au | 4.94h | i(m1+m2+g) | 11.6$`\pm `$ 1.7 |
| 191Pt | 69.6h | c | 40.1$`\pm `$ 4.4 |
| 189Pt | 10.87h | c | 46.8$`\pm `$ 4.8 |
| 188Pt | 10.2d | c | 40.5$`\pm `$ 2.9 |
| 190Ir | 11.78d | c | 0.69$`\pm `$ 0.06 |
| 188Ir | 41.5h | c | 43.2$`\pm `$ 3.2 |
| 188Ir | 41.5h | i(m+g) | 2.93$`\pm `$ 0.69 |
| 186Ir | 16.64h | c | 20.8$`\pm `$ 1.9 |
| 185Ir | 14.4h | c | 34.8$`\pm `$ 2.3 |
| 184Ir | 3.09h | c | 39.5$`\pm `$ 3.0 |
| 185Os | 93.6d | c | 41.8$`\pm `$ 2.8 |
| 183mOs | 9.9h | i(m) | 23.2$`\pm `$ 1.5 |
| 182Os | 22.1h | c | 42.0$`\pm `$ 2.8 |
| 183Re | 70.0d | c | 41.7$`\pm `$ 2.9 |
| 182Re | 12.7h | c | 45.2$`\pm `$ 3.7 |
| 181Re | 19.9h | c | 43.1$`\pm `$ 5.9 |
| 179Re | 19.7m | c | 47.8$`\pm `$ 4.2 |
| 177W | 2.25h | c | 30.1$`\pm `$ 3.5 |
| 176W | 2.30h | c | 30.8$`\pm `$ 4.3 |
| 176Ta | 8.09h | c | 34.5$`\pm `$ 3.6 |
| 173Ta | 3.14h | c | 31.0$`\pm `$ 3.9 |
| 172Ta | 36.8m | c | 17.3$`\pm `$ 2.3 |
| 175Hf | 70.0d | c | 31.3$`\pm `$ 2.3 |
| 173Hf | 23.6h | c | 28.4$`\pm `$ 2.6 |
| 172Hf | 683.017d | c | 24.1$`\pm `$ 1.6 |
| 171Hf | 12.1h | c | 18.2$`\pm `$ 2.8 |
| 170Hf | 16.01h | c | 22.1$`\pm `$ 6.8 |
| 172Lu | 6.7d | c | 23.9$`\pm `$ 1.7 |
| 172Lu | 6.7d | i(m+g) | 0.19$`\pm `$ 0.05 |
| 171Lu | 8.24d | c | 26.1$`\pm `$ 1.8 |
| 170Lu | 48.288h | c | 21.7$`\pm `$ 2.9 |

| Continuation of Table 1 | | | |
| --- | --- | --- | --- |
| Product | T1/2 | Type | Yield (mb) |
| 169Lu | 34.06h | c | 18.6$`\pm `$ 1.2 |
| 169Yb | 32.026d | c | 20.9$`\pm `$ 1.5 |
| 166Yb | 56.7h | c | 16.1$`\pm `$ 1.1 |
| 167Tm | 9.25d | c | 19.4$`\pm `$ 4.0 |
| 165Tm | 30.06h | c | 14.4$`\pm `$ 1.4 |
| 160Er | 28.58h | c | 8.8$`\pm `$ 0.6 |
| 157Dy | 8.14h | c | 5.73$`\pm `$ 0.45 |
| 155Dy | 9.90h | c | 3.66$`\pm `$ 0.27 |
| 155Tb | 5.32d | c | 4.16$`\pm `$ 0.39 |
| 153Tb | 56.16h | c | 2.52$`\pm `$ 0.25 |
| 152Tb | 17.50h | c | 2.10$`\pm `$ 0.17 |
| 153Gd | 241.6d | c | 2.60$`\pm `$ 0.23 |
| 149Gd | 9.28d | c | 2.24$`\pm `$ 0.18 |
| 146Gd | 48.27d | c | 1.26$`\pm `$ 0.09 |
| 147Eu | 24.0d | c | 0.98$`\pm `$ 0.31 |
| 146Eu | 4.59d | c | 1.63$`\pm `$ 0.11 |
| 146Eu | 4.59d | i | 0.37$`\pm `$ 0.05 |
| 143Pm | 265.0d | c | 1.02$`\pm `$ 0.13 |
| 139Ce | 137.64d | c | 0.83$`\pm `$ 0.06 |
| 121mTe | 154.0d | i(m) | 0.44$`\pm `$ 0.04 |
| 121Te | 16.78d | c | 1.11$`\pm `$ 0.11 |
| 119mTe | 4.7d | i(m) | 0.40$`\pm `$ 0.04 |
| 120mSb | 5.76d | i(m) | 0.54$`\pm `$ 0.05 |
| 114mIn | 49.51d | i(m1+m2) | 0.95$`\pm `$ 0.19 |
| 110mAg | 249.79d | i(m) | 1.12$`\pm `$ 0.09 |
| 106mAg | 8.28d | i(m) | 0.89$`\pm `$ 0.08 |
| 105Ag | 41.29d | c | 0.65$`\pm `$ 0.12 |
| 105Rh | 35.36h | c | 4.63$`\pm `$ 0.54 |
| 101mRh | 4.34d | i(m) | 1.29$`\pm `$ 0.16 |
| 103Ru | 39.26d | c | 3.84$`\pm `$ 0.26 |
| 96Tc | 4.28d | i(m+g) | 1.20$`\pm `$ 0.09 |
| 95Tc | 20.0h | c | 1.38$`\pm `$ 0.13 |
| 96Nb | 23.35h | i | 2.31$`\pm `$ 0.19 |
| 95Nb | 34.975d | c | 5.41$`\pm `$ 0.34 |
| 95Nb | 34.975d | i(m+g) | 3.03$`\pm `$ 0.20 |
| 95Zr | 64.02d | c | 2.34$`\pm `$ 0.15 |
| 89Zr | 78.41h | c | 2.30$`\pm `$ 0.16 |
| 88Zr | 83.4d | c | 0.76$`\pm `$ 0.08 |
| 90mY | 3.19h | i(m) | 4.82$`\pm `$ 0.39 |
| 88Y | 106.65d | c | 4.03$`\pm `$ 0.27 |
| 88Y | 106.650d | i(m+g) | 3.41$`\pm `$ 0.25 |
| 87Y | 79.8h | c | 2.94$`\pm `$ 0.23 |
| 85Sr | 64.84d | c | 2.76 $`\pm `$ 0.22 |
| 86Rb | 18.631d | i(m+g) | 5.48 $`\pm `$ 0.66 |
| 83Rb | 86.2d | c | 3.46 $`\pm `$ 0.28 |
| 82mRb | 6.472h | i(m) | 2.73 $`\pm `$ 0.30 |
| 82Br | 35.3h | i(m+g) | 2.17 $`\pm `$ 0.14 |
| 75Se | 119.77d | c | 1.34 $`\pm `$ 0.09 |
| 74As | 17.77d | i | 1.86 $`\pm `$ 0.18 |
| 59Fe | 44.503d | c | 0.91 $`\pm `$ 0.08 |
| 65Zn | 244.26d | c | 0.79 $`\pm `$ 0.19 |
| 46Sc | 83.81d | i(m+g) | 0.35 $`\pm `$ 0.06 |

| Continuation of Table 2 | | | |
| --- | --- | --- | --- |
| Product | T1/2 | Type | Yield (mb) |
| 147Eu | 24.0d | c | 22.4 $`\pm `$ 2.0 |
| 146Eu | 4.59d | c | 23.0 $`\pm `$ 1.9 |
| 146Eu | 4.59d | i | 3.62 $`\pm `$ 0.31 |
| 145Eu | 5.93d | c | 17.8 $`\pm `$ 1.6 |
| 139Nd | 5.5h | c | 2.87 $`\pm `$ 0.43 |
| 139Ce | 137.64d | c | 19.8 $`\pm `$ 1.6 |
| 135Ce | 17.7h | c | 17.8 $`\pm `$ 1.5 |
| 132Ce | 3.51h | c | 16.3 $`\pm `$ 2.7 |
| 132La | 4.8h | c | 14.5 $`\pm `$ 1.6 |
| 131Ba | 11.50d | c | 16.2 $`\pm `$ 1.3 |
| 126Ba | 100m | c | 7.9 $`\pm `$ 1.1 |
| 129Cs | 32.06h | c | 18.7 $`\pm `$ 1.6 |
| 127Xe | 36.4d | c | 15.4 $`\pm `$ 1.3 |
| 125Xe | 16.9h | c | 14.2 $`\pm `$ 1.2 |
| 123Xe | 2.08h | c | 15.6 $`\pm `$ 1.3 |
| 122Xe | 20.1h | c | 11.7 $`\pm `$ 1.0 |
| 121Te | 16.78d | c | 10.7 $`\pm `$ 1.1 |
| 119Te | 16.03h | c | 9.17 $`\pm `$ 0.74 |
| 119mTe | 4.7d | i(m) | 1.97 $`\pm `$ 0.17 |
| 117Te | 62m | c | 8.81 $`\pm `$ 0.77 |
| 118mSb | 5.0h | i(m) | 1.08 $`\pm `$ 0.22 |
| 115Sb | 32.1m |  | 9.85 $`\pm `$ 0.88 |
| 113Sn | 115.09d |  | 7.55 $`\pm `$ 0.67 |
| 111In | 2.8049d |  | 7.44 $`\pm `$ 0.74 |
| 110mIn | 4.9h | i(m) | 3.29 $`\pm `$ 0.29 |
| 109In | 4.2h | c | 5.12 $`\pm `$ 0.43 |
| 106mAg | 8.28d | i(m) | 1.70 $`\pm `$ 0.16 |
| 105Ag | 41.29d | c | 5.33 $`\pm `$ 0.69 |
| 100Pd | 87.12h | c | 1.24 $`\pm `$ 0.27 |
| 100Rh | 20.8h | c | 3.97 $`\pm `$ 0.44 |
| 100Rh | 20.8h | i | 2.68 $`\pm `$ 0.28 |
| 99mRh | 4.7h | c | 2.41 $`\pm `$ 0.28 |
| 97Ru | 69.6h | c | 3.13 $`\pm `$ 0.28 |
| 96Tc | 4.28d | i(m+g) | 1.73 $`\pm `$ 0.20 |
| 93mMo | 6.85h | i(m) | 1.61 $`\pm `$ 0.13 |
| 90Nb | 14.6h | c | 2.58 $`\pm `$ 0.22 |
| 89Zr | 78.41h | c | 3.46 $`\pm `$ 0.28 |
| 88Zr | 83.4d | c | 2.56 $`\pm `$ 0.27 |
| 88Y | 106.65d | c | 3.49 $`\pm `$ 0.34 |
| 88Y | 106.65d | i(m+g) | 1.56 $`\pm `$ 0.22 |
| 87Y | 79.8h | c | 4.13 $`\pm `$ 0.34 |
| 83Sr | 32.41h | c | 1.96 $`\pm `$ 0.93 |
| 84Rb | 32.77d | i (m+g) | 1.31 $`\pm `$ 0.14 |
| 83Rb | 86.2d | c | 3.34 $`\pm `$ 0.58 |
| 82mRb | 6.472h | i(m) | 1.89 $`\pm `$ 0.17 |
| 77Kr | 74.4m | c | 1.71 $`\pm `$ 0.18 |
| 75Se | 119.77d | c | 2.38 $`\pm `$ 0.22 |
| 73Se | 7.15h | c | 1.03 $`\pm `$ 0.11 |
| 74As | 17.77d | c | 1.38 $`\pm `$ 0.16 |
| 69mZn | 13.76h | i(m) | 0.42 $`\pm `$ 0.038 |
| 54Mn | 312.12d | i | 2.51 $`\pm `$ 0.42 |
| 51Cr | 27.704d | c | 4.5 $`\pm `$ 1.4 |

| Continuation of Table 2 | | | |
| --- | --- | --- | --- |
| Product | T1/2 | Type | Yield (mb) |
| 48V | 15.973d | c | 0.557 $`\pm `$ 0.062 |
| 48Sc | 43.67h | i | 0.668 $`\pm `$ 0.091 |
| 43K | 22.3h | c | 0.681 $`\pm `$ 0.084 |
| 28Mg | 20.91h | c | 0.91 $`\pm `$ 0.089 |
| 24Na | 14.959h | c | 4.09 $`\pm `$ 0.34 |
| 7Be | 53.29d | i | 8.7 $`\pm `$ 1.0 |

## Simulation of experimental results

Simulation techniques are of essential importance when forming the set of nuclear constants to be used in designing the ADS facilities, because they are universal and save much time and labour. At the same time, the present-day accuracy and reliability of the simulated results are inferior to experiment. Besides, the simulation codes are of different abilities to work when used to study the reactions that are of practical importance.

Considering the above, the present work is primarily aimed at verifying the simulation codes used most extensively for the above purpose with a view to not only estimating their ability to work when applied to the issues discussed here, but also opening up ways to improve them.

The following eight simulation codes were examined to meet these requirements:

* the CEM95 cascade-exciton code ,
* the CASCADE cascade-evaporation-fission-transport code ,
* the INUCL cascade-preequilibrium-evaporation-fission code ,
* the LAHET (ISABEL and Bertini options) cascade-evaporation-fission code ,
* the YIELDX semi-phenomenological code ,
* the CASCADE/INPE cascade-preequilibrium-evaporation-fission-transport code ,
* the CEM2k cascade-exciton code , a last modification of the CEM95 code,

Contrary to the simulation results, the experimental data include not only the independent, but also (and mainly) cumulative and the supra cumulative yields of residual product nuclei. To get a correct comparison between the experimental and simulation results, theoretical cumulative yields must be calculated on the basis of the simulated independent yields.

Since any branched isobaric chain can be presented to be a superposition of a few linear chains, the simulated cumulative and supra cumulative yields of a n-th nuclide can be calculated as

$$\sigma _n^{cum}=\sigma _n^{ind}+\underset{i=1}{\overset{n1}{}}\sigma _i^{ind}\underset{j=i}{\overset{n1}{}}\nu _j\text{,}$$
(7)
$$\begin{array}{c}\sigma _n^{cum^{}}=\sigma _n^{ind}+\frac{\lambda _{n1}}{\lambda _{n1}\lambda _n}\nu _{n1}\times \left[\sigma _{n1}^{ind}+\underset{i=1}{\overset{n2}{}}\left(\sigma _i^{ind}\underset{j=i}{\overset{n2}{}}\nu _j\right)\right]\hfill \end{array}\text{ .}$$
(8)
The branching ratios of the decay chains were retrieved from . To get a correct comparison between results by different codes, the calculations were renormalized to unified cross sections for proton-nucleus inelastic interactions from .

If an experiment-simulation difference of not above 30% (0.77$`<\sigma _{calc}/\sigma _{exp}<`$1.3) is taken to be the coincidence criterion , the simulation accuracy can be presented to be the ratio of the number of such coincidences to the number of the comparison events. The 30% level meets the accuracy requirements of the cross sections for nuclide production to be used in designing the ADS plants, according to . The mean simulated-to-experimental data ratio can be used as another coincidence criterion:

$$<F>=10^{\sqrt{<\mathrm{log}(\sigma _{cal,i}/\sigma _{exp,i})^2>}},$$
(9)
with its standard deviation

$$S(<F>)=<\left(\mathrm{log}\left(\sigma _{cal,i}/\sigma _{exp,i}\right)\mathrm{log}(<F>)\right)^2>,$$
(10)
where $`<>`$ designates averaging over all NS number of the experimental and simulated results used in a comparison.

The mean ratio $`<F>`$ together with its standard deviation $`S(<F>)`$ defines the interval \[$`<F>:S(<F>)`$ , $`<F>\times S(<F>)`$\] that covers about 2/3 of the simulation-to-experiment ratios.

The two criteria are considered sufficient to derive conclusions about the predictive power of a given code. The default options were committed to practical usage of the simulation codes.

## Comparison of data with simulation results

The results obtained with the above-mentioned codes are presented in:

* Figs. 5 and 6, that show results of a detailed comparison between the simulated and experimental radioactive product yields;
* Figs. 7 and 8, that show the simulated mass distributions of reaction products together with the measured cumulative (and supra cumulative) yields of the products that are at an immediate proximity to the stable isobar of a given mass (the sum of such yields from either sides in case both left- and right-hand branches of the chain are present). Obviously, the displayed simulation results do not contradict the experimental data if calculated values run above the experimental data and follow a general trend of the latter. This is because the direct $`\gamma `$-spectroscopy method used here identifies only radioactive products, that, as a rule, represents a significant fraction of the total mass yield, but, should a stable isobar of the given mass be produced, the $`\gamma `$-spectroscopy data are never equal to the total mass yield;
* In Fig. 9, that shows the experimental and simulated independent yields of reaction products in the form of isotopic mass distributions for several elements.

Table 4 presents the statistics of our comparison between the experimental and simulated reaction product yields in the thin 208Pb and natW samples irradiated by 1.0 GeV and 2.6 GeV protons, respectively. Namely, it shows the total number of measured yields, NT; the number of the measured yields selected to compare with calculations, NG; the number of the product nuclei whose yields were simulated by a particular code, NS; the number of the comparison events when the simulated data differ from the experimental results by not above 30%, NC1.3; the number of the comparison events when the simulated data differ from the experimental results by not more than a factor of 2.0, NC2.0.

Since about 30% of all measured secondary nuclei are not spallation reaction products, an important criterion of the codes is their ability to simulate the high-energy fission and fragmentation processes. Among the codes used here, LAHET, CASCADE, INUCL, CASCADE/INPE, and YIELDX simulate both spallation and fission. The CEM95 and CEM2k codes simulate spallation only, which is explicitly reflected in a smaller number of the products simulated (the parameter NS in Table 4 and in the shapes of the simulation curves in Figs. 5-8.

The following conclusions follow from our analysis of the experiment-to-simulation comparison results presented in Table 4 and in Figs. 5-9:

1. Generally, all codes can quite adequately simulate the weak spallation reactions (the A$``$180 products for 208Pb and the A$``$150 products for natW), with the simulation results differing from experimental data within a factor of 2.
2. In the deep spallation region (150$`<`$A$`<`$180 for 208Pb and 110$`<`$A$`<`$150 for natW), the simulation codes are of very different predictive powers, namely,

* the LAHET (when not shown explicitly as “Bertini”, all results by LAHET are of the ISABEL option), CEM2k, CASCADE/INPE, and YIELDX predictions are actually the same as the experimental data;
    * the CASCADE code simulates the A$`>`$160 product yields adequately. Below A = 160, however, the simulated data get underestimated progressively (up to a factor of 5) compared with experiment (see Fig. 7);
    * the INUCL code underestimates the yields of all the products by a factor of 2-10 in all the above mass ranges (see Figs. 7 and 8).
3. In the mass range characteristic of the fission products (50$`<`$A$`<`$150 for 208Pb and 30$`<`$A$`<`$110 for natW), the INUCL code predictions are in the best agreement with experiment when describing the yields from 208Pb. As a rule, the INUCL-simulated results differ from the data by not above a factor of 1.5. In the case of natW, however, the prediction quality deteriorates substantially. The LAHET-simulated yields are underestimated by a factor of 1.5-10.0 for Pb (Figs. 5 and 7) in the whole fission product mass region and for A$`<`$60 in the case of W (Figs. 6 and 8) but are oversetimated several times for fission fragments with A$`>`$60 from W. The YIELDX-simulated yields are either under- or over-estimated by a factor of up to 30 without showing any physical regularities. The CASCADE/INPE-simulated yields of the 130$`<`$A$`<`$150 reaction products are strongly underestimated (up to 1-2 orders of magnitude), while the simulated 40$`<`$A$`<`$130 product yields agree with the data within a factor of 2, as a rule. Generally, all the codes exhibit the feature noted above for INUCL, namely, the yield prediction quality in the case of natW is much worse compared with 208Pb, probably, because the fission cross sections of high-excited compound nuclei with very low fissility are difficult to calculate.
4. The last version of the improved cascade-exciton model code, CEM2k , shows the best agreement with the 1 GeV Pb-data in the spallation region, especially for the isotopic mass distributions (Fig. 9). At 2.6 GeV (W-target), it overestimates the expected experimental fission cross section of about 41 mb  by a factor of 6. This overestimation of the fission cross section causes an underestimation of the yield of nuclei which are most likely to fission (with a very low fissility) at the evaporation stage of a reaction, after the cascade and preequilibrium stages , i.e., for 147 $`<`$ A $`<`$ 175 (see Fig. 6). Similar disagreement with the 2.6 GeV W-data one can see as well for LAHET and CEM95, that is also related with an overestimation of the fission cross section at 2.6 GeV (see Figs. 6 and 8). The code CEM2k is still under development, its problem with the overestimation of fission cross sections at energies above 1 GeV has yet to be solved, and it has to be complemented with a model of fission fragment production, to be able to describe as well fission products.

## Conclusion

The trends shown by the advances in the nuclear transmutation of radioactive wastes and Spallation Neutron Source (SNS) facilities permit us to expect that the accumulation and analysis study of nuclear data for ADS facilities will have the same rise of academic interest and practical commitments as in the nuclear reactor data during the last five decades. Therefore, the experimental data on the yields of the proton-induced reaction products as applied to the ADS and SNS main targets and structure materials are urgent to accumulate. It should be emphasized that the charge distributions in the isobaric decay chains are important to study as well. The data thus obtained would make it possible, first, to raise the information content of the comparisons between the experimental and simulated results and, second, to lift the uncertainties in experimental determination of the cumulative yields by establishing unambiguous relations between $`\sigma ^{cum}`$ and $`\sigma ^{cum^{}}`$ for many of the reaction product masses.

Regarding the codes benchmarked here, one may conclude that none of them agree well with the data in the whole mass region of product nuclides and all should be improved further. The new CEM2k code developed recently at Los Alamos  agrees with our data in the spallation region the best of the codes tested. But CEM2k has yet to be completed by a model of fission fragmentation, to become applicable in the fission-product region as well.

## Acknowledgement

The authors are grateful to Prof. Vladimir Artisyuk (Tokyo Institute of Technology) for consulting and discussion our results on simulation-to-experimental comparisons.

The work was carried out under the ISTC Project #839 supported by the European Community, Japan (JAERI), Norway and, partially, by the U. S. Department of Energy.