RADIATION-HYDRODYNAMIC MODELS OF THE EVOLVING CIRCUMSTELLAR MEDIUM AROUND MASSIVE STARS

Our simulations are carried out in one and two dimensions. First, the time-dependent, spherically symmetric system of radiation-hydrodynamic equations is solved for the MS and RSG phases, during which the stellar winds are reasonably steady and instabilities are relatively unimportant. Then, the 1D results are remapped to a 2D, cylindrically symmetric (r, z) grid, to continue the evolution of the CSM through the WR stage, where the density gradients and wind–wind interactions lead to instabilities, which severely disrupt the CSM.

The 1D spherically symmetric Eulerian hydrodynamic conservation equations are solved using a second-order, finite volume Godunov-type scheme with outflow-only boundary conditions at the outer boundary. In 2D, the cylindrically symmetric hydrodynamic conservation equations are solved using a hybrid scheme, in which alternate time steps are calculated by a Godunov method and the van Leer flux-vector splitting method (van Leer 1982). This reduces numerical artifacts, which can be produced in Godunov-type schemes when slow shocks are propagating parallel to one of the grid directions. The boundary conditions in the 2D scheme are those of no flux through the axis of symmetry and outflow-only conditions at the outer radial and z boundaries.

The 1D simulations start on a grid of 1000 cells with a physical size of 5 pc, which is large enough to contain the initial Strömgren sphere. We use an expanding grid to avoid large computational domains at the beginning of the simulation: when the outer shock nears the edge of the grid, more (uniform ambient medium) cells are added. Each 1D calculation has the same physical resolution (grid cell size dr = 0.002 pc). By the end of the RSG stage, the stellar wind bubble, H ii region and swept-up neutral shell extend up to 30 pc, meaning that the corresponding grids have of order 6000 cells. The full evolution of the CSM around the massive stars is followed in 1D from the beginning of the MS stage through to the WR stage. However, once the star enters the WR stage, 1D is not adequate to capture the full richness of the structures formed during the wind–wind interactions and so we remap the results onto a 2D grid to continue the calculation in cylindrical r − z coordinates.

The majority of the 2D simulations are carried out on a fixed grid of 500 radial by 1000 z-direction cells with a corresponding physical size of 10 × 20 pc² (i.e., cell size dr × dz = 0.02 pc× 0.02 pc). This size is chosen because the wind–wind interaction takes place within this distance for almost all models, and the radius of observed WR bubbles is R ⩽ 10 pc (Gruendl et al. 2000). One model is run on a slightly larger grid (750 × 1500 cells) because the main wind–wind interaction occurs further from the star, but the spatial resolution is the same.

The transfer of the ionizing radiation is carried out by the method of short characteristics (Raga et al. 1999) adapted for spherical and cylindrical symmetry. We consider only a single point source for the ionizing radiation (the star), which is positioned on the symmetry axis. The hydrodynamics and radiative transfer are coupled through the energy equation, where the source term depends on the photoionization heating and the radiative cooling rate, and through advection equations for the ions and neutrals. The ionization balance equation takes into account both photoionization and collisional ionization together with recombination. This basic code has been extensively tested and applied to the evolution of H ii regions in density gradients (Henney et al. 2005; Arthur & Hoare 2006).

The radiative cooling function takes into account both cooling in collisionally ionized gas (in the hot shocked wind bubble) and in the photoionized gas in the H ii region where collisional ionization of hydrogen is not an important contribution to the cooling rate. Cooling rates are generated using the Cloudy photoionization code (Ferland et al. 1998) for appropriate stellar parameters and used to generate a look-up table. The cooling rates extend down to temperatures of 100 K but photoionization heating ensures that the H ii region has a temperature of ∼10⁴ K.

The stellar wind is input to the grid through the Chevalier & Clegg (1985) thermal wind model, as a volume-distributed source of mass and energy. This avoids the problem of reverse shocks that can occur when there is a decrease in the ram pressure of a stellar wind input in the form of a region of high-density, high-velocity flow near the origin. The volume of the stellar wind injection region is chosen to be sufficiently small so that the stellar wind always reaches its terminal velocity and shocks outside the source region, but sufficiently large so as to be a reasonable geometrical representation of a sphere on a cylindrical grid (in the case of the 2D numerical simulations).

The variation of the stellar mass-loss rate with time is taken from publicly available stellar evolution models (Meynet & Maeder 2003; Eldridge et al. 2006). The stellar wind velocities are calculated using the tabulated stellar parameters (mass, radius, effective temperature, and surface abundances) and the ionizing photon rates are obtained from the appropriate stellar atmosphere models distributed with the Starburst 99 code (Leitherer et al. 1999). All of this information is stored in tables and interpolated as needed using the elapsed simulation time.

We also run a set of models that includes the effects of electron thermal conduction to assess the importance of this physical process for the dynamics of stellar wind bubbles and their observable properties, such as X-ray emission (Weaver et al. 1977). Thermal conduction is taken to be a diffusion process, with the heat flux, $\vec{q}$, given by

$$\begin{equation} \vec{q} = - D \nabla T_{e}, \end{equation} \tag{ 1 }$$

where D is the diffusion coefficient and T_e is the electron temperature (Cowie & McKee 1977). From Spitzer (1962), we have that the electron mean free path can be written as

$$\begin{equation} \lambda _{e} = 2.625 \times 10^{5} T_{e}^{2} /(n_{e} \ln \Lambda), \end{equation} \tag{ 2 }$$

where n_e is the electron number density and ln Λ is the Coulomb logarithm. For a pure hydrogen plasma this can be approximated by

$$\begin{equation} \ln \Lambda = 9.452 + \case{3}{2} \ln T_{e} - \case{1}{2} \ln n_{e}, \end{equation} \tag{ 3 }$$

when T_e ⩽ 4.2 × 10⁵ K, and

$$\begin{equation} \ln \Lambda = 15.90 + \ln T_{e} - \case{1}{2} \ln n_{e}, \end{equation} \tag{ 4 }$$

when T_e ⩾ 4.2 × 10⁵ K.

From Steffen et al. (2008), the diffusion coefficient D is then given by

$$\begin{equation} D = 7.04 \times 10^{-11} \lambda _{e} n_{e} T_{e}^{1/2}, \end{equation} \tag{ 5 }$$

where the numerical coefficient is the result of evaluating physical constants (Spitzer 1962). In spherical symmetry, the nonlinear diffusion equation is

$$\begin{equation} \frac{\partial e}{\partial t} = \rho c_{v} \frac{\partial T_e}{\partial t} = \frac{1}{r^{2}} \frac{\partial }{\partial r} \left(r^{2} D \frac{\partial T_e}{\partial r} \right), \end{equation} \tag{ 6 }$$

and we solve this using a Crank-Nicholson-type scheme (Press et al. 1992). Thermal conduction is also included in the 2D simulations by solving the corresponding diffusion equation in cylindrical symmetry. The effect of conduction is included as an update to the internal energy, e, once the diffusion equation has been solved for the electron temperature, and we assume equipartition between the electrons and ions.

When the electron temperature T_e is high and the electron density n_e low, the electron mean free path λ_e becomes large compared to the size scale of the hot bubble, and the diffusion approximation breaks down. This is the saturated heat flux limit (Cowie & McKee 1977) that can be described as

$$\begin{equation} q_{\mathrm{sat}} = 1.72 \times 10^{-11} T_{e}^{3/2} n_{e}. \end{equation} \tag{ 7 }$$

Saturated conduction is treated by limiting the electron mean free path in the following manner:

$$\begin{equation} \lambda = \min \lbrace 0.244 \,\Delta r, 2.625 \times 10^{5} T_{e}^{2} / n_{e} \ln \Lambda \rbrace, \end{equation} \tag{ 8 }$$

where Δr is the local grid spacing and the numerical constant ensures that the conductive heat flux can never exceed the saturation limit (Steffen et al. 2008). This formulation means that the limiting flux can only be reached at the sharp edge of the hot bubble, where |ΔT_e| ≈ T_e. A drawback of this method of limiting the electron mean free path is that it depends explicitly on the numerical resolution Δr and so our results will be affected by this. However, since all our 1D results are run at the same numerical resolution, comparisons between simulations are perfectly valid.

The temperature found from solving the diffusion equation for every grid cell is used to update the internal energy of that cell once the hydrodynamic and radiation steps have been completed.

4.1.1. Main-sequence Stage

We start from an initial condition of a cold, uniform, neutral medium with density n₀ = 100 cm⁻³ and temperature T = 100 K. The photoionized (H ii) region forms first. The photoionized gas is heated to T ∼ 10⁴ K and begins to expand into the surrounding medium at v_if ≳ 10 km s⁻¹. The conditions at the ionization front require a shock to form in the neutral medium ahead of it, which sweeps up and compresses the neutral gas. Inside the H ii region, the highly supersonic stellar wind forms a two-shock structure: an outer shock sweeps up, heats, and compresses the photoionized gas, and an inner shock brakes and heats the stellar wind. The two regions have the same pressure but different densities and are separated by a contact discontinuity. The outer shocked region can cool down quickly to the temperature of the H ii region due to the high densities in the swept-up material, but its temperature is maintained at 10⁴ K by photoionization. The shocked stellar wind, on the other hand, has a very high temperature (T > 10⁷ K) and low density and so does not cool efficiently. A hot bubble forms inside the H ii region, which is in turn surrounded by a slowly expanding shell of swept-up neutral material. These structures are well known from textbooks (e.g., Dyson & Williams 1997) and the slow time evolution of the stellar parameters does not cause any noticeable modifications.² The evolution of a stellar wind bubble inside an H ii region is different to that when ionization is not taken into account (e.g., van Marle et al. 2004). Once the initial Strömgren region has been formed, the pressure in the H ii region is much higher than that in the surrounding ISM. Moreover, although to begin with the pressure in the hot, shocked, stellar wind bubble is higher than that in the photoionized gas, pressure equilibrium with the H ii is reached after a short time and, thereafter, the pressure across the bubble structure is regulated by both the H ii region and the stellar wind together.

To illustrate the MS evolution, in Figure 5 we show the results from a 40 M_☉ STARS model for simulations with and without thermal conduction after t = 10⁶ yr. The hot, shocked wind bubble, H ii region, and neutral shell are clearly identified. The bubble radius in the simulation with thermal conduction is slightly smaller than that without conduction due to increased radiative losses in the hot bubble. These enhanced radiative losses in the hot, stellar wind bubble reduce the pressure here and enable the H ii region to expand back toward the star, and the photoionized gas has slightly lower density in the thermal conduction simulations that in those without conduction as a result.

Zoom In Zoom Out Reset image size

Table 1. Radii and Masses of the Outer H i Shell at the End of MS Phase

Mass	Model	$R_{\mathrm{H\,{\mathsc {i}}}}$	$M_{\mathrm{{\rm ISM}}}$
(M☉)		(pc)	(M☉)
40	MM2003 N	24.6	2.0 × 105
40	MM2003 R	27.5	2.8 × 105
40	STARS	29.0	3.3 × 105
60	MM2003 N	25.6	2.2 × 105
60	MM2003 R	29.4	3.4 × 105
60	STARS	29.4	3.4 × 105

Download table as:  ASCIITypeset image

For our particular models, we find that by the end of MS stage, the radius of the region affected by the massive star is 25–30 pc (see Table 1) and the mass of swept-up neutral material is >10⁵ M_☉. The exact size of the bubble depends mainly on the length of the MS stage (e.g., models with rotation have longer MSs than those without) and the stellar wind velocity. Higher stellar wind velocities lead to higher pressures in the shocked wind bubbles and therefore faster expansion speeds (e.g., the STARS models have higher MS stellar wind velocities, see Figure 2). We remark that bubbles expanding into denser media will, of course, be smaller, since on dimensional grounds the external radius of a stellar wind bubble in an uniform medium of density n_o is $R \sim (\dot{M_{W}}V_{w}^{2} / n_{0})^{1/5} t^{3/5}$, where $\dot{M_{W}}$ is the mass-loss rate and V_W is the stellar wind velocity (see, e.g., Dyson & Williams 1997), while that of an expanding H ii region follows R_i ∼ (S_*/α_Bn²₀)^3/7c^4/7_IIt^4/7, where α_B and c_II are the recombination coefficient and the sound speed in the ionized gas, respectively, and S_* is the stellar ionizing photon rate. Both of these relations depend inversely on the density (Spitzer 1978).

4.1.2. Post-main-sequence Stage

During the RSG phase the star loses a large amount of its mass in a slow, dense wind (see Figure 1). This wind is subsonic with respect to the hot wind bubble left by the MS stellar wind and so piles up at the bottom of a r⁻² density gradient into a thin dense shell without forming a two-shock pattern.

In Figure 6, we show the density, temperature, and internal energy (i.e., thermal pressure, since e = p/[γ − 1]) as a function of radius at the end of RSG phase for the 40 M_☉ models. We compare the results using the MM2003 stellar evolution model (without rotation) with those of STARS. For the simulations without thermal conduction, we see an inner region of RSG wind, which has a r⁻² density profile and temperature T < 10⁴ K, with a high central density, and bounded by a dense shell of partially ionized material. For the MM2003 model, this region has a radius of ∼6 pc, while for the STARS model the radius is only ∼2 pc, the shell is much narrower and the densities are higher. This is because the RSG phase is much shorter in the STARS models, even though roughly the same amount of stellar mass is lost as in the MM2003 models. The ram pressure of the RSG wind is less than the thermal pressure of the hot, shocked bubble produced during the MS stage. As a result, there is some back filling as the hot diffuse gas pushes back toward the star. The lack of ionizing photons at this stage (see Figure 3), and enhanced opacity due to the RSG wind and neutral shell, leads to recombinations in the outer H ii region at radius R ∼ 27 pc. This is most pronounced for the MM2003 model, since the inner absorbing RSG shell contains more neutral atoms than the shell formed in the STARS case. When thermal conduction is included in the simulations, evaporation of cold material from both the RSG shell and the outer H ii region enhances the density and thus the cooling rate in the hot shocked wind bubble. In the MM2003 case, this results in a much smaller hot bubble (∼19 pc as opposed to ∼27 pc in the simulation without conduction). The CSM in the STARS simulations is less affected by thermal conduction, possibly due to the shorter timescales in the RSG phase.

Zoom In Zoom Out Reset image size

In Figure 7, we show the corresponding circumstellar structures for the 60 M_☉ models at the end of the RSG phase. The MM2003 and STARS results look quite different. This is because the MM2003 model experiences several episodes of enhanced mass loss, whereas the STARS model undergoes a single intense outburst. The ram pressure of the slow wind in the MM2003 case fluctuates and this leads to acoustic waves propagating from the inner edge of the hot bubble outward. In Figure 7, the waves have yet to reach the outer edge of the bubble, but once they do, they will reflect from the dense neutral shell and pass back through the bubble and reflect backward and forward through the hot shocked gas, forming localized heated and compressed regions where oppositely directed waves shock against each other. There is no clearly defined shell of RSG material in the 60 M_☉ models, because the duration of the most intense period of mass loss is very short and just before the onset of the WR phase (see Figure 1), hence the dense material has not had time to move far away from the star before the fast wind starts.

Zoom In Zoom Out Reset image size

We are particularly interested in the radii of the dense, neutral RSG shells formed in these models. These shells will be impacted and broken up by the fast WR winds and so have important consequences for the observables (X-ray and optical emission) of the WR wind-blown bubbles. The shells are most evident in the 40 M_☉ models, reaching radii of ∼6 pc and ∼2 pc for the MM2003 and STARS models, respectively, by the end of the RSG phase. These distances reflect the duration of the RSG stage in each case: 5.6 × 10⁵ yr for the MM2003 model as opposed to 3.1 × 10⁵ yr for the STARS model. Also, the wind velocities in the STARS model are allowed to fall below 30 km s⁻¹ and in general are lower than those of the MM2003 models (see Figure 2). We have not presented results for the MM2003 models with stellar rotation because they are qualitatively similar to the results for the other models and result in a dense neutral shell being formed at ∼2 pc, since the RSG phase in this model only lasts 1.0 × 10⁵ yr. We continue the simulations in 1D until the end of the WR phase but we do not show these models because we are interested in the formation of instabilities during the fast-wind–slow-wind interaction, which do not occur in 1D.

At the end of the RSG phase, we remap our 1D, spherically symmetric distributions of density, momentum, and energy onto a 2D axisymmetric grid using a volume-weighted averaging procedure to ensure conservation of these quantities. The 2D grid is, of computational necessity, coarser than the 1D grid, and we also zoom in on the inner ∼10 pc of the simulation, since this is where the hydrodynamic interaction between the fast WR wind and the slow RSG wind will take place. Also, the observed WR bubbles all have radii 10 pc or less. We continue to follow the radiation-hydrodynamic evolution of the CSM using the stellar wind parameters and the tabulated ionizing photon rate as a function of time shown in Figures 1–3. For all 2D simulations, we present only one quadrant of the calculation, since the other is the same due to symmetry.

Zoom In Zoom Out Reset image size

4.2.1. 40 M_☉ Models

In Figure 8, we present gray-scale images for the wind–wind interaction for the 40 M_☉ STARS model without thermal conduction at four different times corresponding to when the interaction region has reached distances 2, 4, 6, and 8 pc from the central star. These times are, respectively, 39, 800, 65, 900, 86, 000, and 104, 000 yr after the onset of the WR wind. In this figure, we see the disruption of the RSG shell by the fast wind and the formation of clumps. The RSG material consisted initially of a thin, dense shell of partially ionized material located ∼2 pc from the star. The fast WR wind accelerates down the density gradient and slams into the dense shell. The WR wind is shocked and heated by this interaction while the dense shell is shocked and compressed. Vishniac instabilities of the linear thin-shell variety break the shell into dense clumps, which then travel radially outward (Vishniac 1983; Vishniac & Ryu 1989). The clumps have a cometary form, with the dense, cold (10⁴ K photoionized) head pointing toward the central star and tails extending several parsecs radially outward, showing that the heads are traveling more slowly than the diffuse material that surrounds them. Interactions between clump material and the shocked fast wind create a pattern of multiple interacting shocks, which heat the ablated clump gas to X-ray temperatures. For this particular case the RSG shell breaks up at ∼5 pc (see Figure 8).

Figure 9 shows the results obtained when the stellar wind parameters come from the MM2003 40 M_☉ model without rotation and without conduction. The wind–wind interaction region is again shown at 2, 4, 6, and 8 pc from the central star, corresponding to 10, 600, 22, 600, 32, 650, and 36, 600 yr after the onset of the WR wind. For this model, the RSG dense wind material, traveling at ∼30 km s⁻¹, is in the form of a broad shell about 6 pc from the star at the start of the WR stage. The interaction of the WR wind with this shell occurs only a few times 10³ yr after the onset of the fast wind. Compared to the STARS model, the MM2003 models experience Rayleigh–Taylor instabilities, which disrupt the accelerating WR shell before it interacts with the main shell of RSG material, forming small clumps at <2 pc. The main interaction of the WR shell with the RSG shell occurs at 6 pc and compresses the RSG material into a very thin, dense shell, which then corrugates and starts to break up at around 8 pc due to the linear thin-shell instability.

Zoom In Zoom Out Reset image size

Figure 10 shows the results from MM2003 model with stellar rotation but without conduction at 6800, 8800, 12, 800, and 16, 900 yr after the onset of the WR wind. In this model, the interaction takes place much further from the star and the shell of RSG material is not so dense, so the interaction does not produce such violent instabilities and the shell remains essentially intact, being swept up and compressed by the fast wind but not suffering the complete disruption seen in the STARS case discussed above. For this model, the intense mass-loss phase has a relatively short duration, ∼10⁵ yr, and there are several mass-loss episodes before and after the main ejection event. Rather than an RSG stage, the extremely strong mass loss ($\dot{M}_{w} \sim 10^{-3.6}\,M_{\odot }$ yr⁻¹), together with the fact that the stellar effective temperature is relatively high (T_eff ∼ 10⁴ K) for part of this period, indicate an LBV stage, with stellar wind velocities V_w > 200 km s⁻¹. The densest wind material travels quite a large distance from the star before the onset of the WR stage and so the interaction with the fast wind, when it occurs, does not lead to complete breakup of the shell.

Zoom In Zoom Out Reset image size

We find that thermal conduction has little discernible effect either on the density or the temperature structure of the 2D simulations and the results are essentially the same as those shown in Figures 8–10. In Figure 11, we present one set of results including thermal conduction to illustrate this. In this figure, we show the results of the 40 M_☉ MM2003 simulations (without rotation) at 36,600 yr after the onset of the WR wind. There are slight differences between the two sets of results, but the short timescales of the wind–wind interaction stage compared to cooling timescales in the 10⁶ K gas mean that there is no major effect on the hydrodynamics due to thermal conduction during this stage.

We have run further set of STARS simulations to investigate the effects of the radius of the wind injection zone and the grid resolution on the number and distribution of clumps formed during shell disruption. Increasing the grid resolution to 800×1600 cells while maintaining the same number of cells in the wind injection region (radius 30 cells) produced the same results as those presented here. Decreasing the size of the wind injection region (radius 10 cells) did produce fewer clumps but increasing the radius of the wind injection region to 50 cells gave the same results as in Figure 8. The clumps formed in our simulations have diameters between 15 and 50 cells and are photoionized, hence are fully resolved.

Zoom In Zoom Out Reset image size

4.2.2. 60 M_☉ Models

In Figure 12, we show the corresponding results for the 60 M_☉ STARS model without conduction. Again, we show the wind–wind interaction at 2, 4, 6, and 8 pc, corresponding to evolution times of 27, 150, 43, 200, 57, 200, and 71, 250 yr after the onset of the WR wind. The dense RSG (or LBV) material is swept up into a thin shell by the fast WR wind. Linear thin-shell instabilities are induced in the swept-up material and break the shell into many cold (<10⁴ K), dense (>1 cm⁻³) clumps which drift outward, embedded in the hot, shocked WR wind. The morphologies of the 40 M_☉ and 60 M_☉ STARS models are very similar.

Zoom In Zoom Out Reset image size

For the MM2003 60 M_☉ model without rotation, there are several episodes of enhanced mass loss during the RSG/LBV stage, and even in the WR stage there are short-lived outbursts (see Figure 1). This leads to the formation of successive shells of material traveling with different velocities away from the star shown in Figure 13. The most intense mass-loss episode shown in Figure 1 is followed by a reduction in mass-loss rate and increase in wind velocity (see Figure 2), which sweeps up, compresses and accelerates the ejected material. This leads to the formation of Rayleigh–Taylor instabilities, which disrupt the shell once it has reached ∼10 pc from the star. The clumpy shell is traveling quite quickly and soon leaves the region of interest (R < 10 pc). Other, less dense, shells follow, which suffer instabilities but do not contain enough material to form dense, cold clumps. In Figure 13, we show the wind–wind interaction and the formation of different shells and shock patterns as a result of the variability of the mass-loss and wind velocity up to 30,900 yr after the onset of the WR wind.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the MM2003 60 M_☉ model with rotation, from Figure 1 we note that there is no period of particularly enhanced mass loss. In this case, therefore, although the RSG is swept up by the faster WR wind, the swept-up material is not particularly dense and does not form a thin, dense shell nor any instabilities or clumps. In Figure 14, we show the density and temperature in the evolving CSM up to the point when the swept-up material has reached a radius of 8 pc from the star, which occurs 25,500 yr after the onset of the WR stage.

Differences between our results and the work of previous authors can be attributed to several factors: the stellar evolution models, the numerical method, and the initial conditions and physics that is included in the numerical simulations. In this section, we discuss each of these aspects in turn.

Much of the previous numerical work on the evolving CSMs around massive stars has used the same 35 M_☉ and 60 M_☉ stellar evolution models as García-Segura et al. (1996a, 1996b), which are those developed by Langer et al. (1994). These models were used in their full time-dependent form by García-Segura et al. (1996a, 1996b), Freyer et al. (2003, 2006), and Dwarkadas (2007). For both stellar masses considered, the stars evolve to be RSGs (T_eff < 4800 K). One difference between these models and the ones we use in this paper is the stellar wind velocity. In the García-Segura et al. (1996a, 1996b) models, the stellar wind velocity remains high, above 2500 km s⁻¹, throughout the MS stage, whereas in all of our models we see a significant decrease from an initial ∼3000 km s⁻¹ to ∼1000 km s⁻¹. This affects both the pressure and size of the hot, shocked wind bubble at the end of the MS stage. In the case of their 60 M_☉ model, there is an "H-rich WN" stage immediately after the MS stage, in which more than 10 M_☉ of material is lost by the star at high velocity over the space of a million years. None of the models we have studied have a similar evolutionary stage. This mass is moving too fast to remain close to the star during the subsequent evolutionary stages, and thus does not contribute to structure formation in the CSM. Other authors have used the STARS stellar evolution models (Eldridge et al. 2006), where the only differences between that work and ours is that we have corrected the stellar wind velocities (upward) by missing numerical factors, and those of Schaller et al. (1992), for 40 M_☉ and 60 M_☉ (van Marle et al. 2005, 2007), which are very similar to the MM2003 models without rotation. However, van Marle et al. (2005, 2007) approximates the evolution with a three-constant-winds scenario, and so does not follow the episodic ejections of material which are a main characteristic of the post-MS evolution of stars with initial masses ≳ 60 M_☉ (Humphreys 2010). By following the full time variation of the stellar wind parameters, we are able to study the formation and interaction of multiple shells in the CSM in the LBV and WR stages.

All of the numerical simulations by previous authors use well-tested hydrodynamical schemes and so we do not expect major differences to arise from differences between codes. The codes used include ZEUS (Stone & Norman 1992; García-Segura et al. 1996a, 1996b; Ramirez-Ruiz et al. 2005; van Marle et al. 2005; Eldridge et al. 2006; van Marle et al. 2007; Chita et al. 2008; Pérez-Rendón et al. 2009), VH1 (Dwarkadas 2007), and a code due to Yorke & Kaisig (1995) and Yorke & Welz (1996) (Freyer et al. 2003, 2006). Some authors choose spherical coordinates for their 2D models (e.g., García-Segura et al. 1996a, 1996b; van Marle et al. 2005; Dwarkadas 2007; van Marle et al. 2007), while others use cylindrical r − z coordinates, as we do (e.g., Freyer et al. 2003, 2006). Each coordinate system has its advantages and disadvantages. In 2D spherical coordinates, the cell sizes increase with increasing radius, while in r − z coordinates it is more difficult to represent a spherical wind injection zone.

One key ingredient, which can strongly affect the results of the numerical simulations, is the inclusion (or not) of radiative transfer and the modeling of the H ii region that will develop around the star during the MS and WR stages. Of the previous numerical work we have mentioned, only Freyer et al. (2003, 2006), van Marle et al. (2005, 2007), and Chita et al. (2008) (for a 12 M_☉ star) include the effect of ionizing photons on the evolution of the CSM. Both van Marle et al. (2005, 2007) and Chita et al. (2008) use a simple Strömgren radius approach to the radiative transfer, whereby gas is either ionized or neutral depending on whether the ionizing photons reach it or not (it is not clear whether these authors take collisional ionization into account). Freyer et al. (2003, 2006) uses a ray-tracing method and includes collisional ionization. Our radiative transfer method (see Section 2.1) allows us to follow gas at intermediate stages of ionization. For example, once the RSG stage begins, the ionizing photons do not reach the outer radii of the nebula, but the recombination timescale is of order 10⁵/n yr and thus the gas does not recombine immediately. This has important consequences for the pressure across the bubble. In our results, even once the ionizing photons have switched off, we find roughly constant pressure right across the bubble, from the inner edge of the hot, shocked wind through the recombining H ii region and the outer neutral shell, since acoustic waves smooth out the pressure distribution as the H ii region slowly recombines. However, the results presented by van Marle et al. (2004, Figure 4) and van Marle et al. (2005, Figure 1) show a drop in pressure in the ex-H ii region where the number of particles has suddenly been reduced by half due to instantaneous recombination. This will lead to the formation of spurious shells of material as shock waves propagate into the depressurized region from both sides.

The expansion of a stellar wind bubble into an H ii region cannot simply be modeled by a bubble expanding into a uniform medium at 10⁴ K (García-Segura et al. 1996a, 1996b; Ramirez-Ruiz et al. 2005), since the density of the photoionized gas in the H ii region becomes lower with time as the expansion proceeds. That is, not only is the stellar wind bubble expanding, but the H ii region is expanding as well.

In our 2D models we see the Rayleigh–Taylor and linear thin-shell instabilities reported by other authors (García-Segura et al. 1996a, 1996b; Freyer et al. 2003, 2006; Dwarkadas 2007) for the wind–wind interaction stage. The details very much depend on the stellar evolution model being used. Because we do not model the MS stage in 2D we do not see the H ii shadowing instability (Williams 1999) reported by Freyer et al. (2003, 2006) and Arthur & Hoare (2006).

We remark that the initial uniform ambient medium density of our models is n₀ = 100 cm⁻³, which is higher than that used by all the previous authors we have mentioned. The justification for our value is that a massive star will live a large fraction of its lifetime in the vicinity of the molecular cloud where it was born. As the H ii region around the star evolves, the density of the photoionized gas falls with time. For example, electron densities in the Orion nebula (age about one million years) range from 100 cm⁻³ (Felli et al. 1993) to >1000 cm⁻³ (Osterbrock & Flather 1959), while those in the Carina nebula (which harbors several high-mass MS stars and post-MS stars) vary between ∼30 cm⁻³ and ∼350 cm⁻³ (Mizutani et al. 2002). By the end of the MS, the density in our photoionized gas has fallen to <10 cm⁻³. The stellar wind bubble evolves within this changing environment, with the position of the inner wind shock regulated by the pressure in the H ii region. Models in which the external medium is initially less dense (e.g., Freyer et al. 2003, 2006; van Marle et al. 2005, 2007) produce larger, more diffuse H ii regions and consequently larger stellar wind bubbles than in our models. Simulations with a lower density medium and without radiative transfer produce large stellar wind bubbles with thin outer shells of swept-up material (Dwarkadas 2007) instead of broad shells of photoionized gas.

In Figure 15, we show the evolution of neutral and ionized gas kinetic and thermal energy fractions with respect to the total integrated stellar wind mechanical energy, $E_{w}(t) = \int \dot{M}_{w} V_{w}^2/2 \, dt$, resulting from the numerical simulations from the 40 M_☉ STARS models with and without thermal conduction. We do not include the energy of the stellar ionizing photons in this calculation because nearly all of the interstellar UV energy will be radiated away in forbidden-line processes (see Dyson & Williams 1997, chapter 7). At the beginning of the MS stage, the energy fraction profiles show some variations due to the formation and expansion of the initial H ii region around the star (Freyer et al. 2003; Arthur 2007).³ For most of the MS stage, the energy fractions of both kinetic and thermal energy are roughly constant (Freyer et al. 2003, 2006), with the thermal energy being lower in the model with conduction due to enhanced cooling in the hot bubble because the material evaporated from the dense outer shell increases the density, and hence the cooling rate, in the outer part of the bubble. It is evident that the outer neutral shell dominates the kinetic energy fraction of these models, while the hot, shocked bubble dominates the thermal energy.

Zoom In Zoom Out Reset image size

At the end of the MS stage, the neutral shell remains coasting along but the hot bubble suffers a drastic drop in pressure when the (subsonic with respect to the hot bubble) RSG wind starts. This can clearly be seen in Figure 15 where, starting at about 4 × 10⁶ yr, the thermal energy fraction dips while the kinetic energy fraction remains roughly constant. With the onset of the fast WR wind, the thermal energy increases enormously once the high-velocity gas shocks and thermalizes into a new pressure-driven hot bubble. In contrast, the kinetic energy becomes shared between the ionized and neutral components. This is due to photoionization of part of the outer neutral shell and the repressurized hot bubble rejuvenating the expansion of the outer swept-up material.

Between them, the thermal and kinetic energy fractions add up to substantially less than 1, particularly during the MS stage. The remaining energy is lost as radiation both from the hot bubble and from the H ii region and neutral shell which surround the wind bubble.

Our simulated WR wind bubbles bear morphological similarity to optical images of observed WR nebulae, such as RCW 58 and S 308 (Gruendl et al. 2000). Instabilities in the swept-up RSG material shell lead to the formation of dense, cold clumps, which can be compared to the Hα knots, for example, in images of RCW 58. Gruendl et al. (2000) divide a sample of WR bubbles into morphological types depending on the offset of the Hα and [O iii] emission. Type I has no measurable offset between the Hα and [O iii] fronts, Type II morphology has an Hα front trailing close behind an [O iii] front of similar shape, while for Type III the Hα trails far behind a faint [O iii] front and can differ appreciably in shape. Type IV morphology has a faint [O iii] front with no Hα counterpart. Under this classification system, RCW 58 is Type III, while S 308 is Type II. While the Hα is a recombination line produced in gas with temperatures ∼10⁴ K, [O iii] comes from hotter gas >10⁴ K due to either photoionization by harder UV photons or collisional excitation behind shock waves. In our simulations, the dense clumps that result from the instabilities which break up the swept-up shell are cold (T < 10⁴ K) and are photoionized on the side closest to the star. These clumps will therefore be sources of Hα emission and, moreover, lag behind the remainder of the shocked shell, which has temperatures T > 10⁴ K and could be a source of [O iii] emission. The images presented in Section 4.2 were chosen to facilitate comparison with the ring nebula S 308, which has an optical radius of ∼8 pc.

The initial abundances in the stellar evolution models from MM2003 and STARS are solar (Anders & Grevesse 1989) but toward the end of the RSG phase the stellar surface abundances begin to change, becoming enriched with nitrogen in particular. This surface material is expelled from the star in the form of stellar winds during the period of enhanced mass loss and can be expected to enrich the circumstellar nebula. Indeed, spectroscopic abundance determinations of WR ring nebulae (Esteban et al. 1992) show that the nitrogen-to-oxygen ratio, log (N/O), can even become positive in this CSM, while the solar value is log (N/O) = −0.9 (see Table 2).

In Figure 16, we plot the nitrogen-to-oxygen abundance ratio for stellar material expelled since the onset of enhanced mass loss for the different stellar evolution models used in this paper, using the model stellar surface abundances as being representative of the stellar wind chemical composition. These are mass-weighted average values over this time period and do not include material expelled during the MS stage. The MS stellar wind material is spread throughout an extended, diffuse bubble of radius up to some tens of parsecs, while the RSG and WR wind material will form a circumstellar nebula around the WR star. We compare our calculated abundance ratios to those determined for the three WR ring nebulae S 308, NGC 6888, and RCW 58 by Esteban et al. (1992) and listed in Table 2. It is clear from these figures that rotation significantly affects the abundances at the stellar surface, and hence in the stellar wind, throughout the MS stage so that by the start of the enhanced mass-loss period the nitrogen-to-oxygen ratio is already positive. For the evolution models without rotation, the abundances at the start of the RSG stage are still solar, since there has been no mixing with enriched material from the stellar core up to this point.

Zoom In Zoom Out Reset image size

From Figure 16, we see that the abundance ratio for RCW 58 can only come from material expelled in the early stages of enhanced mass loss. The S 308 and NGC 6888 values could be reproduced by various models, except that of the MM2003 40 M_☉ without rotation. For the STARS models, however, these abundances are achieved only once the WR stage has been entered and WR wind material has been added to the circumstellar nebula. For the MM2003 models with rotation, the nitrogen-to-oxygen abundance ratio in the CSM is predicted to be far higher than the observed values throughout the RSG and WR stages. These findings are broadly consistent with the classification of the central star of RCW 58 as late nitrogen-type WR (WNL), while the central stars of S 308 and NGC 6888 are early nitrogen-type WR (WNE).

There are around 200 WR stars in our Galaxy, and only two of them have nebulae detected in diffuse X-rays: S 308 and NGC 6888 (Bochkarev 1988; Wrigge et al. 1994; Wrigge 1999; Chu et al. 2003b; Zhekov & Park 2011). The X-ray emission is very soft and the diffuse emission from S 308 can be fit with a one-temperature component model with (T = 1.1 × 10⁶ K), enhanced nitrogen abundance, and X-ray luminosity L_X ≃ 1.2 × 10³⁴ erg s⁻¹ in the 0.25–1.5 keV energy band, assuming a distance of 1.8 kpc (Chu et al. 2003b). NGC 6888 ASCA and ROSAT X-ray observations are best fit by a two-component model with the dominant component at 1.3 × 10⁶ K and the weaker component at 8 × 10⁶ K and a total X-ray luminosity L_X ≃ 3 × 10³⁴ erg s⁻¹ in the energy range 0.4–2.4 keV, assuming a distance of 1.8 kpc (Wrigge et al. 1994, 1998, 2005). The recent Suzaku observations of NGC 6888 (Zhekov & Park 2011) confirm that 90% of the observed X-ray emission is in the soft 0.3–1.5 keV energy range. In both bubbles, the X-ray emission is internal to the optical [O iii] emission. Such soft X-ray emission indicates that the observations will be strongly affected by the absorbing column density, N_H, along the line of sight. This could be a reason for the non-detection of other wind-blown WR bubbles, such as RCW 58, since these objects are found in the plane of the Galaxy where N_H will be greatest. The two bubbles with detected X-ray emission are relatively close by whereas WR 40, the central star of RCW 58, is more than twice as distant (see Table 2).

5.2.1. X-Ray Simulations

From our numerical simulations, both in one and two dimensions, we can calculate the expected X-ray luminosity for any energy range at any time in the simulation. To do this efficiently, we use the temperature and density in each cell of the hydrodynamic simulation to construct the array of differential emission measures (DEMs), i.e., the emission measure (E.M. = ∫n_en_idV) that corresponds to each bin of a preassigned temperature array for a given numerical simulation. Note that we mask out the unshocked thermal wind in this procedure. The total X-ray spectrum is then found by summing the spectra of each temperature bin and the X-ray luminosity is then found by summing over the required energy range. This method is much more efficient than summing the spectrum (or luminosity) of each individual hydrodynamic cell, since we find that 100 temperature bins are perfectly adequate, while the numerical simulations consist of ∼6000 cells for the 1D models and >10⁵ cells for the 2D models. In Figure 17, we show the corresponding DEM histograms for the 40 M_☉ results presented in Figures 8–10 when the wind–wind interaction is at 8 pc, together with their counterparts that include thermal conduction.

Zoom In Zoom Out Reset image size

The first thing to notice about Figure 17 is that the most of the mass has temperatures T < 10⁶ K. This means that the spectra will be very soft and the absorbing column of neutral hydrogen toward such objects will play a very important role in what is observed. We see no systematic difference between models without thermal conduction and their equivalent models with conduction beyond the differences that can be attributed to slightly different volumes of emitting gas. The main difference we see is between models with complete, if distorted, shells (those of MM2003) and those where the shell has already broken up into clumps (those of STARS). The models with clumps have a much smaller DEM for the hot gas (T > 10^5.5 K) than those where the whole shell is being shocked, because the hot gas in these cases is more diffuse and can flow around the clumps and out of the computational domain.

Once we have the DEM arrays, the X-ray spectra are calculated assuming CIE and the ion fractions of Mazzotta et al. (1998). The line spectrum uses the APED database of CIE line intensities (Smith et al. 2001), while the continuum is calculated with our own code, taking into account free–free, free–bound and two-photon emission. The chemical abundances can be varied according to the requirements of the models.

Typical electron densities in the X-ray emitting gas for our simulations are n_e ∼ 0.1 cm⁻³ and temperatures are in the range (1–4) × 10⁶ K (see Figure 8 through 14 and 17). For these parameters, we can assess the closeness of the different elements (C, N, O, Ne, etc.) to ionization equilibrium using Figure 1 of Smith & Hughes (2010). At 10⁶ K, for nitrogen (the main X-ray spectral characteristic of S 308, Chu et al. 2003b) to be in ionization equilibrium requires n_et > 8 × 10¹¹ cm⁻³s⁻¹, whereas at 4 × 10⁶ K, ionization equilibrium for nitrogen requires n_et > 1 × 10¹¹ cm⁻³s⁻¹. For our typical electron density, these correspond to timescales of 253,000 yr and 32,000 yr, respectively. Thus, it could be argued that a non-equilibrium ionization treatment is required for the X-rays in the shocked stellar wind interaction region, but this is beyond the scope of this paper.

In Figures 18 and 19, we show the X-ray luminosity as a function of time in the WR phase for the 40 M_☉ models, assuming solar abundances (Anders & Grevesse 1989). Figure 18 plots the soft-band (0.25–1.5 keV) and hard-band (1.5–7 keV) X-ray emission for the 1D STARS models with and without thermal conduction, together with the STARS 2D model emission with thermal conduction. From this figure we see that the 1D models produce more intense X-ray emission, by an order of magnitude, than the 2D model. This can be understood by referring to Figure 8 where we see that the RSG shell breaks up into clumps, and the low-density, hot, shocked WR wind flows around them. The X-ray emission starts to diminish after ∼ 70,000 yr because the hot gas flows off the 10 pc 2D grid.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the 1D models, the WR wind is confined by the shell of RSG material, which cannot break up because of the geometrical restriction. The soft X-ray emission of the two 1D models is very similar until 1.5 × 10⁵ yr, when the thermal conduction at the edge of the hot bubble starts to play a dominant role, raising the density and lowering the temperature of the hot gas here and thereby increasing the contribution of the soft X-rays. In the 1D model without conduction, the gradual increase in the stellar wind velocity raises the temperature of the hot, shocked gas and finally results in the hard X-rays being the dominant component. Note that, for the 1D models, we are summing the X-ray contribution from the whole computational grid (except the thermal wind injection region): the contributions from the diffuse, hot, shocked MS wind to the soft and hard X-rays are, respectively, 2.1 × 10³¹ erg s⁻¹ and 2.9 × 10³¹ erg s⁻¹ for the model without conduction and 1.6 × 10³¹ erg s⁻¹ and 2.1 × 10³¹ erg s⁻¹ for the model with thermal conduction (see Figure 6).

Figure 19 shows the complete time evolution of the soft X-ray emission from all of the 40 M_☉ 2D models, assuming solar abundances: MM2003 with and without stellar rotation, and STARS, both with and without thermal conduction. Unfortunately, in all cases, the hot gas moves off the 2D computational grid (radius 10 pc) after a few tens of thousands of years, and this is what causes the reduction in the soft X-ray luminosity. However, we comment that this is a larger size scale than either of the X-ray observed WR wind bubbles. Our results are very varied. The STARS model, in which, as commented above, the shell of RSG material rapidly breaks up into clumps, has the lowest soft X-ray luminosity by an order of magnitude. In the MM2003 models, although the swept-up shells suffer instabilities, they do not break up into clumps before they move off the grid. The differences between the models with and without thermal conduction should be attributed mainly to differences in the position (and volume) of the wind interaction region due to the effect of thermal conduction on the pressure in the MS hot bubble, and only partly due to the local effect of thermal conduction on the interaction region between the hot WR wind and the swept-up RSG material (see Figure 11). Only three of our models attain X-ray luminosities above 10³³ erg s⁻¹: the MM2003 model without stellar rotation for the case without conduction, and the MM2003 models with stellar rotation both with and without conduction. Recall that the estimated total 0.2–1.5 keV band luminosity of S 308 is 1.2 × 10³⁴ erg s⁻¹ (see Table 2) from analysis of the NW quadrant of this object (Chu et al. 2003b), although a recent study of the full set of observations of S 308 suggests that the total luminosity is actually ∼5 × 10³³ erg s⁻¹ (see our companion paper J. A. Toalá et al. 2011, in preparation).

The stellar evolution tracks shown in Figure 4 indicate that none of the MM2003 models, 40 or 60 M_☉, with or without stellar rotation actually become RSGs, since the minimum stellar surface temperature attained by these models is above 4800 K. Rather, they become YSGs. The STARS models depicted in Figure 4 all become RSGs. In Table 3, we list the duration and mass lost in the yellow and RSG phases of evolution, separating the YSGs into those heading to the red and those heading back to the blue. We use the temperature range 4800 K < T_eff < 7500 K to define the YSGs (Drout et al. 2009).

The YSG stage is a particularly interesting stage of stellar evolution, since statistical surveys of nearby galaxies suggest that it should be short lived (Drout et al. 2009). For the STARS stellar evolution models, this is indeed the case, with the YSG phase lasting at most a few thousand years. Furthermore, the most important mass loss occurs during the RSG phase, when the stellar wind velocities are lowest (<30 km s⁻¹, although mass loss in the hotter, bluer phase becomes more important as the initial stellar mass increases). For the MM2003 models, on the other hand, there is no RSG stage, and the YSG stage lasts tens, or even hundreds, of thousands of years. For the MM2003 40 M_☉ model without rotation, the most important mass loss occurs as a YSG heading blueward. However, for the other MM2003 models, the most important mass loss does not even occur in the YSG phase, but in a hotter, bluer evolutionary state.

The different periods of mass loss have consequences for the CSM around the massive star. Intense mass loss in a short-lived RSG phase will result in a slow-moving dense shell of material close to the star; mass loss spread out over hundreds of thousands of years in a YSG phase will lead to an extended region of expanding CSM. Mass lost predominantly in hotter, bluer evolutionary stages may form fast-moving shells of material. When the fast WR wind interacts with the CSM, different structures will be produced depending on where the bulk of the circumstellar mass is and how fast it is moving (García-Segura et al. 1996a, 1996b).

From our results, we see that individual dense clumps are formed when the fast wind interacts with dense, slow-moving material close to the star, such as is the case for the STARS 40 and 60 M_☉ models (see Figures 8 and 12). This is because the linear thin-shell instabilities have time to break up the shell before it has traveled too far from the star. The clump densities and temperatures are such that the clumps will be photoionized and we would expect Hα emission in this case. When the CSM is moving more quickly, then the main fast-wind–slow-wind interaction occurs further from the star and, due to the lower densities (because of the geometrical dilution), the instabilities take longer to break up the swept-up shell. This is what we see in our MM2003 models (Figures 9, 10, 13). In this case, the densities of the clumps and filaments are lower and their temperatures are higher, so we would expect neutral material only at the beginning of the interaction, when the shell is densest.

In Table 2, we present a summary of the observed and derived properties of three WR bubbles S 308, NGC 6888, and RCW 58 and their central stars WR 6, WR 136, and WR 40. In this section, we compare these properties and other details with the results of our numerical simulations.

6.2.1. S 308

The wind-blown bubble S 308 around the WN4 star WR 6 is surprisingly spherical in [O iii] images, apart from a "blowout" region to the northwest. The Hα emission from this object shows a filamentary structure, and Gruendl et al. (2000) classified it as a Type II bubble. The optical radius of S 308 is 9 ± 2 pc at a distance of 1.5 ± 0.3 kpc (Chu et al. 1983, 2003b). Our figures in Section 4.2 were chosen such that the main features of the swept-up shell are roughly at this distance. It is important to point out that we are not trying to fit the wind parameters from the central star to the observed values, we are only choosing hydrodynamic results that approximately reproduce the morphology using the wind parameters obtained from the most appropriate of our stellar models.

The stellar evolution tracks shown in Figure 4 suggest that the central star, WR 6, would be well fit by any one of the MM2003 40 M_☉ models with or without stellar rotation, or by a STARS 45 M_☉ model (Hamann et al. 2006). Our numerical simulations show that the STARS 45 M_☉ model (not shown, but very similar morphologically to the 40 M_☉ model) would produce a Type III bubble (Hα clumps well separated from an [O iii] shell). On the other hand, the MM2003 40 M_☉ model without rotation would appear to produce structures where dense, photoionized gas exists just inside a warmer, collisionally excited (i.e., shocked) region, which could be interpreted as a Type II bubble. Model MM2003 with rotation breaks up the shell too far from the star and at the point depicted in Figure 10 (when the wind–wind interaction is at 8 pc), the Hα emission and [O iii] emission would be coincident.

Spectroscopic analysis of S 308 shows the nebula to have enhanced nitrogen abundances (Esteban et al. 1992). Our simple estimation of the evolution of the nebular abundances with time (see Figure 16), made by time averaging the total ejected mass of the different elements in the post-MS phase, where the relative abundances are given by those at the stellar surface at a given time, suggests that both the STARS model and the MM2003 40 M_☉ model with rotation could reproduce the observed abundances. However, the MM2003 model without rotation would not achieve the nebular abundances until the very end of its evolution.

Another test of the models is to compare the predicted X-ray luminosity and spectrum with observations. As can be seen from Figure 17, our simulations produce large amounts of gas with temperatures T ⩽ 10⁶ K, whereas observations are generally fit by a single temperature component, which is used to estimate the total unabsorbed luminosity. For S 308, Chu et al. (2003b) estimate a total X-ray luminosity in the 0.25–1.5 keV band of L_X = (1.2 ± 0.5) × 10³⁴ erg s⁻¹ by extrapolating their single temperature component (T ∼ 1.1 × 10⁶ K) results from the northwest quadrant assuming a spherical bubble. Chu et al. (2003b) find that the spectrum is well fit using enhanced nitrogen abundances, as determined by Esteban et al. (1992) for an absorbing column density of neutral hydrogen of N_H = 1.1 × 10²¹ cm⁻² (see Table 2). From Figure 19, we see that the STARS 40 M_☉ model completely fails to reproduce the observed X-ray luminosity, while the MM2003 model without rotation achieves this luminosity when thermal conduction is not included, and the MM2003 model with rotation both with and without thermal conduction produces sufficient X-ray luminosity in the soft energy band. We remark that the results in Figure 19 were calculated for solar abundances and that if the enhanced abundances are used, the luminosities are slightly different, depending on the dominant temperature. For example, the soft-band X-ray luminosity corresponding to the 40 M_☉ MM2003 model with rotation and thermal conduction, whose DEM is shown in the lower panel of Figure 17, is 8.2 × 10³³ erg s⁻¹ for solar abundances, but 4.8 × 10³³ erg s⁻¹ when we use the same abundances as Chu et al. (2003b).

We also calculate the predicted X-ray spectra of our models. We note that none of our models have been particularly tailored to S 308 (for example, we have not tried to match the stellar wind velocity to that of the central star WR 6). The main difference between our models and the reported observations is that the simulations contain a range of temperatures that contribute more or less equally to the total spectrum, whereas Chu et al. (2003b) find that the observed spectrum can be well fit by a single temperature component at 1.1 × 10⁶ K. In Figure 20, we show the simulated spectra for the 40 M_☉ MM2003 model without stellar rotation and with thermal conduction, where we have considered solar abundances and also the modified, nitrogen-rich abundances used by Chu et al. (2003b). The assumed distance is 1.5 kpc and the absorption column density is N_H = 1.1 × 10²¹ cm⁻², and we also convolved our model spectra with the instrumental response files appropriate for EPIC/pn on XMM-Newton. From this figure we see that the abundances make a large difference to the appearance of the spectrum. This is principally due to the reduced abundance of iron, which is responsible for the lines at around 0.72 keV. The model spectra are very similar to a single temperature model with T ∼ 3.6 × 10⁶ K, which is hotter than that derived for S 308. It is possible that non-equilibrium ionization effects are important in this nebula. The DEM for this simulation is shown in Figure 17 and we see that most of the gas actually has temperatures below 10⁶ K but the emission from this gas is completely absorbed by neutral hydrogen along the line of sight.

Zoom In Zoom Out Reset image size

6.2.2. NGC 6888

6.2.3. RCW 58

In spherical symmetry, the hydrodynamic equations in conservation form are

$$\begin{equation} \frac{\partial \rho }{\partial t} + \frac{1}{r^2}\frac{\partial r^2 \rho u}{\partial r} = q_{m}, \end{equation} \tag{ A1 }$$

$$\begin{equation} \frac{\partial \rho u}{\partial t} + \frac{1}{r^2}\frac{\partial r^2 (p + \rho u^2)}{\partial r} = \frac{2p}{r}, \end{equation} \tag{ A2 }$$

$$\begin{equation} \frac{\partial e}{\partial t} + \frac{1}{r^2}\frac{\partial r^2 u(e+p)}{\partial r} = G - L + q_{e}, \end{equation} \tag{ A3 }$$

where ρ, u, p are the mass density, radial velocity, and pressure, respectively, e is the total energy, defined by

$$\begin{equation} e = \frac{p}{\gamma - 1} + \frac{1}{2}\rho u^2, \end{equation} \tag{ A4 }$$

G and L are the heating and cooling rates, and γ is the ratio of specific heats. The stellar wind is injected as a thermal wind, with mass-loss rate $q_{m} = 3\dot{M}_{w}/4\pi R_{w}^3$, and thermal energy q_e = 3L_w/4πR³_w, where $\dot{M}_{w}$ is the stellar wind mass-loss rate, $L_{w} = \dot{M}_{w}{V_{w}}^2/2$ is the stellar wind mechanical luminosity, V_w is the stellar wind velocity, and R_w is the radius of the stellar wind injection volume (Chevalier & Clegg 1985). Outside of the stellar wind injection region, i.e., r > R_w, we have q_m = q_e = 0.

In our spherically symmetric simulations, we choose R_w = 30dr, where dr is the grid cell size. In all our 1D simulations, the cell size is set to be dr = 0.002 pc. The stellar wind injection region is small enough such that the stellar wind reaches its terminal velocity well before it shocks. This method of injecting the stellar wind is extremely robust to variations in the stellar wind parameters and does not result in spurious shock waves propagating back toward the center of symmetry.

In addition to the usual hydrodynamic equations, we use advection equations for the neutral and ionized density, which couple the hydrodynamics to the radiative transfer, for example

$$\begin{equation} \frac{\partial \rho n_{X}}{\partial t} + \frac{1}{r^2}\frac{\partial }{\partial r}(r^2 \rho u n_{X}) = 0 \, \end{equation} \tag{ A5 }$$

where n_X is the ionized or neutral number density. The ionized hydrogen fraction is updated using the current photoionization, collisional ionization, and recombination rates in the computational cell:

$$\begin{equation} \frac{\partial y_{i}}{\partial t} = n_e C_{{i-1}} y_{n} + P_\mathrm{i-1} y_{n} - n_e R_{i} y_{i} \, \end{equation} \tag{ A6 }$$

where y_n and y_i are the neutral and ionized hydrogen fractions, n_e is the electron density, C_{i − 1} is the collisional ionization rate of neutral hydrogen and R_i is the radiative recombination rate of ionized hydrogen. Both C_{i − 1} and R_i are functions of temperature only and we use analytic fits for both rates. The photoionization rate in each cell, P_{i − 1}, is obtained from the radiative transfer procedure.

In cylindrical symmetry, the conservation equations in the r − z plane are

$$\begin{equation} \frac{\partial \rho }{\partial t} + \frac{1}{r}\frac{\partial r \rho u_{r}}{\partial r} + \frac{\partial \rho u_{z}}{\partial z}= q_{m}, \end{equation} \tag{ A7 }$$

$$\begin{equation} \frac{\partial \rho u_{r}}{\partial t} + \frac{1}{r}\frac{\partial r (p + \rho u_{r}^2)}{\partial r} + \frac{\partial \rho u_{r} u_{z}}{\partial z} = \frac{p}{r}, \end{equation} \tag{ A8 }$$

$$\begin{equation} \frac{\partial \rho u_{z}}{\partial t} + \frac{1}{r}\frac{\partial r \rho u_{r} u_{z}}{\partial r} + \frac{\partial (p + \rho u_{z}^2)}{\partial z} = 0, \end{equation} \tag{ A9 }$$

$$\begin{equation} \frac{\partial e}{\partial t} + \frac{1}{r}\frac{\partial r u_{r}(e+p)}{\partial r} +\frac{\partial u_{z}(e+p)}{\partial z} = G - L + q_{e} \end{equation} \tag{ A10 }$$

where u_r, u_z are the radial and azimuthal velocities, respectively, p is the pressure, and e is the total energy, defined by

$$\begin{equation} e = \frac{p}{\gamma - 1} + \frac{1}{2}\rho (u_{r}^2 + u_{z}^2). \end{equation} \tag{ A11 }$$

The r and z directions are dealt with using operator splitting.

The energy equation includes radiative cooling and heating due to the absorption of stellar radiation. In photoionized gas, the radiative cooling term should take into account the fact that hydrogen is fully ionized even when the gas temperature is only 10⁴ K, and so collisional ionization of hydrogen is not an important cooling process in an H ii region (see Figure 21). We generate a cooling curve using the Cloudy photoionization code (Ferland et al. 1998) tailored to the ionizing source, which we tabulate and use as a look-up table during the calculation. Cloudy uses the most up-to-date and complete set of atomic data in the astrophysical literature and is constantly updated.

Zoom In Zoom Out Reset image size

The heating term at a given point in the numerical grid depends on the photoionization rate at that point and the stellar effective temperature, where the photoionization rate is calculated using the radiative transfer procedure.

The heating and cooling are source terms in the energy equation, which is solved together with the hydrodynamic equations. If thermal conduction is included in the simulations, this is treated as an update to the energy equation using operator splitting after the hydrodynamic equations have been solved.

In this work, we use the method of short characteristics, described in detail by Raga et al. (1999). This method consists of finding the column density from a source point to any point on the grid by first calculating the column density at all intermediate points, beginning with those closest to the source. The column density at point P is therefore given in terms of the column density at a neighboring point C by the equation

$$\begin{equation} N_{P} = N_{C} + l n_{0,C}, \end{equation} \tag{ A12 }$$

where l is the path length from C to P and n_{0, C} is the neutral hydrogen density at point C. This method can be used both for the direct radiation from point sources and also for the diffuse radiation from a distributed source. In the spherically symmetric case, for direct radiation from a point source, the procedure is trivial. For cylindrical symmetry, the calculation of the path lengths l is slightly more involved, since it is necessary to find the points where the characteristics cross the cell boundaries by linear interpolation.

The photoionization rates are needed to calculate the ionization state and heating rates of the gas. The column densities are therefore calculated at the beginning of every time step and these are used to find the optical depths to each grid cell center and hence the photoionization rates at those positions. The update to the ionized (and neutral) hydrogen fractions is done by operating splitting after the solution of the hydrodynamic equations.