Hydrodynamical models of cometary HII regions
Harrison Steggles, Melvin Hoare, Julian Pittard
TL;DR
This study models the evolution of cometary HII regions driven by massive stars, exploring how stellar mass and core density influence features like limb-brightening and cavity growth.
Contribution
It introduces hydrodynamical models of cometary HII regions with off-centre stars in power-law density cores, highlighting the effects of stellar mass on observable features.
Findings
Stars with mass ≥12 M_sun produce limb-brightening.
Cavity sizes grow and plateau as HII regions evolve.
Spectral energy distributions are similar to uniform density models.
Abstract
We have modelled the evolution of cometary HII regions produced by zero-age main-sequence stars of O and B spectral types, which are driving strong winds and are born off-centre from spherically symmetric cores with power-law () density slopes. A model parameter grid was produced that spans stellar mass, age and core density. Exploring this parameter space we investigated limb-brightening, a feature commonly seen in cometary HII regions. We found that stars with mass produce this feature. Our models have a cavity bounded by a contact discontinuity separating hot shocked wind and ionised ambient gas that is similar in size to the surrounding HII region. Due to early pressure confinement we did not see shocks outside of the contact discontinuity for stars with , but the cavities were found to continue to…
| [] | [] | [] | [] | [] | [] | |||
|---|---|---|---|---|---|---|---|---|
| [] | [] | [] | [] | [] | [] | [] | [] |
|---|---|---|---|---|---|---|---|
| Dist. | [] | |||||
|---|---|---|---|---|---|---|
| [] | [] | |||||
| Dist. | Age [] | |||||
|---|---|---|---|---|---|---|
| [] | [] | |||||
| Dist. | [] | |||||
|---|---|---|---|---|---|---|
| [] | [] | |||||
| Dist. | Age [] | |||||
|---|---|---|---|---|---|---|
| [] | [] | |||||
| [] | |||||
|---|---|---|---|---|---|
| [] | |||||
| Age [] | |||||
|---|---|---|---|---|---|
| [] | |||||
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Hydrodynamical models of cometary Hii regions
H. G. Steggles,1 M. G. Hoare1 and J. M. Pittard1
1School of Physics and Astronomy, University of Leeds, Woodhouse Lane, Leeds LS2 9JT, UK E-mail: [email protected]
(Accepted 2016 December 23. Received 2016 December 22; in original form 2016 August 23)
Abstract
We have modelled the evolution of cometary Hii regions produced by zero-age main-sequence stars of O and B spectral types, which are driving strong winds and are born off-centre from spherically symmetric cores with power-law () density slopes. A model parameter grid was produced that spans stellar mass, age and core density. Exploring this parameter space we investigated limb-brightening, a feature commonly seen in cometary Hii regions. We found that stars with mass 12\text{,}\mathrm{\mathrm{M_{\sun}}} produce this feature. Our models have a cavity bounded by a contact discontinuity separating hot shocked wind and ionised ambient gas that is similar in size to the surrounding Hii region. Due to early pressure confinement we did not see shocks outside of the contact discontinuity for stars with $M_{\star}\leq$40\text{\,}\mathrm{\mathrm{M_{\sun}}}, but the cavities were found to continue to grow. The cavity size in each model plateaus as the Hii region stagnates. The spectral energy distributions of our models are similar to those from identical stars evolving in uniform density fields. The turn-over frequency is slightly lower in our power-law models due to a higher proportion of low density gas covered by the Hii regions.
keywords:
hydrodynamics – radiative transfer – stars: evolution – stars: massive – stars: winds, outflows – Hii regions
††pubyear: 2016††volume: 466††pagerange: 4573–4591
1 Introduction
Massive stars (8\text{,}\mathrm{\mathrm{M_{\sun}}}$$) are an important component of the universe. They are the main sources of stellar feedback that influences the structure, evolution and chemical composition of a galaxy. High-mass stars are born in dense parts of the ISM, erupting with outflows and radiatively driven winds that contribute to the destruction of the natal cloud. Extreme ultra-violet radiation is emitted from high-mass stars, which ionises and heats the surrounding gas. This produces an over-pressure that drives a shock into the ambient medium, sweeping gas into a shell that ultimately disperses the parent molecular cloud and may also trigger further star formation (Hoare et al., 2007). The optical and near-infrared emission from such stars is not able to penetrate the natal clouds due to high column densities of dust. Radio free-free emission can, however, be seen right across the Galaxy (Churchwell, 2002). At wavelengths in the far-infrared, emission coming from warm dust re-emitting the stellar radiation can also be observed from ionised nebulae (Cesaroni et al., 2015).
How high-mass stars are formed is still poorly understood. They have short main-sequence lifetimes (1\text{\times}{10}^{8}\text{,}\mathrm{yr}) which, along with their relatively low formation rates, indicates that only a small proportion ($\simeq$1\text{\times}{10}^{-6}) of stars in the galaxy have high masses. The scarcity of observational data is a significant problem for the testing of theoretical models; some important phases of evolution pass very quickly. The youngest massive stars, however, are associated with UCHii (ultra-compact Hii) regions, a phase that lasts longer than expected (Wood & Churchwell, 1989a). UCHii regions were defined by Wood & Churchwell (1989a) as photo-ionised nebulae with diameters , electron number densities and emission measures that have not yet expanded out of their natal molecular cloud. The observed ionisation from such stars indicates that the ionising photon emission rates lie approximately between corresponding to ZAMS stars with spectral types B2–O5. Wood & Churchwell (1989b) estimated the lifetime of UCHii regions () using the fraction of main-sequence stars observed in the UCHii phase and an adopted value for the main-sequence lifetime (). This lifetime is an order of magnitude longer than that predicted by simple Strömgren sphere expansion ().
UCHii regions are interesting to study as they may provide a lot of information about the early life of massive stars including their properties and, due to the limited number of observed morphological classes, their environments. The Wood & Churchwell (1989a) survey included classifications of UCHii regions into morphological types: cometary (), core-halo (), shell (), irregular or multiple peaked structures () and spherical and unresolved (). Kurtz et al. (1994) made radio-continuum observations of 75 UCHii regions for which the proportions of morphological types agree remarkably well with those in Wood & Churchwell (1989a). Observed at higher spatial resolution, some apparently spherical morphologies uncover more ordered morphologies (Felli et al., 1984). The radio continuum survey of UCHii regions by Walsh et al. (1998) found that most sources were either cometary or irregular (ignoring the unresolved sources). De Pree et al. (2005) revised the morphological scheme by adding the bipolar classification and removing the core-halo morphology. The latter was due to the fact that essentially all UCHii regions were observed to be associated with large-scale diffuse emission that could be separated from the compact cores. The proportions of the morphologies in this survey compare well with the Wood & Churchwell (1989a) and Kurtz et al. (1994) surveys except that of sources were classed as shell-like. Hoare et al. (2007) reviewed the different morphological surveys and noted that a lot of morphologies classified as shell-like could also be classified as cometary.
As UCHii regions have had little time to significantly alter their natal environments, the morphologies that arise are thought to reflect the ambient density field. Cometary types in particular are interesting to study and to test numerical models against because they have a highly regular shape and can overcome the lifetime problem. In the past there were three major models that aimed to explain how UCHii regions form cometary shapes: the champagne flow model (Israel, 1978), the bow shock model (Reid & Ho, 1985; Mac Low et al., 1991; van Buren & Mac Low, 1992) and clumpy/mass-loading models (Dyson et al., 1995; Redman et al., 1996; Williams et al., 1996). A more recent idea says that the situation is probably best described by a combination of these i.e. hybrid models (Arthur & Hoare, 2006).
Cometary Hii regions are mostly observed at the edge of molecular clouds. Along with the shape, this is the reason they were first named “blisters” (Israel, 1978). This led to the champagne flow model, a modification of the simple Strömgren solution in which the ionising star lies in an inhomogeneous (not a uniform) density field. Early modelling by Tenorio-Tagle (1979), Whitworth (1979), Bodenheimer et al. (1979), Bedijn & Tenorio-Tagle (1981) and Tenorio Tagle et al. (1979) included discontinuities in density separating molecular clouds with the ISM. Bodenheimer et al. (1979) found that once the ionised gas reaches the edge of the dense cloud, it accelerates up to . More realistically, Icke (1979a), Icke (1979b), Icke et al. (1980) and Yorke et al. (1983) modelled stars in density gradients, clearly reproducing cometary morphologies and predicting double spectral line profiles corresponding to a velocity splitting of up to along the symmetry axis. Henney et al. (2005) carried out simulations of the photo-evaporation of a cloud with large-scale density gradients. They showed that an ionised flow is set up that has a transient phase with duration , which then becomes approximately stationary with respect to the ionisation front (the quasi-stationary phase) for a large part of its evolution.
Real Hii regions cannot be explained entirely in terms of a density gradient. Not only does the champagne model fail to exhibit limb-brightened morphologies, but also the lifetime problem is not solved due to the unconstrained expansion of the ionisation front. It also does not take into account that OB stars drive stellar winds that can significantly affect the dynamics of nearby gas.
If an ionising star moves supersonically with respect to a uniform density field, the strong stellar wind driven by that star balances the ram pressure produced by the ambient flow, forming a bow shock ahead of the star. This model was introduced by van Buren et al. (1990) as an alternative explanation of cometary Hii regions and later analysed numerically by van Buren & Mac Low (1992). They found that behind the bow shock neutral material is swept up into a very dense, thin shell that traps the expanding ionisation front. This was a good result for the bow shock model as it was known that a confining mechanism was necessary to solve the lifetime problem. Limb-brightening was also found to occur in this model and can be explained as ionised gas that has been compressed by the stellar wind into a thin shell, flowing around the stellar wind cavity. The bow shock model was attractive because the size and shape compared well with observed cometary UCHii regions and it explained velocity gradients commonly seen at the head of cometaries.
Zhu et al. (2015) simulated the evolution of cometary Hii regions in bow shock and champagne flow models and analysed the resulting [Neii] and H2S(2) line profiles. They discovered that the two models can be distinguished using these lines.
Both the bow shock model and the champagne flow model have been found to be inadequate in explaining the recombination line velocity structures of cometary Hii regions. Lumsden & Hoare (1996) investigated the validity of both models in reproducing the velocity structure of the cometary UCHii region G29.96–0.02 and found that individually they poorly describe the object. Gaume et al. (1994) suggested that by including stellar winds and non-uniform ambient density fields more realistic models could be made.
Franco et al. (2007) studied the evolution of Hii regions moving from the centre to the edge of structured molecular clouds with Gaussian cores surrounded by power-law ( and ) halos. They found that UCHii regions remain pressure-confined while they are inside the core but evolve into extended Hii regions as they move into lower density regions of the cloud. A variety of hybrid models were investigated by Arthur & Hoare (2006) including champagne flow plus stellar wind and a combination of bow shock and champagne flow models in which the density gradient and the strength of the stellar wind were varied. With the inclusion of a stellar wind, it was discovered that it is possible for a champagne flow to produce limb-brightened morphologies. This is due to the fact that the stellar wind creates a dense shell that acts as a barrier to photo-evaporated flows, which divert around it. Arthur & Hoare (2006) found that for simple bow shock models the line widths are highest ahead of the massive star and for simple champagne flow models the widths are highest in the tail. They found that for hybrid models, a slow moving star in a steep density gradient has larger velocity widths toward the tail, but a fast moving star in shallow density gradients has its largest velocity widths nearer the star. Making sure to remove the effect of any cloud velocity, a hybrid model that matches the line data of a cometary Hii region could be used to gather information about the density structure of the natal cloud.
Three-dimensional simulations of Hii regions expanding off-centre in turbulent, self-gravitating power-law cores were run by Mac Low et al. (2007). They found that the expanding Hii regions were roughly spherical and noted that this is consistent with the analytical results of Korycansky (1992). These results were also confirmed by Arthur (2007) who performed similar simulations for and power-law density cores.
New Galactic plane surveys are revisiting the lifetime problem and providing large, well selected samples. The Galaxy-wide RMS (Red MSX Source) survey found mid-IR bright compact Hii regions (Lumsden et al., 2013). Mottram et al. (2011) used the results of this survey and determined the lifetime of the compact Hii phase to be or of the source’s main-sequence lifetime. Davies et al. (2011) simulated the RMS results using a particular Galactic gas distribution and different accretion models and compared with the luminosity distribution of the RMS survey. In this work each ionising star in the Galaxy was assumed to be producing Hii regions in a uniform density medium and not blowing a stellar wind (i.e. simple Strömgren expansion).
The aim of the current paper is to produce a grid of more realistic UCHii regions to include in the galaxy model of Davies et al. (2011) by simulating cometary Hii regions. We explore a parameter space spanning stellar mass, the density of the stellar environment and the age of the star and note the behaviour of some observables across this space. In a subsequent paper we will use this model grid to provide more realistic Hii region sizes and fluxes in an improved galaxy model that will be tested against the CORNISH (Co-Ordinated Radio ‘N’ Infrared Survey for High-mass star formation) survey (Hoare et al., 2012; Purcell et al., 2013). The CORNISH survey is an arcsecond resolution radio survey of the northern half of the GLIMPSE region (10\text{\,}\mathrm{\SIUnitSymbolDegree}$<l<$65\text{\,}\mathrm{\SIUnitSymbolDegree}, 1\text{,}\mathrm{\SIUnitSymbolDegree}$$) for compact ionised sources. With such a large unbiased sample of UCHii regions () it will be possible to test the models in the context of high-mass star formation on the Galactic scale.
In section 2 we introduce the numerical scheme and describe the models we simulated. The results of the simulations are presented in section 3 and we discuss the behaviour of the hot stellar wind region, the ionisation front, the emission measures and the spectral indices when stellar mass, age and cloud density are varied. We conclude the paper in section 4 where we summarise our findings.
2 The Model
2.1 Numerical Scheme
Our simulations use TORCH (The Operator-split Radiation, Cooling, Hydro-code), which is a 3D Eulerian fixed grid fluid dynamics code (code tests are given in appendix A). Three coupled physics problems (hydrodynamics, radiative transfer, and heating/cooling) are integrated separately and the result of each is combined to update the solution at each step using the Strang splitting scheme (Strang, 1968).
The fluid dynamics was solved on a two-dimensional axisymmetric grid (Falle, 1991) using a rotated hybrid HLLC-HLL Riemann solver (see appendix B) to calculate numerical fluxes on each inter-cell boundary. Left and right Riemann states are interpolated via a MUSCL scheme (van Leer, 1979) that uses the van Albada slope limiter (van Albada et al., 1982). Rotated hybrid Riemann solvers automatically apply fewer wave solvers normal to shocks (the HLL solver in this code) in order to eliminate carbuncle instabilities. These solvers are very robust and come with an acceptable drop in performance. The minimum and maximum signal velocities in the HLLC solver are approximated as the minimum and maximum eigenvalues of the Roe matrix (Roe, 1981). The hydrodynamics scheme is 2nd order accurate in time and space.
The transfer of ionising radiation is calculated using the C2–Ray method by Mellema et al. (2006). At each step the optical depth to each cell is traced via the method of short-characteristics. In this method the optical depth to the target cell is a weighted average of optical depths in neighbouring cells that lie between the target cell and the cell containing the ionising source. The ionisation fraction in each cell is implicitly solved before the cell’s optical depth can be used in calculations for other cells (which constrains the order in which the cells are iterated over).
Heating from different types of radiation and cooling due to collisional excitation and recombination is calculated using the approximate functions in Henney et al. (2009).
The governing hydrodynamic equations are: conservation of mass,
[TABLE]
conservation of momentum,
[TABLE]
and conservation of energy,
[TABLE]
where is the gas density, u is the fluid velocity, is the thermal pressure, is the total energy density with , and are respectively the heating and cooling rates due to atomic/molecular transitions, is the injection rate of stellar wind energy density and is the injection rate of wind material density as a function of position.
We also have equations describing the advection of the hydrogen ionisation fraction,
[TABLE]
and the rate of hydrogen ionisations and recombinations,
[TABLE]
where is the fraction of hydrogen that is ionised, is the electron number density, is the photo-ionisation rate, is the collisional ionisation coefficient and is the recombination coefficient of hydrogen.
Models were simulated on a numerical grid with square cells that have equal side lengths. The resolution of each grid is given in table 1 along with the physical dimensions. On the edges of the simulation domain we applied Dirichlet boundary conditions. Specifically, we applied reflective conditions on the boundary at and outflow conditions on the other boundaries.
2.2 The Star’s Environment
Our model environment has the same density structure as Model F in Arthur & Hoare (2006). Stars in these models are off-centre in a spherically symmetric density distribution at a distance of from the cloud centre. The density is given by
[TABLE]
where is the distance from the cloud centre, 0.01\text{,}\mathrm{pc}$$ is the cloud core radius, is the density at the cloud centre and the power-law index parameterises the dependence on when . The density at the star’s position is then given by
[TABLE]
where 0.35\text{,}\mathrm{pc}$$ is the distance of the star from the cloud centre, so that .
The power-law index, , has been inferred using a number of different techniques leading to a wide range of values in the literature. Hatchell & van der Tak (2003) found indices between and . Five of the seven dark cloud envelopes that were investigated by Arquilla & Goldsmith (1985) were best characterized by an index of . A range of was found by van der Tak et al. (2000). Pirogov (2009) found an index of that falls more steeply in the outer layers of the dense core. For the current work a value of was adopted.
The initial pressure in a cloud of constant temperature will have the same structure as the density field. Without gravity, gas will move down the pressure gradient leading to bulk motion that can interfere with the dynamics of the star’s stellar wind and ionisation field, especially over the period of each of the simulations here run for. Instead the pressure of the environment was taken to be uniform such that lower temperatures occur towards the cloud core, and higher temperatures occur away from it (the temperature at the position of the star was ).
2.3 Parameters
Using this model environment we explored key parameters, namely the cloud density at the position of the star and the mass of the ZAMS star. Values for these are given in table 1. In table 2 we give the stellar wind parameters, mass-loss rate and terminal velocity (see table 2), that were calculated using the predictions of Vink et al. (2001). These depend on the metallicity (which was assumed to be solar) and also the effective temperature and luminosity of the star. We took values for these (for certain stellar masses) from Davies et al. (2011), which were calculated using the hydrostatic models of Meynet & Maeder (2000). For the same masses, Lyman continuum fluxes were also taken from Davies et al. (2011) who used calculations from Martins et al. (2005) and Lanz & Hubeny (2007). Effective temperature, luminosity and Lyman continuum flux depend only on stellar mass; hence, this is the only free parameter describing the star.
2.4 Stellar Winds
Vink et al. (2001) determined relations for mass-loss rates of O and B stars for a range of metallicities by using a Monte Carlo radiative transfer method (Vink et al., 1999, 2000) on model atmospheres produced by the Improved Sobolev Approximation code (de Koter et al., 1993, 1997). This method allowed them to find the radiative momentum transferred to the wind via photon absorptions and scatterings. Self-consistent solutions, i.e. models for which the radiative momentum was equal to the wind momentum of the model atmosphere, were then used to find mass-loss rates. Strictly speaking, the Vink et al. (2000) model is applicable in the range 15\text{\,}\mathrm{\mathrm{M_{\sun}}}$\leq M_{\star}\leq$120\text{\,}\mathrm{\mathrm{M_{\sun}}}; we are assuming the author’s results can be extended to predict stellar wind parameters for our model stars with 15\text{,}\mathrm{\mathrm{M_{\sun}}}$$.
Mass-loss rates and terminal velocities were calculated for each of our model stars using the values in table 2 and these are plotted in figures 1 and 2, respectively. As can be seen in the table the mass injection rate increases for higher mass stars except between and where there is a decrease. This is due to the bi-stability jump; the star with 6\text{,}\mathrm{\mathrm{M_{\sun}}}$$ is the only one that lies on the cold side of the bi-stability jump where Fe IV recombines to Fe III, an efficient line driver.
In order to simulate the effects of a stellar wind we use the thermal energy injection method developed by Chevalier & Clegg (1985), which has also been used by Comeron (1997) and Arthur & Hoare (2006). Due to the resolution of the numerical grid used to simulate our models we were unable to define a free flowing wind region around the star (as in Rozyczka, 1985) that is sufficiently resolved (low resolution regions seed instabilities that grow enough to render the simulation results unphysical). The stellar wind power is mostly converted to thermal energy at the reverse shock (Castor et al., 1975) so we do not need to reproduce the structure of the unshocked wind region. Instead we can inject the wind luminosity as a rate of mass, , and energy, , into the shocked wind region as it is this power that determines the evolution and structure of the bubble (Weaver et al., 1977). To ensure that the gas in the injection region of our models is adiabatic the gas is not cooled within the injection region in our simulations.
This method also leads to a shocked stellar wind region with the correct temperature. Inside the injection region the temperature is set by the energy injected per unit mass () and the average particle mass. It is nearly uniform in this region, but drops rapidly at its edge as the gas expands supersonically into its surroundings. The supersonic gas at some point passes through a reverse shock, and is re-thermalised to approximately the same temperature as within the wind injection region. The post-shock temperature of the gas is insensitive to the exact location of the shock. If the shock is very near to the edge of the injection region, the pre-shock gas will be only mildly supersonic. The relatively weak shock and consequently relatively small temperature jump in this case is offset by the higher pre-shock temperature of the gas. On the other hand, if the shock occurs at a greater distance the pre-shock temperature of the wind may drop substantially lower, but its Mach number rises commensurably, and the kinetic energy per unit mass of this gas asymptotes to the original thermal energy per unit mass (ready to be thermalised again). If the reverse shock occurs inside the wind injection region the energy and mass within the region builds (keeping the gas at the original temperature) until its pressure overcomes that of the surroundings.
Wind material density and energy density are added within the injection radius at each time-step such that their integrated rates over the volume are and respectively. The rate of injected wind energy density is given by , where is the rate of injected wind material density:
[TABLE]
where is the injection radius, which is given in table 1 for each model.
We tested and confirmed that the evolution of a spherically symmetric wind-blown bubble in a neutral medium with no radiative cooling is as predicted using this method. The size of the injection region ( cells) in all our models was chosen to be large enough that the region is well resolved by the numerical grid and small enough that the injection region has little effect on the evolution of the bubble. An undesirable effect of this is that the injection radius is the radius the wind bubble starts with so that the wind has had a head-start. It does, however, take time to inject enough pressure in the region to blow a wind. This is approximately the time it takes to add pressure into the injection region that is comparable to the ambient pressure i.e.
[TABLE]
where , and are the pressure, hydrogen number density and temperature of the ionised ambient gas, is the Boltzmann constant and is the ratio of heat capacities. The start times, , are given in table 1 for all the models.
2.5 Radiation Field
Ionising radiation from each star was assumed to be monochromatic. Non-ionising FUV (far ultra-violet) radiation was also included to heat the gas in the same way as Henney et al. (2009). It was assumed that recombinations to the ground state of hydrogen are locally reabsorbed i.e. the diffuse field was treated under the on-the-spot approximation. The recombination coefficient was therefore calculated using cubic spline interpolation of the case-B recombination coefficient data in Hummer (1994). Similarly, the collisional ionisation of neutral hydrogen was calculated by cubic spline interpolation of the data in Raga et al. (1997). The cross-section for ionising photons was taken from Osterbrock (1989). Figure 3 shows the Lyman continuum flux plotted against stellar mass for our models.
3 Results and Discussion
Plots of hydrogen number density, pressure, temperature and ionised hydrogen number density are given in figure 4 for a typical cometary Hii region around a star with mass 30\text{,}\mathrm{\mathrm{M_{\sun}}} at an age of $t=$50\text{\,}\mathrm{kyr} that starts the simulation in a local hydrogen number density of 3.2\text{\times}{10}^{4}\text{,}{\mathrm{cm}}^{-3}$$. Nearest the star is the wind injection region of radius from which a wind flows outwards becoming supersonic at its edge. The central temperature is and decreases outside of the injection radius to the wind reverse shock. Gas outside the injection region is accelerated to supersonic speeds. We refer to the gas between the wind injection region and the reverse shock as the “unshocked stellar wind region”. Downstream of the reverse shock is low density shocked stellar wind material, which is at a uniform temperature of and subsonic. Further downstream lies the contact discontinuity separating the shocked stellar wind region and the cooler () ionised ambient gas. Earlier on in the evolution of this Hii region we would expect to see a shock ahead of the contact discontinuity, but pressures either side of the shock have equilibrated so that the shock no longer exists. Finally, the ionised region is bounded by an ionisation front that separates the ionised gas from the neutral gas shock. Photo-ionisation has heated the gas in the ionised region and therefore has created an over-pressure with respect to the surrounding ambient gas, which explains the presence of the shock just outside of the ionisation front.
3.1 Ionisation Fronts
Assuming balance between photo-ionisations and recombinations, for a star in a uniform medium with no stellar wind we have:
[TABLE]
where is the hydrogen number density inside the ionised region, is the case-B recombination coefficient for hydrogen and is the time-dependent radius of the ionised region. The Strömgren radius is the initial radius of the ionised region:
[TABLE]
where is the initial hydrogen number density of the uniform medium. Ionisation heats the gas and sets up an over-pressure; expansion thereafter is driven by this over-pressure.
Raga et al. (2012b) show that if the shock driven by the over-pressure is isothermal and it is assumed that the expanding Hii region is in pressure balance with the shocked neutral gas, the equation of motion of the ionisation front is
[TABLE]
The Spitzer approximation (Spitzer, 1978) can be reproduced by this equation by noting that the second term on the right-hand side of equation is negligible at early times:
[TABLE]
where is the sound crossing time-scale and is the sound speed of the ionised gas.
According to Raga et al. (2012a) equation 12 neglects the inertia of the expanding neutral material and therefore derived solutions should under-estimate the distance to the ionisation front. Including the inertia gives the equation of the shell,
[TABLE]
At early times the solution can be approximated (Hosokawa & Inutsuka, 2006) by
[TABLE]
Simulations of D-type expansion of ionisation fronts (Bisbas et al., 2015) have shown close agreement with equation 15 at early times and with equation 12 at late times.
The sizes of the ionised regions simulated in the current paper are shown in tables 3 and 4, along with the analytical radius of ionised regions that would evolve around stars with the same stellar parameters but in a uniform density medium (equation 12). The isothermal sound speed in the ionised ambient gas was assumed to be , where the ionised gas temperature is 8000\text{,}\mathrm{K}, the average molar mass of ionised hydrogen is $\mu_{\mathrm{H}}=$0.5\text{\,}\mathrm{g}\text{\,}{\mathrm{mol}}^{-1} and is the gas constant, giving a sound speed of 11.5\text{,}\mathrm{km}\text{,}{\mathrm{s}}^{-1}$$. In table 3 the measured ionisation regions show a decrease in size for higher local densities as expected. The sizes also increase for higher masses due to higher Lyman continuum photon fluxes. There is generally close agreement with the analytical radius, with the simulated regions consistently slightly larger. This is due to the contribution of the wind to driving the shock ahead of the ionisation front. Simulations of stars without stellar winds confirm this as they produce Hii regions with radii that are very close to those predicted by equation 12. For higher-mass stars the radius of the simulated ionised region is noticeably larger than the analytical radius, suggesting that the strong winds from these stars are contributing relatively more to the expansion of the Hii region.
As a Hii region expands the pressure inside the region drops until it is equal to the ambient gas pressure, at which point it stagnates. The stagnation radius can be found, by setting in equation 12, to be
[TABLE]
where is the isothermal sound speed of the neutral ambient gas. The temperature of the neutral ambient gas is 300\text{,}\mathrm{K}, so the sound speed is $c_{\mathrm{a}}\simeq$1.58\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}. In our models we therefore have . Setting in equation 13 we can get an approximation for the stagnation time in terms of the sound-crossing time of the initial Strömgren sphere: 60. We can see that for smaller Strömgren radii the stagnation time occurs earlier. Table 4 shows stagnating ionisation fronts that occur at earlier ages for smaller-mass stars.
Most of the simulations of our models did not show the unbounded ionisation front expansion predicted by Franco et al. (1990) for . Only for stars with masses 40\text{,}\mathrm{\mathrm{M_{\sun}}} and local densities $n_{\star}\leq$1.6\text{\times}{10}^{4}\text{\,}{\mathrm{cm}}^{-3} did the ionisation front break free of the cloud and almost always in finger-like structures (a.k.a. the shadowing instability, see Williams, 1999). In all other models the regions are bounded approximately by a sphere, which is consistent with the three-dimensional simulations of off-centre UCHii regions by Mac Low et al. (2007). The regions are at a roughly constant density throughout the simulation and grow to a maximum radius that coincides with when the pressures either side of the ionisation front have equalised.
These results very closely match those in Arthur (2007), in which the condition necessary for Hii regions in power-law density environments to remain bounded during their initial formation stage was given by:
[TABLE]
where i.e. the ratio of the Strömgren radius for a star in a uniform medium to the distance between the star and the cloud centre. Values of this ratio for each set of model parameters were calculated and are given in table 1 by setting in equation 11 equal to the local hydrogen number density, . Our power-law density environments have , so equation 17 reduces to , which is true for all our models. Arthur (2007) also found that in power-law environments, if then pressure balance can halt the breakout of the Hii region during the expansion stage. As mentioned before the pressures do equalise before the Hii regions can become unbounded for most models. However, this is also seen for models with .
3.2 Shocked Stellar Wind Region
Tables 5 and 6 list the sizes of the shocked stellar wind regions as a function of density and time respectively. In table 6 the sizes of the shocked stellar wind regions show a plateau at early times for low-mass stars. The plateau occurs at later times for higher-mass stars. This suggests that the bubbles are becoming pressure-confined (plots of the pressure confirm this).
We use the results of an analysis of stellar wind evolution by Koo & McKee (1992) to see if this behaviour is physical. The terminal wind velocity in all of the models is lower than the critical wind velocity defined in Koo & McKee (1992, see their equation 2.5), indicating that the wind-bubbles are radiative. In this regime the cooling time of the shocked wind is shorter than the time it takes to accumulate a significant mass of gas ahead of the shock. We found that the radii of the shocked stellar wind regions are consistent with pressure-confined fully radiative bubbles.
Numerical diffusion can cause extra cooling, i.e. the cooling length might not be sufficiently resolved and therefore intermediate temperatures are found in a higher volume of the gas. Kelvin-Helmholtz instabilities can also enhance cooling by increasing the surface area of the contact discontinuity. Increasing the resolution of the numerical grid increases the growth of these instabilities to further enhance cooling. These behaviours are competing but overall the level of numerical diffusion has been found to decrease with higher resolution (Parkin & Pittard, 2010). We carried out resolution tests that confirm that the cooling across the contact discontinuity was not enhanced significantly by numerical effects.
The shell radius for a radiative bubble is given by
[TABLE]
where is the mechanical wind luminosity and is the ionised ambient density. For stars with 70\text{,}\mathrm{\mathrm{M_{\sun}}}$$ the wind is driving a shock that follows equation 18 just ahead of the ionisation front, which explains why we are seeing large ionisation front radii in tables 3 and 4 for these stars. The shock ahead of the ionisation front for the lower-mass stars evolves according to equation 12, indicating that the wind-blown bubble has been pressure-confined by the Hii region.
A radiative bubble eventually becomes pressure-confined when the ram pressure of the shell is equal to the ambient pressure. Using equation 18 we can therefore obtain the time at which confinement occurs:
[TABLE]
where is the isothermal sound speed in the ionised ambient medium. The ionised ambient density, , was determined for each snapshot as the density upstream of the shock in the ionised ambient region in the radial direction. Using equation 10 we can get an approximation for this density:
[TABLE]
where is the measured radial distance from star to ionisation front given in tables 3 and 4 ( gives an “average” value for the bubble radius).
For most of the models the transition time is too short to be seen either because the earliest data snapshot was taken at or because the bubble is small and therefore not well resolved. As a result the bubbles quickly stall and consequently are missing the shock outside of the contact discontinuity that separates shocked stellar wind material and the ionised ambient medium. As the shock is missing and the wind-blown bubbles are behaving like radiative bubbles, the radial distance from the star to the contact discontinuity, in tables 5 and 6, should be comparable to the analytical radius of the shock at the confinement time given in equation 19. The approximate final radius of the shocked wind region in the radial direction is
[TABLE]
which is also given for each model in tables 5 and 6. Looking at these tables, the measured transverse distances, , show closer agreement with the confinement radius, , for higher-mass stars. In tables 5 and 6, we can see generally better agreement for stars in lower density environments and at later times respectively. These results suggest that stars with high enough mass have stellar winds that dominate the expansion of the Hii region. For lower-mass stars the shocked stellar wind region is noticeably larger than predicted, suggesting that the expansion of the Hii region is also contributing to the expansion of the shocked stellar wind region.
In all of the models the contact discontinuity eventually stalls. For stars with mass 6\text{\,}\mathrm{\mathrm{M_{\sun}}}$\leq M_{\star}\leq$30\text{\,}\mathrm{\mathrm{M_{\sun}}} in table 6 the stellar wind size in the transverse direction, , plateaus before they reach an age of . This behaviour is easily seen for the star with mass 30\text{,}\mathrm{\mathrm{M_{\sun}}}$$ but is harder to see for lower-mass stars because the contact discontinuity between the shocked stellar wind region and the ionised ambient region is unstable.
3.3 Emission Measures
In order to see what our models would look like if the ionised region is optically thin we calculated the emission measure:
[TABLE]
where is the number density of free electrons and the integration is along a line of sight.
To do this integration we used the cylindrically symmetric ray-tracing scheme presented in the appendix of Dougherty et al. (2003). Grids of the resulting emission measures at a viewing projection angle of are shown spanning two dimensional parameter spaces in figures 5 and 6. In figure 5 the effects of varying stellar mass and the star’s local density were explored at a time of . The evolution over time of emission measures for different stellar masses is plotted in figure 6 at a local number density of 3.2\text{\times}{10}^{4}\text{,}{\mathrm{cm}}^{-3}$$.
The figures clearly show two distinct behaviours. For stars with mass 9\text{,}\mathrm{\mathrm{M_{\sun}}} (and also $M_{\star}=$12\text{\,}\mathrm{\mathrm{M_{\sun}}} for 3.2\text{\times}{10}^{4}\text{,}{\mathrm{cm}}^{-3}) the morphology of the cavity blown out by the stellar wind varies over the simulation time. In these models the stellar wind is weak and is therefore isotropically confined by the Hii region, explaining the non-cometary morphologies. A pressure gradient exists in these Hii regions, which leads to the flow of material down the density gradient. Stars with higher mass have wind bubbles that are not confined down the density gradient, leading to expansion in this direction. These stars have cometary morphologies that do not change appreciably over a period of $200\text{\,}\mathrm{kyr}$ and show limb-brightening. Stars with masses $M_{\star}\geq$70\text{\,}\mathrm{\mathrm{M_{\sun}}} also show limb-brightening if 3.2\text{\times}{10}^{4}\text{,}{\mathrm{cm}}^{-3}$$. The rest have Hii regions that expand past the centre of the dense cloud so that this is the only feature picked up in the emission measure images. For the Hii region in the bottom left of figure 5 we see limb-brightening because the high density core is destroyed by the stellar wind.
For one of the models that does exhibit limb-brightening we show in figure 7 the effects of changing the projection angle on the morphology. The projection angle is defined as
[TABLE]
where is a unit vector directed from star to cloud centre and is a unit vector directed from star to observer. At projection angles closer to limb-brightening is more pronounced i.e. the “limbs” wrap further around the centre of the object. Using the classification scheme of De Pree et al. (2005) the projected morphology is less cometary and more shell-like the closer the axis of symmetry of the Hii region is oriented towards the observer. These results suggest that Hii regions classified as having a shell-like morphology may be cometaries viewed along their axis of symmetry (as noted by Mac Low et al., 1991).
3.4 Spectral Indices
The change in specific intensity, at frequency , across a ray segment in a grid cell is
[TABLE]
where is the specific intensity at the start of the path segment, and are the emission and absorption coefficients respectively, and is the optical depth across the segment given by
[TABLE]
where is the path length through the cell. The emission coefficient (Rybicki & Lightman, 1979) is
[TABLE]
where is the elementary charge, is Coulomb’s constant, is the mass of an electron, is the ionic charge, is Planck’s constant, and is the velocity averaged Gaunt factor (Hummer, 1988). Using Kirchoff’s law,
[TABLE]
where is the Planck function, we have the absorption coefficient:
[TABLE]
We produced synthetic radio continuum maps by accumulating the changes in intensity along rays through the finite cells in our hydrodynamic solutions to pixels on an image plane. The flux density, , in each pixel was obtained by integrating the intensity over the solid angle of the pixel. Spectral indices were approximated using these maps:
[TABLE]
where and are two closely separated frequencies.
A spectral index map for radio continuum emission is plotted in figure 8 for a star in the centre of the parameter space we are exploring for frequencies and . Intermediate spectral indices (between and ) are found across the whole image. At the head and the limbs of the Hii region, where the emission measure is highest, the spectral indices are higher () than those in the cavity blown out by the stellar wind ().
Table 7 shows the spectral indices of models for a range of local densities and masses. These indices were calculated using the total radio continuum flux of each model for frequencies of and . At these frequencies the objects are more optically thick the higher the stellar mass (this is also seen in table 8). Higher-mass stars have more intense ionising radiation and blow stronger winds, so they ionise deeper into the dense cloud, leading to a higher proportion of the total flux from areas of high density. For masses 40\text{,}\mathrm{\mathrm{M_{\sun}}} the spectral indices decrease, which is likely due to the competing effect of larger regions having a higher proportion of flux from optically thin regions. These stars are ionising the regions around the densest part of the cloud (or through it in the case of the star with $M_{\star}=$120\text{\,}\mathrm{\mathrm{M_{\sun}}} in a local density of 12.8\text{\times}{10}^{4}\text{,}{\mathrm{cm}}^{-3}$$).
In table 8 spectral indices are shown to decrease over time. The density in a Hii region decreases as it expands (after its initial Strömgren expansion phase), so the spectral indices should decrease also.
Figure 9 shows the radio spectrum of one of our model Hii regions and also a Hii region that is produced by the same star but in a uniform density environment. Surprisingly, the turn-over frequency is very similar for the power-law case, but it is slightly lower because the Hii region covers more low-density gas down the slope than high-density gas up the slope. The total fluxes have a higher proportion from lower density gas and therefore the integrated spectral indices will lean closer to indices from lower density regions. The close similarity between the SEDs is consistent with the fact that the density slope of ionised gas is much shallower than the initial slope in the neutral ambient gas.
3.5 Model Limitations
If we had included the effects of the entire spectrum of radiation from our model stars the temperature structures of our Hii regions would be different. We would see higher temperatures near the ionisation front due to the hardening of the radiation field (Wood & Mathis, 2004). The extra pressure due to these high temperatures is expected to be higher than in our model Hii regions.
Another limitation of our model is that we took the initial thermal pressure of the neutral gas to be uniform to prevent gas from moving away from the cloud centre (which would happen if the temperature was uniform instead). In a future study the effects of gravity could be included to better model these environments.
Dust is an important factor in the evolution of high density Hii regions that we did not consider. Ionising photons can be absorbed by dust and therefore dusty Hii regions are expected to be smaller. Arthur et al. (2004) simulated dusty compact Hii regions and found that the fraction of ionising photons absorbed by dust decreases as the Hii region expands, and so the effects of dust are less important in the later phases of expansion. They also found that dusty Hii regions stagnate earlier than their dust-free counterparts.
4 Conclusions
We have simulated cometary Hii regions that result from stars that blow stellar winds and are placed off-centre from a power-law density gas cloud. Models with stars of mass 12\text{,}\mathrm{\mathrm{M_{\sun}}}$$ produce limb-brightening. Some are diminished if the star’s Hii region contains the dense centre of the cloud (which is true for stars with high Lyman output near low density clouds). There is a cavity blown out by the stellar wind in each model, which is observed in real cometary Hii regions. The morphological class of a Hii region depends on the viewing projection angle. Morphologies are more shell-like the smaller the angle between the viewing direction and the object’s axis of symmetry. Spectral indices are higher for higher-mass stars, which can be explained as larger Hii regions covering more of the high density gas near the cloud centre. This effect turns over at even higher masses as the volume of ionised low density gas increases more than the volume of ionised high density gas in the cloud centre. Model spectral energy distributions for power-law density models were almost identical to those of the same stars placed in uniform density surroundings, which is consistent with the shallow density structures we found in our simulated ionised regions. The turn-over frequency in the radio spectrum for our power-law density models is slightly lower compared to the models with uniform ambient density. This is because Hii regions in a power-law density medium cover more low density gas than high density gas, which means the proportion of emission from optically thin regions is higher.
In a subsequent paper we will use the grid of simulated models presented in this paper to select Hii regions that correspond to regions generated in a model Milky Way. The ultimate goal is to compare observables from this simulated survey with the CORNISH survey in order to investigate how well our models describe Hii region evolution.
Acknowledgements
We thank the anonymous reviewer for their critical reading of this paper and for their suggestions. This work was supported by the Science and Technology Facilities Council at the University of Leeds.
Appendix A Code Tests
A.1 R-Type Expansion and Shadowing
Two tests were carried out to ensure the accuracy of the radiative transfer scheme. The first test compares the radius of a simulated Strömgren sphere with the analytical solution for the Strömgren radius (Strömgren, 1939):
[TABLE]
where is the recombination time. We assume here that 1\text{\times}{10}^{48}\text{,}\mathrm{phot}\text{,}{\mathrm{s}}^{-1}, $n_{\mathrm{H}}=$100\text{\,}{\mathrm{cm}}^{-3}, and 2.59\text{\times}{10}^{-13}\text{,}{\mathrm{cm}}^{3}\text{,}{\mathrm{s}}^{-1}$$.
The results of this test are shown in figure 10 along with a comparison between first-order and second-order explicit schemes and the implicit scheme TORCH uses. Clearly the implicit scheme out-performs the other schemes, converging to the analytical solution in a time much smaller than the recombination time.
Simulations of UCHii regions with stellar winds should show trapping of the ionisation front behind the swept-up wind bubble. To see if the code works accurately in such a situation, the Strömgren sphere test was repeated but with a dense square clump positioned near the ionising source. The clump effectively blocks photo-ionising radiation from propagating into the “shadow” region. There is, however, some numerical diffusion into the shadow region due to the interpolative nature of the short-characteristics ray-tracing scheme. This will have a negligible effect on the dynamics of any simulated Hii region. Figure 11 shows the evolution of the ionised Hii fraction in the shadow test.
A.2 D-Type Expansion
TORCH was benchmarked against the results of the one-dimensional spherically symmetric D-type expansion tests of Bisbas et al. (2015). The early and late phase simulations, described by this author, of D-type ionisation front expansion into an initially uniform density (5.21\text{\times}{10}^{-21}\text{,}\mathrm{g}\text{,}{\mathrm{cm}}^{-3}) isothermal gas were run. The ionising star in this simulation has a Lyman continuum photon rate of $Q_{\mathrm{Lyc}}=$1\text{\times}{10}^{49}\text{\,}\mathrm{phot}\text{\,}{\mathrm{s}}^{-1}, which produces an initial Strömgren sphere of radius 0.314\text{,}\mathrm{pc} and temperature $T_{\mathrm{HII}}=$10\,000\text{\,}\mathrm{K}. The temperature of the neutral ambient medium was set to 100\text{,}\mathrm{K} for the early phase and $T_{\mathrm{HI}}=$1000\text{\,}\mathrm{K} for the late phase. Coupling between gas and radiation physics was achieved for an isothermal equation of state by setting the temperature of a gas element according to the ionised fraction of hydrogen in that element:
[TABLE]
where is the temperature of fully neutral hydrogen and is the temperature of fully ionised hydrogen.
Early on in the evolution (a few recombination times) a weak-R ionisation front expands according to equation 30. The gas is heated so that the ionised region is over-pressured with respect to the neutral gas. As the ionisation front slows down it becomes R-critical; the pressure wave overtakes the front and steepens into a shock wave, compressing the gas behind it. The R-critical ionisation front transitions into a D-critical front with an isothermal shock ahead of it. D-type expansion is then driven by the over-pressure on the sound crossing time-scale , where is the Strömgren radius given in equation 30 and 12.85\text{,}\mathrm{km}\text{,}{\mathrm{s}}^{-1}$$ is the sound speed of the ionised gas. Figure 12 shows the expansion of the ionisation front in the early phase. The expansion begins by following the Spitzer solution (equation 13) and eventually follows the solution by Hosokawa & Inutsuka (2006) (equation 15). In figure 13, we show the ionisation front evolution in the late phase. The simulated front is a close fit to the StarBench equation (Bisbas et al., 2015), with an error of less than for the simulation duration.
A.3 Wind-Blown Bubble
The wind-blown bubble model of Strickland & Stevens (1998) was simulated using TORCH to test if the cooling instability in the shell is reproduced. A stellar wind with mass loss rate 5\text{\times}{10}^{-5}\text{,}\mathrm{\mathrm{M_{\sun}}}\text{,}{\mathrm{yr}}^{-1}$$ and terminal wind speed was modelled by including source terms and within a radius of 10 grid cells (see section 2.4).
As was found in Strickland & Stevens (1998) the wind sweeps material into a shock-heated shell. The shell starts out thick, but rapidly cools and collapses as it grows. Eventually the shell is so thin that it becomes unstable to thin shell and Vishniac instabilities (Vishniac, 1983). The evolution of the ionised hydrogen number density during the simulation is shown in figure 14 and agrees well with the model in Strickland & Stevens (1998).
A.4 Shadowing Instability
This final test utilises both the radiation and cooling modules in order to reproduce the shadowing instability (Williams, 1999). A wind was set up in the same way as in section 2.4 with wind, star and initial gas parameters taken from the shadowing instability test in Arthur & Hoare (2006). If the cooling time is too short the dense swept-up shell can collapse while the effects of the staggered grid are still appreciable (i.e. the shell has a low resolution). For this reason the energy flux from the heating and cooling module was artificially reduced by a factor of in order to increase the cooling length and so delay the onset of the cooling instability.
Figure 15 shows the evolution of the ionised gas density (similar to figure 2 in Arthur & Hoare, 2006), which behaves as expected. The cooling instability discussed in section A.3 occurs and produces a variation in the optical depth, leading to the tell-tale spokes of the shadowing instability.
Appendix B Rotated Hybrid HLL/HLLC
To solve the Riemann problem across inter-cell boundaries we first implemented a HLLC (Harten-Lax-van Leer-Contact) Riemann solver. After testing it was apparent that this scheme is not robust as unacceptable flow fields can result called carbuncle instabilities (Quirk, 1994). Schemes such as HLL (Harten-Lax-van Leer) and Rusanov do not have this problem as they are highly dissipative. There are a few ways to cure schemes that are prone to the carbuncle phenomenon; the technique we decided to use was to combine HLL and HLLC schemes into a rotated-hybrid Riemann solver (Nishikawa & Kitamura, 2008).
The scheme starts by decomposing the geometric grid cell face normal n into two orthogonal directions:
[TABLE]
and
[TABLE]
where is the velocity difference vector, is a direction tangent to the geometric face and is a small number. The second case in equation 32 ensures that only the HLLC solver is used when streamwise velocity fields are smoothly varying; only at discontinuities will the HLL solver be applied. The flux across this grid cell face is
[TABLE]
Fluxes and are calculated by first finding the velocity in a new coordinate system:
[TABLE]
where v is the velocity in the original coordinate system on either the left or right side of the cell interface and R is the rotation matrix (defined as rotating decomposed Riemann problem directions, or , to align with the geometric cell face normal, n).
Using the Rodrigues’ rotation formula (Rodrigues, 1840) the rotation matrix R, which rotates unit vector onto unit vector n is given by:
[TABLE]
where
[TABLE]
and is the skew-symmetric cross-product matrix of w,
[TABLE]
With rotated left and right Riemann states for each decomposed Riemann direction we can now solve the one dimensional Riemann problem. Solving these problems leaves us with momentum fluxes (calculated using ) that need to be rotated back onto the original coordinate system:
[TABLE]
Along with the rest of the flux components, this gives us and . Using the rotated Riemann solver effectively applies the HLL solver in the direction normal to shocks (suppressing the carbuncle instability) and applies the HLLC solver across shear layers in order to minimise dissipation.
Liou (2000) argues that intermediate cells across a shock exchange information transversely to their neighbours along the shock, which can develop into a carbuncle instability. Considering a face connecting two neighbouring intermediate cells, any significant perturbation travelling in a direction normal to this face will lead to a velocity difference vector oriented in the same direction (Nishikawa & Kitamura, 2008). This means the HLL Riemann solver will introduce dissipation in this direction, suppressing the instability. The resulting scheme pays for its robustness with an acceptable drop in accuracy (using a dissipative Riemann solver) and speed (using two Riemann solvers).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1van Albada et al. (1982) van Albada G. D., van Leer B., Roberts Jr. W. W., 1982, A&A, 108, 76
- 2Arquilla & Goldsmith (1985) Arquilla R., Goldsmith P. F., 1985, Ap J , 297, 436 · doi ↗
- 3Arthur (2007) Arthur S. J., 2007, Ap J , 670, 471 · doi ↗
- 4Arthur & Hoare (2006) Arthur S. J., Hoare M. G., 2006, Ap JS , 165, 283 · doi ↗
- 5Arthur et al. (2004) Arthur S. J., Kurtz S. E., Franco J., Albarrán M. Y., 2004, Ap J , 608, 282 · doi ↗
- 6Bedijn & Tenorio-Tagle (1981) Bedijn P. J., Tenorio-Tagle G., 1981, A&A, 98, 85
- 7Bisbas et al. (2015) Bisbas T. G., et al., 2015, MNRAS , 453, 1324 · doi ↗
- 8Bodenheimer et al. (1979) Bodenheimer P., Tenorio-Tagle G., Yorke H. W., 1979, Ap J , 233, 85 · doi ↗