Accretion spin-up of the massive component in the neutron star stripping model for short gamma-ray bursts

ABSTRACT

In this paper, we use analytical methods to study the last stages of the double neutron star (NS) system evolution. Depending on the initial masses of the components, this evolution can occur either in the framework of the merging scenario or in the NS stripping model. The main new ingredient of this work, compared with previous calculations, is accounting for accretion spin-up of the massive component. This effect leads to a significant decrease in the duration of the stable mass transfer of matter in the stripping mechanism. Within the framework of the Newtonian approximation, we determine the boundary between the merging and stripping scenarios. It is shown that this boundary weakly depends on the total mass of the system and the specific form of the NS equation of state, and is determined mainly by the initial mass ratio of the components. The stripping scenario is realized at M₂/M₁ ≲ 0.8, so it should make a large contribution to the population of gravitational wave events from NS–NS coalescing binaries that are close to us, and their accompanying short gamma-ray bursts. Nevertheless, the value obtained requires further clarification, taking into account relativistic effects, possible non-conservative mass transfer, etc.

accretion, accretion discs, gravitational waves, binaries: close, gamma-ray burst: individual: 170817A, stars: neutron

1 INTRODUCTION

Since the discovery of the pulsar PSR B1913+16 (Hulse & Taylor 1975) as a member of a binary neutron star (NS–NS) system, the interest of astrophysicists in such systems has only increased. This is primarily due to the fact that they are natural laboratories to test general relativity (GR) effects. In particular, the analysis of the radio pulses from PSR B1913+16 has indirectly revealed gravitational wave (GW) emission from two NSs spiralling towards each other (Taylor & Weisberg 1989).

The binary NS inspirals are also of interest as sources of short gamma-ray bursts (GRBs; Blinnikov et al. 1984; Eichler et al. 1989; Narayan, Paczynski & Piran 1992; Nakar 2007). This process is usually described in the merging model (e.g. Faber & Rasio 2012), where two NSs approach each other due to the loss of angular momentum to emit GW and finally merge into a single object, a supramassive neutron star (SMNS) or a black hole (BH). However, there is an alternative to this mechanism proposed by Blinnikov et al. (1984), namely the stripping model. When two NSs spiral inward, the low-mass component first fills its Roche lobe and begins to flow onto its more massive companion. In the end, the low-mass NS (LMNS) reaches the value corresponding to the minimum possible NS mass and explodes, actually producing GRB (see also Blinnikov et al. 1990).

However, for a number of reasons, the stripping model has been neglected for many years. The first reason was connected with the lack of observational data at that time for the NS–NS systems with high mass asymmetry of the components (Thorsett & Chakrabarty 1999), so necessary for the stripping mechanism. There were also doubts about the existence of stable mass transfer during the stripping of the low-mass component (e.g. Lai, Rasio & Shapiro 1994). However, the main reason was the weakness of the GRB: the observed short GRBs with known redshifts had, by orders of magnitude, more energy (Fong et al. 2015; Lien et al. 2016) than that predicted by the stripping model.

Interest in the stripping model was revived on 2017 August 17. After the historical joint detection of the GW event GW170817 and the accompanying GRB170817A, the connection between short GRBs and NS–NS coalescences has been reliably confirmed (Abbott et al. 2017a,b). In our paper Blinnikov et al. (2021), we have shown that many observational properties of GRB170817A, which turned out to be peculiar compared with other short GRBs, are naturally explained in the context of the stripping model. These properties are the anomalously small total isotropic energy of the GRB, the large ejected mass of the red kilonova (Siegel 2019) and its spherical distribution (Sneppen et al. 2023), and the long time delay of |$t_{\mathrm{str}}\approx 1.7 \, \text{s}$| between the GW signal peak and the GRB registration (Clark & Eardley 1977), corresponding to the duration of the stable mass transfer during the stripping of the LMNS. In subsequent work Blinnikov et al. (2022), it was shown that the absence of the stable mass transfer observed in some hydrodynamic simulations may be caused by incorrect selection of the initial conditions [see also Dan et al. 2011 for the white dwarf (WD) systems].

In this paper, we examine the properties of the stripping model using the analytical approach of Clark & Eardley (1977). Compared with the previous calculations (see also Jaranowski & Krolak 1992; Imshennik & Popov 1998; Imshennik 2008), we take into account the effect of accretion spin-up of the massive component obtained earlier in our numerical simulation (Blinnikov et al. 2022). Moreover, we consider the corotation of the LMNS and the tidal spin-down of the massive NS (MNS). This allows us to determine the boundary between the merging and stripping scenarios, as well as the duration of the stable mass transfer t_str (or the stripping time) – the most important dynamical parameter of the NS stripping model.

The outline of this paper is as follows. In Section 2, we describe the problem statement and derive the basic equations that determine the NS–NS system evolution in the stripping model. A discussion of the results is presented in Section 3. Section 4 summarizes our main conclusions and outlines the prospects for further progress.

2 PROBLEM STATEMENT

Let us consider the last stages of the double NS system evolution in more detail, referring to Fig. 1. Two NSs with approximately equal masses approach each other, due to the loss of the total angular momentum of the system J_tot to emit GW, and finally merge (indicated by the arrow to the right). If the binary system has a high enough mass asymmetry of its components (arrow to the left), the LMNS (M₂ in the figure) first fills its Roche lobe and begins to flow onto the MNS (M₁) through the inner Lagrangian point L₁. During accretion, the asymmetry of the system only grows and the components recede each other. The mass transfer lasts on a relatively long time-scale, determined by the rate of the angular momentum loss radiated away by GW. At the same time, the orbital angular momentum of the system J_orb is partially transferred to the total rotational (or spin) angular momentum of the MNS and the accretion disc J₁. The disc forms when the minimum distance R_m of approach of the stream to the centre of the accretor is larger than the equatorial radius of the MNS R₁ (Lubow & Shu 1975; Kramarev & Yudin 2023). Because of the tidal spin-down effect, some part of the rotational angular momentum J₁ of the MNS can be transferred back to the orbital angular momentum J_orb. By increasing the distance between the components, the LMNS, being in the corotation (Yudin, Razinkova & Blinnikov 2020), also spins down and transfers a part of its rotational angular momentum J₂ to the orbital one. If the total mass of the system M_tot is greater than the maximum NS mass M_max (see Appendix  B), then at some point the MNS collapses into the BH. In any case, the result of the stripping mechanism is the same (see Fig. 1): once the donor reaches some very low mass value (∼0.2 M_⊙), the stability of the mass transfer is lost and the remnant M₂ is absorbed by M₁ on the hydrodynamic time-scale. When the LMNS reaches the minimum possible NS mass M_min ∼ 0.1 M_⊙, it loses hydrodynamic stability and explodes.

Figure 1.

The last stages of the double NS system evolution depending on the initial masses of the components.

2.1 General assumptions

Let us move to the mathematical description of the problem statement. As discussed above, the total angular momentum of the system J_tot consists of the orbital angular momentum J_orb and the rotational angular momenta of the components J₁ and J₂. It is lost due to the GW emission and the magneto-dipole radiation of the rotating MNS. The magneto-dipole losses, as we show in Appendix  A, can be neglected, so the following equality takes place (see also Kremer, Sepinsky & Kalogera 2015):

$$\begin{eqnarray} \dot{J}_{\mathrm{GW}}=\dot{J}_{\mathrm{orb}}+\dot{J}_1+\dot{J}_2. \end{eqnarray}$$

(1)

According to Kowalska et al. (2011), the orbital eccentricity of most NS–NS systems at the final stages is negligibly small, which, in particular, corresponds to the results of GW170817 and GW190425 data processing (Lenon, Nitz & Brown 2020). Therefore, we consider the binary to be in a circular Keplerian orbit with the orbital rotational frequency

$$\begin{eqnarray} \Omega _{\mathrm{orb}}=\sqrt{\frac{GM_{\mathrm{tot}}}{a^3}}, \end{eqnarray}$$

(2)

where a is the orbital separation, M_tot = M₁ + M₂ is the total mass of the components. The mass transfer is also assumed to be conservative. The orbital angular momentum of the system is defined by the formula

$$\begin{eqnarray} J_{\mathrm{orb}}=\frac{M_1 M_2}{M_{\mathrm{tot}}}a^2 \Omega _{\mathrm{orb}}. \end{eqnarray}$$

(3)

The rate of the orbital angular momentum loss to emit GW is given by (e.g. Paczyński 1967)

$$\begin{eqnarray} \dot{J}_{\mathrm{GW}}=-\frac{32}{5}\frac{G}{c^5}\frac{M_{1}^2 M_{2}^2}{M_{\mathrm{tot}}^2}a^4 \Omega _{\mathrm{orb}}^5. \end{eqnarray}$$

(4)

In addition to the definition of J₁ and J₂, our main equation (1) must be supplemented by another equation. Prior to the beginning of the mass transfer, this is the condition of constancy of the component masses. Then the low-mass component fills its Roche lobe and the stable mass transfer begins. The closure equation in this case is R₂ = R_R, where the equatorial radius of the LMNS, R₂ = R₂(M₂, J₂), is determined by the NS equation of state (EoS), as well as its spin angular momentum (see Appendix  B). The effective radius of the Roche lobe is parametrized in the standard way (Eggleton 1983):

$$\begin{eqnarray} R_{\mathrm{R}}=a f(q^{\prime }), \,\, f(q^{\prime })=\frac{0.49 (q^{\prime })^{2/3}}{0.6(q^{\prime })^{2/3}+\ln \big [1+(q^{\prime })^{1/3}\big ]}, \end{eqnarray}$$

(5)

where q′ = M₂/M₁. Note that we will concurrently use relation q = M₂/M_tot, the ratio of the donor mass to the total mass of the system.

2.2 Accounting for accretion spin-up of the massive component

We assume that before the mass-transfer process, both NSs are synchronized due to the tidal effects and corotate, i.e. J_{1, 2} = I_{1, 2}Ω_orb, where I_{1, 2} are the moments of inertia of the accretor and donor. After the beginning of the stable mass transfer, the rotational angular momentum of the massive component changes due to the effects of accretion spin-up and tidal spin-down:

$$\begin{eqnarray} \dot{J}_1=\dot{J}_\mathrm{acc}+\dot{J}_\mathrm{tid}. \end{eqnarray}$$

(6)

For the |$\dot{J}_\mathrm{acc}$| term we use the following parametrization:

$$\begin{eqnarray} \dot{J}_\mathrm{acc}=-\dot{M}_2 \mathfrak {j}(q,r_1)a^2\Omega _{\mathrm{orb}}, \end{eqnarray}$$

(7)

where |$\mathfrak {j}$| is the specific angular momentum of accreting matter in orbital units, r₁ is the dimensionless stopping radius (see equation 8 below). As discussed above, during the stripping of the low-mass component, two modes of accretion can take place (Lubow & Shu 1975). In the first case, the accretion stream hits the surface of the accretor with the equatorial radius R₁ and the so-called direct impact accretion takes place. If the minimum distance R_m for which the stream approaches the massive component is larger than the equatorial radius of the accretor R₁, an accretion disc with an outer radius of R_d is formed. The corresponding stopping radius for each mode is equal to

$$\begin{eqnarray} r_1 = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}R_1/a, \,\, R_1\geqslant R_{\mathrm{m}}, \\ R_{\mathrm{d}}/a, \,\, R_1\lt R_{\mathrm{m}}. \end{array}\right. \end{eqnarray}$$

(8)

The detailed calculation and approximation of |$\mathfrak {j}(q,r_1)$|⁠, as well as R_m(q) and R_d(q), are described in our previous paper Kramarev & Yudin (2023). We only emphasize that all our approximations have accuracy better than 1 per cent, but are obtained in the framework of the Newtonian limit and therefore require further refinement to take into account the GR effects.

2.3 Accounting for tidal spin-down of the MNS

Nowadays, several parametric approaches are known to account for the tidal spin-down effect in analytical calculations of the binary system evolution (e.g. Zahn 1977; Eggleton, Kiseleva & Hut 1998). However, the large ambiguity of the input parameters leads to a significant (by orders of magnitude) change in the contribution of this effect (see also Kushnir et al. 2017). To qualitatively account for the tidal spin-down of the MNS, we first consider two limiting cases. In one case, the massive component spins up during accretion according to the formula (7), and there is no tidal spin-down, i.e. |$\dot{J}_\mathrm{tid}=0$|⁠. In the other limit, the MNS, like the LMNS, rotates synchronously both before and after the beginning of the mass transfer, i.e. J₁ = I₁Ω_orb. Accretion spin-up of the massive component is completely absent in this case. A comparison of calculations is given in Fig. 3.

To quantitatively account for the tidal spin-down effect, we exploit the approach used by Marsh, Nelemans & Steeghs (2004), Gokhale, Peng & Frank (2007), and Kremer et al. (2015) to calculate the evolution of WD–WD binaries. In our case, the term expressing tidal spin-down of the MNS can be parametrized as

$$\begin{eqnarray} \dot{J}_\mathrm{tid}=\frac{I_{1}}{\tau _{\mathrm{syn}}}(\Omega _{\mathrm{orb}}-\Omega _1), \end{eqnarray}$$

(9)

where Ω₁ = J₁/I₁ is the massive component spin and τ_syn is its characteristic tidal synchronization time-scale. Although the value of τ_syn can vary within wide limits, we will show below (see Section 3.3) that in this approach the tidal spin-down effect appears to be secondary to the accretion spin-up.

2.4 Corotation of the LMNS

As was shown in Yudin et al. (2020), even if there is an initial angular momentum, the LMNS loses it quickly enough during the mass transfer. Therefore, we assume that the low-mass component always rotates synchronously both before and after the beginning of accretion, i.e. J₂ = I₂(M₂, J₂)Ω_orb. So the third term on the right-hand side of equation (1) can be written as

$$\begin{eqnarray} \dot{J}_2=\beta _2\bigg [I_2 \dot{\Omega }_{\mathrm{orb}}+\Omega _{\mathrm{orb}}\dot{M}_2 \left(\frac{\partial I_2}{\partial M_2}\right)_{\!J_2}\!\!\bigg ], \end{eqnarray}$$

(10)

where we introduce the notation |$\beta _2=\Big [1{-}\Omega _{\mathrm{orb}}\left(\frac{\partial I_2}{\partial J_2}\right)_{\!M_2}\!\!\Big ]^{-1}$| for compactness. The calculation procedure of the moment of inertia I and the equatorial radius R of the rotating NS as functions of M and J is described in Appendix  B.

2.5 The stability criterion for the mass transfer

For stable mass transfer, the size of the Roche lobe must grow faster than the radius of the low-mass component. In accordance with equations (1)–(10), the stability criterion can be written as

$$\begin{eqnarray} \frac{\mathrm{ d} \ln R_2}{\mathrm{ d} \ln M_2}\geqslant \frac{\mathrm{ d} \ln f}{\mathrm{ d} \ln q }- 2\frac{1{-}2q{-}\mathfrak {j}(q,r_1){+}\frac{\beta _2}{ a^2}\left(\frac{\partial I_2}{\partial M_2}\right)_{\!J_2}}{1{-}q{-}\frac{3\beta _2}{q}\frac{I_2}{ M_2 a^2}}. \end{eqnarray}$$

(11)

Analysis of this formula shows that the stability of the mass-transfer process is determined mainly by the NS mass–radius relation in the low-mass range, and by the formula for the specific angular momentum of matter going to spin-up the accretor.

The derived criterion for the stability of the mass transfer determines the boundary between the merging and stripping scenarios, as well as the duration of the stable mass transfer t_str – the most important dynamical parameter of the stripping scenario, corresponding to the time delay between the GW signal peak and the GRB detection (Blinnikov et al. 2021, 2022).

2.6 The gravitational wave and neutrino luminosity

In the light of the GW170817 and GW190425 registration and the further search for the NS–NS inspirals during the fourth observing run (O4) of the LIGO–Virgo–KAGRA GW detectors (Colombo et al. 2022), we are also interested in the GW luminosity in the stripping mechanism. For two point masses moving in circular orbits the GW luminosity is defined by the well-known formula (e.g. Landau & Lifshitz 1975)

$$\begin{eqnarray} L_{\mathrm{GW}}=\frac{32}{5}\frac{G^4}{c^5}\frac{M_{\mathrm{tot}}^{5} q^{2} (1{-}q)^{2}}{a^5}. \end{eqnarray}$$

(12)

With the beginning of the stable mass transfer, some part of the accreting matter energy is radiated, following Clark & Eardley (1977), in the neutrino channel. We estimate the neutrino luminosity using the formula

$$\begin{eqnarray} L_{\mathrm{\nu }}=\frac{G M_1 \dot{M_1}}{R_{\mathrm{\nu }}}, \end{eqnarray}$$

(13)

where R_ν = r₁a is the equatorial radius of the MNS or the radius of the accretion disc (after disc formation). In the case of direct accretion onto the BH, the neutrino luminosity is obviously equal to zero. We note once again that the formula (13) is an estimate. Moreover, some part of the energy is also radiated in the electromagnetic channel, which is important in the context of the search for precursor activity of GRBs (Koshut et al. 1995; Mereminskiy et al. 2022) and the problem of non-conservative mass transfer.

2.7 Consideration of the GR effects

The contribution of the GR effects to our problem can be essential in two hypostases. First, it influences the system parameters: the orbital rotational frequency (Thorne & Hartle 1985) and the effective radius of the Roche lobe of the low-mass component (Ratkovic, Prakash & Lattimer 2005). From this point of view, the impact of the relativistic effects is significant only at the moment of the closest approach of the components |$a\approx 40\, \mathrm{km}$|⁠, when the stripping of matter from the surface of the LMNS begins. For a total mass value of about M_tot = 2–3 M_⊙, we obtain the contribution of the relativistic effects on the order of |$2GM_{\mathrm{tot}}/a c^2 \sim 10 {{\, \rm per\, cent}}$|⁠. It is clear that with an increase in the orbital separation a, the impact of the GR effects on the system parameters decreases.

On the other hand, as discussed in Kramarev & Yudin (2023), the relativistic effects also impact the dynamics of the accretion stream near the surface of the MNS and, consequently, the accretor spin-up process at about 10 per cent. After the collapse of the MNS into the BH, the GR effects must play a major role in the accretion (e.g. Frolov & Novikov 1998). An accurate calculation of the orbital angular momentum transfer to the spin angular momentum of the BH (or the MNS) during accretion would require us to solve the relativistic three-body problem. But instead, we introduce some effective radius of the BH at which the accreting matter stops and transfers its angular momentum, defined by the approximation formula |$\mathfrak {j}(q,r_1)$|⁠, to the BH. As such a radius we can choose, for example, the radius of the equatorial innermost stable circular orbit of the BH,

$$\begin{eqnarray} R_{\mathrm{ISCO}}= \frac{R_\mathrm{ g}}{2}\left(3+Z_2+[(3-Z_1)(3+Z_1+2Z_2)]^{1/2}\right), \end{eqnarray}$$

(14)

including the non-rotating BH radius 3R_g and the limiting rotation radius 4.5R_g, where R_g = 2GM/c² is the Schwarzschild radius of the BH with mass M. Following Frolov & Novikov (1998), we determine some auxiliary functions:

$$\begin{eqnarray} Z_1=1+\left(1-l^2\right)^{1/3}\left[\left(1+l\right)^{1/3}\! \! \!+\left(1-l\right)^{1/3}\right], \end{eqnarray}$$

(15)

$$\begin{eqnarray} Z_2=\left(3 l^2+Z_1^2\right)^{1/2}, \end{eqnarray}$$

(16)

where l = 2J/R_gMc is the specific angular momentum of the BH varying within l ∈ [−1; +1]. One can also choose as such a radius the minimum of the periastrons of all parabolic orbits, which in our notation is written as

$$\begin{eqnarray} R_\mathrm{ b}=R_\mathrm{ g}\left(1+\frac{l}{2}+\sqrt{1+l}\right). \end{eqnarray}$$

(17)

By carrying out a series of NS stripping calculations with different formulae for the effective radius of the massive component after its collapse into the BH, we find that the stripping time t_str – the main parameter of the NS stripping model – weakly depends on the particular type of the formula (see Section 3.4). This means that the approach described above is acceptable for our problem. Unless otherwise specified, we will use 3R_g as the effective radius of the BH hereafter by default.

3 RESULTS

3.1 The boundary between the merging and stripping scenarios

The derived criterion (11) for the stability of the mass transfer, as discussed above, determines the boundary between the merging and stripping scenarios. Before proceeding to the calculations with the specific NS EoS, let us pay attention to the following circumstance. A distinctive feature of the most modern EoSs, satisfying the observational data (e.g. Lattimer & Prakash 2001; Greif et al. 2020; Raaijmakers et al. 2020), is that with NS mass variations from 1 to 2 M_⊙, the NS radius changes by 1–2 km. It turns out that in the moderate mass range, the NS EoS can be described with good accuracy by the polytrope P(ρ) = Kρ^1 + 1/n with index n = 1, for which strictly R(M) = const (e.g. Zel’dovich, Blinnikov & Shakura 1981). It follows that R₁ = R₂ = af(q′). We also make use of the fact that the low-mass component, being in the corotation, spins relatively slowly. Then, in the expression (10) for the derivative of the LMNS spin angular momentum, we can put β₂ = 1.

In view of the above, the stability criterion for the mass transfer (11) can be rewritten as

$$\begin{eqnarray} 0\geqslant \frac{\mathrm{ d} \ln f}{\mathrm{ d} \ln q }- 2\frac{1{-}2q{-}\mathfrak {j}(q,f)+I^{*}_2 f^2}{1{-}q{-}\frac{3 I^{*}_2 f^{2}}{q}}, \end{eqnarray}$$

(18)

where I* = I/MR² is the dimensionless moment of inertia (see Appendix  B). It is easy to notice that the resulting approximation (18) is determined by the mass ratio of the components q (or q′) and weakly depends on the total mass M_tot. This is what we expect from calculations with the different EoSs in agreement with the exact expression for the stability criterion (11).

Fig. 2 shows the boundary between the merging and stripping scenarios for the BSk22, BSk24, BSk25, and BSk26 EoSs (see Pearson et al. 2018). The horizontal axis represents the total mass of the system, and the vertical axis represents the initial mass ratios of the components q and q′. The parameter area M_tot–q above the particular line corresponds to the further merger of the components, and the area below refers to the stripping process of the LMNS. The dotted horizontal lines illustrate the boundary obtained according to the approximate criterion (18) for different values of the dimensionless moment of inertia: I*(M₂ = 1 M_⊙) = 0.25 and I*(M₂ = 1.6 M_⊙) = 0.35. The dotted hyperbolas M₂ = 1 M_⊙ and M₂ = 1.6 M_⊙ determine the limits of the applicability of our approximation. The slight discrepancy between the prediction of the approximate criterion (18) and the exact numerical calculations is due to the contribution of the logarithmic derivative d ln R₂/d ln M₂, which is not strictly equal to zero. Despite this, the predictions of the formula (18) are qualitatively correct: the position of the boundary, in fact, weakly depends on the total mass M_tot and is mainly determined by the mass ratio of the components q′. As can be seen from Fig. 2, the boundary value of the mass ratio is approximately |$q^{\prime }_{\mathrm{str}}(M_{\mathrm{tot}})=0.8\pm 0.05$|⁠. It is reasonable to compare this value with the constraints on the component masses obtained in the LIGO–Virgo post-processing data analysis for the GW170817 (Abbott et al. 2019) and GW190425 (Abbott et al. 2020) events. The initial masses for the case of low component spins, corresponding to our corotation case, are represented by the grey regions in Fig. 2. It can be seen that both detected events could be the result of the NS merging and stripping scenarios.

Figure 2.

The boundary between the merging and stripping scenarios for the different NS EoS. See text for details.

How would the position of the boundary between the scenarios in Fig. 2 change if tidal spin-down of the MNS is taken into account? According to the formula (9), the contribution of this effect to the change of the MNS angular momentum is proportional to the difference between its own and the orbital rotational frequency, Ω₁ and Ω_orb. Therefore, in the case of the MNS corotation, where before the beginning of the mass transfer the frequencies are equal to each other (Ω₁ = Ω_orb), the influence of the mentioned effect must be negligible. However, in the case where the initial rotational frequency of the MNS is more than the orbital one (Ω₁ > Ω_orb), with the beginning of the mass transfer the spin angular momentum of the massive component J₁ must effectively transfer to the orbital one J_orb. Therefore, in the Newtonian approximation, the case of the initial component corotation illustrated in Fig. 2 gives an absolute lower limit on the boundary |$q^{\prime }_{\mathrm{str}}$| (or q_str) between the merging and stripping scenarios. With an increase in the initial rotational frequencies Ω_{1, 2} compared with Ω_orb, the boundary value |$q^{\prime }_{\mathrm{str}}$| between the scenarios also increases, due to the tidal effect and the growth of |$I_2^{*}$| in the approximate criterion (18).

3.2 The importance of accretion spin-up of the massive component

Now let us examine the main element of our work – accounting for accretion spin-up of the massive component. In most previous works (Clark & Eardley 1977; Jaranowski & Krolak 1992; Imshennik & Popov 1998; Imshennik 2008), the calculation of the stripping process was performed without taking into account the spin-up of the accretor, i.e. the authors solved a much simpler equation instead of equation (1):

$$\begin{eqnarray} \dot{J}_{\mathrm{GW}}=\dot{J}_{\mathrm{orb}}. \end{eqnarray}$$

(19)

However, our hydrodynamic simulations (Blinnikov et al. 2022) with Phantom code (Price et al. 2018) showed that the MNS spins up during the stable mass transfer. From general considerations it is clear that accounting for accretion spin-up must lead to a greater loss of the orbital angular momentum of the system and, consequently, to its faster rate of evolution. This is most evident from the comparison of panels (a) and (b) in Fig. 3.

Figure 3.

The stripping time as a function of the initial component masses for the BSk22 EoS. The upper left-hand panel (a) illustrates the results of the simplest calculations, made according to formula (19), and the upper right-hand panel (b) calculations account for the accretion spin-up of the massive component, as well as the corotation of the LMNS, in accordance with equation (1). The lower left-hand panel (c) calculations take into consideration the tidal spin-down of the MNS, assuming its corotation both before and after the beginning of the mass transfer. The colour bars of panels (a) and (c) are consistent with each other. We also add the mass ranges of the GW170817 source: the dark areas correspond to the case of the small initial component spins, and the light areas correspond to the large ones (see Abbott et al. 2019). The blue, red, and green circles correspond to the initial masses M₁ = 1.6 M_⊙ and M₂ = 1.15 M_⊙, for which the evolution of the orbital separation, component masses, and luminosities with time are given in the lower right-hand panel (d). The drop region of the neutrino luminosity for the green and red lines corresponds to the direct impact accretion onto the BH followed by the formation of an accretion disc. See text for a discussion of all details.

Panel (d) in Fig. 3 illustrates the evolution of the basic quantities of the NS–NS system with initial masses M₁ = 1.6 M_⊙ and M₂ = 1.15 M_⊙, corresponding to the mass range of the GW170817 source (Abbott et al. 2019). The zero-point in time corresponds to the closest approach of the components and the beginning of the stripping process. The calculations were carried out up to the end of the stable mass transfer according to the criterion (11). The evolution of all quantities at t ≤ 0.5 s is plotted on a linear scale, and at t > 0.5 s on a logarithmic scale. The red dash–dotted line corresponds to the calculation in the Clark & Eardley (1977) approach, without taking into account accretion spin-up, tidal effects, etc. The blue solid line is the calculation according to the formula (1), where we account for accretion spin-up of the massive component and corotation of the LMNS. One can see how much accounting for the spin-up effect influences the NS stripping process: the stripping time t_str decreases by an order of magnitude! Hence it follows that accretion spin-up of the massive component must necessarily be taken into account in all future calculations of the NS stripping process.

3.3 The contribution of the MNS tidal spin-down effect

During tidal spin-down of the MNS, part of its spin angular momentum transfers back to the orbital angular momentum. Therefore, accounting for the tidal spin-down effect must inevitably increase the stripping time. As discussed above, accounting for this effect correctly is quite difficult. So, we first consider the upper limit corresponding to the tidal synchronization of the MNS, both before and after the beginning of the stripping process. For this purpose, in the angular momentum evolution equation (1) we replace the spin-up formula (6) with (10) (with a corresponding index substitution 2 → 1). The dependence of the stripping time on the initial masses for the case of the MNS corotation is shown in panel (c) in Fig. 3. Inside the painted part, there are two regions corresponding to the compact object type at the MNS site: SMNS (when M_tot < M_max ≈ 2.4 M_⊙) or the BH. This circumstance significantly influences the total stripping time, because after the collapse of the SMNS into the BH (BH formation) the tidal spin-down ceases. That is why at M_tot ≳ M_max the stripping time |$t_{\mathrm{str}}\approx 4 \, \mathrm{s}$|⁠, obtained accounting for the tidal spin-down effect, is only slightly longer than the stripping time corresponding to the calculation based on the simplest equation (19) (compare panels a and c).

The detailed calculation for the initial masses M₁ = 1.6 M_⊙ and M₂ = 1.15 M_⊙ is shown by the green dashed line in panel (d) in Fig. 3. The green and blue lines coincide until the beginning of the mass transfer, as it should be. As noted above, the total stripping time for the case of the MNS corotation turns out to be slightly longer than the corresponding time obtained from the calculation in the Clark & Eardley (1977) approach. This results from the fact that, compared with the case (19), with the growth of the separation a between the components of the system, the orbital angular momentum increases due to the decrease in the spin angular momenta of the components, since J_{1, 2} ∼ Ω_orb ∼ a^−3/2. It is clear that the true line must lie between the blue and green lines corresponding to the two limiting cases.

We also quantitatively account for the tidal spin-down effect according to formula (9). The results of calculations for different values of τ_syn are shown in Fig. 4. It can be seen that even at extremely low values of τ_syn the stripping time t_str changes by only a few tens of per cent. This is due to the fact that most of the stripping time there is accretion onto the BH that is not affected by the considered effect. Thus, within this approach, the tidal spin-down effect is secondary to the accretion spin-up.

Figure 4.

The evolution of the orbital separation, component masses, and luminosities for the NS–NS system after the beginning of the mass transfer for the different tidal synchronization time-scales τ_syn. The initial component masses are M₁ = 1.6 M_⊙ and M₂ = 1.15 M_⊙, and the BSk22 EoS was used. The break in the functions a(t), M₁(t), M₂(t), and log L_GW(t) at small τ_syn corresponds to the derivative discontinuity during the collapse of the MNS into the BH. See text for details.

3.4 The influence of the relativistic effects

As discussed in Section 2.7, the main manifestation of the relativistic effects is the change of the accreting matter stream dynamics near the surface of the MNS. It is supposed that the GR effects must also be of primary importance when considering the dynamics of the accretion matter falling onto the BH. In order to estimate the influence of the GR effects in calculating the specific angular momentum of accreting matter |$\mathfrak {j}(q,r_1)$| going to spin-up the BH in accordance with equation (7), we artificially introduce the effective BH stopping radius R₁, which is 3R_g, 4.5R_g, R_ISCO, or R_b (see formulae 14–17). Fig. 5 shows how the stripping time – the main dynamical parameter of our problem – changes for each radius formula. The small, within hundredths of a second, variations of t_str allow us to argue that the influence of GR effects is not as important as it might seem at first glance. However, it does not eliminate the need to solve the general relativistic problem of the direct impact accretion onto the MNS or BH.

Figure 5.

The stripping time as a function of the initial mass of the MNS for the different effective stopping radii of the BH. The initial mass of the low-mass component is M₂ = 1 M_⊙. The area to the right of M₁ ≈ 1.6 M_⊙ corresponds to the collapse of the MNS into the BH during the LMNS stripping process. The calculations are performed for the BSk22 EoS.

4 CONCLUSIONS

The NS stripping model has experienced a second birth after the joint detection of the GW signal GW170817 and the accompanying GRB170817A (Blinnikov et al. 2021, 2022). In this work, we have calculated the long-term evolution of angular momentum in asymmetric NS–NS binaries undergoing the stripping process. We have advanced the analytical approach of Clark & Eardley (1977), mainly accounting for the accretion spin-up of the massive component (Kramarev & Yudin 2023). It is shown that this effect leads to a significant, more than an order of magnitude, decrease in the stripping time t_str (see Fig. 3), which in the paradigm of the stripping mechanism corresponds to the time delay between the loss of the GW signal and the GRB registration (⁠|$t_{\mathrm{str}} \approx 1.7 \, \mathrm{s}$| for event GW170817–GRB170817A, according to Abbott et al. 2017b). Therefore, the spin-up of the massive component must be necessarily taken into account in all subsequent calculations of the LMNS stripping process.

Moreover, we have considered the MNS tidal spin-down and the LMNS corotation. Because of these effects, some part of the rotational angular momentum of components is transferred back to the orbital angular momentum, which leads to only a slight increase in the stripping time. In this context, we have also derived useful approximate formulae to calculate the moment of inertia and equatorial radius of a rotating NS. As shown in Fig. B1, our formulae produce a good fit to the GR numerical calculations with different EoSs.

Figure B1.

The mass–radius relations of the various rotating configurations for the BSk-type EoSs. The different markers (crosses, asterisks, circles, etc.) illustrate the GR numerical calculations and the lines show the results of our approximation (B2)–(B3). See text for details.

We summarize different effects and their impact on the stripping process in Table 1. In the near future, we intend to examine the influence of the (cold) NS EoS in the low-mass range (Sotani et al. 2014), as we expect it may significantly increase the stripping time. Additionally, we are preparing calculation of the direct impact accretion in various post-Newtonian approaches (Thorne & Hartle 1985; Kaplan, Nichols & Thorne 2009) to account for the GR effects.

Table 1.

Upper part: effects considered in this work and their impact on the stripping time t_str. Lower part: unaccounted effects and their potential influence.

The effect	The impact of considered effect
The MNS spin-down via electromagnetic radiation	This effect is negligible even when considering rapidly rotating magnetar (see Appendix  A).
Tidal spin-down of the MNS	It insignificantly increases tstr even for extremely low values of τsyn (see Fig. 4).
The LMNS corotation	Spin of the LMNS is transferred back to the orbital angular momentum during stripping resulting in small increase in tstr.
Accretion spin-up of the MNS	This effect leads to a significant decrease in tstr (compare panels a and b in Fig. 3).
GR correction to the system parameters	We discuss at the beginning of Section 2.7 that the orbital rotational frequency and radius of the Roche lobe appear to vary by 10 per cent when mass transfer begins. It can drastically change the boundary between scenarios.
Relativistic modification of direct impact accretion	As it was shown, for example, in Appendix  B of Kramarev & Yudin (2023), GR effects can significantly change dynamics of accretion near the MNS surface and, respectively, the accretion spin-up process.
Uncertainty in the (cold) NS EoS	According to the criterion (11), the stability of the mass transfer depends on the logarithmic derivative of the function R(M), so tstr must be very sensitive to the NS EoS in the low-mass range.
Tidal heating of the LMNS crust	Because of the tidal friction of the extended outer LMNS crust during the stripping process, its R(M) relation can drastically change in comparison with the cold configuration.
Non-conservative mass transfer	As we note in Section 2.6, some part of the liberated gravitational potential energy is radiated in the electromagnetic channel. It can lead to a super-Eddington regime of accretion and, correspondingly, non-conservative evolution of a NS–NS system.

Finally, let us discuss the other important result, related to the determination of the boundary between the NS merging and stripping scenarios. We have shown that this boundary weakly depends on the total mass of the system and the concrete form of the NS EoS (see Fig. 2) and is determined mainly by the initial mass ratio of the components, which in the considered approach turned out to be |$q^{\prime }_{\mathrm{str}} = M_2/M_1 \approx 0.8$|⁠. Based on this high value, about a quarter of the observed galactic NS–NS binaries with known masses (Farrow, Zhu & Thrane 2019) must finish their evolution in accordance with the stripping model. Therefore, the stripping mechanism must contribute substantially to the total population of close-to-us (⁠|$d \lesssim 200 \, \mathrm{Mpc}$|⁠) short GRBs, and their corresponding GW sources can be registered during the O4 observing run (Colombo et al. 2022). However, we expect that the subsequent careful accounting for the GR effects, non-conservative mass transfer, and other possible factors (see Table 1) will decrease the resulting value of the critical binary mass ratio. Even if the fraction of the stripping mechanism in the total population of short GRBs is small, we also expect its large contribution to the production of cosmic heavy elements (Panov & Yudin 2020; Yip et al. 2022) due to the relatively large ejected mass (⁠|$M_{\mathrm{min}}\sim 0.1 \, \mathrm{ M}_{\mathrm{\odot }}$|⁠).

ACKNOWLEDGEMENTS

NK is grateful to the RSF 19-12-00229 grant for support. The authors acknowledge Konstantin Manukovsky for calculations of the NS rotating configurations with rns code. We also thank the reviewer for useful comments.

DATA AVAILABILITY

Data generated from computations are reported in the body of the paper. Additional data can be made available upon reasonable request.

APPENDIX A: THE CONTRIBUTION OF THE MNS ELECTROMAGNETIC EMISSION

Let us obtain an upper bound for the effect related to the magneto-dipole angular momentum loss of the MNS. Let |$m=B R_1^3/2$| be the magnetic dipole moment of the NS, where B is the magnetic field at the poles. We estimate the angular momentum loss as

$$\begin{eqnarray} \dot{J}_{\mathrm{m}}\approx \frac{\dot{E}_{\mathrm{m}}}{\Omega _1}, \end{eqnarray}$$

(A1)

where the energy emission is described by the classical formula (e.g. Landau & Lifshitz 1975)

$$\begin{eqnarray} \dot{E}_{\mathrm{m}}=-\frac{2}{3 c^2}m^2 \Omega _1^{4}. \end{eqnarray}$$

(A2)

The rotational frequency of the NS is bounded above by the Keplerian limit (e.g. Tassoul 1978):

$$\begin{eqnarray} \Omega _1 \le \Omega _{\mathrm{K}} \approx \sqrt{\pi G \bar{\rho }}, \end{eqnarray}$$

(A3)

where |$\bar{\rho }=3 M_1/4\pi R_1^{3}$| is the average density of the MNS.

It is reasonable to compare the obtained expression |$\dot{J}_{\mathrm{m}}$| for the magneto-dipole losses with the formula (7) for the angular momentum transferring to the accretion spin-up of the MNS. Taking everywhere |$\Omega _1=\sqrt{\pi G \bar{\rho }}$|⁠, we can write the upper bound for the momenta relation:

$$\begin{eqnarray} \frac{\dot{J}_{\mathrm{m}}}{\dot{J}_{\mathrm{acc}}} \lesssim 0.5 \times {10}^{-9} \! \sqrt{\frac{r_1 (1{-}q)^3}{q^2} r_1} \left(\frac{\mathfrak {j}}{0.5}\right)^{-1} \!\! \frac{t_{\mathrm{str}}}{1 \, \mathrm{s}} \left(\frac{B}{10^{12} \, \mathrm{G}}\right)^{2} \!\!\! \frac{R_1}{10 \, \mathrm{km}}.\\ \end{eqnarray}$$

(A4)

Thus, even in the case of the extremely rotating magnetar with B ∼ 10¹⁵ G, the magneto-dipole angular momentum loss can be neglected with great accuracy, since |$\dot{J}_{\mathrm{m}}/\dot{J}_{\mathrm{acc}} \sim {10}^{-3}$|⁠.

APPENDIX B: THE APPROXIMATE FORMULAE FOR THE MOMENT OF INERTIA AND EQUATORIAL RADIUS OF A ROTATING NS

First of all, let us obtain a consistent description of the R(M) lines of the rotating NSs from the known case of the non-rotating configurations R₀(M₀). It should be taken into account that the case we are interested in – the spin-up accretion of the NS in the stripping mechanism – is extremely complicated. In fact, the mass transfer rate is |$\dot{M}_1 \sim {\rm M}_{\odot }\,\mathrm{s}^{-1}$| (e.g. Clark & Eardley 1977) so the resulting object is likely to be highly hot and unsteady, its rotation should also be non-uniform, etc. So we need a simple and efficient method that will give a good first approximation, and will at least qualitatively take into account the basic trends of the phenomenon.

Let us assume that we have a basic parameter in the problem – the ratio of the spin frequency of the NS Ω to the maximum rotational (Keplerian) frequency,

$$\begin{eqnarray} \Omega _{\mathrm{K}}= \sqrt{\frac{GM}{R^3}}. \end{eqnarray}$$

(B1)

Here we define Ω_K as the spin frequency at which the equatorial mass shedding begins. The expression (B1) in the general case contains a numerical parameter, which for simplicity we further take as equal to unity. For a given EoS we have a mass–radius relation of the static NS R₀(M₀). The results of the numerical calculations (e.g. fig. 4 from Martinon et al. 2014) motivate us to look for the mass M versus equatorial radius R relation of the rotating NS in the form

$$\begin{eqnarray} M=M_0\left[1+\alpha _{\mathrm{M}}\frac{\Omega ^2}{\Omega _{\mathrm{K}}^2}\right]=M_0\left[1+\alpha _{\mathrm{M}}\frac{\Omega ^2 R^3}{GM}\right], \end{eqnarray}$$

(B2)

$$\begin{eqnarray} R=R_0(M_0)\left[1+\alpha _{\mathrm{R}}\frac{\Omega ^2}{\Omega _{\mathrm{K}}^2}\right]=R_0(M_0)\left[1+\alpha _{\mathrm{R}}\frac{\Omega ^2 R^3}{GM}\right]. \end{eqnarray}$$

(B3)

These equations are non-linear and make sense only for Ω ≤ Ω_K. The analysis of data on the Keplerian configurations (Lasota, Haensel & Abramowicz 1996; Haensel et al. 2009; Martinon et al. 2014; Li et al. 2016) gives the next values of the input parameters: α_M ≈ 0.2, α_R ≈ 0.4.

The mass–radius lines obtained by our approximations (B2) and (B3) are shown in Fig. B1. The NS EoSs are taken from Fantina et al. (2013) and Pearson et al. (2018). The green line illustrates the static configuration R₀(M₀), and the red one gives the Keplerian configuration with Ω = Ω_K. The maximum of the red line corresponds to the NS maximum mass M_max. The different coloured lines show the M–R relations with various frequencies ν = Ω/2|$\pi$| (in Hz), where ν = 716 and 1122 Hz correspond to the frequencies of some of the fastest spinning NSs: PSR J1748−2446ad (Hessels et al. 2006) and XTE J1739−285 (Kaaret et al. 2007). For comparison, the numerical calculations (different markers) obtained by the GR code rns (Stergioulas & Friedman 1995; Nozawa et al. 1998) are also given. It is seen that our approximation gives not only qualitative, but also a good quantitative agreement with the GR calculations.

Now we determine the moment of inertia of the rotating NS and replace the variable Ω in the expressions (B2) and (B3) with the new variable J – the spin angular momentum of the NS. As the dimensionless moment of inertia I* = I/MR², we can use, for example, the simplest approximation of Ravenhall & Pethick (1994):

$$\begin{eqnarray} I^{*}=\frac{0.21}{1-x_{\mathrm{GR}}}, \end{eqnarray}$$

(B4)

where x_GR = R_g/R = 2GM/Rc². The main feature of our application of this approximation is that we will use it for the rotating stars, substituting into equation (B4) the equatorial radius R, which implicitly depends on the mass M and the spin angular momentum J of the star. It is motivated by the fact that at a fixed mass M, the moment of inertia I grows with increasing Ω, as it should be, due to a sharp increase of R. However, the dimensionless moment of inertia I* = I/MR² decreases due to the denominator 1 − x_GR, which is also expected from the comparison with the case of the rotating polytropic configurations with index n = 1 (Cook, Shapiro & Teukolsky 1994).

For the particular EoS we can also use more complex approximation formulae for the moment of inertia (e.g. Lattimer & Prakash 2001; Bejger & Haensel 2002). For the BSk22, BSk24, BSk25, and BSk26 EoSs (Pearson et al. 2018), used in the main part of the paper, we derive our own approximation:

$$\begin{eqnarray} I^{*} = \frac{a_1 x^{n_1}(1+a_{2}x^{n_2})}{1+a_{3}x^{n_3}}, \end{eqnarray}$$

(B5)

where x = M/R (in M_⊙ km⁻¹). Approximation coefficients (B5) are shown in the table below.

a1	n1	a2	n2	a3	n3
136.9	1.512	11.04	2.713	214.4	1.261

The presented approximations also satisfy all the properties of the rotating polytropes with index n = 1 described above: the growth of the moment of inertia I and the decrease in the dimensionless moment of inertia I* = I/MR² with increasing Ω or J. Nevertheless, our calculations have shown that the particular type of function I(M, R) has weak influence on the MNS spin-up. Thus, the total stripping time – our main dynamical parameter – varies within hundredths of a second for different approximations of the moment of inertia.

Now let us finally formulate our approach. We set one of the formulae for the moment of inertia I = I(M, R) and determine the spin angular momentum J = IΩ. Then the equations (B2) and (B3) for the mass and radius of the rotating NS can be rewritten as

$$\begin{eqnarray} M=M_0\left[1+\alpha _{\mathrm{M}}\frac{J^2 R^3}{GM I^2}\right], \end{eqnarray}$$

(B6)

$$\begin{eqnarray} R=R_0(M_0)\left[1+\alpha _{\mathrm{R}}\frac{J^2 R^3}{GM I^2}\right]. \end{eqnarray}$$

(B7)

We know the mass M and the angular momentum J of the star. Solving the non-linear equations (B6) and (B7) with respect to M₀ and R, we find the moment of inertia I(M, R) = I(M, J). The static criterion for the maximum mass of an isentropic static configuration discussed by Bisnovatyi-Kogan & Blinnikov (1974) states that it is an extremum (maximum) of the M curve at given J that defines the maximum mass of the star. Exceeding this extremum, the NS collapses into the BH.

In conclusion, we note that the approximation formulae (B6) and (B7) that we obtained to calculate I(M, J) and R(M, J) are not universal and require verification on a larger class of the NS EoSs. The approach described above will be refined and developed in another paper.