Black widow evolution: magnetic braking by an ablated wind
aa r X i v : . [ a s t r o - ph . H E ] J a n MNRAS , 1–11 (2015) Preprint 15 January 2020 Compiled using MNRAS L A TEX style file v3.0
Black widow evolution:magnetic braking by an ablated wind
Sivan Ginzburg ⋆ † and Eliot Quataert Department of Astronomy and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Black widows are close binary systems in which a millisecond pulsar is orbited by acompanion a few per cent the mass of the sun. It has been suggested that the pulsar’srotationally powered γ -ray luminosity gradually evaporates the companion, eventuallyleaving behind an isolated millisecond pulsar. The evaporation efficiency is determinedby the temperature T ch ∝ F / to which the outflow is heated by the flux F on adynamical time-scale. Evaporation is most efficient for companions that fill their Rochelobes. In this case, the outflow is dominated by a cap around the L1 point with an angle θ g ∼ ( T ch / T g ) / , and the evaporation time is t evap = . ( T ch / T g ) − Gyr, where T g > T ch is the companion’s virial temperature. We apply our model to the observed blackwidow population, which has increased substantially over the last decade, consideringeach system’s orbital period, companion mass, and pulsar spin-down power. Whilethe original (Fruchter et al. 1988) black widow evaporates its companion on a few Gyrtime-scale, direct evaporation on its own is too weak to explain the overall population.We propose instead that the evaporative wind couples to the companion’s magneticfield, removes angular momentum from the binary, and maintains stable Roche-lobeoverflow. While a stronger wind carries more mass, it also reduces the Alfv´en radius,making this indirect magnetic braking mechanism less dependent on the flux t mag ∝ t / . This reduces the scatter in evolution times of observed systems, thus betterexplaining the combined black widow and isolated millisecond pulsar populations. Key words: binaries: close – pulsars: general
Black widows are systems with a millisecond pulsar and alow-mass companion, a few per cent the mass of the sun, ona short orbit of several hours. The first such system was de-tected by Fruchter et al. (1988), and a few dozen similar oneshave been found since. Typically, the companion’s orbital pe-riod P orb and its minimum mass m sin i are determined fromthe periodic timing variations in the pulsar’s radio pulse. Upto the last decade, only three black widows were known inthe galactic field, with the rest found in globular clusters.Radio followups of Fermi γ -ray sources, as well as the des-ignated High Time Resolution Universe survey, have sincegreatly increased the number of known field millisecond pul-sars, a significant fraction of which turned out to be blackwidows (Ray et al. 2012; Keith 2013; Roberts 2013). Amongthese discoveries are the extremely low mass companionsto PSR J1719-1438 (Bailes et al. 2011) and to PSR J2322-2650 (Spiewak et al. 2018), with minimum masses similar to ⋆ E-mail: [email protected] †
51 Pegasi b Fellow.
Jupiter—about 20 times lighter than the companion orbitingthe original black widow pulsar.According to the prevailing theory (see Manchester2017, for a review and references), millisecond pulsars werespun up to their fast rotation rates by accreting materialfrom a main-sequence binary companion. In this picture, iso-lated millisecond pulsars are the end products of binary evo-lution, while black widows represent the missing link, withcompanions reduced to a fraction of their original mass.What is the mechanism that sustains mass loss fromblack widow companions, eventually leading to their fulldestruction, leaving behind isolated millisecond pulsars?Gravitational radiation can shrink the binary orbit anddrive it towards Roche-lobe overflow. However, by the timethe companion is reduced to a few percent of the so-lar mass, gravitational waves are too weak (e.g. fig. 1in Romani et al. 2016), leading several authors to con-sider, instead, the evaporation (ablation) of the compan-ion by high-energy photons or particles powered by the pul-sar’s spin-down energy (Kluzniak et al. 1988; Phinney et al.1988; Ruderman et al. 1989a; cf. Eichler & Levinson 1988;Levinson & Eichler 1991). © S. Ginzburg and E. Quataert
The growing number of field black widows discov-ered in the last decade has motivated a renewed theoreti-cal interest in these systems (Benvenuto et al. 2012, 2014,2015; Chen et al. 2013; Jia & Li 2015, 2016; Liu & Li 2017;Ablimit 2019). Much of the theoretical effort has focused onreproducing the observed black widow population by evolv-ing with time binary systems that initially host a main se-quence companion. Despite being a key ingredient in suchevolutionary tracks, the evaporation of the companion by thepulsar’s irradiation is usually parametrized with a simple lin-ear relation between the pulsar’s spin-down power and thecompanion’s mass loss rate, with the coefficient unknown.Here, we calculate the evaporation efficiency by mod-elling the hydrodynamical wind launched off the compan-ion’s surface by the incident pulsar radiation. We find thatthe mass ejection rate does not scale linearly with the pul-sar’s luminosity, but shows a more complex dependence,and is typically lower than previously assumed. We applyour model to the observed black widow population and findthat, for the majority of systems, evaporation by the pul-sar’s radiation on its own is too weak to play a major rolein the companion’s evolution. As an alternative, we suggestthat the evaporative wind may couple to the companion’smagnetic field and remove angular momentum from the sys-tem, thereby maintaining the companion in stable Roche-lobe overflow. Such a mechanism can greatly amplify themass loss rate and explain the observations.The remainder of this paper is organized as follows. InSection 2 we derive general expressions for the mass loss ratein evaporative winds. In Section 3 we apply these expressionsto the observed black widow population and calculate evap-oration efficiencies and time-scales. In Section 4 we considermagnetic braking by interaction of the wind with the com-panion’s magnetic field and reanalyse the observations. Wesummarize and discuss our results in Section 5
In this section we derive the rate of mass ablation off thecompanion’s surface by the pulsar’s irradiation. We omitorder-unity coefficients and focus on the scaling relations.In Section 2.1 we discuss the pressure at which the wind islaunched p . In Sections 2.2 and 2.3 we calculate the outflowrate Û m as a function of p and the heating rate per particle Γ . We discuss Γ in Section 2.4. The equilibrium of optically thin gas subject to ionizing ra-diation has been calculated extensively in the past (see ref-erences below). Balancing cooling with heating and ioniza-tion with recombination yields an equilibrium gas pressure p ( ρ, F ) , where ρ is the density and F is the radiative flux.Heating and ionization processes scale as ∝ ρ F whereas cool-ing and recombination scale as ∝ ρ ; as a result, the equilib-rium ionization state and gas temperature T are functionsof ρ / F for a given irradiation spectrum. Using the ideal gaslaw p = ρ kT / µ , with k denoting Boltzmann’s constant and µ the molecular weight, wee see that p / F is a function of ρ / F .A typical equilibrium curve is given in fig. 2 ofMcCray & Hatchett (1975, see also Basko & Sunyaev 1973; Krolik et al. 1981). At high densities, collisionally excitedHydrogen line cooling, which depends exponentially on thetemperature, balances photoionization heating, and main-tains an almost constant T ∼ K, such that p ∝ ρ .At lower densities, Hydrogen becomes ionized and equi-librium between heating (Compton and photoionization ofheavy elements) and Bremsstrahlung cooling is achievedby raising T . For an approximately constant heating rateper particle, Bremsstrahlung dictates p ∝ ρ T ∝ ρ − (e.g.Rybicki & Lightman 1979), such that p ( ρ ) reaches a localminimum value p . The temperature and pressure keeprising with decreasing ρ until a maximum temperature kT IC = ε / is reached, when inverse Compton scatteringbalances Compton heating ( ε is the typical photon energy;see Rybicki & Lightman 1979). For a more detailed explana-tion of the various cooling and heating mechanisms and theshape of the equilibrium curve, see Buff & McCray (1974),London et al. (1981), and chapter 10 in Krolik (1999).The pressure and density decrease monotonically whengoing out in the companion’s atmosphere. Once the pres-sure drops below the local minimum p , the temperaturemust jump from T ∼ K to T = T IC in order to reachequilibrium. However, for a hard enough spectrum, the as-sociated sound speed c IC = ( kT IC / µ ) / is larger than theescape velocity from the companion v g , launching a hydro-dynamic wind (Basko & Sunyaev 1973; McCray & Hatchett1975; Basko et al. 1977).We follow Begelman et al. (1983) and parametrize thewind launching pressure p with the dimensionless Ξ ′ ≡ Fp c , (1)where c is the speed of light ( F / c is the radiation pressure,which is dynamically unimportant because F is well belowthe companion’s Eddington limit). For a variety of radia-tion spectra, Ξ ′ ∼ (table 2 in London et al. 1981; see alsoKrolik et al. 1981 and Krolik 1999 but notice the differencein definition between Ξ ′ and their ionization parameter Ξ ),and for the remainder of the paper we assume Ξ ′ = . We mark the escaping mass flux with f . Conservation ofmomentum sets an upper limit of f p / v g , since escap-ing gas must accelerate to at least the escape velocity v g by the pressure at the base of the wind p . This is also thelimit derived by Basko et al. (1977) by modelling the out-flow as an isothermal Parker (1958) wind. Ruderman et al.(1989a,b) adopt a similar value for their nominal mass-lossrate. However, as we show below (Fig. 1), black widow sys-tems generally do not quite reach this maximum.We model the outflow as a non-isothermal wind andwrite the equations of mass and momentum conservation inspherical symmetry: ρ v r = const. (2) ρ v d v d r = − d p d r − Gm ρ r , (3)where m is the companion’s mass, r is the distance from itscentre, v is the outflow velocity, and G is the gravitational MNRAS , 1–11 (2015) lack widow evolution constant. We combine equations (2) and (3) with the idealgas law:d v d r − c s v ! = c s r − Gmr − k µ d T d r , (4)with c s ≡ ( kT / µ ) / denoting the isothermal sound speed.The initial conditions for equation (4) are v → , T → ,and p = p at the base of the wind, i.e. on the companion’ssurface r = R . As usual (e.g. Parker 1965), we are seeking fora solution that passes through a sonic point r s , at which v = c s , and accelerates to supersonic velocities. From equation(4), the sonic point is given by (London & Flannery 1982) r s = Gm µ kT + r s d ln T d ln r . (5)The first term on the right hand side of equation (5)gives the isothermal-wind solution. The temperature profiled T / d r in the second term is determined by a competitionbetween the heating and dynamical time-scales of the flow.Begelman et al. (1983) characterised this competition withthe temperature T ch , defined by kT ch = Γ Rc ch , (6)where c ch ≡ ( kT ch / µ ) / , and Γ is the heating rate per parti-cle. Intuitively, this is the temperature reached by the flowon a dynamical time. We solve equation (6): kT ch = µ / ( Γ R ) / ∝ F / . (7)Since the heating rate is proportional to the incident flux(Section 2.4), T ch ∝ F / is a measure of the pulsar’s irradi-ation intensity. Begelman et al. (1983) identified three windregimes, which are distinguished by the characteristic tem-perature T ch and its relation to the Compton temperature T IC and to the virial temperature T g ≡ Gm µ /( kR ) ≪ T IC (no wind is launched if T IC < T g ). The analytical solutionof the equations of motion (including heating) is given byBegelman et al. (1983), while here we provide an intuitivesketch. A summary of the results is provided in Table 1. T g < T IC < T ch In this regime, the heating time kT / Γ is shorter than thedynamical time R /( kT / µ ) / even for the maximum temper-ature T IC . The temperature of the flow therefore rises from alow value at the companion’s surface, reaches T IC and levelsoff at a distance ∆ r ≡ r − R ≪ R . As long as the temperaturerises, d ln T / d ln r ≫ , and there is no solution to equation(5). The sonic point is reached when the asymptotic temper-ature T IC is approached and d ln T / d ln r ≈ , at a distance ∆ r s ≪ R ; the gravitational term in equation (5) is negligible Gm µ /( kT IC ) ≪ R .At the sonic point, the mass flux is given by f = ρ c s = p / c s . By integrating equation (3) we find that p + ρ v isconserved as long as the geometry is planar ( ∆ r ≪ R ) and f = ρ v is constant; the gravitational term in the integrationis negligible by a factor of ( T g / T IC )( ∆ r / R ) ≪ . The pressureat the sonic point is therefore given by p + ρ c s = p , i.e. p = p / . We conclude that the mass flux is this regime is f = p /( c IC ) . T g < T ch < T IC We make the ansatz that in this regime the sonic point islocated at a distance ∆ r s ∼ R from the surface. From thedefinition in equation (6), the temperature there must thenbe ∼ T ch . There is no solution to equation (5) as long as ∆ r ≪ R because d ln T / d ln r ≫ . When ∆ r becomes comparableto R , d ln T / d ln r ≈ and the sonic point is reached; thegravitational term in equation (5) is negligible in this regime Gm µ /( kT ch ) ≪ R . This justifies our assumption that ∆ r s ∼ R .The geometry is roughly planar for ∆ r . R , so p + ρ v is conserved, up to an order-unity correction, by integrationof equation (3) similarly to Section 2.2.1; the gravitationalterm in the integration is negligible by a factor of T g / T ch .The pressure at the sonic point therefore satisfies p + ρ c s ∼ p , i.e. p ∼ p . The mass flux through the sonic point is f = p / c s ∼ p / c ch . T ch < T g < T IC In this regime, our ansatz that the sonic point is at ∆ r s ∼ R —and as a consequence at T ∼ T ch (see Section 2.2.2)—fails.The gravitational term in equation (5) cannot be neglectedin this case because Gm µ /( kT ch ) ≫ R , indicating that thesonic point is reached only at a radius r s ≫ R . Therefore,at ∆ r ∼ R the flow is subsonic ( v ≪ c s ), d ln T / d ln r ∼ asbefore, and the right-hand side of equation (4) is negative: c s ( − d ln T / d ln r ) − v g < . This expression must changesign at the sonic point, and since v g decreases with r , c s . v g at ∆ r ∼ R (similar to an isothermal Parker wind). Accord-ing to equation (4), solutions with c s ≪ v g at ∆ r ∼ R arecharacterised by an exponentially fast acceleration on a scalemuch smaller than R : d ln v / d ln r ∼ v g / c s ≫ . While theseare valid solutions for an isothermal Parker wind, an expo-nentially subsonic velocity at ∆ r ∼ R would give the heatedflow in our case enough time R / v to reach temperatures wellabove T g (i.e. c s ≫ v g ) at ∆ r ∼ R . We conclude that the onlyself consistent solution is the one with c s ∼ v g at ∆ r ∼ R .We find the pressure at ∆ r ∼ R by integrating equation(3) from the surface. The hydrodynamical term is negligibledue to the subsonic velocities, and since for c s ∼ v g the scaleheight is of order R , we find that p ∼ p .The mass flux at ∆ r ∼ R is given by f = ρ v = ( p / c s )( v / c s ) . The Mach number v / c s at ∆ r ∼ R is deter-mined by the condition that the subsonic flow is heatedto T g on a time-scale R / v : kT g = Γ R / v . By comparing thiscondition to equation (6), we find that v / c ch = T ch / T g and v / c s ∼ ( v / c ch )( c ch / v g ) = ( T ch / T g ) / . Finally, the flux is givenby f ∼ ( p / v g )( v / c s ) = ( p / c ch )( T ch / T g ) . As also pointed outby Begelman et al. (1983), the mass flux in the cold regime isnot determined by conditions at the sonic point, but ratherby the Mach number constraint at ∆ r ∼ R .The major difference between Begelman et al. (1983) andour scenario is the heating function. Γ drops as ( r / R ) − forthe accretion-disc winds of Begelman et al. (1983), whereasin our case Γ is constant until r is comparable to the compan-ion’s distance from the pulsar a ≫ R . As discussed above,however, the mass flux in all three wind regimes is deter-mined by conditions at r − R . R , where Γ is approxi-mately constant in both cases. The Begelman et al. (1983) MNRAS000
In this section we derive the rate of mass ablation off thecompanion’s surface by the pulsar’s irradiation. We omitorder-unity coefficients and focus on the scaling relations.In Section 2.1 we discuss the pressure at which the wind islaunched p . In Sections 2.2 and 2.3 we calculate the outflowrate Û m as a function of p and the heating rate per particle Γ . We discuss Γ in Section 2.4. The equilibrium of optically thin gas subject to ionizing ra-diation has been calculated extensively in the past (see ref-erences below). Balancing cooling with heating and ioniza-tion with recombination yields an equilibrium gas pressure p ( ρ, F ) , where ρ is the density and F is the radiative flux.Heating and ionization processes scale as ∝ ρ F whereas cool-ing and recombination scale as ∝ ρ ; as a result, the equilib-rium ionization state and gas temperature T are functionsof ρ / F for a given irradiation spectrum. Using the ideal gaslaw p = ρ kT / µ , with k denoting Boltzmann’s constant and µ the molecular weight, wee see that p / F is a function of ρ / F .A typical equilibrium curve is given in fig. 2 ofMcCray & Hatchett (1975, see also Basko & Sunyaev 1973; Krolik et al. 1981). At high densities, collisionally excitedHydrogen line cooling, which depends exponentially on thetemperature, balances photoionization heating, and main-tains an almost constant T ∼ K, such that p ∝ ρ .At lower densities, Hydrogen becomes ionized and equi-librium between heating (Compton and photoionization ofheavy elements) and Bremsstrahlung cooling is achievedby raising T . For an approximately constant heating rateper particle, Bremsstrahlung dictates p ∝ ρ T ∝ ρ − (e.g.Rybicki & Lightman 1979), such that p ( ρ ) reaches a localminimum value p . The temperature and pressure keeprising with decreasing ρ until a maximum temperature kT IC = ε / is reached, when inverse Compton scatteringbalances Compton heating ( ε is the typical photon energy;see Rybicki & Lightman 1979). For a more detailed explana-tion of the various cooling and heating mechanisms and theshape of the equilibrium curve, see Buff & McCray (1974),London et al. (1981), and chapter 10 in Krolik (1999).The pressure and density decrease monotonically whengoing out in the companion’s atmosphere. Once the pres-sure drops below the local minimum p , the temperaturemust jump from T ∼ K to T = T IC in order to reachequilibrium. However, for a hard enough spectrum, the as-sociated sound speed c IC = ( kT IC / µ ) / is larger than theescape velocity from the companion v g , launching a hydro-dynamic wind (Basko & Sunyaev 1973; McCray & Hatchett1975; Basko et al. 1977).We follow Begelman et al. (1983) and parametrize thewind launching pressure p with the dimensionless Ξ ′ ≡ Fp c , (1)where c is the speed of light ( F / c is the radiation pressure,which is dynamically unimportant because F is well belowthe companion’s Eddington limit). For a variety of radia-tion spectra, Ξ ′ ∼ (table 2 in London et al. 1981; see alsoKrolik et al. 1981 and Krolik 1999 but notice the differencein definition between Ξ ′ and their ionization parameter Ξ ),and for the remainder of the paper we assume Ξ ′ = . We mark the escaping mass flux with f . Conservation ofmomentum sets an upper limit of f p / v g , since escap-ing gas must accelerate to at least the escape velocity v g by the pressure at the base of the wind p . This is also thelimit derived by Basko et al. (1977) by modelling the out-flow as an isothermal Parker (1958) wind. Ruderman et al.(1989a,b) adopt a similar value for their nominal mass-lossrate. However, as we show below (Fig. 1), black widow sys-tems generally do not quite reach this maximum.We model the outflow as a non-isothermal wind andwrite the equations of mass and momentum conservation inspherical symmetry: ρ v r = const. (2) ρ v d v d r = − d p d r − Gm ρ r , (3)where m is the companion’s mass, r is the distance from itscentre, v is the outflow velocity, and G is the gravitational MNRAS , 1–11 (2015) lack widow evolution constant. We combine equations (2) and (3) with the idealgas law:d v d r − c s v ! = c s r − Gmr − k µ d T d r , (4)with c s ≡ ( kT / µ ) / denoting the isothermal sound speed.The initial conditions for equation (4) are v → , T → ,and p = p at the base of the wind, i.e. on the companion’ssurface r = R . As usual (e.g. Parker 1965), we are seeking fora solution that passes through a sonic point r s , at which v = c s , and accelerates to supersonic velocities. From equation(4), the sonic point is given by (London & Flannery 1982) r s = Gm µ kT + r s d ln T d ln r . (5)The first term on the right hand side of equation (5)gives the isothermal-wind solution. The temperature profiled T / d r in the second term is determined by a competitionbetween the heating and dynamical time-scales of the flow.Begelman et al. (1983) characterised this competition withthe temperature T ch , defined by kT ch = Γ Rc ch , (6)where c ch ≡ ( kT ch / µ ) / , and Γ is the heating rate per parti-cle. Intuitively, this is the temperature reached by the flowon a dynamical time. We solve equation (6): kT ch = µ / ( Γ R ) / ∝ F / . (7)Since the heating rate is proportional to the incident flux(Section 2.4), T ch ∝ F / is a measure of the pulsar’s irradi-ation intensity. Begelman et al. (1983) identified three windregimes, which are distinguished by the characteristic tem-perature T ch and its relation to the Compton temperature T IC and to the virial temperature T g ≡ Gm µ /( kR ) ≪ T IC (no wind is launched if T IC < T g ). The analytical solutionof the equations of motion (including heating) is given byBegelman et al. (1983), while here we provide an intuitivesketch. A summary of the results is provided in Table 1. T g < T IC < T ch In this regime, the heating time kT / Γ is shorter than thedynamical time R /( kT / µ ) / even for the maximum temper-ature T IC . The temperature of the flow therefore rises from alow value at the companion’s surface, reaches T IC and levelsoff at a distance ∆ r ≡ r − R ≪ R . As long as the temperaturerises, d ln T / d ln r ≫ , and there is no solution to equation(5). The sonic point is reached when the asymptotic temper-ature T IC is approached and d ln T / d ln r ≈ , at a distance ∆ r s ≪ R ; the gravitational term in equation (5) is negligible Gm µ /( kT IC ) ≪ R .At the sonic point, the mass flux is given by f = ρ c s = p / c s . By integrating equation (3) we find that p + ρ v isconserved as long as the geometry is planar ( ∆ r ≪ R ) and f = ρ v is constant; the gravitational term in the integrationis negligible by a factor of ( T g / T IC )( ∆ r / R ) ≪ . The pressureat the sonic point is therefore given by p + ρ c s = p , i.e. p = p / . We conclude that the mass flux is this regime is f = p /( c IC ) . T g < T ch < T IC We make the ansatz that in this regime the sonic point islocated at a distance ∆ r s ∼ R from the surface. From thedefinition in equation (6), the temperature there must thenbe ∼ T ch . There is no solution to equation (5) as long as ∆ r ≪ R because d ln T / d ln r ≫ . When ∆ r becomes comparableto R , d ln T / d ln r ≈ and the sonic point is reached; thegravitational term in equation (5) is negligible in this regime Gm µ /( kT ch ) ≪ R . This justifies our assumption that ∆ r s ∼ R .The geometry is roughly planar for ∆ r . R , so p + ρ v is conserved, up to an order-unity correction, by integrationof equation (3) similarly to Section 2.2.1; the gravitationalterm in the integration is negligible by a factor of T g / T ch .The pressure at the sonic point therefore satisfies p + ρ c s ∼ p , i.e. p ∼ p . The mass flux through the sonic point is f = p / c s ∼ p / c ch . T ch < T g < T IC In this regime, our ansatz that the sonic point is at ∆ r s ∼ R —and as a consequence at T ∼ T ch (see Section 2.2.2)—fails.The gravitational term in equation (5) cannot be neglectedin this case because Gm µ /( kT ch ) ≫ R , indicating that thesonic point is reached only at a radius r s ≫ R . Therefore,at ∆ r ∼ R the flow is subsonic ( v ≪ c s ), d ln T / d ln r ∼ asbefore, and the right-hand side of equation (4) is negative: c s ( − d ln T / d ln r ) − v g < . This expression must changesign at the sonic point, and since v g decreases with r , c s . v g at ∆ r ∼ R (similar to an isothermal Parker wind). Accord-ing to equation (4), solutions with c s ≪ v g at ∆ r ∼ R arecharacterised by an exponentially fast acceleration on a scalemuch smaller than R : d ln v / d ln r ∼ v g / c s ≫ . While theseare valid solutions for an isothermal Parker wind, an expo-nentially subsonic velocity at ∆ r ∼ R would give the heatedflow in our case enough time R / v to reach temperatures wellabove T g (i.e. c s ≫ v g ) at ∆ r ∼ R . We conclude that the onlyself consistent solution is the one with c s ∼ v g at ∆ r ∼ R .We find the pressure at ∆ r ∼ R by integrating equation(3) from the surface. The hydrodynamical term is negligibledue to the subsonic velocities, and since for c s ∼ v g the scaleheight is of order R , we find that p ∼ p .The mass flux at ∆ r ∼ R is given by f = ρ v = ( p / c s )( v / c s ) . The Mach number v / c s at ∆ r ∼ R is deter-mined by the condition that the subsonic flow is heatedto T g on a time-scale R / v : kT g = Γ R / v . By comparing thiscondition to equation (6), we find that v / c ch = T ch / T g and v / c s ∼ ( v / c ch )( c ch / v g ) = ( T ch / T g ) / . Finally, the flux is givenby f ∼ ( p / v g )( v / c s ) = ( p / c ch )( T ch / T g ) . As also pointed outby Begelman et al. (1983), the mass flux in the cold regime isnot determined by conditions at the sonic point, but ratherby the Mach number constraint at ∆ r ∼ R .The major difference between Begelman et al. (1983) andour scenario is the heating function. Γ drops as ( r / R ) − forthe accretion-disc winds of Begelman et al. (1983), whereasin our case Γ is constant until r is comparable to the compan-ion’s distance from the pulsar a ≫ R . As discussed above,however, the mass flux in all three wind regimes is deter-mined by conditions at r − R . R , where Γ is approxi-mately constant in both cases. The Begelman et al. (1983) MNRAS000 , 1–11 (2015)
S. Ginzburg and E. Quataert
Table 1.
Wind regimes identified by Begelman et al. (1983) fordifferent values of the characteristic temperature T ch which is de-fined in equation (7); the radius of the sonic point is given by r s = R + ∆ r s . Hot Intermediate ColdTemperature T ch > T IC T g < T ch < T IC < T g Sonic point ∆ r s ≪ R ∼ R ≫ R Mass flux f p / c IC p / c ch ( T ch / T g ) p / c ch Figure 1.
Schematic plot of the companion’s mass-loss rate Û m asa function of the characteristic temperature T ch , which is definedin equation (7). T ch is a measure of the flux received from thepulsar F , and its value compared to the virial ( T g ) and Comp-ton ( T IC ) temperatures determines the outflow regime (Table 1).The solid red and orange lines depict the hot and intermediateregimes, respectively. Roche-lobe overflowing (RLO) companionsdiffer from detached systems in the cold regime (dot–dashed bluelines; see Section 2.3). The dashed black line marks the maxi-mum rate assumed by Basko et al. (1977) and Ruderman et al.(1989a,b). This limit, which is linear in the incident luminos-ity ( p ∝ F ; see Section 2.1), is reached only for a single value T ch = T g . disc-wind results are therefore applicable to irradiated com-panions as well. In the cold-wind regime (see Table 1), the mass flux is af-fected by gravity, embodied by the virial temperature T g .We consider two cases with different gravity fields: detachedsystems, in which the companion is far from filling its Rochelobe, and systems in Roche-lobe overflow.In detached systems, the gravitational field is spheri-cally symmetric and equal to g = Gm / r . The cold wind mass-loss rate in this case is simply | Û m | (detached) ∼ π R p c ch (cid:18) T ch T g (cid:19) = π R p v g (cid:18) T ch T g (cid:19) / , (8)where we have substituted the mass flux from Table 1.The gravitational field of companions that fill theirRoche lobe is weakened by the pulsar’s tidal force, especiallynear the L1 and L2 Lagrange points. The effective gravitycomponent normal to the companion’s surface varies approx-imately as g eff ≈ g ( − cos θ ) , where θ is measured relativeto the line connecting the pulsar’s and companion’s centres.Close to the L1 point (the L2 point is not irradiated), theeffective virial temperature is therefore T eff g ≈ T g θ . Inside acritical angle of θ g = (cid:18) T ch T g (cid:19) / (9)the effective gravity is too weak to restrain the flow T eff g < T ch , and the mass flux is given by the intermediate-windsolution f = p / c ch (Table 1). The surface area of the θ . θ g region is ∼ θ g of the total sphere, so the overall mass-lossrate is given by | Û m | (overflowing) ∼ π R θ g p c ch = π R p v g (cid:18) T ch T g (cid:19) / . (10)The cap around L1 dominates the mass flow: while its areais only T ch / T g of the total sphere, the flux there is higher bya factor of ( T g / T ch ) .Fig. 1 shows the three outflow regimes discussed in Ta-ble 1. The hot and intermediate regimes are not sensitive tothe gravity field and the outflow rate in these cases is simply | Û m | ∼ π R f . The cold regime differs for overflowing and de-tached systems, and is given by equations (8) and (10). Fig.1 demonstrates that the commonly assumed relation Û m ∝ F is usually inaccurate. In fact, the limit | Û m | . π R p / v g (Basko et al. 1977; Ruderman et al. 1989a,b) is reached onlyfor systems in which T ch ≈ T g . In order to calculate T ch and determine the relevant regimein Fig. 1, we need to estimate the heating rate per particle Γ , which appears in equation (7). Following the discussionin Section 2.2, we are interested in Γ for gas temperaturesin the range T g . T < T IC (even in the cold regime, thetemperature reaches T g at ∆ r ∼ R ; see Section 2.2.3). Thevirial temperature of black widow companions is typically T g & K, for which Compton scattering is the dominantheating mechanism (e.g. Buff & McCray 1974).The Compton heating rate is given by Γ ∼ σ T F · ( x x ≪ x − ln x x ≫ , (11)where σ T is the Thomson cross section and x ≡ ε /( m e c ) ,with ε denoting the photon energy and m e the electron mass.When x ≪ , each interacting photon deposits a fraction ∼ x of its energy. When x ≫ , photons deposit almostall of their energy, but the cross section is reduced accord-ing to the Klein-Nishina formula (e.g. Rybicki & Lightman1979). Equation (11) implies that the heating rate of a MNRAS , 1–11 (2015) lack widow evolution roughly flat Crab-like γ -ray spectrum (B¨uhler & Blandford2014) is dominated by MeV photons ( ε ∼ m e c ). Further-more, measurements by Compton (Kuiper et al. 2000) andinterpolation of
NuSTAR (Gotthelf & Bogdanov 2017) with
Fermi (Abdo et al. 2013) data suggest that millisecond pul-sar emission might peak around MeV. For these reasons, weestimate that Γ ∼ σ T F .While pair production has a larger cross section thanCompton scattering at the highest photon energies, the re-sulting e ± pairs have to deposit their energy high enoughin the companion’s atmosphere to drive an outflow. A sim-ilar challenge is faced by the TeV e ± wind that poten-tially carries a large fraction of the pulsar’s spin-down power(Ruderman et al. 1989a). According to one scenario, a ∼ G magnetic field may convert the energy of TeV particlesinto MeV photons by synchrotron radiation (Kluzniak et al.1988; Phinney et al. 1988; Ruderman et al. 1989a). Suchsecondary photons can be accounted for by adjusting F ap-propriately. In Section 3 we consider only the direct photonswhen calculating the ablating flux, possibly underestimating F by a factor of a few. Following previous studies (Stevens et al. 1992;Benvenuto et al. 2012, 2014; Chen et al. 2013; Jia & Li2015, 2016; Liu & Li 2017), we define the mass-lossefficiency η as Gm Û mR ≡ − η L γ (cid:18) Ra (cid:19) , (12)where L γ = π a F is the pulsar’s γ -ray luminosity and a isits separation from the companion. η measures the fractionof the incident radiation energy that is invested in overcom-ing the companion’s gravity. As discussed in Section 2.4, weonly consider the γ -ray luminosity that directly couples withthe companion’s upper atmosphere to drive a wind, and notthe total spin-down luminosity. We adopt the same beamingfactor (i.e. no beaming) as Abdo et al. (2013) for consistencywith their L γ measurements.From here onward, we assume that companions arealways close to filling their Roche lobes (Benvenuto et al.2012, 2015; Draghis et al. 2019), so R / a ≃ . ( m / M ) / (Eggleton 1983) where M is the pulsar’s mass. Appar-ently, black widow pulsars occupy the high end of theobserved neutron star mass distribution, with M & M ⊙ (van Kerkwijk et al. 2011; Romani et al. 2015; Linares 2019, M ⊙ denotes the solar mass), presumably due to precedingmass transfer from their companions. Here, we adopt in-stead the canonical M = . M ⊙ , for consistency with pre-vious studies and existing data bases; this difference has aminor effect on our results. The assumption that the com-panion is at its maximal (Roche lobe) size enables us tocalculate an upper limit of the mass-loss rate as a functionof the measured orbital period P orb and mass m . Specifically,using equation (12), we can relate the loss of orbital energyto the pulsar’s luminosity: GM Û ma = − . η L γ . (13)Several earlier studies either assumed that η ∼ . (Stevens et al. 1992; Benvenuto et al. 2012) or left it as afree parameter (Chen et al. 2013; Benvenuto et al. 2014). We now apply the theory developed in Section 2 to evalu-ate η more precisely. First, we determine the relevant windregime by calculating the characteristic temperature T ch . Us-ing equations (7) and (11) kT ch = µ / ( σ T FR ) / , (14)and its ratio to the virial temperature is given by T ch T g = . (cid:18) L γ L ⊙ (cid:19) / (cid:18) m − M ⊙ (cid:19) − / (cid:18) P orb h (cid:19) / , (15)where we replace R with the effective Roche-lobe radius andscale to typical black widow values ( L ⊙ is the solar luminos-ity). Equation (15) indicates that black widow systems aretypically (see Table 2) in the cold wind regime and there-fore (see Fig. 1) evaporate somewhat less efficiently thanassumed by Ruderman et al. (1989a,b). Next, we substitute p from equation (1) into equation (10), which is appropri-ate for overflowing systems in the cold regime. Finally, wefind the efficiency by comparing to equation (12): η ∼ Ξ ′ v g c (cid:18) T ch T g (cid:19) / == . × − (cid:18) L γ L ⊙ (cid:19) / (cid:18) m − M ⊙ (cid:19) / (cid:18) P orb h (cid:19) − / . (16)Black widow evaporation efficiencies are low, so the evapo-ration time-scales are long: t evap = m | Û m | = . η GMmL γ a == Gyr (cid:18) L γ L ⊙ (cid:19) − / (cid:18) m − M ⊙ (cid:19) / (cid:18) P orb h (cid:19) − / . (17)By comparing equations (15) and (17), we find the usefulrelation t evap = . Gyr (cid:18) T ch T g (cid:19) − , (18)which indicates that black widow systems are in the coldregime when t evap > . Gyr.In Fig. 2 we plot P orb ( m , L γ ) lines (dashed black) forwhich the evaporation time t evap = Gyr. Whether or nota pulsar is able to evaporate its companion over a reason-able time-scale depends critically on its γ -ray luminosity. Wepopulate the figure with the observed black widow sample(Table 2), which is limited to pulsars with measured spin-down rates, from which L γ may be estimated. Specifically,Abdo et al. (2013) find that typically ∼ . of a millisecondpulsar’s total spin-down power L PSR is emitted as γ rays inthe 0.1–100 GeV Fermi band (their fig. 10). Following the Originally, van den Heuvel & van Paradijs (1988) introduced η ∼ . by misciting Ruderman et al. (1989b): both papers de-fine an efficiency f , but the definitions differ by a factor of order v g / c , with the more physical Ruderman et al. (1989b) estimatebeing the smaller of the two. In particular, this relation holds, on average, for the 7black widow systems in our sample that are also a subset ofthe Abdo et al. (2013) sample of 40 millisecond pulsars (theirtable 10; see also 3 additional systems in Barr et al. 2013;MNRAS000
Fermi (Abdo et al. 2013) data suggest that millisecond pul-sar emission might peak around MeV. For these reasons, weestimate that Γ ∼ σ T F .While pair production has a larger cross section thanCompton scattering at the highest photon energies, the re-sulting e ± pairs have to deposit their energy high enoughin the companion’s atmosphere to drive an outflow. A sim-ilar challenge is faced by the TeV e ± wind that poten-tially carries a large fraction of the pulsar’s spin-down power(Ruderman et al. 1989a). According to one scenario, a ∼ G magnetic field may convert the energy of TeV particlesinto MeV photons by synchrotron radiation (Kluzniak et al.1988; Phinney et al. 1988; Ruderman et al. 1989a). Suchsecondary photons can be accounted for by adjusting F ap-propriately. In Section 3 we consider only the direct photonswhen calculating the ablating flux, possibly underestimating F by a factor of a few. Following previous studies (Stevens et al. 1992;Benvenuto et al. 2012, 2014; Chen et al. 2013; Jia & Li2015, 2016; Liu & Li 2017), we define the mass-lossefficiency η as Gm Û mR ≡ − η L γ (cid:18) Ra (cid:19) , (12)where L γ = π a F is the pulsar’s γ -ray luminosity and a isits separation from the companion. η measures the fractionof the incident radiation energy that is invested in overcom-ing the companion’s gravity. As discussed in Section 2.4, weonly consider the γ -ray luminosity that directly couples withthe companion’s upper atmosphere to drive a wind, and notthe total spin-down luminosity. We adopt the same beamingfactor (i.e. no beaming) as Abdo et al. (2013) for consistencywith their L γ measurements.From here onward, we assume that companions arealways close to filling their Roche lobes (Benvenuto et al.2012, 2015; Draghis et al. 2019), so R / a ≃ . ( m / M ) / (Eggleton 1983) where M is the pulsar’s mass. Appar-ently, black widow pulsars occupy the high end of theobserved neutron star mass distribution, with M & M ⊙ (van Kerkwijk et al. 2011; Romani et al. 2015; Linares 2019, M ⊙ denotes the solar mass), presumably due to precedingmass transfer from their companions. Here, we adopt in-stead the canonical M = . M ⊙ , for consistency with pre-vious studies and existing data bases; this difference has aminor effect on our results. The assumption that the com-panion is at its maximal (Roche lobe) size enables us tocalculate an upper limit of the mass-loss rate as a functionof the measured orbital period P orb and mass m . Specifically,using equation (12), we can relate the loss of orbital energyto the pulsar’s luminosity: GM Û ma = − . η L γ . (13)Several earlier studies either assumed that η ∼ . (Stevens et al. 1992; Benvenuto et al. 2012) or left it as afree parameter (Chen et al. 2013; Benvenuto et al. 2014). We now apply the theory developed in Section 2 to evalu-ate η more precisely. First, we determine the relevant windregime by calculating the characteristic temperature T ch . Us-ing equations (7) and (11) kT ch = µ / ( σ T FR ) / , (14)and its ratio to the virial temperature is given by T ch T g = . (cid:18) L γ L ⊙ (cid:19) / (cid:18) m − M ⊙ (cid:19) − / (cid:18) P orb h (cid:19) / , (15)where we replace R with the effective Roche-lobe radius andscale to typical black widow values ( L ⊙ is the solar luminos-ity). Equation (15) indicates that black widow systems aretypically (see Table 2) in the cold wind regime and there-fore (see Fig. 1) evaporate somewhat less efficiently thanassumed by Ruderman et al. (1989a,b). Next, we substitute p from equation (1) into equation (10), which is appropri-ate for overflowing systems in the cold regime. Finally, wefind the efficiency by comparing to equation (12): η ∼ Ξ ′ v g c (cid:18) T ch T g (cid:19) / == . × − (cid:18) L γ L ⊙ (cid:19) / (cid:18) m − M ⊙ (cid:19) / (cid:18) P orb h (cid:19) − / . (16)Black widow evaporation efficiencies are low, so the evapo-ration time-scales are long: t evap = m | Û m | = . η GMmL γ a == Gyr (cid:18) L γ L ⊙ (cid:19) − / (cid:18) m − M ⊙ (cid:19) / (cid:18) P orb h (cid:19) − / . (17)By comparing equations (15) and (17), we find the usefulrelation t evap = . Gyr (cid:18) T ch T g (cid:19) − , (18)which indicates that black widow systems are in the coldregime when t evap > . Gyr.In Fig. 2 we plot P orb ( m , L γ ) lines (dashed black) forwhich the evaporation time t evap = Gyr. Whether or nota pulsar is able to evaporate its companion over a reason-able time-scale depends critically on its γ -ray luminosity. Wepopulate the figure with the observed black widow sample(Table 2), which is limited to pulsars with measured spin-down rates, from which L γ may be estimated. Specifically,Abdo et al. (2013) find that typically ∼ . of a millisecondpulsar’s total spin-down power L PSR is emitted as γ rays inthe 0.1–100 GeV Fermi band (their fig. 10). Following the Originally, van den Heuvel & van Paradijs (1988) introduced η ∼ . by misciting Ruderman et al. (1989b): both papers de-fine an efficiency f , but the definitions differ by a factor of order v g / c , with the more physical Ruderman et al. (1989b) estimatebeing the smaller of the two. In particular, this relation holds, on average, for the 7black widow systems in our sample that are also a subset ofthe Abdo et al. (2013) sample of 40 millisecond pulsars (theirtable 10; see also 3 additional systems in Barr et al. 2013;MNRAS000 , 1–11 (2015)
S. Ginzburg and E. Quataert
Table 2.
Black widow sample from the ATNF Pulsar Catalogue (Manchester et al.2005), version 1.61 (September 2019). We include all systems with pulsar periods P PSR < ms (millisecond pulsars are a distinctpopulation, see Lorimer 2008; Manchester 2017), minimum companion masses m sin i < × − M ⊙ (distinguishing black widows from themore massive redbacks; see Chen et al. 2013; Roberts 2013, who discuss this bimodality), and measured spin-down rates Û P PSR from whichspin-down luminosities L PSR are inferred (assuming a pulsar moment of inertia I = g cm ; see Manchester et al. 2005). We excludeJ1737-0811 because it is a wide binary ( P orb = days; Boyles et al. 2013). Each system’s characteristic temperature T ch , evaporationefficiency η , and evaporation time-scale t evap are derived using equations (15)–(17) by assuming a median inclination angle i = ◦ ,and a γ -ray luminosity L γ = . L PSR (fig. 10 in Abdo et al. 2013). The two rightmost columns show each system’s magnetic brakingtime-scale t mag , given by equation (22), for two choices of the companion’s magnetic field B : a constant B = G and a variable B = ( P orb / h ) − . G. The sample is sorted by increasing t evap PSR Name P PSR m sin i P orb L PSR T ch / T g η P PSR / Û P PSR t evap t B = t var B mag (ms) ( − M ⊙ ) (h) ( L ⊙ ) ( − ) (Gyr) (Gyr) (Gyr) (Gyr)J1701-3006F 2.3 2.1 4.9 190 0.88 4.5 0.33 0.58 1.9 2.3J0024-7204P 3.6 1.7 3.5 140 0.73 4.3 0.17 0.85 3.3 2.7J1701-3006E 3.2 3.0 3.8 92 0.44 3.9 0.33 2.4 3.8 3.4J1959+2048 1.6 2.1 9.2 41 0.36 2.4 3.0 3.4 1.7 4.3J2115+5448 2.6 2.2 3.2 44 0.30 3.1 1.1 5.0 6.3 4.7J1513-2550 2.1 1.6 4.3 23 0.24 2.3 3.1 8.1 5.8 5.9J0024-7204R 3.5 2.6 1.6 36 0.21 3.4 0.74 11 17 5.6J1731-1847 2.3 3.3 7.5 20 0.18 2.1 2.9 15 3.2 6.2J1311-3430 2.6 0.82 1.6 13 0.17 2.1 3.9 16 26 8.3J0024-7204O 2.6 2.2 3.3 17 0.16 2.2 2.8 18 9.6 7.2J1446-4701 2.2 1.9 6.7 9.5 0.13 1.5 7.1 25 5.0 8.5J1518+0204C 2.5 3.7 2.1 17 0.11 2.6 3.0 35 17 7.8J2241-5236 2.2 1.2 3.5 6.4 0.11 1.5 10 37 13 10J2214+3000 3.1 1.3 10 4.9 0.11 1.1 6.7 37 3.9 11J2234+0944 3.6 1.5 10 4.4 0.097 1.1 5.7 49 4.1 11J1641+8049 2.0 4.0 2.2 11 0.083 2.3 7.2 66 20 9.5J0023+0923 3.1 1.6 3.3 4.1 0.071 1.3 8.5 91 17 13J1745+1017 2.7 1.4 18 1.5 0.056 0.65 31 140 3.3 17J1544+4937 2.2 1.7 2.9 3.1 0.056 1.3 23 150 23 15J0610-2100 3.9 2.1 6.9 2.2 0.048 0.95 9.9 200 9.3 16J1719-1438 5.8 0.11 2.2 0.41 0.045 0.51 23 220 67 31J0636+5129 2.9 0.69 1.6 1.5 0.045 1.0 26 230 63 21J2322-2650 3.5 0.074 7.8 0.14 0.036 0.26 190 360 21 42J2017-1614 2.3 2.6 2.3 2.0 0.033 1.2 30 420 38 19J2051-0827 4.5 2.7 2.4 1.4 0.026 1.1 11 680 43 23J1836-2354A 3.4 1.7 4.9 0.62 0.021 0.66 46 1000 25 29 discussion in Section 2.4, we estimate that the MeV pho-tons, which drive the outflow from the companion, carry acomparable share of the luminosity.Table 2 and Fig. 2 exhibit a three orders of magnitudediversity in estimated evaporation time-scales. Notwith-standing the order-unity uncertainties in our analysis, it isevident that while some black widows may have significantlyevaporated their companions by γ -ray radiation, many oth-ers have evaporation time-scales that are longer than theage of the universe. Pulsar spin-down times P PSR / Û P PSR ,which are typically several Gyrs (Table 2), may set aneven more stringent constraint. Interestingly, the originalblack widow pulsar, PSR J1959+2048 (a.k.a B1957+20; seeFruchter et al. 1988), whose discovery sparked much of theearly theoretical work on the subject, has one of the short-est estimated time-scales t evap = . Gyr. The two extremelylow mass companions to PSR J1719-1438 (Bailes et al. 2011)and to PSR J2322-2650 (Spiewak et al. 2018), on the other
Bhattacharyya et al. 2013; Ray et al. 2013, which have similar L γ / L PSR fractions). We have chosen to also include black widowswithout direct L γ measurements, increasing our sample to a totalof 26 systems. hand, have two of the longest time-scales, with t evap = Gyr and t evap = Gyr, respectively. Despite the widerange of computed evaporation times t evap , all of our systemsare in the cold wind regime T ch < T g (see Table 2).Even if our analysis underestimates L γ by factors of afew (see Section 2.4), the large scatter in pulsar spin-downluminosities L PSR (Table 2) combined with the strong de-pendence t evap ∝ L − / γ in equation (17) indicates that evap-oration is unlikely to be the sole driver of black widow evo-lution. Any correction to L γ or to the evaporation efficiencythat brings the longest values of t evap into agreement with thesystems’ ages would similarly shorten t evap for the other sys-tems and imply a large population of short-lived (sub Gyr)black widows. Such short evolution times can be ruled outbecause isolated field γ -ray millisecond pulsars are roughlyas common as their presumed black widow progenitors, anddo not vastly outnumber them (Lorimer 2008; Abdo et al.2013). In Section 4 we propose an alternative mechanism forthe evolution of black widows that resolves this puzzle. MNRAS , 1–11 (2015) lack widow evolution -3 -2 -1 Figure 2.
The observed black widow sample (Table 2). The com-panion’s nominal mass m is calculated assuming the median incli-nation angle i = ◦ , with the error bars indicating the minimum( i = ◦ ) and 95 per cent probability ( i = . ◦ ) values. Thecolour of each system indicates its estimated evaporation time t evap , which is also given in Table 2. Dark colours correspond tosystems in which the pulsar has plausibly evaporated the com-panion to its low mass over the system’s lifetime, whereas lightcolours represent systems with evaporation times longer than theage of the universe. The dashed black lines are given by equation(17) with t evap = Gyr and different values of the gamma-rayluminosity L γ . Evaporation is significant for systems that lie leftof the line with the appropriate L γ . In Section 3 we demonstrated that evaporation alone hasdifficulties in explaining how black widow companions losea significant fraction of their mass during the pulsar’s life-time. In principle, any process that removes angular momen-tum from a binary can keep it in stable Roche-lobe overflow(Rappaport et al. 1982), providing an alternative to evapo-ration. As discussed in the introduction, gravitational wave Roche-lobe overflow is stable when / + ξ − α > (e.g.Rappaport et al. 1982; Ginzburg & Sari 2017, and referencestherein), where ξ parametrizes the companion’s mass-radius re-lation R ∝ m ξ , and α is the fraction of the companion’s specificorbital angular momentum that is lost through overflow ( − α presumably returns to the orbit by torques from the disc thatforms around the pulsar). ξ = − / for both fully degenerate coldcompanions and for hot companions that lose mass adiabatically,so we estimate that the overflow is stable as long as at least athird of the angular momentum is conserved. Unstable Roche-lobe overflow would lead to fast destruction of companions ondynamical time-scales; this scenario seems to be at odds with theobserved black widow population (Draghis et al. 2019). -3 -2 -1 Figure 3.
Same as Fig. 2, but showing the magnetic brakingtime-scale t mag , given by equation (22), assuming a companionmagnetic field of B = G ( t mag is also given in the second fromthe right column in Table 2). Compared to Fig. 2, evolution timesare typically shorter, and more importantly, have a much smallerscatter (i.e. similar colour). emission is too slow, so we must look for a different mecha-nism.A well known mechanism to remove angular momentumfrom a magnetized spinning object is magnetic braking by anoutflow which interacts with the co-rotating magnetic field(Weber & Davis 1967). In our case, the companion’s spin istidally locked to its orbit so magnetic braking taps into theorbital angular momentum. Previous studies (Chen et al.2013; Benvenuto et al. 2014) implemented magnetic brak-ing according to the Rappaport et al. (1983) prescription.This prescription, as explicitly stated by Rappaport et al.(1983), was not derived from first principles but rather tai-lored to fit observed rotation velocities of ∼ M ⊙ main se-quence stars. While this braking law might be relevant atearlier stages of their evolution, it cannot be directly ex-trapolated to ∼ − M ⊙ black widow companions. Moreover,whereas main sequence stars power their own winds, the out-flows from black widow companions are driven by incidentpulsar irradiation (Section 2), and therefore remove massand angular momentum at a very different rate. Here we esti-mate the magnetic braking rate by coupling the companion’smass loss Û m (Section 2) to its magnetic field, for which weassume a surface value of B . Our calculation closely followsThompson et al. (2004), where more details can be found.We assume a split monopole geometry (Weber & Davis1967; Thompson et al. 2004) for which the magnetic fielddecays with radius as B ( r ) = B (cid:18) Rr (cid:19) . (19) MNRAS000
Same as Fig. 2, but showing the magnetic brakingtime-scale t mag , given by equation (22), assuming a companionmagnetic field of B = G ( t mag is also given in the second fromthe right column in Table 2). Compared to Fig. 2, evolution timesare typically shorter, and more importantly, have a much smallerscatter (i.e. similar colour). emission is too slow, so we must look for a different mecha-nism.A well known mechanism to remove angular momentumfrom a magnetized spinning object is magnetic braking by anoutflow which interacts with the co-rotating magnetic field(Weber & Davis 1967). In our case, the companion’s spin istidally locked to its orbit so magnetic braking taps into theorbital angular momentum. Previous studies (Chen et al.2013; Benvenuto et al. 2014) implemented magnetic brak-ing according to the Rappaport et al. (1983) prescription.This prescription, as explicitly stated by Rappaport et al.(1983), was not derived from first principles but rather tai-lored to fit observed rotation velocities of ∼ M ⊙ main se-quence stars. While this braking law might be relevant atearlier stages of their evolution, it cannot be directly ex-trapolated to ∼ − M ⊙ black widow companions. Moreover,whereas main sequence stars power their own winds, the out-flows from black widow companions are driven by incidentpulsar irradiation (Section 2), and therefore remove massand angular momentum at a very different rate. Here we esti-mate the magnetic braking rate by coupling the companion’smass loss Û m (Section 2) to its magnetic field, for which weassume a surface value of B . Our calculation closely followsThompson et al. (2004), where more details can be found.We assume a split monopole geometry (Weber & Davis1967; Thompson et al. 2004) for which the magnetic fielddecays with radius as B ( r ) = B (cid:18) Rr (cid:19) . (19) MNRAS000 , 1–11 (2015)
S. Ginzburg and E. Quataert Figure 4.
Each system’s evaporation (grey squares, taken fromFig. 2) and magnetic braking (with a constant B = G; bluecircles, taken from Fig. 3) time-scales as a function of orbital pe-riod. The magnetic braking times and orbital periods are cor-related, with a best fit of t mag ∝ P − . (solid blue line). Avariable magnetic field that increases with the rotation rate as B ∝ Ω . ∝ P − . can remove this correlation (rightmost col-umn in Table 2, best fitted by a horizontal t mag = Gyr).
The escaping mass is forced to co-rotate with the magneticfield until it reaches the Alfv´en radius R A , where the kineticenergy density is comparable to the magnetic energy density ρ v / = B /( π ) . Using equation (19) and | Û m | = π r ρ v , R A is given by B R (cid:18) RR A (cid:19) = | Û m | v ( R A ) . (20)Spherical symmetry is applicable when deriving equation(20) even in the asymmetric Roche-lobe overflowing case(Section 2.3). Although the outflow from the companion’ssurface is dominated by a cap with a polar angle θ g < inthis case, we find below that R A > a ≫ R , allowing the out-flow to symmetrize by the time it reaches R A . On the otherhand, in more realistic geometries, equatorial magnetic fieldlines form closed loops (e.g. fig. 1 in Mestel & Spruit 1987)and are thus unable to shuttle angular momentum to largedistances. In this scenario, only the spherically symmetriccomponent of the wind, which is launched from higher lati-tudes and constitutes only a fraction T ch / T g of the total out-flow, contributes to magnetic braking. Since we find belowthat the magnetic braking time-scale t mag ∝ Û m − / dependsonly weakly on Û m , our results are insensitive to this effect andare given below assuming the simple split monopole geom-etry. The Alfv´en radius marks the transition from a mag-netically dominated flow, which rotates with the companion Disregarding the equatorial L1 cap would slightly alter the scal- and its magnetosphere, to a free radial outflow. Therefore,at the transition v ( R A ) ∼ Ω R A , where Ω = π / P orb is the spin(and through tidal locking also the orbital) angular velocity.As pointed out by Thompson et al. (2004), this relation ap-plies only if Ω R A is larger than the outflow velocity from anon-spinning companion, which is equal to the escape veloc-ity v g in the cold wind regime (Begelman et al. 1983). Usingthe same hierarchy as above Ω R A ≫ Ω R ∼ v g , where the lastsimilarity holds for companions that fill their Roche lobe.Substituting for v ( R A ) in equation (20): (cid:18) R A a (cid:19) = (cid:18) Ra (cid:19) (cid:18) R A R (cid:19) = . π B RP orb M m | Û m | . (21)Since the wind co-rotates with the companion up to R A ,it carries away angular momentum at a rate Û J = Û m Ω R . Thetotal angular momentum of the system is dominated by theorbit J = m Ω a , so the magnetic braking time-scale is t mag = J | Û J | = m | Û m | (cid:18) aR A (cid:19) = ( π ) / . t / MB RP orb ! / == L − / γ P − / (cid:18) m − M ⊙ (cid:19) / (cid:18) B G (cid:19) − / Gyr , (22)where we have substituted t evap from equation (17). In thebottom line of equation (22), L γ is measured in units of L ⊙ and P orb in hours. Equation (22) indicates that in order formagnetic braking to play an important role (i.e. t mag < t evap ),the magnetic field must be strong enough to satisfy R A > a ,justifying our derivation of equation (21).If t mag < t evap , the system evolves, and the compan-ion loses mass, primarily as a result of magnetic brakingwhich keeps the companion in stable Roche-lobe overflow,rather than directly by photoevaporation. Mass is trans-ferred to the pulsar through the L1 point at a very nar-row angle, which is related to the ratio of the companion’sdynamical time to the system’s age (Linial & Sari 2017).The rate of magnetic braking is a function of the evapora-tive outflow, which is launched at a much larger angle θ g and carries mass and angular momentum to a distance of R A . Through this process, the pulsar’s radiation and thewind that it launches affect the system’s evolution indi-rectly. The pulsar may reject—due to its rapid spin or strongwind—the mass transferred to it through Roche-lobe over-flow, and eject it to large distances in a ‘propeller’ mech-anism (Illarionov & Sunyaev 1975). Such mass ejection re-moves angular momentum from the pulsar’s spin, but notfrom the binary orbit, from which the pulsar is tidally decou-pled. Consequently, magnetic braking is powered only by themass lost directly in the evaporative wind, even when thisis outweighed by the Roche-lobe overflow (in stable Roche-lobe overflow, the angular momentum loss rate is, by def-inition, set by external mechanisms, such as gravitationalwave emission or magnetic braking, and not by the Roche-lobe overflow itself). Accordingly, in our notation Û m refersonly to the wind, rather than the total, mass loss rate. ing relations of equation (22) to t mag ∝ L − / γ P − / m / B − / andwould require a somewhat stronger magnetic field B = G toreach the same median t mag = Gyr for our sample.MNRAS , 1–11 (2015) lack widow evolution The magnetic braking time-scale depends weakly on Û m asa result of competing effects: a stronger wind carries awaymore mass, but also constricts the magnetosphere accordingto equation (21), R A ∝ Û m − / , reducing the specific angularmomentum lost. From equation (22), t mag ∝ Û m − R − ∝ Û m − / ,which is a much shallower function of Û m than direct evap-oration t evap ∝ Û m − . This weak dependence on Û m , and byextension on the wind-driving luminosity L γ , may help toreduce the scatter in characteristic time-scales inferred forobserved systems (Fig. 2).In Fig. 3 we re-plot the observed black widow sample,this time coloured according to the magnetic braking time-scale t mag , assuming the same magnetic field of B = G forall companions. This value was chosen so that the median t mag = Gyr. Estimating the magnetic field strength ofblack widow companions from first principles is beyond thescope of this work. Nonetheless, we note that Hot Jupitershave similar radii and exhibit fields similar in strength to ournominal B (Yadav & Thorngren 2017; Cauley et al. 2019).Typical black widow companions are more massive, rotatefaster, and are subject to stronger irradiation—all poten-tially increase B (Christensen 2010). In addition, as brieflymentioned in Section 2.4, G companion fields were his-torically invoked as a mechanism to produce MeV pho-tons through synchrotron radiation (Kluzniak et al. 1988;Ruderman et al. 1989a).As anticipated, Fig. 3 displays a significantly smallerscatter in time-scales compared to Fig. 2. This result indi-cates that magnetic braking, rather than direct evaporation,is a more likely candidate for the dominant evolutionaryprocess that sculpts the observed black widow population. Moreover, whereas the scatter of evaporation times in Fig.2 does not seem to correlate with either mass or orbital pe-riod, Fig. 3 shows some correlation between t mag and the or-bital period: systems on shorter orbits seem to evolve moreslowly. The correlation is easier to notice in Fig. 4, wherewe plot both time-scales as a function of P orb . A possibleexplanation of this trend is that the magnetic field B is notconstant, as we have assumed, but in fact increases with thecompanion’s rotation rate Ω (which is tidally locked to itsorbital motion). From equation (22), a relation B ∝ Ω . can remove the trend in Fig. 4 and infer a population with areasonable scatter, given the order-unity uncertainties of ouranalysis, around a median lifetime of 10 Gyr (rightmost col-umn in Table 2). Such a positive correlation between B and Ω is in accord with many (but not all) dynamo scaling lawsthat have been suggested in the literature (see Christensen2010, for a review; specifically their table 1). The last decade has seen a substantial increase in thepopulation of observed black widow pulsars, with somesystems very different from the original one detected byFruchter et al. (1988). The growing number of black widows Systems evolve on the shorter of the two time-scales. For severalsystems, t evap . t mag by factors of a few (or even less), but in mostcases, t mag < t evap (Table 2). may allow us to statistically check which mechanisms con-trol their formation and evolution. One mechanism, whichwas suggested early on, is the evaporation (ablation) of thelow-mass companion by the high energy radiation of its hostpulsar (Kluzniak et al. 1988; Phinney et al. 1988). In recentyears, evaporation has been implemented in binary evolutioncodes with a simple prescription, which relates the evapo-ration rate to the pulsar’s spin-down luminosity, introduc-ing an unknown efficiency parameter (Benvenuto et al. 2012;Chen et al. 2013).Here, we calculated the evaporation efficiency by study-ing the hydrodynamical Parker (1958) wind launched off thecompanion’s atmosphere. The wind’s structure can be solvedby considering the interplay between thermal and dynamicalprocesses, encapsulated by the characteristic temperature T ch ∝ F / which can be reached in a sound crossing timegiven the incident flux F (Begelman et al. 1983). Black wid-ows are in the ‘cold wind’ regime, characterised by T ch < T g (the companion’s virial temperature). In this regime, grav-ity retards the flow and the mass loss rate does not reachthe maximum value assumed by Basko et al. (1977) andRuderman et al. (1989a,b). The pulsar’s tidal force weakensthe effective surface gravity on companions that fill theirRoche lobe. The outflow in this case is stronger when com-pared to detached companions (but still below the maxi-mum; see Fig. 1), and it is dominated by a cap around theL1 Lagrange point with a polar angle θ g ∼ ( T ch / T g ) / . Thewind launched from this cap is in the ‘intermediate’ regime(Table 1). The time-scale to fully evaporate the companionis t evap = . ( T ch / T g ) − Gyr (equation 18).We applied our model to a sample of 26 black widows,considering each system’s orbital period, companion mass,and pulsar spin-down luminosity. We put a lower limit on t evap by assuming that companions fill their Roche lobe andmaximize their radius, consistent with recent observationalconstraints (Draghis et al. 2019). Even with this assump-tion, most companions evaporate on a time-scale that islonger than either the age of the universe or the pulsar’sspin-down time (Table 2). This conundrum cannot be re-solved by a simple multiplicative correction: the large dis-persion in measured spin-down luminosities, combined withthe strong dependence on the incident flux t evap ∝ F − / ,leads to a three orders of magnitude scatter in evaporationtime-scales, with no apparent trend (Fig. 2). Any arbitraryenhancement of the evaporation efficiency that expedites theevolution of systems with the longest t evap , or even assum-ing that the spin-down luminosity directly couples to thecompanion with 100 per cent efficiency, would also shortenthe shortest t evap , implying a large population of short livedblack widows that quickly destroy their companions. Such apopulation is inconsistent with the relative numbers of ob-served black widows and isolated millisecond pulsars. Weconclude that evaporation on its own is not responsible forthe evolution of most black widows and their eventual trans-formation into isolated millisecond pulsars.Although the pulsar’s irradiation is too weak to directlyevaporate the companion on a relevant time-scale, the evap-orative wind may couple to the companion’s magnetic fieldand remove angular momentum from the binary system.This magnetic braking can keep the companion in stableRoche-lobe overflow and reduce its mass at a much fasterrate. While a stronger wind removes more mass, it also de- MNRAS , 1–11 (2015) S. Ginzburg and E. Quataert creases the Alfv´en radius and by extension the specific an-gular momentum carried away. These competing effects leadto a weak dependence of the magnetic braking time-scale onthe wind strength and on the flux t mag ∝ t / ∝ F − / , sig-nificantly reducing the scatter in evolution times (Fig. 3).Quantitatively, a companion field of about B ∼ G canexplain how the observed population evolves on a time-scaleof ∼ Gyr; the scatter can be further minimized if themagnetic field increases with the companion’s spin rate as B ∝ Ω . . The shorter time-scales and smaller scatter aremore compatible with the combined population of observedblack widows and isolated millisecond pulsars.Previous studies tried to directly constrain theevaporative wind Û m by fitting an ‘intra-binary shock’model to optical observations of black widow companions(Romani & Sanchez 2016). The idea is that the collisionbetween winds emanating from the pulsar and its com-panion forms a shock which reprocesses the pulsar’s spin-down luminosity and heats the companion non-uniformly.In principle, the shape of the shock, and hence the rela-tive strength of the two winds, can be deduced from themodulation in the companion’s light curve. In practice, how-ever, Romani & Sanchez (2016) find a weak dependence on Û m which is degenerate with the other parameters of theirmodel. Recently, Draghis et al. (2019) fit a direct heatingmodel (disregarding the intra-binary shock) to the opticallight curves of nine black widow systems. They find thateight of the systems fill or almost fill their Roche lobes, with(radius) filling factors of 0.7–1 (their table 3, the ninth sys-tem has a filling factor of 0.5). Given the uncertainties of thefitting model, it is plausible that all observed systems com-pletely fill their Roche lobes, as also implied by our combinedmagnetic braking and Roche-lobe overflow scenario. ACKNOWLEDGEMENTS
We thank Jeremy Hare, Adam Jermyn, Sterl Phinney, RogerRomani, Re’em Sari, and Ken Shen for discussions. SG issupported by the Heising-Simons Foundation through a 51Pegasi b Fellowship.
REFERENCES
Abdo A. A., et al., 2013, ApJS, 208, 17Ablimit I., 2019, ApJ, 881, 72Bailes M., et al., 2011, Science, 333, 1717Barr E. D., et al., 2013, MNRAS, 429, 1633Basko M. M., Sunyaev R. A., 1973, Ap&SS, 23, 117Basko M. M., Hatchett S., McCray R., Sunyaev R. A., 1977, ApJ,215, 276Begelman M. C., McKee C. F., Shields G. A., 1983, ApJ, 271, 70Benvenuto O. G., De Vito M. A., Horvath J. E., 2012, ApJ,753, L33Benvenuto O. G., De Vito M. A., Horvath J. E., 2014, ApJ,786, L7Benvenuto O. G., De Vito M. A., Horvath J. E., 2015, ApJ,798, 44Bhattacharyya B., et al., 2013, ApJ, 773, L12Boyles J., et al., 2013, ApJ, 763, 80Buff J., McCray R., 1974, ApJ, 189, 147B¨uhler R., Blandford R., 2014, Reports on Progress in Physics,77, 066901 Cauley P. W., Shkolnik E. L., Llama J., Lanza A. F., 2019,Nature Astronomy, p. 408Chen H.-L., Chen X., Tauris T. M., Han Z., 2013, ApJ, 775, 27Christensen U. R., 2010, Space Sci. Rev., 152, 565Draghis P., Romani R. W., Filippenko A. V., Brink T. G., ZhengW., Halpern J. P., Camilo F., 2019, ApJ, 883, 108Eggleton P. P., 1983, ApJ, 268, 368Eichler D., Levinson A., 1988, ApJ, 335, L67Fruchter A. S., Stinebring D. R., Taylor J. H., 1988, Nature,333, 237Ginzburg S., Sari R., 2017, MNRAS, 469, 278Gotthelf E. V., Bogdanov S., 2017, ApJ, 845, 159Illarionov A. F., Sunyaev R. A., 1975, A&A, 39, 185Jia K., Li X.-D., 2015, ApJ, 814, 74Jia K., Li X.-D., 2016, ApJ, 830, 153Keith M. J., 2013, in van Leeuwen J., ed., IAU SymposiumVol. 291, Neutron Stars and Pulsars: Challenges and Op-portunities after 80 years. pp 29–34 ( arXiv:1210.7868 ),doi:10.1017/S1743921312023095Kluzniak W., Ruderman M., Shaham J., Tavani M., 1988, Nature,334, 225Krolik J. H., 1999, Active galactic nuclei : from the central blackhole to the galactic environment. Princeton University PressKrolik J. H., McKee C. F., Tarter C. B., 1981, ApJ, 249, 422Kuiper L., Hermsen W., Verbunt F., Thompson D. J., Stairs I. H.,Lyne A. G., Strickman M. S., Cusumano G., 2000, A&A,359, 615Levinson A., Eichler D., 1991, ApJ, 379, 359Linares M., 2019, arXiv e-prints, p. arXiv:1910.09572Linial I., Sari R., 2017, MNRAS, 469, 2441Liu W.-M., Li X.-D., 2017, ApJ, 851, 58London R. A., Flannery B. P., 1982, ApJ, 258, 260London R., McCray R., Auer L. H., 1981, ApJ, 243, 970Lorimer D. R., 2008, Living Reviews in Relativity, 11, 8Manchester R. N., 2017, Journal of Astrophysics and Astronomy,38, 42Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ,129, 1993McCray R., Hatchett S., 1975, ApJ, 199, 196Mestel L., Spruit H. C., 1987, MNRAS, 226, 57Parker E. N., 1958, ApJ, 128, 664Parker E. N., 1965, Space Sci. Rev., 4, 666Phinney E. S., Evans C. R., Blandford R. D., Kulkarni S. R.,1988, Nature, 333, 832Rappaport S., Joss P. C., Webbink R. F., 1982, ApJ, 254, 616Rappaport S., Verbunt F., Joss P. C., 1983, ApJ, 275, 713Ray P. S., et al., 2012, arXiv e-prints, p. arXiv:1205.3089Ray P. S., et al., 2013, ApJ, 763, L13Roberts M. S. E., 2013, in van Leeuwen J., ed., IAU Sympo-sium Vol. 291, Neutron Stars and Pulsars: Challenges andOpportunities after 80 years. pp 127–132 ( arXiv:1210.6903 ),doi:10.1017/S174392131202337XRomani R. W., Sanchez N., 2016, ApJ, 828, 7Romani R. W., Filippenko A. V., Cenko S. B., 2015, ApJ, 804, 115Romani R. W., Graham M. L., Filippenko A. V., Zheng W., 2016,ApJ, 833, 138Ruderman M., Shaham J., Tavani M., 1989a, ApJ, 336, 507Ruderman M., Shaham J., Tavani M., Eichler D., 1989b, ApJ,343, 292Rybicki G. B., Lightman A. P., 1979, Radiative processes in as-trophysics. Wiley-VCHSpiewak R., et al., 2018, MNRAS, 475, 469Stevens I. R., Rees M. J., Podsiadlowski P., 1992, MNRAS,254, 19PThompson T. A., Chang P., Quataert E., 2004, ApJ, 611, 380Weber E. J., Davis Leverett J., 1967, ApJ, 148, 217Yadav R. K., Thorngren D. P., 2017, ApJ, 849, L12MNRAS , 1–11 (2015) lack widow evolution van Kerkwijk M. H., Breton R. P., Kulkarni S. R., 2011, ApJ,728, 95van den Heuvel E. P. J., van Paradijs J., 1988, Nature, 334, 227This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000