Peering into the Dark Side: Magnesium Lines Establish a Massive Neutron Star in PSR J2215+5135

We obtained phase-resolved photometric observations of PSR J2215+5135 with the IAC-80 and William Herschel (WHT) telescopes at the Teide and Roque de los Muchachos observatories, respectively, on four different nights (see Table 1). The IAC-80 images were taken on 2014 August 2–3 with the Teide observatory light improved camera (CAMELOT; 0.30 arcsec pixel⁻¹) using the SDSS filters g' and r', an exposure time of either 420 s or 600 s and a binning of 2 × 2 pixels in order to optimize the signal-to-noise ratio (S/N). The WHT images were taken on 2014 August 11 (contemporaneous with our WHT-ISIS spectra) and September 1, with the auxiliary-port camera (ACAM; 0.25 arcsec pixel⁻¹ in imaging mode) using the SDSS filters $g^{\prime} ,r^{\prime}$ and $i^{\prime}$, an exposure time of 60 s and a window of 501 × 501 pixels around the source in order to reduce readout time (readout and filter change resulted in deadtime of 6–9 s per exposure).

Table 1.  Summary of Optical Observations of PSR J2215+5135

Telescope	Instrument	Banda	Date	Time	Exposures	Orbital	Airmass	Seeing
(diameter)	(configuration)	(λ, Å)	(evening)	(UT)	(nr. × duration)	Phase		('')
Photometry
IAC-80	CAMELOT	$g^{\prime} ,r^{\prime}$	2014 Aug 02	22:53-04:01	4 × 420 s	0.4–1.6	1.47–1.09	1.1–2.5
(80 cm)	(bin 2 × 2)	(4639, 6122)			+26 × 600 s
IAC-80	CAMELOT	$g^{\prime} ,r^{\prime}$	2014 Aug 03	23:32-03:52	30 × 420 s	0.3–1.3	1.32–1.09	0.4–0.7
(80 cm)	(bin 2 × 2)	(4639, 6122)
WHT	ACAM	$g^{\prime} ,r^{\prime} ,i^{\prime}$	2014 Aug 11	22:15-00:03	95 × 60 s	0.4–0.8	1.48–1.18	0.7–1.6
(4.2 m)	(win. 501 × 501)	(4639, 6122, 7439)
WHT	ACAM	$g^{\prime} ,r^{\prime} ,i^{\prime}$	2014 Sep 01	23:17-02:06	149 × 60 s	0.4–1.0	1.13–1.09	0.8–1.4
(4.2 m)	(win. 501 × 501)	(4639, 6122, 7439)
Spectroscopy
WHT	ISIS	B:3700–5300	2014 Aug 11	00:23-05:19	17 × 900 s	0.0–1.1	1.17–1.33	0.7–1.6
(4.2 m)	(R600B+R)	R:5500–7200
GTC	OSIRIS	B:4000–5700	2014 Nov 14	19:47-21:06	5 × 935 s	0.6–0.8	1.09–1.11	0.6–1.1
(10.4 m)	(R2000B)
GTC	OSIRIS	B:4000–5700	2014 Nov 15	19:43–23:58	16 × 935 s	0.3–1.3	1.09–1.55	0.7–1.3
(10.4 m)	(R2000B)

Note.

^aEffective wavelength of the photometric filters or approximate wavelength range covered by the spectra, in Angstroms.

Download table as:  ASCIITypeset image

All images were debiased and flat-fielded using standard iraf routines. We then performed aperture photometry using the ultracam pipeline, a variable source extraction radius (set to 1.5–1.7 times the seeing) and a nearby stable nonsaturated reference star. The resulting differential magnitudes are relative to a nearby star in the field and thus insensitive to thin clouds or moderate atmospheric variability (yet the observing conditions were generally good). The absolute flux (apparent magnitude) calibration was done using a nearby AAVSO-APASS star, with uncertainties of 0.05 mag (r') and 0.1 mag (g', i'), and checked against other nearby stars from the USNO-B1 catalog. We also compared these reference star magnitudes with those given by the PANSTARRS catalog and found only a significant difference in the r' band, with a shift of +0.16 mag with respect to the APASS values that we use. We verified that this has no impact on any of the results reported in this work; in particular, the parameters reported in Section 4 are all consistent within the errors when using the PANSTARRS instead of the APASS calibration. Figure 1 (top) shows the J2215 LCs folded at the orbital period.

Zoom In Zoom Out Reset image size

We observed J2215 with the WHT and GTC telescopes on 2014 August 11 and November 14–15, respectively, in order to obtain medium resolution spectra of the companion star and measure its spectral type and velocity along the orbit. For the WHT-ISIS spectra, we chose the R600B (blue arm) and R600R (red arm) gratings centered at 4500 Å and 6400 Å, respectively. The slit width was set to 1'', resulting into a resolution of 105–130 km s⁻¹ (R∼2600) and 65–80 km s⁻¹ (R ∼ 4000) for the blue and red arms, respectively. At the GTC, we used OSIRIS in its long-slit spectroscopy mode, with the R2000B VPH gratings and a slit width of 1'', resulting in a resolution of 145–160 km s⁻¹ (R∼2000). With this campaign we obtained 17 WHT-ISIS and 21 GTC-OSIRIS spectra covering the full 4.14 hr orbit with some redundancy and with exposure times of 900 s and 935 s, respectively (i.e., exposing each spectrum for about 6% of the J2215 orbit; see Table 1).

After applying bias and flat corrections to the trimmed images within iraf, we extracted the spectra and subtracted sky background using the optimal extraction method within starlink/pamela to account for the significant tilt (Marsh 1989). The WHT-ISIS spectra were calibrated in wavelength using interspersed arc spectra (CuNe, CuAr) extracted from the same source extraction regions, taken once every two source spectra. A set of well identified arc lines were satisfactorily fitted with a fourth-order polynomial to produce the wavelength scale (47 and 31 lines in the blue and red arms, resulting in an rms amplitude of residuals of 0.05 Å and 0.02 Å, respectively). We adjusted the same polynomial function to all arc spectra and interpolated in time between adjacent arcs to calibrate the science spectra, thereby accounting for the significant (∼1 Å) drift due to instrument flexure. The GTC-OSIRIS wavelength calibration was done using one set of arcs taken on the second night, fitting the pixel-wavelength relation with a fourth-order polynomial (19 lines giving residuals with an rms amplitude of 0.06 Å). The resulting wavelength calibration was checked and refined using the O i sky line at 5577 Å, which allowed us to correct for residual (≲10 km s⁻¹) shifts in the wavelength solution.

In order to measure the radial velocity and temperature of the companion star in J2215 throughout the orbit, we applied within molly the cross-correlation and optimal subtraction techniques, respectively. As both techniques require comparison spectra (or templates hereafter), we built a library of main-sequence stellar templates with spectral types between O4 and M0 (see the Appendix). The continuum level of all 38 source and 33 template spectra was normalized with a spline fit, and subtracted.

After binning to a same heliocentric velocity scale, excluding telluric lines and broadening the template spectra to the source spectral resolution (Section 2.2), each source and template spectra were cross-correlated to find their relative velocity, allowing shifts between −700 and +700 km s⁻¹. The resulting radial velocity curves (RVCs) from the GTC and WHT spectra were fitted with a sine function V = G +K sin[2π(t – T₀)/P_orb], where V and G are the radial and systemic velocities, K is the semi-amplitude of the RVC, t is the middle time of each spectrum, T₀ is the time of inferior conjunction of the secondary that defines⁴ orbital phase Φ_orb = 0 and P_orb is the orbital period. Having checked that the best-fit period is consistent with (but less precise than) the orbital period from pulsar timing (Abdo et al. 2013), this parameter was subsequently fixed at the pulsar timing value P_orb = 0.1725021049[8] days.

Using the orbital parameters above, we corrected for the systemic velocity and orbital motion and shifted all source spectra to the reference frame of each template. In order to increase the S/N, we averaged four to six source spectra around Φ_orb = 0 ± 0.15 and Φ_orb = 0.5 ± 0.125. We then performed an optimal subtraction using the full GTC spectral range, i.e., subtracted the templates scaled by a factor f_veil from the source averaged spectra, adjusting f_veil to minimize the residual scatter. This is a quantitative way of matching the observed absorption lines from J2215 to a set of templates with known spectral types and temperatures (Marsh et al. 1994).

Because J2215 becomes very faint around Φ_orb = 0 (r ≃ 20 mag; Figure 1, where the cold face of the companion star dominates), the corresponding ISIS spectra have low S/Ns. The tightest constraints on the spectral type and radial velocities come from the higher S/N GTC spectra, so we focus the rest of our analysis on those. Only 1 out of the 21 GTC spectra could not be included in the analysis due to the very low (<100 at peak) number of counts collected.

Thanks to the large GTC collecting area, we were able to measure radial velocities using (i) the full spectral range (4000–5300 Å), (ii) the hydrogen Balmer lines (three 50 Å wide windows centered on Hβ, Hγ and Hδ), and (iii) the Mg i triplet lines present in the J2215 spectra (5152–5199 Å). The corresponding spectral ranges and two representative GTC-OSIRIS spectra are shown in Figure 3. As the absorption spectra of early/hot and late/cold stars are dominated by different sets of lines (see, e.g., the Appendix), the Balmer and Mg i radial velocities allow us to track different parts of the irradiated companion throughout its orbit around the pulsar.

The cross-correlation of two broad lines may yield ambiguous results if the profiles are not exactly the same. The matching wings can give two maxima in the cross-correlation function, at velocities that differ from that measured using the line core ("double peaked" cross-correlation functions). On close inspection of J2215's RVCs, we found that this introduces strong artificial deviations from a sinusoidal function around phase 0.5, only when using Balmer lines and templates of spectral type earlier than F (i.e., when both source and template spectra are dominated by broad lines). For this reason, we include only spectral types later than F0 in our results for the Balmer-line RVCs (Section 3.2).

The optical flux from J2215 varies smoothly along the orbit (no flares are detected in the 60 s cadence data), with one clear maximum and minimum per orbital cycle and a peak-to-peak amplitude of almost two magnitudes (g', B bands). We find that the orbital LCs (Figure 1) are stable over timescales of years, comparing our 2014 observations with the 2010–2011 data presented by Schroeder & Halpern (2014; converting their BVR magnitudes into the SDSS $g^{\prime} r^{\prime} i^{\prime}$ system following Jordi et al. 2006). The ∼0.4 mag change in the [$g^{\prime} -r^{\prime}$] color reveals hotter emission at maximum light (Φ_orb = 0.5). The RVC of the companion shows a large amplitude (K ∼ 400 km s⁻¹) and changes sign near the maximum and minimum light (Φ_orb = 0.5 and 0, respectively), as can be seen already from the WHT spectra (Figure 1) and as first reported by Romani et al. (2015). These properties are indicative of a companion star that is strongly irradiated by the pulsar throughout the compact 4.14 hr orbit.

The J2215 spectra around Φ_orb = 0.5 (Figures 2 and 3) show strong Balmer lines (H1 through H9) consistent with an A5 star, as well as numerous yet weaker metallic (Ca/Fe/Mg) lines. Around Φ_orb = 0, when the companion star presents its cold face to the observer, Balmer lines are much weaker and narrower, while Mg i triplet lines are stronger. Hence the equivalent widths (EW) of Balmer and Mg i lines are anticorrelated along the orbit, as shown in Figure 4. The optimal subtraction analysis (Section 2.3) gives a clean measurement of the temperature and spectral type of the companion star, independent of the measured colors (which may be contaminated by nonstellar light). We find drastic changes between the irradiated and cold sides of the companion star. This is shown qualitatively in Figures 3 and 2, where the J2215 GTC and WHT spectra are compared to A5 and G5 templates degraded to the same resolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows our quantitative results: the reduced chi-squared resulting from the optimal subtraction method (Section 2.3) for templates with a broad range of spectral types (O–M). We thereby measure a spectral type A5±2 for the brightest spectra (i.e., A3–A7 at Φ_orb = 0.5) and G5±5 at the faintest end (G0–K0 at Φ_orb = 0; Romani et al. 2015 report similar yet slightly earlier spectral types of A2 and G0 around phases 0.5 and 0, respectively, but no errors are quoted). These correspond to effective temperatures for the cold ("night") and hot ("day") sides of T_N = 5660${}_{-380}^{+260}$ K and T_D = 8080${}_{-280}^{+470}$ K, respectively (Pecaut & Mamajek 2013). We thus find a drastic temperature contrast between opposite sides of the companion star, where the hot/day side is about 2400 K or 40% hotter than it would be without irradiation. The best matching scale factors are f_veil ≃ 0.8 at both superior and inferior conjunction, suggesting a contribution from nonstellar light (veiling) of about 20% in this GTC-OSIRIS 4000–5300 Å band (which corresponds approximately to filter g').

Zoom In Zoom Out Reset image size

We find that the radial velocities and K values depend systematically on the set of lines or spectral range used to measure them. Namely, as we show in Figure 6, the semi-amplitude of the Mg i line RVC (K_Mg, red circles) is always ∼10% higher than the semi-amplitude of the Balmer-line RVC (K_Balmer, blue squares). Using the same G5 template yields K_Balmer = 382.8 ± 4.7 km s⁻¹ and K_Mg = 420.2 ± 6.2 km s⁻¹.

Zoom In Zoom Out Reset image size

In Figure 7, we present two extreme cases: RVCs of J2215 calculated using Balmer lines (blue symbols) and Mg i-triplet lines (red symbols; see Section 2.3 for the exact wavelength ranges). In both cases the radial velocities were measured by cross-correlating the J2215 spectra with the same G5 template. The Mg i triplet lines yield higher radial velocities than the Balmer lines (in absolute value) around orbital phases 0.25 and 0.75, when the companion star is viewed sideways (in "quadrature"). Because the velocities at quadrature are extreme, this has an important impact on the measured K values and the corresponding neutron star mass (Section 5.1).

Zoom In Zoom Out Reset image size

Most RVCs are reasonably well fitted with a pure sinusoidal function (reduced chi-squared ≲2). There are, however, deviations noticeable in the sine fit residuals (Figure 7), especially in the Balmer-line RVCs (spectral types earlier than F1 were not used in the Balmer-line cross-correlation analysis, see Section 2.3). The fitted K values were verified in a model independent way by measuring the peak to peak semi-amplitude of the RVC: the two maximum and minimum radial velocities were averaged, subtracted and divided by two. The results were always consistent with the K values presented in Figure 6.

Using the full ensemble of lines, on the other hand, yields intermediate values of the RVC semi-amplitude (green triangles in Figure 6; K_All = 398.8 ± 2.3 km s⁻¹ for a G5 template). We also find a clear systematic dependence with the template's spectral type when using the full spectral range: a monotonic increase of K_All toward later (cooler) spectral types. We conclude that, when measuring radial velocities in strongly irradiated systems such as J2215, the spectral range and reference spectra must be chosen carefully in order to measure the pulsar mass accurately (see Section 5 for a further discussion).

The systemic velocities that we find are all in the 40–70 km s⁻¹ range, showing only small changes with spectral range or template spectral type. Averaging the results for spectral types G0–K0, we find G = 49.0 ± 2.5 (stat) ± 8.0 (syst) km s⁻¹ (where the statistical and systematic errors correspond to the standard deviation of the full range and of all three spectral ranges, respectively). From the same sine fits, we find an epoch of zero phase (companion at inferior conjunction) of T0 = 56976.9501 ± 0.0003 (stat) ± 0.0008 (syst) MJD (TDB), which we use together with the radio-pulsar timing P_orb (Section 2.3) in order to compute orbital phases.

xrbcurve includes a nearly Roche-lobe-filling secondary star heated by high-energy photons from the compact object and an accretion disk (not included in this case since no disk emission lines are detected). The binary system's geometry is determined by the orbital inclination i, the mass ratio q = M₂/M₁ (where M₁ and M₂ are the masses of the neutron star and secondary star, respectively), and the Roche-lobe filling factor of the secondary star, f. The orbital period P_orb, the radial velocity semi-amplitude of the secondary star's barycenter K₂ and the distance to the source in kpc (D) set the scale of the system. The light arising from the secondary star depends on its mean, unperturbed effective temperature (T₂), and the gravity darkening exponent β.

The additional light due to irradiation is given by the irradiation efficiency $\eta \equiv$ ${L}_{\mathrm{irr}}/\dot{E}$, which we define in accordance with previous work as the ratio of the heating luminosity (L_irr, assuming isotropic emission by an irradiating point source at the compact object) to the spin-down luminosity of the pulsar ($\dot{E}$ = 5.29 × 10³⁴ erg s⁻¹, Breton et al. 2013). We calculate the resulting increase in the local effective temperature due to the irradiating external source assuming that all irradiating flux is thermalized (Shahbaz et al. 2003). We use nextgen model-atmosphere fluxes (Hauschildt et al. 1999) to determine the intensity distribution on the secondary star and a quadratic limb-darkening law with coefficients taken from Claret (2000), to correct the intensity. Based on the observed mid-G spectral type for the secondary star (Section 3), we fix β at 0.10 (Lucy 1967).

To model the RVCs, we assume that the whole secondary star contributes to both Balmer and Mg i-triplet lines, including the inner/irradiated face. This is based on the corresponding EWs, which follow those expected from an A5–G5 star (Figure 4). We set the Balmer and Mg i-triplet absorption line strengths according to the effective temperature of each surface element of the star. Using our template spectra, we determine the EW versus temperature relationship for the same exact wavelength ranges used in measuring the RVCs (the Appendix; see Figure 13). The line strengths for each surface element on the star are then calculated using its temperature and the EW-temperature relation. Finally, radial velocities are calculated from the model line profile, averaged among all surface elements visible at each orbital phase.

In determining the binary parameters, we use a Markov chain Monte Carlo (MCMC) method convolved with a differential evolution fitting algorithm (see Vrugt 2016; Shahbaz et al. 2017). We use a Bayesian framework to determine our binary model parameters (Gregory 2005). We include the projected semimajor axis of the pulsar orbit measured from radio timing observations (x₁ = a₁ sin i = 0.46814[1] lt-s; Abdo et al. 2013) as an independent constraint on q and K₂ (q = K₂ P_orb / 2πx₁c; where c is the speed of light) using a Gaussian priori. Our MCMC fitting makes use of flat prior probability distributions for the rest of the model parameters. We use 20 individual chains to explore the parameter space and 40,000 iterations per chain. We reject the first 500 iterations and only include every tenth point.

We use a reddening of $E(B-V)$ = 0.38 mag, which we calculate from the X-ray absorbing column density toward J2215, ${N}_{{\rm{H}}}$ = 2.1 × 10²¹ cm⁻² (Gentile et al. 2014) using the conversion from Predehl & Schmitt (1995). Because N_H is in turn estimated from the measured pulsar dispersion measure (He et al. 2013), we verified the accuracy of the reddening toward J2215 in different ways. The total Galactic N_H in the direction of J2215 is 40% higher (Kalberla et al. 2005), although this may include additional absorbing material in the line of sight. Our value of $E(B-V)$ = 0.38 mag is consistent with that measured from IR dust maps in the same direction (0.35±0.02 mag, Schlegel et al. 1998). In order to quantify the possible impact on the measured pulsar mass, we repeated the LC and RV fits with $E(B-V)$ left as a free parameter, with a flat priori. Reassuringly, this yields the same value M₁≃2.3 M_⊙, and a value of $E(B-V)$ fully consistent with (less than 1σ from) the one we use.

We also compute the line-of-sight temperature at orbital phases 0.0 and 0.5, which represent the cold/night (T_N) and hot/day (T_D) temperatures of the secondary star. Using the spectral type measurements explained above, we impose temperature constraints on T_N (5280–5920 K, corresponding to a spectral type G5±5) and T_D (7800–8550 K, corresponding to a spectral type A5±2). Comparing with previous models of J2215 (Table 2), we conclude that our quantitative independent constraints on T_N and T_D are critical in order to find a good solution. From our measured f_veil ≃ 0.8 (Section 3.1), it is clear that there is an extra light component that veils the observed light from the secondary star. To allow for this wavelength-dependent veiling, we include an extra flux component in the model LCs, f_g, f_r, and f_i in the g, r, and i band, respectively. We also allow for possible uncertainties in the absolute flux calibration of the light curves, by including a wavelength-dependant magnitude offset in the same bands.

Table 2.  Comparison of Irradiated Binary Models for J2215

Param.	This Worka	SH14b	R15c	R16d
f	${0.95}_{-0.01}^{+0.01}$	1.00 ± 0.01	≃1	0.905 ± 0.004
q	0.144 ± 0.002	0.155–0.180	0.145	?
i(◦)	${63.9}_{-2.7}^{+2.4}$	49.5–54.3	88.8	83 ± 6
D(kpc)	2.9 ± 0.1	n.a.	3.9?	3.0?
K2(km s−1)	412.3 ± 5.0	329–382e	407.8f	?
η	${2.9}_{-0.2}^{+0.3}$	0.083 ± 0.001	0.97	2.16 ± 0.15
M1(M⊙)	${2.27}_{-0.15}^{+0.17}$	1.9–2.7	1.59	?
M2(M⊙)	${0.33}_{-0.02}^{+0.02}$	0.34–0.44	0.23	?
T2 (K)	${5630}_{-71}^{+52}$	3765–3945	6220	?
TN(K)	5280–5920	n.a.	>6000	6416 ± 58
${T}_{{\rm{D}}}$(K)	7800–8550	4876–5090	9264g	9000
${\chi }^{2}$/dof	164/126 = 1.30	1.5	4.4	1018/263 = 3.9

Notes. Ad-hoc constraints on the models are marked in italics.

^aOur xrbcurve fit to three-band photometric light curves (g, r and i) and two-band spectroscopic radial velocity curves (Mg i and Balmer). See the text for details. ^belc fit to B, V, and R light curves (Schroeder & Halpern 2014, and references therein). Ranges from both their NextGen and PHOENIX models. ^celc fit to the same three-band BVR light curves and one radial velocity curve (Romani et al. 2015). Their preferred "HiT" model (no errors reported). ^dicarus fit to the same BVR LCs and RVC, with irradiation from an extended intrabinary shock and temperature constraint, from Romani & Sanchez (2016). Their preferred "IBS-Td" model. ^ePredicted from the LC model fit. ^fDerived from their best-fit q, i, and M₁. ^gDerived from their best-fit M₁, M₂, and L_irr.

Download table as:  ASCIITypeset image

The model parameters that determine the shape and amplitude of the optical LCs and RVCs are i, T₂, η, D, and K₂ and the extra flux components f_g, f_r, and f_i. There are also a number of extra parameters: the phase shift for the LCs and RVCs (discussed below) as well as the systemic velocities for the Mg i-triplet and the Balmer-line RVCs (which we set to be the same). The $g^{\prime} ,r^{\prime} ,i^{\prime}$-band LCs were phase folded and averaged into 37, 37, and 28 orbital phase bins, respectively. The Mg i-triplet and Balmer-line RVCs contain 20 data points each (see Section 2.3 for details). Given that there are five different data sets with different numbers of data points, to optimize the fitting procedure, we assigned relative weights to them. After our initial search of the parameter space, which resulted in a good solution, we scaled the uncertainties in each data set (i.e., the LCs and the RVCs) so that the total reduced χ² of the fit was ∼1 for each data set separately. The MCMC fitting procedure was repeated to produce the final set of parameters, which were used to determine M₁ and M₂.

Our physical model fits simultaneously the Balmer and Mg i RVCs of J2215 (Figure 7, right) as well as the optical LCs in three bands (Figure 8), with a global reduced χ² of 1.30 for 126 d.o.f. Small residuals are apparent in the LC fits (especially near Φ_orb = 0.8) and in the RVC fits (at ${{\rm{\Phi }}}_{\mathrm{orb}}$ = 0.5). We show in Figure 9 the parameter distributions from our xrbcurve fits to J2215. All model parameters are well constrained, and the overall agreement between the data and model is good. Our best-fit D = 2.9 ± 0.1 kpc is fully consistent with the 3 kpc value found independently from the MSP dispersion measure (Hessels et al. 2011).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2 presents our best-fit values and their 1σ uncertainties, compared to previous studies of J2215. First, as already pointed out by Romani et al. (2015), previous attempts at determining orbital parameters exclusively based on photometric measurements and modeling have been unsuccessful (Schroeder & Halpern 2014; see also Breton et al. 2013). This is clear, e.g., from the K₂ velocities predicted by those photometric fit results, inconsistent with our measured values (Table 2). Thus we conclude that, at least in the presence of strong irradiation, dynamical information is required in order to find a reliable orbital solution.

Second, the orbital inclination depends strongly on the temperatures of both sides of the companion. Indeed, a larger temperature contrast between both sides requires a smaller inclination angle to produce the observed peak to peak magnitude difference. Our temperature constraints on the model are taken from quantitative temperature measurements at Φ_orb = 0.0 and 0.5 (Section 3.1), and lead to a robustly determined i = 639${}_{-2.7}^{+2.4}$. Therefore, independent constraints on the temperature at different orbital phases are also needed to find a robust solution.

Third, we have shown that a point-source irradiation binary model can fit satisfactorily the J2215 data, in clear contrast with previous results (Romani et al. 2015; Romani & Sanchez 2016). This discrepancy might also be due to the different temperature constraints, but a more detailed comparison is warranted. In any case, our results show that irradiation from an extended shock is not required to explain the optical properties of J2215. On the other hand, all four models do agree on the filling factor, showing a nearly Roche-lobe-filling companion in J2215. We find an additional nonstellar flux contribution in the range 0.035–0.07 mJy in all three bands (g, r, and i), between 2 and 10 times fainter than the companion star (at Φ_orb=0 and 0.5, respectively). This extra nonvariable flux component, with a rather flat spectral slope, might be due to synchrotron emission from an intrabinary shock, but at present this interpretation remains tentative. We can also rule out the presence of a quiescent disk (as suggested by Schroeder & Halpern 2014), based on the absence of broad hydrogen and helium emission lines regularly associated with disks.

Our best-fit model predicts a projected rotational velocity for the companion of V sin i =103 ± 1 km s⁻¹. Applying an optimal subtraction (Section 2.3) of the Mg i triplet line region around Φ_orb = 0, with a G5 template broadened in steps of 10 km s⁻¹ up to 210 km s⁻¹, we measure V sin i = 180 ± 20 km s⁻¹ (with a limb-darkening coefficient u = 0.8) and V sin i = 165 ± 15 km s⁻¹ (with no limb darkening, u = 0). Thus, taking into account the uncertainty on the amount of limb darkening, our current observational constraints put V sin i in the range 150–200 km s⁻¹. However, better spectral resolution spectra are required to measure this accurately and compare it to our model prediction (our GTC data have a resolution of 160 km s⁻¹; Section 2.2). Finally, our best-fit value of η implies a very high L_irr = [1.5${}_{-0.1}^{+0.3}$]×10³⁵ erg s⁻¹. An irradiating luminosity three times higher than the spin-down energy budget might be explained by, e.g., beaming/anisotropy of the pulsar wind (e.g., Philippov et al. 2015).

The flux distribution along the surface of the companion star is shown in Figure 10, a by-product of our best-fit model of J2215. These maps illustrate the drastic irradiation or "heating" effects of the pulsar wind and high-energy emission on the temperature distribution of the secondary star. Besides producing a strong temperature gradient between opposite sides (which we measure in Section 3.1), such strong heating shifts the effective center of the secondary star ("center of light") away from its center of mass. Because the strength of the lines varies throughout the surface of the star, this shift is different for different absorption lines. This in turn results in significant distortion of the integrated line profiles as well as the corresponding RVCs (e.g., Phillips et al. 1999; Shahbaz et al. 2000). Thus our model provides a natural explanation for the systematic difference in K velocities that we find and report in Section 3.2.

Zoom In Zoom Out Reset image size

Our model allows for phase shifts in the LC and both RVCs, which may arise physically from, e.g., asymmetric heating of the companion star. The best-fit values for the phase shifts in the Balmer RVC, the Mg i RVC, and the LCs are, respectively, ${\rm{\Delta }}{{\rm{\Phi }}}_{\mathrm{Balmer}}$ = 0.018 ± 0.003, ${\rm{\Delta }}{{\rm{\Phi }}}_{\mathrm{Mg}{\rm{I}}}$ = 0.0002 ± 0.0039 (i.e., consistent with 0 within 1σ), and ${\rm{\Delta }}{{\rm{\Phi }}}_{{\rm{g}}}$ = 0.0082 ± 0.0007 (where the quoted errors are again 1σ statistical). Similar LC phase shifts ∼0.01 have been reported for J2215 (Schroeder & Halpern 2014; Romani et al. 2015; Romani & Sanchez 2016), and interpreted as evidence for asymmetric heating from an intrabinary shock.

However, since the model does not take into account errors on the orbital phases, one must consider carefully the systematic uncertainty on Φ_orb before interpreting such phase shifts. If we used the pulsar timing reference epoch to calculate ${{\rm{\Phi }}}_{\mathrm{orb}}$ (which dates from 2009 December 21, from Abdo et al. 2013), we should include the orbital period derivative ${\dot{P}}_{\mathrm{orb}}$ (not doing so would yield a propagated uncertainty to our 2014 August 02 epoch of 0.02 orbital cycles). But assuming a constant ${\dot{P}}_{\mathrm{orb}}$ over a 5 year timespan is problematic, since redbacks and black widows are known to show large, erratic changes in ${\dot{P}}_{\mathrm{orb}}$ on timescales of years (e.g., Archibald et al. 2013). Instead, we use our own spectroscopic reference time (T0, see Section 3.2) with a systematic uncertainty of 0.0008 days, which corresponds to a systematic uncertainty on Φ_orb of δΦ_orb = 0.0046 cycles.

In light of this, the only phase shift that we deem (marginally) significant is ΔΦ_Balmer, which is larger than 0 at the 3.3σ confidence level (this shift can be seen by close inspection of Figure 7). Thus we conclude that, even if orbital phase shifts are necessary to allow a good statistical fit to the data, their significance and interpretation depend critically on the systematics of the orbital ephemeris. From our data and model of J2215, after quantifying the systematic error on our orbital phases, we find one relatively small phase shift in the Balmer-line RVC and thus only marginal evidence for heating asymmetries.

In a careful study of the strongly irradiated companion star of PSR J2215+5135, we find that the radial velocities depend systematically on the spectral features and spectral range used to measure them (Section 3.2). Thanks to the high-quality GTC optical spectra, we show that magnesium lines yield semiamplitudes of the RVC that are always higher than those inferred from hydrogen Balmer lines or from an undetermined blend of lines. We argue that this new systematic effect arises from the extreme contrast between the the cold and heated sides of the companion star, for which we measure temperatures of T_N = 5660${}_{-380}^{+260}$ K and T_D = 8080${}_{-280}^{+470}$ K, respectively (Section 3.1). Under these circumstances, Balmer lines trace the hot face of the companion star of J2215 in its orbit around the pulsar, while Mg i triplet lines trace its cold, unperturbed face.

We modeled both RVCs and the optical LCs in three bands, and found the center of mass velocity to be K₂ = 412.3 ±5.0 km s⁻¹ (Section 4). Because the center of light of the Mg i lines yields higher velocities (K_Mg = K₂+ΔK, ΔK > 0), it is slightly shifted outwards from the center of mass of the companion. Using the simple expression for the center-of-light displacement ΔR/a = ΔK/(K₂(1+q)) (Wade & Horne 1988) and our best-fit orbital parameters (Table 2), we estimate ${\rm{\Delta }}{\text{}}{R}_{\mathrm{Mg}{\rm{I}}}/{\text{}}{R}_{\mathrm{RL}2}\,\simeq \,+0.11$ (i.e., about 11% of the Roche-lobe radius of the companion). Balmer lines, on the other hand, suffer from a stronger inwards displacement of the center of light, which can yield to systematic errors in dynamical mass measurements if not corrected. We estimate this displacement as above using the measured K_Balmer, and find ${\rm{\Delta }}{\text{}}{R}_{\mathrm{Balmer}}/{\text{}}{R}_{\mathrm{RL}2}\,\simeq \,-0.24$ (i.e., a 24% shift relative to ${R}_{\mathrm{RL}2}$).

Given the large (≃10%) differences in the measured K values and the M₁∝K³ relation, this has important consequences for the measured pulsar masses. For instance, using the different measured K values and our best-fit orbital solution (Section 4) yields inconsistent values for the neutron star mass: M₁ = 1.88 ± 0.16 M_⊙ (from K_Balmer = 382.8 ± 4.7 km s⁻¹) and M₁ = 2.49 ± 0.23 M_⊙ (from K_Mg = 420.2 ± 6.2 km s⁻¹). Thus we find that, in the presence of strong irradiation, the systematic error on K may be equally or more important than the uncertainty on the orbital inclination, i.

To circumvent this and drastically reduce the systematics on K measurements, we have put forward a new method: we measure K velocities using different sets of lines in order to "bracket" the center-of-mass velocity of the companion star. We have deemed as optimal the Balmer and Mg i triplet lines (Section 2.3) that provide a lower and an upper limit on K₂, respectively. These could, in principle, be replaced by other sets of lines that trace the movement of both the cold/dark and irradiated/bright sides of the star. While the traditional "K correction" relies on model assumptions or on the transient nature of irradiation in some systems (Section 5.2), our method provides a direct way of quantifying this correction from the same spectra of the irradiated companion star.

To our knowledge, the only similar studies in the literature involve a combination of emission and absorption lines in white dwarf (WD) binaries (e.g., Parsons et al. 2010; Rodríguez-Gil et al. 2015). Our method relies exclusively on absorption lines from the star's atmosphere, thus eluding the uncertainty sometimes associated with the exact site where Balmer/He/Bowen emission lines are formed. This "empirical K correction" is of particular relevance for the emerging population of MSPs in compact binaries, and it should also be applicable in the broader context of semidetached binaries with strong irradiation/heating effects.

In summary, our findings show that metallic lines, in general, and Mg i lines, in particular, offer a much less distorted view of the center of mass of the secondary star, opening a new way to measure masses in strongly irradiated compact binary MSPs. In the relevant temperature and spectral type range for J2215 and most redback and black widow companions (spectral types A through M), Balmer lines are more sensitive to temperature than Mg i triplet lines (see the Appendix, Figure 13). Indeed, when going from a spectral type A5 to a G5, the EW of Balmer lines decreases by a factor of ≃7 ($48\to 7\,\mathring{\rm A}$), whereas the EW of Mg i triplet lines increases by a smaller factor ≃3 ($1.52\to 4.75\,\mathring{\rm A}$). This may explain why the effects of irradiation are more drastic on Balmer absorption lines than on metallic lines.

Irradiation in compact binaries has been studied in the context of black hole and neutron star LMXBs as well as WD binaries (dwarf novae, DN, post-common envelope binaries, PCEBs, and asynchronous polars, AP). In those systems, irradiation proceeds mostly through X-ray and UV photons from a hot WD or innermost accretion disk. The ratio of the maximum irradiating flux at the companion's surface (near the inner Lagrangian point) over the companion's intrinsic unperturbed flux provides a good way to quantify the importance of irradiation in close binaries: ${f}_{\mathrm{irr}}=\tfrac{{L}_{\mathrm{irr}}\,{R}_{2}^{2}}{{L}_{2}\,{a}^{2}}$. The effects of irradiation on the measured K velocities, on the other hand, are typically parameterized in terms of the so-called K correction, which we define as the ratio between observed and center-of-mass K values: ${f}_{K}=\tfrac{{K}_{\mathrm{obs}}}{{K}_{2}}$. This correction is estimated in the literature in a number of different ways, e.g., by comparing outburst and quiescence K values (Hessman et al. 1984) or by simulating and fitting RVCs with irradiation models (Wade & Horne 1988; Phillips et al. 1999; Muñoz-Darias et al. 2005).

In some cases (compact LMXBs in outburst, very hot WDs), chromospheric/fluorescent emission lines are formed on the inner face of the companion star, leading to a lower limit on the center of mass velocity semi-amplitude K₂ (i.e., f_K < 1). This is the case of Bowen fluorescence lines in LMXBs. It is interesting to compare J2215 with the accreting millisecond pulsars SAX J1808.4–3658 and XTE J1814–338, with Bowen-line outburst K corrections of ${f}_{K}$ = 0.90 and 0.81, respectively (Cornelisse et al. 2009; Wang et al. 2017). These are analogous to (and possibly evolutionary precursors of) black widow and redback MSPs. Assuming an irradiating luminosity in outburst ${L}_{X}\sim {10}^{36}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, we estimate their maximum ratio of irradiating to intrinsic flux and find extreme values, f_irr ∼10³–10⁴. For the persistent atoll sources 4U 1636–536 and 4U 1735–444, with f_K ≃ 0.77 (Casares et al. 2006), we estimate even higher f_irr ∼ 10⁵.

Balmer/He i emission lines in WD binaries are also associated with strong, yet less extreme, irradiation. The AP V1500Cyg, for instance, has an estimated f_K = 0.65 (Balmer and He ii lines, Horne & Schneider 1989) and we find f_irr ∼ 50 (using the stellar parameters in Schmidt et al. 1995). Similarly, we find f_irr ∼ 140 for the PCEB NNSer, with a reported f_K = 0.88 (Balmer lines, Parsons et al. 2010). This suggests that the irradiating flux must be at least 10 times higher than the intrinsic stellar flux in order to form emission lines on the heated face.

In other cases (DNe and BHCs in outburst), absorption lines are partially quenched on the inner face of the companion due to irradiation, so that the effective center of light for these lines is shifted toward the outer face and they provide an upper limit on K₂ (i.e., f_K > 1). Even though stellar atmospheres with external heating are poorly understood, the reduced absorption line strength is often attributed to the reduced vertical temperature gradient in the presence of an external UV/X-ray photon flux. This is the case of the WD binaries ZCha, UGem, and SSCyg, with ${f}_{K}$ estimated at 1.03, 1.04, and 1.26, respectively (Hessman et al. 1984; Wade & Horne 1988; Friend et al. 1990). For these systems, we find mild irradiation, with estimated ${f}_{\mathrm{irr}}$ values in the range 0.1–3. The BHC GRO J1655-40, on the other hand, with f_K ≃ 1.16, has a relatively luminous F6IV companion star that also leads to a mild f_irr ∼ 7 (Orosz & Bailyn 1997; Phillips et al. 1999).

In our case (J2215) and in compact binary MSPs, in general, the relativistic pulsar wind and gamma-ray emission are the dominant sources of irradiation. These feature typical spin-down luminosities $\dot{E}$ = 10³⁴–10³⁵ erg s⁻¹ and gamma-ray luminosities (L_γ) about 10 times lower. Their X-ray luminosities are two to five orders of magnitude lower than $\dot{E}$ and are thus less important in terms of irradiation (L_X = 10³⁰–10³² erg s⁻¹; Gentile et al. 2014; Linares 2014). Indeed, J2215 has L_X = 1.2 × 10³² erg s⁻¹ (0.5–10 keV), L_γ = 1.4 × 10³⁴ erg s⁻¹ (0.1–100 GeV), and $\dot{E}$ =5.3 × 10³⁴ erg s⁻¹ (Linares 2014; Acero et al. 2015; Breton et al. 2013, respectively; for a 2.9 kpc distance). Our measured radial velocity amplitudes K_Balmer and ${K}_{\mathrm{Mg}{\rm{I}}}$, together with the center-of-mass velocity from our best-fit model (K₂) imply K correction factors for J2215 of ${f}_{K}$ = 0.928 and ${f}_{K}$ = 1.019, respectively. In other words, the K values measured using Balmer and Mg i lines are 7.2% lower and 1.9% higher, respectively, than the true center-of-mass K₂.

For comparison, the K correction inferred by van Kerkwijk et al. (2011) for the black widow pulsar PSR B1957+20 was close to 8%: ${f}_{K}$ = 0.918 (where their K_obs was measured using a slightly wider spectral range and a G2 template). Thus, K corrections may be similarly important in black widow and redback binaries. This can be understood qualitatively from the two terms entering f_irr, which have opposite trends: ${R}_{2}^{2}/{a}^{2}$ is lower in black widows (lower irradiating flux due to smaller solid angle), but ${L}_{\mathrm{irr}}/{L}_{2}$ is higher compared to redbacks (as the black widow companions are less luminous).

Using our orbital solution (which gives R_RL2/a = 0.23), our best-fit value of L_irr ≃ 3×$\dot{E}$, and a companion luminosity ${L}_{2}\simeq$ 5.3 × 10³² erg s⁻¹ (from our best-fit T₂ = 5630 K and R₂ =0.39 R_⊙), we estimate ${\text{}}{f}_{\mathrm{irr}}\sim$ 15. In other words, the irradiating flux at the companion's heated face is up to 15 times higher than the intrinsic stellar flux. Compared with UV and X-ray photons, Gamma-ray photons and relativistic particles from the pulsar wind are expected to penetrate deeper into the companion atmosphere. This leads to deep/internal heating of the companion's inner face, so it is not surprising to find no emission lines in J2215, and no quenching of absorption lines either (Figure 4).

We have argued that our "empirical K correction" removes a critical systematic uncertainty in irradiated systems: the difference between center of light and center of mass of the companion. In J2215, our K₂ = 412.3 ± 5.0 km s⁻¹ yields a mass function of 1.2 M_⊙, an absolute lower limit on the pulsar mass. Combined with the tight constraints on q (Section 4; Abdo et al. 2013), this implies a minimum neutron star mass of 1.6 M_⊙. Thus we find that J2215 contains a neutron star more massive than the "canonical" 1.4 M_⊙ double neutron stars, adding to the growing number of such systems (Demorest et al. 2010; Antoniadis et al. 2013; Özel & Freire 2016).

Furthermore, we have shown that J2215 could harbor the most massive pulsar known to date. Our best-fit model yields a very massive neutron star with ${M}_{1}$ = 2.27${}_{-0.15}^{+0.17}$ M_⊙, and a well-constrained inclination of i = 63.9${}_{-2.7}^{+2.4}$ (see Figure 11). This is more massive than the previous well-established record holder (2.01 M_⊙, Antoniadis et al. 2013), at the 97% confidence level. The results of Romani et al. (2015) and Romani & Sanchez (2016) are clearly at odds, although their nearly edge-on model (with a neutron star mass of 1.6 M_⊙) has no uncertainties reported and admittedly fails to describe the data. In a previous study of the original black widow pulsar, van Kerkwijk et al. (2011) found a pulsar mass that was similarly high (2.40 ± 0.12 M_⊙). The main advances of the method we have presented here are (i) an empirical K correction, based on radial velocities measured with two different sets of lines (Section 3.2) and (ii) independent constraints on the temperature imposed on the model, based on absorption line strengths (Section 3.1). We argue that these two new advances make our results more robust compared to previous work (van Kerkwijk et al. 2011; Romani et al. 2015; Romani & Sanchez 2016). There may still be unknown or highly uncertain systematic effects; however, biasing the best-fit model inclination and thus the dynamical mass measurements (see, e.g., the discussion in van Kerkwijk et al. 2011, their Section 4).

Zoom In Zoom Out Reset image size

We therefore conclude that, if confirmed with an independent measurement of the orbital inclination, the massive neutron star in J2215 may place new constraints on the equation of state at supranuclear densities. This would push the limits of the most massive neutron stars in our Galaxy, setting a lower limit of 2.3 M_⊙ on their maximum mass. Since particle interactions in the core provide the pressure necessary to halt its collapse, the maximum mass of a neutron star places independent constraints on how these particles interact (Lattimer & Prakash 2007). For instance, exotic forms of matter such as hyperons or deconfined quarks have been proposed to exist in the central parts of a neutron star, yet they can hardly account for a neutron star as massive as the one we find in J2215 (see also Özel & Freire 2016).

With new Galactic MSPs being currently discovered at a rate of 10–30 per year (Lorimer 2018), the neutron star mass range will be explored further in the next decade, and is likely to continue widening. We have shown here that in the study of strongly irradiated pulsar companions, a controlled measurement of temperatures and velocities throughout the orbit is possible with current instruments and key to finding a robust dynamical solution. Our novel technique, which combines velocity measurements with different absorption lines, temperature measurements, and physical modeling of the binary, should provide a path forward for dynamical mass measurements in this growing population.