X-Ray and Optical Studies of the Redback System PSR J2129–0429

We obtained all the available Swift X-ray Telescope (XRT) (Burrows et al. 2005) data sets taken of J2129 until the end of 2017 from the High Energy Astrophysics Science Archive Research Center (HEASARC) database, amounting to 17 observations in total, listed in Table 2, for a total exposure time of ∼24.7 ks. We used xselect to extract data products, a circular extraction region of 60'' radius for the source region, and an annulus of radius 180'' (excluding the source region of 60'') for the background region. We obtained a total of 106 photons after background subtraction (152 photon counts prior to subtraction), spread across 10 orbital phase bins. We find that the observations sample the orbit roughly evenly across phase bins, which differ in exposure time by a maximum of a factor of two.

The XMM-Newton observatory observed J2129 beginning on 2013 October 28 (Obs-ID 0725070101, PI Roberts) for ∼79 ks with a livetime of ∼74 ks for the MOS detectors and of ∼58 ks for the PN detector. Basic results from this observation were presented in Roberts et al. (2015) (see also Hui et al. 2015). Here, the aim is to see if the XMM-Newton spectra and light curves are consistent with what is seen by NuSTAR and with optical observations.

The observation was analyzed using the XMM-Newton Science Analysis System (SAS) version 15.0.0 (Gabriel et al. 2004), and HEASoft version 6.16. The observation was remarkably free from background flaring. Starting with the pipeline-processed event files, the EPIC event tasks were screened for flaring using the SAS task espfilt, and then barycentered with the task barycen. From these events light curves were created using an 18'' extraction region around the radio timing position of J2129 (a relatively small extraction region was chosen so as to exclude a nearby source, but should still contain ∼75% of the source photons). We created light curves in the 0.3–10 keV energy range, combining the counts from the two MOS detectors. A total of 49 orbital bins were used, with an entire orbit covered by 31 bins. We plot the light curves as a function of orbital phase in Figure 1, with phase 0 defined as the time of ascending node of the pulsar, making superior conjunction of the pulsar phase 0.25. Two successive superior conjunctions were covered by the observation.

Zoom In Zoom Out Reset image size

In all detectors, there appears to be a small dip in the light curve each time at pulsar superior conjunction. If this is a result of the pulsar surface being eclipsed, we would expect this dip to be more pronounced in the soft band than in the hard band (see Section 4). We attempted to verify this by creating folded light curves in the 0.3–1.5 keV and 1.5–10 keV bands, and we present those in Figure 2. These bands were chosen since 1.5 keV is approximately where the emission goes from being mostly thermal to mostly non-thermal at the lower count rates near superior conjunction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since one can generally expect a thermal X-ray component from the NS surface in MSP systems, we fit the spectra in XSPEC (Arnaud 1996) using an absorbed thermal plus power law model, where the thermal component is either a blackbody (BB) or standard XSPEC NS atmosphere (NSA) model. We find a best-fit absorption N_H = 1.8(2.9) × 10²⁰ cm⁻² and temperature kT = 0.21 keV ($\mathrm{log}T=6.1$ K) for the BB (NSA) thermal model for the entire observation, consistent with previously reported results. The light curves suggest that the second peak is somewhat harder than the first peak. We therefore extracted spectra at orbital phases 0.55–0.65 and 0.85–0.95 and fit them separately, keeping the absorption and thermal parameters fixed to the values derived from the entire spectra. Using either a BB or NSA model for the thermal emission, we see a marginally significant hardening of the second peak, with Γ = 1.05 ± 0.07 with a BB thermal component (Γ = 1.07 ± 0.07 with NSA; errors are 90% confidence) for the first peak and Γ = 0.93 ± 0.08 (Γ = 0.94 ± 0.08) for the second peak.

J2129 was observed by NuSTAR (Harrison et al. 2013) twice, between 2015 December 23 and 24 for an exposure of ∼44 ks, and between 2016 April 9 and 11 for an exposure of ∼79 ks (observation IDs 30101057002 and 30101057004, respectively). The data were cleaned, processed, and science products extracted using the standard NuSTAR pipeline NuSTARDAS (v1.5.1) that is part of the HEASoft (v6.17) package. We used the light-curve data produced by the pipeline for our light-curve analysis, choosing to select photons only up to an energy of 40 keV to reduce background noise, since we see little evidence for higher-energy photons in the energy spectrum. We bin the photons in 10 bins based on orbital phase, setting phase 0 to the time of ascending node, measure the cumulative detector exposure time, and keep count of the number of total photons in each bin to calculate count rates and errors, taking care to subtract background photons from our analysis. We use the analytical approximations for 68% confidence regions as given in Gehrels (1986). Because our source and background extraction regions are of unequal sizes, background subtraction inevitably results in fractional counts. Strictly speaking, Poisson statistics only apply to integer values, and as such we admit that the methods we employ are not entirely correct.

Our background regions were carefully selected to avoid any bright background sources that might interfere with the analysis, as we found that the size and location of the background region had a significant impact on our ability to produce reasonable fits to the X-ray spectra. We used a source region filter of radius 80'' centered on the source, and a background region filter of radius 90'' placed on the same detector chip as the source, as close as possible to the source while safely avoiding contamination from the source photons. Due to our sources lying very close to the chip-gaps during both of the observations, on each of the two detectors (FPMA and FPMB), we were unable to use larger extraction regions.

Due to potential systematic effects involved with low-count fitting and background modeling, we employed three different methods of data fitting (see van Dyk et al. 2001 for a summary of statistical issues). (1) We fit the spectra generated by running the pipeline, using the extraction regions detailed above for both the source and the background. (2) We constructed a simple model for the background comprising bremsstrahlung, power -law, and Gaussian line emission, with which we fit the pipeline-extracted background files. We found that our model agrees fairly well with the blank sky as seen by NuSTAR (Wik et al. 2014). We fit the source spectra simultaneously on top of the background model using an absorbed power-law model for the spectrum. (3) We used the cleaned Earth occultation data to account for the scattered photon noise in the data. Using the same extraction regions mentioned earlier, applied to the event files taken during the occultation observing mode, we manually extracted spectra using xselect. We then subtracted the occultation data from the science (pipeline-extracted) source and background data, adjusting for count rate variation between occultation and science observing modes. We then fit a simple absorbed power-law to the occultation-subtracted spectrum. We used the C-statistic for fitting when using methods 1 and 2, and we used χ² for fits using method 3.

Considering the large amount of uncertainty involved in fitting low-count data with substantial noise, we feel confident with our X-ray spectrum model as our best fits using the three different methods are in rough agreement. We used XSPEC version 12.9.0 for the spectral fitting presented in this paper.

The companion of J2129 was first identified through its UV counterpart in Swift data. The various Swift UVOT (Roming et al. 2005) observations used a variety of UV filters, including u, uvw1, uwv2, and uwm2. Exposure and sky coordinate images for these filters were obtained from NASA's HEASARC. These level 2 products have been previously processed through the standard pipeline for bias, bad pixel adjustment, and coordinate transformation. Packages from HEASoft's Swift software were used for data processing. Specifically, uvotdetect was used to detect the source on the level 2 products, using the provided calibration database for the SWIFT UVOT instrument on HEASARC. Photometry was done with the uvotsource package using a standard sky substraction region and the suggested source region by uvotdetect.

For uw1, there were nine exposures totaling an effective exposure of 7448 s. Only one of these observations was made in 2016, while the rest were all taken for the same epoch in 2013. U exposures were more distributed, with a couple of observations each year between 2013 and 2016. While the light curve was constrained on these filters, systematic and statistical errors inherent to the processing through the UVOT pipeline and associated HEASoft mission-specific packages precluded any long-term trend analysis. As for uvm2 and uw2, the source was too faint on these filters for analysis or light-curve constraints.

Fitting parameters for uw1 return an average AB magnitude of 21.30 ± 0.07. The fit suggests the light curve is dominated by ellipsoidal modulation, in correspondence with g', r', i' data (see Figure 3). In comparison, the standard processing for XMM-Newton Optical Monitor (OM) uvw1 data for an observation with a total exposure of 67,400 s in 2014 leads to an average AB magnitude of 20.70 ± 0.02, a difference of 0.6 from the Swift value. However, a systematic difference of 0.83 ± 0.23 is found with other sources in the field of view between the XMM-Newton OM and Swift sources. The XMM-Newton value is significantly brighter than expected from an extrapolation of g', r', i', and z' values using the Bellm et al. (2016) fits for the temperature. While this could be interpreted as an UV excess during the XMM-Newton observation, the consistency of Swift uw1 values before and after the XMM-Newton observation with the extrapolation from longer wavelengths suggests the XMM-Newton OM value is in error.

For this study, J2129 was monitored across six epochs, starting on 2014 November and ending on 2017 October. These observations were made with 1 and 2 m telescopes from the LCOGT network, using exposures which were typically 200 s in duration. Data were obtained at each epoch using the standard i', r', and g' SDSS filters. Images obtained in 2016 and 2017 were processed using the LCOGT BANZAI reduction pipeline while earlier observations used the ORAC pipeline (Jenness et al. 2013). Table 3 summarizes the start and end dates for each epoch, and the number of measurements taken during each. Some measurements had to be discarded due to poor image quality. A total of 10 exposures were discarded, accounting for about 4% of all observations taken. A few were discarded for shutter issues, where only a fraction of the field of view was observed; in these cases, the source would be within the unobserved region. A few other exposures were discarded due to poor image quality where the source was undetectable with the g' filter or the measured errors were so large that the measurement gave no meaningful constraints on the observed trends. A few others were discarded due to visible patterns present in the image, either as a result of lens condensation or bad flat-field, that interfered directly with the area of extraction and affected the integrity of the image.

Image Reduction and Analysis Facility v2.16 (IRAF; Tody 1986) and SAOImage DS9 v7.5 (Joye & Mandel 2003) were used to carry out aperture photometry analysis on the reduced images. Aperture radii were optimized for every image with respect to the source only. In order to measure trends in the average orbital flux, the procedure outlined in Bellm et al. (2016) was followed to fit the light curves with a simple model of sines and cosines. Equation (1) shows this model where ϕ is the orbital phase and coefficients a₁ through a₅ account for the shape of the light curve.

$$\begin{eqnarray}m(\phi ) & = & {m}_{0}+{a}_{1}\cos (2\pi (\phi -0.25))\\ & & +{a}_{2}\cos (4\pi (\phi -0.25))+{a}_{3}\cos (6\pi (\phi -0.25))\\ & & +{a}_{4}\sin (2\pi (\phi -0.25))+{a}_{5}\sin (4\pi (\phi -0.25)).\end{eqnarray} \tag{ 1 }$$

The Swift-XRT observations were taken over the course of approximately six years. Although the observations are short and collectively do not sample the entire orbit, the light curve we are able to produce is consistent with the XMM-Newton observation in terms of orbital variation as a function of phase, as can be seen in Figure 4, indicating that the prominent double-peaked feature noted by Roberts et al. (2015) persists over several years. We use the online PIMMS tool (v4.8d) (Mukai 1993) provided by HEASARC to convert the Swift count-rates into their XMM-Newton equivalent, setting the Galactic N_H = 4 × 10²⁰ cm⁻², which we find is the best fit to our data, and the source model to a power law of Γ = 1.0.

Zoom In Zoom Out Reset image size

The NuSTAR observations exhibit the same double-peaked phenomenon, as shown in Figure 5. This double-peaked feature has been hinted at in other spiders (Bogdanov et al. 2005, 2014a, 2014b; Gentile et al. 2014; Romani 2015; Salvetti et al. 2015); however, J2129's features are the most prominent seen so far.

Zoom In Zoom Out Reset image size

We find that the X-ray spectrum extends to around ∼40 keV, is well fit with a single absorbed power law with a very hard photon index of Γ ∼ 1.05–1.3, and shows no evidence for breaking. It overlaps well with the XMM-Newton data in the soft regime, as can be seen in Figure 6. We fit the data with an absorbed power law, with either a blackbody or an NS atmosphere component (Romani 1987; Zavlin et al. 1996). The best-fit parameters we find for the spectrum using the three methods detailed in Section 2.3 are shown in Table 4. There does not appear to be a significant change in the overall spectrum between the two NuSTAR observations, so we only show parameters for simultaneous fitting of both.

Zoom In Zoom Out Reset image size

Table 4.  Best-fit Parameters for the Full X-Ray Spectrum (XMM-Newton and NuSTAR) of J2129

	Method 1	Method 2	Method 3
NH (×1020 cm−2)	5.38 ± 3.2	7.60 ± 5.3	6.05 ± 5.3
$\mathrm{log}{T}_{\mathrm{eff}}({\rm{K}})$	6.07 ± 0.107	5.89 ± 0.115	5.85 ± 0.164
Γ	1.03 ± 0.102	1.18 ± 0.063	1.35 ± 0.060
reduced χ2 (DOF)	1.27 (250)	2.75 (284)	1.49 (250)
NH (×1020 cm−2)	3.59 ± 2.2	2.93 ± 3.1	4.88 ± 3.8
kT (keV)	0.180 ± 0.028	0.170 ± 0.040	0.132 ± 0.035
Γ	1.09 ± 0.056	1.16 ± 0.058	1.36 ± 0.053
reduced χ2 (DOF)	1.28 (250)	2.75 (284)	1.49 (250)

Note. We try fitting to an absorbed power law with either an NSA or a BB component. For the NSA model we use the canonical values M_ns = 1.4 M_⊙ and R_ns = 10 km. The values do not change significantly if we instead use M_ns = 1.7 M_⊙ and R_ns = 12 km.

Download table as:  ASCIITypeset image

We also note that emission from the second peak appears to have significantly hardened between the first and second observations, which are less than four months apart, as can be seen in the hardness ratios in the bottom left plot of Figure 7. To compute hardness ratios, we follow the same procedure used to produce our light curves, and then divide the soft by the hard count rates. Due to the statistical biases this method introduces, we additionally use the two Bayesian-based methods developed by Park et al. (2006) to produce better estimates of the hardness ratio. Their method employs Bayesian analysis on a Poisson posterior distribution, integrating using either a Markov chain Monte Carlo (MCMC) technique or numerically through Gaussian quadrature (see also van Dyk et al. 2001). Hardness ratios computed using both those methods, using the BEHR program (Park et al. 2006), are shown in the plots on the right in Figure 7. We used the mode of each distribution as data points, seeing as it provides the most robust hardness ratio estimate.

Zoom In Zoom Out Reset image size

Fitting only NuSTAR photons that fall between phase 0.5 and 0.75, to an absorbed power-law model, we find a photon index 1.10 ± 0.13 (fixed N_H = 4.4 × 10²⁰ cm⁻²), with reduced χ² = 1.05. We find that the spectrum of the first peak does not change significantly between observations. However, the spectra we obtain for the second peak vary significantly between the two observations, as reflected in Figure 7. The best-fit parameters for the spectrum of the second peak (0.85–0.95 phase) are shown in Table 6 for our three analysis methods. The variation in spectral index is considerably less obvious when we instead use the photons between phases 0.75–1.0, indicating that the mechanism that causes emission at the peaks is highly geometrically constrained.

In comparison, the spectrum of the first peak (0.55–0.65 phase) shows little evidence for variation, with a best-fit photon index of Γ ∼ 0.9 using the first two methods of analysis (Table 6). While employing the occultation subtraction method, we found difficulty in producing results that agree between the two detectors. For this reason we do not use this method for studying the first peak.

Simultaneously fitting the individual peaks (0.55–0.65 and 0.85–1.0 in phase for the first and second peaks respectively) with the XMM-Newton data, assuming an absorbed power-law with a blackbody component model, it is clear that the hard emission (Γ < 1.2) is a persistent phenomenon. We find a best-fit photon index of Γ = 0.83 ± 0.132 (χ² = 0.607 with 59 DOF) corresponding to the first peak, a value consistent with the NuSTAR spectrum at that phase. Emission from the second peak, however, appears to be somewhat erratic. We find a best-fit photon index of Γ = 1.08 ± 0.126 (χ² = 1.460 with 40 DOF), almost halfway between the spectra from the first and second NuSTAR observations.

The orbital variations in the LCOGT light curves were dominated by the a₂ and a₅ sinusoidal components, which roughly account for ellipsoidal modulation and residual distortions, respectively (Bellm et al. 2016). Figure 8 shows the fitted light curves for the i' filter, spanning all epochs from 2014 November to 2017 October. From 2014 November to 2016 June a dimming was observed, i.e., an increase in average apparent magnitude, whereas from 2016 June to 2017 October the apparent magnitude of the system was observed to decrease. Figure 9 shows these trends in all filters. The magnitude values we use for these analyses are the fit m₀ values. In Figure 10, we plot the first and last epoch for each filter. There is a clear decrease in the g' intensity and an increase in the i' intensity, indicating cooling. There appears to a systematic shift in the light curve, something that has been seen in other systems (Li et al. 2014); however, in this case there are not enough data to conclusively say that the shift is real and not due to fluctuations in the data. In Figure 11 we plot a color–color diagram which suggests continuous cooling since 2014.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The two-peaked orbital light curve of J2129 provides a puzzle for the shock geometry. As pointed out by Roberts et al. (2015), if the peaks are due to Doppler boosting along the shock front, then their orbital phases at ∼0.6 and ∼0.9 would suggest the shock is directed away from the companion rather than wrapping around it, as would be expected if a powerful pulsar wind is interacting with a much weaker stellar wind. However, with only one observation of the orbit, the peaks could have represented random flaring rather than be associated with particular orbital phases. Here we have shown that, in two observations with NuSTAR and the average light-curve of the Swift observations, these peaks are persistent phenomena and so must reflect something intrinsic and relatively stable about the shock geometry. The large variance seen also suggests that the system must have been observed at a high inclination angle. This is consistent with the photometric results of the companion and the interpretation of the X-ray dips seen at orbital phase 0.25 as an eclipse of the pulsar's surface emission by a nearly Roche lobe filling companion and a high inclination angle of the system (Bellm et al. 2016). The greater X-ray dips in the lower energy band (0.3–3 keV, see Figure 2) are due to the obscuration of thermal emission from the pulsar surface.

In addition, if the shock region wrapped closely around the companion, then the size of the emission region is roughly that of the companion. The optical observations reported in Bellm et al. (2016) suggest that the companion fills ∼95% of the Roche lobe in the period from 2012 to 2014, from which we calculate that it takes up about 1% of the pulsar sky. From this estimate, assuming an isotropic pulsar wind, then we would estimate a spin-down conversion efficiency of about 20% into 0.3–40 keV X-rays. If the wind is fully confined to the equator, we would estimate the companion to intercept ∼7% of the wind, implying instead a conversion efficiency of ∼2%–3%. These very high efficiencies suggest that the pulsar wind may be interacting with something much more extended than the companion star itself.

Two models have been proposed to explain how heating from an intra-binary shock can produce a double-peaked light curve such as that of J2129. Romani & Sanchez (2016, hereafter RS16) construct a model where the pulsar wind flux is concentrated equatorially, modeling it with sin⁴θ dependence, and shocks against an isotropic companion wind. They find that the ratio of the momentum flux from the two winds is

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{\dot{M}}_{w}{v}_{w}}{\dot{E}}\end{eqnarray} \tag{ 2 }$$

where ${\dot{M}}_{w}$ and v_w are the companion wind mass-loss rate and velocity respectively, and largely determines the geometry of the intra-binary shock. According to their model, J2129 fits nicely with the case where η > 1, although η may have to be as high as 10 or greater to be in agreement with the model.

Another model, presented by Wadiasingh et al. (2017, hereafter W17), finds similar conclusions. Their model follows from the thin-shell approximation for the collision of two isotropic winds (Canto et al. 1996), with a bow shock forming around the pulsar. Doppler-boosting along the shock would correspond to the orbital modulation we observe. Constraining the X-ray model using the fraction of radio eclipses observed in the system, the synthetic light curves they produce are dependent mostly on the shock radius R (from the MSP) and the inclination angle i, with greater modulation as the bulk velocity at the shock β increases. Based on their model, given J2129's inclination angle of i ∼ 80°, we would expect R ≤ 0.1 r₀, where r₀ is the orbital separation between the pulsar and the companion, and β > 0.9. Both models imply a shock dominated by the companion's wind, η ≫ 1 in RS16 and R ≤ 0.1 r₀ in W17, a scenario that we find difficult to justify on physical grounds. Wadiasingh et al. (2015) argue that an intra-binary shock surrounding the pulsar in a transitional redback system is not entirely unexpected as it may have recently been in an accretion-powered state. In such a case the wind may shock against the remnants of the accretion disk, although Bellm et al. (2016) find no evidence to support such a scenario. The companion "wind" directed toward the pulsar may also be enhanced by some material overflowing its Roche lobe and falling into the pulsar's gravity well (q-RLOF) with the formation of a disk being inhibited by the pulsar wind. Whether or not this would result in a fairly stable shock front of the required shape would need to be investigated further. Chernyakova et al. (2017) find flux and spectral variations on hour-long timescales in the X-ray binary system LSI +61 303, which they attribute to the the stellar wind probing different structures as it moves around the system. Such a scenario may be able to account for some of the short-term modulations we see if any evidence for an accretion disk can be found.

Estimating the magnetic field at the intra-binary shock based on these two models, we find a lower limit of 13 G based on a shock distance of ∼0.3 r₀ (RS16). However W17 predict that the shock would occur within 800 light-cylinder radii of the pulsar, with a magnetic field of ≥42 G at the shock. Sanchez & Romani (2017) derive an expression for the shock radius as a function of the momentum flux ratio, where for η = 15 the shock is as close to the pulsar as 0.01 r₀, whereas for η = 3 it occurs at 0.23 r₀. According to these values the magnetic field strength at the shock may be 18 G  < B < 420 G.

The second peak falls off more slowly than the first peak rises, which may be due to instabilities in the trailing wind at the shock, as is seen most clearly in Figures 1 and 2. The simulations of Bosch-Ramon et al. (2012) show turbulent mixing at the contact discontinuity between the pulsar and the companion's winds (see also Bosch-Ramon et al. 2015). The trailing edge appears to be more turbulent than the leading side, and is likely much more sensitive to the geometry.

The photon index we find, Γ ∼ 1.1–1.2 (and the even harder photon indices at the peaks), is inconsistent with Doppler-boosted emission from shock-accelerated particles, which tends to be much softer with Γ ≥ 1.5 (Sironi et al. 2015). A much more likely mechanism for the production of such hard X-ray power-law emission is driven magnetic reconnection, for which the expected photon index can be as hard as Γ ∼ 1 (Zenitani & Hoshino 2001; Larrabee et al. 2003). Some numerical calculations (Kagan et al. 2016; Werner et al. 2016) suggest that the photon index is dependent on the plasma magnetization parameter, defined as the ratio of Poynting to kinetic energy flux,

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{{B}^{2}}{4\pi \gamma {{nm}}_{e}{c}^{2}},\end{eqnarray} \tag{ 3 }$$

where n is the plasma number density and γ is the average bulk Lorentz factor, such that Γ hardens from ∼1.4 (for σ = 3) to ∼1.2 (for σ = 1000) as the magnetization increases. Magnetic reconnection is believed to contribute to particle acceleration in the pulsar wind, as was first discussed by Coroniti (1990) and Michel (1994). In this scenario, the "striped" pulsar wind forms a shock against the companion's wind, which drives reconnection in the alternatingly magnetized stripes of plasma in the pulsar wind. It is thought that the magnetically striped wind annihilates fully through magnetic reconnection at the termination shock of a pulsar wind nebula (PWN), a possible solution to the "σ problem" in pulsar winds (Lyubarsky 2003). Particle-in-cell 2D and 3D simulations come to similar conclusions regarding the hardness of the power law (Lyubarsky & Liverts 2008; Sironi & Spitkovsky 2011; Werner & Uzdensky 2017). This makes it a much more viable acceleration mechanism for the observed emission, considering the hardness of J2129's spectrum.

We may be able to constrain the cascade pair production multiplicity factor κ, if we interpret the variation in the spectrum of the second peak as resulting from partial, rather than full, magnetic dissipation at the shock, where the pulsar wind still maintains alternating sheets of magnetized plasma. Pétri & Lyubarsky (2007, hereafter PL07) argue similarly in an attempt to explain the modulation observed in the double pulsar system PSR J0737−3039 (McLaughlin et al. 2004). We note that 2D and 3D simulations (e.g., Sironi & Spitkovsky 2011) have found results contrary to those of PL07, in that they conclude that full dissipation always occurs. Although it is widely accepted that the striped wind is completely annihilated by the time it reaches the termination shock in a PWN, that may not necessarily be the case for intra-binary shocks seen in redbacks or BWs (Hoshino & Lyubarsky 2012). If we follow the condition PL07 find for full magnetic dissipation at the shock:

$$\begin{eqnarray}&&\displaystyle \frac{R}{{R}_{\mathrm{lc}}}\geqslant \displaystyle \frac{2\pi }{3}\kappa \end{eqnarray} \tag{ 4 }$$

where the light cylinder radius, corresponding to J2129, and the separation between the pulsar and the companion's surface are R_lc = 363 km and R ∼ 2.1 × 10⁶ km respectively, we estimate a lower bound for the pair multiplicity factor, κ ∼ 2760, assuming the shock occurs exactly at the companion's surface where the pulsar's magnetic field would be ∼6 G. Lyubarsky & Kirk (2001) calculate a theoretical maximum radius for the striped wind, which in this case would be much greater than ∼100 R, hence the partial reconnection interpretation may be plausible. Although the value of the pair multiplicity is uncertain, our calculations fall safely within κ ∼ 10³–10⁴ believed to be the pair multiplicity for MSPs (Harding & Muslimov 2011). If the intra-binary shock lies very close to the pulsar as predicted by RS16 and W17, however, the pair multiplicity may have to be several orders of magnitude greater. Constraints on the pair multiplicity factor may help inform theories on the nature of pulsar radio emission. This interpretation, however, leaves open the question of why the spectrum at the first peak is left undisturbed.

We calculate estimates of the magnetization, σ, that agree well with the simulations of Werner et al. (2016) and Kagan et al. (2016). High magnetization (σ ≫ 1) in the wind is a requirement for relativistic magnetic reconnection. Assuming a bulk plasma with Lorentz factor γ = 100 (Lyutikov et al. 1999), and using the plasma density approximation scaled from the Goldreich–Julian density (Goldreich & Julian 1969; Lyubarsky & Kirk 2001),

$$\begin{eqnarray}&&n=\displaystyle \frac{\kappa {B}_{\mathrm{lc}}{\rm{\Omega }}}{2\pi {ec}}{\left(\displaystyle \frac{{R}_{\mathrm{lc}}}{R}\right)}^{2}\end{eqnarray} \tag{ 5 }$$

where B_lc is the magnetic field strength at the light cylinder, e is the electron charge, Ω is the angular velocity of the pulsar, and κ = 2760 from the lower limit calculated above. We find that for a shock at a distance of 0.1 r₀ with B = 42 G, as would be the case in W17, we find a magnetization factor σ ≲ 1300. For a shock at 0.3 r₀ and B = 13 G, as in RS16, we calculate σ ≲ 1100. If the shock instead occurs near the surface of the companion we find σ ≲ 1400, values that are all consistent with what is needed for our measured photon index and the simulations of Werner et al. (2016) and Kagan et al. (2016). Interestingly, Lyubarsky (2003) finds an estimate of the plasma magnetization given based on the pulsar spin-down power, which gives roughly double the estimates we find above, σ ∼ 2300 at 0.1 r₀, for J2129.

Based on the flux estimates, we calculate J2129's power-law X-ray luminosity in the 0.3–40 keV band to be around (5–9) × 10³¹ $(\tfrac{d}{1800\,\mathrm{pc}})$ erg⁻¹ s⁻¹, depending on the analysis method and model used. We choose the optically derived distance rather than the dispersion measure distance as it is more constrained (see Table 5). This corresponds to roughly 0.12% to 0.23% of the spin-down power, $\dot{E}=3.9\,\times$ 10³⁴ erg s⁻¹ (Roberts et al. 2014). Given this energy budget, we would expect the power-law emission to cut off or begin to curve below 10 MeV.

Sironi & Spitkovsky (2014) calculate a maximum photon energy that depends on the magnetization and the spectral index, γ_max ∼ [σ(3−2Γ)/2(Γ−1)]^1/(3−2Γ), which corresponds to synchrotron emission, in the case of J2129, of about ∼12–280 keV for a photon index Γ ∼ 1.2. For Γ ∼ 1.1 the maximum photon energy is ∼0.3–7.4 keV. For greater Γ, the emission would fall well within the range detectable by the Fermi Large Area Telescope. For an emission cut-off ∼40 keV, as our data suggest, we would expect Γ ∼ 1.17, consistent with the observed photon index. We adopt a magnetization of σ = 1000 for these calculations.

Recent work by Sanchez & Romani (2017) suggests that the intra-binary shock reprocesses the pulsar wind, with particles being ducted to the companion's magnetic poles, resulting in additional heating when compared to direct heating, and leads to better agreement with BW light curves. The model seems to apply to redbacks as well. Such a model would necessitate the presence of a strong companion magnetic field. In this case, the magnetic field may play a role in the position and shape of the shock, as well as the inferred magnetic field in the emitting region. A strong companion magnetic field may also contribute toward observed stochastic orbital variations in tight binaries (Applegate 1992; Applegate & Shaham 1994).

Wolff et al. (2009) suggest that magnetic cycling in the companion may account for orbital variation in the binary period, a plausible scenario in the case of J2129 considering the rapid rotation of its companion. It may even be possible, depending on the degree of magnetic cycling in the companion, to see changes in the X-ray spectra that correlate with the magnetic cycles of the star.

It has been shown that heating plays a role in the dynamics of redback systems, although the contribution is much stronger for BWs. Optical studies by Breton et al. (2013) show that eclipsing MSP companions exhibit ∼10%–30% efficiency in recycling the incoming pulsar wind energy flux for surface heating. However in the case of J2129 the companion is very mildly irradiated, with a irradiation efficiency η ∼ 2%. The X-ray flux we observe would at most contribute to 10% of the heating. However, the companion has been observed to dim and then brighten up, but has an overall cooling trend (Bellm et al. 2016, this paper). This suggests some alternative heating source in the past, perhaps due to increased X-rays during a recent accretion episode or tidal heating.

Using the temperature scale of Pinsonneault et al. (2013) to convert colors to effective temperatures, and correcting for extinction, (A_v = 0.14 ± 0.09, Bellm et al. 2016), we estimate the change in the temperature of the companion over the years of observation. The values in parentheses are errors corresponding to measurement errors on the color. In the g − r colors, the temperature went down from 5222 (5199–5245) K to 5104 (5082–5127) K within the three years between 2014 November and 2017 October. The temperatures derived from the g − i values dropped from 5104 (5025–5157) K to 5001 (4971–5031) K. We convert the changes in magnitudes in the g − r and g − i colors over the three years to a ratio of intensities (I₁₄/I₁₇ = 10^−Δm/2.5) and, assuming that the ratio of the companion radii squared is directly proportional to the ratio of intensities, we can estimate the change in the companion's size by accounting for the blackbody temperatures of the measured colors, from which we find that 19 R_Nov14 ≈ 20 R_Oct17. Considering that J2129 was, at the time of Bellm et al. (2016), 95% Roche lobe filling, this estimate implies that the system is currently, or almost, filling its Roche lobe and may enter an accreting state very soon. We plan to monitor the system in radio, optical, and X-rays in anticipation of such an event.