Properties of Hot Nuclear Matter

Abstract

A fully quantitative description of the equilibrium and dynamical properties of hot nuclear matter will be needed for the interpretation of the available and forthcoming astrophysical data, providing information on the post-merger phase of a neutron star coalescence. We discuss the results of a recently developed theoretical model, based on a phenomenological nuclear Hamiltonian including two- and three-nucleon potentials, to study the temperature dependence of average and single-particle properties of nuclear matter relevant to astrophysical applications. The potential of the proposed approach for describing dissipative processes leading to the appearance of bulk viscosity in neutron star matter is also outlined.

1. Introduction

The interpretation of the presently available and future astronomical data providing information on the post-merger phase of coalescing binary neutron stars will require an accurate description of the properties of dense nuclear matter at temperatures as high as 100 MeV [1,2,3,4,5]. Of great importance, in this context, will be the development of a consistent framework suitable for modelling both the equilibrium configurations—determining the equation of state (EOS) of neutron star matter—and dissipative processes, involving mechanisms that lead to the appearance of bulk viscosity [6] and neutrino emission [7].

The EOS of hot nuclear matter is often derived from dynamical models based on the independent-particle approximation, using Skyrme-type effective interactions [8] or the formalism of quantum field theory and the relativistic mean field (RMF) approximation [9]. More comprehensive studies have been performed within the framework of Nuclear Many-Body Theory, in which the description of nuclear dynamics is based on a phenomenological Hamiltonian strongly constrained by the observed properties of the two- and three-nucleon systems. Calculations along this line have been carried out using both G-matrix perturbation theory and the variational approach based on the formalism of correlated wave functions and the cluster expansion technique; see, e.g., Refs. [10,11].

The authors of Ref. [12] have developed a procedure to renormalise the coordinate-space nuclear Hamiltonian by introducing screening effects arising from short-range nucleon–nucleon correlations. The resulting density-dependent effective potential—which includes the contributions of both two- and three-nucleon forces—is well behaved and can be employed to carry out perturbative calculations on the basis of the eigenstates of the noninteracting system. The extension of this formalism to the case of a nonzero temperature—involving a proper definition of the grand canonical potential needed to achieve thermodynamic consistency—is based on the assumption that at temperature $T\ll m_{\pi }$, $m_{\pi }\approx 138$ MeV, which is the isospin averaged mass of the $\pi$-meson, thermal effects do not significantly affect strong interaction dynamics [13].

In this paper, we discuss the main features of the approach of Refs. [12,13], as well as its application to a variety of equilibrium and dynamical properties of nuclear matter [14]. The nuclear Hamiltonian and the derivation of the effective interaction are described in Section 2, while the perturbative calculation of the nuclear matter EOS at a finite temperature is outlined in Section 3. Thermal effects on the single-particle properties of charge-neutral $\beta$-stable matter are discussed in Section 5. Finally, Section 6 is devoted to the calculation of the bulk viscosity coefficient, driving the damping of the density oscillations of neutron stars.

2. Nuclear Hamiltonian

Nuclear Many-Body Theory (NMBT) is based on the hypothesis that all nucleon systems—from the deuteron to neutron stars—can be described in terms of point-like protons and neutrons, the dynamics of which are dictated by the Hamiltonian

$$H=\sum _i\frac{p_i^2}{2m}+\sum _{i<j}{v}_{ij}+\sum _{i<j<k}{V}_{ijk},$$

(1)

with m and ${\mathit{p}}_i$ being the mass and momentum of the i-th particle1.

Nucleon–nucleon (NN) potentials that are local or semi-local in coordinate space are usually written in the form

$$v_{ij}=\sum _p{v}^p\left(r_{ij}\right)O_{ij}^p,$$

(2)

where $r_{ij}=|{\mathbf{r}}_i-{\mathbf{r}}_j|$ is the distance between the interacting particles. They are designed to reproduce the measured properties of the two-nucleon system, in both bound and scattering states, and reduce to the Yukawa one-pion-exchange potential at large distances. The sum in Equation (2) includes up to eighteen terms, with the corresponding operators, $O^p$, being required to describe the strong spin–isospin dependence and noncentral nature of nuclear forces ($i=1,\dots ,6$), as well as the occurrence of the spin–orbit and other angular-momentum-dependent interactions ($i=7,\dots 14$). Highly accurate phenomenological potentials, such as the Argonne $v_{18}$ (AV18) model, also feature additional terms accounting for small violations of charge symmetry and charge independence ($i=15,\dots ,18$) [15].

The addition of the three-nucleon (NNN) potential $V_{ijk}$ is needed to model the effects of irreducible three-body interactions, reflecting the appearance of processes involving the internal structure of the nucleons. Nuclear Hamiltonians comprising an AV18 NN potential and a phenomenological NNN potential, designed to explain the binding energies of ³He and ⁴He and the empirical equilibrium density of isospin-symmetric matter—such as the widely used Urbana IX (UIX) model [16,17]—have been shown to possess a remarkable predictive power. The results of Quantum Monte Carlo (QMC) calculations, extensively reviewed in Ref. [18], demonstrate that the AV18 + UIX Hamiltonian is capable of describing the energies of the ground and low-lying excited states of nuclei with mass number $A\le 8$ with an accuracy within a few percent.

2.1. Renormalisation of the Nucleon–Nucleon Interaction

Owing to the presence of a strong repulsive core, the matrix elements of the NN potential between eigenstates of the noninteracting system are large, and standard perturbation theory cannot be used to carry out calculations of nuclear-matter properties.

The renormalisation scheme based on the formalism of correlated basis functions (CBFs) and the cluster expansion technique [19,20] allows one to determine an effective interaction suitable to carry out perturbative calculations in coordinate space. This approach, originally proposed by Cowell and Pandharipande in the early 2000s [21,22] and further developed by the authors of Refs. [23,24,25], has been extensively employed to study nuclear matter and neutron stars [12,13,14]. The resulting potential can be written as in Equation (2), including terms with $i\le 6$, associated with the operators

$$\begin{array}{c}\hfill O_{ij}^{p\le 6}=[1,({\mathsf{\sigma }}_i·{\mathsf{\sigma }}_j),S_{ij}]\otimes [1,({\mathsf{\tau }}_i·{\mathsf{\tau }}_j)].\end{array}$$

(3)

Here, ${\mathsf{\sigma }}_i$ and ${\mathsf{\tau }}_i$ are Pauli matrices acting in spin and isospin space, respectively, while the angular dependence is described by the tensor operator $S_{ij}$, defined as

$$\begin{array}{c}\hfill S_{ij}=\frac{3}{r_{ij}^2}({\mathsf{\sigma }}_i·{\mathbf{r}}_{ij})({\mathsf{\sigma }}_j·{\mathbf{r}}_{ij})-({\mathsf{\sigma }}_i·{\mathsf{\sigma }}_j).\end{array}$$

(4)

Note that the one-pion-exchange potential can also be written in terms of the $O_{ij}^{p\le 6}$ defined by Equations (3) and (4).

The CBF effective interaction is defined by the equation

$$\langle H\rangle = \langle {\Psi }_0|H|{\Psi }_0 \rangle =T_F + \langle {\Phi }_0|\sum _{i<j}{v}_{ij}^{\mathrm{eff}}|{\Phi }_0\rangle ,$$

(5)

where $|{\Phi }_0\rangle$ and $T_F$ denote the ground state of the noninteracting Fermi gas at density $\varrho$ and the corresponding energy, respectively, while H is the nuclear Hamiltonian of Equation (1). The correlated ground state, $|{\Psi }_0\rangle$, is obtained from the corresponding Fermi gas state $|{\Phi }_0\rangle$ through the transformation

$$|{\Psi }_0\rangle \equiv \frac{F|{\Phi }_0\rangle }{\langle {\Phi }_0|F^{†}F{|{\Phi }_0\rangle }^{1/2}},$$

(6)

where the operator F is a symmetrised product of two-body correlation operators, whose structure is chosen in such a way as to reflect the complexity of NN forces.

The effective interaction employed to obtain the results discussed in this paper has been derived following the procedure described in Ref. [24], which allows for the inclusion of the contribution of three-nucleon clusters to the ground-state expectation value $\langle H\rangle$ appearing in the left-hand side of Equation (5). This feature is essential to take into account three-nucleon interactions, which play a dominant role in the high-density regime relevant to astrophysical applications. It should be also pointed out that the nucleon–nucleon correlation functions are determined in such a way as to simultaneously reproduce the ground-state energies of pure neutron matter (PNM) and isospin-symmetric matter (SNM) obtained from highly accurate many-body calculations, carried out using the Auxiliary Field Diffusion Monte Carlo (AFDMC) technique or the Fermi Hyper-Netted Chain/Single-Operator Chain (FHNC/SOC) [12] variational approach. Therefore, the potential $v_{ij}^{\mathrm{eff}}$ defined by Equation (5) can be used to describe nuclear dynamics in nuclear matter with any neutron excess.

2.2. The CBF Effective Interaction

The nuclear Hamiltonian employed to obtain the CBF effective interaction consists of the Argonne $v_6^{\prime }$ (AV6P) NN potential—determined by projecting the full AV18 potential onto the operator basis of Equation (3) [26]—and the UIX NNN potential [17]. The AV6P potential predicts the binding energy and electric quadrupole moment of the deuteron with accuracies of 1%, and 4%, respectively, and provides an excellent fit of the elastic NN scattering phase shifts in the ${}^1{\mathrm{S}}_0$ channel—which is dominant in pure neutron matter—up to lab energy ∼600 MeV, well above the pion production threshold.

The UIX potential is written in the form

$$\begin{array}{c}\hfill V_{ijk}=V_{ijk}^{2\pi }+V_{ijk}^R,\end{array}$$

(7)

where the first term is the attractive Fujita–Miyazawa potential [27]—describing two-pion-exchange NNN interactions with excitation of a $\Delta$-resonance in the intermediate state—while the purely phenomenological $V_{ijk}^R$ is an isoscalar repulsive term, the strength of which is fixed in such a way as to reproduce the saturation density of isospin-symmetric matter inferred from nuclear systematics [16,17].

Note that the CBF effective interaction depends on density through both dynamical correlations, described by the operator F of Equation (6), and statistical correlations, arising from the antisymmetric nature of the state $|{\Phi }_0\rangle$. In Figure 1, the radial dependence of $v^{\mathrm{eff}}$ in the spin–isospin channel $S=0$ and $T=1$ is displayed for baryon densities $\varrho =$ 0.04, 0.32, and 0.48 fm${}^{-3}$. The effect of renormalisation clearly emerges from the comparison with the bare AV6P potential.

Recent studies of the EOS of cold neutron matter—performed by the authors of Ref. [28] using accurate computational techniques—show that the predictions of the somewhat simplified AV6P + UIX Hamiltonian are very close to those obtained from the full AV18 + UIX model employed by Akmal, Pandharipande, and Ravenhall [29].

3. Many-Body Perturbation Theory at Finite Temperature

Let us consider, for simplicity, a one-component Fermi system. The derivation of perturbation theory at a finite temperature is based on the solution of the Bloch equation

$$\begin{array}{c}\hfill -\frac{\partial \Phi }{\partial \beta }=(H-\mu N)\Phi ,\end{array}$$

(8)

where

$$\begin{array}{c}\hfill \Phi \left(\beta \right)=e^{-\beta (H-\mu N)},\end{array}$$

(9)

with the initial condition $\Phi \left(0\right)=1$; see, e.g., Ref. [30]. In the above equations, $\beta =1/T$, while H and $\mu$ denote the Hamiltonian and the chemical potential, respectively.

The perturbative expansion of the grand canonical partition function $Z=\mathrm{Tr}\Phi$ is easily obtained by exploiting the formal similarity between Equation (8) and the time-dependent Schrödinger equation of quantum mechanics and rewriting the Hamiltonian in the form

$$\begin{array}{cc}\hfill H& =H_0+H_I.\hfill \end{array}$$

(10)

The substitution of Equation (10) into the right-hand side of the Bloch equation, leading to

$$\begin{array}{c}\hfill -\frac{\partial \Phi }{\partial \beta }=[(H_0-\mu N)+H_I]\Phi =(H_0^{\prime }+H_I)\Phi ,\end{array}$$

shows that the formalism of time-dependent perturbation theory can be readily generalised by replacing $t\to -i\beta$ and using $H_0^{\prime }$ to define operators in the interaction picture.

The fundamental relation

$$\begin{array}{c}\hfill \Omega =-\frac{1}{\beta }lnZ=-PV=F-\mu N=E-TS-\mu N,\end{array}$$

(11)

where V is the normalisation volume, provides the link between the grand canonical potential $\Omega$, the pressure P, and the free energy $F=E-TS$, with E and S being the energy and entropy of the system, respectively; see, e.g., Ref. [31]. From Equation (11), it follows that

$$\begin{array}{c}\hfill P=-\frac{\Omega }{V},S=-\frac{\partial \Omega }{\partial T},N=-\frac{\partial \Omega }{\partial \mu }.\end{array}$$

(12)

Here, we will discuss the application of the above results to a system described by the Hamiltonian (10), with

$$\begin{array}{c}\hfill H_0=\sum _k{e}_k{a}_k^{†}{a}_k,\end{array}$$

(13)

where, in general

$$\begin{array}{c}\hfill e_k=\frac{{\mathbf{k}}^2}{2m}+U_k=t_k+U_k,\end{array}$$

(14)

and

$$\begin{array}{c}\hfill H_I=\sum _{k,k^{\prime },q,q^{\prime }}\langle k^{\prime }{q}^{\prime }|v|kq\rangle a_{k^{\prime }}^{†}{a}_{q^{\prime }}^{†}{a}_q{a}_k-\sum _k{U}_k{a}_k^{†}{a}_k.\end{array}$$

(15)

Here, the label k specifies both the particle momentum and the discrete quantum numbers corresponding to one-particle states, $a_k^{†}$ and $a_k$ denote creation and annihilation operators, respectively, and v is the potential describing interparticle interactions. The single-particle potential $U_k$, which in principle does not affect the results of calculations of physical quantities, is chosen in such a way as to improve the convergence of the perturbative expansion or fulfil specific conditions; see, e.g., Ref. [32].

It has to be pointed out that, according to Equation (11), the pressure can be written in the form

$$\begin{array}{c}\hfill P=\varrho \left(\mu -\frac{F}{N}\right),\end{array}$$

(16)

with $\varrho =N/V$, implying that at equilibrium, that is, for $P=0$, $\mu =F/N$. This result is the generalisation of the Hugenholtz–Van Hove theorem [33] to the case of a nonvanishing temperature.

It should be emphasised that, when used in conjunction with the CBF effective interaction discussed in Section 2.2, the perturbative approach based on the Hamiltonian defined by Equations (10) and (13)–(15) allows two- and three-nucleon interactions to be taken into account in a fully consistent fashion.

3.1. Perturbative Expansion

At the first order in $H_I$, the grand canonical potential is given by [34]

$$\begin{array}{c}\hfill \Omega ={\Omega }_0+{\Omega }_1,\end{array}$$

(17)

with

$$\begin{array}{cc}\hfill {\Omega }_0& =-\frac{1}{\beta }\sum _{k}ln\left\{1+e^{-\left[\beta (e_k-\mu )\right]}\right\},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\Omega }_1& =\frac{1}{2}\sum _{k{k}^{\prime }}\langle k{k}^{\prime }|v|k{k}^{\prime }{\rangle }_A{n}_k{n}_{k^{\prime }}-\sum _k{U}_k{n}_k,\hfill \end{array}$$

(19)

where $|k{k}^{\prime }{\rangle }_A=|k{k}^{\prime }\rangle -|k^{\prime }k\rangle$ denotes an antisymmetrised two-particle state, and $n_k$ is the Fermi distribution, defined as

$$\begin{array}{c}\hfill n_k={\left[1+e^{\beta (e_k-\mu )}\right]}^{-1}.\end{array}$$

(20)

From Equations (18) and (19) it follows that the free energy per particle,

$$\begin{array}{c}\hfill \frac{F}{N}=\frac{1}{N}({\Omega }_0+{\Omega }_1)+\mu ,\end{array}$$

(21)

can be cast in the form

$$\begin{array}{cc}\hfill \frac{F}{N}& =\frac{1}{N}\{\sum _k{t}_k{n}_k+\frac{1}{2}\sum _{k,k^{\prime }}\langle k{k}^{\prime }|v|k{k}^{\prime }{\rangle }_A{n}_k{n}_{k^{\prime }}\hfill \\ \hfill & +\frac{1}{\beta }\sum _k\left[n_{k}ln{n}_k+(1-n_k)ln(1-n_k)\right]\}+\mu \left(1-\frac{1}{N}\sum _k{n}_k\right).\hfill \end{array}$$

(22)

In principle, for any assigned values of temperature and chemical potential, the above equations provide a scheme for the determination of the equation of state of nuclear matter at a finite temperature, $P=P(\mu ,T)$. However, because the baryon number is conserved by all known interactions, in nuclear matter it is more convenient to use the baryon density as an independent variable and determine the chemical potential from the relation

$$\begin{array}{c}\hfill \varrho =-\frac{1}{V}\frac{\partial }{\partial \mu }\left({\Omega }_0+{\Omega }_1\right).\end{array}$$

(23)

In the $T\to 0$ limit, the momentum distribution reduces to the Heaviside step function $\theta (\mu -e_k)$, with the chemical potential being given by $\mu =e_{k_F}$. The Fermi momentum is related to the particle number density $\varrho$ through $k_F={\left(6{\pi }^2\varrho /\nu \right)}^{1/3}$, where $\nu$ denotes the degeneracy of momentum eigenstates.

3.2. Thermodynamic Consistency

For $T\ne 0$ and density-dependent potentials, thermodynamic consistency is not trivially achieved at any given order of perturbation theory. A clear manifestation of this difficulty is the mismatch between the value of pressure obtained from Equation (16) and the one resulting from the alternative—although, in principle, equivalent—thermodynamic expression

$$\begin{array}{c}\hfill P=-\frac{\partial F}{\partial V}={\varrho }^2\frac{\partial }{\partial \varrho }\frac{F}{N}.\end{array}$$

(24)

A procedure fulfilling the requirement of thermodynamic consistency by construction can be derived from a variational approach, based on the minimisation of the trial grand canonical potential [35]

$$\begin{array}{cc}\hfill \tilde{\Omega }=\sum _k{t}_k{n}_k& +\frac{1}{2}\sum _{k,k^{\prime }}\langle k{k}^{\prime }|v|k{k}^{\prime }{\rangle }_A{n}_k{n}_{k^{\prime }}\hfill \\ \hfill & +\frac{1}{\beta }\sum _k\left[n_{k}ln{n}_k+(1-n_k)ln(1-n_k)\right]-\mu N,\hfill \end{array}$$

(25)

with respect to the form of distribution $n_k$. Note that the above expression—the use of which is fully legitimate in the variational context—can also be obtained in first-order perturbation theory neglecting terms involving $\partial {\Omega }_1/\partial T$ and $\partial {\Omega }_1/\partial \mu$ [34].

The condition,

$$\begin{array}{c}\hfill \frac{\delta \tilde{\Omega }}{\delta n_k}=0,\end{array}$$

(26)

turns out to be satisfied by the distribution function

$$\begin{array}{c}\hfill n_k={\left\{1+e^{\beta [(t_k+U_k+\delta e)-\mu ]}\right\}}^{-1},\end{array}$$

(27)

with

$$\begin{array}{c}\hfill U_k=\sum _{k^{\prime }}\langle k{k}^{\prime }|v|k{k}^{\prime }{\rangle }_A{n}_{k^{\prime }},\end{array}$$

(28)

and

$$\begin{array}{c}\hfill \delta e=\frac{1}{2}\sum _{k,k^{\prime }}\langle k{k}^{\prime }|\frac{\partial v}{\partial \varrho }{|k{k}^{\prime }\rangle }_A{n}_k{n}_{k^{\prime }}.\end{array}$$

(29)

Within the above scheme, which reduces to the standard Hartee–Fock approximation in the case of density-independent potentials, all thermodynamic functions at a given temperature and baryon density can be consistently obtained using the distribution $n_k$ of Equation (27). Note, however, that, because both $U_k$ and $\delta e$ depend on $n_k$ (see Equations (28) and (29)), calculations must be carried out self-consistently, applying an iterative procedure.

4. Equilibrium Properties of Hot Nuclear Matter

Consider now a nucleon system at temperature T, baryon number density, $\varrho$ and proton number density ${\varrho }_p=Y_p\varrho$. At the first order in the CBF effective interaction $v^{\mathrm{eff}}$, the internal energy per nucleon can be written in the form [14]

$$\begin{array}{c}\hfill \frac{E}{N}=\frac{1}{N}\left\{\sum _{\alpha \mathbf{k}}\frac{{\mathbf{k}}^2}{2m}{n}_{\alpha }(k,T)+\frac{1}{2}\sum _{\alpha \mathbf{k}}\sum _{{\alpha }^{\prime }{\mathbf{k}}^{\prime }}\langle \alpha k,{\alpha }^{\prime }{k}^{\prime }|v^{\mathrm{eff}}{|\alpha k,{\alpha }^{\prime }{k}^{\prime }\rangle }_A{n}_{\alpha }(k,T)n_{{\alpha }^{\prime }}(k^{\prime },T)\right\}.\end{array}$$

(30)

In the above equations, the index $\alpha =n,p$ labels neutrons and protons, respectively, $\mathbf{k}$ is the nucleon momentum, $k=|\mathbf{k}|$, and $|\alpha k,{\alpha }^{\prime }{k}^{\prime }{\rangle }_A$ denotes an antisymmetrised two-nucleon state. The conservation of the baryon number obviously implies that $Y_n=1-Y_p$.

The temperature dependence is described by the Fermi distributions

$$\begin{array}{c}\hfill n_{\alpha }(k,T)={\left\{1+exp\left[\beta (e_{\alpha k}-{\mu }_{\alpha })\right]\right\}}^{-1},\end{array}$$

(31)

where the single-particle energies are defined as

$$\begin{array}{c}\hfill e_{\alpha k}=e_{\alpha k}^{\mathrm{HF}}+\delta e,\end{array}$$

(32)

with

$$\begin{array}{cc}\hfill e_{\alpha k}^{\mathrm{HF}}=\frac{{\mathbf{k}}^2}{2m}& +\sum _{{\alpha }^{\prime }{\mathbf{k}}^{\prime }}\langle \alpha k,{\alpha }^{\prime }{k}^{\prime }|v^{\mathrm{eff}}{|\alpha k,{\alpha }^{\prime }{k}^{\prime }\rangle }_A{n}_{\alpha }(k^{\prime },T),\hfill \end{array}$$

(33)

and

$$\begin{array}{c}\hfill \delta e=\frac{\varrho }{2}\sum _{\alpha \mathbf{k}}\sum _{{\alpha }^{\prime }{\mathbf{k}}^{\prime }}\langle \alpha k,{\alpha }^{\prime }{k}^{\prime }|\frac{\partial v^{\mathrm{eff}}}{\partial \varrho }{|\alpha k,{\alpha }^{\prime }{k}^{\prime }\rangle }_A{n}_{\alpha }(k,T)n_{{\alpha }^{\prime }}(k^{\prime },T),\end{array}$$

(34)

The correction to the Hartree–Fock (HF) spectrum is needed to satisfy the requirement of thermodynamic consistency and vanishes in the case of a density-independent potential [13]. The chemical potentials ${\mu }_{\alpha }$ are determined by the normalisation conditions

$$\begin{array}{c}\hfill 2\sum _{\alpha \mathbf{k}}{n}_{\alpha }(k,T)=N_{\alpha },\end{array}$$

(35)

where $N_{\alpha }$ denotes the number of particles of species $\alpha$, the fractional density of which is $Y_{\alpha }=N_{\alpha }/N_B={\varrho }_{\alpha }/\varrho$. Note that the above definitions imply that both the single-nucleon energies and the chemical potentials depend on temperature through the Fermi distribution.

The free energy per nucleon is obtained from

$$\begin{array}{c}\hfill \frac{F}{N}=\frac{1}{N}\left(E-TS\right),\end{array}$$

(36)

with the internal energy per nucleon of Equation (30) and the entropy per nucleon defined as

$$\begin{array}{c}\hfill \frac{S}{N}=-\sum _{\alpha \mathbf{k}}\left\{n_{\alpha }(k,T)ln{n}_{\alpha }(k,T)+\left[1-n_{\alpha }(k,T)\right]ln\left[1-n_{\alpha }(k,T)\right]\right\}.\end{array}$$

(37)

Figure 2 shows the density and temperature dependence of the free energy per nucleon of SNM and PNM, corresponding to proton fraction $Y_p=0.5$ and 0, respectively, obtained from the procedure described above using the CBF effective interaction.

5. Thermal Effects on Nuclear Matter Properties

In the temperature regime discussed in this paper, thermal modifications of nuclear matter properties arise primarily from the Fermi distributions, defined by Equation (31). Comparison to the $T\to 0$ limit shows that the probability distribution $n_{\alpha }(k,T>0)$ is reduced from unity in the region corresponding to ${\mu }_{\alpha }-T\lesssim e{}_{\alpha k}\lesssim {\mu }_{\alpha }$ and acquires nonvanishing positive values for ${\mu }_{\alpha }\lesssim e{}_{\alpha k}\lesssim {\mu }_{\alpha }+T$. It follows that, for any given temperature T, the extent of thermal modifications to the Fermi distribution is driven by the ratio $2T/{\mu }_{\alpha }$. This observation in turn implies that, because the chemical potential is a monotonically increasing function of the particle density ${\varrho }_{\alpha }$ over a broad range of temperatures, for any given T thermal effects turn out to be more significant at lower ${\varrho }_{\alpha }$. However, they become vanishingly small in the high-density regime, in which degeneracy dominates. For example, in SNM the nucleon chemical potential grows from 32 to 144 MeV as the baryon number density increases from 0.32 to 0.48 fm${}^{-3}$ [13].

Charge-Neutral $\beta$-Stable Matter at Finite Temperature

In charge-neutral matter consisting of neutrons, protons, and leptons in equilibrium with respect to the weak interaction processes

$$\begin{array}{c}\hfill n\to p+\mathsf{\ell }+{\overline{\nu }}_{\mathsf{\ell }},p+{\mathsf{\ell }}^{-}\to n+{\nu }_{\mathsf{\ell }},\end{array}$$

(38)

where $\mathsf{\ell }=e,\mu$ labels the lepton flavour, the proton fraction $Y_p$ is uniquely determined by the equations

$$\begin{array}{c}\hfill {\mu }_n-{\mu }_p={\mu }_{\mathsf{\ell }},Y_p=\sum _{\mathsf{\ell }}{Y}_{\mathsf{\ell }}.\end{array}$$

(39)

At densities such that the electron chemical potential does not exceed the rest mass of the muon, $m_{\mu }=105.7$ MeV, the sum appearing in the above equation includes electrons only. At higher densities—typically at $\varrho \gtrsim {\varrho }_0$, with ${\varrho }_0=0.16{\mathrm{fm}}^{-3}$ being the baryon number density of isospin-symmetric matter in thermodynamic equlibrium at $T=0$—the appearance of muons becomes energetically favoured and should be taken into account.

The solid lines of Figure 3 show the density dependence of the proton fractions corresponding to the $\beta$-equilibrium in matter consisting of protons, neutrons, electrons, and muons, or $npe\mu$ matter, at $T=$ 0 (triangles) and 50 MeV (circles) [14]. The results have been obtained assuming that neutrinos do not interact with matter and therefore have vanishing chemical potential. For comparison, the same quantities in $npe$ matter, in which the muon contribution is not included, are displayed by the dashed lines. At temperatures as high as 50 MeV, the assumption that neutrinos are nondegenerate may be questionable, because the neutrino mean free path in matter, $\lambda$, is known to decrease with increasing temperatures. The present analysis should be seen as the first step in a systematic study of neutrino trapping, based on a fully consistent evaluation of $\lambda$, which is currently underway.

The most prominent thermal effect emerging from Figure 3 is a significant departure from the monotonic behaviour observed in cold matter. The emergence of a minimum in the density dependence of the proton fraction—which turns out to be largely unaffected by the appearance om muons—results from the balance between the thermal and degeneracy contributions to the chemical potentials appearing in Equation (39). For $T\gtrsim 20$ MeV and a low density, typically $\varrho \lesssim {\varrho }_0$, the thermal contribution—whose leading order term can be written in the form $\delta {\mu }_{\alpha }\propto T^2/{\varrho }_{\alpha }^{1/3}$—turns out to be much larger for protons than for neutrons, and the $\beta$-equilibrium requires large proton fractions. The results displayed in Figure 3, showing that $Y_p$ does not exceed 25% for ${\varrho }_0/2\le \varrho \le 4{\varrho }_0$, imply that thermal effects in $\beta$-stable matter mainly affect the proton distributions.

The description of single-particle dynamics in interacting many-body systems is largely based on the use of the effective mass $m_{\alpha }^{★}$, defined as

$$\begin{array}{c}\hfill \frac{1}{m_{\alpha }^{★}}={\left(\frac{1}{k}\frac{d{e}_{\alpha k}}{dk}\right)}_{k=k_{F_{\alpha }}},\end{array}$$

(40)

with $e_{\alpha k}$ given by Equation (32). The effective mass dictates the nucleon dispersion relation in matter, which plays a critical role in determining the rates of many processes relevant to neutron star properties, such as the bulk viscosity to be discussed in Section 6. The momentum dependence of the nucleon spectra in cold nuclear matter is often parametrised according to [7]

$$\begin{array}{c}\hfill e_{\alpha k}=\frac{k^2}{2m_0^{★}}+U_{\alpha },\end{array}$$

(41)

where $m_0^{★}$ denotes the value of $m_{\alpha }^{★}$ at $T=0$, while the offset $U_{\alpha }$ from the free-nucleon spectrum—the value of which depends on both the temperature and density—is determined by the requirement that the above approximation reproduce the results of the full microscopic calculation at $k=0$.

The solid lines of Figure 4 show the momentum dependence of the proton spectra in charge-neutral $\beta$-stable matter at baryon density $\varrho =2{\varrho }_0$ and temperature $T=0$ and 50 MeV. A comparison with the dashed lines indicates that at $T=50$ MeV the apprroximation of Equation (41) fails to provide an accurate representation of $e_{\alpha k}$ for $k>0.5{\mathrm{fm}}^{-1}$. However, the quadratic approximation turns out to be in remarkable agreement with the microscopic result if $m_0^{★}$ is replaced with the effective mass computed at $T=50$ MeV. At $\varrho =2{\varrho }_0$ ($3{\varrho }_0$ ), the values of $U_{\alpha }$ turn out to increase from −98 (−50) MeV to −62 (−10) MeV as the temperature grows from 0 to 50 MeV.

6. Bulk Viscosity of Neutron Star Matter

Bulk viscosity in charge-neutral $\beta$-stable matter appears in the aftermath of an instantaneous departure from $\beta$-equilibrium, resulting from the change in density and pressure induced by some perturbation. The role of bulk viscosity in determining the maximum rotation rate of neutron stars—the value of which is limited by the onset of the Chandrasekahr–Friedman–Shutz (CFS) instability—driven by gravitational wave emission [36,37,38]—has been discussed by many authors; for a review, see, e.g., Ref. [39]. More recent studies, motivated by the observation of the gravitational wave event GW170817 [40] focused on the effect of bulk viscosity in neutron star mergers, leading to a damping of density oscillations [41].

6.1. Dissipative Processes in Fluids

The description of fluids is based of the continuity equation expressing conservation of mass

$$\begin{array}{c}\hfill \frac{\partial \varrho }{\partial t}=-\mathbf{\nabla }·\left(\varrho \mathrm{v}\right),\end{array}$$

(42)

where $\mathbf{v}=\mathbf{v}(\mathbf{r},t)$ and $\varrho =\varrho (\mathbf{r},t)$ denote the velocity field and matter density, respectively. For ideal fluids, that is, in the absence of dissipation, the force acting on a fluid element is simply related to the pressure of the surrounding medium through

$$\begin{array}{c}\hfill \mathbf{F}=-\mathbf{\nabla }P,\end{array}$$

(43)

and the equation of motion reduces to Euler’s equation [42]

$$\begin{array}{c}\hfill \frac{\partial \mathbf{v}}{\partial t}+(\mathbf{v}·\mathbf{\nabla })\mathbf{v}=-\frac{1}{\varrho }\mathbf{\nabla }P.\end{array}$$

(44)

The above equations can be combined to obtain

$$\begin{array}{c}\hfill \frac{\partial \left(\varrho {\mathrm{v}}_i\right)}{\partial t}=-{\nabla }_j{\Pi }_{ij},\end{array}$$

(45)

where a sum on the index j is implicit and ${\Pi }_{ij}$ is the momentum flux tensor, defined as

$$\begin{array}{c}\hfill {\Pi }_{ij}=P{\delta }_{ij}+\varrho {\mathrm{v}}_i{\mathrm{v}}_j.\end{array}$$

(46)

In viscous fluids, the form of the continuity equation does not change, but the momentum flux tensor needs to be modified by adding a term describing the deviations from the ideal fluid behaviour, due the occurrence of processes involving irreversible momentum transfer. The resulting expression can be written in the form

$$\begin{array}{c}\hfill {\Pi }_{ij}=P{\delta }_{ij}+\varrho {\mathrm{v}}_i{\mathrm{v}}_j+\delta {\Pi }_{ij},\end{array}$$

(47)

where

$$\begin{array}{c}\hfill \delta {\Pi }_{ij}=-\eta \left[{\nabla }_j{\mathrm{v}}_i+{\nabla }_i{\mathrm{v}}_j+\frac{2}{3}{\delta }_{ij}(\nabla ·\mathbf{v})\right]-\zeta {\delta }_{ij}(\nabla ·\mathbf{v}),\end{array}$$

(48)

with the two velocity-independent quantities $\eta$ and $\zeta$ being referred to as shear and bulk viscosity coefficients, respectively. In the case of incompressible fluids, in which $\varrho$ does not depend on either $\mathbf{r}$ or t, Equation (42) implies

$$\begin{array}{c}\hfill \mathbf{\nabla }·\mathrm{v}=0,\end{array}$$

(49)

and the contribution of bulk viscosity vanishes.

A pulsation of frequency $\omega$ induces a variation of the fluid density described by the equation

$$\begin{array}{c}\hfill \varrho \left(t\right)={\varrho }_{\mathrm{eq}}+\delta \varrho cos\omega t,\end{array}$$

(50)

where ${\varrho }_{\mathrm{eq}}$ is the density corresponding to chemical equilibrium and $(\delta \varrho /{\varrho }_{\mathrm{eq}})\ll 1$. The rate of energy density dissipation due to bulk viscosity, averaged over the pulsation period $\tau =2\pi /\omega$, is given by [43]

$$\begin{array}{c}\hfill 〈\frac{d{ϵ}_{\mathrm{diss}}}{dt}〉=-\frac{1}{\tau }{\int }_0^{\tau }dt{\nabla }_i\left[-v_j\zeta {\delta }_{ij}(\mathbf{\nabla }·\mathbf{v})\right]=\frac{1}{\tau }{\int }_0^{\tau }dt\zeta {(\mathbf{\nabla }·\mathbf{v})}^2.\end{array}$$

(51)

6.2. Bulk Viscosity of $\beta$-Stable Matter

The authors of Refs. [44,45] have derived the expression of the bulk viscosity coefficient of neutrino-transparent $\beta$-stable $npe$ matter under the assumption that neutrinos and antineutrinos are produced through the processes

$$\begin{array}{cc}\hfill p+e& \to n+{\nu }_e,\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill n& \to p+e+{\overline{\nu }}_e,\hfill \end{array}$$

(53)

referred to as direct Urca reactions. They followed the scheme originally proposed in Ref. [46] and employed a simple parametrisation of the nuclear matter EOS assuming a quadratic dependence on the neutron excess $\alpha =1-2Y_p$.

In direct Urca processes leading to the emission of nondegenerate neutrinos, the requirement of momentum conservation entails a relation involving the momenta of the degenerate fermions, which can be exploited to define a threshold by applying the Fermi surface approximation. While being admittedly inconsistent at $T>0$—and certainly not expected to be accurate at $T=50$ MeV—in our work this definition is used as a baseline to analyse thermal modifications of the onset of the Urca process originating from the modifications of matter composition illustrated in Figure 3.

A comprehensive and detailed analysys of the appearance of bulk viscosity associated with the occurrence of Urca processes in $npe$ and $npe\mu$ matter—including both the neutrino-transparent and neutrino-trapped regimes—has recently been carried out by the authors of Ref. [47].

Bulk viscosity is associated with deviations from beta equilibrium, signalled by a nonvanishing difference between the neutrino and antineutrino production rates

$$\begin{array}{c}\hfill \Delta \Gamma ={\Gamma }_{{\nu }_e}-{\Gamma }_{{\overline{\nu }}_e},\end{array}$$

(54)

with $\Delta \Gamma$ being a function of the variable $\delta \mu$, describing the departure from chemical equilibrium. In the case of nondegenerate neutrinos,

$$\begin{array}{c}\hfill \delta \mu ={\mu }_n-{\mu }_p-{\mu }_e.\end{array}$$

(55)

Under the assumption that $\delta \mu \ll T\ll {\mu }_i$, $\Delta \Gamma$ can be expanded in powers of ${\delta }_{\mu }$. At the leading order one finds

$$\begin{array}{c}\hfill \Delta \Gamma =\lambda {\delta }_{\mu },\end{array}$$

(56)

with

$$\begin{array}{c}\hfill \lambda =2{\left(\frac{\partial {\Gamma }_{\nu }}{\partial \delta \mu }\right)}_{\delta \mu =0}.\end{array}$$

(57)

The energy density dissipation rate—averaged over the period of the density fluctuation of Equation (50)—can be written in terms of the change in pressure associated with the departure from chemical equilibrium using

$$\begin{array}{cc}\hfill 〈\frac{d{ϵ}_{\mathrm{diss}}}{dt}〉& =\frac{1}{\varrho }\frac{1}{\tau }{\int }_0^{\tau }dt\delta P_{\mathrm{chem}}\frac{d\varrho }{dt}\hfill \\ \hfill & =\frac{\delta \varrho }{\varrho }\frac{\omega }{\tau }{\int }_0^{\tau }dt\delta P_{\mathrm{chem}}sin\omega t,\hfill \end{array}$$

(58)

with [46]

$$\begin{array}{c}\hfill \delta P_{\mathrm{chem}}=-\lambda C^2\frac{\delta \varrho }{\varrho }\frac{\omega }{{\omega }^2+{(2\lambda B/\varrho )}^2}sin\omega t+\dots \end{array}$$

(59)

where the ellipsis refers to additional terms giving vanishing contributions to the integral. Here, $\lambda$ is given by Equation (57), and the constants B and C are defined as

$$\begin{array}{c}\hfill B={\left(\frac{\partial \delta \mu }{\partial \alpha }\right)}_{\delta \varrho =0},C={\left[\varrho \left(\frac{\partial \delta \mu }{\partial \varrho }\right)\right]}_{\alpha ={\alpha }_{\mathrm{eq}}},\end{array}$$

(60)

where ${\alpha }_{\mathrm{eq}}$ is the value of neutron excess corresponding to the $\beta$-equilibrium. The substitution of Equation (59) into Equation (58) yields

$$\begin{array}{c}\hfill 〈\frac{d{ϵ}_{\mathrm{diss}}}{dt}〉=-\lambda \frac{{\omega }^2}{2}{\left(\frac{\delta \varrho }{\varrho }\right)}^2\frac{C^2}{{\omega }^2+{(2B\lambda /\varrho )}^2}.\end{array}$$

The bulk viscosity coefficient is obtained combining the above equation with Equation (51), which can be rewritten using the continuity equation associated with the conservation of the baryon number, implying

$$\begin{array}{c}\hfill \mathbf{\nabla }·\mathbf{v}=-\frac{1}{\varrho }\frac{\partial \varrho }{\partial t}=\omega \frac{\delta \varrho }{\varrho }sin\omega t.\end{array}$$

(61)

Following this procedure, one finally arrives at

$$\begin{array}{c}\hfill \zeta =-\lambda \frac{C^2}{{\omega }^2+{(2B\lambda /\varrho )}^2}.\end{array}$$

(62)

The above result shows that the calculation of the bulk viscosity coefficient involves $\lambda$, which is obtained from the neutrino and antineutrino production rates, and the constants B and C defined by Equation (60). We have studied the density and temperature dependence of $\zeta$ using the nuclear matter model described in the previous sections, which allows us to take into account thermal modifications of the chemical potentials and effective masses within a fully consistent framework.

6.3. Calculation of the Bulk Viscosity Coefficient

The bulk viscosity coefficient $\zeta$ has been obtained from Equation (62), with B and C given by Equation (60). The calculation of the coefficient $\lambda$, defined by Equations (56) and (57), involves the rate of production of antineutrinos of flavour ℓ in neutron decay processes, which can be written in the form

$$\begin{array}{cc}\hfill {\Gamma }_{{\overline{\nu }}_{\mathsf{\ell }}}=\int \frac{d^3{p}_n}{{\left(2\pi \right)}^3}& \frac{d^3{p}_p}{{\left(2\pi \right)}^3}\frac{d^3{p}_{\mathsf{\ell }}}{{\left(2\pi \right)}^3}\frac{d^3{p}_{\nu }}{{\left(2\pi \right)}^3}{f}_n(1-f_p)(1-f_{\mathsf{\ell }})\hfill \\ \hfill & \times {\left(2\pi \right)}^4\delta (E_n-E_p-E_{\mathsf{\ell }}-E_{\nu }){\delta }^{\left(3\right)}({\mathit{p}}_n-{\mathit{p}}_p-{\mathit{p}}_{\mathsf{\ell }}-{\mathit{p}}_{\nu })\sum _{\mathrm{spins}}{|M|}^2,\hfill \end{array}$$

(63)

as well as the rate of neutrino production associated with lepton capture by protons. In the above equation, ${\mathit{p}}_j$—with $j=n,p,\mathsf{\ell },\nu$—denotes the particle momentum, while the $f_j$ are the Fermi distributions of degenerate particles, defined by Equation (20). Treating the nucleons as nonrelativistic and neglecting ${\mathit{p}}_{\nu }$ in the argument of the momentum-conserving $\delta$-function, as justified for thermal neutrinos, the spin-summed squared transition amplitude can be written in the form

$$\begin{array}{c}\hfill \sum _{\mathrm{spins}}{|M|}^2=2G^2(g_V^2-3g_A^2).\end{array}$$

(64)

Here, $G=G_{F}cos{\theta }_c$, with $G_F$ and ${\theta }_c$ being the Fermi constant and the Cabibbo mixing angle, respectively, while $g_V=1$ and $g_A=-1.26$ denote the vector and axial-vector coupling constants.

The above transition probability is a constant, which can be moved out of the integral appearing in the right-hand side of Equation (63). The integrations can then be performed either numerically, as described by the authors of Ref. [48], or applying the widely used approximation known as phase space decomposition, thoroughly discusssed in, e.g., Ref. [49]. Within this scheme, the magnitudes of the momenta of the strongly degenerate fermions—that is, the nucleons and the charged lepton—are replaced by the corresponding Fermi momenta, and the effective masses, defined as in Equation (40), are used to rewrite the integration measures in the form

$$\begin{array}{c}\hfill d^3{p}_j=d{\Omega }_j{p}_{Fj}{m}_j^{★}d{E}_j,\end{array}$$

(65)

where $j=n,p,\mathsf{\ell }$, and ${\Omega }_j$ is the solid angle specifying the direction of the vector ${\mathit{p}}_j$.

In this work, the calculation of the rates of neutrino and antineutrino production in the Urca processes (52) and (53) have been performed using an improved version of the phase space decomposition. The important new feature of this procedure lies in the inclusion of thermal modifications of the effective masses, which were found to be sizeable by the authors of Ref. [14]. Note that this entails a temperature dependence of the integration measures. Thermal modifications of the chemical potentials, determining both the Fermi distributions and matter composition, have also been taken into account.

Equation (62) can be written in a somewhat more transparent form in terms of the equilibration rate $\gamma =-2B\lambda /\varrho$. The resulting expression,

$$\begin{array}{c}\hfill \zeta =\varrho \frac{C^2}{B}\frac{1}{2}\frac{\gamma }{{\omega }^2+{\gamma }^2},\end{array}$$

(66)

exhibits a resonant maximum located at $\gamma =\omega$. This property, which depends on both the density and temperature, has been thoroughly analysed by the authors of Ref. [41], who also discussed the differences between the isothermal and adiabatic treatment of thermodynamic quantities.

Under the assumption that during a neutron star merger the heat flow between adjoining fluid elements is negligible, the derivatives involved in the calculations of the quantities B and C defined by Equation (60) must be computed keeping the entropy per baryon constant. The alternative procedure based on isothermal derivation, while being equivalent in the $T\to 0$ limit, leads to sizeably different results in the high-temperature regime.

Figure 5 shows the behaviour of $\zeta$ as a function of density for temperatures 10, 30, and 50 MeV, with the frequency of the density oscillation driving the appearance of viscosity being set to $\omega =2\pi \times 1\mathrm{kHz}$, a value typical of oscillations occurring in neutron star mergers [6]. The results were obtained by performing the derivatives of Equation (60) using both the isothermal (solid lines and open circles) and adiabatic (dashed lines) definitions.

The peculiar density dependence featured by the solid line corresponding to the highest temperature, $T=50$ MeV, results from the occurrence of a minimum in the proton fraction $Y_p\left(\varrho \right)$—clearly visible in Figure 3—and from the fact that the value of $Y_p$ exceeds the threshold for the onset of direct Urca processes at all densities. A comparison to the corresponding dashed line shows that the minimum of $\zeta \left(\varrho \right)$ is entirely washed out when the isothermal derivative is replaced by the adiabatic one. Because the proton fraction is still above the Urca threshold, however, one finds $\zeta \left(\varrho \right)\ne 0$ over the whole density range of the figure.

A similar pattern is displayed by the solid line and the open circles representing the results obtained at $T=30\mathrm{MeV}$, although in this case the minimum of $Y_p\left(\varrho \right)$ is not reflected by a minimum in $\zeta \left(\varrho \right)$. One only observes an isolated point with $\zeta \ne 0$ located at $\rho =0.25{\rho }_0$ and a density region extending up to $\rho >2.25{\varrho }_0$ in which $\zeta =0$. At larger densities, $Y_p$ is always above the threshold of Urca processes, and $\zeta \ne 0$. At $T=10$ MeV, the solid and dashed lines turn out to lie on top of one another.

It should be noted that, although the Urca process is always allowed at high temperatures, the corresponding values of $\zeta$ are very small, and the bulk viscosity is unlikely to be distinguished from other dissipative phenomena active in a neutron star merger.

The temperature dependence of the bulk viscosity coefficient at $\omega =2\pi \times 1\mathrm{kHz}$ and densities $\varrho =$ 2.5, 3, and 4 ${\varrho }_0$ is illustrated in Figure 6. The maximum at $T\approx 2\mathrm{MeV}$, whose value monotonically increases with density, is clearly apparent. Owing to the low-temperature range spanned by the figure, in this instance the results obtained using isothermal and adiabatic derivatives turn out to be nearly identical. For this reason, only isothermal results are shown.

7. Discussion

The description of neutron star mergers, which will be needed for the interpretation of current and future astronomical observations, requires a quantitative understanding of both the equilibrium and dynamical properties of hot and dense matter. In this context, the availability of a dynamical model strongly constrained by phenomenology and suitable for use in finite-temperature perturbation theory will be crucial.

The approach described in this paper has been employed to obtain the EOS of matter with arbitrary neutron excess in the density region extending up to 4 ${\varrho }_0$—in which the applicability of the description in terms of nucleons is supported by electron-nucleus scattering data [50]—and temperatures up to 50 MeV. In this density regime—relevant to the description of neutron stars of mass compatible with the observational upper bound—the nonrelativistic approximation is expected to be justified, because the speed of sound in matter predicted by our approach never exceeds the speed of light.

Single-nucleon properties, such as the quasipartlcle spectra and effective masses, have also been computed within the same theoretical framework, in which thermal effects are consistently taken into account. Exploratory studies of dissipative processes leading to the damping of density oscillations in neutron star mergers suggest that a detailed treatment of thermal effects is, in fact, needed to clarify important issues, such as the onset of the Urca process.

The capability of consistently including muons is an important feature of our description of neutron star matter. However, in view of the fact that the most striking thermal effect emerging from the present study—that is, the breakdown of the monotonic density dependence of the proton fraction—is not significantly affected by the presence of muons (see Figure 3), the present analysis limited the discussion of bulk viscosity in the simpler case of $npe$ matter. The role of muons in neutron star mergers is thoroughly discussed in Ref. [51].

The extension of the approach outlined in this paper to the description of neutrino emission processes and the neutrino mean free path along the line described in Refs. [24,25] does not involve any conceptual issues and is currently under way. In principle, the description of Urca processes involving hyperons [52] may also be possible. However, it would require a generalisation of the effective Hamiltonian allowing for a consistent description of all baryonic interactions. The development of such a model is beyond the scope of our work.