Circumstellar medium around rotating massive stars at solar metallicity
Key Words.:ISM: bubbles – ISM: evolution – ISM: kinematics and dynamics – Stars: circumstellar matter – Stars: mass-loss
11email: email@example.com 22institutetext: École Normale Supérieure, Lyon, CRAL, UMR CNRS 5574, Université de Lyon, France 33institutetext: Ioffe Physical Technical Institute of the Russian Academy of Sciences, Saint Petersburg, Russia 44institutetext: Laboratoire Univers et Particules de Montpellier, Université Montpellier 2, CNRS/IN2P3, CC 72,
Place Eugène Bataillon, F- 34095 Montpellier Cedex 5, France 55institutetext: CSCS - Swiss National Supercomputing Centre, Via Trevano 131, 6900 Lugano, Switzerland
Aims:Observations show nebulae around some massive stars but not around others. If observed, their chemical composition is far from homogeneous. Our goal is to put these observational features into the context of the evolution of massive stars and their circumstellar medium (CSM) and, more generally, to quantify the role of massive stars for the chemical and dynamical evolution of the ISM.
Methods:Using the A-MAZE code, we perform 2d-axisymmetric hydrodynamical simulations of the evolution of the CSM, shaped by stellar winds, for a whole grid of massive stellar models from to and following the stellar evolution from the zero-age main-sequence to the time of supernova explosion. In addition to the usual quantities, we also follow five chemical species: H, He, C, N, and O.
Results:We show how various quantities evolve as a function of time: size of the bubble, position of the wind termination shock, chemical composition of the bubble, etc. The chemical composition of the bubble changes considerably compared to the initial composition, particularly during the red-supergiant (RSG) and Wolf-Rayet (WR) phases. In some extreme cases, the inner region of the bubble can be completely depleted in hydrogen and nitrogen, and is mainly composed of carbon, helium, and oxygen. We argue why the bubble typically expands at a lower rate than predicted by self-similarity theory. In particular, the size of the bubble is very sensitive to the density of the ISM, decreasing by a factor of for each additional dex in ISM density. The bubble size also decreases with the metallicity of the central star, because low-metallicity stars have weaker winds. Our models qualitatively fit the observations of WR ejecta nebulae.
Massive stars are the main drivers of the chemical enrichment of and kinetic energy deposition in the interstellar medium (ISM), first through their winds and then during the explosion of the supernova (SN). Since the pioneering analytical and spherically symmetric work of Weaver et al. (1977, 1978), the development of even more efficient and refined hydrodynamical codes has allowed better and better understanding of the evolution of the circumstellar medium (CSM) surrounding massive stars. Garcia-Segura et al. (1996a, b) performed one of the first multi-D simulation of massive star environments, following different mass-loss events: main-sequence (MS) wind, red-supergiant (RSG) wind, Wolf-Rayet (WR) wind, etc. More recently, Freyer et al. (2003, 2006) have studied the energetic impact of stellar winds on the ISM in more detail. The effect of stellar rotation has been discussed in various works (Chita et al., 2007, 2008; van Marle et al., 2008), as has the inclusion of ionisation in the computations (Toalá & Arthur, 2011; Dwarkadas & Rosenberg, 2013). The chemical enrichment of the bubble created by the winds has also been studied by Kröger et al. (2006). More extensive studies of the evolution of the CSM around whole grids of stellar models have been done by Eldridge et al. (2006).
Recently, the Geneva stellar evolution group has released a new model grid at solar metallicity (), with their most recent code (Ekström et al., 2012), including updated opacities and nuclear reaction rates, as well as the effect of internal rotation. These models were calibrated to reproduce a variety of observed characteristics: the characteristics of the Sun at its present age, the width of the MS in the Hertzsprung-Russell diagram (HRD), the position of red (super-) giants in the HRD, and the observed chemical composition and averaged surface velocities of B-type dwarf stars. They also succeed in partially reproducing the observed population of WR stars (and subtypes) at solar metallicity (Georgy et al., 2012). The goals of this paper are to simulate wind blown bubbles for the whole set of rotating models of the Geneva grid and to discuss both common and specific characteristics. We address, in particular, the size, temperature, and density of the bubble, its chemical composition, and potential (non-) appearance as a nebula, all as a function of time and initial mass of the central star and with special attention to the effects of various successive wind episodes. We do not aim to perform a very-high resolution study of a peculiar case of stellar evolution, as is usually done (see recently, e.g., van Marle & Keppens, 2012), but want to cover a wide variety of different cases instead.
The paper is organised as follows. Section 2 briefly recalls some important trends concerning the evolution of the central stars and explains the extraction of the wind properties (velocity, density) from the stellar evolution models. Section 3 presents the hydrodynamical code used in this work, as well as the set-up of the simulations. Our results are presented and discussed in Section 4, conclusions are drawn in Section 5.
2 The input model
2.1 Stellar models
We use as input the most recent grid of rotating stellar models provided by the Geneva group (Ekström et al., 2012). We only consider the following massive star models (see Table 1): , , , , , , , and . Models of lower mass stars have radiative winds that are too weak to be of interest in this framework. The considered high-mass models were computed from the zero-age main sequence up to the end of central carbon burning and cover the main stages of stellar evolution: main sequence, RSG, or WR phase, if any. After carbon is exhausted in the centre of the star, the core evolution becomes so fast that the surface is insensitive to what happens in the core, and thus evolves no more.
Among the available data, the following ultimately enter our simulations of the CSM: the luminosity and effective temperature of the star, its mass-loss rate, and its surface chemical composition, all of them being provided at each time step of the stellar evolution computation. The mass-loss rates are determined according to several prescriptions: Vink et al. (2001), if applicable, or de Jager et al. (1988) otherwise for stars hotter than ; de Jager et al. (1988) for red supergiants with ; Crowther (2001) when ; Nugis & Lamers (2000) and Gräfener & Hamann (2008) for WR stars (for more details, see Ekström et al., 2012).
For the massive star models, we have the following trends for the evolution(Georgy et al., 2012):
The follows a “classical” evolution, spending its MS in the blue part of the HRD, then crossing it up to the red supergiant branch, and ending its life on the red side of the HRD.
The and models, after spending some time as a red supergiant, cross back the HRD up to , and enter a WR phase at the very end of their lives, as WNL stars.
The and never reach the red supergiant branch, with a minimal . Then, they become WR stars, being successively WNL, WNE, and finally WC stars, with very high .
Models more massive than already become WR during the MS, completely avoiding the region of the HRD with .
|For these models already reaching the WR phase during the MS, the MS mass loss includes the amount of mass lost before the WR phase.||–|
|For these models already reaching the WR phase during the MS, the MS mass loss includes the amount of mass lost before the WR phase.||–|
|For these models already reaching the WR phase during the MS, the MS mass loss includes the amount of mass lost before the WR phase.||–|
A summary of the amount of mass lost during each stage is given in Table 1.
2.2 Extraction of the wind parameters: velocity and density
In this section, we describe how we reconstruct the physical parameters of the stellar winds from the surface characteristics of the stellar models. First, the velocity of the wind far from the star is obtained using the relations by Kudritzki & Puls (2000):
with the escape velocity at the stellar surface. To determine , we follow the procedure described in Lamers & Cassinelli (1996, see their eq. 16), and note that for hot massive stars it is crucial to account for the radiation pressure. The above formulation of is suitable for hot stars. For cool stars (), we set the velocity to an average value for red supergiant winds of .
The density of the wind at a distance (far enough) from the star of radius is given by
with the mass flux (per unit time and surface). This mass flux is computed as a function of the flux at the surface with . Even if the effect remains small for stars rotating far from their critical velocity, the anisotropy of the stellar wind is accounted for. The mass flux at the stellar surface is computed as in Georgy et al. (2011). For the initial velocity considered in this work, the maximal contrast between the polar and equatorial mass flux is . We thus have and ( being the colatitude) instead of a purely spherical symmetry. Upon mapping of the stellar wind parameters to radius angular momentum conservation of the wind is enforced.
2.3 Typical mass-loss histories
Figure 1 shows the time evolution of the terminal wind velocity and density at a distance from the central star, for the four different mass-loss histories discussed above. For the star, we see two distinct mass-loss episodes. The first one occurs during the MS with fast radiative winds (), but with low density . After the MS, the star becomes a RSG, with very slow winds (set here to ), but with a higher density of around .
The behaviour of the model is at first similar to that of the . However, after the RSG phase, the star evolves back towards the blue part of the HRD, becoming a WR star. During this last phase, the winds again speed up, while density remains high compared to its value during the MS.
At higher mass, the phase of a slow and dense wind progressively disappears. The star exhibits a very fast and relatively dense wind during the last few years before the final collapse. The different and time-dependent mass-loss histories have a strong impact on the power injected in the ISM by the winds. In Fig. 2, the power of the winds of our models is shown as a function of time. During the MS, the wind power is roughly constant for the models between and . It is followed by a more or less short period of low-power wind (corresponding to the phase when the star has left the MS, but is still not a WR star), before again entering a strong power wind phase (WR phase, except for the model, which avoids this phase). For the more massive models (), the power of the wind increases slightly during the MS (with a possible slight decrease at the very end of the MS), before becoming still stronger during the WR phase.
3 Details on the hydro-simulations
3.1 The A-MAZE numerical code and tools
The simulations were performed with the A-MAZE computational tool box (Walder & Folini, 2000; Folini et al., 2003), comprising massively parallel adaptive 3D MHD and radiative transfer codes, together with tools of data analysis and visualisation. This toolbox is validated well in that it has been used for a variety of different astrophysical flow and radiative transfer problems: colliding winds in binaries (Nussbaumer & Walder, 1993; Walder, 1998) and their spectral response (Folini & Walder, 2000), wind accretion (Walder et al., 2008), or supersonic turbulence (Folini & Walder, 2000, 2006). For this study, we have added some new tools to the box.
The first new tool is the extension to meshes on general curved manifolds for the case of ideal hydrodynamics (Euler equations). The grid is constructed on the basis of a Cartesian parameter space and a generic, not necessarily differentiable map from this parameter space to the physical space. The advantage of this procedure is that the adaptive mesh algorithm implemented in A-MAZE works also for a mesh on the curved manifold and the same Riemann-solvers as on a Cartesian mesh can be used on the manifold. The implementation closely follows the one described in Calhoun et al. (2008) and Colella (1990). This new tool of A-MAZE will be presented in detail in a forthcoming paper.
For this study, we used a specific map which results in a mesh identical to a spherical mesh with an equidistant angle discretisation and a logarithmic radial discretisation. With this mesh, the self-similar evolution phase of the bubble is well recovered, demonstrating the suitability of the chosen approach for the astrophysical problem under consideration (see Fig. 13 and the discussion in Sect. 4.6).
The second new tool is an elaborate, object-oriented implementation of the description of time-dependent stellar boundary conditions for entire stellar systems. For the present study, it allows us to accommodate in an easy and transparent way the time-dependent boundary conditions of our single star. Although not used in the present study, the implementation can equally handle boundary conditions for entire stellar systems, where each individual star has its own stellar evolution/time-dependent boundary conditions. Also possible here are, again in a transparent way for the user, changes of the outflow characteristics, e.g. stellar wind, jet, or supernova explosion. Finally, this new tool can account for proper motions of the stars relative to each other or relative to the computational grid.
From a technical code point of view, a stellar system consists of a given number of stars having different attributes that may be time-dependent: radius, shape, temperature, luminosity, rotation velocity, magnetic fields, etc. Different attributes of (time-dependent) interaction from the stars with their environment can be linked to each star. Among the interaction attributes implemented so far are isotropic and axisymmetric winds, accelerated or launched with the final velocity, jets, nova and supernova explosions. Stellar systems are implemented in three fortran90 modules: one, which defines a stellar system, a star with its attributes and its interactions, by means of fortran90 types; one, which treats the time-dependent boundaries for the different interaction types, and one that links these modules to the general structure of the code. Henceforth, the code is able to automatically change from one star type to another and from one interaction type to another. For instance, if a star explodes as supernova, the star type changes from an ordinary, evolved star to a neutron star or black hole, whereas the interaction type changes from wind to explosion and afterwards to accretion. The object-oriented implementation easily allows to extend the list of stellar or interaction attributes.
The stars may have fixed locations or may move on prescribed orbits or may self-consistently move under the influence of the gravitational forces between them. At the beginning of the simulation, a star track model is loaded. If one wants to follow the evolution of the stellar system after the first explosion of a star, one also loads the explosion characteristics (energy, asphericity, ejected yields, nature of compact remnant, etc.).
For this first study based on A-MAZE, the system only consists of one star at a fixed position, and the interaction type is always the same: aspherical winds shed with their final velocity. The radius, shape, temperature, luminosity, rotation speed, mass loss, and wind velocity however vary in time. The full power of this new A-MAZE tool can be grasp from Fig. 14 in the discussion section.
3.2 Initial and boundary conditions, numerical settings
The simulations presented in this paper were performed on a spherical grid, with a logarithmic mesh along the radial direction. The number of cells in the -direction is , and in the -direction it varies between 150 and 250, depending on the spatial extension of the physical domain. The radial extension of the domain lies between , where is chosen to be slightly smaller than times the minimal stellar radius during its evolution, and is chosen large enough to be more than the maximal spatial extension of the bubble. In the direction, the simulations cover the range , with the stellar equator being located at . The parameters are summarised in Table 2.
The winds are injected at a distance equal to times the current stellar radius if this value is less than , or else at , with a velocity and a density as discussed in Sect. 2.2. Boundary conditions are reflecting along the symmetry axis (i.e. for or ) and free flow at the outer boundary.
The interstellar medium is generally assumed to have a uniform density of , with some initial random perturbations (). Individual models with higher ISM density are discussed in Sect. 4.3. As the UV-flux of massive stars is very strong, the Strömgren radius around our stars is larger than the size of the created bubbles. We thus assume that the ISM is ionised and has a temperature of . If during its evolution the effective temperature of the star becomes lower than , the ISM and wind can cool down to a temperature of . The chemical composition of the ISM is the same as the initial composition of the stellar wind, i.e. solar: , , and , with solar mixture of metals (Asplund et al., 2005). We therefore do not consider the structure of the molecular cloud out of which the star is born. Also, all our simulations assume a single star. The effects of binarity or dense clusters will not be considered.
The cooling function used in this work is the one described in Wiersma et al. (2009). It allows small modifications in the chemical composition of the medium, which is adapted well to our problem, where the chemical composition of the winds evolves according to the stellar surface. However, the possible range of metallicities allowed by this cooling function is not sufficient to properly treat the most advanced stages of our WR stars. In this case, we use the cooling function with the minimum allowed hydrogen fraction.
4.1 Circumstellar medium around massive stars for different mass-loss histories
In this Section, we discuss the broad evolution of the CSM in the various mass-loss history described in Section 2.3. A more detailed discussion can be found in the next sections.
The main sequence phase:
This phase is common and qualitatively similar for all the models computed in the framework of this paper. During this phase, the quick stellar wind (characteristic of massive O-type stars) progressively pushes the ISM away. The size of the bubble increases progressively during the MS. Unlike other similar studies (Garcia-Segura et al., 1996b, a), where a strong over-density (compared to the density of the ISM) is observed at the outer edge of the bubble, we do not obtain such a structure in our simulations. We come back to this point in more detail in Sect. 4.4 and only stress here that our result is coherent with the recent study by Toalá & Arthur (2011), who accounted for the effects of ionisation and found an over-density only further away from the star at the location of the ionisation front. The typical aspect of the CSM at the end of the MS is shown in Fig. 3 (top-left panel). Near the star, the density of the stellar wind decreases with increasing radius, up to the wind-termination shock, where the density is minimal. In the bubble, the density is then roughly constant, up to the edge with the ISM.
The O-star – RSG scenario:
This evolutionary path is followed by models with an initial mass . After the MS phase, these stars enter the RSG phase, characterised by a slow and dense wind. This wind creates a dense shell around the star, which progressively extends through the bubble left after the MS (see Fig. 3, top-right and bottom-left panel). In the transition region from the slow RSG wind to the fast MS wind, the velocity difference of the two winds results in a rarefaction wave, possibly accompanied by a shock. Instabilities are generated from numerical seeds as these waves hit the previously existing wind termination shock (e.g. Walder & Folini, 1998). The associated slight density enhancements are further amplified and, to some degree, preserved by cooling effects. Rayleigh-Taylor instabilities of the RSG shell are also likely to be present (Dwarkadas, 2007). Details of the instabilities are washed out by the coarse resolution of the simulation, so what remains are the over-dense rays apparent in Fig. 3. Also at the onset of the slow wind, the pressure inside the bubble decreases slightly, causing the bubble to shrink. The high density eases the cooling, and the dense matter shell is relatively cold compared to the rest of the bubble, with temperatures of the order of . At the end of its life, the star is enshrouded in an asymmetrical and dense region, surrounded by the remnant of the bubble created during the MS (Fig. 4 left panel). The high-density structure seen in the interior of the bubble is the result of two large roll-ups that ultimately gain over other, higher order modes that are initially present. A full 3D simulation would certainly show a different structure.
The O-star – RSG – WN scenario:
In the framework of the new grids of stellar models from the Geneva group (Ekström et al., 2012), the models with an initial mass again cross the HRD after becoming a RSG as described above, and end their stellar lives as WN stars. In that case, after a short episode of slow wind, the stellar wind becomes much more rapid. As shown in Fig. 3 (bottom-right panel), the WR wind progressively pushes the dense shell produced during the RSG wind episode away. This shell progressively dilutes in the previously built bubble, increasing its density and softening the edge of the interface between the bubble and the ISM. Near the central star, a structure similar to the one existing during the MS appears, with a shock braking the wind of the WR star.
The O-star – YSG – WN – WC scenario:
This kind of evolution is typical of stars with initial mass in the range to . In that case, at the end of the MS, the star’s effective temperature decreases briefly down typically to (and thus completely avoiding the RSG branch). During this short phase (during which the star is a yellow supergiant, YSG), the wind density increases, and the wind velocity decreases, but not as much as in the RSG case. However, the effect is relatively similar, with the ejection of a dense and cold shell around the star, which is progressively expelled away by the upcoming very fast wind of the WR star. We do not observe instabilities leading to the formation of a thin and dense shell, as usually reported by other authors (for example, Garcia-Segura et al., 1996a). As mentioned in van Marle & Keppens (2012), this is certainly because our spatial resolution is not high enough to capture this behaviour. The final result is that the mean density of the bubble is slightly increased by the dilution of the slow wind and that the edge of the cavity is not as sharp as during the MS.
The O-star – WN – WC scenario:
For the most massive stars, the WR phase begins before the end of the MS. In that case, the star remains at very high during its entire evolution. Particularly, there is no ejection of a dense and slow wind before the onset of the WR wind. The evolution after the MS is thus in the continuity of the first time evolution, except that the bubble growth velocity increases slightly, because the WR wind deposits more momentum. The typical aspect of the bubble at the end of the star’s life is shown in Fig. 4 (right panel).
4.2 Size and chemical composition of the bubble
The left-hand panel of Fig. 5 shows the size of the bubble, defined as the region where the medium consists of at least stellar material and the position of the wind termination shock (right panel). Thanks to the increasing mass-loss rate and wind velocity with increasing stellar mass, the size of the bubble is larger for higher initial mass. The final size of the bubble spans from for the bubble around the model to more than around the model. Near the end of the tracks of the models, a decrease in the size of the bubble is apparent, which occurs during the RSG phase. During this phase, the pressure in the bubble decreases, causing it to shrink. For the models ending as a WR star, the growth of the bubble increases at the onset of the fast and dense WR wind.
The position of the wind termination shock follows the same trend. Because it depends strongly on the wind parameters (mass-loss rates and velocity), which can change slightly during the stellar evolution, it is not perfectly smooth during the MS. It disappears during the RSG phase for models having such a stage. During the WR phase, owing to the complexity of the CSM near the star, it is sometimes difficult to attribute a position to this shock. This explains the somewhat chaotic curves during this phase. For the most massive models, the position of this shock is mostly farther than from the star.
An interesting quantity is the ratio between the amount of mass lost by the star during its entire life, , and the amount of interstellar medium mass pushed away by the stellar wind during this time, . This ratio, shown in Fig. 6 for all models, is representative of the density contrast between the mean density in the bubble and the ISM density since most of the matter present in the bubble originates in the stellar winds. As can be seen, the mean density in the bubble is the lowest for the most massive models. The mean density in the bubble around the model is roughly one tenth of the initial ISM density. For the model, this number drops to less than . Although is highest for the most massive models, the increase in bubble size with stellar mass compensates for and even dominates it, explaining this contrast.
Owing to the action of mass loss and mixing, the chemical composition of the stellar wind changes progressively as a function of time (see Georgy et al., 2012, their Fig. 6). This produces a progressively inhomogeneous modification of the chemical composition of the CSM in the bubble. Observationally, several measurements show that most of the WR bubbles are enriched in He and N (Esteban et al., 1992), typically corresponding to matter processed by the CNO cycle. More recently, Fernández-Martín et al. (2012) have measured various abundances at different locations in the WR nebula NGC 6888 (the central star is a WN6 star) and found different compositions in different regions. They particularly found that the WR nebula can be divided into three zones, with different chemical compositions: an inner shell showing a strong nitrogen enrichment, an O/H ratio slightly deficient (compared to the solar value), and a strong overabundance of He. This shell is surrounded by an outer shell, where He/H abundance is higher than the solar value, but where both N/H and O/H do not seem to be enriched. Finally, there is a surrounding skin, with composition similar to the one expected for the ISM at this galactic radius.
In our simulations, the distribution of chemical species around the star during its life can vary a lot and lead to very inhomogeneous structures. In Fig. 7, we show the situation at the end of the MS (for a star, top panel) and immediately before entering the WR phase (for a star, bottom panel). The compositions shown are averages over constant radii (x-axis). Both cases are quite similar, the chemical composition in the bubble being almost constant (except near the edge of the bubble) and showing the products of CNO-cycle burning at low temperature (increase in He and N, decrease in H and C, O almost constant) for the model, and at high temperature for the (decrease in H, C, and O, increase in He and N).
For stars having a RSG phase after the MS, the situation evolves slowly during this phase with slow and dense winds. The development of a large external convective zone below the stellar surface means that more CNO-cycle elements are rapidly brought to the surface and carried away by the winds. This is illustrated in Fig. 8, where we see that the chemical composition in the external zone of the bubble remains similar to the one at the end of the MS, while near the star, where the shell of dense matter produced by the RSG winds is located, we observe a medium enriched in He and N and depleted in H, C, and O compared to the external zone.
Above , the stars have a WR phase at the end of their lives. Figure 9 shows the abundance profile in the bubble for the and models. The first one evolves into a WR star after a RSG stage, while the second already becomes a WR star during the MS, avoiding the RSG stage altogether. In the case of the star, we see that the abundance profile at the end of the WR phase has changed since the end of the RSG stage (see the bottom panel of Fig. 8). The stellar winds have progressively filled the bubble with material strongly modified by the CNO-cycle burning: a strong H depletion, with He becoming the dominant element in a large fraction of the bubble. In the meantime, C and O abundances decrease, and are replaced by N. A complete view of the evolution of chemical species in the bubble of the model is shown on the left-hand panel of Fig. 10.
The case of the is even more extreme. The very strong mass loss enriches the wind not only in the burning products of hydrogen, but also in the products of He-burning (mostly C and O) at the end of the evolution. This creates two distinct regions in the bubble. In the external region (between and ), we find the H-burning products. Particularly, the He and N abundances are maximal; however, this area is also polluted in C and O from the internal region, increasing the abundance of these elements. In the inner region, C becomes the dominant species. Oxygen is also strongly enhanced, while He is depleted. Finally, H and N are completely absent in this region. The time evolution as a function of the radius for each element illustrating these various effects is shown in the right-hand panel of Fig. 10. Looking at both panels of the same figure, we also see that the stellar winds are efficient to quickly modify the chemical composition inside the bubble, the abundances variation following the chemical composition of the stellar surface quite closely.
The case of NGC 6888 and its central star WR 136 (Johnson & Hogg, 1965) qualitatively fits our results. According to the stellar parameters given in Hamann et al. (2006) and to the HRD tracks by Ekström et al. (2012), the central WR star (a WN6 star) corresponds to a star with initial mass around . During the WN phase of this model, the abundance profile is similar to the one (top panel of Fig. 9), with a central region strongly modified by the products of the H burning (He and N increased, H, C, and O depleted). Near the edge of the bubble, we find a transition region towards the composition of the ISM.
4.3 ISM density and metallicity effect
The final radius of the bubble, hence the size of the region where the stellar winds substantially modify the chemical composition of the medium, mainly relies on the assumed density of the ISM and on the efficiency of the stellar winds at injecting momentum into the circumstellar medium. This last point itself depends on the metallicity of the central star, since the radiative mass-loss rate depends on the metal content of the stellar atmosphere (, Vink et al., 2001)222Changing the metallicity of the star not only changes the total amount of mass lost, but also affects substantially the evolutionary tracks in the HRD..
To quantify these effects, we computed six additional models: four and models at , with (10 times our standard value) and (100 times the standard value), and two models at the metallicity of the Small Magellanic Cloud (SMC, ), one of and one of (Georgy et al., 2013, in press). The final bubble radii of these models are summarised in Table 3. Increasing the ISM density decreases the size of the bubble. The multiplying factor is dependent on the initial mass ( for the model, for the ). A more complete discussion can be found in Sect. 4.6. Decreasing the metallicity by a factor of 7 (between and ) decreases the size of the bubble by a factor of .
These results show the extreme dependence of this kind of simulation on the various initial and boundary conditions. Peculiar results obtained for one stellar model of a given mass and metallicity are hardly applicable to other situations. This demonstrates the value of grids of models as presented here.
Another interesting question in this context is what the evolution of the bubble in a very inhomogeneous medium would be, which is a typical situation for star-forming regions where most of the massive stars are found. The simulations presented in this work are valid for a single massive star, and do not take possible binarity or the presence of one or more other massive stars in the surrounding few 10 pc into account, as would be the case in a cluster, for example. In the latter case, the situation would be much more complex, with interactions between the stellar winds. Some exploratory work into this direction is shown in Fig. 14, where we see complex turbulent interaction regions of the individual stellar winds. Mixing of chemical species is probably more efficient in such interaction regions, as is the outward transport of matter and momentum: the interaction regions in Fig. 14 extend farther outward than do the (partial) bubbles of the individual stars (see also van Marle et al., 2012). The possibility that the star is moving with respect to the ISM also changes the behaviour of the CSM (see for example the recent work by Decin et al., 2012; Mackey et al., 2012).
4.4 Nebulae around O- and WR-stars
Unlike other studies (e.g. Garcia-Segura et al., 1996b, a) where a double-shock structure appears in the simulation, our computations show only one shock (the wind termination shock) and no shock at the outer edge of the bubble in most cases. In our simulations, such a structure only appears in the very early evolution and in some cases in the WR state. This is due to the expansion velocity of the bubble through the ISM in our simulations, which is significantly subsonic for most of our models, except the most massive ones (above -, see Fig. 11). An outer shock, able to significantly compress the ISM matter, thus can only develop for these very massive star models above approximately - or during the first few hundred thousand years of the star’s lifetime. For stars with a lower mass and subsonic expansion, no outer shock can form and no shock-compressed shell can establish. Instead, ISM matter piled up by the stellar wind is simply pushed away by a pressure wave at approximately the sound speed.
Given the critical role of the bubble’s expansion velocity with respect to the sound speed of the ISM, one may ask about the role of the temperature of the ISM, which determines its sound speed. The ISM temperature itself depends on the degree of ionisation. Taking ionisation effects into account in our simulations would be desirable, but is not feasible for as many models as we study here, owing to computational limitations. Looking at the literature, results from radiation-hydrodynamical models for single-star models that include ionisation effects show somewhat conflicting results: for a star, Toalá & Arthur (2011) find a completely ionised ISM and swept-up shell (), whereas for a similar star Dwarkadas & Rosenberg (2013) find that the swept up shell is only partly ionised (). However, even in this latter case, the temperature of the swept-up shell is always higher than .
Returning to our own study, we take from Fig. 11 that even for such a low temperature, the expansion velocity of the bubble is subsonic for MS O-type stars with a mass lower than about . Nebulosities based on a compressed swept-up shell around such stars thus seem unlikely to exist in the framework of the hypothesis assumed in our computations (fully ionised and homogeneous ISM, single star). By contrast, for MS O-type stars with masses above about , such nebulosities may exist. As mentioned above, taking ionisation into account could change this picture, particularly during the first stages. The simulation of the formation of HII regions in star-forming region requires physical input that is not included in our simulations (see e.g. Tremblin et al., 2012b, a). Most probably, more detailed simulations than presented in this work are required to explain the observed bubbles around young OB-stars by the GLIMPSE survey (Churchwell et al., 2006).
The case is different for WR stars. Several WR stars exhibit ejecta nebulae around them, which are observable in visible wavelength (H, Stock & Barlow, 2010). These ejecta nebulae are interpreted as resulting from the interaction of slow and dense matter ejected during a previous mass-loss episode, with the fast wind of the WR star. Most of the nebulae analysed by Stock & Barlow (2010) are found around WN type WR stars (10 out of 13 in the Milky-Way sample). This is coherent with the picture that emerges from our simulations, showing that the dense region around the central star created during the RSG phase is progressively diluted during the WN phase, thereby decreasing its density, hence its observability. Moreover, in the framework of single-star evolution, the stars that have a RSG phase and then a WR phase finish their lives as WN stars (Georgy et al., 2012). We thus expect most of the WR ejecta nebulae to be found around WN stars. The relatively small number of WN stars showing such a nebula (, according to Stock & Barlow, 2010) is also understandable because the nebula is only visible close to the WR star (and thus only at the beginning of the WR stage), and also because several WN stars originate in more massive progenitors that completely avoid the RSG phase (Georgy et al., 2012).
4.5 Impact of the grid resolution
The simulations presented in this paper have a lower resolution than other recent work (e.g. van Marle & Keppens, 2012; Krause et al., 2013). The choice is deliberate, for it allows us to compute the ISM time evolution for a whole grid of stellar models, an impossible endeavour to tackle at much higher resolution owing to CPU time limitations. The price to pay is a restricted ability to catch the development of small-scale instabilities. However, since our aim is to describe the broad chemical features and bubble sizes, this limitation is not a severe one as we demonstrate in the following.
Doubling the spatial resolution of our grid, we recomputed four simulations: the and models, since they are representative of the evolution of the stellar winds studied here, against ISM densities of and . The evolution of the was followed up to , and the up to . Results from the resolution study are compared in Fig. 12 for the low-density ISM case. The high-density ISM case behave similarly. As can be seen, the bubble is slightly larger in the better resolved simulation. The contact interface is steeper, and its inner edge is shifted outward by about 9%. The resolution effect is smaller on the outer shock. The steepening of the contact interface results because its width in terms of grid cells remains about constant, independent of numerical resolution (about 3 to 5 cells for the high-resolution integrator used here). Its smaller spatial extension is crucial since the interface is the main cooling region of the bubble. Here, the density increases rapidly from the low values in the bubble interior towards the high value found in the dense region of cold ISM material and, at the same time, the temperature decreases from the very high value of the shocked wind. The intermediate temperature and density values are the most favourable for radiative cooling. Thus, in summary, higher numerical resolution results in a thinner contact interface and reduced radiative losses, thus more energy remains to drive the bubble expansion and the bubble evolves slightly faster than in the simulation with lower numerical resolution (see also the blue and red curves in Fig. 13). This effect is well known and described, for instance, in Folini & Walder (1995), Walder & Folini (1996), or Parkin & Pittard (2010).
Even if the size of the bubble shows some sensitivity to the resolution, results are not too different with regard to latitudinal averages of the density, temperature, and chemical composition (see Fig. 12). We see that except for the small shift produced by the slightly different size of the bubble at any given time, the averaged quantities we consider in this work are remarkably similar.
From the above results and considerations we conclude that resolution effects are not very important for our main conclusions, which are targeted at large-scale features. This is not to say that much higher resolutions or 3D simulation, both potentially allowing for more and different instabilities, would not add to the picture arising from the present study.
4.6 Non-self-similar evolution
Weaver et al. (1977) have shown that a bubble evolves in a self-similar way under many circumstances. For instance, for a purely adiabatic evolution (polytrophic EOS with ), a constant wind power, , a constant density of the ISM , and a negligible pressure of the ISM, the position of the shock towards the ISM (index s) and the position of the contact discontinuity (index c) evolve as , where is a parameter of order unity. The velocities of these waves scale as . However, as also discussed by Weaver et al. (1977), the bubble does not always evolve in a self-similar way, in particular when radiative cooling affects its dynamics or when the self-similarity is broken by the presence of a substantial pressure in the ISM.
Figure 11 and Table 3 indicate that in the presented simulations a self-similar evolution is the exception. This is confirmed by a closer analysis, as summarised in Fig. 13. Indeed, this can be expected, since radiative cooling is included in our models, the wind power is not constant over the evolution (Fig. 2), and the pressure of the ISM (in particular in the models with a dense ISM) cannot be neglected. To understand the behaviour of our models in this context, we computed additional models: some with an increased resolution (factor of 2 in each direction, see Sect. 4.5), and some without any radiative cooling (adiabatic case).
In a first phase, present in all our simulations, the size of the bubble increases linearly, , approximately until . The duration of this linear phase depends somewhat on the stellar mass and ISM density. The reason for the existence of this linear phase is –similar to the case of a supernova remnant (see e.g. Chevalier, 1982)– that sufficient ISM material has to be collected in front of the snow plough before the interaction zone slows down and adopts self-similar behaviour.
In a next phase, we find clear differences between the adiabatic models (black in Fig. 13) and the models including radiative cooling and using different numerical resolutions (blue and red lines in Fig. 13). While some models show self-similar evolution of the form , others show a shallower dependence on time. A clear law is seen for the adiabatic model against . This model essentially meets the underlying assumptions (see above) of the self-similarity law. The self-similar phase lasts until the sharp rise in wind power close to the end of the star’s life (Fig. 2). Small deviations from the self-similar evolution we ascribe to the non-constant wind-power even during the MS evolution.
Reduction of the wind power tends to reduce the slope below . The adiabatic model with already deviates towards a power law coefficient that is clearly smaller than . The deviation starts when the pressure of the ISM becomes similar to the pressure of the shocked wind. Similarly, a denser ISM (, dashed lines in Fig. 13) also leads to an earlier and more pronounced deviation from . For the model, the slope is close to .
Radiative cooling also tends to reduce the slope of the power law evolution, as part of the wind energy is radiated and cannot be used to drive the bubble. While the radiatively cooling model against still has a slope close to , the model against has a shallower slope close to , and the model against has a much shallower slope of only about . The much more pronounced dependency on the ISM density, as compared to the adiabatic case, is expected since cooling is proportional to density squared. Nevertheless, the ISM density does not translate directly into slope changes, because the main cooling region of the bubble is the contact interface where temperature and density take intermediate values between the hot bubble and the cold ISM. This is detailed in Sect. 4.5 and in the literature (Folini & Walder, 1995; Walder & Folini, 1996; Parkin & Pittard, 2010), along with an explanation of why higher spatial resolution results in reduced radiative losses and a slightly faster bubble expansion, i.e. a slightly steeper slope for the high-resolution simulations, as is apparent in Fig. 13 (red lines are steeper than blue lines).
We conclude that (1) bubbles do not generally evolve self-similarly over their entire lifetimes and (2) that changes in the ISM density or the metallicity of the star cannot be taken into account by a simple analytical formula. Instead, the values given in Table 3 mirror these differences much better.
5 Discussion and conclusions
Returning once more to the complex flow shown in Fig 14, we briefly want to highlight the potential impact of such flow conditions for high-energy astrophysics. Future studies in this direction are envisaged, profiting from our enlarged toolbox, especially the elaborate implementation of stellar systems as described in Sect. 2. Powerful stellar winds of young massive stars and core-collapse supernova explosions with strong shock waves can convert a sizeable part of the kinetic energy release into fluctuating magnetic fields and relativistic particles. The star-forming regions and compact young star cluster are considered as favourable sites for energetic particle acceleration and could be seen as bright sources of non-thermal emission with the upcoming high-energy instruments, either at hard X-rays or at TeV gamma-rays by the Cerenkov Telescope Array (CTA).
Rich associations of OB-stars, particularly Cygnus OB, have recently been detected at gamma-rays. The Fermi Large Area Telescope has detected photon emission from a 50-parsec-wide cocoon-like structure in the active Cygnus X star-forming region (Ackermann et al., 2011). The authors proposed that the cocoon is filled with freshly accelerated cosmic rays that flood the cavities carved by the stellar winds from young stellar clusters. The gamma-ray luminosity estimated is about , which is a small fraction of the kinetic power of the stellar winds. However, to produce the observed gamma-ray luminosity, one should maintain much greater power in cosmic rays, and this requires an efficient conversion of the kinetic power of the stellar winds and supernovae. An extended source of very high energy emission was also detected with the High Energy Stereoscopic System (H.E.S.S.) from Cyg OB2 region (Aharonian et al., 2002) and a very massive young compact cluster, Westerlund I (Abramowski et al., 2012). Starburst galaxies NGC 253, M82, and some others (see for a recent review Ohm, 2012) have demonstrated the high-energy emission spectra that are harder than that of the Milky Way or M31, where the global diffuse TeV regime emission has not been reported so far.
Complex supersonic flows carrying magnetic fluctuations of a broad dynamical range of scales can efficiently interact with relativistic particles. The interaction results in particle acceleration, which in turn modifies the plasma flows and affects specific mechanisms of magnetic field amplification in the vicinity of collisionless shocks (van Marle et al., 2012). Relativistic particles are subject to Fermi acceleration at strong astrophysical collisionless shocks (see for a review Blandford & Eichler, 1987). Ensembles of multiple shocks accompanied with large-scale MHD motions are very efficient particle accelerators (e.g. Bykov, 2001; Bykov et al., 2013; Parizot et al., 2004; Ferrand & Marcowith, 2010), and they can be the sources of the galactic cosmic rays accelerated beyond the spectral knee (e.g. Bykov & Toptygin, 2001).
In this paper, we present hydrodynamical simulations of the time evolution of the CSM around massive stars ( to ), based on the recent grid at solar metallicity provided by Ekström et al. (2012). We show the differences between various evolutionary scenarios, by extracting from our data various averaged quantities as a function of (stellar evolution-) time and/or the distance to the central star: radius of the bubble, position of the wind termination shock, density, temperature, and chemical composition of the gas. The differences found, which cannot be overcome by simple scaling laws, demonstrate the need for grids of models as a complement to very detailed (and computationally expensive) studies of one particular evolutionary track.
For the range of models considered, we find bubble radii right before the SN explosion between roughly 10 and 100 pc. Much smaller bubbles are obtained for higher ISM densities (factor 3 for ten times higher density) and/or lower metallicities of the massive star. Average densities within the bubble are generally (much) lower than ISM densities, except very close to the star. The chemical composition in the bubble can be very inhomogeneous and very different from the ISM composition, owing to the progressive modification of the chemical composition of the star. This finding is in line with observations of the WR nebula NGC 6888.
Concerning, more generally, the (non-) existence of observable nebulae around massive stars, our models allow for WR ejection nebulae (WR wind against RSG wind) around young WN type WR stars, in line with observations. On the scale of the entire bubble, our models suggest that the wind blown bubble is confined by an (observable) outer shock wave only for the most massive stars and then only during their WR phase. During MS, a shock at the outer rim of the bubble forms only at the very early phases or for the most massive models. Later on, no such shock is present, because the bubble expands subsonically. A shock front may, however, still develop from the action of the photoionisation front (Toalá & Arthur, 2011) or if hydrogen can recombine and the gas can cool to very low temperatures (Dwarkadas & Rosenberg, 2013). A complete view of the formation of nebulosity around MS massive stars thus requires taking physical effects into account that are not yet included in our simulations.
Leaving the concept of single stars in favour of stellar clusters, wind collision zones very likely alter the above results: mixing of chemical species is likely more efficient, and observable signatures more pronounced.
At the end of their lives, the models presented in this work explode as a SN. The question of how these explosions disperse the chemical species expanding through the previously formed bubble is important in the matter of the chemical evolution of the galaxies. This point will be studied in a forthcoming paper.
Acknowledgements.The authors thank the anonymous referee for his pertinent comments that have significantly improved the quality of the present work. The authors thank the Pôle Scientifique de Modélisation Numérique in Lyon for the access to their computing facilities. They are also grateful to R. Hirschi, for allowing them to use the Shyne cluster to perform part of the simulations presented in this paper. CG acknowledges support from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 306901, as well as from the Swiss National Science Foundation. CG, RW, and DF acknowledge support from the Programme National de Physique Stellaire. AB was supported by Russian Academy of Sciences Presidium and OFN-17 programme.
- Abramowski et al. (2012) Abramowski, A., Acero, F., Aharonian, F., et al. 2012, A&A, 537, A114
- Ackermann et al. (2011) Ackermann, M., Ajello, M., Allafort, A., et al. 2011, Science, 334, 1103
- Aharonian et al. (2002) Aharonian, F., Akhperjanian, A., Beilicke, M., et al. 2002, A&A, 393, L37
- Asplund et al. (2005) Asplund, M., Grevesse, N., & Sauval, A. J. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 336, Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis, ed. T. G. Barnes, III & F. N. Bash (San Francisco: ASP), 25–+
- Blandford & Eichler (1987) Blandford, R. & Eichler, D. 1987, Phys. Rep, 154, 1
- Bykov (2001) Bykov, A. M. 2001, Space Sci. Rev., 99, 317
- Bykov et al. (2013) Bykov, A. M., Gladilin, P. E., & Osipov, S. M. 2013, MNRAS, 429, 2755
- Bykov & Toptygin (2001) Bykov, A. M. & Toptygin, I. N. 2001, Astronomy Letters, 27, 625
- Calhoun et al. (2008) Calhoun, D. A., Helzel, C., & Leveque, R. J. 2008, SIAM Review, 50, 723
- Chevalier (1982) Chevalier, R. A. 1982, ApJ, 258, 790
- Chita et al. (2008) Chita, S. M., Langer, N., van Marle, A. J., García-Segura, G., & Heger, A. 2008, A&A, 488, L37
- Chita et al. (2007) Chita, S. M., van Marle, A. J., Langer, N., & García-Segura, G. 2007, in Revista Mexicana de Astronomia y Astrofisica, vol. 27, Vol. 30, Revista Mexicana de Astronomia y Astrofisica Conference Series, 80–83
- Churchwell et al. (2006) Churchwell, E., Povich, M. S., Allen, D., et al. 2006, ApJ, 649, 759
- Colella (1990) Colella, P. 1990, Journal of Computational Physics, 87, 171
- Crowther (2001) Crowther, P. A. 2001, in Astrophysics and Space Science Library, Vol. 264, The Influence of Binaries on Stellar Population Studies, ed. D. Vanbeveren (Dordrecht: Kluwer Academic Publishers), 215
- de Jager et al. (1988) de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259
- Decin et al. (2012) Decin, L., Cox, N. L. J., Royer, P., et al. 2012, A&A, 548, A113
- Dwarkadas (2007) Dwarkadas, V. V. 2007, ApJ, 667, 226
- Dwarkadas & Rosenberg (2013) Dwarkadas, V. V. & Rosenberg, D. L. 2013, High Energy Density Physics, 9, 226
- Ekström et al. (2012) Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146
- Eldridge et al. (2006) Eldridge, J. J., Genet, F., Daigne, F., & Mochkovitch, R. 2006, MNRAS, 367, 186
- Esteban et al. (1992) Esteban, C., Vilchez, J. M., Smith, L. J., & Clegg, R. E. S. 1992, A&A, 259, 629
- Fernández-Martín et al. (2012) Fernández-Martín, A., Martín-Gordón, D., Vílchez, J. M., et al. 2012, A&A, 541, A119
- Ferrand & Marcowith (2010) Ferrand, G. & Marcowith, A. 2010, A&A, 510, A101
- Folini & Walder (1995) Folini, D. & Walder, R. 1995, in Annals of the Israel Physical Society, Vol. 11, Asymmetrical Planetary Nebulae, ed. A. Harpaz & N. Soker, 253
- Folini & Walder (2000) Folini, D. & Walder, R. 2000, Ap&SS, 274, 189
- Folini & Walder (2006) Folini, D. & Walder, R. 2006, A&A, 459, 1
- Folini et al. (2003) Folini, D., Walder, R., Psarros, M., & Desboeufs, A. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 288, Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, 433
- Freyer et al. (2003) Freyer, T., Hensler, G., & Yorke, H. W. 2003, ApJ, 594, 888
- Freyer et al. (2006) Freyer, T., Hensler, G., & Yorke, H. W. 2006, ApJ, 638, 262
- Garcia-Segura et al. (1996a) Garcia-Segura, G., Langer, N., & Mac Low, M.-M. 1996a, A&A, 316, 133
- Garcia-Segura et al. (1996b) Garcia-Segura, G., Mac Low, M.-M., & Langer, N. 1996b, A&A, 305, 229
- Georgy et al. (2013) Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, ArXiv e-prints, 1308.2914
- Georgy et al. (2012) Georgy, C., Ekström, S., Meynet, G., et al. 2012, A&A, 542, A29
- Georgy et al. (2011) Georgy, C., Meynet, G., & Maeder, A. 2011, A&A, 527, A52
- Gräfener & Hamann (2008) Gräfener, G. & Hamann, W.-R. 2008, A&A, 482, 945
- Hamann et al. (2006) Hamann, W.-R., Gräfener, G., & Liermann, A. 2006, A&A, 457, 1015
- Johnson & Hogg (1965) Johnson, H. M. & Hogg, D. E. 1965, ApJ, 142, 1033
- Krause et al. (2013) Krause, M., Fierlinger, K., Diehl, R., et al. 2013, A&A, 550, A49
- Kröger et al. (2006) Kröger, D., Hensler, G., & Freyer, T. 2006, A&A, 450, L5
- Kudritzki & Puls (2000) Kudritzki, R.-P. & Puls, J. 2000, ARA&A, 38, 613
- Lamers & Cassinelli (1996) Lamers, H. J. G. L. M. & Cassinelli, I. P. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 98, From Stars to Galaxies: the Impact of Stellar Physics on Galaxy Evolution, ed. C. Leitherer, U. Fritze-von-Alvensleben, & J. Huchra, 162–+
- Mackey et al. (2012) Mackey, J., Mohamed, S., Neilson, H. R., Langer, N., & Meyer, D. M.-A. 2012, ApJ, 751, L10
- Nugis & Lamers (2000) Nugis, T. & Lamers, H. J. G. L. M. 2000, A&A, 360, 227
- Nussbaumer & Walder (1993) Nussbaumer, H. & Walder, R. 1993, A&A, 278, 209
- Ohm (2012) Ohm, S. 2012, in American Institute of Physics Conference Series, Vol. 1505, American Institute of Physics Conference Series, ed. F. A. Aharonian, W. Hofmann, & F. M. Rieger, 64–71
- Parizot et al. (2004) Parizot, E., Marcowith, A., van der Swaluw, E., Bykov, A. M., & Tatischeff, V. 2004, A&A, 424, 747
- Parkin & Pittard (2010) Parkin, E. R. & Pittard, J. M. 2010, MNRAS, 406, 2373
- Stock & Barlow (2010) Stock, D. J. & Barlow, M. J. 2010, MNRAS, 409, 1429
- Toalá & Arthur (2011) Toalá, J. A. & Arthur, S. J. 2011, ApJ, 737, 100
- Tremblin et al. (2012a) Tremblin, P., Audit, E., Minier, V., Schmidt, W., & Schneider, N. 2012a, A&A, 546, A33
- Tremblin et al. (2012b) Tremblin, P., Audit, E., Minier, V., & Schneider, N. 2012b, A&A, 538, A31
- van Marle & Keppens (2012) van Marle, A. J. & Keppens, R. 2012, A&A, 547, A3
- van Marle et al. (2008) van Marle, A. J., Langer, N., Yoon, S.-C., & García-Segura, G. 2008, A&A, 478, 769
- van Marle et al. (2012) van Marle, A. J., Meliani, Z., & Marcowith, A. 2012, A&A, 541, L8
- Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574
- Walder (1998) Walder, R. 1998, Ap&SS, 260, 243
- Walder & Folini (1996) Walder, R. & Folini, D. 1996, A&A, 315, 265
- Walder & Folini (1998) Walder, R. & Folini, D. 1998, Ap&SS, 260, 215
- Walder & Folini (2000) Walder, R. & Folini, D. 2000, in Astronomical Society of the Pacific Conference Series, Vol. 204, Thermal and Ionization Aspects of Flows from Hot Stars, ed. H. Lamers & A. Sapar, 281
- Walder et al. (2008) Walder, R., Folini, D., & Shore, S. N. 2008, A&A, 484, L9
- Weaver et al. (1977) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377
- Weaver et al. (1978) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1978, ApJ, 220, 742
- Wiersma et al. (2009) Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99