EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 665 (September 2022) A&A, 665 (2022) A35 Full HTML
Open Access
Issue
A&A
Volume 665, September 2022
Article Number A35
Number of page(s) 18
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202243531
Published online 09 September 2022

© S. Green et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1 Introduction

Stellar-wind bubbles form around massive stars because the ram pressure of the radially expanding wind pushes the interstellar medium (ISM) outwards (Weaver et al. 1977). Studying these bubbles provides a constraint on the mass-loss rate, M, of the driving star(s) by the simple argument that stronger winds drive larger stellar-wind bubbles (Gvaramadze et al. 2012). For stars moving with respect to the surrounding ISM the constraint becomes stronger because the bubble relaxes quickly to a stationary state (e.g. Mac Low et al. 1991; Arthur & Hoare 2006; Comerόn & Kaper 1998), even for subsonic motion (Mackey et al. 2016).

Runaway massive stars make up ~25% of all OB stars (Gies 1987; Blaauw 1993). When these runaway massive stars move supersonically through the ISM they develop an extensive bow shock, produced by the compression of gas and dust from the parent star and the ISM (Gull & Sofia 1979; van Buren & McCray 1988). Intense radiation from the star in question heats and ionizes this gas and dust. This makes these large-scale structures easily observable in optical nebular lines (e.g. Hα) and mid-infrared (mid-IR) emission (van Buren & McCray 1988; Meyer et al. 2014; Acreman et al. 2016). Well-known examples are are ζ Ophiuchi (Gull & Sofia 1979; Gvaramadze et al. 2012), Vela X-1 (Kaper et al. 1997; Gvaramadze et al. 2018), and BD+43° 3654 (van Buren et al. 1995; Comerόn & Pasquali 2007).

ζ Ophiuchi is one of the closest O-type stars to Earth at a distance of 135±3 pc (Gaia Collaboration 2021). Its bow shock was discovered in optical narrow-band images (Gull & Sofia 1979) and infrared (van Buren & McCray 1988) observations. It is one of the best studied cases of a bow shock around a runaway massive star, and it became one of the most iconic images from the Spitzer Space Telescope (Werner et al. 2004). Even with its proximity to Earth and its brightness, so far there have been no detailed 2D or 3D simulations of this wind bubble and bow shock.

ζ Ophiuchi is thought to have originated in a binary system and have been ejected from this system when its binary companion exploded as a supernova (Blaauw 1993; Tetzlaff et al. 2010; Neuhäuser et al. 2020). During the binary interaction before the supernova event, mass transfer occurs from the primary to the secondary star, leading to the observed rapid rotation and enrichment of He and N. The significant space velocity of ζ Ophiuchi (Neuhäuser et al. 2020) was imparted when the binary system became unbound after the primary exploded.

Toalá et al. (2016) detected diffuse X-ray emission in the vicinity of ζ Ophiuchi using Chandra observations (Weisskopf et al. 1996). They concluded that their detection was thermal X-ray emission and its morphology is consistent with expectations from radiation-hydrodynamic models of wind bubbles produced by moving stars (Mackey et al. 2015). This is to date the only detection of thermal X-ray emission from a stellar wind bubble around an isolated runaway O star.

Green et al. (2019) investigated the Bubble Nebula, NGC 7635, with two-dimensional (2D) hydrodynamic simulations, arguing that the nebula is a bow shock produced by wind-ISM interaction driven by the runaway star BD+60°2522. They made predictions of diffuse thermal X-ray emission, finding that this could be detectable with current instruments. Toalá et al. (2020) have since shown that the predicted flux is a factor of ≳ 10 too high, via a non-detection in observations with XMM-Newton. To investigate this discrepancy further we here apply our numerical model to the bow shock of ζ Ophiuchi to see if our method also overpredicts the thermal X-rays for a bubble with detected emission.

Until now most simulations of bow shocks (Mac Low et al. 1991; Comerόn & Kaper 1998; Meyer et al. 2014; Mackey et al. 2015; Green et al. 2019) and wind bubbles (Garcia-Segura et al. 1996; Freyer et al. 2003; Toalá & Arthur 2011) have been 2D because of computational requirements. Only in the last decade have 3D hydrodynamical (HD) and magnetohydrodynamical (MHD) simulations become feasible. These were first used to model bow shocks around cool stars with slow winds, where computational requirements are less severe (Wareing et al. 2007b,a; Mohamed et al. 2012; Meyer et al. 2021). Very recently this has also been applied to hot stars with fast winds (Scherer et al. 2020; Mackey et al. 2021; Baalmann et al. 2021), also for wind bubbles expanding around young stars in molecular clouds (Geen et al. 2021; Walch et al. 2015). In Mackey et al. (2021) we developed and tested a method for 3D, nested-grid, MHD simulations of bow shocks around runaway stars that are rotating and magnetised, based on established methods from heliospheric modelling (Pogorelov et al. 2004).

Building on our results from 2D hydrodynamical simulations of NGC 7635, here we investigate the bow shock of ζ Ophiuchi with 3D MHD simulations. Our aim is to determine whether models can produce X-ray emission comparable with the detected diffuse emission around ζ Ophiuchi, while simultaneously matching the infrared emission. In Sect. 2, we review the infrared observational data on the bow shock formed by ζ Ophiuchi and reanalyse Chandra archival observations to investigate the thermal diffuse X-ray emission from the nebula formed by ζ Ophiuchi. In Sects. 3-3.1, we present the numerical methods and simulation set-up. We describe our three simulations in Sects. 3.1.1-3.1.3 and describe the features of each simulation. In Sect. 3.2, the global properties of the simulation X-ray emission are described. In Sect. 4, we use the three simulations to produce synthetic infrared, Hα, soft X-ray, emission measure, and radio 6 GHz emission maps and compare them with the relevant observations. We discuss our results in Sect. 5 and conclude in Sect. 6.

Table 1

Summary of the parameters of ζ Ophiuchi.

2 Observational data

In Table 1, we show a summary of some relevant physical parameters of ζ Ophiuchi. The velocity with which ζ Ophiuchi is moving through the ISM is somewhat uncertain, and new measurements have led to some revision (Neuhäuser et al. 2020) compared with the peculiar velocity estimate, v* = 26.5 km s−1, calculated by Gvaramadze et al. (2012). To calculate the peculiar transverse velocity (vtr) in the plane of the sky, and the peculiar radial velocity (vr) we followed the method described in Green et al. (2019) but using the Gaia EDR3 data. We used an updated distance to the Galactic Centre of 8.15 ± 0.11 kpc (Bobylev & Bajkova 2021) and the circular Galactic rotation velocity of 240 km s−1 (Reid et al. 2009). The solar peculiar velocity (Schönrich et al. 2010) was taken to be (U⊙, V⊙, W⊙) = (11.1,12.2,7.3) km s−1. Using this we obtained the values of vtr and vr listed in Table 1. For vr we used the heliocentric radial velocity quoted in Gvaramadze et al. (2012) and obtain vr = −2.5 km s−1. Recently, however, Zehe et al. (2018) reviewed previous radial velocity measurements and provide an updated value corresponding to vr = 24.7 km s−1 (i.e. the vector of the peculiar velocity is inclined by an angle of about 40 degrees to the plane of sky). If correct, this radial velocity would imply a factor of 1.3 higher peculiar velocity, .

All of the parameters shown in Table 1 are relevant because they are used as initial conditions in our simulations, but there is disagreement in the literature on the values of and vr. Table 2 from Gvaramadze et al. (2012) gives estimates for the M for ζ Ophiuchi using a variety of established methods in the literature, and the values range from 10−7 to 10−9 M⊙ yr−1. Table 1 from Zehe et al. (2018) gives estimates for the vr for ζ Ophiuchi using a variety of established methods in the literature, and the values range from −35.0 to +30.0 km s−1. The two parameters and vr are important, and affect our modelling of the ζ Ophiuchi bow shock; both are able to affect the size and brightness. Until these values are constrained better, accurately modelling a bow shock is difficult.

thumbnail Fig. 1

Spitzer Space Telescope 24 μm observational data of ζ Ophiuchi in units of Mjy ster−1. The star is at the origin with coordinates in parsecs relative to the position of the star, for a distance of 135 pc.

2.1 Infrared observations

Figure 1 shows 24 μm observational infrared data from the Spitzer Space Telescope. This 24 μm image was obtained from the NASA/IPAC infrared science archive1. The bow shock was observed on 2007 April 11 as part of Program ID 30088 (PI: A. Noriega-Crespo), using the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004). At 24 μm the angular resolution is 6 arcsec. The maximum brightness of the 24 μm data is 170 MJy ster−1. We estimated the background emission observed in the 24 μm Spitzer data of 47 MJy ster−1.

The standoff distance, R0, is given by (Baranov et al. 1970) (1)

where is the mass-loss rate of the stellar wind, v∞ is the wind terminal velocity, ρISM is the density of the ISM, v* is the star’s space velocity, and cs is the sound speed in the ISM. By using the visible bow shock size in observational data it is possible to measure a value for R0 and in turn estimate the value of other variables in equation 5.1. In Gvaramadze et al. (2012) this R0 is measured from IR observations to be R0 = 0.16 pc for d = 112 pc. From this Gvaramadze et al. (2012) derived = 2.2 × 10−8 M⊙ yr−1 by using extra constraints on the local gas density from the size of the H II region around ζ Oph. Gull & Sofia (1979) derived a very similar value = 2.3 × 10−8 M⊙ yr−1 by using radio and Hα estimates of the ISM gas density.

Using the observational data presented in Fig. 1 we obtain RCD = 202 arcsec, and RFS = 405 arcsec. Using the distance of 112 pc as in Gvaramadze et al. (2012) this gives RCD = 0.11 pc and RFS = 0.22 pc. With the updated Gaia EDR3 distance of 135 pc, we obtain RCD = 0.13 pc and RFS = 0.26 pc. The relative thickness of the bow shock, , is 0.5. For all figures, including Fig. 1, we use d = 135 pc to convert angular to linear scales.

thumbnail Fig. 2

Chandra diffuse X-ray emission around ζ Ophiuchi in units of counts s−1 pixel−1 with a 3σ contour line (0.3–2 KeV) overlayed. Coordinates are in parsecs relative to the position of the star assuming the distance to the star is 135 pc. The star is at the centre of the white circle and its emission has been subtracted from the image (see Sect. 2.2 for details).

2.2 X-ray observations

To compare the observed diffuse X-rays around ζ Ophiuchi to our simulations, we used the Chandra observations presented in Toalá et al. (2016). This Chandra observation was performed on 2013 July 3 (Obs. ID 14540, PI L. M. Oskinova) with the Advanced CCD Imaging Spectrometer I-array (ACIS-I, Garmire et al. 2003) as the primary instrument. We reduced and analysed the Chandra observations using the CIAO v4.12.12 (Fruscione et al. 2006) software package with CALDB v4.9.33. The dataset was reduced using the contributed script chandra_repro, resulting in a flare-filtered and dead time-corrected exposure time of 72.1 ks.

As noted by Toalá et al. (2016), ζ Ophiuchi is such a bright X-ray emitter that it resulted in photon pile-up at the source location and also a readout streak along the ACIS-I detector column (their Fig. 2, left). A readout streak can have adverse effects on both source detection and diffuse emission analysis, so we removed the readout streak before proceeding with our analysis. We did this by identifying the streak region using the acis_streak_map CIAO task, extracting a streak background spectrum from a streak-parallel background region free of point sources using dmextract, and correcting the streak using acisreadcorr with the streak region, background spectrum and the sky X and Y position of ζ Ophiuchi.

With the readout streak removed, we produced an image of the diffuse X-ray emission by broadly following the CIAO analysis thread ‘An Image of Diffuse Emission’4, with some minor differences. We identified point sources in the field using the CIAO task wavdetect with a point spread function (PSF) map generated using mkpsfmap with an encircled counts fraction (ECF) of 0.95 in the ACIS energy range of 0.3–8 keV. The resulting source regions can be used to excise the point sources from the image and synthetically fill the source regions with realistic background using the dmfilth task. However, given the brightness of ζ Ophiuchi, we treated this source differently as counts in the wings of its PSF could easily contaminate any faint diffuse emission in its surroundings. We replaced the ζ Ophiuchi source region in the wavdetect source list with a conservative region determined for an ECF of 0.99 and multiplied by a factor of 2. Unfortunately, the large ‘hole’ left when excising ζ Ophiuchi and the counts still present in its extended PSF wings caused problems for dmfilth; the background used to fill the hole comprised mostly source counts in its extended PSF wings. Therefore, any apparent structure at or inside the ζ Ophiuchi excision region should be interpreted as nothing more than an artefact of the dmfilth task.

We then generated a counts image of the diffuse emission with 2″ bins in the 0.3–2 keV range, before creating the final exposure corrected and adaptively smoothed image using the dmimgadapt task using a minimum smoothing scale of 10″. The image region around ζ Ophiuchi is shown in Fig. 2. It is clear from the image that the counts in the wings of the ζ Ophiuchi PSF are dominating the apparent diffuse emission close to the star. However, this appears to be immersed in a larger scale faint emission extending in a roughly south-east-north-west direction. We interpret this as diffuse emission associated with the mixing region as detected by Toalá et al. (2016). We defined its extent as regions where the 0.3-2 keV count rate per pixel rises to 3σ above the average background, determined from nearby source-free regions. This is shown by the contour in Fig. 2.

To correctly derive a clean diffuse surface brightness image to compare to our synthetic X-ray emission maps would require the separation of photons from the ζ Ophiuchi PSF and diffuse emission components, which is not possible using an imaging analysis. Therefore, we performed a spectral analysis where the ζ Ophiuchi and diffuse emission could be fitted as separate components. This allowed us to quantify the diffuse flux contribution, and hence the average surface brightness within the extended emission contour.

The extended and background spectra, and weighted response files were created using the CIAo task specextract. These were fitted using XSPEC (Arnaud 1996) version 12.11.1 with abundance tables set to those of Wilms et al. (2000), photoelectric absorption cross-sections set to those of Balucinska-Church & McCammon (1992), and atomic data from ATOMDB 3.0.95. Detected point sources were masked. Additional masking was applied to the region around ζ Ophiuchi to reduce the contribution of the photons in the PSF wings to the spectra, shown by the white exclusion region in Fig. 2.

When fitting X-ray spectra of faint extended X-ray emission, careful consideration and treatment of the background is necessary. The background components present in all Chandra spectra can be split into two categories: the astrophysical X-ray background (AXB) from photons passing through the telescope optical system, and the detector non-X-ray background (NXB). To constrain the NXB, which includes high-energy particles and fluorescent emission, we used the ‘stowed’ background data and analytical model for the ACIS-I CCDs (see Bartalucci et al. 2014). We reprocessed and reprojected the stowed event file to match the ζ Ophiuchi observation using the CIAO tasks acis_process_events and reproject_events, respectively. Stowed spectra were extracted from the same extended emission and background regions on the detector that were used to constrain the NXB contributions in each.

The AXB typically comprises up to four components (Snowden et al. 2008; Kuntz & Snowden 2010), namely the unabsorbed thermal emission from the Local Hot Bubble (LHB), absorbed cool and hot thermal emission from the Galactic halo (CGH, HGH, respectively), and an absorbed power law representing unresolved background active galactic nuclei (AGN). Detailed descriptions of the components and their model parameters are given in Kavanagh et al. (2020), among others. The halo and AGN components are also absorbed by material in the Galaxy. The foreground Galactic absorption component was fixed at a column density of NH = 6.0 × 1021 cm−2 based on the Dickey & Lockman (1990) HI maps, determined using the HEASARC NH Tool6. To constrain the AXB in the region of ζ Ophiuchi, we fitted our source and background spectra simultaneously in XSPEC, linking and scaling the AXB model parameters as necessary, and using the stowed spectra fits to constrain the NXB in each. We found that the CGH component was not statistically required, so we omitted this component from the final model fits. This is unsurprising as the high foreground absorption through the Galaxy effectively absorbs the cool halo emission.

For the extended emission region, we assumed that the spectrum comprises real diffuse emission as well as a contribution from photons in the ζ Ophiuchi PSF wings. The X-ray spectral parameters were already well constrained in the analysis of Toalá et al. (2016) using both Chandra and Suzaku data, so we adopted the same model: a thermal plasma plus power law (apec+powerlaw in XSPEC) for ζ Ophiuchi with the ratio of the normalisations set to 0.93, and a thermal plasma (apec) for the diffuse emission with the temperature and photon index fixed.

The best-fit results are given in Table 2 and the spectra are shown in Fig. 3. It is clear from our results and Fig. 3 that in spite of our efforts to reduce the contamination by the ζ Ophiuchi PSF photons, they still contribute significantly to the spectrum (blue dashed lines in Fig. 3). Nevertheless, the diffuse emission component is observed (red line in Fig. 3). We determined its de-absorbed flux in the 0.3–2 keV range to be erg cm−2 s−1 which is equivalent to a luminosity of erg s−1 for d = 135pc. Taking into account the sky area of the spectral extraction region, this corresponds to an average surface brightness inside the contour in Fig. 2 of erg cm−2 s−1 arcsec−2.

3 Computational methods and simulation properties

We use the PION radiation-MHD code (Mackey et al. 2021) to model ζ Ophiuchi as an O star moving through the interstellar medium and emitting a stellar wind. PION solves the MHD equations on a computational grid in the (x, y, z) 3D space. The mass, energy, momentum densities, and magnetic vector B are defined at the centre of each computational cell, and evolved with time according to the ideal MHD equations. The ∇ · B = 0 constraint is maintained using a modified version of the Dedner et al. (2002) divergence-cleaning method with some improvements from Derigs et al. (2018), plus the Powell et al. (1999) source terms. Further details of the MHD implementation can be found in Mackey et al. (2021).

For numerical convenience, a reference frame is chosen in which the star is stationary and located at the origin (x, y, z) = (0,0,0) of a rectangular cuboid. The ISM flows past the star in the negative x-direction, interacting with the stellar wind as it does so. A passive scalar variable is used to distinguish between the ISM and wind gas. For the sake of simplicity, the ISM is assumed to be homogeneous. A split-monopole magnetic field is also set up at the star where its rotation winds the magnetic field into a Parker spiral aligned with the z-axis.

Static mesh-refinement, in the form of a multiple nested grid, was used in our simulations. This allows us to model the apex of the bow shock at high resolution, and regions further from the star at lower-resolution (e.g. the tail of the bow-shock and wake).

We used the Monte Carlo radiation transport and hydrodynamics code TORUS (Harries et al. 2019) to calculate a synthetic infrared emission map (at wavelength 24 μm, in MJy ster−1) from our simulations. The procedure used is described in Mackey et al. (2016) and Green et al. (2019), but modified to work for 3D nested-grid simulations as described in Mackey et al. (2021).

To calculate the infrared emission maps we assume a gas-to-dust mass ratio of 160 (Zubko et al. 2004), which is comprised of 70% silicate grains (Draine 2003) and 30% carbonaceous grains (Zubko et al. 1996). We assume that no polycyclic aromatic hydrocarbons survive within the H II region. For both the silicate and carbonaceous grains, we assume minimum and maximum grain sizes of 0.005 and 0.1 μm respectively. The size distribution between these limits is a power law dn/da ∝ a−q (Mathis et al. 1977), where we take q = 3.3. For the stellar spectrum we use a Kurucz (1993) spectral model with the same temperature and mass as ζ Ophiuchi. The dust is also removed wherever the temperature is more than 106 K; this is to ensure a dust-free wind.

An independent ray-tracing method to calculate synthetic images from PION simulations was developed to allow us to compare simulations with real observational data extending the uniform-grid methods used in Mackey et al. (2013). For a 3D Cartesian grid, straight lines are projected through the 3D simulation at an angle θ to one of the coordinate axes, and emission from each cell is added along the ray according to the local quantities in that cell (e.g. density, temperature). The intensity along the rays are assigned to pixels on a 2D image to compare with observational data. Local absorption within the simulation domain is not included.

Table 2

X-ray spectral fit results (see Sect. 2.2 for details).

thumbnail Fig. 3

Chandra X-ray spectra of the extended emission of ζ Ophiuchi. Left: Chandra X-ray spectra of the extended emission region (black) and the background (red). The best-fit models are shown as solid lines and the fit residuals are shown in the lower panel. Right: extended emission spectrum only with the additive model components shown. The green dash-dot-dot line represents the combined NXB components constrained using the stowed observations, the dashed blue lines show the ζ Ophiuchi emission components, the magenta dash-dot line is the AXB component, and the red solid line shows the diffuse emission component.
Table 3

Stellar and interstellar medium parameters used for simulations Z01, Z02, and Z03.

thumbnail Fig. 4

3D bow shock of ζ Ophiuchi produced by pion. Top: 3D volumetric density image from the Z01 ζ Ophiuchi simulation in Sect. 3.1.1, with the colour scale showing gas density. Bottom: Transfer function overplotted on the colour scale showing the alpha value assigned to each density value for the top image. The alpha parameter describes the opacity of the gas as a function of density.

3.1 Simulations Z01, Z02, and Z03

We used features from the infrared observations (Fig. 1) to constrain the simulation parameters, specifically the brightness and the measured distances RCD and RFS. To create a model that best matched the size and infrared brightness of ζ Ophiuchi, several different 3D simulations were run, varying the uncertain parameters v*, , and ρISM while still producing a bow shock with approximately the right size and shape. Table 3 shows the parameters used to run the three simulations presented in the following sections. All three simulations have a box length of Lbox = 7 × 1018 cm in each dimension (or 2.273 pc3), resolved with 2563 grid cells per level with three refinement levels. The focus of the mesh refinement is (x, y, z) = (1.5 × 1018, 0,0) cm. On the finest level the cell size is ∆x = 6.84 × 1015 cm, or 2.2 × 10−3 pc.

For the three simulations we impose a stellar surface field of 1G in a split-monopole configuration. This is well below the non-detection of 〈Bz〉 = 118 ± 61 G found by Bagnulo et al. (2015, Table 5), and also low enough that the star does not have a magnetosphere. For the ISM magnetic field we take |B0| = 10 μG, oriented partly in the direction of motion of the star and partly in the perpendicular plane: B0 = (5.0,8.66,1.0) μG. This is stronger than the typical magnetic field in the Galactic plane, but the ISM density that we use is also a few times above the mean density in the Plane. The motivation for certain parameters in each of the three simulations are discussed further in Sect. 4.

A projected image of the gas density is shown in Fig. 4 with rotations of 45° to the x- and the y-axis. The transfer function is chosen to highlight the contact discontinuity and the apex of the bow shock, and is plotted on the colour scale. This is not intended to mimic an observable quantity such as spectral line or continuum emission, only to show the 3D shape of the bow shock.

thumbnail Fig. 5

Simulation Z01. Left: plots of the Log10 gas density (g cm−3), slice through the (x–z) plane at y = 0. Right: plots of the Log10 gas temperature (K), slice through the (x–y) plane at z = 0. The borders of the nested-grid levels are shown in black. The star is at the origin.

3.1.1 Simulation Z01

The Z01 simulation uses v* = 26.5 km s−1, which is the case where the line-of-sight radial velocity is very small and the bow shock is seen almost edge-on. We use an ISM number density of ions, ni = 8 cm−3 and a mass-loss rate of = 1.46 × 10−8 M⊙ yr−1.

The Z01 simulation was run for 0.13492 Myr which is just over 1.5 crossing times (Lbox/v*), and 50 dynamical timescales of the bow shock (R0/V*). Numerical results from simulation Z01 are shown in Fig. 5, after 0.13492 Myr, when the bow shock has reached equilibrium (i.e. the overall shape of the bow-shock/bubble does not change and has reached a stationary state). The left panel shows log10 of the gas density (g cm−3) in a slice at y = 0 in the (x, z) plane, and the right panel shows log10 of the gas temperature for a slice at z = 0 in the (x, y) plane. The boundaries of the two refined grid-levels are shown in the right panel of the figure. Each level of refinement is a factor of two more refined than the previous level. The wind injection region is a sphere of radius 20 cells.

At the apex of the bow shock, the peak density is ρ = 5.6 × 10−23 g cm−3. The density remains at approximately this value throughout the simulation because the bow shock is dynamically stable. The temperature at the apex of the bow shock is also in equilibrium with the rest of the ISM (−8.2 × 103 K) (i.e. the forward shock is approximately isothermal). Some shock heating is evident from the faint outline of the forward shock visible in the temperature plot.

The contact discontinuity can clearly be seen in the density and temperature plots as both quantities change by a factor of −103 from the shocked ISM to the shocked wind. Inside the shocked wind region, the gas density is as low as 6.5 × 10−27 g cm−3 and as high as 3 × 10−25 g cm−3. From this simulation, we measure RTS = 0.07 pc, RCD = 0.11 pc, and RFS = 0.22 pc.

The ‘tortilla chip’ shape of the shocked wind region is due to the orientation of the ISM’s magnetic field which is primarily in the ŷ direction (B0 = [5, 8.66,1] μG). In the (x, z) plane the magnetic field lines run almost perpendicular to the outward pressure gradient of the bubble and the magnetic pressure compresses the shocked wind. On the other hand, in the (x, y) plane the field lines run parallel to the outward pressure of the bubble, which allows the shocked wind to expand into a wider ‘tortilla chip’ shape. This is shown in Fig. 6 where slices in the (y, z) plane at x = −0.3 pc (top panel) and x = −1 pc (bottom panel) of the log10 of magnitude of the magnetic field (|B|/G) are plotted. The streamlines overlaid in the two panels show the direction of the magnetic field in this plane. The thickness of the streamlines represent the strength of the field with thin lines showing weak magnetic fields and thicker lines showing stronger magnetic fields. Both panels show how the pressure from the magnetic field is counteracting the pressure from the expansion of the bubble in the z-direction, but does not restrict the bubble expansion as much in the y-direction.

3.1.2 Sim Z02

The Z02 simulation was run for 0.08967 Myr, again approximately 1.5 crossing times and 50 dynamical timescales of the bow shock. The final snapshot is shown in Fig. 7, showing the (x, y) plane at z = 0. From this simulation we can measure RTS = 0.09 pc, RCD = 0.13 pc, and RFS = 0.26 pc. The simulation differs from Z01 in that we used v* = 40 km s−1, corresponding to the case where there is a significant radial component to the star’s space velocity. We use a similar mass-loss rate to Z01, but a lower ISM number density ni = 4 cm−3 so that the ram pressure of the ISM is similar, and the bow shock has the correct size. This simulation has a bow shock with approximately the same density as Z01 at the apex, and a correspondingly similar infrared brightness (see Sect. 4.1).

3.1.3 Sim Z03

The Z03 simulation (see Table 1 for parameters), was run for 0.08918 Myr, which is just over 1.5 crossing times and 50 dynamical timescales of the bow shock. The final snapshot from the simulation Z03 is shown in Fig. 7, for which we measure RTS = 0.09 pc, RCD = 0.13 pc, and RFS = 0.26 pc. Again we show the (x, y) plane at z = 0. The bow shock is approximately the same size and shape as for Z02 because we again used v* = 40 km s−1, but here the ISM density and mass-loss rate are both twice as high as for Z02. This means that the total pressure throughout all parts of the bow shock is doubled, and the density at the apex of the bow shock is also twice that of Z02. This means that the infrared brightness (see Sect. 4.1) is significantly higher than for Z01 and Z02, and the observational data. The simulation was run because the higher pressure should give rise to higher X-ray emission in the shocked stellar wind (see Sect. 3.2).

thumbnail Fig. 6

Log10 magnetic field magnitude (|B|/G). Both plots are of Z01 in the (y–z) plane at x = −0.3 pc (top) and x = −1.0 pc (bottom). The streamlines overlaid in the panels show the direction of the magnetic field in this plane.

3.2 X-ray emission

The emissivity as a function of temperature for different X-ray bands was calculated using XSPEC v12.9.1 (Arnaud 1996) and tabulated. Solar abundances from Asplund et al. (2009) as implemented in XSPEC were used. This was used in the same way as in Green et al. (2019) to calculate the X-ray luminosity from each grid cell in simulation snapshots. Summing over all cells gives the total predicted X-ray luminosity, and ray-tracing through the domain gives synthetic X-ray surface-brightness maps. Absorption within the simulation was neglected.

Figure 8 shows the predicted soft (0.3–1 keV), medium (1–2 keV), and hard (2–10 keV) X-ray luminosity of the plasma as a function of time for the simulations Z01 (top), Z02 (middle), and Z03 (bottom). For each energy band we show an upper bound on the X-ray luminosity, corresponding to all of the X-ray emission from the whole simulation domain, and a lower bound corresponding to emission within a cube of side 0.5 × 0.5 × 0.5 pc centred on the star. This smaller cube is similar to the region observed with Chandra in Fig. 2. This allows us to best compare the values predicted from the simulation with observational values. Figure 8 plots the unabsorbed flux (erg cm−2 s−1) for a distance of 135 pc and the luminosity (erg s−1) versus time. The observed luminosity with uncertainties is shown as the grey shaded region.

For comparison, the mechanical luminosity of the wind, , is Lw ≈ (1–2) × 1034erg s−1 using the values in Table 1. The measured X-ray flux corresponds to a luminosity ≈ 105x less than the mechanical luminosity of the stellar wind. The Z01 simulation has a shocked-wind region that produces a soft X-ray luminosity (0.3-1 keV) of ~1028–1029 erg s−1 (flux 10−14–10−13 erg cm−2 s−1). For this simulation the stellar motion is almost entirely in the plane of the sky, so the significant emission from the wake behind the star is not measured by the Chandra observations. The correct comparison is therefore to the lower flux limit, which is about 10x lower than the observed flux. For Z02 we find a similar soft X-ray luminosity, with lower limit ~1028 erg s−1 and upper limit 4 × 1028 erg s−1. In this case the star has a significant line-of-sight velocity and some of the emission from the wake could be projected onto the Chandra field of view, but still the flux is too low by at least a factor of three. The Z03 simulation has the highest pressure bubble and therefore the largest soft X-ray luminosity of the three simulations: 4 × 1028 to 1029 erg s−1. Here again the star has a significant line-of-sight velocity and so much of the emission from behind the star is projected onto the Chandra field of view. The X-ray flux from this simulation is consistent with the flux measured by Chandra, although still on the low side. The luminosity and flux in the region around the apex of the bow shock stays constant over the whole life of the simulated nebula, because the flow is not subject to instability. The total luminosity throughout the simulation does have some fluctuation on long timescales, similar to what was found by Green et al. (2019), although the level of the oscillations appears to be smaller.

The X-ray emission from the hot gas in the Z01 simulation is further analysed by calculating the differential emission measure (DEM) as a function of T, defined by (2)

where ne is the electron number density in cell k and ∆Vk is the volume of cell k (Toalá & Arthur 2018). Figure 9 shows the DEM of the simulated nebula (unabsorbed) from the Z01 simulation, after 0.1349 Myr of evolution. The DEM shows a profile strongly skewed towards lower temperatures, with a power-law behaviour similar to that shown by Toalá & Arthur (2018) for stellar-wind bubbles with turbulent mixing layers and a power-law exponent of approximately −1.75 (the black line in the figure).

We can also use the X-ray emissivity in a given energy band together with the DEM profile, to calculate TA for the simulated wind bubble Toalá & Arthur (2018). Figure 10 shows the evolution of TA as a function of time for different X-ray energy bands. The soft X-ray emission shown (0.3–2 keV) has a mean temperature of about 106.4 K. Whereas, the hard X-ray emission shown (2–10 keV) has a mean temperature of about 107.25 K. The figure also shows that the medium X-ray emission between 1 and 2 keV has a mean temperature of 106.8 K. The values of TA are almost constant for the duration of the simulation, and consistent with previous work using 2D simulations (Toalá & Arthur 2018; Green et al. 2019).

thumbnail Fig. 7

Simulations Z02 and Z03. Top: plots of the Log10 gas density (g cm−3) of simulation Z02. Bottom: plots of the Log10 gas density (g cm−3) of simulation Z03. Slices through the (x–y) plane at z = 0. The star is at the origin.
thumbnail Fig. 8

Synthetic X-ray unabsorbed flux (erg cm−2 s−1) and luminosity (erg s−1) plot of ζ Ophiuchi as it evolves in time (Myr). Top: Z01. Middle: Z02. Bottom: Z03. Three X-ray bands: soft (0.3–1 keV), medium (1–2 keV), and hard (2–10 keV). The grey shaded box shows the observed Chandra luminosity and flux values. The solid lines represent the Chandra field of view and the dashed lines are the whole simulation grid.
thumbnail Fig. 9

Differential emission measure profile of the simulated nebula (unabsorbed) from the Z01 simulation, after 0.13492 Myr of evolution. The black line represents the slope of the DEM profile between 104.2 K and 107.2 K.
thumbnail Fig. 10

Log10 dominant temperature, TA (K), of the simulated nebula (unabsorbed) as it evolves in time (Myr).

4 Synthetic emission

4.1 Infrared synthetic data

Synthetic infrared emission maps of the thermal dust emission at 24 μm (MJy ster−1), calculated from the Z01, Z02, and Z03 simulations, are shown in the top, middle, and bottom panels of Fig. 11, respectively. These were calculated using TORUS as described in Sect. 3. It was estimated that the background emission in the 24 μm Spitzer image was 47 MJy ster−1, and so this constant background value was added to the emission maps. Because of the high Galactic latitude of the source, it is likely that this is foreground emission from the Local Bubble. Possible sources are thermal dust emission, or spectral lines of [O IV] and [Ne V] in the MIPS 24 μm band (Chu et al. 2009). We made synthetic images at angles from 0° to 90° (with 15° steps) between the line of sight and the velocity vector of the star.

We show the 90° image for the Z01 simulation in the top panel Fig. 11 because if vr = −2.5 km s−1 then v* is mostly in the plane of the sky. In the middle and bottom panel we show the 40° projection for the Z02 and Z03 simulations because if vr = +24.7 km s−1 then v* is no longer in the plane of the sky. The other angles from the Z01 simulation can be seen in Fig. A.1, where we show synthetic emission maps in 24 μm, 70 μm, Hα, radio 6 GHz, emission measure, and soft X-rays at angles of 0°–90°.

In the top panel of Fig. 11 it can be seen that the bulk of the infrared emission corresponds to the region between the contact discontinuity and the forward shock, at 0.11–0.22 pc. This is consistent with the position and width of the infrared arc seen in Spitzer data in Fig. 1 for the previously used distance of d = 112 pc (Gvaramadze et al. 2012), but for the new Gaia distance of d = 135 pc the bow shock is too small. The opening angle of the bow shock in the synthetic image is smaller than in the observational data as well. In the other two panels, we can see that the infrared emission is at 0.13–0.26 pc. This is consistent with the position and width of the infrared arc seen in the Spitzer data in Fig. 1. The opening angle of the bow shock in these synthetic images is similar to the observational data. The bow shock in all panels, however, is smooth, whereas the observations show some density structure. This is probably because we use a uniform ISM without any turbulent density structure.

The maximum brightness of the 24μm synthetic snapshot in the top (184 MJy ste−-1) and middle (171 MJy ster−1) panels, matches the maximum brightness of the Spitzer image: 170 MJy ster−1. The bottom panel (383 MJy ster−1) does not match the maximum brightness of the Spitzer image. The overall agreement between the Z01 and Z02 synthetic infrared emission maps and the observational infrared data is good. Due to the simplifications in our model (i.e. uniform ISM density) it is not possible to match every aspect of the observations.

thumbnail Fig. 11

Synthetic infrared emission maps of the bow shock in units of MJy ster−1 at 24 μm. Top: Z01. Middle: Z02. Bottom: Z03. The star is at the origin with coordinates in parsecs relative to the position of the star. The image is rotated to the same orientation as the Spitzer image in Fig. 1.

4.2 Hα synthetic data

Our 3D ray-tracing method described in Mackey et al. (2013) and Green et al. (2019) has been used to produce synthetic Hα emission maps of the Z01, Z02, and Z03 simulations, shown in Figs. 12 and A.1. The apex of the bow shock contains the brightest Hα emission with intensity 8 × 10−16 erg cm−2 sv1 arcsec−2.

The H II region surrounding ζ Ophiuchi is less dense than the bow shock, but is also significantly larger. The H II region has a diameter of ≈10 pc and a mean density of ni ≈ 3 cm−3 (Gvaramadze et al. 2012), for an emission measure of ≈90 cm−6 pc. This background emission has not been added to the synthetic images, and so the background emission level is determined by the ISM number density and the size of the domain (2.27 pc), which is only one-fourth of the true diameter of the H II region. As such it is the brightness of the bow shock above the background emission that should be compared with observations.

The line of sight passing through the bow shock has a length of ≈0.4 pc (Fig. 12), and therefore the bow shock should have in order to have a comparable emission measure to the H II region, and be easily visible above the background. In all three simulations this condition is fulfilled, and so the bow shock should be visible in Hα. Gull & Sofia (1979) detected the bow shock in narrow-band images, more clearly in [OIII] than in Hα, but a digital version of their image is not available for more detailed comparison. New narrow-band observations of this bow shock would be valuable to further constrain its physical properties.

4.3 X-ray synthetic data

Using the 3D ray-tracing method we produced synthetic soft X-ray emission maps, shown in Fig. 13, using the same snapshot as in Figs. 5 and 7. The panels show the final snapshot from simulations Z01, Z02, and Z03 from top to bottom. For Z01 we plot a projection along the z-axis so that the image x- and y-axes are also the simulation x- and y-axes. For Z02 and Z03, because the model includes a significant radial velocity along the line of sight, we plot a projection along a line 45° between the z- and x-axes.

The surface brightness seen in Fig. 13 should be compared with the average surface brightness of the Chandra observations in Fig. 2, of erg cm−2 s−1 arcsec−2. The X-ray emission in the top panel (Z01) near the apex of the bow shock is about 10 times too faint, in agreement with the global luminosity calculation above. Simulation Z02 (middle panel) has somewhat brighter emission, but still much fainter than the observed flux. As expected, simulation Z03 (bottom panel) is much brighter in soft X-rays than the other two simulations (by design), but still is a factor of three fainter than the observed emission.

Figure 14 shows the soft X-ray emission presented in Fig. 13 set to the same axis scales and rotated to the same direction as the Chandra image in Fig. 2. The top panel is from the Z01 simulation, the middle panel from the Z02 simulation, and the bottom panel from the Z03 simulation. The majority of the X-ray emission is coming from the apex of the bow shock at the contact discontinuity where the shocked ISM and shocked wind meet, and also following this boundary downstream from the star. This is the region of mixed wind and ISM gas, containing intermediate temperatures that cool quickly and radiate efficiently in soft X-rays. The emission from the synthetic map at the apex of the bow shock can be compared to the outer contour from the observational image, where the outline of the emitting region is also an arc that traces the inner edge of the infrared bow shock. However, the observational data shows that some of the emission may also be occurring closer to the star which is not seen in the synthetic image. From the uncertainty due to pixel pile-up from stellar emission within the white circle in the observational image, this fainter emission observed closer to the star could also have a contribution from stellar emission, although we have gone to great lengths to remove the stellar contamination as much as possible.

thumbnail Fig. 12

Synthetic Hα emission maps of the simulated nebula around ζ Ophiuchi on a linear scale in units of erg cm−2 s−1 arcsec−2. Top: Z01 with the line of sight being the y-axis. Middle: Z02. Bottom: Z03. Coordinates in parsecs relative to the position of the star (white cross).
thumbnail Fig. 13

Synthetic soft X-ray (0.3–2 KeV) emission maps of the simulated nebula around ζ Ophiuchi (unabsorbed). Top: Z01 at an angle of 90° with respect to the direction of stellar motion. Middle: Z02 at an angle of 45°. Bottom: Z03 at an angle of 45°. Coordinates are in parsecs relative to the position of the star (white cross) and the colour scale is from zero to maximum in erg cm−2 s−1 arcsec−2. The image normal direction is the z-axis of the simulation.
thumbnail Fig. 14

Synthetic soft X-ray (0.3–2 KeV) emission maps of the simulated nebula around ζ Ophiuchi (unabsorbed) rotated to the direction of motion of ζ Ophiuchi in Fig. 1 and set to the same scale. Top: Z01 at an angle of 90° with respect to the direction of stellar motion. Middle: Z02 at an angle of 45°. Bottom: Z03 at an angle of 45°. Coordinates are in parsecs relative to the position of the star (white cross) and the colour scale is from zero to maximum in erg cm−2 s−1 arcsec−2. The image normal direction is the z-axis of the simulation.
thumbnail Fig. 15

Synthetic radio 6 GHz and emission measure maps. Top: synthetic radio 6 GHz emission maps of the simulated nebula around ζ Ophiuchi on a linear scale with units MJy ster−1. Bottom: synthetic emission measure maps of the simulated nebula around ζ Ophiuchi on a linear scale in units of cm−6 pc. Both are the Z01 simulation. Coordinates are in parsecs relative to the position of the star (white cross).

4.4 Radio and emission measure synthetic data

Using our 3D ray-tracing method again, we produced synthetic 6 GHz radio emission (assuming thermal bremsstrahlung) and synthetic emission measure maps of the Z01 simulation, shown in Fig. 15. The emission measure is defined by (3)

where ne is the electron number density as a function of position along a line of sight and emission measure is traditionally measured in units of cm−6 pc. Both panels show a projection where the line of sight is the y-axis and the image coordinate axes are x and z. The maximum brightness of the radio emission is 0.035 MJy ster−1, and 300 cm−6 pc for the emission measure. In both instances the maximum brightness occurs at the apex of the bow shock.

As discussed in Sect. 4.2, ζ Ophiuchi is located within a large H II region that itself has significant emission measure (and also thermal radio emission) at the level of −100 cm−6 pc. Figure 15 shows that the background-subtracted emission measure of the bow shock is ≈150 cm−6 pc, which should be detectable. The fact that the background emission is almost as bright as the bow shock, combined with the large size of the bow shock on the sky, makes this a challenging observation, unless the bow shock emits bright non-thermal emission as well as thermal. Downstream from the apex, the background-subtracted emission measure is only ~(50–100) cm−6 pc, and would be much more challenging to detect against the non-uniform background emission from the H II region.

Until recently only one bow shock (of BD+43°3654) had been detected at radio frequencies (Benaglia et al. 2010), but in the past year a number of others have been detected: Vela X-1 (van den Eijnden et al. 2022), bow shocks in the NGC 6357 and RCW 49 regions (Van den Eijnden et al. 2022), and also non-thermal emission from the Bubble Nebula, NGC 7635 (Moutzouri et al. 2022). These detections, together with our simulation results, suggest it could be worthwhile attempting to detect the bow shock of ζ Ophiuchi at radio frequencies.

5 Discussion

5.1 Simulation motivation

Three simulations, Z01, Z02, and Z03, were run to create synthetic observations that are comparable to ζ Ophiuchi’s features. Here we describe the motivation for running each simulation, and summarise the results obtained.

Z01. the Z01 simulation was modelled using parameters based on the system parameters calculated by Gvaramadze et al. (2012), where the distance (d) to ζ Ophiuchi was estimated to be 112 pc and the stellar space velocity was 26.5 km s−1. As mentioned in Sect. 4.1, the synthetic 24 μm infrared emission map from Z01 had a comparable peak intensity and bow shock stand-off distance to the observations, but the opening angle of the bow shock was too small. The bow shock is now too small when compared to the updated distance, d = 135 pc, from Gaia EDR3 (Gaia Collaboration 2021). Figure 2 shows diffuse X-ray emission with luminosity erg s−1 from Chandra observations. This value is at least eight times more than that predicted with our Z01 simulation (≈1028 erg s−1). In addition, the majority of X-ray emission was seen to occur around the apex of the bow shock, whereas the emission is seen to be closer to the star in observations.

We discuss in Sect. 2 that the radial velocity of ζ Ophiuchi could be significantly higher than previously estimated (Zehe et al. 2018): vr = 24.7 km s−1 as oppose to vr = −2.5 km s−1. This implies a larger space velocity v* and different viewing angle of the bow shock, which requires simulations with different parameters, motivating simulations Z02 and Z03.

Z02. The Z02 simulation was set up using a new distance to the star of 135 pc from Gaia Collaboration (2021), a new peculiar transverse velocity (29 km s−1) and a peculiar radial velocity (24.7 km s−1) from Zehe et al. (2018). We approximate the total space velocity as v* = 40 km s−1. To compensate for the higher space velocity, we reduced the ISM number density so that the standoff distance would be correct. Z02 was found to have a comparable infrared intensity to the Spitzer observations. The bow shock stand-off distance and opening angle were also comparable to observations. We obtained RTS = 0.09 pc, RCD = 0.13 pc, and RFS = 0.26 pc which are consistent with the position and width of the infrared arc seen in Spitzer and WISE data. The synthetic X-ray emission maps from this simulation reduced the amount of emission from the bow-shock, potentially bringing the emission closer to the star. The calculated thermal X-ray emission from the simulated wind bubble (luminosity ≈(1-4) × 1028 erg s−1) is fainter than the Chandra diffuse X-ray observations ( erg s−1), although the simulation upper limit is only a factor of 2 below the observational lower limit.

Z03. Because the X-ray emission of Z02 was still too faint, we set up Z03 with a mass-loss rate from the star and ISM number density that were both two times higher. We expected that the two times higher density within the wind bubble would produce X-ray luminosity approximately four times higher (because it scales with ), but that the size and shape of the bow shock would be unchanged. These expectations were proven correct: Z03 was found to have a comparable bow shock stand-off distance and opening angle to observations, but its infrared intensity was about two times higher than in the Spitzer observations because of the higher gas density. The calculated thermal X-ray emission from the simulated wind bubble shows a more comparable luminosity (4 × 1028 to 1 × 1029 erg s−1) to the Chandra diffuse X-ray observations ( erg s−1).

Overall none of the simulations provides a perfect match to the available observational data. In particular the morphology of the observed X-ray emission is not matched by the simulation results, which always show a bright arc at the location of the contact discontinuity, and the overall X-ray luminosity of our simulations is a factor of a few too low except for Z03. A priori, we expected the opposite to be the case, because similar models overpredicted the thermal X-ray emission from NGC 7635 (Green et al. 2019; Toalá et al. 2020). It is tempting to argue that the inclusion of thermal conduction in our simulations would increase the X-ray emission, and perhaps move the peak of the emission closer to the wind-driving star (cf. Meyer et al. 2014), but this would exacerbate the problems with NGC 7635. Nevertheless, it would be valuable to run simulations similar to those presented here, but including anisotropic thermal conduction as implemented, for example, in PLUTO (Mignone et al. 2012; Meyer et al. 2017) to compare the predictions for synthetic X-ray emission. From our simulations we found an emission-weighted temperature of log10(TA/K) ≈ 6.4, or 0.22 keV, comparable to the temperature derived from observations (0.2 keV) by Toalá et al. (2016). Thermal conduction is expected to lower the emission-weighted temperature even further (Toalá & Arthur 2011), and so we do not expect it to have a strong effect on our results.

Another explanation for the discrepancy between our simulations and observations could be that the separation of stellar and diffuse X-ray photons in the Chandra data reduction is not perfect, and that the true diffuse emission could be fainter than our measured value. While we made every effort to quantify the contribution of the stellar photons to the spectrum from the extended emission region, it is challenging to disentangle this from the intrinsically diffuse emission. Further observations would be very valuable, for example with XMM-Newton (Jansen et al. 2001). We used the spectral parameters and flux of the diffuse emission determined in this work with the WebPIMMs tool7 to determine that only a modest exposure time is required with EPIC (Strüder et al. 2001; Turner et al. 2001) to better constrain the extent and parameters of the diffuse emission. However, contamination from the PSF wings of the ζ Ophiuchi emission will still make the analysis demanding. In the longer term, Athena (Nandra et al. 2013) and Lynx (Gaskin et al. 2019) should easily allow us to trace and characterise the diffuse emission to even lower flux levels, with the high-sensitivity, subarcsecond imaging of Lynx being particularly suited to the case of ζ Ophiuchi (see Kavanagh 2020, for a comparison of the capabilities of future X-ray observatories).

5.2 Effects of resolution and magnetic fields

In the 2D simulations of Green et al. (2019), the X-ray emission was mainly produced by the mixing of wind and ISM gas in the wake downstream from the star, but here this is not the case and the brightest X-ray emission is from the apex of the bow shock. There are three differences between the simulations presented here and those in Green et al. (2019): 3D instead of 2D with somewhat lower resolution, MHD instead of hydrodynamics, and ISM density significantly lower in this work than in the previous paper. Any of these three effects could play a role in the qualitatively different X-ray emission maps that we obtain here.

To investigate this we ran the simulation Z02 with a lower resolution of 1283 cells and three refinement levels, and also without a magnetic field (i.e. hydrodynamics, at both resolutions). The results of the four simulations are shown in Fig. 16 at approximately the same evolutionary stage, for a slice through the simulation at z = 0. The hydrodynamic simulations are prone to Kelvin-Helmholtz (KH) instability at the contact discontinuity; these can barely develop at 1283 resolution, but are much clearer at 2563 resolution. This instability is suppressed in the MHD simulations; this appears not to be a resolution effect and is expected on physical grounds (Frank et al. 1996).

Figure 17 shows the time evolution of X-ray luminosity from the four simulations from the (0.5 pc)3 region around the star (solid lines) and the full simulation domain (dashed lines). For the full domain both MHD simulations are more luminous than the two hydrodynamic simulations, despite the lack of KH instability-driven mixing. This may be because the HLL MHD solver is diffusive and generates a broad mixing region at the wind-ISM interface. For both hydrodynamics and MHD, the higher-resolution simulations have lower X-ray emission than the lower-resolution simulations, the difference being about a factor of two in luminosity. This indicates that the results have not yet converged, a sign that the numerical diffusivity of the contact discontinuity is impacting the derived emission. For the emission from the (0.5 pc)3 region around the star, the hydrodynamical and MHD simulations give very similar results at the same resolution, apart from the modulation of the hydrodynamical results by the KH instability. Figure 18 shows synthetic X-ray emission maps from the four simulation snapshots; in all cases the brightest emission is from the apex of the bow shock, although emission from the wake shows some differences.

None of the simulations shows the strong mixing of dense clumps into the shocked wind that was found for simulations of NGC 7635 (Green et al. 2019). We can speculate that this is so because the ISM here is lower density, with less inertia, easily pushed aside and not able to penetrate into the low-density wake. This would need to be tested with a systematic study of bow shocks in different density environments. Alternatively, it could simply be a lack of resolution, but we did not have the computational resources for this project to run higher resolution simulations. The issue of energy dissipation in stellar-wind bubbles and bow shocks, through turbulent mixing or thermal conduction, is currently being actively investigated from theoretical (Green et al. 2019; Geen et al. 2021; Lancaster et al. 2021a,b) and observational (Rosen et al. 2014; Olivier et al. 2021) perspectives. It is important to constrain this because the level of energy dissipation in wind bubbles crucially determines whether or not stellar-wind feedback is important for the dynamics of the ISM in galaxies.

thumbnail Fig. 16

log10 of gas density (gcm−3) of simulation Z02 for four different cases: top left with MHD and resolution 1283 cells per level, top right with MHD and 2563, bottom left with hydrodynamics and 1283, and bottom right with hydrodynamics and 2563. Slices through the (x–y) plane at z = 0 are shown, with the star at the origin.
thumbnail Fig. 17

Synthetic soft (0.3–1 keV) X-ray unabsorbed flux (erg cm−2 s−1) and luminosity (erg s−1) plot of emission from the bow shocks shown in Fig. 16 as they evolve over time (Myr). The solid lines represent the Chandra field of view and the dashed lines are the whole simulation grid, for the simulations indicated in the legend.

6 Conclusions

This paper is the continuation of a project to investigate thermal emission from stellar wind bubbles surrounding runaway O stars. We followed up our 2D study of NGC 7635 with a detailed 3D study of the bow shock of ζ Ophiuchi. To date, this has the only detection of diffuse thermal X-ray emission from a wind bubble of an isolated massive star, and ours is the first investigation of this emission with multi-dimensional simulations. Three-dimensional MHD simulations have been run to model the interaction of the star’s wind with the interstellar medium, using a range of stellar and ISM parameters appropriate for comparison with ζ Ophiuchi.

We reanalysed the Chandra archival observations to produce a diffuse emission map and estimate the total diffuse thermal X-ray flux from the shocked wind region of the bow shock. Overall our results are consistent with the published analysis by Toalá et al. (2016), and in addition we present a new map plotting the diffuse emission and a bounding contour enclosing the emission above 3σ. We found a total unabsorbed X-ray flux in the 0.32 keV band corresponding to diffuse X-ray luminosity of .

Three simulations (Z01, Z02, and Z03) were presented and compared with observational data. The Monte Carlo radiative-transfer code TORUS was used to post-process the simulations to generate synthetic 24μm emission-map predictions to compare with observational Spitzer MIPS data. A ray-tracing projection code was also used to produce synthetic emission measure, radio bremsstrahlung, Hα, and soft (0.3–2 keV) X-ray emission maps to compare with the relevant observational data.

The Z01 simulation, for the case where the star’s space velocity is entirely in the plane of the sky, was found to have a comparable mid-IR intensity and bow shock stand-off distance to observations, although the opening angle of the bow shock was too small. The maximum brightness of our 24 μm synthetic emission maps from Z01 are also comparable with the corresponding observational data. The synthetic X-ray emission maps show that the majority of X-ray emission occurs at the apex of the bow shock at the contact discontinuity. This arc-shaped emission region is not apparent in the observational data. Calculated thermal X-ray emission from the simulated wind bubble is significantly fainter (luminosity ~1028–1029erg s−1) than the Chandra diffuse X-ray observations ( erg s−1).

Recent spectroscopic observations suggest that ζ Ophiuchi could have a significant radial velocity, however, comparable to its plane-of-sky velocity. This motivated simulations Z02 and Z03, with a faster moving star (40 km s−1 compared with 26.5 km s−1 for Z01). Simulation Z02 has many similar properties to Z01, and fits the observational data somewhat better: the opening angle of the bow shock is closer to that inferred from mid-IR observations, and the morphology of the X-ray emission is more of a filled bubble rather than an arc because of the different angle to the line of sight. The total diffuse X-ray flux remains well below observations.

Simulation Z03, with a higher mass-loss rate and ISM density, has a higher total pressure and density in all parts of the bow shock, and hence more intense emission at all wavelengths. The predicted X-ray emission is much closer to the observational values (luminosity ~[0.4–1] × 1029 erg s−1), although the predicted mid-IR emission is two times brighter than observed. We have made no attempt to fine-tune the dust properties for this calculation, and so the factor of 2 disagreement in mid-IR intensity is perhaps not so significant.

Comparison between hydrodynamic and MHD simulations showed that the Kelvin–Helmholtz instability is apparent at the contact discontinuity in the former case, but absent (or too weak to see at these resolutions) when magnetic fields are included. Our results for X-ray luminosity are also somewhat resolution-dependent and not converged, in the sense that higher resolution has lower luminosity. Further work is required to assess whether instabilities would arise in even higher-resolution MHD simulations.

The shocked-wind region around ζ Ophiuchi is the closest object to Earth where bubble energetics and dissipative processes for the wind of a massive star can be investigated, and as such it is an ideal laboratory for constraining the relevant physical processes. This first numerical study of the bow shock and wind bubble around ζ Ophiuchi does not give simple answers to the important questions, but our work can be used as a basis for building more complicated models including inhomogeneous and turbulent ISM, anisotropic thermal conduction, particle acceleration and transport, and more detailed wind models. Better observational data would also be very helpful because the existing X-ray dataset has significant contamination of the diffuse emission by stellar emission.

thumbnail Fig. 18

Synthetic soft X-ray (0.3-2 KeV) emission maps of the simulated nebula around ζ Ophiuchi (unabsorbed) for the snapshots shown in Fig. 16. Top left: Simulation Z02 with MHD and at 1283. Top right: simulation Z02 with MHD and at 2563. Bottom left: simulation Z02 with HD and at 1283. Bottom right: Simulation Z02 with HD and at 2563. All images are projected at an angle of 45° to the x-axis. Coordinates are in parsecs relative to the position of the star (white cross) and the colour scale is from zero to maximum in units of erg cm-2 s-1 arcsec-2.

Acknowledgements

S.G. acknowledges funding from the Dublin Institute for Advanced Studies. We acknowledge the SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities and support (project dsast023b). J.M. acknowledges funding from a Royal Society-Science Foundation Ireland University Research Fellowship (UF140375, 14/RS-URF/3219) and a University Research Fellowship Renewal (URFR201015, 20/RS-URF-R/3712). T.J.H. is funded by a Royal Society Dorothy Hodgkin Fellowship. M.M. acknowledges funding from a Royal Society Research Fellows Enhancement Award (RGF\EA\180214, 17/RS-EA/3468). We are grateful to the referee for a detailed report that substantially improved the presentation of our work. This work has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration, and the SIMBAD database, operated at CDS, Strasbourg, France. This work also has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018), Numpy (Harris et al. 2020), matplotlib (Hunter 2007) and yt (Turk et al. 2011).

Appendix A Projected emission at all angles

Figure A.1 shows synthetic emission maps from the Z01 simulation projected at angles 0°–90° with respect to the velocity vector of the star (from left to right), and for the observational tracers (from top to bottom) 24 μm thermal dust emission, 70 μm thermal dust emission, Hα, emission measure, 6 GHz bremsstrahlung flux, and soft X-rays (0.3–2 keV).

thumbnail Fig. A.1

Snapshot from the Z01 simulation used to generate synthetic images of ζ Ophiuchi at angles from 0° – 90° with respect to the direction of stellar motion. The first row shows the 24μm emission, the second row shows the 70μm emission, the third row shows the Hα emission, the fourth row shows the emission measure, the fifth row shows the 6Ghz Radio, and the sixth row shows the soft X-ray emission. From left to right, the projection angles are 0°, 15°, 30°, 45°, 60°, 75°, and 90° with respect to the x-axis. Each image is centred on the star (white cross) and is roughly 1.0×1.0 parsec.

References

  1. Acreman, D. M., Stevens, I. R., & Harries, T. J. 2016, MNRAS, 456, 136 [NASA ADS] [CrossRef] [Google Scholar]
  2. Arnaud, K. A. 1996, in Astronomical Data Analysis Software and Systems V, eds. G. H. Jacoby, & J. Barnes, Astronomical Society of the Pacific Conference Series, 101, 17 [NASA ADS] [Google Scholar]
  3. Arthur, S. J., & Hoare, M. G. 2006, ApJS, 165, 283 [NASA ADS] [CrossRef] [Google Scholar]
  4. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  5. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  7. Baalmann, L. R., Scherer, K., Kleimann, J., et al. 2021, A&A, 650, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Bagnulo, S., Fossati, L., Landstreet, J. D., & Izzo, C. 2015, A&A, 583, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Balucinska-Church, M., & McCammon, D. 1992, ApJ, 400, 699 [Google Scholar]
  10. Baranov, V. B., Krasnobaev, K. V., & Kulikovskii, A. G. 1970, Akad. Nauk. SSSR Dokl., 194, 41 [NASA ADS] [Google Scholar]
  11. Bartalucci, I., Mazzotta, P., Bourdin, H., & Vikhlinin, A. 2014, A&A, 566, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Benaglia, P., Romero, G. E., Martí, J., Peri, C. S., & Araudo, A. T. 2010, A&A, 517, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Blaauw, A. 1993, in Massive Stars: Their Lives in the Interstellar Medium, eds. J. P. Cassinelli, & E. B. Churchwell, Astronomical Society of the Pacific Conference Series, 35, 207 [NASA ADS] [Google Scholar]
  14. Bobylev, V. V., & Bajkova, A. T. 2021, Astron. Rep., 65, 498 [NASA ADS] [CrossRef] [Google Scholar]
  15. Chu, Y.-H., Gruendl, R. A., Guerrero, M. A., et al. 2009, AJ, 138, 691 [NASA ADS] [CrossRef] [Google Scholar]
  16. Comerόn, F., & Kaper, L. 1998, A&A, 338, 273 [Google Scholar]
  17. Comerόn, F., & Pasquali, A. 2007, A&A, 467, L23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Dedner, A., Kemm, F., Kröner, D., et al. 2002, J. Comput. Phys., 175, 645 [Google Scholar]
  19. Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., & Bohm, M. 2018, J. Comput. Phys., 364, 420 [NASA ADS] [CrossRef] [Google Scholar]
  20. Dickey, J. M., & Lockman, F. J. 1990, ARA&A, 28, 215 [Google Scholar]
  21. Draine, B. T. 2003, ApJ, 598, 1026 [Google Scholar]
  22. Frank, A., Jones, T. W., Ryu, D., & Gaalaas, J. B. 1996, ApJ, 460, 777 [Google Scholar]
  23. Freyer, T., Hensler, G., & Yorke, H. W. 2003, ApJ, 594, 888 [NASA ADS] [CrossRef] [Google Scholar]
  24. Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, SPIE Conf. Ser., 6270 [Google Scholar]
  25. Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Garcia-Segura, G., Mac Low, M.-M., & Langer, N. 1996, A&A, 305, 229 [NASA ADS] [Google Scholar]
  27. Garmire, G. P., Bautz, M. W., Ford, P. G., Nousek, J. A., & Ricker, G. R., Jr. 2003, SPIE Conf. Ser., 4851 28 [NASA ADS] [Google Scholar]
  28. Gaskin, J. A., Swartz, D., Vikhlinin, A. A., et al. 2019, J. Astron. Telescopes Instrum. Syst., 5, 1 [CrossRef] [Google Scholar]
  29. Geen, S., Bieri, R., Rosdahl, J., & de Koter, A. 2021, MNRAS, 501, 1352 [Google Scholar]
  30. Gies, D. R. 1987, ApJS, 64, 545 [NASA ADS] [CrossRef] [Google Scholar]
  31. Green, S., Mackey, J., Haworth, T. J., Gvaramadze, V. V., & Duffy, P. 2019, A&A, 625, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Gull, T. R., & Sofia, S. 1979, ApJ, 230, 782 [CrossRef] [Google Scholar]
  33. Gvaramadze, V. V., Langer, N., & Mackey, J. 2012, MNRAS, 427, L50 [NASA ADS] [CrossRef] [Google Scholar]
  34. Gvaramadze, V. V., Alexashov, D. B., Katushkina, O. A., & Kniazev, A. Y. 2018, MNRAS, 474, 4421 [Google Scholar]
  35. Harries, T. J., Haworth, T. J., Acreman, D., Ali, A., & Douglas, T. 2019, Astron. Comput., 27, 63 [NASA ADS] [CrossRef] [Google Scholar]
  36. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
  37. Howarth, I. D., & Prinja, R. K. 1989, ApJS, 69, 527 [NASA ADS] [CrossRef] [Google Scholar]
  38. Howarth, I. D., & Smith, K. C. 2001, MNRAS, 327, 353 [NASA ADS] [CrossRef] [Google Scholar]
  39. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  40. Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Kaper, L., van Loon, J. T., Augusteijn, T., et al. 1997, ApJ, 475, L37 [NASA ADS] [CrossRef] [Google Scholar]
  42. Kavanagh, P. J. 2020, Ap&SS, 365, 6 [NASA ADS] [CrossRef] [Google Scholar]
  43. Kavanagh, P. J., Sasaki, M., Breitschwerdt, D., et al. 2020, A&A, 637, A12 [EDP Sciences] [Google Scholar]
  44. Kuntz, K. D., & Snowden, S. L. 2010, ApJS, 188, 46 [NASA ADS] [CrossRef] [Google Scholar]
  45. Kurucz, R. L. 1993, Phys. Scr., 47, 110 [CrossRef] [Google Scholar]
  46. Lancaster, L., Ostriker, E. C., Kim, J.-G., & Kim, C.-G. 2021a, ApJ, 914, 89 [NASA ADS] [CrossRef] [Google Scholar]
  47. Lancaster, L., Ostriker, E. C., Kim, J.-G., & Kim, C.-G. 2021b, ApJ, 914, 90 [NASA ADS] [CrossRef] [Google Scholar]
  48. Mac Low, M.-M., van Buren, D., Wood, D. O. S., & Churchwell, E. 1991, ApJ, 369, 395 [NASA ADS] [CrossRef] [Google Scholar]
  49. Mackey, J., Langer, N., & Gvaramadze, V. V. 2013, MNRAS, 436, 859 [NASA ADS] [CrossRef] [Google Scholar]
  50. Mackey, J., Gvaramadze, V. V., Mohamed, S., & Langer, N. 2015, A&A, 573, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Mackey, J., Haworth, T. J., Gvaramadze, V. V., et al. 2016, A&A, 586, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Mackey, J., Green, S., Moutzouri, M., et al. 2021, MNRAS, 504, 983 [NASA ADS] [CrossRef] [Google Scholar]
  53. Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [Google Scholar]
  54. Meyer, D. M.-A., Mackey, J., Langer, N., et al. 2014, MNRAS, 444, 2754 [NASA ADS] [CrossRef] [Google Scholar]
  55. Meyer, D. M. A., Mignone, A., Kuiper, R., Raga, A. C., & Kley, W. 2017, MNRAS, 464, 3229 [Google Scholar]
  56. Meyer, D. M. A., Mignone, A., Petrov, M., et al. 2021, MNRAS, 506, 5170 [NASA ADS] [CrossRef] [Google Scholar]
  57. Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7 [Google Scholar]
  58. Mohamed, S., Mackey, J., & Langer, N. 2012, A&A, 541, A1 [Google Scholar]
  59. Moutzouri, M., Mackey, J., Carrasco González, C., et al. 2022, A&A 663, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Nandra, K., Barret, D., Barcons, X., et al. 2013, ArXiv e-prints [arXiv:1306.2307] [Google Scholar]
  61. Neuhäuser, R., Gießler, F., & Hambaryan, V. V. 2020, MNRAS, 498, 899 [CrossRef] [Google Scholar]
  62. Olivier, G. M., Lopez, L. A., Rosen, A. L., et al. 2021, ApJ, 908, 68 [NASA ADS] [CrossRef] [Google Scholar]
  63. Pogorelov, N. V., Zank, G. P., & Ogino, T. 2004, ApJ, 614, 1007 [NASA ADS] [CrossRef] [Google Scholar]
  64. Powell, K., Roe, P., Linde, T., Gombosi, T., & de Zeeuw, D. 1999, J. Comput. Phys., 154, 284 [NASA ADS] [CrossRef] [Google Scholar]
  65. Reid, M. J., Menten, K. M., Zheng, X. W., Brunthaler, A., & Xu, Y. 2009, ApJ, 705, 1548 [NASA ADS] [CrossRef] [Google Scholar]
  66. Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25 [Google Scholar]
  67. Rosen, A. L., Lopez, L. A., Krumholz, M. R., & Ramirez-Ruiz, E. 2014, MNRAS, 442, 2701 [NASA ADS] [CrossRef] [Google Scholar]
  68. Scherer, K., Baalmann, L. R., Fichtner, H., et al. 2020, MNRAS, 493, 4172 [Google Scholar]
  69. Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 [Google Scholar]
  70. Snowden, S. L., Mushotzky, R. F., Kuntz, K. D., & Davis, D. S. 2008, A&A, 478, 615 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18 [Google Scholar]
  72. Tetzlaff, N., Neuhäuser, R., Hohle, M. M., & Maciejewski, G. 2010, MNRAS, 402, 2369 [CrossRef] [Google Scholar]
  73. Toalá, J. A., & Arthur, S. J. 2011, ApJ, 737, 100 [CrossRef] [Google Scholar]
  74. Toalá, J. A., & Arthur, S. J. 2018, MNRAS, 478, 1218 [CrossRef] [Google Scholar]
  75. Toalá, J. A., Oskinova, L. M., González-Galán, A., et al. 2016, ApJ, 821, 79 [CrossRef] [Google Scholar]
  76. Toalá, J. A., Guerrero, M. A., Todt, H., et al. 2020, MNRAS, 495, 3041 [CrossRef] [Google Scholar]
  77. Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9 [NASA ADS] [CrossRef] [Google Scholar]
  78. Turner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27 [CrossRef] [EDP Sciences] [Google Scholar]
  79. van Buren, D., & McCray, R. 1988, ApJ, 329, L93 [NASA ADS] [CrossRef] [Google Scholar]
  80. van Buren, D., Noriega-Crespo, A., & Dgani, R. 1995, AJ, 110, 2914 [NASA ADS] [CrossRef] [Google Scholar]
  81. van den Eijnden, J., Heywood, I., Fender, R., et al. 2022, MNRAS, 510, 515 [Google Scholar]
  82. Van den Eijnden, J., Saikia, P., & Mohamed, S. 2022, MNRAS, 512, 5374 [CrossRef] [Google Scholar]
  83. Walch, S., Girichidis, P., Naab, T., et al. 2015, MNRAS, 454, 238 [Google Scholar]
  84. Wareing, C. J., Zijlstra, A. A., & O’Brien, T. J. 2007a, MNRAS, 382, 1233 [NASA ADS] [CrossRef] [Google Scholar]
  85. Wareing, C. J., Zijlstra, A. A., & O’Brien, T. J. 2007b, ApJ, 660, L129 [NASA ADS] [CrossRef] [Google Scholar]
  86. Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377 [Google Scholar]
  87. Weisskopf, M. C., O’dell, S. L., & van Speybroeck, L. P. 1996, SPIE Conf. Ser., 2805, 2 [NASA ADS] [Google Scholar]
  88. Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1 [NASA ADS] [CrossRef] [Google Scholar]
  89. Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914 [Google Scholar]
  90. Zehe, T., Mugrauer, M., Neuhäuser, R., et al. 2018, Astron. Nachr., 339, 46 [NASA ADS] [CrossRef] [Google Scholar]
  91. Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321 [Google Scholar]
  92. Zubko, V., Dwek, E., & Arendt, R. G. 2004, ApJS, 152, 211 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Summary of the parameters of ζ Ophiuchi.

In the text
Table 2

X-ray spectral fit results (see Sect. 2.2 for details).

In the text
Table 3

Stellar and interstellar medium parameters used for simulations Z01, Z02, and Z03.

In the text

All Figures

thumbnail Fig. 1

Spitzer Space Telescope 24 μm observational data of ζ Ophiuchi in units of Mjy ster−1. The star is at the origin with coordinates in parsecs relative to the position of the star, for a distance of 135 pc.
In the text
thumbnail Fig. 2

Chandra diffuse X-ray emission around ζ Ophiuchi in units of counts s−1 pixel−1 with a 3σ contour line (0.3–2 KeV) overlayed. Coordinates are in parsecs relative to the position of the star assuming the distance to the star is 135 pc. The star is at the centre of the white circle and its emission has been subtracted from the image (see Sect. 2.2 for details).
In the text
thumbnail Fig. 3

Chandra X-ray spectra of the extended emission of ζ Ophiuchi. Left: Chandra X-ray spectra of the extended emission region (black) and the background (red). The best-fit models are shown as solid lines and the fit residuals are shown in the lower panel. Right: extended emission spectrum only with the additive model components shown. The green dash-dot-dot line represents the combined NXB components constrained using the stowed observations, the dashed blue lines show the ζ Ophiuchi emission components, the magenta dash-dot line is the AXB component, and the red solid line shows the diffuse emission component.
In the text
thumbnail Fig. 4

3D bow shock of ζ Ophiuchi produced by pion. Top: 3D volumetric density image from the Z01 ζ Ophiuchi simulation in Sect. 3.1.1, with the colour scale showing gas density. Bottom: Transfer function overplotted on the colour scale showing the alpha value assigned to each density value for the top image. The alpha parameter describes the opacity of the gas as a function of density.
In the text
thumbnail Fig. 5

Simulation Z01. Left: plots of the Log10 gas density (g cm−3), slice through the (x–z) plane at y = 0. Right: plots of the Log10 gas temperature (K), slice through the (x–y) plane at z = 0. The borders of the nested-grid levels are shown in black. The star is at the origin.
In the text
thumbnail Fig. 6

Log10 magnetic field magnitude (|B|/G). Both plots are of Z01 in the (y–z) plane at x = −0.3 pc (top) and x = −1.0 pc (bottom). The streamlines overlaid in the panels show the direction of the magnetic field in this plane.
In the text
thumbnail Fig. 7

Simulations Z02 and Z03. Top: plots of the Log10 gas density (g cm−3) of simulation Z02. Bottom: plots of the Log10 gas density (g cm−3) of simulation Z03. Slices through the (x–y) plane at z = 0. The star is at the origin.
In the text
thumbnail Fig. 8

Synthetic X-ray unabsorbed flux (erg cm−2 s−1) and luminosity (erg s−1) plot of ζ Ophiuchi as it evolves in time (Myr). Top: Z01. Middle: Z02. Bottom: Z03. Three X-ray bands: soft (0.3–1 keV), medium (1–2 keV), and hard (2–10 keV). The grey shaded box shows the observed Chandra luminosity and flux values. The solid lines represent the Chandra field of view and the dashed lines are the whole simulation grid.
In the text
thumbnail Fig. 9

Differential emission measure profile of the simulated nebula (unabsorbed) from the Z01 simulation, after 0.13492 Myr of evolution. The black line represents the slope of the DEM profile between 104.2 K and 107.2 K.
In the text
thumbnail Fig. 10

Log10 dominant temperature, TA (K), of the simulated nebula (unabsorbed) as it evolves in time (Myr).
In the text
thumbnail Fig. 11

Synthetic infrared emission maps of the bow shock in units of MJy ster−1 at 24 μm. Top: Z01. Middle: Z02. Bottom: Z03. The star is at the origin with coordinates in parsecs relative to the position of the star. The image is rotated to the same orientation as the Spitzer image in Fig. 1.
In the text
thumbnail Fig. 12

Synthetic Hα emission maps of the simulated nebula around ζ Ophiuchi on a linear scale in units of erg cm−2 s−1 arcsec−2. Top: Z01 with the line of sight being the y-axis. Middle: Z02. Bottom: Z03. Coordinates in parsecs relative to the position of the star (white cross).
In the text
thumbnail Fig. 13

Synthetic soft X-ray (0.3–2 KeV) emission maps of the simulated nebula around ζ Ophiuchi (unabsorbed). Top: Z01 at an angle of 90° with respect to the direction of stellar motion. Middle: Z02 at an angle of 45°. Bottom: Z03 at an angle of 45°. Coordinates are in parsecs relative to the position of the star (white cross) and the colour scale is from zero to maximum in erg cm−2 s−1 arcsec−2. The image normal direction is the z-axis of the simulation.
In the text
thumbnail Fig. 14

Synthetic soft X-ray (0.3–2 KeV) emission maps of the simulated nebula around ζ Ophiuchi (unabsorbed) rotated to the direction of motion of ζ Ophiuchi in Fig. 1 and set to the same scale. Top: Z01 at an angle of 90° with respect to the direction of stellar motion. Middle: Z02 at an angle of 45°. Bottom: Z03 at an angle of 45°. Coordinates are in parsecs relative to the position of the star (white cross) and the colour scale is from zero to maximum in erg cm−2 s−1 arcsec−2. The image normal direction is the z-axis of the simulation.
In the text
thumbnail Fig. 15

Synthetic radio 6 GHz and emission measure maps. Top: synthetic radio 6 GHz emission maps of the simulated nebula around ζ Ophiuchi on a linear scale with units MJy ster−1. Bottom: synthetic emission measure maps of the simulated nebula around ζ Ophiuchi on a linear scale in units of cm−6 pc. Both are the Z01 simulation. Coordinates are in parsecs relative to the position of the star (white cross).
In the text
thumbnail Fig. 16

log10 of gas density (gcm−3) of simulation Z02 for four different cases: top left with MHD and resolution 1283 cells per level, top right with MHD and 2563, bottom left with hydrodynamics and 1283, and bottom right with hydrodynamics and 2563. Slices through the (x–y) plane at z = 0 are shown, with the star at the origin.
In the text
thumbnail Fig. 17

Synthetic soft (0.3–1 keV) X-ray unabsorbed flux (erg cm−2 s−1) and luminosity (erg s−1) plot of emission from the bow shocks shown in Fig. 16 as they evolve over time (Myr). The solid lines represent the Chandra field of view and the dashed lines are the whole simulation grid, for the simulations indicated in the legend.
In the text
thumbnail Fig. 18

Synthetic soft X-ray (0.3-2 KeV) emission maps of the simulated nebula around ζ Ophiuchi (unabsorbed) for the snapshots shown in Fig. 16. Top left: Simulation Z02 with MHD and at 1283. Top right: simulation Z02 with MHD and at 2563. Bottom left: simulation Z02 with HD and at 1283. Bottom right: Simulation Z02 with HD and at 2563. All images are projected at an angle of 45° to the x-axis. Coordinates are in parsecs relative to the position of the star (white cross) and the colour scale is from zero to maximum in units of erg cm-2 s-1 arcsec-2.
In the text
thumbnail Fig. A.1

Snapshot from the Z01 simulation used to generate synthetic images of ζ Ophiuchi at angles from 0° – 90° with respect to the direction of stellar motion. The first row shows the 24μm emission, the second row shows the 70μm emission, the third row shows the Hα emission, the fourth row shows the emission measure, the fifth row shows the 6Ghz Radio, and the sixth row shows the soft X-ray emission. From left to right, the projection angles are 0°, 15°, 30°, 45°, 60°, 75°, and 90° with respect to the x-axis. Each image is centred on the star (white cross) and is roughly 1.0×1.0 parsec.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.