THE FORMATION OF THE ECCENTRIC-ORBIT MILLISECOND PULSAR J1903+0327 AND THE ORIGIN OF SINGLE MILLISECOND PULSARS

The orbits of the triple stars were integrated using a regularized version of the Burlish-Stoer integrator (Aarseth & Zare 1974a, 1974b), keeping the numerical error at machine precision, and allowing a maximum relative energy error of O (10⁻¹⁴) per integration of the outer orbit.

During the dynamical evolution, we resolved the mass transfer and consequential change in orbital parameters of the LMXB that is orbited by the outer star.

The rate of mass transfer in the inner binary is calculated from the empirical relation based on the orbital period at the onset of RLOF in the LMXB: $\langle \dot{m} \rangle = 6 \times 10^{-10} (\mbox{${P}_{\rm orb}$}\,({\rm initial})/1\, {\rm day})$ M_☉ yr⁻¹ (Shore et al. 1994).

The NS was allowed to accrete at most at the Eddington rate and any surplus mass is assumed to leave the binary with the angular momentum of the accreting NS, but lost adiabatically from the triple.

We perform the mass transfer in the inner binary every time the outer orbit has had 10 revolutions, after which the numerical orbit integration was updated using the newly calculated orbital parameters which resulted from carrying out the mass transfer. We varied the interval between which mass transfer in the inner binary was conducted between every 1–1000 outer orbits, but this choice did not significantly affect the results. (In Section 5.3, we present the results of a series of simulations where we decoupled the mass-transfer process from the gravitational evolution after the triple has become dynamically unstable.)

We initialized 10³ binaries and calculated the evolution for each up to an age of at most 10 Gyr or until the mass of the donor star drops to the mass of the degenerate helium core for a population-II star (≲ 0.4 M_☉; Tauris & Savonije 1999).

Each triple was initialized by randomly selecting the eccentricity of the outer orbit from the thermal distribution and the initial separations of the inner and outer orbits for each triple-evolution calculation are selected using the iterative procedure described in Section 3 (see the thick solid and dashed curves in Figure 1). For each simulated binary, we randomly selected the inclination of the inner orbit with respect to the outer orbit, the longitude of the ascending node, the argument of periastron, and the phases of the two orbits. The preference in the orbital elements introduced by the supernova explosion may affect the rate of ejected MSPs relative to those that stay in a binary, but we ignore that complication. We found no significant correlations between the final outcome of the simulations and longitude of the ascending node, the argument of periastron, or the phases of the two orbits. The anisotropic velocity caused by the mass loss in the supernova of the inner binary is therefore not expected to have a significant effect on the surviveability of the triple. During the evolution of the triples, small eccentricities (≲ 0.1) are commonly induced in the inner orbit, particularly when the triple approaches the regime where it becomes dynamically unstable.

We stop a simulation when the eccentricity of the inner (LMXB) orbit exceeds 0.3 for more than 10⁵ years. Such high eccentricities can be induced shortly by a strong interaction with the outer star, which typically results in the break-up of the triple, or by a long term secular perturbations of the inner orbit by the outer star (Kozai 1979). In the latter case, orbital variations result naturally from the secular evolution of the triple and do not directly lead to a dynamical instability. We still decided to terminate such simulations because it becomes hard to follow the mass-transfer process within the inner binary. It would require extensive hydrodynamical simulations to study the consequences of a Kozai resonances in a triple with a Roche lobe-filling inner binary. About 10% of our simulations were stopped as a consequence of this effect. Mass transfer in eccentric orbits is generally ill understood, although courageous attempts are underway to tackle this problem (Lajoie & Sills 2011a, 2011b).

The majority (917) of the binaries become dynamically unstable long before the other stopping criteria apply, in which cases we continue to resolve the dynamics by integrating the equations of motion and resolving the internal mass transfer until one of the stars escapes. The orbital separation of the LMXB at which the triple is expected to become unstable is indicated by the dotted curve in Figure 1.

An illustrative example of the evolution of the period of the inner and outer orbits is presented in Figure 2. The evolution of the inner LMXB also drives the expansion of the outer orbit, in particular by the dynamical coupling between both orbits and in a lesser extend by the mass lost from the inner LMXB in those cases that mass transfer proceeds non-conservatively. This relative softening causes the final MSP binaries to be somewhat wider than observed in J1903+0327, and as a consequence the triple remains stable for somewhat longer than expected based on our analytic energy balance (see Section 3), and eventually results in the MSP being somewhat more massive (by about 0.2 M_☉) than observed in J1903+0327. We can compensate for this by adopting a slightly (∼20 %) smaller outer orbital separation at the onset of the LMXB phase.

Zoom In Zoom Out Reset image size

The resulting parameters of the MSP with a G-star binary are presented in Figure 3; they straddle the observed orbital separation and eccentricity of J1903+0327. The phase of mass transfer for these binaries lasted for about 23.4 ± 9.3 Myr, after which the donor was ejected in about one-fifth ∼20% of the cases. The resulting binaries, for those with a < 10⁴ R_☉, had an average orbital separation of 175 ± 145 R_☉, and an eccentricity of 0.65 ± 0.19. During the LMXB phase, the NSs were able to grow from 1.30 M_☉ to 1.64 ± 0.12 M_☉. The mean mass of the ejected donor was 0.72 ± 0.13 M_☉.

Zoom In Zoom Out Reset image size

Our success in reproducing the observed parameters of the binary MSP J1903+0327 demonstrates that its progenitor may well have been born as a triple star. However, the pin-pointed search of parameter space described here makes it impossible to derive the birthrate for J1903+0327, which we will calculate in the next section.

According to the Hipparcos database, the ratio of hierarchical, higher order, multiple stellar systems (404) to binaries (1438) is 404/1842 ∼ 0.22 and the outer stars have a rather flat mass distribution (Eggleton & Tokovinin 2008, 2010). We assume that the same ratio holds for (inner) binaries that are progenitors to LMXBs, even though there are none in this catalogue. We further require that the mass of the outer stars is about ten times lower than that of the inner binary. We then find a ratio of suitable triples to LMXB progenitor binaries of a few percent. We present in Section 6.2 observational evidence that such triples do exist.

We will not dwell on the details of the common-envelope evolution, but assume that it leads to the spiral-in of the inner two components without significantly affecting the outer star. The consequences of the common envelope do not qualitatively affect our result, but have a strong effect on the derived birthrate (see Section 5.4). In Section 6.2 we discuss that the orbit of the tertiary star that orbits the LMXB 4U 2129+47 may have experienced a considerable reduction during the common envelope, but we consider this insufficient evidence to draw general conclusions regarding the effect of the common envelope on the outer orbit.

The effect of the supernova explosion in the inner binary on the orbital parameters of the triple can be calculated relatively straightforwardly, and has a profound effect on the survivability of the triple because the weakly bound outer orbit is easily disrupted.

The mass loss in the supernova explosion causes the inner binary to be ejected (Blaauw 1961; Boersma 1961). In order to keep the triple bound, a small asymmetric velocity kick imparted to the newly formed NS is required (Hills 1983; Brandt & Podsiadlowski 1995; Tauris & Takens 1998). In Figure 4, we show the survival probability of several combinations of inner and outer orbits as function of the asymmetric kick magnitude. In these calculations, the effect of the Blaauw–Boersma kick was self-consistently taken into account to calculate the survival probability. This fraction ranges from ∼40% for small kicks to zero for kicks significantly above $100\,\mbox{${\rm km\, s}^{-1}$}$.

Zoom In Zoom Out Reset image size

To study the future evolution of the triple through the LMXB phase, we need to derive the orbital parameters as a function of the most probable kick velocity and the amount of mass lost in the supernova. We consider several combinations of both, and the distinction between them can be made by the mass of the initial primary star and the moment it filled its Roche lobe.

In Section 3, we discussed the most likely range of parameters at birth and adopt these same numbers here to study the survivability of the triple. Here we make a distinction between two types of supernovae; in one case we adopt that the inner primary at the moment of the supernova is a 1.9 M_☉ helium star with a degenerate ONeMg core that is the leftover of the initial 9–12 M_☉ primary star. If the initial primary star was born somewhat more massive, 10–13 M_☉, and was stripped from its hydrogen envelope at a later stage of its evolution, its helium core may have grown to ∼2.7 M_☉. The ONeMg star that experiences an electron-capture supernova loses ∼0.6 M_☉ and if the exploding star has a more massive helium core, the mass loss is ∼1.4 M_☉. In both cases, we adopt an NS mass of 1.30 M_☉.

Apart from the smaller mass lost in the explosion of the ONeMg core, the NS is also expected to receive a smaller velocity kick upon birth. For these electron capture supernovae, we adopted a Gaussian distribution for the velocity kick with a dispersion of 20 km s⁻¹ (Podsiadlowski et al. 2004; Scheck et al. 2004) in a random direction. From studies of single radio pulsars in the solar neighborhood, several kick velocity distributions have been constructed (Lyne & Lorimer 1994; Hansen & Phinney 1997; Cordes & Chernoff 1998; Arzoumanian et al. 2002), most of which yield considerably higher velocity than for the ONeMg supernovae. These latter kick velocity distributions are expected to be more suitable for NSs formed from an isolated star or a non-interacting binary, in which case the NS is formed by the collapse of an iron core (van den Heuvel 2004; Podsiadlowski et al. 2004). These kicks range from single (zero-centered) Gaussian distributions with a velocity dispersion of 265 km s⁻¹ (Hobbs et al. 2005) to more complicated distributions like the proposed two Gaussians with dispersions of 175 km s⁻¹ and 700 km s⁻¹ and relative probability of 0.86 and 0.14, respectively (Cordes & Chernoff 1998).

We study the probability that a triple survives the supernova by means of Monte Carlo simulations. The simulated triples had the following characteristics: the inner binary consists of a 1.9 M_☉ helium star with a degenerate ONeMg core or a ∼2.7 M_☉ helium star, and a 0.8–2.0 M_☉ secondary. The latter was selected randomly with equal probability within the interval. The orbital separation of the circular inner binary was taken to be flat in log  space between 1 R_☉ and 100 R_☉. We adopted these parameters from the population synthesis calculation of LMXBs (Willems & Kolb 2003), in particular using the results of their models KM25 to KM100 of (Willems & Kolb 2003).

The mass of the tertiary star was selected to be less massive than the secondary, but not smaller than 0.8 M_☉ (in theory there is no lower limit for the mass of the tertiary star). The outer orbital separation was chosen with a probability distribution flat in log  space with a maximum of 10⁴ R_☉. The minimum separation was chosen to be consistent with a dynamically stable initial triple (adopting a 10 M_☉ primary and an inner orbital separation of 200 R_☉) and the earlier selected secondary and tertiary masses. The latter two stars were assumed not to accrete any material throughout the common envelope and supernova explosion. The eccentricity of the outer orbit before the supernova was selected at random from the thermal distribution between a circular orbit and a maximum which was chosen such that the triple is dynamical stable. The selection of the minimum orbital separation of the outer star and its eccentricity therewith becomes an iterative procedure. The other orbital elements are selected randomly, as we described in Section 4.2.

For each system we calculate the effect of the combined Blaauw–Boersma and intrinsic velocity kick on the inner and the outer orbits; the latter kick was assumed to be isotropic. In our simulations, we varied the mass loss in the supernova and the velocity distribution of the asymmetric kick. For electron-capture supernovae (Podsiadlowski et al. 2004) the fraction of surviving triples is one-third. For higher kick velocities, we adopted that the exploding star was 2.7 M_☉ and as a consequence the fraction of survivors drops to 1/25 (Arzoumanian et al. 2002), 1/28 (Hobbs et al. 2005), and 1/50 (Cordes & Chernoff 1998), where the quoted literature refers to the adopted kick velocity distribution. Note here that the smaller amount of mass lost in the electron-capture supernova explosions helps considerably in preserving more triples, as opposed to the more violent kicks.

We conclude that for every three inner binaries that survive the electron-capture supernova the outer tertiary star remains in orbit around the inner binary, but that this fraction may drop considerably (to 1/50) when more mass is ejected in the supernova shell and the kick velocity is higher. Varying the amount of mass lost in the supernova explosion has a profound effect on the survivability of the triple, in particular since the Blaauw–Boersma kick imparted on the inner binary is proportional to this mass loss, and to the relative orbital velocity of the inner binary. It is interesting to note that the vast majority of the triples that survive the supernova explosion are dynamically stable, but their orbital eccentricities tend to be considerably higher than that indicated by the thermal distribution.

We synthesize the Galactic population of binaries like J1903+0327 by randomly selecting 10⁴ triples that survived the supernova explosion, and continue their evolution in AMUSE (see Section 4.1). Instead of performing a self-consistent evolution as adopted in Section 4.2 to validate the proposed scenario, we tentatively decoupled the mass-transfer process from the orbit integration. The population synthesis simulations start by resolving the mass transfer in the inner binary until the triple becomes dynamically unstable (Mardling & Aarseth 2001), after which we continue the simulation by resolving the dynamics of the three-body system until one star is ejected. During this latter part we ignore the mass transfer.

Note that this decoupled approach, though computationally cheaper by about a factor of 10³, is not a priori less reliable than the self-consistent simulations in Section 4.2, because our numerical methods provide no self-consistent way to resolve the mass-transfer process in eccentric orbits.

The majority (≳ 90%) of our simulations lead to the disruption of the triple, which results in the ejection of the initial primary (MSP), the secondary (partially stripped sub-giant or white dwarf), or the tertiary (main-sequence F, G, or K dwarf) star, leaving the other two stars in a binary. The respective ratios at which these occur are 0.013, 0.683, and 0.304 for circular outer orbits (which is unlikely after the supernova) and 0.03, 0.489, and 0.481 for highly eccentric (e ∼ 0.9) orbits. For extremely high eccentricities (e ∼ 0.99), these fractions become 0.05, 0.75, and 0.20, respectively. The outer orbits tend to have high eccentricities because they often barely survive the supernova kick, and as a consequence, the total fraction of systems in which the MSP remains in a binary with the original outer component is ≲ 0.5; the fraction of single MSPs is ∼0.05.

The two methods through which we resolved the triples (self-consistent, see Section 4.1, as opposed the here adopted decoupled approach) give slightly different results, but these can be related to the variation of the implementations. The largest effect is illustrated in Figure 2, where we demonstrate how the slow changes in the inner orbit drive the secular evolution of the outer orbit. The main consequences of this coupling are the higher mass of the NS by the time the triple becomes dynamically unstable and the small fraction of triples that survive the entire evolution because the adiabatic expansion of the outer orbit prevents the triple from becoming dynamically unstable. The consequential loss of systems, however, is well compensated by the increased probability that the donor in the LMXB is ejected when the triple becomes dynamically unstable. The differences between the two numerical implementations affect our estimates for the birthrate (see Section 5.4) on a ≲ 20% level. For clarity, the statistical uncertainties between the two different numerical approaches depend on Poissonian arguments rather than on the lack of our understanding of parts of the physical process, such as the common-envelop evolution, the process of mass transfer in non-circular orbits, and the nonlinear effects in the dynamically unstable configuration just before the triple is resolved. We control the latter by adopting a high order and numerically extremely precise algorithm (see Section 4.1). However, the tipping-point physics of ejecting the MSP, its close white dwarf companion or the main-sequence outer star remains elusive, even when integrating near machine precision.

The ratio of the birthrates of ordinary LMXBs to systems like J1903+0327 is a combination of the factors derived above: a few percent for the number of suitable triples, a reduction of a factor of ∼3 owing to the effect of an electron-capture supernova, and then a fraction ≲ 0.5 of systems in which the inner secondary is ejected. For larger average kick velocities (Cordes & Chernoff 1998; Arzoumanian et al. 2002), the fraction of triples that survive the supernova drops from 1/3 to 1/50 (see Section 5.2). The higher kicks cause an even stronger reduction in the number of single MSPs, because outer orbits tend to be highly eccentric, which leads to a reduction of the probability that the MSP is ejected once the triple becomes dynamically unstable. The birthrate for systems like J1903+0327 is then in the range of 2 × 10⁻⁴ to 3 × 10⁻³ times the birthrate of ordinary LMXBs. Estimates of the total number of LMXBs in the Galaxy are in the range 10⁴–10⁵ (Cote & Pylyser 1989). The typical lifetime of an LMXB is of the order of a Gyr, whereas the lifetime of systems like J1903+0327 is limited by the lifetime of the MSP and the difference between that of the inner secondary star and the outer tertiary star. Both stars have a mass of the order of 1 M_☉ and their lifetimes exceed several Gyr, which is a factor of a few longer than the lifetime of the LMXB phase. The number of J1903+0327-like systems in the Galaxy then is at least 30–300 for electron-capture supernovae but drops to 3–30 when we adopt considerably more mass to be lost in the supernova and higher velocity kicks (Cordes & Chernoff 1998; Arzoumanian et al. 2002).

In the triple scenario for the formation of J1903+0327 proposed by Champion et al. (2008), the observed ∼95 day orbital period is that of an unobserved massive 0.9–1.1 M_☉ white dwarf, while the observed G dwarf has a considerably wider orbit. The observed eccentricity of the inner orbit would in this case be driven by the secular evolution via the Kozai mechanism (Kozai 1962). This scenario has a number of serious shortcomings, one of which is the unusually high white dwarf mass, which leaves little room to produce a massive 1.67 M_☉ NS. Since its discovery, the observed change in eccentricity, ${\dot{e}} \sim 10^{-16} \,{\rm s}^{-1}$, is three orders of magnitude smaller than predictions based on the Kozai mechanism (Gopakumar et al. 2009), which excludes the presence of a highly inclined tertiary star. We exclude this scenario therefore from further consideration.

An alternative to the previous model would be the direct formation of a MSP in a supernova explosion, or via the fallback of material in a circum NS disk (Liu & Li 2009). A variety of arguments against the direct formation of rapidly spinning pulsar with a small surface magnetic field (2 × 10⁸ G) was provided by Champion et al. (2008). These include the following. First, of the 50 NSs in young supernova remnants, none are fast-spinning low-field pulsars (Kaspi & Helfand 2002). Second, a "born-fast" scenario for J1903+0327 would likely account for the 15 isolated MSPs detected in the galactic disk, but the spin distribution, space velocities, and energetics of these single MSPs are indistinguishable from those of recycled, not "born fast" binary MSPs (Archibald et al. 2009), while their space velocities and scale heights do not match those of non-recycled single pulsars (Lorimer et al. 2007). Third, magnetic fields in young pulsars likely originate either from dynamo action in the proto-NS (Thompson & Duncan 1993) or through compression of "frozen-in" fields of the progenitor star during collapse (Ferrario & Wickramasinghe 2006), and no young pulsars with B ≲ 10¹⁰ Gauss, like J1903+0327, are known. We add a fourth argument: before the direct collapse, the NS progenitor ascends the giant branch star and common-envelope evolution with the 95 day-orbit G dwarf would have dramatically reduced the orbital period, in the same way as normal LMXBs form. For these four reasons, a core-collapse "born-fast" MSP can be ruled out.

The above arguments hold equally well against the accretion-induced collapse of a massive and rapidly rotating white dwarf to an NS. While this may lead to millisecond rotation periods, the NS formed here will be just as hot and differentially rotating during the early liquid phases as an NS formed by core collapse. In both cases, 10⁵³ erg in gravitational energy is released, which will erase all memory of the violent mechanism in which the NS was formed. Therefore, dynamo action in both cases will not depend on the formation mechanism, and there is no fundamental reason for expecting a weak magnetic field in NSs formed by accretion-induced collapse. Furthermore, the direct accretion-induced collapse of a white dwarf to an MSP in population studies produces binaries with a orbital period ≲ 20 days, which is considerably shorter than observed in J1903+0327 (Chen et al. 2011). In addition, this scenario requires the white dwarf to accrete from an evolved companion star, which is inconsistent with the observed main-sequence companion G star. This model thus fails to explain J1903+0327.

Finally, a formation scenario where J1903+0327 is formed when the donor star in the inner binary is ablated and destroyed, such as the case for the "black-widow" system PSR B1957+20 (Fruchter et al. 1988), has thee problems. First, the observed timescales for straightforward evaporation of the donor star are too long (Champion et al. 2008). Second, formation of such a system likely involves an exchange interaction (King et al. 2003), which would be greatly impeded by the outer G dwarf companion. Third, even if no exchange took place, the slow evaporation of the donor in the inner binary requires the triple to be dynamically stable with respect to the G star in its current orbit. Since the supernova explosion can at most account for a reduction of a factor of two in the orbital separation, the common envelope should in that case be responsible for a further reduction from the initial orbital separation of ≳ 560 R_☉ to the currently observed ∼123 R_☉, which requires considerable fine tuning. We therefore conclude that J1903+0327 is unlikely to have originated through a black-widow-like scenario.

Each of the above models has serious shortcomings, and we conclude that none of the scenarios discussed above give a satisfactory explanation for the formation of J1903+0327.

There are currently no known triples with parameters suitable for evolving into systems like J1903+0327. However, there is observational evidence that such triples exist, as it is proposed that 4U 2129+47 (V1727 Cyg), a 5.24 hr LMXB, is accompanied by a spectral type F dwarf in an eccentric orbit of about 175 days (Garcia et al. 1989; Bothwell et al. 2008; Lin et al. 2009). This triple could have formed in the same way as J1903+0327, with the exception that after the common envelope and the subsequent supernova explosion, the inner binary period was smaller than the bifurcation period of about one day (Pylyser & Savonije 1988, 1989). The consequence of such a short orbital period is that the LMXB evolves to an even shorter orbital period. Interestingly the 175 day orbit of the outer star is too small to have been dynamically stable at the birth of the triple (Mardling & Aarseth 2001). We argue that the common envelope and/or the supernova may have reduced the orbital separation of the outer star.

We validated the probability that a binary orbit shrinks as a result of the supernova by means of population synthesis and conclude that in a fraction of 0.511 of the triples that survive the supernova the separation of the circularized outer orbit is smaller than the initial orbit. However, if the progenitor of 4U 2129+47 had parameters comparable to what we derived for J1903+0327 in Section 2, the common-envelope phase of the inner binary must have resulted in a reduction of the outer orbit as well, by ≳ 210 R_☉. Such a reduction in the separation of the outer orbit by the common envelope enormously boosts the survivability of the triple in the supernova explosion (see Figure 4), and therefore dramatically increases the birthrate of binaries like J1903+0327.

Following the same scenario as for forming 4U 2129+47 but with a less massive (≲ 8 M_☉) initial inner primary star, the inner binary could evolve into a cataclysmic variable such as EC 19314−5915, with an orbital period of 4.75 hr. The observed radial velocity of ∼9 km s⁻¹ has been attributed to a G8 dwarf in orbit around the cataclysmic variable (Buckley et al. 1992), which is consistent with a semimajor axis of ∼2400 R_☉.

The real missing link would be the discovery of a relatively wide (≳ 30 R_☉) LMXB that is orbited by a tertiary low-mass main-sequence star. The mass-transfer phase in our simulations averaged about 23 Myr (see Section 4.2), while the lifetime of the MSP in J1903+0327 is at least 1 Gyr. We therefore expect that the Galaxy contains at most seven such wide triple LMXBs.

To study the effect of the adopted mass-transfer rate in the LMXB phase, we performed several simulations using the time-variable mass-transfer rates from detailed simulations of LMXBs (Tauris & Savonije 1999). The consequence is that shortly after the onset of Roche lobe contact the rate of mass transfer exceeds our adopted constant rate of $\langle \dot{m} \rangle = 6 \times 10^{-10} (\mbox{${P}_{\rm orb}$}({\rm initial})/1\, {\rm day})$ M_☉yr⁻¹ to drop below that rate at later time (see Section 4.1). This more realistic mass-transfer rate causes the inner orbits to expand more dramatically shortly after the start of the mass-transfer process and results on average in slightly lower NS masses and somewhat wider orbits. Though we did not further study the consequences of this more realistic approach, we realize that the more dramatic widening of the inner orbit in the early LMXB phase causes an increase of the birthrate of systems like J1903+0327.