Bulk Viscous Damping of Density Oscillations in Neutron Star Mergers

In this paper, we discuss the damping of density oscillations in dense nuclear matter in the temperature range relevant to neutron star mergers. This damping is due to bulk viscosity arising from the weak interaction ``Urca'' processes of neutron decay and electron capture. The nuclear matter is modelled in the relativistic density functional approach. The bulk viscosity reaches a resonant maximum close to the neutrino trapping temperature, then drops rapidly as temperature rises into the range where neutrinos are trapped in neutron stars. We investigate the bulk viscous dissipation timescales in a post-merger object and identify regimes where these timescales are as short as the characteristic timescale $\sim$10 ms, and, therefore, might affect the evolution of the post-merger object. Our analysis indicates that bulk viscous damping would be important at not too high temperatures of the order of a few MeV and densities up to a few times saturation density.


Introduction
The recent detections of gravitational waves by the LIGO-Virgo collaboration, in particular, the multimessenger binary-neutron star (BNS) merger event GW170817 [1], motivate studies of the transport properties of dense nuclear matter at temperatures and densities relevant to BNS mergers [2][3][4][5][6][7]. The mass of the post-merger object typically would exceed the maximum mass of a neutron star and, as a consequence, it would collapse to a black hole on the timescales ranging from tens of milliseconds up to seconds depending on the mass of the post-merger object [8][9][10]. While gravitational waves in the post-merger phase have not been observed in the GW170817 event due to lack of detector sensitivity at high frequencies, numerical relativity studies of BNS mergers in their highly non-linear regime predict intense emission of gravitational waves in the kHz frequency range during the initial phase of post-merger phase lasting typically 10 ms (for recent simulations see, for example, [11][12][13][14]). Improvements at least by a factor of 2 compared to advanced LIGO design sensitivity are necessary to measure the dominant frequency component of the signal of GW170817-like event and by factors of 4-5 to observe sub-dominant features of post-merger signal [15]. It is expected that dissipation or damping of matter flows in the merged stars could influence the gravitational waves emitted during the post-merger phase. Recent estimates of the role of the thermal conduction and shear and bulk viscosities indicate that, of these, damping of density oscillations via bulk viscosity has the strongest influence [2,6,7]. After a brief introduction to the problem of computation of the bulk viscosity, we present an extension of our recent work [7] on the bulk viscosity of nucleonic matter which includes the estimates of timescales of the damping of the oscillations by the bulk viscosity. The relevance of the bulk viscosity will be assessed by comparing the damping timescales of density oscillations to the characteristic timescales of the initial phase of post-merger ∼ 10 ms (over which the post-merger object is expected to emit intense gravitational waves) as well as the longer-term phase ∼ 1 s.
The bulk viscosity of nuclear matter at temperatures up to about 1 MeV has been studied extensively [16][17][18][19][20][21][22][23][24][25][26] in the context of oscillations of neutron stars and, in particular, as a source of damping of (unstable) r-mode oscillations; for reviews see [27,28]. More recently, interest in neutron star mergers has motivated studies of bulk viscosity of nucleonic matter at temperatures up to several tens of MeV, covering both the regime where neutrinos escape from neutron stars and the regime where they are trapped [6,7]. Such high temperatures significantly affect the phase space occupation of the fermions and, therefore, the rates of the weak-interaction processes [29]. Once the temperature becomes high enough for neutrinos to be trapped, they affect the composition of matter and ensure that direct Urca processes are always kinematically possible, and modified Urca is a subleading correction.
In this contribution, we extend recent work [7] to compute the timescales associated with bulk-viscous damping of oscillations. We use two models of the equation of state (EoS) and associated composition of dense matter which are based on the relativistic density functional theory of nuclear matter. The overall picture is that, at densities from around n 0 (nuclear saturation density) to around 3n 0 , the bulk viscosity reaches a maximum at temperature T 2 ÷ 6 MeV, which is in the regime where neutrinos are not (or not completely) trapped. At these temperatures, therefore, the damping timescale is at a minimum, with values ranging from about 10 ms at low densities n B ≤ n 0 down to milliseconds (or even tenths of milliseconds depending on the EoS) at n B 3n 0 . This means that bulk viscous damping can have noticeable effects during the ∼ 10 ms of initial (gravitational-wave-emission) phase of the post-merger. At higher temperatures where the neutrino-trapping occurs the bulk viscosity falls by orders of magnitude, which implies much longer damping timescales which are larger than the characteristic timescales involved. This paper is organized as follows. In Section 2 we discuss the formalism for computing the bulk viscosity and the approximations involved in such a computation. Our focus is on the beta equilibration processes of neutron decay and electron capture and the microscopic relaxation rates associated with these processes. Section 3 starts with a brief discussion of the thermodynamic properties of nuclear matter derived from density functional theory in Section 3.1. Bulk viscosity and the oscillation damping timescale are discussed in Sections 3.2 and 3.3. Our main results are summarized in Section 4. We use natural units (h = c = k B = 1) and the metric signature g µν = diag(1, −1, −1, −1).

Urca processes and bulk viscosity
We start with a brief reminder of the bulk viscosity of nuclear matter composed of neutrons, protons, electrons, muons; for more details see [6,7]. For simplicity, we will neglect the muonic contribution to equilibration rates, although we include their contribution to static thermodynamic quantities such as susceptibilities.
Above the trapping temperature T tr 5 MeV, the neutrino mean-free-path is smaller than the size of a neutron star, so neutrinos are trapped in the merger region. Under these conditions beta equilibrium is established via neutron decay and electron capture, and their inverse processes In β-equilibrium the chemical potentials of particles obey the relation The particle fractions of baryonic matter for any given temperature T, baryon number density n B = n n + n p and lepton densities n L l = n l + n ν l = Y L l n B (since we ignore muon reactions, the lepton fractions Y L l need to be fixed for each flavor separately) are found by imposing the beta-equilibrium condition (3), a similar condition for muons µ p + µ µ = µ n + µ ν µ and the charge neutrality condition n p = n e + n µ . At lower temperatures T T tr the neutrino mean free path is larger than the size of a neutron star so the system is neutrino-transparent. Neutrinos cannot occur in initial states, therefore the reactions (1) and (2) proceed only in one direction (from left to right). To determine the composition of matter in this case we use the ordinary zero-temperature β-equilibrium conditions µ n = µ p + µ e and µ µ = µ e . Reference [29] found that there are significant corrections to these conditions at T 1 MeV; nevertheless the bulk viscosity in the neutrino transparent regime is not affected significantly [6].
If the matter is driven out of β-equilibrium, for example by compression and rarefaction, the left and right-hand sides of Equation (3) do not balance anymore. The deviation from β-equilibrium is then measured by a quantity As a result, the Urca processes (1) and (2) will go faster in one direction than in the other until the beta equilibrium of matter is restored. Consider now small-amplitude density oscillations in nuclear matter with a frequency ω. The baryon and lepton conservation implies for periodic perturbations δn where n B0 = n n0 + n p0 and n L0 = n e0 + n ν0 are the unperturbed background densities of baryons and leptons, and θ is the divergence of fluid velocity. The compression and rarefaction of matter implies perturbations in particle densities which can be separated into instantaneous equilibrium and non-equilibrium parts where n j0 are the static values of particle densities. The variations δn eq j (t) stand for the shifts of the equilibrium state for the instantaneous values of n B (t) and n L (t), whereas δn j (t) are the deviations of the particle densities from those equilibrium values. There exist two choices of the instantaneous equilibrium state. Below, we follow our recent work [7], for the alternative see Ref. [30]. We compare below these two approaches and explain why they give the same result for the bulk viscosity.
The non-equilibrium perturbations δn j (t) drag matter out of chemical equilibrium by leading to a small chemical potential shift (4) which can be written in terms of particle densities as where The off-diagonal elements A np and A pn are non-zero because of the cross-species strong interaction between neutrons and protons. For small amplitude density oscillations we only need to evaluate the derivatives in Equation (8) at µ ∆ = 0. If there were no flavor-changing weak processes, then the particle densities would just oscillate around their static equilibrium values according to The weak interactions lead to an imbalance between the rates of direct and inverse Urca processes which in the "subthermal regime" (µ ∆ T) can be written as [18,19,30] with Γ p and Γ n being the production rates of protons and neutrons, respectively. The production rate (10) should be added to the right hand sides of Equation (9) with a plus sign for p, e and a minus sign for n, ν, for example, Substituting here Equation (7), exploiting the conditions δn B = δn n + δn p , δn p = δn e + δn µ and δn L = δn e + δn ν , and using Equation (5) for n B , n L and an analogous equation for n µ (as muons are assumed not to participate in any reactions, their fraction is conserved) one finds with a "beta-disequilibrium-proton-fraction" susceptibility given by Equations (12)- (15) are the extensions of Eqs. (37)-(39) of Ref. [7] as we included non-zero muon density n µ here, which was previously neglected. The final formula for the bulk viscosity, however, remains the same after this addition.
In the next step we find δn eq j using the definition of the instantaneous β-equilibrium state: A n δn eq n + A ν δn eq ν − A p δn eq p − A e δn eq e = 0, which gives Note that these expressions are the solutions of the balance equations (11) in the limit of infinite relaxation rate λ → ∞, which implies necessarily µ ∆ → 0. One can check this also by pushing λ to infinity directly in general solutions (12)-(15) which will reduce then to Equations (17)- (20). Now the non-equilibrium density perturbations can be found according to Equation (6) δn where γ = λA has a dimension of frequency and measures the relaxation rate of particle densities to their equilibrium values, and is the "beta-disequilibrium-baryon-density" susceptibility, with Y n = n n /n B being the neutron fraction. The non-equilibrium part of the pressure -the so-called bulk viscous pressure, can be now computed as where we used the Gibbs-Duhem relation dp = sdT + ∑ l n l dµ l . Substituting the solutions (21) in Equation (23) one finds the bulk viscosity from the definition Π = −ζθ The susceptibility prefactor C 2 /A is a pure thermodynamic quantity and depends only on the EoS, whereas the relaxation rate γ = λA [Equation (10)] depends on the microscopic scattering amplitudes of weak interactions. It is seen from Equations (10) and (16) that γ is actually the derivative Thus, the quantity γ measures how the proton net production rate increases when the neutron fraction increases at fixed baryon density, i.e., γ measures how fast the system reacts to a change in the chemical composition of matter. The quantity γ −1 has a dimension of time and can be interpreted as a relaxation time of the system to its beta-equilibrium state. The bulk viscosity in Equation (24) has the classic resonant form and for density oscillations of a given frequency, ω attains a resonant maximum at the temperature where the relaxation rate matches the oscillation frequency, γ(T) = ω.
The value of the bulk viscosity at that maximum is which is independent of the microscopic interaction rates.
In the regime of slow equilibration, where γ ω, the bulk viscosity takes the form In the fast equilibration regime γ ω, and the bulk viscosity, in this case, reduces to The physical reason for this resonant maximum is easy to understand. In the limit where the relaxation rate is much smaller than the oscillation frequency there would effectively be an additional conserved quantity since both proton number and neutron number would be conserved. The proton fraction would be independent of density, so density oscillations would not drive the system out of chemical equilibrium. There would then be no bulk viscosity.
In the opposite limit of fast equilibration or slow density oscillations, γ ω, weak interactions are able to restore the chemical equilibrium of matter on timescales much smaller than the oscillation period. This means that the matter is practically always beta-equilibrated while undergoing compression and rarefaction and, therefore, will not experience any bulk viscosity.
Above we computed the bulk viscous pressure using its standard definition where δP is the shift of the pressure from its static equilibrium value P 0 = P(n j0 ) for arbitrary perturbations δn j , and δP eq is the instantaneous shift of the equilibrium pressure, which depends on δn eq j and does not contribute to the bulk viscous pressure.
According to the two limiting cases of vanishing bulk viscosity discussed above the bulk viscous pressure can be defined also in an alternative way where δP 0 is the instantaneous shift of the pressure in a certain non-equilibrium state with conserved particle fractions which corresponds to the limit λ → 0. According to Equation (30) we can take as alternatives of the beta-equilibrium shifts δn eq j (t) the shifts δn 0 j (t) given by Equation (9), as these are the solutions of exact balance equations (11) in the limit λ → 0. This was just the choice of the equilibrium state in Ref. [30]. Then instead of Equation (21) we will have Computing the bulk viscous pressure from Equation (23) we obtain the same result for ζ, as expected.

Urca process rates
As it was shown in Ref. [7] the equilibration rate λ of the neutrino-trapped matter is dominated by the electron capture process because the neutron decay rate is exponentially damped at low temperatures as ∼ exp(−µ ν /T). The microscopic β-equilibration rate of the electron capture process is given by [7] GeV −2 is the Fermi coupling constant, θ c is the Cabibbo angle with cos θ c = 0.974, and g A = 1.26 is the axial-vector coupling constant; the index j runs over the four participating particles j = {n, p, e, ν}, f (p j ) are the Fermi distributions of particles andf (p j ) = 1 − f (p j ). The inverse process Γ p will be given by an analogous expression by replacing all f j functions withf j .
For small departures from β-equilibrium µ ∆ T the imbalance between the direct and inverse rates can be linearized in µ ∆ with the coefficient of the linear expansion given by [7] where m * is the effective nucleon mass, α j = µ * j /T, µ * j are the effective chemical potentials of particles, f (x) = (e x + 1) −1 and g(x) = (e x − 1) −1 are the Fermi and the Bose distributions, respectively, In the case of neutrino-transparent matter the equilibration rate λ is given by a similar expression, see Refs. [6,7] for details. The low-temperature limit of λ in neutrino-trapped matter is given by [7] and in the neutrino-transparent matter [6,7,20] where p Fj are the Fermi-momenta of the particles. The θ-function in Equation (36) blocks the direct Urca processes at low densities where the proton and electron Fermi momenta are not sufficiently large to guarantee the momentum conservation. In the case where neutrinos are trapped in matter the momentum conservation can be always satisfied for certain particle momenta and the rate is always finite.

Numerical results
We start with the thermodynamics of nuclear matter, which is derived from the relativistic density functional theory based on phenomenological baryon-meson Lagrangians, for reviews see [31,32]. We use the density-dependent baryon-meson coupling model DDME2 [33] as applied to finite-temperature nucleonic matter, see for details Ref. [34][35][36].

Thermodynamics of nuclear matter
The Lagrangian density of the model considered reads where N sums over nucleons, λ -over the leptons, and ψ i are the fermionic Dirac fields with masses m i . The meson fields σ, ω µ , and ρ µ are the effective mediators of strong interaction between baryons, ω µν = ∂ µ ω ν − ∂ ν ω µ and ρ µν = ∂ µ ρ ν − ∂ ν ρ µ are the field strength tensors of ω µ and ρ µ mesons, respectively, m σ , m ω , and m ρ are meson masses, and g i are the baryon-meson coupling constants with i = σ, ω, ρ.
The density dependence of the particle fractions Y j = n j /n B is shown in Fig. 1. The left panel is for a low temperature T = 1 MeV where the system is neutrino-transparent, and the right panel shows the results for neutrino-trapped matter at temperature T = 50 MeV for lepton fraction fixed at Y L = 0.1. In the first case muons appear only above a certain baryon density n B n 0 , where the condition µ e ≥ m µ 106 MeV is satisfied, whereas in the neutrino-trapped case the muons' threshold disappears.
Within the framework of the model above the susceptibilities A and C in the non-relativistic limit for nucleons are given by (see Ref. [7] for details) ν -t r a n s . where σ is the σ meson mean field and ρ 03 is the ρ meson mean field which is non-zero in asymmetric nuclear matter. These expressions for susceptibilities are derived for isothermal density perturbations, which is the case only if the thermal conduction is fast enough to smoothen the temperature gradients during one period of oscillation [2,7]. This might happen, e.g., in the presence of turbulent flows in the merger region which could generate temperature and density variations on distance scales of the order of a few hundred meters. An order-of-magnitude estimate of thermal relaxation timescale is given in Ref. [2], which is 1 sec ×(z typ /km) 2 (T/10 MeV) 2 , where z typ is the typical length scale of thermal gradients. Assuming z typ 100 m and temperatures 1 ÷ 10 MeV (this is the temperature range where the bulk viscosity is relevant to mergers, see below), the thermal relaxation time will lie in the interval 0.1-10 ms, which is below the characteristic timescale of binary neutron star mergers. Thus, for thermal gradients on this distance scale the assumption of isothermal matter is the relevant one. On the scales over which thermal conduction is inefficient the matter should be treated as iso-entropic. The isothermal and adiabatic susceptibilities differ at most by a factor of 2 in the relevant density and temperature range, see Ref. [6] for further details. Figure 2 shows the ratio C 2 /A of susceptibilities as a function of density for three values of the temperature. The susceptibility C is an increasing function of density. At sufficiently high temperatures T 30 MeV, it is negative at low density and crosses zero at a temperature-dependent critical density where the proton fraction reaches a minimum as a function of the density. At that point, the system becomes scale-invariant, as the compression does not drive the matter out of beta equilibrium, and the bulk viscosity vanishes at that critical point.
The ratio C 2 /A grows rapidly with the density in both cases of neutrino-transparent and neutrino-trapped matter and is sensitive to the temperature only close to the point where C = 0. We see that C 2 /A is approximately an order of magnitude smaller in the neutrino-trapped matter the reason being much larger values of A dominated by the contribution of neutrinos as compared to neutrino-transparent matter. Figure 3 shows the relaxation rate γ = λA [Equations (10) and (25)] as a function of temperature for various densities. The equilibration rate λ of the neutrino-trapped matter is dominated by the electron capture process because the neutron decay rate is exponentially damped at low temperatures as λ 1 ∼ exp(−µ ν /T), whereas λ 2 has approximately a quadratic increase with the temperature as suggested by Equation (35). It is seen from Fig. 3(b) that the relaxation rate γ of the neutrino-trapped matter is several orders of magnitude larger than the oscillation frequencies f = ω/2π 1 kHz which are typical to neutron star mergers. This means that the neutrino-trapped matter is always in the fast equilibration regime where the bulk viscosity is independent of the oscillation frequency and is given by Equation (28).

Beta relaxation rates and bulk viscosity
In the neutrino-transparent regime, in contrast, the relaxation rate is much slower, reaching values in the kHz range where it can resonate with typical density oscillations in mergers. We see this in Figure 3(a) where the relaxation rate γ crosses the ω = 2π kHz (corresponding to f = 1 kHz) line at temperatures 4 ÷ 5 MeV indicating that the neutrino-transparent matter possesses a resonant maximum at those temperatures, as it was found also in Refs. [2,6].
The density dependence of the bulk viscosity is shown in Fig. 4. The oscillation frequency is fixed at f = 1 kHz in the case of neutrino-transparent matter whereas the neutrino-trapped matter features a frequency-independent bulk viscosity, as discussed above. The bulk viscosity of neutrino-transparent matter mainly increases with the density, the increase being faster at low temperatures T 3 MeV where ω γ, as seen from Fig. 3. The density dependence of the bulk viscosity of neutrino-trapped matter mainly follows that of the susceptibility C 2 /A because γ depends weakly on the baryon density in this case. At sufficiently  the uncertainty in the exact value of neutrino trapping temperature T tr which is supposed to lie in that range. Figure 5 plots the dependence of the bulk viscosity on the temperature. As already discussed above, the bulk viscosity of neutrino-transparent matter attains its maximum at the temperature where γ(T) = ω, whereas the temperature dependence of the bulk viscosity in the neutrino-trapped matter is mainly decreasing (up to point where the matter becomes scale-invariant) because the relaxation rate is already too fast: the resonant maximum would be at lower temperatures. Since the relaxation rate γ rises as T 2 [Equation (35)] we expect [from Equation (28)] that ζ ∝ T −2 in this regime.
In Fig. 6, panel (a) we combine the results obtained for neutrino-transparent and neutrino-trapped matter by interpolating the results between these two regimes in the temperature range 5 ≤ T ≤ 10 MeV which is regarded as the transition region. Close to the transition temperature, the bulk viscosity is much larger in the neutrino transparent regime because of much lower beta relaxation rate and larger susceptibility C 2 /A as well. As a result, the resonant peak of the bulk viscosity occurs always at or below the neutrino-trapping temperature. Hence we can already anticipate that the bulk viscosity is going to play an important role in the dynamics of neutron star mergers in the regime of neutrino-transparent rather than neutrino-trapped matter.
For comparison, we show also the bulk viscosity of nuclear matter for an alternative model NL3 [37] in panel (b) of Fig. 6. This model has density-independent meson-nucleon couplings but contains non-linear terms in the σ-meson field. The results obtained within the two models differ mainly in the low-temperature regime, where the model NL3 features much higher viscosities than the model DDME2 above the saturation density. The reason for this is that NL3 model has a threshold of direct Urca opening at around n B 2.5n 0 whereas the model DDME2 does not have a threshold up to densities n B = 5n 0 . Because of the threshold, the model NL3 has much faster relaxation rates at n B 2n 0 than the model DDME2, and, as a consequence, the maximum of bulk viscosity is shifted to lower values of the temperature as compared to the DDME2 model.

Estimating damping timescales
In this subsection we examine the characteristic timescales of damping of density oscillations by the bulk viscosity in neutron star mergers. The energy density stored in baryonic oscillations with amplitude δn B is given by where is the isothermal incompressibility of nuclear matter. The energy dissipation rate by the bulk viscosity per unit volume is The characteristic timescale required for damping of oscillations τ ζ = /(d /dt) is then given by The incompressibility of nuclear matter is shown in Fig. 7 for the two parametrizations discussed above. It is an increasing function of the density and at low densities decreases from its value at zero temperature as the temperature is increased.  Figures 8a and 8b show, for two EoSs, the damping timescales of oscillations with frequency f = 1 kHz. We use the interpolation of bulk viscosity between neutrino-transparent and neutrino-trapped regimes shown in Fig. 6. As nuclear incompressibility depends weakly on the temperature, the temperature dependence of the damping timescale closely follows that of the bulk viscosity. As a result, the damping timescale attains a minimum at the temperature where the bulk viscosity has a maximum for a fixed value of the density. The minimal value of the damping timescale is In the limits of slow and fast equilibration the damping timescale is given by Thus, in contrast to the bulk viscosity, the damping timescale becomes frequency-independent in the low-temperature (slow equilibration) regime, and decreases with the frequency at high temperature (fast equilibration) regime. The density dependence of τ ζ reflects the density dependence of the ratio of nuclear incompressibility and the bulk viscosity. The density dependence of these two quantities almost compensates each other in the neutrino-transparent regime of the DDME2 model. In the case of NL3 model the density dependence of the bulk viscosity dominates and the damping timescale mainly decreases with density. The reason for this is the increase of the reaction rates with density as a result of fast opening of phase space for direct Urca reactions in the NL3 model.
The gray shaded areas in Fig. 8 show the temperature regions where the damping timescale is smaller than the characteristic timescales for the early ( 10 ms, dark shaded areas) and long-term ( 1 s, light shaded areas) post-merger evolution, respectively. It is seen that in massive neutron star mergers, bulk viscosity strongly damps density oscillations at densities n B n 0 in the temperature range 3 T 6 MeV for model DDME2 and 1 T 6 MeV for model NL3. On the timescale of long-term evolution the damping is efficient also at lower densities, and the whole range of temperatures where the bulk viscosity is relevant is 2 T 7 MeV for model DDME2 and 1 ≤ T ≤ 7 MeV for model NL3.
Comparing the two panels of Fig. 8 we observe that at low temperatures T ≤ 4 MeV and at densities n B ≥ 2n 0 the damping timescales are much shorter for the EoS model NL3 which has a direct Urca threshold. Thus, the dynamics and observational signatures of post-merger objects potentially contain information on whether the direct Urca processes are operative in the high-density domain of the neutron stars.
The temperature where the damping timescale reaches its minimum for f = 1 kHz is located around T 5 MeV at densities below the direct Urca threshold and around T 3 MeV above the threshold. These results agree well with the results of Ref. [6] obtained within the Fermi-surface approximation. The exact computation, however, obtained using exact beta-equilibrium condition for the neutrino-transparent matter at finite temperatures [29] suggests that the minimum of τ ζ is shifted to lower temperatures at densities which are below the direct Urca threshold, and, as a result, the minimum always appears around T 3 MeV [6] (note that the authors of Ref. [6] included also the modified Urca processes in their calculations, which, however, do not change the location of the maximum of bulk viscosity as their contribution is subdominant above T = 2 MeV).
The density-dependence of the damping timescale found here differs from that of Ref. [6] where τ ζ was found to reach its minimum at low densities n B n 0 . Apart from this, we find much lower values for the damping timescale at the minimum. This discrepancy arises because of the non-relativistic approximation for nucleon susceptibilities used in this work. This approximation works well at low densities, but strongly overestimates the susceptibility C 2 /A at higher densities n B ≥ 2n 0 .
The triangles in Fig. 2 show the values of the susceptibility C 2 /A obtained in Ref. [6] for the models DD2 and IUFSU. We see that, although the relativistic corrections to the spectrum of nucleonic excitations are about 20% at n B = 2n 0 , they need to be included in the susceptibilities. This will require also a fully relativistic study of the beta-equilibration rates which is relegated to a future work.  Figure 11. Same as Fig. 10 but for frequency f = 10 kHz.
In Fig. 9 we show the damping timescale for 10 kHz density oscillations. The minimum value of τ ζ , in this case, is smaller than in the case of f = 1 kHz by factors between 2 and 10, and the values of τ ζ in the neutrino-trapped regime are smaller by two orders of magnitude. However, the damping timescales of neutrino-trapped matter always remain larger than a second, since the bulk viscosity is not high enough to affect the evolution of mergers in this regime. Figures 10 and 11 show the dependence of the bulk viscous damping timescale on the density and temperature colormap for the oscillation frequency fixed at f = 1 kHz and f = 10 kHz, respectively. The white solid and blue dashed lines show where the damping timescale becomes equal to the timescales of 10 ms and 1 s, respectively. In the areas shaded in dark red the bulk viscous damping timescale is τ ζ ≤ 10 ms, therefore, the damping of density oscillations by the bulk viscosity is very efficient in those regimes. In the regions shaded in blue the role of the bulk viscosity in damping of oscillations is negligible, as the damping timescale τ ζ ≥ 1 s there.
For completeness we comment also on how our results will change if larger lepton fractions are considered. The case Y L = 0.4 was studied in our previous work Ref. [7], where the bulk viscosity was shown to be reduced by factors from 1 to 3 as compared to the Y L = 0.1 case. The pressure, and, therefore, also the nuclear incompressibility is only slightly sensitive to the lepton fraction. As a consequence, the damping timescales in the Y L = 0.4 case will be larger than in the Y L = 0.1 case by factors of a few, but the overall quantitative picture will remain the same.
In closing, we stress again that at densities n B ≥ 2n 0 the relativistic corrections to the spectrum of nucleonic excitations become important for the computation of the bulk viscosity, and our results at high densities need to be improved accordingly. Also the appearance of hyperons and other heavy baryons needs to be taken into account. Finally, we note that in the case of hybrid stars with quark cores, the bulk viscosity of quark matter can be important for damping of density oscillations (for computations in the case of cold compact stars see [30,[38][39][40][41][42][43][44][45][46]).

Summary
We have reviewed the computation, ingredients, and approximations involved in computations of bulk viscosity of nuclear matter at finite temperatures relevant to binary neutron star mergers. The bulk viscosity arises from the direct Urca β-equilibration reactions. A novel ingredient relative to the studies of cold neutron stars is the trapped neutrino component coexisting with the nuclear matter at temperatures T 5 MeV. The concrete computations were carried out with the relativistic density functional approach to the EoS of nuclear matter with two different parametrizations.
At a given value of oscillation frequency ω ≡ 2π f the bulk viscosity shows the standard resonant form (24), with a maximum where the beta relaxation rate γ matches ω. This resonant maximum is achieved in the temperature range where neutrinos escape from the merger region, since the relaxation rate at temperatures of a few MeV is sufficiently low to match the density oscillation frequency. The reason for lower relaxation rates as compared to the neutrino-trapped case is the suppression of the direct Urca processes at the relevant temperatures and densities.
When the temperature rises to the threshold for neutrino trapping (T ∼ 5 MeV) the bulk viscosity experiences a sharp fall by several orders of magnitude as the relaxation rate rises and the material enters the fast beta-equilibration regime with γ ω. In this regime the bulk viscosity is independent of the frequency and decreases with the temperature approximately as ζ ∝ T −2 . At temperatures of about 30 MeV a new feature appears: the bulk viscosity drops to zero at the temperature where the beta-disequilibrium-baryon-density susceptibility C vanishes, and then rises again at higher temperatures. The susceptibility vanishes because the particle fractions become independent of the density and the material becomes scale-invariant.
The main new result of this work concerns the timescales of damping of density oscillations in neutron star mergers by the bulk viscous dissipation. As an input we used the results for the bulk viscosity in Ref. [7]. We find that the damping timescale has a minimum as a function of temperature, which is located at temperatures in the range 3÷6 MeV for various densities. Assuming oscillation frequency of 1 kHz we find that the damping timescale at its minimum is of the order of ms, i.e., much shorter in the entire density range considered than the characteristic timescales of initial (∼ 10ms) and long-term (∼ 1s) post-merger evolution. We further find that the timescales of damping of density oscillations are shorter at the higher densities. If the temperature is above the neutrino trapping temperature, the damping timescales are much longer as the bulk viscosity is strongly suppressed. Finally, we note that bulk viscous dissipation could be of interest in the context of hydrodynamics simulations of supernovas, where electron capture rates on protons and nuclei could be out of equilibrium (for recent numerical simulations, see [47][48][49][50]).