The essential structure of modeling nuclear matter using RBHF theory is outlined in this section. A more detailed explanation of the approach can be found in Poschenrieder and Weigel [24]; Weber [14], with finite temperature extensions given in our previous work [30, 34]. Nuclear matter at supranuclear densities can be described as a complex, many-body system whose dynamics are governed by the Lagrangian density: L = L N + ∑ M L M + L M N , where L N denotes the Lagrangian of non-interacting nucleons, L M is the Lagrangian density of different free meson fields, and L M N describes the interaction between nucleons and mesons. In the relativistic framework, nucleons are treated as effective Dirac particles, which are described by the relativistic Dirac equation. The equations of motion for the various particle fields within the many-body system are derived from the Euler-Lagrange equation and solved using the Martin-Schwinger hierarchy of coupled Green’s functions [35].

The formal structure of the RBHF approach is to solve a system of highly nonlinear, coupled equations, which include the Dyson equation for the two-body Green’s function G 1 , a Bethe-Salpeter type integral equation for the scattering matrix T , and the equation for the self-energy Σ . In this formalism, the Dyson equation is used to determine the two-body Green’s function G 1 for all nucleons as: p − m − Σ p G 1 p = 1 , where p = γ μ p μ = γ μ g μ ν p ν = γ 0 p 0 − ∑ i = 1 3 γ i p i and γ 0 and γ i are the Dirac matrices. Following the determination of the G 1 , the in-medium scattering matrix T is determined with Bethe-Salpeter type integral equation: T P ; p , p ′ = V p − p ′ + ∫ d 4 p ′ ′ 2 π 4 V p − p ′ ′ Λ P 2 + p ′ ′ , P 2 − k ′ ′ T P ; p ′ ′ , p ′ , (1) where Λ is the intermediate nucleon-nucleon propagator and V represents the repeated sums of two-particle interactions given by a one-boson-exchange (OBE) potential, which describes the interaction among two nucleons in terms of the exchange of scalar, pseudo-scalar, and vector mesons. In this work, we use the so-called Bonn B potential [36] to describe the OBE interaction, which employs a pseudoscalar type of pion-nucleon coupling. This potential is widely used in high-density nuclear matter studies due to its reliability and numerical stability, particularly in RBHF theory. While newer relativistic OBE potentials have been proposed (e.g., in [37]), the Bonn-B potential has been extensively tested in dense nuclear environments, making it a robust and reliable choice for modeling the equation of state under extreme conditions.

The final coupled equation in the formal scheme is for the self-energy Σ , written in momentum space as: Σ p = i 2 π 4 ∫ d 4 p ′ p − p ′ 2 T p + p ′ p − p ′ 2 − p − p ′ 2 T p + p ′ p ′ − p 2 G 1 p ′ . The self-consistent calculations are carried out using a complete basis of particles ( Φ λ ) and antiparticles ( θ λ ) to decouple the integral equations and make the two-body propagator Λ diagonal, where λ = ± 1 / 2 are the helicity eigenvalues.

The particle propagator Λ in Equation 1 takes the form of the Brueckner propagator, which is defined as the following at finite temperatures: Λ p ⃗ , p ⃗ ′ ; P 0 = 2 π f | p ⃗ | f | p ⃗ ′ | P 0 − ω 1 | p ⃗ | − ω 1 | p ⃗ ′ | , where f signifies the Fermi-Dirac distribution functions, given by: f 1 p ⃗ = 1 e β ω 1 p ⃗ − μ + 1 , f 2 p ⃗ = 1 e β − ω 2 p ⃗ + μ + 1 , where “1” indicates the positive energy states and “2” indicates the thermally-excited negative energy states. We recall here that at finite temperatures, the behavior of nuclear matter undergoes important modification, attributed to thermal baryonic excitations surpassing the Fermi surface. As T → 0 , the Fermi-Dirac distribution for positive energy states becomes: f 1 p ⃗ → Θ μ − ω 1 p ⃗ , and for negative energy states, f 2 → 0 .

An elegant technique used to make the many-body equations numerically tractable and to calculate the key quantities of many-body systems is to utilize the spectral representation of the G 1 function [24]. G 1 can then be defined in Fourier space at finite temperatures as [14, 38]: G 1 p 0 , p ⃗ = ∫ d ω Ξ ω , p ⃗ ω − p 0 − μ 1 + i η − 2 i π sign p 0 − μ Ξ p 0 − μ , p ⃗ e β | p 0 − μ | + 1 where Ξ represents the spectral function which is dependent on the single-particle energy ω , μ is the nucleon chemical potential, and η is used to circumvent a singularity occurring as integrals are carried out in the complex plane. The temperature inclusion arises through β = 1 / k B T , where k B is the Stefan-Boltzmann constant and T is the temperature.

Once a self-consistent solution to the coupled system of equations is found, the self-energy Σ and spectral function A are used to determine the EOS. The number density ρ of the system follows from: ρ = 4 2 π 3 ∫ d 3 p Ξ 0 ω , p ⃗ f p , where Ξ 0 is the time-like component of the spectral function and f denotes the Fermi-Dirac distributions. The pressure of the system at finite temperatures is determined from the free energy per nucleon, denoted as F , which is defined as: F ρ , T = U ρ , T − T S ρ , T , where U is the internal energy, T is the temperature, and S is the entropy ter Haar and Malfliet [39]. In this approach, T is held constant and S is determined by: S ρ , T = − 1 2 π 3 ρ ∫ d 3 p 1 − f ln 1 − f + f ⁡ ln f . (2) Both particles and antiparticles contribute to S , but the antiparticle contribution is very small ( f 2 ≪ f 1 ) . Therefore, Equation 2 can be approximated with only the particle contribution, f = f 1 . Once the entropy S , and subsequently the free energy F , are calculated, the pressure is derived as: P ρ , T = ρ 2 ∂ F ρ , T ∂ ρ .

Using the outlined theory, two models for the EOS of neutron star matter are constructed at temperatures T = 10 and 50 MeV, shown visually in Figure 1. As shown in previous work Farrell and Weber [30], the maximum mass of each EOS for non-rotating and uniformly rotating stellar sequences at their mass-shedding limit is over 2 M ⊙ , as required by observational constraints. These models will be used as input to determine equilibrium models for differentially rotating objects, of which the formalism is discussed below.

The theoretical framework for modeling differential rotation in neutron stars described in this work follows from the framework laid out by Komatsu et al. [40], which was then modified in Cook et al. [8] (referred to as CST throughout the text). The equations shown in this section directly follow the modifications introduced in CST.

To model differentially rotating neutron stars, we begin with the definition of the line element [8]: d s 2 = − e γ − ρ d t 2 + e 2 α d r 2 + r 2 d θ + e γ − ρ r 2 ⁡ sin 2 ⁡ θ d ϕ − ω d t 2 , (3) where the metric potentials ρ , γ , α , and ω are dependent on both the radial r and polar θ coordinates. For both uniform and differential rotation, Equation 3 models neutron stars as stationary, axisymmetric configurations of a (self-gravitating) perfect fluid [8]. Under the assumption of neutron star matter as a perfect fluid, sources of non-isotropic stresses such as magnetic fields or heat transport are ignored [41]. This assumption also allows neutron star matter to be described by the energy-momentum (or stress-energy) tensor given by: T κ σ = ϵ + P u κ u σ + g κ σ ⁡ P , where u is the fluid’s 4-velocity, κ and σ are indices ranging from 0 to 3, and ϵ and P are given by the underlying EOS. Equilibrium models for neutron stars must obey the equation of hydrostatic equilibrium as Einstein’s field equation, given as: R κ σ − 1 2 R g κ σ = 8 π T κ σ , where R κ σ is the Ricci tensor, R is the curvature scalar, and g κ σ is the metric tensor.

Equilibrium models for neutron stars must obey the equation of hydrostatic equilibrium, which has the form: h P − h p = ∫ P p P d P ϵ + P = ln u t − ln ⁡ u p t − ∫ Ω c Ω F Ω d Ω , where h ( P ) is enthalpy as a function of the pressure, u t is the time-like component of the 4-velocity u μ , and the subscripts p and c denote the variable’s value at the pole or center, respectively. The integrand of the final integral term, F ( Ω ) , is the function that defines the rotation law of the matter in the case of differential rotation. Following CST, we define F ( Ω ) as a linear rotation law: F Ω ̄ = A 2 Ω ̄ c − Ω ̄ , (4) where A is a parameter that dictates the degree of differential rotation within the star. In the case of uniform rotation, Equation 4 disappears as the value for the frequency at the center of the star is constant throughout (i.e., Ω c = Ω ). It is important to note that the choice of rotation law directly impacts the maximum mass for a given EOS. While many implementations of CST use the linear rotation law in Equation 4 (see, for example, [42, 43]), other studies of differential rotation have explored the impact of using either modified versions of the linear law or non-linear rotation laws, as shown in Galeazzi et al. [10]; Hanauske et al. [11]; Zhou et al. [44].

Using the linear rotation law in Equation 4, the equation of hydrostatic equilibrium can be integrated to give: h P − h p = 1 2 γ p + ρ p − γ − ρ − ln 1 − v 2 + A ̂ 2 Ω ̂ − Ω ̂ c 2 , (5) where A ̂ is the rotation parameter scaled as A ̂ = A / r ̄ e . The matrix of angular frequency Ω ̂ can be derived using the following equation: Ω ̂ c − Ω ̂ = 1 A ̂ 2 Ω ̂ − ω ̂ s 2 1 − μ 2 e − 2 ρ 1 − s 2 − Ω ̂ − ω ̂ 2 s 2 1 − μ 2 e − 2 ρ . As the rotation parameter A appears throughout the numerical scheme scaled and inverted, we follow the lead of previous work which parameterized calculated sequences of differentially rotating stars by values of A ̂ − 1 = 0.3, 0.5, 0.7, and 1.0 [8, 10]. Uniform rotation is obtained in the limit A ̂ − 1 → 0 , and an upper bound of the scaled rotation parameter in this study is A ̂ − 1 = 1.0 .

In the numerical scheme, the metric potentials ( ρ , γ , ω ̂ , and α ) are used to determine a value for the radius at the equator, r e of the star: r ̄ e 2 = 2 h P ̄ ϵ ̄ c − h p γ ̂ + p + ρ ̂ p − γ ̂ m − ρ ̂ m , (6) which is equivalent to Equation 5 evaluated at the location of the maximum (denoted by subscript m ) density of the star. In CST, this location is assumed to be at the star’s center. The equation for updating r e changes when the maximum density within the star is not in the center, which is the case in very deformed configurations. We modify the original CST algorithm by instead fixing the maximum interior density, thus redefining Equation 6 as: r e 2 = 2 h P m − h p + ln 1 − v m 2 − A ̂ 2 Ω ̂ m − Ω ̂ c 2 γ ̂ p + ρ ̂ p − γ ̂ m − ρ ̂ m , where the subscript m denotes the quantity’s value at the coordinates of the maximum density.

In this section, stellar sequences are constructed over a range of constant central densities for the two EOS models at temperatures of 10 and 50 MeV. These sequences are calculated with a fixed value for the rotation parameter set the degree of differential rotation; as mentioned in Section 2.2, sequences are parameterized by fixing A ̂ − 1 to be 0.3, 0.5, 0.7, or 1.0. When assuming the maximum density of the star is no longer in the center, more extreme configurations with lower values of r r a t i o can be calculated. Sequences presented in this section are also calculated over a range of r r a t i o values to demonstrate the parameter’s impact on the maximum mass for a given EOS.

For each EOS, sequences over a range of constant central densities are computed for varying values of r r a t i o . These sequences are shown for the EOS model with temperature T = 10 MeV in Figure 2 and for T = 50 MeV in Figure 3. In both figures, a sequence of stars rotating uniformly at their mass-shedding limit is plotted for comparison. As echoed in previous work, more extreme differential rotation ( A ̂ − 1 → 1.0 ) paired with more extreme structural deformation ( r r a t i o → 0.0 ) results in higher masses at lower central densities, which then taper off to the Kepler limit at higher central densities. This trend is also encountered in Morrison et al. [43], who make similar modifications to the CST algorithm.

Rotating neutron stars formed from a core-collapse suprenova (CCSN) or binary stellar mergers may experience nonaxisymmetric instabilities that directly impact their rotation rates and overall stability. Previous studies (see [45, 46]) in Newtonian gravity have shown rotational instabilities arise from non-radial toroidal modes, i.e., e i m ϕ ( m = ± 1,2 , … ) , which result in the stability parameter, or the ratio of rotational ( T ) to gravitational ( W ) energy β = T / | W | , exceeding some critical value β c . These rotational instabilities are likely to impact a star’s gravitational radiation signal, making the study of such instabilities an important topic in the wake of new gravitational wave detectors. In this section, we will focus on determining if calculated stellar models are subject to the so-called bar-mode instability, where m = ± 2 , which is expected to be the fastest-growing mode and the subject of many instability studies for both uniformly and differentially rotating neutron stars.

Two mechanisms cause rotating stars to be unstable to bar-mode deformation: secular and dynamical instabilities. In Newtonian theory, uniformly rotating incompressible neutron stars become secularly unstable to bar-mode deformation at a critical value of β c ≥ 0.14 with similar findings in studies of post-Newtonian theories [47]. In general, the secular instability grows only in the presence of a small dissipative mechanism like viscosity or gravitational radiation at lower rotation rates [48]. The secular instability usually has a longer growth time when compared to the dynamic timescale of the system. A similar critical value of β c has been observed in numerical studies of relativistic stars but is also dependent on the compaction (M/R) of the star and the dissipative mechanism. For example, viscosity-driven secular instability has been shown to occur at β c > 0.14 in more compact configurations with higher rotation rates, but gravitational radiation-driven instabilities occur at β c < 0.14 at lower rotation rates. A more in-depth review of this topic is given by Paschalidis and Stergioulas [49].

The dynamical bar-mode instability occurs independent of any dissipative mechanism and with a growth rate determined by the dynamical timescale of the system, which is generally faster than the timescale of growth for secular instabilities. Therefore, numerical simulations of hydrodynamical equations are necessary to determine the onset threshold of the dynamical bar-mode instability. Many simulations have been carried out in Newtonian theory, the consensus of which gives the critical value β c ≥ 0.27 [4]. Simulations of the dynamical bar-mode instability in general relativity are less common, as solving the nonlinear hydrodynamical equations in full relativity is more complex. However, there have been reliable studies carried out by Shibata and Uryu (2000); Saijo et al. [50] for uniformly rotating stars and by Bodenheimer and Ostriker [51]; Shibata et al. [52, 53]; Camarda et al. [54]; Di Giovanni et al. [55] for differentially rotating stars.

The relativistic simulation of differentially rotating stars carried out by Shibata et al. [52] finds that the critical value of the stability parameter for the dynamic bar-mode instability is β c ≈ 0.24 − 0.25 , slightly lower than the Newtonian limit. These simulations were carried out using the same linear rotation law as in this chapter (and the proceeding one) for similar values of r r a t i o . Therefore, we adhere to their specified threshold of β to determine whether differentially rotating stars are stable or unstable to the dynamical bar-mode deformation. Specifically, we use the upper limit for the critical threshold, so stars deemed unstable in this section will have a stability parameter β c ≥ 0.25 .

A visual representation of stable and unstable models using this criterion at T = 50 MeV for two degrees of differential rotation is shown in Figure 4. The two panels in each row show the same information, mass vs. central density for a range of r r a t i o values, where the figures in the right column categorize each stellar model as stable with a blue square or unstable with a red dot based on the calculated stability parameter for each stellar model. The top row presents a lesser degree of differential rotation ( A ̂ − 1 = 0.3 ) , where 81.6% of calculated stellar models are considered stable against dynamical bar-mode deformation. In contrast, the bottom row has the highest degree of differential rotation ( A ̂ − 1 = 1.0 ) , and the increase in differential rotation decreases the percentage of stable models to 67.2%. For both cases of differential rotation, the majority of the unstable configurations have small values for r r a t i o , a trend that will be echoed in the next section.

Table 1 gives the percentage of unstable models and the average r r a t i o of those unstable models for four degrees of differential rotation for both temperatures of 10 and 50 MeV. The general trend shows the percentage of unstable models decreases as temperature increases. For both temperatures, r ̄ r a t i o generally tends to fall between 0.4 and 0.5. Stars with r r a t i o ≤ 0.5 are extremely deformed and thus likely vulnerable to dynamical instabilities. These findings imply temperature plays less of a role when compared to the star’s deformation on its stability against dynamical bar-mode excitation.

In this section, we examine density and frequency maps of individual stellar models. For each temperature, two models are computed: the first with a lesser degree of differential rotation ( A ̂ − 1 = 0.3) and the second with a larger degree of differential rotation ( A ̂ − 1 = 1.0). All four models are computed with an r r a t i o = 0.4, a significant degree of structural deformation. Information on the stars’ masses, radii, and stability parameters are given in Table 2.

The density and frequency maps for T = 10 MeV are shown in Figure 5, where the top row shows results for A ̂ − 1 = 0.3 and the bottom for A ̂ − 1 = 1.0 . The lesser degree of differential rotation is shown very clearly in the frequency map, where the difference of frequency values in the star only spans 320 Hz. In contrast, the star with a higher degree of differential rotation sees a span of 5,600 Hz from the highest to lowest frequency values. The density maps show an interesting depiction of the overall structural deformation. For A ̂ − 1 = 0.3 , the star takes an ellipsoid-like shape, but for A ̂ − 1 = 1.0 , the star instead takes a quasi-toroidal shape. While each star was initialized with the same central density of 400 MeV/fm³, the star with A ̂ − 1 = 1.0 experiences much higher densities outside of the center of the star. A similar story is seen for T = 50 MeV in Figure 6 but with even higher densities observed at the higher temperature.

As shown in Section 3.2, the stability of the star depends not only on the deformation characterized by r r a t i o but also on the degree of differential rotation. In Table 2, all four models have the same value of r r a t i o . For both temperatures, the stars computed with A ̂ − 1 = 0.3 have stability parameters T / | W | well below the critical limit of 0.25. However, the stars with A ̂ − 1 = 1.0 are well above the critical limit for β , where the stability parameter is 0.64 for T = 10 MeV and 0.75 for T = 50 MeV. These stars are well beyond the threshold for the dynamical bar-mode instability and thus likely unphysical.

The timescale over which the dynamical bar-mode instability develops, also known as the dynamical timescale, is proportional to R 3 / 2 M 1 / 2 . While this is not a definite indication of the growth time of the instability, which would instead require a full simulation in both time and space (see Ref. [52] for a good example), for extremely unstable configurations the dynamical timescale can provide some idea of how long these stars may exist with the bar-mode instability excited before collapse. For the unstable configuration at T = 10 MeV, the dynamical timescale is ∝ 35.8 s, and for the unstable configuration at T = 50 MeV, the dynamical timescale is ∝ 34.9 s.

The preceding sections present results dependent on the the underlying theoretical frameworks discussed in Section 2, which employ important approximations, discussed below.

The finite temperature EOS models at T = 10 and 50 MeV are derived from calculations utilizing the RBHF approximation, described in Section 2.1. While this method offers significant improvements over other approaches to modeling relativistic numerical matter, it still involves key approximations [36]. For instance, nucleon-nucleon interactions are modeled using OBE potentials, with the Bonn-B potential selected in this paper for stability at high densities. Although RBHF theory does not require adjustable parameters, which is an advantage, it omits higher-order quantum corrections [56, 57]. One of the most important higher-order effects missing in RBHF is the inclusion of three-nucleon forces, which are known to significantly influence high-density matter [58]. Other potentially influential effects excluded in the RBHF approximation include higher-order relativistic many-body corrections, ring diagrams, vertex corrections, self-energy insertions, and medium polarization effects; a more complete discussion of these effects can be found in [14, 59, 60].

Variations in the choice of OBE potential and omission of quantum corrections, such as three-body forces, can introduce errors in the RBHF EOS of up to 10%, as demonstrated in Brockmann and Machleidt [36]. Using chiral effective field theory, Hebeler and Schwenk [61] and Tews et al. [62] concluded that including three-body forces affects neutron star radius and mass predictions by 5%–10%.

In addition, the EOS models include finite temperatures (10 MeV and 50 MeV), which are essential for modeling neutron stars formed in extreme events like supernovae or binary neutron star mergers Steiner et al. [63]. However, as mentioned before, the neutron star crust is assumed to be at zero temperature, while the core is modeled at finite temperatures. This approximation may introduce inconsistencies in the determination of properties like neutron star mass and radius, especially in higher temperature regimes. In particular, the radius increases due to thermal pressure in the crust, potentially making the star slightly larger.

When modeling differential rotation in neutron stars, the rotation profile, given in Equation 4, is parameterized by A ̂ − 1 , representing the degree of differential rotation. This profile, while consistent with the original CST algorithm discussed in Section 2.2, is chosen for mathematical convenience rather than based on astrophysical observations, potentially limiting its accuracy for real neutron stars. Uncertainties in the choice of the rotation profile may cause errors of 5%–15% in star mass and radius predictions; for more on uncertainties related to differential rotation profiles and their impact on neutron star stability and mass limits, see Baumgarte and Shapiro [64]. Additionally, the assumption of axial symmetry is a simplification. In reality, neutron stars may exhibit more complex geometries, especially under differential rotation and dynamical instabilities.

For a comprehensive review discussing EOS uncertainties, neutron star mass-radius relationships, and the challenges of matching theoretical models with astrophysical data, see Lattimer and Prakash [65].