Kompaneets equation for neutrinos: Application to neutrino heating in supernova explosions

We derive a `Kompaneets equation' for neutrinos, which describes how the distribution function of neutrinos interacting with matter deviates from a Fermi-Dirac distribution with zero chemical potential. To this end, we expand the collision integral in the Boltzmann equation of neutrinos up to second order in energy transfer between matter and neutrinos. The distortion of the neutrino distribution function changes the rate at which neutrinos heat matter, as the rate is proportional to the mean square energy of neutrinos, $E_\nu^2$. For electron-type neutrinos the enhancement in $E_\nu^2$ over its thermal value is given approximately by $E_\nu^2/E_{\nu,\rm thermal}^2=1+0.086(V/0.1)^2$ where $V$ is the bulk velocity of nucleons, while for the other neutrino species the enhancement is $(1+\delta_v)^3$, where $\delta_v=mV^2/3k_BT$ is the kinetic energy of nucleons divided by the thermal energy. This enhancement has a significant implication for supernova explosions, as it would aid neutrino-driven explosions.


Introduction
The Boltzmann equation is ubiquitous in physics. In a system in which matter and radiation interact, a useful approximation can be obtained by expanding the collision integral in the Boltzmann equation up to second order in energy transfer between matter and radiation.
One example is Kompaneets's equation [1]. This equation describes how the distribution function of photons evolves via interaction with thermal electrons in the non-relativistic limit. A solution to the Kompaneets equation in the optically-thin limit is known as the thermal Sunyaev-Zel'dovich (tSZ) effect [2], which describes a distortion of the black-body spectrum of the cosmic microwave background photons by inverse-Compton scattering off hot electrons in galaxy clusters. The tSZ effect has been routinely detected towards O(10 3 ) galaxy clusters [3]. Another example is the Fermi acceleration [4]. The so-called second-order Fermi acceleration distorts the distribution of charged particles by stochastic acceleration due to time-dependent electromagnetic fields. Both examples can be formulated in the same form; namely, a diffusion equation for the distribution function of photons or charged particles in momentum space.
In this paper, we apply the same approximation to the Boltzmann equation describing neutrinos interacting with matter. Specifically, we consider a system in which isotropic neutrinos interact with nucleons, and expand the collision integral in the Boltzmann equation up to second order in energy transfer between neutrinos and nucleons. We do not assume that the system is optically thin. By solving this Kompaneets-like equation for neutrinos, we obtain distortions of the neutrino distribution function from a Fermi-Dirac distribution with zero chemical potential as a function of the temperature and bulk velocity of nucleons. Our result has a significant implication for neutrino-driven supernova explosions, as the distortion of the distribution function changes the rate at which neutrinos heat nucleons.
The rest of this paper is organized as follows. In Section 2, we derive the Kompaneets-like equation for neutrinos interacting with nucleons. Nucleon's motion includes both thermal and bulk motion. We present solutions of this equation in Section 3, including the effect of opacity of electron-type neutrinos. We summarize our results and discuss their implications in Section 4. In Appendix A we review the matrix element of neutrino-nucleon scattering. In Appendix B we provide an alternative derivation of the main result of this paper following the argument of [5]. Throughout this paper we shall set the speed of light to be unity, c = 1.

Kompaneets equation for neutrinos
We follow derivations of the Kompaneets equation for photons interacting with electrons given in ref. [6,7], and derive a similar equation for neutrinos interacting with nucleons. To this end, we expand the collision integral of the Boltzmann equation up to second order in energy transfer. We shall ignore the mass of neutrinos throughout this paper, as the typical neutrino temperature that we consider in this paper (e.g., that in the supernova engine) is much greater than the current upper bound on the mass of neutrinos on the order of 1 eV.
The Boltzmann equation is where N denotes nucleons (neutrons and protons), d/dt is the Lagrangian time derivative along with the trajectory of a phase-space volume element, |M N | 2 is the spin-averaged matrix element of neutrino-nucleon scattering, f (p, t) and g N (q) are the distribution functions of neutrinos and nucleons with 3-momenta p and q, respectively. The 4-dimensional Dirac delta function δ (4) (p + q − p ′ − q ′ ) ensures energy and momentum conservation.
We assume that nucleons are non-relativistic, i.e., E N (q) ≡ m 2 N + |q| 2 ≈ m N + |q| 2 /(2m N ), with the distribution function given by 2 /14 where n N , m N , T N and v N are the number density, mass, temperature and bulk velocity of nucleons, respectively, and k B is the Boltzmann constant. Performing integration over q ′ , we obtain The energy transfer is on the order of The next step is to expand |M N | 2 , E(q + p − p ′ ), the delta function and the distribution functions up to second order in energy transfer. In this paper we consider elastic neutrinonucleon scattering, i.e., ν + N → ν + N , and ignore other processes.
There are three types of the so-called "neutrinosphere", which are determined by different micro processes [8]. Namely, • Number sphere: The optical depth of emission and absorption of neutrinos is of order unity. For electron-type neutrinos ν e andν e , the electron/positron capture and its inverse processes (i.e., ν e + n ↔ p + e − andν e + p ↔ n + e + ) are important. For other two types of neutrinos ν µ , ν τ ,ν µ ,ν τ , the pair production/annihilation processes (i.e., νν ↔ γγ, νν ↔ e + e − , N N ↔ N N νν) determine the opacity. • Energy sphere: Inelastic scattering by electrons is important here. Electrons receive energy from neutrinos, as the electron rest mass energy (511 keV) is much smaller than the typical neutrino energy (∼ 10 MeV), which is determined by the matter temperature at the number sphere. Inside the energy sphere the neutrinos are thermalized due to energy transfer with electrons, which are tightly coupled with baryons. • Transport sphere: Beyond the energy sphere, elastic scattering by nucleons and nuclei determines the opacity. As the rest mass energy of these particles is much larger than the typical neutrino energy, scattering can be treated as elastic 1 .
For ν e andν e , all neutrinospheres are almost coincident, and thus the spectrum is almost thermal. For ν µ , ν τ ,ν µ ,ν τ , on the other hand, the radii of neutrinospheres are separated from each other [8], and thus neutrinos may establish a nonthermal component, which would be produced between the energy and transport spheres. This is the region we consider in this paper.
1 Note that ref. [8] investigated how the recoil term affects the neutrino spectrum.

3/14
With the matrix element for neutrino-nucleon scattering given in Appendix A, the right hand side of Eq. (3) becomes where ∓ takes − for neutrinos and + for anti-neutrinos, and f (0) (p) is the zeroth-order distribution function of neutrinos, which is a Fermi-Dirac distribution with zero chemical potential. We shall ignore the feedback of distorted neutrino spectrum in this paper (f (1) and f (2) in [7]). Here,p is a unit vector, cos θ =p ·p ′ , G F is the Fermi coupling constant and F N 1 , F N 2 and F N A are the form factors appearing in |M N | 2 (see Appendix A). The derivative of the Dirac delta function will be handled by integration by parts.
Eq. (8) is the main result of this paper. Comparing this to the result for photon-electron scattering (see Eq. (15) of [6]), we find two differences. First, the power of p is different by two because the neutrinonucleon cross section is proportional to p 2 whereas the Thomson scattering cross section σ T is independent of photon monenta. Second, we have 1 − f (0) in the last term in Eq. (8) instead of 1 + f (0) in Eq. (11), because of Fermi statistics.

Chemical potential distortion
From now on we shall assume that protons and neutrons share the same temperature T and bulk velocity V . We shall approximate that masses of protons and neutrons are the same and denote them as m, i.e., m N → m, and drop the superscript (0) on f . Let us define dimensionless variables x ≡ p/(k B T ) and where We then write Eq. (8) as where the prime denotes derivative with respect to x, and δ v ≡ mV 2 /(3k B T ) is the ratio of nucleon's bulk kinetic energy to thermal energy.

5/14
An equilibrium solution, df /dy = 0, is given by where µ ν is a chemical potential. The bulk velocity effectively enlarges the temperature by a factor of 1 + δ v . As the number of neutrinos is conserved in our set up, we can obtain the emergent spectrum after stochastic scattering by bulk fluid motion as follows. For initially a thermal spectrum (i.e., µ ν = 0), we get the same number density of neutrinos if we have µ ν = −0.341, −4.34, and −80.3k B T for δ v = 0.1, 1, and 10, respectively. Therefore, scattering distorts the neutrino thermal spectrum by giving a non-zero chemical potential, which is well known as "µ-distortion" in the cosmic microwave background research [10].
With the spectrum of Eq. (15), we can calculate the neutrino-annihilation rate which is one of the important heating processes in supernova explosion. The energy deposition rate via neutrino annihilation (ν +ν → e + + e − ) is given by [11,12] where The subscriptν indicates the corresponding quantities for anti-neutrinos. The overall factor C includes the weak interaction coefficients and information of the angular distribution of neutrinos. To calculate this factor we need to determine geometry of the neutrino-emitting source. Since this factor is expected not to change significantly by including the neutrino acceleration process, we concentrate on the effect of the spectral change here. For simplicity, we assume that the spectra of ν andν are identical. Then, we obtainĖ Here we use F 3,ν ∝ ǫ ν , as F 2,ν does not change by scattering processes alone. Roughly speaking, ǫ ν = (1 + δ v ) ǫ ν thermal and ǫ 2 ν = (1 + δ v ) 2 ǫ 2 ν thermal , where ǫ ν thermal and ǫ 2 ν thermal are mean energy and mean square energy based on purely-thermal distribution function. Therefore, we findĖ

Including opacity
The calculation of the energy deposition rate given above includes only scattering; thus, it can describe only µ-and τ -type neutrinos. However, there are also terms related to emission and absorption. Absorption is important for electron-type neutrinos, which leads to thermalization of the spectrum of ν e andν e . In the following we include all relevant terms. The revised transfer equation in which both O(v) and O(v 2 ) terms are taken into account is written as follows [see also 13, 14, for photon cases], where κ a and κ sc are opacities (inverse of the mean free path) for absorption and scattering, respectively, κ t ≡ κ a + κ sc , and f eq is the neutrino distribution function in thermal equilibrium. When we omit the second line in the right hand side of this equation, we get exactly the same equation for neutrino transfer based on diffusion approximation, which is solved in a number of numerical simulations; see Eq. (A27) of [15] for a flux-limited diffusion approximation and Eqs. (5) and (6) of [16] for an isotropic diffusion source approximation. Both of them are originally derived from the Boltzmann equation including velocity dependent terms up to O(v) and are commonly used in neutrino-radiation hydrodynamics simulations. The effect of the second term proportional to ∇ · V in the right hand side is analogous to the bulk Comptonization process in photon cases and the first-order Fermi acceleration for charged particles. This term modifies the neutrino spectrum as well when there is a compressional flow, such as a shock or an accretion flow onto a black hole. In core-collapse supernovae, a neutron star is formed and compression is almost negligible in the optically thick region for neutrinos so that this first-order acceleration does not work at all [17].

Numerical solution
Assuming one-zone (i.e. ∇f = 0) and incompressible bulk flow (i.e. ∇ · V = 0), Eq. (19) reads The scattering opacity is given from Eqs. (9), (10) and (13) as where Y n = n n /n b and Y p = n p /n b are the number fractions of free neutrons and protons with n b being the number density of baryons, and m = 1.67 × 10 −24 g. Since all nuclei are photodisintegrated into neutrons and protons, Y n + Y p = 1 here. The absorption opacity (ν e + n → p + e − for ν e andν e + p → n + e + forν e ) is given by [9] κ a = 3g 2 where g A = 1.2723 is the axial-vector coupling constant, and σ 0 = 1.705 × 10 −44 cm 2 is a reference neutrino cross section. The small corrections due to the mass difference between a neutron and a proton, and to weak magnetism and recoil are neglected here. For characteristic values Y ≡ Y n = Y p = 0.5 is used in Eq. (22). For constant temperature and chemical potential, V 2 = 0, and f eq ν = (e (ǫν −µν )/kBT + 1) −1 , we have a steady-state solution (∂f /∂t = 0) of Eq.  The initial spectrum and an analytic solution without absorption are shown by the blueand green-dashed lines, respectively. For numerical solutions, two different values of y are used (y = 1 and 10 −5 for the thick and thin lines, respectively). (Bottom) Difference between numerical solutions at y = 1 and the initial spectrum. The low energy part (x 7) is reduced, while the high energy part is increased due to up-scattering.
Writing the Boltzmann equation in a dimensionless form, we have where Θ T ≡ (κ sc /κ a )(k B T /m). In Figure 1, numerical solutions of Eq. (23) are shown by the red solid lines. For this calculation, we use k B T = 6 MeV, δ v = 1, and Y p = 0.1. Since we use κ a ∝ Y p , solutions shown in this figure are valid forν e . The initial condition at y = 0 is given by (e x + 1) −1 , which is shown by the blue-dashed line in this figure. The green-dashed line is an analytic solution of Eq. (15), which is realized when the absorption is not taken into account. We find that, for higher energy (x 10), the numerical solutions deviate from the thermal equilibrium (blue-dashed line), and the spectra become slightly non-thermal toward the analytic solution (green-dashed line). This is because the ratio of the two terms in the right hand side of Eq. (20) is ∝ p/m, such that the absorption term dominates in the low energy regime. We find that the numerical solutions are very similar for y 0.001. 8/14

Implications for neutrino heating rate
Since the neutrino interaction rate is a strong function of neutrino energy, the neutrino heating rate is affected when the neutrino spectrum changes. The neutrino heating rate is given as [18] where With the numerical solutions for the neutrino spectrum (Figure 1), we find the following results: They are valid for δ v 1. For larger δ v , non-linear terms in δ v contribute. Combining them gives an enhancement in the square mean energy over its thermal value as Thus, the neutrino heating rate Q + ν is amplified inside the gain region, where the neutrino heating overwhelms neutrino cooling between shock and neutrinospheres (see Eqs. (24) and (27)). The neutrino cooling rate is reduced inside the protoneutron star (see Eqs. (26) and (28)). Both of them make shock revival of supernova easier.

Summary and discussion
In this paper, we have derived the Kompaneets-like equation for neutrinos by expanding the collision integral of the Boltzmann equation up to second order in energy transfer, or O(v 2 ), including thermal and bulk velocities of nucleons. We also include absorption and emission of neutrinos, to arrive at the full neutrino transport equation given in Eq. (19). The dimensionless form of the equation suitable for numerical calculations is given in Eq. (23), and the numerical solutions are presented in Figure 1.
We find that the distortion of the neutrino spectrum due to interaction with nucleons leads to a larger neutrino heating rate in the gain region and a smaller neutrino cooling rate in the protoneutron star, which provides a better condition of supernova explosions than solutions without the effects we find in this paper.
The formulation given in this paper is a natural extension of the transfer equation that is solved in neutrino-radiation hydrodynamics simulations for core-collapse supernovae, in which velocity dependent terms in the collision integral are included only up to O(v). Note that the current equation is derived in the non-relativistic limit, i.e., k B T ≪ m and V ≪ 1. Relativistic corrections are known to amplify spectral distortion of the SZ effect [e.g. 19], which could even enhance the neutrino heating rates. 9/14 The numerical solutions presented in Section 3.3 give the emergent spectrum of electrontype neutrino (ν e ) and anti-neutrinos (ν e ), as we include the absorption and emission term (the last term of the right hand side in Eq. 23), which is relevant only for charged-current reactions. For other heavier-leptonic flavors, i.e., muon-type and tauon-type, this term does not appear in the transfer equation so that their spectrum would be like Eq. (15). The neutrino annihilation rate of these heavier-leptonic neutrinos [11,20] can also be amplified by a factor (1 + δ v ) 3 , which would make supernova explosion easier as well.
Our finding motivates self-consistent hydrodynamics simulations including the spectral distortion of neutrinos, which will be needed to calculate the quantitative impact on supernova explosions.

A. Neutrino-nucleon scattering A.1. Matrix element
In this section, we follow [21,22] to write the expression for the spin-averaged matrix element of neutrino-nucleon scattering. It is given by where ∓ takes − for neutrinos and + for anti-neutrinos, and where all momenta are 4-vectors, and F N 1 , F N 2 and F N A are the so-called neutral-current Dirac, Pauli and axial form factors, respectively, which depend on Q 2 ≡ −(p − p ′ ) 2 .
When nucleons are non-relativistic, we have where the momenta on the left hand side are 4-vectors, whereas those on the middle and right hand side are 3-vectors and their magnitudes. p c.m. and θ c.m. are the energy and scattering 10/14 angle in the center-of-mass frame. Eliminating p c.m. and θ c.m. , we obtain . (A8) Using this relationship to expand the matrix element in powers of p/m N and q/m N up to their linear order, we obtain wherep ·p ′ ≡ cos θ.

A.2. Total cross section
In the center-of-mass frame, the differential cross section is given by Integrating this over angles, we obtain the total cross section in the limit of non-relativistic nucleons (s ≈ m 2 N ) as For the energy scale we are interested in, the only relevant form factors are the ones with Q 2 = 0. Thus, Here, g A = 1.2723 is the axial-vector coupling constant, sin 2 θ W = 0.23122 is the weak angle, G F = 1.1663787 × 10 −5 GeV −2 is the Fermi coupling constant, and µ p and µ n are magnetic moments of protons and neutrons, respectively. As F N 2 does not contribute to the angleaveraged collision integral (8), we do not need to evaluate them in this paper. We thus find (A15)

B. Alternative derivationà la Rybicki and Lightman
In this section, we provide an alternative derivation of the main result of this paper, Eq. (8), following the argument given in Section 7.6 of [5]. For ease of comparison, let us change our notation of momenta to follow closely that of [5]: q → p, q ′ → p 1 , p → ωn, and p ′ → ω 1n1 .

11/14
The Boltzmann equation is where dσ N /dΩ = (σ N /4π)(1 + δ N cos θ) (with N = p and n) is the differential cross section of scattering by free nucleons, g N (p) the distribution function of nucleons with 3-momenta p, ω and ω 1 the neutrino energies,n andn 1 the unit vectors of the neutrino 3-momenta, and cos θ =n ·n 1 .

B.1. Thermal nucleons
Since the typical energy transfer is small compared to nucleon's kinetic energy, we define and expand the right hand side of Eq. (B1) up to second order in ∆. In a frame in which the initial nucleon velocity is β = p/E with E = |p| 2 /(2m N ), we have Note that β ≡ |β| ≈ k B T /m N ≪ 1. Comparing Eqs. (B3) and (B2), we find where x ≡ ω/(k B T ). Expanding the neutrino distribution function with ω 1 up to second order in ∆, we obtain where the primes denote derivatives with respect to ω. We also expand the distribution function of non-relativistic nucleons up to second order in ∆ = ( Here, n N is the number density of nucleons. Then Eq. (B1) becomes To check the calculation, note that a Fermi-Dirac distribution f FD (x) = 1/(e x+η + 1) is a steady-state solution to Eq. (B6). Here, η is an integration constant corresponding to a chemical potential term. 12/14 Let us evaluate the integral in the second term in Eq. (B6). Writing (n 1 −n) · p = |n 1 − n|p cos ξ andn ·n 1 = cos θ, we have Using d 3 p = p 2 dpd(cos ξ)dϕ and dΩ = d(cos θ)dφ, the integral evaluates to The integral in the first term in Eq. (B6) is more complicated. Therefore, we follow [5] and use the conservation law to calculate it. Since scattering does not change the number of neutrinos, dxx 2 ∂f /∂t = 0. This equation is satisfied by the continuity equation, where j is a flux of neutrinos in momentum space. Next we write the unknown integral as The Boltzmann equation then becomes where α N ≡ 1 − δ N /3. Requiring this equation be equal to Eq. (B9), we write an ansatz: with two unknown functions g and h, so that j ′ has a term in f ′′ but no higher derivatives. We then determine the function h from the equilibrium form, f FD , which has the following relation; f ′ FD + f FD (1 − f FD ) = 0. If j vanishes in equilibrium, then h = f (1 + f ). We find Comparing the f ′′ terms in Eqs. (B11) and (B13), we obtain Comparing the f (1 − f ) terms gives which yields I = α N (6 − x)x. Note that this makes f ′ terms consistent as well. The final result is the following equation: By changing the variable from t to y via dy ≡ N α N kBT mN n N σ N (x = 1)dt, we obtain df dy which agrees with Eq. (14) without the bulk velocity effect δ v . 13/14

B.2. Including bulk motion
In the above derivation, we have used the Maxwell distribution for nucleons. Next, let us include the bulk motion. The distribution function of nucleons including bulk motion is given by g N (p) = n N (2πm N k B T ) −3/2 e −(p−mN vN ) 2 /2mN kBT . Then, the bulk velocity of a fluid, v N , is given by where . . . denotes the average over the momentum distribution of nucleons. The second moment of p/m N is where v ≡ p/m N − V is the velocity of a nucleon in the rest frame of a fluid element.
Assuming that nucleon's momentum distribution in the rest frame of a fluid element is a Maxwell distribution, we have v 2 = 3k B T /m N . Including the bulk velocity in Eq. (14), we obtain df dy which agrees with Eq. (14).