Damping of density oscillations in neutrino-transparent nuclear matter

TL;DR

This paper calculates the damping times of density oscillations in nuclear matter relevant to neutron star mergers, considering temperature and density effects on bulk viscosity due to weak interactions.

Contribution

It provides new estimates of bulk viscous damping times at various densities and temperatures, extending previous low-temperature models to more realistic merger conditions.

Findings

Damping times range from 5 to 20 ms for 1 kHz oscillations.

Bulk viscosity can significantly influence neutron star merger dynamics.

Higher frequency oscillations damp faster, affecting merger evolution.

Abstract

We calculate the bulk-viscous dissipation time for adiabatic density oscillations in nuclear matter at densities of 1-7 times nuclear saturation density and at temperatures ranging from 1 MeV, where corrections to previous low-temperature calculations become important, up to 10 MeV, where the assumption of neutrino transparency is no longer valid. Under these conditions, which are expected to occur in neutron star mergers, damping of density oscillations arises from beta equilibration via weak interactions. We find that for 1 kHz oscillations the shortest dissipation times are in the 5 to 20 ms range, depending on the equation of state, which means that bulk viscous damping could affect the dynamics of a neutron star merger. For higher frequencies the dissipation time can be even shorter.

Figures

Equations

$$\begin{array}[]{rcl@{\quad}l}n&\rightarrow&p+e^{-}+\bar{\nu}_{e}&\mbox{direct Urca}\\ n+N&\rightarrow&N+p+e^{-}+\bar{\nu}_{e}&\mbox{modified Urca}\end{array}$$

$$\begin{array}[]{rcl@{\quad}l}e^{-}+p&\rightarrow&n+\nu_{e}&\mbox{direct Urca}\\ N+e^{-}+p&\rightarrow&N+n+\nu_{e}&\mbox{modified Urca}\end{array}$$

$${\Gamma }_{dU,nd}$$

$$\times {\left(2\pi \right)}^4{\delta }^4(p_n-p_p-p_e-p_{\nu })f_n\left(1-f_p\right)\left(1-f_e\right)$$

$${\Gamma }_{dU,ec}$$

$$\times {\left(2\pi \right)}^4{\delta }^4(p_n-p_p-p_e+p_{\nu })\left(1-f_n\right)f_p{f}_e,$$

$$\sum _{\text{spins}}|{\mathcal{M}}_{\text{dU}}{|}^2$$

$$\times \left(1+3g_A^2+\left(1-g_A^2\right)\frac{{\mathbf{p}}_e\cdot {\mathbf{p}}_{\nu }}{E_e{E}_{\nu }}\right),$$

$${\mu }_n={\mu }_p+{\mu }_e\text{(Fermi Surface approx)}.$$

$$\sum _{\text{spins}}|{\mathcal{M}}_{\text{dU}}{|}^2=32G^2(1+3g_A^2)E_n^{\ast }{E}_p^{\ast }{E}_e{E}_{\nu }$$

$${\Gamma }_{dU,nd}$$

$${\vartheta }_{dU}$$

$$A_{dU}$$

$${\Gamma }_{\text{mU,n}}=A_{mU}{G}^2{f}_{\pi NN}^4{g}_A^2\frac{m_n^3{m}_p}{m_{\pi }^4}\frac{p_{Fn}^4{p}_{Fp}}{{\left(p_{Fn}^2+m_{\pi }^2\right)}^2}{\vartheta }_n{T}^7,$$

$${\vartheta }_n$$

$${\Gamma }_{\text{mU,p}}$$

$$\displaystyle\times\dfrac{p_{Fn}(p_{Fn}\!-\!p_{Fp})^{4}}{\bigl{(}(p_{Fn}\!-\!p_{Fp})^{2}+m_{\pi}^{2}\bigr{)}^{2}}\,\vartheta_{p}\,T^{7}$$

$${\vartheta }_p$$

$$n_B(t)=n_B+\Delta n\sin \omega t,$$

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

$$\zeta =\frac{C^2}{B}\frac{\gamma }{{\omega }^2+{\gamma }^2},$$

$$B$$

$$C$$

$$\lambda=\frac{\partial(\Gamma_{n\rightarrow p}-\Gamma_{p\rightarrow n})}{\partial\mu_{\Delta}}\bigg{\rvert}_{\mu_{\Delta}=0}.$$

$$\lambda ={\lambda }_{dU}+{\lambda }_{mU,n}+{\lambda }_{mU,p},$$

$${\lambda }_{dU}$$

$${\lambda }_{mU,n}$$

$$\times \frac{p_{Fn}^4{p}_{Fp}}{(p_{Fn}^2+m_{\pi }^2{)}^2}{\vartheta }_n{T}^6,$$

$${\lambda }_{mU,p}$$

$$\times \frac{p_{Fn}(p_{Fn}-p_{Fp}{)}^4}{((p_{Fn}-p_{Fp}{)}^2+m_{\pi }^2{)}^2}{\vartheta }_p{T}^6.$$

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

$${\mu }_{\Delta }={\mu }_n-{\mu }_p-{\mu }_e-{\mu }_{\delta }.$$

$$E_i=U_i+m_i+\frac{p_i^2}{2m_i}\text{mbox}for\text{leavevmode}i=(p,n)$$

$${\zeta }_{\text{max}}=\frac{C^2}{2B\omega }.$$

$$\varepsilon=\frac{1}{2}(\Delta n)^{2}\frac{\partial^{2}\varepsilon}{\partial n_{B}^{2}}\bigg{|}_{x_{p},s/n_{B}}=\frac{\kappa_{S}^{-1}}{2}\left(\frac{\Delta n}{n_{B}}\right)^{2}\ ,$$

$$\kappa_{S}^{-1}=n_{B}\frac{\partial P}{\partial n_{B}}\bigg{|}_{x_{p},s/n_{B}}\ .$$

$$\frac{\mathrm{d\epsilon }}{\mathrm{dt}}=\frac{{\omega }^2}{2}{\left(\frac{\Delta n}{n_B}\right)}^2\zeta .$$

Full text

I Introduction

The observation of gravitational waves from neutron star mergers TheLIGOScientific:2017qsa has drawn attention to the importance of understanding the properties of nuclear matter at the densities and temperatures attained in mergers. Recent studies Abbott:2018exr ; Landry:2018prl of the gravitational wave signals estimate that the two neutron stars involved in GW170817 reached central densities of at least twice nuclear saturation density ($n_{\text{sat}}\equiv 0.16\text{ fm}^{-3}$). Numerical simulations Perego:2019adq ; Hanauske:2019qgs ; Hanauske:2017oxo ; Baiotti:2016qnr ; Kastaun:2016elu ; Bernuzzi:2015opx ; Foucart:2015gaa ; kiuchi:2012mk ; sekiguchi:2011zd find that during the first 10 to 20 ms after the neutron stars make contact the material at densities of up to several times $n_{\text{sat}}$ can reach temperatures of many tens of MeV, perhaps up to 80 Hanauske:2019qgs ; Hanauske:2017oxo or even 100 MeV Perego:2019adq , and co-moving fluid elements are subjected to strong density oscillations at a typical frequency of 1 kHz Alford:2017rxf . This raises the possibility of bulk viscosity in nuclear matter playing an important role if it is strong enough to damp those oscillations on a short enough timescale to affect the dynamics of the merger111We use the term “merger” to refer to the late stages of the inspiral as well as the process that begins when the stars touch..

An initial estimate for $npe\mu$ matter Alford:2017rxf suggested that the bulk viscous dissipation time could be extremely short: in the range of a few milliseconds. In this work we focus on $npe$ matter in the neutrino-transparent regime ($T\lesssim 5$ to 10 MeV Alford:2018lhf ; Roberts:2016mwj ; Haensel:1987zz ; 1979ApJ…230..859S ; Sawyer:1975js ), where neutrinos escape from the merger region, and make a detailed study. We calculate the dependence of the bulk viscous dissipation time on density and temperature, and we explore its sensitivity to the equation of state. We find that the lowest dissipation times occur at temperatures of about 3 MeV, which lies in the neutrino-transparent regime, since neutrinos with energies of 3 MeV have mean free paths on the order of a few kilometersAlford:2018lhf ; Sawyer:1975js ; 1979ApJ…230..859S . For very low densities ($n_{B}\approx 0.5n_{\text{sat}}$), the dissipation time of 1 kHz oscillations is as low as 5 ms. At densities above $n_{\text{sat}}$, the dissipation times can be as low as 20 ms. These times are similar for both equations of state that we study.

Bulk viscosity for oscillations in the kHz range arises from beta equilibration via weak interactions. Neutrino transparency means there is no Fermi sea of neutrinos, so in the relevant Urca processes neutrinos only occur in final states Yakovlev:2000jp ; Shapiro:1983du . This deviation from detailed balance leads to corrections to the standard Fermi Surface (FS) approximation when the temperature rises above about 1 MeV Alford:2018lhf . Those corrections are included in our calculations.

We work in natural units, where $\hbar=c=k_{B}=1$. All data presented in our figures can be found in the Supplemental Material supp .

II Nuclear matter and the Urca process

Beta equilibration establishes the stable proton fraction via Urca processes. In this work we will focus on processes involving electrons, leaving muon contributions for future investigation. If there is a deficit of protons, protons are created via the neutron decay processes

[TABLE]

Here, $N$ denotes a spectator neutron or proton. If there is a deficit of neutrons, neutrons are created via the electron capture processes

[TABLE]

The direct Urca process is in general faster than the modified Urca process, since it involves two fewer particles.

The rate of the direct Urca neutron decay and electron capture processes are given by the twelve dimensional phase space integrals Yakovlev:2000jp

[TABLE]

where $f_{i}=1/\{1+\exp[(E_{i}-\mu_{i})/T]\}$ are the Fermi-Dirac distributions for species $i=n,p,\text{ or }e$ with chemical potential $\mu_{i}$ and the matrix element is Yakovlev:2000jp

[TABLE]

where $G^{2}=G_{F}^{2}\cos^{2}{\theta_{c}}=1.29\times 10^{-22}{\rm MeV}^{-4}$ and $g_{A}=1.26$. We will describe the nucleon energy dispersion relations and discuss the significance of using $E^{*}$ instead of $E$ for the nucleons in Sec. V.1. We note here that the matrix elements in Alford:2018lhf denoted $\langle|\mathcal{M}|^{2}\rangle$ are in fact $\sum_{\text{spins}}|\mathcal{M}|^{2}/(2E_{1}\times 2E_{2}\times...\times 2E_{n})$ for a process with $n$ particles in the initial and final states combined.

The rate of modified Urca neutron decay and electron capture with a nucleon spectator is given by the eighteen dimensional phase space integrals given in Haensel:2001mw .

III Bulk viscosity in the Fermi Surface approximation

III.1 Urca rates

At low temperatures $T\lesssim 1\,{\rm MeV}$ Alford:2018lhf , where the Fermi surfaces are sharply defined, we can make the approximation that only particles near the Fermi surface can participate in Urca processes. We call this the “Fermi Surface (FS) approximation”. At these temperatures the beta equilibrium condition is Yakovlev:2000jp ; Shapiro:1983du

[TABLE]

This condition enforces the equality of Urca rates for proton production and proton capture, yielding a proton fraction that is constant in time Yuan:2005tj . In the Fermi Surface approximation the phase space integrals (Sec. II) can be simplified by fixing the momentum magnitudes to the corresponding Fermi momenta.

In the Fermi surface approximation, the direct Urca matrix element (Eq. II) simplifies to

[TABLE]

under the assumption of non-relativistic nucleons Yakovlev:2000jp . In the Fermi Surface approximation, in (low-temperature) beta-equilibrium (Eq. 6) the rate of direct Urca neutron decay and direct Urca electron capture are equal and are given by Alford:2018lhf ; Yakovlev:2000jp ; Haensel:2000vz

[TABLE]

We see that in this approximation, the direct Urca process has a density threshold, above which the process is kinematically allowed.

In the Fermi surface approximation, the neutron-spectator modified Urca rates (both neutron decay and electron capture, which are equal in beta equilibrium (Eq. 6)) are given by

[TABLE]

and the proton-spectator modified Urca rates are given by Alford:2018lhf ; Yakovlev:2000jp

[TABLE]

[TABLE]

III.2 Bulk viscosity

Consider a fluid element of nuclear matter subjected to a small-amplitude, periodic baryon density oscillation

[TABLE]

where $\Delta n\ll n_{B}$. Since the equilibrium proton fraction varies with density, changing the density will temporarily push the nuclear matter out of beta equilibrium by an amount Yakovlev:2000jp ; 1992A&A…262..131H ; fi09000a

[TABLE]

We will consider only “subthermal” density oscillations where $\mu_{\Delta}\ll T$. In response to the density change, the nuclear matter will try to reestablish beta equilibrium via the Urca process which has a characteristic rate $\gamma(n_{B},T)$ Yakovlev:2000jp . The bulk viscosity of neutrino-transparent nuclear matter is given by Alford:2010gw ; Schmitt:2017efp ; Haensel:1992zz

[TABLE]

where $B$ and $C$ are susceptibilities of the nuclear matter which depend on the equation of state but not on the weak-interaction equilibration rate,

[TABLE]

In this paper we assume that there is negligible heat flow between adjacent fluid elements during the merger. This is valid as long as the thermal equilibration time in the absence of neutrinos is much longer than about 10 ms. From Eq. (1) of Ref. Alford:2017rxf this will be true as long as density oscillations (and the resultant thermal gradients) have wavelengths longer than about a meter. This criterion is obeyed in current simulations, whose spatial resolution is tens of meters at best. Since heat does not flow we use adiabatic susceptibilities, evaluating the derivatives Eq. 30 at constant entropy per baryon $S/N_{B}$, or equivalently, at constant entropy density per baryon density $s/n_{B}$. See App. A for a comparison of adiabatic and isothermal thermodynamic quantities.

The equilibration rate $\gamma=B\lambda$, where

[TABLE]

In the Fermi Surface approximation, we can analytically compute $\lambda$ from Eqs. 8,12, and 16, Yakovlev:2000jp ; Haensel:1992zz ; Reisenegger:1994be ; Friman:1978zq ; Sawyer:1989dp

[TABLE]

where

[TABLE]

From Eq. 29 it follows that, for an oscillation of fixed frequency $\omega$, the bulk viscosity has a resonant maximum when the equilibration rate (which varies as a function of density and temperature) coincides with the oscillation frequency, i.e., when $\gamma(n_{B},T)=\omega$. A major goal of this paper will be to map the regions in density and temperature where this maximum is achieved.

IV Intermediate temperatures

IV.1 Urca process at intermediate temperatures

Below the direct Urca threshold, neutrons, protons, and electrons on their Fermi surfaces cannot participate in the direct Urca process while conserving energy and momentum. Direct Urca still occurs, but it involves particles away from their Fermi surfaces, and so the rate is Boltzmann suppressed. As temperatures rise above 1 MeV, up to the neutrino trapping temperature, the Boltzmann suppression lessens, and the direct Urca rate becomes comparable to and then, eventually, greater than the modified Urca rate, thus broadening the direct Urca threshold. In contrast, the presence of a spectator nucleon means that neutrons, protons, and electrons close to their Fermi surfaces can participate in a modified Urca process, and so the Fermi surface approximation is still appropriate for modified Urca, even at the moderately high temperatures that we consider Alford:2018lhf .

The beta equilibrium condition becomes

[TABLE]

where the additional chemical potential $\mu_{\delta}$ arises from the absence of detailed balance: neutron decay and electron capture processes are not exact inverses of each other. $\mu_{\delta}$ is a function of temperature and baryon density, and its value is determined by the requirement that the rates for neutron decay and electron capture rates must balance Alford:2018lhf . These rates include direct Urca contributions, calculated by integration of the full phase space, and modified Urca contributions for which we can use the Fermi Surface approximation as described in the appendix of Ref. Alford:2018lhf .

To illustrate the breakdown of the Fermi Surface approximation at temperatures relevant for neutron star mergers, we show in Fig. 1 the comparison between rates calculated in the Fermi Surface approximation (using the low-temperature beta equilibrium condition Eq. 6) and the full phase space integration (with the general beta equilibrium condition Eq. 36). We show the neutron decay rate, which is equal to the electron capture rate in each case, calculated for neutrino-transparent nuclear matter described by the IUF equation of state, described below. We see that at $T=4\,{\rm MeV}$ the Fermi surface approximation makes the direct Urca threshold seem unphysically sharp, underestimating the below-threshold rates by an order of magnitude and overestimating the above-threshold rates by a factor of 3 to 5.

IV.2 Bulk viscosity at intermediate temperatures

To calculate bulk viscosity at temperatures where the Fermi Surface approximation is not valid, we use the appropriate characterization of the deviation from equilibrium,

[TABLE]

When $\mu_{\Delta}=0$ the system is in true beta equilibrium (36). In this regime, the equilibration rate $\lambda$ (Eq. 32) no longer has a simple analytic form, as the direct Urca rates (3) and (4) can only be simplified to three-dimensional numerical integrals Alford:2018lhf . We obtain $\lambda$ by calculating the difference $\Gamma_{n\rightarrow p}-\Gamma_{p\rightarrow n}$ at proton fractions around the beta equilibrium value where the rates are equal. The slope of the difference as a function of $\mu_{\Delta}$ at $\mu_{\Delta}=0$ is $\lambda$.

V Results

V.1 Models of nuclear matter

To gauge the sensitivity of our results to the equation of state of nuclear matter, we will use two representative equations of state, one stiffer (DD2Hempel:2009mc ; Typel:2009sy ) and one softer (IUF Hempel:2009mc ; Fattoyev:2010mx ; RocaMaza:2008ja ). Both are tabulated at the CompOSE website (https://compose.obspm.fr/eos/18/ and https://compose.obspm.fr/eos/22/). Both are relativistic mean field theories, where nucleons interact strongly by exchanging sigma, omega, and rho mesons. The couplings between the mesons and nuclei are chosen to reproduce nuclear observables like the nuclear binding energy, saturation density, symmetry energy, incompressibility, among others glendenning2000compact ; Dutra:2014qga ; Oertel:2016bki . Aside from the standard nucleon-meson linear couplings, the IUF model has self-interactions in the sigma and omega fields, plus interactions among the mesons (see Eq. 7 in Dutra:2014qga ). DD2 has the standard linear nucleon-meson couplings, but instead of nonlinear self-interactions or meson-meson interactions, it promotes the nucleon-meson couplings to density-dependent functions Dutra:2014qga . DD2 has no direct Urca threshold (ie, direct Urca is kinematically forbidden at all densities in the Fermi surface approximation) (see Fig. 2 in Alford:2017rxf ), while IUF has a direct Urca threshold near $4n_{\text{sat}}$. The maximum mass neutron star for DD2 is $2.42\,M_{\odot}$ and for IUF it is $1.96\,M_{\odot}$.

In our calculations we assume the neutrinos and electrons are ultrarelativistic free particles, and nucleonic excitations have the dispersion relation Roberts:2016mwj ; Roberts:2012um

[TABLE]

where $U_{i}$, the nuclear mean field, is chosen as a function of density and temperature so that the Fermi energy $E_{F,i}\equiv E_{i}(p_{F,i})$ matches the chemical potential $\mu_{i}$ for the given EoS. For the mass $m_{i}$, we use the rest mass in vacuum. The microscopic origin of the nuclear mean field can be understood through the framework of relativistic mean field theories, where the mean field $U$ is a function of the vacuum expectation values of the $\omega$ and $\rho$ mesons which are the strong force carriers between nucleons Fu:2008zzg .

As alluded to in Sec. II, the neutron and proton energies can be written as $E_{i}=U_{i}+E_{i}^{*}$, where $E_{i}^{*}=m_{i}+p^{2}/(2m_{i})$. In the rate calculations (3) and (4), $E^{*}$ should be used for the energies in the matrix element and in the energy factors in the denominator, while $E$ should be used in the energy delta function and the Fermi Dirac factors Leinson:2002bw ; Roberts:2016mwj ; Fu:2008zzg . However, in the approximation we used for the direct Urca matrix element (7), the $E^{*}$ factors cancel out.

V.2 Bulk viscosity

In Fig. 2, we show the bulk viscosity of nuclear matter with the DD2 and IUF EoS, when subjected to a 1 kHz density oscillation, which is a typical frequency for neutron star mergers Alford:2017rxf . The dashed lines are the bulk viscosity with Urca rates calculated in the Fermi Surface approximation while the solid lines use the exact Urca rates.

In Fig. 2(a), corresponding to the DD2 EoS, the exact bulk viscosity peaks at a temperature that is 1-2 MeV lower than would be predicted by the Fermi Surface approximation. This is because DD2 never allows direct Urca (the threshold is at infinite density), and we know (see Fig. 1) that the Fermi Surface approximation underestimates the below-threshold Urca rate. This means that in the Fermi Surface approximation the temperature must be pushed up to a higher value in order for the equilibration rate to match the oscillation frequency, which is where the resonant peak occurs (Sec. III.2).

In Fig. 2(b), corresponding to the IUF EoS, which has a direct Urca threshold near $4n_{0}$. Here we see two distinct behaviors. For densities $n_{0}$ and $3n_{0}$, which are below threshold, the behavior is similar to that seen for DD2: the Fermi Surface approximation only includes modified Urca processes, but the exact calculation includes below-threshold direct Urca processes which increase the total rate, moving the resonant peak to lower temperatures. Above the threshold density, the Fermi Surface approximation for direct Urca overestimates the total Urca rate, since the exact phase space integration leads to only a gradual opening of the phase space around the direct Urca threshold, hence the resonant peak moves to higher temperatures than predicted by the Fermi Surface approximation.

As can be seen from Eq. 29, the maximum value of bulk viscosity at a frequency $\omega$ is

[TABLE]

In Fig. 3 we plot $C^{2}/B=2\omega\zeta_{\text{max}}$ for a representative range of densities and temperatures for which nuclear matter is likely neutrino-transparent. We see that for a given frequency, the maximum value of bulk viscosity varies by 1-2 orders of magnitude, and depends more strongly on density than on temperature. Most notably, we can see that $C^{2}/B$ rises rapidly at low densities, then levels off at $n\sim 2n_{\rm sat}$ to a value about an order of magnitude larger then its value at $n=n_{\rm sat}$. (This could already be seen in Fig. 2).

In Fig. 4, we plot the bulk viscosity as a function of density and temperature (the curves in Fig. 2 are cross-sections through Fig. 4). For a fixed density, as the temperature rises, the beta equilibration rate $\gamma$ rises rapidly because of the increase in available phase space. At temperatures of a few MeV, the reequilibration rate closely matches the oscillation frequency of 1 kHz, then bulk viscosity reaches a maximum. At higher temperatures, the reequilbration is too fast and the bulk viscosity drops.

We see that for the DD2 EoS, the bulk viscosity peak is at a temperature of about 3 MeV for all densities, which is a lower temperature than predicted by the Fermi surface approximation. For IUF, the FS approximation would suggest two different peaks in bulk viscosity: one below the direct Urca threshold corresponding to the near-equality of the modified Urca rate and the density oscillation frequency, and one above the threshold, corresponding to the near-equality of the direct Urca rate and the density oscillation frequency. However, the gradual opening of the direct Urca threshold coming from the exact direct Urca calculation melds these two peaks into one broad peak. At low density, the peak is at 3-4 MeV, but as density increases it moves down to 2 MeV.

V.3 Energy dissipation time

The most direct indicator of the importance of bulk viscous damping is the dissipation time $\tau_{\rm diss}$ for density oscillations. Since the merging stars settle down into a massive remnant in tens of milliseconds, bulk viscous damping will be important if $\tau_{\rm diss}$ is tens of milliseconds or less. To calculate the dissipation time, we need the energy of an oscillation and the rate at which that energy is dissipated by bulk viscosity. The energy density of an adiabatic baryon density oscillation $n_{B}(t)=n_{B}+(\Delta n)\sin{(\omega t)}$ is Alford:2017rxf

[TABLE]

where $\kappa_{S}$ is the adiabatic compressibilityschroeder1999introduction ; compose_manual

[TABLE]

See Sec. III.2 for a discussion of the assumption of adiabaticity.

We note that the adiabatic compressibility depends on the EoS. However, it is a common feature of all nucleonic EoSs that nuclear matter becomes more incompressible at high densities, so the inverse compressibility $1/\kappa_{S}$ rises with density, as shown in Fig. 5 for a range of EoSs including those used in this work. This means that at higher density, oscillations in the density store more energy.

To facilitate comparison with previous work (for example, Alford:2017rxf ), we mention that the “stiffness” of nuclear matter is often described via the nuclear incompressibility $K$ Shapiro:1983du ; glendenning2000compact ; Schmitt:2010pn . $K$ is conventionally defined at saturation density, zero temperature, and for symmetric nuclear matter, and is approximately 250 MeV Oertel:2016bki . Some works have extended the definition of the nuclear incompressiblity to densities above nuclear saturation Dexheimer:2007mt . At zero temperature, $n_{\text{sat}}$, and for symmetric nuclear matter, the adiabatic $\kappa_{S}$ can be related to the nuclear incompressibility $K$ by $K=9/(\kappa_{S}n_{0})$ compose_manual ; Schmitt:2010pn .

The rate of energy density dissipation is given by Alford:2010gw ; Sawyer:1980wp

[TABLE]

Using Eq. 40, the energy dissipation time is

[TABLE]

Note that one can also define 1990ApJ…363..603C a decay time for the amplitude, which would be longer by a factor of two since the energy of an oscillation goes as the square of the amplitude.

In Fig. 6, we plot the dissipation time of a 1 kHz density oscillation as a function of density and temperature for two different EoSs, using the exact Urca rates. We first discuss the physical content and implications of the exact results (Fig. 6), then compare them to the Fermi Surface approximation, shown in Fig. 7.

Temperature dependence. The adiabatic compressibility is relatively independent of temperature, so the bulk viscosity dominates the temperature dependence of the dissipation time. As discussed in Sec. V.2, for a given density, the bulk viscosity increases, reaches a resonant maximum when the beta reequilibration rate $\gamma$ matches the oscillation frequency $\omega$, and then decreases as temperature increases. This leads to minimum dissipation time at approximately the temperature at which the bulk viscosity reaches its maximum, for a given density.

Density dependence. The adiabatic inverse compressibility strongly increases as a function of density, as seen in Fig. 5. While the bulk viscosity was weakly dependent on density, the dissipation time at high density is strongly increased due to the several order-of-magnitude rise of the adiabatic inverse compressibility. Physically, oscillations in high density nuclear matter have a lot of energy due to the high incompressiblity of dense nuclear matter (see Eq. 40). Thus, it takes correspondingly longer time for those high energy oscillations to damp. As a result of the behavior of the compressibility of nuclear matter, the minimum of dissipation time is likely to be located at a low density.

It is worth noting that the bulk viscosity varies non-monotonically with density. It rises as density increases from $0.5n_{\text{sat}}$, reaches a peak at several times $n_{\text{sat}}$, and then falls off at high density. This can be seen by noting that the maximum bulk viscosity is $\zeta_{\rm max}=(1/2\omega)C^{2}/B$ (Eq. 39), which is plotted in Fig. 3. It is clear that the particular features of the rise and fall in bulk viscosity as a function of density depend on the EoS. Throughout the range of densities that we consider, the bulk viscosity prefactor $C^{2}/B$ varies by 1-2 orders of magnitude. However, the inverse compressibility rises by three orders of magnitude over that density range, so it has a more substantial effect on the density dependence of the dissipation time.

For the DD2 EoS, as seen in Fig. 6(a), the minimum dissipation times lie around temperatures of 3 MeV for all densities, indicating that the reequilibration rate doesn’t change strongly with density, which is expected since only modified Urca and below-threshold direct Urca are acting. As a function of density, the dissipation times get longer as density increases. This behavior comes from the dramatic monotonic rise of the inverse compressibility as a function of density. The bulk viscosity prefactor $C^{2}/B$ rises by one order of magnitude from $0.5n_{\text{sat}}$ to 3 or 4 $n_{\text{sat}}$, and then slightly decreases at higher densities, but it doesn’t vary rapidly enough to compete with the rise of the inverse compressibility, and thus the dissipation time rises monotonically with density. DD2 has a minimum dissipation time of about 6 ms, which occurs only at low density ($0.5n_{\text{sat}}$) at temperatures of just under 3 MeV. Only fluid elements with densities under twice saturation density would dissipate energy on timescales relevant for mergers.

As seen in Fig. 6(b), the behavior of the dissipation time scale for the IUF EoS is more complicated. The lowest dissipation times do occur at temperatures of around 3 MeV, since the resonant peak of bulk viscosity is around that temperature. However, the nonmonotonic behavior of $C^{2}/B$ as a function of density is more dramatic for the IUF EoS than for DD2, so it competes with the rapidly rising inverse compressibility as density increases, leading to two minima in the dissipation time. The first is at low density, where the nuclear inverse compressibility is decreasing rapidly as the density decreases to the lowest value for which we trust our equation of state, $n=0.5n_{\text{sat}}$. There, energy dissipation can occur in as little as 5 ms. There is also a local minimum around $n=2n_{\text{sat}}$, where the bulk viscosity prefactor $C^{2}/B$ has a local maximum (see Fig. 3(b)) and dissipation times reach down to 19 ms. For the IUF EoS, dissipation occurs on merger timescales in fluid elements up to four times saturation density, in contrast to the behavior of DD2.

It is interesting to compare the Fermi Surface approximate results (Fig. 7) and the exact results (Fig. 6) for each EoS. For DD2, the use of the exact Urca rates just increases the total Urca rate and thus the bulk viscosity is maximized at a lower temperature than would be predicted by the Fermi Surface approximation. For IUF, the Fermi Surface approximate result would predict a sharp change in the behavior of the bulk viscosity at the direct Urca threshold, $n=4n_{\text{sat}}$ (for a generic example of this behavior, see Figure 1 in Haensel:2001mw ). However, at the temperatures of interest to us the exact Urca rates show a gradual increase with density and thus the bulk viscosity does not change suddenly at the threshold density.

V.4 Higher frequency oscillations

There is evidence from simulations Hotokezaka:2013iia ; Chaurasia:2018zhg ; Clark:2015zxa ; Rezzolla:2016nxn ; Dietrich:2016hky ; Maione:2017aux ; Bauswein:2018bma (see also the review Baiotti:2019sew ) that eccentric binary neutron star mergers excite oscillations at frequencies above 1 kHz. We plot the dissipation times for 3 kHz and 5 kHz oscillations in Figs 8 and 9. We see that at these higher frequencies, bulk viscosity plays a bigger role, and density oscillations can be damped in as little as 1 ms, and for a broad range of temperatures and densities, oscillations can be damped in under 25 ms.

We note that, at a given density, a higher temperature is required to make the reequilibration rate $\gamma$ match a density oscillation which has a frequency above 1 kHz, and thus the region of maximum bulk viscosity is moved to higher temperatures. For example, a 5 kHz density oscillation has maximum bulk viscosity (and thus minimum damping time) at about $T=4\text{ MeV}$ (see Fig. 9), while a 1 kHz density oscillation has maximum bulk viscosity at around $T=3\text{ MeV}$ (see Fig. 6).

VI Conclusions

We have calculated the bulk-viscous dissipation time in nuclear matter at temperatures and densities relevant to neutron star mergers. We assumed the material was transparent to neutrinos, which should be valid for temperatures up to about 5 MeV, and we studied the damping of oscillations with frequencies in the 1 kHz range, which are seen in simulations of mergers. The main uncertainty in our result is the form of the nuclear matter equation of state at supranuclear densities, so we performed calculations for two different equations of state, one stiffer, DD2, and one softer, IUF, which differ in their treatment of the nucleon-meson interaction (see Sec. V.1). Our main results are displayed in Fig 6.

Bulk viscous damping will play a significant role at densities and temperatures where the dissipation time is comparable to or less than the typical timescale of the merger, which is in the range of tens of milliseconds. Both equations of state show a similar overall pattern: bulk viscosity damps oscillations on timescales comparable to a merger for nuclear matter at temperatures of 2-4 MeV and for densities between $0.5n_{\text{sat}}$ to $2n_{\text{sat}}$, with IUF also exhibiting fast damping for densities up to $4n_{\text{sat}}$. Both EoSs have minimum dissipation times of about 5 ms, occurring at $0.5n_{\text{sat}}$, while IUF has another local minimum of dissipation time, about 20 ms, occurring at $2n_{\text{sat}}$. The occurrence of dissipation times in the 10 ms range leads us to conclude that bulk viscous damping should be seriously considered for inclusion in future simulations. The strong dissipation that we see at low density may be relevant to the density oscillations that were found to be associated with mass ejection in the outer regions of the merger (See Fig. 9 of Ref. DePietri:2019khb ).

There are several directions in which this topic could be further developed. In parallel with this work, an analogous calculation of the bulk viscosity for nuclear matter with trapped neutrinos was conducted Alford:2019kdw . This is appropriate for nuclear matter at temperatures well above 5 MeV, and is relevant to mergers because temperatures up to 80 or even 100 MeV are predicted by simulations. It should be noted that the neutrino-transparent and neutrino-trapped regimes are the simplified extremes of a continuum, the whole of which is probably realized at different regions and stages of a merger. Between these extremes lies the regime where the spectrum of neutrinos includes a low energy population that escapes, a high energy tail that is trapped, and an intermediate energy range where the mean free path is comparable to the distance scale of the fluid flows, requiring explicit inclusion of neutrinos in the dynamics of the nuclear fluid Ardevol-Pulpillo:2018btx ; Perego:2014qda ; Galeazzi:2013mia ; Sekiguchi:2012uc ; Rosswog:2003rv ; 1999JCoAM.109..281M .

Another limitation of our calculation is the assumption of low-amplitude density oscillations. We calculated the “subthermal” bulk viscosity, but simulations show high amplitude density oscillations Alford:2017rxf for which the suprathermal bulk viscosity Alford:2010gw is relevant. This could extend the region of large bulk viscosity down to lower temperatures, since suprathermal effects allow high-amplitude oscillations to experience the maximum bulk viscosity $\zeta_{\rm max}$ (Eq. 39) at lower temperatures Alford:2010gw .

Our discovery of short bulk-viscous dissipation times at densities below nuclear saturation density, primarily due to the low inverse compressibility, underscores the need for a detailed understanding of the structure of nuclear matter below saturation density. The DD2 and IUF EoSs predict uniform nuclear matter down to densities of 0.25 to 0.4 $n_{\text{sat}}$ respectively, which is why we restricted our calculations to densities above $0.5n_{\text{sat}}$. However, a sequence of mixed “pasta” phases has been predicted at densities between $0.2$ and $0.7$ $n_{\text{sat}}$ Grill:2014aea ; Fattoyev:2017zhb ; Pais:2015lma ; Pais:2014hoa ; Oyamatsu:1993zz . It has been noted Yakovlev:2018jia ; Gusakov:2004mj that the appearance of free protons in certain pasta phases would open up the direct Urca process, albeit with such a reduced rate that it would take temperatures of tens of MeV—which is well above the pasta melting temperature of a few MeV Roggero:2017pag —to reach the resonant peak of bulk viscosity. Thus, it is important to know how and at what densities and temperatures nuclear matter transforms from a uniform phase to a mixed phase. Based on our findings above, we expect subthermal bulk viscosity to be large for these low densities, down to the density at which uniform nuclear matter transitions to a pasta or spherical nuclei phase.

As mentioned in section II, we did not consider Urca processes involving muons, and did not include muons in the EoSs. The presence of muon Urca processes would increase the equilibration rate $\gamma$ for densities at which muons are present. In addition, muon-electron conversion would give rise to a separate contribution to the bulk viscosity Alford:2010jf . The calculation of bulk viscous damping time in Ref. Alford:2017rxf uses EoSs that contain muons. Above the onset density for muons the nuclear matter susceptibilities are larger, which would lead to larger bulk viscosity and thus shorter dissipation times compared to the muonless EoSs considered in this work. We are therefore planning to perform a full study of bulk viscous dissipation in EoSs that include muons.

There is evidence that properly including in-medium effects in the nucleon propagator can lead to a large increase in the modified Urca rate just below the direct Urca threshold Shternin:2018dcn . We have not included this in our analysis, but it could potentially lead to a shift of the resonant peak of bulk viscosity to lower temperatures for a range of densities near the direct Urca threshold.

VII Acknowledgments

We thank Armen Sedrakian, Kai Schwenzer, Alex Haber, Reed Essick, Andreas Reisenegger, Alessandro Drago, and Jürgen Schaffner-Bielich for discussions. We also thank the Frankfurt Institute for Advanced Study where part of this work was performed. This research was partly supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Award No. #DE-FG02-05ER41375.

Appendix A Adiabatic and isothermal oscillations

Most previous works use the isothermal susceptibilities

[TABLE]

often only considering the zero temperature caseYakovlev:2018jia ; Alford:2010gw ; Haensel:2000vz ; Haensel:1992zz ; Sawyer:1989dp . As discussed in section III.2, because thermal equilibration is so slow in neutrino-transparent nuclear matter in merger conditions, we must use the adiabatic susceptibilities

[TABLE]

Note that at zero temperature, adiabatic and isothermal quantities become equivalentschroeder1999introduction ; swendsen2012introduction . Often, it is convenient to work with thermodynamic derivatives at constant temperature $T$, or baryon density $n_{B}$, or proton fraction $x_{p}$. In particular, these three variables are the degrees of freedom in the CompOSE database of EoSs compose_manual . Using a Jacobian coordinate transformation swendsen2012introduction , we can relate adiabatic derivatives (derivatives at constant entropy per baryon) to isothermal derivatives. The adiabatic susceptibility derivatives are related to the isothermal susceptibility derivatives by

[TABLE]

The isothermal compressibility is

[TABLE]

and the adiabatic compressibility is given by

[TABLE]

The adiabatic derivative can be obtained from the isothermal derivative by

[TABLE]

Above $n_{\text{sat}}$, the adiabatic and isothermal derivatives are within 25% of each other for all temperatures considered here. For nuclear matter that is below $n_{\text{sat}}$ with $T>5\,\text{MeV}$, there are noticeable differences between the isothermal and adiabatic suseptibility $C$ and the compressibility $\kappa$. The susceptibility $B$ is not sensitive to differences between adiabaticity and isothermality (the differences are below 10%).

Below $n_{\text{sat}}$, the adiabatic $C$ is greater than the isothermal $C$ by as much as a factor of 2.5 (DD2) or 5.5 (IUF). These large differences are at temperatures above 5 MeV. Thus, the adiabatic $C^{2}/B$ is larger than the isothermal version by factors of up to 6 (DD2) or 30 (IUF). However, these large differences occur at low densities ($\approx 0.5n_{\text{sat}}$) and high temperatures ($T\approx 10\,\text{MeV}$) where the bulk viscosity is small anyway, since the equilibration rate $\gamma$ is much faster than a 1 kHz density oscillation. In the regions where bulk viscosity is large, the difference between adiabatic and isothermal susceptibilities is at most a factor of 2 in the quantity $C^{2}/B$.

At the densities and temperatures where bulk viscosity is large, the isothermal compressibility is at most 20% larger than the adiabatic compressibility, which means that adiabatic density oscillations would lose energy slightly more slowly than isothermal density oscillations. At densities below $n_{\text{sat}}$ and temperatures above 5 MeV, the isothermal compressibility can be up to 40% (DD2) or 80% (IUF) larger than the adiabatic value, but the bulk viscosity is too small for fluid elements under these conditions for this to matter.