- Split View
-
Views
-
Cite
Cite
Jonathan Mackey, Norbert Langer, Vasilii V. Gvaramadze, Dynamics of H ii regions around exiled O stars, Monthly Notices of the Royal Astronomical Society, Volume 436, Issue 1, 21 November 2013, Pages 859–880, https://doi.org/10.1093/mnras/stt1621
- Share Icon Share
Abstract
At least 25 per cent of massive stars are ejected from their parent cluster, becoming runaways or exiles, travelling with often-supersonic space velocities through the interstellar medium (ISM). Their overpressurized H ii regions impart kinetic energy and momentum to the ISM, compress and/or evaporate dense clouds, and can constrain properties of both the star and the ISM. Here, we present one-, two- and (the first) three-dimensional simulations of the H ii region around a massive star moving supersonically through a uniform, magnetized ISM, with properties appropriate for the nearby O star ζ Oph. The H ii region leaves an expanding overdense shell behind the star and, inside this, an underdense wake that should be filled with hot gas from the shocked stellar wind. The gas column density in the shell is strongly influenced by the ISM magnetic field strength and orientation. Hα emission maps show that H ii region remains roughly circular, although the star is displaced somewhat from the centre of emission. For our model parameters, the kinetic energy feedback from the H ii region is comparable to the mechanical luminosity of the stellar wind, and the momentum feedback rate is >100 times larger than that from the wind and ≈10 times larger than the total momentum input rate available from radiation pressure. Compared to the star's eventual supernova explosion, the kinetic energy feedback from the H ii region over the star's main-sequence lifetime is >100 times less, but the momentum feedback is up to 4 times larger. H ii region dynamics are found to have only a small effect on the ISM conditions that a bow shock close to the star would encounter.
INTRODUCTION
About 25 per cent of massive stars are classified as isolated (Gies 1987; Blaauw 1993), meaning they are not currently part of a star cluster or association. A number of recent studies (e.g. de Wit et al. 2005; Gvaramadze & Bomans 2008; Schilbach & Röser 2008; Gvaramadze et al. 2011, 2012a) have shown that all but a handful of nearby O stars are very likely to be either classical runaway stars (with space velocity v⋆ ≥ 30 km s−1) or to at least have a significant peculiar velocity directed away from a star cluster that could plausibly have been their birth place. We will refer to all such stars as exiles: stars that have been ejected from their place of birth. Nearby stars such as ζ Oph and α Ori (Betelgeuse) belong in this category; indeed by definition all massive stars closer to us than the nearest region of massive star formation must be exiles.
Bow shocks from these stars have been studied for many years and in some depth; see Villaver, Manchado & García-Segura (2012), Mohamed, Mackey & Langer (2012) and references therein for bow shock modelling around cool stars, and e.g. van Marle et al. (2006), Arthur & Hoare (2006), Comeron & Kaper (1998) and references therein for bow shocks around hot stars. H ii regions around exiled massive stars have received comparatively less attention, probably because at increasingly hypersonic velocities the H ii region produces ever-weaker shocks whereas the bow shock gets stronger.
Analytic and one-dimensional (1D) numerical predictions were made by Rasiwala (1969) for the shapes and properties of H ii regions around moving stars, with the simplifying assumption that recombinations do not occur. More detailed numerical models (including recombinations but no hydrodynamics) by Thuan (1975) predicted a somewhat cometary shape for the H ii region with a broad recombination front in the wake behind the star, and noted that the H ii region reaches a steady state in a time shorter than the lifetime of a massive star. An analytic model by Raga (1986) showed that the H ii region should become more cometary as the stellar space velocity increases. Raga et al. (1997a) made the first two-dimensional axisymmetric radiation-hydrodynamics simulations of an H ii region around a hypersonically moving star. The H ii region shape was found to agree broadly with predictions, being measurably (but not strongly) aspherical for the parameters chosen in the simulation. In addition, the hydrodynamic expansion of the H ii region was found to drive a weak, outward-moving, conical shell of overdense gas trailing behind the star, which could be observable in infrared dust emission. Other research has focused on the H ii regions produced by stars encountering a dense molecular interstellar medium (ISM; Tenorio Tagle, Yorke & Bodenheimer 1979; Mac Low et al. 1991; Arthur & Hoare 2006), or stars moving from a dense cloud into lower density ISM at velocities v⋆ ≤ 12 km s−1 (Franco et al. 2007), or slowly moving young stars in proto-star clusters (Peters et al. 2010; Dale & Bonnell 2011).
Chiang & Rappaport (1996) studied H ii regions around stationary and moving supersoft X-ray sources. They assumed isothermal ionized gas, and also ignored hydrodynamics, but included a non-equilibrium ionization calculation with an approximate treatment of the relative motion between the ISM and the radiation source. They found that recombination-line emission-maps of H ii regions are slightly aspherical for moving stars with v⋆ = 30 km s−1, having a sharper upstream edge and more extended downstream edge. This distortion becomes more dramatic with v⋆ = 100 and 300 km s−1. The general features of their results should also be found for H ii regions around moving main-sequence O stars, although the softer radiation spectrum of O stars means that ionization fronts (I-fronts) will be much thinner. The distortion they find for v⋆ = 100 km s−1 is much larger than was found by Raga et al. (1997a), a consequence of the different radiation spectra as well as the very different Strömgren radii of the modelled H ii regions.
The velocity range between transonic and hypersonic stellar motion (10 km s−1 < v⋆ < 50 km s−1) has not yet been explored for exiled stars moving through the diffuse ISM. The distribution of stellar space velocities for exiles is, however, weighted strongly towards this range (Eldridge, Langer & Tout 2011). In addition, many O stars currently in clusters will become exiles later in life (e.g. because of a binary supernova explosion; Blaauw 1961), so the observed 25 per cent fraction of isolated O stars is a lower limit to the total fraction that become exiles (cf. Eldridge et al. 2011). It is shown here that in such cases the energy and momentum imparted to the ISM from the H ii region can be comparable to, or larger than, that from stellar winds, at least for moderate-velocity stars such as ζ Oph with v⋆ = 26.5 km s−1 (Gvaramadze, Langer & Mackey 2012b).
The aims of this work are as follows:
to investigate the gas dynamics of an H ii region produced by an exiled star moving with supersonic (but not hypersonic) velocities through the warm neutral medium (WNM), taking parameters similar to ζ Oph as an example case;
to predict the effects of the H ii region on the ISM, in terms of kinetic energy and momentum feedback and
to assess the effects of the H ii region gas dynamics on the ISM near the star, to deduce the properties of the ISM that a stellar wind bow shock encounters.
We address these questions by modelling H ii regions with one-, two- and three-dimensional (1D, 2D and 3D) magnetohydrodynamic (MHD) simulations including non-equilibrium photoionization. These are the first 3D simulations of H ii regions around supersonically moving exiled massive stars, and also the first to include an ISM magnetic field. We make simulated observations to investigate the shape of the H ii region, the structure of the shell around it, its emission properties, the effects of an ISM magnetic field and the stability of the I-front to perturbations. In addition, we investigate how the density and velocity of the ISM at the star are changed by the presence of the H ii region, and discuss the potential consequences for the stellar wind bow shock (not modelled here). Our motivation is partly to test the analytic model of Gvaramadze et al. (2012b) for constraining the mass-loss rates of massive stars by simultaneous observation of their H ii region and bow shock. This model ignored possible (magneto)hydrodynamic complications, for example the dynamic response of gas to photoionization and density inhomogeneity in the ISM.
Section 2 describes the simulation code and suite of simulations we have run. Section 3 describes 1D simulations of H ii regions around stars moving with different velocities from sonic to hypersonic, comparing the position of the upstream I-front to analytic predictions. Results from 3D simulations are presented in Section 4, where the general morphology of the H ii regions is discussed, as well as the kinetic energy and momentum feedback, and the properties of the ISM near the star. The results are discussed further in Section 5, and our conclusions are presented in Section 6. Effects of spatial and temporal numerical resolution are discussed for a range of stellar space velocities in Appendix A. Effects of limited spatial resolution in 3D simulations are studied in Appendix B with higher resolution 2D simulations. Equations for heating and cooling rates are listed in Appendix C.
NUMERICAL METHODS AND SIMULATION SETUP
We use the radiation-MHD code pion (Mackey & Lim 2010, 2011; Mackey 2012) for the simulations presented here, solving either the Euler or ideal MHD equations on a uniform rectilinear grid, coupled to a microphysics integrator to solve for the non-equilibrium neutral fraction of hydrogen, yn, and the ionization-dependent heating/cooling rates. Simulations are run on 1D, 2D and 3D Cartesian grids, with spatial derivatives set to zero for the dimensions not calculated (slab symmetry). The finite-volume integration scheme (including heating and cooling source terms) is algorithm A3 in Mackey (2012) (based on Falle, Komissarov & Joarder 1998) and is second-order-accurate in space and time and dimensionally unsplit. Radiation transfer is solved in the on-the-spot approximation, solving only for direct radiation from point sources with a raytracer that calculates the column density of neutral hydrogen and total gas density from the source to every cell. The radiation flux obeys an inverse-square law regardless of the grid dimensionality. For the 1D simulations, we model a line along the star's direction of motion, passing through the star. In two dimensions this is a plane containing the star and its velocity vector, and in three dimensions the full space is modelled. The motivation for this in 1D is that it allows us to model the flow through the upstream I-front realistically (i.e. with correct boundary conditions).
Microphysics
Near the star there is strong photoionizing radiation, and also an elevated FUV radiation field that we model following Henney et al. (2009) for the extra FUV heating rate. In photoionized gas there is strong cooling from forbidden-line emission from ionized C and O. We do not include stellar winds so the maximum gas temperature is T ≈ 12 000 K.
To calculate τν, rays are traced using the short characteristics method with the interpolation scheme proposed in Mellema et al. (2006). Two raytracing are performed per timestep, one of which uses time-centred values of |$n_{\small h}$| and yn to ensure the overall scheme is second-order-accurate in time (Mackey 2012). Algorithm A3 in Mackey (2012) is an explicit finite-volume integration scheme, so the timestep must be limited by the velocity of any I-fronts in the simulation (in addition to the hydrodynamic timestep limit). It has been shown that a sufficient criterion to use with A3 is that the timestep, Δt, should satisfy |$\Delta t \le 0.25/\skew3\dot{y}_n$| (criterion dt02 in Mackey 2012), and this is used throughout this paper (for further discussion on timestep criteria see Appendix A).
These two equations are coupled in that most of the rates depend on T and xe. They are integrated together using high-order backwards differencing and adaptive substepping to fixed relative and absolute error tolerances with the cvode numerical integration library (Cohen & Hindmarsh 1996). It has been demonstrated (Bovino et al. 2013) that high-order integration methods which integrate the energy equation together with the rate equations are significantly better than uncoupled schemes.
Volumetric heating/cooling rates as a function of T are plotted in Fig. 1 for different gas densities, with and without an ionizing radiation field. The curve for photoionized gas is largely independent of density and distance from the star, except that the equilibrium temperature increases with source attenuation because of spectral hardening. For neutral gas far from any ionizing sources, the heating rate scales with density whereas cooling scales with density squared. This means the equilibrium temperature decreases from about 7000 K at |$n_{{\small h}}=0.1\,\mathrm{cm}^{-3}$| to 55 K at |$n_{{\small h}}=100\,\mathrm{cm}^{-3}$|. The main coolants at these temperatures are far-infrared fine-structure metal lines. At T ≳ 104 K, collisional excitation of neutral H becomes dominant in neutral gas, and C and O optical forbidden lines in ionized gas.
Simulation setup
The simulations here are designed to be relevant for the nearby runaway O9.5 V star ζ Oph (Martins, Schaerer & Hillier 2005), discussed in Gvaramadze et al. (2012b). A mean number density of |$n \approx 1.1 n_{{\small h}} = 3\,\mathrm{cm}^{-3}$| was derived for the ISM in its vicinity; here we use a similar value |$n_{{\small h}}=2.5\,\mathrm{cm}^{-3}$|, or ρ0 = 5.845 × 10−24 g cm−3. The initial gas pressure is pg = 3.795 × 10−13 dyne cm−2, corresponding to a temperature of T = 1000 K, and the initial H+ fraction (1 − yn) = 0.0021 is the equilibrium value at this temperature and density. A random adiabatic perturbation with maximum amplitude of 25 per cent in pressure is applied to each cell in multidimensional simulations. A point source of ionizing and FUV photons (hereafter ‘the star’) is placed at the origin, and the ISM flows past this at a velocity of v⋆ = 26.5 km s−1, again motivated by the case of ζ Oph.
The star has an ionizing photon luminosity Q0 = 3.63 × 1047 s−1, distributed according to a blackbody spectrum for an effective temperature T⋆ = 30 500 K, appropriate for an O9.5 V star. In addition it has an FUV photon luminosity Lfuv = 3 × 1047 s−1, implemented as in Henney et al. (2009); this acts in addition to the background interstellar radiation field, which is assumed to have the standard parameters (taken from Wolfire et al. 2003).
Simulations have been run at a number of different spatial resolutions, with zone-size decreasing by factors of 2 from Δx = 0.32 to 0.005 pc (identified by name as r1 to r7, with r1 having the coarsest resolution). In 3D only r1 and r2 were possible with available computing resources. To estimate resolution effects, models with resolution r3 and r4 were run in 2D (see Appendix B), and up to resolution r7 in 1D (see Section 3). The highest resolutions (r6 and r7) have grid zones that are optically thin to ionizing photons even when fully neutral (the mean free path for a 13.6 eV photon in neutral gas at |$n_{{\small h}}=2.5\,\mathrm{cm}^{-3}$| is ≈0.021 pc, comparable to the zone size for r5). Spatially resolving the I-front in 3D will be difficult to achieve without adaptive spatial resolution (cf. Cantalupo & Porciani 2011), but the results presented here show this is not necessary to obtain meaningful results, at least for v⋆ < 100 km s−1.
For 2D and 3D calculations, the supersonic motion of gas across the grid generates anisotropic numerical diffusivity in grid-aligned flows, and this can produce spurious numerical instability (Quirk 1994; Sutherland, Bicknell & Dopita 2003). If the flow is grid-aligned, the I-front instability is much stronger along the grid-axis and this eventually destroys the solution. A number of corrections were attempted but none completely removes this effect, so the best solution was to rotate the bulk flow by 40° so that the numerical diffusivity in the x and y directions are comparable. The flow velocity is rotated in 3D simulations for the same reason, and is set at 50° to the z-axis and 40° to the x-axis.
ONE-DIMENSIONAL SIMULATIONS
The primary parameter determining the properties of an I-front is its Mach number relative to the isothermal sound speed in photoionized gas |$\mathcal {M}\equiv v/c_i$|, and this depends only weakly on the stellar radiation field, ISM number density and ISM metallicity. For this reason, only the star's space velocity is varied here; the other properties of the star and the ISM are kept constant. Stellar space velocities from 10 to 200 km s−1 have been simulated, specifically v⋆ ∈ {10, 20, 25, 30, 40, 50, 75, 100, 150, 200} km s−1. The models have been run for spatial resolutions r1 to r7 and for different timestep criteria.
1D simulations allow us to model both the upstream I-front and the downstream recombination front, although multidimensional effects (not to mention the shocked stellar wind) are expected to alter the recombination front. As well as studying the physical properties of the fronts, this also allows a very clean test of the numerical error of our solution as a function of spatial and temporal resolution. The simulations relax to a stationary state in all cases where v⋆ is greater than Mach 2 in ionized gas. The effects of spatial and temporal resolution on the 1D simulations are discussed in Appendix A; here we note that it is not necessary to spatially resolve the I-front (cell optical depths could be τ > 1 without loss of accuracy) at least for v⋆ < 100 km s−1. It is necessary, however, to use a smaller-than-usual Courant–Friedrichs–Lewy (CFL) number of Ccfl = 0.1 for v⋆ < 100 km s−1, and even smaller for larger v⋆. When the I-front is spatially resolved the timestep criterion is not so important.
Profiles of the ISM density, H+ fraction (1 − yn), and temperature are shown in Fig. 2 for all models with v⋆ ≥ 25 km s−1, with a profile for a static star in a static medium (i.e. the original case considered by Strömgren 1939) overplotted for reference. We have not plotted lower velocities where the I-front is D-type because the H ii region expands further than the Strömgren radius Rs, and the dense shell that forms remains unstable for many grid advection times (the time for the background flow to cross the grid).
At larger velocities the I-front retreats with increasing v⋆. To quantify this, the position of the I-front (defined as the furthest upstream grid zone with yn < 0.5) is plotted as a function of v⋆ in Fig. 3 for simulations with multifrequency radiation (the default case) and also for models with monochromatic radiation and the same ionizing photon luminosity (by number of photons). The model with monochromatic radiation has a hotter interior part of the H ii region, hence a smaller recombination rate and a larger H ii region.
The normalization of the curve depends on tr, which in turn depends on the temperature-dependent recombination rate αb. The monochromatic radiation model has a roughly constant temperature in the H ii region, whereas the multifrequency radiation model has a cooler interior and a hotter border because of spectral hardening. This implies a higher recombination rate and therefore a smaller H ii region. The analytic curves were normalized to a H ii region temperature of T = 7600 K for monochromatic radiation and 6300 K for multifrequency radiation.
These 1D results and their comparison with analytic theory provide the basic understanding of the global asymmetry of moving H ii regions, which will be confirmed and further explored by multidimensional models in the following sections. They also allow verification and testing of the algorithms in terms of their convergence properties and timestepping requirements (see Appendix A).
Comparison to previous work
An R-type I-front with v⋆ = 50 km s−1 was modelled by Henney et al. (2005) as part of a study of steady I-fronts that included advection self-consistently. They found that the supersonic I-front has a strong temperature peak at the I-front that is absent in the static case, and that the location of the I-front is closer to the radiation source by about 5 per cent. The temperature peak in our case is perhaps a bit stronger and increases in amplitude with v⋆, but is qualitatively similar. The position of the I-front in our calculations is modified by much more than 5 per cent compared to the static case, although Henney et al. (2005) use plane–parallel radiation so this cannot be compared directly to our results.
Rasiwala (1969) studied H ii regions around moving stars assuming isothermal gas and neglecting recombinations, so agreement with our results is not expected. The best agreement is in the position and properties of the upstream I-front, but the downstream recombination front is of course not captured at all by the Rasiwala (1969) calculation. Subsequent calculations including recombination (but without hydrodynamics) were made by Thuan (1975) for the H ii region around a v⋆ = 50 km s−1 O star in a diffuse ISM with hydrogen number density |$n_{{\small h}}=0.1\,\mathrm{cm}^{-3}$|. The upstream and downstream I-front radii were predicted to be Rup = 140 pc and Rdn = 265 pc, respectively, compared to Rs = 190 pc. In our case for v⋆ = 50 km s−1, we find Rup/Rs = 0.85, significantly larger than the value of 0.74 from Thuan (1975); this is simply explained by the different values for tr and Rs in the two calculations. Raga et al. (1997a) calculated Rup/Rs = 0.93 for an O5 star with v⋆ = 100 km s−1 and |$n_{{\small h}}=1\,\mathrm{cm}^{-3}$|; this again differs from the values obtained here because of the differing tr and Rs (we consider an O9.5V star with a much smaller ionizing photon luminosity). The ratio Rup/Rs decreases as the pre-factor in the quadratic term of equation (9) increases; as long as this is taken into account our results agree well with these previous calculations. The 1D profiles in Chiang & Rappaport (1996) show much broader I-fronts as a consequence of the much harder radiation spectrum they consider, and so they are not (and are not expected to be) comparable to our results.
3D SIMULATIONS
The 3D simulations run are listed in Table 1. They consist of a 3D domain with {x, y, z} ∈ [ − 32.64, 18.56] pc with a Cartesian coordinate system, and with the star at the origin. There are simulations where the ISM magnetic field is aligned with (BA3r1, BA3r2) and perpendicular to (BT3r1, BT3r2) the direction of motion and simulations with no magnetic field (HD3r1, HD3r2). All MHD models have a field strength B = 7 μG, corresponding to a plasma parameter β ≡ 8πpg/B2 = 0.2 in the neutral ISM, and β = 2.7 in fully ionized gas at T = 7000 K. The neutral medium is therefore magnetically dominated, whereas the ionized gas is gas pressure dominated. Models with a perpendicular field have the field vector in the x–y plane. The simulations have inflow boundary conditions on the upstream boundaries, and only-outflow conditions downstream, and it is assumed that the star ‘switches on’ instantly at t = 0. As a result, the initial evolution is not very meaningful because no stars are born in the WNM. The time for the ISM to advect across the diameter of the H ii region (≈20 pc) is about 0.75 Myr, so we expect the effects of initial conditions to disappear soon after this time. This is also about the time that gas affected by the star's radiation first leaves the downstream boundary. A time-stationary state takes ≈1.5 Myr to be established.
ID . | (Nx, Ny, Nz) . | Δx . | |$\boldsymbol {B}$|-field (μG) . |
---|---|---|---|
HD3r1 | (160, 160, 160) | 0.32 | (0, 0, 0) |
HD3r2 | (320, 320, 320) | 0.16 | (0, 0, 0) |
BA3r1 | (160, 160, 160) | 0.32 | (7, 0, 0) |
BA3r2 | (320, 320, 320) | 0.16 | (7, 0, 0) |
BT3r1 | (160, 160, 160) | 0.32 | (0, 7, 0) |
BT3r2 | (320, 320, 320) | 0.16 | (0, 7, 0) |
ID . | (Nx, Ny, Nz) . | Δx . | |$\boldsymbol {B}$|-field (μG) . |
---|---|---|---|
HD3r1 | (160, 160, 160) | 0.32 | (0, 0, 0) |
HD3r2 | (320, 320, 320) | 0.16 | (0, 0, 0) |
BA3r1 | (160, 160, 160) | 0.32 | (7, 0, 0) |
BA3r2 | (320, 320, 320) | 0.16 | (7, 0, 0) |
BT3r1 | (160, 160, 160) | 0.32 | (0, 7, 0) |
BT3r2 | (320, 320, 320) | 0.16 | (0, 7, 0) |
ID . | (Nx, Ny, Nz) . | Δx . | |$\boldsymbol {B}$|-field (μG) . |
---|---|---|---|
HD3r1 | (160, 160, 160) | 0.32 | (0, 0, 0) |
HD3r2 | (320, 320, 320) | 0.16 | (0, 0, 0) |
BA3r1 | (160, 160, 160) | 0.32 | (7, 0, 0) |
BA3r2 | (320, 320, 320) | 0.16 | (7, 0, 0) |
BT3r1 | (160, 160, 160) | 0.32 | (0, 7, 0) |
BT3r2 | (320, 320, 320) | 0.16 | (0, 7, 0) |
ID . | (Nx, Ny, Nz) . | Δx . | |$\boldsymbol {B}$|-field (μG) . |
---|---|---|---|
HD3r1 | (160, 160, 160) | 0.32 | (0, 0, 0) |
HD3r2 | (320, 320, 320) | 0.16 | (0, 0, 0) |
BA3r1 | (160, 160, 160) | 0.32 | (7, 0, 0) |
BA3r2 | (320, 320, 320) | 0.16 | (7, 0, 0) |
BT3r1 | (160, 160, 160) | 0.32 | (0, 7, 0) |
BT3r2 | (320, 320, 320) | 0.16 | (0, 7, 0) |
Observable properties of simulations
The images show some of the same features of previous axisymmetric models of hypersonic stars (Raga et al. 1997a), notably the almost circular H ii region and the tail of overdense gas expanding from its lateral edges. Here, because the Mach number of the flow is much lower, there is noticeable expansion of this overdense shell in a cone shape, and it forms at a smaller angle from the velocity vector of the star. The |$N_{\mathrm{H\,{\small {I}}}}$| contours show the overdense shell, and also the underdense H ii region and the wake behind it. Underdensity in the H ii region is because the gas is ionized (and so does not show up in H i) but in the wake where the gas has recombined, the gas is underdense because the H ii region has displaced an outward-flowing expanding shell.
The left-hand plot with zone-size Δx = 0.32 pc (r1) has relaxed to a steady state with a stable I-front, but this is only because of the low numerical resolution. The right-hand plot with Δx = 0.16 pc (r2) shows instability in the I-front, resulting in substructure in both Hα emission and |$N_{\mathrm{H\,{\small {I}}}}$|. Still higher resolution 2D simulations in Appendix B show that the instability is only marginally resolved with resolution r2. The unstable I-front generates knots of dense gas on small scales; for r2 there are a few of them present in the I-front at any time, and they move slowly away from the symmetry axis around to the sides of the H ii region where they join the expanding shell. This instability was not found by Raga et al. (1997a), probably because our model has a much lower Mach number and the H ii region is in a higher density environment where shocks are more dissipative. The evolution of these knots is similar to that depicted in García-Segura & Franco (1996) for the expansion of D-type I-fronts. In our case, the I-front is moving at a constant velocity that is (coincidentally) close to the limit between R-type and D-type (Mach 2 in ionized gas, Mihalas & Mihalas 1984). In regions where the I-front is expanding there is no upstream shock, but in overdense regions where it is receding, a radiative shock forms that enhances the density and mass of the overdense clumps. The underdense regions that push neutral gas towards the overdense regions initially resemble the spear-shaped features in previous work (García-Segura & Franco 1996; Franco et al. 2007; Whalen & Norman 2008), but they soon saturate and develop a bubble shape. The lifetime of the knots is much longer than the length of time a parcel of gas spends in the knot, so while they are formed from dynamical instability they seem to be quite stable structures.
The dense knots are seeds of strong photoevaporation flows that expand conically into the H ii region when they are near the symmetry axis. The photoevaporation flow would expand spherically for a static ionizing source (cf. the Eagle Nebula pillars; Hester et al. 1996) but here the supersonic advection drags the photoevaporation flow downstream, resulting in a cone. As the dense knots move off-axis to larger angles on the I-front, the photoevaporation flow is no longer symmetric because one side of the flow tries to expand upstream and fails, leading to the Hα arcs seen in Fig. 4 that trail downstream from the bottom to the centre of the H ii region. Features like this should be observable in H ii regions around moving stars if this I-front instability is present, although ISM inhomogeneities may drive similar photoevaporation flows.
The magnitude of the overdensity and underdensity in neutral H (compared to the mean) is |$N_{\mathrm{H\,{\small {I}}}}\approx 3\!-\!4\times 10^{20}\,\mathrm{cm}^{-2}$| and |$N_{\mathrm{H\,{\small {I}}}}\approx -2\times 10^{20}\,\mathrm{cm}^{-2}$|. The undisturbed background gas density provides |$N_{\mathrm{H\,{\small {I}}}}\approx 1.8\times 10^{20}\,\mathrm{cm}^{-2}$| through the 24 pc diameter of the H ii region, so the contrast is only a factor of 2-3. This may be difficult to measure observationally, especially in the Galactic plane where there is strong H i background emission. It should, however, have similar radial velocity to the mean Hα emission, so this could help to identify the conical shell.
The upstream I-front is much sharper than the downstream recombination front because of their different character; the recombination length ℓr ≈ v⋆tr (recall tr = (αbne)−1 is the recombination time) sets the recombination front thickness; the I-front is very thin in comparison. This is an obvious feature of the Hα emissivity plots, and can also be seen in the Hα map of ζ Oph's H ii region (Gvaramadze et al. 2012b, see also Section 5). In addition the Hα emission in our simulations has a clear gradient, being brighter upstream from the star than downstream. This is because the H ii region expands as gas flows downstream, so the recombination rate decreases and the emissivity along with it. This is just as clear in the HD3r2 plot (see right-hand panel of Fig. 4) even though the gas is more disturbed.
Equivalent results are shown in Fig. 5 for the MHD simulations with a magnetic field aligned with the stellar velocity, BA3r1 and BA3r2. In this case, the higher resolution model BA3r2 also has a stable I-front and there is no significant difference between the two models. It has been found for H ii regions around static stars (Krumholz, Stone & Gardiner 2007; Arthur et al. 2011; Mackey & Lim 2011) that D-type I-front expansion along field lines is largely stable, whereas expansion perpendicular to field lines is unstable, producing ‘beads’ of dense gas tied to field lines (see also the 2D studies of Williams 2002, 2007). Similar results are seen here for R-type I-fronts, although we note that 2D simulations at higher resolution do become unstable even when velocity and magnetic field vectors are aligned (see Appendix B). This agrees with the finding (Newman & Axford 1967; Williams 1999) that R-type I-fronts should be generically unstable to the shadowing instability for finite density perturbations, suggesting that the stability of BA3r1 and BA3r2 is at least partly a resolution effect.
Another feature of simulations BA3r1 and BA3r2 is that the expanding neutral shell is much less dense than in the hydrodynamic models, because magnetic pressure resists compression and expansion perpendicular to field lines. The fast magnetosonic speed is larger than the sound speed so the effective Mach number of the flow is lower, hence the outward-moving shock makes a larger angle with the direction of gas flow. The fast magnetosonic shock is also less compressive than its hydrodynamical equivalent. Both of these effects mean that the shell is more diffuse and its column density is lower.
Fig. 6 shows projections through model BT3r2, this time at t = 4.5 Myr and from two different projection angles. Both projections are perpendicular to the bulk gas flow as before, but the left-hand plot is projected such that the background magnetic field is along the LOS, whereas in the right-hand plot it is vertical and in the image plane. The left-hand plot is very similar to the results in Fig. 5 for BA3r2, in that the shell column density is low, the angle it makes with the flow velocity vector is quite large, and the upstream I-front is more stable than the hydrodynamic model HD3r2. This relates directly to the magnetic field orientation, which in both cases inhibits gas flow in the image plane. The right-hand plot is much more similar to the hydrodynamic model HD3r2, in that there is instability in the I-front and a strongly compressed shell downstream.
The kink in the upstream Hα emission in the left-hand plot seems to cross most of the H ii region, and its counterpart in the right-hand image shows the kink end-on as a v-shape notch in the upstream I-front. It appears that the 3D I-front instability seen in model HD3r2 has reverted to an almost 2D instability in BT3r2, because of the large-scale ordered magnetic field. It is not clear why the ripple in the I-front would form in this way since its connectivity is across the field lines, not along them. Gas can only flow freely along field lines, so these unstable clumps must preferentially form in a 2D fashion, but it is not clear why they should be connected across the full length of the I-front. Higher resolution simulations are required to investigate if this feature is numerical or physical.
Shape of H ii region
Plots of the H ii region radius are shown for the different 3D simulations in Fig. 7. The range of values for the upstream and downstream I-front radius are plotted with error bars, the maximum and minimum I-front radius over its full surface is shown with the blue points, and the mean I-front radius is plotted as the blue crosses. The upstream I-front radius is about 9.5 pc in all models (compared to 9.83 pc in the 1D simulations), whereas the Strömgren radius is Rs ≈ 10.6 pc.1 The downstream radius is not well defined because the recombination length is significant, so yn = 0.5 is chosen to define the radius. In this case it is about 13.5-14 pc, so the diameter is 23-23.5 pc, or 8.5-11 per cent larger than 2Rs. The mean radius calculated from the spherically averaged ionized volume is about 12 pc in all models (1.13Rs). This larger mean radius is expected because the mean density within the H ii region is less than the ISM density. The maximum radius is always in the range of directions 120° < θ < 150° (where θ measures angles from the star's velocity vector, with the star at the origin), implying that the H ii region is flattened in the downstream direction.
Simulations BA3r1 and BA3r2 are basically the same because both are stable and in a steady state at late times. BT3r2 and HD3r2 have I-front instability, so here the upstream I-front has a range of radii: about 8−10 pc for HD3r2 and 8.8−9.8 pc for BT3r2. The maximum and minimum radius also fluctuate because of this instability.
Despite the clear asymmetry in the H ii region radius relative to the star, it is much more spherical in Hα emission with respect to its own geometric centre. It is also noteworthy that the star is located downstream of the peak Hα intensity, but upstream of the H ii region centre.
H ii region feedback on the ISM
Figs 4–6 show that the H ii region has a measurable impact on the ISM that persists once the star has passed. The quantity of compressed gas as a function of time and the density distribution of gas after 4 Myr are plotted in Fig. 8 for each of the 3D simulations. Note that only gas still on the simulation domain is counted, so this measures the mass of dense gas generated in the time from when the gas is compressed to when it leaves the simulation domain (≈0.8−1 Myr). The density distribution is shown in terms of gas mass with a density greater than |$n_{{\small h}}$|. After about 1.5 Myr all models approach a stationary state.
Simulations BA3r1 and BA3r2 have the weakest density enhancement with almost no gas compressed to ρ > 2ρ0. This is because all of the pressure gradient is across field lines, so magnetic pressure and tension resist gas compression. BT3r1 and BT3r2 have more dense gas because field lines only resist compression in one of the two directions perpendicular to motion, and HD3r1 and HD3r2 have the densest gas because there is no magnetic field to resist gas flows. In HD3r2, there is 100 M⊙ of gas that has been compressed by about a factor of 10, and 2D results show that at higher resolution even stronger compression is obtained (Appendix B). The lower panel of Fig. 8 demonstrates that the simulations have converged for weak compression, but that the mass of strongly compressed gas is very resolution-dependent, as would be expected. These results suggest that the passage of massive stars through the ISM could compress WNM gas sufficiently for it to make the transition to cold neutral medium. They also demonstrate the importance of the ISM magnetic field in determining the compression ratio in these weak radiative shocks. Fig. 8(b) also shows that more than half of the gas in the simulation is underdense compared to the undisturbed upstream gas; the star creates an underdense cylindrical volume in the ISM with cross-section of about 10 pc radius.
Despite these results, it is not obvious that the H ii region feedback effects are significant compared to those expected from stellar winds. The main stellar feedback mechanism in terms of energetics is photoheating, but this energy is quickly dissipated once the star has passed and the gas has recombined. For static stars, the conversion efficiency of ionizing photon luminosity to kinetic energy was found to be ≲0.1 per cent (Freyer, Hensler & Yorke 2003) (although they included the ionization energy of H as well as the excess heating energy in this estimate), whereas for stellar winds ≈10 per cent of the input kinetic energy survives dissipative processes and drives motion in the ISM (Weaver et al. 1977; García-Segura, Mac Low & Langer 1996; Krause et al. 2013).
Kinetic energy feedback
The only persistent feedback from the H ii region is the part of the photoheating energy that goes into driving kinetic energy in the surroundings, from H ii region expansion and I-front instability. To measure this, we calculate the total gas momentum and kinetic energy on the simulation domain as a function of time, with the bulk flow subtracted off, plotted in Fig. 9. The slope of this function gives the feedback rate at early times (before disturbed gas leaves the downstream boundary), and the stationary state value at late times divided by the grid-crossing time gives another approximate feedback rate. The kinetic energy increases somewhat with spatial resolution because of waves/shocks generated by I-front instability, which is more vigorous at higher resolution. A stronger trend is that the models with aligned magnetic field have lower kinetic energy than those with a transverse field, which in turn lie below the zero field (hydrodynamic) models. This is because the tension of the aligned magnetic field resists H ii region dynamical expansion in all directions perpendicular to the direction of motion. The perpendicular field models allow expansion along field lines (i.e. in one of the two perpendicular directions), and hydrodynamic models have no magnetic tension to inhibit expansion.
As a comparison, the kinetic energy plot also shows the mechanical luminosity of the stellar wind driven by ζ Oph, |$L_w = 0.5\dot{M}v_w^2$|, using a mass-loss rate of |$\dot{M} \approx 2.2\times 10^{-8} \,\mathrm{M}_{{\odot }}\,\mathrm{yr}^{-1}$| and wind velocity of vw ≈ 1500 km s−1 (Gvaramadze et al. 2012b). With the stellar luminosity and spectrum used here, the mean energy per ionizing photon is 17.2 eV, so the mean heating per ionization is 3.6 eV. If we consider the kinetic energy of the ISM as a percentage of this heating input power, we obtain an efficiency of conversion from photoheating to kinetic energy of ϵ ≈ 0.5–1 per cent. Even with this low efficiency the kinetic energy gained from the H ii region is comparable to that from the stellar wind and in reality is probably larger because much of the wind kinetic energy gets lost through dissipative processes.
Momentum feedback
Perhaps a truer measure of the feedback is the rate of momentum input to the ISM, because momentum is conserved regardless of the dissipative properties of the gas. Here the resolution dependence is not apparent and the differences between the models are also less pronounced, but a similar dependence on magnetic field orientation is seen. For comparison, we also show the momentum input we would obtain from the wind of ζ Oph, with parameters given above. Its momentum input rate is multiplied by 100 to show it on the same scale, so evidently it is insignificant compared to the H ii region's dynamical expansion. The rate of momentum input from H ii region expansion (measured as the slope of the line for t < 0.8 Myr) is also about 10 times larger than the total momentum output rate of the radiation from ζ Oph, (L/c)t, plotted as the green dot–dashed line in Fig. 9 for a luminosity of L = 6.4 × 104 L⊙ (Howarth & Smith 2001, scaled to a distance of 112 pc). Each photon would therefore need to have ≥10 scatterings before escaping the H ii region in order to significantly affect the dynamics presented here. This is unlikely given the low ISM density.
Comparison to supernova feedback
By both measures (kinetic energy and momentum) the H ii region expansion has stronger feedback than the stellar wind for an O9.5 V star, even though it is moving supersonically through the ISM. This feedback is not very violent, and may not be as easily observed as a stellar wind bow shock, but it is more important in terms of energetics. Compared to the kinetic energy of the eventual supernova explosion of ζ Oph, Esn ≈ 1051 erg, both feedback effects are not very important. If we take the main-sequence lifetime of ζ Oph to be τms ≈ 7 Myr, then the total kinetic energy imparted to the ISM from the H ii region is ≈3.8 × 10−3Esn. The momentum input to the ISM is, however, similar to that of a supernova, at 8.8 × 104 M⊙ km s−1. If we take the supernova to have 4 M⊙ of ejecta with velocity 5000 km s−1 (giving a kinetic energy of 1051 erg), then we see that the H ii region imparts about 4× as much momentum to the ISM during the star's main-sequence lifetime.
The calculations here have only considered a single stellar mass, ISM density and stellar space velocity. It is not trivial to estimate how the relative importance of the different feedback processes scale with all of these parameters, particularly stellar mass, but it should be explored in future work.
ISM near the star
Models of bow shocks around hot stars usually assume (for simplicity) that the H ii region does not affect the ISM density or dynamics significantly (e.g. Comeron & Kaper 1998). This is clearly true in the hypersonic limit, and false in the subsonic limit, but the effects for intermediate v⋆ have not so far been explored. We analysed the mean gas properties near the star for all simulation snapshots to investigate their spatial and temporal variations. Two averaging volumes have been used: r < 1.0 pc and r < 0.32 pc, corresponding to typical bow shock sizes. Within these volumes the mean density |$\langle \rho \rangle = (1/N)\sum _{i=1}^{N} \rho _i$|, root-mean-squared (rms) density variance |$\sigma = \sqrt{(\langle \rho ^2\rangle -\langle \rho \rangle ^2)}/\langle \rho \rangle$|, and mean velocity (in the star's rest frame) |$\langle v\rangle = (1/N)\sum _{i=1}^{N} |v_i|$| of the ISM were calculated. Here N is the number of grid zones in the averaging volume.
Plots of these three quantities as a function of time are shown in Fig. 10. At early times (t < 1 Myr), the H ii region undergoes large-scale changes in shape leading to a decrease in density and subsequent increase. This is a strongly damped oscillation driven by the initial out-of-equilibrium state. The low-resolution simulations, HD3r1, BA3r1 and BT3r1, show no later time variation because they have reached a steady state. The higher resolution aligned magnetic field model BA3r2 also reaches a steady state, with values very similar to BA3r1. Models BT3r2 and HD3r2, however, show persistent fluctuations in gas properties near the star, because of the I-front instability and the waves it creates in the H ii region. The mean gas density 〈ρ〉 is ≈10 per cent below the ISM mean density because the H ii region has expanded slightly. The mean value is independent of the averaging volume for low-resolution models, indicating that a linear density gradient through the volume is a good approximation. At higher resolution, the temporal fluctuations in 〈ρ〉 are stronger on smaller scales, as expected if there are isolated planar shocks passing through.
The rms density variance is larger for a larger averaging volume because there are more cells and so more independent volume elements. For the small averaging volume, the density values are strongly correlated because of numerical diffusion. The rms fluctuations in density are about 1-5 per cent for this spatial resolution, although this is expected to rise with higher resolution simulations that have stronger I-front instability (see Appendix B).
Interestingly, the velocity flowing past the star is 7.5–10 per cent smaller than the ISM bulk velocity. This is caused by the I-front jump conditions for an R-type I-front, which predict a small density increase and velocity decrease according to equation (6). The H ii region temperature is T1 ≈ 6300 K, so the isothermal sound speed is about c1 ≈ 9 km s−1 giving |$\mathcal {M}_1\approx 3$|, and we expect ρ1/ρ0 ≈ 1.1 and v1/v0 ≈ 0.9. The velocity decrease is seen in Fig. 10 but the density increase is not. Inspection of slices through the simulations shows that at the upstream I-front this density increase is found, but that the H ii region expansion has reduced the density at the position of the star. This multidimensional effect was of course not seen in the 1D simulations (see Fig. 2).
These changes to the ISM properties generated by the H ii region dynamics are relatively small perturbations to the bulk flow properties, and are not expected to disrupt or greatly perturb the stellar wind bow shock that should be present. It remains to be seen, however, what effect an inhomogeneous ISM will have on this picture. In the cloud-zapping regime (Bertoldi 1989) where clumps are photoionized by R-type I-fronts, the H ii region should homogenize the ISM (e.g. McKee, van Buren & Lazareff 1984), but if clumps are large enough to trap the I-front then a photoevaporation flow can be set up with velocity up to v ≲ 30 km s−1, which could strongly affect the bow shock.
DISCUSSION
Comparison to the H ii region of ζ Oph
The nearest example of an H ii region around an exiled O star is that of ζ Oph, Sh 2-27. In Fig. 11, we show the Hα image of this H ii region originating from the Southern Hα Sky Survey Atlas (SHASSA; Gaustad et al. 2001).2 A detailed comparison to our simulations is beyond the scope of this paper, and will be pursued in future work; here we note some of the morphological similarities to (and differences from) the simulated H ii regions. The most striking difference is the clumpy nature of the Hα emission in the observational image, most likely related to underlying ISM inhomogeneity and/or foreground patchy extinction. The ridge of strong absorption towards the bottom of the H ii region is correlated with a molecular cloud (Complex 4 in de Geus, Bronfman & Thaddeus 1990), indicating significant extinction due to foreground clouds that may or may not be physically associated with the H ii region. Complex 1(E) from de Geus et al. (1990) also provides some extinction around (l, b) = (5°, 23°), downstream from and just below the star. Even allowing for patchy extinction, it still seems clear that underlying ISM inhomogeneity is a significant component of the ζ Oph H ii region (cf. Wood et al. 2005).
Nevertheless, if we consider only the upper half of the H ii region where there is less foreground extinction, there are definite similarities to our simulated Hα maps. The arc of the upstream I-front is sharp, whereas the downstream recombination front has a much more gradual decrease in intensity with distance from the star. In addition, the downstream quadrant is generally of lower intensity than the upstream quadrant. The distance to the upstream I-front is smaller than that to the downstream recombination front (Rup/Rdn ≈ 3/4), and the H ii region has its largest radius at about 120-150° from the star's velocity vector. All of these qualitative features are also found in the simulation results.
The LOS magnetic field through Sh 2-27 has been measured using Faraday rotation by Harvey-Smith, Madsen & Gaensler (2011), and found to be reasonably strong at 6.1 μG (with some uncertainty arising from assumptions about gas clumping and distance). If the field is this strong then the compression factor of the neutral gas shell found in our simulations would be significantly reduced by the magnetic pressure, and an H i column density map of Sh 2-27 would resemble the left-hand panel of Fig. 6 more than the right-hand panel. If the plane-of-sky magnetic field is much weaker than 6.1 μG, then strongest shell compression is along the LOS, in which case the shell may be observable in position–velocity H i data.
In Gvaramadze et al. (2012b), where the mass-loss rate of ζ Oph was estimated based on the sizes of its H ii region and bow shock, it was assumed that the mean H ii region gas density is the same as the mean density at the bow shock, and furthermore that the ISM upstream from the bow shock remains at rest. We have shown in Fig. 10 that these assumptions are reasonably well justified, in that the ISM density near the star (ρi) is only about 10 per cent lower than that of the undisturbed ISM (ρ0), and the ISM has been accelerated only slightly in passing through the R-type I-front (to vi = 2−2.4 km s−1 in the upstream direction). These are both small modifications of the upstream ram pressure in the star's rest frame, reducing |$\rho _i v_i^2$| by about 25 per cent. The bow shock standoff distance, |$R_{\mathrm{{\small so}}}\equiv \sqrt{\dot{M}v_w/4\pi \rho _iv_i^2}$|, is then increased by about 12 per cent. The analytic model in Gvaramadze et al. (2012b) can therefore be applied with only minor corrections even for the relatively slowly moving exile ζ Oph.
ISM inhomogeneity
We have not addressed ISM clumping in this work, assuming as a first step that the ISM is smooth. A large clumping factor could affect the Gvaramadze et al. (2012b) analysis because Rs is determined by recombination which scales with ρ2, whereas |$R_{\mathrm{{\small so}}}$| is determined by ram pressure which is linear in ρ. A very clumpy medium would therefore have a smaller ratio |$R_{\rm s}/R_{\mathrm{{\small so}}}$| than the maximal value obtained for a smooth medium. Inspection of the observed Hα maps in Fig. 11 strongly suggests some degree of underlying inhomogeneity in the ISM. This remains a caveat to the Gvaramadze et al. (2012b) analysis, and will be studied further in future work.
This question has also been approached from another angle by Wood et al. (2005), who used Hα, [N ii] and [S ii] emission maps from the H ii region around ζ Oph to constrain the porosity (or equivalently clumpiness) of the ISM through which it is passing. They used a radiative transfer code to predict the emission properties of a static density field consisting of clumps of various sizes and filling factors embedded in a low-density medium. They compared the circularly averaged predictions to observations, for Hα brightness and the [N ii]/Hα and [S ii]/Hα line ratios, as a function of distance from the star. It was found that some clumping of the ISM was required to explain the spatial distribution of the H ii region surface brightness in these lines, given the modelling assumptions.
Figs 4–6 show, however, that circular symmetry is a bad approximation for Hα emission because the upstream I-front has very different properties to the downstream recombination front. This should also hold for other emission lines (cf. Henney et al. 2005). It would be very interesting to repeat the Wood et al. (2005) analysis, but with a 2D spatially resolved comparison using a method such as that proposed by Wood et al. (2013). Even with this refinement, however, it is not clear that the supersonic relative motion between the star and the ISM can be neglected in the radiative transfer modelling because ionization equilibrium is never attained (although a stationary state is). Supersonic advection compresses the upstream I-front significantly (Raga et al. 1997a; Henney et al. 2005) and significantly changes the temperature and ionization structure of the downstream H ii region border from an I-front to a recombination front (Newman & Axford 1968; Williams & Dyson 1996). No static model will be able to capture the downstream recombination front correctly, or the spatial compression of the upstream I-front. The degree to which this is important should be investigated in future work where we will include more ions in the microphysics network, allowing us to predict line ratios directly.
Isolated O and B stars in the Small Magellanic Cloud
Recently Oey et al. (2013) presented observations of 14 single-star H ii regions in the Small Magellanic Cloud (SMC), for each of which the ionizing source is a late-O or early-B star located at least 28 pc from its nearest O- or B-type neighbour. These were interpreted as candidate static stars that formed in isolation, largely on the basis of a lack of stellar wind bow shock, the circular shape of the H ii region, and the approximately central location of the star within the H ii region.
Regarding the shape of the H ii region, we have demonstrated that, at least for v⋆ = 26.5 km s−1, the H ii region remains approximately circular in projected Hα emission. The main difference for an exile is that the upstream edge is sharper and has higher surface brightness than the downstream edge (see Fig. 12 which shows an Hα emission map of simulation HD3r2 overlaid with linearly spaced contours from 10 to 90 per cent of the peak emission). A number of the H ii regions in the Oey et al. (2013) sample have similar asymmetries in their Hα emission, although we have no way to constrain whether the asymmetries are correlated with their unknown proper motion. Similar results to ours were also obtained by Raga et al. (1997a) for a much higher velocity exile (v⋆ = 100 km s−1) in a somewhat lower density medium (n = 1 cm−3). Both of these results show that the shape of the H ii region in Hα emission is not a strong indicator of stellar motion unless we are discussing properties at the 10-25 per cent level (but see Chiang & Rappaport 1996, for results from much harder spectrum X-ray sources in a denser medium). Furthermore, there is a well-known degeneracy between H ii region asymmetry produced by ISM density gradients and by stellar motion (see Arthur & Hoare 2006, for a detailed investigation). When the likely presence of ISM inhomogeneity is added to this, it may be difficult to draw strong conclusions regarding stellar motion from H ii region shapes in Hα emission maps. A better probe of stellar motion is probably the neutral shell that forms around the H ii region. For v⋆ ≳ 20–30 km s−1 there will be no upstream or downstream dense shell, whereas a static star should have a shell of swept-up ISM in nearly all directions.
The offset of the star from the centre of the H ii region depends on how the centre is defined. Identifying the H ii region border by eye, many of the stars in the Oey et al. (2013) sample are offset from the H ii region centre by a similar degree to the star in our simulations. This is a weak statement, however, and should be made more quantitative. From Fig. 12, the star is significantly offset downstream from the 90 per cent intensity contour, and significantly offset upstream from the 10 per cent contour, but is very close to the centre of the (almost circular) 50 per cent contour.
The Oey et al. (2013) H ii region sample is clearly an important data set for probing the formation process and environment of isolated massive stars, but we argue that further quantitative predictions for the properties of H ii regions around moving stars (with a range of spatial velocities) are required before conclusions can be drawn about whether these stars are exiles or formed in situ.
Implications of our results
We have seen that the dynamical response of the ISM to a passing exiled star is that an overdense expanding conical shell is swept up, leaving an underdense wake behind it. Stellar winds are not included in our simulations because the range of scales is too large, but hot shocked-wind material will also gather in the underdense wake behind the star. It is therefore plausible that exiled massive stars create underdense ‘tunnels’ through the ISM which can then be maintained by energy from supernova explosions (cf. Cox & Smith 1974), since the tunnels always originate at star clusters. The momentum feedback of the H ii region expansion is at least as strong as that of a supernova for the parameters considered here, and will be even stronger for earlier O stars. This suggests that exiled massive stars can make a significant contribution to the dynamics of gas in galaxy discs. Furthermore, work quantifying the dependence of this feedback mechanism on stellar mass, ISM density and metallicity are warranted to deduce its importance at a global level.
An important point, first made by Raga et al. (1997a) and emphasized here, is that the H ii regions around exiled stars remain roughly circular, and the star stays near the centre of the H ii region unless it is moving with a really extreme velocity (also for H ii regions around X-ray sources; Chiang & Rappaport 1996). The H ii region produced by a moving star remains quite spherical because it is always dynamically young, in contrast to the often-complex shapes of old H ii regions around static stars. This is only true for H ii regions in which the stellar wind bow shock does not absorb a significant fraction of the star's ionizing photon output. In a dense medium, the H ii region can be trapped by the bow shock (Mac Low et al. 1991; Arthur & Hoare 2006) for two reasons: the bow shock standoff distance decreases with increasing density less rapidly (|$R_{\mathrm{{\small so}}}\propto \rho ^{-1/2}$|) than the H ii region radius (Rs∝ρ− 2/3), and gas compression in the bow shock increases with density because of the shorter cooling time.
Conroy & Kratter (2012) suggested that runaway massive stars could have contributed significantly to the reionization of the Universe because they move out of their small protogalaxies within their main-sequence lifetime, and so all of their ionizing radiation escapes their parent galaxy. Dynamical feedback effects from H ii regions should become much stronger at low metallicity because the H ii region is hotter and less thermal energy is dissipated through radiative cooling. We suggest that, in addition to the process considered by Conroy & Kratter (2012), the exiled stars also created lower density tunnels from the centres of protogalaxies to their outskirts, increasing the escape fraction of ionizing photons and potentially providing channels through which supernova-enriched gas can more easily escape to the intergalactic medium (IGM). Of course the exiled stars themselves also explode in the IGM, leading to in situ enrichment and the generation of intergalactic shockwaves. Gvaramadze, Gualandris & Portegies Zwart (2009) suggested that high-velocity (v⋆ > 200 km s−1) O and B stars ejected from galactic discs could be the ionizing sources of extraplanar H ii regions observed ≈0.5–1.5 kpc above some galactic discs (e.g. Tüllmann et al. 2003). Such stars could travel outside the virial radius of a 108 M⊙ galaxy at redshift z ∼ 10 (see Conroy & Kratter 2012). The small size and lower escape velocities of high-redshift galaxies make exiled stars much more important at high redshift than in the local universe.
CONCLUSIONS
We have studied the dynamics of H ii regions around supersonically moving hot stars with multidimensional MHD simulations including non-equilibrium photoionization. While the response of the ISM to the H ii region is not violent, it is significant and leaves a long-lasting imprint on the ISM once the star has passed. Only one stellar velocity has been considered, v⋆ = 26.5 km s−1, matching the peculiar velocity of the nearby O9.5 V star ζ Oph. We have also considered only one ionizing photon luminosity, again matching observations of this star.
The H ii region expansion leaves a cone-shaped shell in the wake behind the star, truncated on the H ii region edge at the angle where the normal velocity of gas through the I-front reaches Mach 2 (with respect to ionized gas). This shell should have column density of a few ×1020 cm−3 in neutral H and may be observable with kinematic H i data (or dust emission; Raga et al. 1997a). Directly downstream from the star is an underdense wake that should be filled by hot gas from the stellar wind (not modelled here). The gas expansion driven by the H ii region should make the stellar wind's wake longer-lived than would be the case without this expansion, because this creates a lower density and pressure environment.
The shocked shell is affected by the presence and orientation of a large-scale ISM magnetic field. Both H ii region expansion and shock compression are strongly inhibited perpendicular to magnetic field lines, leading to a less overdense shell with a much lower H i column density. The larger wave speeds in a magnetized plasma also give the shell a wider opening angle. The underdense wake, by contrast, is not significantly affected by the ISM magnetic field.
Kinetic energy is generated in the ISM from the H ii region's expansion at a rate comparable to the mechanical luminosity of the stellar wind from ζ Oph. When we take into account that (at least for static stars) about 90 per cent of the wind energy is lost to dissipation (e.g. García-Segura et al. 1996), the H ii region expansion has stronger feedback energetically than the stellar wind. This is also true of the momentum feedback, where the H ii region expansion generates more than 100 times the momentum in the ISM than the stellar wind. The momentum input rate is also about 10 times the total radiation momentum from the stellar luminosity. The kinetic energy feedback rate is affected somewhat (∼50 per cent level) by the ISM magnetic field, but momentum feedback is very little affected. Compared to the star's eventual supernova explosion, the total kinetic energy feedback from the H ii region over the star's main-sequence lifetime is >100 times less, but the momentum feedback is up to 4 times larger.
The unstable upstream I-front perturbs gas within the H ii region, but the perturbations have an amplitude of only ≈ 5-10 per cent near the star and are unlikely to affect the bow shock to any extent (although further simulations are needed to ensure that the perturbations do not grow to non-linear amplitudes with higher resolution).
A density gradient is created by the H ii region expansion, such that the ionized gas is densest near the upstream I-front and decreases downstream. This results in Hα emission maps that are brighter upstream from the star than downstream. In addition, the upstream I-front is very thin whereas the downstream recombination front is broad, and this is also be reflected in Hα emission maps. Similar features can be seen in Hα emission from the H ii region around ζ Oph. From this we conclude that gas dynamics and non-equilibrium ionization are important ingredients in the observational properties of H ii regions from supersonically moving exiles. Radiative transfer models that assume ionization equilibrium or ignore the relative motion between the star and the ISM may have only limited applicability to stars such as ζ Oph.
We have not studied the effects of ISM inhomogeneity and turbulent motions in this work. Rather we have provided a baseline study to see what the dynamics generated purely by the photoionization process are. In future work, we will include clumpy and turbulent ISM structure (cf. Mellema et al. 2006; Arthur et al. 2011) to investigate the degree to which this substructure is destroyed or enhanced by the passage of a hot star.
It does not seem to have been appreciated until now that H ii regions around exiled O stars can be very good laboratories for studying the physics of I-fronts, in particular their stability. The early R-type expansion of H ii regions around newly formed stars is deeply embedded in molecular clouds and so is difficult to observe, and both the I-front velocity and the stellar radiation spectrum are time-varying. Runaway stars, by contrast, are typically in low-extinction environments and have a constant velocity R-type I-front that (at least in regions of constant ISM density) should relax to a steady state. Some of them, for example ζ Oph, are also much closer to us than the nearest massive star-forming regions, so we can observe their I-fronts in much more detail and strongly test theoretical predictions and numerical simulations.
JM is funded by a fellowship from the Alexander von Humboldt Foundation. This work was supported by the Deutsche Forschungsgemeinschaft priority programme 1573, ‘Physics of the Interstellar Medium’. The authors acknowledge the John von Neumann Institute for Computing for a grant of computing time on the JUROPA supercomputer at Jülich Supercomputing Centre. We thank M.S. Oey for useful comments on a draft version of this paper. JM acknowledges the Nordita programme on Photo-Evaporation in Astrophysical Systems (June 2013), which enabled useful discussions about the results obtained in this work. We thank the referee, Alex Raga, for useful suggestions which improved the presentation of the paper.
Alexander von Humboldt Professor.
A temperature-dependent recombination rate is used, and the H ii region is not isothermal, but for a mean temperature T ≈ 6300 K we find Rs ≈ 10.6 pc.
The Southern H-Alpha Sky Survey Atlas (SHASSA) is supported by the National Science Foundation.
REFERENCES
APPENDIX A: RESOLUTION EFFECTS IN 1D SIMULATIONS
As discussed in Section 3, a CFL number of Ccfl = 0.1 was required for the simulations presented in the text. This choice is justified here, by varying the temporal and spatial resolution of the 1D simulations. The problem is that for |$\mathcal {M}>2$|, the H ii region structure relaxes to a stationary state with respect to the computational grid (even though gas is continuously recombining and being ionized), so the microphysics time-scales go to infinity. In this case, the CFL condition is the only active timestep restriction on the simulation. The effects of varying Ccfl on the upstream I-front position are shown in temperature plots for two different spatial resolutions (r1 and r3) in Fig. A1, for four different values of v⋆. The errors increase with v⋆ for a given Ccfl and Δx. The thin dotted line in the higher velocity figures shows the highest resolution model (r7) with Ccfl = 0.1; in this case the solution has converged, but we have had to spatially resolve the I-front (cf. Cantalupo & Porciani 2011). Even for v⋆ = 25 km s−1, the lowest resolution simulations have an incorrect I-front position with Ccfl = 0.3.
Fig. A2 shows the gas temperature through the upstream I-front at different spatial resolutions for Ccfl = 0.1 and v⋆ = 25 km s−1 (above) and 200 km s−1 (below). It shows that the temperature peak is well modelled for Ccfl = 0.1 for the low-velocity model, but for v⋆ = 200 km s−1 the solution lags behind the true solution at low spatial resolution. Once the I-front is spatially resolved (with Δτ ≤ 1 per zone) the solution has converged. The positions of the I-front are plotted for different Ccfl in Fig. A3, showing that even for v⋆ = 100 km s−1 the I-front position is modelled accurately with Ccfl = 0.1, but larger velocities are not. It is concluded that Ccfl = 0.1 provides nearly resolution-independent I-front positions and temperatures, at least for v⋆ < 100 km s−1.
Such constraints were not necessary for the test calculations in Mackey (2012) that did not include any hydrodynamics. In that case, the accuracy of the solution was independent of the I-front velocity, but for a static ISM and for an I-front that was moving rapidly with respect to the grid. It seems that the strong advection requires tight coupling to the microphysics integration to obtain an accurate solution (at least for this algorithm).
APPENDIX B: RESOLUTION EFFECTS IN 2D SIMULATIONS
A suite of 2D simulations have been run to test the effects of grid resolution on the H ii region feedback effects. They are analogous to the 3D simulations listed in Table 1 but with the same and higher resolutions, and are listed in Table B1. They consist of a 2D domain with x, y ∈ [ − 32.64, 18.56] pc with Cartesian (i.e. infinite slab) geometry. Fig. B1 shows simulations BT2r2 and BA2r2, having the same resolution and magnetic field configuration as the 3D models BT3r2 and BA3r2 in Section 4. Hydrogen number density |$n_{{\small h}}$| is plotted in colour, with contours showing the ion fraction (1 − yn), and solid lines crossing the domain showing the magnetic field orientation (in black). They show most of the same features at BT3r2 and BA3r2. The model with aligned spatial velocity and magnetic field is stable, but the model where they are perpendicular shows some instability, closer to the hydrodynamic result. The lateral expanding shell is more compressive in BT2r2 than BA2r2, and the opening angle of the shell downstream is narrower because the shock velocity is lower. Both models have similar minimum density, although the structure of the downstream wake is different because of the magnetic field orientation.
ID . | (Nx, Ny, Nz) . | Δx . | |$\boldsymbol {B}$|-field (μG) . |
---|---|---|---|
HD2r1 | (160, 160, 1) | 0.32 | (0, 0, 0) |
HD2r2 | (320, 320, 1) | 0.16 | (0, 0, 0) |
HD2r3 | (640, 640, 1) | 0.08 | (0, 0, 0) |
HD2r4 | (1280, 1280, 1) | 0.04 | (0, 0, 0) |
BA2r1 | (160, 160, 1) | 0.32 | (7, 0, 0) |
BA2r2 | (320, 320, 1) | 0.16 | (7, 0, 0) |
BA2r3 | (640, 640, 1) | 0.08 | (7, 0, 0) |
BA2r4 | (1280, 1280, 1) | 0.04 | (7, 0, 0) |
BT2r1 | (160, 160, 1) | 0.32 | (0, 7, 0) |
BT2r2 | (320, 320, 1) | 0.16 | (0, 7, 0) |
BT2r3 | (640, 640, 1) | 0.08 | (0, 7, 0) |
BT2r4 | (1280, 1280, 1) | 0.04 | (0, 7, 0) |
ID . | (Nx, Ny, Nz) . | Δx . | |$\boldsymbol {B}$|-field (μG) . |
---|---|---|---|
HD2r1 | (160, 160, 1) | 0.32 | (0, 0, 0) |
HD2r2 | (320, 320, 1) | 0.16 | (0, 0, 0) |
HD2r3 | (640, 640, 1) | 0.08 | (0, 0, 0) |
HD2r4 | (1280, 1280, 1) | 0.04 | (0, 0, 0) |
BA2r1 | (160, 160, 1) | 0.32 | (7, 0, 0) |
BA2r2 | (320, 320, 1) | 0.16 | (7, 0, 0) |
BA2r3 | (640, 640, 1) | 0.08 | (7, 0, 0) |
BA2r4 | (1280, 1280, 1) | 0.04 | (7, 0, 0) |
BT2r1 | (160, 160, 1) | 0.32 | (0, 7, 0) |
BT2r2 | (320, 320, 1) | 0.16 | (0, 7, 0) |
BT2r3 | (640, 640, 1) | 0.08 | (0, 7, 0) |
BT2r4 | (1280, 1280, 1) | 0.04 | (0, 7, 0) |
ID . | (Nx, Ny, Nz) . | Δx . | |$\boldsymbol {B}$|-field (μG) . |
---|---|---|---|
HD2r1 | (160, 160, 1) | 0.32 | (0, 0, 0) |
HD2r2 | (320, 320, 1) | 0.16 | (0, 0, 0) |
HD2r3 | (640, 640, 1) | 0.08 | (0, 0, 0) |
HD2r4 | (1280, 1280, 1) | 0.04 | (0, 0, 0) |
BA2r1 | (160, 160, 1) | 0.32 | (7, 0, 0) |
BA2r2 | (320, 320, 1) | 0.16 | (7, 0, 0) |
BA2r3 | (640, 640, 1) | 0.08 | (7, 0, 0) |
BA2r4 | (1280, 1280, 1) | 0.04 | (7, 0, 0) |
BT2r1 | (160, 160, 1) | 0.32 | (0, 7, 0) |
BT2r2 | (320, 320, 1) | 0.16 | (0, 7, 0) |
BT2r3 | (640, 640, 1) | 0.08 | (0, 7, 0) |
BT2r4 | (1280, 1280, 1) | 0.04 | (0, 7, 0) |
ID . | (Nx, Ny, Nz) . | Δx . | |$\boldsymbol {B}$|-field (μG) . |
---|---|---|---|
HD2r1 | (160, 160, 1) | 0.32 | (0, 0, 0) |
HD2r2 | (320, 320, 1) | 0.16 | (0, 0, 0) |
HD2r3 | (640, 640, 1) | 0.08 | (0, 0, 0) |
HD2r4 | (1280, 1280, 1) | 0.04 | (0, 0, 0) |
BA2r1 | (160, 160, 1) | 0.32 | (7, 0, 0) |
BA2r2 | (320, 320, 1) | 0.16 | (7, 0, 0) |
BA2r3 | (640, 640, 1) | 0.08 | (7, 0, 0) |
BA2r4 | (1280, 1280, 1) | 0.04 | (7, 0, 0) |
BT2r1 | (160, 160, 1) | 0.32 | (0, 7, 0) |
BT2r2 | (320, 320, 1) | 0.16 | (0, 7, 0) |
BT2r3 | (640, 640, 1) | 0.08 | (0, 7, 0) |
BT2r4 | (1280, 1280, 1) | 0.04 | (0, 7, 0) |
Results from simulations BT2r3 and BA2r3 are shown in Fig. B2, showing that with 2 times higher resolution both models have unstable I-fronts upstream from the star, and there are many more dense knots than were found in the 3D simulation. Fig. B3 has still higher resolution models BT2r4 and BA2r4, showing even greater degree of instability in the upstream I-front than BT2r3 and BA2r3, and also stronger compression in clumps and shell. R-type I-fronts have been shown (Newman & Axford 1967; Williams 1999) to be unstable to weak density perturbations via a shadowing instability. As long as the magnetic field is not sufficiently strong to make the gas flow entirely 1D, this instability should also be present in magnetized R-type I-fronts. The stability of the simulation BA2r2, and the similar 3D simulations BA3r1 and BA3r2, is at least partly an artefact of the increased stiffness of magnetic field lines at lower resolution.
The kinetic energy and momentum in these simulations is plotted in Fig. B4 in the same way as the 3D results in Fig. 9. The absolute values of the results are not normalized physically because the slab-symmetry implies a value per unit depth (here cm−1), but they are useful as a resolution study. The results are similar to those from 3D simulations; here the kinetic energy also increases slowly with spatial resolution and is not converging to a maximum value, because the I-front instability gets stronger with higher resolution.
The ISM density distribution function is plotted in Fig. B5, where by comparison to Fig. 8 we see that the higher resolution 2D simulations have stronger compression than their lower resolution counterparts. We expect a similar trend in 3D, if we could run such high-resolution models in 3D.
The properties of the ISM near the star are plotted in Fig. B6, which again should be compared to the 3D results in Fig. 10, although symmetry restrictions mean that the results are only indicative. There is a clear trend of increasing spatial and temporal variability with simulation resolution, also obvious from the simulation snapshots. The lowest resolution models are not shown because all relax to a steady state with only small variations near the star. For the next lowest resolution, BA2r2 relaxes to an almost-stationary state, whereas stronger I-front instability in BT2r2 and HD2r2 leads to continuing fluctuations for the duration of the simulation. These fluctuations are stronger again for models with resolutions r3 and r4. Only one averaging volume is plotted (r < 1 pc) but the trends at smaller volumes are the same as for 3D results.
In 2D, the mean density is unchanged near the star after the initial evolutionary phase (differing from 3D), and the mean velocity of gas drops by about 10 per cent (similar to 3D). At high resolution, the root-mean-squared (rms) density fluctuations are also ≈10 per cent, and the temporal fluctuations in mean density and velocity are up to 15 per cent. The H ii region does measurably affect gas properties near the star, but it is still at a level that can be treated as a small perturbation.
From the trends with resolution in the 2D simulations, it is clear that many quantities have not fully converged numerically, so we should be cautious in drawing strong conclusions regarding the compression of gas to high densities and the properties of the ISM near the star. The rate of momentum feedback seems quite robust, being barely changed when the spatial resolution changes by a factor of 8. The kinetic energy input also only changes by about 20 per cent over the same range of resolutions so, while it has not converged, it is not changing dramatically.
APPENDIX C: HEATING AND COOLING RATES FOR ATOMIC AND PHOTOIONIZED GAS
The heating and cooling functions used in the rate equation for energy (equation 5) should comprise the main coolants in WNM and photoionized gas. They are described in more detail here, term by term as they appear in equation (5). All terms have units erg s−1 and should be multiplied by |$n_{\small h}$| to get volumetric rates.
Heating processes
Cooling processes
Radiative recombination and bremsstrahlung cooling from H+ are interpolated from the tables for the total cooling coefficient in Hummer (1994), although an analytic fit would have been significantly more efficient. We assume that He is singly ionized when H is ionized, so we also add in a bremsstrahlung cooling term for He+ (its recombination cooling is neglected). These comprise the terms (Crr + Cff)ne(1 − yn) in equation (5), but neither term is important in photoionized gas at solar metallicity.
Collisional excitation cooling from H0, Ccxynne, is a strong coolant in partially ionized gas, and its rate is obtained by interpolating data from table 11 in Raga, Mellema & Lundqvist (1997b).
SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article:
Figure 4. Projections through the 3D simulations HD3r1 (left) and the higher resolution HD3r2 (right) at t = 4 Myr.
Figure 5. As Fig. 4, but for the simulations with a flow-aligned magnetic field, BA3r1 and BA3r2.
Figure 6. As Fig. 4, but for simulation BT3r2 with a magnetic field perpendicular to the gas flow, at t = 4.5 Myr, and from two different orientations (Supplementary Data).
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.