Gravitational Waves from Accreting Neutron Stars Undergoing Common-envelope Inspiral

A star in a close binary can expand over the course of its stellar evolution to the point where it envelops its companion star, thus initiating a CE phase. It is thought that the CE phase involves a contact phase, followed by a plunge-in phase in which the companion undergoes a fast inspiral, and finally a self-regulated phase (Ivanova et al. 2013). If enough orbital energy is dissipated in the envelope, the system transitions to a post-CE phase in which the remainder of the envelope is ejected. In this paper, we ignore the initial plunge, which typically lasts one to two orbits in 3D simulations, and treat the entire orbital evolution as self-regulated and driven by local gravitational drag in the vicinity of the NS companion. The adequacy of this approximation remains to be shown for the binary system parameters we consider here. We are undertaking this task using 3D simulations that will be presented in a future paper.

The energy formalism is commonly used to parameterize CE evolution (e.g., Webbink 1984) using the efficiency α_CE with which the orbital energy release unbinds the envelope. The binding energy of the envelope, E_bind, is then related to the orbital energy dissipated over the course of the inspiral, ΔE_orb, via

$$\begin{eqnarray}&&{E}_{\mathrm{bind}}={\alpha }_{\mathrm{CE}}{\rm{\Delta }}{E}_{\mathrm{orb}}={\alpha }_{\mathrm{CE}}\left(\displaystyle \frac{{{GM}}_{\mathrm{core}}{M}_{\mathrm{NS}}}{2{a}_{{\rm{f}}}}-\displaystyle \frac{{{GM}}_{1}{M}_{\mathrm{NS}}}{2{a}_{{\rm{i}}}}\right),\end{eqnarray} \tag{ 1 }$$

where M₁ is the mass of the primary, M_core is the mass of the primary's core, M_NS is the mass of the NS companion, a_i is the initial orbital separation (taken here to be the radius of the primary), and a_f is the final separation. The core–envelope split for the primary is an uncertainty, such that M_core usually includes both the mass of the primary's core and the leftover envelope that remains bound to it. Some observed systems require α_CE > 1, implying that additional sources of energy contribute to unbinding the envelope.⁴

Neutron stars in low-mass X-ray binaries can accrete from their companion stars at rates as high as the Eddington accretion rate, given by

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{Edd}}=\displaystyle \frac{4\pi {{GM}}_{\mathrm{NS}}}{\zeta c{\kappa }_{{\rm{T}}}},\end{eqnarray} \tag{ 2 }$$

where κ_T is the Thomson scattering opacity and ζ is the efficiency with which the mass accretion rate powers the luminosity, typically taken to be ζ ∼ 0.1. We re-express the Eddington accretion limit in terms of typical NS parameters as

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{Edd}}\approx 3.5\times {10}^{-8}\left(\displaystyle \frac{{M}_{\mathrm{NS}}}{1.33{M}_{\odot }}\right){\left(\displaystyle \frac{{\kappa }_{{\rm{T}}}}{0.34{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}}\right)}^{-1}\ {M}_{\odot }\,{\mathrm{yr}}^{-1}.\end{eqnarray} \tag{ 3 }$$

An NS accreting mass can also accrete angular momentum and be spun up. The accretion torque on an NS can be approximated as (e.g., Bildsten 1998)

$$\begin{eqnarray}&&{N}_{\mathrm{acc}}\approx \dot{M}\sqrt{{{GM}}_{\mathrm{NS}}{R}_{\mathrm{eq}}},\end{eqnarray} \tag{ 4 }$$

where M_NS is the mass of the NS and R_eq is the radius of the NS at its equator. If we assume torque balance, so that GW radiation losses equal the accretion torque (discussed in more detail in Section 4), the estimated strain amplitude from the accreting neutron star is

$$\begin{eqnarray}&&{h}_{0}\approx 1.8\times {10}^{-26}{\left(\displaystyle \frac{D}{2.8\mathrm{kpc}}\right)}^{-1}{\left(\displaystyle \frac{\dot{M}}{{\dot{M}}_{\mathrm{Edd}}}\right)}^{1/2}{\left(\displaystyle \frac{{f}_{\mathrm{GW}}}{600\mathrm{Hz}}\right)}^{-1/2},\end{eqnarray} \tag{ 5 }$$

where we have scaled quantities to those appropriate for Sco X-1 and assumed M_NS = 1.33M_⊙ and R_eq = 10 km. Here, D is the distance to the LMXB and f_GW is the GW frequency.

The accretion flow during CE inspiral can be parametrized using the theory of BHL accretion (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Edgar 2004). The characteristic length scale for accretion onto an object of mass M moving at a speed ${v}_{\infty }$ relative to a uniform medium with density ${\rho }_{\infty }$ and sound speed c_s is given by

$$\begin{eqnarray}&&{R}_{\mathrm{BHL}}\equiv \displaystyle \frac{2{GM}}{{v}_{\infty }^{2}+{c}_{{\rm{s}}}^{2}}.\end{eqnarray} \tag{ 6 }$$

If, following MR15, we ignore the sound speed, then the characteristic length scale is the Hoyle–Lyttleton (HL) accretion radius, R_HL (Hoyle & Lyttleton 1939). We define the following parameters:

$$\begin{eqnarray}&&{R}_{{\rm{a}}}\equiv {R}_{\mathrm{HL}}=\displaystyle \frac{2{{GM}}_{\mathrm{NS}}}{{v}_{\infty }^{2}},\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&{\tau }_{{\rm{a}}}\equiv \displaystyle \frac{{R}_{{\rm{a}}}}{{v}_{\infty }},\end{eqnarray} \tag{ 7b }$$

$$\begin{eqnarray}&&{{ \mathcal M }}_{\infty }\equiv \displaystyle \frac{{v}_{\infty }}{{c}_{{\rm{s}}}}.\end{eqnarray} \tag{ 7c }$$

Here, R_a is the accretion radius of the NS, τ_a is the inspiral accretion timescale, and ${{ \mathcal M }}_{\infty }$ is the local upstream Mach number. The HL accretion rate is then

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{HL}}=\pi {R}_{{\rm{a}}}^{2}{\rho }_{\infty }{v}_{\infty }.\end{eqnarray} \tag{ 8 }$$

Assuming this form for the accretion rate, the gravitational drag force (e.g., Iben & Livio 1993) and corresponding energy loss rate experienced by the inspiraling companion can be estimated using

$$\begin{eqnarray}&&{F}_{{\rm{d}},\mathrm{HL}}=\pi {R}_{{\rm{a}}}^{2}{\rho }_{\infty }{v}_{\infty }^{2}={\dot{M}}_{\mathrm{HL}}{v}_{\infty },\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{HL}}=\pi {R}_{{\rm{a}}}^{2}{\rho }_{\infty }{v}_{\infty }^{3}={\dot{M}}_{\mathrm{HL}}{v}_{\infty }^{2}={F}_{{\rm{d}},\mathrm{HL}}{v}_{\infty }.\end{eqnarray} \tag{ 9b }$$

The actual accretion rate, drag force, and energy loss rate are reduced relative to the HL rates during CE inspiral as discussed in Section 5.2.

In hypercritical accretion, thermal neutrino cooling is sufficient to allow the accretion rate to exceed the Eddington limit (Houck & Chevalier 1991). Once the NS surface temperature exceeds ∼1 MeV, photons can be trapped in the flow and a neutrino-cooling layer is expected to form outside the NS surface. The photon-trapping radius is approximately

$$\begin{eqnarray}{R}_{\mathrm{tr}} & = & \displaystyle \frac{\dot{M}{\kappa }_{{\rm{T}}}}{4\pi c}\\ & \approx & 5.7\times {10}^{13}\left(\displaystyle \frac{\dot{M}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)\left(\displaystyle \frac{{\kappa }_{{\rm{T}}}}{0.34\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}}\right)\ \mathrm{cm}.\end{eqnarray} \tag{ 10 }$$

The shock radius is approximately

$$\begin{eqnarray}&&{R}_{\mathrm{sh}}\approx 1.6\times {10}^{8}{\left(\displaystyle \frac{\dot{M}}{{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{-0.37}\ \mathrm{cm}.\end{eqnarray} \tag{ 11 }$$

The general criterion for super-Eddington accretion to occur via efficient thermal neutrino cooling is R_tr > R_sh, which implies

$$\begin{eqnarray}\dot{M}\gtrsim {\dot{M}}_{\mathrm{hyper}} & \equiv & 8.9\times {10}^{-5}{\left(\displaystyle \frac{{\kappa }_{{\rm{T}}}}{0.34{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}}\right)}^{-0.73}\\ & & \sim {10}^{4}{\dot{M}}_{\mathrm{Edd}},\end{eqnarray} \tag{ 12 }$$

where the last expression holds for typical NS masses. In our model treatment, the transition from Eddington to hypercritical accretion occurs instantaneously. Physically, the density of the gas falling into the accretion radius of the NS is high enough such that temperatures rise to the point where thermal neutrino cooling will start to occur and allow for hypercritical accretion.

Farther away from the NS, transport of thermal radiation can efficiently cool the material that is gravitationally captured by the NS. As an NS inspirals within the CE, it creates density wakes that provide additional gravitational drag. Photon bubbles can travel faster through the less dense regions (e.g., Gammie 1998) and transport thermal energy away from the accreting material more efficiently than advection and conduction. Regions near the NS may also be unstable to the neutrino-bubble (Socrates et al. 2005) and photon-bubble instability, such that perturbations in the flow are seeded, grow exponentially, and saturate into shock trains (e.g., Turner et al. 2005; Begelman 2006; Turner et al. 2007).

As emphasized by MR15, a finite density gradient in the stellar envelope breaks the symmetry of HL accretion. Material gravitationally captured within the accretion radius may not be accreted onto the inspiralling companion due to the flow's angular momentum. The background density gradient can be parametrized using the nondimensional quantity

$$\begin{eqnarray}&&{\epsilon }_{\rho }\equiv \displaystyle \frac{{R}_{{\rm{a}}}}{{H}_{\rho }}=-\displaystyle \frac{{R}_{{\rm{a}}}}{\rho }\displaystyle \frac{d\rho }{{dr}},\end{eqnarray} \tag{ 13 }$$

where ρ is the density of the background medium and H_ρ is the density scale height. MR15 found that the accretion rate decreases relative to ${\dot{M}}_{\mathrm{HL}}$ with steeper density gradients. (We describe the MR15 model in more detail in Section 5.) This result suggests that NSs can survive the CE phase since they accrete only a small fraction of their initial mass during CE inspiral. The actual accretion rate $\dot{M}$ during inspiral appears to satisfy ${\dot{M}}_{\mathrm{Edd}}\lt {\dot{M}}_{\mathrm{hyper}}\lt \dot{M}\lt {\dot{M}}_{\mathrm{BHL}}\lt {\dot{M}}_{\mathrm{HL}}$. Hence, using ${\dot{M}}_{\mathrm{HL}}$ for the accretion rate should provide an upper limit on the GW emission due to accretion during CE inspiral, assuming torque balance. Since the accretion rate suppression in the MR15 model can be significant, we use MR15 instead of the HL rate in what follows.

NS surface deformations away from axisymmetry, i.e., mountains, produce a time-varying quadrupole moment and thus GW emission. There are two types of NS mountains that may be present as the NS accretes during CE evolution: thermal mountains and magnetic mountains. Bildsten (1998) first described how electron captures in the crust of Eddington-accreting NSs may generate density asymmetries on the NS surface if a transverse temperature gradient is present, i.e., either from compositional variation, asymmetric local accretion flow, buried magnetic fields, etc., and thus lead to GW emission. The typical temperature of the crust of Eddington-accreting NSs is of order 10⁸ K such that the crust remains solid. The maximum quadrupole moment of such thermal mountains depends on the properties of the NS crust where the crust ellipticity can be approximated as (Glampedakis & Gualtieri 2017)

$$\begin{eqnarray}&&{\varepsilon }_{\mathrm{th}}\lesssim \displaystyle \frac{{\mu }_{\mathrm{cr}}{\sigma }_{\mathrm{br}}{V}_{\mathrm{cr}}}{{{GM}}_{\mathrm{NS}}^{2}/{R}_{\mathrm{NS}}}\sim {10}^{-5}\left(\displaystyle \frac{{\sigma }_{\mathrm{br}}}{0.1}\right),\end{eqnarray} \tag{ 18 }$$

where μ_cr is the shear modulus of the crust, V_cr is the volume of the crust, and σ_br is the crustal breaking strain. For magnetic mountains on NSs with magnetar-level field strengths, B ∼ 10¹⁵ G, the ellipticity may take on values of ε_mag ∼ 10⁻⁶–10⁻⁵. Our work is agnostic to the mechanism by which mountains are produced and only requires that such quadrupole moments are achievable. Several previous works have estimated what values of the maximum quadrupole moment may be sustained by the NS crust. Haskell et al. (2006) applied the formalism presented by Ushomirsky et al. (2000) to both accreted and non-accreted NS crusts, and determined the maximum quadrupole moment to be Q_max ≈ 2 × 10³⁹ g cm². Horowitz & Kadau (2009) performed 3D molecular dynamics simulations of local regions of NS crust to obtain estimates for the breaking strain. With their estimate for the breaking strain, Horowitz & Kadau (2009) obtained a maximum NS crust ellipticity of ε ≈ 4 × 10⁻⁶. With this bound on the ellipticity, we thus use a maximum quadrupole moment of Q_max = 4 × 10³⁹ g cm² for this work. With this upper limit on the maximum quadrupole moment, the strain amplitudes we estimate in Section 6 should be regarded as upper limits.

When the accretion torque approximately balances the GW radiation torque (N_GW + N_acc ≈ 0), the quadrupole moment can be expressed as a function of the accretion rate by using Equations (4) and (15) to obtain

$$\begin{eqnarray}&&{Q}_{\mathrm{tb}}={\left(\displaystyle \frac{5}{32}\displaystyle \frac{{c}^{5}}{G{{\rm{\Omega }}}^{5}}\dot{M}\sqrt{{{GM}}_{\mathrm{NS}}{R}_{\mathrm{eq}}}\right)}^{1/2}.\end{eqnarray} \tag{ 19 }$$

Since the mass and accretion rate change slowly compared to the GW period, the implicit time dependence they introduce into Q modulates the GW strain amplitude as the NS inspirals within the CE.

We treat our model NS undergoing CE evolution as beginning in torque balance and spinning up if the required quadrupole moment to maintain torque balance is greater than the maximum quadrupole moment that the NS crust can sustain. At hypercritical accretion rates, temperatures are of order MeV such that the the solid NS crust melts into liquid and erases surface deformations. As long as there exists a solid NS crust, thermal mountains and/or magnetic mountains may still be formed before they melt away. Our treatment of the crust-melting process is described in Section 5.6. For the magnetar scenario, the crust remains deformed by the global magnetic field as it is melting. For the thermal-mountain scenario, the electron-capture rates just need to vary enough with temperature in order to produce the required quadrupole moment. We include a derivation of the required electron-capture rate temperature dependence in Appendix. At accretion rates relevant to the fallback scenario, the accretion torque is large enough such that the NS is spun up to secular instability, which excites r-modes that produce GWs. In Section 6, we will show that accreting NSs in the CE case do not spin up to secular instability.

We compute the single-star evolution of a 12 M_⊙ and a 20 M_⊙ primary with solar metallicity using the MESA stellar evolution code version 8845 (Paxton et al. 2011, 2013, 2015). We consider CE phases occurring at two different evolutionary stages: at the base and the tip of the red giant branch (RGB). Due to mass loss, the precise donor star masses used are slightly smaller than their zero-age main-sequence (ZAMS) values.⁵ During the inspiral, the primary's structure is taken to be fixed. As pointed out by MR15, this approximation is most reasonable for large binary mass ratios. The density profiles and nondimensional density gradients of our stellar models are shown in Figure 1. A density inversion appears near the surface of the stellar models because of the interaction of the iron opacity bump with the nearly Eddington radiative flux in the envelope. This produces a spike in _ρ. Three-dimensional radiation hydrodynamics simulations (Jiang et al. 2015) show that in low-density regions, such inversions develop a porous structure that allows radiative and convective energy transport to compete, expanding the envelope and partially eliminating the inversion. In our calculations, the presence of the inversion does not matter significantly for the GW prediction. Given the present uncertainty regarding how to properly treat such inversions in 1D stellar evolution codes, we have chosen to leave the inversion unmodified in this paper.

Zoom In Zoom Out Reset image size

To estimate the NS accretion rate at a given radius within the primary, we use the HL rate expressions (Equations (8) and (9)) modified using the model of MR15. MR15 performed simulations of accretion onto a spherical sink region with a background density gradient as a model for the local accretion dynamics during CE evolution. They approximate the gravitational drag obtained in their simulations using

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{{\rm{d}},\mathrm{MR}15}({\epsilon }_{\rho })}{{F}_{{\rm{d}},\mathrm{HL}}}\approx {f}_{1}+{f}_{2}{\epsilon }_{\rho }+{f}_{3}{\epsilon }_{\rho }^{2},\end{eqnarray} \tag{ 20 }$$

where the coefficients f_i (i = 1, 2, 3) are given by

$$\begin{eqnarray}&&{f}_{1}=1.91791946,\end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray}&&{f}_{2}=-1.52814698,\end{eqnarray} \tag{ 21b }$$

$$\begin{eqnarray}&&{f}_{3}=0.75992092.\end{eqnarray} \tag{ 21c }$$

The accretion rate in their simulations is approximately

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{\dot{M}}_{\mathrm{MR}15}({\epsilon }_{\rho })}{{\dot{M}}_{\mathrm{HL}}}\right)\approx {\mu }_{1}+\displaystyle \frac{{\mu }_{2}}{1+{\mu }_{3}{\epsilon }_{\rho }+{\mu }_{4}{\epsilon }_{\rho }^{2}},\end{eqnarray} \tag{ 22 }$$

with accretion rate coefficients μ_i (i = 1, 2, 3, 4) given by

$$\begin{eqnarray}&&{\mu }_{1}=-2.14034214,\end{eqnarray} \tag{ 23a }$$

$$\begin{eqnarray}&&{\mu }_{2}=1.94694764,\end{eqnarray} \tag{ 23b }$$

$$\begin{eqnarray}&&{\mu }_{3}=1.19007536,\end{eqnarray} \tag{ 23c }$$

$$\begin{eqnarray}&&{\mu }_{4}=1.05762477.\end{eqnarray} \tag{ 23d }$$

We plot the accretion and drag ratios computed using Equations (20) and (22) and the density gradients from our MESA stellar models in Figure 2. Note that in the region of the density inversion the drag is significantly enhanced while the accretion rate is significantly diminished relative to the HL prediction. Hence, including the inversion in the primary model causes the NS to pass through this region rapidly without emitting much GW energy. In some respects, this mimics the effect of the initial dynamical plunge seen in 3D CE simulations, which we do not explicitly include.

Zoom In Zoom Out Reset image size

Accretion onto the NS during the CE phase is a fairly uncertain process, and the wind-tunnel model for the orbital evolution may not work well in our case since the primary star structure may respond to the orbiting NS. To consider the possibility that the accretion rate is diminished even further below the MR15 model, we parameterize the accretion rate using an efficiency factor η ≤ 1. Since unstable cooling solutions exist for ${\dot{M}}_{\mathrm{Edd}}\lt \dot{M}\lt {\dot{M}}_{\mathrm{hyper}}$, the accretion rate we use is

$$\begin{eqnarray}\dot{M}=\left\{\begin{array}{ll}\eta {\dot{M}}_{\mathrm{MR}15} & \eta {\dot{M}}_{\mathrm{MR}15}\geqslant {\dot{M}}_{\mathrm{hyper}}\ {\rm{or}}\ \leqslant {\dot{M}}_{\mathrm{Edd}},\\ {\dot{M}}_{\mathrm{Edd}} & {\dot{M}}_{\mathrm{Edd}}\lt \eta {\dot{M}}_{\mathrm{MR}15}\lt {\dot{M}}_{\mathrm{hyper}}.\end{array}\right.\end{eqnarray} \tag{ 24 }$$

The case $\dot{M}\lt {\dot{M}}_{\mathrm{Edd}}$ is implicitly allowed in our model, but in our integrations we find that the accretion rate is always $\geqslant {\dot{M}}_{\mathrm{Edd}}$. The corresponding drag force is then

$$\begin{eqnarray}&&{F}_{{\rm{d}}}=\eta {F}_{{\rm{d}},\mathrm{MR}15},\end{eqnarray} \tag{ 25 }$$

where, even for Eddington-limited accretion rates, momentum dissipation occurs within the accretion radius of the NS, which contributes to the drag force.

We consider the same orbital evolution model as MR15. In this model, the inspiral is driven by local gravitational drag within the envelope such that the orbital decay rate, $\dot{a}$, obeys

$$\begin{eqnarray}&&\dot{a}=\dot{E}\displaystyle \frac{{da}}{{{dE}}_{\mathrm{orb}}},\end{eqnarray} \tag{ 26 }$$

where $\dot{E}$ is the energy dissipated by gravitational drag,

$$\begin{eqnarray}&&\dot{E}=-{F}_{{\rm{d}}}({\epsilon }_{\rho }){v}_{\infty }.\end{eqnarray} \tag{ 27 }$$

The gravitational drag force F_d(_ρ) is obtained from Equation (25). Assuming a nearly circular orbit, the inspiral velocity ${v}_{\infty }$ obeys

$$\begin{eqnarray}&&{v}_{\infty }^{2}=\displaystyle \frac{G({M}_{\mathrm{NS}}+m(a))}{a},\end{eqnarray} \tag{ 28 }$$

where $m(a)={\int }_{0}^{a}4\pi \rho (r){r}^{2}\,{dr}$ is the envelope mass enclosed within the orbital separation. The orbital energy E_orb is

$$\begin{eqnarray}&&{E}_{\mathrm{orb}}=-\displaystyle \frac{{{GM}}_{\mathrm{NS}}m(a)}{2a}.\end{eqnarray} \tag{ 29 }$$

We follow the evolution of the accretion rate, $\dot{M}$, and the orbital separation, a, by integrating Equation (26). Following MR15, we terminate the orbital evolution when the dissipated energy begins to exceed the binding energy of the envelope,

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{orb}}\geqslant {E}_{\mathrm{bind}},\end{eqnarray} \tag{ 30 }$$

where we take α_CE = 1 and where the binding energy is

$$\begin{eqnarray}&&{E}_{\mathrm{bind}}={\int }_{m(a)}^{M}\left(u-\displaystyle \frac{{Gm}(a)}{r}\right)\,{dm}.\end{eqnarray} \tag{ 31 }$$

Here, u is the specific internal energy of the envelope gas. This criterion may not correspond to the true end of the CE phase.

From the strain amplitude and/or its time derivative, it may be possible to constrain the orbital decay rate of the NS. We thus define the quantity

$$\begin{eqnarray}&&\xi \equiv \displaystyle \frac{{\tau }_{\mathrm{ins}}}{{\tau }_{\mathrm{orb}}},\end{eqnarray} \tag{ 32 }$$

where the inspiral timescale τ_ins and the orbital timescale τ_orb respectively obey

$$\begin{eqnarray}&&{\tau }_{\mathrm{ins}}=\displaystyle \frac{a}{| \dot{a}| },\end{eqnarray} \tag{ 33a }$$

$$\begin{eqnarray}&&{\tau }_{\mathrm{orb}}=\displaystyle \frac{2\pi }{\sqrt{G({M}_{\mathrm{NS}}+m(a))}}{a}^{3/2},\end{eqnarray} \tag{ 33b }$$

and $\dot{a}$ obeys Equation (26). We use the ξ parameter to characterize how the orbital decay of the NS evolves within the CE, and describe its features as follows. The limit of ξ ≫ 1 corresponds to a slow migration, where the inspiral occurs over many orbital periods. The range ξ ∼ 1 corresponds to a merger within an orbital period, i.e., a moderate inspiral. The limit of ξ ≪ 1 corresponds to a merger on a timescale much shorter than the orbital period, i.e., a rapid inspiral. In terms of CE evolution, the plunge-in phase may be characterized by ξ ≪ 1 increasing to ξ ∼ 1, followed by the self-regulated phase, where ξ ≫ 1. Since we are treating the CE inspiral as quasi-circular, cases where ξ ≪ 1 indicate a breakdown in the validity of this assumption.

To self-consistently determine the equatorial radius given the NS spin, we adopt the same approach as PT12. We treat the NS as an axisymmetric Maclaurin spheroid and take any internal viscous dissipation to be negligible since it operates on timescales several orders of magnitude longer than dissipation via GW radiation. The GW emission then reacts immediately to changes in the accretion rate.

The NS spin is parametrized using the spin parameter $\beta =T/| W|$, where T is the rotational kinetic energy and W is the gravitational potential energy. From the Maclaurin spheroid model, the spin parameter can be expressed in terms of the NS ellipticity as

$$\begin{eqnarray}&&\beta =\displaystyle \frac{3}{2{e}^{2}}\left[1-\displaystyle \frac{e{(1-{e}^{2})}^{1/2}}{{\sin }^{-1}e}\right]-1.\end{eqnarray} \tag{ 34 }$$

The NS spin obeys

$$\begin{eqnarray}&&{{\rm{\Omega }}}^{2}=\displaystyle \frac{2\pi G\bar{\rho }}{{q}_{n}}\left[\displaystyle \frac{{(1-{e}^{2})}^{1/2}}{{e}^{3}}(3-2{e}^{2}){\sin }^{-1}e-\displaystyle \frac{3(1-{e}^{2})}{{e}^{2}}\right],\end{eqnarray} \tag{ 35 }$$

where $\bar{\rho }=3M/4\pi {R}_{0}^{3}$ is the average density, R₀ is the radius of the NS with zero rotation, and

$$\begin{eqnarray}&&{q}_{n}=(1-n/5){\kappa }_{n},\end{eqnarray} \tag{ 36 }$$

with n as the polytropic index and κ_n as a constant of order unity (see Table 1 in Lai et al. 1993). We take n = 0.5, which is stable against mass shedding. From the ellipticity, the NS's equatorial radius is

$$\begin{eqnarray}&&{R}_{\mathrm{eq}}=\displaystyle \frac{\bar{R}}{{(1-{e}^{2})}^{1/6}},\end{eqnarray} \tag{ 37 }$$

where $\bar{R}$ is the mean radius of the NS given by

$$\begin{eqnarray}&&\bar{R}={R}_{0}{\left[\displaystyle \frac{{\sin }^{-1}e}{e}{(1-{e}^{2})}^{1/6}(1-\beta )\right]}^{-n/(3-n)}.\end{eqnarray} \tag{ 38 }$$

We take the spin parameter to be bounded by β_sec, so that GW emission prevents the NS from becoming secularly unstable. We plot the ellipticity and spin frequency, f_s = Ω/(2π), as a function of the spin parameter in Figure 3.

Zoom In Zoom Out Reset image size

The NS crust ellipticity is the source of GWs in our model of NS–CE evolution. The temperature of the NS crust while accreting at the Eddington rate is $\sim {10}^{8}\ {\rm{K}}$ (e.g., Bildsten 1998). At MeV temperatures relevant to hypercritical accretion, the outer crust melts into ocean, which erases surface deformations that were previously present. The lower layers may still have surface deformations and continue to generate GWs before melting away. Once the entire crust melts, then there is no GW emission from the NS for the rest of the CE evolution. In this scenario, we expect there to be GW emission as long as the crust-melting timescale is longer than the CE inspiral timescale. Physically, the local crust-melting rate depends on a number of factors including the composition, density, etc., which are uncertain. As a conservative assumption, we take the accretion rate as a bound on the crust-melting rate, which sets the crust-melting timescale to

$$\begin{eqnarray}&&{\tau }_{\mathrm{melt}}\sim \displaystyle \frac{{M}_{\mathrm{cr}}}{\dot{M}},\end{eqnarray} \tag{ 39 }$$

A lower crust-melting rate would be advantageous as it would allow for the NS crust to remain present on a longer timescale, thus allowing for more GW emission to occur during the CE inspiral. We approximate the crust to be close to constant density such that the ellipticity of the crust is a function of the crust mass that is present, thus setting the maximum quadrupole moment to

$$\begin{eqnarray}{Q}_{\max } & = & 4\times {10}^{39}\left(\displaystyle \frac{{M}_{\mathrm{acc}}}{{M}_{\mathrm{cr}}}\right)\ {\rm{g}}\,{\mathrm{cm}}^{2}\\ & = & 4\times {10}^{39}\left(\displaystyle \frac{{M}_{\mathrm{acc}}}{0.05{M}_{\odot }}\right)\ {\rm{g}}\,{\mathrm{cm}}^{2},\end{eqnarray} \tag{ 40 }$$

where M_acc is the amount of mass the NS accretes during CE inspiral and where we have taken the initial crust mass to be M_cr = 0.05M_⊙ (Barkov & Komissarov 2011). The model for the quadrupole moment is summarized as

$$\begin{eqnarray}&&Q=\min ({Q}_{\max },{Q}_{\mathrm{tb}}),\end{eqnarray} \tag{ 41 }$$

where Q_tb is the quadrupole moment obtained from torque balance, Equation (19), as described in Section 4.1. With the above prescription, the NS crust sustains a maximum quadrupole moment and spins up according to Equation (17) if Q_tb < Q_max. If Q_tb > Q_max, then the NS maintains torque balance. With the spin frequency Ω and the quadrupole moment Q, the strain amplitude h₀ is computed using Equation (16).

The parameter values we use for the spinning NS are M_NS = 1.33 M_⊙ and R_NS = R₀ = 10 km. We take the NS to have a spin frequency of Ω/(2π) = 250 Hz with the corresponding ellipticity and spin period being e₀ = 0.177 and P₀ = 4 ms, respectively. In order to factor out the distance to the source, we will present the strain amplitude relative to its Eddington-accreting NS in torque balance spinning at the initial rotation rate, P₀. The initial separation is taken to be a_i = R₁, where R₁ is the radius of the primary. We summarize the parameter values in Table 1.

Given the NS parameters and the primary structure, we determine the orbital separation, accretion rate, and strain amplitude as functions of time by integrating Equation (26). The resulting evolution for the 12 M_⊙ and 20 M_⊙ models is shown in Figures 4 and 5.

Zoom In Zoom Out Reset image size

In all cases we find that the accreted mass (panel (E) in Figures 4 and 5) is less than ∼0.1 M_⊙ at the end of the integration, consistent with MR15. Thus, the NS does not accrete enough mass to collapse into a BH. However, for both initial masses, the inspiral behavior (panel (A)) is strikingly different for stars at different evolutionary stages. For donors at the tip of the RGB, the NS appears to end the inspiral at a finite separation, with ξ increasing. However, for donors at the base of the RGB, the NS undergoes a rapid plunge upon reaching the termination criterion, suggesting that these systems will merge and form Thorne–Żytkow objects. In these cases, the accretion rate increases rapidly at the end of the inspiral. We have tested the sensitivity of this result to numerical effects by carrying out the integrations with 1/10 the time step and integrating past the termination criterion, finding that the behavior is unchanged.

Zoom In Zoom Out Reset image size

Careful comparison with the model profiles (Figure 1) suggests that the more evolved stars avoid merger because they reach the termination criterion before the NS reaches the hydrogen-burning shell, where _ρ displays a discontinuity. Although reaching a burning shell should violate the assumption that the donor profile is static (and the final plunge certainly violates the quasi-circular orbit assumption), the response of the star at this point might be to expand and thus further increase the drag. However, note from panel (D) that except for the more evolved cases with η = 1.0 and the very end of the merger cases (if run past the termination criterion), generally ξ > 1. This suggests that the quasi-circular approximation is at least self-consistent for these cases.

The strain amplitude (panel (H)) is typically louder than an Eddington-limited NS, particularly for the merger cases, implying that NS–CE systems can potentially be louder sources of GWs than LMXBs. This also implies that NS–CE systems may potentially be detected out to greater distances than Galactic LMXBs, which we discuss in more detail in Section 7.3. Since the accretion rate monotonically increases over the course of the inspiral (see panel (C) of Figures 4 and 5), the strain amplitude also increases with time, but experiences a turnover as a significant amount of crust has melted. For the more evolved stars, it levels off at ∼3×–10× its initial value, while for the less evolved stars, which merge, the strain amplitude increases by another factor of ∼3–10 over a timescale of less than a year. The loudest NS–CE systems should therefore be the ones with more compact, less evolved donors. Merging and non-merging cases can be distinguished based on the sign of the second time derivative of the strain amplitude.

Reducing the accretion efficiency η below 1.0 can cause the accretion rate to fall within the region of unstable cooling solutions found by Houck & Chevalier (1991). This is seen in the η = 0.1 models. The accretion rates during those periods are thus modulated by the Eddington limit, making the system no louder than an LMXB. The gravitational drag decreases when the accretion rate is reduced, resulting in a longer inspiral timescale. Thus, the η = 0.1 models have an inspiral duration ∼10 times longer than for the η = 1 models. We tabulate the duration of the inspiral, the factor by which the initial orbital separation decreases, and the factor by which the strain amplitude increases in our models in Table 2.

The spin evolution is shown in panel (F) of Figures 4 and 5, and spin ratios and values are tabulated in Tables 2 and 3, respectively. The spin-up in our models suggests that the GW waveform is that of slow chirp, where both the GW amplitude and GW frequency increase over the course of the inspiral as long as there exists a significant amount of crust. Smaller accretion torques produce a slower spin-up such that the waveform may be relatively continuous over an integration time of ${ \mathcal O }(10\,\mathrm{hr})$. Our models do not reach secular instability.

Table 3.  Model Quantities at the Final Time Step

Model	h0,f/h0,Edd	Ωf/(2π) (Hz)	${\dot{M}}_{{\rm{f}}}\ ({M}_{\odot }\,{\mathrm{yr}}^{-1})$	ρ,f	af (au)	${\dot{a}}_{{\rm{f}}}\ (\mathrm{au}\ {\mathrm{yr}}^{-1})$	ξf
M1	7.94	402	2.50	1.29	2.26 × 10−2	−7.20	2.32
M2	7.76	402	0.261	1.43	2.22 × 10−2	−0.874	19.5
M3	13.6	313	5.91 × 10−2	0.972	0.444	−2.26	1.71
M4	13.6	313	5.91 × 10−2	0.972	0.444	−0.226	17.1
M5	0.00	675	4.74	0.817	1.78 × 10−2	−4.65	5.22
M6	0.00	669	0.485	0.819	1.75 × 10−2	−0.457	54.1
M7	13.4	324	2.11 × 10−2	0.521	1.93	−1.59	1.59
M8	13.4	324	2.11 × 10−2	0.521	1.93	−0.159	15.9

Download table as:  ASCIITypeset image

The orbital separation at any point during CE evolution would be difficult to measure with electromagnetic observations due to the high optical depth of the envelope. Direct detection of a GW signal from an accreting NS within the CE may provide the only way to constrain the orbital separation between the NS and the primary's core during the inspiral. However, since this constraint relies on knowledge of the primary mass and evolutionary state, the NS mass, and the distance to the system, pre-CE electromagnetic observations would be needed to obtain useful constraints.

The orbital separation can be constrained as follows. We estimate the accretion rate from the GW strain amplitude using Equations (16) and (19). From the accretion rate, the background density gradient can be computed via Equations (22)–(41). The accretion rate as a function of the strain amplitude is plotted in panel (A) of Figure 6. During the earlier portions of the inspiral when a majority of the crust remains, higher strain amplitudes imply higher accretion rates. While the functional form of the MR15 model admits multiple solutions for _ρ given $\dot{M}$, our forward integration with an assumed stellar model selects the correct one. Using the data from panels (B) and (H) of Figures 4 and 5, we plot in the top panel of Figure 6 the nondimensional density gradient, _ρ, as a function of the strain amplitude relative to the strain amplitude for the NS at its Eddington limit. (The NS mass and the distance are used to compute h_0,Edd here.) Time is implicitly a parameter in this plot, and since the strain amplitude increases monotonically with time, both increase from left to right. Because of the early steep drop in _ρ, for a given fixed uncertainty in h₀/h_0,Edd, the best constraints should be obtained later in the evolution, where _ρ is smaller.

Zoom In Zoom Out Reset image size

Assuming a stellar model for the primary, the orbital separation can be constrained from the dimensionless density gradient. Figure 7 shows the dependence of a on _ρ; clearly, a is most sensitive to _ρ where _ρ is better determined from h₀/h_0,Edd. Note that while the connection between h₀/h_0,Edd and _ρ depends on the value of η, the dependence of a on _ρ comes from the stellar model only. We also plot a as a function of the strain amplitude in the middle panel of Figure 6, where higher strain amplitudes correspond to lower orbital separations for all of our models.

Zoom In Zoom Out Reset image size

We list final values of the strain amplitude, density gradient, and separation for our models in Table 3. Since the stars at the base of the RGB appear to undergo merger, the final separation in these cases is relatively insensitive to the mass of the primary and η. For more evolved stars, the dependence on η is essentially undetectable, but changing the mass from 12 to 20 M_⊙ increases the final separation by a factor of 4.3 and decreases the strain amplitude by about 40%. Hence, while this technique could be used to measure the amount of CE inspiral, successful application will require some type of constraint on the primary mass.

As we have seen, the inspiral parameter ξ is a measure of how rapid the inspiral is compared to the orbital period. A measurement of this parameter would provide information about the CE drag mechanism as well as the envelope structure. This is another area in which electromagnetic observations can give us very little information during the CE phase, but combining pre-CE electromagnetic constraints with the GW emission during CE can potentially distinguish between plunging and nonplunging scenarios.

Estimating ξ requires an estimate of a (Section 7.1) to determine the orbital period τ_orb and an estimate of $\dot{a}$ to determine the inspiral timescale τ_ins. In both cases, we first obtain _ρ from h₀/h_0,Edd as described in the previous section. We then use Equation (25) to determine the drag force and Equations (26)–(28) to compute $\dot{a}$. (As with the accretion rate, the functional form of the MR15 drag force model allows multiple solutions for _ρ, but an assumed stellar model lets us choose the correct one.) The inspiral parameter is plotted as a function of the strain amplitude in panel (D) of Figure 6, as a function of _ρ in the bottom panel of Figure 7, and as a function of orbital separation in Figure 8. We also tabulate values at the end of the inspiral in Table 3. As expected from the orbital evolution, the final orbital decay rates for the models at the base of the RGB are greater in magnitude than those for the more evolved models. However, the final inspiral parameter values are smaller at the tip of the RGB, though they are still greater than 1. At face value, this seems inconsistent with the apparent merging behavior of the less evolved primaries. However, ξ is more rapidly increasing for the more evolved cases. If the integration is carried beyond the termination criterion, ξ continues to increase, while for the stars at the base of the RGB, it drops sharply.

Zoom In Zoom Out Reset image size

As previously mentioned, reducing the accretion rate increases the inspiral timescale, which increases ξ, as seen in the η = 0.1 models relative to the η = 1 models. Reducing the accretion rate also reduces the strain amplitude, which shifts the model curves to the left in the ξ–h₀ plane. Early in the orbital evolution, ξ ≪ 1, and the quasi-circular approximation is invalid; we have noted that this is due to uncertainties in the outer envelope structure. During the middle stages of evolution where these uncertainties begin to affect _ρ less, the inspiral parameter depends weakly on the strain amplitude. However, at the end where the evolution resolves into stabilization as a close binary or merger, ξ undergoes a sharp uptick. In the nonmerger cases, this is particularly pronounced, while in the merger cases it is followed by a sharp drop when integrating past the termination criterion. This behavior is nearly independent of η. This abrupt change in behavior offers a (model-dependent) way to distinguish the two scenarios.

Any accreting, spinning NS with an asymmetry with respect to its rotation axis is a GW source that may be detectable in the aLIGO band. In Figure 9, we compare the GW emission properties of NS–CE systems from models M3 and M7 with LMXBs accreting at the Eddington limit, NSs undergoing fallback accretion in supernovae, a subsample of strain ratio upper limits from Meadors et al. (2017), and expected strain ratio upper limits for aLIGO, assuming a source distance of 2.8 kpc. As we have seen, even when the density gradient is steep and we make conservative assumptions about the accretion efficiency, NS–CE systems can be about 10 times louder than LMXBs like Sco X-1. However, fallback-accretion rates have been estimated to be ${\dot{M}}_{\mathrm{fall}}\sim {10}^{-4}\mbox{--}{10}^{-2}{M}_{\odot }\,{{\rm{s}}}^{-1}$ over a timescale of ${ \mathcal O }(\sim {10}^{2}\ {\rm{s}})$ (MacFadyen et al. 2001; Zhang et al. 2008; PT12). Such accretion rates can exert an accretion torque large enough to spin up the NS to the point of secular instability, β ∼ β_sec. Thus, fallback-accreting systems should be up to  ∼ 10³–10⁴ times louder than NS–CE systems. Their characteristic frequencies should reach higher values than those for X-ray binaries, but may overlap the binary range before reaching secular instability. However, since the NS–CE systems have a GW modulation timescale ranging from ∼1 month—100 yr, it should be easy to distinguish them from fallback-accretion events.

Zoom In Zoom Out Reset image size

For the aLIGO design sensitivity curve, we have taken a coherent integration time of ${\tau }_{\mathrm{coh}}=30\,\mathrm{hr}$, similar to what the Einstein@Home uses for their continuous GW searches (e.g., LIGO Scientific Collaboration et al. 2009; Abbott et al. 2016a). The coherence time is for an individual frequency such that N continuous-waveform searches within a given interval of the aLIGO band results in a total integration time of N × τ_coh. If we use its value in the range of frequencies relevant to X-ray binaries, hyperaccreting NS–CE systems should almost be detectable at 2.8 kpc once aLIGO achieves design sensitivity and the loudest merging NS–CE systems may already be detectable at 10 times this distance. If the NS has achieved the maximum quadrupole moment that it can sustain, the maximum distance at which NS–CE systems would have strain amplitudes comparable to Sco X-1 can be estimated as

$$\begin{eqnarray}&&D\approx 62\,\mathrm{kpc}{\left(\displaystyle \frac{{h}_{0}}{3.5\times {10}^{-26}}\right)}^{-1}{\left(\displaystyle \frac{P}{1.4\mathrm{ms}}\right)}^{-2}\left(\displaystyle \frac{{Q}_{\max }}{2\times {10}^{39}\ {\rm{g}}\,{\mathrm{cm}}^{2}}\right).\end{eqnarray} \tag{ 42 }$$

Thus, by the time Sco X-1 becomes detectable via gravitational waves, hyperaccreting NS–CE systems with ${ \mathcal O }(1\ \mathrm{ms})$ NS spin periods should be detectable as far away as the Magellanic Clouds. NS–CE systems with larger NS spin periods may be detectable only within the Galaxy. For example, an NS spin period of ∼4 ms would only be detectable to ∼8 kpc, i.e., the Galactic center.