A&A 369, 263-268 (2001)
DOI: 10.1051/0004-6361:20010070
K. G. Pavlakis ^{1,2} - R. J. R. Williams^{3} - J. E. Dyson ^{2} - S. A. E. G. Falle^{4} - T. W. Hartquist ^{2}
1 - Physics Department, University of Crete,
714 09 Heraklion, Crete, Greece
2 -
Department of Physics and Astronomy,
University of Leeds, Leeds LS2 9JT, UK
3 -
Department of Physics and Astronomy,
University of Cardiff, PO Box 913,
Cardiff CF24 3YB, UK
4 -
Department of Applied Mathematics,
University of Leeds, Leeds LS2 9JT, UK
Received 9 October 2000 / Accepted 4 January 2001
Abstract
We study the evolution of a globule of neutral material immersed in
the more tenuous hotter plasma of an H II region surrounding
newly born OB stars. The neutral globule is illuminated by the direct
ionizing radiation of OB stars, and by diffuse radiation emitted by
recombination in the surrounding ionized gas. We perform 2D, time
dependent axisymmetric hydrodynamic simulations, and find that, for
values of the diffuse field of the order of 10% of the direct field,
the evolution of the globule is completely different to its evolution
when the diffuse field is neglected.
Key words: interstellar medium: clouds - interstellar medium: H II regions - kinematics and dynamics -
hydrodynamics
The formation of cometary tails behind globules in some planetary nebulae (e.g. Meaburn et al. 1992; O'Dell & Handron 1996) has been attributed to the shadowing of gas behind the globules from the stellar sources of ionizing radiation (e.g. Cantó et al. 1998). Similar explanations for the formation of tails behind globules in some regions of high-mass star formation (e.g. Bally et al. 1998) have been proposed (Richling & Yorke 1998; Störzer & Hollenbach 1999, and references therein). An alternative possibility is that a tail arises due to evaporation driven by the flow of a subsonic hot shocked stellar wind around a globule (Dyson et al. 1993).
The literature on the formation and development of cometary tails due to shadowing of radiation is very large. The most recent analytic and semi-analytic approaches to the problem include those of Bertoldi (1989); Bertoldi & McKee (1990); Henney et al. (1996, 1997); Johnstone et al. (1998); Mellema et al. (1998), and Störzer & Hollenbach (1999), while recent relevant numeric simulations include those of Lefloch & Lazareff (1994, hereafter LL94), Cantó et al. (1998); Mellema et al. (1998), and Richling & Yorke (1998). Of those who have conducted numerical simulations only Cantó et al. (1998) and Richling & Yorke (1998) took account of the effects of a diffuse component of the ionizing radiation field. Cantó et al. (1998) study the development of tails in cylindrical symmetry without considering the variation of the ionizing radiation with distance from the source, while the lifetime of the tail was extended in the work of Richling & Yorke (1998) as the angular momentum of the evaporated material prevented it from collapsing onto the symmetry axis.
In this paper we present results of simulations of the evolution of a non-gravitating, non-magnetic globule in response to the direct ionizing radiation from a star that turns on instantaneously and simultaneously with a diffuse component of the ionizing radiation field. In each simulation the globule was assumed to be initially of uniform density and temperature and of spherical shape. Section 2 contains the basic equations and descriptions of the assumptions, initial conditions, and the code used. We present our results in Sect. 3.
We followed the approach of LL94 as closely as possible but also included the effects of a diffuse component of the radiation field.
The equations that we solved include the continuity equation and
the momentum equation:
(1) |
(2) |
We assume that
(4) |
Figure 1: Approximate analytic solutions to the radiative transfer equation. a) Uniform density ionized gas b) 1/r density ionized gas (appropriate for an evaporative wind, however note that local ionization equilibrium has been assumed in both approximations). The solid lines show the two-beam approximation of Cantó et al., while the dashed lines show the results for our one-beam approximation | |
Open with DEXTER |
The calculation of the diffuse field is more tricky. Ideally, a full
radiative transfer solution would be required to provide an accurate
determination of the radiation field, but this is computational
expensive in a 2D, time-dependent hydrodynamic treatment. Instead,
we base the treatment in the present exploratory work on the
work of Cantó et al. (1998). These authors developed an iterative
2-beam treatment of the diffuse radiation field in cylindrical
symmetry, which has since been verified against a more complete
treatment of the radiative transfer equations by Razoumov & Scott
(1999). We further simplify their treatment by including only the
inward-going part of the radiation field. This is a reasonable
further approximation for the centrally-peaked density distributions
we find in our simulations, for which the radiation field will indeed
be predominantly inwards for the radii at which most of the
recombinations occur. We determine the variation of
^{} with distance from the axis, r, by integrating
In Fig. 1, we compare the results of the analytic
approximation to the two-beam treatment developed by Cantó et al. (for the limit of local ionization equilibrium) with analogous results
for our one-beam treatment. The neutral tail (if present) is assumed
to lie inside the radius, ,
where the inward-going
radiation flux reaches zero. This radius is shown as a function of
the parameters
(7) |
(8) |
We see from both panels of Fig. 1 that the one-beam approximation somewhat overestimates the radius of the column predicted for large diffuse fields (small or ). Since we are interested in the conditions under which tails are inhibited by hydrodynamic effects, the conservative nature of this approximation is satisfactory for exploratory work.
To calculate the pressure, we assumed the same relationship
between temperature, T and
taken by LL94.
That is, for densities greater than 20 per cent of the
initial cloud density:
To enable us to compare our results as closely with those of LL94 as possible, we assumed, as they did, that , which is smaller than the threshold value. A decreased value of this cross-section may be appropriate, because of the hardening of the radiation field close to the ionization front. However, the ionization fronts are unresolved in our simulations so the precise value of this constant is unimportant for our results.
Following LL94, we take the value of the recombination coefficient,
,
to be
independent of the gas temperature. However, because the IF
structures are not resolved, this may in fact be more accurate than
including the variation with temperature. In particular, it
circumvents problems with the advection speed of the ionization front
analogous to those found for underresolved detonation waves, as
described by LeVeque (1998). At local fractional ionization ,
the gas has temperature T given by Eq. (9). Within a
volume V, there are
neutral hydrogen atoms and
hydrogen ions, respectively, where n is the hydrogen nucleon
density. Now interpret
as determining the partition of the gas
into two phases in pressure equilibrium, one fully ionized and one
fully neutral. The cool, neutral phase has temperature ,
density
and occupies volume (1-f_{x})V, while the hot
ionized phase has temperature ,
density
and
occupies volume f_{x} V. Conservation of mass demands
and
.
Since there is
pressure equilibrium
.
Hence the total recombination rate within V is then
In the initial configuration, a uniform sphere of gas with a radius of 0.5 pc centred at , and cm^{-3} is surrounded by initially neutral gas with cm^{-3}. At t = 0, the radiation field switches on, and thereafter the flux of Lyman continuum photons is . Thus, our initial conditions are the same as those of LL94 (model 2) except for our inclusion of which we took to have a constant value at . For our model A, , and for our model B, .
For these values, the initial values of
are 0.92 (model A)
and 0.31 (model B), so no shadowed tail would be expected to form
(cf. Fig. 1a). In fact, it is difficult to find
circumstances in which the value of
would be substantially
greater than unity for a long thin tail, since
(11) |
The simulations were carried out with a two-dimensional hydrodynamics code based on an exact Riemann solver, using the algorithm described by Falle (1991).
We solved the equations in cylindrical coordinates, and we performed
the calculations on a
grid with
cm and
cm. Free flow boundary
conditions were imposed for all the variables on the upper and lower
z-boundary and the outer r-boundary. At the inner r-boundary,
mirror conditions were imposed for the r-component of the velocity
while free boundary conditions were assumed for the other variables.
Figure 2: Evolution of a photoionized globule, case A with diffuse flux 0.05 times the direct flux. Plots show the gas density and selected velocity vectors at a) , b) , c) and d) after the start of the simulation. The density scale is defined by the colour bar on each plot: the contours shown are at densities 10, 100, 1000, ... . The scale of the velocity vectors is such that would corresponds to a linear scale of | |
Open with DEXTER |
The equations for the ionization structure and radiative transfer are
treated as follows. The ionization fraction is included in the
hydrodynamic scheme as a passively advected variable, which treats the
advective term in Eq. (3). The direct and diffuse
fields are propagated across a cell
(12) |
(13) |
Then the ionization fraction is evolved via
(14) |
(15) |
(16) |
To test our scheme against previous results, we first took and followed the evolution of a globule. Our results are very similar to those given in Fig. 4 of LL94. Specifically, at time , the early stages of the formation of neutral "ears'' are clearly seen. As the "ears'' develop, they converge towards the axis of symmetry in the region shadowed from the direct radiation field. Until the "ears'' have converged on the axis at a time , the neutral gas distribution is narrowest in its most upstream direction. After convergence, the neutral gas is broadest in its most upstream direction and remains so for a period of several tenths of a My. The convergence of the "ears'' produces a globule with a dense core but which is extended in the positive (downstream) z direction with a length: width ratio of 4-5:1 at maximum compression. Re-expansion of the globule followed by recompression leads to a transient place where the globule has roughly constant radial dimensions; eventually a continuous structure with a cometary lead and long tail is produced which lasts for about before it becomes completely ionized.
In complete contrast, our results (Figs. 2 and 3) show a total suppression of the formation of "ears''. Although that region in the downstream direction where "ears'' form in the case is still shielded from the direct radiation, even a modest diffuse radiation field ionizes that small amount of gas that otherwise would have produced "ears''. Consequently the formation of a long relatively thin globule by the convergence of neutral "ears'' no longer occurs.
At the rear of the computational grid, the ionized gas is nearly static both in the simulations presented in LL94 and in the present paper. This gas is maintained close to the density in the initial conditions by the zero-gradient outer boundary condition on the grid. The shock at its edge can be taken to represent part of a termination shock in the ionized wind from the globule, which marks the transition between the wind and the intermediate pressure ionized gas which fills the overall HII region. However, since this gas remains behind the clump at all times and is separated from it by the supersonic ionized wind, it has no effect on the dynamics of the globule.
In our case A model, the compressed layer driven by ionized gas
pressure advances more quickly into the globule in surface regions
illuminated mainly by the direct radiation field than those
illuminated mainly by the diffuse field. A globule with a roughly
triangular cross-section is produced (Fig. 2a,
)
which lasts for about one tenth of a My. At a later time,
the compressed layer reflects from the z axis and neutral gas
velocities in the positive z direction are generated and the globule
elongates in this direction (Fig. 2b,
). The converging shock driven by the ionization
front has a structure analogous to a shaped charge (Batchelor 1968).
This generates a flow of neutral gas along the axis away from the
ionization source with supersonic velocities (Fig. 2c,
). The globule expands to a maximum radius at
and is then recompressed. The internal
velocities are predominantly in the positive z direction and the
globule becomes long and thin (Fig. 2d,
). The globule then undergoes further compression
and rarefaction, as most of the mass is lost by photoevaporation. In
Fig. 3 we show the state of the globule at
.
To reach this time, we extended the computational
grid in the z direction: results at earlier timesteps with this
larger grid were indistinguishable from those shown in
Fig. 2. The globule has a length about its original
size, but its width is now decreased by a factor of 5. It has been
displaced about 3 pc from its original position and little mass
remains.
Figure 3: Gas density in case A after , cf. Fig. 2 | |
Open with DEXTER |
Figure 4: Variation of neutral mass in globule with time. The upper curve is Model A, the lower curve Model B | |
Open with DEXTER |
For case B, in which the diffuse field was increased, the overall development of the structures is similar to that seen in case A. The decrease in the mass of the neutral gas in each case is shown in Fig. 4. The mass loss curves are remarkably similar, with a plateau at 0.1 My corresponding to the greatest compression of the globule, seen in Fig. 2b. In both cases, almost all the neutral gas has disappeared in the first 1 My after the radiation field switches on. The globule in case A accelerates somewhat more rapidly than that in case B, perhaps as a result of the increased diffuse field in case B decreasing the area of the globule as seen from upstream.
We conclude that ordered cometary structures lasting up to times of Mys are not produced if even a modest diffuse radiation field exists, at least for the densities in the global HII region studied in the paper, as found in the simpler geometry studied by Cantó et al. (1998). Large globules dissolve into small globules which contain very little of the original neutral mass. It remains, however, to study the stability of the neutral tails predicted by Cantó et al. (1998) for higher densities in the global HII region.
In future work we will address the photo-evaporation of magnetically supported globules. An added complication here is the modification to ionization front jump conditions produced by the magnetic fields (cf. Redman et al. 1998; Williams et al. 2000).
Acknowledgements
We thank the referee, Dr. S. Richling, for a very helpful report. RJRW thanks PPARC for support through an Advanced Fellowship.