The Role of Strong Gravity and the Nuclear Equation of State on Neutron-star Common-envelope Accretion

Even for the highest spinning pulsars observed to date, such NSs can be considered as slowly rotating objects, meaning that rotation can thus be treated as a small perturbation to the Tolman–Oppenheimer–Volkoff (TOV) solution for nonrotating NSs (Oppenheimer & Volkoff 1939; Tolman 1939). Here, ≡ Ω/Ω_k is a dimensionless spin parameter, where Ω is the angular spin frequency of the star, and Ω_k is the Keplerian angular spin frequency ${{\rm{\Omega }}}_{{\rm{k}}}=\sqrt{{{GM}}_{\mathrm{TOV}}/{R}_{\mathrm{TOV}}^{3}}$, with M_TOV and R_TOV the mass and radius of our NS if it were not rotating. We solve for the structure of slowly rotating NSs to second-order in ≪ 1 using the Hartle–Thorne approximation (Hartle 1967; Hartle & Thorne 1968) with the same set of 46 hadronic EoSs from Silva et al. (2020, Appendix A). This set of EoSs is simultaneously consistent with the LIGO-Virgo observations of GW170817 (Abbott et al. 2017) and the NICER observation of PSR J0030 + 0451 (Miller et al. 2019a; Riley et al. 2019).

For a given EoS, and a chosen value of the central density and spin frequency Ω, the second-order solution to the Einstein equations in the Hartle–Thorne approximation allows us to calculate macroscopic properties of the star (e.g., Berti et al. 2005). These properties include the spin-corrected mass M_NS = M_TOV + ² δ M, the spin-corrected equatorial radius R_NS = R_TOV + ² δ R, the leading-order-in-spin moment of inertia I_NS and the spin-corrected dimensionless gravitochemical potential Φ_NS (Alécian & Morsink 2004).

In the context of accreting NSs, the gravitochemical potential can be interpreted as a susceptibility to changes in baryon number or the fraction of baryon mass that gets converted into gravitational mass. In the nonrotating limit, Φ_NS simplifies to

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{{\rm{\Omega }}\to 0}{{\rm{\Phi }}}_{\mathrm{NS}}={{\rm{\Phi }}}_{\mathrm{TOV}}=\sqrt{1-2{{ \mathcal C }}_{\mathrm{TOV}}}=\sqrt{1-\displaystyle \frac{2{{GM}}_{\mathrm{TOV}}}{{c}^{2}{R}_{\mathrm{TOV}}}},\end{eqnarray} \tag{ 1 }$$

where ${{ \mathcal C }}_{\mathrm{TOV}}$ is the compactness of a given nonrotating NS. We will use Φ_NS for our calculations, where we elaborate in Appendix B on how this is calculated with our EoS catalog. Equation (1) provides a fast approximation for population synthesis or as a subgrid prescription for global hydrodynamic simulations. We later compare in Section 3 how well this approximation compares to using Φ_NS, with a more detailed quantification shown in Appendix D.

We plot in Figure 1 the gravitochemical potential Φ_TOV versus the gravitational mass M_TOV for nonrotating NSs, as predicted from our EoS catalog. Each curve represents Φ_TOV for a different EoS, with different colors corresponding to different dimensionless NS binding energy $| {{ \mathcal B }}_{\mathrm{TOV}}| /({M}_{\mathrm{TOV}}{c}^{2})$. Observe that the gravitochemical potential decreases as the gravitational mass increases, and also it decreases as the NS compactness increases. Since the gravitochemical potential is inversely related to the binding energy, this figure tells us that binding energy conversion is enhanced for higher mass or higher compactness NSs.

Zoom In Zoom Out Reset image size

As the NS mass and spin increase during the CE inspiral, we are then able to track the temporal evolution of all of NS macroscopic quantities. For example, as the NS accretes, its gravitational mass responds to the baryon mass accretion rate ${\dot{M}}_{{\rm{b}}}$ as well as the angular momentum that the accreted mass carries. The resulting NS gravitational-mass accretion rate is

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{NS}}=\displaystyle \frac{\partial {M}_{\mathrm{NS}}}{\partial {M}_{{\rm{b}}}}{\dot{M}}_{{\rm{b}}}+\displaystyle \frac{\partial {M}_{\mathrm{NS}}}{\partial {J}_{\mathrm{NS}}}{\dot{J}}_{\mathrm{NS}}={{\rm{\Phi }}}_{\mathrm{NS}}{\dot{M}}_{{\rm{b}}}+\displaystyle \frac{{\rm{\Omega }}}{{c}^{2}}{\dot{J}}_{\mathrm{NS}},\end{eqnarray} \tag{ 2 }$$

where c is the speed of light, and J_NS = I_NSΩ is the NS spin angular momentum. Similarly, as the NS accretes, its spin angular momentum will also change, as given by

$$\begin{eqnarray}&&{\dot{J}}_{\mathrm{NS}}={\dot{I}}_{\mathrm{NS}}{\rm{\Omega }}+{I}_{\mathrm{NS}}\dot{{\rm{\Omega }}}=\displaystyle \frac{{{dI}}_{\mathrm{NS}}}{{{dM}}_{\mathrm{NS}}}{\dot{M}}_{\mathrm{NS}}{\rm{\Omega }}+{I}_{\mathrm{NS}}\dot{{\rm{\Omega }}}.\end{eqnarray} \tag{ 3 }$$

With this at hand, we can now solve for the temporal evolution of the angular frequency and the gravitational mass. We assume the NS accretes from a Keplerian accretion disk, where matter captured within the NS accretion radius carries angular momentum and spirals several orders of magnitude down to the scale of several NS radii. Approximating the total torque as the accretion torque, ${\dot{J}}_{\mathrm{NS}}\approx {\dot{M}}_{\mathrm{NS}}\sqrt{{{GM}}_{\mathrm{NS}}{R}_{\mathrm{NS}}}$ (e.g., Brown et al. 2000), in Equations (2) and (3), and for now ignoring other external torques on the NS, we then find

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{NS}}=\displaystyle \frac{{{\rm{\Phi }}}_{\mathrm{NS}}{\dot{M}}_{{\rm{b}}}}{1-\sqrt{{{GM}}_{\mathrm{NS}}{R}_{\mathrm{NS}}}\,({\rm{\Omega }}/{c}^{2})},\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&\dot{{\rm{\Omega }}}=\displaystyle \frac{{{\rm{\Phi }}}_{\mathrm{NS}}{\dot{M}}_{{\rm{b}}}}{{I}_{\mathrm{NS}}}\displaystyle \frac{\sqrt{{{GM}}_{\mathrm{NS}}{R}_{\mathrm{NS}}}-{\rm{\Omega }}\,({{dI}}_{\mathrm{NS}}/{{dM}}_{\mathrm{NS}})}{1-\sqrt{{{GM}}_{\mathrm{NS}}{R}_{\mathrm{NS}}}\,({\rm{\Omega }}/{c}^{2})}.\end{eqnarray} \tag{ 5 }$$

The evolution Equations (4) and (5) are generic for any slowly rotating NS accreting from a Keplerian disk. In general, however, the accretion and NS's angular-momentum evolution may be more complex. Such complications can arise if the NS's magnetic field pressure is comparable to the pressure of the radiation and accreting plasma, or if there is feedback from the accretion itself (e.g., Grichener & Soker 2019; Soker et al. 2019; López-Cámara et al. 2020). For a given EoS, we can then find the right-hand sides of the above equations as a function of M_NS and Ω, which leads to a closed system of ordinary differential equations, once ${\dot{M}}_{{\rm{b}}}$ is prescribed. In the CE inspiral context, the baryon mass accretion rate depends on the primary star's envelope structure, which we discuss in the following subsection.

We evolve single massive stars with the MESA (v12778) stellar-evolution code (Paxton et al. 2010, 2013, 2015, 2018, 2019) to obtain their interior structure. We consider a total of six primary red-giant stars with masses of M_⋆/M_⊙ = (12, 12, 16, 16, 20, 20) with respective radii R_⋆/R_⊙ = (173, 594, 322, 672, 872, 1247). Here, the smaller radii at a given mass corresponds to the RGB base, while the larger radii at a given mass corresponds to the RGB tip. For our CE inspiral calculations, we take the envelope structure to be constant in time.

As an NS inspirals in the CE, the envelope plasma supersonically flows past the NS and may be captured within the NS's accretion radius R_a = 2GM_NS/v², where v is the upstream flow velocity, which, in the NS's rest frame is the orbital velocity. If the upstream flow is homogeneous, then from Hoyle–Lyttleton (HL) theory (Hoyle & Lyttleton 1939), the accretion rate and local drag force obey

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{HL}}=\pi {R}_{{\rm{a}}}^{2}\rho v=\displaystyle \frac{4\pi \rho {G}^{2}{M}_{\mathrm{NS}}^{2}}{{v}^{3}},\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{F}_{{\rm{d}},\mathrm{HL}}={\dot{M}}_{\mathrm{HL}}v=\pi {R}_{{\rm{a}}}^{2}\rho {v}^{2}=\displaystyle \frac{4\pi \rho {G}^{2}{M}_{\mathrm{NS}}^{2}}{{v}^{2}},\end{eqnarray} \tag{ 6b }$$

where ρ is the upstream mass density. For NS accretion in stellar-envelope environments, the density and temperature may be high enough for neutrino cooling, such that the accretion rate exceeds the Eddington limit (Houck & Chevalier 1991). The envelope's local density scale height may be comparable in size to the NS accretion radius, which breaks the symmetry that HL theory assumes and thus requires a treatment of this effect.

To model the accretion and drag, we use the fitting formulae from De et al. (2020, see their Appendix A). The accretion and drag coefficients, C_a and C_d are defined such that the baryonic mass accretion rate and the local drag force are

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{b}}}={C}_{{\rm{a}}}{\dot{M}}_{\mathrm{HL}},\quad {C}_{{\rm{a}}}={C}_{{\rm{a}}}({ \mathcal M },q,{R}_{\mathrm{sink}}),\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&{F}_{{\rm{d}}}={C}_{{\rm{d}}}{F}_{{\rm{d}},\mathrm{HL}},\quad {C}_{{\rm{d}}}={C}_{{\rm{d}}}({ \mathcal M },q).\end{eqnarray} \tag{ 7b }$$

These coefficients are both functions of the upstream Mach number ${ \mathcal M }$ and the mass ratio q between the compact object and the enclosed mass within the orbit. The accretion coefficient also depends on the sink radius R_sink, given that the wind-tunnel simulations only resolve the accretion flow up to a sphere with radius 0.05R_a surrounding the point-mass accretor. Thus, some fraction of matter that flows into the region within 0.05R_a ultimately ends up accreting onto the NS. For each EoS, we use R_NS as the sink radius. In Appendix C, we describe in more detail how we compute these accretion and drag coefficients.

We plot in Figure 2 the stellar profiles of the density, upstream Mach number, polytropic exponent, and envelope binding energy (panels A, B, C, and D, respectively) for the primary masses of M_⋆/M_⊙ ∈ (12, 16, 20) that we consider here. For a given evolutionary stage and for most of the primary's radii, the δ_ρ parameter decreases as the primary's mass increases, such that the accretion rate will be greater and result in higher accreted mass. Primary stars that are smaller in size will have higher envelope binding energy, which thus requires more energy dissipation during the CE phase in order for successful envelope ejection. In Section 3, we quantify how much more NSs accrete when in envelopes with higher binding energies compared to less bound envelopes.

Zoom In Zoom Out Reset image size

Given the primary's envelope structure, we can now model the CE inspiral as follows. We approximate the orbital inspiral with Newtonian gravity (Blanchet 2014), given that on the scales of CE evolution, gravity is weak and the orbital velocities are nonrelativistic, i.e., v_orb/c ≪ 1. With this in mind, the orbital energy throughout the inspiral is E = − GM_NS m_⋆/(2a), where ${m}_{\star }={m}_{\star }(a)={\int }_{0}^{a}4\pi \rho {r}^{2}\,{dr}$ is the mass enclosed within the NS's separation from the primary's center. The orbital velocity at any given time obeys ${v}^{2}=G\left[{M}_{\mathrm{NS}}+{m}_{\star }(a)\right]/a$, since we consider the inspiral to be quasi-circular. The change in the binary orbital energy as the NS inspirals thus obeys

$$\begin{eqnarray}&&{dE}=\displaystyle \frac{\partial E}{\partial {m}_{\star }}{{dm}}_{\star }+\displaystyle \frac{\partial E}{\partial {M}_{\mathrm{NS}}}{{dM}}_{\mathrm{NS}}+\displaystyle \frac{\partial E}{\partial a}{da},\end{eqnarray} \tag{ 8 }$$

such that we can then solve for da/dt

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{da}}{{dt}} & = & a\left(-\displaystyle \frac{\dot{E}}{E}+\displaystyle \frac{{\dot{m}}_{\star }}{{m}_{\star }}+\displaystyle \frac{{\dot{M}}_{\mathrm{NS}}}{{M}_{\mathrm{NS}}}\right)\\ & = & \displaystyle \frac{a}{1-4\pi \rho {a}^{3}/{m}_{\star }}\left(-\displaystyle \frac{\dot{E}}{E}+\displaystyle \frac{{\dot{M}}_{\mathrm{NS}}}{{M}_{\mathrm{NS}}}\right),\end{array}\end{eqnarray} \tag{ 9 }$$

where ${\dot{m}}_{\star }=4\pi \rho {a}^{2}\dot{a}$ since we assume a static envelope. We take the energy decay rate to be the drag luminosity $\dot{E}=-{F}_{{\rm{d}}}({\delta }_{\rho })v$, which dominates over the gravitational-wave luminosity from the orbital motion, and which can be obtained using both Equations 6(b) and 7(b).

We summarize our integration procedure as follows. We first precompute the NS properties shown in Equations (4) and (5) as well as Appendix B for each EoS in our catalog, which are then stored as tables to interpolate from at each timestep of an orbital integration. We then explicitly integrate Equations (4), (5), and (9) to obtain M_NS, Ω, and a throughout the CE inspiral. The NS properties at each point in the NS's evolution correspond to a Hartle–Thorne NS with gravitational mass M_NS and spin Ω that we obtain from our precomputed tables. Our orbital integrations are carried out for each of our six primary stellar models and for each EoS in our catalog, varying the initial NS gravitational mass M_NS,0 and initial NS spin Ω₀. We terminate these orbital integrations when the dissipated orbital energy ΔE_orb = E(a) − E(a₀) is equal to the primary envelope's binding energy E_env,bind given by

$$\begin{eqnarray}&&{E}_{\mathrm{env},\mathrm{bind}}={\int }_{{m}_{\star }(r)}^{{m}_{\star }({R}_{\star })}\left(u-\displaystyle \frac{{{GM}}_{\star }(r)}{r}\right)\,{dM},\end{eqnarray} \tag{ 10 }$$

where u is the stellar fluid's internal energy and where the integration coordinate is the primary's mass coordinate. This amounts to assuming a CE efficiency parameter (e.g., Webbink 1984) of α_CE = 1.

We plot the NS evolution for the often fiducial case of a pre-CE NS mass of 1.4M_⊙ and a primary of 12M_⊙ in the top panel of Figure 3. For this case, the primary is taken to be at the RGB base such that a₀ = 173R_⊙ and we take the initial NS spin to be Ω₀/(2π) = 50 Hz. The black curves correspond to Φ = 1, i.e., not accounting for NS binding energy. Each curve corresponds to a different EoS in our catalog.

Zoom In Zoom Out Reset image size

In all cases, the NS accretes no more than a few percent of its pre-CE mass, due to the suppressed accretion rate from the envelope density gradient. The gravitational-mass gain as well as the spin-up further decreases, since some of the accreted baryon mass-energy is converted into binding energy. In the bottom panel of Figure 3, we plot the final spin-up ΔΩ_f/(2π) = (Ω_f − Ω₀)/(2π) and the final gravitational-mass gain for each EoS model and for varying initial NS spins of Ω₀/(2π) = (10, 50, 100, 200, 500) Hz as blue, orange, green, red, and purple points, respectively. Higher initial NS spins increase the NS binding energy, such that less gravitational mass is gained and the spin-up decreases.

With different EoSs, there is an anticorrelation between the mass gain and spin-up, where an EoS that allows for higher gravitational-mass gain results in a lower spin-up when starting with the same initial NS spin. A larger increase in ΔM_NS is a result of a larger Φ_NS, i.e., higher baryon mass converted to gravitational mass. The gravitochemical potential Φ_NS is proportional to the inertia I_NS, such that less compact NSs are harder to spin-up because they have higher Φ_NS and higher I_NS. Conversely, more compact NSs will gain less gravitational mass and spin-up more because they have lower Φ_NS and lower I_NS.

Given that the relativistic mass deficit is greater for more massive NSs (see Figure 1), we then vary the pre-CE NS mass. We run inspirals for the following set of pre-CE NS masses M₀/M_⊙ ∈ [1.2, 1.8] with a step size of 0.1M_⊙. We plot in the top panels of Figure 4 the mean accreted masses when varying the EoS as solid lines with the shaded region corresponding to the ±2σ deviation. The dashed lines correspond to Φ = 1, i.e., taking the accreted gravitational mass to be equivalent to the accreted baryonic mass. The width of the dashed line encompasses the ±2σ region. In the bottom panels of Figure 4, we plot the corresponding spin-up ΔΩ_f = (Ω_f,0 − Ω₀)/(2π).

Zoom In Zoom Out Reset image size

An increasing pre-CE NS mass results in a systematically decreasing accreted NS mass across all of our models. This is because at constant α_CE, having a more massive NS results in a larger dissipated orbital energy, such that envelope ejection is achieved at wider separations and such that the accreted baryonic mass is reduced compared to lower-mass NSs. It remains to be seen whether or not this trend will hold in global 3D hydrodynamic CE simulations when the initial NS mass is varied. Models at the RGB base result in higher accreted mass and spin-up, which is due to the larger envelope binding energy from their smaller sizes as compared to the RGB tip (bottom panel of Figure 2).

As previously shown in Figure 1, NSs with higher gravitational mass will convert a larger fraction of the accreted mass into binding energy. We plot in Figure 5 the distributions of the ratio of the accreted gravitational mass to the accreted baryonic mass from our RGB base models. In each panel, each distribution from bottom ascending to top is for initial NS gravitational masses of M_NS,0/M_⊙ = (1.2, 1.4, 1.6.,1.8), respectively. The green and magenta curves correspond to initial NS spins of 50 Hz and 200 Hz, respectively. We also plot a black dashed curve that corresponds to using Φ_TOV,0 as a fast approximation, i.e., the gravitochemical potential of the nonrotating NS at its initial properties.

Zoom In Zoom Out Reset image size

Higher initial spins tend to decrease ΔM_NS/ΔM_b, though the model with M_NS,0 = 1.8M_⊙ and M_⋆ = 20M_⊙ at the RGB base exhibits opposite behavior, albeit slight. Since lower-mass NSs accrete more gravitational mass compared to the higher-mass NSs in our models, the differences between the ratio distributions at various spins is also higher as well. Distributions for the RGB tip case will be similar, though the separation between distributions at the same initial NS mass will be smaller since the RGB tip cases resulted in less mass gain (Figure 4).