PULSAR BINARY BIRTHRATES WITH SPIN-OPENING ANGLE CORRELATIONS

When inverting relation (3) to calculate a birthrate ${\cal R}_i$, we determine N_psr following the procedure originally described in KKL. In order to examine the effects of newly discovered pulsars and estimated pulsar beaming fractions, we considered the same surveys we used in PSC. This includes all surveys listed in KKL and three more surveys, i.e., the Swinburne intermediate-latitude survey using the Parkes multibeam system (Edwards et al. 2001), the Parkes high-latitude survey (Burgay et al. 2003), and the mid-latitude drift-scan survey with the Arecibo telescope (Champion et al. 2004). Explicitly, N_psr is the most likely number of a given pulsar binary population in our Galaxy. We estimate N_psr via a synthetic survey of 10⁶ similar pulsars pointing toward us, as described in KKL. If N_det synthetic pulsars are found in our virtual survey, we estimate N_psr = 10⁶/N_det; this ratio agrees with the slope labeled α in KKL and defined in their Equation (8). Within Poisson error of a few % (=$1/\sqrt{N_{\rm det}}$), our results are consistent with a reanalysis of previous simulations (cf. Table 1 in PSC, from KKL06 and KKL).

In the relation (3), the effective lifetime τ_i encapsulates any and all factors needed to convert between the present-day number $\bar{N}_i$ and a birthrate $f_b^{-1}{\cal R}$ of pulsars that emit along our line of sight. For example, pulsars with a long visible lifetime τ do not reach their equilibrium number (${\cal R} \tau$), instead accumulating steadily with time. Assuming a steady birthrate, the effective lifetime for a binary containing a pulsar i is the smaller of the Milky Way's age and the pulsars' lifetime:

$$\begin{equation} \tau _i = \mbox{\rm min}(\tau _{{\rm mw}},\tau _{\rm age}+\mbox{\rm min}(\tau _{{\rm mrg}},\tau _{d})). \end{equation} \tag{ 4 }$$

The lifetime of a pulsar binary is defined as the sum of the current age of the system (τ_age) and an estimate of the remaining detectable lifetime. For the present age of PSRs B1913+16 and B1534+12, we use the spin-down age of a pulsar, an estimate based on extrapolating pure magnetic-dipole spin-down backward from present to some high initial frequency (Arzoumanian et al. 1999).

For all other pulsars, we adopt the proper-motion corrected ages presented in Kiziltan & Thorsett (2009). Because the pulsar spends most of its time at or near its current age, the error in the current age should be relatively small for any recycled pulsar, unless its current state is very close to the endpoint of recycling. For pulsar binaries considered in this work, even the most extreme possibility—a true age τ_age = 0—changes lifetimes by less than 30%. The only exception is PSR J0737−3039A (see Section 3.1). For sufficiently tight binaries (orbital periods less than ∼10 hr), the remaining lifetime is usually limited by inspiral through gravitational radiation (Peters 1964). When the second-born pulsars are observed, however, the binary's detectable lifetime is instead limited by how long it radiates, as measured by the time until it reaches the pulsar "death-line" (τ_d) as shown in Figure 1 (Chen & Ruderman 1993; Zhang et al. 2000; Harding et al. 2002; Contopoulos & Spitkovsky 2006). We note that our results are fairly robust to the uncertainties in the death-line: (1) merging timescales are more important in birthrate estimates for recycled pulsars (τ_mrg>τ_d), and (2) the death-line is relatively well determined for the non-recycled pulsars in our sample, see, e.g., the discussion in Section 3.3 of PSRs J1141−6545 and in Section 3.1 of J1906+0746 as well as Figure 1.⁴

Zoom In Zoom Out Reset image size

The reconstructed birthrate for wide binaries (τ_i ≃ O(τ_mw)) is sensitive to variations in the star formation history of the Milky Way. Small volumes of the Milky Way, such as the Hipparchos-scale volume of stars, can have O(1) relative changes in the star formation history fluctuation on O(0.3–3 Gyr) timescales (see, for example, Figure 4 in Gilmore 2001). On the larger scales over which these radio pulsar surveys are sensitive, however, these fluctuations average out; see, for example, observations of open clusters and well-mixed dwarf stars in de la Fuente Marcos & de la Fuente Marcos (2004), Hernández et al. (2001), and references therein. In particular, for lifetimes τ_i < 0.2 Gyr relevant to the most significant tight PSR–NS binaries, the star formation rate is constant to within tens of percent, from Figure 15 in de la Fuente Marcos & de la Fuente Marcos (2004). On longer timescales, observations of other disk galaxies and phenomenological models for galaxy assembly also support a nearly constant star formation rate (see Naab & Ostriker 2006; Schoenrich & Binney (2009); Fuchs et al. (2009) and references therein). For pulsar binaries with τ_i>3 Gyr, observations and models suggest the star formation rate trends weakly upward with time (see, e.g., Aumer & Binney 2009 and Fuchs et al. 2009).

For simplicity and to facilitate comparison with previous results, we adopt a constant star formation rate in most of this paper and figures; our final best estimates (Figure 11) include a small correction for exponentially decaying disk star formation, based on the candidate star formation history

$$\begin{eqnarray} {\cal G} = \frac{\dot{\Sigma }_{*}(t)}{\dot{\Sigma }_*(0)}&\propto & e^{- 0.09t/\,{\rm Gyr}}, \end{eqnarray} \tag{ 5 }$$

where t = 0 is the present, drawn from Aumer & Binney (2009); other proposals, such as profiles expected from a Kennicutt–Schmidt relation (Fuchs et al. 2009), are easily substituted. Assuming negligible delay between star formation and compact binary formation (i.e., ${\cal R}\propto {\cal G}$), Equation (3) for the average number of pulsars seen $\bar{N}$ at present in terms of the present-day birthrate ${\cal R}$ generalizes to

$$\begin{eqnarray} &&\bar{N}_{X, {\rm with\; time}} = f_{\rm b}^{-1} R_X(0) \int _{-\min (\tau _X,\tau _{mw})}^o\!\! {\cal G} \ dt = N_{X} \tau _{i,X} \langle{\cal G}\rangle.\nonumber\\ \end{eqnarray} \tag{ 6 }$$

Based on the candidate star formation history of Equation (5), $\left<\cal G\right>\approx 1.6$ for the longest averaging time. Thus, because there was more star formation available to form the widest, long-lived PSR–NS binaries than if stars formed at a steady rate, the present-day formation rate for these wide PSR–NS binaries is roughly 1/1.6 ≈ 0.65 times smaller than that shown in Figure 8.

Except for two pulsars, the empirical beaming geometry of any pulsar in a binary is not tightly constrained. Further, as discussed below, observations of single pulsars suggest that, even restricting to pulsars with similar evolutionary state (e.g., spin), the pulsar beam is randomly aligned relative to its spin axis. Nonetheless, because of the Poisson statistics of pulsar detection (KKL), only one property of the intrinsic pulsar geometry distribution matters: the fraction F(P_s)[ ≡ f⁻¹_b,eff] of all randomly selected pulsars whose beam crosses our line of sight. Assuming the pulsar beam drops off rapidly—typical models involve Gaussian cones—the fraction of all pulsars F(P_s) of a given spin period and luminosity that emit toward us is well-defined and essentially independent of distance. The main ingredients necessary to calculate F(P_s) are the half-opening angle ρ of the radio beam and the misalignment angle between pulsar's spin and magnetic axes α.

All of our analysis leading to the posterior rate prediction, Equations (1)–(3) generalize to an arbitrary beam geometry distribution upon substituting 1/f_b → F.

To calculate F, for simplicity and following historical convention we assume the bipolar pulsar beam subtends a hard-edged cone with opening angle ρ, misaligned by α < π/2 from the rotation axis. As the pulsar rotates, this beam subtends a solid angle

$$\begin{eqnarray} \frac{4\pi }{f_{\rm b}}& =& {2\pi }\times 2 \int _{\rm {\rm min}(0, \alpha -\rho)}^{\rm {max}(\alpha +\rho,\pi /2)} d\cos \theta, \end{eqnarray} \tag{ 7 }$$

where the beaming correction factor f_b is the inverse of the fraction of solid angle the beam subtends, given opening and inclination angles (ρ, α). This simple model for beam geometry and its correlation with pulsar spin have been explored ever since pulsar polarization data and the rotating vector model made alignment constraints possible (see, e.g., Rankin 1993; Gil & Han 1996; Kramer et al. 1998; Weltevrede & Johnston 2008). In one form or another, this theoretically motivated semi-empirical correlation has been adopted as an ingredient in most models for synthetic pulsar populations (see, e.g., Arzoumanian et al. 1999; Gonthier et al. 2004; Gonthier et al. 2006; Faucher-Giguère & Kaspi 2006; Story et al. (2007); and references therein).⁵

Conversely, fully self-consistent population models which account for space distribution and kicks, luminosity evolution, beam structure and shape evolution, accretion during binary evolution, and spin-down are required, in order to extract all the degenerate parameters that enter into a model by comparing its predictions with observations.

In this paper, we focus on quantifying how much the addition of a spin-dependent opening angle ρ(P_s) combined with a misalignment angle α distribution influences birthrate estimates. As a simple approximation, we adopt α and ρ(P_s) distributions that are consistent with the observed pulse widths at 10% intensity level; see Zhang et al. (2003) and Kolonko et al. (2004). Specifically, we assume the misalignment angle is uniformly distributed between [0,2/π]:

$$\begin{eqnarray} &&{\cal P}(\alpha)d \alpha = 2/\pi d\alpha. \end{eqnarray} \tag{ 8 }$$

Other distributions for α have been proposed, from a random vector ${\cal P}(\alpha)=\sin \alpha$ (strongly disfavored by the high frequency with which low-α pulsars have been detected; most detected pulsars would be orthogonal) to more tightly aligned distributions; see for example Figure 3 in Kolonko et al. (2004). While limited current observations cannot tightly constrain the intrinsic misalignment angle distribution, they strongly suggest tight alignments (low α) should be attained at least as often as a flat distribution implies. In Table 3, we compare our reference model with several alternative misalignment distributions; only for exceptional misalignment assumptions (e.g., p(α) ∝ 1/α) will our final results change significantly.

For pulsars with spin period larger than 10 ms, we calculate ρ(P_s) applying a model consistent with classical observations of isolated pulsars' beam geometry (see, e.g., Weltevrede & Johnston 2008, and references therein), as shown in Figure 2:

$$\begin{eqnarray} \rho (P_{\rm s})&=& 5{^\circ}\!\!.4 P_s^{-0.5} \qquad P_s>10\,{\rm ms.} \end{eqnarray} \tag{ 9a }$$

The available single-pulsar data do not support a compelling model for rapidly spinning young and recycled pulsars. We adopt an ad hoc power-law form as our fiducial choice:⁶

$$\begin{eqnarray} &=& 54^\circ (P_s/10\,{\rm ms}) \qquad P_s<10 \,{\rm ms}. \end{eqnarray} \tag{ 9b }$$

Zoom In Zoom Out Reset image size

In both regions, we allow for a small Gaussian error in ln ρ to allow for model uncertainty; we conservatively adopt σ_ln ρ = ln 1.3 (i.e., "30%" error; Figure 2 shows a 2σ interval about our fiducial choice).⁷ Based on these two independent distributions, we estimate the fraction of pulsars with spin period P that point toward us as

$$\begin{eqnarray} F(P_{\rm s})&\equiv & \left< f_{\rm b}^{-1} \right>\equiv 1/f_{\rm b,eff}. \end{eqnarray} \tag{ 10 }$$

The trend of f_b,eff(P_s) shown in Figure 3 largely agrees with previous estimates of the beaming fraction for nonrecycled pulsars (see, e.g., Tauris & Manchester 1998).

Zoom In Zoom Out Reset image size

In addition to our fiducial choice, we have explored several other short-period beaming models. As a benchmark for comparison, two extreme cases are provided in Figure 2: a "very narrow" (solid) and "very wide" (dashed) opening angle model, where the P^−1/2 truncates at 14° (149 ms) and 54° (10 ms), respectively. For relevant spin periods, these extreme alternatives imply noticeably different amounts of beaming correction from each other and our fiducial model, up to O(×2) (see Figure 3). Nonetheless, our best semi-empirical estimates for binary pulsar birthrates are fairly or highly insensitive to the precise opening angle model ρ(P) adopted for rapidly spinning pulsars.⁸ As described below, for pulsars with 10 ms < P_s < 100 ms, we anchor our theoretical bias with observed geometries of comparable binary pulsars. In this period interval, the choice for ρ(P) effectively serves as an upper bound on plausible ρ (and sets a lower bound on f_b).

Any specific choice of opening angle and misalignment distribution ${\cal P}(\alpha,\rho |P_{\rm s})$ determines that model's probability ${\cal P}(f_b|P_{\rm s})$ that a randomly selected pulsar with spin period P_s has beaming correction factor f_b

$$\begin{eqnarray} {\cal P}(f_b|P_{\rm s})d f_b &=&\int \delta (f_b-f_b(\alpha,\rho)) {\cal P}(\rho |P_{\rm s})d\rho {\cal P}(\alpha)d \alpha.\nonumber\\ \end{eqnarray} \tag{ 11 }$$

Further, because 1/f_b is the probability that a given randomly selected pulsar points toward us, the distribution p_e of f_b among the fraction F of all pulsars aligned with our line of sight is

$$\begin{eqnarray} {\cal P}_{e}(f_b|P_{\rm s})df_b &=& \frac{{\cal P}(f_b|P_{\rm s})/f_b}{F(P_{\rm s})}. \end{eqnarray} \tag{ 12 }$$

For our fiducial beaming model, Figure 4 compares the median value of f_b for the two distributions as a function of ρ(P_s). It shows that the beaming correction factors f_b for the detected population of pulsars can differ substantially from the underlying population, if only because those pulsars with narrow beam coverage (i.e., f_b ≫ 1) are less likely to be detected. Only extremely narrow pulsar beams ρ < 10° will lead to a typical detected pulsar with f_b ≃ 6, comparable to the measured value for known PSR–NS binaries that was adopted in previous analyses as the canonical value (see Figure 3).

Zoom In Zoom Out Reset image size

In a few cases, observations ${\cal O}_i$ of pulsar i provide some information about that pulsar's geometry, such as confidence intervals on α and ρ or even a posterior distribution function ${\cal P}(\alpha,\rho |{\cal O}_i)$. When a unique, superior constraint exists, we could have simply factored this information into the average that defines F = 1/f_b,eff (Equation 9; e.g., α constraints for PSR J1141−6545; the preferred geometries of PSRs B1913+16 and B1534+12).

However, in cases where competing proposals exist (e.g., the two models for PSR J0737-3039A suggested by Demorest et al. 2004 and Ferdman et al. 2008, also reviewed in Kramer & Stairs 2008; see Section 3.1), or where our prediction is extremely sensitive to unavailable priors, such as the mean misalignment angle of PSR J1141−6545 (Manchester et al. 2010; M. Kramer 2008, private communication; Section 3.3), we explicitly provide multiple solutions in the text. When multiple competing predictions are available, our final prediction averages over a range of predictions between them. Explicitly, if ${\cal P}_f(\log f_{\rm b})$ is a distribution in X ≡ log(f_b/f_b,eff) reflecting our a priori model uncertainty in f_b,eff, and ${\cal P}(\log {\cal R}|f_{\rm b})= {\cal R} {\cal P}({\cal R}|f_{\rm b}) \ln 10$ is our posterior estimate for the birthrate ${\cal R}$ given known beaming f_b (cf. Equation (1)), then a posterior estimate that reflects uncertain beaming is

$$\begin{eqnarray} {\cal P}(\log {\cal R})d \log {\cal R}&=&\int dX {\cal P}_f(X+\log f_{\rm b,eff})\nonumber\\ && \times {\cal P}(\log {\cal R} - X |f_{\rm b,eff}). \end{eqnarray} \tag{ 13 }$$

For the many PSR–NS binaries which have spin period between 10 ms and 100 ms, we calculate P(log R) using Equation (11). Specifically, we average P(log R) over the beaming correction factor between f_b,eff (listed in Table 1) and the canonical value of 6, adopting a systematic error distribution P_sys(log f_b) which is uniform in log f_b between these limits. For pulsars with spin period outside of this range, we use Equation (1) with f_b,eff listed in Table 1 or 2.

Table 1. Properties of PSR–NS Binaries Considered in this Work

PSR Name	Ps	$\dot{P_{\rm s}}$	Mpsr	Mc	Porb	e	$f_{\rm {b,\rm obs}}$	$f_{\rm {b,eff}}$	τaage	τmgr	τd	Npsr	C	Refb
	(ms)	10−18 (ss−1)	(M☉)	(M☉)	(hr)				(Gyr)	(Gyr)	(Gyr)		(kyr)
Tight binaries
B1913+16	59.	8.63	1.44	1.39	7.75	0.617	5.72	2.26	0.0653	0.301	4.31	576	111	1,2
B1534+12	37.9	2.43	1.33	1.35	10.1	0.274	6.04	1.89	0.200	2.73	9.48	429	1130	3,4
J0737-3039A	22.7	1.74	1.34	1.25	2.45	0.088		1.55	0.142	0.086	14.2	1403	105	5
J0737−3039B	2770.	892.			2.45	0.088		14.	0.0493		0.039			6
J1756−2251	28.5	1.02	1.4	1.18	7.67	0.181		1.68	0.382	1.65	16.1	664	1821	7
J1906+0746	144.	20300.	1.25	1.37	3.98	0.085		3.37	0.000112	0.308	0.082	192	126	8,9
Wide binaries
J1518+4904	40.94	0.028	1.56	1.05	206.4	0.249		1.94	29.2	>τH	51.0	276	18,700	10,11
J1811−1736	104.18	0.901	1.60	1.00	451.2	0.828		2.92	1.75	>τH	7.9	584	5860	12,13
J1829+2456	41.01	0.053	1.14	1.36	28.3	0.139		1.94	12.3	>τH	43.0	271	19,000	14
J1753−2240c	95.14	0.97	1.25	1.25	327.3	0.303		2.80	1.4	>τH	8.2	270	13,900	15

Notes. For most pulsars, f_b,eff averages over the half-opening angle ρ and misalignment angle α. For PSR J0737−3039B, we adopt the preferred choice for α ≃ 90° and average only over the stated uncertainties in ρ(P_s). For PSRs B1913+16 and B1534+12, where both α, ρ measurements are available, we adopt the values of f_b, obs from Kalogera et al. (2001). The final column is C = τ_i/N_psrf_b, also see, Equations (1) and (4). When numbers are uncertain, this table shows self-consistent fiducial choices. Significant uncertainties are included by explicit convolutions described in the text. Small uncertainties are ignored; for example, our Monte Carlo estimates for N_psr have Poisson sample-size errors of roughly $1/\sqrt{N_{\rm det}}\simeq O(2$%–5%), where N_det = 10⁶/N_psr. ^aWhenever available, we use the spin-down ages corrected for the Shklovskii effects given in Kiziltan & Thorsett (2009). As for PSRs B1913+16 and B1534+12, we adapt the results from Arzoumanian et al. (1999). For PSRs J1906+0746, J1811−1736, J1829+2456, J1753−2240, which are not mentioned in Kiziltan & Thorsett (2009), we adopt the characteristic age as the current age of the pulsar. ^bReferences. (1) Hulse & Taylor 1975; (2) Wex et al. 2000; (3) Wolszczan 1991; (4) Stairs et al. 2002; (5) Burgay et al. 2003; (6) Lyne et al. 2004; (7) Faulkner et al. 2005; (8) Lorimer et al. 2006; (9) Kasian and PALFA consortium 2008; (10) Nice et al. 1996; (11) Janssen et al. 2008; (12) Lyne et al. 2000; (13) Kramer et al. 2003; (14) Champion et al. 2004; (15) Keith et al. 2009. ^cThe nature of the companion of PSR J1753-2240 is not yet clear, and it can be either a WD or NS (Keith et al. 2009). In this work, we assume PSR J1753−2240 is another wide NS–NS binary. Given that its small contribution to the total rate estimates, we note that the nature of the companion would not change the main results shown in this work. The masses shown for PSR J1753−2240 are half the total binary mass. All plausible mass pair choices lead to a merger time >10 Gyr; the masses otherwise do not influence our results.

Download table as:  ASCIITypeset image

Table 2. Properties of Tight PSR–WD Binaries Considered in this Work

PSR Name	Ps	$\dot{P_{\rm s}}$	Mpsr	Mc	Porb	e	$f_{\rm {b,eff}}$	τaage	τmgr	τd	Npsr	C	Refb
	(ms)	10−18 (ss−1)	(M☉)	(M☉)	(hr)			(Gyr)	(Gyr)	(Gyr)		(kyr)
J0751+1807	3.48	0.00779	1.26	0.12	6.32	<10−7	2.62	6.66	9.48	>τH	2404	1588	1,2
J1757−5322	8.87	0.0278	1.35	0.67	10.9	<10−6	1.26	7.16	8.0	145	1082	7335	3
J1141−6545	393.9	4295.	1.3	0.986	4.74	0.172	5.46	0.00145	0.60	0.10	346	53	4,5
J1738+0333	5.85	0.0241	1.7	0.2	8.5	4 × 10−6	1.69	3.71	10.8	>τH	609	9716	6

Notes. For all binaries here, $f_{\rm {b,eff}}$ averages over the half-opening angle ρ and misalignment angle α. Monte Carlo sampling uncertainty in N_psr is roughly $1/\sqrt{N_{\rm det}}\simeq O(3\%\hbox{--}5\%)$. ^aFor PSR J0751+1807, we use the spin-down ages corrected for the Shklovskii effects (Kiziltan & Thorsett 2009). For other pulsars, we use the characteristic age. ^bReferences: (1) Lundgren et al. 1995; (2) Nice et al. 2008; (3) Edwards & Bailes 2001; (4) Kaspi et al. 2000; (5) Bailes et al. 2003; (6) Jacoby 2005.

Download table as:  ASCIITypeset image

Similarly, if the lifetime is uncertain, one can marginalize $P(\log {\cal R})$ over the lifetime τ; if the relative likelihood of different lifetimes ${\cal P}_{\tau }(\log \tau _i)$ is known, then defining Y = log(τ/τ_i),

$$\begin{eqnarray} {\cal P}(\log {\cal R})d \log {\cal R}&=&\int dY {\cal P}_\tau (Y + \log \tau _i) {\cal P}(\log {\cal R} + Y |\tau _i),\nonumber\\ \end{eqnarray} \tag{ 14 }$$

Though usually the lifetime is relatively well determined, being dominated by relatively well-determined merger or death timescale, we do use this expression to marginalize over the considerable uncertainty in the current age of PSR J0737−3039A, based on the proposed range of lifetimes presented by Lorimer et al. (2007; see Figure 5).

Zoom In Zoom Out Reset image size

The pulsar binaries we consider in this work and the estimated f_b,eff based on our standard model are summarized in Tables 1 and 2. Tight PSR–NS binaries contribute to our estimate of the overall PSR–NS birthrate (which, for these short-lived systems, is equivalent to their merger rate). Among this set, both PSRs B1913+16 and B1534+12 have both ρ and α measurements. For PSR B1913+16, we adopt $f_{b,\rm obs}=5.72$ based on ρ = 124 (Weisberg & Taylor 2002)⁹ and α = 156°. For PSR B1534+12, we use $f_{b,\rm obs}=6.04$ based on ρ = 4°.87 and α = 114° (Arzoumanian et al. 1996). (The canonical value of f_b ∼ 6 is obtained from the average of the $f_{b,\rm obs}$ for these pulsars.) Alternative choices for α and ρ have been proposed for both pulsars (Figure 3); combining the most extreme of these options leads to values comparable to our model's preferred values: f_b,eff= 2.2 and 1.8 for PSRs B1913+16 and B1534+12, respectively. As tightly beamed pulsars are unlikely to be seen in our fiducial or tightly beamed model (the solid and dashed curves in Figure 3),¹⁰ we adopt an alternative approach for pulsars with spin periods in the otherwise poorly constrained region 10 and 100 ms, the interval containing most PSR–NS binaries.

Pulsars in tight PSR–NS binaries may have wide opening angles, given their spin periods (see Figure 2). If so, these binaries likely have f_b < 6. However, the observations suggest that PSRs B1913+16 and B1534+12 have narrower beams (see Figure 2, region defined by dotted lines). If these narrower opening angles are characteristic of tight, recycled pulsars in PSR–NS binaries, the appropriate f_b could be closer to 6. In particular, given the similarity in spin period of PSR J0737−3039A to those pulsars and the lack of other constraints in that period interval (see Figure 2), we assume log f_b could take on any value between log 1.5 (the value we estimate in our spin model) and log 6.

Given posterior likelihoods, we could explicitly and systematically include observational constraints on the beaming geometry of PSR J0737−3039A as described earlier; see, e.g., the posterior constraints in Demorest et al. (2004) and Ferdman et al. (2008). Observations support two alternate scenarios. In one, the pulse is interpreted as from a single highly aligned pole (α < 4°). Because of its tight alignment, in this model the beaming correction factor should be large: at least as large as those for binary pulsars (f_b ≃ 6, assuming ρ ∼ 30°, from ρ(P_s)), and potentially larger (f_b ≃ 30 assuming ρ = 10°, based on observed opening angles for PSRs B1913+16 and B1534+12). In the other scenario, favored by recent observations (Ferdman et al. 2008), the pulse profile is interpreted as a double pole orthogonal rotator α ≃ π/2 with a fairly wide beam (ρ ∼ 60°–90°, consistent with ρ(P_s)). This latter case is consistent with our canonical model and leads to a comparable f_b. Comparing with the assumptions presented earlier, so long as we ignore the possibility of tight alignment and narrow beams, our preferred model and uncertainties for PSR J0737−3039A already roughly incorporate its most significant modeling uncertainties. Considering that the contribution from PSR J1906+0746 is comparable with that of the PSR J0737−3039A, our best estimate for the birthrate of merging PSR–NS binaries is not very sensitive to changes in a nearly orthogonal-rotator geometry model for PSR J0737−3039A. However, because we cannot rule out the most extreme scenarios for PSR J0737−3039A, for completeness we also describe implications of a unipolar model: the beam shape constraints summarized by Figure 6 translate to a prior on log f_b that is roughly uniform between log 6 and log 30.

Zoom In Zoom Out Reset image size

In Figure 7, we show that the probability distribution of PSR–NS merger rate with our best estimates for the beaming correction f_b,eff, assuming the reference model of KKL06. $P({\cal R}$)'s follow from convolving together birthrate distributions based on each individual pulsar binary, where those birthrate distributions are calculated as described in previous sections. Including PSR J1906+0746, we found the median PSR–NS merger rate is ≃89 Myr⁻¹, which is smaller than what we predicted in KKL06 (∼123 Myr⁻¹, cf. their peak value is 118 Myr⁻¹), assuming the same τ_age and τ_mrg listed in Table 2, but used f_b,J0737 = 5.9, due entirely to the smaller beaming correction factors for PSR J0737-3039A allowed for in this work. Though our best estimate for the merger rate is slightly smaller than previous analyses, the difference is comparable to the Poisson-limited birthrate uncertainty and much smaller than the luminosity model uncertainty described in the appendix. Being nearly unchanged, our study has astrophysical implications in agreement with prior work such as PSC and O'Shaughnessy et al. (2010).

Zoom In Zoom Out Reset image size

Although the merging timescale is relatively well-defined, we note that our estimate does not include O(30%) uncertainties in the current binary age. In this work, for example, we fix the total age of PSR J0737−3039A to be ∼230 Myr. Figure 5 shows a marginalized P(${\cal R}$) using the age constraints from (Lorimer et al. 2007) between 50 and 180  Myr (These models take into account interaction between the two pulsars; we omit the two most extreme models 2 and 3 with rapid magnetic field decay).

In addition to those in tight orbits, some PSR–NS binaries have wide orbits (orbital period is larger than 10 hr). These wide binaries would never merge through gravitational radiation within a Hubble time. Their estimated death timescales imply that most of the known pulsars in wide binaries will remain visible for time comparable to or in excess of 10 Gyr (see Table 1) and that binaries have accumulated over the Milky Way to their present number. In this work, we assume that PSR J1753−2240 (Keith et al. 2009) is another wide PSR–NS binary. Motivated by the fact that spin periods of pulsars found in wide orbits are comparable to those relevant to tight PSR–NS binaries, we adopt the same techniques for estimating f_b,eff. For PSRs J1811−1736 and J1753−2240, we adopt the f_b,eff values expected from their spin periods. On the other hand, for PSRs J1518+4904 and J1829+2456, we average ${\cal P}({\cal R}|f_{\rm b})$ over a range of ln f_b between the value predicted (ln f_b,eff listed in Table 1) and ln 6. Our estimate for the birthrate of "wide" PSR–NS binaries, shown in Figure 8, is slightly higher than previously published estimates by PSC. The birthrate of wide PSR–NS binaries is not significantly changed due to the discovery of PSR J1753−2240, because of its resemblance to PSR J1811−1736. Contrary to Section 2, the previous analyses had assumed that the pulsar population reached number equilibrium, e.g., O'Shaughnessy et al. (2005); our effective lifetime is 2 times shorter. At the same time, our beaming correction is roughly three times smaller, so the relevant factors mostly cancel; our median prediction is almost exactly 1.5 times lower than the previous estimate.

Zoom In Zoom Out Reset image size

The reconstructed $P({\cal R}$) assuming steady-state star formation has a median at 0.84 Myr⁻¹ and is narrower than the previous work due to the discovery of PSR J1753−2240, with a ratio between the upper and lower rate estimates at 95% confidence interval, ${\cal R}_{95\%}/{\cal R}_{5\%}$, is 8.8 (compared to 13.8 without J1753−2240). As in PSC, the birthrate estimate of "tight" PSR–NS binaries is still roughly 100 times greater than that of "wide" binaries, even though both estimates have been modified through lifetime and f_b,eff corrections. Interestingly, known "wide" PSR–NS binaries were discovered through recycled pulsars exclusively. Following PSC, we interpret the difference in birthrate as evidence for a selection bias: only a few progenitors of PSR–NS binaries undergo enough mass transfer to spin up the PSR, but not enough to bring the orbit close enough to merge.

With the exception of PSR J1141−6545, the lifetimes of tight PSR–WD binaries are limited by their gravitational wave inspiral time (Peters 1964). The net visible lifetimes of these binaries are significantly in excess of the age of the Milky Way, and are strongly recycled with P < 10 ms. Based on comparison with measured isolated pulsars, a likely range for their opening angles and thus f_b,eff can be estimated (see Table 2 for a summary). Combining this information, birthrates for the three millisecond PSR–WD binaries, PSRs J0751+1807, J1757−5322, and the newly discovered J1738+0333, are easily estimated (Figure 9). The birthrate for tight PSR–WD binaries, however, is dominated by the non-recycled pulsar PSR J1141−6545. For this unusual binary, the visible pulsar lifetime is dictated by its death timescale (τ_d), rather than its gravitational-wave merger time (τ_mrg), e.g., Kim et al. (2004).

Zoom In Zoom Out Reset image size

Based on different models for the death line, we expect the total lifetime of PSR J1141−6545 to be between 100and170 Myr; for the purposes of a rate estimate, we allow τ to be logarithmically distributed within these limits. Additionally, recent observations imply its misalignment angle to be 160°⁺⁸₋₁₆ (or equivalently, 12° < α < 36°, which is adapted in this work) at a 68% confidence level (Manchester et al. 2010; M. Kramer 2008, private communication; see Figure 10). Without this alignment information we would already predict fairly substantial beaming: f_b,eff ≈ 5.46 solely based on ρ(P_s) ≈ 86 (Table 2). With a nearly polar beam, however, the beaming could be 10 times tighter. We therefore present two estimates for the birthrate of PSR J1141−6545, contingent on α: our current a priori model, shown as a solid curve in Figure 10; and an estimate that includes tight beaming (dotted curve).

Zoom In Zoom Out Reset image size

Pulse profiles are an essential ingredient in empirical pulsar population modeling. The pulse profile of PSR J1141−6545 has significantly broadened, by a factor 8, since the pulsar was discovered in 2000 (M. Kramer 2008, private communication). This is mainly attributed to the effects of the geodetic precession. Recalculating N_psr for PSR J1141-6545 with the current, broader pulse, the most likely value of N_psr is now ∼1000 instead of 350 that is estimated with the pulse profile in 2000 (see Kim et al. 2004 for the details).¹¹ We expect that an analysis that treats both geodetic precession and the duration of pulsar radio surveys imply a slightly different birthrate for the PSR J1141−6545 than what has been shown in Kim et al. (2004). In this work, however, we use the pulse profile presented in the discovery paper (Kaspi et al. 2000) and assume no evolution in the pulse profile, as we only focus on the effects of the beaming correction factor with various beam geometry models. The effects of the pulse profile evolution of PSR J1141−6545 to the Galactic birthrate of PSR–WD binaries, taking into account more detailed corrections for observational biases such as the Doppler smearing, will be discussed in a separate paper (C. Kim et al. 2010, in preparation). Our preferred birthrate estimate agrees with previously published estimates based on a "fiducial" beaming factor f_b = 6 (PSC; compare C = 53 kyr here with A = 47 kyr in their Table 1 for PSR J1141−6545), as well as with early estimates that adopt f_b = 1 (Kim et al. 2004 predict a peak-probability galactic birthrate of 4f_b Myr⁻¹, based on an N_psr = 400).¹² By coincidence the fiducial beaming factor agrees with our best beaming estimates based on ρ(P) data for comparable-spin pulsars.

Alternate misalignment models. In order to calculate f_b,eff, we adopt a flat distribution as our standard model for the α distribution defined in Equation (8). Several alternative distributions for pulsar misalignment have been proposed; some, however, do not self-consistently account for detection bias. For example, a randomly aligned beam vector on the sphere (p(α)dα = dcos α) implies that nearly all detected pulsars should be nearly orthogonal rotators. The observation of many pulsars with smaller α, e.g., Kolonko et al. (2004), is not consistent with that distribution. On the other hand, a flat distribution and even more centrally concentrated ones (e.g., p(α) = cos αdα) are plausible, considering observational biases can explain the lack of detected pulsars with small α (say, <45° from Figure 3 in Kolonko et al. 2004). But because the expectation defining f_b,eff (Equation (9)) is governed by the smallest typical beaming fraction of detected pulsars, the predictions for f_b,eff are fairly independent of misalignment model, assuming it is not too concentrated near the equator or poles. For this reason, and without more information to quantify uncertainties in a reconstructed α distribution with which to base a comparison (e.g., the selection bias of having an α measurement), we choose a flat distribution for simplicity.

What beaming to use? No one beaming correction factor applies for all circumstances: even for a fixed beam geometry distribution, different "natural" choices f_b make sense for different questions (Figure 4). Beaming depends strongly on spin. Finally, the canonical choice f_b = 6 implicitly requires exceptionally tight beaming, strong alignment, or both (Figure 4 and Table 3). If one number must be adopted, use the spin-dependent relation proposed by Tauris & Manchester (1998) for pulsars with P>100 ms (Figure 3).

For clarity, in the text we discussed the impact of spin-dependent pulsar beaming on the birthrate given a fixed pulsar population model. The results shown in this work are based on our reference model (a power-law distribution for a pseudoluminosity dlog N/dlog L = −1, with 0.3 mJy kpc² as the minimum intrinsic luminosity at 400 MHz, Gaussian distribution in radial direction where the radial scale length is assumed to be R_o = 4 kpc, and exponential function in z, with a scale height of Z_o = 1.5 kpc). See KKL, KKL06 for further details on the pulsar population model. Although the rate estimates are sensitive to model parameters, or the fraction of faint pulsar assumed in a model, the qualitative feature described here are robust with different model assumptions made. As discussed in KKL06, however, the reconstructed pulsar population depends sensitively on the luminosity function model. In this appendix, we quantify the uncertainties attributed to a pulsar luminosity function following KKL06. For example, if the cumulative probability of a luminosity greater than L is modeled by P(>L) = (L_min/L)^1−p, a function with two parameters p and L_min, then empirically we have found the number of pulsar binaries N_psr implied by one detection scales as

$$\begin{eqnarray} N_{\rm psr}&\propto & L_{\rm min} 10^{-1.6 p}. \end{eqnarray} \tag{ A1a }$$

Assuming that a pulsar luminosity function is similar for both millisecond pulsars (P_s < 20 ms) and pulsars found in compact binaries, KKL06 adapt the results shown in Cordes & Chernoff (1997) and estimate that the uncertainty in L_min and p can be described by the two uncorrelated distributions

$$\begin{eqnarray} {\cal P}(L_{\rm min})&=& 1.22 L_{\rm min}(1.7-L_{\rm min}) \quad L_{\rm min}\in [0, 1.7]\,{\rm mJy}\,{\rm kpc}^2\nonumber\\ \end{eqnarray} \tag{ A1b }$$

$$\begin{eqnarray} {\cal P}(p)&=& \frac{1}{\sqrt{2\pi \sigma _p^2}}e^{-(p-2)^2/2\sigma _p^2} \quad \sigma _p =0.12 \quad p \in [1.4, 2.6].\nonumber\\ \end{eqnarray} \tag{ A1c }$$

The above scaling relation and PDF allow us to generalize any result for ${\cal P}(\log {\cal R}|L_{\rm min},p)$ presented in the main text, which assumed a specific pulsar luminosity model L_min, ref = 0.3 mJy kpc² and p_ref = 2.0, to a "global" result $p_g(\log {\cal R})$ that fully marginalizes over pulsar model uncertainties:

$$\begin{eqnarray} p_g(\log {\cal R})&=& \int {\cal P}(\log {\cal R}|\log L_{\rm min},p){\cal P}(\log L_{\rm min})\nonumber\\ && \times{\cal P} (p) d\log L_{\rm min} dp \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} \log Q&\equiv & {-}\log (L_{\rm min}/L_{\rm min,ref}) + 1.6(p-p_{{\rm ref}}) \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} {\cal P}(\log {\cal R}|L_{\rm min},p)&=& {\cal P}(\log {\cal R} - Q | L_{\rm min,ref}, p_{{\rm ref}}). \end{eqnarray} \tag{ A4 }$$

Therefore, changing variables to z = 1.6(p − p_ref), we find the fully marginalized PDF ${\cal P}_g(\log {\cal R})$ follows from the results presented here via convolution with an ambiguity function ${\cal A}$:

$$\begin{eqnarray} {\cal P}_g(\log {\cal R})&=& {\cal P}(\log {\cal R}|L_{\rm min,ref}, p_{\rm ref}) * {\cal P}(-\log L_{\rm min}) * {\cal P}(z)\nonumber\\ \end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray} {\cal A}&\equiv & {\cal P}(-\log L_{\rm min}) * {\cal P}(z). \end{eqnarray} \tag{ A6 }$$

Figure 11 shows our estimate of the global ambiguity function, given the model of Equation (A1). This distribution is roughly as wide as the distributions shown in the text, whose width is limited by Poisson statistics of the number of observed binaries. Thus, without a better resolved luminosity model, even several new binary pulsars would not reduce the overall uncertainty in the rate estimate.

Though simple, a global power-law luminosity distribution leads to predictions that depend sensitively on the low-luminosity cutoff; see Equation (A1). Recent detailed pulsar population synthesis studies by Faucher-Giguère & Kaspi (2006, henceforth FK) suggest a simpler, lognormal luminosity function fits the whole pulsar population

$$\begin{eqnarray} p(\log L)d \log L &=& \frac{1}{\sqrt{2\pi \sigma _{\log L}^2}}e^{-(\log L -\langle \log L\rangle)^2/2\sigma _{\log L}^2} d\log L,\nonumber\\ \end{eqnarray} \tag{ A7 }$$

where σ_log L = 0.9; see their Figure 15. As shown in Figure 12, however, a naive lognormal distribution differs substantially from our canonical model: even adopting <log L> = log 0.3( mJy kpc⁻²), substantially larger than the best-fit model in FK with <log L> = −1.1, a lognormal distribution predicts more bright and fewer faint pulsars than our fiducial power-law luminosity model. Translating to N_psr and birthrates, compared to our fiducial power-law distribution roughly four to six times more binaries must be present in the fiducial lognormal (<log L> = −1.1) for one to be detected, depending on the pulsar spin and width being simulated. For comparison, for each pulsar binary our "extreme" lognormal distribution (<log L> = log 0.3) predicts roughly the same numbers N_psr of pulsars as the fiducial power-law distribution.¹³ Further, the best-fit lognormal distribution from FK is really a spin-dependent luminosity $L\propto P^{\epsilon _P}\dot{P}^{\epsilon _{\dot{P}}}$, with weakly constrained exponents $\epsilon _{\dot{P}}$ and _P. That model implicitly introduces time-dependent selection effects, beyond the scope of this paper. The lognormal model described above is presented as a convenient summary, not a fundamental distribution; FK thus do not describe how reliable its parameters are. Lacking control over the model, we defer a detailed discussion of lognormal luminosity models to a future paper.

Zoom In Zoom Out Reset image size