# The bow shock and kinematics of PSR J1959+2048

Roger W. Romani, Adam Deller, Lucas Guillemot, Hao Ding, Martijn de Vries, Chase Parker, Robert T. Zavala, Aurélien Chalumeau, Ismaël Cognard

► To cite this version:

Roger W. Romani, Adam Deller, Lucas Guillemot, Hao Ding, Martijn de Vries, et al.. The bow shock and kinematics of PSR J1959+2048. The Astrophysical Journal, 2022, 930, pp.101. 10.3847/1538-4357/ac6263. insu-03672003

HAL Id: insu-03672003

https://insu.hal.science/insu-03672003

Submitted on 4 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# The Bow Shock and Kinematics of PSR J1959+2048

Roger W. Romani¹, Adam Deller², Lucas Guillemot³,⁴, Hao Ding², Martijn de Vries¹, Chase Parker¹, Robert T. Zavala⁵, Aurélien Chalumeau³,⁴,⁶, and Ismaël Cognard³,⁴,⁵

¹Department of Physics/KIPAC, Stanford University, Stanford, CA 94305-4060, USA; rwr@astro.stanford.edu

²Center for Astrophysics & Computing, Swinburne University of Technology, Mail Number H74, PO box 218, Hawthorn, VIC 3122, Australia

³Laboratoire de Physique et Chimie de l'Environnement et de l'Espace LPC2E UMR7328, Université d'Orléans, CNRS, F-45071 Orléans, France

⁴Station de Radioastronomie de Nançay, Observatoire de Paris, PSL University, CNRS, Université d'Orléans, F-18330 Nançay, France

⁵U.S. Naval Observatory, Flagstaff Station, 10391 West Naval Observatory Road, Flagstaff, AZ 86005-8521, USA

⁶Université de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France

Received 2022 February 28; revised 2022 March 28; accepted 2022 March 28; published 2022 May 9

## Abstract

Pulsar Hα bow shocks provide rare opportunities to constrain the energetics and kinematics of the relativistic pulsar wind. We have acquired optical imaging and integral field unit spectroscopy of the bow shock of the millisecond pulsar PSR J1959+2048, measuring the shock symmetry axis at a position angle = 213.2 ± 0.2 and showing that this slow nonradiative shock has a broad-to-narrow line component ratio $I_b/I_n = 4$. The data show that the pulsar's velocity lies 2.2 out of the plane of the sky. Coupled with a new fit for its timing proper motion, giving $\mu_{tot} = 30.05$ mas yr$^{-1}$ and a Very Long Baseline Array (VLBA) interferometric parallax measurement giving $d = 2.57^{+1.84}_{-0.77}$ kpc (90% range), we have unusually complete information on the pulsar kinematics. The bow shock constraints on the wind momentum flux imply that, at the best-fit parallax distance, the pulsar moment of inertia must be very large and/or the Hα efficiency at its modest shock velocity must be very high.

Unified Astronomy Thesaurus concepts: Pulsar wind nebulae (2215)

## 1. Introduction

The millisecond pulsar PSR J1959+2048 (hereafter J1959) was the first “black widow” companion-evaporating pulsar to be discovered (Fruchter et al. 1988), with a spin period 1.607 ms and an observed period derivative $\dot{P} = 1.685 \times 10^{-20}$ s s$^{-1}$. The short spin period and low inferred dipole field suggest appreciable recycling via mass accretion. Indeed it has been proposed that the pulsar has accreted sufficient mass to reach $2.4M_\odot$ (van Kerkwijk et al. 2011), although more recent observations find γ-ray (Clark et al. 2022) and X-ray (Kandel et al. 2021) eclipses, implying an edge-on view of the system and a more modest $M_{PSR} \approx 1.8M_\odot$. The pulsar mass, and hence moment of inertia, affect the spindown energy loss rate. In turn this affects the possibility that J1959 may complete evaporation of its companion before spindown, thus providing a path to the otherwise puzzling isolated millisecond pulsars.

The spindown power can be probed by study of the pulsar's relativistic wind. Most commonly these winds are observed as X-ray synchrotron nebulae. These are continuum emission sources and as such cannot be used to measure the full spindown energetics. However, J1959 is also the first known example of a pulsar wind nebula enclosed in an Hα bow shock (Kulkarni & Hester 1988). These rare objects (only eight others are known; Brownsberger & Romani 2014, hereafter BR14) provide additional opportunities to study the relativistic pulsar wind as it shocks against the interstellar gas. In particular, the Hα line kinematics provide information on the pulsar's radial velocity and the wind's momentum flux, as well as the state of the surrounding interstellar medium. We describe here new optical and radio observations that provide unusually complete

information on the 3D pulsar velocity and its relativistic wind. We find that the pulsar spindown power may require a remarkably large moment of inertia, although uncertainties in the system distance and the Hα shock efficiency temper this conclusion.

## 2. Observations and Data Reduction

### 2.1. New SAMI Hα Image

J1959 was observed on 2018 August 6 at the SOAR 4.1 m telescope using the SAMI camera on the SOAR Adaptive Module, which uses a UV laser to provide ground-layer adaptive optics (AO) corrections (Tokovinin et al. 2016). We obtained 3 × 1200 s in Hα ($\lambda = 6569.9 \pm 9.3Å$) and added continuum images with 2 × 30 s in the Sloan Digital Sky Survey (SDSS) griz filter bands. Images were binned 2 × 2 for a scale of 0.09/binned pixel and the AO system was used in all exposures. Conditions were fair with the individual continuum frames delivering (corrected) image FWHM ≈ 0.8 (i) to ≈ 1.1 (g) and the integrated Hα image FWHM = 0.7. The data were processed with the SAMI pipeline for basic calibration and cosmic ray cleaning. To improve the coordinate system we matched a set of stars surrounding the pulsar to GAIA positions. Camera distortion left a scatter in the absolute coordinates of 0.1, but allowed definition of the field orientation to 0.02. Calibration exposures of planetary nebulae PN G018.1+20.1 (Na 1) and PN G043.3+11.6 (M3-27) were used to establish a flux scale (Frew et al. 2013).

For Figure 1 we combine the Hα frames (red) with r (green) and g (blue). The binary was at orbital phase $\phi_B \approx 0.05-0.2$, nearing optical minimum, so the pulsar optical counterpart (the leading, southwest member of the pair near the bow shock focus is relatively faint in the combined images. At maximum it has magnitude $r = 19.6$, similar to that of the star to the northeast, 1.5 away in 2018.

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** SOAR SAMI image of the PSR B1959+2048 Hα bow shock. This color image combines Hα (red) with SDSS r (green) and SDSS g (blue). Note that the pulsar companion (leading star of the pair near the bow shock focus) is faint, near optical minimum. Note also the scalloping at the Hα limb, especially on the nebula's western edge. Image spans 2.4 × 2.2.</figcaption></figure>

## 2.2. KOALA Hα Data Cube

We observed J1959 on 2014 October 19 and 20 with the KOALA Integral Field Unit (IFU) on the Anglo-Australian Telescope (AAT) feeding the AAOmega spectrograph (Ellis et al. 2012). With the AAT seeing of 1.5–2.5 we could not access the finer details of the bow shock, but KOALA provides 1000 fibers at 25 pitch covering $50'' \times 27''$, which allowed us to take spectra across the face of the nebula. We set up the AAOmega double spectrograph with the 2000R grating in the red channel, covering 6273–6737 Å at 0.48 Å (22 km s$^{-1}$) resolution. This sufficed to resolve the velocity structure of the Hα nebula. The blue arm employed the 1500V grating, with ~3920/5420 Å providing access to the higher Balmer lines. However, with the lower ~58 km s$^{-1}$ resolution and low S/N, these data did not constrain the line kinematics, so we focus here on the Hα results.

We obtained 6 × 1800 s exposures on J1959, with the pointing dithered several arcseconds between integrations. These were interspersed with flat field exposures and spectrophotometric flux standards. All data were processed with the instrument's 2dfdr software. This included spatial flattening via the (non-Hα) sky lines. The six individual frames were combined into a final flux-calibrated position–velocity data cube. Pulsar bow shocks are nonradiative, so the optical emission line flux is dominated by Balmer lines, with Hα by far the strongest. To isolate this Hα emission, we generated a median filtered continuum image from 6530 to 6600 Å and subtracted this from each velocity channel. This largely removes the field stars, although of course the Hα photospheric absorption results in negative stellar features (Figure 2). These are most prominent for the hotter stars; the differing radial velocities for these field stars are apparent with the absorption maximizing in different channels.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** KOALA IFU Hα position–velocity cube slices, with wavelengths in Å in blue above each slice and velocity in magenta below. The oval region in the central velocity slice provides the line spectrum shown in Figure 3. Dark spots represent the Hα absorption line of the field stars, whose varied depth and velocity prevent perfect continuum subtraction.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** KOALA IFU Hα profiles. Right: synthesized long-slit spectrum along the velocity axis. Note the oversubtraction of a field star with strong Hα at $x \approx -18$. Middle: synthesized 2D spectrum for an effective slit transverse to the velocity axis. Right: velocity profile from the elliptical region $\sim 3r_a$ behind the pulsar (see Figure 2), with narrow and broad component fits.</figcaption></figure>

Nonradiative shocks should have two velocity components: a narrow line at the ambient ISM velocity, formed by excitation of neutral H drifting into the heated postshock gas, and a broad component formed by charge exchange onto $H^+$ ions accelerated in the postshock flow (Heng 2010). Thus across the face of the nebula we expect a triple Hα profile: a narrow line at the ambient medium radial velocity flanked by the near-side blueshifted flow and the far-side redshifted flow. The narrow component can be seen in Figure 3, both as a faint, unresolved horizontal line in the synthesized long-slit spectra and as a low-velocity peak from the spectrum extracted behind the pulsar.

Although the IFU data does not approach the spatial resolution of Figure 1, the velocity structure provides important constraints on the radial velocity of the ambient medium and on the bow shock inclination. The shock apex, at $x > 0$ in the left panel of Figure 3, is brightest at negative velocities, while the shock body ($x < 0$) has emission extending farther to the red. In the cross-cut spectrum (middle panel) we can see the limb spreads to larger positive velocities. These both suggest a symmetry axis inclined *out* of the plane of the sky, as we quantitatively fit in Section 5. Here, with the red and blue components well separated, we can measure the narrow component, seen here for the first time in a pulsar Hα bow shock.

### 3. Long Baseline Timing Proper Motion

Over the past 16 yr, J1959 has been regularly monitored with the Nançay Radio Telescope (NRT) in France, primarily at the frequency of 1.4 GHz. This monitoring provided a precise

pulse ephemeris for use in the gated Very Long Baseline Array (VLBA) observations (Section 4). These data can also be used to measure the proper motion of the pulsar after careful long-term filtering to model and remove residuals due to changes in the pulsar's orbit and/or dispersion measure (DM). The radio pulsar timing analysis was carried out by analyzing data recorded with the latest pulsar backend in operation at Nançay, NUPPI (a version of the Green Bank Ultimate Pulsar Processing Instrument designed for the NRT; see Guillemot et al. 2016, for details on this backend and on the data analysis procedure). The radio timing data consisted of times of arrival (TOAs) recorded between MJD 55903 (2011 December 8) and 59361 (2021 May 27). For each NRT observation we extracted one TOA per 128 MHz of bandwidth to be sensitive to potential DM changes, and per 10 minutes of observation to get sufficient signal-to-noise while having short integration times compared to the orbital period. In order to exclude data recorded during eclipses or ingress/egress phases we conservatively discarded TOAs with orbital phases ranging from 0.17 to 0.38, with 0.25 corresponding to superior conjunction.

The TOA data were then analyzed using the TEMPO2 (Hobbs et al. 2006) and ENTERPRISE (Ellis et al. 2019) pulsar timing libraries. After a first step of the analysis where we modeled the rapid orbital period variations with nine orbital frequency derivatives with TEMPO2, we then fixed the latter parameters as well as the longitude of periastron passage and the orbital eccentricity, and used ENTERPRISE to fit for the astrometric, rotation, and other binary parameters. We also fitted for white noise in the TOAs, achromatic red noise representing potential pulsar spin noise or noise from imperfect orbital
modeling, and chromatic red noise representing potential DM noise. This analysis yielded $\mu_{R.A.} = -15.81 \pm 0.05\text{mas y}^{-1}$ and $\mu_{\text{decl.}} = -25.54 \pm 0.08\text{mas y}^{-1}$, i.e., about $10 \times$ the accuracy of the proper motion reported in Arzoumanian et al. (1994).

## 4. Very Long Baseline Interferometry Astrometry

J1959's dispersion is $29.1 \text{ cm}^{-3} \text{ pc}$, giving a model distance of 1.73 kpc based on the YMW16 model (Yao et al. 2017), or 2.5 kpc based on the NE2001 model (Cordes & Lazio 2002). Given that uncertainties in the global electron density distribution models limit the accuracy of DM-based distance estimates, and that much of the interpretation is sensitive to the distance, we embarked on a new campaign to better constrain the distance of J1959 via a measurement of annual geometric parallax.

J1959 was observed twice by the VLBA under project code BD179 in 2015 during the target selection phase of the MSPSRπ project (Vigeland et al. 2018; Ding et al. 2020, 2021). These observations resulted in one detection and one nondetection, and J1959 was not selected for further MSPSRπ observations. It was later recognized that the nondetection observation took place near the time of pulsar superior conjunction, where the radio signal is absorbed by the ionized companion wind.

We subsequently observed J1959 a total of 13 further times under a dedicated VLBA program (project code UZ002) over the period 2017 October to 2019 April. These additional observations were timed to avoid superior conjunction, when the pulsar emission was predicted to be eclipsed. All observations were made using the same recording setup, which sampled 256 MHz of bandwidth in dual polarization, selecting regions of lower radio frequency interference in the range 1392–1744 MHz, and correlated with the DiFX software correlator in Socorro (Deller et al. 2011). Around 75 minutes of on-target time was acquired in each epoch, with the exception of the initial BD179 epoch, which was around half the duration.

The VLBA data were reduced using a ParselTongue (Kettenis et al. 2006) pipeline that has previously been used for numerous VLBA pulsar data sets (e.g., Deller et al. 2019) and can be obtained from https://github.com/dingswin/psrvlbireduce/. The source J200352+232800, separated by ~3° from J1959, was used as the primary phase and delay reference calibrator, while J1925+2106 was used to calibrate the instrumental bandpass. For the measurement of parallax and proper motion at ~1.5 GHz as attempted here, the precision of the relative astrometry of the target source is typically dominated by the residual direction-dependent ionospheric propagation effects between the phase calibrator source and the target. Accordingly, and as been used in most recent pulsar astrometry campaigns, we made use of a secondary “in-beam” phase calibrator source located closer to J1959 on the sky, which was identified in preliminary MSPSRπ observations. This source, J200022.8+2059338, is located just 16' from J1959, enabling an order-of-magnitude reduction in the systematic astrometric offsets compared to the primary phase calibrator. The positions of all calibrator sources are shown in Table 1; the positions obtained for J1959 are directly tied to our estimated position for J200022.8+205938, such that any constant offset in the position of J200022.8+205938 would lead to an equivalent (time-independent) translation in the reference position for J1959. Following standard practise, we concatenated the data from all epochs in order to build a high-

<table><caption>**Table 1**
Reference Positions for VLBA Calibrator Sources</caption><thead><tr><td>Calibrator</td><td>R.A. (J2000)</td><td>decl. (J2000)</td></tr></thead><tbody><tr><td>J2003+2328</td><td>20<sup>b</sup>03 <sup>m</sup>52.56148</td><td>23.28<sup>s</sup>00.3281</td></tr><tr><td>J200022.8+205938</td><td>20<sup>b</sup>00<sup>m</sup>22.82142</td><td>20.59<sup>s</sup>38.2975</td></tr></tbody></table>

fidelity model of each calibrator source, and rereduced the data using these final models and after excising radio frequency interference. This yielded a data set with 14 images of J1959, with S/N values ranging from 6 to 20, from which positions and statistical uncertainties were extracted with an image-plane fit. The resultant position time series for the pulsar can be obtained from https://github.com/dingswin/publication_related_materials/tree/main/J1959+2048.

We followed the approach taken by Deller et al. (2019) and Ding et al. (2021) to fit the parallax and proper motion of J1959 from this set of positions, namely, (1) adding an estimated systematic astrometric uncertainty to each epoch based on the empirical estimator described by Deller et al. (2019), and then (2) using a bootstrap process (Efron & Tibshirani 1991) to infer the parallax and proper motion (and their uncertainties). After applying the estimated systematic uncertainty, the reduced $\chi^2$ of a simple least squares fit to the full data set of 14 epochs was 2.1, indicating that the per-epoch astrometric uncertainty may have been somewhat underestimated; however, the use of a bootstrap sampling process provides robustness against this. We performed two separate fits, one with the proper motion fixed to the timing values, and one where it is also fit from the VLBA data—no significant effect on the fitted parallax is seen, and we use the timing-constrained values from this point on. With this, the final parallax is $\pi = 0.39_{-0.10}^{+0.09}$ mas; with the substantial uncertainties, the resulting distance distribution is distinctly non-Gaussian. The bootstrap-inferred probability distribution function (PDF) of d is shown in Section 6.

Armed with this new distance estimate, we can correct the observed pulsar proper motion to the rest frame at the pulsar position, to facilitate a comparison with the observed bow shock geometry. Using the Oort constants from Bovy (2017) we compute a differential rotation correction $\Delta\mu_{R.A.} = -1.18$ mas yr$^{-1}$, $\Delta\mu_{\text{decl.}} = -2.99$ mas yr$^{-1}$. Alternatively, assuming a flat Galactic rotation curve with $v = 220$ km s$^{-1}$ at the Solar circle, we obtain $\Delta\mu_{R.A.} = -2.57$ mas yr$^{-1}$, $\Delta\mu_{\text{decl.}} = -4.23$ mas yr$^{-1}$, i.e., significantly different. These reduce the total proper motion from 30.03 mas yr$^{-1}$ to 26.9 mas yr$^{-1}$ or 25.1 mas yr$^{-1}$, respectively.

The proper motion and distance measurements allow us to correct the observed spindown rate for the Shklovskii (1970) effect

$$ \dot{P} = \dot{P}_{\text{obs}} - 2.43 \times 10^{-21} P \mu_{\text{mas y}^{-1}}^2 d_{\text{kpc}}, \quad (1) $$

which, using the best-fit parallax d, gives J1959's intrinsic $\dot{P}$ as $7.7 \times 10^{-21}$ s s$^{-1}$, $>2 \times$ smaller than the observed value.

## 5. Modeling

Wilkin 2000 (hereafter W00) has given an analytic model for the shape of the contact discontinuity (CD) of a thin (momentum-balance) bow shock

$$ r(\theta) = r_0[3(1 - \theta \cot \theta)(1 - c_2/12) + 3(c_2 \sin^2 \theta)/4]^{1/2} / \sin \theta, \quad (2) $$
where $\theta$ is measured from the bow shock nose and the pulsar momentum flux is axisymmetric with a distribution $\propto 1 - c_2/3 + c_2 \cos^2 \theta_*$ about a wind axis. In general the wind axis will be misaligned with the pulsar velocity (e.g., as for PSR J2124-3358; see Romani et al. 2017), but for J1959 we see no such asymmetry and so restrict our analysis to the aligned case $\theta = \theta_*$, measured from the direction of motion. A spherical wind has anisotropy parameter $c_2 = 0$, an equatorial wind $c_2 < 0$. The characteristic standoff distance is related to the pulsar spindown power $\dot{E}$ and velocity $v_p$ and the ISM hydrogen number density $n_H$ and mean mass per particle $\gamma_H m_p$ (giving a mass density $\rho = \gamma_H m_p n_H$) by

$$r_0 = (\dot{E}/4\pi c \rho v_p^2)^{1/2}. \qquad (3)$$

For a “thin” shock, the CD with apex at $r_0$ marks the H$\alpha$ front. In practice, the H$\alpha$ emission comes from near the forward shock and the observed H$\alpha$ standoff should be $r_a \sim 1.3r_0$ at the apex (Aldcroft et al. 1992; Bucciantini 2002). Also the shock should spread somewhat faster than the CD, as postshock pressure widens the standoff downstream. A simple approximation to match numerical models increases the transverse scale by $1.25 \times$ (which we refer to here as a “wide” shock model; BR14). Another source of widening can be due to anisotropy in the pulsar wind, if we take $c_2 < 0$; $c_2 \approx -3/4$ is close to our “wide” model. A combination of the two effects can occur with anisotropy broadening the shock toward the wind equator and postshock pressure broadening the shock past its apex.

To match a W00-style model to the data, we need at least three parameters: the CD standoff scale $r_0$, the symmetry axis' orientation $PA_{BS}$, and its inclination to the line of sight $i_{BS}$. Since the bow shock models have rounded fronts, the projected image is only weakly dependent on $i_{BS}$ for values near $\pi/2$. The H$\alpha$ velocity structure is substantially more sensitive. Ideally one would use the mass density and velocity expressions of W00 to fit the position-velocity cube for the bow shock parameters. However, these assume a thin pure momentum-balance flow all at the CD. In practice we find that while such models reproduce the general structure of observed bow shock H$\alpha$ position-velocity cubes, the match is insufficient for a meaningful fit. This is likely due to the varying width of the postshock flow, and density/ionization variations in the impinging ISM. Instead we focus on the zone behind the highly curved apex, where the shock limb presents a simple wedge. If we assume azimuthal symmetry we can use the wedge opening angle $\theta_{op}$ and the velocity splitting in this zone to get a nearly model-independent measurement of the inclination, for a shock axis near the plane of the sky.

Behind the pulsar, at the $-3r_a$ location of our face spectrum, we measure $\theta_{op} = \pm 18^\circ$ (between the analytic values $16^\circ$ for the thin shock and $20^\circ$ for the wide shock models). The flow along this shock provides the red and blue line shifts measured in Figure 3 (values are on the plot; the small $\pm 0.2 \text{ km s}^{-1}$ line center statistical fit errors are likely dominated by systematic effects). For $\theta_{op} = \pm 18^\circ$ this requires a velocity axis inclination $i_v = 2^\circ \pm 18^\circ$ out of the plane of the sky and effective velocity parallel to the H$\alpha$ shock front $v_\parallel = 150 \text{ km s}^{-1}$ (Figure 4). While $\theta_{op}$ is difficult to measure in the presence of the shock crenelation and other detailed structure, $i_v$ is quite insensitive to

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Measured blue line velocity shift $\Delta\lambda_B$ and the red/blue line shift velocity ratio are shown by the shaded horizontal bands, as measured from the H$\alpha$ spectrum fit in Figure 3. Diagonal lines give the inclination constraints for various bow shock opening angles (estimated at $\theta_{op} = \pm 18^\circ$ from the H$\alpha$ image). This implies $i_v = 2^\circ \pm 18^\circ$ out of the plane of the sky and a shock front $v_\parallel \approx 150 \text{ km s}^{-1}$ (solid lines). Estimates for $\theta_{op} = \pm 17^\circ$ and $\pm 19^\circ$ are shown by the dotted and dashed lines, respectively.</figcaption></figure>

<table><caption>**Table 2**
PSR J1959+2048 Bow Shock Geometry</caption><thead><tr><th>Model<sup>a</sup></th><th>c<sub>2</sub></th><th>θ<sub>0</sub><br>(")</th><th>PA<br>(Deg)</th><th>χ<sup>2</sup>/DoF<br>32093 DoF</th></tr></thead><tbody><tr><td>Thin</td><td>0</td><td>3.624<sup>+0.019</sup><sub>-0.013</sub></td><td>214.1 ± 0.2</td><td>1.239</td></tr><tr><td>Thin</td><td>-3/4</td><td>3.903<sup>+0.013</sup><sub>-0.011</sub></td><td>213.3 ± 0.2</td><td>1.167</td></tr><tr><td>Wide</td><td>0</td><td>2.836<sup>+0.010</sup><sub>-0.006</sub></td><td>213.2 ± 0.2</td><td>1.165</td></tr><tr><td>Wide</td><td>-3/4</td><td>2.986<sup>+0.023</sup><sub>-0.012</sub></td><td>212.9 ± 0.2</td><td>1.241</td></tr></tbody></table>

**Note.**

this value. However $v_\parallel$ varies by $\sim 8 \text{ km s}^{-1}$ for a $1^\circ$ opening angle variation.

Armed with $i_v$, we can fit the bow shock image for the shock scale and position angle (PA). We fit to a $23'' \times 13''$ region around the apex, masking the regions around bright stars. We used the affine invariant Markov Chain Monte Carlo algorithm of Goodman & Weare (2010), implemented through the python package EMCEE, to sample the likelihood function and obtain posterior distributions for the fit parameters. Our metric is the rms residual and we record the geometrical parameters and the goodness of fit in Table 2 (the nebula flux and background normalization are also free). The nebula PA differs slightly between models, but is quite consistent between the best-fit narrow/$c_2 = -3/4$ and wide/$c_2 = 0$ models, which we adopt as fiducial. We also varied the mask regions and fit to the continuum-subtracted H$\alpha$ frames, getting consistent fit parameters for all models.

We can now discuss the 3D orientation of the pulsar motion and bow shock. The $\gamma$-ray eclipse requires occultation of the

<footnote><sup>a</sup> Velocity angle to plane of sky fixed at $i_v = 2.18^\circ$.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left: continuum-subtracted Hα image of J1959+2048. The timing proper motion is shown by the magenta arrow extending from the pulsar counterpart position, with a 10σ error ellipse at the tip. The proper motion corrected for differential Galactic rotation (red) has reduced velocity, but is very close to the nebula symmetry axis (green arrow); proper motion arrows show the 1000 yr shift. The agreement limits the transverse motion of the local ISM (see the text). Right: side view of our inferred geometry for the J1959 system, with the contact discontinuity (CD) and Hα-emitting forward shock, here for the wide/c2 = −3/4 model. The pulsar motion heads out of the plane of the sky, resulting in a larger bow shock red component shift and smaller blue component shift. The binary orbit, within 7° of edge-on as viewed from Earth, likely has an orbit normal aligned with the bow shock axis (red line), but this normal can be offset, while still in the plane of the sky (red dotted ellipse).</figcaption></figure>

pulsar by the companion photosphere; with the photometrically estimated companion size, this requires $i = 90 \pm 7^{\circ}$ for the orbit. Note that this places the orbit normal within $7^{\circ}$ of the plane of the sky, consistent with our $i_v = 2^{\circ}18$. Of course, these measurements do not themselves require alignment; both could lie in the plane of the sky at a large angle to each other. However, if part of the bow shock spreading is due to pulsar wind $c_2 < 0$ anisotropy, this should be at the pulsar spin equator. We see no evidence for this spin being misaligned with $v_{PSR}$. Finally if, as expected, the spin axis of this recycled pulsar coincides with the orbit normal (common angular momentum source), we can then infer that the orbit normal and space velocity are aligned. This might be directly checked with pulsar scintillation studies, which can be sensitive to the orbital velocity of the pulsar, checking its projection on the plane of the sky (Ord et al. 2002).

Our narrow component velocity implies that the local ISM has a radial velocity $\sim -7.8$ km s$^{-1}$. The pulsar motion out of the plane of the sky adds another $v_{PSR} \sin(i_v) \approx 15$ km s$^{-1}$ toward Earth (see Figure 5). The good agreement between the nebula axis and the proper motion vector suggests a transverse motion $\sim 3$ km s$^{-1}$ or less. However, the correction for differential rotation at the $\sim 2.6$ kpc pulsar distance dominates this estimate. If we use a flat $V = 220$ km s$^{-1}$ rotation curve to correct instead of the locally measured Oort constants, we get an angle offset of $\sim 1^{\circ}3$, implying a transverse motion of $\sim 7$ km s$^{-1}$. In any case transverse motion appears to be small. Note that the differential rotation contributes $3-5$ mas y$^{-1}$ of the pulsar proper motion; the local space velocity is thus decreased.

## 6. Discussion and Conclusions

In BR14 it was shown how with known kinematics the Hα flux of a bow shock can be used to derive the impinging momentum flux of the ISM (in the pulsar rest frame) and compare with the momentum balance from the pulsar wind at the CD (Equation (3)). With our planetary nebula flux calibrations we establish the sensitivity of the SOAR Hα exposures at $1.57 \pm 0.08 \times 10^{-16}$ erg$^{-1}$cm$^{-2}$count$^{-1}$. Measuring the flux of the apex region from $+\theta_a$ in front of the pulsar to $-2\theta_a$ behind, we find an apex flux of $1.7 \times 10^{-3}$Hα cm$^{-2}$s$^{-1}$, in good agreement with the independent measurement in BR14. For a wide shock (or for $c_2 = -3/4$; Equation (2)), integrating back to $-2\theta_a$, the nebula sweeps up

$$ \dot{N} = \pi(3.55 r_a)^2 v n_H f_{HI} \quad (4) $$

neutral H atoms, where $r_a = 1.3(1 + 2c_2/3)^{1/2}\theta_0 d$ and $f_{HI}$ is the neutral fraction. These neutrals pass into the shocked ISM, where they suffer collisional excitation and charge exchange, radiating Hα before ionization. We next need an efficiency for incoming protons to emit Hα. Unfortunately, most studies of nonradiative shocks are aimed at supernova remnant velocities $v > 10^3$ km s$^{-1}$, and the low-velocity case appropriate to J1959 and many other pulsar bow shocks is poorly explored.

Indeed, with a broad-to-narrow line ratio $I_b/I_n$ fit as $\approx 4$ (averaged over the red and blue components in the postapex spectrum in Figure 3) the narrow component is much weaker than the $I_b/I_n \sim 0.7-1.1$ seen for nonradiative shocks in supernova remnants. This is almost certainly a result of the very low effective shock velocity, but also implies that the
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Hα production efficiency (left scale) and broad-to-narrow line ratio (right log scale), computed from modeling in Heng & McCray (2007). Multiple transfers and excitations increase the broad line component at low shock velocity.</figcaption></figure>

postshock electrons and ions are not fully equilibrated in the emission zone. The narrow component is sufficiently detected to constrain the upstream radial velocity, but is not resolved here. With neutral medium temperatures of <10⁴ K, we expect thermal velocity widths of <10 km s⁻¹, so full exploitation of the narrow component requires much higher S/N and spectral resolution. This would be quite valuable, as the narrow line profile can reveal not only the state of the upstream medium, but may also show (e.g., though superthermal core, non-Gaussian wings) evidence for preheating by bow shock-generated cosmic rays escaping upstream.

To estimate the Hα efficiency we can use Figure 6 of Heng & McCray (2007), which gives, for each incoming H, the number of initial excitations N₀ (which produce the narrow line) and the number of charge transfers N_T and subsequent excitations N_En (which produce the broad lines) as a function of the shock velocity. With their assumption that $\epsilon_{tr} \approx 0.03$ of the transfers produce an Hα photon, while the narrow line emits in Case B ($\epsilon_B \approx 0.27$ Hα/H) and the broad line emits in Case A ($\epsilon_A \approx 0.048$ Hα/H), we can write the total Hα efficiencies as the sum of the narrow line $\epsilon_n$ and broad line $\epsilon_b$ contributions:

$$ \epsilon_{H\alpha} = \epsilon_n + \epsilon_b = \epsilon_B N_{E0} + (\epsilon_A N_{En} + \epsilon_{tr} N_T). \quad (5) $$

These efficiencies are shown in Figure 6 for the case where the postshock electrons and ions do not reach full equilibrium, along with the approximate $\epsilon_{H\alpha} \approx 0.6v_7^{-1/2}$ used in BR14, which was matched at 10³ km s⁻¹. We can further compute $I_b/I_n = \epsilon_b/\epsilon_b$, also plotted in Figure 6. We see that the $I_b/I_n \approx 4$ that we find for J1959 can be produced in such nonequilibrated shocks for low shock velocity. For shocks with full electron-ion equilibration, ionization is enhanced and $I_b/I_n$ is substantially smaller and inconsistent with the J1959 data. Note that, unlike the v in Equation (4), the velocity used in $\epsilon_{H\alpha}$ is the effective shock velocity, a fraction of the pulsar velocity with respect to the local ISM for an oblique shock.

This incoming momentum flux balances against the pulsar momentum flux $\dot{E}/c = I\Omega\dot{\Omega}/c = 4\pi^2 P/c P^3$ at the CD. Correcting as above for the Shklovskii (1970) effect, we use a spindown luminosity $\dot{E} = 7.32 \times 10^{34} I_{45} \text{erg s}^{-1}$.

Now, noting that the observed Hα flux $f_{H\alpha}$ gives the sweep up rate as

$$ \dot{N} = 4\pi d^2 f_{H\alpha} / \epsilon_{H\alpha} \quad (6) $$

we can combine Equations (3), (4), (5), and (6) to solve for the moment of inertia

$$ I = 16\pi c \gamma_H m_p v d^2 f_{H\alpha} / [34.5 f_{HI} \epsilon_{H\alpha} \Omega \dot{\Omega}], \quad (7) $$

where the numerical factor in the denominator is for the $c_2 = -3/4$ shock case, $\dot{\Omega}$ is corrected for the Shklovskii effect, and the velocities are inferred from the (differential rotation and projection corrected) proper motions. With $v \propto d$ our moment of inertia estimate is very sensitive to the pulsar distance, with $I \propto d^{7/2}$. Note that as in BR14, the upstream density $n_H$ (and hence the angular size of the CD apex angle $\theta_0 = r_0/d$) drops out of the I estimate. In principle the I estimate is directly proportional to the observed flux. In practice, it also depends on the efficiency assumption, projection effects (retired here via our iBS measurement), and the upstream neutral fraction $f_{HI}$.

Even beyond its d sensitivity the effective shock velocity strongly affects the emission efficiencies for small v. Interestingly the large observed $I_b/I_n = 4$ behind the shock indicates a very low effective velocity in this region. Of course, the shock is highly oblique in this region, with opening angle $\theta_{op} = 18^\circ$, so we expect an effective local shock speed $v_p \sin(18^\circ) = 0.3v_p$. The predicted $I_b/I_n$ matches quite well with the observed value at a speed $v_{eff} \approx 150$ km s⁻¹. At this speed the very large increase in $I_b$ is due to the low ionization rate, which allows many excitations.

In Figure 7 we plot the log of the neutron star moment of inertia (in units of $10^{45}$ g cm⁻²) inferred from Equation (7) for various effective shock velocities as a function of distance, using the Heng & McCray (2007) model efficiency. The leftmost curve shows the value for the normal shock at the apex where the full pulsar velocity $v_p$ is appropriate, while the curves to the right show I for 0.5$v_p$ (e.g., for a shock obliquity 60°) and for 0.3$v_p$ (for the effective velocity at $\theta_{op} = 18^\circ$ at $-3r_a$). Horizontal lines show the moment of inertia values from Lattimer & Schutz (2005; at J1959's 1.8 M☉ mass) for stiff equations of state allowing 2 M☉ neutron stars. This plot assumes a neutral upstream medium $f_{HI} = 1$.

We see that our bow shock-derived $I_{45}$ values are quite large at our most probable parallax distance estimate. This is due to the large momentum flux required to form J1959's bright extended bow shock at large distances. These large values would require a very stiff equation of state. This is mitigated for $\le 2$ kpc distances and effective velocities below the normal shock apex value. Improved parallax measurements and improved modeling of low-velocity nonradiative shocks can help probe this tension. Higher $\epsilon_{H\alpha}$ than assumed here could also help. Finally our model assumes momentum flux $\dot{E}/c = I\Omega\dot{\Omega}/c$ for a relativistic pulsar wind. While appropriate for almost all pulsars, we note that J1959 presents the only bow shock surrounding a “spider” (companion-evaporating) pulsar. Since the companion wind shocks against the pulsar wind (Kandel et al. 2021), and the postshock pressure does work against the massive companion wind, the total wind
<figure><img src="image_10.png" /></figure>

**Figure 7.** Bow shock Hα flux estimates of the neutron star moment of inertia. The red curves show Log($I_{45}$) for three effective shock velocities, using the HM07 Hα efficiencies, as a function of pulsar distance. All assume a neutral upstream medium and a wide spherical wind bow shock. For comparison Log [$I_{45}$($1.8 M_⊙$)] from three stiff equations of state allowing 2 M⊙ neutron stars are shown (dotted lines). The dashed line shows the distance PDF (arbitrary normalization) from the parallax data.

momentum flux should be increased beyond that of the direct relativistic pulsar wind. Further modeling is needed to see if such an increase is significant, helping to explain J1959's unusually bright bow shock.

Our study has extracted unusually complete geometrical and kinematic information on this black widow binary pulsar. We now have good knowledge of the 3D motion of the pulsar and can relate this to the binary (and implicitly spin) axis. Their apparent near-alignment, which might be checked with an Ord et al.-style (2002) scintillation velocity study, argues for a physical connection or strong selection effects.

We conclude that the quest to use the bow shock to measure a precise equation of state is frustrated by the very strong distance dependence coupled with $ε_{Hα}$'s extreme velocity sensitivity at the low effective speeds relevant for J1959's bow shock. Thus, improved modeling of low-velocity nonradiative shocks, including the highly oblique shocks appropriate to much of the face of the nebula, is needed before strong constraints can be made on the moment of inertia. Of course, since the dense matter equation of state is universal (while presently unknown), it may make more sense to turn this argument around: with the known pulsar mass providing (via the equation of state) a known *I*, one might use the Hα flux modeling to extract a precise distance.

We wish to thank Andrei Tokovinin and Ceasar Briceno at SOAR for assistance with the SAMI observations. At the AAT,

Andy Green assisted with the KOALA data acquisition; we also thank Newton Cheng for assistance with their reduction.

R.W.R. was supported in part by NASA grant 80NSSC21K0896.

The Nançay radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS), and partially supported by the Region Centre in France. We acknowledge financial support from the “Programme National Hautes Energies” (PNHE) of CNRS/INSU, France.

The authors acknowledge use of the Very Long Baseline Array under the US Naval Observatory's time allocation. The Very Long Baseline Array (VLBA) is operated by the National Radio Astronomy Observatory (NRAO). The NRAO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

## ORCID iDs

Roger W. Romani https://orcid.org/0000-0001-6711-3286

Adam Deller https://orcid.org/0000-0001-9434-3837

Lucas Guillemot https://orcid.org/0000-0002-9049-8716

Robert T. Zavala https://orcid.org/0000-0002-9402-2870

Ismaël Cognard https://orcid.org/0000-0002-1775-9692

## References

Aldcroft, T. L., Romani, R. W., & Cordes, J. M. 1992, ApJ, 400, 638

Arzoumanian, Z., Fruchter, A. S., & Taylor, J. H. 1994, ApJL, 426, L85

Bovy, J. 2017, MNRAS, 468, L63

Brownsberger, S., & Romani, R. W. 2014, ApJ, 784, 154

Bucciantini, N. 2002, A&A, 387, 1066

Clark, K., Kerr, M., & Breton, R. 2022, NatAs, submitted

Cordes, J. M., & Lazio, T. J. W. 2002, arXiv:astro-ph/0207156

Deller, A. T., Brisken, W. F., Phillips, C. J., et al. 2011, PASP, 123, 275

Deller, A. T., Goss, W. M., Brisken, W. F., et al. 2019, ApJ, 875, 100

Ding, H., Deller, A. T., Fonseca, E., et al. 2021, ApJL, 921, L19

Ding, H., Deller, A. T., Freire, P., et al. 2020, ApJ, 896, 85

Efron, B., & Tibshirani, R. 1991, Sci, 253, 390

Ellis, J. A., Vallisneri, M., Taylor, S. R., & Baker, P. T. 2019, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust Pulsar Inference SuitE, Astrophysics Source Code Library, ascl:1912.015

Ellis, S. C., Ireland, M., Lawrence, J. S., et al. 2012, Proc. SPIE, 8446, 84460V

Frew, D. J., Bojičić, I. S., & Parker, Q. A. 2013, MNRAS, 431, 2

Fruchter, A. S., Stinebring, D. R., & Taylor, J. H. 1988, Natur, 333, 237

Goodman, J., & Weare, J. 2010, CAMCOS, 5, 65

Guillemot, L., Smith, D. A., Laffon, H., et al. 2016, A&A, 587, A109

Heng, K. 2010, PASA, 27, 23

Heng, K., & McCray, R. 2007, ApJ, 654, 923

Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655

Kandel, D., Romani, R. W., & An, H. 2021, ApJL, 917, L13

Kettenis, M., van Langevelde, H. J., Reynolds, C., & Cotton, B. 2006, in ASP Conf. Ser. 351, Astronomical Data Analysis Software and Systems XV, ed. C. Gabriel et al. (San Francisco, CA: ASP), 497

Kulkarni, S. R., & Hester, J. J. 1988, Natur, 335, 801

Lattimer, J. M., & Schutz, B. F. 2005, ApJ, 629, 979

Ord, S. M., Bailes, M., & van Straten, W. 2002, ApJL, 574, L75

Romani, R. W., Slane, P., & Green, A. W. 2017, ApJ, 851, 61

Shklovskii, I. S. 1970, SvA, 13, 562

Tokovinin, A., Cantarutti, R., Tighe, R., et al. 2016, PASP, 128, 125003

van Kerkwijk, M. H., Breton, R. P., & Kulkarni, S. R. 2011, ApJ, 728, 95

Vigeland, S. J., Deller, A. T., Kaplan, D. L., et al. 2018, ApJ, 855, 122

Wilkin, F. P. 2000, ApJ, 532, 400

Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29