Searching for Faint X-ray Emission from Galactic Stellar Wind Bow Shocks

We present a stacking analysis of 2.61 Msec of archival Chandra observations of stellar wind bow shocks. We place an upper limit on the X-ray luminosity of IR-detected bow shocks of $<2\times10^{29}$ erg s$^{-1}$, a more stringent constraint than has been found in previous archival studies and dedicated observing campaigns of nearby bow shocks. We compare the X-ray luminosities and $L_X/L_{\rm bol}$ ratios of bow shock driving stars to those of other OB stars within the Chandra field of view. Driving stars are, on average, of later spectral type than the"field of view"OB stars, and we do not observe any unambiguously high $L_X/L_{\rm bol}$ ratios indicative of magnetic stars in our sample. We additionally asses the feasibility of detecting X-rays from stellar wind bow shocks with the proposed Lynx X-ray Observatory. If the X-ray flux originating from the bow shocks is just below our Chandra detection limit, the nearest bow shock in our sample (at $\sim$0.4 kpc with an absorbing column of $\sim10^{21}$ cm$^{-2}$) should be observable with Lynx in exposure times on the order of $\sim$100 kiloseconds.


INTRODUCTION
Stellar wind bow shocks produced by runaway OB stars (M > 8 M ) are believed to be a major source of high-energy emission in the Milky Way (Mohamed et al. 2012;del Valle et al. 2015;Toalá et al. 2016;del Valle & Pohl 2018). Non-thermal emission arises from relativistic particles (mainly electrons) being accelerated by a magnetic field at the shock front (del Valle & Romero 2012;Meyer et al. 2017) via first-order Fermi acceleration; models of this process typically use the magnetic field strength of the ambient ISM (∼few µG; Meyer et al. 2017) or equipartition arguments with respect to the kinetic power from the driving star's stellar wind (∼tens of µG; del Valle & Romero 2012). The electrons then participate in inverse Compton scattering with IR photons emitted by swept-up dust in the bow shock, scattering the IR photons up to X-ray and gamma ray energies (Peri et al. 2012;del Valle & Pohl 2018).
These models, however, predict only faint X-ray luminosities, ∼ (1 − 10) × 10 29 erg s −1 (e.g., Pereira et al. 2016;Meyer et al. 2017;De Becker et al. 2017;del Valle & Pohl 2018), making direct X-ray observations of bow shocks difficult. Although López-Santiago et al. (2012) reported the first high energy detection of the AE Aurigae bow shock, it was later shown by Toalá et al. (2017) that the X-ray emission was not spatially coincident with the IR bow shock. Toalá et al. (2016) report a detection of diffuse, thermal X-ray emission surrounding ζ Oph using Chandra, although pile-up during the observation makes direct interpretation of the Xray spectrum difficult. Non-thermal radio emission has been observed in BD +43 • 3654 (Benaglia et al. 2010), although observations with Suzaku and XMM-Newton have only yielded upper limits on the X-ray luminosity (3σ upper limit of 1.1×10 32 erg s −1 ; Terada et al. 2012;Toalá et al. 2016). Recent 2D hydrodynamical models by del Valle & Pohl (2018) suggest that only a very small fraction of the total stellar wind kinetic power (∼ 10 −5 ) is converted into nonthermal emission, and that such emission is more likely to be detected at radio wavelengths. At these faint luminosities, any hope of directly observing X-rays from bow shocks >0.5 kpc from the Sun requires prohibitively long exposure times.
In contrast, the massive stars driving the bow shocks may appear as X-ray point sources that are detectable in modest exposure times. Roughly 10% of OB stars exhibit strong (often dipolar) magnetic fields (Petit et al. 2013;ud-Doula & Nazé 2016). In these systems, charged particles from the stellar wind become trapped and channeled along magnetic field lines, leading to magnetically confined wind shocks (Babel & Montmerle 1997a,b). These stars have systematically higher X-ray luminosities and L X /L bol ratios than non-magnetic OB stars. The prototypical example is θ 1 Ori C (Gagné et al. 2005), and magnetic fields have also been detected in Tr 16-22 and Tr 16-13 in the Carina Nebula (Nazé et al. 2012). Stellar wind bow shock nebulae formation, however, is likely not affected by the magnetic surface properties of driving stars, as the radial dipolar component of the stellar B-field falls off as ∼ 1/r 2 . The stellar B-field is thus negligible at the standoff radius R 0 (typically at ∼ 0.1 − 1 pc).
We have leveraged 91 archival Chandra X-ray observations (with a total exposure time of 2.61 Msec) containing 60 IR-bright Galactic stellar-wind bow shocks. The majority (∼70%) of these bow shocks are located in apparently isolated environments, making them good candidate runaway stars. By stacking the X-ray images, we have created the deepest X-ray exposure of IR stellar wind bow shocks to date. We additionally performed a comparison of the X-ray properties of bow shock driving stars to other massive stars detected within the Chandra field of view. In Section 2 we describe the IR bow shock sample, the Chandra archival observations we use, and describe our data reduction and processing procedures (including X-ray point source detection). In Section 3 we describe our bow shock stacking analysis, and Section 4 presents a comparison of the X-ray properties of driving stars to other OB stars. We discuss our results in Section 5 and prospects for observations with future X-ray facilities, and summarize our findings in Section 6.

SAMPLE SELECTION
The largest published catalog of mid-IR stellar bow shocks contains 709 sources identified in Spitzer and WISE images of the Galactic Plane (Kobulnicky et al. 2016), and the latest data release from the citizen science initiative the Milky Way Project (MWP; Jayasinghe et al. in prep) has found an additional 282 IR stellar wind bow shocks that we include in our initial Galactic bow shocks sample.
The massive OB driving stars of these stellar wind bow shocks may also appear as X-ray point sources. The superior angular resolution of the Chandra X-ray Observatory is therefore necessary to separate potential X-ray emission from the bow shock driving star from the bow shock itself. We searched the Chandra archive for ACIS observations containing the position of at least one IR-identified bow shock within 5 of the nominal aim point. This radius was chosen to ensure reasonably small ( 2 ) and symmetric PSFs at the position of the bow shock. The search returned 91 archival observations (containing 60 unique bow shocks) with a total of 2.61 Msec exposure time. In Table 1, we summarize the observation identification numbers (hereafter referred to as "ObsIDs"), the date of the observation, which ACIS instrument was used (S or I), the nominal aim-point, the effective exposure time (see next section), and the number of bow shocks contained within 5 of the nominal aim point. Most archival observations were originally used to study massive stars, star clusters, and young stellar objects in H II regions, such as the Carina Nebula and M17. Other observations targeted pulsars and their resulting wind nebulae, supernova remnants, or follow-up observations of Swift or Fermi sources.

Data Reduction and Reprocessing
We used CIAO v4.10 and CALDB 4.7.8 to reprocess all the observations using the chandra repro task using standard reduction procedures 1 . Point source detection was performed using the CIAO task wavdetect (Freeman et al. 2002) with scales of 1 , 2 , and 4 on the 0.5-7 keV image. The absolute astrometric accuracy of Chandra observations 2 is typically ∼0.8 . We compare the Xray point source positions to 2MASS positions and use CIAO task reproject aspect to refine the positional accuracy of our X-ray images. Figure 1 shows the positional difference between Chandra and 2MASS sources as a function of off-axis angle for Obs ID 3501, a particularly crowded field; only sources within the dashed-line box (defined as off-axis angles < 5 and positional difference < 2 ) are retained in our analysis.
All X-ray sources in our preliminary list were then masked out of the X-ray images. We inspected each observation for background flares using the lc clean script; background light curves were clipped at 5σ to create good-time intervals. Each observation was then filtered using the new good time intervals (e.g., the "effective" exposure time listed in Table 1), and we restricted the energy range for each observation to 0.5-7 keV. We then use the CIAO task fluximage to create exposure maps and exposure-corrected images of each IR bow shock region. The cleaned, energy-restricted event files were used for the remainder of our analysis. Note- Table 1 is published in its entirety in machine-readable format available from the journal. A portion is shown here for guidance regarding its form and content. To create our final X-ray source lists, only sources with >5 counts within 5 of the nominal aim point were preserved. The sources were then visually examined for possible false detections. The positions of the detected X-ray sources were checked against the position of the bow shock driving stars and other known OB stars using SIMBAD. We derived distances to the bow shock driving stars and other known OB stars in the field of view (hereafter referred to as "FOV stars") using Gaia DR2 parallaxes (Gaia Collaboration et al. 2018).   Figure 2. Histogram of distances of X-ray detected FOV massive stars (black). The distances of the 5 bow shock driving stars with reliable parallaxes are indicated with downward blue arrows. FOV stars span a much larger range in distances (4.5±1.0 kpc) than driving stars (2.3±0.5 kpc).
shows the distribution of distances for driving stars and FOV stars with reasonably small distance uncertainties (i.e., d/σ > 3). We find that bow shock driving stars are systematically closer (2.3±0.5 kpc) than FOV stars (4.5±1.0 kpc). Table 3 provides a summary of all the Xray detected OB stars in our study, in order of decreasing 0.5-7 keV counts; our procedure for estimating the star's X-ray luminosity is described in Section 4. Only 5 out of the 65 unique bow shocks (∼8%) contained within a Chandra observation had their driving star detected with >5 counts.

Bow Shock "Postage Stamps"
We created "postage stamp" Chandra images for each of the 60 unique bow shocks from the exposure-corrected "fluxed" images. We masked the location of the driving star for each bow shock with a circular region that was approximately twice the radius of the PSF full width at half maximum at the detector location, so that our resulting stacked image did not contain low-level X-ray emission from the driving stars. Generally, the size of the Chandra PSF is less than the stand-off radius R 0 between the driving star and the bow shock (see, e.g., Primini et al. 2011, their Appendix A for details on modeling the Chandra PSF), so inadvertent masking of Xray counts that are genuinely associated with the IR bow shock location is unlikely. The size of the box was set to three times the stand-off radius (R 0 ) of the bow shock. Bow shocks in the Kobulnicky et al. (2016) catalog have measured position angles giving the orientation of the IR bow shock relative to Galactic north; we use the same methodology to measure the position angles of the bow shocks in the Milky Way Project using the 24 µm images.
For each image, we define a vector located at the coordinates of the driving star, with a magnitude (length) R 0 (the standoff radius) that is pointed at the position angle of the IR bow shock. Figure 3 show two example IR bow shocks from the Milky Way Project with deep (>100 ks) corresponding X-ray image and bow shock vectors.

BOW SHOCK STACKING ANALYSIS
To create the final stacked X-ray image of the bow shocks, we used the CIAO task dmregrid to rotate and resize each image vector so that the IR bow shock is located at an (x, y) position (0.5n,0.67n), where n is the size of the rescaled image. We tested several final image sizes (50, 70, and 100 pixels on a side) but our choice of final image size did not affect our results.
Our stacked X-ray image reveals an average photon flux at the location of the IR bow shock of ∼1.1×10 −7 ph s −1 . We use the CIAO task modelflux to convert this photon flux to an energy flux, assuming an absorbed power law spectral model (Γ = 2) and a typical Milky Way absorbing column of ∼ 6 × 10 21 cm −1 . The predicted 0.5-7 keV energy flux is 4.4×10 −16 erg s −1 cm −2 (unabsorbed). Assuming an average distance to a Galactic bow shock of ∼2 kpc, this flux corresponds to a total unabsorbed 0.5-7 keV luminosity of ∼2.2×10 29 erg s −1 . For comparison, the predicted "background" flux for this image (selected from a region in the stacked image that is not expected to contain any X-ray sources) is 3.9×10 −16 erg s −1 cm −2 (corresponding to a flux of 1.9×10 29 erg s −1 at a distance of 2 kpc).
Although ∼70% of the IR bow shocks in the Kobulnicky et al. (2016) catalog are in "isolated" environments, ∼34% of the total X-ray exposure time in our sample was targeted at star-forming regions that likely contain significant diffuse X-ray emission. Spectral fitting of Chandra observations of five massive starforming regions  found the diffuse emission to be dominated by at least two thermal components, one at kT ∼0.2-0.6 keV and one with kT ∼0.5-0.9 keV, with an integrated surface brightness of ∼ (3−500)×10 30 erg s −1 pc −2 . The physical size of a "postage stamp" that is 70 pixels on a side (with a plate scale of 0.492 per pixel) at an average distance of 2 kpc is ∼0.32 pc on a side (∼0.10 pc 2 ). The total X-ray luminosity in each "postage stamp" is therefore expected to be (0.3-50)×10 29 erg s −1 , in agreement with our background estimate above. We therefore do not find any evidence of excess X-ray emission in our stacked image that can be attributed to stellar wind bow shocks.
We perform an additional stacking analysis using only the 42 bow shocks from the Kobulnicky et al. (2016) in "isolated" environments and the 19 bow shocks from the Milky Way Project that were not associated with a Galactic bubble. Using only the "isolated" bow shocks (∼67% of our sample by number) yields a total exposure time of 1.73 Msec. We follow the same stacking analysis as for our full sample; at the location of the IR bow shocks the average photon flux is 1.2×10 −7 , yielding an energy flux of 4.9×10 −16 erg s −1 cm −2 (∼2.4×10 29 erg s −1 at 2 kpc). The predicted background flux is consistent with the full bow shock sample; we again find no evidence for excess X-ray emission at the location of the IR bow shock. Figure 4 show the stacked X-ray image of the isolated bow shocks (both in "raw" photon counts, and an image that has been smoothed for display purposes); the full sample images look similar.

X-RAY CHARACTERISTICS OF BOW SHOCK DRIVING STARS
We additionally investigated whether the characteristics of X-ray detected bow shock driving stars differed significantly from X-ray detected FOV stars. The X-ray spectra of many of the massive stars listed in Table 3 have been previously analyzed in the literature; here, we only aim to approximate their X-ray luminosities assuming a simplified X-ray spectral model.
We assume a typical X-ray spectrum composed of an absorbed thermal plasma (tbabs*apec in XSPEC), with  Note-Known X-ray binaries were excluded from our sample. a Spectral types were taken from SIMBAD.
b Distances were derived using the Gaia second data release (Gaia Collaboration et al. 2018).
c For sources detected in multiple ObsIDs, we report the largest off-axis angle. N H ∼ 6 × 10 21 cm −2 and a plasma temperature of ∼0.9 keV (see, e.g. Nazé 2009, for a detailed study of Xray spectra of massive stars with XMM-Newton), and use PIMMS 3 to convert the observed count rate for each source into a 0.5-7 keV flux. This flux is then converted into a luminosity for stars with Gaia distance estimates.
In Table 4, we compare our estimated luminosities for four well-studied massive stars to the values available in the literature. The most significant deviations occur where the Gaia distances disagree with those assumed in earlier studies; we therefore scale the literature values to the Gaia distances for a more direct comparison.
In general, our approximate luminosities agree with the literature values within a factor of a few. We therefore do not anticipate systematic differences in luminosity across our full massive star sample. We find no evidence of systematic differences in Xray luminosity between the FOV stars and the five bow shock driving stars. The average 0.5-7 keV luminosity of stars in our sample is ∼ 5 × 10 31 erg s −1 , while the average bow shock driving star luminosity is ∼ 3 × 10 31 erg s −1 . We also investigated whether the stars in our sample follow the well-known L X ∼ 10 −7 L bol relationship (Chlebowski 1989;Berghoefer et al. 1997;Nazé et al. 2011). We obtain the spectral type of each star from SIMBAD. We adopted model luminosities of massive stars (see Crowther 2007, their Table 2) to infer L bol for different spectral types of stars in our sample. In , constructed from 61 individual Chandra exposures totaling 1.73 Msec for bow shocks in isolated environments. The raw figure is shown on the left, while the right panel has been smoothed for display purposes only. The red "X" shows the location of the central driving stars and the red arrow indicates the location of the IR bow shocks. No excess X-ray emission is detected at the location of the IR bow shocks. Notea Literature luminosity has been scaled to the distance reported in Table 3.
cases where SIMBAD lists multiple spectral type for a star, we assume the earliest type. Figure 5 shows the distributions of L X and L X /L bol ratios for our sample. We find an average L X /L bol ∼ 1.1 × 10 −7 , in excellent agreement with previous studies. None of the stars in our sample exhibit unusually high L X /L bol ratios expected from magnetic OB stars (Petit et al. 2013;ud-Doula & Nazé 2016;Gagné et al. 2005;Nazé et al. 2012).
We note that all of the X-ray detected bow shocks driving stars are of "late-types" (e.g., later than O6), while the FOV star sample contains both early-(e.g., earlier than O5.5) and evolved Wolf-Rayet stars. The number of X-ray detected Wolf-Rayet and early-O type stars is small, and exhibit L X /L bol ratios that are statistically consistent with our sample average.

DISCUSSION
Previous attempts to study the X-ray emission from stellar wind bow shocks have yielded only upper limits on the order of L X 10 30 erg s −1 using archival XMM-Newton observations (Toalá et al. 2017(Toalá et al. , 2016 and dedicated observations (De Becker et al. 2017) of individual bow shocks. Our stacking analysis adds a new, more stringent upper limit, L X 2 × 10 29 erg s −1 .
The non-thermal radiation model developed by De Becker et al. (2017) assumes X-rays are generated via inverse Compton scattering at a magnetically-confined shock front. The swept-up dust in the bow shock is heated by photons from the driving star, and emits the reprocessed energy in the IR. These IR photons then interact with the energetic electrons at the shock front and are scattered into the X-ray and γ-ray regime. The IR luminosity therefore some fraction χ IR of the bolometric. luminosity of the dust at the shock front (e.g., van Buren & McCray 1988). The predicted X-ray luminosity additionally depends on the magnetic field strength B at the shock front and the power in relativistic electrons (expressed as a fraction χ rel of the total power available in the bow shock). When upper limits of χ IR = χ rel = 1 and B = 10 −4 G are assumed (measurements of the ambient ISM magnetic field are consistent with ∼ 2 − 6µG, Meyer et al. 2017;Harvey-Smith et al. 2011;Fiedler & Mouschovias 1993;Troland & Heiles 1986), their model predicts an integrated X-ray luminosity of ∼ 4 × 10 29 erg s −1 -a factor of ∼2 above our upper limit. Our non-detection indicates that, if X-ray photons are indeed being produced in IR bow shocks, the production mechanism is much less efficient than the models derived by De Becker et al. (2017).
State-of-the art hydrodynamical codes now include particle acceleration, allowing X-ray production in bow shock nebulae to be self-consistently modeled (van Marle et al. 2018). For example, the detailed, 2D hydrodynamical treatment by del Valle & Pohl (2018), which includes diffusion of particles and advection of energy out of the bow shock region, predicts X-ray luminosities on the order of a few ×10 28 erg s −1 consistent with our upper limit, with potentially higher luminosities (∼ 10 30 erg s −1 ) produced in the γ-ray regime. Although these models predict higher luminosities in the γ-ray than the X-ray, systematic searches for γ-ray emission from bow shocks have similarly yielded only upper limits (Schulz et al. 2014;H. E. S. S. Collaboration et al. 2018).
It is unlikely that X-ray emission from stellar wind bow shocks, if it is indeed produced through the mechanisms previously proposed, is detectable by present Xray telescopes. Even the future Athena X-ray Observatory (Nandra et al. 2013) will not possess the sensitivity and angular resolution required to detect such faint Xray emission from IR bow shocks.
The proposed NASA X-ray flagship mission Lynx, however, will have 50× greater sensitivity than Chan-dra with similar angular resolution over a significantly wider field of view 4 . We use the Simulated Observations of X-ray Sources 5 (SOXS) package in Python to determine the feasibility of detecting faint X-ray emission from bow shocks with Lynx.
For consistency with our assumptions in Section 3, we assume the non-thermal bow shock spectrum follows a power law with Γ = 2 over an energy range of 0.1-10 keV, subject to absorption due to intervening gas and dust in the Milky Way. We used the nearest bow shock in our sample: driving star HD 34078, at 0.41 kpc. The bow shock nebula is formed at R 0 = 9.9 from its driving star, with a column density along the line of sight N H = 1.2 × 10 21 cm −2 . The spectrum is then renormalized so that the 0.5-7 keV flux is just below our Chandra detection limit at the distance of HD 34078.
Photons are drawn from the underlying spectrum and spatially distributed according to a β-model surface brightness profile. The X-ray surface brightness S as a function of radius r is given by where S 0 is the core surface brightness, r C is the core radius, and β is the slope parameter. We assume r c ∼ 3 and β = 1 (with an ellipticity of 0.3) to approximately match the typical appearance of the 24 µm emission. After the photons are assigned a sky location, our simulated bow shock is then "observed" in the 0.2-7 keV energy range with the Lynx high-definition X-ray imager (HDXI) using the SOXS instrument simulator for a given exposure time. The instrument simulator additionally adds a Galactic foreground and instrument noise to our image; we refer the reader to the SOXS User's Guide for further details. Figure 6 shows the results of a 100 ks Lynx exposure of the HD 34078 bow shock just below our Chandra detection limit. An X-ray excess clearly evident in the radial surface brightness profile; the bow shock is detectable with a signal-to-noise ratio of ∼10. The same bow shock placed farther away (∼2 kpc) is detectable with a signal-to-noise ratio of ∼3 in ∼500 ks.

CONCLUSIONS
We have performed a stacking analysis leveraging 2.61 Msec of archival Chandra observations to search for faint X-ray emission from Galactic stellar wind bow shocks. Our stacked image shows no evidence for excess emission at the expected bow shock location to a flux limit of ∼ 4.4 × 10 −16 erg s −1 cm −2 which corresponds to a luminosity of ∼ 2 × 10 29 erg s −1 assuming a distance of ∼2 kpc (typical of the bow shocks in our sample). This provides the most stringent observational upper limit in the X-ray regime to date.
We assess the plausibility of detecting faint X-ray emission from bow shocks with the proposed future Xray mission Lynx. If stellar wind bow shocks are indeed producing X-rays at just below our Chandra detection limit, the least obscured cases (e.g., the nearest, or runaways at higher Galactic latitudes) should be detectable by Lynx with modest exposure times.  Figure 6. Simulated 100 ks Lynx observations (0.2-7 keV) of the HD 34078 (at a distance of 0.41 kpc) stellar wind bow shock just below the flux detection limit set by Chandra. The left panel shows the smoothed "raw" counts image. The color bar shows the pixel value, ranging from 0 (black) to ∼ 3 (white). The right panel shows the background-subtracted 0.2-7 keV surface brightness profile extracted from the simulated image, showing the excess X-ray emission is clearly detected.