Bulk viscosity of baryonic matter with trapped neutrinos

TL;DR

This paper investigates the bulk viscosity of dense baryonic matter with trapped neutrinos, revealing how it varies with temperature and composition, and assessing its impact on neutron star mergers and gravitational wave damping.

Contribution

It models the bulk viscosity in neutrino-trapped dense matter using a relativistic density functional approach, highlighting its temperature dependence and the effects of lepton fraction.

Findings

Bulk viscosity peaks near the neutrino trapping temperature.

Bulk viscosity decreases as T^{-2} in the neutrino-trapped regime.

Bulk viscosity is significantly lower than in the neutrino-transparent regime.

Abstract

We study bulk viscosity arising from weak current Urca processes in dense baryonic matter at and beyond nuclear saturation density. We consider the temperature regime where neutrinos are trapped and therefore have nonzero chemical potential. We model the nuclear matter in a relativistic density functional approach, taking into account the trapped neutrino component. We find that the resonant maximum of the bulk viscosity would occur at or below the neutrino trapping temperature, so in the neutrino trapped regime the bulk viscosity decreases with temperature as $T^{-2}$, this decrease being interrupted by a drop to zero at a special temperature where the proton fraction becomes density-independent and the material scale-invariant. The bulk viscosity is larger for matter with lower lepton fraction, i.e., larger isospin asymmetry. We find that bulk viscosity in the neutrino-trapped regime is smaller by several orders of magnitude than in the neutrino-transparent regime, which implies that bulk viscosity in neutrino-trapped matter is probably not strong enough to affect significantly the evolution of neutron star mergers. This also implies weak damping of gravitational waves emitted by the oscillations of the postmerger remnant in the high-temperature, neutrino-trapped phase of evolution.

Tables

Table 1. Table 1: Meson masses m i subscript 𝑚 𝑖 m_{i} and their couplings g i ​ N subscript 𝑔 𝑖 𝑁 g_{iN} to the baryons in DD-ME2 parametrization. The remaining parameters specify the density dependence of the couplings.

	$\sigma$	$\omega$	$\rho$
$m_i$ [MeV]	550.1238	783.0000	763.0000
$g_{iN}\left(n_0\right)$	10.5396	13.0189	3.6836
$a_i$	1.3881	1.3892	0.5647
$b_i$	1.0943	0.9240	—
$c_i$	1.7057	1.4620	—
$d_i$	0.4421	0.4775	—

Equations

$$n\rightleftarrows p+e^{-}+{\bar{\nu }}_e,$$

$$p+e^{-}\rightleftarrows n+{\nu }_e.$$

$${\mu }_n={\mu }_p+{\mu }_e+{\mu }_{\bar{\nu }},$$

$${\mu }_p+{\mu }_e={\mu }_n+{\mu }_{\nu },$$

$${\mu }_{\Delta }\equiv {\mu }_n+{\mu }_{\nu }-{\mu }_p-{\mu }_e.$$

$$\sum |{\mathcal{M}}_{Urca}{|}^2=32{\tilde{G}}^2(k\cdot p^{\prime })(p\cdot k^{\prime }),$$

$$|{\mathcal{M}}_{Urca}{|}^2=32{\tilde{G}}^2{p}_0{p}_0^{\prime }{k}_0{k}_0^{\prime },$$

$${\Gamma }_{1p}({\mu }_{\Delta })$$

$${\Gamma }_{1n}({\mu }_{\Delta })$$

$${\Gamma }_{1p}({\mu }_{\Delta })$$

$${\Gamma }_{1n}({\mu }_{\Delta })$$

$$\displaystyle y_{0}=\frac{m^{*}}{2Tz^{2}}\bigg{(}\alpha_{n}-\alpha_{p}+y-z^{2}\frac{T}{2m^{*}}-\frac{\mu_{\Delta}}{T}\bigg{)}^{2}-\alpha_{p},$$

$${\lambda }_1$$

$${\Gamma }_1$$

$${\lambda }_1$$

$${\Gamma }_{1p}=\alpha m^{\ast 2}{\tilde{G}}^2{T}^5{p}_{Fe}\theta (p_{Fp}+p_{Fe}-p_{Fn}).$$

$${\lambda }_1=\frac{17}{480\pi }{m}^{\ast 2}{\tilde{G}}^2{T}^4{p}_{Fe}\theta (p_{Fp}+p_{Fe}-p_{Fn}).$$

$${\Gamma }_{2p}({\mu }_{\Delta })$$

$${\Gamma }_{2n}({\mu }_{\Delta })$$

$${\Gamma }_{2p}({\mu }_{\Delta })$$

$${\Gamma }_{2n}({\mu }_{\Delta })$$

$${\lambda }_2$$

$${\Gamma }_2$$

$${\Gamma }_2=\frac{m^{\ast 2}{\tilde{G}}^2}{12{\pi }^3}{T}^3{p}_{Fe}{p}_{F\nu }(p_{Fe}+p_{F\nu }-|p_{Fn}-p_{Fp}|).$$

$${\lambda }_2$$

$${\Gamma }_{2n}=\alpha m^{\ast 2}{\tilde{G}}^2{T}^5{p}_{Fe}\theta (p_{Fp}+p_{Fe}-p_{Fn}).$$

$${\lambda }_2=\frac{17}{480\pi }{m}^{\ast 2}\widetilde{G^2}{T}^4{p}_{Fe}\theta (p_{Fp}+p_{Fe}-p_{Fn}),$$

$$\lambda =\frac{17}{240\pi }{m}^{\ast 2}\widetilde{G^2}{T}^4{p}_{Fe}\theta (p_{Fp}+p_{Fe}-p_{Fn}),$$

$$\delta n_i(t)=-\frac{\theta }{i\omega }{n}_{i0},i=\{B,L\},$$

$$n_j(t)=n_{j0}+\delta n_j(t),$$

$${\mu }_{\Delta }(t)=(A_{nn}-A_{pn})\delta n_n(t)+A_{\nu \nu }\delta n_{\nu }(t)-(A_{pp}-A_{np})\delta n_p(t)-A_{ee}\delta n_e(t)\equiv \sum _i{s}_i{A}_i\delta n_i(t)$$

$$A_{ij}={\left(\frac{\partial {\mu }_i}{\partial n_j}\right)}_0,$$

$$\frac{\partial }{\partial t}\delta n_j(t)+\theta n_{j0}=0,\delta n_j(t)=-\frac{\theta }{i\omega }{n}_{j0}.$$

Full text

I Introduction

The recent detection of gravitational waves by the LIGO-Virgo collaboration, in coincidence with electromagnetic counterparts, has brought into focus the study of the physics of binary neutron star mergers Abbott2017 . In these events, a postmerger object is formed which either evolves into a stable neutron star or collapses to a black hole, once it cannot be supported by the differential rotation. As seen in numerical simulations Perego:2019adq ; Hanauske:2019qgs ; Hanauske:2017oxo ; Kastaun:2016elu ; Bernuzzi:2015opx ; Foucart:2015gaa ; kiuchi:2012mk ; sekiguchi:2011zd ; Ruiz:2016 ; East:2016 ; Most2019 ; Bauswein2019 ; Baiotti:2016qnr ; Faber2012:lrr there are significant density oscillations in the postmerger remnant, which can generate observable gravitational waves.

These oscillations will be damped eventually by dissipative processes on characteristic secular timescales controlled by thermodynamics of background matter and the kinetics of the relevant dissipative process. Bulk viscosity is known to be one of the dissipative processes that could efficiently damp certain classes of oscillations of general relativistic equilibria.

Studies of damping mechanisms in the context of binary neutron star mergers are still at an embryonic stage. Reference Alford2017 has suggested that modified Urca processes can produce significant bulk viscous dissipation on timescales of order a few milliseconds, i.e., on timescales that are relevant to neutron star merger and postmerger evolution. This might then affect the emitted gravitational signal. It remains less clear whether shear viscosity and thermal conduction could play a significant role Alford2017 . In parallel, the electrical conductivity was computed in the relevant regime to assess its impact on the evolution of the electromagnetic fields. It was shown that the Hall effect could be important on characteristic timescales of the merger and postmerger evolution Harutyunyan:2016a ; Harutyunyan:2018 .

To quantify the amount of bulk viscous dissipation a more detailed analysis is required that takes into account the realistic temperature and density conditions encountered in this context. The aim of this work is to obtain the bulk viscosity of dense matter created in neutron star mergers in the temperature and density regime characteristic for such events. Specifically, we will consider the dominant weak-interaction processes of Urca type at temperatures that are above the trapping temperature $T_{\rm tr}\simeq 5$ MeV Roberts:2012um ; Alford:2018lhf . In this regime, the neutrinos have a mean free path that is significantly shorter than the stellar size, and consequently a nonzero chemical potential. This affects the composition of the background baryonic matter. Thus, compared to the extensively studied case of cold neutron stars, the key new features that arise in the neutron-star merger and postmerger context is the higher temperatures at which the weak reactions take place and the significantly different background composition of baryonic matter.

The bulk viscosity of baryonic matter has been studied extensively in the low-temperature limit for purely nucleonic matter Sawyer1979ApJ ; Sawyer1980ApJ ; Sawyer1989 ; Haensel1992PhRvD ; Dong2007 ; Alford2010JPhG ; Kolomeitsev2015 , including possible leptonic contributions Alford:2010jf , for hyperonic matter Jones2001PhRvD ; Lindblom2002 ; Dalen2004PhRvC ; Haensel2002 ; Chatterjee2006 ; Chatterjee2008ApJ , and for quark matter Madsen1992 ; Drago2005 ; Alford:2006gy ; Alford:2007rw ; Manuel:2007pz ; SaD2007a ; SaD2007b ; Alford:2008pb ; Sinha2009 ; Wang2010PhRvD ; Huang2010 including the effects of an interface with the nuclear matter envelope AlfordHan2015 . For a review see Schmitt2017 . At temperatures $T\lesssim 1$ MeV the physics of bulk viscosity is affected by the superfluidity of baryons Haensel2000 ; Haensel2001 ; Gusakov2007 , but we will assume that the temperatures are always larger than the critical temperatures for the pairing of various baryons, as our focus is on the regime where neutrinos are trapped, i.e., temperatures $T\geq T_{\rm tr}$. The neutrino transparent case is discussed in Ref. Alford:2019zzz and Appendix D.

Physical conditions that are similar to those we are aiming to study are encountered also in proto-neutron stars born in supernova explosions. In this case, the matter is much more isospin symmetrical, but the rest of the physics is quite analogous. We will cover this case as well having in mind the possibility of observations of oscillations of a proto-neutron star, should a supernova explosion occur within the detectable range.

In the present study we will assume that thermal conduction is efficient enough to keep matter isothermal while it is being compressed and uncompressed. A rough estimate (Ref. Alford2017 , Eq. (2)) gives the timescale for relaxation of thermal gradients order of 1 sec $\times(z_{\rm typ}/{\rm km})^{2}(T/10\,{\rm MeV})^{2}$, where $z_{\rm typ}$ is the typical scale of thermal gradients. Thus, for thermal gradients on a distance scale of about 30m or less, the relaxation time would be 1 ms or less, i.e., the characteristic time-scale of binary neutron star merger. So, for density fluctuations on this distance scale the assumption of isothermal matter is the relevant one. In particular, if turbulent flow arises in the merger then this could give flows and density variations on such distance scale. On the scales over which thermal conduction is inefficient, the matter should be treated as isoentropic. The formalism presented below can be simply adapted to this case. We anticipate that for an adiabatic calculation the bulk viscosity will be of the order as found here. We will return to this problem in a separate work.

This work is organized as follows. In Sec. II we discuss the rates of the two relevant processes. Section III derives the corresponding formulas for the bulk viscosity of matter. In Sec. IV we first describe the properties of the background matter derived on the basis of the density functional theory at a finite temperature which accounts for a neutrino component with nonzero chemical potential (Sec. IV.1). This is followed by a discussion of perturbed quantities and bulk viscosity in Sec. IV.2. Our conclusions are collected in Sec. V.

In this work we use the natural (Gaussian) units with $\hbar=c=k_{B}=1$, and the metric $g_{\mu\nu}=\textrm{diag}(1,-1,-1,-1)$.

II Urca process rates

We will consider the simplest composition of baryonic matter consisting of neutrons ($n$), protons ($p$), electrons ($e$), muons ($\mu$), and neutrinos at densities in the range $0.5n_{0}$ to $3n_{0}$ ($n_{0}\simeq 0.16$ fm*-3*) and temperatures in the range $T_{\rm tr}\simeq 5$ to 50 MeV. Other constituents and forms of matter have been proposed, but we will focus on the standard scenario for this regime, which can serve as a starting point for future explorations of more complex phases of baryonic matter. Note that positrons do not appear in matter in substantial amounts because the electron chemical potential is of the order of 100 MeV, see Sec. IV.1. The weak processes involving positrons will be suppressed roughly by a factor $\exp(-\mu_{e}/T)\simeq 0.1$ at $T=50$ and 0.01 at $T=30$ MeV.

In the dynamically evolving environment of a neutron star merger, fluid elements undergo rhythmic cycles of compression and decompression, which can lead to bulk viscous dissipation if the rate at which the proton fraction relaxes toward its equilibrium value (“beta equilibrium”) is comparable to the frequency of the compression cycles. To analyze this, we consider the simplest beta equilibration processes, the Urca processes:

[TABLE]

The first process is the $\beta$-decay of a neutron and the second one is the electron capture on a proton. If the matter is in $\beta$-equilibrium, then the chemical potentials of particles obey the relation

[TABLE]

where the neutrino and antineutrino chemical potentials are related by $\mu_{\bar{\nu}}=-\mu_{\nu}$, which leaves us with a single relation. As noted above, the matter can be driven out of $\beta$-equilibrium by a cycle of compression and rarefaction, and this can be characterized via a nonzero value of the chemical potential that measures the deviation from $\beta$-equilibrium

[TABLE]

The rate at which $\mu_{\Delta}$ relaxes to zero is a measure of the speed at which the chemical constitution of the matter adjusts to a change in pressure. We start with the computation of the $\beta$-equilibration rate assuming a given value of $\mu_{\Delta}\neq 0$.

The squared matrix element of processes (1) and (2) is given by the well-known expression Greiner2000gauge

[TABLE]

where $p^{\prime}$ and $p$ refer to the four-momenta of the neutron and proton, $k$ and $k^{\prime}$ to the four-momenta of neutrino (antineutrino) and electron, respectively, and the sum includes summation over the spins of neutron, proton, and electron. Note that each of baryon four-momenta is dotted into a four-lepton momentum. We consider only the Standard Model neutrinos (antineutrinos) which are left-handed (right-handed) only, therefore they have only one projection of helicity that has to be counted. In the following, we will use the nonrelativistic limit of the matrix element (6) because in the temperature and density range that we consider the baryons are nonrelativistic to an accuracy of about 10%.

Thus we keep only the contribution of timelike parts of the scalar products in the matrix element,

[TABLE]

where index 0 refers to the timelike component of a four-vector, $\tilde{G}^{2}\equiv G_{F}^{2}\cos^{2}\theta_{c}(1+3g_{A}^{2})$, where $G_{F}=1.166\cdot 10^{-5}$ GeV*-2* is the Fermi coupling constant, $\theta_{c}$ is the Cabibbo angle ($\cos\theta_{c}=0.974$) and $g_{A}=1.26$ is the axial-vector coupling constant.

II.1 The rates of the processes $n\rightleftarrows p+e^{-}+\bar{\nu}_{e}$

The $\beta$-equilibration rate for the neutron decay $n\rightarrow p+e^{-}+\bar{\nu}_{e}$ is given by

[TABLE]

where $f(p)=\{\exp[(E_{p}-\mu)]+1\}^{-1}$ etc. are the Fermi distributions of particles, with energies $E_{p}$ and chemical potential $\mu$. Similarly, the rate of the inverse process, i.e., $p+e^{-}+\bar{\nu}_{e}\rightarrow n$ is given by

[TABLE]

Some of the phase-space integrals in Eqs. (8) and (9) can be carried out analytically; the details are relegated to Appendix A. We find

[TABLE]

where $\alpha_{j}\equiv\mu_{j}^{*}/T$ and index $j=n,p,e,\nu$ labels the constituents of matter with $\mu_{j}^{*}$ being the effective chemical potentials of particles (see Sec. IV.1), $m^{\ast}$ stands for the nonrelativistic effective mass of a nucleon 111Here we do not distinguish the neutron and proton effective masses, but a straightforward generalization would give $m^{*2}\to m_{n}^{*}m_{p}^{*}$, where the indices $n$ and $p$ refer to neutrons and protons. , the Fermi and Bose functions of nondimensional arguments have the form $f(x)=(e^{x}+1)^{-1}$ and $g(x)=(e^{x}-1)^{-1}$,

[TABLE]

and the integration limits $x_{\rm min}$ and $x_{\rm max}$ are given by $x_{\rm min/max}=(y+\alpha_{e}+\alpha_{\nu}\mp z)/2$. The integration variables $y$ and $z$ are the transferred energy and momentum, respectively, normalized by the temperature, and the variable $x$ is the normalized-by-temperature antineutrino energy computed from the chemical potential, i.e., $x=(\epsilon_{\bar{\nu}}+\mu_{\nu})/T$ (recall that antineutrino chemical potential is $-\mu_{\nu}$). Note that in our rate calculations we numerically evaluate the full three-dimensional integral. We do not use the Fermi surface approximation (assuming that all momenta lie close to the Fermi momentum) because it is no longer valid at the temperatures of interest to us.

When the system is in beta equilibrium, $\mu_{\Delta}=0$ and the rates of the direct and inverse processes are equal $\Gamma_{1n}=\Gamma_{1p}\equiv\Gamma_{1}$ in exact $\beta$-equilibrium. For small departures from $\beta$-equilibrium, $\mu_{\Delta}\ll T$, only the terms that are linear in the departure $\mu_{\Delta}$ are of interest; the coefficients of the expansion involve the derivatives

[TABLE]

where $\Gamma_{1}$ is the rate in $\beta$-equilibrium (as defined above) and is given explicitly by

[TABLE]

Note that if in some density-temperature range neutrinos are trapped and are degenerate, i.e., $\mu_{\nu}\gg T$, then one can approximate $\epsilon_{\nu}\simeq\mu_{\nu}$, $y\simeq 1$, therefore $x-y\simeq 2\mu_{\nu}/T\gg 1$, and the electron Fermi function $f(x-y)$ in Eq. (14) vanishes. If one were to extrapolate the neutrino-trapped rate to low temperature, one would find that $\Gamma_{1}=\lambda_{1}=0$ in this limit.

In the case of neutrino-transparent matter, one should drop the antineutrino distribution $f(x)$ and substitute $\mu_{\nu}=0$ in the neutron-decay rate $\Gamma_{1p}$, while the rate of the inverse process $\Gamma_{1n}$ vanishes. In this case, the $\lambda_{1}$ parameter is given by

[TABLE]

In the limit of strongly degenerate matter ($T\lesssim 1\,{\rm MeV}$ Alford:2018lhf ) we find the following limits for $\Gamma_{1p}$ and $\lambda_{1}$ for the neutrino-transparent case

[TABLE]

where $\alpha=3\left[\pi^{2}\zeta(3)+15\zeta(5)\right]/16\pi^{5}=0.0168$, and $p_{Fi}$ are the Fermi-momenta of the particles, and

[TABLE]

II.2 The rates of the processes $n+\nu_{e}\rightleftarrows p+e^{-}$

The computation of the rates of the processes (2) is carried out in an analogous manner. The rates of the direct and the inverse processes are given, respectively, by

[TABLE]

These rates can be reduced to the following (see Appendix A)

[TABLE]

with $\bar{x}_{\rm min}=(z-y-\alpha_{e}-\alpha_{\nu})/2$. One can verify that $\Gamma_{2p}=\Gamma_{2n}$ when $\mu_{\Delta}=0$. In this case, the $\lambda$-parameter is given by

[TABLE]

where

[TABLE]

If we extrapolate the neutrino-trapped result into the low-temperature regime (where in reality neutrinos are no longer trapped) of highly degenerate matter where the fermionic chemical potentials satisfy the condition $\mu_{i}\gg T$ ($i\in n,p,e,\nu$) the rate $\Gamma_{2}$ given by Eq. (23) reduces to

[TABLE]

Thus, the computation of the parameters $\lambda_{1,2}$, which determine the nonequilibrium relaxation rate of Urca processes at arbitrary degeneracy of the involved fermions, nonzero temperature and in the presence of neutrino trapping reduces to an evaluation of three-dimensional integrals given by Eqs. (14) and (23). These results are essentially exact, the only approximation being the neglect of the terms $O(m^{*}/E)$, where $E=\sqrt{p^{2}+m^{*2}}$, in the tree-level weak-interaction matrix element. Note, however, that many-body correlations in the baryonic matter, which arise from a resummation of particle-hole diagrams are not included yet. In other words, our results correspond to the evaluation of the polarization tensor of baryonic matter in the one-loop approximation.

We consider also the neutrino-transparent case, where we have $\Gamma_{2p}=0$, and for $\lambda_{2}$ we find

[TABLE]

In the limit of strongly degenerate matter ($T\lesssim 1\,{\rm MeV}$ Alford:2018lhf ), the rate $\Gamma_{2n}$ is the same as in the case of $\Gamma_{1p}$

[TABLE]

Similarly, the low-temperature limit for $\lambda_{2}$ reads

[TABLE]

therefore, the total rate $\lambda\equiv\lambda_{1}+\lambda_{2}$ is given by

[TABLE]

which agrees with the corresponding expression given in Ref. Haensel2000 .

III Bulk viscosity

The purpose of this section is to derive a microscopic formula for the bulk viscosity. We will assume that the matter is composed of neutrons, protons, electrons, muons, and neutrinos. Although there are large-amplitude density oscillations in a merger Alford2017 ; Perego:2019adq , we will restrict our analysis to the “subthermal” case where the matter is only slightly perturbed from equilibrium by density oscillations of some characteristic frequency $\omega$. We will not study the “suprathermal” case of high amplitude oscillations, but it can only lead to an enhancement of the bulk viscosity when the beta equilibration rate is slower than the density oscillation frequency, and, as we will show, in the neutrino-trapped temperature range we are in the opposite regime of fast equilibration. In our analysis we take into account the contribution of muons to the thermodynamic quantities, but neglect their contribution to the bulk viscosity, as it is subdominant to the processes involving electrons.

In the case where neutrinos are trapped, the equilibrium with respect to weak interactions implies the conditions (3). The charge neutrality condition implies $n_{p}=n_{e}+n_{\mu}$. These two conditions are sufficient to fix the number densities of the constituents for any given temperature $T$, the value of the baryon number density $n_{B}=n_{n}+n_{p}$ and the lepton number density $n_{L}=n_{e}+n_{\nu}$. Here the neutrino net density is given by ${n}_{\nu}\equiv\tilde{n}_{\nu}-\tilde{n}_{\bar{\nu}}$ where $\tilde{n}_{\nu}$ and $\tilde{n}_{\bar{\nu}}$ are the neutrino and antineutrino number densities, respectively.

Consider now small-amplitude density oscillations in the matter, with characteristic timescales that are long compared to the strong interaction timescale $\sim 10^{-23}$ s. Since the strong interactions establish thermal equilibrium, the particle distributions are always thermal (Fermi-Dirac or Bose-Einstein); the only deviation from equilibrium that is induced by the oscillations is a departure from beta equilibrium, which can be expressed in terms of a single chemical potential $\mu_{\Delta}$ (5).

The perturbed densities are written as $n_{B}(t)=n_{B0}+\delta n_{B}(t)$, and $n_{L}(t)=n_{L0}+\delta n_{L}(t)$, with $n_{B0}$ and $n_{L0}$ being the unperturbed background densities of baryons and leptons. The time-dependence of the density perturbations is taken as $\delta n_{B}(t),\delta n_{L}(t)\sim e^{i\omega t}$. The continuity equation $\partial n_{i}/\partial t+{\rm div}\,n_{i}\bm{v}=0$ then implies

[TABLE]

where $\bm{v}$ is the bulk (hydrodynamic) velocity of matter and $\theta={\rm div}\,\bm{v}$. (Note that we consider only linear perturbation in densities.)

The density perturbations above imply density perturbation of particle number which can be separated into equilibrium and nonequilibrium parts

[TABLE]

where $j=\{n,p,e,\nu\}$ labels the particles. The variations $\delta n^{\rm eq}_{j}(t)$ denote the shift of the equilibrium state for the instantaneous values of the baryon and lepton densities $n_{B}(t)$ and $n_{L}(t)$, whereas $\delta n^{\prime}_{j}(t)$ denote the deviations of the corresponding densities from their equilibrium values. Due to the nonequilibrium shifts $\delta n^{\prime}_{j}(t)$ the composition balance of matter is disturbed leading to a small shift $\mu_{\Delta}(t)=\delta\mu_{n}(t)+\delta\mu_{\nu}(t)-\delta\mu_{p}(t)-\delta\mu_{e}(t)$, which can be written as

[TABLE]

where $s_{i}=+1$ for $n,\,\nu$ and $-1$ for $p,\,e$; $A_{n}=A_{nn}-A_{pn}$, $A_{p}=A_{pp}-A_{np}$, and $A_{e}=A_{ee}$, $A_{\nu}=A_{\nu\nu}$ with

[TABLE]

and index 0 denotes the static equilibrium state. The off-diagonal elements $A_{np}$ and $A_{pn}$ are nonzero because of the cross-species strong interaction between neutrons and protons. Since we treat the electrons and neutrinos as ultrarelativistic noninteracting gas, we have kept only the terms that are diagonal in indices $i,j$, which we will further denote simply as $A_{e}$ and $A_{\nu}$. The computation of susceptibilities $A_{i}$ is performed in Appendix B.

If the weak processes are turned off, then a perturbation conserves all particle numbers, therefore

[TABLE]

Once the weak reactions are turned on, there is an imbalance between the rates of weak processes given by Eqs. (14) and (23). To linear order in $\mu_{\Delta}$ the imbalance is given by Sawyer1989 ; Haensel1992PhRvD ; Huang2010

[TABLE]

where $\Gamma_{j}=\Gamma_{1j}+\Gamma_{2j}$ and $\lambda=\lambda_{1}+\lambda_{2}$ see Eqs. (13) and (22) (in the case of neutrino-transparent matter one should use Eqs. (17) and (27) for $\lambda_{j}$). Then instead of Eq. (33) we will have the following rate equations which take into account the loss and gain of particles via the weak interactions

[TABLE]

We next substitute $\delta n_{j}(t)\sim e^{i\omega t}$, eliminate $\mu_{\Delta}$ using Eq. (31), and employ the relations $\delta n_{p}=\delta n_{B}-\delta n_{n}$, $\delta n_{e}=\delta n_{p}$ and $\delta n_{L}=\delta n_{e}+\delta n_{\nu}$, and Eq. (29) to find

[TABLE]

where

[TABLE]

so $A$ is the “beta-disequilibrium–proton-fraction” susceptibility: it measures how the out-of-beta-equilibrium chemical potential is related to a change in the proton fraction. In order to separate the nonequilibrium parts of $\delta n_{j}$ we need to find also the equilibrium shifts $\delta n_{j}^{\rm eq}$. According to the definition of the $\beta$-equilibrium state we have $\mu^{\rm eq}_{n}(t)+\mu^{\rm eq}_{\nu}(t)-\mu^{\rm eq}_{p}(t)-\mu^{\rm eq}_{e}(t)=0$, therefore

[TABLE]

Using the relations $\delta n^{\rm eq}_{n}+\delta n^{\rm eq}_{p}=\delta n_{B}$, $\delta n^{\rm eq}_{e}=\delta n^{\rm eq}_{p}$, $\delta n^{\rm eq}_{e}+\delta n^{\rm eq}_{\nu}=\delta n_{L}$, and substituting also $\delta n_{B}$ and $\delta n_{L}$ from Eq. (29) we find

[TABLE]

Then, according to Eq. (30), we find for $\delta n^{\prime}_{j}$

[TABLE]

with

[TABLE]

so $C$ is the “beta-disequilibrium–baryon-density” susceptibility: it measures how the out-of-beta-equilibrium chemical potential is related to a change in the baryon density at fixed proton fraction. Now we are in a position to compute the full nonequilibrium pressure which is given by

[TABLE]

where the nonequilibrium part of the pressure, referred to as bulk viscous pressure, is given by

[TABLE]

Using the Gibbs-Duhem relation $dp=sdT+\sum_{j}n_{j}d\mu_{j}$, which is valid also out of equilibrium, we can write 222Note that the temperature is assumed to be constant because, as argued in Sec. I, we assume that the thermal equilibration rate is much larger than the chemical equilibration rate.

[TABLE]

from which we can identify

[TABLE]

where we used the symmetry relation $A_{np}=A_{pn}$. Collecting the results (45), (46), (51) we find that the bulk viscous pressure (49) is given by

[TABLE]

The bulk viscosity is the real part of $-\Pi/\theta$,

[TABLE]

which has the classic resonant form depending on two quantities: the prefactor $C^{2}/A$ which is a ratio of susceptibilities (40),(47), depending only on the EoS, and the relaxation rate $\lambda A$ which depends on the weak interaction rate $\lambda=\lambda_{1}+\lambda_{2}$ (13),(22) and the susceptibility that relates $\mu_{\Delta}$ to the proton fraction. Note that if we extrapolate the neutrino-trapped calculation to the low-temperature, degenerate limit, we can compute the susceptibility $A$ analytically, see Appendix B.

IV Numerical results

To quantify the amount of dissipation through bulk viscosity in the present context we need first to specify the properties of $\beta$-equilibrated nuclear matter. We choose to do so using the density functional theory (DFT) approach to the nuclear matter, which is based on phenomenological baryon-meson Lagrangians of the type proposed by Walecka, Boguta-Bodmer and others glendenning2000compact ; weber_book ; Sedrakian2007 . We will use the parametrization of such a Lagrangian with density-dependent meson-nucleon coupling Lalazissis2005 and will apply the DFT to nuclear matter with trapped neutrinos, see also Colucci2013 .

IV.1 Beta-equilibrated nuclear matter

The Lagrangian density of matter can be written as ${\cal L}={\cal L}_{N}+{\cal L}_{l}$, where the baryonic contribution is given by

[TABLE]

where $N$ sums over nucleons, $\psi_{N}$ are the nucleonic Dirac fields with masses $m_{N}$. The meson fields $\sigma,\omega_{\mu}$, and ${\bm{\rho}}_{\mu}$ mediate the interaction among baryon fields, $\omega_{\mu\nu}$ and ${\bm{\rho}}_{\mu\nu}$ represent the field strength tensors of vector mesons and $m_{\sigma}$, $m_{\omega}$, and $m_{\rho}$ are their masses. The baryon-meson coupling constants are denoted by $g_{iN}$ with $i=\sigma,\omega,\rho$. The leptonic contribution is given by

[TABLE]

where $\lambda$ sums over the leptons $e^{-},\mu^{-}$, $\nu_{e}$ and $\nu_{\mu}$, which are treated as free Dirac fields with masses $m_{\lambda}$; the mass of electron neutrino is negligible and is set to zero. We do not consider electromagnetic fields, therefore their contribution is dropped. The coupling constants in the nucleonic Lagrangian are density-dependent and are parametrized according to the relation $g_{iN}(n_{B})=g_{iN}(n_{0})h_{i}(x)$, for $i=\sigma,\omega$, and $g_{\rho N}(n_{B})=g_{\rho N}(n_{0})\exp[-a_{\rho}(x-1)]$ for the ${\bm{\rho}}_{\mu}$-meson, where $n_{B}$ is the baryon density, $n_{0}$ is the saturation density and $x=n_{B}/n_{0}$. The density dependence of the couplings is encoded in the functions

[TABLE]

This parametrization has in total eight parameters, which are adjusted to reproduce the properties of symmetric and asymmetric nuclear matter, binding energies, charge radii, and neutron radii of spherical nuclei, see Table 1. We recall that the Lagrangian of this model has only linear in meson field interaction terms and the coupling nucleon-meson constants are density-dependent. We will also employ below the NL3 model Lalazissis1997 as an alternative, which has density-independent meson-nucleon couplings but contains nonlinear in meson fields terms.

From the Lagrangian densities (IV.1) and (55) we obtain the pressure of the nucleonic component

[TABLE]

where $m_{N}^{*}=m_{N}-\sigma_{0}g_{\sigma N}$ is the relativistic (Dirac) effective nucleon mass, $\mu^{*}_{N}=\mu_{N}-g_{\omega N}\omega_{0}-g_{\rho N}I_{3}\rho^{0}_{3}$ is the nucleon chemical potential including the time component of the fermion self-energy, $I_{3}$ is the third component of nucleon isospin and $\sigma_{0}$, $\omega_{0}$ and $\rho^{0}_{3}$ are the mean values of the meson fields, $E^{N}_{k}=\sqrt{k^{2}+m^{*2}_{N}}$ are the single particle energies of nucleons. The first three terms in this expression are associated with mean values of the mesonic fields, whereas the last term is the fermionic contribution which is temperature-dependent.

The leptonic contribution to the pressure is given by

[TABLE]

where $g_{\lambda}$ is the leptonic degeneracy factor and $E_{k}^{\lambda}=\sqrt{k^{2}+m_{\lambda}^{2}}$ are the single particle energies of leptons. At nonzero temperature, the net entropy of the matter is the sum of the nucleon contribution

[TABLE]

and the lepton contribution

[TABLE]

where $\bar{f}(y)=1-f(y)$. The energy density of the system and other thermodynamical parameters, for example, the free-energy can be computed in an analogous manner.

In the case of the neutrino-transparent medium, the composition of hadronic matter includes neutrons, protons, electrons and muons. Neutrinos are assumed to escape. We will use the approximation that the beta equilibrium conditions become $\mu_{n}=\mu_{p}+\mu_{e}$ and $\mu_{\mu}=\mu_{e}$. At the temperatures of interest to us there are corrections to these expressions Alford:2018lhf but we will neglect them.

If the neutrinos are trapped in the matter they contribute to the energy density and entropy of matter. We take into account the two lightest flavors of neutrinos, electron and muon neutrinos, and their antineutrinos 333Because we restrict our consideration to the Standard Model, where there are no neutrino oscillations, $\tau$ neutrinos can be ignored because $\tau$ lepton is too heavy to appear in matter.. In this case, the chemical equilibrium conditions read $\mu_{n}+\mu_{\nu_{l}}=\mu_{p}+\mu_{l}$ with $l=\{e,\mu\}$, and the lepton number $n_{L}$ conservation implies $n_{l}+n_{\nu_{l}}=n_{L_{l}}=Y_{L_{l}}n_{B}$, where the lepton fractions $Y_{L_{l}}$ should be fixed for each flavor separately. In our numerical calculations we will consider two cases: (i) $Y_{L}=0.1$ for both flavors, which is typical for matter in binary neutron star mergers; (ii) $Y_{L}\equiv Y_{L_{e}}=0.4$ and $Y_{L_{\mu}}=0$ which are typical for matter in supernovae and proto-neutron stars Prakash1997 ; Malfatti2019 ; Weber2019 .

The neutrino transparent case $Y_{\nu}=0$ which applies for temperatures below the transparency temperature $T_{\rm tr}\simeq 5$ MeV is discussed in Ref. Alford:2019zzz and Appendix D.

Before presenting our results on the bulk viscosity in Sec. IV.2 we first discuss the thermodynamics of the underlying relativistic density functional model of nuclear matter, which will be used as the background equilibrium for our subsequent perturbation analysis. We will employ the DD-ME2 parametrization of the density functional given in Ref. Lalazissis2005 .

Figures 1 and 2 show the particle fractions $Y_{j}=n_{j}/n_{B}$ as functions of the baryon density normalized to the nuclear saturation density, which is $n_{0}=0.152$ fm*-3* in the DD-ME2 model. Figures 1 and 2 refer to the cases of neutron star mergers and supernovae, respectively. The results are shown for two temperatures $T=5$ MeV (which is of the order of $T_{\rm tr}$) [panels (a)] and $T=50$ MeV [panels (b)], which is close to the upper limit of the temperature range achieved in these events. Comparing the panels (a) and (b) in Figs. 1 and 2 we see that the particle fractions are generally not sensitive to the temperature for the given value of $Y_{L}$. The only exception is for low lepton-fraction matter, where at low density and high temperature (Fig. 1(b)), the net electron neutrino density becomes negative, indicating that there are more (electron) antineutrinos than neutrinos. As one would expect, merger (low lepton fraction) matter has much smaller electron neutrino fraction; the electron and muon fractions are $\sim 10\%$, so by charge neutrality $Y_{p}=Y_{e}+Y_{\mu}$ the proton fraction is $\sim 20\%$.

Since the bulk viscosity is related to departure from beta equilibrium, it is instructive to examine the chemical potentials of particles as functions of density and temperature. These are shown as functions of the baryon density in Figs. 3 and 4.

There are two different chemical potentials for baryons: the thermodynamic chemical potentials $\mu_{n}$ and $\mu_{p}$, which enter into the thermodynamic relations and the $\beta$-equilibrium condition (3), and the effective chemical potentials $\mu_{n}^{*}$ and $\mu_{p}^{*}$, which enter into the baryon distribution functions and are defined after Eq. (57). We show also the effective nucleon mass $m^{*}$ with dotted lines.

In low lepton-fraction matter at low density and high temperature (Fig. 3(b)), we see that the neutrino chemical potentials become negative, as expected from the particle fraction results (Fig. 1(b)) which showed that there are more antineutrinos than neutrinos in this regime.

For completeness and reference, we show the equation of state (EoS) of the DD-ME2 model at two temperatures in Fig. 5. For $T=50$ MeV the EoS is shown for two values of lepton fractions corresponding to a supernova and binary-neutron-star merger settings. The increase in the pressure at larger temperatures is due to the additional thermal contribution from baryons and the contribution from trapped neutrinos which is absent in the case when $T=0.1$ MeV.

IV.2 Perturbations and bulk viscosity

Now we turn to the discussion of perturbations on the background equilibrium of matter presented above and concentrate on the $\beta$-relaxation rates and the bulk viscosity.

Our numerical calculations show that the main contribution to the beta-disequilibrium–proton-fraction susceptibility $A$ (40) comes from neutrinos whereas electron contribution is minor. The baryon contributions are much smaller than those of leptons because of their finite mass, see Eq. (133). The contribution from $\rho$ meson is negligible in the whole regime of interest. The susceptibility $A$ does not depend strongly on the density and the temperature and has roughly the same order of magnitude $A\sim 10^{-3}$ MeV*-2* in the relevant portion of the phase diagram.

The beta-disequilibrium–proton-fraction susceptibility $C$ is an increasing function of density and has the same order of magnitude in the cases of neutron star mergers and supernovas. At sufficiently high temperatures $T\gtrsim 30$ MeV $C$ crosses zero at a temperature-dependent critical value of the density, close to saturation density. The vanishing of $C$ arises when the proton fraction in beta-equilibrated matter is independent of the density (passing through a minimum in this case). At the critical density the system is scale-invariant: it can be compressed and remain in beta equilibrium. Thus the bulk viscosity vanishes at the critical density.

Figure 6 shows the ratio $C^{2}/A$ as a function of density. As seen from the figure, this ratio is temperature-sensitive only close to the point were $C$ crosses zero.

In Appendix C we discuss the rates of neutron decay and electron capture that combine to establish beta equilibrium. Electron capture dominates because the neutron decay process involves antineutrinos the population of which is damped by a factor of $\exp(-\mu{{}_{\nu}}/T)$.

In Fig. 7 we show the beta equilibrium relaxation rate $\lambda A$, which determines where the bulk viscosity reaches its resonant maximum [Eq. (53)]. For comparison here we show also the case of neutrino-transparent matter (solid line). The relaxation rate is slowest in the neutrino-transparent case and increases with the lepton fraction in the neutrino-trapped case.

The relaxation rate $\lambda A$ of the neutrino-trapped matter is several orders of magnitude larger than the oscillation frequencies $f=\omega/2\pi\lesssim 10$ kHz typical to neutron star mergers and supernovas. In Fig. 7 the horizontal lines for different oscillations frequencies intersect the $\lambda A$ curves at low temperatures $T\lesssim 0.1$ MeV, indicating that the resonant maximum occurs at low temperatures where the assumption of neutrino trapping is no longer valid. The neutrino-trapped regime lies at higher temperatures, where the bulk viscosity is independent of the oscillation frequency and takes the form $\zeta\approx C^{2}/(\lambda A^{2})$.

In contrast, the neutrino-transparent matter features bulk viscosity which strongly depends on the oscillation frequency, see Ref. Alford:2019zzz and Appendix D for details.

We also see from Fig. 7, that $\lambda A$ is almost independent of the baryon density in the range $0.5\leq n_{B}/n_{0}\leq 3$ for neutrino-trapped matter. At moderately temperatures $\lambda$ scales as $\lambda\propto T^{2}$ for temperatures $T\leq 10$ MeV: this scaling is clearly seen from Eqs. (22) and (24), which are applicable as long as the fermions are semidegenerate, i.e., $T\leq 10$ MeV in the relevant density range.

Figure 8 shows the density dependence of the bulk viscosity for various values of the temperature. Because the beta relaxation rate $\lambda A$ is almost independent of the baryon density, the density dependence of the bulk viscosity follows that of the susceptibility prefactor $C^{2}/A$, and, therefore, as noted above, may drop to zero at a critical density where the system becomes scale-invariant. The critical density is present only at sufficiently high temperatures $T\geq 30$ MeV.

The temperature dependence of the bulk viscosity is shown in Fig. 9. The temperature dependence of $\zeta$ arises mainly from the temperature dependence of the beta relaxation rate $\lambda A\propto T^{2}$ (see Eqs. (22) and (24)), so the bulk viscosity decreases as $\zeta\propto T^{-2}$ in the neutrino-trapped regime, as can be seen also from Fig. 9. This scaling breaks down at special temperatures where the bulk viscosity has zeros when the matter becomes scale-invariant (see our discussion of Fig. 6).

In order to check whether the high-temperature behavior found in Fig. 9 is a universal behavior or is specific to the EoS we used, we compared our results with those obtained in the framework of NL3 model, see Fig. 10. The figure shows that the minimums arise independently of the equation of state, and, therefore, are typical to the high-temperature regime of dense nuclear matter.

Comparing the results shown in panels (a) and (b) of Figs. 8 and 9 we see that the bulk viscosity is generally by a factor of few smaller for larger lepton fractions, which is a consequence of larger values of $\lambda A$ for higher $Y_{L}$, as was seen from Fig. 7. However, the order of magnitude of the bulk viscosity is the same in both cases.

In Figs. 11 and 12 we combine and compare our results for the neutrino-trapped matter with the results for neutrino-transparent matter (see also Alford:2019zzz ). In the interval $5\leqslant T\leqslant 10$ MeV we interpolate the numerical results for the bulk viscosity between the two regimes. We see that the bulk viscosity in the neutrino transparent regime is larger, for two reasons. First, the beta relaxation rate is slower, so the resonant peak of the bulk viscosity occurs within its regime of validity, whereas for neutrino trapped matter the regime of validity starts at temperatures well above the resonant maximum. Second, the prefactor $C^{2}/A$ is larger in the neutrino-transparent matter, so the bulk viscosity reaches a higher value at its resonant maximum.

It, therefore, seems likely that bulk viscosity will have its greatest impact on neutron star mergers in regions of the merger that are neutrino transparent rather than neutrino trapped Alford:2019zzz .

V Conclusions

In this work, we have studied the bulk viscosity of hadronic component of neutron stars composed of neutrons, protons, electrons, muons, and neutrinos. The main new ingredient of our study is the trapped neutrinos, which modify significantly the composition of the background equilibrium matter. We have derived semianalytical expressions for the weak interaction rates in the case of trapped neutrinos and corresponding expressions for the bulk viscosity. Our numerical study of the relevant quantities displays the following features:

(a) Electron capture dominates. In neutrino-trapped matter, beta equilibration, and hence bulk viscosity, is dominated by the electron capture process and its inverse (2). Neutron decay and its inverse (1) involve antineutrinos and are therefore suppressed by factors of $\exp(-\mu_{\nu}/T)$.

(b) Role of susceptibilities. The beta-disequilibrium–baryon-density susceptibility $C$ (47) plays an essential role in the bulk viscosity since it measures the degree to which matter is driven out of beta equilibrium when compressed at constant proton fraction (i.e., without any weak interactions). We find that $C$ vanishes, i.e., the proton fraction becomes independent of density, at a critical density that is in the vicinity of saturation density at high temperatures $T\gtrsim 30$ MeV. The subthermal bulk viscosity vanishes at this critical density because the equilibrium value of the proton fraction is density-independent, so compression does not drive the system out of beta equilibrium.

(c) Temperature dependence. The bulk viscosity as a function of temperature at fixed oscillation frequency $\omega$ shows the standard resonant form (53), with a maximum when the beta relaxation rate $\lambda A$ matches $\omega$. At the temperatures where our assumption of neutrino trapping is valid ($T\gtrsim 5\,{\rm MeV}$), we are always in the regime where beta relaxation is fast ($\lambda A$ is greater than typical frequencies for neutron star density oscillations, which are in the kHz range) so $\zeta\approx C^{2}/(\lambda A^{2})$. Since the relaxation rate rises with temperature (due to increasing phase space), the bulk viscosity drops with increasing temperature as $\zeta\propto T^{-2}$. This scaling can be understood, by noting that the factor $\lambda A$ scales approximately as $T^{-2}$. Given that the remaining factor $C^{2}/A$ mildly depends on the temperature (except the points where it goes to zero, see Fig. 6), we recover the $\zeta\propto T^{-2}$ scaling. At some temperatures and densities the material becomes scale-invariant, so $C$ goes through zero, driving the bulk viscosity to zero at those points.

(d) Dependence on lepton fraction. The dependence of the bulk viscosity on lepton fraction can be inferred from Figs. 8, 9, and 10. It is seen that the bulk viscosity is smaller for larger lepton fraction, specifically comparing the cases of supernova matter with $Y_{L}=0.4$ and neutron-star/binary-merger matter with $Y_{L}=0.1$ one finds that the bulk viscosity of supernova matter is a few times smaller compared to that of neutron-star/binary-merger matter.

(e) Effect of neutrino trapping. In astrophysical scenarios of supernova and binary-mergers, the object formed in the aftermath of these events will cool eventually below the neutrino trapping temperature $T_{\rm tr}$. Therefore, one may ask how the bulk viscosity of matter changes as it passes through $T_{\rm tr}$. Figure 11 shows the temperature dependence of the bulk viscosity, where we combine the results obtained in the neutrino-trapped ($T>T_{\rm tr}$) and neutrino-transparent matter ($T<T_{\rm tr}$). It is seen that the bulk viscosity of matter with neutrinos is several orders of magnitude smaller than that for the matter which is transparent to neutrinos.

(f) Dependence on the density-functional model. All computations have been carried out for two alternative relativistic density functional models one being the density-dependent DD-ME2 parametrization the other being the NL3 nonlinear parametrization. The key results show insignificant dependence on the chosen model, therefore we conclude that our results are largely independent of this input (see also Ref Yakovlev2018 ).

(g) Relevance for mergers. It seems likely that the presence of neutrino-trapped matter will not be an important source of bulk viscosity in mergers. Ref. Alford:2019zzz found that the bulk viscosity in neutrino transparent matter was just enough to yield dissipation times in the 20 ms range. We find that bulk viscosity in neutrino-trapped matter is thousands of times smaller, so the corresponding dissipation times are likely to be too long to affect a merger, provided other input in the analysis of Ref. Alford:2019zzz does not change by orders of magnitude.

In the future it would be interesting to extend our discussion to more complicated compositions of dense matter which would include heavy baryons such as hyperons and delta-isobars Li2018_PLB ; Li2019ApJ . We also hope that this work will help to clarify which dissipative processes are important enough to be worth including in hydrodynamic simulations of neutron star mergers.

Acknowledgments

We thank Steven Harris, Kai Schwenzer, Alex Haber for discussions and Tobias Fischer for output tables of supernova simulations and discussions. M. A. is supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Award No. DE-FG02-05ER41375. A. H. and A. S. acknowledge the partial support of the European COST Action “PHAROS” (CA16214) and the State of Hesse LOEWE-Program in Helmholtz International Center for FAIR. A.S. acknowledges the support by the DFG (Grant No. SE 1836/4-1).

Appendix A Phase space integrals

For further computations it is convenient to write the energy conservation in the form $\delta(p_{0}\pm k_{0}+k^{\prime}_{0}-p^{\prime}_{0})=\delta(\epsilon_{p}\pm\epsilon_{\bar{\nu}/\nu}+\epsilon_{e}-\epsilon_{n}-\mu_{\Delta})$, where we added and subtracted $\mu_{\Delta}$ in the argument of the $\delta$-function, and denoted by $\epsilon_{i}$ the energies of the particles computed from their chemical potentials, e.g., $\epsilon_{p}=p_{0}-\mu_{p}$. To simplify the calculation of the $\beta$-equilibration rates it is useful to introduce a so-called “dummy” integration, so that, by substituting Eq. (7) into the rates (8) and (18) we obtain

[TABLE]

where

[TABLE]

with $\delta^{(4)}(p-p^{\prime}+q)=\delta({\bm{p}}-{\bm{p}}^{\prime}+{\bm{q}})\delta(\epsilon_{p}-\epsilon_{p^{\prime}}+\omega-\mu_{\Delta})$, $\delta^{(4)}(k^{\prime}\pm k-q)=\delta({\bm{k}}^{\prime}\pm{\bm{k}}-{\bm{q}})\delta(\epsilon_{k^{\prime}}\pm\epsilon_{k}-\omega)$. The rates of the inverse processes (9) and (19) can be obtained from Eqs. (61) and (62) by replacing $f(p_{i})\to 1-f(p_{i})$ for all particles. Thus, the problem reduces to the computation of three $q$-dependent integrals $I_{1}(q)$, $I_{2}(q)$ and $I_{3}(q)$ in Eq. (63)-(65).

To compute the integral $I_{1}(q)$ we use the following identity between the Fermi $f(p)$ and Bose $g(p)$ functions

[TABLE]

Then, integrating over neutron momentum and separating the angular part of the remaining integral we obtain ($\tilde{\omega}=\omega-\mu_{\Delta}$)

[TABLE]

where $x$ is the cosine of the angle between $\bm{p}$ and $\bm{q}$. Using the nonrelativistic spectrum for the nucleons $\epsilon_{i}=p^{2}/2m^{*}-\mu_{i}^{*}$, $i=\{p,n\}$, we obtain for $I_{1}(q)$

[TABLE]

where $x_{0}$ is the zero of the argument of the $\delta$-function, i.e.,

[TABLE]

The step-function sets the following limit on the momentum of a particle

[TABLE]

Taking now the momentum integral we finally obtain (note that $g(q)$ depends only on $\omega$)

[TABLE]

where the lower limit $\epsilon_{\rm min}$ follows from $p_{\rm min}$, i.e.,

[TABLE]

We now transform the second integral (64)

[TABLE]

where we neglected the electron mass and introduced the short-hand notation $\omega^{\prime}=\omega+\mu_{e}-\mu_{\nu}$. The argument of the $\delta$-function is zero at

[TABLE]

therefore for the $\delta$-function, we obtain (note that the $\delta$-function also implies $\omega^{\prime}-k\geq 0$)

[TABLE]

and

[TABLE]

The condition $|y_{0}|\leq 1$ sets the following limits on the momentum $k$

[TABLE]

therefore the following two inequalities must be satisfied simultaneously

[TABLE]

Note that because of the condition $\omega^{\prime}\geq k\geq 0$ we have $\omega^{\prime}+q\geq 0$. Then the first condition in Eq. (78) gives $k\geq(\omega^{\prime}-q)/2$. Next, if $\omega^{\prime}\geq q$, the second condition implies $k\leq(\omega^{\prime}+q)/2\leq\omega^{\prime}$. If $\omega^{\prime}\leq q$ instead, then the first condition is satisfied automatically, and the second one gives $k\geq(\omega^{\prime}+q)/2\geq\omega^{\prime}$, which is not allowed. Substituting these results into Eq. (73) we obtain

[TABLE]

The integral $I_{3}(q)$ given by Eq. (65) is transformed as follows

[TABLE]

The argument of the $\delta$-function is zero at

[TABLE]

and for the $\delta$-function we obtain (note that the $\delta$-function also implies $\omega^{\prime}+k\geq 0$)

[TABLE]

therefore

[TABLE]

The condition $|y_{0}|\leq 1$ sets the following limits on the momentum $k$

[TABLE]

therefore the following two inequalities must be satisfied simultaneously

[TABLE]

Because of the condition $\omega^{\prime}+k\geq 0$ we have $\omega^{\prime}+q+2k\geq 0$. Then the second condition in Eq. (85) gives $q\geq\omega^{\prime}$. Next, if $\omega^{\prime}\geq-q$, the first condition implies $k\geq(q-\omega^{\prime})/2\geq-\omega^{\prime}$. If $\omega^{\prime}\leq-q$ instead, then the first condition implies $k\leq(q-\omega^{\prime})/2\leq-\omega^{\prime}$, which is not allowed. Substituting these results into Eq. (80) we obtain

[TABLE]

Combining now Eqs. (61), (62), (71), (79) and (86) we obtain the final formulas for $\Gamma_{1p}$ and $\Gamma_{2p}$ given in the main text by Eqs. (10) and (20). The rates of the inverse processes $\Gamma_{1n}$ and $\Gamma_{2n}$ given by Eqs. (11) and (21) can be obtained in an analogous manner.

For the derivatives of the rates $\Gamma_{1p}$ and $\Gamma_{1n}$ we obtain

[TABLE]

In the same way we obtain for $\Gamma_{2p}$ and $\Gamma_{2n}$

[TABLE]

From these expressions, it is straightforward to obtain Eqs. (13), (15), (22) and (25) of the main text.

A.1 Low-temperature Urca rates

In the case of highly degenerate matter we have $\mu_{i}/T\to\infty$, therefore $\epsilon_{\rm min}/T\to\pm\infty$. Thus we find from Eq. (71) (for $\mu_{\Delta}=0$)

[TABLE]

In terms of momenta the condition $\epsilon_{\rm min}\leq 0$ can be written as $|p_{Fn}^{2}-p_{Fp}^{2}-q^{2}|\leq 2qp_{Fp}$, where we neglected $\omega\sim T$ terms, therefore

[TABLE]

In the case of neutrino-transparent matter, we can also set $q=p_{Fe}$, therefore for Eq. (91) we obtain

[TABLE]

To obtain the low-temperature limit of integral $I_{2}$ in the case of neutrino-transparent matter we drop the neutrino momentum and neutrino distribution and approximate $|\bm{k}^{\prime}|=p_{Fe}$ in the first equation of (73), which gives

[TABLE]

Performing the integrations over $k^{\prime}$ and $\Omega_{k}$ we obtain

[TABLE]

Substituting Eqs. (93) and (95) into the neutron decay rate (61) we find

[TABLE]

In the same manner we can obtain the low-$T$ result for $\Gamma_{2n}$ in the neutrino-transparent matter

[TABLE]

The integrals appearing in Eqs. (96) and (97) can be computed successively

[TABLE]

and

[TABLE]

from which we obtain the low-$T$ results (16) and (26) of the main text.

Next, we find the low-temperature limit of Eqs. (15) and (25). By replacing $\omega\to\omega-\mu_{\Delta}$ in Eq. (93) we find

[TABLE]

Therefore, the low-temperature limit for $\lambda_{1}$ can be found by replacing $y\to\{y[1+g(y)]-1\}/T$ in Eq. (96). Thus

[TABLE]

Similarly, $\lambda_{2}$ can be obtained from Eq. (97) by replacing $y\to-[yg(y)-1]/T$

[TABLE]

The double integrals in Eqs. (101) and (102) are identical and are equal to $17\pi^{4}/120$. Substituting this value we obtain the final results given by Eqs. (17) and (27).

In the case where neutrinos remain trapped in the degenerate regime, one finds $I_{2}(q)=0$, because in this limit the distribution function of antineutrinos is exponentially suppressed. To obtain the low-temperature limit of $I_{3}(q)$ we substitute $k=p_{F\nu}$ and $\omega^{\prime}=\mu_{e}-\mu_{\nu}=p_{Fe}-p_{F\nu}$ in Eq. (86) ($\epsilon_{k}\equiv k-p_{F\nu}$) to find

[TABLE]

The inequalities (85) in this approximation are independent of $\epsilon_{k}$ and imply

[TABLE]

Combining now Eqs. (62), (91) and (103) and approximating $p_{F\nu}/T\to\infty$ for $\Gamma_{2}$ we obtain

[TABLE]

The last two integrals give $2\pi^{2}/3$, and for the $q$-integral we have

[TABLE]

In neutron star matter we have typically $p_{Fn}+p_{Fp}\geq p_{Fe}+p_{F\nu}\geq|p_{Fn}-p_{Fp}|\geq|p_{Fe}-p_{F\nu}|$, and we obtain the final result given by Eq. (24).

Appendix B Computation of susceptibilities $A_{j}$

To compute the susceptibilities $A_{ij}$ given by Eq. (32) we use the following formula for the particle densities

[TABLE]

where $g_{i}$ is the degeneracy factor, and $f(p)$ and $\bar{f}(p)$ are the distribution functions for particles and antiparticles, respectively. For neutrons, protons, and electrons we have $g_{i}=2$, and for neutrinos $g_{\nu}=1$.

Differentiating the left and right sides of Eq. (107) with respect to $n_{j}$ and exploiting the expressions

[TABLE]

in the case of baryons we obtain

[TABLE]

where

[TABLE]

The average values of the meson fields are given by Chatterjee2007

[TABLE]

which gives (recall that $\mu^{*}_{i}=\mu_{i}-g_{\omega}\omega_{0}-g_{\rho}\rho_{03}I_{3i}$)

[TABLE]

The scalar field is given by

[TABLE]

with $U(\sigma)$ being the self-interaction potential of the scalar field, therefore

[TABLE]

The last term is suppressed in the nonrelativistic limit and can be neglected, after which we obtain

[TABLE]

Substituting this into Eq. (109) we obtain the following equations for coefficients $B_{ij}$

[TABLE]

where

[TABLE]

In the case of $i\neq j$ we find from Eq. (116)

[TABLE]

Substituting these expressions into Eq. (116) for $i=j$ we obtain

[TABLE]

and

[TABLE]

In the nonrelativistic limit we will use the expansion $m/\sqrt{m^{2}+p^{2}}\simeq 1-p^{2}/2m^{2}$ in the integrals (110). We will further drop the contribution of antiparticles because it is not important for the regime of interest. Then ${I}_{0i}^{+}={I}_{0i}^{-}\simeq\tilde{I}_{2i}$, ${I}_{1i}^{+}={I}_{1i}^{-}\simeq\tilde{I}_{2i}-\tilde{I}_{4i}/2m^{*2}$, and ${I}_{2i}^{+}={I}_{2i}^{-}\simeq\tilde{I}_{2i}-\tilde{I}_{4i}/m^{*2}$, where

[TABLE]

Then in the nonrelativistic limit we find for Eq. (117)

[TABLE]

and

[TABLE]

Then

[TABLE]

We next obtain the two combinations relevant to the bulk viscosity

[TABLE]

Substituting the expression for $\beta$ from Eq. (117) and recalling Eq. (112) we obtain

[TABLE]

For leptons we have simply

[TABLE]

where the lepton energies in the integral are taken as $\epsilon_{p}=p-\mu_{L}$. Then

[TABLE]

and

[TABLE]

Taking into account also Eq. (111) and the nonrelativistic limit of Eq. (113) we obtain

[TABLE]

The terms containing $U(\sigma)$ vanish in the case of DD-ME2 model and are numerically very small in the case of NL3 model. We find also, that the terms $\propto g_{\rho}$ in Eqs. (130) and (132) are negligible in comparison to the first four terms. The last term in Eq. (132) is comparable to the rest of the terms.

In the case of degenerate matter the susceptibilities can be computed analytically

[TABLE]

which agree with the results of Refs. Chatterjee2007 ; Chatterjee2008ApJ .

Appendix C Beta equilibration rates

In this appendix we discuss the numerical results for the $\beta$-equilibration rates $\Gamma_{1}$ and $\Gamma_{2}$ given by Eqs. (14) and (23). Figures 13 and 14 show the temperature dependence of the neutron decay and the electron capture rates, respectively, for various values of the density and lepton fraction.

Figures 13 and 14 demonstrate that the electron capture process and its inverse (2) dominate the beta equilibration and hence the bulk viscosity. The electron capture rate $\Gamma_{2}$ is always many orders of magnitude larger than the neutron decay rate $\Gamma_{1}$. This is because we are studying neutrino-trapped matter, where the trapped species is typically neutrinos rather than antineutrinos (Figs. 1, 2) so any process involving antineutrinos will be suppressed by factors of $\exp(-\mu{{}_{\nu}}/T)$.

As seen from the upper panels of the figures, $\Gamma_{1}$ is a rapidly increasing function of the temperature and exponentially vanishes at low temperatures because of vanishing antineutrino density in the degenerate matter. The threshold of temperature below which $\Gamma_{1}$ practically vanishes is located at higher temperatures for higher lepton fractions, because the suppression of the antineutrino density is stronger for larger $Y_{L}$.

The temperature dependence of the electron capture rate differs significantly from that of $\Gamma_{1}$. Indeed, $\Gamma_{2}$ increases with the temperature according to a power-law $\Gamma_{2}\propto T^{3}$ up to the temperatures $T\simeq 10$ MeV, as seen from the low-temperature limit given by Eq. (24). This scaling breaks down at higher temperatures $T\geq 10$ MeV and sufficiently low densities $n_{B}\leq n_{0}$, where the finite temperature effects become important in the evaluation of the integral (23). We have checked numerically that the exact result for $\Gamma_{2}$ (23) tends to its low-temperature limit (24) as $T\leq 1$ MeV.

Comparing the upper and lower panels in Figs. 13 and 14, we see that $\Gamma_{1}$ is smaller for larger lepton fraction, whereas $\Gamma_{2}$ shows the opposite behavior. The reason for this behavior is clear: the neutron decay rate $\Gamma_{1}$ is proportional to the antineutrino (number) density, whereas the electron capture rate $\Gamma_{2}$ is proportional to the neutrino density. Because the electron neutrino fraction increases with the increase of $Y_{L}$, as was seen from Figs. 1 and 2, the antineutrino population becomes more suppressed at higher $Y_{L}$, thus leading to smaller $\Gamma_{1}$ and larger $\Gamma_{2}$.

Appendix D Bulk viscosity of neutrino-transparent matter

In this appendix we are interested in the domain of temperatures $T\leq T_{\rm tr}$ where the matter is neutrino transparent ($Y_{\nu}=0$). Below the temperature $T\simeq 1\,{\rm MeV}$, the modified Urca process becomes important Alford:2018lhf therefore our results are strictly relevant for the temperature domain $1\leq T\leq 5$ MeV. To account for uncertainty in the value of $T_{\rm tr}$ the numerical results will be shown up to $T=10$ MeV. The particle fractions are shown in Fig. 15. In this case muons appear only above a certain baryon density $n_{B}\gtrsim n_{0}$, where the condition $\mu_{e}\geq m_{\mu}\simeq 106$ MeV is satisfied. Below this threshold, the proton and electron fractions are equal, as required by the charge neutrality condition, whereas above the threshold the condition $Y_{p}=Y_{e}+Y_{\mu}$ is satisfied.

The relevant $\beta$-equilibration rates $\Gamma_{1p}$ and $\Gamma_{2n}$ are shown in Fig. 16 and as functions of the temperature. Both quantities rapidly increase with the temperature and are exponentially damped in the low-temperature limit. The reason is that the condition $p_{Fp}+p_{Fe}\geq p_{Fn}$ is never satisfied for the given model of hadronic matter because of very small proton fraction $Y_{p}\leq 7$ %, see Fig. 15. As a consequence, the direct Urca processes for neutrino-transparent matter are always blocked at low temperatures, therefore in that regime, the modified Urca processes should be accurately taken into account Alford:2018lhf .

Quantitatively, the neutron decay rate $\Gamma_{1p}$ in the case of $Y_{\nu}=0$ is larger than in the case of $Y_{L}=0.1$, whereas the electron capture rate is smaller. This result can be anticipated from Eqs. (14) and (23), where one should substitute $\alpha_{\nu}=0$ in the neutrino-transparent case. However, the rate $\Gamma_{1p}$ again remains much smaller than $\Gamma_{2n}$, and only at sufficiently high densities approaches $\Gamma_{2n}$, as seen in Fig. 17, since both quantities have the same low-temperature limit given by Eqs. (16) and (26).

We remark also, that the behavior of $\lambda_{1}$ and $\lambda_{2}$ is very similar to that of $\Gamma_{1p}$ and $\Gamma_{2n}$, respectively. At densities $n_{B}\geq n_{0}$ we have approximately $\lambda_{1}\simeq\Gamma_{1p}/T$ and $\lambda_{2}\simeq\Gamma_{2n}/T$, as it was the case of neutrino-trapped matter, see Eqs. (13) and (22).

Next, we show the density dependence of the bulk viscosity in Fig. 18. As in the case of $Y_{\nu}\neq 0$, the frequency dependence of $\zeta$ can be neglected at sufficiently high temperatures $T\simeq 10$ MeV (see Fig. 7), and, because the product $\lambda A$ is almost density-independent, the bulk viscosity as a function of density increases as $C^{2}$. For smaller temperatures $T\lesssim 5$ MeV the frequency dependence of $\zeta$ becomes important, and the bulk viscosity drops rapidly with $\omega$. We show the bulk viscosity for three values of the frequency in Fig. 18 for $T=5$ MeV. The general increase of $\zeta$ with the density, in this case, is again caused by the factor $C^{2}$, but the hump structures of $\zeta$ in the range arise from the weak density-dependence of the product $\lambda A$. In the case of $T=1$ MeV we have already the opposite limit $\omega\gg\lambda A$ (see Fig. 7), therefore the bulk viscosity scales as $\zeta\propto\lambda C^{2}$, which rapidly increases with the density and drops to zero at sufficiently low densities $n_{B}\lesssim 2n_{0}$. We have checked that our low-temperature result for $T=0.1$ MeV agrees with the result of the bulk viscosity shown in Fig. 2 of Ref. Haensel2000 obtained for direct Urca processes for the same state of hadronic matter.

For comparison, we show also the bulk viscosity of the neutrino-transparent matter for the nuclear model NL3 in Fig. 19. The results obtained within two models DD-ME2 and NL3 differ mainly in the low-temperature region, where the bulk viscosity is strongly suppressed at low densities because of blocking of direct Urca processes. The suppression sets in at lower densities in the case of the NL3 model because this model predicts larger proton and electron fractions than the DD-ME2 model.

Figures 20 and 21 show the temperature dependence of the bulk viscosity for models DD-ME2 and NL3, respectively. It has a maximum at the temperature where $\omega=\lambda A$, and the slope of the curve on the right side of the maximum is larger than in the case of $Y_{\nu}\neq 0$ because of stronger $\lambda=\lambda(T)$ dependence. Using the approximate scaling $\lambda\propto T^{4}$, we find that $T_{\rm max}=(\omega/\lambda_{0}A)^{1/4}$, where $\lambda_{0}$ is the value of $\lambda(T)$ at $T=1$ MeV, i.e., the maximum shifts to higher temperatures for higher frequencies. As seen from a comparison of Figs. 8 and 9 with Figs. 18 and 20, the bulk viscosity is larger in the case of neutrino-transparent matter, as expected.

We also compare our results with the high-frequency (low-temperature) result of Ref. Haensel1992PhRvD obtained from direct Urca processes for the neutrino-transparent matter. Choosing the parameters, e.g., $T=1$ MeV, $f=1$ kHz, $n_{B}=2\,n_{0}$, we find from Figs. 1 and 2 and Eq. (13) of Ref. Haensel1992PhRvD the values $\zeta\simeq 3.5\cdot 10^{28}$ g cm*-1* s*-1* for the model I, and $\zeta\simeq 10^{29}$ g cm*-1* s*-1* for model II, whereas our calculations give much lower result $\zeta\simeq 10^{21}$ g cm*-1* s*-1*. The reason for this is the lower proton fraction in our model, which leads to blocking of direct Urca processes at low temperatures and, therefore, to lower bulk viscosity.