A&A 369, 263-268 (2001)
DOI: 10.1051/0004-6361:20010070

The modification by diffuse radiation of "cometary tail'' formation behind globules

K. G. Pavlakis 1,2 - R. J. R. Williams3 - J. E. Dyson 2 - S. A. E. G. Falle4 - 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


1 Introduction

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.

2 The model

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:

\begin{displaymath}%
\frac{ \partial \rho }{ \partial t } + \nabla\cdot\left(\rho {\vec v}
\right) = 0
\end{displaymath} (1)

and

\begin{displaymath}%
\frac{\partial \left(\rho {\vec v} \right)} {\partial t}
+ \nabla\cdot\left(
\rho {\vec v} {\vec v} \right) =
-\nabla p
\end{displaymath} (2)

where $\rho, {\vec v}, p$, and t are the mass density, velocity, pressure, and time respectively. We calculate the fractional ionization, $x_{\rm e}$ by solving the ionization equation

 \begin{displaymath}\frac{\partial \left( \rho x_{\rm e} \right)}{\partial t}
+ ...
...m e} \right)
-\alpha_{\rm B} n_{\rm H} x^{2}_{\rm e} \right]
\end{displaymath} (3)

where $n_{\rm H} = \rho/m_{\rm H}$ is the density of hydrogen nucleons, $m_{\rm H}$ is the proton mass, a is the mean photoionization cross section, c is the speed of light, E is the local Lyman continuum photon density and $\alpha_{\rm B}$ is the case B recombination coefficient (although we include the diffuse continuum incident from outside the grid explicitly, we treat recombinations within the numerical grid in the on-the-spot approximation).

We assume that

\begin{displaymath}%
E = E_{\rm dir} + E_{\rm dif},
\end{displaymath} (4)

where $E_{\rm dir}$ and $E _{\rm dif}$ are the direct and diffuse components of E. We assume that the direct radiation field is plane parallel, so $E_{\rm dir} = F_{\rm dir}/c$ and the direct radiation flux $F_{\rm dir}$ satisfies

 \begin{displaymath}%
\frac{{\rm d} F_{\rm dir}}{{\rm d} z} =
- a n_{\rm H} \left( 1 - x_{\rm e} \right) F_{\rm dir},
\end{displaymath} (5)

where z measures the distance of a point (in the direction of the globule) away from the plane perpendicular to the symmetry axis and containing the star.
  \begin{figure}
\par\subfigure[]{\includegraphics[angle=-90,width=8cm,clip]{h2485...
...ure[]{\includegraphics[angle=-90,width=8cm,clip]{h2485f1b.ps} }
\par\end{figure} 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 $F_{\rm dif} =
\pi I_{\rm dif} = cE_{\rm dif}/2$[*] with distance from the axis, r, by integrating

 \begin{displaymath}%
\frac{{\rm d} \left(F_{\rm dif} \right)} {{\rm d}r} =
2a n_{\rm H} \left(1 - x_{\rm e} \right)F_{\rm dif},
\end{displaymath} (6)

from the edge of the grid towards the symmetry axis. The doubling of the optical depth compared to the plane parallel case results from the larger path-length through which the more oblique photons must travel.

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, $r_{\rm i}$, where the inward-going radiation flux reaches zero. This radius is shown as a function of the parameters

\begin{displaymath}%
\xi = {r_{\rm s} n_{\rm i}^2 \alpha_{\rm B}\over F_{\rm dif}}
\end{displaymath} (7)

for the uniform density case calculated by Cantó et al., and

\begin{displaymath}%
\eta = {r_{\rm s} n_{\rm i}^2 r_{\rm i}^4 \alpha_{\rm B}\over F_{\rm dif}}
\end{displaymath} (8)

for a 1/r density distribution more characteristic of a cylindrical evaporative flow. In these expressions, $r_{\rm s}$ is the shadow radius; $r_{\rm i}$, the radius of the neutral core and $n_{\rm i}$ is the density in the ionized gas at the edge of the tail. The additional factor of $r_{\rm i}^4$ in $\eta$ results from the fundamental parameter for a photoevaporating column being its mass per unit length of neutral material ( $\propto n_{\rm i} r_{\rm i}^2$). For fixed neutral mass per unit length, the radiation field formally cannot destroy the tail, but its lifetime while being photoevaporated would decrease in proportion to $r_{\rm i}$.

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 $\xi$ or $\eta$). 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 $x_{\rm e}$ taken by LL94. That is, for densities greater than 20 per cent of the initial cloud density:

 \begin{displaymath}%
T^{-1} \!=\! \left[ \left( 2 T_{\rm i} \right)^{-1}
\!+\! ...
...rt{x_{\rm e}} \right)^{2} \right]
\left[1 + x_{\rm e} \right]
\end{displaymath} (9)

with $T_{\rm i}$ the temperature of fully ionized gas and $T_{\rm n}$that of fully neutral gas. For densities less than 20 per cent of the initial density of the cloud, the gas was assumed to be warm neutral material, at a temperature $T_{\rm i}$ independent of its ionization state. We used $T_{\rm i} = 10^{4}$ K and $T_{\rm n} =
10^{2}$ K.

To enable us to compare our results as closely with those of LL94 as possible, we assumed, as they did, that $a = 3~10^{-18}{\rm\,cm^{2}}$, 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, $\alpha_{\rm B} = 2.7~10^{-13}{\rm\,cm^3\,s^{-1}}$, 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 $x_{\rm e}$, the gas has temperature T given by Eq. (9). Within a volume V, there are $(1-x_{\rm e})nV$ neutral hydrogen atoms and $x_{\rm e}nV$hydrogen ions, respectively, where n is the hydrogen nucleon density. Now interpret $x_{\rm e}$ 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 $T_{\rm n}$, density $n_{\rm n}$ and occupies volume (1-fx)V, while the hot ionized phase has temperature $T_{\rm i}$, density $n_{\rm i}$ and occupies volume fx V. Conservation of mass demands $(1-x_{\rm e})n =
(1-f_x) n_{\rm n}$ and $x_{\rm e} n = f_x n_{\rm i}$. Since there is pressure equilibrium $n_{\rm i}/n_{\rm n} \equiv \epsilon \ll 1$. Hence the total recombination rate within V is then

 \begin{displaymath}%
\alpha_{\rm B} n_{\rm i}^2 f_x V = \alpha_{\rm B} n^2 x_{\rm e}^2 V \left(1+
\epsilon {1-x_{\rm e}\over x_{\rm e}}\right),
\end{displaymath} (10)

where the left-hand side is explicitly the recombination rate in the fully-ionized part of the gas, and $\alpha_{\rm B}$ is the recombination coefficient at the equilibrium temperature of the ionized gas. If, as in the present paper, the temperature variation of the recombination coefficient is ignored then assuming a recombination rate $\alpha_{\rm B} n^2 x_{\rm e}^2$ gives accurate results near poorly resolved IFs (except, as a result of the term in brackets on the r.h.s. of Eq. (10), close to $x_{\rm e}=0$ where the recombination rate is small), while if the temperature variation of $\alpha_{\rm B}$ is included, the recombination rate may be substantially overestimated.

In the initial configuration, a uniform sphere of gas with a radius of 0.5 pc centred at $z = 0{\rm\,pc}$, $T = 100{\rm\,K}$ and $n_{\rm H} = 1.55~10^{3}$ cm-3 is surrounded by initially neutral gas with $n_{\rm H} = 15.5$ cm-3. At t = 0, the radiation field switches on, and thereafter the flux of Lyman continuum photons is $F_{\rm dir} \left(z = -0.75{\rm\,pc} \right) =
2.18~10^{9}{\rm\,cm^{-2}\,s^{-1}}$. Thus, our initial conditions are the same as those of LL94 (model 2) except for our inclusion of $F_{\rm dif}$ which we took to have a constant value at $r = 0.6{\rm\,pc}$. For our model A, $F_{\rm dif} \left( r = 0.6
{\rm\,pc} \right) = 0.05 F_{\rm dir} \left( z = -0.75 {\rm\,pc}
\right)$, and for our model B, $F_{\rm dif} \left(r = 0.6 {\rm\,pc}
\right) = 0.15 F_{\rm dir} \left( z = -0.75 {\rm\,pc} \right)$.

For these values, the initial values of $\xi$ 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 $\xi$ would be substantially greater than unity for a long thin tail, since

\begin{displaymath}%
\xi = {r_{\rm s} n_{\rm i}^2 \alpha_{\rm B}\over F_{\rm dif...
...{F_{\rm dir}\over F_{\rm dif}} {r_{\rm s}\over\ell_{\rm Str}},
\end{displaymath} (11)

where $\ell_{\rm str}$ is the Strömgren length due to the direct field (certainly an upper limit on the length of the tail). For a long, thin tail to form $\ell_{\rm Str}\gg r_{\rm s}$, and so the direct flux must be very much greater than the diffuse flux for $\xi\mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil
$\displaystyle...
...r{\offinterlineskip\halign{\hfil$\scriptscriptstyle .... This is not the case on average through an ionized nebula, although it may occur close to the ionizing star (Rubin 1968).

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 $64 \times 190$ grid with $\Delta z = 4.06\
10^{16}$ cm and $\Delta r = 2.89 \ 10^{16}$ 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.

  \begin{figure}
\par\includegraphics[width=16.5cm,clip]{h2485f2.ps}\end{figure} 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) $t=0.088 {\rm\,My}$, b) $t=0.158{\rm\,My}$, c) $t=0.200 {\rm\,My}$ and d) $t=0.327{\rm\,My}$ 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, ... ${\rm\,cm^{-3}}$. The scale of the velocity vectors is such that $45{\rm\,km\,s^{-1}}$ would corresponds to a linear scale of $0.4{\rm\,pc}$
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

\begin{displaymath}%
F_{\rm dir} \left( z + \Delta z \right) =
F_{\rm dir} \lef...
...t) {\rm e}^{ -a n_{\rm H} \left(
1-x_{\rm e}\right) \Delta z}
\end{displaymath} (12)

and

\begin{displaymath}%
F_{\rm dif} \left( r - \Delta r \right) =
F_{\rm dif} \lef...
...) {\rm e}^{ -2a n_{\rm H}
\left(1-x_{\rm e}\right) \Delta r},
\end{displaymath} (13)

cf. Eqs. (5) and (6). The diffuse field is switched on only when the main ionization front has reached the large-z edge of the numerical grid, as no diffuse radiation would be expected until this time.

Then the ionization fraction is evolved via

\begin{displaymath}%
x_{\rm e} \left( t + \Delta t \right) = f x_{\rm e}(t) + (1-f) x_{\rm eq}
\end{displaymath} (14)

with

\begin{displaymath}%
f = \left[1 +
{
\left\{acE\left[1-x_{\rm e}(t)\right]-\alp...
...\over x_{\rm eq} - x_{\rm e} \left( t \right)
}
\right]^{-1}.
\end{displaymath} (15)

This agrees with the ionization part of Eq. (3) to first order in $\Delta t$; for large $\Delta t$ the ionization reaches the equilibrium value, $x_{\rm e} = x_{\rm eq}$, which satisfies

\begin{displaymath}%
a c E(1-x_{\rm eq}) = \alpha_B n_{\rm H} x_{\rm eq}^2.
\end{displaymath} (16)

3 Results and discussion

To test our scheme against previous results, we first took $F_{\rm
dif} = 0$ and followed the evolution of a globule. Our results are very similar to those given in Fig. 4 of LL94. Specifically, at time $t \simeq 0.1{\rm\,My}$, 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 $t
\approx 0.2{\rm\, My}$, 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 $1.4{\rm\,My}$ 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 $F_{\rm
dif} = 0$ 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, $t=0.088 {\rm\,My}$) 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, $t=0.158{\rm\,My}$). 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, $t=0.200 {\rm\,My}$). The globule expands to a maximum radius at $t
\approx 0.23{\rm\,My}$ and is then recompressed. The internal velocities are predominantly in the positive z direction and the globule becomes long and thin (Fig. 2d, $t=0.327{\rm\,My}$). 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 $t=0.790{\rm\,My}$. 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.

  \begin{figure}
\par\includegraphics[width=4.5cm,clip]{h2485f3.ps}\par\end{figure} Figure 3: Gas density in case A after $t=0.790{\rm\,My}$, cf. Fig. 2
Open with DEXTER


  \begin{figure}
\par\includegraphics[angle=-90,width=8cm,clip]{h2485f4.ps}\par\end{figure} 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.

References

 
Copyright ESO 2001