format: text

text:
# 1 Introduction

## 1 Introduction

When he conjectured the existence of the neutrino, Pauli imposed very stringent bounds on its electrical neutrality.

Nevertheless, the direct (=weak) interaction of neutrino with electrons and positrons shows a principle possibility of the interaction of neutrino with the electromagnetic field. In vacuum, however, this interaction is negligible. The situation is extremely changed in media with free carriers of the electric charge such as a dense plasma of metals, stars, the lepton plasma of the early Universe, etc.

It was shown in  that the electrodynamics of neutrinos is changed in such media comparing with the electrodynamics of neutrinos in vacuum in such a way that the contribution of the neutrino electromagnetic vertex in medium cannot be regarded as a small correction like the case of radiative corrections in vacuum. The main distinction is the appearance of the induced electric charge of neutrino which is proportional to the density of free carriers of electricity 111There appears also induced magnetic moment of neutrino , cross-sections of the neutrino scattering off nuclei are modified in plasmas especially in low energy region  and so on.. Namely this leads to the new qualitative effect: the appearance of the long-distance forces inevitably leads to the collective interactions of neutrino with an ensemble of charged particles via electromagnetic field.

The approximation of the neutrino propagation in an isotropic plasma considered in  and then in  has some natural bounds for applications in astrophysics. Therefore one needs to have the self-consistent system of kinetic equations of the most general kind applicable for other astrophysical plasmas including magnetized stars.

Really this turns out to be a new problem unsolved before in the neutrino kinetics while in the one-particle Schrödinger equation approach one can show the importance of the axial vector potential for the neutrino propagation in a magnetized plasma ($`V_A<\overline{\psi }_e\gamma _i\gamma _5\psi _e>B_i`$) that provides the interaction of neutrino with the magnetic field B without any neutrino magnetic moment  and could explain, e.g. the neutron star kick .

In the present work we take into account the spin interactions of electrons that are important for neutrino propagation in a dense degenerate electron gas polarized by the external magnetic field.

In subsection 2.1. we discuss the equilibrium state of an electron gas polarized by the external magnetic field $`𝐁_00`$.

Then in subsection 2.2 we present the full set of the coupled Relativistic Kinetic Equations (RKE) for electrons and neutrinos in the lepton plasma with an external magnetic field including the self-consistent electromagnetic field contribution. These RKE’s are derived using the Bogolyubov method analogously to the way described in details for an isotropic plasma in .

In subsection 2.3 we complete derivation of RKE’s recasting them for the gauge invariant distribution functions and prove the lepton current conservation,
$`j_\mu ^{(a)}(𝐱,t)/x_\mu =0`$, in a magnetized plasma both for neutrinos and electrons.

We conclude this section deriving in subsection 2.4 the dispersion equation for the lepton density perturbations in an isotropic plasma that coincides with results .

In section 3 we generalize the Bloch equation for the evolution of the magnetization density perturbation in the presence of neutrino fluxes. In section 4, neglecting spatial dispersion we obtain the increment of spin waves excited by the neutrino beam and compare it with the case of the neutrino driven instability of plasma waves in an isotropic plasma.

In section 5 we discuss results and give some estimates of relevant parameters in a polarized medium.

## 2 Relativistic Kinetic Equations for lepton plasma in Vlasov approximation

The kinetic equations in the Standard Model (SM) are derived from the quantum Liouville equation for the non-equilibrium statistical operator $`\widehat{\rho }(t)`$ using the Bogolyubov method with the total interaction Hamiltonian given by the Feynman diagrams for the neutrino scattering off electrons and the usual ee-scattering in QED . For a magnetized plasma we do not consider $`\nu \nu `$-scattering and neglect also weak $`ee`$-scattering comparing with the electromagnetic interaction of charged particles through the photon exchange.

The one-particle density matrix $`\widehat{f}_{r^{}r}^{(a)}(𝐩,𝐱,t)=_𝐤e^{i\mathrm{𝐤𝐱}}Tr(\widehat{\rho }(t)\widehat{b}_{𝐩𝐤/2,r}^{+(a)}\widehat{b}_{𝐩+𝐤/2,r^{}}^{(a)})`$ is determined as the statistical average of the product of creation ($`\widehat{b}_{𝐩r}^{+(a)}`$) and annihilation ($`\widehat{b}_{𝐩^{}r^{}}^{(a)}`$) operators for the lepton $`a=e,\nu `$. In the case of electrons it takes the form

$$\widehat{f}_{\lambda ^{}\lambda }^{(e)}(𝐩_e,𝐱,t)=f^{(e)}(𝐩_e,𝐱,t)\frac{\delta _{\lambda ^{}\lambda }}{2}+S_i^{(e)}(𝐩_e,𝐱,t)\frac{(\sigma _i)_{\lambda ^{}\lambda }}{2},$$

where $`\lambda =\pm 1`$ is the eigenvalue of the spin projection on an external magnetic field. E.g. for the magnetic field $`𝐁=(0,0,B_0)`$ only the spin component $`\mathrm{\Sigma }_z`$ commutes with the Hamiltonian for electrons, $`[\mathrm{\Sigma }_z,\widehat{}^{(e)}]_=0`$, hence $`(\sigma _z)_{\lambda ^{}\lambda }=\lambda \delta _{\lambda ^{}\lambda }`$. In the absence of neutrino magnetic moment the one-particle density matrix for massless Dirac neutrinos corresponds to the active (left-handed, $`r=1`$) neutrino distribution function only,

$$\widehat{f}_{r^{}r}^{(\nu )}(𝐪,𝐱,t)=\frac{(1r)\delta _{r^{}r}}{2}f^{(\nu )}(𝐪,𝐱,t),$$

where in an uniform medium the helicity $`r=\pm 1`$ conserves, $`(\sigma _iq_i)_{r^{}r}=qr\delta _{r^{}r}`$.

The number density distribution functions are related with the lepton densities $`n^e(𝐱,t)=(d^3p/(2\pi )^3)f^{(e)}(𝐩,𝐱,t)`$ for electrons and $`n^\nu (𝐱,t)=(d^3q/(2\pi )^3)f^{(e)}(𝐪,𝐱,t)`$ for neutrinos.

### 2.1 Equilibrium state in a polarized electron gas

In the linear approximation for a slightly inhomogeneous medium the distribution functions have the form

$`f^{(\nu )}(𝐪,𝐱,t)=f_0^{(\nu )}(𝐪)+\delta f^{(\nu )}(𝐪,𝐱,t),`$
$`f^{(e)}(𝐩,𝐱,t)=f_0^{(e)}(\epsilon _p)+\delta f^{(e)}(𝐩,𝐱,t),`$
$`𝐒^{(e)}(𝐩,𝐱,t)=𝐒_0^{(e)}(\epsilon _p)+\delta 𝐒^{(e)}(𝐩,𝐱,t),`$ (1)
where the neutrino background is given by $`f_0^{(\nu )}(𝐪)`$ not being in equilibrium with plasma environment (e.g. for the powerful neutrino flux outside of the SN neutrinosphere), and we consider the uniform polarized equilibrium electron gas for which the one-particle Wigner density matrix takes the form

$$\widehat{f}_{\lambda ^{}\lambda }^{(0e)}(\epsilon _p)=\frac{\delta _{\lambda ^{}\lambda }}{2}f^{(0)}(\epsilon _p)+\frac{(\sigma _j)_{\lambda ^{}\lambda }}{2}S_j^{(0)}(\epsilon _p).$$
(2)
Here $`f^{(0)}(\epsilon _p)=g_e[\mathrm{exp}(\epsilon _p\zeta )/T)+1]^1`$ is the equilibrium Fermi function; $`g_e=2`$ is the Lande factor;

$$S_j^{(0)}(\epsilon _p)=\frac{\mu _BB_j^{(0)}}{\gamma }\frac{df^{(0)}(\epsilon _p)}{d\epsilon _p}$$
(3)
is the equilibrium spin distribution; $`\gamma =\epsilon _p/m_e`$ is the electron gamma factor, and we use hereafter units $`\mathrm{}=c=1`$.

Note that in (3) we assumed the quasi-classical electron spectrum in a realistic external magnetic field $`B_0`$ obeying the inequality $`eB_0T^2`$, when the Landau spectrum

$$\epsilon (\lambda ,p_z,n)=(m_e^2+p_z^2+eB_0(2n+1\lambda ))^{1/2}$$
(4)
reduces to

$$\epsilon _p(\lambda )=\epsilon _p\frac{\lambda eB_0}{2\epsilon _p},$$
(5)
where $`\lambda =\pm 1`$ and we changed $`eB_0(2n+1)=p_{}^2`$ for large Landau numbers $`n1`$. Thus, the full spectrum (5) contains the small paramagnetic (spin) correction to the continuous spectrum $`\epsilon _p=\sqrt{p^2+m_e^2}.`$

One can easily check that in this quasi-classical limit the exact expression for magnetization of the electron gas,

$`M_j^{(0)}=`$ $`\mu _B<\overline{\psi }_e\gamma _j\gamma _5\psi >_0=\mu _B{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{eB_0}{(2\pi )^2}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑p_zTr[\sigma _j\widehat{f}^{(0e)}(\epsilon (\lambda ,p_z,n))]=`$ (6)
$`=\mu _B{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{eB_0\delta _{jz}}{(2\pi )^2}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑p_z{\displaystyle \underset{\lambda }{}}\lambda f_{\lambda \lambda }^{(0e)}(\epsilon (\lambda ,p_z,n)),`$
takes the quasi-classical form

$$M_j^{(0)}=\mu _B\frac{d^3p}{(2\pi )^3}S_j^{(0)}(\epsilon _p)=2\mu _B^2B_{0j}D(\epsilon _p)\frac{df^{(0)}(\epsilon _p)}{d\epsilon _p}𝑑\epsilon _p,$$
(7)
where the spin distribution $`S_j^{(0)}(\epsilon _p)`$ is given by eq. (3) and $`D(\epsilon _p)=pm_e/(2\pi )^2`$. In particular, in the non-relativistic (NR) limit this background magnetization corresponds to the spin paramagnetism of the free electron gas in metal .

Note that for the degenerate electron gas from the quasi-classical Eq. (7) one obtains the same value $`M_z^{(0)}=\mu _BeB_0p_{F_e}/2\pi ^2`$ as it follows from the exact quantum eq. (6) for an arbitrary strong magnetic field for which electrons populating the main Landau level (n=0) contribute only.

Note also that this magnetization determines the axial vector potential of a probing neutrino in a magnetized plasma too ,

$$V_A=G_F\sqrt{2}c_A\frac{q_j}{q}\frac{M_j^{(0)}}{\mu _B},$$

that changes the spectrum of the ultrarelativistic neutrino, $`\epsilon _q=q+V+V_A`$, comparing with the standard one in an isotropic medium (e.g. in the Sun), $`\epsilon _q=q+V`$, and modifies the neutrino oscillations in a magnetized SN described by the one-particle Schrödinger equation for two neutrino species ,

$$i\left(\begin{array}{c}\dot{\nu }_a\hfill \\ \dot{\nu }_b\hfill \end{array}\right)=\left(\begin{array}{cc}V+V_Ac_2\delta & s_2\delta \\ s_2\delta & 0\end{array}\right)\left(\begin{array}{c}\nu _a\\ \nu _b\end{array}\right),$$
(8)
where $`c_2=\mathrm{cos}2\theta `$, $`s_2=\mathrm{sin}2\theta `$, $`\delta =\mathrm{\Delta }m^2/4E`$ are the neutrino mixing parameters; $`\nu _b`$ is the sterile neutrino wave function.

In our kinetic approach for the SM lepton plasma we do not consider oscillations ($`s_2=0`$) assuming the massless spectrum $`\epsilon _q=q`$ for active neutrino of the given kind $`\nu _a`$, $`a=e,\mu ,\tau `$.

### 2.2 The master kinetic equations for perturbations in plasma

Within the linear approximation (2.1) the RKE for the Lorentz-invariant Wigner number density distribution functions $`\delta f^{(a)}(𝐩,𝐱,t)=Tr\left[\widehat{\delta f}^{(a)}(𝐩,𝐱,t)\right]`$ take the covariant forms

$`q_\mu {\displaystyle \frac{\delta f^{(\nu )}(𝐪,𝐱,t)}{x_\mu }}+G_F\sqrt{2}c_V{\displaystyle \frac{d^3p_e}{(2\pi )^3}\frac{(p_e^\mu q_\mu )}{\epsilon _{p_e}}\frac{\delta f^{(e)}(𝐩_e,𝐱,t)}{x_j}\frac{f_0^{(\nu )}(𝐪)}{q_j}}+`$
$`+{\displaystyle \frac{G_Fc_Am_eq^\mu }{\sqrt{2}}}{\displaystyle \frac{d^3p_e}{(2\pi )^3\epsilon _{p_e}}\frac{\delta a_\mu ^{(e)}(𝐩_e,𝐱,t)}{x_j}\frac{f_0^{(\nu )}(𝐪)}{q_j}}=0,`$ (9)
for neutrinos and

$`p_\mu {\displaystyle \frac{\delta f^{(e)}(𝐩,𝐱,t)}{x_\mu }}+e\delta F_{j\mu }(𝐱,t)p^\mu {\displaystyle \frac{f_0^{(e)}(\epsilon _p)}{p_j}}+eF_{j\mu }^{(0)}p^\mu {\displaystyle \frac{\delta f^{(e)}(𝐩,𝐱,t)}{p_j}}+`$ (10)
$`+G_Fc_V\sqrt{2}{\displaystyle \frac{d^3q}{(2\pi )^3}\frac{(p^\mu q_\mu )}{\epsilon _q}\frac{\delta f^{(\nu )}(𝐪,𝐱,t)}{x_j}\frac{f_0^{(e)}(\epsilon _p)}{p_j}}+`$
$`+G_Fc_A\sqrt{2}m_e{\displaystyle \frac{d^3q}{(2\pi )^3}\frac{(q_\mu a_i^\mu (𝐩))}{\epsilon _q}\frac{S_i^{(0e)}(\epsilon _p)}{p_j}\frac{\delta f^{(\nu )}(𝐪,𝐱,t)}{x_j}}=0,`$
for electrons where $`c_V=2\xi \pm 0.5`$, $`c_A=0.5`$ are the vector and axial weak couplings with upper (lower) sign for electron (muon or tau) neutrinos correspondingly, $`\xi =\mathrm{sin}^2\theta _W`$ is the Weinberg parameter; the tensor $`F_{jk}^{(0)}=e_{jkl}B_{0l}`$ corresponds to the external magnetic field term.

Here we have introduced in eq. (2.2) the 4-vector $`a_\mu ^{(e)}(𝐩,𝐱,t)a_\mu ^i(𝐩)S^i(𝐩,𝐱,t)`$ that is the statistical generalization  of the Pauli-Lubański 4-vector

$$a^\mu (𝐩)=(a_i^\mu (𝐩)\zeta _i)=[\frac{𝐩\stackrel{}{\zeta }}{m_e};\stackrel{}{\zeta }+\frac{𝐩(𝐩\stackrel{}{\zeta })}{m_e(\epsilon _p+m_e)}]$$

and has the components

$$a_\mu ^{(e)}(𝐩,𝐱,t)=[\frac{\mathrm{𝐩𝐒}^{(e)}(𝐩,𝐱,t)}{m_e};𝐒^{(e)}(𝐩,𝐱,t)+\frac{𝐩(\mathrm{𝐩𝐒}^{(e)}(𝐩,𝐱,t))}{m_e(\epsilon _e+m_e)}].$$
(11)
Note that in the presence of the external magnetic field $`𝐁_0`$ the neutrino parts of these Vlasov equations differ from the analogous ones used in  due to the inclusion of axial vector interactions ($`c_A`$) in the total $`\nu `$e-scattering amplitude that allows us to account for collective interactions of $`\nu _\mu `$ and $`\nu _\tau `$-neutrinos too for which the vector coupling $`c_V`$ is small, $`c_V=2\xi 0.50`$.

To check general symmetries we temporarily refuse the linearization Eq. (2.1) and write down the RKE for the total electron spin distribution function $`𝐒^{(e)}(𝐩,𝐱,t)`$ accounting for the weak interaction of electrons with neutrinos:

$`{\displaystyle \frac{S_i^{(e)}(𝐩,𝐱,t)}{t}}+𝐯{\displaystyle \frac{S_i^{(e)}(𝐩,𝐱,t)}{𝐱}}+eF_{j\mu }(𝐱,t){\displaystyle \frac{p^\mu }{\epsilon _p}}{\displaystyle \frac{S_i^{(e)}(𝐩,𝐱,t)}{p_j}}+`$
$`+2\mu _B\left[{\displaystyle \frac{E_i(𝐱,t)(\mathrm{𝐯𝐒}^{(e)}(𝐩,𝐱,t))v_i(𝐄(𝐱,t)𝐒^{(e)}(𝐩,𝐱,t))}{1+\gamma }}+{\displaystyle \frac{[𝐒^{(e)}(𝐩,𝐱,t)\times 𝐁(𝐱,t)]}{\gamma }}\right]+`$
$`+{\displaystyle \frac{2G_F\sqrt{2}c_Am_e}{\epsilon _p}}e_{ikl}{\displaystyle \frac{d^3q}{(2\pi )^3}\frac{(q_\mu a_k^\mu (𝐩))}{\epsilon _q}f^{(\nu )}(𝐪,𝐱,t)S_l^{(e)}(𝐩,𝐱,t)}+`$
$`+{\displaystyle \frac{G_F\sqrt{2}c_Am_e}{\epsilon _p}}{\displaystyle \frac{d^3q}{(2\pi )^3}\frac{(q_\mu a_i^\mu (𝐩))}{\epsilon _q}\frac{f^{(e)}(𝐩,𝐱,t)}{p_j}\frac{f^{(\nu )}(𝐪,𝐱,t)}{x_j}}+`$
$`+{\displaystyle \frac{G_F\sqrt{2}c_V}{\epsilon _p}}{\displaystyle \frac{d^3q}{(2\pi )^3}\frac{(p^\mu q_\mu )}{\epsilon _q}\frac{S_i^{(e)}(𝐩,𝐱,t)}{p_j}\frac{f^{(\nu )}(𝐪,𝐱,t)}{x_j}}+`$
$`+{\displaystyle \frac{G_F\sqrt{2}c_V}{\epsilon _pm_e}}[e_{kmn}{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}{\displaystyle \frac{p_mq_k}{\epsilon _q}}\left(v_s{\displaystyle \frac{f^{(\nu )}(𝐪,𝐱,t)}{x_s}}\right)`$
$`e_{\mu m\rho \sigma }{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}{\displaystyle \frac{q^\mu p^\rho a_n^\sigma (𝐩)}{\epsilon _q}}{\displaystyle \frac{f^{(\nu )}(𝐪,𝐱,t)}{x_m}}]e_{ijn}S_j^{(e)}(𝐩,𝐱,t)=0.`$ (12)
This is the relativistic generalization in SM of the kinetic equation for spin waves in nonferromagnetic metals . But it should be kept in mind that we consider here the Fermi gas of free electrons in contrast with the quasi-particle approach for metals and also neglected exchange interactions (neglecting exchange Feynman diagrams for the ee-scattering means here that the long-range forces are dominant).

The system of RKE’s is completed by the Maxwell equations for the electromagnetic field $`F_{\mu \nu }(𝐱,t)=\delta F_{\mu \nu }(𝐱,t)+F_{\mu \nu }^{(0)}`$ accounting for the spin wave contribution,

$$\frac{1}{c}\frac{𝐁(𝐱,t)}{t}=[\times 𝐄(𝐱,t)],$$

and

$`[\times [\times 𝐄(𝐱,t)]]+{\displaystyle \frac{1}{c^2}}{\displaystyle \frac{^2𝐄(𝐱,t)}{t^2}}={\displaystyle \frac{4\pi }{c^2}}{\displaystyle \frac{}{t}}(e{\displaystyle }{\displaystyle \frac{d^3p}{(2\pi )^3}}𝐯f^{(e)}(𝐩,𝐱,t)+`$
$`+\mu _B[{\displaystyle }{\displaystyle \frac{d^3p}{(2\pi )^3}}{\displaystyle \frac{[\times 𝐒^{(e)}(𝐩,𝐱,t)]}{\gamma }}{\displaystyle }{\displaystyle \frac{d^3p}{(2\pi )^3}}(𝐯){\displaystyle \frac{[𝐯\times 𝐒^{(e)}(𝐩,𝐱,t)]}{\gamma +1}}]).`$ (13)
In eqs. (2.2)–(2.2) $`e=e`$ is the electric charge of the electron; $`G_F`$ is the Fermi constant; $`\mu _B=e\mathrm{}/2m_ec`$ is the Bohr magneton; $`𝐯=𝐩/\epsilon _p`$ is the electron velocity; the latin indices run $`i,j=1,2,3`$ and the greek ones are $`\mu =0,1,2,3`$ for scalars written in the Feynman metrics, $`A_\mu B^\mu =A_0B_0A_iB_i`$.

In the non-relativistic (NR) limit $`𝐯1`$ ($`\gamma 1`$) the last term drops out and eq. (2.2) coincides with the Maxwell equation in a magnetized medium:

$$[\times (𝐁4\pi 𝐌)]=4\pi 𝐣+\frac{𝐄}{t},$$
(14)
where $`𝐌(𝐱,t)=\mu _Bd^3p𝐒(𝐩,𝐱,t)/(2\pi )^3`$ is the magnetization density of NR plasma.

Multiplying eq. (2.2) by the energy $`\epsilon _p`$ and the spin distribution function $`S_i(𝐩,𝐱,t)`$ we obtain the covariant RKE for the Lorentz-invariant product $`a_\mu (𝐩,𝐱,t)a^\mu (𝐩,𝐱,t)=(𝐒^{(e)}(𝐩,𝐱,t))^2`$ where the 4-vector $`a_\mu ^{(e)}(𝐩,𝐱)`$ is given by (11).

Note that for the uniform electron beam ($`S_i^{(e)}(𝐩,t)=S_i(t)(2\pi )^3n_{e0}\delta ^{(3)}(𝐩𝐩_0)`$)
omitting the neutrino term and integrating the spin RKE (2.2) over $`d^3p`$ we obtain the one-particle spin evolution equation

$$\frac{d𝐒(t)}{dt}=\frac{2\mu _B}{\gamma _0}\left[𝐒\times 𝐁\right]+\frac{2\mu _B}{\gamma _0+1}\left[𝐒\times \left[𝐄\times 𝐯_0\right]\right],$$

that turns out to be exactly the Bargman-Mishel-Telegdi equation for the electron spin motion with the normal magnetic moment $`\mu _B`$ . Here $`𝐯_0=𝐩/\epsilon _{p_0}`$ is the velocity of the electron beam, $`\gamma _0=\epsilon _{p_0}/m_e`$ is its $`\gamma `$-factor.

In general, it is possible to generalize eq. (2.2) adding the electromagnetic scattering of electrons through anomalous magnetic moment (Schwinger correction $`\mu ^{}=(\alpha /2\pi )\mu _B`$) that could lead to additional terms similar to electromagnetic terms in the neutrino spin RKE .

Below we consider some important particular cases of the master equations checking their consistency with the known results .

### 2.3 Lepton current conservation $`j_\mu ^{(a)}(𝐱,t)/x_\mu =0`$

Let us rewrite eq. (2.2) as the classical equation

$$\frac{\delta f^{(\nu )}(𝐪,𝐱,t)}{t}+\dot{𝐱}\frac{\delta f^{(\nu )}(𝐪,𝐱,t)}{𝐱}+\dot{𝐪}\frac{f_0^{(\nu )}(𝐪)}{𝐪}=0,$$
(15)
where $`\dot{𝐱}=𝐧=𝐪/q`$ is the velocity of massless neutrino, the derivative $`\dot{𝐪}`$ is given by

$`\dot{𝐪}=+`$ $`\left[G_F\sqrt{2}c_V\delta n^{(e)}(𝐱,t)+{\displaystyle \frac{G_F}{\sqrt{2}}}c_A\delta A_0(𝐱,t)\right]`$ (16)
$`G_F\left[\sqrt{2}c_V[𝐧\delta 𝐣^{(e)}(𝐱,t)]+{\displaystyle \frac{c_A}{\sqrt{2}}}[𝐧\delta 𝐀^{(e)}(𝐱,t)]\right].`$
Here $`\delta j_\mu ^{(e)}(𝐱,t)=(\delta n^{(e)}(𝐱,t);\delta 𝐣^{(e)}(𝐱,t))=(d^3p/(2\pi )^3)(p_\mu /\epsilon _p)\delta f^{(e)}(𝐩,𝐱,t)`$ is the four-vector of the electron current density perturbation ;
$`\delta A_\mu ^{(e)}(𝐱,t)=m_e(d^3p/(2\pi )^3)\delta a_\mu ^{(e)}(𝐩,𝐱,t)/\epsilon _p`$ is the axial four-vector of the spin density perturbation222In NR plasma the spin density is given by the 3-vector component only, $`\delta A_\mu (𝐱,t)=(0,\delta 𝐒(𝐱,t))`$ .; the four-vector spin distribution $`\delta a_\mu ^{(e)}(𝐩,𝐱,t)`$ is given by (11);$`_i=/x^i=^i`$.

Analogously the electron RKE (10) can be rewritten as the quasi-classical one,

$`{\displaystyle \frac{\delta f^{(e)}(𝐩,𝐱,t)}{t}}+\dot{𝐱}{\displaystyle \frac{\delta f^{(e)}(𝐩,𝐱,t)}{𝐱}}+Tr\left(\dot{\widehat{𝐩}}{\displaystyle \frac{\widehat{f}^{(0e)}(\epsilon _p)}{𝐩}}\right)=0,`$ (17)
where the particle number distribution function is obtained via the summing over spin quantum numbers, $`\delta f^{(e)}(𝐩,𝐱,t)=Tr(\delta \widehat{f}^{(e)}(𝐩,𝐱,t)`$; $`\dot{𝐱}=𝐯=𝐩/\epsilon _p`$ is the electron velocity;

$`\dot{\widehat{𝐩}}=`$ $`\left(e\left(\delta 𝐄(𝐱,t)+[𝐯\times 𝐁(𝐱,t)]\right)+G_F\sqrt{2}c_V\left[{\displaystyle \frac{p^\mu }{\epsilon _p}}\delta j_\mu ^{(\nu )}(𝐱,t)\right]\right)\delta _{\lambda \lambda ^{^{}}}+`$ (18)
$`+G_F\sqrt{2}c_A{\displaystyle \frac{m_e}{\epsilon _p}}\delta j_\mu ^{(\nu )}(𝐱,t)a_k^\mu (𝐩)(\sigma _k)_{\lambda \lambda ^{^{}}}`$
is the force matrix which accounts for the Lorentz force with electromagnetic fields $`\delta 𝐄`$, $`𝐁=𝐁_0+\delta 𝐁`$ as well as the weak interaction terms $`G_F`$ 333 Such tems can be also obtained from the weak interaction Hamiltonian for a probe electron moving in a neutrino medium
$$\widehat{H}_{weak}=\widehat{L}_{weak}=G_F\sqrt{2}\frac{\overline{u}_\lambda (p)\gamma ^\mu u_\lambda ^{^{}}(p)}{2\epsilon _p}\delta j_\mu ^{(\nu )}(𝐱,t)+G_F\sqrt{2}c_A\frac{\overline{u}_\lambda (p)\gamma ^\mu \gamma _5u_\lambda ^{^{}}(p)}{2\epsilon _p}\delta j_\mu ^{(\nu )}(𝐱,t).$$
;

$$\widehat{f}^{(0e)}(\epsilon _p)=\frac{\delta _{\lambda ^{^{}}\lambda }}{2}f^{0e}(\epsilon _p)+\frac{(\stackrel{}{\sigma }\widehat{\stackrel{}{b}}^{(0)})_{\lambda ^{^{}}\lambda }}{2}S^{(0e)}(\epsilon _p)$$
(19)
is the equilibrium density matrix as given in (2) with the ort $`\widehat{𝐛}^{(0)}=𝐁_0/B_0`$ separated from the isotropic distribution function $`S^{(0e)}(\epsilon _p)=(\mu _BB_0/\gamma )\text{d}f^{(0e)}/\text{d}\epsilon _p`$ given by (3).

Substituting (18), (19) into (17) one can easily check that the electron RKE (10) takes the form

$`{\displaystyle \frac{\delta f^{(e)}(𝐩,𝐱,t)}{t}}+𝐯{\displaystyle \frac{\delta f^{(e)}(𝐩,𝐱,t)}{𝐱}}+e\delta F_{j\mu }(𝐱,t){\displaystyle \frac{p^\mu }{\epsilon _p}}{\displaystyle \frac{f_0^{(e)}(\epsilon _p)}{p_j}}+eF_{j\mu }^{(0)}{\displaystyle \frac{p^\mu }{\epsilon _p}}{\displaystyle \frac{\delta f^{(e)}(𝐩,𝐱,t)}{p_j}}+`$
$`+G_F\sqrt{2}c_V\left({\displaystyle \frac{\delta n^{(\nu )}(𝐱,t)}{x_j}}v_k{\displaystyle \frac{\delta j_k^{(\nu )}(𝐱,t)}{x_j}}\right){\displaystyle \frac{f_0^{(e)}(\epsilon _p)}{p_j}}+`$
$`+G_F\sqrt{2}c_Am_e\left({\displaystyle \frac{\delta n^{(\nu )}(𝐱,t)}{x_j}}{\displaystyle \frac{a_0(𝐩)}{\epsilon _p}}{\displaystyle \frac{a_k(𝐩)}{\epsilon _p}}{\displaystyle \frac{\delta j_k^{(\nu )}(𝐱,t)}{x_j}}\right){\displaystyle \frac{S^{(0e)}(\epsilon _p)}{p_j}}=0,`$ (20)
where we introduced the unit polarization four-vector $`a^\mu (𝐩)=a_i^\mu (𝐩)\widehat{b}_i^{(0)}`$, $`a_\mu a^\mu =1`$.

On first glance both the neutrino (15) and electron (17) RKE’s do not obey the lepton number conservation law, $`j_\mu ^{(a)}(𝐱,t)/x_\mu 0`$, (due to the last terms in (16) and the last three terms in (2.3) correspondingly).

This non-conservation is true for the gauge non-invariant Wigner distributions $`f^{(a)}(𝐩,𝐱,t)\stackrel{~}{f}^{(a)}(𝐩,𝐱,t)=e^{i\mathrm{𝐩𝐲}}d^3y\stackrel{~}{f}^{(a)}(𝐱𝐲/2,𝐱+𝐲/2,t)`$ we used above, for which the distribution function in the coordinate representation (within integrand),

$$\stackrel{~}{f}^{(a)}(𝐱_1,𝐱_2,t)=Tr\left(\widehat{\rho }(t)\widehat{\mathrm{\Psi }}^{(a)+}(𝐱_2)\widehat{\mathrm{\Psi }}^{(a)}(𝐱_1)\right),$$

is not invariant with respect to gauge transformation of the corresponding interaction fields.

For example, in QED plasma the kinetic equation for the gauge non-invariant distribution of charged particles $`\stackrel{~}{f}^{(e)}(𝐱_1,𝐱_2,t)`$ derived from the quantum Liouville equation by the same Bogolyubov method does not obey the electric current conservation law since the force term depends on the electromagnetic potentials $`A_\mu (𝐱,t)`$ which do not enter as combinations expressed via field strengths, $`𝐄`$, $`𝐁`$ .

The recasting of such RKE for the gauge-invariant Wigner distribution

$$f^{(e)}(𝐩,𝐱,t)=\stackrel{~}{f}^{(e)}(𝐩+e𝐀(𝐱,t),𝐱,t)=d^3ye^{i\mathrm{𝐩𝐲}}f^{(e)}(𝐱𝐲/2,𝐱+𝐲/2,t),$$
(21)
where the gauge invariant distribution in the coordinate representation $`f^{(e)}(𝐱_1,𝐱_2,t)`$ is connected with the gauge non-invariant one,
$`\stackrel{~}{f}^{(e)}(𝐱_1,𝐱_2,t)=Tr\left(\widehat{\rho }(t)\widehat{\mathrm{\Psi }}^{(e)+}(𝐱_2)\widehat{\mathrm{\Psi }}^{(e)}(𝐱_1)\right)`$, by the important phase factor :

$$f^{(e)}(𝐱_1,𝐱_2,t)=\mathrm{exp}\left[ie(𝐱_2𝐱_1)_0^1𝑑\xi 𝐀(𝐱_2+\xi (𝐱_1𝐱_2),t)\right]\stackrel{~}{f}^{(e)}(𝐱_1,𝐱_2,t),$$
(22)
allows to obtain the usual Lorentz force term in the standard Boltzman equation for charged particles :

$$\frac{f^{(e)}(𝐩,𝐱,t)}{t}+𝐯\frac{f^{(e)}(𝐩,𝐱,t)}{𝐱}+e\left(𝐄+[𝐯\times 𝐁]\right)\frac{f^{(e)}(𝐩,𝐱,t)}{𝐩}=0,$$

for which, of course, the electric current $`j_\mu ^{(e)}(𝐱,t)=d^3p(p_\mu /\epsilon _p)f^{(e)}(𝐩,𝐱,t)/(2\pi )^3`$ is conserved, $`j_\mu ^{(e)}/x_\mu =0`$.

Such conservation is stipulated by the presence of the phase factor in the distribution (22) which in turn is invariant under the standard gauge transformation (with an arbitrary gauge function $`\chi (𝐱,t)`$ obeying the d’Alambert equation),

$`\widehat{\mathrm{\Psi }}^{(e)}(𝐱_1)e^{ie\chi (𝐱_1,t)}\widehat{\mathrm{\Psi }}^{(e)}(𝐱_1),`$
$`\widehat{\mathrm{\Psi }}^{(e)+}(𝐱_2)e^{+ie\chi (𝐱_2,t)}\widehat{\mathrm{\Psi }}^{(e)}(𝐱_2),`$
$`𝐀(𝐱_2+\xi (𝐱_1𝐱_2),t)𝐀(𝐱_2+\xi (𝐱_1𝐱_2),t)`$
$`{\displaystyle \frac{\chi (𝐱_2+\xi (𝐱_1𝐱_2),t)}{𝐱_2}},`$ (23)
or, equivalently, this arbitrary phase $`\chi (𝐱,t)`$ cancels in (22). Such invariance is crucial for macroscopic physics since it provides the physical sense of the Wigner function (21) and the conservation of the macroscopic electric current.

Deriving electron RKE (17) we used such gauge transformation while for the weak interaction terms we did not.

Thus, a recipe of the gauge invariance restoration in RKE’s (15), (17) should be similar to the change of arguments in (21): we should substitute the generalized (congugate) neutrino momentum,

$$𝐐=𝐪+G_F\sqrt{2}c_V𝐣^{(e)}(𝐱,t)+\frac{G_F}{\sqrt{2}}c_A𝐀^{(e)}(𝐱,t),$$
(24)
for the neutrino gauge invariant Wigner distribution

$$f^{(\nu )}(𝐪,𝐱,t)=\stackrel{~}{f}^{(\nu )}(𝐐,𝐱,t),$$
(25)
and the matrix of the generalized momentum

$$\widehat{𝐏}=[𝐩+G_F\sqrt{2}c_V\delta 𝐣^{(\nu )}(𝐱,t)]\delta _{\lambda \lambda ^{^{}}}G_F\sqrt{2}c_A\delta n^{(\nu )}(𝐱,t)(\stackrel{}{\sigma })_{\lambda \lambda ^{^{}}},$$
(26)
for the electron gauge invariant Wigner distribution

$$f^{(e)}(𝐩,𝐱,t)=\stackrel{~}{f}^{(e)}(𝐏,𝐱,t),$$
(27)
where $`𝐏=2^1Tr\widehat{𝐏}`$.

The RKE’s for gauge invariant (25),(27) are obtained from (15), (17) as

$$\frac{\delta \stackrel{~}{f}^{(\nu )}(𝐐,𝐱,t)}{t}+\dot{𝐱}\frac{\delta \stackrel{~}{f}^{(\nu )}(𝐐,𝐱,t)}{𝐱}+\dot{𝐐}\frac{\stackrel{~}{f}_0^{(\nu )}(𝐐)}{𝐐}=0,$$

and

$$\frac{\delta \stackrel{~}{f}^{(e)}(𝐏,𝐱,t)}{t}+\dot{𝐱}\frac{\delta \stackrel{~}{f}^{(e)}(𝐏,𝐱,t)}{𝐱}+Tr\left[\dot{\widehat{𝐏}}\frac{\widehat{f}_0^{(e)}(\epsilon _P)}{𝐏}\right]=0,$$

where generalized momenta $`𝐐`$, $`\widehat{𝐏}`$ are given by Eqs. (24) and (26) correspondingly; the background distributions $`\stackrel{~}{f}_0^{(a)}`$ can be changed to $`f_0^{(\nu )}(𝐪)`$ for neutrinos and $`\widehat{f}_0^{(e)}(\epsilon _p)`$ (19) for electrons since we retain the lowest order $`G_F`$ in the corresponding weak interaction terms.

Substituting in last RKE’s the momenta (24) and (26), for which we take into account both the derivatives (16), (18) and the total time derivatives $`\mathrm{d}\delta 𝐣^{(\mathrm{e})}(𝐱,\mathrm{t})/\mathrm{dt}=\delta 𝐣^{(\mathrm{e})}(𝐱,\mathrm{t})/\mathrm{t}+(𝐧)\delta 𝐣^{(\mathrm{e})}(𝐱,\mathrm{t})`$, $`\mathrm{d}\delta 𝐀^{(\mathrm{e})}(𝐱,\mathrm{t})/\mathrm{dt}=\delta 𝐀^{(\mathrm{e})}(𝐱,\mathrm{t})/\mathrm{t}+(𝐧)\delta 𝐀^{(\mathrm{e})}(𝐱,\mathrm{t})`$ in the neutrino RKE, $`\mathrm{d}\delta \mathrm{j}_\mu ^{(\nu )}(𝐱,\mathrm{t})/\mathrm{dt}=\delta \mathrm{j}_\mu ^{(\nu )}(𝐱,\mathrm{t})/\mathrm{t}+(𝐯)\delta \mathrm{j}_\mu ^{(\nu )}(𝐱,\mathrm{t})`$ in the electron one, and then using the identities

$$(𝐧)\delta 𝐣^{(e)}(𝐱,t)(𝐧\delta 𝐣^{(e)}(𝐱,t))[𝐧\times \times \delta 𝐣^{(e)}(𝐱,t)],$$

$$(𝐧)\delta 𝐀^{(e)}(𝐱,t)(𝐧\delta 𝐀^{(e)}(𝐱,t))[𝐧\times \times \delta 𝐀^{(e)}(𝐱,t)],$$

$$(𝐯)\delta 𝐣^{(\nu )}(𝐱,t)(𝐯\delta 𝐣^{(\nu )}(𝐱,t))[𝐯\times \times \delta 𝐣^{(\nu )}(𝐱,t)],$$

$$(𝐯)\delta n^{(\nu )}(𝐱,t)\widehat{𝐛}^{(0)}(𝐯\widehat{𝐛}^{(0)}\delta n^{(\nu )}(𝐱,t))(𝐯\times \times \widehat{𝐛}^{(0)}\delta n^{(\nu )}(𝐱,t)),$$

we arrive to the final forms of the lepton RKE’s.

Namely, accounting for the definition of the neutrino gauge invariant distribution (25) the neutrino RKE (15) which is equivalent to the master (2.2) takes finally in the Vlasov approximation the form

$`{\displaystyle \frac{\delta f^{(\nu )}(𝐪,𝐱,t)}{t}}+`$ $`𝐧{\displaystyle \frac{\delta f^{(\nu )}(𝐪,𝐱,t)}{𝐱}}+\delta F_{j\mu }^{(V)}(𝐱,t){\displaystyle \frac{q^\mu }{\epsilon _q}}{\displaystyle \frac{f_0^{(\nu )}(𝐪)}{q_j}}+`$ (28)
$`+\delta F_{j\mu }^{(A)}(𝐱,t){\displaystyle \frac{q^\mu }{\epsilon _q}}{\displaystyle \frac{f_0^{(\nu )}(𝐪)}{q_j}}=0,`$
where the antisymmetric tensors $`\delta F_{jk}^{(V,A)}(𝐱,t)`$ entering effective Lorentz force terms are given by the vector and axial vector currents correspondigly,

$`\delta F_{j0}^{(V)}(𝐱,t)/G_F\sqrt{2}c_V=_j\delta n^{(e)}(𝐱,t){\displaystyle \frac{\delta j_j^{(e)}(𝐱,t)}{t}},`$
$`\delta F_{jk}^{(V)}(𝐱,t)/G_F\sqrt{2}c_V=e_{jkl}(\times \delta 𝐣^{(e)}(𝐱,t))_l,`$
$`\sqrt{2}\delta F_{j0}^{(A)}(𝐱,t)/G_Fc_A=_j\delta A_0^{(e)}(𝐱,t){\displaystyle \frac{\delta A_j^{(e)}(𝐱,t)}{t}},`$
$`\sqrt{2}\delta F_{jk}^{(A)}(𝐱,t)/G_Fc_A=e_{jkl}(\times \delta 𝐀^{(e)}(𝐱,t))_l;`$ (29)
and $`/𝐐=/𝐪`$ from the relation (24).

First three terms in RKE (28) were derived in . In the paper  analogous vector coupling terms ($`c_V`$ for electron neutrinos) were obtained in cold hydrodynamics.

Analogously for the electron gauge invariant distribution (27) we obtain finally from (2.3) (equivalent to the master equation (10)) the electron RKE

$`{\displaystyle \frac{\delta f^{(e)}(𝐩,𝐱,t)}{t}}+𝐯{\displaystyle \frac{\delta f^{(e)}(𝐩,𝐱,t)}{𝐱}}+e\delta F_{j\mu }(𝐱,t){\displaystyle \frac{p^\mu }{\epsilon _p}}{\displaystyle \frac{f_0^{(e)}(\epsilon _p)}{p_j}}+eF_{j\mu }^{(0)}{\displaystyle \frac{p^\mu }{\epsilon _p}}{\displaystyle \frac{\delta f^{(e)}(𝐩,𝐱,t)}{p_j}}+`$
$`+F_{j\mu }^{(weak)}(𝐱,t){\displaystyle \frac{p^\mu }{\epsilon _p}}{\displaystyle \frac{f_0^{(e)}(\epsilon _p)}{p_j}}G_F\sqrt{2}c_A[{\displaystyle \frac{\delta n^{(\nu )}(𝐱,t)\widehat{𝐛}^{(0)}}{t}}(𝐯\times \times \delta n^{(\nu )}(𝐱,t)\widehat{𝐛}^{(0)})+`$
$`+{\displaystyle \frac{m_e}{\epsilon _p}}(𝐚(𝐩)\delta 𝐣^{(e)}(𝐱,t))]{\displaystyle \frac{S^{(0e)}(\epsilon _p)}{p_j}}=0,`$ (30)
where the weak vector term ($`c_V`$) has the Lorentz structure with the tensor components

$`F_{j0}^{(weak)}(𝐱,t)/G_F\sqrt{2}c_V=_j\delta n^{(\nu )}(𝐱,t){\displaystyle \frac{\delta j_j^{(\nu )}(𝐱,t)}{t}},`$
$`F_{jk}^{(weak)}(𝐱,t)/G_F\sqrt{2}c_V=e_{jkl}(\times \delta 𝐣^{(\nu )}(𝐱,t))_l;`$ (31)
$`/𝐏=/𝐩`$ from the relation (26); and the axial vector term does not contribute to the continuity equation since the last term with the 3-vector component of the four-vector $`a_\mu (𝐩)`$ introduced after (2.3) ,

$$𝐚(𝐩)=\widehat{𝐛}^{(0)}+\frac{𝐩(𝐩\widehat{𝐛}^{(0)})}{m_e(\epsilon _p+m_e)},$$

vanishes under integration $`d^3p(a_k(𝐩)/\epsilon _p)S^{(0e)}(\epsilon _p)/p_j`$ due to the isotropy of the background and the odd power of momenta $`𝐩`$ under the integral.

Obviously, RKE’s (28), (2.3) obey the lepton current conservation law

$$\frac{j_\mu ^{(a)}(𝐱,t)}{x_\mu }=0,$$
(32)
or the lepton currents are conserved in a magnetized plasma.

### 2.4 Dispersion equation for lepton density perturbations in isotropic plasma

Here checking master equations (2.2), (10) against known results for an isotropic plasma in the absence of an external magnetic field, $`B_0=0`$, hence neglecting spin waves, we show that accounting for the lepton current conservation (32) such observables like the spectra of lepton density perturbations do not depend whether we apply initial RKE’s for gauge non-invariant distribution functions (2.2), (10), or same equations written in the completed forms, (28),(2.3).

The situation is similar to the case of standard plasma where the initial RKE for the gauge non-invariant Wigner distribution $`\stackrel{~}{f}^{(e)}(𝐩,𝐱,t)`$ is often more suitable for the obtaining of concrete results than Boltzman equation with the Lorentz force . Nevertheless, the electric current should be written in the gauge invariant form obeying (32).

As well as authors  we consider the electron plasma waves (EPW) driven by the intense neutrino flux, neglecting also the transverse wave contribution, i.e. retaining only the electric field in the Lorentz force for electrons.

Using the Fouirer representation in the RKE’s (2.2), (10) one can easily obtain the algebraic system for the current perturbations
$`\delta j_\mu ^{(\nu )}(\omega ,𝐤)=d^3q(q_\mu /q)\delta f^{(\nu )}(𝐪,𝐤,\omega )/(2\pi )^3`$,
$`\delta j_\mu ^{(e)}(\omega ,𝐤)=d^3p(p_\mu /\epsilon _p)\delta f^{(e)}(𝐩,𝐤,\omega )/(2\pi )^3`$,

$`\delta n^{(\nu )}(\omega ,𝐤)+G_F\sqrt{2}c_V[\delta n^{(e)}(\omega ,𝐤){\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}{\displaystyle \frac{k_jf_0^{(\nu )}(q)/q_j}{\omega \mathrm{𝐤𝐧}}}`$
$`\delta j_k^{(e)}(\omega ,𝐤)k_j{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}{\displaystyle \frac{n_kf_0^{(\nu )}(q)/q_j}{\omega \mathrm{𝐤𝐧}}}]=0,`$
$`\delta n^{(e)}(\omega ,𝐤)[1+\chi _e(\omega ,𝐤)]+G_F\sqrt{2}c_V[\delta n^{(\nu )}(\omega ,𝐤){\displaystyle }{\displaystyle \frac{d^3p}{(2\pi )^3}}{\displaystyle \frac{k_jf_0^{(e)}(\epsilon _p)/p_j}{\omega \mathrm{𝐤𝐯}}}`$
$`\delta j_k^{(\nu )}(\omega ,𝐤)k_j{\displaystyle }{\displaystyle \frac{d^3p}{(2\pi )^3}}{\displaystyle \frac{v_kf_0^{(e)}(\epsilon _p)/p_j}{\omega \mathrm{𝐤𝐯}}}]=0,`$ (33)
where $`𝐧=𝐪/q`$ is the velocity of the massless neutrino and the term $`\chi _e(\omega ,𝐤)`$ is connected with the longitudinal permittivity,

$$\chi _e(\omega ,𝐤)=\epsilon _l(\omega ,𝐤)1=\frac{4\pi e^2}{k^2}\frac{d^3p}{(2\pi )^3}\frac{k_jf_0^{(e)}(\epsilon _p)/p_j}{\omega \mathrm{𝐤𝐯}}.$$
(34)
Note that starting instead of (2.2), (10) from the corresponding RKE’s (28), (2.3) derived in the previous subsection for the gauge invariant Wigner distribution functions $`f^{(\nu )}(𝐐,𝐱,t)`$, $`f^{(e)}(𝐏,𝐱,t)`$ we obtain the same algebraic system (2.4) omitting axial terms ($`c_A`$) for $`𝐁_0=0`$.

Really additional terms coming from the total time derivatives,

$$\frac{\delta 𝐣^{(e)}(𝐱,t)}{t}+(𝐧)\delta 𝐣^{(e)}(𝐱,t)$$

for $`\dot{𝐐}`$ via the second term in (24), and

$$\frac{\delta 𝐣^{(\nu )}(𝐱,t)}{t}+(𝐯)\delta 𝐣^{(\nu )}(𝐱,t)$$

for $`\dot{𝐏}`$ via the second term in (26), are proportional in the Fourier representation to the $`\stackrel{ˇ}{C}`$erenkov factors $`\omega \mathrm{𝐤𝐧}`$, $`\omega \mathrm{𝐤𝐯}`$ correspondingly that exactly cancels with the poles within integrands in (2.4) leading to zero contributions in brackets for both equations (2.4): $`\delta j^{(e)}(\omega ,𝐤)(d^3q/(2\pi )^3)f_0^{(\nu )}(q)/q_j=0`$, and
$`\delta j^{(\nu )}(\omega ,𝐤)(d^3p/(2\pi )^3)f_0^{(e)}(\epsilon _p)/p_j=0`$ .

Now using the symmetrical tensor entering the second line in (2.4) (and the analogous one in the fourth line)

$$\frac{d^3q}{(2\pi )^3}\frac{n_kf_0^{(\nu )}(q)/q_j}{\omega \mathrm{𝐤𝐧}}=A\delta _{kj}+B\frac{k_kk_j}{k^2},$$

where the scalars $`A`$, $`B`$ are given by

$$A=\frac{1}{2}\left[\frac{d^3q}{(2\pi )^3(\omega \mathrm{𝐤𝐧})}\left(n_i\frac{f_0^{(\nu )}(q)}{q_i}\frac{\omega }{k^2}k_i\frac{f_0^{(\nu )}(q)}{q_i}\right)\right]$$

$$B=\frac{1}{2}\left[\frac{d^3q}{(2\pi )^3(\omega \mathrm{𝐤𝐧})}\left(\frac{3\omega }{k^2}k_i\frac{f_0^{(\nu )}(q)}{q_i}n_i\frac{f_0^{(\nu )}(q)}{q_i}\right)\right],$$

we obtain the factor multiplying the electron current $`\delta j_k^{(e)}`$ in the second line (and the analogous one in the fourth line in (2.4))

$$k_j\frac{d^3q}{(2\pi )^3}\frac{n_kf_0^{(\nu )}(q)/q_j}{\omega \mathrm{𝐤𝐧}}=k_k(A+B)=\frac{k_k\omega }{k^2}\frac{d^3q}{(2\pi )^3(\omega \mathrm{𝐤𝐧})}k_i\frac{f_0^{(\nu )}(q)}{q_i}.$$
(35)
Substituting eqs. (34), (35) we rewrite the system (2.4) as

$`\delta n^{(\nu )}(\omega ,𝐤)+G_F\sqrt{2}c_V[\delta n^{(e)}(\omega ,𝐤){\displaystyle \frac{\omega }{k^2}}k_k\delta j_k^{(e)}(\omega ,𝐤)]\times `$
$`\times {\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3(\omega \mathrm{𝐤𝐧})}}k_i{\displaystyle \frac{f_0^{(\nu )}(q)}{q_i}}=0,`$
$`G_F\sqrt{2}c_V\left[\delta n^{(\nu )}(\omega ,𝐤){\displaystyle \frac{\omega }{k^2}}k_k\delta j_k^{(\nu )}(\omega ,𝐤)\right]{\displaystyle \frac{k^2\chi _e(\omega ,𝐤)}{4\pi e^2}}+`$
$`+\delta n^{(e)}(\omega ,𝐤)\left(1+\chi _e(\omega ,𝐤)\right)=0.`$ (36)
Using the lepton current conservation (32), $`k_k\delta j_k^{(a)}(\omega ,𝐤)=\omega \delta n^{(a)}(\omega ,𝐤)`$, we obtain from (2.4) the dispersion equation in the lowest order over the Fermi constant ($`G_F^2`$) :

$$1+\chi _e(\omega ,𝐤)\frac{\mathrm{\Delta }_\nu k^2c_V^2}{\omega _{pe}^2}\left(1\frac{\omega ^2}{k^2}\right)^2\chi _e(\omega ,𝐤)\frac{d^3q}{(2\pi )^3(\omega \mathrm{𝐤𝐧})}k_i\frac{\widehat{f}_0^{(\nu )}(q)}{q_i}=0,$$
(37)
where $`\mathrm{\Delta }_\nu =2G_F^2n_{e0}n_{\nu 0}/m_e`$ is the weak parameter introduced in ; $`\widehat{f}_0^{(\nu )}=f_0^{(\nu )}(q)/n_{\nu 0}`$ is the normalized neutrino background distribution function; $`\omega _{pe}=\sqrt{4\pi \alpha n_{e0}/m_e}`$ is the non-relativistic expression for the plasma frequency.

Let us stress the importance of additional terms in RKE’s (28), (2.3) discussed after eq. (34) which did not contribute to (2.4) while they become important when we claim the lepton current conservation deriving dispersion equation (37). We can conclude that the master RKE’s (2.2), (10) plus the lepton current conservation laws provided by corresponding eqs. (28), (2.3) lead to the correct issues known in literature, in particular, to the dispersion equation (37).

The prediction of the streaming instability driven by the neutrino beam in a dense plasma (i.e., outside the neutrinosphere of a supernova)  was criticized in . In this section, however, we do not touch these issues following from eq. (37).

Let us turn to the case of magnetized plasma.

## 3 Bloch equation in the presence of neutrino beam

Substituting the linear decomposition (2.1) into eq. (2.2), then integrating the latter over $`d^3p`$ and multiplying by $`\mu _B`$ we have generalized Bloch equation for the magnetization density perturbations $`𝐦(𝐱,t)=\mu _Bd^3p\delta 𝐒(𝐩,𝐱,t)/(2\pi )^3`$ in a polarized NR electron gas that is the base of the theory of the electron paramagnetic resonance in the absence of weak interactions  while in SM with neutrinos such equation takes the form

$`{\displaystyle \frac{m_j(𝐱,t)}{t}}+`$ $`2\mu _B\left[\left(𝐦(𝐱,t)\chi _0𝐛(𝐱,t)\right)\times 𝐁_0\right]_j+{\displaystyle \frac{\mathrm{\Sigma }_{kj}(𝐱,t)}{x_k}}+`$ (38)
$`+{\displaystyle \frac{G_F\sqrt{2}c_A\mu _Bn_{0e}}{m_e}}{\displaystyle \frac{d^3q}{(2\pi )^3}\frac{\delta f^{(\nu )}(𝐪,𝐱,t)}{x_j}}=0.`$
Here $`\chi _0=2\mu _B^2d^3p(df^{(0)}/d\epsilon _p)/(2\pi )^3`$ is the static susceptibility of the polarized electron gas 444 The static susceptibility is small in a degenerate electron gas , $`\chi _0=\alpha v_{F_e}/4\pi ^21`$ , in contrast with the varying one $`\chi (t)`$ (see below eq (41)). This is the reason why the static magnetic induction $`B_0=(1+4\pi \chi _0)H_0`$ and the magnetic field strength $`H_0`$ practically coincide there.; $`𝐛(𝐱,t)`$ is the magnetic field perturbation in the total field $`𝐁(𝐱,t)=𝐁_0+𝐛(𝐱,t)`$; $`n_{e0}=d^3pf^{(0)}(\epsilon _p)/(2\pi )^3`$ is the density of the electron gas. The pseudotensor $`\mathrm{\Sigma }_{kj}(𝐱,t)=\mu _Bd^3pv_k\delta S_j(𝐩,𝐱,t)/(2\pi )^3`$ describes spatial inhomogeneity of the magnetization and the last new vector term ($`c_A`$) corresponds to the parity violation for the evolution of the macroscopic axial vector $`m_j`$ in SM.

Note that last two lines in the complicated spin RKE (2.2) are the relativistic corrections and do not contribute to (38) in NR plasma. Moreover, we omitted the small term $`G_F\sqrt{2}c_A\chi _0[𝐣^{(\nu )}(𝐱,t)\times 𝐁_0]_j`$ originated from the third line (without derivatives) in (2.2) when comparing it with the third term in the Bloch equation (38) stipulated by the standard spin precession in QED plasma (arising due to the last term in the second line in (2.2)).

We neglected also the fifth term in (2.2) that is proportional to the vector coupling $`c_V`$ since the background polarization is small in our WKB approximation, $`S^{(0e)}(\mu _BB_0/\gamma )\text{d}f^{(0e)}/\text{d}\epsilon _pf^{(0e)}(\epsilon _p)`$, and this term is much less than the previous one in (2.2) which is proportional to the axial coupling $`c_A`$ retained in (38). This is in agreement with estimates of the small polarization $`<\lambda >=n_{e0}/n_e0.01`$ made in  (Nunokawa et al) even for strong magnetic field in SN where $`n_{0e}`$ is the electron density on the main Landau level, $`n_e`$ is the total electron density.

In order to break the chain of hydrodynamical equations for magnetization moments, $`m_j`$, $`\mathrm{\Sigma }_{ij}`$, etc, one can consider long-wave spin waves with the wave lengths $`\lambda _s`$ that are much bigger than the electron Larmour radius, $`\lambda _sl_{eH}`$. In such case the spin perturbation function is given by the first two moments:

$$\mu _B\delta S_j(𝐩,𝐱,t)=\frac{df^{(0)}}{d\epsilon _p}\left[\frac{d^3p}{(2\pi )^3}\frac{df^{(0)}}{d\epsilon _p}\right]^1\times \left(m_j(𝐱,t)+\frac{3v_i}{v^2}\mathrm{\Sigma }_{ij}(𝐱,t)\right),$$

that allows us to complete the hydrodynamical system by the equation for the pseudotensor $`\mathrm{\Sigma }_{ij}`$,

$`{\displaystyle \frac{\mathrm{\Sigma }_{ij}(𝐱,t)}{t}}+`$ $`\mathrm{\Omega }_e\left(e_{jtl}\widehat{n}_l^B\mathrm{\Sigma }_{it}(𝐱,t)e_{itl}\widehat{n}_l^B\mathrm{\Sigma }_{tj}(𝐱,t)\right)+`$ (39)
$`+{\displaystyle \frac{<v^2>}{3}}{\displaystyle \frac{}{x_i}}\left(m_j(𝐱,t)\chi _0b_j(𝐱,t)\right)`$
$`{\displaystyle \frac{G_Fc_A\mu _Bn_{e0}}{\sqrt{2}m_e}}{\displaystyle \frac{d^3q}{(2\pi )^3}n_j(q)\frac{\delta f^{(\nu )}(𝐪,𝐱,t)}{x_j}}=0.`$
Here $`\mathrm{\Omega }_e=eB_0/m_e`$ is the electron cyclotron frequency;
$`<v^2>=[d^3p𝑑f^{(0)}/𝑑\epsilon _p]^1d^3pv^2𝑑f^{(0)}/𝑑\epsilon _p`$ is the average of the velocity squared; $`\widehat{n}^B=𝐁_0/B_0`$ is the unit pseudovector and last true tensor term describes the parity violation through weak interactions. Note that eq. (39) would be important accounting for the large spatial dispersion $`k<v>\omega `$ outside the region (40) and claims an inclusion of exchange interactions (from the ee-scattering exchange Feynman diagram) we omitted here.

## 4 Neutrino driven streaming instability of spin waves

In this section we derive from generalized Bloch equation (38) the dispersion equation in magnetized plasma in the presence of neutrino beam.

Apparently master RKE’s for number density distribution functions (2.2), (10) should be consistent with the spin RKE (2.2) and its consequence (38). Therefore we are checking below whether more general (28), (2.3) are necessary. As we find below there is no difference between the use of these equations for the spin wave propagation in plasma.

Let us consider for simplicity the case of long-wave perturbations with the spectrum $`\omega (𝐤)`$ obeying the inequalities

$$\frac{k^2}{\omega }\omega k<v>.$$
(40)
As the mean electron velocity $`<v>`$ is small in NR plasma, $`<v>1`$, and keeping in mind the Maxwell equation written in the Fourier representation $`𝐛=[𝐤\times \delta 𝐄]/\omega `$ one finds from the condition $`k/\omega <v>`$ that the electric field contribution in the spin RKE (2.2) ($`\delta 𝐄`$) is negligible comparing with the magnetic one even for the maximum frequency in whole space-like region $`\omega k`$ relevant to the $`\stackrel{ˇ}{C}`$erenkov resonance with neutrinos, $`\omega =\mathrm{𝐤𝐧}`$. Without electric fields the transversal components of the permeability and the susceptibility tensors ($`i,j=x,y`$) appearing in the second term of (38) are diagonal in plasma, $`\mu _{ij}(\omega ,𝐤)=\mu (\omega ,𝐤)\delta _{ij}`$, $`\chi _{ij}(\omega ,𝐤)=\chi (\omega ,𝐤)\delta _{ij}`$, that differs this medium from ferromagnets.

Moreover, the latter inequality in (40), $`\omega k<v>`$, means the high-frequency approximation, for which the spatial uniform susceptibility $`\chi (t)`$ and the permeability $`\mu (t)=1+4\pi \chi (t)`$ can be considered instead of more general ones, $`\chi (𝐱,t)`$ and $`\mu (𝐱,t)`$, that allows us to neglect the complicated pseudotensor term $`\mathrm{\Sigma }_{ij}`$ in (38) as well as the spatial dispersion in the perturbations $`m_\pm (𝐱,t)=𝑑t^{}\chi (tt^{})h_\pm (𝐱,t^{})`$ and $`b_\pm (𝐱,t)=𝑑t^{}\mu (tt^{})h_\pm (𝐱,t^{})`$ correspondingly where $`h_\pm (𝐱,t)=h_x(𝐱,t)\pm ih_y(𝐱,t)`$ is the magnetic field strength perturbation.

Under such conditions we find from the linearized spin RKE (2.2) the susceptibility in QED plasma neglecting neutrinos,

$$\chi (\omega )=\frac{\pm \mathrm{\Omega }_e\chi _0}{\omega \pm \mathrm{\Omega }_e(14\pi \chi _0)},$$
(41)
that may be finite at the paramagnetic resonance $`\omega =\mathrm{\Omega }_e`$ given by the equation $`1+4\pi \chi (\omega )=0`$. Latter follows from the shortened Maxwell equation $`𝐤\times (𝐛4\pi 𝐦)=0`$ when one neglects electric field terms.

Then substituting into (38) the solution of the neutrino RKE (2.2) obtained in the same linear approximation (2.1),

$`\delta f^{(\nu )}(𝐪,𝐤,\omega )=G_F\sqrt{2}c_V{\displaystyle \frac{k_k}{(\omega \mathrm{𝐤𝐧})}}{\displaystyle \frac{f_0^{(\nu )}(q)}{q_k}}{\displaystyle \frac{d^3p_e}{(2\pi )^3}(1𝐯_e𝐧)\delta f^{(e)}(𝐩_e,𝐤,\omega )}+`$
$`+{\displaystyle \frac{G_Fc_A}{\sqrt{2}}}{\displaystyle \frac{k_kn_l(q)}{(\omega \mathrm{𝐤𝐧})}}{\displaystyle \frac{f_0^{(\nu )}(q)}{q_k}}{\displaystyle \frac{d^3p_e}{(2\pi )^3}\delta S_l(𝐩_e,𝐤,\omega )},`$ (42)
we can rewrite the generalized Bloch equation (38) in the Fourier representation as

$`i\omega m_j(𝐤,\omega )+`$ $`2\mu _B[(𝐦(𝐤,\omega )\chi _0𝐛(𝐤,\omega ))\times 𝐁_0]_j+{\displaystyle \frac{c_A^2\mathrm{\Delta }_\nu }{2}}A_l^{(\nu )}(\omega ,𝐤)m_l(𝐤,\omega )ik_j+`$ (43)
$`+{\displaystyle \frac{c_Ac_V\mathrm{\Delta }_\nu k_j}{8\pi m_e}}(A_i^{(\nu )}(\omega ,𝐤){\displaystyle \frac{B^{(\nu )}(\omega ,𝐤)}{\omega }}k_i)\times `$
$`\times (e_{ikl}k_k[b_l(𝐤,\omega )4\pi m_l(𝐤,\omega )]+\omega \delta 𝐄_i)=0.`$
Here the factor $`\mathrm{\Delta }^\nu `$ is given after eq. (37); the vector $`A_l^{(\nu )}A_l^{(\nu )}(\omega ,𝐤)`$ and the scalar $`B^{(\nu )}B^{(\nu )}(\omega ,𝐤)`$ depend on the neutrino background distribution $`\widehat{f}_0^{(\nu )}(𝐪)`$,

$$A_l^{(\nu )}=\frac{d^3q}{(2\pi )^3}\frac{\widehat{f}_0^{(\nu )}(𝐪)}{q}\left(\frac{[(\mathrm{𝐤𝐧})^2k^2]n_l(q)}{(\omega \mathrm{𝐤𝐧})^2}+\frac{(\mathrm{𝐤𝐧})n_lk_l}{\omega \mathrm{𝐤𝐧}}\right),$$
(44)
$$B^{(\nu )}=\frac{d^3q}{(2\pi )^3}\frac{\widehat{f}_0^{(\nu )}(𝐪)}{q}\left(\frac{[(\mathrm{𝐤𝐧})^2k^2]}{(\omega \mathrm{𝐤𝐧})^2}\right).$$
(45)
Obtaining last term in eq. (43) ($`c_Ac_V`$) we used the exact Maxwell equation (14) when we substituted the convection current $`\delta j_i(𝐤,\omega )=d^3pv_i\delta f^{(e)}(𝐩,𝐤,\omega )`$ in the first line of eq. (4).

Note that the electromagnetic current $`j_\mu ^{(e)}(𝐱,t)`$ conserves automatically since we can neglect weak interactions in the electron RKE (10) while retaining them in the neutrino RKE (2.2). In turn the neutrino influence the electron spin evolution coming from the solution (4) is not changed if we would substitute the solution of the complete equation (28) instead of the initial one (2.2).

Really in the last case there appear the two additional terms in (4),

$$G_F\sqrt{2}c_V\frac{f_0^{(\nu )}(q)}{q_k}\delta j_k^{(e)}(\omega ,𝐤)\frac{G_F}{\sqrt{2}}c_A\frac{f_0^{(\nu )}(q)}{q_k}\delta m_k(\omega ,𝐤)$$
(46)
which do not contribute to the generalized Bloch equation (38).

One can easily see that in NR plasma with a low electron density $`n_{0e}m_e^32\times 10^{31}cm^3`$ the general condition of the macroscopic description $`2\pi /k(n_0)^{1/3}`$ means that long wave lengths exceeding the Compton one, $`2\pi /km_e^1`$ are possible only, and due to (40) frequencies obey the inequality $`\omega km_e`$.

Accounting for low $`\omega `$ and $`k`$ and comparing the term $`c_Vc_A`$ with the previous one in eq. (43) ($`c_A^2`$) one finds that the convection current gives a negligible contribution. Note that for $`\nu _{\mu ,\tau }`$-neutrinos with $`c_V1`$ such correction becomes even less than for $`\nu _e`$. Omitting also the small term $`\chi _0`$ we arrive to the shortened form of (43),

$$i\omega m_j(𝐤,\omega )+2\mu _B[𝐦(𝐤,\omega )\times 𝐁_0]_j+\frac{c_A^2\mathrm{\Delta }_\nu }{2}A_l^{(\nu )}(\omega ,𝐤)m_l(𝐤,\omega )ik_j=0.$$
(47)
Thus, from the generalized Bloch equation (47) we have derived finally the dispersion equation in SM,

$$(\omega ^2\mathrm{\Omega }_e^2)\left(\omega \frac{c_A^2\mathrm{\Delta }^{(\nu )}}{2}A_z^{(\nu )}k_z\right)\frac{\omega ^2c_A^2\mathrm{\Delta }^{(\nu )}}{4}(A_{}^{(\nu )}k_++A_+^{(\nu )}k_{})=0,$$
(48)
with the vector $`A_i^{(\nu )}`$ as given in eq. (44).

In the particular case of the neutrino beam, $`\widehat{f}_0^{(\nu )}(𝐪)=(2\pi )^3\delta ^{(3)}(𝐪𝐪_0)`$, $`𝐧_0=𝐪_0/q_0`$, for the transversal neutrino propagation $`𝐧_0𝐁_0`$ with the beam direction along x-axis, $`𝐧_0=(1,0,0)`$, and the magnetic field $`𝐁_0=(0,0,B_0)`$ we obtain from eq. (48) the resonant excitation of spin waves $`\omega =\mathrm{𝐤𝐧}_0+i\delta =\mathrm{\Omega }_e+i\delta `$ with the increment $`\delta `$,

$$\delta \mathrm{\Omega }_e^{1/3}\frac{\sqrt{3}}{4}\left(\frac{\mathrm{\Delta }^{(\nu )}}{q_0}\right)^{1/3}(\sqrt{2}c_Ak\mathrm{sin}\theta _{q_0})^{2/3}\mathrm{\Omega }_e\frac{\sqrt{3}}{4}\left(\frac{\mathrm{\Delta }^{(\nu )}}{q_0}\right)^{1/3}(\sqrt{2}c_A\mathrm{sin}\theta _{q_0})^{2/3},$$
(49)
where we substituted the scalar product $`\mathrm{𝐤𝐧}_0=k\mathrm{cos}\theta _{q_0}=k\mathrm{sin}\theta \mathrm{cos}\varphi `$ denoting $`\theta _{q_0}`$ as the angle between the neutrino beam direction and the wave vector $`𝐤`$; $`\theta `$ is the angle of $`𝐤`$ with respect to the magnetic field $`𝐁_0`$. Note that we relied in the last inequality on the frequency approximation $`k\omega \mathrm{\Omega }_e`$ assumed above.

Let us compare this increment with the fastest one in the case of isotropic plasma following from the dispersion equation (37) ,

$$\delta _{weak}=\frac{\sqrt{3}}{2}\omega _{pe}\left(\frac{\mathrm{\Delta }_\nu }{q_0}\frac{\mathrm{sin}^6\theta _{q_0}}{\mathrm{cos}^4\theta _{q_0}}\right)^{1/3}.$$
(50)
One can easily see the advantage of the excitation of collective modes by the intense neutrino flux in a polarized electron gas. For a strong magnetic field in a dense plasma obeying $`\mathrm{\Omega }_e>\omega _{pe}`$ the angular dependence in eq. (49) $`(\mathrm{sin}\theta _{q_0})^{2/3}`$ gives a less suppression of the increment for the small angles $`\theta _{q_0}\mathrm{arccos}(<v>/c)`$ for which Landau damping for growing modes is absent. This is due to the absence of the factor $`(1\omega ^2/k^2)^2`$ in the new dispersion equation (48). It is obvious that such angular dependence would be especially dangerous and important for the relativistic plasma case.

There is the second advantage of spin waves enhanced via the weak axial vector currents ($`c_A=0.5`$) instead of the case of plasma waves excited via the weak vector currents with the small vector coupling in the case of muon and tau neutrinos (choosing the lower sign in $`c_V=2\xi \pm 0.5`$ where $`\xi 0.23`$ is the Weinberg parameter ). This is the reason why authors  considered the case of electron neutrinos only and put for them $`c_V1`$. Note that during the main neutrino burst in SN all neutrino species are produced in the hot SN core via the pair annihilation $`e^+e^{}\nu _a\stackrel{~}{\nu }_a,`$ $`a=e,\mu ,\tau `$.

On the other hand, there are other arguments against streaming instability driven by neutrino beams in SN (e.g. not collimated beam)  to be relevant for neutrino propagation in a magnetized medium. Nevertheless, we have just showed that in a polarized electron gas the dispersion equations are quite different from the case of the isotropic plasma  that stimulates a future exploration of collective plasma phenomena in the presence of intense neutrino fluxes.

## 5 Discussion and conclusions

The above derivation shows that the Bogolyubov method starting from QFT Feynman diagrams remains a power and straightforward tool to obtain RKE’s. Being supported with an additional gauge restoring transformation similar to  it leads to the same results obtained within Hamiltonian approach of  for the case of neutrinos propagating in an isotropic plasma.

In the case of a polarized electron gas the master RKE’s derived above by the Bogolyubov method allow to analyse a new phenomenon- spin wave propagation in NR plasma enhanced by the neutrino beam. Note that complete RKE’s (28), (2.3) are not necessary to derive the dispersion equation (48) for spin waves as we showed in (46).

The violation of parity in the SM lepton plasma given by axial vector currents ($`c_A`$) leads to the growth of spin eigen modes through their excitation by the intense neutrino flux. These spin waves are generally coupled to the magnetosonic ones analogously to the case of spin waves in ferromagnets  or can transfer their energy to electromagnetic and plasma waves at the cross of spectra that finally could lead to the heating of ions and the background plasma.

The possible explanation of the shock revival in SN by different collective mechanisms including the neutrino driven streaming instability seems to be very perspective goal for future studies in the case of polarized electron gas.

The application of these mechanisms in the magnetized plasma behind the shock and outside the SN neutrinosphere is self-consistent with plasma parameters expected there. Really, we do not consider neutrino collisions with matter within this region using the collisionless neutrino RKE (2.2) and may also neglect the electron-ion collision frequency $`\nu _{ei}`$ comparing with the cyclotron frequency at the paramagnetic resonance $`\omega =\mathrm{\Omega }_e`$. Hence collisionless Vlasov approximation should be valid for the spin equation (2.2) as well as for the electron RKE (10) since the Debye number is large, $`N_D=n_{0e}r_D^31`$. Indeed, in the field $`B_0=10^{12}Gauss`$ the cyclotron frequency reaches $`\mathrm{\Omega }_e=1.7\times 10^7B_01.7\times 10^{19}sec^1`$. This frequency is comparable with the plasma one at the density $`n_{e0}10^{29}cm^3`$, $`\omega _{pe}=1.8\times 10^{19}sec^1`$, and turns out to be larger than e.g. the electron-proton collision frequency $`\nu _{ep}=50n_{e0}(cm^3)/(T_e(K))^{3/2}sec^11.6\times 10^{17}sec^1`$ in the surrounding NR plasma with the temperature $`T_e10^9K`$. For these parameters the Debye radius $`r_D10^9cm`$ corresponds to the large plasma parameter $`N_D1001`$.

Under such conditions for the mean neutrino energy $`q_010MeV`$ the increment (49) reaches the maximum value $`\delta 10^{10}sec^1`$ as in previous optimistic estimates. Note that this means too sharp collimated neutrino beam with the spread of directions $`𝐧_0`$ for the fixed wave number $`k`$ not exceeding very small value $`\mathrm{\Delta }(\mathrm{𝐤𝐧}_0)/\mathrm{\Omega }_e<\delta /\mathrm{\Omega }_e10^9`$. If neutrinos move along radii beyond the neutrinosphere $`r>R_\nu `$ this estimate is too optimistic since there is a spread of the angular distribution of neutrino trajectories that damages simple model with parallel rays assumed here .

Nevertheless, we think the simple dispersion equation (48) based on the enhancement of pure magnetic field perturbations is the only particular case of the general kinetic equations derived here.

Note that even in this approximation there are some advantages of our model comparing with the isotropic plasma case discussed in  as we have just shown in previous section (after eq. (50)).

A more general case with the overlap of electromagnetic eigen modes in a polarized medium and accounting for the spatial dispersion in such plasma seems to be a more realistic model for a magnetized SN while it is beyond of the scope of the present work .

We would like to acknowledge discussions with A.S. Volokitin and A.I. Rez. This work was supported for V.S. by the RFBR grant 00-02-16271.


text:
# Meromorphic Solutions to a Differential–Difference Equation Describing Certain Self-Similar Potentials

## 1. Introduction

Let $`L=_x^2+u(x)`$ be a Schrodinger operator that we factorize as $`L=A^+A^{}+\lambda `$, where $`A^\pm =_x+f(x)`$ and $`x,\lambda `$. Then the function $`f(x)`$ satisfies the Riccati equation $`f^2(x)f^{}(x)+\lambda =u(x)`$. If $`\stackrel{~}{L}=A^{}A^++\lambda `$ denotes a new Schroedinger operator, obtained from $`L`$ by permuting the operator factors $`A^\pm `$, then the potential $`\stackrel{~}{u}(x)`$ of $`\stackrel{~}{L}`$ is given by $`\stackrel{~}{u}(x)=f^2(x)+f^{}(x)+\lambda `$. The differential-difference equation (DDE)

(1)
$$[f(x)+f(x+a)]^{}+f^2(x)f^2(x+a)=\mu ,$$

where $`a`$ and $`\mu `$ are some complex constants, was derived in \[S1, S2\] to describe potentials of the Schrodinger operator $`L`$ satisfying the self-similarity constraint $`\stackrel{~}{L}=ULU^1+\mu `$, where $`U`$ is the translation operator $`Uf(x)=f(x+a)`$. In the case $`\mu =0`$ equation (1) is satisfied by ( \[S1, S2\])

(2)
$$f(x)=\frac{1}{2}\frac{𝒫^{}(xx_0)𝒫^{}(a)}{𝒫(xx_0)𝒫(a)},$$

where $`𝒫`$ denotes the Weierstrass elliptic function and $`x_0`$ is an arbitrary constant. The function (2) is a meromorphic function with only simple poles.

The aim of the present paper is to prove existence of meromorphic solutions with simple poles to (1) in the case $`\mu 0`$. Without any loss of generality, we assume $`a`$ to be a positive real number. Indeed, if $`a=re^{i\varphi }`$, where $`\varphi `$ and $`r>0`$, then the transformation $`xe^{i\varphi }x`$, $`fe^{i\varphi }f`$ reduces (1) to an equation of the same type with the step $`r^+`$. Our approach, which utilizes some ideas of \[To\], consists of essentially three statements: 1) there exist two different formal power series solutions (in powers of $`x^{\frac{1}{2}}`$) to (1); 2) for any formal solution $`\widehat{f}(x)`$ there exists an actual solution $`f(x)`$, analytic (when $`|x|`$ is sufficiently large) in a sector $`S`$ on the complex $`x`$-plane of opening greater than $`\pi `$ and having the asymptotic expansion $`\widehat{f}(x)`$ in $`S`$; 3) any such solution can be meromorphically continued onto $``$ so that $`f(x)`$ may have only first order poles. These statements are proven in Sections 3 and 2 respectively. The author wants express his gratitude to V. Spiridonov for interesting discussions held at the NATO ASI “Special functions 2000”, Tempe, Arizona and for the following correspondence. The author also want to use this opportunity to thank the organizers of the NATO ASI “Special functions 2000”.

## 2. Meromorphic continuation

Let $`R`$ be a simply-connected domain, bounded by piece-wise smooth curves $`\eta _1(\xi )`$ and $`\eta _2(\xi )`$, where $`x=\xi +i\eta `$ is a complex number and $`\eta _2(\xi )>\eta _1(\xi )`$ for all $`\xi (\mathrm{},+\mathrm{})`$. The values $`\eta _2=+\mathrm{}`$ and $`\eta _1=\mathrm{}`$ are allowed. We say that the domain $`R`$ is $`a`$-wide on the interval $`I`$ if $`\eta _2(\xi )\eta _1(\xi )>a`$ for all $`\xi I`$.

###### Theorem 2.1.

A solution $`f(x)`$ to the equation (1) that is analytic in some domain $`R`$ that is $`a`$-wide on $``$ admits a meromorphic continuation on the whole complex plane and all possible singularities of $`f(x)`$ are first order poles.

###### Proof.

Let $`R_1=R`$ and let $`R_n`$ denote the domain bounded by the curves $`\eta _1(\xi )`$ and $`\eta _2(\xi )+(n1)a`$, where $`n=2,3,\mathrm{}`$. We first prove that $`f(x)`$ can be meromorphically continued to the domain $`R_{\mathrm{}}=_1^{\mathrm{}}R_n`$, i.e., to the right of the domain $`R`$, by considering (1) as the Riccati equation

(3)
$$f^{}(x)=f^2(x)+h(x),$$

where $`h(x)=\mu f(xa)^{}f^2(xa)`$. The latter equation is equivalent to to the second order linear differential equation

(4)
$$u^{\prime \prime }(x)+h(x)u(x)=0$$

through the standard transformations

(5)
$$u(x)=\mathrm{exp}_{\stackrel{~}{x}}^xf(t)𝑑t\mathrm{and}f(x)=\frac{u^{}(x)}{u(x)},$$

where $`\stackrel{~}{x}R`$ is a point such that $`\stackrel{~}{x}aR`$. Note that, according to (5), the function $`u(x)`$ is analytic in $`R`$.

Consider (4) in the region $`Q_1`$, where $`Q_n=R_{n+1}\backslash R_n`$, $`n^+`$. According to the assumption of the theorem, the function $`h(x)`$ is analytic in $`Q_1`$, so that solutions of the linear equation (4) are analytic in $`Q_1`$. Thus, we can analytically continue $`u(x)`$ onto $`R_2`$. If $`u(x)`$ does not attain zero value in $`Q_1`$, we get an analytic continuation of $`f(x)`$ onto $`R_2`$ by (5). However, if $`u(x_0)=0`$ for some $`x_0Q_1`$, then $`f(x)`$ has a first order pole in $`x_0`$, which is an isolated singularity of $`f(x)`$. Thus, we obtained the required meromorphic continuation of $`f(x)`$ onto $`R_2`$. This process can be continued to the domains $`R_3,R_4`$, etc. in the same fashion. However, now we have to consider a possibility that $`h(x)`$ has a singularity at $`x_0+a`$. Then, $`u(x)`$ and, correspondingly, $`f(x)`$ may have singularities at $`x_n=x_0+na`$, where $`n^+`$. We need to show that these possible singularities of $`f(x)`$ are first order poles only.

Let

$`h(x)`$ $`=h_0+h_1(xx_0)+h_2(xx_0)^2+\mathrm{}\mathrm{and}`$
(6) $`u(x)`$ $`=(xx_0)+u_2(xx_0)^2+u_3(xx_0)^3+\mathrm{}`$
be the Taylor expansions of $`h(x)`$ and $`u(x)`$ near $`x=x_0`$ in (4). Comparing the like powers of $`xx_0`$ in (4), we obtain

(7)
$$u_2=0,u_3=\frac{h_0}{6},u_4=\frac{h_1}{12},\mathrm{},$$

so that the principle part of $`f(x)`$ at $`x=x_0`$ is $`(xx_0)^1`$.

Combining the expression for $`h(x)`$ with (3), we obtain

(8)
$$h(x+a)=\mu h(x)2f^2(x)=\mu h(x)2[(\mathrm{ln}u(x))^{}]^2.$$

Direct computations show that

(9)
$$[\mathrm{ln}u(x)]^{}=\frac{1}{xx_0}\frac{h_0}{3}(xx_0)\frac{h_1}{4}(xx_0)^2+O(xx_0)^3,$$

and

(10)
$$2([\mathrm{ln}u(x)]^{})^2=\frac{2}{(xx_0)^2}+\frac{4h_0}{3}+h_1(xx_0)+O(xx_0)^2,$$

so that

(11)
$$h_1(x+a)=h(x+a)=\frac{2}{(xx_0)^2}+(\mu +\frac{h_0}{3})+O(xx_0)^2.$$

Considering now this equation near the point $`x_1Q_2`$, we obtain

(12)
$$h_1(x)=\frac{2}{(xx_1)^2}+(\mu +\frac{h_0}{3})+O(xx_1)^2.$$

The proof that all the points $`x_n`$ are either regular points or first order poles of $`f(x)`$ follows by induction from the following two lemmas.

###### Lemma 2.2.

The differential equation

(13)
$$u^{\prime \prime }(x)+h_n(x)u(x)=0,$$

where the coefficient $`h_n(x)`$ has the form

(14)
$$h_n(x)=\frac{N(N1)}{(xx_n)^2}+A_0+A_2(xx_n)^2+\mathrm{}+A_{2N2}(xx_n)^{2N2}+O(xx_n)^{2N1},$$

posseses two linearly independent solutions

(15)
$$u(x)=(xx_n)^N\left[1+u_2(xx_n)^2+\mathrm{}+u_{2N}(xx_n)^{2N}+O(xx_n)^{2N+1}\right]$$

(16)
$$v(x)=(xx_n)^{N+1}\left[1+v_2(xx_n)^2+\mathrm{}+u_{2N}(xx_n)^{2N}+O(xx_n)^{2N+1}\right]$$

that are analytic in a neighborhood of $`x=x_n`$. Here $`A_k,u_k,v_k`$ and $`n,N`$ with $`n+1N2`$.

###### Proof.

Frobenius multipliers of (13) at $`x=x_n`$ are $`N`$ and $`N+1`$, so the first terms of $`u(x)`$ and $`v(x)`$ are $`(xx_n)^N`$ and $`(xx_n)^{N+1}`$ respectively. Suppose,

(17)
$$u(x)=(xx_n)^N+\underset{k=1}{\overset{\mathrm{}}{}}u_k(xx_n)^{N+k}.$$

Computing $`u_k`$, we see that the odd coefficients $`u_1=u_3=\mathrm{}=u_{2N1}=0`$, so that $`u(x)`$ is in the form (15). Similar arguments work for the second solution $`v(x)`$. ∎

Arguments of Lemma 2.2 show that solutions to (13) have no branching at $`x=x_n`$ since Frobenius multipliers are integer and there are no logarithms. Note that the general solution to (13) has the form

(18)
$$u(x)=(xx_n)^{N+1}\left[u_0+u_2(xx_n)^2+\mathrm{}+u_{2N2}(xx_n)^{2N2}+O(xx_n)^{2N1}\right]$$

where $`u_00`$. The only nontrivial special solution to (13) is proportional to $`u(x)`$ given by (15). In any case, the point $`x_n`$ is a simple pole of the solution $`f(x)=[\mathrm{ln}u(x)]^{}`$ to (3). The following Lemma shows that the new coefficient

(19)
$$h_{n+1}(x)=\mu h_n(xa)2[(\mathrm{ln}u(xa))^{}]^2,$$

where $`u(x)`$ is a solution to (13), has the form (14) with another $`N`$. Thus, the point $`x_{n+1}`$ is also a simple pole or a regular point of the solution $`f(x)`$. (Note than $`x_{n+1}`$ is a regular point of $`h_{n+1}`$ if the new $`N+1`$.)

###### Lemma 2.3.

The new coefficient $`h_{n+1}(x)`$ has the form (14) with the new $`N`$ equal to $`N1`$ if $`u(x)`$ is given by (18) or equal to $`N+1`$ if $`u(x)`$ is given by (15).

###### Proof.

Consider, for example, the case when $`u(x)`$ is given by (15). Using expansion (14) for $`h_n(x)`$, we obtain the first odd coeficient $`u_{2N+1}`$ in (15) is

(20)
$$u_{2N+1}=\frac{A_{2N+1}}{4N(2N+1)}.$$

On the other hand, we see that $`\mathrm{ln}u(x)=N\mathrm{ln}(xx_n)+u_2(xx_n)^2+\mathrm{}`$, where the leading odd term of $`\mathrm{ln}u(x)`$ is $`u_{2N+1}(xx_n)^{2N+1}`$. Then

$$[\mathrm{ln}u(x)]^{}=(xx_n)^1\left[N+2u_2(xx_n)^2+\mathrm{}\right]$$

where the leading odd term in the square brackets is $`(2N+1)u_{2N+1}(xx_n)^{2N+1}`$. Finally, we get the leading odd term of $`2([\mathrm{ln}u(x)]^{})^2`$ as $`4N(2N+1)u_{2N+1}(xx_n)^{2N1}`$. So, according to (19) and (20), the leading odd term of $`h_{n+1}(x)`$ at $`x_{n+1}`$ is of the order $`(xx_{n+1})^{2N+1}`$. The leading term of $`h_{n+1}(x)`$ is $`\frac{N(N1)2N^2}{(xx_{n+1})^2}=\frac{(N+1)N}{(xx_{n+1})^2}`$. So, $`h_{n+1}(x)`$ has the form (14) with the new $`N`$ equal to $`N+1`$. The case when $`u(x)`$ is given by (18) can be considered in a similar way. ∎

Consider, for example, the singular point $`x_1`$. According to (12), we have $`N=2`$ at this point. Solving the corresponding initial value problem for (13), we get the solution $`u(x)`$ of either type (15) or (18). In any case, $`f(x)`$ has a first order pole. The corresponding function $`h_2(x)`$, according to Lemma 2.3, has a second order pole at $`x=x_2`$ with the principal part $`\frac{6}{(xx_2)^2}`$ if $`u(x)`$ is proportional to (15), and is regular at $`x=x_2`$ if $`u(x)`$ is given by (18). We can continue these arguments to show that all points $`x_n`$ are either first order poles or regular points. Thus, we proved meromorphic continuation of $`f(x)`$ on $`R_{\mathrm{}}`$.

To prove the meromorphic continuation to the left of the domain $`R`$ let us note that the transformation $`x=ta,g(t)=f(t)`$ reduces the equation (1) to the equation of the same type

(21)
$$[g(t)+g(t+a)]^{}+g^2(t)g^2(t+a)=\mu ,$$

which have an analytic solution $`g(t)`$ on $`tR_{\mathrm{}}`$. We can now use the previous arguments to continue $`g(t)`$ to the right on the whole complex plane. ∎

###### Corollary 2.4.

Let $`I`$ be an interval and let $`f(x)`$ be a solution to the equation (1) that is analytic in a domain $`R`$ that is $`a`$-wide on $`I`$. Then $`f(x)`$ admits a meromorphic continuation on the strip $`\mathrm{}xI`$ and all possible singularities of $`f(x)`$ are first order poles.

## 3. Asymptotic solutions

In this section we rewrite (1) as

(22)
$$\mathrm{}f^{}(z)=\mu +\mathrm{\Delta }f^2(z),$$

where $`z=x+b,b=\frac{a}{2}`$ and the operators $`\mathrm{},\mathrm{\Delta }`$ act on a function $`g(z)`$ as

(23)
$$\mathrm{}g(z)=g(z+b)+g(zb)\mathrm{\Delta }g(z)=g(z+b)g(zb).$$

Some important for us properties of $`\mathrm{},\mathrm{\Delta }`$ are given by identities

(24)
$$2\mathrm{}[fg]=\mathrm{}f\mathrm{}g+\mathrm{\Delta }f\mathrm{\Delta }g2\mathrm{\Delta }[fg]=\mathrm{\Delta }f\mathrm{}g+\mathrm{}f\mathrm{\Delta }g,$$

where $`f,g`$ are given functions.

In this section we construct two different formal power series solutions $`\widehat{f}^\pm (z)`$ to (1) and prove existence of corresponding actual solutions $`f_n^\pm (z)`$. These solutions are analytic in the corresponding sectors $`S_n^\pm ,n=0,1,2`$, specified below, when $`|z|`$ is sufficiently large and admit asymptotic expansions

(25)
$$f_n^\pm (z)\widehat{f}^\pm (z),z\mathrm{},zS_n^\pm .$$

This fact together with Theorem 2.1 prove existence of non-trivial meromorphic solutions to (1).

###### Proposition 3.1.

Equation (22) possesses two formal power series solutions

(26)
$$\widehat{f}_\pm (z)=\pm \lambda \sqrt{z}+\frac{1}{2b}+\underset{k=2}{\overset{\mathrm{}}{}}y_k^\pm z^{\frac{k}{2}},$$

where $`\lambda =\sqrt{\frac{\mu }{2b}}`$ and the coefficients $`y_k^\pm `$ are defined uniquely.

###### Proof.

The substitution $`f=\lambda z^{\frac{1}{2}}+c+\widehat{y}(z)`$ reduces the equation (22) to

(27)
$$\frac{\lambda }{2}\mathrm{}z^{\frac{1}{2}}+\mathrm{}\widehat{y}^{}=\mu +(\lambda \mathrm{\Delta }z^{\frac{1}{2}}+\mathrm{\Delta }\widehat{y})(\lambda \mathrm{}z^{\frac{1}{2}}+2c+\mathrm{}\widehat{y}).$$

Comparing leading coefficients and taking into account $`\mathrm{\Delta }z^{\frac{1}{2}}\mathrm{}z^{\frac{1}{2}}=2b`$, we get $`\lambda =\pm \sqrt{\frac{\mu }{2b}}`$. Taking into account (24), the latter equation can be now rewritten as

(28)
$$2\lambda \mathrm{\Delta }[z^{\frac{1}{2}}\widehat{y}]=\mathrm{}\widehat{y}^{}\mathrm{\Delta }\widehat{y}^22c\mathrm{\Delta }\widehat{y}+\lambda (\frac{1}{2}\mathrm{}z^{\frac{1}{2}}2c\mathrm{\Delta }z^{\frac{1}{2}}).$$

To expand $`\widehat{y}`$ in powers of $`z^{\frac{1}{2}}`$ we need the free term of (28) to be of order $`O(z^1)`$ or less. Thus, we obtain $`c=\frac{1}{2b}`$, so that the free term is of the order $`O(z^{\frac{3}{2}})`$. Equation (28) can be now rewritten as

(29)
$$\mathrm{\Delta }[z^{\frac{1}{2}}\widehat{y}]=\frac{1}{2\lambda }\left(\mathrm{}\widehat{y}^{}\mathrm{\Delta }\widehat{y}^2\frac{1}{b}\mathrm{\Delta }\widehat{y}\right)+\frac{\mathrm{}z^{\frac{1}{2}}}{4}\frac{\mathrm{\Delta }z^{\frac{1}{2}}}{2b}.$$

It is clear that the expression in the left hand side is the dominant term of the latter equation and that the substitution $`\widehat{y}=_{k=2}^{\mathrm{}}y_k^\pm z^{\frac{k}{2}}`$, satisfying (29), defines the coefficients $`y_k`$ uniquely. ∎

Let $`\widehat{f}(z)=\lambda \sqrt{z}+\frac{1}{2b}+\mathrm{}`$ denote one of the formal solutions $`\widehat{f}^\pm (z)`$. Given $`\lambda `$, we define angles

(30)
$$\beta _n=\frac{2}{3}\pi (1+2n)\frac{2}{3}\mathrm{arg}\lambda ,$$

where $`n=0,1,2`$ and $`\mathrm{arg}\lambda [0,2\pi )`$. If $`\beta _n+\pi `$ is not a multiple of $`2\pi `$, we define the sector $`S_n^+`$ on the Riemann surface of $`z^{1/3}`$ by extending a small sector bisected by $`\mathrm{arg}z=\beta _n`$ independently in both positive and negative directions until it either hits the negative real direction or the ray $`\mathrm{arg}z=\beta _n\pm \pi `$ respectively. If $`\beta _n`$ is not a multiple of $`2\pi `$, the sector $`S_n^{}`$ on the Riemann surface of $`z^{1/3}`$ is defined by extending a small sector bisected by $`\mathrm{arg}z=\beta _n`$ independently in both positive and negative directions until it either hits the positive real direction or the ray $`\mathrm{arg}z=\beta _n\pm \pi `$ respectively. In the case $`\beta _n=2\pi k,k`$ or $`\beta _n=(2k+1)\pi ,k`$, the sectors $`S_n^{}`$ or $`S_n^+`$ respectively are considered to be empty.

###### Theorem 3.2.

If $`\widehat{f}(z)`$ is a formal solution to the equation (1) and $`S_n^+`$ is a nonempty sector described above, then there exists and actual solution $`f(z)`$ to (1) that is analytic in sufficiently remote part of any proper subsector $`S`$ of $`S_n^+`$ and

(31)
$$f(z)\widehat{f}(z),z\mathrm{},zS.$$

###### Proof.

Our main idea is to reduce the considered DDE to an integro-differential equation (IDE) and to show that the latter equation can be solved by successive iterations in a proper sectorial neighborhood of infiniti. Substituting $`f(z)=\lambda \sqrt{z}+\frac{1}{2b}+y(z)`$, we obtain equation (29) for $`y(z)`$. The inverse Laplace transform $`^1`$, applied to (29), yields (see, for example, \[Ob\])

$`2\mathrm{sinh}(pb)^1[z^{\frac{1}{2}}y](p)`$ $`={\displaystyle \frac{1}{2\lambda }}\left(2p\mathrm{cosh}(pb)Y(p)+2\mathrm{sinh}(bp)\left[Y^2(p)+{\displaystyle \frac{1}{b}}Y(p)\right]\right)`$
(32) $`+{\displaystyle \frac{\mathrm{sinh}bp}{2\sqrt{\pi p}}}(\mathrm{coth}bp{\displaystyle \frac{1}{bp}}),`$
where $`Y(p)=^1[y](p)`$ and $`F(p)^2=_0^pF(pq)F(q)𝑑q`$. After a simple algebra, the latter equation becomes

(33)
$$^1[z^{\frac{1}{2}}y](p)=\frac{\mathrm{\Xi }(bp)}{2\lambda b}Y(p)\frac{1}{2\lambda }Y^2(p)\frac{\mathrm{\Xi }(bp)}{2\sqrt{\pi }bp^{3/2}},$$

where $`\mathrm{\Xi }(x)=x\mathrm{coth}x1`$. Note that $`\mathrm{\Xi }(x)`$ is a meromorphic function with simple poles at the points $`i\pi k`$, where $`k\backslash \{0\}`$, and that $`\mathrm{\Xi }(x)`$ has not more than linear growth in any non-vertical direction $`\mathrm{arg}x=const`$.

Separating the linear part $`bp`$ of $`\mathrm{\Xi }(bp)`$ along the positive real axis, we can rewrite (33) as

(34)
$$^1[z^{\frac{1}{2}}y](p)\frac{p}{2\lambda }Y(p)=\frac{\mathrm{\Xi }(bp)bp}{2\lambda b}Y(p)\frac{1}{2\lambda }Y^2(p)\frac{\mathrm{\Xi }(bp)}{2\sqrt{\pi }bp^{3/2}},$$

where the function $`\mathrm{\Xi }(bp)bp`$ is bounded on $`[0,\mathrm{})`$. Applying now the Laplace transform $``$ to (34), we reduce (29) to the IDE

(35)
$$y^{}(z)+2\lambda z^{\frac{1}{2}}y(z)=\frac{1}{b}\xi (z)y(z)y^2(z)\frac{\lambda }{b\sqrt{\pi }}r(z),$$

where $`\xi (z)=[\mathrm{\Xi }(bp)bp](z)`$, $`r(z)=[\mathrm{\Xi }(bp)p^{3/2}](z)`$ and the convolution is defined by

(36)
$$f(z)g(z)=\frac{1}{2\pi i}_{Ai\mathrm{}}^{A+i\mathrm{}}f(s)g(zs)𝑑s$$

with a sufficiently large $`A>0`$. Considering (35) as a perturbed linear equation $`y^{}(z)+2\lambda z^{\frac{1}{2}}y(z)=0`$, we rewrite the former as

(37)
$$y(z)=e^{\frac{4}{3}\lambda z^{3/2}}_{\gamma (z)}e^{\frac{4}{3}\lambda t^{3/2}}W[y(t)]𝑑t,$$

where the nonlinear operator $`W[y]`$ denotes the right hand side of (35) and the contour of integration $`\gamma (z)`$ is to be specified below. Equation (37) can be rewritten in the operator form as $`y=W[y]`$.

Let us assume for a while that $`\frac{\pi }{4}<\mathrm{arg}\lambda <\frac{7\pi }{4}`$ and consider sector $`S_0^+`$. According to (30), this choice of $`\lambda `$ allows us to find a proper closed subsector $`SS_0^+`$ that contains the right half-plane. We want to solve the IDE (35) in a sufficiently remote part of the sector $`S`$ by successive approximations. In order to formulate the statement more precisely we need to introduce the following notations.

Let $`\mathrm{\Sigma }`$ be the image of the sector $`S`$ under the transformation $`\zeta =z^{3/2}`$. Let $`z_0`$ be a sufficiently remote point on the ray $`\mathrm{arg}z=\beta _0`$ and let $`\mathrm{\Sigma }_{\zeta _0}`$, where $`\zeta _0=z_0^{3/2}`$ denote the parallel shift of $`\mathrm{\Sigma }`$ so that the vertex of $`\mathrm{\Sigma }`$ is shifted to $`\zeta _0`$. For every $`\zeta \mathrm{\Sigma }_{\zeta _0}`$ we define a contour $`\mathrm{\Gamma }(\zeta )`$ as a ray eminating from $`\zeta `$ and such that: $`\mathrm{\Gamma }(\zeta )=\{\tau :\mathrm{arg}(\tau \zeta )=\frac{3}{2}\beta _0\}`$ if $`|\mathrm{arg}\zeta \frac{3}{2}\beta _0|<\frac{\pi }{2}`$; $`\mathrm{\Gamma }(\zeta )=\{\tau :\mathrm{arg}(\tau \zeta )=\mathrm{arg}\zeta +\frac{\pi }{2}\}`$ if $`\mathrm{arg}\zeta \frac{3}{2}\beta _0<\frac{\pi }{2}`$ and; $`\mathrm{\Gamma }(\zeta )=\{\tau :\mathrm{arg}(\tau \zeta )=\mathrm{arg}\zeta \frac{\pi }{2}\}`$ if $`\mathrm{arg}\zeta \frac{3}{2}\beta _0>\frac{\pi }{2}`$. Then the region $`S_{z_0}`$ and the contour $`\gamma (z)`$ are the images of $`\mathrm{\Sigma }_{\zeta _0}`$ and $`\mathrm{\Gamma }(\zeta )`$ under the transformation $`z=\zeta ^{2/3}`$.

Let $`y_0(z)0`$, $`y_k(z)=W[y_{k1}](z),k=1,2,\mathrm{}`$ and $`\delta y_k=y_ky_{k1}`$. We will show that the solution to (37) is given by

(38)
$$y(z)=\underset{k=1}{\overset{\mathrm{}}{}}\delta y_k,$$

where the series converges absolutely and uniformly in $`S_{z_0}`$ for sufficiently large $`|z_0|`$. This can be done by introducing the Banach space $``$ of functions $`h(z)`$, such that $`h`$ is analytic on $`S_{z_0}`$ and satisfy $`|h(z)|B|z^2|`$ there with some constant $`B>0`$. (Note that $`B`$ depends on $`h`$.) The norm of $`h`$ is $`sup_{zS_{z_0}}|z^2h(z)|`$. According to Lemma 14.2 from \[Wa\], the integral operator $`:`$ is a bounded linear operator, where $``$ is proportional to $`|z_0^{\frac{1}{2}}|`$. We start to study the nonlinear operator $`𝒲[y]`$ by considering the convolution

$`\xi (z)y(z)`$ $`=[(\mathrm{\Xi }(bp)bp)Y(p)](z)=[\{bp(\mathrm{coth}bp1)e^{2bp}\}Y(p)](z)`$
(39) $`+[(e^{2bp}1)Y(p)](z)=\mu (z)y(z)+y(z+2b)y(z),`$
where, according to \[OB\],

(40)
$$\mu (z)=[\{bp(\mathrm{coth}bp1)e^{2bp}]=\frac{1}{2b}\mathrm{\Psi }^{}(1+\frac{z}{2b})\frac{1}{z+2b}.$$

Here $`\mathrm{\Psi }(x)=\frac{\mathrm{\Gamma }^{}(x)}{\mathrm{\Gamma }(x)}`$ denotes the logarithmic derivative of the Gamma-function. It is clear that $`\mu (z)`$ is a meromorphic function with double poles at $`z=2kb,k^+`$. The asymptotic expansion

(41)
$$\mathrm{\Psi }^{}(x)\frac{1}{z}+\frac{1}{2z^2}+\underset{k=1}{\overset{\mathrm{}}{}}\frac{B_{2k}}{z^{2k+1}},z\mathrm{},|\mathrm{arg}z|<\pi ,$$

where $`B_{2k}`$ are Bernoulli numbers, follows from the Stirling formula (see \[GR\]). Combining the latter facts, we obtain the estimate

(42)
$$|\mu (z)|\frac{b\mathrm{\Psi }_0}{|z+2b|^2},zS,$$

where the constant $`\mathrm{\Psi }_0`$ depends only on the sector $`S`$.

###### Proposition 3.3.

If $`y`$ then $`\mu (z)y(z)`$ and

(43)
$$\frac{1}{b}\mu yMy,$$

where $`M>0`$ does not depend on $`y`$ and $`z_0`$.

###### Proof.

Consider first the case $`z`$ belongs to the right half-plane $`\widehat{S}=\{z:\mathrm{}zmz_0\}`$, where the number $`m>0`$ is choosen so that $`\widehat{S}S_{z_0}`$. Note that $`m`$ does not depend on $`|z_0|`$. Setting $`A=mz_0`$ in (36), we obtain

(44)
$$\mu (z)y(z)=\frac{1}{2\pi i}_{Ai\mathrm{}}^{A+i\mathrm{}}y(s)\mu (zs)𝑑s$$

Note that we are integrating over the boundary of $`\widehat{S}`$ and thus $`\mathrm{}(zs)0`$. The fact that $`\mu (z)y(z)`$ is analytic in $`\widehat{S}`$ follows from the properties of $`\mu (z)`$ and $`y(z)`$ immediately. Utilizing (42) and the fact that $`y`$, we obtain

$`{\displaystyle \frac{1}{by}}|\mu (z)y(z)|`$ $`{\displaystyle \frac{\mathrm{\Psi }_0}{2\pi }}{\displaystyle _{Ai\mathrm{}}^{A+i\mathrm{}}}{\displaystyle \frac{|ds|}{|s|^2|z+2bs|^2}}`$
(45) $`={\displaystyle \frac{\mathrm{\Psi }_0}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\eta }{(A^2+\eta ^2)[(uA)^2+(v\eta )^2]}},`$
where $`s=A+i\eta `$ and $`z+2b=u+iv`$ with $`\eta ,u,v`$.

Let $`I`$ be the latter integral and $`g(\eta )`$ denote its integrand, which has simple poles at the points $`\eta =iA`$ and $`\eta =v+i(uA)`$ in the upper half-plane. Computing $`I`$ via the residues of $`g(\eta )`$ in the upper half-plane, we obtain

$`I=2\pi i{\displaystyle \text{Res }g(\eta )}`$ $`={\displaystyle \frac{\pi }{v+i(u2a)}}\left[{\displaystyle \frac{1}{A(viu)}}+{\displaystyle \frac{1}{(uA)(v+iu)}}\right]`$
(46) $`={\displaystyle \frac{\pi }{A|z+2b|^2}}\left(1+{\displaystyle \frac{A}{uA}}\right)`$
Thus

(47)
$$\frac{1}{by}|\mu (z)y(z)|\frac{\mathrm{\Psi }_0}{2b|z|^2}$$

since $`A`$ can be taken greater than $`2b`$.

Consider now the half-plane $`\widehat{S}_\varphi `$ that is obtained by rotating $`\widehat{S}`$ on the angle $`\varphi `$, where $`\varphi `$ is choosen in such a way that $`\widehat{S}_\varphi S_{z_0}`$. Define now another convolution $`m_\varphi y`$ by (44), where the contour of integration is the boundary of $`\widehat{S}_\varphi `$. Using the same arguments as above, we can obtain the estimate

(48)
$$\frac{1}{by}|\mu (z)_\varphi y(z)|\frac{M_\varphi }{|z|^2},$$

where $`z\widehat{S}_\varphi `$ and $`M_\varphi `$ continuously depends on $`\varphi `$. However, $`m(z)_\varphi y(z)`$ is an analytic continuation of $`m(z)y(z)`$, since the functions coincide on $`\widehat{S}_\varphi \widehat{S}`$. Taking $`M=\mathrm{max}M_\varphi `$, where $`M_0=\frac{\mathrm{\Psi }_0}{2b}`$, we complete the proof of the proposition. ∎

To complete the proof of the theorem we use the standard technique to show the convergence of iterations (38). Using properties of the Laplace transform and of the operator $``$, one can show that $`y_1`$. Let $`K=y_1`$, and let us prove by induction that $`\delta y_n2^{1n}K`$ if $`|z_0|`$ is sufficiently large. Indeed, according to the estimate of the convolution,

(49)
$$\delta y_n\left(M+2+\frac{4K}{|z|^2}\right)\delta y_{n1},$$

where the induction assumption

(50)
$$y_{n1}^2y_{n2}^2y_{n1}+y_{n_2}\frac{|\delta y_{n1}}{|z|^2}4K\frac{|\delta y_{n1}}{|z|^2}$$

was used to estimate the nonlinear term of $`𝒲`$. It remains to choose $`z_0`$ so that $`\left(M+2+\frac{4K}{|z|^2}\right)<\frac{1}{2}`$ for $`zS_{z_0}`$ to complete the proof of the theorem for $`\frac{\pi }{4}<\mathrm{arg}\lambda <\frac{7\pi }{4}`$ the sector $`S_0^+`$.

Let us now consider the general case $`0\mathrm{arg}\lambda <2\pi `$. The sector $`S_0^+`$ is given by $`\frac{\pi +2\mathrm{arg}\lambda }{3}<\mathrm{arg}z<\pi `$ if $`\mathrm{arg}\lambda [0,\pi ]`$ and by $`\pi <\mathrm{arg}z<\frac{\pi 2\mathrm{arg}\lambda }{3}`$ if $`\mathrm{arg}\lambda [\pi ,2\pi ]`$. Note that the opening of the sector $`S_0^+`$ is greater than $`\pi `$ for any $`\mathrm{arg}\lambda `$ and that $`S_0^+`$ contains the right hlaf-plane if $`\frac{\pi }{4}<\mathrm{arg}\lambda <\frac{7\pi }{4}`$. Let us now choose a proper closed subsector $`SS_0^+`$ of opening greater than $`\pi `$, and let $`\mathrm{arg}z=\alpha `$ be the bisector of $`S`$.Clearly, $`|\alpha |<\frac{\pi }{2}`$. Let

(51)
$$_\alpha [Y](z)=_0^{e^{i\alpha }\mathrm{}}e^{zp}Y(p)𝑑p$$

define the Laplace transform along the ray $`\mathrm{arg}p=\alpha `$. The contour for the corresponding inverse Laplace transform as well as for the corresponding convolution $`_\alpha `$ is a straight line perpendicular to $`\mathrm{arg}z=\alpha `$. Therefore, we can use our previous arguments to show the uniform and absolute convergence of iterations (38) in a properly constructed $`S_{z_0}S`$. In the same fashion the theorem can be proven for any nonempty sector $`S_n^+`$, $`n`$.

Recall that the function $`\mathrm{\Xi }(bp)`$ has poles on the imaginary axis and, therefore, $`_\alpha \mathrm{\Xi }`$ is not defined for $`\alpha =\pm \frac{\pi }{2}`$. However, we can define $`_\alpha \mathrm{\Xi }`$ for $`\alpha `$ such that $`|\alpha \pi |<\frac{\pi }{2}`$ and repeat the previous arguments for sectors $`S_n^{}`$, $`n`$.

###### Corollary 3.4.

Let $`f_n^\pm ,n=0,1,2`$ denote a solution of (1), analytic in some remote part of $`S_n^\pm `$ as described in 3.2. Then, according to 2.1, $`f_n^\pm `$ is a meromorphic function on $``$ that can have only simple poles.


text:
# Giant strongly connected component of directed networks

## Abstract

We describe how to calculate the sizes of all giant connected components of a directed graph, including the strongly connected one. Just to the class of directed networks, in particular, belongs the World Wide Web. The results are obtained for graphs with statistically uncorrelated vertices and an arbitrary joint in,out-degree distribution $`P(k_i,k_o)`$. We show that if $`P(k_i,k_o)`$ does not factorize, the relative size of the giant strongly connected component deviates from the product of the relative sizes of the giant in- and out-components. The calculations of the relative sizes of all the giant components are demonstrated using the simplest examples. We explain that the giant strongly connected component may be less resilient to random damage than the giant weakly connected one.

05.10.-a, 05-40.-a, 05-50.+q, 87.18.Sn

The giant components of a network are components which relative sizes are finite (nonzero) in the large network limit. The knowledge of these sizes provides the basic information about the global topology of a network. The understanding of the topological structure of networks and its change under external action is the central problem of the statistical physics of random networks . Actually, this is the natural generalization of the general percolation theory.

The most interesting networks in Nature, including the World Wide Web, are directed graphs, i.e., their vertices are connected by directed edges . In the general case, the structure of the directed graph looks as it is shown in Fig. 1 (all the notions are introduced and explained in the figure caption). In particular, the World Wide Web has such a structure .

In Refs. , the previous strong results of mathematicians  were developed, and it was proposed the general theory of percolation phenomena in networks with arbitrary degree distributions and statistically uncorrelated (randomly connected) vertices. Of course, the last assumption is not true for most of growing nets in Nature. Nevertheless, the direct conclusions of such an approach proved to explain the behavior of real networks .

In paper , it was shown how to find the relative sizes of the following giant components of directed graphs with statistically uncorrelated vertices: (i) the giant weakly connected component, $`W`$; (ii) the giant in-component, $`I`$; and (iii) the giant out-component, $`O`$. We emphasize that, for brevity and consistency, we use the definitions of the giant in- and out-components other than in Refs.  (see the caption of Fig. 1).

Here we demonstrate how to calculate the relative size $`S`$ of, perhaps, the most important part of the directed graph, of its giant strongly connected component (iv). (In the GSSC, every pair of vertices is connected in both directions, i.e., from one of the vertices, one can approach the other by moving either along or against the edge directions.) This allows us to completely describe the total structure of directed graphs with arbitrary degree distributions and statistically uncorrelated vertices. For the demonstration, we use the networks with the simplest degree distributions providing non-trivial results.

We have to briefly remind a very usefull approach of Ref. . The Z-transforms (or generating functions) are used. For the undirected graph, $`\mathrm{\Phi }(x)_kP(k)x^k`$, and, for the directed one, $`\mathrm{\Phi }(x,y)_{k_i,k_o}P(k_i,k_o)x^{k_i}y^{k_o}`$ . Here, $`P(k)P^{(w)}(k)=_{k_i}P(k_i,kk_i)`$ is the degree distribution ($`k=k_i+k_o`$ is the total number of connections of a node) and $`P(k_i,k_o)`$ is the joint distribution of in- and out-degrees. When all the connections are inside the network, the average in- and out-degrees are equal: $`_x\mathrm{\Phi }(x,1)|_{x=1}=_y\mathrm{\Phi }(1,y)|_{y=1}z^{(d)}`$. Therefore the average degree is $`z=2z^{(d)}`$. If one ignores the directedness of edges, the degree distribution of the directed network, in the Z-representation, takes the form $`\mathrm{\Phi }^{(w)}(x)=\mathrm{\Phi }(x,x)`$. In this case, the distribution of the number of connections minus one of any of the end vertices of a randomly chosen edge corresponds to $`\mathrm{\Phi }_1^{(w)}(x)\mathrm{\Phi }^{(w)}(x)/z`$.

The giant weakly connected component exists if $`\mathrm{\Phi }_1^{(w)}(1)>1`$, that corresponds to the well known criterium of Molloy and Reed

$$\underset{k}{}k(k2)P(k)>0.$$
(1)
The size of the GWCC, $`W`$, can be easily obtained from the relations

$$W=1\mathrm{\Phi }^{(w)}(t_c),x_c=\mathrm{\Phi }_1^{(w)}(t_c).$$
(2)
From Eq. (1), one sees that the existence of the GWCC crucially depends on the size of the fraction of dead ends in the network. Indeed, $`P(1)`$ is the only term in Eq. (1) that prevents the GWCC. In Fig. 2,a, the evolution of the giant connected component of the undirected graph induced by the change of some control parameter is schematically shown. From Eq. (1), it is also clear that the divergency of the second moment of the degree distribution makes the GWCC extremely stable . If the exponent $`\gamma `$ of the power-law degree distribution $`P(k)k^\gamma `$ is smaller than $`3`$ or equal to it, one has to remove at random almost all the vertices or edges of the network to eliminate the GWCC .

In a similar way, it is easy to study the GIN and GOUT componets of the directed network . One introduces the Z-transform of the out-degree distribution of the vertex approachable by following a randomly chosen edge when one moves along the edge direction, $`\mathrm{\Phi }_1^{(o)}(y)_x\mathrm{\Phi }(x,y)|_{x=1}/z^{(d)}`$. Also, $`\mathrm{\Phi }_1^{(i)}(x)_y\mathrm{\Phi }(x,y)|_{y=1}/z^{(d)}`$ corresponds to the in-degree distribution of the vertex which one can approach moving against the edge direction. The GIN and GOUT are present if $`\mathrm{\Phi }_1^{(i)}(1)=\mathrm{\Phi }_1^{(o)}(1)=_{xy}^2\mathrm{\Phi }(x,y)|_{x=1,y=1}/z^{(d)}>1`$, that is

$`{\displaystyle \underset{k_i,k_o}{}}(2k_ik_ok_ik_o)P(k_i,k_o)=`$ (3)
$`2{\displaystyle \underset{k_i,k_o}{}}k_i(k_o1)P(k_i,k_o)=2{\displaystyle \underset{k_i,k_o}{}}k_o(k_i1)P(k_i,k_o)>0.`$ (4)
(5)
In this case, there exist the non-trivial solutions of the equations

$$x_c=\mathrm{\Phi }_1^{(i)}(x_c),y_c=\mathrm{\Phi }_1^{(o)}(y_c).$$
(6)
They have the following meanings. $`x_c<1`$ is the probability that the connected component, obtained by moving against the edge directions starting from a randomly chosen edge, is finite. $`y_c<1`$ is the probability that the connected component, obtained by moving along the edge directions starting from a randomly chosen edge, is finite. Then $`P(k_i,k_o)x_c^{k_i}`$ and $`P(k_i,k_o)y_c^{k_o}`$ are the probabilities that a vertex with $`k_i`$ incoming and $`k_o`$ outgoing edges has finite in- and out-components, respectively. The in- and out-components of a vertex are sets of vertices which are approachable from it moving against and along edges, respectively, plus the vertex itself. Summation of these expressions over $`(k_i,k_o)`$ yields the total probabilities that the in- and out-components of a randomly chosen vertex are finite, respectively. Therefore, they are equal to $`\mathrm{\Phi }(x_c,1)`$ and $`\mathrm{\Phi }(1,y_c)`$, respectively. Thus, the relative sizes of the GIN and GOUT are

$$I=1\mathrm{\Phi }(x_c,1),O=1\mathrm{\Phi }(1,y_c).$$
(7)
Here we show that from Eq. (6), it is possible to find not only $`I`$ and $`O`$ but also the relative size $`S`$ of the GSCC using the considerations similar to Ref. . Suppose that a vertex has $`k_i`$ incoming and $`k_o`$ outgoing edges. They are assumed to be statistically independent. Then the probability that all the incoming edges come from finite in-components is $`x_c^{k_i}`$. The probability that this vertex has the infinite in-component is equal to $`1x_c^{k_i}`$, that is, at least one of the $`k_i`$ incoming edges has to come from the GIN. Similarly, $`1y_c^{k_o}`$ is the probability that the vertex has the infinite out-component. The vertex belongs to the GSCC if its in- and out-components are both infinite; the corresponding probability is equal to $`(1x_c^{k_i})(1y_c^{k_o})`$. Then the total probability that a vertex belongs to the GSCC is equal to $`_{k_i,k_o}P(k_i,k_o)(1x_c^{k_i})(1y_c^{k_o})`$. Finally, the relative size of the GSCC takes the form

$`S`$ $`=`$ $`{\displaystyle \underset{k_i,k_o}{}}P(k_i,k_o)(1x_c^{k_i})(1y_c^{k_o})=`$ (9)
$`1\mathrm{\Phi }(x_c,1)\mathrm{\Phi }(1,y_c)+\mathrm{\Phi }(x_c,y_c).`$
Therefore, $`\mathrm{\Phi }(x_c,y_c)`$ is equal to the probability that both the in- and out-components of a vertex are finite. One can write $`\mathrm{\Phi }(x_c,y_c)=1W+T`$. Knowing $`W`$, $`S`$, $`I`$, and $`O`$, it is easy to obtain the relative size of tendrils,

$$T=W+SIO.$$
(10)
Eqs. (2), (6)–(10) allow us to obtain all the giant components of the directed networks with arbitrary joint in,out-degree distributions and statistically uncorrelated vertices. It is usefull to rewrite the main Eq. (9) in the form

$$S=IO+\mathrm{\Phi }(x_c,y_c)\mathrm{\Phi }(x_c,1)\mathrm{\Phi }(1,y_c).$$
(11)
If the the joint distribution of in- and out-degrees factorizes, $`P(k_i,k_o)=P^{(i)}(k_i)P^{(o)}(k_o)`$, Eq. (11) takes the simple form $`S=IO`$. Otherwise, such factorization of $`S`$ is impossible. At the threshold, $`x_c=y_c=1`$, and $`I`$, $`O`$, and $`S`$ simultaneously approach zero.

We have no intention to calculate the sizes of the giant components for real networks, for instance, for the WWW, for the following reasons. There are some correlations between their vertices, and their joint degree distributions are unknown yet (nevertheless, see the attempt of the calculation of $`I`$ and $`O`$ for the WWW in Ref. ). Instead of this, for demonstration, we consider two simplest non-trivial nets. In the first of them, the joint in,out-degree distribution factorizes: $`P(k_i,k_o)=[p(\delta _{k_i,0}+\delta _{k_i,1})+(12p)\delta _{k_i,3}]\times `$ $`[p(\delta _{k_o,0}+\delta _{k_o,1})+(12p)\delta _{k_o,3}]`$. The results, the dependence of the sizes of the giant components on $`p`$ are shown in Fig. 3 and, schematically, in Fig. 2,b. The curves $`I(p)=O(p)`$ approach the threshold linearly, and $`S(p)`$ – quadratically but the range of the quadratic dependence is narrow. Over a wide range of $`p`$, $`S(p)1/2p`$ (also see the next case, Fig. 4).

In real growing nets, the joint in,out-degree distributions do not factorize just because of their growth . Therefore, for comparison, we calculate the sizes of the giant components for the network with the distribution $`P(k_i,k_o)=p(\delta _{k_i,0}\delta _{k_o,1}+\delta _{k_i,1}\delta _{k_o,0})+(12p)\delta _{k_i,3}\delta _{k_o,3}`$. This means that large in- and out-degrees correlate as well as small in- and out-degrees. The results are plotted in Fig. 4 (also see the schematic plot 2,b) . One sees that, in this case, the size of the GSCC noticeably differs from the product $`IO`$. We should note that a similar deviation is present in the WWW. From the data of Ref.  for the WWW, $`I0.490`$, $`O0.489`$, so $`IO0.240`$, that is less than the measured value $`S0.277`$ but is not far from it.

Figures 24 demonstrate that, in the wide enough range of parameters, the following situation may be realized. The directed graph may have the GWCC and have no the GSCC. Only the stability of the GWCC to random damage was discussed yet . Nevertheless, just the stability of the GSCC is the most important, e.g., for the WWW. Is it possible that the GSCC are less resilient to failures than the GWCC? Let us briefly discuss this problem.

One can see from Eq. (5), that the GSCC is extremely resilient if the average $`k_ik_o`$ diverges. Let us consider two limiting situations. In the first one, the joint in,out-degree distribution factorizes, so $`k_ik_o=k_i^2`$. In this case, if the distributions are of a power-law form then, for the robustness of the GSCC, the corresponding exponent $`\gamma _i`$ or $`\gamma _o`$ should be smaller than $`2`$ or equal it. This is a very strong requirement. Here, the smallest exponent of $`\gamma _i`$ and $`\gamma _o`$ is also equal to the exponent $`\gamma `$ of the degree distribution. For the resilience of the GWCC to random damage, it should be $`\gamma 3`$  that ensures the divergence of $`k^2`$. Therefore, for such distributions, when $`2<\mathrm{min}(\gamma _i,\gamma _o)=\gamma 3`$, random damage can destroy the GSCC but can not eliminate the GWCC.

In the other limiting case, $`P(k_i,k_o)=P(k_i)\delta _{k_i,k_o}`$, the correlations between in- and out-degrees are very strong. This form more resembles the joint degree distributions of the real growing networks. In such an event, $`k_ik_o=k_i^2`$, and the conditions for the resilience of the GSCC and GWCC, $`\gamma 3`$, coincide. One should note that the real distributions are between the considered limiting cases.

In summary, we have shown how to obtain the size of the giant strongly connected component of the directed network with the arbitrary degree distribution and statistically uncorelated vertices. This allows us to find all the giant components of such a graph and to describe its basic structure. Using the simplest examples and the general considerations we have demonstrated that the correlations between in- and out-degrees subsequently influence the global topology of the network.

SND thanks PRAXIS XXI (Portugal) for a research grant PRAXIS XXI/BCC/16418/98. SND and JFFM were partially supported by the project POCTI/1999/FIS/33141. We also thank P.L. Krapivsky for useful discussions.

Electronic address: sdorogov@fc.up.pt
Electronic address: jfmendes@fc.up.pt
Electronic address: alnis@samaln.ioffe.rssi.ru