EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 695 (March 2025) A&A, 695 (2025) A16 Full HTML
Open Access
Issue
A&A
Volume 695, March 2025
Article Number A16
Number of page(s) 12
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/202452259
Published online 27 February 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Millisecond pulsars (MSPs) are widely considered to be old, spun-up neutron stars (NSs) characterized by high rotational frequencies, 𝒪(102 Hz), and weak surface magnetic fields, 𝒪(108 D). The formation of MSPs can be traced back to a recycling phase in a low-mass X-ray binary (LMXB; see Tauris & van den Heuvel 2023). During the recycling process, mass and angular momentum are accreted onto the NS from the companion star when the latter fills its Roche lobe (e.g., Bhattacharya & van den Heuvel 1991; Tauris & Savonije 1999; Tauris & van den Heuvel 2023). During mass transfer, tidal interactions circularize the orbit on a timescale that is much shorter (∼ 104 yr) than the mass accretion phase. Therefore, MSPs are expected to populate binary systems with very small eccentricities (Phinney 1992; Verbunt & Phinney 1995).

Although the vast majority of MSPs are indeed situated in systems with highly circularized orbits, a considerable fraction appear to be solitary (∼27% of known MSPs have spin periods ≤30 ms; see the ATNF pulsar catalog; Manchester et al. 2005)1. These isolated MSPs (iMSPs) raise questions about the recycling process, particularly about possible mechanisms that could eject, consume, or disrupt the MSP companion stars (Verbunt et al. 1987; Sun et al. 2019; Nurmamat et al. 2019; Jiang et al. 2020; Antoniadis 2021).

In addition, recent pulsar surveys have revealed the existence of binary MSPs with eccentricities of up to e ≃ 1. With the exception of MSPs in globular clusters and PSR J1903+0327 (Champion et al. 2008), which is a Galactic field MSP with a main-sequence companion and an eccentricity of e ≃ 0.4, eccentric MSPs (eMSPs) appear to have orbital periods between 20 and 50 d and eccentricities of 𝒪(0.1), (Deneva et al. 2013; Barr et al. 2013; Knispel et al. 2015; Camilo et al. 2015; Antoniadis 2014; Antoniadis et al. 2016; Stovall et al. 2019). Although chaotic gravitational interactions, such as the ejection of the least massive member of a hierarchical triple system progenitor, have the potential to explain the formation of eMSPs in clusters (e.g., Freire et al. 2011; Portegies Zwart et al. 2011), they are unlikely to produce systems that have similar binary parameters (see Deneva et al. 2013; Barr et al. 2013; Knispel et al. 2015; Camilo et al. 2015; Octau et al. 2018), such as those found in the Galactic field. However, several attempts have been made to investigate alternative formation channels that can accommodate the observed properties of eMSPs (e.g., see Freire & Tauris 2014; Antoniadis 2014; Jiang et al. 2015, 2021; Ginzburg & Chiang 2022, for an overview of various eMSP formation scenarios).

The formation of both iMSPs and eMSPs could be related to abrupt changes in the interior of the NS (Freire & Tauris 2014; Jiang et al. 2015; Alvarez-Castillo et al. 2019). For instance, the dependence of the NS structure on the equation of state (EoS) highlights the possibility of exotic phases and phase transitions in their interiors. Jiang et al. (2015) find that eccentricities of ∼0.11 − 0.15 will be induced by instantaneous gravitational mass loss of the NS due to a phase transition in its interior, even without invoking an explicit kick mechanism. However, for the abovementioned scenarios, it can be troublesome to explain the “instantaneity” of the transition when the NS spin is included. As pointed out by Glendenning et al. (1997) in the context of their scenario for an accretion-induced phase transition to a hybrid star branch that is connected with the NS one, the reconfiguration of the mass distribution in the NS interior is accompanied by a change in the moment of inertia. This would entail a “pirouette effect” (a spin-up), which in turn would inhibit the completion of the transition due to the sensitive dependence of the density profile inside the star on its spin state caused by centrifugal forces. Glendenning et al. (1997) obtain in their scenario a timescale of 105 years for the completion of the phase transition, which is in stark contrast with the requirement that the mass defect occur on a timescale much shorter than the orbital period.

A phase transition scenario that circumvents this problem assumes the existence of a third family of compact stars (henceforth “twin stars”; see Gerlach 1968) beyond the conventional white dwarf (WD) and NS classifications. Such a third family could emerge as a result of a sufficiently strong phase transition in the NS core, leading to the formation of hybrid stars with an outer core of hadronic matter and an inner core of de-confined quark matter (e.g., Schertler et al. 2000; Benic et al. 2015; Alvarez-Castillo et al. 2019; Alvarez-Castillo 2021). Given that baryon number and angular momenta can be conserved simultaneously during this process, the transition can occur almost instantaneously, on a dynamic timescale, which is insignificant compared to the orbital period.

Within the context of MSP formation scenarios, such an instantaneous transition to the third-family branch could be triggered by mass accretion during, or after, the LMXB phase (see, e.g., Zdunik et al. 2006; Bejger et al. 2017; Alvarez-Castillo et al. 2019). This transition could in turn induce instantaneous mass loss and/or a kick that can cause significant changes to the orbital dynamics (Jiang et al. 2015). More specifically, the sudden decrease in gravitational mass due to quark de-confinement in the core (owing to the “catastrophic rearrangement” and strong compactification of matter when transitioning to the third-family branch, which increases the gravitational binding energy; see Mishustin et al. 2003) may significantly increase eccentricity or even disrupt the binary (see also Jiang et al. 2021). The main aim of this work is to explore this possibility in more detail.

The investigation of such scenarios is important for advancing our understanding of dense matter under extreme conditions. Recent advances in observational techniques, such as gravitational wave (GW) detections from NS mergers (e.g., Abbott et al. 2017) and high-precision pulsar timing (e.g., Miller 2019; Riley 2019; Miller 2021; Bogdanov 2021), offer unprecedented opportunities to test and refine theoretical models that predict phase transitions (e.g., Bauswein et al. 2022). For instance, GWs from unstable fundamental f-mode oscillations can reveal the internal stellar structure and matter distribution characteristic of a first-order phase transition, as shown in Pradhan et al. (2024), where estimations for the same class of model used in this work are presented.

Khosravi Largani et al. (2024) simulated NS evolution with an accretion-induced phase transition in a neutrino radiation hydrodynamics treatment, implementing a three-flavor Boltzmann neutrino transport and a microscopic quark-hadron hybrid model EoS. This approach results in the production of gravitational radiation from f-mode excitation and the associated neutrino emission. The phase transition triggers an expansion of the accreted envelope with supersonic velocities resulting in a mass loss of ∼ 10−3 M⊙. It is worth noting that finite temperatures are considered in this approach, which correspond to an entropy per baryon of 0.1 − 0.5. Since finite temperatures soften the quark matter EoS, the onset densities for de-confinement are lowered and the latent heat is increased so that twin stars of identical mass but different radii and compactness (Gerlach 1968; Schertler et al. 2000) become possible even when at the vanishing temperature no gravitational instability is associated with the de-confinement transition (Hempel et al. 2016; Carlomagno et al. 2024b). The associated mass defect of such an accretion-induced thermal star quake can be as much as 3% of the solar mass (Carlomagno et al. 2024a). In fact, with the precision of the next generation of GW detectors, we expect to be able to study these scenarios. In general, the timely synergy between theoretical developments and observational capabilities places us at the forefront of uncovering the complex NS physics, where phase transitions emerge as a major component (Bauswein et al. 2022).

In this work we combined detailed binary evolution calculations with a simple analytic model for the response of the NS spin and radius to accretion to investigate whether twin stars can form in LMXBs. We also investigated whether such accretion-induced phase transitions could help explain the observed eMSP and iMSP populations. The text is organized as follows: In Sect. 2 we describe the methodology and physical assumptions. We present the simulation results in Sect. 3 and conclude with a summary and discussion on the implications of our calculations in Sect. 4.

2. Methodology

2.1. Initial setup and numerical calculations

Our main goal was to investigate whether first-order phase transitions can be triggered in LMXBs and recycled MSPs under realistic mass accretion and spin evolution scenarios. To this end, our working setup was divided into three distinct components:

  1. An EoS for dense nuclear matter that allows for a first-order phase transition paired with a mass-radius relation for spinning NSs;

  2. A stellar evolution model for the binary system, providing the mass-accretion profile as a function of time for the NS;

  3. An accretion model that allows us to calculate the spin, mass, and radius evolution of the NS during and after the LMXB phase.

They are described in detail in the remainder of the section.

2.1.1. Equation of state and mass-radius relation

For the EoS, we selected a multi-polytrope model suitable for describing both nuclear and quark matter within NSs. This EoS was chosen for its ability to (a) account for a strong first-order phase transition from hadronic to quark matter above 1.4 M⊙ (depending on the angular momentum) and (b) satisfy current constraints derived from multi-messenger astronomy (Antoniadis et al. 2013; Abbott 2018; Miller 2019; Riley 2019; Cromartie et al. 2020).

The outer layers, characterized by densities n ≤ 0.5n0, where n0 denotes the nuclear saturation density, are modeled using the BPS EoS (Baym et al. 1971). An intermediate density regime with 0.5n0 < n ≤ 1.1n0 is assumed to consist of homogeneous matter in β-equilibrium. For regions with densities n > 1.1n0, we distinguished four regimes, each described by a polytropic EoS of the form

P ( n ) = κ i ( n n 0 ) Γ i , $$ \begin{aligned} P(n) = \kappa _i \left( \frac{n}{n_0} \right)^{\Gamma _i}, \end{aligned} $$(1)

henceforth referred to as the ACB5 EoS as in Paschalidis et al. (2018).

The first polytrope, (i = 1), is a fit to the stiffest EoS version provided in Hebeler et al. (2013). The second polytrope, (i = 2), describes the first-order phase transition. The polytropic index Γ allows us to distinguish between a sharp transition (Maxwell construction, Γ = 0) and a non-sharp transition (Gibbs construction, Γ ≠ 0). Since we required a strong transition to produce a large jump in energy density, the second polytrope was defined as a region of constant pressure Pc = κ2, where Γ2 = 0 is imposed by the Maxwell construction. The remaining two polytropes (i = 3, 4) describe regions with densities that exceed the critical value for a first-order transition and correspond to stiff quark matter that should be able to support a NS with masses up to ∼2.0 M⊙ (e.g., Cromartie et al. 2020). The different density regimes are thermodynamically joined in a consistent manner.

To infer the corresponding mass-radius relation for rigidly spinning NSs, we used the numerical integrator described in Alvarez-Castillo et al. (2019). In this code, the integration of the structure equations proceeds radially outward from the center of the NS, where the density is highest, to its surface, where the pressure falls to zero, thus determining the total mass and radius for various central densities and angular momenta, J.

The resulting mass-radius relations are illustrated in Fig. 1 where the baryonic mass (left) and gravitational mass (right) mass are plotted against the equatorial radius of the NS. The lower right section of each panel depicts the mass-radius relation for NSs within the hadronic branch, corresponding to the first ACB5 polytrope and indicating a composition of purely hadronic matter. Moving toward the upper left segment of each panel, we shift into the twin (hybrid) star branch, representing NSs that have undergone a phase transition to include a quark matter core, described by the ACB5’s third and fourth polytropes. This transition is shown by a dashed black curve, representing the second polytrope of the ACB5 EoS and indicating the highest mass that a stable NS composed entirely of hadronic matter can sustain. Beyond this threshold, additional mass triggers a phase transition to quark matter, moving the star to the hybrid branch.

thumbnail Fig. 1.

Left Panel: Baryonic mass vs. equatorial radius. Dotted and dot-dashed curves indicate lines of constant angular momentum and constant frequency, respectively. The dashed black line highlights the points of maximum stability in the hadronic branches, signaling the onset of a phase transition (PT). Thin, horizontal magenta lines show trajectories where angular momentum (J) and baryonic mass (Mb) are conserved. The thick magenta curve connects points on the hybrid star branches that can be reached through a collapse conserving both J and Mb. Right Panel: Similar to the left panel, but with the y-axis representing gravitational mass. Star markers denote the endpoint trajectories for direct (M2) and delayed (M3) collapse models. The arrows are inclined because gravitational mass, unlike baryonic mass, is not conserved during the PT.

The role of spin and angular momentum in dictating the precise moment of transition is crucial; higher angular momentum allows a NS to support more mass in the hadronic phase. This effect is illustrated by the dotted curves, which denote varying levels of angular momenta, normalized by J0 = G M⊙2/c. Here, the innermost (dark blue) dotted line represents a static NS, and the outermost (light yellow) line corresponds to J/J0 = 0.45. Beyond this point, the hadronic branch extends past the shedding limit, where equatorial matter attains Keplerian escape velocity and is ejected.

For the calculations that follow, we assumed that a NS reaching the transition line undergoes a phase transition on the dynamic timescale, under conservation of the angular momentum and baryonic mass. The transition moves the NS from the maximal stability threshold on the hadronic branch (dashed line), to a corresponding position (in terms of angular momentum and baryonic mass) on the third family. This principle is visually represented by the magenta curve, linking points on the hybrid branch that could be attained through the transition and defining the birth radius of the twin star.

2.1.2. Mass transfer and binary evolution models

To obtain realistic mass-transfer models, we performed detailed numerical binary calculations using the implicit one-dimensional code Modules for Experiments in Stellar Astrophysics (MESA v12778; Paxton et al. 2011, 2013, 2015, 2018). Our binary systems consisted of a zero-age main-sequence (ZAMS) donor with an initial metallicity Z = 0.02 and an NS accretor, which was treated as a point mass. We calculated three models chosen as representative of three distinct cases: an LMXB in which a phase transition never occurs (M1), one in which the transition occurs during the mass-transfer phase (M2), and a model for which the transition occurs after the accretion phase, due to the spin-down of the NS (M3). The initial binary parameters for each model are summarized in Table 1.

Table 1.

Initial binary parameters for our stellar evolution models.

For the donor-star evolution, we used the type-2 OPAL Rosseland mean opacity tables (Iglesias & Rogers 1996). Convection was modeled using the modified mixing length theory prescription of Henyey et al. (1965) with a mixing length parameter of αML = 1.8. In our models, we enabled convective premixing while trying to avoid increases in the abundance of species that were burned, where possible. Stability against convection was determined according to the Ledoux criterion (Ledoux 1947). Furthermore, we used semi-convection, which we treated as a diffusive process, with an efficiency parameter of αSEM = 1.0 (Langer 1991). Lastly, we accounted for the overshooting of convective material beyond the convective layers by adopting an exponential overshooting efficiency of fov, core = 0.016 and fov, env = 0.0174 in the core and envelope, respectively.

Rotational mixing was taken into consideration, including the effects of Eddington-Sweet circulations, secular and dynamical instability, and the Goldreich-Schubert-Fricke instability (see Heger et al. 2000, for details). The contribution of these instabilities to the total diffusion coefficient was reduced by a mixing efficiency factor, fc = 1/30 (Chaboyer & Zahn 1992; Heger et al. 2000). Moreover, we employed the Spruit-Tayler dynamo to compute the internal magnetic field strength and the corresponding transport of angular momentum, as described in Spruit (2002), Heger et al. (2005).

To calculate the mass loss rate due to stellar wind, we implemented the cool red giant branch wind scheme from Reimers (1975):

M ˙ 1 , wind = − 4 × 10 − 13 η ( R 1 R ⊙ ) ( L 1 L ⊙ ) ( M ⊙ M 1 ) [ M ⊙ yr − 1 ] , $$ \begin{aligned} \dot{M}_{1, \text{ wind}} = -4 \times 10^{-13} \eta \left(\frac{R_1}{R_\odot }\right) \left(\frac{L_1}{L_\odot }\right)\left(\frac{M_\odot }{M_1}\right)\quad [M_\odot \,\text{ yr}^{-1}], \end{aligned} $$(2)

using a scaling factor of η = 0.1. In cases of hydrogen-poor donor stars whose hydrogen mantles have been stripped (surface hydrogen mass fraction XS < 0.4), we applied the prescription of Nugis & Lamers (2000) with a scaling factor of η = 1.0. This factor is used in MESA under the Dutch hot wind scheme (Glebbeek et al. 2009).

We assumed that our binary models were initially tidally synchronized to the orbit at the beginning of the evolution (ZAMS stage) and that the initial eccentricity was negligible. These assumptions are justified as tidal forces would circularize (and synchronize) the orbit on a much shorter timescale compared to the main-sequence lifetime of a low-mass donor star (e.g., Verbunt & Phinney 1995).

To compute the mass transfer rates during Roche lobe overflow (RLOF), we used the prescription suggested by Kolb & Ritter (1990). In addition, an isotropic reemission scenario of mass transfer was adopted (see Tauris & van den Heuvel 2006, for a review). In this scenario, we assumed that mass flows conservatively from the donor to the NS accretor via RLOF, and a fraction of the transferred material, β, is lost (reemitted) from the vicinity of the NS as an isotropic fast wind, carrying away the specific angular momentum of the NS. Hence, the NS’s efficiency of accretion, ϵ, is defined as

ϵ = 1 − α − β − δ , $$ \begin{aligned} \epsilon = 1 - \alpha - \beta - \delta , \end{aligned} $$(3)

where α is the fraction of mass lost directly from the donor via winds, and δ is the fraction of mass lost from a circumbinary coplanar toroid. Here, we omitted any losses from winds or circumbinary toroids (α = δ = 0) and assumed that 50% of the transferred mass is reemitted (β = 0.5) resulting in an accretion efficiency of ϵ = 0.5. Furthermore, we assumed that the accretion rate of the NS is limited by the Eddington mass accretion rate, which for accreted material composed of pure ionized hydrogen takes the form

M ˙ Edd = 1.5 × 10 − 8 ( M 2 , i 1.3 M ⊙ ) [ M ⊙ yr − 1 ] , $$ \begin{aligned} \dot{M}_{\text{ Edd}} = 1.5 \times 10^{-8} \left(\frac{M_{2,i}}{1.3 M_\odot }\right)\quad [M_\odot \,\text{ yr}^{-1}], \end{aligned} $$(4)

where M2, i is the initial mass of the accretor (Misra et al. 2020). Here, the mass accreted would not change the Eddington limit significantly, and thus we fixed the Eddington limit to ṀEdd = 1.5 × 10−8 M⊙ yr−1. Excess material, exceeding the Eddington accretion rate, was assumed to carry the specific angular momentum of the NS.

For orbital angular momentum losses, we considered mechanisms such as GW radiation, mass lost from the system, spin-orbit coupling due to tidal effects, and magnetic braking for donors with convective envelopes. The orbital angular momentum loss due to GW radiation were calculated using the formula

d J gw dt = − 32 5 G 7 / 2 c 5 M 1 2 M 2 2 ( M 1 + M 2 ) 1 / 2 α 7 / 2 , $$ \begin{aligned} \frac{dJ_{\text{ gw}}}{dt} = - \frac{32}{5}\frac{G^{7/2}}{c^5}\frac{M_1^2 M_2^2 (M_1 + M_2)^{1/2}}{\alpha ^{7/2}}, \end{aligned} $$(5)

where G is the gravitational constant, c is the speed of light in vacuum, α is the binary separation, and M1, M2 are the masses of the donor and the accretor, respectively (Tauris & van den Heuvel 2006). Angular momentum loss due to magnetic braking was computed following the prescription from Rappaport et al. (1983):

d J mb dt = − 3.8 × 10 − 30 M 1 R ⊙ 4 ( R 1 R ⊙ ) γ Ω 1 3 [dyn cm] , $$ \begin{aligned} \frac{dJ_{\text{ mb}}}{dt} = -3.8 \times 10^{-30} M_1 R_\odot ^4 \left(\frac{R_1}{R_\odot }\right)^{\gamma }\Omega _1^3\quad \text{[dyn} \text{ cm]}, \end{aligned} $$(6)

where M1, R1, and Ω1 are respectively the mass, radius, and angular velocity of the donor star. Here, the magnetic braking index was fixed to γ = 4.

We let our models evolve until the donor star forms a WD or the number of computational steps needed for convergence exceeds a limit of 3 × 105. In all cases, the spatial and temporal resolution were set to their default values. Nonetheless, to assess whether rapid torque fluctuations (see Sect. 2.2) during the equilibrium phase are physical or numerical in nature, and whether they impact our conclusions, we performed a time-resolution sensitivity study, by running extra MESA tracks. To this end, we varied the time resolution during the mass-transfer phase of our models across a range of values, from 8 × 10−5 to 8 × 10−3, to test the sensitivity of the spin evolution and check the robustness of our conclusions.

2.1.3. Evolution of NS spin and structure

To simulate the response of the interior of the NS to mass accretion, we adopted the accretion model of Tauris (2012), which allowed us to calculate the accretion torque as a function of time, N(t). We summarize this model here for the sake of completeness.

The torque model is based on two primary assumptions: (a) spherical accretion, which simplifies the inflow dynamics to a scenario where the flow of ionized gas is predominantly governed by gravity until it reaches the magnetosphere, and (b) the magnetic dipole assumption, which posits that the NS surface magnetic field can be approximated as a dipole. As matter falls on the surface of the NS, the ram pressure of the infalling material is counteracted by the magnetic pressure (e.g., Corbet 1996). Assuming a spherical accretion scenario, these two pressures can be equated to determine the so-called magnetospheric radius,

R mag ≈ 7.8 ( B 10 8 G ) 4 / 7 ( R ns 10 km ) 12 / 7 × ( M ns 1.4 M ⊙ ) − 1 / 7 ( M ˙ M ˙ Edd ) − 2 / 7 [ km ] , $$ \begin{aligned} R_\mathrm{mag}&\approx 7.8 \left(\frac{B}{10^8\,\mathrm{G} }\right)^{4/7} \left(\frac{R_\mathrm{ns} }{10\,\mathrm{km} }\right)^{12/7} \nonumber \\&\times \left(\frac{M_\mathrm{ns} }{1.4\,\mathrm{M} _{\odot }}\right)^{-1/7} \left(\frac{\dot{M}}{\dot{M}_\mathrm{Edd} }\right)^{-2/7}\,[\mathrm{km} ], \end{aligned} $$(7)

inside of which the plasma flow is dictated by the magnetic field corotating with the NS (Andersson et al. 2005). This radius provides an estimate of the location of the inner edge of the accretion disk. The outer boundary of the magnetosphere is indicated by the light-cylinder radius:

R lc = c P spin 2 π , $$ \begin{aligned} R_\mathrm{lc} = \frac{cP_\mathrm{spin} }{2\pi }, \end{aligned} $$(8)

which signifies the location where the corotation velocity equals the speed of light (e.g., Rappaport et al. 2004; Tauris 2012; Bhattacharyya & Chakrabarty 2017, and references therein).

Another important length scale is the corotation radius,

R co = ( G M ns ) 1 / 3 ( P spin 2 π ) 2 / 3 , $$ \begin{aligned} R_\mathrm{co} = (G M_\mathrm{ns} )^{1/3}\left(\frac{P_\mathrm{spin} }{2\pi }\right)^{2/3}, \end{aligned} $$(9)

which is defined as the radial distance at which material corotating with the NS attains the Keplerian angular frequency (e.g., Corbet 1996; Bhattacharyya & Chakrabarty 2017; Bhattacharyya 2023). If Rmag < Rco, the accreted material follows the magnetic field lines and is deposited on the surface, carrying angular momentum and exerting a positive torque that tends to spin up the NS. However, if Rmag > Rco a centrifugal barrier arises, inhibiting the accretion flow as the material accelerates to super-Keplerian velocities (e.g., Bhattacharyya & Chakrabarty 2017; Bhattacharyya 2023, and references therein). Within this so-called propeller regime, the accelerated plasma may reach the escape speed and be ejected from the system. The expulsion of matter exerts a negative (braking) torque extracting angular momentum from the NS and causing it to spin down. Therefore, it is expected that accretion cannot spin up the NS beyond the point where Rmag = Rco. The fastest spin a NS can acquire from accretion is the so-called equilibrium spin period:

P eq ≈ 0.26 ( B 10 8 G ) 6 / 7 ( R 10 km ) 18 / 7 × ( M 1.4 M ⊙ ) − 5 / 7 ( M ˙ M ˙ Edd ) − 3 / 7 [ ms ] $$ \begin{aligned} P_\mathrm{eq}&\approx 0.26 \left(\frac{B}{10^8\,\mathrm{G} }\right)^{6/7} \left(\frac{R}{10\,\mathrm{km} }\right)^{18/7} \nonumber \\&\times \left(\frac{M}{1.4\,\mathrm{M} _{\odot }}\right)^{-5/7} \left(\frac{\dot{M}}{\dot{M}_\mathrm{Edd} }\right)^{-3/7}\,\mathrm{[ms]} \end{aligned} $$(10)

(Andersson et al. 2005).

To calculate the resulting accretion torque as a function of time, we used the formulation suggested by Tauris (2012):

N ( t ) = n ( ω ) [ M ˙ ( t ) G M ns R mag ( t ) ξ + μ 2 9 R mag 3 ( t ) ] − E ˙ dipole ( t ) Ω ( t ) , $$ \begin{aligned} N(t) = n(\omega ) \left[\dot{M}(t) \sqrt{GM_\mathrm{ns} R_\mathrm{mag} (t)}\xi + \frac{\mu ^2}{9R^3_\mathrm{mag} (t)}\right] - \frac{\dot{E}_\mathrm{dipole} (t)}{\Omega (t)}, \end{aligned} $$(11)

Equation (11) describes the net torque governing the spin evolution of the NS. The first term represents the accretion torque. The second term corresponds to the magnetic torque arising from the interaction between the NS’s dipole field and the magnetosphere. The third term represents the spin-down torque caused by magnetic dipole radiation. Here, n(ω) is a dimensionless quantity defined as n ( ω ) = tanh ( 1 − ω δ ω ) $ n(\omega) = \tanh \left( \frac{1 - \omega}{\delta \omega} \right) $, where ω is the fastness parameter given by ω = ( R mag R co ) 3 $ \omega = \sqrt{\left( \frac{R_{\mathrm{mag}}}{R_{\mathrm{co}}} \right)^3} $. The parameter ξ is a factor of order unity and is fixed to ξ = 1 in our calculations. The parameter μ denotes the magnetic dipole moment, Ėdipole(t) represents the time-dependent dipole energy loss, and Ω is the spin frequency. The mass-transfer rate, Ṁ, was taken directly from our MESA binary stellar evolution models. The terms G, Mns, and Rmag represent the gravitational constant, NS mass, and magnetic radius, respectively. The width of the transition zone near the magnetospheric boundary was assumed to be small enough (δω = 0.02) so that the quantity n(ω) corresponds to a step function, n(ω) = ± 1. Finally, for a recycled pulsar, it is anticipated that the magnetic field decreases significantly during the early accretion phase, dropping by several orders of magnitude from its original strength (Shibazaki et al. 1989). Because this reduction in the magnetic field is believed to occur over a timescale shorter than the duration of the mass-transfer phase (Bhattacharya & van den Heuvel 1991; Bhattacharyya & Chakrabarty 2017), we considered the surface magnetic field strength to be constant with a value of B = 1.0 × 108 G.

In practice, we used the temporal evolution of mass loss from the donor star Ṁd(t) provided by MESA, together with the mass-radius relations described above to calculate the evolution of the NS radius, spin and mass as a function of time. More specifically, the angular momentum Jns, was calculated as the discrete summation:

J ns ( t n ) = J ns ( t 0 ) + ∑ i = 1 n Δ t i − 1 N i − 1 , $$ \begin{aligned} J_\mathrm{ns} (t_n) = J_\mathrm{ns} (t_0) + \sum _{i=1}^n \Delta t_{i-1} N_{i-1}, \end{aligned} $$(12)

where the timestep Δt was taken from our MESA simulations, and the torque was evaluated using Eq. (11). Since the NS is considered old, we initialized the angular momentum as Jns(0) = 0. At each step, the radius, Rns(Mns, Jns), and spin frequency, fns(Mns, Jns), were interpolated from our tabulated EoS data, using a linear radial basis function.

Although the temporal evolution of the NS mass was already provided by MESA assuming that the NS accretes half of the mass lost by the donor (see Sect. 2.1.2), we reevaluated this parameter a posteriori using the torque model described in this section. More specifically, we assumed that mass accretion is limited by the Eddington mass accretion rate, and that the NS accumulates mass only when the fastness parameter is negative. Due to this treatment, the NS mass evolves slightly differently than the accreting point mass as simulated with MESA, with M3 showing the largest difference, showing an increase in the average accretion efficiency of approximately 10%.

2.2. Accretion torque fluctuations

The torque equation, as expressed in Eq. (11), leads to fluctuations around the equilibrium frequency during the spin equilibrium phase, a phenomenon extensively reported in existing literature (e.g., Bildsten et al. 1997; Chandra 2025). For example, Tauris (2012) highlights how minor variations in the mass-transfer rate can impact the fastness parameter, resulting in rapid oscillations in the direction of the accretion torque. Similarly, Elsner et al. (1980) discuss the occurrence of alternating periods of spin-up and spin-down in pulsating X-ray sources, attributing these phases to variations in the luminosities of disk-fed NSs, which in turn affect the accretion torque. To assess the impact of these torque fluctuations on our findings, we conducted further analysis using additional MESA models to explore the sensitivity of spin evolution to the timing resolution. We find that larger time intervals fail to accurately capture the Roche lobe decoupling phase, which interrupts the spin equilibrium. However, the evolution trajectories for spin, mass, and radius remained consistent across both larger and smaller time intervals. Consequently, even when adjusting the time resolution by an order of magnitude, the pre-transition mass is altered by less than 3%. Therefore, we conclude that our findings are sufficiently robust against naturally occurring torque fluctuations. However, studies of spin variations in numerous X-ray pulsar and accreting MSPs have raised a number of issues regarding the standard accretion model. Observational evidence of unexpected spin reversals and spin-down phases during outbursts in X-ray pulsars (e.g., Bildsten et al. 1997; Galloway et al. 2002) or accretion during the propeller phase challenge the conventional model used in this study and call for more detailed investigations on the accretion torque exerted on magnetized NSs.

2.3. Orbital reconfiguration

The evolution of the orbit and the NS spin were followed until the NS reached the phase transition threshold (see Fig. 1). It was then assumed that the phase transition occurs instantaneously (i.e., on a timescale that is much shorter than the orbital period). To investigate the impact of the phase transition on the orbit, we used the prescriptions of Hills (1983) and Tauris et al. (2017) in which the post-transition semimajor axis is given as

a f a i = 1 − Δ M / M 1 − 2 Δ M / M − ( w / v rel ) 2 − 2 cos θ ( w / v rel ) , $$ \begin{aligned} \frac{a_\mathrm{f} }{a_\mathrm{i} } = \frac{1 - \Delta M/M}{1 - 2\Delta M/M - (w/v_\mathrm{rel} )^2 - 2\cos \theta (w/v_\mathrm{rel} )}, \end{aligned} $$(13)

where ΔM is the instantaneous mass defect corresponding to the released gravitational binding energy during the phase transition, M is the total mass of the pre-transition system, vrel is the relative velocity between the two stars ( v rel = G M / a i $ v_{\mathrm{rel}} = \sqrt{GM/a_i} $), w is the magnitude of the kick velocity, and θ is the kick angle between the kick velocity vector, w, and the pre-transition orbital velocity vector. The eccentricity of the post-transition binary system is given as

e = 1 + 2 E orb , f L orb , f 2 μ f G 2 M f , 1 2 M f , 2 2 , $$ \begin{aligned} e = \sqrt{1 + \frac{2E_\mathrm{orb,f} L_\mathrm{orb,f} ^2}{\mu _\mathrm{f} G^2 M_\mathrm{f,1} ^2 M_\mathrm{f,2} ^2}}, \end{aligned} $$(14)

where L orb , f = a i μ f ( v rel + w cos θ ) 2 + ( w sin θ sin ϕ ) 2 $ L_{\mathrm{orb,f}} = a_{\mathrm{i}} \mu_{\mathrm{f}} \sqrt{(v_{\mathrm{rel}} + w\cos\theta)^2 + (w\sin\theta \sin\phi)^2} $ is the post-transition orbital angular momentum, with ϕ being the kick angle on the plane perpendicular to the pre-transition velocity vector, μf is the post-transition reduced mass, and Eorb, f = −GMf, 1Mf, 2/2af is the post-transition orbital energy.

To investigate the post-transition orbital configurations for our models, we considered a mass defect of 0.016 M⊙ and secondary kicks with magnitudes (w) up to 100 km s−1 (see Sect. 2.4 for the justification). Furthermore, we assumed that the kick lacks a preferential orientation, leading us to model the kick angles θ and ϕ as uniformly distributed variables within a 4π solid angle.

2.4. Mass defect and secondary kicks

To provide a reliable estimate for the mass defect parameter in our simulations, we discuss here the estimates from the literature for the case of a “catastrophic rearrangement” due to a mass twin transition. According to Mishustin et al. (2003), the energy defect in the gravitational field of about ΔE ∼ 7 × 1051 erg can heat the star to an initial temperature of T ≃ 40 MeV. At this temperature, neutrinos are copiously produced and trapped inside the star because their mean free path is much smaller than the size of the star. This leads to the establishment of a neutrino chemical potential μν in the direct and inverse β-decay processes, which establish the β-equilibrium including neutrinos. In a microscopic model for the quark de-confinement transition, the neutrino chemical potential influences on the critical density of the transition. Estimating the mass defect, one compares the gravitational masses of NSs in the initial and final states (before and after the transition) at a fixed baryon number. Without considering neutrino trapping, the mass defect for the EoS adopted here was calculated by Alvarez-Castillo (2021) to be in the range of 0.5–1.0% for NS masses between 1.2 and 1.6 M⊙. The mass defect due to the neutrino trapping–untrapping transition, modeled by a drop of the neutrino chemical potential from 150 MeV to zero, has been found by Aguilera et al. (2004) to be ΔM ≈ 0.01 M⊙. Due to a nonlinear dependence of ΔM on μν one can estimate the mass defect to be doubled for μν ≈ 200 MeV.

In Aguilera et al. (2004) a “normal” de-confinement transition without catastrophic rearrangement to a mass twin configuration has been considered. A scenario of neutrino untrapping in the cooling of a color superconducting quark star from a two-flavor, superconducting (2SC) phase to a three-flavor (Color-Flavor-Locked; CFL phase) mass-twin configuration has been examined in great detail for μν = 200 MeV by Sandin & Blaschke (2007). The untrapping transition led to a mass defect of 0.05 M⊙ and the subsequent transition to the twin star with a CFL quark matter core produced a further mass defect of 0.03 M⊙. It has recently been shown that the mass defect associated with an accretion-induced thermal star quake can reach 3% of a solar mass (Carlomagno et al. 2024a).

Figure 2 shows the baryonic mass (dot-dashed line) and gravitational mass (solid line) as functions of the central energy density for a fixed angular momentum, J/J0 = 0.3. From this plot, the mass defect expected from a phase transition from hadronic to quark matter can be constructed by analyzing changes in these mass functions and identifying the critical points of the transition. By locating the central energy density at the onset of the transition and noting that the baryonic mass remains constant, one can determine the energy jump and identify the central energy density marking the completion of the transition. Tracing the corresponding gravitational mass at these points allowed us to compare the gravitational mass before and after the transition. The difference between these values represents the mass defect, which quantifies the mass converted into energy during the phase transition.

thumbnail Fig. 2.

Baryonic mass (Mb) and gravitational mass (Mg) as functions of the central energy density (εc) for a fixed angular momentum (J). The solid black line represents the gravitational mass, while the dashed black line represents the baryonic mass. Vertical gray lines mark the central energy densities at the onset (εc = 275 MeV/fm3) and completion (εc = 595 MeV/fm3) of the phase transition (PT). Black dotted lines and the blue-shaded region denote unstable configurations. The inset zooms into the PT region, illustrating the mass defect (ΔMg) of 0.016 M⊙ associated with the transition from hadronic to quark matter. The horizontal gray lines indicate constant values of Mb and Mg during the PT. The star marker denotes the birth mass of the twin star configuration at Mg = 1.499 M⊙.

Our present EoS in the form of a multi-polytrope model is too simple to enable a more precise quantitative estimate of the mass defect for a realistic dense-matter model with neutrino untrapping and (three-flavor) color superconductivity. This task is left to separate work. In any case, the estimates provided in this subsection justify a variation of the mass defect parameter in a range 0.1% up to 5% assumed in Jiang et al. (2015, see also Jiang et al. 2021).

Besides orbital rearrangement due to a mass defect, we also considered the possibility of a secondary kick during the transition (Hobbs et al. 2005; Bombaci & Popov 2004). For the range of mass defects considered here, anisotropic emission could give rise to secondary kicks with magnitudes of up to w ≃ 20 000 km s−1. One possible mechanism that can give rise to such anisotropies is the channeling of the neutrino emission along a preferential axis formed by magnetic vortex lines together with parity violation in the strong internal magnetic field of the NS (see Berdermann et al. 2006; Kaminski et al. 2016; Fukushima & Yu 2024, for recent examples with typical B-field strengths of B ∼ 1012 Gauss).

The electromagnetic rocket effect, which was introduced in Harrison & Tademaru (1975) and revisited by Agalianou & Gourgouliatos (2023), is based on a displacement of a similarly strong dipolar magnetic field from the center of the NS and its sufficiently fast rotation. As it was discussed in Agalianou & Gourgouliatos (2023), observational evidence provided by the Neutron Star Interior Composition Explorer (NICER) X-ray observatory, particularly from investigating the location of hot spots on the surfaces of MSPs (Miller 2019; Riley 2019), indicates that the magnetic field diverges from a centered dipole configuration. Based on the results for PSR J0030+0451, Kalapotharakos et al. (2021) estimated the magnetic field structure and found an off-centered magnetic field consisting of dipole and quadruple components.

Since both the neutrino and electromagnetic rocket mechanisms require strong magnetic fields of at least the young pulsar field strength B ∼ 1012 Gauss, such effects seem at first glance unlikely for “old” MSPs with low magnetic fields of B ∼ 108 Gauss. However, in systems with mass accretion from a companion star like in the case of LMXBs, the argument of Chevalier (1989) can be applied, that the strong magnetic field of the interior gets buried in the crust of the pulsar (e.g., Bhattacharya & van den Heuvel 1991; Cumming et al. 2001; Konar & Choudhuri 2004, and references therein). Thus, a small surface magnetic field is compatible with a super-strong magnetic field in the interior. Due to such a magnetic field profile, kick velocities of a few dozen km s−1 can be produced in MSPs, as it has been demonstrated by Agalianou & Gourgouliatos (2023) by applying the electromagnetic rocket effect to the case of PSR J0030+0451.

A gravitational recoil effect arising from the anisotropic redistribution of matter following the phase transition, could also produce a secondary kick. This effect is independent of the magnetic field strength, making it a viable secondary kick mechanism for pulsars across a wide range of magnetic field intensities, including old MSPs with weak magnetic fields. Considering the above, we included a possible secondary kick mechanism in our study by parametrically varying the kick velocity parameter in the range w = 0, …, 100 km s−1.

3. Results

3.1. Mass and spin accretion

Figure 3 illustrates the evolution of our binary models from approximately the onset of the LMXB phase. The top and middle panels showcase that, due to their comparable initial setups, M1 and M2 have similar characteristics. Initially, the donor star in each model spends most of its life in the main sequence. After approximately 12.2 Gyr, it reaches the tip of the red giant branch and it fills its Roche lobe, initiating a mass flow toward the NS with a mass-transfer rate of Ṁ ≃ 10−9 Ṁ yr−1. As detailed in Sect. 2, the interaction of the accreted matter with the magnetic field of the NS generates an accretion torque (Tauris 2012). We find that the exerted torque significantly decreases the NS spin to nearly 1.1 ms shortly after the RLOF is initiated. Over the next circa 56 Myr the mass-transfer rate gradually decreases, allowing the spin of the NS to approach a stable equilibrium period, exhibiting only minor fluctuations due to rapid torque reversals. At this point, the binary orbit widens to a period of Porb ≃ 16.5 d with the donor star separating from its Roche lobe after shedding nearly half of its total mass. In contrast, the NS has gained approximately 0.2 M⊙. As the mass-transfer rate decreases, so does the ram pressure, allowing the magnetospheric boundary to extend outward. Initially, during the early stages of Roche lobe decoupling, the NS spin decreases while remaining in equilibrium. However, the reduction in the mass-transfer rate over time disrupts the equilibrium, allowing the NS to enter a relatively brief phase, lasting approximately 28 Myr, where its rotation rate further decreases. Eventually, the pulsar’s spin relaxes to an average value of ⟨Pns⟩ = 2.8 ms. Following this, the donor star refills its Roche lobe, initiating a second episode of mass transfer that lasts for approximately another ≃63 Myr. In this stage, the NS mass increases, accompanied by a steady contraction in radius, ultimately reaching approximate values of 1.3 M⊙ and 14 km, respectively, for M1. Simultaneously, the donor star transitions into a low-mass helium WD with Md ≃ 0.29 M⊙. In all cases, during the late stages of the evolution, the magnetospheric radius surpasses the light-cylinder radius. When this happens, the evolution of the NS is only driven by the magnetic dipole radiation torque (i.e., the third term of Eq. (11))

thumbnail Fig. 3.

Time evolution of M1 (top panel), M2 (middle panel), and M3 (bottom panel). The first column shows the evolution of mass-transfer rate (left y-axis) and NS radius (right y-axis). The second column shows the accretion torque. The third column shows the NS mass (left y-axis, dashed black line) and various radii (right y-axis: magnetic, corotation, light-cylinder, and NS radius). The fourth column shows the spin period (left y-axis) and orbital period (right y-axis). Gray regions indicate post-transition evolution, with non-solid tracks for M2 and M3. For a detailed description of the three spin evolution scenarios, see the main text.

In contrast, toward the end of this second mass-transfer phase, M2 surpasses the maximum mass that can be supported by pure hadronic matter. This event is visually represented in Fig. 1 when the green curve intersects the phase transition boundary. At this point, the NS mass is approximately 1.5 M⊙, and it has a radius of about 14.6 km, while the orbital period of the binary system is 45.2 d. To model the NS’s core transitioning from hadronic to quark matter due to accretion, we assumed a minor loss in gravitational mass of ΔM = 0.016 M⊙. This assumption requires the core’s instantaneous collapse, with both angular momentum (J) and baryonic mass (Mb) conserved, yielding a denser object, now with a reduced radius of approximately 13.7 km. This transition is indicated by the black arrow in Fig. 1, where its position on the third-family branch for compact objects is marked with a green star. For M2, the twin star is 6.2% more compact than its purely hadronic counterpart.

The structure of the donor star at the time of the phase transition is shown on the left panel of Fig. 4. The bottom half shows the helium core mass, which is approximately equal to 0.27 M⊙, and the perceived color of the star (as seen by an observer) based on its surface effective temperature. The inline plot illustrates the evolution of the donor star on the Hertzsprung-Russell diagram with the yellow circle marking its current position. The top half depicts the energy generation that takes place on a thin layer above the helium core with the color bar indicating the energy generation rate divided by the losses due to neutrino emission. The dominant mixing processes can also be discerned; at this point, the donor star has a fully convective envelope with a mass of approximately 0.1 M⊙ and a radius of ∼17 R⊙.

thumbnail Fig. 4.

Structure of the donor star. The top half of each figure shows the energy generation and interior mixing of the donor star at the time of the phase transition for M2 (left) and M3 (right). The bottom half shows the helium core mass and the color of the star based on its surface temperature. The yellow circle in the inline plot shows their position on a standard Hertzsprung-Russell diagram; the x-axis is the logarithm of the effective temperature, and the y-axis is the logarithm of surface luminosity. Dotted areas depict regions where convection is the dominant mixing mechanism.

Binary systems harboring rapidly rotating NSs in wider orbits, such as in M3, can experience a delayed accretion-induced phase transition as a result of the spin-down phase following the cessation of mass transfer (see the blue curve in Fig. 1). Owing to the longer initial orbital period, M3 experiences a single mass-transfer episode due to a reduced gravitational influence from the NS, allowing the donor star to maintain a more stable evolutionary path. As a consequence, the donor star does not expand significantly until it reaches a later stage in its evolution, where it fills its Roche lobe only once. The mass transfer commences approximately 104 Myr later compared to the first RLOF for M1 and M2 (bottom panel in Fig. 3). The duration of the mass transfer is close to 33 Myr and, at the end of it, the donor and the NS possess a mass of approximately 0.33 M⊙ and 1.56 M⊙, respectively. Subsequent to Roche lobe detachment, the NS spin rate decreases for approximately 117 kyr at which point becomes unstable and undergoes a phase transition. Similarly to M2, its trajectory and final position on the third-family branch are marked with a black arrow and a blue star, respectively. The resulting twin star is ∼8.4% more compact than the initial NS, settling on a radius of approximately 13.9 km.

The structure of the donor star at the time of the phase transition can be seen in Fig. 4. Although similar to the structure of the donor star in M2, the envelope mass is lower (∼0.04 M⊙) and its total radius almost twice that of M2 (∼34 R⊙) as it enters the asymptotic giant branch phase.

3.2. Orbital reconfiguration and post-transition evolution

Figure 5 illustrates the post-transition binary system configurations for our progenitor models, M2 and M3, for kick velocities with fixed magnitudes and random orientations imparted on the compact star. For both M2 and M3, even in the absence of a kick (w = 0 km s−1; depicted by the black markers), the binary systems exhibit an increase in their eccentricity that is directly proportional to the mass defect during the mass transition, e = ΔM/M ≃ 7 × 10−2, for ΔM = 0.016 M⊙. For M3, these orbital changes are preserved, as the system does not experience additional mass transfer after the phase transition, and thus the tidal forces are negligible (Freire & Tauris 2014; Antoniadis 2014; Jiang et al. 2015). When factoring in a secondary kick mechanism, as the kick velocity increases, there is a progressively larger spread in the distribution of both eccentricity and binary separation, leading to a more diverse range of possible post-kick orbits. For a kick of certain magnitude, wider systems such as M3 are more likely to produce more eccentric systems or eject the compact star, resulting in the formation of an eMSP or iMSP, respectively.

thumbnail Fig. 5.

Monte Carlo simulation for the post-transition distribution of orbital parameters from M2 (left) and M3 (right). A transition-triggered mass loss of ΔM = 0.016 M⊙ and a fixed magnitude kick velocity were assumed. The kick orientation was randomly sampled from a uniform distribution.

Evolutionary calculations following the phase transition, provided by MESA, in M2 may be unreliable as they do not take into account the effect of the orbital reconfiguration on the mass transfer. In contrast, the simulation for M3, where the phase transition occurs during Roche lobe decoupling, remains generally reliable. However, it is important to note that the final orbital period in M3 may still be inaccurate due to the limitations in capturing the orbital reconfiguration. In addition, the radius of the hybrid star decreases as a result of the phase transition. Due to the conservation of angular momentum, this reduction leads to a 15% increase in the spin frequency of the compact star. Moreover, because of the magnetic flux conservation during the phase transition, the surface magnetic field of the hybrid star also increases. The new magnetic field strength can be calculated as

B hybrid = B ns ( R ns R hybrid ) 2 ≃ 1.14 × 10 8 G , $$ \begin{aligned} B_\mathrm{hybrid} = B_\mathrm{ns} \left(\frac{R_\mathrm{ns} }{R_\mathrm{hybrid} }\right)^2 \simeq 1.14 \times 10^8\,\mathrm{G} , \end{aligned} $$(15)

where Bns and Rns are the magnetic field and the radius of the NS before the transition, respectively, and Rhybrid is the radius of the hybrid star after the transition.

For M2 where the phase transition takes place during the LMXB phase, the magnetospheric radius outgrows the corotation radius (i.e., the fastness parameter becomes greater than one) forcing the hybrid star to enter a propeller phase and inhibit further accretion. Hence, it is conceivable that the NS mass stops increasing despite the ongoing mass transfer. Moreover, if the magnetospheric radius outgrows the light cylinder radius, the magnetic field is strong enough to prevent accretion altogether. This may lead to the formation of a binary system that consists of a rotation-powered pulsar and an evolved “spider”-like companion (see Strader et al. 2019; Papitto & de Martino 2022, and references therein).

Although we did not perform detailed calculations with MESA for the post-transition evolution of M2 and M3, we used the same accretion model to compute the mass, radius, and spin evolution of the hybrid stars (dash double dotted lines in Fig. 3). We find that M2 continues to go through sporadic episodes of mass accretion reaching a final mass of approximately 1.51 M⊙.

In contrast, M3, which had already decoupled from its Roche lobe prior to the strong phase transition, retained its post-transition mass of approximately 1.55 M⊙. Following the cessation of mass transfer, both models are further compactified as a result of the spin-down and angular momentum loss (see also Fig. 1). By the end of the evolution, the hybrid star of M2 has a radius of approximately 12.4 km compared to its pure hadronic counterpart with a radius of 14.6 km at the time of transition. Similarly, in M3, the hybrid star has a radius of 12.2 km whereas the pure hadronic star has a radius of 15.2 km.

4. Discussion

4.1. Summary

In this work we investigated the potential role of rapid first-order phase transitions from ordinary baryonic to quark matter in the evolution of LMXB and post-LMXB systems. We performed detailed binary evolution modeling using MESA and Monte Carlo simulations to explore the influence of mass transfer, accretion- and spin-down-induced phase transitions, and secondary kick mechanisms on the formation and evolution of these systems. Our results demonstrate that:

  1. Twin stars can form as a result of mass transfer in LMXBs

  2. In systems with sufficiently long orbital periods, strong phase transitions can occur after the LMXB phase, when the NS is a rotation-powered MSP. In this case, the eccentricity of the resulting system can be several orders of magnitude larger than what is expected for recycled MSPs (Phinney 1992).

  3. If the transition is accompanied by a secondary kick (see Sect. 2.3 for the justification), then a considerable number of systems can be disrupted, even if the kick magnitude is moderate. Therefore, this mechanism may be a viable formation path for iMSPs.

  4. Secondary kicks also provide a mechanism for producing binary MSPs with very long orbital periods (> 1000 d) from systems that are initially more compact.

4.2. Formation of isolated and eccentric millisecond pulsars

As discussed in Sect. 3.2, the minimum post-transition eccentricity is directly related to the mass defect, ΔM. The EoS model adopted in this work predicts a mass defect of ΔM = 0.016 M⊙, which in turn produces eccentricities of e ≤ ΔM/M ≃ 7 × 10−2.

For M2 this residual eccentricity would rapidly dissipate as a result of tidal forces acting on the envelope of the donor star (Zahn 1977). In contrast, for M3 the strong phase transition occurs during the spin-down phase (i.e., after the LMXB phase), when the donor star is already evolving toward the WD cooling track. For this reason, the resulting orbital changes would persist (Freire & Tauris 2014; Antoniadis 2014), providing a potential mechanism for explaining eMSPs. Although this particular model has an orbital period that is too long compared to those of known eMSPs with WD companions, it shows that strong phase transitions provide a possible formation mechanism for eMSPs.

Because high-precision pulsar timing can constrain extremely small eccentricities (≤10−8), even models with no neutrino trapping and with small mass defects that would produce smaller deviations from the Phinney (1992) relation could potentially be observed. Our work suggests that the MSP eccentricity distribution, especially in binary MSPs with long orbital periods, can be used to constrain not only stellar physics (Phinney 1992; Antoniadis 2014), but also nuclear matter models.

When considering secondary kick mechanisms, higher kick velocities produce a broader distribution of eccentricities and orbital periods, suggesting a diverse range of possible post-transition configurations. Figures 5 and 6 illustrate that both M2 and M3 have a high probability of producing iMSPs or binary MSPs with extremely long orbital periods. More specifically, the red and blue curves in Fig. 6 (representing iMSPs from M2 and M3, respectively) illustrate the expected fraction of binary systems that, due to the imparted kick velocities, achieve eccentricities greater than unity – indicating the complete disruption of the binary and the formation of an iMSP. As the kick velocity increases, a larger fraction of systems become unbound. In contrast, the black curve on the left y-axis, representing the fraction of eMSPs in M3, shows an inverse relationship: as the kick velocity increases, the fraction of eMSPs declines. This suggests that higher kick velocities are more likely to disrupt the binary systems entirely, thereby preventing the formation of eMSPs.

thumbnail Fig. 6.

Expected fraction of eMSPs (black curve, left y-axis) and iMSPs (red and blue curves, right y-axis) as a function of kick velocity, assuming a mass defect of ΔM = 0.016 M⊙. For iMSPs, we consider the systems with eccentricities greater than unity. M2 cannot create eccentric systems due to the rapid re-circularization of the orbit.

In our study we considered three different mechanisms that could introduce such a kick: the neutrino rocket effect, the electromagnetic rocket effect, and the GW recoil scenario. In the last case, even if a small percentage (∼0.5%) of the converted energy is emitted in the form of GWs, it can account for the desired kick velocities. However, a fully relativistic approach to modeling the GW emission is necessary to accurately quantify the results of such a recoil. Nevertheless, the possibility that these three mechanisms could coexist, or sequentially influence the binary during various evolutionary stages, is intriguing since each contributes differently to the system dynamics. Future pulsar surveys that are more sensitive to these systems will therefore play an important role in constraining the potential birth signatures of twin stars (Watts et al. 2015; Antoniadis 2021).

Acknowledgments

This work was supported by the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of “Science and Society” Action Always strive for excellence – “Theodoros Papazoglou” (Project Number: 01431). S.Ch and J.A. also acknowledge support from the European Commission under project ARGOS-CDS (Grant Agreement number: 101094354). D.B. acknowledges support from the Polish National Science Center under grant No. 2019/33/B/ST9/03059 (before 10/2023) and under grant No. 2021/43/P/ST2/03319 (after 04/2024). D.E.A-C. is supported by the program Excellence Initiative-Research University of the University of Wroclaw of the Ministry of Education and Science. This work made extensive use of the following software: Astropy (Astropy Collaboration 2013, 2018), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Pandas (The pandas development team 2020), Matplotlib (Hunter 2007), SciencePlots (Garrett 2021), MesaReader (Bill & Josiah 2017), Jupyter (Kluyver et al. 2016). Figure 4 was created with TULIPS (https://astro-tulips.readthedocs.io/en/latest/) (Laplace 2022).

References

  1. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Phys. Rev. Lett., 119, 161101 [Google Scholar]
  2. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101 [Google Scholar]
  3. Agalianou, V., & Gourgouliatos, K. N. 2023, MNRAS, 522, 5879 [NASA ADS] [CrossRef] [Google Scholar]
  4. Aguilera, D. N., Blaschke, D., & Grigorian, H. 2004, A&A, 416, 991 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Alvarez-Castillo, D. E. 2021, Astron. Nachr., 342, 234 [NASA ADS] [CrossRef] [Google Scholar]
  6. Alvarez-Castillo, D. E., Antoniadis, J., Ayriyan, A., et al. 2019, Astron. Nachr., 340, 878 [NASA ADS] [CrossRef] [Google Scholar]
  7. Andersson, N., Glampedakis, K., Haskell, B., & Watts, A. L. 2005, MNRAS, 361, 1153 [NASA ADS] [CrossRef] [Google Scholar]
  8. Antoniadis, J. 2014, ApJ, 797, L24 [NASA ADS] [CrossRef] [Google Scholar]
  9. Antoniadis, J. 2021, MNRAS, 501, 1116 [Google Scholar]
  10. Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448 [Google Scholar]
  11. Antoniadis, J., Kaplan, D. L., Stovall, K., et al. 2016, ApJ, 830, 36 [NASA ADS] [CrossRef] [Google Scholar]
  12. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  14. Barr, E. D., Champion, D. J., Kramer, M., et al. 2013, MNRAS, 435, 2234 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bauswein, A., Blaschke, D., & Fischer, T. 2022, Astrophysics in the XXI Century with Compact Stars (Singapore: World Scientific), 283 [CrossRef] [Google Scholar]
  16. Baym, G., Pethick, C., & Sutherland, P. 1971, ApJ, 170, 299 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bejger, M., Blaschke, D., Haensel, P., Zdunik, J. L., & Fortin, M. 2017, A&A, 600, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Benic, S., Blaschke, D., Alvarez-Castillo, D. E., Fischer, T., & Typel, S. 2015, A&A, 577, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Berdermann, J., Blaschke, D., Grigorian, H., & Voskresensky, D. N. 2006, Prog. Part. Nucl. Phys., 57, 334 [NASA ADS] [CrossRef] [Google Scholar]
  20. Bhattacharya, D., & van den Heuvel, E. P. J. 1991, Phys. Rep., 203, 1 [Google Scholar]
  21. Bhattacharyya, S. 2023, Galaxies, 11, 103 [NASA ADS] [CrossRef] [Google Scholar]
  22. Bhattacharyya, S., & Chakrabarty, D. 2017, ApJ, 835, 4 [NASA ADS] [CrossRef] [Google Scholar]
  23. Bildsten, L., Chakrabarty, D., Chiu, J., et al. 1997, ApJS, 113, 367 [CrossRef] [Google Scholar]
  24. Bill, W., & Josiah, S. 2017, https://doi.org/10.5281/zenodo.826958 [Google Scholar]
  25. Bogdanov, S., et al. 2021, ApJ, 914, L15 [NASA ADS] [CrossRef] [Google Scholar]
  26. Bombaci, I., & Popov, S. B. 2004, A&A, 424, 627 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Camilo, F., Kerr, M., Ray, P. S., et al. 2015, ApJ, 810, 85 [NASA ADS] [CrossRef] [Google Scholar]
  28. Carlomagno, J. P., Contrera, G. A., Grunfeld, A. G., & Blaschke, D. 2024a, Universe, 10, 336 [NASA ADS] [CrossRef] [Google Scholar]
  29. Carlomagno, J. P., Contrera, G. A., Grunfeld, A. G., & Blaschke, D. 2024b, Phys. Rev. D, 109, 043050 [CrossRef] [Google Scholar]
  30. Chaboyer, B., & Zahn, J. P. 1992, A&A, 253, 173 [NASA ADS] [Google Scholar]
  31. Champion, D. J., Ransom, S. M., Lazarus, P., et al. 2008, Science, 320, 1309 [NASA ADS] [CrossRef] [Google Scholar]
  32. Chandra, A. D. 2025, PASA, accepted [arXiv:2311.13303] [Google Scholar]
  33. Chevalier, R. A. 1989, ApJ, 346, 847 [Google Scholar]
  34. Corbet, R. H. D. 1996, ApJ, 457, L31 [NASA ADS] [CrossRef] [Google Scholar]
  35. Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020, Nat. Astron., 4, 72 [NASA ADS] [CrossRef] [Google Scholar]
  36. Cumming, A., Zweibel, E., & Bildsten, L. 2001, ApJ, 557, 958 [NASA ADS] [CrossRef] [Google Scholar]
  37. Deneva, J. S., Stovall, K., McLaughlin, M. A., et al. 2013, ApJ, 775, 51 [NASA ADS] [CrossRef] [Google Scholar]
  38. Elsner, R. F., Ghosh, P., & Lamb, F. K. 1980, ApJ, 241, L155 [NASA ADS] [CrossRef] [Google Scholar]
  39. Freire, P. C. C., & Tauris, T. M. 2014, MNRAS, 438, L86 [NASA ADS] [CrossRef] [Google Scholar]
  40. Freire, P. C. C., Bassa, C. G., Wex, N., et al. 2011, MNRAS, 412, 2763 [NASA ADS] [CrossRef] [Google Scholar]
  41. Fukushima, K., & Yu, C. 2024, arXiv e-prints [arXiv:2401.04568] [Google Scholar]
  42. Galloway, D. K., Chakrabarty, D., Morgan, E. H., & Remillard, R. A. 2002, ApJ, 576, L137 [NASA ADS] [CrossRef] [Google Scholar]
  43. Garrett, J. D. 2021, http://doi.org/10.5281/zenodo.4106649 [Google Scholar]
  44. Gerlach, U. H. 1968, Phys. Rev., 172, 1325 [NASA ADS] [CrossRef] [Google Scholar]
  45. Ginzburg, S., & Chiang, E. 2022, MNRAS, 509, L1 [Google Scholar]
  46. Glebbeek, E., Gaburov, E., de Mink, S. E., Pols, O. R., & Portegies Zwart, S. F. 2009, A&A, 497, 255 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Glendenning, N. K., Pei, S., & Weber, F. 1997, Phys. Rev. Lett., 79, 1603 [NASA ADS] [CrossRef] [Google Scholar]
  48. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  49. Harrison, E. R., & Tademaru, E. 1975, ApJ, 201, 447 [NASA ADS] [CrossRef] [Google Scholar]
  50. Hebeler, K., Lattimer, J. M., Pethick, C. J., & Schwenk, A. 2013, ApJ, 773, 11 [CrossRef] [Google Scholar]
  51. Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
  52. Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 [Google Scholar]
  53. Hempel, M., Heinimann, O., Yudin, A., et al. 2016, Phys. Rev. D, 94, 103001 [NASA ADS] [CrossRef] [Google Scholar]
  54. Henyey, L., Vardya, M. S., & Bodenheimer, P. 1965, ApJ, 142, 841 [Google Scholar]
  55. Hills, J. G. 1983, ApJ, 267, 322 [NASA ADS] [CrossRef] [Google Scholar]
  56. Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974 [Google Scholar]
  57. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  58. Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
  59. Jiang, L., Li, X.-D., Dey, J., & Dey, M. 2015, ApJ, 807, 41 [NASA ADS] [CrossRef] [Google Scholar]
  60. Jiang, L., Wang, N., Chen, W.-C., et al. 2020, A&A, 633, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Jiang, L., Wang, N., Chen, W.-C., et al. 2021, RAA, 21, 231 [NASA ADS] [Google Scholar]
  62. Kalapotharakos, C., Wadiasingh, Z., Harding, A. K., & Kazanas, D. 2021, ApJ, 907, 63 [NASA ADS] [CrossRef] [Google Scholar]
  63. Kaminski, M., Uhlemann, C. F., Bleicher, M., & Schaffner-Bielich, J. 2016, Phys. Lett. B, 760, 170 [CrossRef] [Google Scholar]
  64. Khosravi Largani, N., Fischer, T., Shibagaki, S., Cerdá-Durán, P., & Torres-Forné, A. 2024, A&A, 687, A245 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, eds. F. Loizides, & B. Schmidt, (IOS Press), 87 [Google Scholar]
  66. Knispel, B., Lyne, A. G., Stappers, B. W., et al. 2015, ApJ, 806, 140 [NASA ADS] [CrossRef] [Google Scholar]
  67. Kolb, U., & Ritter, H. 1990, A&A, 236, 385 [NASA ADS] [Google Scholar]
  68. Konar, S., & Choudhuri, A. R. 2004, MNRAS, 348, 661 [NASA ADS] [CrossRef] [Google Scholar]
  69. Langer, N. 1991, A&A, 252, 669 [NASA ADS] [Google Scholar]
  70. Laplace, E. 2022, Astron. Comput., 38, 100516 [NASA ADS] [CrossRef] [Google Scholar]
  71. Ledoux, P. 1947, ApJ, 105, 305 [NASA ADS] [CrossRef] [Google Scholar]
  72. Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993 [Google Scholar]
  73. Miller, M. C., et al. 2019, ApJ, 887, L24 [CrossRef] [Google Scholar]
  74. Miller, M. C., et al. 2021, ApJ, 918, L28 [NASA ADS] [CrossRef] [Google Scholar]
  75. Mishustin, I. N., Hanauske, M., Bhattacharyya, A., et al. 2003, Phys. Lett. B, 552, 1 [CrossRef] [Google Scholar]
  76. Misra, D., Fragos, T., Tauris, T. M., Zapartas, E., & Aguilera-Dena, D. R. 2020, A&A, 642, A174 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Nugis, T., & Lamers, H. J. G. L. M. 2000, A&A, 360, 227 [NASA ADS] [Google Scholar]
  78. Nurmamat, N., Zhu, C., Lü, G., et al. 2019, JA&A, 40, 32 [Google Scholar]
  79. Octau, F., Cognard, I., Guillemot, L., et al. 2018, A&A, 612, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Papitto, A., & de Martino, D. 2022, in Astrophysics and Space Science Library, eds. S. Bhattacharyya, A. Papitto, & D. Bhattacharya, Astrophys. Space Sci. Lib., 465, 157 [NASA ADS] [CrossRef] [Google Scholar]
  81. Paschalidis, V., Yagi, K., Alvarez-Castillo, D., Blaschke, D. B., & Sedrakian, A. 2018, Phys. Rev. D, 97, 084038 [NASA ADS] [CrossRef] [Google Scholar]
  82. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  83. Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
  84. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  85. Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
  86. Phinney, E. S. 1992, Phil. Trans. R. Soc. London Ser. A, 341, 39 [CrossRef] [Google Scholar]
  87. Portegies Zwart, S., van den Heuvel, E. P. J., van Leeuwen, J., & Nelemans, G. 2011, ApJ, 734, 55 [NASA ADS] [CrossRef] [Google Scholar]
  88. Pradhan, B. K., Chatterjee, D., & Alvarez-Castillo, D. E. 2024, MNRAS, 531, 4640 [NASA ADS] [CrossRef] [Google Scholar]
  89. Rappaport, S., Verbunt, F., & Joss, P. C. 1983, ApJ, 275, 713 [Google Scholar]
  90. Rappaport, S. A., Fregeau, J. M., & Spruit, H. 2004, ApJ, 606, 436 [NASA ADS] [CrossRef] [Google Scholar]
  91. Reimers, D. 1975, Problems in Stellar Atmospheres and Envelopes (Springer-Verlag, New York, Inc.), 229 [CrossRef] [Google Scholar]
  92. Riley, T. E., et al. 2019, ApJ, 887, L21 [CrossRef] [Google Scholar]
  93. Sandin, F., & Blaschke, D. 2007, Phys. Rev. D, 75, 125013 [CrossRef] [Google Scholar]
  94. Schertler, K., Greiner, C., Schaffner-Bielich, J., & Thoma, M. H. 2000, Nucl. Phys. A, 677, 463 [CrossRef] [Google Scholar]
  95. Shibazaki, N., Murakami, T., Shaham, J., & Nomoto, K. 1989, Nature, 342, 656 [NASA ADS] [CrossRef] [Google Scholar]
  96. Spruit, H. C. 2002, A&A, 381, 923 [CrossRef] [EDP Sciences] [Google Scholar]
  97. Stovall, K., Freire, P. C. C., Antoniadis, J., et al. 2019, ApJ, 870, 74 [NASA ADS] [CrossRef] [Google Scholar]
  98. Strader, J., Swihart, S., Chomiuk, L., et al. 2019, ApJ, 872, 42 [Google Scholar]
  99. Sun, S., Li, L., Liu, H., et al. 2019, PASA, 36, e005 [NASA ADS] [CrossRef] [Google Scholar]
  100. Tauris, T. M. 2012, Science, 335, 561 [NASA ADS] [CrossRef] [Google Scholar]
  101. Tauris, T. M., & Savonije, G. J. 1999, A&A, 350, 928 [NASA ADS] [Google Scholar]
  102. Tauris, T. M., & van den Heuvel, E. P. J. 2006, Compact stellar X-ray sources, 39 (Cambridge, UK: Cambridge University Press), 623 [NASA ADS] [CrossRef] [Google Scholar]
  103. Tauris, T. M., & van den Heuvel, E. P. J. 2023, Physics of Binary Star Evolution. From Stars to X-ray Binaries and Gravitational Wave Sources (Princeton University Press) [Google Scholar]
  104. Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170 [Google Scholar]
  105. The pandas development team 2020, https://doi.org/10.5281/zenodo.3509134 [Google Scholar]
  106. Verbunt, F., & Phinney, E. S. 1995, A&A, 296, 709 [NASA ADS] [Google Scholar]
  107. Verbunt, F., van den Heuvel, E. P. J., van Paradijs, J., & Rappaport, S. A. 1987, Nature, 329, 312 [NASA ADS] [CrossRef] [Google Scholar]
  108. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  109. Watts, A., Espinoza, C. M., Xu, R., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array, AASKA14, 43 [Google Scholar]
  110. Zahn, J. P. 1977, A&A, 57, 383 [Google Scholar]
  111. Zdunik, J. L., Bejger, M., Haensel, P., & Gourgoulhon, E. 2006, A&A, 450, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1.

Initial binary parameters for our stellar evolution models.

In the text

All Figures

thumbnail Fig. 1.

Left Panel: Baryonic mass vs. equatorial radius. Dotted and dot-dashed curves indicate lines of constant angular momentum and constant frequency, respectively. The dashed black line highlights the points of maximum stability in the hadronic branches, signaling the onset of a phase transition (PT). Thin, horizontal magenta lines show trajectories where angular momentum (J) and baryonic mass (Mb) are conserved. The thick magenta curve connects points on the hybrid star branches that can be reached through a collapse conserving both J and Mb. Right Panel: Similar to the left panel, but with the y-axis representing gravitational mass. Star markers denote the endpoint trajectories for direct (M2) and delayed (M3) collapse models. The arrows are inclined because gravitational mass, unlike baryonic mass, is not conserved during the PT.
In the text
thumbnail Fig. 2.

Baryonic mass (Mb) and gravitational mass (Mg) as functions of the central energy density (εc) for a fixed angular momentum (J). The solid black line represents the gravitational mass, while the dashed black line represents the baryonic mass. Vertical gray lines mark the central energy densities at the onset (εc = 275 MeV/fm3) and completion (εc = 595 MeV/fm3) of the phase transition (PT). Black dotted lines and the blue-shaded region denote unstable configurations. The inset zooms into the PT region, illustrating the mass defect (ΔMg) of 0.016 M⊙ associated with the transition from hadronic to quark matter. The horizontal gray lines indicate constant values of Mb and Mg during the PT. The star marker denotes the birth mass of the twin star configuration at Mg = 1.499 M⊙.
In the text
thumbnail Fig. 3.

Time evolution of M1 (top panel), M2 (middle panel), and M3 (bottom panel). The first column shows the evolution of mass-transfer rate (left y-axis) and NS radius (right y-axis). The second column shows the accretion torque. The third column shows the NS mass (left y-axis, dashed black line) and various radii (right y-axis: magnetic, corotation, light-cylinder, and NS radius). The fourth column shows the spin period (left y-axis) and orbital period (right y-axis). Gray regions indicate post-transition evolution, with non-solid tracks for M2 and M3. For a detailed description of the three spin evolution scenarios, see the main text.
In the text
thumbnail Fig. 4.

Structure of the donor star. The top half of each figure shows the energy generation and interior mixing of the donor star at the time of the phase transition for M2 (left) and M3 (right). The bottom half shows the helium core mass and the color of the star based on its surface temperature. The yellow circle in the inline plot shows their position on a standard Hertzsprung-Russell diagram; the x-axis is the logarithm of the effective temperature, and the y-axis is the logarithm of surface luminosity. Dotted areas depict regions where convection is the dominant mixing mechanism.
In the text
thumbnail Fig. 5.

Monte Carlo simulation for the post-transition distribution of orbital parameters from M2 (left) and M3 (right). A transition-triggered mass loss of ΔM = 0.016 M⊙ and a fixed magnitude kick velocity were assumed. The kick orientation was randomly sampled from a uniform distribution.
In the text
thumbnail Fig. 6.

Expected fraction of eMSPs (black curve, left y-axis) and iMSPs (red and blue curves, right y-axis) as a function of kick velocity, assuming a mass defect of ΔM = 0.016 M⊙. For iMSPs, we consider the systems with eccentricities greater than unity. M2 cannot create eccentric systems due to the rapid re-circularization of the orbit.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.