Nonradiative shocks accelerate particles and heat astrophysical plasmas. While supernova remnants are the most wellstudied example neutron star NS bow shocks are also nonradiative and Balmer dominated. NS bow shocks are likely ubiquitous in the interstellar medium due to their large speeds imparted at birth and they are thought to be a discrete source population contributing to the Galactic cosmicray spectrum. To date nine NS bow shocks have been directly observed in H images. Most of these shocks have been characterized using narrowband H imaging and slit spectroscopy which do not resolve the multicomponent velocity structure of the shocks and their spatial geometry. Here we present integral field spectroscopy of three NS bow shocks J07422822 J17412054 and J22256535 the Guitar Nebula. We measure the shock properties simultaneously in four dimensions the 2D projected shock morphology the radial velocity structure and the H flux. The broadtonarrow line ratio I SUBbSUBI SUBnSUB is inferred from radial velocity profiles and for J17412054 the narrow line is detected in multiple regions of the shock. The inferred line ratios and widths suggest that NS bow shocks represent a lowshock velocity regime V  200 km sSUP1SUP in which I SUBbSUBI SUBnSUB is high distinct from the shock regime probed by supernova remnants. Our results illustrate a need for nonradiative shock models at velocities lower than previously considered which will reveal the electronion temperature ratios and particle acceleration efficiencies of these bow shocks.
['Neutron stars', 'Stellar bow shocks', 'Shocks', 'Supernova remnants', 'Cosmic rays', '1108', '1586', '2086', '1667', '329', 'Astrophysics - High Energy Astrophysical Phenomena']
Probing the Lowvelocity Regime of Nonradiative Shocks with Neutron Star Bow Shocks
{'Probing the Low-Velocity Regime of Non-Radiative Shocks with Neutron Star Bow Shocks': '<!-- image -->  \n1  \nCahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA 2 Carnegie Science Observatories, Pasadena, CA 91101, USA', 'ABSTRACT': 'Non-radiative shocks accelerate particles and heat astrophysical plasmas. While supernova remnants are the most well-studied example, neutron star (NS) bow shocks are also non-radiative and Balmer-dominated. NS bow shocks are likely ubiquitous in the interstellar medium due to their large speeds imparted at birth, and they are thought to be a discrete source population contributing to the Galactic cosmic ray spectrum. To date, nine NS bow shocks have been directly observed in H α images. Most of these shocks have been characterized using narrowband H α imaging and slit spectroscopy, which do not resolve the multi-component velocity structure of the shocks and their spatial geometry. Here we present integral field spectroscopy of three NS bow shocks: J0742 -2822, J1741 -2054, and J2225 + 6535 (the Guitar Nebula). We measure the shock properties simultaneously in four dimensions: the 2D projected shock morphology, the radial velocity structure, and the H α flux. The broad-to-narrow line ratio ( I b /I n ) is inferred from radial velocity profiles, and for J1741 -2054 the narrow line is detected in multiple regions of the shock. The inferred line ratios and widths suggest that NS bow shocks represent a low shock velocity regime ( V ≲ 200 km/s) in which I b /I n is high, distinct from the shock regime probed by supernova remnants. Our results illustrate a need for non-radiative shock models at velocities lower than previously considered, which will reveal the electron-ion temperature ratios and particle acceleration efficiencies of these bow shocks.  \nKeywords: Neutron stars - Bow shocks - Cosmic rays - Supernova remnants - Interstellar medium', '1. INTRODUCTION': "Optical imaging is the main method by which neutron star bow shocks are detected and characterized, as these shocks are non-radiative and chiefly emit in H α (Kulkarni & Hester 1988). The H α emission arises from collisional excitation and charge exchange of the interstellar medium (ISM) where it confronts the pulsar wind. Only nine H α -emitting neutron star bow shocks are known, likely because most of the neutron stars targeted in bow shock searches reside in regions of low neutral gas fraction (Chatterjee & Cordes 2002; Brownsberger & Romani 2014), and because optical surveys are severely limited by extinction for neutron stars residing deep in the Galactic plane. Nevertheless, many neutron stars are expected to produce bow shocks because the majority have transverse velocities larger than the fast magnetosonic speed of the ISM (Verbunt et al. 2017).  \nNeutron star bow shocks are one of the few observed analogs of non-radiative, Balmer-dominated shocks in super-  \nCorresponding author: Stella Koch Ocker  \nsocker@caltech.edu  \nnova remnants (SNRs; Heng 2010). Similar to SNRs, these shocks are expected to have a narrow H α line component (width ∼ 10 s km s -1 ) at the radial velocity of the ambient ISM arising from collisional excitation of neutral hydrogen crossing the shock interface, and a broad H α line component (width ∼ 100 s km s -1 ) arising from charge exchange between pre- and post-shock particles (Heng 2010). For bow shocks, these two emission components give rise to a triplepeaked radial velocity profile consisting of the narrow H α line flanked by the blue and red-shifted components of the broad line from the near and far sides of the shock (Romani et al. 2022). The broad-to-narrow line intensity ratio ( I b /I n ) and the broad line width ( W b ) are sensitive diagnostics of the shock physics, including electron-ion temperature equilibration (Raymond et al. 2023) and particle energy distributions (Nikoli'c et al. 2013). Moreover, I b /I n and W b constrain the physical velocity of the shock, which combined with proper motion measurements yield an independent distance estimate (van Adelsberg et al. 2008). While I b /I n and W b are routinely constrained for SNRs (e.g. Raymond et al. 2023 and refs. therein), there is only one published measurement of I b /I n for a neutron star bow shock to date (Romani et al. 2022).  \nRadial velocity profiles are also critical to resolving the physical geometry of the shock, as 2D images only constrain the shock in projection. The shock morphology is determined by conditions for ram pressure balance between the pulsar wind and the ISM (e.g. Wilkin 1996). Velocity profiles are sensitive to the shock's inclination relative to the plane of the sky, which combined with a model for the shock structure yields the physical shock stand-off radius and the radial velocity of the neutron star, which cannot be inferred from pulsar timing or radio imaging (Romani et al. 2017; de Vries & Romani 2020; Romani et al. 2022). The ram pressure conditions required to explain observed neutron star bow shocks have been exploited to characterize the energetics of pulsar winds (Kulkarni & Hester 1988; Cordes et al. 1993; van Kerkwijk & Kulkarni 2001) and their anisotropy (Romani et al. 2010), pulsar distances and velocities (Jones et al. 2002; de Vries & Romani 2020), the neutron star moment of inertia (Romani et al. 2017), and ISM density fluctuations (Chatterjee & Cordes 2004; Ocker et al. 2021). Most of these prior studies were based on narrowband H α imaging and slit spectroscopy, necessitating model-based assumptions of the shock inclination angle and subsequent constraints.  \nIn this study, we leverage integral field spectroscopy (IFS) to spatially resolve the H α velocity structure of three neutron star bow shocks: J0742 -2822, J1741 -2054, and J2225 + 6535 (the Guitar Nebula). The observations are described in Section 2. The radial velocity profiles are used to constrain I b /I n , W b , and the shock inclinations (Section 3). We interpret our results in the context of analogous shocks in SNRs in Section 4.", '2. OBSERVATIONS': "Weobserved J0742 -2822, J1741 -2054, and J2225 + 6535 with the Keck Cosmic Web Imager (KCWI) on the Keck II Telescope on UT 2024-02-01 (J0742 -2822) and UT 202407-12 (J1741 -2054 and J2225 + 6535). All sources were observed with the RH1 grating and medium slicer, covering a field-of-view of 16 '' . 5 by 20 '' . 4 at resolutions of 0 . 25 ˚ A and 0 '' . 7 , although seeing was typically 1.5 to 3 times larger (the latter, worst seeing was on 2024-07-12). Nine 600 s exposures were obtained for J0742 -2822 and J1741 -2054 and eighteen 600 s exposures for J2225 + 6535. Offset sky fields were observed throughout each session in addition to flux standards GD50 and G138 -31.  \nData were reduced using the KCWI Data Reduction Pipeline (https://pypi.org/project/kcwidrp/; Neill et al. 2023), described in Morrissey et al. (2018). Bias frames, internal flats, and dome flats were used to correct the science images, and cosmic ray rejection was performed using the Laplacian edge detection algorithm described in van Dokkum (2001). Continuum bars and arc lamps were used to perform the geometry and wavelength calibrations. A sky model was con-  \nstructed from the offset fields and sky subtraction quality was verified using off-target regions. Flux standards GD50 and G138 -31 were used to perform absolute flux calibration at an error level < 4% . Individual flux and wavelength calibrated exposures were combined using drizzle (Gonzaga et al. 2012).", '3. RESULTS': 'Figure 1 shows our KCWI observations of the three bow shocks, illustrating the evolution of the shock emission over ≈ 6559 ˚ A -6568 ˚ A , equivalent to radial velocities ≈ -250 to +150 km/s (the exact range shown differs by shock). The spatial structure of the emission has a strong radial velocity dependence: limb brightening, lobes, and rings are apparent at different velocities and are indicative of changing conditions across the apex, body, and tail of each shock. We interpret these features using spectra averaged over different portions of the shocks. Figure 2 shows the aperture-averaged velocity profiles and the regions of the shocks they correspond to. These regions were chosen to illustrate the highest signal-to-noise (S/N) narrow line detections and evolution of velocity features from shock apex to tail. The narrow line is detected in all three shocks. For J0742 -2822 and J2225 + 6535, the narrow line is barely detected just behind the shock apexes, whereas for J1741 -2054, the narrow line is detected across multiple parts of the shock.', '3.1. Interpretation of Radial Velocity Profiles': "We model the velocity profiles as the sum of N Gaussian components, where N is determined by the number of distinct peaks in the velocity profile and validated through comparison of the reduced χ 2 to that obtained for a model with one less component. The number of components required to explain each profile depends on spatial location, and in turn whether the faint narrow line is detected, which requires accurate sky subtraction. When the narrow line is detected, we estimate I b /I n using the best-fit amplitudes of the narrow and broad components. We quote lower limits when the narrow line is covered by fewer than four resolution elements, i.e., it is not well separated from the broad line emission. The interpretation of I b and I n for bow shocks differs from the interpretation of shocks in SNRs. In the case of bow shocks, which are essentially conical shells, the velocity profile traces the line-of-sight through the shock viewed in projection. The observed narrow line intensity is thus the sum of the I n contributions from both the near and far sides of the shock, both of which are at the same ambient ISM velocity. We therefore estimate I b /I n by taking the sum of the blue and red-shifted broad line amplitudes ( I blue b , I red b ) and dividing by the amplitude of the narrow line:  \nI b I n ≡ I blue b + I red b I n . (1)  \nFigure 1. KCWI observations of three pulsar bow shocks, J0742 -2822 (top), J1741 -2054 (middle), and J2225 + 6535 (the Guitar Nebula; bottom). Images are shown for 12 radial velocity bins centered around the peak shock emission; the spectra have been smoothed from the native resolution of 0.25 ˚ Ato show the large velocity range (larger bins are used for the Guitar due to its wider extension to negative velocities). These slices are continuum subtracted using off-band portions of the data cube; photospheric H α absorption causes many stars to show negative residuals, centered on different low velocity channels. Over-subtracted stars are shown by the grey dashed circles and are masked in subsequent analysis. The white tick marks are 5 '' apart; each image is about 13 '' on a side ( ≈ 0.1 pc, 0.02 pc, and 0.04 pc on a side, from top to bottom). The dynamic range of each image is normalized to its peak flux level in order to illustrate shock structure at the lowest and highest velocities.  \n<!-- image -->  \nInterpretation of the broad-line width W b is less straightforward because the apparent line width is a function of inclination angle. When the shock is close to the plane of the sky ( i ≈ 0 · ), the widths of the blue and redshifted broad lines ( W blue b , W red b ) should be roughly equal. For inclination angles out of the plane of the sky, W blue b > W red b , while for inclination angles into the plane of the sky, W blue b < W red b . The exact dependence of W blue b and W red b on inclination will depend on the 3D spatial morphology of the shock, and hence will differ depending on the shock stand-off radius, the properties of the pulsar wind (e.g. whether it is equatorial or not), and the ISM density. Producing such 3D models for each shock is beyond the scope of this paper, but we nonetheless report empirical constraints on the apparent line width to quantitatively illustrate the velocity range of each shock.  \nEven in the absence of detailed shock models, we can place constraints on the plausible range of inclination angles based on the difference in velocity extrema of the blue and redshifted broad lines ( V blue , V red ). If we assume an axisymmetric shock, then | V blue | ≈ | V red | when the inclination relative to the plane of the sky is i ≈ 0 · . When i < 0 · (out of the plane of the sky), | V blue /V red | > 1 , whereas for i > 0 · (into the plane of the sky), | V blue /V red | < 1 . We use the thinshell approximation for an isotropic stellar wind provided by Wilkin (1996) as a starting point, noting caveats below. In this scenario, the velocity tangent to the shock at the contact discontinuity is given by Wilkin (1996) as  \nv t = v ∗ √ ( θ -sin θ cos θ ) 2 +(˜w 2 -sin 2 θ ) 2 ) 2 α (1 -cos θ ) + ˜w 2 , (2)  \nwhere θ is the position angle along the shock measured relative to the shock apex, v ∗ is the ISM velocity, α = v ∗ /v w is the ratio of the ISM velocity to the stellar wind velocity ( ≪ 1 for relativistic pulsar winds), and ˜ w 2 is given by  \n˜ w 2 = 3(1 -θ cot θ ) . (3)  \nWe evaluate v t for a range of inclination angles and equivalent θ , and estimate the resulting ratio of V blue /V red . Here we implicitly assume that the observed radial velocity depends on θ in the same fashion as v t , normalized such that negative inclination i corresponds to inclination out of the plane of the sky. Note that the dependence of v t on v ∗ does not matter in this context because we are only interested in the velocity ratio between the near and far sides of the shock. Figure 3 shows | V blue /V red | vs. i for different values of α ; for relativistic pulsar winds, we use the limit α → 0 .  \nIn general, the velocity of the H α emission will not exactly follow that of the contact discontinuity; e.g., Bucciantini (2002) find that while v t in the external layer where H α emission is produced (between the contact discontinuity and forward shock) tends to track the velocity of the contact discontinuity near the apex, v t in the external layer increases more rapidly in the tail of the shock. The Wilkin model assumes conservation of momentum (isothermal shocks), whereas non-radiative shocks conserve energy and gas accelerates down the shock, leading to v t larger than in the isothermal case. Anisotropy in the pulsar wind and density variations in the ISM will further modify the velocity structure of the shock from the simple form in Equation 2 (Wilkin  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 2. Heliocentric radial velocity profiles for J0742 -2822 (left), J1741 -2054 (middle), and J2225 + 6535 (right). The upper panels illustrate the regions of each shock used to calculate the aperture-averaged velocity profiles shown in the bottom panels, labeled (a,b,c,d). The images in the upper panels are averaged over the full velocity range of each shock, and the grey dashed circles indicate stars that are subsequently masked. The velocity profiles are each normalized to unity peak flux and shown with an arbitrary offset so that specific velocity features can be compared. Error bars correspond to the standard deviation of flux within each region at a given wavelength. Solid lines indicate the best-fit model for the velocity profile, which is treated as a sum of Gaussian components (§ 3).  \n<!-- image -->  \n2000; Vigelius et al. 2007). The inclination angles given below should thus be regarded as fiducial estimates, to inform more detailed future modeling.", '3.2. Individual Pulsar Results': 'Table 1 summarizes the quantitative constraints obtained from the radial velocity profiles. Results relevant to each pulsar are discussed below, along with information relevant to the interpretation of each measurement.', '3.3. J0742 -2822': 'The shock emission of this pulsar is dominated by the blueshifted broad line, which has a low-amplitude extension to large negative velocities, suggesting that the shock is viewed at substantial inclination out of the plane of the sky. A high inclination angle would also explain the ring feature that appears to move in velocity space to the right of the images shown in Figure 1 (e.g. Barkov et al. 2020). The narrow line is weakly detected in Region (b). Due to its low S/N ( ≈ 5 ), we report a lower limit on the broad-to-narrow line  \nratio, I b /I n ≥ 3 . 0 ± 0 . 5 . Since we needed to average over a large region to detect the narrow line, the resulting spectrum covers a range of azimuthal angles, which broadens the apparent line width.  \nGiven the large asymmetry in the velocity profile, we place an upper limit on W b using the full-width-at-half-maximum (FWHM) of the blueshifted broad line, W b ≲ 160 km/s. To estimate the velocity ratio between the blue and redshifted broad lines, we first correct for the ambient ISM radial velocity measured from the narrow line. If we use the bestfit peaks of the broad lines, then we find | V blue /V red | ≈ 1 . However, this approach ignores the higher flux density of the blueshifted broad line and its extension to large negative velocities. We therefore calculate the integrated fluxes of the best-fit blue and red broad lines as functions of wavelength. To capture the velocity extrema, which are the most sensitive indicators of the shock inclination, we then evaluate the wavelengths at which the integrated fluxes of the blue and redshifted lines reach 15% and 85% of their total values, re-  \nTable 1. Bow Shock Properties  \nNOTE-The broad-to-narrow line ratio I b /I n is quoted as a lower limit when the narrow line S/N is low. The narrow line radial velocity V n is reported in the heliocentric frame; for J0742 -2822 and J1741 -2054, the apparently substantial V n are broadly consistent with Galactic differential rotation after correcting to the Local Standard of Rest. The broad-line velocity ratio | V blue /V red | is obtained by calculating the velocity extrema of the blue and redshifted broad lines using the integrated flux of each broad line (§ 3). The broad line width W b is reported as an upper limit when the blue and redshifted broad lines are highly asymmetric and not well separated; in all cases, W b is based on the apparent full-width-at-half-maximum of the broad line and does not account for the effect of inclination. Quantities on the right half of the table are drawn from the literature. Distances indicated by † are estimated from NE2001 (Cordes & Lazio 2002); pulsar transverse velocities ( V psr ⊥ ) are estimated from the proper motion and distance. The projected angular stand-off radii θ BS correspond to the most recently published values, and are based on differing shock models. References: (1) Brownsberger & Romani (2014); (2) Bailes et al. (1990); (3) Auchettl et al. (2015); (4) de Vries et al. (2022); (5) Deller et al. (2019).  \nFigure 3. Ratio of blue to redshifted broad line radial velocity | V blue /V red | vs. shock inclination angle. Negative angles correspond to inclination out of the sky and positive angles into the sky. The horizontal lines and shaded regions show inferred values and 1 σ uncertainties of | V blue /V red | for J0742 -2822 (blue), J1741 -2054 (green), and J2225 + 6535 (purple). The expected velocity ratio is estimated for an axisymmetric shock from an isotropic stellar wind (Wilkin 1996), for different values of α = V ∗ /V w , the ratio of the ISM velocity to stellar wind velocity. We constrain approximate ranges of inclination angle for each shock using the limit α → 0 (solid black curve), as appropriate for a relativistic wind.  \n<!-- image -->  \nspectively. Using this approach, which accounts for the actual partitioning of flux between the blue and redshifted lines, we find V blue = -141 ± 20 km/s and V red = 82 ± 13 km/s, yielding a velocity ratio | V blue /V red | = 1 . 7 ± 0 . 4 . For an axisymmetric shock and isotropic wind, this ratio suggests  \nan inclination angle -11 · ≲ i ≲ -2 · (see Figure 3). This inclination is consistent with the previously published upper limit of | i | ≲ 25 · (Jones et al. 2002), which was based on the projected surface brightness required to see the shock in narrowband imaging.', '3.4. J1741 -2054': 'The overall morphology of this bow shock is highly complex. At large radial velocities, the shock emission stretches across two main lobes that are oriented perpendicular to the pulsar proper motion, which cuts diagonally from the upper left to bottom right of the images in Figure 1 (Romani et al. 2010). This two-lobed structure may support previous interpretations of narrowband H α images, which found that the flatness of the shock apex was consistent with an equatorial pulsar wind (Romani et al. 2010).  \nThe narrow line is detected in three regions behind the shock (a, b, and c). Using the best-fit amplitudes of the broad and narrow components, we find I b /I n = 2 . 0 ± 0 . 1 , 2 . 3 ± 0 . 2 , and 2 . 2 ± 0 . 1 in Regions (a), (b) and (c), respectively. The radial velocity of the narrow line shifts by < 2 km/s between each region (Table 1). Given these variations at the 1 σ level, we find that I b /I n is statistically consistent with being constant across the shock.  \nWe constrain the broad line width using the FWHM averaged across both the blue and redshifted broad line components. We find W b ≈ 60 km/s and is consistent with being constant across the three velocity profiles behind the shock apex (see Table 1). These widths are also consistent at the 1 σ level with the FWHM of the single broad component detected at the shock apex in Region (d), which gives a width of 63 . 90 ± 0 . 01 km/s.  \nThe high S/N and clear separation between the blue and red sides of the velocity profiles shown in Figure 2 yield precise constraints on | V blue /V red | shown in Table 1. For an axisym-  \nmetric shock, the mean value | V blue /V red | = 1 . 34 ± 0 . 04 suggests an inclination angle -6 · ≲ i ≲ -2 · . An inclination angle out of the plane of the sky is supported by the blueshifted broad line width, which is up to 20 km/s wider than the redshifted broad line. This estimate is smaller but technically consistent with Romani et al. (2010), who found i = -15 ± 10 · based on the flatness of the shock apex and velocity shifts in split spectra.', '3.5. Guitar Nebula': "Our observation of the Guitar Nebula was substantially impacted by seeing, which blurs the shock's undulating morphology seen best in space-based images (Chatterjee & Cordes 2002; Chatterjee et al. 2004; Ocker et al. 2021; de Vries et al. 2022). The narrow line is only detected in Region (a) close to the shock apex. Further down the shock body, the narrow line blends with the redshifted broad line, likely due to the shock inclination and the broad line width decreasing downstream from the shock apex. Given that seeing was substantially larger than the shock stand-off radius (Table 1), the narrow line may include flux from preshock ionization. We thus report a lower limit I b /I n ≥ 3 . 1 ± 0 . 2 .  \nThe radial velocity profile is asymmetrically skewed negative, again suggesting an inclination out of the plane of the sky. Given this large skew, we report an upper limit on the broad line width from Region (a), W b ≲ 150 km/s, based on the FWHM of the blueshifted broad line. We constrain | V blue /V red | in the same manner as for J0742 -2822 and J1741 -2054, finding | V blue /V red | = 1 . 9 ± 0 . 1 and -11 · ≲ i ≲ -7 · . This inclination angle is larger than previous estimates because our data reveal a more significant blueshift than previous IFS observations (de Vries et al. 2022). Higher spatial-resolution IFS data (e.g. at similar resolution to de Vries et al. 2022) covering a larger portion of the shock could reveal concentric rings of H α emission where the shock is pinched, presumably due to ISM density fluctuations (Barkov et al. 2020).", '4. DISCUSSION': "The broad-line widths inferred above are all below the smallest values included in state-of-the-art simulations of non-radiative, Balmer-dominated shocks (e.g. Ghavamian et al. 2001; van Adelsberg et al. 2008), implying that the bow shocks fall into a low shock velocity regime. This conclusion is further compounded by the inferred inclination angles, which suggest that the widths reported in Table 1 are overestimates of the true, de-projected W b .  \nFollowing Romani et al. (2022), we use the relationship between H α production efficiencies and I b /I n modeled by Heng & McCray (2007) to estimate the effective shock velocity v s . Figure 4 shows I b /I n vs. v s for the three neutron star bow shocks observed with KCWI and the one other neutron star bow shock with a detection of both broad and narrow  \nlines, J1959 + 2048 (Romani et al. 2022). For J1741 -2054, I b /I n ≈ 2 , and we find v s ≈ 200 km/s. The pulsar proper motion and estimated distance give a transverse velocity of 197 km/s, which is broadly consistent with the effective shock velocity. The limits on I b /I n for J0742 -2822 and the Guitar Nebula similarly yield low velocities ( v s ≲ 160 km/s). Published non-radiative models do not treat this lowvelocity regime because SNR shocks become radiative at a transition velocity V trans ≈ 260 E 1 / 17 51 n 2 / 17 km / s , where E 51 is energy in units of 10 51 ergs and n is density in units of cm -3 (Blondin et al. 1998). Below this transition velocity, hydrogen is ionized before it reaches the SNR shock, and the broad-to-narrow H α line diagnostic cannot be measured. The closest available shock model is Ghavamian et al. (2001), which models the dependence of of I b /I n on electron-ion equilibration for a lowest shock velocity v s ≈ 250 km/s. Naively applying this model to the I b /I n inferred here would yield electron-ion temperature ratios β = T e /T p ≲ 0 . 05 , dramatically smaller than values expected for SNRs, which follow a well-established trend in which β declines with shock velocity (Ghavamian et al. 2013; Raymond et al. 2023). This discrepancy may be reconciled in part if bow shocks correspond to lower temperatures than SNRs, because I b /I n depends on the relative rates of charge transfer (which varies slowly with shock speed) and collisional excitation (which will be smaller for lower temperatures). Additionally, Ghavamian et al. (2001) assume that no Ly β is converted to H α , which may not be true for small broad line widths. Allowing for Ly β conversion could further increase I b /I n from the Ghavamian et al. (2001) predictions. Accounting for all of these effects, in addition to the potential role of Coulomb collisions, will be critical to accurately determine β for the bow shock properties we infer.  \nOur results illustrate a strong need for both non-radiative shock simulations at low velocities ( ≲ 200 km/s) and highresolution IFS observations of the full neutron star bow shock sample, in order to establish whether these sources represent a distinct regime of β and v s from SNRs. The electronion temperature ratio is a fundamental diagnostic of nonthermal electron acceleration, and of the degree to which plasma turbulence transfers energy between ions and electrons. In SNRs, the decrease in β with higher v s has been interpreted in the context of diffusive shock acceleration as a decline in electron injection efficiency with increasing shock speed (Morlino & Celli 2021), which can ultimately contribute to the differing slopes of observed cosmic ray electron and proton spectra (Aguilar et al. 2019a,b). Pulsar wind nebulae (PWNe) play a critical role in cosmic ray acceleration across a broad energy range (10s Gev - TeV). They are commonly invoked to explain the local positron excess (Abeysekara et al. 2017; Petrov et al. 2020; Di Mauro et al. 2021; Orusa et al. 2021). The role of discrete sources in ob-  \nFigure 4. Broad-to-narrow line ratio ( I b /I n ) vs. effective shock velocity for both SNRs (blue) and neutron star bow shocks (orange). SNR constraints are from the recent literature compilation of Raymond et al. (2023), and orange points indicate values inferred for the three neutron star bow shocks in this study and one bow shock from Romani et al. (2022). All four bow shocks appear to have velocities below the SNR radiative transition, which is shown by the grey dashed line for a nominal energy E = 10 51 erg and density n = 1 cm -3 (Blondin et al. 1998).  \n<!-- image -->  \ned cosmic ray spectra continues to be a topic of debate (e.g. Butsky et al. 2024), with neutron stars serving as one of the main source populations of interest (e.g. Thaler et al. 2023). PWNe with bow shocks represent an important subset of neutron stars' contribution to the cosmic ray spectrum, as reacceleration of ultrarelativistic particles at the bow shock interface can contribute a significant fraction of positrons detected at ≳ 10 s GeV (Bykov et al. 2017). While the escape of these particles can be directly traced through UV and Xray observations (Rangelov et al. 2017; de Vries & Romani 2022), the H α velocity structure offers a critical additional handle on the underlying physics, not only by resolving the  \n3D geometry of the wind and shock, but also by constraining the conditions for electron-proton equilibration and, with sufficient resolution, the underlying particle distributions. Further expanding IFS observations to the tails of these shocks may additionally constrain the role of mass-loading in shaping these shocks' larger-scale, undulating morphologies, the source of which remains debated (Morlino et al. 2015).  \nFacility: Keck:II (KCWI)  \nSoftware: KCWI DRP (Neill et al. 2023), drizzle  \n(Gonzaga et al. 2012)", 'ACKNOWLEDGMENTS': 'The authors thank John Raymond, Roger Romani, and Leon Garcia for their invaluable input on this work. SKO and MC are supported by the Brinson Foundation through the Brinson Prize Fellowship Program. SKO is a member of the NANOGrav Physics Frontiers Center (NSF award PHY2020265). Part of this work was performed at the Aspen Center for Physics, which is supprted by National Science Foundation grant PHY-2210452. Some of the data presented herein were obtained at Keck Observatory (Keck Program C422, PI: Ocker), which is a private 501(c)3 non-profit organization operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. Rosalie McGurk and Michael McDonald gave advice on observing and data reduction. Caltech and Carnegie Observatories are located on the traditional and unceded lands of the Tongva people. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the Native Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.', 'REFERENCES': "Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, Science, 358,  \n911, doi: 10.1126/science.aan4880 Aguilar, M., Ali Cavasonza, L., Ambrosi, G., et al. 2019a, PRL, 122, 041102, doi: 10.1103/PhysRevLett.122.041102 Aguilar, M., Ali Cavasonza, L., Alpat, B., et al. 2019b, PRL, 122, 101101, doi: 10.1103/PhysRevLett.122.101101 Auchettl, K., Slane, P., Romani, R. W., et al. 2015, ApJ, 802, 68, doi: 10.1088/0004-637X/802/1/68 Bailes, M., Manchester, R. N., Kesteven, M. J., Norris, R. P., & Reynolds, J. E. 1990, MNRAS, 247, 322 Barkov, M. V., Lyutikov, M., & Khangulyan, D. 2020, MNRAS, 497, 2605, doi: 10.1093/mnras/staa1601  \n```\nBlondin, J. M., Wright, E. B., Borkowski, K. J., & Reynolds, S. P. 1998, ApJ, 500, 342, doi: 10.1086/305708 Brownsberger, S., & Romani, R. W. 2014, ApJ, 784, 154, doi: 10.1088/0004-637X/784/2/154 Bucciantini, N. 2002, A&A, 387, 1066, doi: 10.1051/0004-6361:20020495 Butsky, I. S., Hopkins, P. F., Kempski, P., et al. 2024, MNRAS, 528, 4245, doi: 10.1093/mnras/stae276 Bykov, A. M., Amato, E., Petrov, A. E., Krassilchtchikov, A. M., &Levenfish, K. P. 2017, SSR, 207, 235, doi: 10.1007/s11214-017-0371-7 Chatterjee, S., & Cordes, J. M. 2002, ApJ, 575, 407, doi: 10.1086/341139\n```  \n- -. 2004, ApJL, 600, L51, doi: 10.1086/381498\n- Chatterjee, S., Cordes, J. M., Vlemmings, W. H. T., et al. 2004, ApJ, 604, 339, doi: 10.1086/381748  \nCordes, J. M., & Lazio, T. J. W. 2002, arXiv e-prints, astro.  \nhttps://arxiv.org/abs/astro-ph/0207156  \n- Cordes, J. M., Romani, R. W., & Lundgren, S. C. 1993, Nature, 362, 133, doi: 10.1038/362133a0\n- de Vries, M., & Romani, R. W. 2020, ApJL, 896, L7, doi: 10.3847/2041-8213/ab9640\n- -. 2022, ApJ, 928, 39, doi: 10.3847/1538-4357/ac5739\n- de Vries, M., Romani, R. W., Kargaltsev, O., et al. 2022, ApJ, 939, 70, doi: 10.3847/1538-4357/ac9794\n- Deller, A. T., Goss, W. M., Brisken, W. F., et al. 2019, ApJ, 875, 100, doi: 10.3847/1538-4357/ab11c7  \nDi Mauro, M., Donato, F., & Manconi, S. 2021, PRD, 104,  \n083012, doi: 10.1103/PhysRevD.104.083012  \n- Ghavamian, P., Raymond, J., Smith, R. C., & Hartigan, P. 2001, ApJ, 547, 995, doi: 10.1086/318408\n- Ghavamian, P., Schwartz, S. J., Mitchell, J., Masters, A., & Laming, J. M. 2013, SSR, 178, 633, doi: 10.1007/s11214-013-9999-0\n- Gonzaga, S., Hack, W., Fruchter, A., & Mack, J. 2012, The DrizzlePac Handbook (STSci)\n- Heng, K. 2010, PASA, 27, 23, doi: 10.1071/AS09057\n- Heng, K., & McCray, R. 2007, ApJ, 654, 923, doi: 10.1086/509601\n- Jones, D. H., Stappers, B. W., & Gaensler, B. M. 2002, A&A, 389, L1, doi: 10.1051/0004-6361:20020651\n- Kulkarni, S. R., & Hester, J. J. 1988, Nature, 335, 801, doi: 10.1038/335801a0  \nMorlino, G., & Celli, S. 2021, MNRAS, 508, 6142,  \ndoi: 10.1093/mnras/stab2972  \n- Morlino, G., Lyutikov, M., & Vorster, M. 2015, MNRAS, 454, 3886, doi: 10.1093/mnras/stv2189\n- Morrissey, P., Matuszewski, M., Martin, D. C., et al. 2018, ApJ, 864, 93, doi: 10.3847/1538-4357/aad597\n- Neill, D., Matuszewski, M., Martin, C., Brodheim, M., & Rizzi, L. 2023, KCWI DRP: Keck Cosmic Web Imager Data Reduction Pipeline in Python, Astrophysics Source Code Library, record ascl:2301.019\n- Nikoli'c, S., van de Ven, G., Heng, K., et al. 2013, Science, 340, 45, doi: 10.1126/science.1228297\n- Ocker, S. K., Cordes, J. M., Chatterjee, S., & Dolch, T. 2021, ApJ, 922, 233, doi: 10.3847/1538-4357/ac2b28\n- Orusa, L., Manconi, S., Donato, F., & Di Mauro, M. 2021, JCAP, 2021, 014, doi: 10.1088/1475-7516/2021/12/014\n- Petrov, A. E., Bykov, A. M., & Osipov, S. M. 2020, in Journal of Physics Conference Series, Vol. 1697, Journal of Physics Conference Series (IOP), 012002,  \ndoi: 10.1088/1742-6596/1697/1/012002  \n- Rangelov, B., Pavlov, G. G., Kargaltsev, O., et al. 2017, ApJ, 835, 264, doi: 10.3847/1538-4357/835/2/264\n- Raymond, J. C., Ghavamian, P., Bohdan, A., et al. 2023, ApJ, 949, 50, doi: 10.3847/1538-4357/acc528\n- Romani, R. W., Shaw, M. S., Camilo, F., Cotter, G., & Sivakoff, G. R. 2010, ApJ, 724, 908, doi: 10.1088/0004-637X/724/2/908 Romani, R. W., Slane, P., & Green, A. W. 2017, ApJ, 851, 61, doi: 10.3847/1538-4357/aa9890\n- Romani, R. W., Deller, A., Guillemot, L., et al. 2022, ApJ, 930, 101, doi: 10.3847/1538-4357/ac6263\n- Thaler, J., Kissmann, R., & Reimer, O. 2023, Astroparticle Physics, 144, 102776, doi: 10.1016/j.astropartphys.2022.102776 van Adelsberg, M., Heng, K., McCray, R., & Raymond, J. C. 2008, ApJ, 689, 1089, doi: 10.1086/592680\n- van Dokkum, P. G. 2001, PASP, 113, 1420, doi: 10.1086/323894 van Kerkwijk, M. H., & Kulkarni, S. R. 2001, A&A, 380, 221, doi: 10.1051/0004-6361:20011386\n- Verbunt, F., Igoshev, A., & Cator, E. 2017, A&A, 608, A57, doi: 10.1051/0004-6361/201731518\n- Vigelius, M., Melatos, A., Chatterjee, S., Gaensler, B. M., & Ghavamian, P. 2007, MNRAS, 374, 793,  \ndoi: 10.1111/j.1365-2966.2006.11193.x  \n- Wilkin, F. P. 1996, ApJL, 459, L31, doi: 10.1086/309939 -. 2000, ApJ, 532, 400, doi: 10.1086/308576"}