A 1.05 M_☉ COMPANION TO PSR J2222−0137: THE COOLEST KNOWN WHITE DWARF?

Radio observations of PSR J2222 to measure the Shapiro delay occurred in the last week of 2011 May with the 100 m Robert C. Byrd GBT.¹⁰ We had a 6 hr observation taken around superior conjunction of the binary system augmented by five 2 hr observations at each of the other five Shapiro extrema, all using the Green Bank Ultimate Pulsar Processing Instrument (GUPPI; DuPlain et al. 2008). The 800 MHz of bandwidth centered at 1500 MHz in two orthogonal polarizations was separated into 512 Nyquist-sampled frequency channels of width 1.5625 MHz via a polyphase filter bank. These channels, sampled at 8 bits, provided full polarization information and an effective time resolution of 0.64 μs. Each channel was coherently dedispersed at the nominal dispersion measure of the pulsar (3.27761 pc cm⁻³ at the time, although we later refined this measurement). Each observing session was broken into 30 minute observations of PSR J2222 separated by 60 s calibration scans of the extragalactic radio source 3C 190. The calibration scans were taken in the same mode as the pulsar observations, but also included a 25 Hz noise diode inserted into the receiver.

Data reduction was performed using the PSRCHIVE package (Hotan et al. 2004). Flux calibration used the on- and off-source scans of 3C 190. This was followed by removal of radio frequency interference by the psrzap utility. The calibrated pulse profile determined from the long observation covering conjunction is given in Figure 1. The data were aligned in time using the best ephemeris (below), divided into 16 frequency channels, and refit for dispersion measure and rotation measure using a bootstrap error analysis. We found that the period-averaged flux density varied by a factor of a few over the course of long observations due to scintillation, with an average of 1–2 mJy at 1500 MHz. Individual times of arrival (TOAs) were measured from the folded total-intensity profiles using the frequency domain algorithm in PSRCHIVE (Taylor 1992). A template was created by fitting three Gaussians to the summed pulse profile. From these Gaussian components, we created a noise-free template with the phase of the fundamental component in the frequency domain rotated to zero. The observations were divided into two minute segments, with one TOA measured for each segment. Note that since interstellar scintillation caused the flux to vary considerably, there was a proportional change in the TOA precision that varied over the data set.

Zoom In Zoom Out Reset image size

These data were combined with previous data taken for the discovery observations of PSR J2222 (Boyles et al. 2013) to produce a timing model. We used the "DD" model (Damour & Deruelle 1985, 1986) in TEMPO,¹¹ which incorporates the Shapiro delay. The astrometric data for this model were taken from Deller et al. (2013), and we used the DE421 JPL ephemeris (Folkner et al. 2009). Timing fits with no Shapiro delay were statistically unacceptable, with an rms residual of 9.3 μs (χ² = 4539.4 for 931 degrees of freedom), and a clear Shapiro delay signature was obvious in the residuals (Figure 2). With the Shapiro delay included in the fit the rms residual was 4.2 μs (χ² = 930 for 929 degrees of freedom), with no obvious remaining structure in the residuals (varying the astrometric parameters within the uncertainties from Deller et al. 2013 changed the timing results by ≪1σ). The Shapiro delay determines the inclination of the orbit and the companion mass; this is then combined with the binary mass function to determine the pulsar's mass. Due to the combination of several different and much less precise observing modes from earlier monitoring with the high-precision Shapiro delay campaign, we estimated the timing parameters with a bootstrap error analysis. We give the full timing results, with 1σ error estimates from the bootstrap analysis, in Table 1.

Zoom In Zoom Out Reset image size

Table 1. Fitted and Derived Timing Parameters for PSR J2222−0137

Parameters	Value
Timing Parameters
Spin period (s)	0.032817859053065(3)
Period derivative (s s−1)	5.865(7) × 10−20
Dispersion measure (pc cm−3)	3.2842(6)
Rotation measure (rad m−2)	+2.6(1)
Reference epoch (MJD)	55743
Right ascensionb (J2000)	22:22:05.969101(1)
Declinationb (J2000)	−01:37:15.72441(4)
R.A. proper motionb (mas yr−1)	44.73(4)
Decl. proper motionb (mas yr−1)	−5.68(6)
Parallaxb (mas)	$3.742^{+0.013}_{-0.016}$
Position epochb (MJD)	55743
Span of timing sata (MJD)	55005–55922
Number of TOAsa	943
rms residual (μs)	4.2
Binary parametersd
Orbital period (days)	2.4457599929(3)
Projected semi-major axis (lt-s)	10.8480276(12)
Epoch of periastron (MJD)	55742.13242(0)
Orbital eccentricity	3.8086(15) × 10−4
Longitude of periastron (deg)	119.778(12)
Mass function (M☉)	0.22907971(8)
sin i	0.9985(3)
Companion mass (M☉)	1.05(6)
Derived parameters
Distanceb (pc)	267.3$^{+1.2}_{-0.9}$
Transverse velocityb (km s−1)	57.1$^{+0.3}_{-0.2}$
Orbital inclination i (deg)	86.8(4)
Shklovskii period derivative (s s−1)	4.33(5) × 10−20
Intrinsic period derivativec (s s−1)	1.54(5) × 10−20
Surface magnetic fieldc (109 G)	0.719
Spin-down luminosityc (1031 erg s−1)	1.72
Characteristic agec (Gyr)	33.8
Pulsar mass (M☉)	1.20(14)
Flux density at 1500 MHz (mJy)	1–2

Notes. Values in parentheses are uncertainties on the last digit. For the timing data derived here, the uncertainties were derived from a bootstrap analysis and are quoted at the 1σ level. ^aDuring the initial timing observations we calculated a TOA every 10 minutes. During the new observations described here we calculated a TOA every two minutes. ^bValues are from Deller et al. (2013) and were held fixed for the timing fit. ^cValues are corrected for Shklovskii effect. ^dWe used the "DD" model (Damour & Deruelle 1985, 1986).

Download table as:  ASCIITypeset image

Our data consist of high-quality coherently dedispersed data from an intensive one week campaign and a few other epochs. The remainder of the data were both less precise and less uniform, with a wider range of observation frequency and instrumental setup. This makes it difficult (if not impossible) to robustly constrain long-term secular changes like periastron precession ($\dot{\omega }$; Lorimer & Kramer 2012). Nonetheless, we tried a fit with $\dot{\omega }$ fixed to the value predicted by general relativity (≈008 yr⁻¹). The resulting fit was good, with the rms decreasing to 3.8 μs. The pulsar and companion masses each increased by about 1σ compared to the values in Table 1. Given the small eccentricity and inhomogeneous data set with large gaps we do not believe that fitting for $\dot{\omega }$ is viable at this time, but encourage further long-term monitoring of this system to establish its secular behavior.

We observed the position of PSR J2222 at optical and near-infrared wavelengths, as listed in Table 2. The deepest Keck observations used the red side of the Low-Resolution Imaging Spectrometer (LRIS; Oke et al. 1995) on the 10 m Keck I telescope. The data were reduced using standard procedures in IRAF, subtracting the bias, dividing by flatfields, and combining the individual exposures. The seeing was about 08 in the combined R image, and 07 in the combined I image. We computed an astrometric solution fitting for a shift and separate scales and rotations along each axis (i.e., a six parameter fit) using 100 non-saturated stars identified from the Sloan Digital Sky Survey (SDSS) Data Release 10 (DR10; Ahn et al. 2014), giving rms residuals of 02 in each coordinate. We did photometric calibration relative to SDSS photometry, identifying 23 well-detected, well-separated, non-saturated stars, and transforming from the SDSS filter set to Johnson–Cousins using the appropriate transformation equations.¹² The zero-point uncertainty was <0.01 mag, although there are systematic uncertainties coming from our filter transformations. We see no object at the position of the pulsar (Figure 3); the closest object is about 2'' from the position of the pulsar (about 10σ away) and appears extended (R = 23.1 ± 0.1 and statistical position uncertainties of ±03 in each coordinate). We determined the 3σ upper limits using sextractor (Bertin & Arnouts 1996) to determine the magnitude that gave a 0.3 mag uncertainty (verified with fake-star tests), which we give in Table 2.

Zoom In Zoom Out Reset image size

We observed PSR J2222 in r-band with the Goodman Spectrograph (Clemens et al. 2004) on the 4.1 m Southern Astrophysical Research (SOAR) telescope over two nights in 2013 July. All exposures were dithered and binned by a factor of two in both dimensions. The frames were bias-subtracted and flattened with a dome flat. We then used a median of the data (having masked the scattered-light halos of three saturated stars) from the second night constructed without registration to create a sky flat, which we smoothed with a 20 × 20 pixel boxcar filter. This corrects for larger-scale brightness variations. Cosmic rays were interpolated on individual exposures using the lacosmic routine (van Dokkum 2001). The seeing varied considerably over the course of the observations, going from 11 to 2''. We then shifted each exposure by an integer number of pixels for registration and summed them. The final summed image has an effective seeing of 13 and a total exposure time of 2.6 hr. The photometric zero-point was again computed relative to the SDSS DR10 data, using 31 stars. The astrometric solution was done using six 30 s exposures through http://astrometry.net (Lang et al. 2010). As with the Keck data, we see no object at the position of the pulsar (Figure 3) and give a 3σ upper limit in Table 2.

While they were taken through different filters and with very different instruments/resolutions, we tried combining the Keck R-band and SOAR r-band images using swarp (Bertin et al. 2002). We still see no source at the position of the pulsar. The data are sufficiently different that a limiting flux is difficult to compute, but it could be as much as 0.3 mag fainter than the limits in Table 2.

The near-infrared observations come from the NIRC2 camera¹³ on the 10 m Keck II telescope, and used the Laser Guide Star Adaptive Optics (AO) system (van Dam et al. 2006). The data were taken through thin clouds and the AO corrections were not optimal, resulting in a delivered image quality of 02 FWHM. The images were reduced using a custom pipeline implemented with python and pyraf using dark frames and dome-flats. A sky fringe frame was created by combining dithered images of multiple targets with the bright stars masked. We used SExtractor (Bertin & Arnouts 1996) for the preliminary detection and masking of stars. The fringe frame was subtracted from the flat-fielded data after being scaled to the appropriate sky background level. Before coadding the frames, each frame was corrected for optical distortion using a distortion solution measured for NIRC2.¹⁴ A faint glare has been visible in the lower right (southwest) corner of the NIRC2 wide camera images starting in 2009 August. The shape and amplitude of the glare vary with telescope orientation, resisting correction through surface fitting or modeling. Instead we masked the glare using a triangular region. There was no independent photometric calibration that night, and only a single star is visible on the co-added image. To determine a photometric zero-point, we used photometry for that star from the SDSS DR10. We then employed the empirical main-sequence color relations from Covey et al. (2007), inferring the z − K_s color from the observed g − i color (we ignore differences between K_s and K' filters). For this star (SDSS J222204.76−013658.9) we infer a spectral type of K2.5 and predict K_s = 16.9. We expect zero-point uncertainties of ±0.2 mag or so based on comparison of the other SDSS colors to those predicted using Covey et al. (2007). Again, we see no object at the position of the pulsar, and give 3σ upper limits in Table 2.

Since we do not detect the optical counterpart of the companion, the first inference is that the companion could be a low-mass NS. It would be the lowest mass NS known (Lattimer 2012; Özel et al. 2012; Kiziltan et al. 2013), although it is only a roughly 2–3σ excursion from the mean of the companions in DNS systems (Özel et al. 2012; Kiziltan et al. 2013): rare, given the ≈10 DNS systems, but not impossible.

In that case, its eccentricity of 3.8 × 10⁻⁴ would be a factor of 200 lower than any other DNS system (PSR J1906+0746 has the lowest eccentricity of e = 0.085, although this may be an NS–WD system; Kasian 2012; van Leeuwen et al. 2014). In Figure 4, we show the eccentricity versus component masses for all DNS and NS–WD systems with well-determined masses. In fact, there are three NS–WD systems with higher eccentricities: PSR J1141−6545, which was likely not recycled (Kaspi et al. 2000); PSR J0337+1715, which has had its eccentricity increased by dynamical interactions (Ransom et al. 2014); and PSR J0621+1002 (likely an IMBP, with the eccentricity the result of unstable mass transfer; Phinney & Kulkarni 1994; Camilo et al. 2001).

Zoom In Zoom Out Reset image size

The normal formation scenario for a DNS involves two core-collapse supernova explosions, with the eccentricity the result of the second explosion and its kick, and no final mass-transfer phase to circularize the orbit (e.g., Tauris & van den Heuvel 2006). In contrast, formation via an electron-capture supernova (ECS; Miyaji et al. 1980) could result in a significantly lower NS mass (Schwab et al. 2010; Ferdman et al. 2013) along with a lower supernova kick (Podsiadlowski et al. 2004; van den Heuvel 2004). PSR J2222 has a low transverse velocity (58 km s⁻¹), although higher than some systems thought to be the products of ECSs (given the age of the system, this velocity may be more related to motion in the Galactic potential than birth conditions). This may reflect the velocity dispersion of the progenitor systems. However, the contrast between PSR J2222 and other systems thought to be the results of ECSs (e.g., PSR J1906+0746 or PSR J0737−3039; Ferdman et al. 2013) is extreme, with the ratio of eccentricities above 200 as mentioned previously. In a scenario without a kick we can place an upper limit on the amount of material that could have been ejected by the explosion to (M_PSR + M_c)e = 8 × 10⁻⁴ M_☉ (with M_psr the pulsar mass and M_c the current companion mass; e.g., Bhattacharya & van den Heuvel 1991). This is a much tighter bound than in any of the other systems proposed for this mechanism, and difficult to reconcile with the change in binding energy needed to collapse to an NS $\sim G M_c^2/R_c c^2\approx 0.1\,M_\odot$ (with R_c ≈ 15 km for an NS), presumably released as neutrinos (e.g., Freire & Tauris 2014): this leads to the horizontal line in Figure 4, above which all confirmed DNS systems are found. In order to have a DNS system with such a low eccentricity, we need to invoke increasingly exotic (and perhaps implausible) evolutionary scenarios. For instance, if the system began as a hierarchical triple (Champion et al. 2008; Ransom et al. 2014; Tauris & van den Heuvel 2014), then the inner components could have formed a standard eccentric DNS system early on. Later evolution of the outer member could have led to a circum-binary accretion disk that would have worked to circularize the inner system, after which the outer object would have exploded or otherwise been ejected from the system.

The other possible scenario is that the companion could be a massive WD, making the system an IMBP. Its orbital eccentricity is somewhat high compared to most low-mass binary pulsars of similar periods (based on Phinney 1992), but not nearly as high as a DNS, consistent with an IMBP classification (Camilo et al. 2001). It falls in the locus of other CO WDs in the "Corbet" (binary period versus spin period) diagram in Tauris et al. (2012). The pulsar mass is lower than most pulsar–WD binaries, but is consistent with the short orbital-period IMBP discussed by Ferdman et al. (2010) which may indicate a similar formation mechanism involving a common envelope (Tauris & van den Heuvel 2006).

However, as a WD it would be extremely faint: far fainter than any of the optical companions to IMBPs currently known (van Kerkwijk et al. 2005; Jacoby et al. 2006; Pallanca et al. 2013) or indeed any WD companion to a millisecond pulsar (MSP) with a similar mass (Antoniadis et al. 2011); it is perhaps the faintest WD ever observed. With the apparent magnitude limits from Table 2, we can compute absolute magnitude limits in each band. We use the distance 267 ± 1 pc (Deller et al. 2013), and we estimate the extinction to be A_V = 0.12 mag from Drimmel et al. (2003). In terms of bolometric luminosity, the most constraining limit ends up coming from the R-band data, where we limit M_R > 19.1 (the r-band limit of M_r > 19.2 is very similar, given slight differences in bolometric correction). For comparison, the companion to PSR J1022+1001 with a median companion mass of 0.85 M_☉ has M_R ≈ 14 (Lundgren et al. 1996). In Figure 5, we plot the absolute magnitude against mass for pulsar+WD systems as well as select cool WDs with parallax distances: even compared to the observed truncation of the cooling sequence in old halo globular clusters like NGC 6397 (Richer et al. 2013) or M4 (Bedin et al. 2009), the putative companion is far fainter: at the distance of NGC 6397, our limit of M_R > 19.1 translates to an apparent magnitude of R > 31.6, compared to R ≈ 29, or M_R ≈ 16 for the coolest WDs seen in NGC 6397. Some of the difference comes from the change in radius: a 1.0 M_☉ WD has a radius about 65% of that of a typical 0.6 M_☉ WD, leading to a 1 mag change in brightness at the same effective temperature. However, the difference in Figure 5 is more like 2.5 mag, so the companion to PSR J2222 must also be cooler than the known thick disk/halo WDs.

Zoom In Zoom Out Reset image size

Beyond the absolute magnitude, which is directly computable from observable quantities, we can limit the radius/temperature of a putative WD by using our R-band absolute magnitude limit to constrain the bolometric luminosity. This is more complicated, as it involves atmosphere calculations in an uncertain and poorly tested regime, but it should be reasonably reliable. We use the synthetic photometry and evolutionary models from Tremblay et al. (2011) and Bergeron et al. (2011) for H and He atmospheres, respectively.¹⁵ For isolated WDs pure He atmospheres can be largely excluded because of Bondi-Hoyle accretion from the ISM (Bergeron 2001), and even small amounts of hydrogen mixed into the helium can cause near-infrared flux deficiencies like pure hydrogen (see below; Bergeron & Leggett 2002). However, the binary orbit and MSP wind in this system could have inhibited such accretion and therefore a He atmosphere is possible. In any case, a pure He atmosphere will serve as a limiting case compared to the H models. These models are used to convert the absolute magnitude limits into temperature limits, so for simplicity we use the 1.0 M_☉ models (differences in bolometric corrections as a function of mass are small, <0.05 mag).

The most constraining limit is again from the R-band data, where M_R > 19.1 implies T_eff < 1700 K (see Figure 6) for an H atmosphere. The He-atmosphere models do not extend to sufficiently cool temperatures but stop at T_eff = 3500 K with M_R = 17.7. At lower temperatures the details of the atmospheric physics are rather uncertain, but a blackbody is likely an acceptable approximation (P. Bergeron 2014, private communication). With an He atmosphere an effective temperature <3000 K would be required (Figure 6). The H limits are more constraining since more of the flux appears in the optical regime rather than the near-infrared—a consequence of collisionally induced absorption by molecular H₂ (Bergeron et al. 1995; Hansen 1998). These limits change slightly with mass given the small but finite mass uncertainties, since the radius would change with mass: going to the 1.2 M_☉ H model we can constrain T_eff < 2100 K (at our nominal mass of 1.05 M_☉ the radius of a C/O WD is about 0.0073 R_☉, and it scales as R∝M^−1.6). As inferred from Figure 5, the companion to PSR J2222 would be far cooler than any known WD from other surveys (e.g., Kilic et al. 2012; Catalán et al. 2012, 2013), where the coolest objects tend to have T_eff ≈ 3800 K.

Zoom In Zoom Out Reset image size

However, we cannot exclude such a very cool WD on age grounds. WD cooling curves, which start out having more massive objects warmer at the same age, eventually cross to have more massive objects cooler at the same age (Figure 7; this is also visible in Figure 5). This is because massive WDs crystallize earlier, at a higher T_eff (but at a similar internal temperature), at which point the faster Debye cooling takes over (Mestel & Ruderman 1967; van Horn 1968; Chabrier et al. 2000). Cooling ages for these models may not be reliable, as the impacts of state changes, sedimentation, and chemical processes are not precisely known, and the atmospheres are not trivial to calculate (Montgomery et al. 1999; Chabrier et al. 2000; Althaus et al. 2007; Salaris et al. 2010; Tremblay et al. 2011). But we believe conservatively that the cooling age is close to 10 Gyr, almost certainly >8 Gyr. In Figure 7, we show example cooling curves, computed for thin and thick DA atmospheres and C/O WDs (likely irradiation is a negligible perturbation to the WDs surface temperature, given the measured spin-down luminosity of the pulsar). For the model closest to the best-fit mass of PSR J2222 we would infer that the true age is near 9 Gyr, with the possible range from 6–12 Gyr. The upper limit provided by the pulsar's characteristic spin-down age (34 Gyr after correction for the Shklovskii effect; Shklovskii 1970) is not constraining; the assumption that the pulsar's initial spin period is much shorter than the current spin period is clearly not valid. Instead, we take as our upper limit to the age that of the Milky Way's halo (11.4 ± 0.7 Gyr; Kalirai 2012) minus the ≈70 Myr required for the main-sequence lifetime of a ≈6 M_☉ progenitor (Koester & Reimers 1996; Williams et al. 2009), although this does not really exclude any models. Such an age would, however, imply a lower limit to the (re-)birth period of about 25 ms, assuming spin-down with a braking index n = 3 (magnetic dipole radiation). We note that the cooling models in Figure 7 may not be the only solution for this progenitor: changing the WD composition (likely it is below the transition to O/Ne/Mg WDs based on Nomoto 1984; Iben & Tutukov 1985, although binary evolution could change that; also see Lazarus et al. 2014) or atmosphere (helium, carbon, etc) could lead to different solutions, and to draw robust conclusions we need to explore a wider range of models with better observational constraints. There are also considerable complications and uncertainties in models for these temperatures: for instance, the models of Salaris et al. (2010; the BaSTI database) give rather different ages as T_eff never drops below 4000 K for 1.0 M_☉ models, even for ages of >14 Gyr, while Althaus et al. (2007) and Chabrier et al. (2000) do have ≈1.0 M_☉ models go below 4000 K (note that the models in Althaus et al. 2007 are primarily O/Ne rather than C/O). However, we believe the T_eff upper limits to be more robust, as they do tend to agree between different calculations.

Zoom In Zoom Out Reset image size

While extreme, the companion to PSR J2222 may not be especially unique. Similar ultra-cool WDs are presumably present in globular clusters and in the field even if they are often too faint to identify on their own. Individual ultra-cool WDs can be identified but only if very nearby, like the two objects in Kilic et al. (2012) at ≈30 pc. If we correct roughly for the different progenitor masses between the Kilic et al. (2012) systems and PSR J2222 (Kalirai et al. 2008) and use a Salpeter (1955) initial mass function, we would estimate ≈200 massive WDs of a similar age within 300 pc, which is of the same order as the luminosity function from Rowell (2013, also see Giammichele et al. 2012) extrapolated¹⁶ to M_bol > 19.

Instead, binary systems are the best way to identify cold WDs (e.g., Parsons et al. 2012), which is effectively the technique used here. However, even in binary systems where we know that a source is present, the systems will often be too distant for good constraints (i.e., PSR J1454−5846 in Figure 5; Jacoby et al. 2006). We still require a fortuitously nearby system for useful observations. The occurrence of a nearby massive WD like the companion to PSR J2222 is reasonably consistent with expectations based on the observed binary population: there are five pulsar binaries from the ATNF Pulsar Catalog¹⁷ (Manchester et al. 2005) within 300 pc, and the other four have low-mass He WD companions. This 1/5 ratio is similar to that for CO WD compared to He WD companions in the whole ATNF catalog (also see Tauris et al. 2012), and the pulsars' spin-down ages appear to have similar distributions for both companion types.

Finally, we can ask whether an NS is the most likely companion to an ultra-cool WD. Most binaries are assumed to have mass ratios near one (Pinsonneault & Stanek 2006, but see Sana et al. 2012), but a binary composed of two ultra-cool WDs would be just as hard to detect optically as a single object. If the companion were a lower-mass WD or a main-sequence star the binary could be visible, although it would require spectroscopic follow-up to identify the companion and in the absence of GAIA this has not been done for the majority of stars within a few hundred pc. So the situation of PSR J2222, with an NS companion, is reasonably plausible as the initial mass ratio would have been close to one and the chances of companion follow-up and identification after discovery of the pulsar are high.