1. Introduction

Runaway O stars pass through the Galactic disk at velocities of ≈ 100 km s⁻¹ (see, e.g., ^{Cruz-González et al. 1974}). These stars interact with the interstellar medium of the disk, forming an inner stellar wind bow shock, and an outer “cometary” HII region. Observations of these structures are presented, e.g., by ^{Van Buren et al. (1995)} and ^{Noriega-Crespo et al. (1997)}. We note that these HII regions have a somewhat elongated “tear-drop” shape, and that the “cometary” label that is applied to them is somewhat misleading.

The problem of the cometary HII region for a star in hypersonic motion was solved analytically by ^{Rasiwala (1969)} and numerically by ^{Thuan (1975)}. Axisymmetric gasdynamic simulations of this flow were presented by ^{Raga et al. (1997)}.

The standard, ram-pressure balance stellar wind bow shock has the well used analytic solution of ^{Dyson (1975)}. Also, at a somewhat surprisingly later time, an analytic solution to the full, thin shell stellar wind bow shock problem (including the “centrifugal pressure”) was found (^{Wilkin 1996}; ^{Cantó et al. 1996}). The stellar wind bow shock problem has been studied extensively with gasdynamic numerical simulations (see, e.g., ^{Stevens et al. 1992}).

Stellar wind bow shock models have been used to model compact/ultracompact HII regions formed by young O stars within molecular clouds (e.g., ^{Van Buren 1990}; ^{Mac Low 1991}), moving at velocities of a few km/s through the surrounding cloud. The models for these objects generally assume that the ionization front is trapped within the stellar wind bow shock, so that the flow differs from the detached bow shock/HII region found in runaway O stars. Gasdynamic simulations of bow shock compact HII regions have been presented by ^{Arthur & Hoare (2006)}.

Simulations of bow shocks around runaway stars have been computed including thermal conduction (^{Meyer et al. 2014}) and magnetic fields (^{Meyer et al. 2017}). The cometary HII regions surrounding these bow shocks have been modeled in detail (in axisymmetry and 3D, including an environmental magnetic field) by ^{Mackey et al. (2013)}. Interestingly, even in the work of ^{Raga et al. (1997)}, which modeled both the bow shock and the cometary HII region of runaway O stars, these two structures of the flow have been considered separately. In the present paper we explore the regime in which the bow shock, while not trapping the ionization front, produces a substantial absorption of the stellar ionizing photons, therefore modifying the surrounding HII region.

We describe the flow by considering the thin shell analytic solution (of ^{Wilkin 1996}; ^{Cantó et al. 1996} and ^{Cantó et al. 2005}) for the stellar wind bow shock. A review of this solution is presented in § 2. We then write a modified form for the equation of ^{Rasiwala (1969)} and ^{Raga et al. (1997)}, including the absorption of ionizing photons by the bow shock (§ 3) and integrate it numerically to obtain the selfconsistent HII region+bow shock solution for a set of parameters (§ 4). Finally, in the summary (§ 5) we diScuss the possible application of the new model for the case of compact/ultracompact HII regions, and the further steps that will be necessary to study the full parameter space of these objects.

2. The Bow Shock

A star with an isotropic wind of mass loss rate Ṁ_W and (terminal) wind velocity v_w moving hypersonically at a velocity v_a through a uniform medium of density ρ_a produces a double bow shock (with one shock stopping the wind and the second shock pushing the environment). This is shown schematically in Figure 1, as seen in a reference system moving with the star.

Fig. 1 Schematic diagram showing a star with an isotropic wind interacting with a uniform, streaming environment (of density ρa and velocity va). A two-shock stellar wind bow shock with a locus R_b (θ) (two thick, solid curves) and a detached ionization front of locus R_IF (θ) (thick, dashed curve) are formed. The angles θ, α and γ (used in the derivation of the model equations) are shown. 

^{Wilkin (1996)} and ^{Cantó et al. (1996)} found an analytic solution for the steady, thin shell stellar bow shock problem:

with

with and where the polar angle θ and the spherical radius are shown in Figure 1.

Interestingly, the axisymmetric gasdynamic simulations of ^{Raga et al. (1997)} show that even though the bow shock driven by a runaway O star has a complex, time-dependent structure and a separation of ≈ 0.3R₀ between the two bow shocks, it still approximately follows the thin shell solution of equation (1). We therefore adopt this thin shell solution as a valid description of the stellar wind bow shock.

We also use the emission measure:

i.e., the square of the ion number density (n=ρ/m¯ for fully ionized H and singly ionized He) integrated across the width h of the thin shell. This integral was calculated by ^{Cantó et al. (2005)} for a bow shock travelling into an environment with a density gradient, and their results can be used (by setting the density gradient to zero) to obtain the emission measure for a bow shock in a uniform environment:

with

where m¯=1.3 m_H for gas with 0.9 H and 0.1 He number abundance, c₀ ≈ 10 km s⁻¹ is the sound speed of the photoionized gas, and f_a (θ), f_w (θ) and g_σ (θ) are given in Appendix A.

3. The Equation for The HII Region

For a star that emits S_* ionizing photons per unit time, if the bow shock shell does not trap the ionization front, we will have a cometary HII region surrounding the stellar wind bow shock (see the schematic diagram of Figure 1). For a hypersonic flow, the ionization front limiting this HII region will not produce a hydrodynamic perturbation, and the environmental material will therefore flow with density ρ_a and velocity ν_a until it intercepts the bow shock (see Figure 1).

The equation determining the locus R_IF (θ) of the ionization front is:

where the terms in the equation represent:

per unit time within the solid angle between θ and θ + dθ. In this equation, r_IF = R_IF sin θ is the cylindrical radius of the ionization front and α_H is the H recombination coefficient (we set α_H = 2.56 × 10⁻¹³cm³ s⁻¹).

Combining equations (3-6) we obtain the dimensionless equation:

where R_b/R₀ is given by equation (1),

is the Mach number of the impinging environment with respect to the isothermal sound speed of the photoionized gas,

is the Strömgren radius and

where na=ρa/m¯ is ion+atom environmental number density. The G(θ) function (describing the absorption of the bow shock shell as a function of θ) is given in Appendix A.

Our equation (7) is a generalization of equation (12) of ^{Raga et al. (1997)} (which has a dimensionless parameter ξ = χR₀/R_s) to the case of a cometary HII region with an inner stellar wind bow shock producing an appreciable absorption of the stellar ionizing photons.

For θ = 0 from equation (7) we have:

where we have used the fact that G(0) = 3/2 (see equation A17). This is a cubic equation which can be inverted (analytically or numerically) to obtain the on-axis radius R_IF,0 of the ionization front. This value can the be used to start a numerical integration of the ordinary differential equation (7) which then gives the shape R_IF (θ) of the ionization front.

Depending on the values of the parameters χ, R_S /R₀ and M_a, we find that equation (11) has either:

This latter situation corresponds to the case in which the head of the bow shock traps the ionization front,

4. Numerical Results

As there is a considerable number of free parameters, a general exploration of the possible parameter space is lengthy (and also not very enlightening). We therefore choose to explore a limited range of parameters, which is appropriate for the case of a runaway O star passing through the plane of the Galaxy.

For the runaway O star, we choose a O7V star with wind velocity v_w = 2500 km s⁻¹, mass loss rate Ṁw = 4.5 × 10⁻⁷M_⊙yr⁻¹ and ionizing photon rate S_∗ = 10⁴⁹s⁻¹. For the streaming environment we choose three possible velocities (v_a = 20, 100 and 300 km s⁻¹) and four possible number densities (n_a = 1, 100, 10⁴ and 10⁶ cm⁻³). The values of R₀, R_S/R₀ and χ corresponding to the chosen models are given in Table 1 (we do not list M_a = v_a /10 km s⁻¹).

For the combinations of these environmental velocities and densities we first calculate the axial standoff distance R_IF,0 of the ionization front from equation (11). The results of this exercise are given in Table 1. It is clear that for all of the chosen parameters the bow shock does not trap the ionization front, and that an external HII region is produced. The smaller HII region (relative to the size of the bow shock) is found for the v_a = 20 km s⁻¹, n_a = 10⁶ cm⁻³ model, in which the on-axis radius of the ionization front is 3 times the bow shock radius (see Table 1).

With the R_IF (θ = 0) = R_IF,0 initial condition, we integrate equation (7) to obtain the shape of the ionization front. The results for the v_a = 20 km s⁻¹ models are shown in Figure 2, the ones for the v_a = 100 km s⁻¹ models in Figure 3, and the results for the v_a = 300 km s⁻¹ models in Figure 4,

Figures 2-4 show the stellar wind bow shock (equation 1), the ionization front R_IF obtained from equation (7), and the ionization front R_I obtained from equation (12) of ^{Raga et al, (1997)}. R_I is the ionization front that would be obtained for a star producing a S_∗ ionizing photon rate, but with no wind (and therefore no stellar wind bow shock). Interestingly, in all but the v_a = 20 and 100 km s⁻¹, n_a = 10⁶cm⁻³ model, no difference can be seen in the leading hemisphere of the HII region (i.e., for z < 0).

Fig. 2 The ionization front (thick, solid outer curve) and the stellar wind bow shock (thin, solid inner curve) obtained for models with v_a = 20 km s⁻¹ and n_a = 1 (top frame), 100, 10⁴ and 10⁶ cm⁻³ (bottom). The dashed curves show the ionization front that would be obtained for a star with no wind. The (z, r) axes are given in units of R₀. 

Fig. 3 The same as Figure 2, but for the v_a = 100 km s⁻¹ models (see Table 1). 

Fig. 4 The same as Figure 2, but for the v_a = 200 km s⁻¹ models (see Table 1). 

Substantial differences between R_IF and R_I can be seen in the trailing hemisphere. The absorption of ionizing photons by the bow shock wings and the lack of absorption within the stellar wind region lead to the formation of either:

Except for the trailing region close to the bow shock wings, and except for the v_a = 300 km s⁻¹, n_a = 1 cm⁻³ model (top frame of Figure 4, which has a clearly non-spherical HII region), the ionization front always has an almost spherical shape, with a radius close to the Strömgren radius R_S for a star with S_∗ in a stationary medium of density n_a (given in the fourth column of Table 1). The sphere, however, is not centered on the position of the star, but on a position downstream of the star. This effect is described by the “χ << 1 solution” of equation (13) of ^{Raga et al. (1997)}, which to first order in χ corresponds to a sphere of radius R_S centered on a point with a downstream offset of χR₀/3 with respect to the position of the star.

5. Discussion

Our exploration of the parameter space for the interaction of a runaway O7 star with the surrounding ISM shows that there is:

It is clear that we have not explored parameters for which the ionization front is trapped by the head of the stellar wind bow shock. Which parameters would be necessary for doing this? The condition for ionization front capture at the head of the bow shock can be obtained by setting R_IF,0 = R₀ in equation (11), which gives:

where R_S is given by equation (9), R₀ by equation (2), M_a by equation (8) and χ by equation (10) as functions of the physical parameters of the flow. If we fix the stellar parameters for an O7V star, the remaining free parameters are the density and velocity of the streaming environment. Equation (12) then gives us the minimum environmental density n_a,c as a function of velocity v_a necessary for the ionization front to be captured by the head of the bow shock. The results of this exercise are shown in Figure 5.

We see that for the v_a = 1 → 40 km s⁻¹ range shown in Figure 5, the critical density is in the n_a,c ∝ v_a^1/2 regime obtained from equation (12) for << χ 1. It is clear that for v_a 1 km s⁻¹, large environmental densities are required in order for the bow shock to capture the ionization front!

Fig. 5 Minimum environmental density n_a,c (for capturing the ionization front) as a function of velocity v_a for the parameters of an O7V star (see the text). 

For v_a >10 km s⁻¹ stellar motions, n_a,c >10⁷ cm⁻³ environmental densities are required, so that it is unlikely that one will ever find such objects with ionization fronts captured by the stellar wind bow shocks. For lower values of v_a, one could have captured ionization fronts provided that the objects travel within a dense molecular cloud (of densities ≈ 105 → 10⁷ cm⁻³, see Figure 5).

We are, of course, now talking about compact/ultracompact HII regions within dense molecular clouds, produced by massive stars with velocities of a few km/s with respect to the surrounding cloud. Interestingly, as we see in Figure 5, we could also have O stars with v_a = 1 → 10 km s⁻¹ with bow shocks that do not capture the ionization fronts. This parameter regime is interesting because it is not correctly described by our present model.

For v_a < 10 km s⁻¹, the steady bow shock+HII region configuration (for a case in which the IF is not captured by the bow shock) will have a motion through the neutral environment that is subsonic with respect to the c₀ ≈ 10 km s⁻¹ isothermal sound speed of the ionized gas. Such an ionization front will have a strong hydrodynamical effect, driving a shock wave ahead of it into the neutral environment. The post-shock material will be almost at rest with respect to the star, and therefore will not shock again before reaching the contact discontinuity against the stellar wind.

The flow configuration will then have an inner shock involving the stellar wind, and an outer shock surrounding an extended HII region. This outer bow shock will be somewhat peculiar, since all of the post-shock flow (even far along the bow shock wings) will be subsonic (because the pre-shock flow is subsonic with respect to c₀). Therefore, the post-bow shock flow will not have the typical division between a compact axial subsonic region and an extended supersonic region farther away from the flow axis.

In a future paper we will continue with a study of this interesting regime (already explored by ^{Arthur & Hoare 2006}) with axisymmetric gasdynamic+radiative transfer simulations. Clearly, the codes needed to do this have already been fully developed (see, e.g., ^{Mackey et al. 2021}).