text:
# Interplay of linear and nonlinear impurities in the formation of stationary localized states

## I introduction

The discrete nonlinear Schrödinger equation (DNLSE) describes a number of phenomena in condensed matter physics, nonlinear optics and other fields of physics. The DNLSE in one dimension is written as

$$i\frac{dC_n}{dt}=ϵ_nC_n+V(C_{n+1}+C_{n1})+\chi _nf_n(|C_n|)C_n,$$
(1)
where $`n=1,2,3,\mathrm{},N`$. The $`ϵ_n`$ is the linear site energy at the site $`n`$, $`V`$ is the nearest-neighbor hopping matrix element, and $`\chi _n`$ is the nonlinear parameter associated with the $`n`$th site of the lattice. Since $`_n|C_n|^2`$ can be made unity by choosing appropriate initial conditions, $`|C_n|^2`$ may be considered to be the probability of finding a particle at the $`n`$th site in the lattice. Consequently, $`C_n`$ may be interpreted as the probability amplitude of a particle at the $`n`$th site. In the context of the condensed matter physics, the nonlinear term $`f_n(|C_n|)`$ arises from the coupling of vibrations of masses at the lattice point to the motion of a quasiparticle in the limit of high frequency lattice vibrations. If the lattice oscillators are purely harmonic in nature and are independent of one another, it follows that $`f_n(|C_n|)=|C_n|^2`$. As an application of DNLSE particularly in condensed matter physics, we cite among others the exciton propagation in Holstein molecular chain. In general, the exciton propagation in quasi one dimension having short range electron-phonon interaction can be adequately modeled by the DNLSE. Other examples include nonlinear optical responses in superlattices formed by dielectric or magnetic slabs and the mean field theory of aperiodic array of twinning planes in the high $`T_c`$ superconductors.

One important feature of the DNLSE is that it can produce stationary localized (SL) states. The SL states might play an important role in the nonlinear DNA dynamics and also in energy localization in nonlinear lattices. It has been shown that the presence of impurities can produce SL states in one, two and three dimensions. The formation of SL states in one-dimensional chain with a single nonlinear impurity and that with a dimeric nonlinear impurity have been studied in detail. The same problem has also been studied starting from an appropriate Hamiltonian. The fixed point of the Hamiltonian which generates the appropriate DNLSE can also produce the correct equation, giving the formation of SL states. The formation of SL states in a perfect nonlinear chain containing a single nonlinear impurity and also in a perfect nonlinear chain with a dimeric nonlinear impurity has been studied. The formation of SL states in a perfectly nonlinear Cayley tree due to various nonlinear impurities has also been considered earlier.

On the other hand, there has been a detailed study on the formation of localized states due to linear impurities in one-dimensional chain as well as in the Cayley tree. The number of localized states do not exceed the number of impurities in case of linear systems. We note that the stationary localized states have appeared due to the presence of either linear impurities or nonlinear impurities. Thus, one may naturally ask the question of what happens when both the linear and nonlinear impurities are present at the same time in the system. In order to investigate the mixed role of linear and nonlinear impurities in the formation of SL states, we consider a simple one-dimensional chain as well as a Cayley tree consisting of a single linear impurity and a single nonlinear impurity. Furthermore, we consider a one-dimensional chain with linear and nonlinear site energies at alternate sites of the lattice.

This paper is organized as follows. The one-dimensional chain with a linear impurity and a nonlinear impurity is considered in Sec. II. In Sec. III, we consider the one-dimensional chain with linear and nonlinear site energies at alternate sites. The Cayley tree with a linear impurity and a nonlinear impurity is considered in Sec. IV. We finally summarize our observations in Sec. V.

## II One-dimensional chain with linear and nonlinear impurities

Consider a one-dimensional chain consisting of a linear impurity of strength $`ϵ`$ placed at the zeroth site and a nonlinear impurity of strength $`\chi `$ at the first site. The relevant Hamiltonian for such a system can be written as

$$H=\frac{1}{2}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}(C_nC_{n+1}^{}+C_n^{}C_{n+1})+\frac{ϵ}{2}|C_0|^2+\frac{\chi }{4}|C_1|^4,$$
(2)
where the nearest-neighbor hopping matrix element is assumed to be unity. As a possible solution for the stationary localized states, we consider

$$C_n=\varphi _ne^{iEt},$$
(3)
where $`E`$ is the stationary state energy. It is also well known that impurity states in one-dimensional system are exponentially localized. Therefore, the presence of two impurities at the two consecutive sites (site 0 and site 1) in the chain, suggests us to consider the following forms of $`\varphi _n`$:

$`\varphi _n`$ $`=`$ $`\varphi _1\left[\mathrm{sgn}(E)\eta \right]^{n1}n1`$ (4)
$`\varphi _{|n|}`$ $`=`$ $`\varphi _0\left[\mathrm{sgn}(E)\eta \right]^{|n|}n0.`$ (5)
The form of $`\varphi _n`$ given above are exact and may be derived from the Green’s function analysis. Here $`\eta `$ lies between 0 and 1 and is given as

$$\eta =\frac{1}{2}\left[|E|\sqrt{E^24}\right].$$
(6)
Since we are dealing with the localized states, $`\varphi _0`$ and $`\varphi _1`$ can be assumed to be real without forsaking the mathematical rigor. Three possible cases may be considered, i.e., $`\varphi _1=\varphi _0`$, $`\varphi _1=\varphi _0`$, and $`\varphi _1\varphi _0`$, the solutions for which are called, respectively, symmetric, antisymmetric, and asymmetric solutions. To take into account of all three possibilities, we introduce a variable $`\beta =\varphi _0/\varphi _1`$. Consequently, $`\beta =1`$ corresponds to the symmetric solutions, $`\beta =1`$ to the antisymmetric solutions, and $`\beta \pm 1`$ to the asymmetric solutions. Substituting the ansatz of Eq. (4) and the definition of $`\beta `$ to the normalization condition $`|C_n|^2=1`$, we obtain

$$|\varphi _0|^2=\frac{1\eta ^2}{1+\beta ^2}.$$
(7)
Using Eqs. (3), (4) and (6) in Eq. (2), we obtain the effective Hamiltonian for the reduced dynamical system:

$$H_{\mathrm{eff}}=\mathrm{sgn}(E)\eta +\frac{(1\eta ^2)}{(1+\beta ^2)}\left(\beta +\frac{ϵ}{2}+\frac{\chi }{4}\frac{\beta ^4(1\eta ^2)}{(1+\beta ^2)}\right).$$
(8)
We first consider the case of $`\beta =\pm 1`$. The effective Hamiltonian for this case becomes

$$H_{\mathrm{eff}}^\pm =\mathrm{sgn}(E)\eta \pm \frac{(1\eta ^2)}{2}+\frac{ϵ(1\eta ^2)}{4}+\frac{\chi }{16}\beta ^4(1\eta ^2)^2,$$
(9)
where “$`+`$” refers to the symmetric states and “$``$” corresponds to the anti-symmetric states. For a particular system, $`ϵ`$ and $`\chi `$ are constants, while $`\eta `$ is a variable which determines the energy of the localized state. The value of $`\eta `$ should be determined self-consistently and, for this reason, we treat $`\eta `$ as the dynamic variable in the effective Hamiltonian of the reduced dynamical system. The fixed-point solutions corresponding to the localized states can be found, from the variabtional principle, by the condition $`H_{\mathrm{eff}}^\pm /\eta =0`$ with $`\eta [0,1]`$. Since the form of the probability profile is rigorous, we obtain the correct values of $`\eta `$, satisfying the equation

$$\frac{4}{\chi ^\pm }=\frac{\eta (1\eta ^2)}{sgn(E)(\frac{ϵ}{2}\pm 1)\eta }=f_\pm (\eta ).$$
(10)
¿From Eq. (9) it is clear that the states will appear above the band for $`\chi ^\pm >0`$ and below the band for $`\chi ^\pm <0`$. \[Note that the band edges are at $`E/V=\pm 2`$.\] We will consider positive values of $`\chi ^\pm `$ and, hence, the sgn($`E`$) may be taken to be positive. For symmetric states, the value of $`f_+(\eta )`$ exhibits the extreme values: $`f_+(\eta )0`$ as $`\eta 0`$ and $`f_+(\eta )\mathrm{}`$ as $`\eta 1/(ϵ/2+1)`$. Thus, $`f_+(\eta )`$ increases as $`\eta `$ increases from $`\eta =0`$ and diverges at $`\eta =1/(ϵ/2+1)`$, and it becomes negative for $`\eta >1/(ϵ/2+1)`$. This means that there will always be only one solution (symmetric state) for all positive values of $`ϵ`$ and $`\chi `$. For antisymmetric states, we found that $`f_{}(\eta )0`$ as $`\eta 0`$ and $`f_{}(\eta )\mathrm{}`$ as $`\eta 1/(ϵ/21)`$. The divergence of $`f_{}(\eta )`$ is, however, confined between 0 and 1 if and only if $`ϵ>4`$. Therefore, for $`ϵ>4`$, $`f_{}(\eta )`$ increases as $`\eta `$ increases from $`\eta =0`$, diverges at $`\eta =1/(ϵ/21)`$, and finally becomes negative for $`\eta >1/(ϵ/21)`$. This implies that there will always be one antisymmetric state for $`ϵ>4`$ for any value of $`\chi >0`$. However, the situation is different for $`ϵ<4`$; the value of $`f_{}(\eta )`$ approaches 0 as $`\eta 0`$ and $`\eta 1`$ and remains positive and finite for $`0<\eta <1`$. This implies that $`f_{}(\eta )`$ has at least one maximum in the range of $`0<\eta <1`$. We confirm by a graphical analysis that there is only one maximum. Therefore, there will be a critical value of $`\chi `$ separating two antisymmetric states from the region of no antisymmetric state for $`ϵ<4`$. The solid line in Fig. 1 is such a critical line and is denoted by $`\chi _{cr}^a`$. Thus, there are three regions in the entire plane of $`ϵ`$ versus $`\chi `$, having different number of symmetric and antisymmetric SL states. The region of $`ϵ>4`$ has one symmetric and one antisymmetric SL states, while the region of $`ϵ<4`$ is divided into two regions by the critical line $`\chi _{cr}^a`$, below of which has only one symmetric state and above of which has one symmetric state and two antisymmetric states. Therefore, the maximum number of SL states contributed by the symmetric and anti-symmetric solutions is three.

Next, we consider the asymmetric solutions, i.e., for $`\beta \pm 1`$. The effective Hamiltonian in Eq. (6) should be considered as a function of two variables, $`\eta `$ and $`\beta `$, while $`ϵ`$ and $`\chi `$ remain constant for a particular system. The fixed-point solutions corresponding to the SL states can be obtained from the following two extreme conditions

$`H_{\mathrm{eff}}/\eta `$ $`=`$ $`0`$ (11)
$`H_{\mathrm{eff}}/\beta `$ $`=`$ $`0.`$ (12)
After a trite algebra, the fixed-point equations reduce to

$`(1+\beta ^2)(1+\beta ^22\eta \beta ϵ\eta )\chi \eta \beta ^4(1\eta ^2)`$ $`=`$ $`0`$ (13)
$`(1\beta ^4)ϵ\beta (1+\beta ^2)+\chi \beta ^3(1\eta ^2)`$ $`=`$ $`0.`$ (14)
These coupled equations can not be solved analytically and, accordingly, we analyze them numerically to find the number of asymmetric SL states (number of solutions with $`\eta [0,1]`$).

It is found that there are three critical lines separating different number of states in the plane of $`ϵ`$ vs $`\chi `$ denoted by the dashed curves in Fig. 1. Maximum number of asymmetric SL states obtained is two. Considering all possible states we find that maximum number of SL states is five. The “s”, “a” and “as” stand for symmetric, antisymmetric and asymmetric states, respectively. For example, the region marked as “1s+2a+2as” has one symmetric SL state, two antisymmetric SL states, and two asymmetric SL states. It should be noted that the number of states changes one by one as we go from one region to another in the phase space. As an example, as we cross from the region “1s+1as” to the region “1s+2a+1as”, the number of states increases by two; however, the number of states on the critical line which separates two regions is “1s+1a+1as”, yielding the increment of the number of states one by one. The maximum number of SL states with two nonlinear impurities was also found to be five. We thus conclude that one linear impurity and one nonlinear impurity play a similar role to the two nonlinear impurities as far as the maximum number of SL states are concerned. However, the phase diagrams for two cases are different.

Figure 2 shows the variation of energy of the states as a function of $`\chi `$ for a fixed value of $`ϵ=2`$. As $`\chi `$ increases, the energy of the state increases (goes away from the band edge) and, beyond $`\chi 10.39`$, two more states appear as expected from the phase diagram in Fig. 1. The energy of one of the two new states increases and that of the rest decreases and approaches to the band edge as $`\chi `$ increases. Therefore, for $`\chi <10.39`$, the localization length of the state, i.e., the characteristic length of the profile, decreases as $`\chi `$ increases, whereas, for $`\chi >10.39`$, the localization length of one of the three states increases as $`\chi `$ increases. Thus, for very large $`\chi `$, two of the states are strongly localized and the rest is weakly localized.

## III Linear and nonlinear site energies at alternate sites

In this section, we consider a one-dimensional chain with alternate site energies; the even sites have a linear energy $`ϵ`$ and the odd sites have a nonlinear energy $`\chi |C_n|^2`$. Such a system can be described by the Hamiltonian given as

$`H`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}(C_nC_{n+1}^{}+C_n^{}C_{n+1})`$ (16)
$`+{\displaystyle \frac{ϵ}{2}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}|C_{2n}|^2+{\displaystyle \frac{\chi }{4}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}|C_{2n+1}|^4.`$
For $`\chi =0`$, the system reduces to a periodically modulated linear system. The translational symmetry of the systems are preserved for both $`\chi =0`$ and $`\chi 0`$. Using $`C_n=\varphi _n\mathrm{exp}(iEt)`$, we obtain the Hamiltonian in terms of $`\varphi _n`$. It is, however, not possible to find the exact ansatz for the localized states in this system, although a certain rational ansatz can be considered. For example, the onsite peaked solution, intersite peaked or intersite dipped solutions are possible. We will consider each of these in the subsequent subsections.

### A Onsite peaked solution

Since the system has linear and nonlinear sites at regular intervals, the solution may have a peak either at the linear site or at the nonlinear site. We first consider the solution peaked at the linear site, i.e. the solution is peaked at the central site of the system. We therefore consider the ansatz $`\varphi _n=\varphi _0\eta ^{|n|}`$. Using this and the normalization condition, we obtain the effective Hamiltonian as

$$H_{\mathrm{eff}}=\frac{2\eta }{1+\eta ^2}+\frac{ϵ}{2}\frac{(1+\eta ^4)}{(1+\eta ^2)^2}+\frac{\chi }{2}\frac{\eta ^4(1\eta ^2)^2}{(1+\eta ^2)^2(1\eta ^8)}.$$
(17)
The fixed-point equation, $`H_{\mathrm{eff}}/\eta =0`$, reduces to

$`2`$ $`(1\eta ^8)^2(1\eta ^4ϵ\eta +ϵ\eta ^3)2\chi \eta ^5(1\eta ^2)^2(1\eta ^8)`$ (18)
$`+`$ $`2\chi \eta ^3(1\eta ^4)\left((12\eta ^2)(1\eta ^8)+2\eta ^8(1\eta ^2)\right)=0`$ (19)
Solving Eq. (14), we found that there are two separate regions in the plane of $`ϵ`$ versus $`\chi `$: one with three SL states another with only one SL state, as shown in Fig. 3.

The localized solution peaked at the nonlinear site may be obtained considering that each site is shifted by one lattice site along either positive or negative direction, and the resulting fixed-point equation is give as

$`2`$ $`(1\eta ^2)(1\eta ^8)^2(1+\eta ^2+ϵ\eta )`$ (21)
$``$ $`\chi \eta (1\eta ^4)(1\eta ^{16})+\chi \eta (1\eta ^4)(4\eta ^6(1\eta ^2)`$ (22)
$``$ $`(1\eta ^2)^2(1+\eta ^8)(1+\eta ^4))=0.`$ (23)
The phase diagram for the SL states is presented in Fig. 4. We again find that there are two separate regions, similar to the case peaked at the linear site. Thus, at least two onsite peaked SL states (one of them is peaked at the linear site and another peaked at the nonlinear site) may be obtained for any finite values of $`\chi `$ and $`ϵ`$. However, the regions of three SL states peaked at the nonlinear site appears from much lower values of $`\chi `$ and $`ϵ`$, compared to the case peaked at the linear site.

### B intersite peaked and dipped solutions

For the intersite peaked and dipped solutions we use the dimeric ansatz in Eq. (4). After short calculatons, we obtain the effective Hamiltonian for the reduced dynamical system as

$`H_{eff}`$ $`=`$ $`sgn(E)\eta +{\displaystyle \frac{\beta (1\eta ^2)}{1+\beta ^2}}`$ (24)
$`+`$ $`{\displaystyle \frac{ϵ}{2}}{\displaystyle \frac{(1+\beta ^2\eta ^2)}{(1+\beta ^2)(1+\eta ^2)}}+{\displaystyle \frac{\chi }{4}}{\displaystyle \frac{(1\eta ^2)^2(\eta ^4+\beta ^4)}{(1+\beta ^2)^2(1\eta ^8)}}`$ (25)
with $`\beta `$ and $`\eta `$ as defined earlier. We first consider the symmetric and antisymmetric states, i.e., $`\beta =\pm 1`$. The fixed-point equation yields

$$\frac{4}{\chi ^\pm }=\frac{\eta }{(1\eta )(1+\eta ^2)^2}=f^\pm (\eta ),$$
(26)
where “+” referes to the symmetric states and “-” to the antisymmetric states. It should be noted that $`ϵ`$ has no role in the formation of symmetric and antisymmetric SL states in the system. ¿From Eq. (17) it is clear that only one symmetric state appears for any finite value of $`\chi `$. However, for antisymmetric states, the situation is different. The value of $`f^{}(\eta )`$ in Eq. (17) initially increases as $`\eta `$ increases and finally decreases as $`\eta `$ approaches to 1, leaving a maximum at $`\eta [0,1]`$. Thus, there are two critical values of $`\chi `$ such that there is no anti-symmetric SL state for $`\chi <18.72`$, two anti-symmetric SL states for $`18.72<\chi <32`$ and one anti-symmetric SL state for $`\chi 32`$, shown as the two vertical lines in Fig. 5.

We next consider the asymmetric states of $`\beta 1`$. In this case $`ϵ`$ plays its role, unlike the case of symmetric and antisymmetric states. The fixed-point equations are given as

$`(1+\eta ^2)(1+\beta ^2)(1\eta ^8)^2(1+\beta ^22\beta \eta )`$ (27)
$``$ $`ϵ\eta (1\beta ^4)(1\eta ^8)^2\chi \eta (1\eta ^4)^2`$ (28)
$`\times `$ $`\left((\eta ^4+\beta ^4)(1+\eta ^2+\eta ^4\eta ^6)\eta ^2(1\eta ^8)\right)=0`$ (29)
and

$`(1\eta ^4)(1+\beta ^2)(1\eta ^8)(1\beta ^2)`$ (30)
$``$ $`ϵ\beta (1+\beta ^2)(1\eta ^8)(1\eta ^2)`$ (31)
$`+`$ $`\chi \beta (1\eta ^2)(1\eta ^4)(\beta ^2\eta ^4)=0.`$ (32)
Solving Eqs. (18) and (19) we obtain various critical lines separating the $`\chi ϵ`$ plane into various regions with various number of asymmetric SL states as shown in Fig. 5. Complete phase diagram will be obtained by superposing the Figs. 3, 4, and 5. The figures are shown separately because otherwise it will appear very complicated. The maximum number of states obtained in this case is 9. We note that the presence of both linear and nonlinear sites increases the number of SL states compared to the fully nonlinear case. No SL state is obtained for $`\chi =0`$. This is expected because the system with $`\chi =0`$ reduces to a periodically modulated linear system which again may be renormalized to a perfect linear system.

## IV Cayley tree with a linear impurity and a nonlinear impurity

Here we consider a Cayley tree with a linear impurity and a nonlinear impurity embedded in it. The purpose of this work is to observe how the SL states vary according to the geometry of the underlying lattice structure. The structure of the Cayley tree for the coordination number $`Z=3`$ and the connectivity $`K=2`$ is shown in Fig. 6.

We pick a bond and assign the numbers 0 and 1 on its two ends. The neighboring sites of the site 1 are numbered with an increasing order and those of the site 0 are numbered with a decreasing order. Thus, all points with the same number are assumed to be in the same generation. The number of points in the $`n`$th generation is $`K^{n1}`$ if $`n1`$ and $`K^{|n|}`$ if $`n0`$. We note that all sites in a given generation have the same probability amplitude.

We now consider the motion of a particle on a Cayley tree of the connectivity $`K`$ with two impurities ($`ϵ_0`$ and $`ϵ_1`$) embedded at sites 0 and 1, respectively. In the tight-binding formalism with nearest-neighbor hopping only, the equations governing the motion of a particle are

$`i{\displaystyle \frac{dC_n}{dt}}`$ $`=`$ $`KC_{n+1}+C_{n1},n>1`$ (33)
$`i{\displaystyle \frac{dC_n}{dt}}`$ $`=`$ $`KC_{|n|1}+C_{|n|+1},n<0`$ (34)
$`i{\displaystyle \frac{dC_1}{dt}}`$ $`=`$ $`KC_2+C_0+ϵ_1C_1`$ (35)
$`i{\displaystyle \frac{dC_0}{dt}}`$ $`=`$ $`KC_1+C_1+ϵ_0C_0,`$ (36)
where $`C_n`$ denotes the probability amplitude at any point in the $`n`$th generation. We assume that the matrix element of the nearest-neighbor hopping is unity and that all points in a given generation have the same site energy. The normalization condition for the site amplitudes becomes

$$\underset{n=\mathrm{}}{\overset{0}{}}K^{|n|}|C_n|^2+\frac{1}{K}\underset{n=1}{\overset{\mathrm{}}{}}K^n|C_n|^2=1.$$
(37)
Solving the Eq. (20) appears to be nontrivial; however, by a simple transformation of variables (i) $`\tau =\sqrt{K}t`$, (ii) $`ϵ_0=\stackrel{~}{ϵ}_0/\sqrt{K}`$ and $`ϵ_1=\stackrel{~}{ϵ}_1/\sqrt{K}`$, and (iii) $`C_n=K^{(n1)/2}\stackrel{~}{C}_n`$ for $`n1`$ and $`C_{|n|}=K^{|n|/2}\stackrel{~}{C}_n`$ for $`n0`$ in Eq. (20), we obtain

$`i{\displaystyle \frac{d\stackrel{~}{C}_n}{d\tau }}`$ $`=`$ $`\stackrel{~}{C}_{n+1}+\stackrel{~}{C}_{n1},\text{for}n>1\text{and}n<0.`$ (38)
$`i{\displaystyle \frac{d\stackrel{~}{C}_1}{d\tau }}`$ $`=`$ $`\stackrel{~}{C}_2+{\displaystyle \frac{1}{\sqrt{K}}}\stackrel{~}{C}_0+\stackrel{~}{ϵ}_1\stackrel{~}{C}_1,`$ (39)
$`i{\displaystyle \frac{d\stackrel{~}{C}_0}{d\tau }}`$ $`=`$ $`\stackrel{~}{C}_1+{\displaystyle \frac{1}{\sqrt{K}}}\stackrel{~}{C}_1+\stackrel{~}{ϵ}_0\stackrel{~}{C}_0.`$ (40)
which are precisely the equations for a one-dimensional chain, with the normalized probability amplitude $`\stackrel{~}{C}_n`$. The normalization condition in Eq. (21) reduces to $`_{\mathrm{}}^{\mathrm{}}|\stackrel{~}{C}_n|^2=1`$. Therefore, the motion of a particle on a Cayley tree is mapped onto that on a one-dimensional chain with a bond defect between the impurity sites. In other words, the nearest-neighbor hopping matrix element between the zeroth and the first site is reduced from unity to $`1/\sqrt{K}`$. It can be shown that the Green’s function $`G_{0,0}(E)`$ calculated from Eq. (7) would yield $`\stackrel{~}{G}_{0,0}(\stackrel{~}{E}=E\sqrt{K})`$ for a Cayley tree of the connectivity $`K`$. Here, we consider a linear impurity at the zeroth site and a nonlinear impurity at the first site. Therefore, the impurities are defined as $`\stackrel{~}{ϵ}_0=\stackrel{~}{ϵ}`$ and $`\stackrel{~}{ϵ}_1=\stackrel{~}{\chi }|C_1|^2`$ with $`\stackrel{~}{\chi }=\chi \sqrt{K}`$ and $`\stackrel{~}{ϵ}=ϵ\sqrt{K}`$.

Now, following the same approach as in the case of a one-dimensional chain, we obtain all possible symmetric, antisymmetric and asymmetric SL states for various values of the impurity strengths. The full phase diagram of SL states in the plane of $`ϵ`$ versus $`\chi `$ for the Cayley tree ($`K=4`$) with a linear impurity and a nonlinear impurity embedded in it is shown in Fig. 7. We note that the maximum number of SL states in this case is 7. Thus, a linear impurity and a nonlinear impurity play the same role as that of the two nonlinear impurities as long as the maximum number of SL states is concerned. However, the phase diagram is completely different from that for a Cayley tree with two nonlinear impurities. The one-by-one increment of the number of SL states from one region to another is retained.

## V summary

Formation of SL states due to the presence of linear and nonlinear impurities is studied. Three different cases are considered. First, a one-dimensionl chain with a linear impurity and a nonlinear impurity is considered, and the maximum number of SL states is found to be five. If one considers that the maximum number of SL states due to two linear impurities was found to be two and that due to two nonlinear impurities was five, one may conclude that the presence of one linear impurity and one nonlinear impurity plays the same role as that of two nonlinear impurities as long as the maximum number of SL states is concered. Variation of the energy of SL states as a function of nonlinear strength is studied and it is found that when a pair of states of the same type appear at the same value of $`\chi `$, one of them goes away from the band edge and another approaches the band edge as $`\chi `$ increases. Second, a linear chain with linear and nonlinear site energies at alternate sites is considered and a rich phase diagrams for SL states are obtained. The maximum number of SL states are found to be 9, which is larger than that for any system studied so far in this direction. These SL states may be called self-localized states because they appear even though the translational symmetry in the system is preserved. Third, as for the case of a different lattice structure, the Cayley tree with a linear impurity and a nonlinear impurity is considered and a rich phase diagram of SL states is obtained. The maximum number of SL states is found to be 7.

## ACKNOWLEDGMENTS

This work was supported by the Korea Research Foundation Grant (KRF-2000-015-DP0101). The authors are grateful for the support.


format: text

text:
# Constraints on diffuse neutrino background from primordial black holes

## I Introduction

Some recent inflation models (e.g., the hybrid inflationary scenario ) predict the ”blue” power-spectrum of primordial density fluctuations. In turn, as is well known, the significant abundance of primordial black holes (PBHs) is possible just in the case when the density fluctuations have an $`n>1`$ spectrum ($`n`$ is the spectral index of the initial density fluctuations, $`n>1`$ spectrum is, by definition, the ”blue perturbation spectrum”).

Particle emission from PBHs due to the evaporation process predicted by Hawking  may lead to observable effects. Up to now, PBHs have not been detected, so the observations have set limits on the initial PBH abundance or on characteristics of a spectrum of the primordial density fluctuations. In particular, PBH evaporations contribute to the extragalactic neutrino background. The constraints on an intensity of this background (and, correspondingly, on an PBH abundance) can be obtained from the existing experiments with atmospheric and solar neutrinos. The obtaining of such constraints is a main task of the present paper.

The spectrum and the intensity of the evaporated neutrinos depend heavily on the PBH’s mass. Therefore, the great attention should be paid to the calculation of the initial mass spectrum of PBHs. We use in this paper the following assumptions leading to a prediction of the PBH’s mass spectrum.

1. The formation of PBHs begins only after an inflation phase when the universe returns to the ordinary radiation-dominated era. The reheating process is such that an equation of state of the universe changes almost instantaneously into the radiation type (e.g., due to the parametric resonance ) after the inflation.

2. It is assumed, in accordance with analytic calculations  that a critical size of the density contrast needed for the PBH formation, $`\delta _c`$, is about $`1/3`$. Further, it is assumed that all PBHs have mass roughly equal to the horizon mass at a moment of the formation, independently of the perturbation size.

3. Summation over all epochs of the PBH formation can be done using the Press-Schechter formalism . This formalism is widely used in the standard hierarchial model of the structure formation for calculations of the mass distribution functions (see, e.g., ).

It was shown recently that near the threshold of a black hole formation the gravitational collapse behaves as a critical phenomenon . In this case the initial mass function will be quite different from that which follows from the standard calculations of refs. .The first calculations of the PBH initial mass function (e.g., in ref.) have been done under the assumption that all PBHs form at the same horizon mass (by other words, that all PBHs form at the smallest horizon scale immediately after reheating). The initial PBH mass spectrum for the case of the critical collapse, based on the Press-Schechter formalism, was obtained in the recent work .

The calculations in the present paper are based on the standard  picture of the gravitational collapse leading to a PBH formation. The case of the critical collapse will be considered in a separate work.

The plan of the paper is as follows.

In Sec.II we give, for completeness, the brief derivation of a general formula for the initial PBH mass spectrum. The final expression is presented in a form, which is valid for an arbitrary relation between three physical values: the initial PBH mass, the fluctuation mass (i.e., the mass of the perturbed region) at a moment of the collapse, and the density contrast in the perturbed region (also at a moment of the collapse). This expression contains the corresponding results obtained in refs.  as particular cases. As in refs. , the derivation is based on the linear perturbation theory and on the assumption that a power spectrum of the primordial fluctuations can be described by a power law.

In Sec.III we derive the approximate formula for a calculation of the extragalactic neutrino background from PBH evaporations. We do not use cosmological models of the inflation and of the spectrum of primordial fluctuations, so there are two free parameters: the reheating temperature and the spectral index. At the end of the section, some examples of instantaneous neutrino spectra from the evaporation of an individual black hole are presented.

In Sec.IV we perform the normalization of the primordial power spectrum of density perturbations by COBE data on the angular power spectrum of CMB temperature fluctuations. We consider the case of flat cosmological models with nonzero cosmological constant. We show, that the result of this normalization rather weakly depends on the matter content of the present Universe, $`\mathrm{\Omega }_m=1\mathrm{\Omega }_\mathrm{\Lambda }`$.

In Sec.V we present the results of numerical calculations of the neutrino background spectra from PBH evaporations, with taking into account the effects of neutrino absorption in the space. For background neutrino energies $`10100MeV`$ we studied relative contributions to the background intensity of different cosmological redshifts (and it is shown that at high reheating temperatures the characteristic values of the redshifts are very large).

In Sec.VI we consider the possibilities of constraining the spectral index, using experiments with neutrinos of natural origin.

In Sec.VII the spectral index constraints, followed from our calculations and from available data of the neutrino experiments, are given.

## II The initial mass spectrum of PBHs

As is pointed out in the Introduction, we use the Press-Schechter formalism which allows to carry out the summation over all epochs of PBHs formation. According to this formalism, the mass spectrum of density fluctuations ( i.e., the number density of regions with mass between $`M`$ and $`M+dM`$ ) is calculated by the formula

$$n(M,\delta _c)dM=\frac{\rho _i}{M}\left|\frac{\beta }{M}(M,\delta _c)\right|dM.$$
(1)
Here, $`\beta (M,\delta _c)`$ is the fraction of regions having sizes larger than $`R`$ and density contrast larger than $`\delta _c`$ ,

$$\beta (M,\delta _c)=2_{\delta _c}^{\mathrm{}}P(M,\delta )𝑑\delta ,$$
(2)
$$P(M,\delta )=\frac{1}{\sqrt{2\pi }\sigma _R(M)}\mathrm{exp}\left(\frac{\delta ^2}{2\sigma _R^2(M)}\right).$$
(3)
In these formulas, $`\delta `$ is the initial density contrast, $`\sigma _R`$ is the standard deviation of the density contrast of the regions with a comoving size $`R`$ and mass $`M`$. The factor of $`2`$ in Eq.(2) ensures the correct normalization of $`n(M,\delta _c)`$, integrated over all $`M`$. This factor is usually added to take into account (approximately) the contribution of underdense regions to gravitationally bound objects .

It is convenient to introduce the double differential distribution $`n(M,\delta )`$, the integral over which gives the total number of fluctuated regions:

$$n=n(M,\delta )dMd\delta ,n(M,\delta _c)=n(M,\delta )d\delta .$$
(4)
Using Eqs.(1-3), for $`n(M,\delta )`$ one has the expression

$$n(M,\delta )=\sqrt{\frac{2}{\pi }}\frac{\rho _i}{M}\frac{1}{\sigma _R^2(M)}\frac{\sigma _R}{M}\left|\left(\frac{\delta ^2}{\sigma _R^2(M)}1\right)\right|\mathrm{exp}\left(\frac{\delta ^2}{2\sigma _R^2(M)}\right).$$
(5)
To obtain the mass spectrum of PBHs one must introduce the $`M_{BH}`$ variable. Besides, we will use the variable $`\delta ^{}`$, connected with $`\delta `$ by the relation

$$\delta ^{}=\delta \left(\frac{M}{M_i}\right)^{2/3}.$$
(6)
Here, $`M_i`$ is the horizon mass at the moment of a beginning of the growth of density fluctuations. This new variable is the density contrast at the moment of the collapse. The new distribution function is

$$n_{BH}(M_{BH},\delta ^{})=n(M_{BH},\delta )\frac{d\delta }{d\delta ^{}}\frac{dM}{dM_{BH}}.$$
(7)
Further, we assume that there is some functional connection between $`M_{BH}`$, $`M`$ and $`\delta ^{}`$:

$$M_{BH}=f(M,\delta ^{}).$$
(8)
In this case one can rewrite Eq.(7) in the form

$$n_{BH}(M_{BH},\delta ^{})=n(M,\delta ^{}\left(\frac{M}{M_i}\right)^{2/3})\left(\frac{M}{M_i}\right)^{2/3}\frac{1}{df(M,\delta ^{})/dM}$$
(9)
and the PBH mass spectrum is given by the integral

$$n_{BH}(M_{BH})=n_{BH}(M_{BH},\delta ^{})𝑑\delta ^{}.$$
(10)
Now, to connect the mass spectrum with the spectral index one can use the relations

$$\sigma _R=\sigma _H(M)\left(\frac{M}{M_i}\right)^{2/3},\sigma _H(M)M^{\frac{1n}{6}}.$$
(11)
In these relations $`\sigma _H(M)`$ is the horizon crossing amplitude, i.e., the standard deviation of the density contrast in the perturbed region at the moment when this fluctuation crosses horizon. The derivation of Eqs.(11) is given, for completeness, in Sec.IV.

Using Eqs.(11) in the expression (5) for $`n(M,\delta )`$ one obtains

$$n(M,\delta )=\frac{n+3}{6}\sqrt{\frac{2}{\pi }}\frac{\rho _i}{M_i^2}\frac{1}{\sigma _H}\left(\frac{M}{M_i}\right)^{4/3}\left|\frac{\delta ^2}{\sigma _R^2}1\right|e^{\delta ^2/2\sigma _R^2}.$$
(12)
Substituting this expression in the r.h.s. of the Eq.(9) and using Eq(10) one has

$$n_{BH}(M_{BH})=\frac{n+3}{6}\sqrt{\frac{2}{\pi }}\rho _i\frac{1}{M^2\sigma _H}\left|\frac{\delta _{}^{}{}_{}{}^{2}}{\sigma _H^2}1\right|e^{\frac{\delta _{}^{}{}_{}{}^{2}}{2\sigma _H^2}}\frac{1}{df(M,\delta ^{})/dM}𝑑\delta ^{}.$$
(13)
Values of $`M`$ in r.h.s of Eq.(13) are expressed through $`M_{BH}`$, $`\delta ^{}`$ by the relation (8).

Formula (13) is the final expression for the PBH mass spectrum, the main result of this Section. It is valid for any relation between PBH mass $`M_{BH}`$, mass of the original overdense region $`M`$ and the density contrast $`\delta ^{}`$.

In the particular case of a near critical collapse

$$f(M,\delta ^{})=M_i^{1/3}M^{2/3}k(\delta ^{}\delta _c)^{\gamma _k}\xi M_i^{1/3}M^{2/3},$$
(14)
therefore

$$M=M_{BH}^{3/2}M_i^{1/2}\xi ^{3/2};df/dM=\frac{2}{3}M_i^{1/2}M_{BH}^{1/2}\xi ^{3/2}.$$
(15)
Eq.(14) can be rewritten in the form, derived in the ref. , using the connection between $`M`$ and the horizon mass $`M_h`$ (which is equal to a fluctuation mass at the moment when the fluctuation crosses horizon):

$$M_h=M_i^{1/3}M^{2/3}.$$
(16)
The resulting formula for the spectrum,

$$n_{BH}(M_{BH})=\frac{n+3}{4}\sqrt{\frac{2}{\pi }}\rho _i\sqrt{M_i}M_{BH}^{5/2}_{\delta _c}^1\frac{1}{\sigma _H}\left|\frac{\delta ^{}^2}{\sigma _H^2}1\right|e^{\frac{\delta ^2}{2\sigma _H^2}}\xi ^{3/2}𝑑\delta ^{},$$
(17)
coincides, as can be easily proved, with the expression derived in ref. .

The mass spectrum formula for the Carr-Hawking collapse can be obtained analogously, with the relation $`f(M,\delta ^{})=\gamma ^{1/2}M_i^{1/3}M^{2/3}`$ (see Eq.(32) below), or from Eq.(17) using the substitutions

$$\gamma _k0,k\gamma ^{1/2},\delta _c\gamma $$
(18)
and the approximate relation

$$_\gamma ^1𝑑\delta ^{}\left(\frac{\delta ^{}^2}{\sigma _H^2}1\right)e^{\frac{\delta ^{}^2}{2\sigma _H^2}}\gamma e^{\frac{\gamma ^2}{2\sigma _H^2}}.$$
(19)
In this case one obtains the expression derived in ref. :

$$n_{BH}(M_{BH})=\frac{n+3}{4}\sqrt{\frac{2}{\pi }}\gamma ^{7/4}\rho _iM_i^{1/2}M_{BH}^{5/2}\sigma _H^1\mathrm{exp}\left(\frac{\gamma ^2}{2\sigma _H^2}\right).$$
(20)
One can see from Eqs.(17) and (20) that the PBH mass spectrum following from the Press-Shechter formalism has quasi power form in both considered cases $`\left(M_{BH}^{5/2}\right)`$. In Carr-Hawking case $`M_{BH}^{min}M_h`$ , in contrast with this in the critical collapse case it is possible that $`M_{BH}M_h`$ .

## III Formula for neutrino diffuse background from PBHs

The starting formula for a calculation of the cosmological background from PBH evaporations is

$$S(E)=n_{com}\frac{1}{4\pi a_0^2\rho _{}^2}f\left(E(1+z)\right)𝑑V_{com}.$$
(21)
Here , $`n_{com}`$ is the comoving number density of the sources (in our case the source is an evaporating PBH of the definite mass $`m`$ ), $`a_0`$ is the scale factor at present time, $`t=t_0`$, $`f(E)`$ is a differential energy spectrum of the source radiation. $`V_{com}`$ is a comoving volume of the space filled by sources, therefore

$$dV_{com}=a_0^3\frac{\rho ^2d\rho }{\sqrt{1k\rho ^2}}d\mathrm{\Omega }.$$
(22)
Here, $`k`$ is the curvature coefficient, and $`\rho `$ is the radial comoving coordinate. Using the change of the variable,

$$\frac{d\rho }{\sqrt{1k\rho ^2}}=\frac{dt}{a},$$
(23)
the comoving number density can be expressed via the initial density $`n_i`$,

$$n_{com}=n_{phys}(t_0)=n_i\left(\frac{a_i}{a}\right)^3\left(\frac{a}{a_0}\right)^3=n_i\left(\frac{a_i}{a_0}\right)^3.$$
(24)
Substituting Eqs. (22) - (24) in Eq. (21) one obtains

$$S(E)=n_i𝑑t\frac{a_0}{a}\left(\frac{a_i}{a_0}\right)^3f\left(E(1+z)\right).$$
(25)
In our concrete case the source of the radiation is a Hawking evaporation:

$$n_if\left(E(1+z)\right)=𝑑mn_{BH}(m,t)f_H(E(1+z),m).$$
(26)
Here , $`n_{BH}(m,t)`$ is the PBH mass spectrum at any moment of time, $`f_H(E,m)`$ is the Hawking function ,

$$f_H(E,m)=\frac{1}{2\pi \mathrm{}}\frac{\mathrm{\Gamma }_s(E,m)}{exp\left(\frac{8\pi GEm}{\mathrm{}c^3}\right)(1)^{2s}}.$$
(27)
Here, $`E`$ is the energy of an evaporated particle (it lies in the interval \[$`0`$,$`\mathrm{}`$\] for massless particles), $`\mathrm{\Gamma }_s(E,m)`$ is the coefficient of the absorption by a black hole of a mass $`m`$, for an particle having spin s and energy $`E`$.

In our calculations (we considered only the simplest case of a nonrotating uncharged hole) we used the approximation when $`\mathrm{\Gamma }_s`$ depends only on $`mE`$ (as is the case for massless particles), i.e. we ignore the decrease of $`\mathrm{\Gamma }_s`$ at nonrelativistic energies of evaporated particles. But, working in this massless limit, we nevertheless used in a low $`mE`$ region the exact forms of $`\mathrm{\Gamma }_s`$, which are obtained numerically (concretely, we used the plots for absorption cross section, $`\sigma \mathrm{\Gamma }_s/E^2`$, as a function of $`mE`$, for $`s=1/2`$ and $`s=1`$, from the Fig.1 of work ), rather than the low-energy limits obtained analytically in works .

The initial spectrum of PBHs is given by Eq. (20). The minimum value of PBH mass in the specrtum can be obtained from the following considerations. It is known  that the condition of a recollapse of the perturbed region is the following: radius of a maximum expansion of this region, $`R_c`$, must obey the inequality

$$R_c>\gamma ^{1/2}R_h=\gamma ^{1/2}(ct_c).$$
(28)
Here, $`R_h`$ is the horizon radius at $`t_c`$, the moment of time, when the radius of the perturbation region is equal to maximum one. It is known also  that $`t_c`$ and $`t_i`$, the initial moment of time, are connected by the relation

$$\frac{t_i}{t_c}\delta _i,$$
(29)
where $`\delta _i`$ is the initial density contrast for which there is the condition following from Eq.(28):

$$\delta _i\gamma \left(\frac{M}{M_i}\right)^{2/3}.$$
(30)
Using Eq.(28) and the formula for horizon mass at $`t_i`$,

$$M_i=\frac{4}{3}\pi (ct_i)^3\rho (t_i),$$
(31)
one obtains the connection between $`M_{BH}`$ and $`t_c`$:

$$M_{BH}=\frac{4}{3}\pi R_c^3\rho (t_c)=\gamma ^{3/2}M_i\frac{t_c}{t_i}.$$
(32)
The use of Eqs.(29,30) in Eq.(32) leads to the relation

$$M_{BH}=\gamma ^{1/2}M_i^{1/3}M^{2/3}.$$
(33)
From Eq.(33) one follows that the minimum value of a PBH mass in the PBH mass spectrum is given by the simple formula

$$M_{BH}^{min}=\gamma ^{1/2}M_i.$$
(34)
Strictly speaking, the minimum value of $`\frac{t_c}{t_i}`$ is equal to $`1`$ (corresponding to $`\delta _i^{max}=1`$). Therefore, the absolute minimum value of PBH mass, following from Eq.(32), is equal to $`\gamma ^{3/2}M_i`$. However, the use of such a value as a border value of the PBH mass spectrum is clearly inconsistent. It contradicts with the relation (33) and, in general, with the fact that, due to the exponential damping , the bulk of the PBH intensity is determined just by the minimum value of $`\delta _i`$ (given by Eq.(30)).

To take into account the existence of the minimum we must add to the initial spectrum expression the step factor $`\mathrm{\Theta }(M_{BH}\gamma ^{1/2}M_i)`$. The connection of the initial mass value $`M_{BH}`$ and the value at any moment $`t`$ is determined by the solution of the equation

$$\frac{dm}{dt}=\frac{\alpha (m)}{m^2}.$$
(35)
The function $`\alpha (m)`$ accounts for the degrees of freedom of evaporated particles and determines the lifetime of a black hole. In the approximation $`\alpha =const`$ the solution of Eq.(35) is

$$M_{BH}\left(3\alpha t+m^3\right)^{1/3}.$$
(36)
This decrease of PBH mass leads to the corresponding evolution of a form of the PBH mass spectrum. At any moment one has

$$n_{BH}(m,t)dm=\frac{m^2}{\left(3\alpha t+m^3\right)^{2/3}}n_{BH}\left((3\alpha t+m^3)^{1/3}\right)\times \mathrm{\Theta }\left[m\left((\gamma ^{1/2}M_i)^33\alpha t\right)^{1/3}\right]dm.$$
(37)
Substituting Eqs. (27), (37) in the integral in Eq.(26) , we obtain the final expression for the spectrum of the background radiation:

$`S(E)={\displaystyle \frac{c}{4\pi }}{\displaystyle }dt{\displaystyle \frac{a_0}{a}}\left({\displaystyle \frac{a_i}{a_0}}\right)^3{\displaystyle }dm{\displaystyle \frac{m^2}{(3\alpha t+m^3)^{2/3}}}n_{BH}\left[(3\alpha t+m^3)^{1/3}\right]\times `$ (38)
(39)
$`\mathrm{\Theta }[(m((\gamma ^{1/2}M_i)^33\alpha t)^{1/3}]g_\nu f_H(E(1+z),m).`$ (40)
Here, $`g_\nu `$ is a number of degrees of freedom. In the following we will interest in ($`\nu _e+\stackrel{~}{\nu }_e`$) - spectrum, so $`g_\nu =2`$.

One should note that the corresponding expressions for the spectrum in refs.  contain the factor $`\left(\frac{a_i}{a}\right)^3`$ instead of the correct factor $`\left(\frac{a_i}{a_0}\right)^3`$. It leads to a strong overestimation of large $`z`$ contributions in $`S(E)`$ (see below, Fig.5).

It is convenient to use in Eq.(38) the variable $`z`$ instead of $`t`$. In the case of flat models with nonzero cosmological constant one has

$$\frac{dt}{dz}=\frac{1}{H_0(1+z)}\left(\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_r(1+z)^4+\mathrm{\Omega }_\mathrm{\Lambda }\right)^{1/2},$$
(41)
$$\mathrm{\Omega }_r=(2.410^4h^2)^1,\mathrm{\Omega }_m=1\mathrm{\Omega }_r\mathrm{\Omega }_\mathrm{\Lambda }.$$
(42)
The factor $`(\frac{a_i}{a_0})^3`$ can be expressed through the value of $`t_{eq}`$, the moment of matter-radiation density equality:

$$\left(\frac{a_i}{a_0}\right)^3(1+z_{eq})^3\left(\frac{t_i}{t_{eq}}\right)^{3/2}=\left(\frac{2}{3}(2\sqrt{2})\right)^{3/2}H_0^{3/2}\left(2.410^4h^2\right)^{3/4}t_i^{3/2}.$$
(43)
The constant $`H_0`$, entering the Eqs.(41,43) is connected with $`h`$ by definition:

$$H_0=100h\text{ km}\text{ s}^1\text{ Mpc}^1.$$
(44)
Throughout the paper we suppose that $`h=0.67`$. One can see from Eqs.(43,44) that the factor $`\left(\frac{a_i}{a_0}\right)^3`$ does not depend on $`h`$ and $`\mathrm{\Omega }_m`$, as it must be.

Integrating over PBH’s mass in Eq.(38), one obtains finally, after the change of the variable $`t`$ on $`z`$, the integral over $`z`$:

$$S(E)=F(E,z)d\mathrm{log}_{10}(z+1).$$
(45)
In analogous calculations of the photon diffuse background integral over $`z`$ in the expression for $`S(E)`$ is cut off at $`z=z_0700`$ because for larger $`z`$ the photon optical depth will be larger than unity . In contrast with this , interactions of neutrinos with the matter can be neglected up to very high values of $`z`$. Therefore the neutrino diffuse background from PBH evaporations is much more abundant. The neutrino absorption effects are estimated below, in Sec.VI.

The evaporation process of a black hole with not too small initial mass is almost an explosion. So, for a calculation of spectra of evaporated particles with acceptable accuracy it is enough to know the value of $`\alpha `$ for an initial value of the PBH mass only. Taking this into account and having in mind the steepness of the PBH mass spectrum, we use the approximation

$$\alpha (m)=\alpha (M_{BH}^{min})=\alpha (\gamma ^{1/2}M_i),$$
(46)
and just this value of $`\alpha `$ is meant in the expressions (36)-(38).

The very detailed calculation of the function $`\alpha (m)`$ was carried out in the works . Here we use the simplified approach in which $`\alpha (m)`$ is represented by the dependence

$$\alpha =\alpha _0+\underset{i}{}a_i10^{25}\mathrm{\Theta }(b_ilog_{10}(m)).$$
(47)
Here, $`\mathrm{\Theta }(x)`$ is the step function. Coefficient $`\alpha _0`$ gives the summary contribution of $`e^+`$, $`e^{}`$, $`\nu `$ and $`\gamma `$ and is equal to $`8.4210^{25}g^3sec^1`$ . All other coefficients are collected in the Table I.

The coefficients $`b_i`$ determine the value of the PBH mass beginning from which particles of $`i`$-type can be evaporated. This value, $`M_{BH}^{b_i}`$, is obtained from the relations

$$\left(\frac{10^{13}g}{M_{BH}^{b_i}}\right)\text{ GeV}kT_{BH}^{b_i}\frac{m_i}{\kappa (s_i)},$$
(48)
$$\kappa (\frac{1}{2})4,\kappa (1)6.$$
(49)
So,

$$b_i=log_{10}M_{BH}^{b_i}.$$
(50)
Writing Eqs.(48,49) we took into account that the evaporation of massive particles becomes essential once the peak value in the energy distribution exceeds the particle rest mass. The peak value depends on the particle spin (see, e.g., ):

$$E_{peak}(s=1/2)4kT_{BH};E_{peak}(s=1)6kT_{BH}.$$
(51)
The values of masses of (u,d,s,c,b,g,$`\tau `$)-particles used in our calculations of $`b_i`$ are the same as in work , top quark mass was taken equal to $`170\text{GeV}`$. The resulting function $`\alpha (m)`$ is shown on Fig.1. For comparison, the corresponding dependence from the work  (drawn using the Eq.7 of ) is also shown.

According to Eq.(34) the minimum value of PBH mass in the mass spectrum is about the horizon mass $`M_i`$. In turn, the value of $`M_i`$ is connected with $`t_i`$, the moment of time just after reheating, from which the ”normal” radiation era and the process of PBHs formation start,

$$M_i=\frac{1}{8}\frac{M_{pl}}{t_{pl}}t_i.$$
(52)
The corresponding reheating temperature is given, approximately, by the relation (see Eq.(113))

$$kT_{RH}\left(\frac{t_i}{1s}\right)^{1/2}\text{MeV}.$$
(53)
If $`kT_{RH}`$ values lie in the interval $`(10^810^{10})\text{ GeV}`$ (just for this interval we obtain the spectral index constraints in the present paper), the corresponding $`M_i`$ values lie in the mass interval $`(10^{11}10^{15})\text{ g}`$. So, as one can see, we really need for our aims the all information about $`\alpha `$-parameter, that contains in Fig.1, including the highest values of $`\alpha `$ at $`m10^{11}\text{ g}`$.

The formula (38), as it stands, takes into account the contribution to the neutrino background solely from a direct process of the neutrino evaporation. If we suppose that particles evaporated by a black hole propagate freely, the calculation of other contributions to the neutrino background can be performed using our knowledge of particle physics .

On Fig.2 the typical result of our calculation of instantaneous neutrino spectra from evaporating black hole is shown. The spectrum of the straightforward (direct) $`(\nu _e+\stackrel{~}{\nu }_e)`$ \- emission is given by the Hawking function $`f_H(E,m)`$, Eq.(27). The $`(\nu _e+\stackrel{~}{\nu }_e)`$ spectrum arising from decays of $`(\mu ^++\mu ^{})`$ evaporated directly is calculated by the formula

$$f^{(\mu )}(E,m)=\mathrm{\Theta }(b_\mu \mathrm{log}_{10}m)g_\mu f_H(E_\mu ,m)\frac{dn_\nu (E_\mu ,E)}{dE}𝑑E_\mu ,$$
(54)
where $`dn_\nu ^\mu /dE`$ is the neutrino spectrum in a $`\mu `$ \- decay, $`g_\mu =4`$. To evaluate the electron neutrino spectrum resulted from fragmentations of evaporated quarks, we used the simplest chain

$$(u,d)\text{quarks}\pi \mu \nu _e,$$
(55)
and the formula

$$f^{(q)}(E,m)=\mathrm{\Theta }(b_q\mathrm{log}_{10}m)g_qf_H(E_q,m)\frac{dn_\pi ^q(\xi )}{d\xi }\frac{d\xi }{dE_q}\frac{dn_\mu ^\pi (E_\mu ,E_\pi )}{dE_\mu }\frac{dn_\nu ^\mu (E_\mu ,E)}{dE}𝑑E_q𝑑E_\pi 𝑑E_\mu .$$
(56)
Here, $`g_q=24`$ for a summary contribution from $`u,\stackrel{~}{u},d,\stackrel{~}{d}`$ \- quarks, $`dn_\pi ^q/d\xi `$ is the $`q\pi `$ fragmentation function, for which the simple parametrization was taken:

$$\frac{dn_\pi ^q}{d\xi }=\frac{15}{16}(\xi 1)^2\xi ^{3/2},\xi =\frac{E_\pi }{E_q},$$
(57)
and $`dn_\mu ^\pi /dE_\mu `$ is the $`\mu `$ \- spectrum in a decay $`\pi \mu +\nu _\mu `$.

Analogous formulas are used for calculations of the neutrino spectrum from other channels of the neutrino production, for instance from decays of evaporated $`W`$-bosons ($`We+\nu _e,`$ $`W\mu \nu _e`$).

The formula (57) for the fragmentation function was suggested in the paper . It leads to a simple multiplicity growth law of $`\sqrt{E}`$ (such a law follows from a naive statistical model of an jet fragmentation ). It had been shown in ref.  that Eq.(57) fits rather well the PETRA data and can be used instead of the more ingenious formula based on the leading logarithm approximation of QCD. In turn, it was shown recently  using the programme HERWIG  that the parametrizations of ref. , being rather good at low $`\xi `$, lead to an significant overestimation of the yield of large $`\xi `$ final states. This fact as well as the simplicity of Eq.(57) were decisive for our choice of the fragmentation function because, as we will show, we need just the upper estimate of the yield of quark fragmentations.

In all calculations in the present paper we neglected the contribution to neutrino fluxes from neutron decays in the jets. These decays give the distinct bumps in instantaneous $`\nu \stackrel{~}{\nu }`$-spectra at $`E1\text{ MeV}`$ . In general, neutrinos of such low energies have too small cross sections for neutrino-nucleus and neutrino-nucleon interactions (see sec.VI) and are inessential for our main aim: the constraining of the spectral index (neutrinos with energies $`10100\text{ MeV}`$ are mainly responsible for this). Besides, we will see at Sec.V that at large values of one of the parameters of our approach (reheating temperature $`T_{RH}`$) these neutrinos are not noticeable in the summary neutrino background even at $`E1\text{ MeV}`$ because of strong redshift effects.

In order to check the accuracy of our calculations of instantaneous neutrino spectra we calculated for one value of a black hole temperature the summary neutrino spectrum (excluding secondary $`\nu _\tau `$ only). The result is shown on Fig.3 together with the corresponding curve from ref.. One can see from this figure that in the practically interesting region of neutrino energies ($`E>100\text{ MeV}`$) our result differs from the Monte-Carlo result of ref. not more than on a factor of 2.

The relative contribution of different channels to the total $`\nu _e`$ spectrum strongly depends on the black hole temperature. One can see from Fig.2, that at high temperature decays of massive particles evaporated by the black hole become rather important at high energy tail of the spectrum (if the corresponding branching ratios are not too small).

The total instantaneous neutrino spectrum from a black hole evaporation is given by the sum

$$f(E,m)=g_\nu f_H(E,m)+f^{(\mu )}(E,m)+f^{(q)}(E,m)+\mathrm{},$$
(58)
and the total background neutrino spectrum is given by the same Eq.(38), except the change $`g_\nu f_H(E,m)f(E,m)`$.

## IV Normalization of the perturbation amplitude

For numerical calculations of the neutrino background we need normalization of the dispersion $`\sigma _R`$ entering the PBH mass spectrum formula. For this aim we must consider the time evolution of this dispersion. The connection of $`\sigma _R(t)`$ with a power spectrum $`P(k,t)`$ of the density perturbations is given by the expression

$`\sigma _R^2(t)={\displaystyle \underset{0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{k^3}{2\pi ^2}}P(k,t)W^2(kR){\displaystyle \frac{dk}{k}}`$ (59)
(60)
$`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{k^4}{(aH)^4}}\delta _H^2(k,t)W^2(kR){\displaystyle \frac{dk}{k}}.`$ (61)
Here, $`W(kR)`$ is the smoothing window function. The second line of this equation defines the (nonsmoothed) horizon crossing amplitude $`\delta _H(k,t)`$. In general, it depends on the time and we need the initial amplitude $`\delta _H(k,t_i)`$, whereas the normalization on COBE data (which had been performed, in particular, in the work ) determines this amplitude at present time, $`t=t_0`$. Only in the particular case of a flat cosmological model with $`\mathrm{\Lambda }=0`$ the horizon crossing amplitude is time-independent and the normalization on COBE data gives $`\delta _H(k,t_i)`$ straightforwardly.

Assuming a power law of primordial density perturbations,

$$P(k,t)k^n,$$
(62)
and choosing for $`W(kR)`$ the top-hat form,

$$W(kR)=\frac{3j_1(kR)}{kR},$$
(63)
one obtains from Eq.(59)

$$\sigma _R^2(t)=C_n^2\frac{k_{fl}^4}{(aH)^4}\delta _H^2(k_{fl},t)\frac{k_{fl}^4}{(aH)^4}\sigma _H^2(k_{fl},t).$$
(64)
In Eq.(64), $`k_{fl}`$ is the comoving wave number, characterizing the perturbed region,

$$k_{fl}=\frac{1}{R}$$
(65)
and $`C_n`$ is a constant which depends on the spectral index $`n`$. For $`n1`$ one has $`C_n4.8`$ . The fluctuation power per logarithmic interval, or variance, at present time is given by the expression

$$\frac{k^3}{2\pi ^2}P(k,t_0)=\frac{1}{H_0^4}k^4A_S^2(k)\frac{g^2(\mathrm{\Omega }_m)}{\mathrm{\Omega }_m^2}.$$
(66)
Here, $`A_S(k)`$ is an amplitude of the scalar perturbations, k-dependence of which arises from a deviation from scale invariance ($`A_S^2k^{n1}`$). The factor $`g^2(\mathrm{\Omega }_m)`$ takes into account the suppression of the growth of density perturbations due to nonzero $`\mathrm{\Omega }_\mathrm{\Lambda }`$ (density perturbations stop growing once the universe becomes cosmological constant dominated). The factor $`\mathrm{\Omega }_m^2`$ in the denominator of Eq.(66) arises because the potential fluctuations are reduced by $`\mathrm{\Omega }_m`$: due to this, for fixed COBE normalization, the power spectrum of density perturbations must be, aside from other factors, enhanced by $`\mathrm{\Omega }_m^2`$ .

For arbitrary moment of time in matter era the expression for the power spectrum is evident generalization of Eq.(66):

$$\frac{k^3}{2\pi ^2}P(k,t)=\frac{k^4}{H_0^4}A_S^2(k)\left(\frac{g(a,\mathrm{\Omega }_m)}{\mathrm{\Omega }_m}\right)^2.$$
(67)
Here, the function $`g(a,\mathrm{\Omega }_m)`$ accounting for the growth of density perturbations is given by the expression

$$g(a,\mathrm{\Omega }_m)=\frac{5\mathrm{\Omega }_m}{2a}\frac{da}{d\tau }\underset{0}{\overset{a}{}}\left(\frac{da^{}}{d\tau }\right)^3𝑑a^{},$$
(68)
where $`\tau H_0t`$. In our case

$$\frac{da}{d\tau }=\mathrm{\Omega }_m\frac{1}{a}+\mathrm{\Omega }_\mathrm{\Lambda }a^2$$
(69)
and

$$g(a,\mathrm{\Omega }_m)=\frac{5}{2}\mathrm{\Omega }_ma^{3/2}(\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }a^3)^{1/2}\underset{0}{\overset{a}{}}\frac{x^{3/2}dx}{(\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }x^3)^{3/2}}.$$
(70)
The normalization of this function is such that

$$g(a,1)=a.$$
(71)
For the present time ($`a=a_0=1`$) one has

$$g(1,\mathrm{\Omega }_m)=g(\mathrm{\Omega }_m),$$
(72)
where $`g(\mathrm{\Omega }_m)`$ is the function entering Eq.(66), as it must be. The horizon crossing amplitude $`\delta _H(k,t)`$ defined by the Eq.(59) is determined by the relation

$$\delta _H^2(k,t)=\frac{(aH)^4}{H_0^4}A_S^2(k)\left(\frac{g(a,\mathrm{\Omega }_m)}{\mathrm{\Omega }_m}\right)^2\stackrel{tt_0}{}A_S^2(k)\left(\frac{g(\mathrm{\Omega }_m)}{\mathrm{\Omega }_m}\right)^2.$$
(73)
One can see that the time dependence of $`\delta _H(k,t)`$ arises just due to nonzero value of $`\mathrm{\Lambda }`$. If $`\mathrm{\Lambda }=0`$, one has from Eq.(73), using Eq.(71),

$$\delta _H^2(k,t)|_{\mathrm{\Lambda }=0}=A_S^2(k).$$
(74)
Our normalization of perturbation amplitude on COBE data is based on two inputs.

1. Connection between CMB fluctuations and scalar density fluctuations is described  by the following relations (derived for a lowest-order reconstruction of the inflationary potential):

$$S\frac{5C_S^2}{4\pi }=0.104f_S^{(0)}(\mathrm{\Omega }_\mathrm{\Lambda })A_S^2(k_{}),$$
(75)
$$f_S^{(0)}(\mathrm{\Omega }_\mathrm{\Lambda })=1.040.82\mathrm{\Omega }_\mathrm{\Lambda }+2\mathrm{\Omega }_\mathrm{\Lambda }^2,$$
(76)
$$0.0\mathrm{\Omega }_\mathrm{\Lambda }0.8.$$
(77)
Here, $`C_2^S`$ is the scalar contribution to the angular power spectrum of CMB temperature fluctuations for $`l=2`$, $`k_{}a_0H_0`$ is the comoving wave number at the present horizon scale. The dependence $`f_S^{(0)}`$ on $`\mathrm{\Omega }_\mathrm{\Lambda }`$ is due to the integrated Sachs-Wolfe effect, i.e., due to the evolution of the potentials from last scattering surface till the present time.

2. From 4-year COBE data  one has (assuming that the scalar contribution dominates over the tensor on the large scales) the following results for small $`l`$ values :

$$\left(\frac{l(l+1)C_l^S}{2\pi }\right)^{1/2}(1.03\pm 0.07)10^5,$$
(78)
$$n=1.02\pm 0.24.$$
(79)
From Eqs.(73-78) one obtains, finally, the normalization of the horizon crossing amplitude:

$$\delta _H(k_{},t_0)\delta _H^{(COBE)}210^5\frac{1}{\sqrt{f_S^{(0)}(\mathrm{\Omega }_\mathrm{\Lambda })}}\left(\frac{g(\mathrm{\Omega }_m)}{\mathrm{\Omega }_m}\right),$$
(80)
$$A_S(k_{})=210^5\frac{1}{\sqrt{f_S^{(0)}(\mathrm{\Omega }_\mathrm{\Lambda })}}.$$
(81)
All the dependence of $`\delta _H(k,t)`$ on $`k`$ is contained in an amplitude $`A_S(k)`$ and, according to our assumption, Eq.(62), is very simple:

$$A_S(k)k^{\frac{n1}{2}}.$$
(82)
Expressing the variable $`k`$ through $`M_h`$, horizon mass at the moment when the scale $`\lambda \frac{1}{k}`$ crosses the Hubble radius, one obtains

$$A_S(M_h)\{\begin{array}{ccc}\hfill M_h^{\frac{1n}{4}}& ,\hfill & \hfill \text{radiation era ,}\\ & & \\ \hfill M_h^{\frac{1n}{6}}& ,\hfill & \hfill \text{matter era .}\end{array}$$
(83)
Writing Eqs.(83) we took into account that the comoving mass density in radiation era decreases with time , while in matter era it is conserved. Now, using the Eq.(73) one can calculate the amplitude $`\delta _H(k)`$ at the moment of matter - radiation equality:

$$\delta _H(M_{eq},t_{eq})=\frac{(a_{eq}H_{eq})^2}{H_0^2}A_S(M_{eq})\frac{g(a_{eq},\mathrm{\Omega }_m)}{\mathrm{\Omega }_m}.$$
(84)
Here, $`M_{eq}`$ is the horizon mass at $`t_{eq}`$. Since $`a_{eq}1`$, one has

$$H_{eq}^2H_0^2\frac{\mathrm{\Omega }_m}{a_{eq}^3},g(a_{eq},\mathrm{\Omega }_m)a_{eq}.$$
(85)
Therefore, according to Eq.(80),

$$\delta _H(M_{eq},t_{eq})=A_S(M_h)=2.010^5\frac{1}{\sqrt{f_S^{(0)}(\mathrm{\Omega }_\mathrm{\Lambda })}}\left(\frac{M_{eq}}{M_{h0}}\right)^{\frac{1n}{6}}.$$
(86)
Here, $`M_{h0}`$ is the present horizon mass. Using Eq.(83) one obtains the time-independent amplitude $`\delta _H(M_h)`$ in radiation era:

$$\delta _H(M_h)=210^5\frac{1}{\sqrt{f_S^{(0)}(\mathrm{\Omega }_\mathrm{\Lambda })}}\left(\frac{M_{eq}}{M_{h0}}\right)^{\frac{1n}{6}}\left(\frac{M_h}{M_{eq}}\right)^{\frac{1n}{4}}.$$
(87)
The values of horizon masses $`M_{eq}`$ and $`M_{h0}`$ in Eq.(87) depend on $`\mathrm{\Omega }_m`$:

$$M_{eq}=\mathrm{\Omega }_m^2h^41.510^{49}g;M_{h0}=\mathrm{\Omega }_mh^1610^{55}g.$$
(88)
The expression for $`\sigma _H(M_h)`$ differs from Eq.(87) only by the numerical coefficient (see Eq.(64)):

$$\sigma _H(M_h)=C_n\delta _H(M_h),$$
(89)
and $`M_h`$ now is a function of $`k_{fl}`$ .

One can see from Eqs.(87,88) that the horizon crossing amplitude $`\sigma _H(M_h)`$ , normalized on COBE data, depends on $`\mathrm{\Omega }_m`$ only through the $`\mathrm{\Omega }_\mathrm{\Lambda }`$-dependence of the $`f_S^{(0)}`$ function:

$$\sigma _H(M_h)\frac{1}{\sqrt{f_S^{(0)}(\mathrm{\Omega }_\mathrm{\Lambda })}}M_h^{\frac{1n}{4}}.$$
(90)
The relations (87) and (90) are main results of this Section. The final expression for $`\sigma _R(t)`$ is, according to Eq.(64),

$$\sigma _R(t)=\left(\frac{M_{hor}(t)}{M_{fl}(t)}\right)^{2/3}\sigma _H(M_h).$$
(91)
Here, $`M_{fl}(t)`$ is the fluctuation mass,

$$M_{fl}(t)=\frac{4}{3}\pi \left(\frac{a(t)}{k_{fl}}\right)^3\rho (t)$$
(92)
and $`M_{hor}(t)`$ is the horizon mass at a moment t. At the initial moment of time one has ($`M_{hor}(t_i)=M_i,M_{fl}(t_i)=M`$)

$$\sigma _R(t_i)\sigma _R(M)=\left(\frac{M_i}{M}\right)^{2/3}\sigma _H(M_h).$$
(93)
The expressions (11) of Sec.II can be obtained from Eqs.(90 , 93) using the connection between $`M`$ , $`M_i`$ and $`M_h`$ for radiation era (Eq.(16)).

## V Neutrino background spectra

A precise calculation of the PBH neutrino background must include also a taking into account the neutrino absorption during a travelling in the space. In the case of the photon background the most important absorption process is a pair production on neutral matter  due to which the photons from PBHs evaporated earlier than $`z700`$ are absent today. In our case, the analog of an optical depth of the universe for the neutrino emitted at a redshift z and having today an energy $`E`$ is given by the integral<sup>*</sup><sup>*</sup>* In formulas of Secs.V, VI we use the convention $`\mathrm{}=c=k=1`$.

$$\tau (z,E)=\underset{0}{\overset{z}{}}\sigma \left(E(1+z^{})\right)n(z^{})\frac{dt}{dz^{}}𝑑z^{}.$$
(94)
Here, $`\sigma (E)`$ is the neutrino interaction cross section, $`n(z)`$ is a number density of the target particles.

Two processes are potentially ”dangerous”: neutrino-nucleon inelastic scattering growing linearly with an energy, and annihilations with neutrinos of the relic background

$$\nu _e+Ne^{}+anything,$$
(95)
$$\nu _e+\stackrel{~}{\nu }_e(relic)\underset{i}{}(f_i+\stackrel{~}{f}_i).$$
(96)
Here, $`f_i`$ are charged fermions (leptons and quarks). As is known, the relic neutrino background exists, with a Planck distribution, beginning from the epoch of neutrino decoupling ($`z10^{10}`$).

The cross section of the neutrino-nucleon inelastic scattering can be approximated (in the neutrino energy interval from $`E1\text{GeV}`$ to $`E10^5\text{GeV}`$) by the simple formula:

$$\sigma _{\nu N}(E)0.510^{38}\left(\frac{E}{\text{GeV}}\right)cm^2.$$
(97)
At $`E10^5\text{ GeV}`$ the cross section grows with energy more slowly, so the use of Eq.(97) in this region leads to the overestimation of the cross section value. But we will use this formula for the entire interval of neutrino energy, having in mind that practically we need only the upper estimate of this cross section (because, as we will see, the contribution of this channel to the total absorption factor is small).

The number density of the target nucleons can be estimated using the existing constraints on fraction of critical density in baryons :

$$0.004<\mathrm{\Omega }_Bh^2<0.021.$$
(98)
From these inequalities and for $`h=0.67`$ one obtains for the present nucleon number density :

$$n_N^0=\frac{n_N(z)}{(1+z)^3}=(410^8÷210^7)cm^3.$$
(99)
In numerical calculation we used the value $`10^7`$ for $`n_N^0`$. The $`\tau _{\nu N}`$-function is given by the simple integral:

$$\tau _{\nu N}(E,z)0.510^{38}\left(\frac{E}{\text{GeV}}\right)10^7\underset{0}{\overset{z}{}}(1+z)^4\frac{dt}{dz^{}}𝑑z^{}.$$
(100)
The expression (94), as it stands, is valid only in a case when a target particle is at rest. In the case of an absorption through annihilations of evaporated neutrinos with relic neutrinos the current expression for the $`\tau _{\nu \stackrel{~}{\nu }}`$ \- function is (see, e.g., ref. )

$$\tau _{\nu \stackrel{~}{\nu }}(E,z)=\underset{0}{\overset{z}{}}(1\mathrm{cos}\theta _{\nu \stackrel{~}{\nu }})\sigma _{\nu \stackrel{~}{\nu }}(E(1+z^{}),E_\nu (1+z^{}),\theta _{\nu \stackrel{~}{\nu }})n_{\stackrel{~}{\nu }}(z^{})\frac{dt}{dz^{}}𝑑z^{}.$$
(101)
The angular brackets indicate an average over the relic antineutrino energy distribution and over the angle $`\theta _{\nu \stackrel{~}{\nu }}`$ between the two colliding particles in the cosmic frame.

The cross section of the $`\nu \stackrel{~}{\nu }`$ \- annihilation is expressed by the formula

$$\sigma _{\nu \stackrel{~}{\nu }}(E,E_{\stackrel{~}{\nu }},\theta _{\nu \stackrel{~}{\nu }})=\frac{G_F^2}{4\pi }\left[N_{eff}^{NC}(s)+N_{eff}^{CC}(s)\right]s,$$
(102)
$$s=2EE_{\stackrel{~}{\nu }}(1\mathrm{cos}\theta _{\nu \stackrel{~}{\nu }}).$$
(103)
Here, $`s`$ is a square of the summary energy of colliding particles in the center of mass system, the coefficients $`N_{eff}`$ are the effective numbers of annihilation channels (for neutral and charged currents).

Formula (102) is valid in the energy region before the $`Z`$ boson pole, i.e., when

$$sM_Z^2.$$
(104)
As one can see from Eq.(103), maximum values of $`s`$, for which we must know the cross section, are about $`ET_{\stackrel{~}{\nu }}^0(1+z)^2`$ (where $`T_{\stackrel{~}{\nu }}^0`$ is the present temperature of the relic antineutrino gas, $`z`$ is the upper limit of the integral in Eq.(101), and $`E`$ is the neutrino energy today). So, the cross section formula (102) can be used if

$$z\left(\frac{M_Z^2}{ET_{\stackrel{~}{\nu }}}\right)^{1/2}310^8\left(\frac{E}{\text{GeV}}\right)^{1/2}.$$
(105)
We will see from the numerical results that this condition is satisfied in all the practically important cases.

The coefficients $`N_{eff}`$ are given by the following expressions:

$$N_{eff}^{NC}(s)=\underset{f}{}\mathrm{\Theta }(\frac{s}{4}m_f^2)\frac{2}{3}n_f\left(18t_{3f}q_f\mathrm{sin}^2\theta _W+8q_f^2\mathrm{sin}^4\theta _W\right),$$
(106)
$$N_{eff}^{CC}=\mathrm{\Theta }(\frac{s}{4}m_e^2)\frac{16}{3}\mathrm{sin}^2\theta _W.$$
(107)
Here, $`n_f`$ is the number of colours (1 for leptons, 3 for quarks ), $`t_{3f}`$ and $`q_f`$ are the third component of the weak isospin and the electric charge ( in units of the positron charge ), respectively, $`\theta _W`$ is the electroweak mixing angle ($`\mathrm{sin}^2\theta _W0.23`$), $`\mathrm{\Theta }(x)`$ is the step function.

Using the relations

$$(1\mathrm{cos}\theta _{\nu \stackrel{~}{\nu }})^2E_{\stackrel{~}{\nu }}=\frac{4}{3}\rho _{\stackrel{~}{\nu }}\frac{1}{n_{\stackrel{~}{\nu }}},$$
(108)
$$\frac{s}{4}m_f^2ET_{\stackrel{~}{\nu }}^0m_f^2,$$
(109)
where $`\rho _{\stackrel{~}{\nu }}`$ is the relic antineutrino energy density, one obtains the final expression for $`\tau _{\nu \stackrel{~}{\nu }}`$:

$$\tau _{\nu \stackrel{~}{\nu }}(E,z)\frac{G^2}{4\pi }\frac{8}{3}\rho _{\stackrel{~}{\nu }}^0E\underset{0}{\overset{z}{}}(1+z^{})^5\frac{dt}{dz^{}}N_{eff}\left(ET_{\stackrel{~}{\nu }}^0(1+z^{})^2\right)𝑑z^{},$$
(110)
$$N_{eff}=N_{eff}^{NC}+N_{eff}^{CC}.$$
(111)
Here, $`\rho _{\stackrel{~}{\nu }}^0`$ is the present electron antineutrino energy density,

$$\rho _{\stackrel{~}{\nu }}^0=\frac{7}{8}\frac{\pi ^2}{30}(T_{\stackrel{~}{\nu }}^0)^4.$$
(112)
Results of the calculations (with $`h=0.67`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$) of $`\tau _{\nu N}`$ and $`\tau _{\nu \stackrel{~}{\nu }}`$ are shown on Fig.4 for several values of neutrino energy. One can see, that the contribution of $`\nu N`$ channel to the total $`\tau `$ is negligibly small everywhere. For typical neutrino energy $`100\text{ MeV}`$ the absorption is essential beginning from $`z310^6`$. A growth of the absorption factor with $`z`$ is very fast, therefore the redshifts for which the condition (105) is violated are never really important in a calculation of the absorption.

Our calculation of neutrino spectra from evaporating PBHs contains two parameters: a spectral index $`n`$ and a time of an end of the inflation $`t_i`$ (which, by assumption, is a time when density fluctuations develop). We assume that at $`t_i`$ the universe has as a result of the reheating the equilibrium temperature $`T_{RH}`$. The connection of $`T_{RH}`$ and $`t_i`$ is given by the standard model:

$$t_i=0.301g_{}^{1/2}\frac{M_{pl}}{T_{RH}^2}\frac{0.24}{(T_{RH}/1MeV)^2}s$$
(113)
($`g_{}100`$ is the number of the degrees of freedom in the early universe).

On Fig.5 the integrand of the background spectrum integral, considered in Sec.III is shown as a function of redshift $`z`$ for several values of the parameter $`T_{RH}`$. The inclusion of the absorption effects leads to an appearing of the additional factor in this integral:

$$S(E)=e^{\tau (E,z)}F(E,z)d\mathrm{log}_{10}(z+1),$$
(114)
$$\tau (E,z)=\tau _{\nu \stackrel{~}{\nu }}(E,z)+\tau _{\nu N}(E,z).$$
(115)
Each curve in Fig.5 has a strong cut off on some redshift value. This feature is connected with the existence in our model the minimum value of PBH mass. The PBH mass spectrum is steeply falling function of the mass, so the masses near minimum give a largest contribution to the neutrino background. The moment of their evaporation is, approximately,

$$t_{ev}\frac{\left(M_{PBH}^{min}\right)^3}{3\alpha },$$
(116)
and the corresponding redshift is determined by the relation

$$z_{ev}+1(z_{eq}+1)\left(\frac{t_{eq}}{t_{ev}}\right)^{1/2}.$$
(117)
Larger masses evaporate at larger times and smaller redshifts. If, for instance, $`T_{RH}=10^{10}\text{GeV}`$, one has $`t_i=0.2410^{26}\text{s}`$ ; $`M_i=7,510^{10}\text{g}`$ , $`z_{ev}10^7`$. So, at $`T_{RH}=10^{10}GeV`$ the redshifts of order of $`10^7`$ give a largest contribution to the neutrino background (if there is no neutrino absorption; the absorption effects are essential just near $`z10^7`$, and, as one can see from the figure, the actual $`z`$-dependences of the integrand have peak values at $`z(2÷3)10^6`$, depending on a neutrino energy).

The cut off value $`z_{ev}`$ strongly depends on $`T_{RH}`$:

$$z_{ev}\frac{1}{\sqrt{t_{ev}}}\left(M_{BH}^{min}\right)^{3/2}t_i^{3/2}T_{RH}^3.$$
(118)
Again, this estimate does not take into account the absorption effects. The real situation is seen on the figure. In particular, at $`T_{RH}>10^{10}\text{ GeV}`$ this cut off is the same for all $`T_{RH}`$ values, as a result of a strong absorption at $`z310^6`$.

On Figs.6,7 spectra of electron neutrino background are shown separately for several channels of the neutrino production. One can see from Fig.6 that 1) the contribution to the summary background from the direct $`\nu _e`$ emission is, at $`T_{RH}10^9\text{GeV}`$, absolutely dominant for $`E1\text{ MeV}`$ and 2) the relative importance of different channels changes with an increase of $`T_{RH}`$. It is seen also that at $`T_{RH}10^9\text{GeV}`$ the contribution of the quark fragmentation channel is very small at $`E1\text{ MeV}`$ $`(5\%)`$, so the evident underestimation of this channel in our calculation, connected with the neglect of the contribution of heavy quark fragmentations, has (at large $`T_{RH}`$ values) no particular importance.

Situation with the contribution from quark fragmentations is slightly different at smaller $`T_{RH}`$ values (see Fig.7). Due to the flattening out of the curve for direct $`\nu _e`$ emission the relative contribution of the fragmentation channel at $`E1\text{MeV}`$ grows with an decrease of $`T_{RH}`$. The highest relative contribution is at $`T_{RH}10^{8.5}\text{GeV}`$. At this $`T_{RH}`$ the minimum value of mass in PBH mass spectrum ($`M_{BH}^{min}`$) coincides with the threshold of quark evaporations ($`m_{thr}10^{14.1}\text{g}`$, according to Fig.1). At further decrease of $`T_{RH}`$ the relative contribution of quark fragmentations at $`E1\text{MeV}`$ fastly decreases and, at some $`T_{RH}`$ value, it is strictly zero (when, due to correspondingly large $`M_{BH}^{min}`$ value, the age of the universe is not enough for ”warming-up” the PBHs).

It follows from Fig.7 that at neutrino energies larger than $`(5÷10)\text{MeV}`$ the contribution of fragmentation channels is negligibly small ( $`5\%`$). As was mentioned above, the region of smaller neutrino energies is not important for the calculation of spectral index constraints, so some possible underestimation of the quark fragmentation contribution in this region (due to the simplified treatment of them in the present paper) is not sufficient.

Some typical results of calculations of summary electron neutrino background spectra are shown on Figs.8,9. The main features of such spectra have been revealed in previous works : $`E^3`$ dependence above $`100\text{ MeV}`$ and the flattening out of the spectra at lower energies. The form of background spectra depends on the form of an initial PBH mass distribution. In refs. the simple power mass spectrum of PBHs containing no cut off (minimum mass value) were used. The behavior of background spectra at low energies in this case is, aside from jet fragmentation contribution, $`E^1`$ and turnover energy is $`100\text{MeV}`$ . In the case of PBH mass spectrum with the cut off the turnover energy depends on $`T_{RH}`$  (it decreases as $`T_{RH}`$ increases, as is seen at Fig.8). According to Figs.8,9 the turnover energy is about $`1\text{ MeV}`$ at $`T_{RH}=10^9\text{ GeV}`$. Besides, it was shown in ref. that at low energies the photon background spectra are almost flat and,as is seen in Figs.8,9, the same is also true in the case of the neutrino background spectra.

## VI Constraints on the spectral index

For an obtaining of the constraints on the spectral index we use three types of neutrino experiments.

1. Radiochemical experiments for the detection of solar neutrinos. There are data on solar neutrino fluxes from the famous Davis experiment  and the $`GaGe`$ experiment . The basic process used in radiochemical experiments is a charged-current neutrino capture reaction

$$\nu _e+A_Z^Ne^{}+A_{Z+1}^{N1}.$$
(119)
In general, the cross section for the neutrino absorption via a bound-bound transition can be calculated using the approximate formula

$$\sigma =\{\begin{array}{cccc}& \frac{G_F^2}{\pi }\left(1^2+\left(\frac{g_a}{g_v}\right)^2\sigma ^2\right)p_eE_e\hfill & \hfill ,& E<100MeV\hfill \\ & \text{const. }\hfill & \hfill ,& E>100MeV.\hfill \end{array}$$
(120)
Here, $`p_e`$ and $`E_e`$ are the momentum and the energy of an electron produced in a neutrino capture reaction. Neutrino and electron energies are connected by the relation

$$E_e=EE_{thr}+m_e.$$
(121)
Stopping of the rise of a neutrino capture cross section (for a bound-bound transition) with neutrino energy is due to influence of a weak nuclear formfactor of the transition. The approximate value of the turnover energy in Eq.(120) is determined from the relations

$$Eq_{max}\frac{1}{R_A}m_\pi .$$
(122)
Here, $`q_{max}`$ is a maximum momentum transfer from a neutrino to the nucleus (in a model with step-like nuclear formfactor), $`R_A`$ is a size of the nucleus.

In the case of the Davis experiment the neutrino capture reaction is

$$\nu +^{37}Cl^{37}Ar+e^{}.$$
(123)
Here we took into account the super- allowed transition only, i.e., the transition to the isotopic analog state in $`{}_{}{}^{37}Ar`$, for which

$$1^2=3,\sigma ^2=0.2,E_{thr}5MeV.$$
(124)
This transition gives dominant contribution to the effect in the detector, in spite of the rather high energy threshold, due to the large value of the Fermi strength $`1^2`$.

In the case of the $`GaGe`$ reaction,

$$\nu +^{71}Ga^{71}Ge+e^{},$$
(125)
the main contribution gives the ground state - ground state transition ($`E_{thr}=0.242MeV`$). In $`{}_{}{}^{71}Ge`$ also there is an isotopic analog state, but its energy is too high, $`9\text{ MeV}`$ (it is slightly higher than the particle emission threshold, so this transition is irrelevant). The cross section of the ground state - ground state transition in the reaction (125) is given by the expression  (for $`E<100\text{ MeV}`$)

$$\sigma \frac{G^2}{\pi }p_eE_e\frac{2J_A+1}{2J_B+1}\frac{6163.4}{ft}.$$
(126)
Here, $`J_A`$ and $`J_B`$ are angular momenta of $`{}_{}{}^{71}Ga`$ and $`{}_{}{}^{71}Ge`$, respectively, $`ft`$ value is determined by the half-life time of $`{}_{}{}^{71}Ge`$. According to ref. one has

$$J_A=\frac{3}{2},J_B=\frac{1}{2},\mathrm{lg}ft=4.37.$$
(127)
The number of target atoms $`N^T`$ in the detectors is about $`2.210^{30}`$ for the $`ClAr`$ experiment  and $`10^{29}`$ for the $`GaGe`$ experiments . The average statistics is $`1.5\text{ captures/day}`$ ($`ClAr`$) and $`1\text{ capture/day}`$ ($`GaGe`$). The constraints are calculated using the relation

$$4\pi N^T10^5S(E)\sigma (E)𝑑E<1.$$
(128)
2. The experiment on a search of an antineutrino flux from the Sun . In some theoretical schemes (e.g., in the model of a spin - flavor precession in a magnetic field) the Sun can emit rather large flux of antineutrinos. LSD experiment  sets the upper limit on this flux, $`\mathrm{\Phi }_{\stackrel{~}{\nu }}/\mathrm{\Phi }_\nu 1.7\%`$. In this experiment the neutrino detection is carried out using the reaction

$$\stackrel{~}{\nu }_e+pn+e^+.$$
(129)
The number of target protons is $`8.610^{28}`$ per 1 ton of the scintillation detector, and the obtained upper limit is $`0.28`$ antineutrino events per year per ton . The cross section of the reaction (129) is well known (see, e.g.,). In particular, in the low energy region ($`Em_p`$) it is given by the formula

$$\sigma (\stackrel{~}{\nu }+p+e^++n)=\frac{G_F^2}{\pi }\left(1+3\left(\frac{g_a}{g_v}\right)^2\right)p_eE_e,$$
(130)
$$E_e=E(m_nm_p).$$
(131)
The threshold of this reaction is given by the relation

$$E_{thr}=m_nm_p+m_e1.8\text{ MeV}.$$
(132)
The cross section of the reaction (129) grows with the neutrino energy up to $`E2GeV`$, and is a constant ($`0.510^{38}cm^2`$) at larger energies.The product of this cross section and a PBH antineutrino background spectrum has the broad maximum at $`E100MeV`$. The constraint is determined from the condition that the calculated effect in the LSD detector is smaller than the upper limit obtained in ref..

3. The Kamiokande experiment on a detection of atmospheric electron neutrinos . In this experiment the electrons arising in the reaction

$$\nu _e^{atm}+np+e^{}$$
(133)
in the large water Cherenkov detector were detected and, moreover, their energy spectrum was measured. This spectrum has a maximum at the energy about $`300MeV`$. The spectrum of the atmospheric electron neutrinos is calculated(see, e.g.,) with a very large accuracy (assuming an absence of the neutrino oscillations) and the experimentally measured electron spectrum coincides, more or less, with the theoretical prediction. The observed electron excess at $`E100MeV`$ (which is a possible consequence of the oscillations) is not too large. We use the following condition for an obtaining the our constraint : the absolute differential intensity of the PBH neutrino background at the neutrino energy $`E0.3GeV`$ cannot exceed the theoretical differential intensity of the atmospheric electron neutrinos at the same energy (otherwise the total number of electrons in the detector and the energy spectrum of these electrons must be different from the observed ones).

## VII Results and discussions

Fig.10,11 shows our results for the spectral index constraints. It is seen that the best constraints are obtained using the Kamiokande atmospheric neutrino data and the LSD upper limit on an antineutrino flux from the Sun.

It appears that the spectral index constraints obtained in the present work are not stronger than the corresponding constraint following from COBE data (see Eq.(79)). One should remember, however, that our approach directly constrains the perturbation amplitude $`\sigma _H(M_h)`$ rather than the spectral index. Comoving length scales of perturbations, for which we constrain the amplitude, are very small in comparison with the size of present horizon. For example, if $`T_{RH}310^8\text{ GeV}`$ the minimum PBH mass value is about $`10^{14}\text{ g}`$. In this case the time of PBH formation is

$$t_i\frac{M_{BH}^{min}}{M_{pl}}t_{pl}10^{24}\text{ s},$$
(134)
the corresponding scale factor is, according to Eq.(43), about$`10^{22}`$, and for the comoving length scale one obtains

$$\lambda _i\frac{ct_i}{a}10^9\text{ pc}.$$
(135)
It is clear that the normalization on COBE data has sense only if it is assumed that the spectral index is independent on scale. The curves of Fig.10,11 were obtained using just this assumption; it is very easy to calculate from them the corresponding constraints on $`\sigma _H(M_h)`$.

For a demonstration of the sensitivity of these results to a chosen $`\mathrm{\Omega }_\mathrm{\Lambda }`$-value we calculate the constraints for two extreme $`\mathrm{\Omega }_\mathrm{\Lambda }`$-values, for which the parametrization for $`f_S^{(0)}`$ (Eq.(76)) is still valid. It is clear that some weakening of the constraint with an increase of $`\mathrm{\Omega }_\mathrm{\Lambda }`$ is due to the corresponding decrease of values of the function $`\sigma _H(M_h)`$ entering the initial PBH mass spectrum (see Eq.(20)). This decrease of $`\sigma _H(M_h)`$ takes place solely due to the dependence of $`f_S^{(0)}`$ on $`\mathrm{\Omega }_\mathrm{\Lambda }`$ (see Eq.(76)).

The behavior of the constraint curve $`n(T_{RH})`$ is sharply different from that was obtained in ref.  (authors of ref.  used the initial PBH mass function following from the near critical collapse scenario with a domination of the earliest epoch of PBHs formation, and diffuse extragalactic photon data).

The slight bend at $`T_{RH}10^{10}\text{GeV}`$ of the constraint curves (Fig.10) is an effect of a neutrino absorption in the space. Constraint based on the atmospheric neutrino experiment is especially sensitive to the absorption because neutrinos of larger energies (in average) are responsible for it.

Flattening out of the constraint curve $`n(T_{RH})`$ with increase of $`T_{RH}`$ is connected with the dependence $`M_{BH}^{min}(T_{RH})`$. According to the relations (34,52,113), the minimum value of PBH mass in the initial mass spectrum is inversely proportional to $`T_{RH}^2`$. More exactly, the relation between $`M_{BH}^{min}`$ and $`T_{RH}`$ is

$$M_{BH}^{min}(g)=\frac{7.210^{30}}{(T_{RH}/1GeV)^2}.$$
(136)
At small $`T_{RH}`$ the minimum value of PBH mass in the initial PBH mass spectrum is large and the neutrino background is dominated by the evaporations at recent epochs. In opposite, at large $`T_{RH}`$, when $`M_{BH}^{min}`$ is small, the neutrino background is dominated by the evaporations at earlier times. Thus, using the relation

$$z_{eq}2.410^4\mathrm{\Omega }_mh^2$$
(137)
and the curves, corresponding to $`T_{RH}=10^9\text{GeV}`$ in Fig.5, one can see that in cases when $`T_{RH}>10^9\text{GeV}`$ the neutrino background is dominated by PBH evaporations during the radiation era.

Introducing, as usual , the function $`M^{}(M_{BH}^{min})`$ by the formula

$$M^{}(M_{BH}^{min})=\left[3\alpha (M_{BH}^{min})t_0\right]^{1/3},$$
(138)
where $`t_0`$ is the age of the universe, one can easily see that at low $`T_{RH}`$ (before the flattening out of our $`n(T_{RH})`$ curve) the inequality $`M_{BH}^{min}>M^{}(M_{BH}^{min})`$ takes place while at large $`T_{RH}`$ one has $`M_{BH}^{min}<M^{}(M_{BH}^{min})`$. We remind that $`M^{}`$ is the mass of a black hole, formed in the early universe and evaporating today. The turning-point (i.e., the point of a change of the slope of the $`n(T_{RH})`$ curve) corresponds to the equality $`M_{BH}^{min}=M^{}(M_{BH}^{min})`$ and in the case $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$ it is at $`T_{RH}10^8\text{GeV}`$. The shift of this turning-point with an increase of $`\mathrm{\Omega }_\mathrm{\Lambda }`$ is explained by the corresponding increase of an age of the universe:

$$t_0=\frac{2}{3H_0\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }}}\mathrm{ln}\frac{1+\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }}}{\sqrt{1\mathrm{\Omega }_\mathrm{\Lambda }}}\stackrel{\mathrm{\Omega }_\mathrm{\Lambda }0}{}\frac{2}{3H_0}.$$
(139)
One can compare our spectral index constraints with the corresponding results of ref. , where the same initial mass spectrum of PBHs had been used. Some difference in resulting constraints ($`\mathrm{\Delta }n0.02`$) may be connected simply with the fact that authors of ref.  use slightly different formula for $`\sigma _H(M_h)`$, namely,

$$\sigma _H(M_h)=\sigma _H(M_{h,0})\left(\frac{M_h}{M_{h,0}}\right)^{\frac{1n}{4}}.$$
(140)
In other respects the constraints are quite similar although in ref.  they were based on diffuse extragalactic photon background data.

One should note that usually the calculations of these constraints are accompanied by the calculation of the bounds based on requirement that the energy density in PBHs does not overclose the universe at any epoch ($`\mathrm{\Omega }_{BH}<1`$). For a setting of such bounds one must consider the cosmological evolution of the system PBHs + radiation. We intend to carry out these calculations in a separate paper. Estimates show that at large $`T_{RH}`$ values ($`T_{RH}10^{10}\text{GeV}`$) the constraints based on $`\mathrm{\Omega }_{BH}<1`$ are stronger than those based on the neutrino experiments (it is the reason why we did not explore the region $`T_{RH}>10^{10}\text{GeV}`$ in the present paper).

The neutrino background spectra obtained in the present paper may be useful somewhere else, irrespective of the spectral index constraint problem. The main feature of such spectra at $`E1MeV`$ and at $`T_{RH}10^910^{10}\text{ GeV}`$ is the very small contribution ($`5\%`$) of all secondary channels of the neutrino production at evaporation (except the decay of directly evaporated muons) (see Fig.6). There are two reasons for this: the steepness of the initial PBH mass spectrum and, most essentially, the existence in our model the cut off in this spectrum, i.e., the minimum PBH mass value. Evaporation of PBHs having masses near $`M_{BH}^{min}`$ gives the dominant contribution in the background, so there is a strong maximum in $`z`$ distributions at low energies (Fig.5). For example, at $`T_{RH}=10^{10}\text{GeV}`$ the peak value of $`z`$ is $`10^6`$. Neutrino energy $`1\text{MeV}`$ corresponds in this case to energy at evaporation $`1\text{TeV}`$. For an evaporation of such high energy neutrinos the PBH mass must be $`10^{10}\text{g}`$ in the case of direct neutrinos and much smaller in the case of secondary neutrinos (because their energies are much lower in average (see Fig.2)). At the same time, the minimum value of PBH mass in the mass spectrum at $`T_{RH}=10^{10}\text{GeV}`$ is $`10^{11}\text{g}`$. So, it should be clear that, in our example, the secondary neutrinos with $`E1\text{MeV}`$ can appear only after the essential evolution of the initial PBH mass spectrum (see Eq.(37) and the detailed discussion of this evolution in ref.). This evolution leads to an increase of the total mass interval and to the corresponding decrease of a PBH number density (per unit mass interval) and, as a result, to a relative smallness of the secondary neutrino flux.

One can compare our results for the neutrino background spectra with those of ref.. According to the conclusions of ref., the secondary neutrino fluxes are larger relatively and even dominant at $`E10\text{MeV}`$. Similar results were obtained in earlier work . In our model the dominance of secondary neutrino fluxes at $`E<10\text{ MeV}`$ takes place at some values of $`T_{RH}`$ only (see Fig.7).

The results for neutrino diffuse background spectra and the constraints on the spectral index obtained here do not take into account the possible existence of the photospheres  around a black hole. In studies of photon background spectra two kinds of photospheres should be considered: QCD photosphere (quark-gluon plasma arising from interactions of quarks and gluons before hadronization) and QED photosphere ($`e^+e^{}\gamma `$ plasma, in which directly evaporated photons are processed). As for the latter, in our case one should consider instead the ”electroweak (or neutrino) photosphere”  which may form only at extremely high black hole temperatures and, probably, has no influence on the background spectra. QCD photosphere may be more important for the problem. Thus, the recent calculations  of instantaneous photon spectra from the model QCD photosphere show that some increase of photon fluxes (in comparison with the usual case when QCD photosphere does not form and photons appear simply through the direct quark fragmentations with subsequent $`\pi ^0`$-decays) definitely takes place. For example, according to these calculations, at $`T_{BH}=50\text{GeV}`$ the peak value of instantaneous photon spectrum increases on a factor $`4`$. If it is true and if this increase is the same at all black hole temperatures, the contribution to neutrino background spectra from the quark fragmentations obtained in our calculations must be increased approximately on the same factor.

###### Acknowledgements.

We wish to thank G. V. Domogatsky for valuable discussions and comments. We are grateful also to H.I.Kim for informing us about his work on the same problem and for useful remarks.