PSR J2234+0611: A New Laboratory for Stellar Evolution

PSR J2234+0611 was discovered in the Arecibo Observatory 327 MHz Drift Scan Survey in December 2012 (Deneva et al. 2013). After discovery of the pulsar, initial follow-up observations were performed, also with the Arecibo 305 m telescope, using the "L-wide" receiver at a center frequency of 1.5 GHz and recorded with the Puerto Rican Ultimate Pulsar Processing Instrument (PUPPI) in search mode, allowing for offline folding of each observation to get the observed period at each epoch. The preliminary orbital parameters resulting from these observations were already reported in Deneva et al. (2013).

We then folded the data using the new orbit and began to refine the timing solution by generating pulse times-of-arrival (ToAs) and performing pulsar timing analysis using Tempo.²³ Subsequent data was recorded using PUPPI in coherent dedispersion and online folding mode. Figure 2 shows the profile for PSR J2234+0611 at 430 and 1.5 GHz from roughly 30 minute duration coherent fold mode observations.

Zoom In Zoom Out Reset image size

PSR J2234+0611 was immediately found to have excellent timing precision and therefore was added to the pulsar timing array's (PTA's) efforts to detect low frequency gravitational waves, in particular to the North American Nanohertz Observatory for Gravitational Waves (NANOGrav, Demorest et al. 2013) PTA. Observations of the pulsar have continued under that project, using the Arecibo 305 m radio telescope with the "L-wide" receiver (with frequency coverage between 1130 and 1730 MHz) and the 430 MHz receiver with a cadence of about 3 weeks. For both types of observations, the PUPPI back-end was used, with coherent dedispersion and folding mode, as for other PTA pulsars; these observations are described in detail by Arzoumanian et al. (2018) but extend later in time than the data presented in that paper. Current timing solution parameters from data spanning 5 yr are given in Tables 1 and 2.

Table 1.  Non-binary Parameters for PSR J2234+0611

Observation and data reduction parameters
Reference Epoch (MJD)	56794.093186
Span of timing data (MJD)	56347–58291
Number of ToAs	5882
Solar wind parameter, n0 (cm−3)	6
Overall individual ToA rms residual (μs)	0.58
rms residual for incoherent L-band (μs)	0.35
rms residual for coherent 430 MHz (μs)	1.51
rms residual for coherent L-band (μs)	0.59
χ2	5891.76
Reduced χ2	1.013
Spectral parameters
Mean flux density at 430 MHz, S430 (mJy)	1.3
Mean flux density at 1400 MHz, S1400 (mJy)	0.6
Spectral Index, α	−0.7
Astrometric and spin parameters
R.A., α (J2000)	22:34:23.073090(2)
Decl., δ (J2000)	06:11:28.68633(7)
Proper motion in α, μα (mas yr−1)	25.30(2)
Proper motion in δ, μδ (mas yr−1)	9.71(5)
Parallax, ϖ (mas)	1.03(4)
Spin frequency, ν (Hz)	279.5965821510426(5)
Spin frequency derivative, $\dot{\nu }$ (10−16 Hz s−1)	−9.3920(1)
Dispersion measure, DM (pc cm−3)	10.778
Derived parameters
Galactic longitude, l	+72.99
Galactic latitude, b	−43.01
Magnitude of proper motion, μ (mas yr−1)	27.10(2)
Position angle of proper motion, Θμ (deg, J2000)	69.0(1)
Position angle of proper motion, Θμ (deg, Galactic)	111.5(1)
DM-derived distance, d1 (kpc)	0.68
DM-derived distance, d2 (kpc)	0.86
Parallax-derived distance, d (kpc)	0.97(4)
Galactic height, z (kpc)	−0.651(26)
Transverse velocity, vT (km s−1)	123(5)
Spin period, P (ms)	3.576581631673107(6)
Spin period derivative, $\dot{P}$ (10−20 s s−1)	1.20142(1)
Intrinsic spin period derivative, ${\dot{P}}_{\mathrm{int}}$ (10−20 s s−1)	${0.647}_{-0.025}^{+0.023}$
Surface magnetic flux density, B0 (108 Gauss)	1.5
Characteristic age, τc (Gyr)	8.8
Spin-down power, $\dot{E}$ (1033 erg s−1)	5.6

Note. Timing parameters and 1σ uncertainties derived using Tempo in Barycentric Dynamical Time (TDB), using the DE421 Solar System ephemeris, the NIST UTC time timescale and the DDGR orbital model. d₁ is derived using the NE2001 (Cordes & Lazio 2002) Galactic model, d₂ using the YMW16 (Yao et al. 2017) Galactic model. Estimate of v_T, ${\dot{P}}_{\mathrm{int}}$, and derived parameters assume distance from the measured parallax and its uncertainty.

Download table as:  ASCIITypeset image

The timing analysis of the PUPPI data is similar to that described by Arzoumanian et al. (2018). The ToAs are derived from the integrated pulse profiles using the standard PSRCHIVE routines. The ToA analysis is made using Tempo. To convert the telescope ToAs (corrected to the International Bureau of Weights and Measures version of Terrestrial Time, TT) to the Solar System barycentre, we used the Jet Propulsion Laboratory's DE436 solar system ephemeris; the resulting timing parameters are presented in Barycentric Dynamical Time (TDB). We used the same method used by NANOGrav to estimate variations of the dispersion measure (DM), but with the ToAs grouped in intervals of 32 days (the orbital period), instead of 6 days as is the norm for the NANOGrav pulsars. DM values are reported as offsets relative to an arbitrary fiducial value of 10.778 pc cm⁻³.

We used three orbital models to analyze the data, all based on the description of Damour & Deruelle (1985, 1986). The first is the "DDGR" model, which assumes the validity of general relativity (GR) and where we fit directly for the total mass of the system (M_tot) and the companion mass (M_c). The second model is basically the theory-independent DD model, but with the orthometric parameterization of the Shapiro delay described by Freire & Wex (2010); this was implemented in Tempo by Weisberg & Huang (2016), where it is designated as the "DDFWHE" orbital model. The third model is again based on the DD model but takes into account the kinematic effects described by Kopeikin (1995, 1996); this was implemented in Tempo2 by Edwards et al. (2006), where it is designated as the "T2" model; it was implemented in Tempo by one of us (I.H.S.), where it is designated as the "DDK" model.

The reason for the usage of these three orbital models is that, as we will show, no single model alone fully captures all the constraints on the masses and orbital orientation of this system. In the DDFWHE and DDK solutions, we used the Einstein delay calculated in the DDGR solution; the reason for this is because it cannot be determined independently with our data, and because it is strongly correlated with $\dot{x}$ (see A. Ridolfi et al. 2018, in preparation). Furthermore, the orthometric ratio of the Shapiro delay (ς) in the DDFWHE solution and the orbital inclination (i) in the DDK solution are derived from the $s\equiv \sin i$ parameter calculated by the DDGR solution; the reason being the extremely small signature of the Shapiro delay. In Section 3 we discuss the significance of these parameters.

As part of the NANOGrav data analysis procedures, the data have been flux and polarization calibrated, allowing straightforward measurements of the polarization profile (Figure 2) and mean flux density. We have taken flux density values from a preliminary analysis of the upcoming 12.5 yr data release (Z. Arzoumanian et al. 2018, in preparation). The data in this preliminary release was polarization and flux calibrated using the same methods as the NANOGrav 9 yr data release (The NANOGrav Collaboration et al. 2015). However, the observed flux density for PSR J2234+0611 varies over a fairly wide range due to scintillation by the interstellar medium. Using psrflux from the PSRCHIVE pulsar suite, we calculated the mean value from 43 observations at 430 MHz, ranging from 0.2 to 5.3 mJy and 50 observations at 1.5 GHz, ranging from 0.03 to 3.3 mJy, to get an estimate for the mean flux density at these frequencies. The resulting mean values and spectral index are given in Table 1.

For this pulsar, we obtain a highly significant measurement of the parallax, 1.05(4) mas (all uncertainties are 68.3% confidence limits) resulting in a pulsar distance d of 0.95(4) kpc. This distance can be compared with the prediction of the DM models. The NE2001 model (Cordes & Lazio 2002) predicts a distance of 0.68 kpc, while the YMW16 model (Yao et al. 2017) predicts a distance of 0.86 kpc. A relative uncertainty of about 20% is generally assigned to these estimates. Our parallax measurement is certainly in better agreement with the YMW16 model.

This measurement, together with the measurement of the proper motion, allows a relatively accurate measurement of the Heliocentric transverse velocity, 123(5) km s⁻¹. Combining this with the systemic radial velocity of −20(34) km s⁻¹ measured by Antoniadis et al. (2016a), we obtain a 3D heliocentric velocity of ${124}_{-5}^{+10}\,\mathrm{km}\,{{\rm{s}}}^{-1}$. This velocity is smaller than that used in the detailed analysis of the Galactic motion of PSR J2234+0611 made by Antoniadis et al. (2016a), mostly because they were using a preliminary value of the parallax that yielded a larger distance, however the qualitative conclusions obtained by Antoniadis et al. (2016a) remain valid: the 3D velocity of this system is similar to what has been observed for other nearby recycled pulsars (e.g., Gonzalez et al. 2011). We will return to this topic in Section 5, particularly in the discussion on the formation of the system.

For any assumed distance we can estimate the magnitude of the kinematic effects on the variation of the Doppler shift factor (D) using the simple expressions provided by Shklovskii (1970) for the effect of the centrifugal acceleration (proportional to the square of the total proper motion, μ) and Damour & Taylor (1991) for the effect of the difference in the Galactic accelerations of the pulsar's system and the Solar System projected along the direction from the pulsar to the Earth, a_l:

$$\begin{eqnarray}&&\displaystyle \frac{\dot{D}}{D}\,\equiv -\displaystyle \frac{{\mu }^{2}d+{a}_{l}}{c}\end{eqnarray} \tag{ 1 }$$

where c is the speed of light. In order to estimate a_l, we use the expressions presented by Lazaridis et al. (2009), where the equation for the vertical acceleration should be valid to a Galactic height of ∼±1.5 kpc (the Galactic height of PSR J2234+0611 is −0.651(26) kpc). In those expressions we use the distance to the center of the Galaxy measured by the GRAVITY experiment (Gravity Collaboration et al. 2018), r₀ = 8.122(31) kpc and a revised value for the rotational velocity of the Galaxy derived using the latter r₀ (McGaugh 2018), v_Gal = 233.3 km s⁻¹. We obtain a_l/c = −1.53 × 10⁻¹⁹ s⁻¹ (for a comparison, we can use the Galactic model presented by McMillan 2017 to obtain a_l/c = −1.76 × 10⁻¹⁹ s⁻¹, which is a similar number). For the proper motion term we obtain μ²d/c = 1.702 × 10⁻¹⁸ s⁻¹, an order of magnitude larger. Adding both terms, we obtain $\dot{D}/D\,=-1.550\,\times {10}^{-18}\,{{\rm{s}}}^{-1}$.

The contribution of this effect to the spin period derivative is given by ${\dot{P}}_{\mathrm{kin}}=-P\,\dot{D}/D={5.54}_{-0.25}^{+0.23}\,\times \,{10}^{-21}\,{\rm{s}}\,{{\rm{s}}}^{-1}$. Subtracting this from the observed $\dot{P}$ in Table 1 we obtain the intrinsic spin period derivative (${\dot{P}}_{\mathrm{int}}={6.47}_{-0.25}^{+0.23}\,\times {10}^{-21}\,{\rm{s}}\,{{\rm{s}}}^{-1}$), which is about half of the observed $\dot{P}$. From this and the spin period P, we derive a surface magnetic flux density B₀ ≃ 1.5 × 10⁸ G, the rate of loss of rotational energy $\dot{E}\,\simeq \,5.6\,\times {10}^{33}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, and a characteristic age τ_c ≃ 8.8 Gyr using the standard equations summarized by Lorimer & Kramer (2004). The cooling age for the WD companion is 1.5 Gyr, which according to Antoniadis et al. (2016a) is comparable to the age of the system. This is compatible with τ_c since the latter represents an upper limit for the age that assumes that the initial spin period P_init was much smaller than the currently observed P. Assuming an n = 3 braking index and an age of 1.5 Gyr, we obtain P_init = 3.25 ms.

This rate of change of the Doppler shift factor will also be a dominant contributor to the observed variation of the orbital period, ${\dot{P}}_{{\rm{b}},\mathrm{obs}}$. According to Lorimer & Kramer (2004):

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\dot{P}}_{{\rm{b}}}}{{P}_{{\rm{b}}}}\right)}^{\mathrm{obs}}=-\,\displaystyle \frac{\dot{D}}{D}\,+{\left(\displaystyle \frac{{\dot{P}}_{{\rm{b}}}}{{P}_{{\rm{b}}}}\right)}^{\mathrm{GW}}\,+{\left(\displaystyle \frac{{\dot{P}}_{{\rm{b}}}}{{P}_{{\rm{b}}}}\right)}^{\dot{{\rm{m}}}}\,+{\left(\displaystyle \frac{{\dot{P}}_{{\rm{b}}}}{{P}_{{\rm{b}}}}\right)}^{{\rm{T}}},\end{eqnarray} \tag{ 2 }$$

the first term, the kinematic contribution to ${\dot{P}}_{{\rm{b}},\mathrm{obs}}$ is given by ${\dot{P}}_{{\rm{b}},\mathrm{kin}}=-{P}_{{\rm{b}}}\,\dot{D}/D={4.28}_{-0.18}^{+0.19}\times {10}^{-12}\,\,{{\rm{s}}}^{-1}$. The second term in Equation (2) is due to the loss of orbital energy caused by the emission of gravitational waves. For PSR J2234+0611, this term is, assuming the validity of GR, given by ${\dot{P}}_{{\rm{b}},\mathrm{GR}}\,=2.62\times {10}^{-17}\,{\rm{s}}\,{{\rm{s}}}^{-1}$ (this is the estimate provided by the DDGR model for the masses derived by that model). This is about 5 orders of magnitude smaller than ${\dot{P}}_{{\rm{b}},\mathrm{obs}}$ and its uncertainty. The third term is caused by radiative mass loss from the system. Assuming that this is dominated by the loss of rotational energy for the pulsar, it is given by Damour & Taylor (1991):

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\dot{P}}_{{\rm{b}}}}{{P}_{{\rm{b}}}}\right)}^{\dot{{\rm{m}}}}=\displaystyle \frac{8\pi G}{{T}_{\odot }{c}^{5}}\displaystyle \frac{I}{{M}_{\mathrm{tot}}}\displaystyle \frac{{\dot{P}}_{\mathrm{int}}}{{P}^{3}}\,\sim \,3.8\,\times \,{10}^{-21}{{\rm{s}}}^{-1}\end{eqnarray} \tag{ 3 }$$

where ${T}_{\odot }={{GM}}_{\odot }{c}^{-3}=4.925490947\,\mu {\rm{s}}$ is a solar mass (M_⊙) in time units, c is the speed of light, and G is Newton's gravitational constant, I is the moment of inertia of the pulsar, I ≃ 10³⁸ kg m². Thus, ${\dot{P}}_{{\rm{b}}}^{\dot{m}}=1.05\,\times \,{10}^{-14}\,{\rm{s}}\,{{\rm{s}}}^{-1}$, which is about 40 times smaller than ${\dot{P}}_{{\rm{b}},\mathrm{kin}}$. Finally, the last term in Equation (2) is caused by tidal dissipation. For PSR J2234+0611, this term should be negligible: the WD mass and atmospheric parameters indicate that the star is well within its Roche lobe and no mass loss occurs. Consequently, the tidal dissipation timescale (Zahn 1977) is of the order 20 Gyr, well above the characteristic age of the pulsar; τ_c ≃ 1.5 Gyr.

Thus, the only relevant term appears to be ${\dot{P}}_{{\rm{b}},\mathrm{kin}}$. This matches the observation (${\dot{P}}_{\mathrm{obs}}=3.3\,\pm \,2.5\,\times \,{10}^{-12}\,{\rm{s}}\,{{\rm{s}}}^{-1}$ for the DDGR and DDFWHE solutions, ${\dot{P}}_{\mathrm{obs}}=4.9\,\pm \,2.5\,\times {10}^{-12}\,{\rm{s}}\,{{\rm{s}}}^{-1}$ for the DDK solution, see Table 2); while for the DDK solution we have a 2σ "detection" of this effect.

The post-Keplerian effect measured to the highest significance for PSR J2234+0611 is the rate of advance of periastron, $\dot{\omega }$. According to Lorimer & Kramer (2004), the observed effect is given, in the absence of a third component in the system, by:

$$\begin{eqnarray}&&{\dot{\omega }}_{\mathrm{obs}}={\dot{\omega }}_{\mathrm{rel}}+{\dot{\omega }}_{{\rm{k}}}+{\dot{\omega }}_{\mathrm{SO}}.\end{eqnarray} \tag{ 4 }$$

The third term is caused by spin–orbit coupling, a result of the finite size of the companion white dwarf; for wide systems like PSR J2234+0611 this effect is negligible.

The first term is caused by relativistic effects. Assuming GR, we can estimate the total mass of the binary, M_tot (Robertson 1938) from ${\dot{\omega }}_{\mathrm{rel}}$ by inverting the well-known expression derived by Taylor & Weisberg (1982):

$$\begin{eqnarray}&&{M}_{\mathrm{tot}}=\displaystyle \frac{1}{{T}_{\odot }}{\left[\displaystyle \frac{{\dot{\omega }}_{\mathrm{Rel}}}{3}(1-{e}^{2})\right]}^{\tfrac{3}{2}}{\left(\displaystyle \frac{{P}_{{\rm{b}}}}{2\pi }\right)}^{\tfrac{5}{2}}.\end{eqnarray} \tag{ 5 }$$

The DDGR model assumes that ${\dot{\omega }}_{\mathrm{rel}}={\dot{\omega }}_{\mathrm{obs}}$, i.e., all other terms are assumed to be negligible. As we see below this assumption cannot be made for PSR J2234+0611. From this assumption, the DDGR model obtains M_tot = 1.6798(29) M_⊙. The $\dot{\omega }$ provided by the DDFWHE solution yields, assuming GR, an identical M_tot. This constraint is represented by the solid red line in Figure 4.

Zoom In Zoom Out Reset image size

However, the $\dot{\omega }$ measured by the DDK solution is smaller than that measured by the DDFWHE model by a small (${\rm{\Delta }}\dot{\omega }=9.60\times {10}^{-6}\,\deg \,{\mathrm{yr}}^{-1}$) but highly significant (9.3σ) amount. The reason is that, for PSR J2234+0611, as for another wide, precisely timed system, PSR J1903+0327 (Freire et al. 2011) the second term in Equation (4), ${\dot{\omega }}_{{\rm{k}}}$, is larger than the measurement uncertainty. This term is given by Kopeikin (1995), here re-arranged as in Freire et al. (2011):

$$\begin{eqnarray}&&{\dot{\omega }}_{{\rm{k}}}=\displaystyle \frac{\mu }{\sin i}\cos \left({{\rm{\Theta }}}_{\mu }-{\rm{\Omega }}\right),\end{eqnarray} \tag{ 6 }$$

where Θ_μ is the position angle (PA) of the proper motion and Ω is the PA for the line of nodes. In the DDK orbital model, the PAs are measured in Equatorial (J2000) coordinates, starting from North through East and an inclination smaller than 90° corresponds to a system where the line-of-sight component of the angular momentum points toward the Earth.

Although Θ_μ is measured directly from the proper motion (see Table 1), the orientation of the line of nodes Ω is generally harder to determine. In Figure 4, we display with the dashed red lines the total masses assuming minimal or maximal contributions of ${\dot{\omega }}_{{\rm{k}}}$ to ${\dot{\omega }}_{\mathrm{obs}}$, M (estimated from Equation (6) by setting $\cos \left({{\rm{\Theta }}}_{\mu }-{\rm{\Omega }}\right)=\pm 1$). This shows clearly that ${\dot{\omega }}_{{\rm{k}}}$ is potentially much larger than the uncertainty in the measurement of ${\dot{\omega }}_{\mathrm{obs}}$.

However, we can estimate the total mass more accurately since, in the DDK model, we can determine i and Ω with high precision (see details in Section 3.6, and DDK solution in Table 2). Using these values the model internally estimates ${\dot{\omega }}_{{\rm{k}}}$ and automatically subtracts it from the "measured" ${\dot{\omega }}_{\mathrm{obs}}$, reporting only the part (presumably) caused by the relativistic effects, ${\dot{\omega }}_{\mathrm{rel}}$. Assuming GR, this yields a lower binary mass (M_tot = 1.6526(29) M_⊙) than estimated by the DDGR model. We consider this to be accurate since it takes the proper motion into account.

In the DDGR model, we not only obtain a precise (but in this case inaccurate) estimate for M_tot, but also a precise estimate for the companion mass (M_c = 0.300(13) M_⊙). Given the mass function of the system, the estimated M_c implies sin i ∼ 0.668; this implies either $i\,\sim 42\,\deg$ or i ∼ 138 deg. The measurement is possible because of the presence of the Shapiro delay, however, the fact that the Shapiro delay is detected at all for a system with an orbital inclination so far from edge-on (90°) is unusual. The detection in this case is the result of two factors: one is the high timing precision of the Arecibo observations of this pulsar, the second is the large eccentricity of the orbit; the latter helps separate the Shapiro delay from the normal "Roemer" delays caused by the geometry of the orbital motion relative to the line of sight.

The far from edge-on inclination means that the Shapiro delay is not easy to measure. When using the DDFWHE model to fit for both Shapiro delay parameters, h₃ and ς, both values are measured to relatively low confidence. In order to better quantify the Shapiro delay, we use the best-fit value of $s\equiv \sin i$ that corresponds to the masses determined by the DDGR model (s = 0.6677654 ...) to derive (Freire & Wex 2010):

$$\begin{eqnarray}&&\varsigma =\displaystyle \frac{s}{1+\sqrt{1-{s}^{2}}}=0.382811...,\end{eqnarray} \tag{ 7 }$$

this is represented by the blue dashed line in Figures 4 and 5. Fixing this in the DDFWHE model, we obtain a significant h₃ = 82 ± 14 ns; an unusually small value that is a consequence of the inclination of the system. The mass and inclination constraints introduced by this measurement and its ±1σ uncertainties are shown by the solid blue curves in Figure 4. The region where these h₃ lines cross the $\dot{\omega }$ lines provides a good explanation of the DDGR estimate for M_c and its related uncertainty.

Zoom In Zoom Out Reset image size

As seen in Table 2, both the DDGR and DDFWHE timing solutions contain a precise measurement of a change in the projected semimajor axis ($x=a\ \sin \ i$ where a is the semimajor axis and i is the orbital inclination) of the pulsar's orbit, $\dot{x}=-2.79(7)\times {10}^{-14}$ $\mathrm{lt}-{\rm{s}}\,{{\rm{s}}}^{-1}$, thus ${(\dot{x}/x)}^{\mathrm{obs}}=-1.99(5)\,\times {10}^{-15}\,{{\rm{s}}}^{-1}$. Following Lorimer & Kramer (2004), the observed change in $\dot{x}$ can, in the absence of a third object in the system, be written in terms of various contributions as:

$$\begin{eqnarray}\begin{array}{rcl}{\left(\displaystyle \frac{\dot{x}}{x}\right)}^{\mathrm{obs}} & = & {\left(\displaystyle \frac{\dot{x}}{x}\right)}^{{\rm{k}}}\,+{\left(\displaystyle \frac{\dot{x}}{x}\right)}^{\mathrm{GW}}+\displaystyle \frac{d{\epsilon }_{A}}{{dt}}\\ & & -\ \displaystyle \frac{\dot{D}}{D}\,+{\left(\displaystyle \frac{\dot{x}}{x}\right)}^{\dot{m}}\,+{\left(\displaystyle \frac{\dot{x}}{x}\right)}^{\mathrm{SO}}.\end{array}\end{eqnarray} \tag{ 8 }$$

The first term is caused by the changing geometry due to the motion of the system relative to the Earth, it is given by (Kopeikin 1995):

$$\begin{eqnarray}&&{\dot{x}}_{{\rm{k}}}=x\mu \cot i\sin ({{\rm{\Theta }}}_{\mu }-{\rm{\Omega }}),\end{eqnarray} \tag{ 9 }$$

where we have, again, re-written the terms as in Freire et al. (2011), except for the latter's negative sign, so that we are in the right-handed convention being used in the DDK model. As we will see below, this is the only term that can account for the observations.

The second term is from the decrease of the size of the orbit caused by gravitational wave emission; this is given by

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\dot{x}}{x}\right)}^{\mathrm{GW}}=\displaystyle \frac{2}{3}\displaystyle \frac{{\dot{P}}_{{\rm{b}},\mathrm{GW}}}{{P}_{{\rm{b}}}}=-6.31\times \,{10}^{-24}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 10 }$$

where we used the predicted ${\dot{P}}_{{\rm{b}},\mathrm{GW}}$ from the DDGR solution. This is more than eight orders of magnitude smaller than the measured value.

The third term, caused by aberration, is proportional to the geodetic precession rate for the pulsar, this is given by Barker & O'Connell (1975) as:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{geod}}={\left(\displaystyle \frac{2\pi }{{P}_{{\rm{b}}}}\right)}^{5/3}{T}_{\odot }^{2/3}\displaystyle \frac{1}{1-{e}^{2}}\displaystyle \frac{{M}_{c}(4{M}_{\mathrm{tot}}-{M}_{c})}{2{M}_{\mathrm{tot}}^{4/3}},\end{eqnarray} \tag{ 11 }$$

the result for PSR J2234+0611 is Ω^geod = 5.6 × 10⁻¹⁴ rad s⁻¹, i.e., the geodetic precession cycle has a length of 3.6 Myr. The aberration contribution to $\dot{x}$ is then given by (Damour & Taylor 1992):

$$\begin{eqnarray}&&\displaystyle \frac{d{\epsilon }_{A}}{{dt}}=\displaystyle \frac{P}{{P}_{b}}\displaystyle \frac{{{\rm{\Omega }}}^{\mathrm{geod}}}{\sqrt{1-{e}^{2}}}\displaystyle \frac{\cot \lambda \sin 2\eta +\cot i\cos \eta }{\sin \lambda },\end{eqnarray} \tag{ 12 }$$

where η and λ are the polar coordinates of the pulsar's spin. For PSR J2234+0611, the non-geometric factors (the first two fractions in the equation above) amount to −7.28 × 10⁻²³ s⁻¹, making this term about eight orders of magnitude smaller than the observed value.

The fourth term, $-\dot{D}/D=1.541\times {10}^{-18}\,{{\rm{s}}}^{-1}$, is three orders of magnitude smaller than the observed effect.

The fifth term can be derived from Equation (10), with the orbital variability given by Equation (3), from this we obtain ${(\dot{x}/x)}^{\dot{m}}=2.53\,\times \,{10}^{-21}\,{{\rm{s}}}^{-1}$, which is also extremely small.

Finally, the sixth term is due to changes in the orbital plane of the system from spin–orbit coupling; these are extremely small in such a wide system.

Since the first term in Equation (8) is the only measurable contribution to ${\dot{x}}_{\mathrm{obs}}$, we will now assume that the latter is described by Equation (9). Using that equation, we can combine ${\dot{x}}_{\mathrm{obs}}$ with the Shapiro delay to constrain the system geometry as shown in Figure 5. The orange lines show the $\cos i$ and Ω that are consistent with the measured $\dot{x}$ while the dashed blue lines show the $\cos i$ compatible with the assumed ς. These cross in four locations, listed in Table 3. They represent the four possible orbital orientations of the system according to the DDK and DDFWHE timing solutions.

Table 3.  Details for the Grid Regions

Region	$\cos i$	Ω	Best $\cos (i)$	Best Ω	Best Mtot	Min χ2
1	−0.92 to −0.52	340 to 540	−0.748	43.8	1.6512	5872.9
2	0.52 to 0.92	900 to 1100	0.748	94.4	1.6518	5881.7
3	0.52 to 0.92	2100 to 2300	0.716	220.6	1.7058	5926.0
4	−0.92 to −0.52	2700 to 2900	−0.704	278.4	1.7058	5929.6

Download table as:  ASCIITypeset image

For most pulsars where these constraints are available, we cannot eliminate the degeneracy implied by these four i-Ω solutions. However, if the binary system is relatively nearby and has a high timing precision, then apart from the secular variation of ω (${\dot{\omega }}_{\mathrm{obs}}$) and x (${\dot{x}}_{\mathrm{obs}}$) there are yearly cyclical variations in these parameters caused by Earth's orbit around the Sun (Kopeikin 1996). These are taken into account in the DDK model.

In Table 3, we can see that the quality of the local ${\chi }^{2}$ minima are clearly not identical, being significantly better for solution 1 (this is the DDK solution presented in Table 2); the latter solution is significantly better than either the DDGR or the DDFWHE solutions. A possible reason for this is that we have detected the yearly cyclical variations of x or ω or both. We quantify this statement in the next section.

Given all the correlations and caveats related to the different orbital models, and in order to better determine M_tot, M_c, M_p, i, Ω, their uncertainties and correlations, we have made a self-consistent χ² map of the Ω–$\cos i$–M_tot space using the DDK orbital solution with the assumption that GR is the correct theory of gravity. These parameters are chosen because they have, a priori, a constant probability density for randomly aligned orbits.

For each point in the grid of i, Ω, and M_tot values, we hold Ω and i fixed in the DDK model (from this it estimates all kinematic effects) and derive other relevant post-Keplerian parameters from the known mass function, i and M_tot (M2; OMDOT and GAMMA) using the GR equations. All these parameters are fixed inputs to the DDK model used to do the timing analysis for that grid point. The Einstein delay (GAMMA) must be calculated and used in the fit because, for wide binary pulsars like PSR J2234+0611, it is strongly correlated with $\dot{x}$ in the DDFWHE model and with Ω in the DDK model (see A. Ridolfi et al. 2018, in preparation). We then run Tempo, fitting for all other timing parameters not mentioned above, recording the value of χ² for each combination of Ω, $\cos i$, and M_tot.

Given the computational expense, our χ² map does not cover the full space; it consists instead of four disconnected regions around the four local χ² minima listed in Table 3; the cos i and Ω bounds sampled around these minima are also listed there. These variables are sampled with step sizes of 0.004 and 0.2 deg, respectively. For each grid section, we mapped the third variable, M_tot, from $1.641\,{M}_{\odot }$ to 1.731 M_⊙ with a step size of $0.0006\,{M}_{\odot }$. The quality of the fits in the regions outside these bounds are extremely low, for that reason those regions were not sampled.

The resulting 3D grids of χ² values are then used to calculate a three-dimensional probability density function (pdf) for ${\rm{\Omega }},\cos i,{M}_{\mathrm{tot}}$, as discussed by Splaver et al. (2002):

$$\begin{eqnarray}&&p({\rm{\Omega }},\cos i,{M}_{\mathrm{tot}})\propto {e}^{({\chi }_{\min }^{2}-{\chi }^{2})/2},\end{eqnarray} \tag{ 13 }$$

where ${\chi }_{\min }^{2}$ is the lowest χ² of the whole grid.

This 3D pdf is then projected onto two planes: the $\cos i\mbox{--}{M}_{c}$ plane and derived M_c–M_p plane (see the contours in the main panels in Figure 4) and the $\cos i$-Ω plane (see the contours in the main panel of Figure 5). It is also projected along three axes, $\cos i$ (top left side panels in Figures 4 and 5), M_c (and derived M_p, see top right and right side panels in Figure 4), and Ω (depicted in the right side panel of Figure 5).

The resulting 1D pdfs show that, as hinted by the ${\chi }^{2}$ values in Table 3, the probabilities for the four i-Ω solutions are far from identical. The solution with the lowest χ² (number 1) is preferred, with a total probability of 98.786%. The second most likely solution (number 2), has a total probability of 1.214%, it is still visible in the side panels of Figure 5 as a separate peak with very small amplitude. Solutions 3 and 4 have probabilities that are too small for our numerical precision, being thus definitively excluded. The discrimination between solutions 1 and 2, i.e., between the two possible values of i and Ω does not yet reach a statistical significance equivalent to 3σ, but they imply that the absolute orbital orientation of the system will be precisely known in the near future. Despite the fact that we cannot yet point out a single solution to equivalent 3σ significance, the exclusion of two of the solutions to high significance implies a significant detection of the annual orbital parallax.

The values derived for the quantities we set out to determine are: $i={138.7}_{-2.2}^{+2.5}\,\deg$ (68.27% C. L.), ${138.7}_{-4.2}^{+5.1}\,\deg$ (95.45% C. L.), ${\rm{\Omega }}={43.7}_{-2.2}^{+2.3}\,\deg$ (68.27% C. L.), and ${43.7}_{-4.4}^{+5.2}\,\deg$ (95.45% C. L.). In Figure 5, we see a fine correlation between i and Ω, which is a direct consequence of the precisely measured $\dot{x}$.

For the component masses, the situation is very clear: both solutions with measurable probability have the same total mass, ${M}_{\mathrm{tot}}={1.6518}_{-0.0035}^{+0.0033}\,{M}_{\odot }$ (68.27% C. L.), ${1.6518}_{-0.0070}^{+0.0066}\,{M}_{\odot }$ (95.45% C. L.). For the component masses we obtain ${M}_{c}={0.298}_{-0.012}^{+0.015}\,{M}_{\odot }$ (68.27% C. L.), ${0.298}_{-0.021}^{+0.034}\,{M}_{\odot }$ (95.45% C. L.), ${M}_{p}={1.353}_{-0.017}^{+0.014}\,{M}_{\odot }$ (68.27% C. L.), and ${1.353}_{-0.040}^{+0.025}\,{M}_{\odot }$ (95.45% C. L.). The fine Ω-i correlation also affects the mass measurements: for values of Ω closer to Θ_μ, the more face-on inclinations result in a more massive companion and a less massive pulsar.

The measurements made by the Bayesian analysis are in good agreement with the values inferred by the results in Section 3. For example, the total mass is well described by the $\dot{\omega }$ of the DDK solution, as it must since we used the latter model to map the masses assuming GR. The individual masses are well described by the intersection of the latter constraint with the h₃ of the DDFWHE solution. The constraints these impose on the range of inclinations plus the constraints imposed by the detection of $\dot{x}$ in the DDGR/DDFWHE solutions provide a good description of the range of Ω, plus its strong correlation with i near the best DDK solution.

Our analysis of PSR J2234+0611 is informed by previous work on two other eccentric binary systems, PSRs J1946+3417 and J1950+2414. Mass measurements of these three pulsars can be combined to constrain theories of formation for eccentric binary MSPs. The rotation-delayed accretion induced collapse (RD-AIC) hypothesis presented by Freire & Tauris (2014) for the formation of the eccentric MSPs has been excluded already by the mass measurement for PSR J1946+3417 presented in Barr et al. (2017): the mass of that pulsar (M_p = 1.828(22) M_⊙) is too large to have resulted from the collapse of a massive WD.

The RD-AIC theory could, in principle, generate larger MSP masses if we allow for differential rotation of the massive WD progenitor to the MSP: with differential rotation WDs can be much more massive than the ∼1.48 M_⊙ upper mass limit for rigidly rotating WDs. However, even in such a case the systems produced by RD-AIC would still have small peculiar velocities, otherwise the range of observed orbital eccentricities would not be as small as the observed range (see details in Freire & Tauris 2014). Such a possibility is difficult to reconcile with the large peculiar velocity measured for PSR J1946+3417 (in particular, its large vertical velocity relative to the Galaxy, see Barr et al. 2017) and the observed velocity of PSR J2234+0611 (Antoniadis et al. 2016a).

The large mass of PSR J1946+3417 is consistent with the hypothesis proposed by Jiang et al. (2015), which is also based on an instantaneous loss of binding energy of the more massive component. However, in this hypothesis the more massive component starts as a massive MSP. As it spins down, the centrifugal support is steadily reduced, causing a slow but steady increase in the central pressure with time, until a critical threshold is reached and the phase transition happens, presumably forming a quark star or some other type of exotic object that is still observable as a MSP. The sudden decrease in mass (owing to the larger binding energy of the new exotic remnant) results in the large orbital eccentricity. Other properties of PSR J1946+3417, like the large vertical velocity relative to the Galactic disk, are also consistent with this hypothesis, since the original system already shared the large velocity (relative to typical stars in the Galactic disk) typical of MSP-WD binaries. However, if there is a single pressure threshold for this phase transition, the masses observed for the MSPs in these eccentric systems should lie in a relatively narrow range (which is nevertheless finite because of differences in the spin periods, which would result in different NS masses for the same central pressure at which the phase transition occurs).

The masses measured for PSR J2234+0611 (1.35 M_⊙) and for PSR J1950+2414 (M_p = 1.495(24) M_⊙, W. W. Zhu et al. 2018, in preparation) are inconsistent with this hypothesis, since they are much smaller than the mass observed for PSR J1946+3417—clearly, a single pressure threshold for a phase transition does not provide a good description of these systems. Indeed, the observed MSP masses within this class appear to be as broad as those observed for the general MSP population (Antoniadis et al. 2016b; Özel & Freire 2016).

All measurements thus far are consistent with the expectations of the hypothesis proposed by Antoniadis (2014). This proposes that the orbital eccentricity is caused by material ejected from the companion due to unstable hydrogen shell burning. This hypothesis predicts that the MSPs in these systems should have a range of masses and Galactic velocities similar to those of the general MSP population; the observations are thus far consistent with this prediction.

Regarding the companions to the MSPs in these systems, all hypotheses advanced to date predict that they should be Helium white dwarfs with masses similar to what should be expected from the Tauris & Savonije (1999) relation. For PSRs J1946+3417, J1950+2414, and J2234+0611, the mass ranges predicted by this relation are 0.275–0.303, 0.268–0.296, and 0.281–0.310, respectively. In the case of PSR J2234+0611, our measured WD mass is in agreement with that prediction. For PSR J1946+3417, the companion mass M_c = 0.2556(19) M_⊙ is marginally consistent with this expectation, being lighter than expected. The companion of PSR J1950+2414 has a mass (${M}_{c}={0.280}_{-0.004}^{+0.006}\,{M}_{\odot }$) that is also well within the range expected by Tauris & Savonije (1999) for its orbital period. We note that within the context of the CB disk scenario, depending on the lifetime and mass of the disk, there could be a significant (up to ∼10%) reduction of the orbital separation. This effect would result in somewhat larger masses for a given orbital period, compared to the Tauris & Savonije (1999) relation—the opposite of what is observed for PSR J1946+3417.

The distance to PSR J2234+0611 is very well measured through the detection of timing parallax, ϖ = 1.03(4) mas. This corresponds to a distance d = 0.97(4) kpc. The uncertainty of 40 pc for J2234+0611 places it among the best measured pulsar distances.

The distance estimate will improve further in the near future. As the timing baseline T increases, the precision of ${\dot{P}}_{{\rm{b}},\mathrm{obs}}$ will also improve quickly. This will result in an additional precise distance estimate from the inversion of Equation (2), which will only be limited by knowledge of the Galactic potential. Measurements of this distance can be corroborated by VLBI campaigns. The component masses will also improve significantly, particularly the total mass; for the individual masses significant improvements will depend on advances in timing precision.

The precise distance and mass estimates presented here, together with the spectroscopic constraints on the WD atmospheric properties, transform the system into a laboratory for testing WD physics. As discussed in detail in Antoniadis et al. (2016a), the aforementioned measurements yield a radius estimate of ${R}_{\mathrm{WD}}={0.024}_{-0.002}^{+0.004}$ R_⊙ and a surface gravity of $\mathrm{log}g={7.11}_{-0.16}^{0.08}$ dex, both of which are model-independent. This is important for two reasons: first, PSR J2234+0611 is only the second system after PSR J1909−3744 for which independent atmospheric parameters can be obtained (Antoniadis et al. 2016a). Second, the surface temperature of T_eff ≃ 8600 K, obtained from atmospheric modeling, indicates that the WD envelope is convective. Spectroscopic 1D models for cool convective atmospheres are suspected to produce spurious results, but quantitative estimates and empirical corrections are difficult to obtain due to the lack of measurements. For both these reasons, PSR J2234+0611 becomes particularly important for calibrating atmospheric models. Currently, the precision of such tests is severely limited by the poor quality of the optical spectra, but could be improved significantly with further optical observations.

In addition, PSR J2234+0611 can also be used to test the predictions of WD mass–radius relations. One of the main remaining uncertainties in low-mass WD cooling models is the size of the hydrogen envelope that surrounds the degenerate Helium core. The latter can significantly affect the stellar radius, as well as the main energy source (residual hydrogen shell burning versus thermal cooling) and, consequently, the cooling age. Here again, our estimates are broadly consistent with the predictions for thin-envelope models, but a detailed test is limited by measurement uncertainties of the WD atmospheric parameters (see Antoniadis et al. 2016a, for details). For PSR J2234+0611, a future precision measurement of its envelope size is also important for probing its formation history, since the thin-shell instabilities on the proto-WD required for creating a CB disk are also expected to reduce significantly the size of the WD envelope (see Istrate et al. 2014, 2016; Antoniadis et al. 2016a, and references therein).

Last but not least, PSR J2234+0611 is within a few 100 K from the ZZ-Ceti instability strip for low-mass WDs. Kilic et al. (2018) recently reported on photometric observations of the system and found no pulsations. Consequently, the improved mass estimate reported in this work can further constrain the instability mechanism and the structure of WD envelopes (see Figure 5 in Kilic et al. 2018, and references therein for details).