A METAL-RICH LOW-GRAVITY COMPANION TO A MASSIVE MILLISECOND PULSAR

The number and variety of millisecond pulsars (MSPs) has expanded dramatically in recent years, due to the combination of significant radio surveys (e.g., Bailes et al. 2011; Keith et al. 2012; Boyles et al. 2013; Lynch et al. 2013) and the advent of the Fermi satellite (see Ray et al. 2012 for a recent review). These millisecond pulsars offer unprecedented laboratories to study particle acceleration (Abdo et al. 2009), binary evolution (Benvenuto et al. 2012), and the masses of both white dwarfs (Corongiu et al. 2012) and neutron stars (Demorest et al. 2010; Romani et al. 2012).

Quite a number of discoveries have in fact been the result of the confluence of these trends, picking out unidentified γ-ray sources as priority targets for ongoing radio surveys. This technique resulted in the discovery of the eclipsing MSP PSR J1816+4510 (Kaplan et al. 2012) as part of the ongoing Green Bank North Celestial Cap (GBNCC) survey (K. Stovall et al. 2013, in preparation). This source has a spin period of 3.2 ms, an orbital period of 8.66 hr, and shows both γ-ray pulsations and radio eclipses for ≈10% of the orbit.

While the radio and γ-ray properties resembled a black widow or redback system (Roberts 2011), with a non-degenerate companion being ablated by the pulsar's wind, the blue colors and lack of variability suggested instead that the companion to PSR J1816+4510 was more similar to a helium-core white dwarf (van Kerkwijk et al. 2005), even though the inferred radius seemed too large and the radio eclipses were unprecedented among those systems.

Here we present optical spectroscopy of the companion to PSR J1816+4510 (Section 3.1), taken with the Keck I, Palomar 5 m, and William Herschel telescopes. We confirm through radial velocity variations that the source identified by Kaplan et al. (2012) is in fact the companion (Section 3.2), and that, at first blush, it resembles a hot, low-gravity white dwarf (Section 3.1). However, detailed analysis shows that the spectrum has strong lines of helium and metals (calcium, silicon, magnesium, and possibly iron). We discuss the implications of the spectrum and radial velocity measurements in Section 4, along with possible origins for the unique features of this system and ways they can be exploited to learn more about neutron stars and their companions. We conclude in Section 5.

We observed the likely counterpart of PSR J1816+4510 with a number of instruments over 2012, as summarized in Table 1. We primarily report data from the blue sides of the spectrographs, although the red sides were analyzed as well. Flux calibrators were observed with the same setups (GD 153 for the WHT data, Feige 34 for the 2012 April P200 data, BD+33°2642 for the Keck data and the 2012 May P200 data). The slit width varied between 07 and 15 depending on the conditions. For all observations of PSR J1816+4510 except for the last two Keck observations, we did not have the slit at the parallactic angle as we instead included both the pulsar and a reference star on the slit to provide a velocity standard (for the 2012 June Keck observations the reference star was instead observed right after the pulsar). We used two different stars, as listed in Table 1 and shown in Figure 1.

Zoom In Zoom Out Reset image size

For all of the observations our reduction followed standard procedures, although we used MIDAS for the William Herschel Telescope (WHT) data and IRAF for the other data. We subtracted biases and overscans then flattened the data using dome flats, where the shape of the lamp was removed with a high-order polynomial. We then extracted the spectra of the pulsar and the reference star using optimal weighting, subtracting background as measured on both sides of the object spectrum. Wavelength calibration was done using arc spectra, and fluxes were calibrated using the spectrophotometric standards, correcting for differences in airmass using extinction curves appropriate for each site.

3.1. The Spectrum of PSR J1816+4510

To start, we focus on the highest-quality spectra which are the 2012 March Keck data. We show a set of blue and red-side spectra in Figure 2. One immediately notices the broad Balmer lines, consistent with our tentative identification of the optical counterpart as a low-mass white dwarf (Kaplan et al. 2012). However, it is also clear that Balmer lines are not the only ones present. We also see strong but narrower lines from He i, with at least 20 lines present between 3000 Å and 7000 Å, as well as narrow lines from Ca ii, Si ii, and Mg ii, and possibly Fe i (see Figure 4 for a more detailed view).

Zoom In Zoom Out Reset image size

We fit the data to white-dwarf model atmospheres to constrain the atmospheric parameters log₁₀(g) and T_eff as well as the He and metal abundances. While we obtain reasonable results, we note that the validity of our models has not really been tested for gravities log₁₀(g) < 5.50 where our best fits are. Instead those gravities are more typical of main-sequence/sdB stars. With this caveat, and acknowledging the inhomogeneous data set, our fitting procedure was as follows.

Our model grid covered gravities log₁₀(g) = 4–6 (in steps of 0.25 dex) and temperatures T_eff = 14, 000–20, 000 K (in steps of 250 K). We also included helium, with relative abundances log₁₀He/H = −1 to +1 by number. For each spectrum we interpolated the model spectrum to the observed wavelength grid for a range of trial velocities, convolving the model with the pixel width and a Gaussian determined by the seeing, and truncating by the slit; we explored also convolving by the expected radial velocity change during each observation (see below), but that did not change the results noticeably. We normalized the continuum using some line-free regions and a third-order polynomial fit to each spectrum. The χ²(log₁₀(g), T_eff, v) was found by summing over the 3750–4455 Å region, excluding metal lines, and done separately for each value of He/H.

After optimizing for velocity (see below), the fit was covariant among remaining parameters, with the helium abundance altering not just the line depths but also the continuum shape. We find reasonable fits for log₁₀He/H = −0.5 to +0.5: going to log₁₀He/H = −1 (i.e., solar) increases χ² by 8670 compared to log₁₀He/H = 0, for 9372 degrees-of-freedom. For each value of the helium abundance, the best-fit locus in (log₁₀(g), T_eff) space was also covariant, with a positive correlation between the fitted parameters (Figure 3). Overall more helium resulted in a lower effective temperature and surface gravity. Our best-fit was log₁₀He/H = 0.0 ± 0.5, log₁₀(g) = 4.9 ± 0.3, and T_eff = 16, 000 ± 500 K; we believe these uncertainties to be conservative. We see no evidence for variations in effective temperature with orbital phase, down to a limit of about ±500 K (limited by the non-uniformity of the data and the uneven phase sampling).

Zoom In Zoom Out Reset image size

With the background model determined we examined the metal lines. We fitted for the abundances of the individual species listed in Table 2 (helium was fit along with the spectrum, as discussed above). Other elements were present in our models but they did not produce significant lines in the 3700–7000 Å range. While we could see lines from sodium in the spectrum, they had different velocities from the rest and therefore are presumably mostly interstellar; hence we quote only an upper limit. For iron, no individual line is significantly detected, but the large number of lines makes the presence of iron very likely (Figure 4). There may be other species present at low levels, and we could also determine a number of upper limits, but given the model uncertainties and our difficulty in interpreting just the current data we leave the analysis here until the models and data have improved. Overall, we find that the metals have abundances consistent with roughly 10 times solar abundance relative to hydrogen, or roughly solar relative to helium.

Zoom In Zoom Out Reset image size

Table 2. Elemental Abundances

Species	log10(N/NH)	log10(N/NH)solar	[X/H]
He	0.0 ± 0.5	−1.0	1.0 ± 0.5
Naa	<−3.0	−5.7	<2.7
Mg	−3.8 ± 0.2	−4.4	0.6 ± 0.2
Si	−3.8 ± 0.2	−4.5	0.7 ± 0.2
Ca	−4.5 ± 0.3	−5.6	1.1 ± 0.3
Feb	−2.8 ± 0.3	−4.5	1.7 ± 0.3

Notes. Solar abundances are from Anders & Grevesse (1989). ^aMost of the Na seen in the spectrum was interstellar, with significant wavelength shifts compared to the other species. ^bDetection of Fe consisted of numerous lines near the noise threshold. Individual lines were not significant, but the ensemble was.

Download table as:  ASCIITypeset image

3.2. The Radial Velocity of PSR J1816+4510

We determined radial velocities of the companion to determine its orbit. Unfortunately, in the blue part of the spectrum where the majority of the spectral lines are there are not many sky lines that can be used as absolute velocity references. This deficiency is made more severe by the two-armed nature of the spectrographs used, where the red- and blue-side wavelength calibrations are separate. Instead, we oriented the slit to include both the pulsar's companion and another reference star. As discussed above, the observations used two separate reference stars (Table 1).

Using the models of Munari et al. (2005), we fit for the spectral type of our reference stars over the wavelength range 3700 Å–5400 Å, allowing for an overall scaling to account for errors in flux calibration, slit losses, etc. Reference star 1 (S1) appears to be of spectral type mid- to late-F (effective temperature 6250 K), while reference star 2 (S2) is a bit cooler, roughly G2 (5750 K), as seen in Figure 1. This is reasonably consistent with the observed colors (from Two Micron All Sky Survey (Skrutskie et al. 2006) and other archival data), where S2 is slightly redder than S1.

To determine the radial velocities of the pulsar's companion and the reference star, we used appropriate model atmospheres, minimizing the χ² between the model and the data for a range of trial velocities. At each velocity we convolved the model as described above, and allowed for a low-order polynomial difference between the model and the data. We restricted the wavelength range to 3700 Å–5400 Å (the maximum wavelength was 4450 Å for the WHT data) so as to focus on a well-calibrated part of the data above the highest-order Balmer lines and Balmer break but below the dichroic cutoffs.

For the pulsar, the fits were good, with reduced χ² typically near 1 for the WHT and P200 data. The higher-quality Keck data had reduced χ² of ≈4, so we scaled the velocity uncertainties such that the reduced χ² was 1. The fits were also good for the reference stars, although not quite as good as for the pulsar. Instead χ² was typically about 1.5 for the P200 and WHT data and up to 6 for the Keck data. This may be from a slight mismatch in resolution (although we explored changing our resolution), but more likely is from having an imperfect model grid that is coarse and that may not reflect the true metallicity or rotation. In any case, the best-fit velocities did not depend much on the exact model or the choices of the other parameters.

The use of two reference stars complicates our analysis slightly, as each star can of course have its own peculiar velocity. We therefore obtained two spectra where S1 and S2 were on the slit together without the pulsar. These exposures were right before the observations of the pulsar on 2012 August 18 and used the same configuration and reduction. We found a consistent radial velocity offset of −54 ± 1 km s⁻¹ for S2 relative to S1. In what follows, we apply this offset to determine the velocities of the pulsar relative to S1.

We show the best-fit radial velocities in Figure 5. The systemic velocity is relative to S1 (in no case were the finite durations of the exposures significant compared to the orbital period, with corrections of <1 km s⁻¹). We have included a ±12 km s⁻¹ systematic velocity uncertainty in the data to account for limitations in our reduction and differences between the instruments. This will also account for mis-alignment of the slit intended to cover two objects (an offset of approximately one half of the slit width corresponds to a velocity shift of approximately 10 km s⁻¹). With that, we were able to fit a circular orbit using the phase established by the radio ephemeris (Stovall et al.), which has since been confirmed in γ-rays (Kaplan et al. 2012) and is more accurate than an optical-only ephemeris. With all of the data we find a radial velocity amplitude K_c = 343 ± 7 km s⁻¹ for the companion, with χ² = 14.5 for 14 degrees of freedom and a systemic velocity of −99 ± 8 km s⁻¹ relative to S1. If we only use the P200 data including S2—a somewhat more homogeneous data set—we find a consistent K_c = 339 ± 11 km s⁻¹ (χ² = 3.5 for 4 degrees of freedom). We can also reference all of the velocities to the 5577 Å sky line. Here we might worry that the line is at the extreme edge of our spectra where the wavelength solutions are not as reliable, but the result is consistent: K_c = 347 ± 11 km s⁻¹ and a systemic velocity of −108 ± 12 km s⁻¹ (χ² = 20.9 for 12 degrees of freedom).

Zoom In Zoom Out Reset image size

Combining the best-fit radial velocity amplitude of the companion, K_c = 343 ± 7 km s⁻¹, with the measured projected semi-major axis of the pulsar x_PSR = 0.595402 ± 0.000002 s, we can infer the mass ratio:

$$q\equiv \frac{M_{\rm PSR}}{M_{\rm c}}=\frac{K_{\rm c}P_{\rm B}}{2\pi x_{\rm PSR} c}=9.54\pm 0.21.$$

The minimum masses for the pulsar and companion are then: M_psrsin ³i = 1.84 ± 0.11 M_☉ and M_csin ³i = 0.193 ± 0.012 M_☉, where i is the inclination of the orbit. Note that, unlike black widow systems that are strongly affected by heating such that the photocenter and center of mass are not the same (van Kerkwijk et al. 2011; Romani et al. 2012), the minimal modulation seen for PSR J1816+4510 (Kaplan et al. 2012; D. L. Kaplan et al. 2013, in preparation), suggests that we are measuring the true center-of-mass motion.

Combining this mass determination with our spectroscopically-determined gravity, we can determine the radius of the companion and the distance to the system. In what follows, we give the values for the minimum mass, equivalent to that for an edge-on orbit, which is not unreasonable given the eclipses (although since they are radio eclipses from ionized material, constraints on the inclination are weak). The minimum radius is 0.26  ±  0.08 R_☉, implying an average density of 15 g cm⁻³. To get the distance, we use the bolometric corrections appropriate for this effective temperature from Tremblay et al. (2011, using the 0.2 M_☉ models),⁷ with M_bol − R = −1.52 ± 0.09 mag; while this is for a larger gravity (log₁₀(g) = 6.6) than what we measure, the dependence of the bolometric correction on gravity is quite small. We assume an extinction A_V = 0.22 mag, which comes from fitting our photometry to the same synthetic photometry as we used for the bolometric corrections.⁸ Overall we find d = 4.5 ± 1.7 kpc using the R-band photometry from Kaplan et al. (2012), where the uncertainty is dominated by the uncertainty on the measured gravity and temperature ($d\propto T_{\rm eff}^{0.75}10^{-\log _{10}(g)/2}$), and we get similar results integrating our model atmosphere directly. The distance is almost a factor of two larger than the dispersion-measure distance discussed in Kaplan et al. (2012), suggesting that, as is found for pulsars out of the Galactic plane generally (Roberts 2011), the dispersion-measure distance is an underestimate. The nominal inferred radius is correspondingly larger, indicating that this source is even further from equilibrium than previously thought with a thermal timescale of ≈1 Myr (the radii of 0.2 M_☉ WDs are ≲ 0.03 R_☉), although it is still well within the Roche-lobe radius of 0.56 R_☉. The distance also implies a larger γ-ray luminosity, L_γ = 5 × 10³⁴ erg s⁻¹, roughly equal to the total spin-down luminosity, $\dot{E}=5 \times 10^{34}\,{\rm erg\,s}^{-1}$ (ignoring possible corrections for secular acceleration and assuming the standard moment of inertia I = 10⁴⁵ g cm²), so $L_\gamma / \dot{E}\approx 1.0_{-0.7}^{+0.9}$. Note that the $L_\gamma < \dot{E}$ requirement is somewhat rough, with observed values of $L_\gamma /\dot{E}$ (based on poorly-known dispersion measure distances, and assuming beaming into 4π sr) ranging from 10% to >200% (Ransom et al. 2011); given the uncertainty in our photometric distance we see no problems.

For the nominal $\dot{E}$, we expect ≈1.1 × 10¹¹ erg s⁻¹ cm⁻² of incident flux at the surface of the companion. This compares with σT⁴_eff = 3.7 × 10¹² erg s⁻¹ cm⁻² coming from the companion itself, so we would expect variations of at most ±100 K coming from the day and night sides of the companion (assuming an albedo of 100%). Our spectra cannot yet constrain such a variation, but this should be possible in the future with high signal-to-noise spectroscopy or photometry.

Our very low gravity brings up a question of nomenclature. Initially we called the companion a white dwarf, and indeed we fit the spectra with white dwarf model atmospheres. But this gravity is really more consistent with a sdB star or a proto-WD, such as GALEX J171708.5+675712 (Vennes et al. 2011) or HD 188112 (Heber et al. 2003). As shown in Figure 6, PSR J1816+4510 is even more extreme in the gravity axis than those sources. It almost certainly is not a standard sdB star, as these typically have effective temperatures >25,000 K and masses of ≈0.5 M_☉ (able to sustain core He fusion) compared to 16,000 K and 0.2 M_☉ here.⁹ It may still be contracting and cooling from a recent hydrogen shell flash, as discussed in Kaplan et al. (2012). We note that the approximate mass of 0.2 M_☉ is at the limit where He-core WDs are expected to undergo shell flashes instead of steady burning (Panei et al. 2007): if we could determine a more precise mass and establish whether or not unstable flashes occurred, it will help our understanding of the evolution of low-mass WDs.

Zoom In Zoom Out Reset image size

Another question is the origin of the metals.¹⁰ While planetesimal accretion has been suggested for other low-gravity WDs (Figure 6; Gänsicke et al. 2012), this seems very unlikely for PSR J1816+4510, where the relativistic wind of the pulsar and the ionized wind of the companion will both suppress accretion. The solar Ca/He ratio also supports a primordial origin for the metals (cf. Dufour et al. 2010). At the other mass extreme, comparison with sdB stars also shows differences. For instance, Edelmann et al. (2003) show that helium abundance increases with effective temperature, becoming close to solar only for T_eff > 30,000 K. Geier (2013) finds no correlation of metal abundances with T_eff or the helium abundance, although they are generally sub-solar. The ionized wind that likely causes the radio eclipses may be part of the answer here, as it might alter diffuse equilibrium (much as inferred for sdB stars winds; Unglaub & Bues 1998) and hence the onset of shell burning. It also may remove enough angular momentum to create ultra-compact binaries such as PSR J1719−1438 and its planetary-mass companion (van Haaften et al. 2012). Alternatively, if a shell flash occurred recently, this might not only explain the low surface gravity but also the abundances, as it might have mixed up the interior and the present low gravity would delay sedimentation compared to an ordinary white dwarf.

Another interesting comparison is the black widow system PSR J1311−3430 (Romani et al. 2012). Its companion has a lower gravity, a much lower mass, and is likely less degenerate than the companion of PSR J1816+4510, yet has a comparable maximum temperature and a spectrum dominated by lines of He i, with small amounts of metals (Ca ii, Mg ii) visible and no detectable hydrogen. Could PSR J1816+4510 evolve into a system like PSR J1311−3430? Current suggestions for the origin of PSR J1311−3430 involve substantial radiative stripping (Benvenuto et al. 2012), rather than unstable hydrogen burning, so perhaps not. Indeed, there are indications for a dense wind from PSR J1311−3430 that appears driven by irradiation; irradiation appears less important for PSR J1816+4510, but detailed calculations are necessary to determine this quantitatively.

The mass of the neutron star is considerably above the canonical value of 1.4 M_☉. This is consistent with the accretion necessary to have spun it up to its current period (Tauris et al. 2012) and with the masses of a number of other pulsars with low-mass He white dwarf companions (van Kerkwijk et al. 2005; Corongiu et al. 2012), almost to the level of black-widow systems (van Kerkwijk et al. 2011; Romani et al. 2012). While our current measurement is not as extreme as some (Demorest et al. 2010), the brightness of the companion and its narrow spectral lines mean that we can potentially obtain a precision of <1% on K_c and hence on the mass ratio without uncertainty about the photocenter (cf. van Kerkwijk et al. 2011). At this point our uncertainty in the masses will be entirely dominated by uncertainties in the inclination. A Shapiro delay measurement is likely impossible because of the eclipses, so we may have to rely instead on combining spectroscopy with modeling the photometric light curve (Reynolds et al. 2007; Romani et al. 2012); we have begun detailed time-resolved photometry to search for irradiation or ellipsoidal variations (expected to be roughly 2%). Given the masses we infer, we can expect the system to merge due to gravitational radiation in about 1 Gyr (Shapiro & Teukolsky 1983); the actual merger time might be shorter if energy is dissipated through tidal interactions, and the bloated nature of the companion could make this more significant than otherwise inferred (Burkart et al. 2012).

A related possibility is that we might be able to determine the rotation of the companion, which would help understand the history of mass transfer in the system (e.g., Pablo et al. 2012). If rotating synchronously with the orbit, we would expect the equatorial velocity to be ≈35 km s⁻¹. This high velocity combined with the low surface gravity and the presence of narrow He lines means that we might separate the pressure broadening responsible for the Balmer lines from rotation. Initial estimates suggest this might be possible, but again it requires considerably better data than we have right now along with improved atmospheric modeling.

We have described an initial spectroscopic study of the companion to the millisecond pulsar PSR J1816+4510. Our findings largely agree with the initial speculation from Kaplan et al. (2012): that the companion is hot and large, consistent with neither black widow companions nor with standard white dwarfs. However, we now find that it is even more distinct, in showing strong lines of helium and metals in its spectrum, together with the abnormally low gravity and ionized-gas eclipses.

We are left with a number of questions. (1) What is the lifetime of such an object? It cannot be too small compared to typical MSP ages, since otherwise it would be too unlikely to find it among the ∼100 known MSP systems. (2) Why are there strong metals and helium? Are they related? Is the cause of metals similar to that of normal low-mass WDs, or is it unique to this system? (3) What are the detailed metal abundances? Are they all the same relative to solar? (4) Why is it eclipsing? A temperature of 16,000 K is generally thought to be too low for a wind. Is that related to the pulsar (is it ablating)? Or is it a more general property of the companion? Finally, (5) does the presence of the pulsar significantly affect (through its relativistic wind) the evolution of the companion? These questions have important consequences for our understanding of white dwarf and neutron star evolution, as well as for the study of interacting binary systems.

Partially based on observations made with the William Herschel Telescope operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and NASA; the Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. We thank David Levitan, Assaf Horesh, and Yi Cao for assistance with observations; research at Caltech was funded by a grant from the NSF. We thank an anonymous referee for helpful suggestions.

Facilities: Keck:I (LRIS) - KECK I Telescope, Hale (Double Spectrograph) - Palomar Observatory's 5.1m Hale Telescope, ING:Herschel (ISIS) - Isaac Newton Group's William Herschel Telescope