Beta equilibrium under neutron star merger conditions

We calculate the nonzero-temperature correction to the beta equilibrium condition in nuclear matter under neutron star merger conditions, in the temperature range $1\,$MeV$<T \lesssim 5\,$MeV. We improve on previous work by using a consistent description of nuclear matter based on the IUF and SFHo relativistic mean field models. This includes using relativistic dispersion relations for the nucleons, which we show is essential in these models. We find that the nonzero-temperature correction can be of order $10$ to $20\,$MeV, and plays an important role in the correct calculation of Urca rates, which can be wrong by factors of $10$ or more if it is neglected.


I. INTRODUCTION
Nuclear matter in neutron stars settles into beta equilibrium, meaning that the proton fraction is in equilibrium with respect to the weak interactions.In this paper we will study the conditions for beta equilibrium in ordinary nuclear matter (where all the baryon number is contributed by neutrons (n) and protons (p)) in the temperature range 1 MeV ≲ T ≲ 5 MeV.This regime, which arises in neutron star mergers [1][2][3][4], is cool enough so that neutrinos are not trapped, but warm enough so that there are corrections to the low-temperature equilibrium condition.It has previously been shown [5] that in this regime the full beta equilibrium condition is where ∆µ is a correction that arises from the violation of detailed balance (neutrino transparency) and the breakdown of the Fermi surface approximation (see Sec. II).In nuclear matter in the temperature regime discussed here, the proton fraction will equilibrate towards the value given by Eq. (1).Even if equilibrium is not reached on the timescale of a merger, one needs to know the correct equilibration condition in order to analyze phenomena associated with this relaxation process, such as bulk viscosity and neutrino emission.At low temperatures (T ≪ 1 MeV) ∆µ is negligible, but in the temperature regime under consideration here it has been estimated to be up to tens of MeV [5].The calculation in Ref. [5] went beyond the Fermi-surface approximation by performing the phase space integral for the equilibration rate over the entire momentum space.However, it used a very crude model of the in-medium nucleons, assigning them their vacuum mass and assuming that their kinematics remained nonrelativistic at all densities.
In this paper we improve on the analysis of Ref. [5].We treat nuclear matter consistently using relativistic mean field models [6,7] with fully relativistic dispersion relations for the nucleons.We show that this makes a considerable difference to the beta equilibration rates because in these models the nucleons at the Fermi surface become relativistic at densities of a few times nuclear saturation density n 0 .We calculate the direct Urca rate using the entire weak-interaction matrix element rather than its non-relativistic limit, and evaluate the full phase space integral.
Other authors have evaluated direct Urca phase space integrals in calculations of the direct Urca rate, the neutrino emissivity, or the neutrino mean free path.Fully relativistic computations of direct Urca phase space integrals are uncommon in the literature, but they do appear.Refs.[8][9][10][11] calculate the neutrino mean free path using a fully relativistic formalism, while integrating over the full phase space.Ref. [10] calculates the direct Urca electron capture rate using a fully relativistic formalism and performs the full phase space integration.While these calculations perform the full integration over phase space, they focus on high temperatures (T ≳ 5 MeV) where neutrinos are trapped and where the direct Urca threshold is blurred over a wide density range.In this temperature regime, which can be reached in mergers as well [1,[12][13][14], beta equilibrium is given by with µ ν being the neutrino chemical potential.As discussed in more detail in Sec.II, the neutrino trapped beta equilibration condition does not require an additional finite temperature correction.This paper will examine the phase space integral at lower temperatures where the direct Urca threshold is apparent and a key feature in the physics of beta equilibration or neutrino emission.
Other works use the relativistic formalism, but assume the nuclear matter is strongly degenerate (using the Fermi surface approximation, described below), and thus their results have a sharp direct Urca threshold density [15][16][17].Ref. [18] uses the Fermi surface approximation, but develops a way to incorporate the finite 3-momentum of the neutrino, slightly blurring the threshold at nonzero temperature.Some works do the full phase space integration, but use nonrelativistic approximations for the matrix element and nucleon dispersion relations [5,[19][20][21].The vast majority of calculations use non-relativistic approximations of the matrix element and the nucleon dispersion relations, together with the Fermi surface approximation [22][23][24][25][26][27][28][29][30][31][32][33][34].All of these calculations are approximations of the full phase space integration using the fully relativistic formalism.Under certain conditions, the approximations match well with the full calculation, and have the advantage of being simple.
In Sec.III we introduce the two relativistic mean field models, IUF and SFHo, that we use.Sec.IV describes our calculation of the rate of direct Urca processes, where we integrate over the entire phase space in order to include contributions from the region that would be kinematically forbidden in the low temperature limit.Sec.V describes our calculation of the modified Urca contribution to the rate, where we use the Fermi surface approximation since there is no kinematically forbidden region for those processes in the density range that we consider.Sec.VI presents our results, and Sec.VII provides our conclusions.
We work in natural units, where ℏ = c = k B = 1.

II. BETA EQUILIBRATION
Beta equilibration in npe − matter is established by the Urca processes [35].The modified Urca processes N + p + e − → N + n + ν, (here, N represents a "spectator" neutron or proton) operate at all densities in the core of the neutron star1 .The direct Urca processes p + e − → n + ν are exponentially suppressed when the temperature is much less than the Fermi energies and the density is in the range where k F n > k F p + k F e .In nuclear matter, the proton fraction rises as the density rises above n 0 and eventually may reach a "direct Urca threshold" where Above this threshold density beta equilibration is dominated by direct Urca, since (when kinematically allowed) it is faster than modified Urca.
In nuclear matter at temperatures greater than, say, 10 MeV, the neutrino mean free path is short and the nuclear matter system (for example, a protoneutron star) is neutrino-trapped and has conserved lepton number Y L = (n e + n ν )/n B .In this case, the Urca processes (3) and ( 4) can proceed forward and backward, as the nuclear matter contains a population of neutrinos (or antineutrinos).In beta equilibrum, the forward and reverse processes have equal rates (detailed balance) and the beta equilibrium condition is given by balancing the chemical potentials of the participants in the equlibration reactions [6,37] In cooler nuclear matter, at the temperatures considered in this work, the neutrino mean free path is comparable to or longer than the system size and therefore neutrinos are not in thermodynamic equilibrium: they escape from the star.Neutrinos can then occur in the final state but not the initial state of the Urca processes.Beta equilibrium is still achieved, but now by a balance of the neutron decay and the electron capture processes.However, the principle of detailed balance is not applicable because electron capture is not the time-reverse of neutron decay.There is then no obvious equilibrium condition that can be written down a priori.
In the limit of low temperature (T ≪ 1 MeV) the Fermi surface approximation becomes valid: the particles participating in the Urca processes are close to their Fermi surfaces, and the neutrino carries negligible energy ∼ T .The beta equilibrium condition can then be obtained by neglecting the neutrino, so that neutron decay and electron capture are just different time orderings of the same process n ↔ p e − , and detailed balance gives2 µ n = µ p +µ e (low temperature, ν-transparent).(6) At temperatures T ≳ 1 MeV corrections to the Fermi surface approximation start to become significant, particularly for the protons whose Fermi energy is in the 10 MeV range.Then one cannot neglect the nonzerotemperature correction to (6) The correction ∆µ is a function of density and temperature, and its value in beta equilibrium is found by explicitly calculating the neutron decay and electron capture rates and adjusting ∆µ so that they balance [5] (see also [39], where a similar calculation was done in the context of a hot plasma).In this paper we perform that calculation.
For weak interactions we use the Fermi effective theory, which is an excellent approximation at nuclear energy scales.The main approximations arise in our treatment of the strong interaction.To describe nuclear matter and the nucleon excitations we use two different relativistic mean field models, both consistent with known phenomenology and chosen to illustrate a plausible range of behaviors.We describe these models in Sec.III.For the modified Urca process we model the nucleon-nucleon interaction with one-pion exchange [31,40].

III. NUCLEAR MATTER MODELS
We will use two different equations of state, IUF [41] and SFHo [42], to calculate the Urca rates and the nonzero-temperature correction ∆µ.These are both consistent at the 2σ level with observational constraints on the maximum mass and the radius of neutron stars.
For the radius of a star of mass 2.06 M ⊙ , SFHo predicts R = 10.3 km, consistent with R = 12.39 +1.30 −0.98 km from NICER and XMM analysis of PSR J0740+6620 [43].For the radius of a 1.4 M ⊙ neutron star, IUF predicts R = 12.7 km and SFHo predicts R = 11.9 km, consistent with R = 11.94 +0.76  −0.87 km obtained by a combined analysis of X-ray and gravitational wave measurements of PSR J0740+6620 in Ref. [46].Both models reproduce properties of nuclear matter around saturation within known experimental constraints.Although both models predict a slightly different value for the saturation density n 0 , we will use n 0 = 0.16 fm −3 in all our plots.
It is still not determined whether there is a direct Urca threshold or not in nuclear matter at neutron star densities [47][48][49][50][51], so we choose one equation of state (IUF) with a threshold slightly above 4 n 0 and one (SFHo) with no threshold, as shown in Fig. 1.Our approach could be applied to any equation of state where the beta process rates can be calculated.As we will see in Sec.VI C, the density dependence of the momentum surplus k F p + k F e − k F n is an important factor in the behavior of the direct Urca rates at low temperature, but the density dependence of the nucleon effective masses and Fermi momenta has a noticeable impact as well.
The coupling constants for SFHo are shown in Appendix A. Notice that the constants are taken from the online CompOSE database (https://compose.obspm.fr/), and are different from the values provided in Ref. [42].
A key feature of our calculation is that we use the full relativistic dispersion relations for the nucleons.In Figs.
2 and 3 we illustrate the importance of this in relativistic mean field theories, where the nucleon effective mass drops rapidly with density 3 .We plot the Dirac effective mass [56] and the Fermi momentum of the neutrons and protons in these two EoSs.While around nuclear saturation density n 0 , the nucleons are nonrelativistic, as the density rises to several times n 0 , the nucleon effective mass has dropped significantly below its vacuum value.Neutrons on their Fermi surface become relativistic at 2 − 3n 0 , while protons on their Fermi surface remain nonrelativistic until the density rises to 3 − 6n 0 .In Figs. 13 and 14, we show that using a non-relativistic approximation would lead to Urca rates that are incorrect by about an order of magnitude, although for direct Urca neutron decay the discrepancy can be many orders of magnitude.

IV. BETA EQUILIBRATION VIA DIRECT URCA
We calculate the in-medium direct Urca rates for neutron decay and electron capture using the relativistic weak-interaction matrix element and the relativistic dispersion relations for the nucleons and electrons.We also integrate over the full momentum phase space, not relying on the Fermi surface approximation.This is impor- 3 While the precipitous drop in the nucleon Dirac effective mass with increasing density is a common feature in relativistic mean field theories [52,53], we note that in two recent treatments that go beyond the mean field approximation, the drop in the effective mass was not as dramatic [54,55].tant because in the "dUrca-forbidden" density range the Fermi surface approximation would say the direct Urca rate is zero, so nonzero-temperature corrections are the leading contribution.These become significant (comparable to modified Urca) at the temperatures of interest here, T ≳ 1 MeV [5].
In relativistic mean field models the dispersion rela-tions for the neutrons, protons, and electrons are where the nucleons' effective mass m * i and energy shift U i depend on density and temperature [10].The unshifted energies E * i arise in the phase space normalization and the Dirac traces [9].

A. Neutron decay
The direct Urca neutron decay rate is [31,57] For a more detailed explanation of this expression and its evaluation, see Appendix B. As described there, it can be reduced to 5-dimensional momentum integral (B34) where R, S, and M ϕ0 are defined in Eqs.(B15), (B16), and (B17).The antineutrino energy E ν is given by which becomes a function of the remaining integration variables, k n , k p , and k e .Note that there are Fermi-Dirac distributions for the neutrons, proton vacancies, and electron vacancies, but none for the neutrinos because we work in the neutrino-transparent regime where neutrinos escape from the star and do not form a Fermi gas.We evaluate this integral numerically using a Monte-Carlo algorithm.

B. Electron capture
The expression for the electron capture rate can be obtained from that for neutron decay (B1) by making the following changes: (1) the energy-momentum delta function now corresponds the the process p e − → n ν, and (2) there are Fermi-Dirac distributions for proton and electron particles, and neutron vacancies, Evaluating this expression takes us through the same steps as for neutron decay, except that the neutrino energy is now and the requirement that this be positive leads to different limits on the momentum integrals, We calculate the rate of the modified Urca processes (3) using the relativistic dispersion relations of the nucleons in the phase space integration, but unlike the direct Urca rate we do not perform the phase space integration exactly, which would be difficult because the involvement of the spectator particles would lead to an 11-dimensional numerical integral over momentum.Instead we use the Fermi Surface approximation.This is reasonable for modified Urca as long as the Fermi surfaces are not too thermally blurred, i.e. when the temperature is below the lowest Fermi kinetic energy, which is that of the proton.The modified Urca processes do not have a density threshold in the range of densities we consider here (see footnote 1), so the Fermi surface approximation never predicts a vanishing rate.In this work we explore the temperature range 1 MeV < T < 5 MeV, and the proton's Fermi kinetic energy is at least 10 MeV in the density range n > n 0 , so the Fermi surface approximation is justified for modified Urca rates.The first paragraph of Sec.IV contains a discussion of why we need to go beyond the Fermi surface approximation in our direct Urca rate calculations.For the matrix elements that arise in modified Urca (C1) and (C16), we use the standard results (see, e.g., [31]), which were calculated assuming non-relativistic nucleons.It has been pointed out [58] that the standard calculation of the modified Urca matrix element [40], which we use here, is based on a very crude approximation for the propagator of the internal off-shell nucleon.A more accurate treatment would lead to different modified Urca rates and shift our predicted values of ∆µ; we defer such a calculation to future work.

A. Neutron Decay
Modified Urca can proceed with either a neutron spectator or a proton spectator.From Fermi's Golden rule, we have the rate for the neutron decay process Here, s = 1/2 because of the identical particles appearing in the process.N 1 and N 2 are neutrons in the nspectator process and for the p-spectator neutron decay process, N 1 and N 2 are protons.The matrix element is different for each process (see Eqs. C1 and C16).The detailed derivation of the modified Urca rates is in Appendix C. For n-spectator neutron decay, allowing the system to deviate from the low-temperature beta equilibrium condition ( 6) by amount we obtain where f ≈ 1 is the N -π coupling [31], and The functions Li n (x) are polylogarithms of order n [59].
For p-spectator neutron decay, we obtain where

B. Electron Capture
The electron capture modified Urca rate can be obtained in a similar way to neutron decay, by changing the sign of the neutrino 4-momentum in the energymomentum delta function and interchanging the particle and hole Fermi-Dirac factors, Through a similar calculation, we find that the modified Urca neutron decay and electron capture rates in the Fermi surface approximation are related by and VI. RESULTS A. Beta equilibrium at nonzero temperature Figs. 4 and 5 show our final results for the nonzerotemperature correction ∆µ required to achieve beta equilibrium, for the IUF and SFHo equations of state respectively.The key features are • At low temperatures T ≲ 1 MeV, the Fermi Surface approximation is valid and beta equilibrium is achieved with a negligible correction ∆µ (see Sec. II).
• As the temperature rises through the neutrinotransparent regime, the value of ∆µ rises.Nonzero-temperature correction ∆µ required for beta equilibrium (Eq.7) with the SFHo EoS.

T=1 MeV
• We only provide results for temperatures up to 5 MeV because at temperatures of around 5 to 10 MeV the neutrino mean free path will become smaller than the star, invalidating our assumption of neutrino transparency.
• The figures indicate that the nonzero-temperature correction reaches values of 10 to 20 MeV before neutrino trapping sets in.
• The density dependence of ∆µ appears very different for different EoSs.For IUF the largest values are reached at moderate densities, near the direct Urca threshold.For SFHo, ∆µ has a minimum at those densities.
In the rest of this section we will explain these features of our results.] Urca rates at T=3MeV FIG. 6. Urca (direct plus modified) rates for IUF and SFHo EoSs at T = 3 MeV.When ∆µ = 0 (dashed lines) the rates for neutron decay (nd) and electron capture (ec) do not balance.With the correct choice of ∆µ (Figs. 4, 5) the neutron decay and electron capture rates (solid lines) become equal and the system is in beta equilibrium.
The temperature dependence follows from the breakdown of the Fermi surface approximation.At T ≲ 1 MeV the Urca processes are dominated by modes close to the Fermi surfaces of the neutron, protons, and electrons.The energy of the emitted neutrino is of order T which is negligible, so the direct Urca process is effectively n ↔ p e − , for which the equilibrium condition is µ n = µ p + µ e , i.e. ∆µ = 0.As the temperature approaches the Fermi energy of the protons, the Fermi surface approximation breaks down.Modes far from the proton and electron Fermi surfaces begin to play a role, and the energy of the emitted neutrino becomes important.The processes that establish beta equilibrium, n → p e − νe and p e − → n ν e , are not related by time reversal, so the principle of detailed balance does not apply.This means that even below the direct Urca threshold density, direct Urca processes can be fast enough and sufficiently different in their rates to require a correction ∆µ to bring them into balance.As we will explain below, at ∆µ = 0 electron capture is much less suppressed than neutron decay, requiring a positive value of ∆µ to decrease the proton fraction and equalize the rates.
The density dependence of the correction ∆µ is more complicated, depending on specific features of the equations of state.We will discuss this in more detail below.

B. Urca Rates
Fig. 6 illustrates how, without the nonzerotemperature correction ∆µ (dashed lines), the neutron decay (nd) and electron capture (ec) rates become very different when the temperature rises to 3 MeV.For both EoSs, electron capture is significantly faster than neutron decay, so a positive ∆µ will be required to balance the rates and establish beta equilibrium (solid lines).This is because a positive ∆µ reduces the proton fraction.The resultant change in the phase space near the neutron and proton Fermi surfaces enhances the neutron decay rate and suppresses electron capture, bringing the two processes into balance with each other.
For IUF, the mismatch between electron capture and neutron decay is greatest just below the IUF direct Urca threshold density of 4 n 0 , which explains why for IUF ∆µ reaches its highest value there (Fig. 4).For SFHo, the mismatch is smallest at that density, which explains why for SFHo ∆µ reaches a local minimum there (Fig. 5).] Urca rate (IUF), T=3MeV, Δμ=0 FIG. 7. Urca rates calculated using the IUF EoS at T = 3 MeV.Because ∆µ = 0 there is a large mismatch between the direct Urca rates for neutron decay and electron capture.Modified Urca (with neutron spectator (n) and proton spectator (p)) rates are calculated in the Fermi surface approximation and therefore match automatically.
Figs. 7 and 8 give further insight in to the density dependence of the rates by showing the separate contributions from direct and modified Urca.
For IUF (Fig. 7), in the dUrca-forbidden density range one would expect that the direct Urca rates should be exponentially suppressed at low temperature, leaving the modified Urca rates which automatically balance when ∆µ = 0 because they are calculated in the Fermi Surface approximation.We see that the direct Urca neutron decay rate is indeed strongly suppressed, but the direct Urca electron capture rate only shows a slight reduction below the threshold, and remains well above the modified Urca rates.This mismatch is what leads to a positive correction ∆µ in beta equilibrium.We will explain below why this is the case.
For SFHo (Fig. 8), the analysis is similar: neutron decay is heavily suppressed as expected in the dUrcaforbidden region (up to infinite density), but electron capture is much less suppressed.In the middle density range (3 to 5 n 0 ) where mUrca is dominant there is no need for a correction, since the mUrca rates balance at ∆µ = 0.However, at lower or higher densities the di- Urca rates calculated using the SFHo EoS at T = 3 MeV.Because ∆µ = 0 there is a large mismatch between the direct Urca rates for neutron decay and electron capture.Modified Urca (with neutron spectator (n) and proton spectator (p)) rates are calculated in the Fermi surface approximation and therefore match automatically.
rect Urca electron capture rate becomes large enough to dominate, so a positive ∆µ will be required to pull it down and establish equilibrium between neutron decay and electron capture.
In the next subsection we analyze the imbalance between electron capture and neutron decay rates in the dUrca-forbidden density range.This imbalance is the reason why a nonzero ∆µ is required in beta equilibrium.We can understand the difference in the rates, and their density dependence, by looking at which parts of the phase space dominate the rate integrals.This is largely determined by the Fermi-Dirac factors in the rate integrals, since the matrix element depends only weakly on the magnitudes of the momenta.

C. Direct Urca suppression factors
The density and temperature dependence of the direct Urca rates is dominated by the Fermi-Dirac factors.Below the dUrca threshold density, at zero temperature all direct Urca processes would be forbidden, but at nonzero temperature the Fermi surfaces are blurred, so there is some non-zero occupation of particle and hole states in regions of momentum space where the direct Urca process is kinematically allowed.The rate is governed by the Fermi-Dirac suppression factors for those momentum states.
At each density and temperature we search for the combination of momenta that is least suppressed, i.e. that maximizes the product of Fermi-Dirac factors in the rate integral while maintaining energy-momentum conservation.The magnitude of that product of Fermi-Dirac factors tells us how suppressed the whole process will be, at that density and temperature.
Below the direct Urca threshold density, considering particles near their Fermi surfaces, the neutron has a momentum larger than the sum of proton and electron momenta, even if the proton and electron are coaligned (see Fig. 1).In this regime, the direct Urca kinematics will become essentially one-dimensional, as this is how the electron and proton momenta can come closest to adding up to the large neutron momentum.We take the neutron momentum to be positive, so a negative momentum indicates motion in the direction opposite of the neutron.For momentum conservation to hold, the electron and proton will have to be away from their Fermi surfaces.In the assumption of one-dimensional kinematics, we determine the optimal momenta {k opt n , k opt p , k opt e , k opt ν } as follows.For neutron decay, we maximize f n (1 − f p )(1 − f e ) and for electron capture we maximize (1 − f n )f p f e .Energy and (one-dimensional) momentum conservation impose two constraints on the momentum, leaving two independent momenta over which to maximize.
The results of this maximization exercise are shown for the IUF EoS in Figs. 9 and 10, and for SFHo in Figs.11  and 12.
The left panels show how far from their Fermi surfaces the particles are in the least Fermi-Dirac-suppressed kinematic configuration.For each particle i we show which is the extra energy the particle with its optimal momentum has relative to its Fermi energy.The curves only exist in the dUrca-forbidden region, which for IUF ends slightly above 4 n 0 .(In the dUrca-allowed region all particles can be on their Fermi surfaces, so the curves would be trivially zero and are not shown).The right panels show the maximum value of the Fermi-Dirac factor, which gives the overall suppression of the process.Neutron decay Direct Urca neutron decay is suppressed because the neutrons at their Fermi surface have just enough energy to make a proton and electron near their Fermi surfaces (this is a consequence of the beta equilibrium condition (6)), but too much momentum (Fig. 1).The process can still proceed (with an exponential suppression factor) by exploiting the thermal blurring of the Fermi surfaces.Figs. 9 (IUF) and 11 (SFHo) show that the best option is to create a proton at energy γ p above its Fermi surface and an electron at energy γ e = −γ p which is below its Fermi surface.The co-linear proton and electron now have more momentum then when they were both on their Fermi surfaces because the proton's momentum rises rapidly with γ p because the proton is less relativistic, whereas the electron's momentum drops more slowly as γ e becomes more negative, because the electron is ultrarelativistic.The creation of the proton incurs no Fermi-Dirac suppression because states above the Fermi surface are mostly empty, but the creation of the electron is suppressed by a Fermi-Dirac factor of e −|γe|/T reflecting the scarcity of electron holes available to take such an electron.The net suppression of the rate, e −|γe|/T e −|γn|Θ(γn)/T , is shown in the right panels of Figs. 9 (IUF) and 11 (SFHo).For IUF we see the strongest suppression at around 2 n 0 , which explains the density dependence of the IUF neutron decay rate shown in Fig. 7.For SFHo, since the dUrca-forbidden region extends up to infinite density, and the momentum deficit remains large across the density range surveyed, we see stronger suppression that does not relent at the upper end of the density range, explaining the almost total suppression seen in Fig. 8.
We can understand the density dependence of γ p in terms of the one-dimensional model within which the maximization was performed.
We assume that, as seen in Figs. 9 (IUF) and 11 (SFHo), the neutron remains on its Fermi surface, and the neutrino takes negligible energy/momentum, since lack of momentum to build the final state is the main obstacle.Conservation of energy and momentum then tells us that Using the dispersion relations (8) we can solve for k opt p and k opt e and, after using that E F n = E F p + E F e (since we have assumed ∆µ = 0), we find where ∆k ≡ k F n −k F p −k F e is the momentum deficit (we plotted the surplus −∆k in Fig. 1).From this analysis, we learn that the density dependence of γ p , and therefore the rate, not only depends on the momentum deficit ∆k, but on the relative behavior of the neutron and proton Fermi momenta and their effective masses.
Although the momentum deficit ∆k in IUF monotonically shrinks with density, γ p shows a slight increase at low densities due to the fast drop of the effective proton mass m * p (see Fig. 3).This fast decrease counterintuitively leads E * F p to drop with density, while the real Fermi energy, which includes the nuclear mean field, U p , rises with density as expected.Closer to the threshold density, the momentum deficit dominates the behavior of γ p and the rate, so that γ p goes to zero at the threshold as ∆k approaches zero, leading to k opt p = k F p as expected.
For the SFHo EoS, the direct Urca momentum deficit is only varying weakly with density (see again Fig. 1).Although the momentum deficit is slowly falling, γ p continues to rise with density as shown in Fig. 11.This is due to the neutron Fermi momentum which rises fast enough that the denominator in Eq. (27) decreases by more than a factor of five in the studied density range while the momentum surplus stays nearly constant in comparison.

Electron capture
In the dUrca-forbidden density range, using the onedimensional kinematics described above, we find that the optimal kinematics for electron capture has a proton above its Fermi surface and an electron close to its Fermi surface combining to make a neutron slightly below its Fermi surface and a neutrino.The Fermi-Dirac suppression factor is e −γp/T e −|γe|Θ(γe)/T e −|γn|Θ(−γn)/T , reflecting the scarcity of protons and electrons above their Fermi surfaces, and of neutron holes below the neutron Fermi surface.
Figs. 10 and 12 show the corresponding energy excesses γ i and Fermi-Dirac suppression factors.In the right panels we see that in the dUrca-forbidden region, electron capture is somewhat suppressed but not nearly as suppressed as neutron decay.This is because, as we explain below, it is able to proceed using a proton that is much closer to its Fermi surface than is possible for neutron decay, and there is correspondingly less Fermi-Dirac suppression (compare the left panels of Figs. 9 vs. 10, and Figs.11 vs. 12).
The special feature of electron capture is that there is a very efficient way to exploit the thermal blurring of the Fermi surfaces.Given a momentum shortfall ∆k ≡ k F n − k F p − k F e , we can start with a proton whose momentum is less than ∆k above the Fermi surface.The rarity of finding such a proton leads to a Fermi-Dirac suppression factor of e −|γp|/T .This proton captures an electron near its Fermi surface with momentum parallel to the proton's.At this point their combined momentum is not enough to make a neutron on its Fermi surface, and there is excess energy.But we can use that excess energy to create, along with a neutron on its Fermi surface, a neutrino whose momentum partly cancels the neutron momentum, so the combined momentum of the proton and electron is enough to create that final state.
Because of the "help" from the neutrino, the proton does not need to be as far above its Fermi surface as the proton in neutron decay, so the electron capture rate is suppressed by a smaller Fermi-Dirac factor, The density dependence of the suppression factors (right panels of Fig. 10 for IUF and Fig. 12 for SFHo) explain the density dependence of the direct Urca electron capture rates shown in Figs.7 and Figs. 8. To understand the density dependence of γ p , we can perform a similar analysis as for neutron decay.We now assume neutron and electron to be on their Fermi surfaces, as shown in Figs. 10 (IUF) and 12 (SFHo), which is not as good as an assumption compared to the neutron decay analysis, but still helps us to gain insight into the behavior of the rates.Energy-momentum conservation again allows us to deduce that which leads, following the same procedure as in the neutron decay case, to For IUF at low densities, we can neglect the proton Fermi momentum compared to the effective mass.The behavior of γ p is then again dominated by the effective proton mass, whose rapid decrease overcomes the rising neutron Fermi momentum at low densities.This pushes the proton further away from its Fermi surface at low densities, before the momentum surplus dominates the behavior of γ p as the threshold is approached.As for neutron decay, ∆k = 0 at the threshold, therefore the rate is again dominated by particles on their respective Fermi surfaces.
For SFHo, the momentum surplus is becoming smaller from n 0 to 3 n 0 while the combination of the effective masses and Fermi momenta in (30) varies slowly with density.This allows the behavior of the momentum surplus ∆k to dominate the behavior of γ p at low densities, so both are increasing and therefore pushing the proton further away from its Fermi surface initially.At higher densities, SFHo is seemingly approaching asymptotically a direct Urca threshold.Both the momentum surplus and the Fermi momenta and effective masses in Eq. ( 30) are pushing the ideal proton momentum back closer to the Fermi surface.Overall, the behavior of the electron capture rate in SFHo can therefore largely be explained by the density dependence of the momentum surplus.

D. Non-relativistic Rate vs. Relativistic Rate
In Sec.III we emphasized that as the density rises above about 2n 0 relativistic corrections become important in the nucleon dispersion relations.In this section we illustrate the importance of relativistic corrections in the neutron decay rate.

Direct Urca neutron decay
Fig. 13 shows various approximations to the direct Urca neutron decay rate at T = 3 MeV (with ∆µ = 0).We show the rate calculated with fully relativistic dispersion relations, with the non-relativistic dispersion relation and with the "vacuum dispersion relation" used in [5], where m N = 940 MeV, and m eff,N is chosen such that For the non-relativistic curves, we use a corresponding non-relativistic approximation of the rescaled matrix element (B3), see Refs.[5,31], and the derivation in Appendix C of [60].We see that relativistic corrections make an enormous difference to the rate.The non-relativistic approximation is reasonably accurate at low density (where the nucleons are indeed non-relativistic) but overestimates the rate by up to eight orders of magnitude (at T = 3 MeV) between 2 n 0 and the direct Urca threshold at 4.1 n 0 .Due to the breakdown of the non-relativistic approximation, the direct Urca threshold condition is incorrectly already fulfilled below two times saturation density, which explains the steep increase of the non-relativistic rate around this density.For a detailed discussion of the density dependence of the relativistic rate, see Sec.VI C.
The thermal blurring of the Fermi energy, which is proportional to the temperature T , translates to a blurring in momentum space of order T /v F , where v F is the Fermi velocity.In the correct relativistic treatment, v F has an upper bound of 1, whereas for the non-relativistic dispersion relation, the Fermi velocity grows without a limit.This leads to a suppression of the non-relativistic rate at higher densities which partially cancels the effects of the earlier threshold.
The "vacuum dispersion relation" gives a rate that is one to eight orders of magnitude too large (at T = 3 MeV), and is less suppressed at higher densities since the corresponding Fermi velocity stays comparatively small in the plotted density range.

Modified Urca neutron decay
Fig. 14 shows the importance of using relativistic dispersion relations in calculating modified Urca.The rates are calculated for the IUF equation of state in the Fermi surface approximation at T = 3 MeV.The relativistic rate is about 1 to 2 orders of magnitude smaller than the non-relativistic rate.The modified Urca rates are not sensitive to the direct Urca threshold because of the spectator providing extra momenta.Much of the difference between the non-relativistic calculation and the relativistic calculation comes from the prefactors, as shown in Sec.V and Eqs.(C14), (C17), (C15) and (C18).The relativistic rates are suppressed by i m * i /E * i , where i is the index for each of the nucleons participating the interaction.Notice that the proton-spectator modified Urca rate is always less than the neutron-spectator rate because the proton Fermi surface, and its accompanying phase space, is smaller.

VII. CONCLUSIONS
We have investigated the conditions for beta equilibrium in nuclear matter in neutron stars, focusing on the temperature range where the material is cool enough that neutrinos escape (T ≲ 5 MeV) but warm enough so that nonzero temperature corrections to the Fermi surface approximation play an important role (T ≳ 1 MeV).
Previous work [5] found that a nonzero-temperature correction ∆µ to the traditional beta equilibrium condition (Eq.( 7)) was required to balance the rate of neutron decay against the rate of electron capture.We have improved on that calculation by using a consistent description of nuclear matter, based on two relativistic mean field models, IUF and SFHo.
We find that when using relativistic mean field models it is important to use the full relativistic dispersion relations of the nucleons.In these theories the effective masses drop quickly with density, so the neutrons become relativistic at densities of 2 to 3 n 0 .Using non-relativistic nucleon dispersion relations can make the modified Urca rates wrong by an order of magnitude and the direct Urca rates wrong by many orders of magnitude.
Our results for the nonzero-temperature correction ∆µ are shown in Figs. 4 and 5.We find that it rises with the temperature, and can be of order 10 to 20 MeV for temperatures in the 3 to 5 MeV range.The density dependence is quite different for the two EoSs that we studied, and we showed in detail how it depends on specific properties of the EoS.
We find that the nonzero-temperature correction plays an important role in the correct calculation of Urca rates.Using the naive (low-temperature) beta equilibrium con-dition µ n = µ p + µ e at T = 3 MeV would yield electron capture rates that are too large by an order of magnitude, and neutron decay rates that are too small by an order of magnitude (Fig. 6).This would significantly affect calculations of neutrino emissivity in the cooler regions of a neutron star merger, and therefore the estimated energy loss due to neutrinos.Currently used neutrino leakage schemes (e.g.Ref. [61] and references therein), which often treat the temperature range T ≲ 5 MeV as neutrino free streaming, need to be adapted to the corrected beta equilibrium.Additionally, the bulk viscosity of nuclear matter [62] depends on the rate of the Urca process which restores the system to beta equilibrium.The improved calculation of the Urca rates presented here will modify the temperatures and densities at which bulk viscosity reaches its maximum strength.Using the correct beta equilibrium condition also affects the equation of state: a recent study estimated its impact to be at the 5% level [63], and it would be interesting to evaluate the impact by performing a merger simulation using an EoS that incorporates the the finite-temperature correction described in this paper.[42].
There is no neutrino Fermi-Dirac factor because we assume the medium is neutrino-transparent, i.e., neutrinos escape the star.The spin-summed matrix element is [11] |M where G = G F cos θ c , G F = 1.166 × 10 −11 MeV −2 is the Fermi constant and θ c = 13.04 • is the Cabbibo angle.As they originate from spin summations (see Appendix B of [9]), the 4-vector dot products in the matrix element (B2) are k µ = (E * , k).
It is convenient to define the rescaled dimensionless matrix element In the nonrelativistic limit, since g A ≈ 1, M ≈ (1 + 3g 2 A ) ∼ 4 [11,20,31,34,65,66].The neutron decay rate can now be written The 12-dimensional integral can be reduced to a 5-dimensional integral as follows.Integrating over the 3-momentum conservation delta functions reduces the integral to 9 dimensions (compare (E.1) in Ref. [60]) The remaining delta function imposes energy conservation in the creation of the neutrino: Each momentum integral can be written in polar co-ordinates as d 3 k = k 2 dkdzdϕ where z = cos θ.Setting up the following coordinate system (see Appendix E in [60]) allows us to integrate over z n and ϕ n yielding a factor of 4π and over ϕ p yielding a factor of 2π, which eliminates three angular integrals, so that (compare (E.5) in [60]) where Note that for simplicity we label the electron azimuthal angle as ϕ (rather than ϕ e ).The factor of Θ(E ν ) restricts the integral to the region of momentum space where the neutrino energy E ν (k n , k p .k e ) is positive, which is a requirement for the emission of a neutrino.This condition leads to the upper limits on the proton and electron momenta.If we perform the integrals in the order shown in (B10) then the electron momentum integral is the inner integral, so it is performed for known values of k n and k p , so the constraint E ν > 0 corresponds to E e < E n − E p .Similarly, the k p integral is performed for a known value of k n , so its range is constrained by requiring that there is enough energy to create an electron (of unknown momentum) and a neutrino, E p < E n − m e .This leads to upper limits on the proton and electron integral, In the delta function in Eq.(B11), where Since g(ϕ) depends on ϕ only via cos ϕ there will be either zero or two solutions to g(ϕ) = 0, so where M ϕ0 is the dimensionless rescaled matrix element (B3) evaluated at ϕ 0 , which can be either of the two solutions of g(ϕ) = 0, It does not matter which solution we use for ϕ 0 because g is a function of cos ϕ and M depends only on cos ϕ and sin 2 ϕ, so the integrand has the same value for both the solutions.The theta function Θ(S − |E 2 ν − R|) imposes the condition that there are two solutions (rather than none), by limiting the integral to the domain where −1 < cos ϕ 0 < 1.
We now use (B14) and (B18) to evaluate the integrand in (B17).Firstly, the Jacobian of the delta function is Using (B19) in (B17), Secondly, substituting (B18) in to (B2) gives the matrix element where Limits of angular integration To speed up the numerical evaluation of (B20) we implement the theta function as limits on the range of integration over z p and z e .The condition S > |E 2 ν − R| can be written (using (B15), (B16)) as where The inequality (B25) is obeyed for z − e < z e < z + e where Note that if the roots are real then they are always within the physical range z e ∈ [−1, 1].We can therefore put bounds on z p by requiring that (B29) has real roots, ) This means that z − p < z p < z + p , where In this case, however, these bounds are not necessarily within the physical range z p ∈ [−1, 1], so the true bounds on the z p integral are We can now write the angular integral as Using this in (B10) we obtain The second line corresponds to the I integral (B11), (B33).It is natural to group a factor of E ν with M ϕ0 to cancel the factor of E ν in the denominator (B21) which can cause numerical problems at the edge of the kinematically allowed momentum range where E ν → 0. The neutron decay rate can therefore be computed as a 5-dimensional momentum integral (B34), obtaining the integration ranges from (B12), (B13), (B32), and (B29), the matrix element from (B21), and the Jacobian (square root denominator) from (B15), (B16).

Appendix C: Modified Urca neutron decay rate
The matrix element is (4.16) in [31,60] where f ≈ 1 is the N-π coupling and s = 1/2 for the identical particles.The conventional way of doing the integral is to divide the integral into an energy integral and an angular integral (termed "phase space decomposition" [35]) We use relativistic dispersion relations for nucleons where U is the mean field contribution to the energy.We define E * ≡ √ k 2 + m * 2 , then dE * = kdk/E * .We use ultra-relativistic dispersion relations for electron and neutrino, then dE = dk (the electron mass m e = 0.511 MeV is negligible compared to its momentum).Therefore, we can convert the momentum integral to an energy integral, and the rate integral becomes Γ mU,nd(n) = 42G 2 g 2 A f 4 (2π) 14 m 4 π dΩ n dΩ p dΩ e dΩ ν dΩ N1 dΩ N2 × δ (3) Notice that it is most common to set ⃗ k ν = 0 in the momentum conserving delta function but keep E ν in the energy delta function.
In the Fermi surface approximation, we set all momenta to Fermi momenta and we will have E e = k e = k F e , k ν = E ν .Now, the rate integral becomes × dΩ n dΩ p dΩ e dΩ ν dΩ N1 dΩ N2 δ (3) For the energy integral, we do a change of variable, where the approximation is valid because µ * ≫ T .For neutrino, µ = 0 and m = 0, so the lower bound is 0.Then, the energy integral, which we denote as I, becomes =T 7 dx n dx p dx e dx ν dx N1 dx N2 For the angular integral, we can look up [25], which calculated the n-dimensional angular integral for n=3,4,5, and obtain where Γ mU,ec(p) (ξ) = Γ mU,nd(p) (−ξ) , where (C19) FIG.2.Density dependence of the neutron's (Dirac) effective mass and Fermi momentum for the IUF and SFHo EoSs, showing that neutrons at the Fermi surface become relativistic at densities above 2 to 3 n0.

14 density n B /n 0 Γ [MeV 4 ]
FIG. 8.Urca rates calculated using the SFHo EoS at T = 3 MeV.Because ∆µ = 0 there is a large mismatch between the direct Urca rates for neutron decay and electron capture.Modified Urca (with neutron spectator (n) and proton spectator (p)) rates are calculated in the Fermi surface approximation and therefore match automatically.

FIG. 9 . 1 density n B /n 0
FIG.9.The optimal kinematics for neutron decay for the IUF EoS.Left panel: the least suppressed kinematic arrangement, showing the energy distance γ of each particle from its Fermi surface.Right panel: the Fermi-Dirac suppression factor, e −|γe|/T e −|γn|Θ(γn)/T which is dominated by the difficulty of finding an electron hole at energy γe below its Fermi surface.

FIG. 11 .FIG. 12 .
FIG.11.The optimal kinematics for neutron decay at T = 3 MeV for SFHo, obtained by maximizing the Fermi-Dirac products.The suppression factor, e −|γe|/T e −|γn|Θ(γn)/T is dominated by the difficulty of finding an electron hole below its Fermi surface.

12 density n B /n 0 Γ [MeV 4 ]direct
FIG.13.Direct Urca neutron decay rate calculated using relativistic, non-relativistic and the vacuum dispersion relations at T = 3 MeV for IUF.