# Three-dimensional simulations of globule and pillar formation around HII regions: turbulence and shock curvature

P. Tremblin, E. Audit, V. Minier, W. Schmidt, N. Schneider

► To cite this version:

P. Tremblin, E. Audit, V. Minier, W. Schmidt, N. Schneider. Three-dimensional simulations of globule and pillar formation around HII regions: turbulence and shock curvature. Astronomy and Astrophysics - A&A, 2012, 546, 10.1051/0004-6361/201219224. insu-03622047

HAL Id: insu-03622047

https://insu.hal.science/insu-03622047

Submitted on 28 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Three-dimensional simulations of globule and pillar formation around HII regions: turbulence and shock curvature

P. Tremblin<fnref n="1" />, E. Audit<fnref n="1" />,<fnref n="2" />, V. Minier<fnref n="1" />, W. Schmidt<fnref n="3" />, and N. Schneider<fnref n="1" />

Received 14 March 2012 / Accepted 26 July 2012

## ABSTRACT

**Aims.** We investigate the interplay between the ionization radiation from massive stars and the turbulence inside the surrounding molecular gas using three-dimensional (3D) numerical simulations.

**Methods.** We used the 3D hydrodynamical code HERACLES to model an initial turbulent medium that is ionized and heated by an ionizing source. Three different simulations were performed with different mean Mach numbers (1, 2, and 4). A non-equilibrium model for the ionization and the associated thermal processes was chosen. This turned out to be crucial when turbulent ram pressure is on the same order as the ionized-gas pressure.

**Results.** The density structures initiated by the turbulence cause local curvatures of the dense shell formed by the ionization compression. When the curvature of the shell is sufficient, the shell collapses in on itself to form a pillar, while a smaller curvature leads to the formation of dense clumps that are accelerated with the shell and therefore remain in the shell during the simulation. When the turbulent ram pressure of the cold gas is sufficient to balance the ionized-gas pressure, some dense-gas bubbles have enough kinetic energy to penetrate the ionized medium, forming cometary globules. This suggests that there is a direct relation in the observations between the presence of globules and the relative significance of the turbulence compared to the ionized-gas pressure. The probability density functions present a double peak structure when the turbulence is low relative to the ionized-gas pressure. This could be interpreted in observations as an indication of the turbulence inside molecular clouds.

**Key words.** stars: formation – HII regions – ISM: structure – turbulence – methods: numerical

## 1. Introduction

Although massive stars have a great impact on their environment, the importance of their radiative feedback on star-formation rates is still a matter of discussion. For example, Dale et al. (2005) found a slight enhancement of the star-formation activity when the radiative feedback is taken into account whereas they presented simulations on a scale of a giant molecular cloud with almost no impact (see Dale & Bonnell 2011). A closer look on small scales helps us to address this question by identifying the detailed mechanisms that form the dense structures eventually leading to star formation. These structures consist mainly of pillars of gas pointing toward the ionizing source, globules detached from the parent molecular cloud, and dense clumps at the interface between the HII region and the cloud. Optical, infrared (IR), and far-IR observations, using the Hubble Space Telescope and the Spitzer and Herschel satellites impressively revealed these features, mainly in high-mass star-forming regions (e.g. Hester et al. 1996; Gerin et al. 2009; Deharveng et al. 2010; Schneider et al. 2010; Zavagno et al. 2010). Several models have been proposed to explain their formation. The collect and collapse scenario described by Elmegreen & Lada (1977) or shadowing instabilities in the ionization front (e.g. Williams 1999) concentrate on the formation of dense clumps at the interface. Bertoldi (1989) and Lefloch & Lazareff (1994) with the radiation-driven implosion of clumps studied the formation of

globules and Mackey & Lim (2010) proposed the shadowing effect to explain the formation of pillars.

In a previous paper (Tremblin et al. 2012, Paper I hereafter), we presented a new approach emphasizing the importance of the curvature of the dense shell at the edge of the HII region to explain the formation of pillars and dense clumps in the shell. When the shell is sufficiently curved, a pillar will form as the result of the shell collapsing in on itself, while a smaller curvature will trigger lateral flows forming dense clumps and dips inside the shell. However, this scenario and the previous ones presented above are rather idealized set-ups and have to be validated in more realistic situations, e.g. by taking into account a turbulent molecular cloud.

Recent studies (see Mellema et al. 2006; Gritschneder et al. 2010; Arthur et al. 2011) started to concentrate on the interplay between the ionization coming from the massive stars and the turbulence inside the molecular cloud. They found that pillars, dense clumps, and globules emerge naturally in their models, although the detailed processes forming these structures are difficult to identify with the turbulence.

In the present study, we concentrate on the interplay between turbulence and ionization in the light of the detailed processes presented in Paper I. We first present the numerical method used and some tests of the out-of-equilibrium ionization module, and then three ionized turbulent simulations at different mean Mach numbers (1, 2, and 4). We then compare our results with the

<footnote n="1">¹ Laboratoire AIM Paris-Saclay, CEA/Irfu – Univ. Paris Diderot – CNRS/INSU, Centre d'Études de Saclay, 91191 Gif-Sur-Yvette, France
e-mail: pascal.tremblin@cea.fr</footnote>

<footnote n="2">² Maison de la Simulation, CEA-CNRS-INRIA-UPS-UVSQ, USR 3441, Centre d'Études de Saclay, 91191 Gif-Sur-Yvette, France</footnote>

<footnote n="3">³ Institut für Astrophysik der Universität Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany</footnote>
previous idealized situations presented in Paper I and provide a new approach to explaining the formation of globules. We finally show that the probability density functions of the gas can be used in observations to determine the importance of turbulence relative to the ionized-gas pressure.

## 2. Numerical methods

We used the HERACLES code<fnref n="1" /> (González et al. 2007) with the ionization module described in paper I and the cooling module described in Audit & Hennebelle (2010, based on Wolfire et al. 2003) to model the molecular cloud and the ionization coming from the OB cluster. In addition to the physics of Paper I, the turbulence in the cloud was modeled using the turbulence module described in Schmidt et al. (2006) and Schmidt et al. (2009). Gravity is an important ingredient for an ultra-compact HII region (e.g. Peters et al. 2010), although we are interested in large-scale HII regions (e.g. the Rosette molecular cloud, Schneider et al. 2010). Self-gravity and the gravity from the ionizing source were considered in Paper I without any noticeable change to the formation of the structures. We therefore do not consider gravity in the present study.

We first describe the model we used for ionization. It is an out-of-equilibrium model solving the equations

$$ 
\begin{aligned}
d(F_{\gamma})/dx &= -n_{H}(1-X)\sigma_{\gamma}F_{\gamma} + \delta_{0}(x)S_{*} \\
d(Xn_{H})/dt &= I - R = n_{H}(1-X)\sigma_{\gamma}F_{\gamma} - \beta X^{2}n_{H}^{2} \\
de/dt &= I \times e_{\gamma} - R \times k_{b}T/(\gamma - 1),
\end{aligned}
\quad (1) $$

in which X is the ionization fraction (advected with the hydrodynamic), F_γ is the photon flux, S* the source term, n_H the total hydrogen density, σ_γ the cross-section for ionization, I the ionization rate, R the recombination rate, β is given by $2 \times 10^{-10} T^{-0.75}$ cm³/s, e is the internal energy of the gas, e_γ the energy transferred from the ionizing photons to the electrons, and T the temperature of thermal equilibrium between all the species. The first equation is the ray tracing of the ionizing photons coming from the source. The second one is the ionization/recombination photo-chemistry evolution and the third one is the associated thermodynamic evolution. One-dimensional (1D) and three-dimensional (3D) tests of the out-of-equilibrium ionization and recombination processes are performed in Sect. 3.

The turbulent forcing is modeled by an Ornstein-Uhlenbeck process in Fourier space. For each mode of the force field, random increments with a Gaussian distribution are added after each time step. The wave-numbers of the forcing modes are in the interval [0, 2k₀], where k₀ is given by the typical length scale of the forcing L = 2π/k₀ and L = 4 pc is the size of the box in our case. By projecting the modes in transversal and longitudinal directions with respect to the wave vectors, divergence-free and rotation-free components are produced in physical space. The resulting force field is statistically homogeneous, isotropic, and stationary. The rms magnitude is on the order of V²/L, where V is the characteristic velocity of the flow. We perform three simulations with V resulting in a stationary transonic and supersonic velocity of Mach 1, 2, and 4. The turbulence used is solenoidal with a ratio of compressive forcing power to the total forcing power of 10%.

Three turbulent scenarios are investigated in Sect. 4, a transonic turbulent molecular cloud (mean Mach number around 1) and two supersonic cases (with mean Mach numbers of 2 and 4). The box is a cube of 4 parsec at a resolution of 0.01 pc (400³),

<figure><img src="image_2.png" /><figcaption>Fig. 1. Cut of the density field and velocity field in the radiation-driven implosion of a clump inside an ionized medium. The gas behind the clump in the white box is out of equilibrium and recombines slowly. The gas at the front of the clump escapes by the means of the rocket effect and is cooled by the strong expansion.</figcaption></figure>

the mean density is 500 cm⁻³, and the temperature is initially at 25 K. During the turbulent evolution, the boundary conditions are periodic. The simulations are run until a statistically steady state is obtained, i.e. a constant mean Mach number in the box. We then turn on the ionization with a flux of 10⁹ photons per second per squared centimeter, which is the typical flux of an O4 star at 30 parsecs. The ionization originates from the top of the box and this face is changed according to the reflexive boundary condition for the hydrodynamic, while the opposite side is set to be a free flow. The turbulent stirring is maintained during the ionization phase. In some runs, the turbulence was not maintained and there was no noticeable change, the large-scale modes do not have time to influence the small-scale ones during the ionization phase. In Sect. 4, we investigate the interplay between turbulence and ionization for the three different simulations.

## 3. Tests of the out-of-equilibrium ionization and recombination processes

When turbulence is included in an ionization simulation, the main difference is the possibility that some ionized gas gets into the shadow of the dense unionized gas because of the mixing. This is typically an out-of-equilibrium state for the ionized gas that will begin to recombine. To study this state, we investigate simplified set-ups in which hot ionized gas is put in the shadow of dense cold gas and will therefore be out of equilibrium. We simulate the impact of ionization on a clump ($D_{\text{clump}} = 0.5$ parsec) inside a homogeneous ionized medium (see Fig. 1). We also perform an “equivalent” 1D simulation where the clump is an infinite slab, to identify the effect of the geometry. The gas in the shadow of the clump/slab starts to recombine since it is out of equilibrium. The typical recombination time for the ionized gas is given by $1/\beta n$ that is on the order of $10^4-10^5$ years for ionized gas at a density of $1-10/\text{cm}^3$. It is on the order of the gas dynamical time $D_{\text{clump}}/c_{\text{io}} = 10^5$ years, where $c_{\text{io}}$ is the sound speed of the ionized gas. Therefore, an out-of-equilibrium model is needed to treat this state of the gas. We first studied the profile of the ionization fraction and the effect of the limited-variation time-stepping and then the structure of the temperature in these tests.

Our method solves the two first equations in Eq. (1) implicitly, while the thermal balance is solved explicitly and sub-cycled

<footnote n="1">¹ http://irfu.cea.fr/Projets/Site_heracles/index.html</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Results of our 1D and 3D simulations of an isolated clump inside a completely ionized medium exposed to an ionization flux of $10^9$ photons per second per squared centimeter. Note that in the 1D simulation the clump is in fact an infinite slab. **Top:** ionization profile through the middle of the box at $t = 500$ ky. **Bottom:** temperature profile at $T = 500$ ky. The ionization source is at $x = 0$ parsec. In the 1D simulation, the clump is located at $x = 1.5$ parsec, while in the 3D simulations it is at 2.75 parsec. The two 3D simulations show that we do not needed to limit too strongly the variation in the ionization fraction. The gas behind the clump is in the shadow and therefore recombine. Its recombination time is long and the gas is not in equilibrium. These simulations also show the importance of the 3D effects on the temperature structure.</figcaption></figure>

to limit the variation in $T$. The variation in $X$ is usually limited to 10% in the implicit step. However, we often do not need to limit the variations that much. The ionization fraction and the temperature profiles are plotted in Fig. 2, after 500 ky of evolution and for limited variations of 10% and 500%. There is almost no difference between the two. At the front, the equilibrium for ionization is reached even with long time-steps using the implicit method, the ionization fraction being almost at one. At the back, the recombination time of a hot ionized plasma is quite long. For example, in our case the recombination time is on the order of 100 ky which also explains why the ionization fraction is still around 0.2 after 500 ky. It is therefore possible to allow for long time-steps.

We now investigate the temperature structure. The 3D profiles in Fig. 2 are taken through the middle of the box (the white dashed line in Fig. 1), thus through the middle of the clump, which is located between 2.5 and 3 parsecs. In the 1D case, the clump is a slab located after 500 ky between 1 and 2 parsecs. The difference in the position of the two is easy to understand: in the 1D case, the ionized material cannot be evacuated on the side and the column density in front of the ionization front is higher.

There are only a few ionization processes occurring at the ionization front, where the dynamic is dominated by the expansion of the gas evaporating from the slab. In the 3D case, the ionized gas escapes on the side of the clumps, therefore ionization can penetrate further within the clump and deposit more energy. As a consequence, the clump is pushed further away by rocket motion and the energy deposit at the surface of the clump leads to the expansion of the ionized gas and the shock surface that can be seen in Fig. 1. The ionized gas is cooled to 3000–6000 K by the expansion at the front.

The temperature structure is also quite different behind the clump (white box in Fig. 1) because of the shadowed recombination. In the 1D simulation, the thermal equilibrium is reached at the front, within the HII region. Behind the slab, the gas is not in equilibrium and recombines, the ionization fraction decreasing from 1 to $\approx 0.2$. Moreover, most of the energy emitted during the recombination process is radiated away by the cooling, so that the temperature of the gas does not vary much, increasing only slightly, from $\approx 8000$ K to $\approx 9000$ K. The temperature structure is quite different in 3D as shown in Fig. 2. The cooling by recombination behind the clump should decrease the pressure of the gas as in the 1D case. However in 3D, it leads to lateral gas flows from the hot ionized gas exposed to the radiation that compress the hot recombined gas in the shadow. Therefore the recombination behind the clump is done more or less at constant pressure leading to an increase in temperature behind the clump, which is at first sight surprising.

This simple test shows that the equilibrium assumption for ionization is invalid when hot ionized gas gets in the shadow of dense unionized gas (an ionization fraction of 20% instead of 0). The equilibrium assumption is usually justified by the recombination time being much shorter than the dynamical time. This is usually true in the HII region, although this is not the case at the interface between the HII region and the cold gas when turbulence is included. At a high level of turbulence, hot ionized gas in particular will be mixed with cold dense gas and therefore be in the shadow of the cold gas, exactly as we found in this simplified test. The ionization equilibrium is not reached for a non-negligible part of the ionized gas that can get into the shadow of the dense unionized gas. This is especially the case when the level of turbulence is sufficient to balance the ionized-gas pressure and mix the unionized and ionized phases as we see in the Sect. 4.

# 4. Turbulent simulation

## 4.1. Transonic turbulence

The first line of Fig. 3 presents three snapshots of the column density for the transonic simulation, when ionization is turned on, 500 ky and 1 My after ionization has started. The hot ionized gas expands in the cold turbulent medium and triggers the propagation of a shock ahead of the ionization front. This is the collect part of the classical collect and collapse scenario (Elmegreen & Lada 1977) that is observed to be operating in many regions (see Zavagno et al. 2010; Thompson et al. 2012; Deharveng et al. 2010). The 500-ky snapshot clearly shows a three-phase medium of the hot ionized gas (at the top), the shocked dense region, and the unperturbed cold cloud. The shocked region is apparently 0.5–1 parsec wide at 500 ky and gets wider in time, up to 1.5 parsec at 1 My. This width is a lot larger than one would expect from the collect process in a homogeneous medium.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Column density snapshots of the three turbulent simulations (row 1 Mach 1, row 2 Mach 2, row 3 Mach 4) at three different times (Col. 1: before ionization starts; Col. 2: 500 ky after ionization is turned on; Col. 3: 1 My).</figcaption></figure>

The typical width of the dense shell is 0.1 parsec in this case. However, Fig. 4 clearly shows that the real thickness of the shell is indeed of order 0.1 parsec but in contrast to the collect process in a homogeneous medium, the shell is not flat but highly curved. Therefore, the column density plots show a wide shocked region owing to the projection effects of this highly turbulent shocked region.

The properties of these phases can be studied using the distribution of the mass fraction in either the density-pressure plane or the density-Mach number plane. Figures 5 and 6 show the 2D plots of the mass fraction before and after ionization is turned on. Before ionization is turned on, the distribution of the gas is at thermal equilibrium and the average Mach number is of order 1 with a dispersion of 0.48. After 1 My, the three-phase medium is formed and the three phases are clearly separate in the density-Mach number plot. The hot ionized gas can be seen on the density-pressure plot, on the isothermal equilibrium curve for the ionized gas, and on the density-Mach number plot in the

sub-sonic area. The unperturbed medium remains at the same location in the graph with respect to Fig. 5, while the shocked material is pushed to high Mach number and high densities. The distribution of the mass fraction is quite similar to the one found in Paper I, the only difference is the initial distribution of the gas, which was previously homogeneous at $500 \text{ cm}^{-3}$ with a density or interface modulation. In both cases, the shocked region contains points with densities higher than the one corresponding to the collect and collapse scenario (red dashed line in Figs. 5 and 6), i.e. the density after the collect phase in a homogeneous medium. These points represent 3% of the shell in mass fraction and help us to reach densities high enough to trigger the gravitational collapse.

## 4.2. Supersonic turbulence and comparison

The simulations with mean Mach numbers of two and four contain the same three phase structures: the ionized region, the
<figure><img src="image_5.png" /><figcaption>Fig. 4. Density cuts through the three different simulations at t = 500 ky. Left Mach 1, middle Mach 2, and right Mach 4.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Two-dimensional plot of the mass fraction at given density and pressure (top) and at given density and Mach number (bottom), before ionization is turned on in the simulation with a mean Mach number of one. Blue to black contours are increasing mass-fraction contours (blue 1E-6, green 1E-5, yellow 5E-5, orange 1E-4, red 5E-4, black 2E-3). The dashed-red line at $4.8 \times 10^4 \text{ cm}^{-3}$ is the maximum density achieved in a plane-parallel 1D simulation for a homogeneous density of 500 cm$^{-3}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Two-dimensional plot of the mass fraction at given density and pressure (top) and at given density and Mach number (bottom), 1 My after ionization is turned on. Blue to black contours are increasing mass-fraction contours (blue 1E-6, green 1E-5, yellow 5E-5, orange 1E-4, red 5E-4, black 2E-3). The dashed-red line at $4.8 \times 10^4 \text{ cm}^{-3}$ is the maximum density achieved in a plane-parallel 1D simulation for a homogeneous density of 500 cm$^{-3}$.</figcaption></figure>

shocked region, and the unperturbed region. However, in these cases the initial gas has more turbulent ram pressure to resist the ionized-gas expansion. The pressure of the ionized gas and the ram pressure caused by turbulence can be estimated from

$$p_{\text{II}} = 2\bar{n}_{\text{II}}k_{\text{b}}T_{\text{II}} \\ p_{\text{turb}} = \overline{\rho_{0}v_{0}^{2}}, \qquad (2)$$

where $n_{\text{II}}$ and $T_{\text{II}} \approx 7700$ K are the density and the temperature in the ionized gas and $\rho_0$ and $v_0$ are the density and the mean velocity in the initial turbulent medium. The density of the ionized gas is nearly constant with time in the three simulations and does not depend on the initial turbulence. It has a value of around $10 \text{ cm}^{-3}$

and the corresponding pressure at t = 500 ky is given in Table 1 with the initial turbulent ram pressure. The ionized-gas pressure dominates the turbulent ram pressure for a mean Mach number of one and two. The consequent progression of the ionized-gas expansion is clearly visible in the column density plots of Fig. 3. However, the plots also show that this progression nearly halts at a mean Mach number of four because the turbulent ram pressure is then higher than the ionized-gas pressure (see Table 1). Another interesting difference between the Mach-four case and the others is the presence of many globules. The presence of these globules in our simulations appears to be linked to the turbulent motions dominating over the ionized-gas expansion, as we investigate in more detail in Sect. 5.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Positions of the 0.95-, 0.5-, and 0.05-transition of the *xy*-average ionization fraction. Between 0 and the lower bound, the ionization fraction of the gas is between 1 and 0.95, between the lower and middle bounds, the fraction is between 0.95 and 0.5, and between the middle and upper bounds the ionization fraction is between 0.5 and 0.05.</figcaption></figure>

<table><caption>**Table 1.** Initial turbulent ram pressure and ionized-gas pressure at *t* = 500 ky.</caption><thead><tr><th>Mean Mach number</th><th>1</th><th>2</th><th>4</th></tr></thead><tbody><tr><td>p<sub>turb</sub>/kb at t = 0 ky [K/cm<sup>3</sup>]</td><td>3.2 × 10<sup>4</sup></td><td>1.2 × 10<sup>5</sup></td><td>4.5 × 10<sup>5</sup></td></tr><tr><td>p<sub>II</sub>/kb at t = 500 ky [K/cm<sup>3</sup>]</td><td>1.7 × 10<sup>5</sup></td><td>1.9 × 10<sup>5</sup></td><td>1.8 × 10<sup>5</sup></td></tr></tbody></table>

The ionization fraction changes locally from one in the ionized region to zero in the unionized region across a very small zone (smaller than 0.01 pc) when turbulence can be neglected. Therefore, the position at which the fraction is less than 0.5 gives the position of the ionization front. When turbulence is high, the position of the front is difficult to find since there is mixing between ionized gas and cold gas (see Fig. 4). However, global parameters such as the mean position and apparent width of the ionizing front can still be studied as part of the whole ionization fraction field. We considered the average of the vertical profile of the ionization fraction. It varies between one in the completely ionized region and zero in the cold region, slowly decreases in between because of the apparent width caused by the projection. We define the mean position of the ionizing front $Z_X$ as that where the average profile is at 0.5 and the width $\Delta Z_X$ as the distance between the 0.95 and 0.05 positions. This corresponds to the equations

$$ 
\begin{align} 
Z_X &= \int_0^{L_{\text{box}}} H(\langle X(x, y, z) \rangle_{x,y} - 0.5) dz, \\ 
\Delta Z_X &= \int_0^{L_{\text{box}}} H(\langle X \rangle_{x,y} - 0.05) - H(\langle X \rangle_{x,y} - 0.95) dz, \\ 
H(x-y) &= 1, \quad \text{if } x \ge y \text{ or} \nonumber \\ 
&= 0, \quad \text{if } x < y, \tag{3} 
\end{align} 
$$

where $X(x, y, z)$ is the ionization fraction field and $H$ the Heavy-side function. Both $Z_X$ and $\Delta Z_X$ are shown in Fig. 7. The front progression is decelerates at time and its width increases as the Mach number increases. This is caused by the greater density contrast with increasing turbulence. In the three cases, the width of the shell is mainly caused by the projection of the column density. The real shell is thin (see Fig. 4), but is very spatially dispersed by the inner structures of the cloud. Therefore, it appears to have a wider thickness.

<table><caption>**Table 2.** Percentage of cells out of ionization equilibrium.</caption><thead><tr><th>Time</th><th>Mach 1</th><th>Mach 2</th><th>Mach 4</th></tr></thead><tbody><tr><td>250 ky</td><td>1%</td><td>7%</td><td>8%</td></tr><tr><td>500 ky</td><td>3%</td><td>9%</td><td>12%</td></tr><tr><td>750 ky</td><td>5%</td><td>10%</td><td>13%</td></tr><tr><td>1 My</td><td>6%</td><td>12%</td><td>17%</td></tr></tbody></table>

### 4.3. Ionization and thermal equilibrium

The scheme that we used (as described in Sect. 2) is able to treat the ionization physics out of equilibrium. An interesting quantity to follow is therefore the fraction of gas that is out of equilibrium. In equilibrium, the ionization rate compensates for the recombination rate. From Eq. (1), the ionization fraction in equilibrium is $X = (\sqrt{1 + 4/y} - 1)y/2$ where $y = \sigma_{\gamma}F_{\gamma}/n_{\text{H}}\beta$. In the molecular cloud, $F_{\gamma}$ is 0 hence $X$ is 0 at equilibrium. In the HII region, $y \gg 1$, so the ionization fraction tends to 1. At the ionization front, the ionization fraction can be in equilibrium and in-between 0 and 1, although this region is very small and can be neglected in the simulation (being smaller than 1% of the box). The thickness of the ionization front is given by $1/n_{\text{H}}\sigma_{\gamma}$ which is on the order of $10^{-4}$ pc for $n_{\text{H}} = 500 \text{ cm}^{-3}$. We therefore define the fraction of gas out of equilibrium as the percentage of cells that have a ionization fraction between 0.05 and 0.95, the values at different times for the three simulations being given in Table 2. The more turbulent the medium, the more out of equilibrium. We note that this is similar to what was observed for the thermal instability in the interstellar medium (see Sánchez-Salcedo et al. 2002; Gazol et al. 2005; Audit & Hennebelle 2005). Almost one fifth of the simulation box at $t = 1$ My is out of equilibrium for a mean Mach number of four. This means that the hypothesis of an ionization equilibrium is quite inappropriate when studying the impact of ionization on a turbulent cloud. This is mainly caused by the ionized gas getting into a shadowed region. In these regions, the recombination time is long, between 1 ky and 100 ky for densities around 1.−10. cm⁻³ and a temperature around 8000 K, as we have explained in Sect. 3.

It is also important to monitor the temperature in our scheme. We solve the thermal balance between ionization and recombination which allows the gas to be out of equilibrium, as shown in Sect. 2. This behavior can be identified in Fig. 8, where the gas ahead of the dense parts facing the ionization flux is ionized but has a temperature between 2000 K and 6000 K. This is caused by the 3D expansion of the dense ionized gas at the top of the structures facing the ionization radiation. In addition, the regions that are in the shadow of the dense structures tend to be overheated at almost 10 000 K while the equilibrium temperature is around 7750 K. This is because the recombination process in the shadow operates at constant pressure, the pressure being imposed by the hot ionized and exposed gas surrounding the shadowed gas. Therefore, as we have demonstrated in Sect. 2 for the simple situation of a clump shadowing parts of the hot ionized gas, the temperature behind the clump will increase. In the turbulent simulation, the cooling processes for the unionized gas also play a role (see Paper I for details), cooling the recombined gas so that the temperature does not reach 14 000 K but equilibrates around 10 000 K. This result is counter-intuitive since one would expect the gas to stay at 7750 K during the recombination step, the kinetic energy of the electrons being radiatively lost from the system. However, the 3D versus 1D study done in Sect. 2 clearly show the importance of the lateral motions from the exposed gas
<figure><img src="image_9.png" /></figure>

**Fig. 8.** Temperature cut at *t* ≈ 710 ky after ionization is turned on. The 3D temperature effects studied in the clump simulation (Sect. 2) are clearly visible through the snapshot.

to the shadowed gas in achieving a recombination at constant
pressure behind the structures.

**5. Structures and observational signatures**

**5.1. Pillars**

It is unclear how to distinguish between the structures in the col-
umn density maps of Fig. 3 because of the confusion along the
line of sight. We therefore investigate some density cuts to help
us identify pillars and globules. Figure 9 presents the formation
of a pillar in the Mach 1 simulation between 200 ky and 700 ky.
The pillar formed is one parsec long and has a complex struc-
ture, with successive dense parts along its length. Figure 10 is a
map along the pillar of the density and the direction of the veloc-
ity field. The structure of the field is similar to the one already
obtained in Paper I for non-turbulent set-ups. The pillar has two
dense heads at vertical positions of 1.4 pc and 1.6 parsec, a dense
base between 1.8 pc and 2 pc and two holes on the side at 2.4 pc.
The velocity field shows that the pillar forms at its base by the
convergence of the dense parts of the shell, while the holes are
formed by the divergence of the dense parts created in a concave
zone. This is the same as the configuration previously seen in
Paper I, without any turbulence. This suggests that the mecha-
nism at the origin of the formation of pillars does not depend
too much on the turbulence. Here the turbulence gives only the
initial conditions in terms of the density contrasts and the initial
structures on which pillars will form. The holes on the side of the
pillars in the non-turbulent set-ups of Paper I were in fact part of
a complete annulus around the pillars caused by the symmetry
of our initial conditions. The major difference with a turbulent
scenario is that there is no symmetry in the initial conditions.
This leads to holes that are local and do not extend to form a full
annulus.

The curvature of the shell can also be seen in spectra of the
line-of-sight velocity. Figures 11 and 12 show on the left the
column density maps of the pillar and on the right the mass
histograms as a function of the line-of-sight velocity in squares
of 0.05 parsec. Figure 11 is a snapshot at *t* ≈ 240 ky when the
pillar did not have yet formed and Fig. 12 at *t* ≈ 710 ky for
the whole pillar. The line-of-sight velocity spectra are wide (be-
tween −2 and +2 km s⁻¹) before the pillar is formed. This was al-
ready identified in Paper I as the signature of the dense shell that
is curved in on itself with the two components, blue-shifted and

<figure><img src="image_10.png" /></figure>

**Fig. 9.** Density cut at *t* ≈ 240 ky and *t* ≈ 710 ky after ionization is turned on. The white box indicates the areas of the column density and line-of-sight velocity spectra made in Figs. 11 and 12. The arrow indicates the direction of the line of sight; for the cross, the line of sight is perpendicular to the density cut plane.

red-shifted, that are going to collide to form the pillar. When the
collision occurred, the wide spectra were no longer visible and
the line-of-sight velocity spectra peaked around a null velocity
(see Fig. 12). This demonstrates that the scenario considered in
Paper I for the study of non-turbulent set-ups is also applicable
to our turbulent simulation.

On the basis of this analysis, we can infer that molecular
tracers that are optically thin should reflect the same spectral
structure as the line-of-sight velocity spectra of the simulations.
Therefore, wide spectra could be an indication of nascent pillars,
whereas evolved pillars should have a spectra peaked around
the velocity of the expanding bubble projected along the line of
sight. This velocity is zero in our case since the line of sight is
perpendicular to the direction of the expansion.

**5.2. Dense clumps at the edge of the ionized gas**

The difference between the formation mechanism of pillars and
clumps inside the shell is also clearly visible in Fig. 9. In the first
snapshot, the position at which the pillar is forming can be iden-
tified as the position at which the curvature of the shell is the
highest. In the other parts of the shell, the curvature is insuffi-
cient to form pillars but instead triggers the formation of clumps
inside the shell. This phenomenon was already seen in Paper I,
where the higher the curvature, the longer the pillar. If the shell
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Cut of the density field and velocity-field orientation along the pillar at *t* ≈ 710 ky.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Left: column density in a 0.6 × 0.5 × 0.5 pc³ box around the area in which the pillar will form at *t* ≈ 240 ky. Right: mass-weighted histogram of the line-of-sight velocity in the same box (similar to optically thin observational line spectra). Each spectrum is made on a square of 0.1 × 0.1 pc² drawn on the column density map. The spectra are drawn between −4 and 4 km s⁻¹ in 80 bins (horizontal axis) and for mass between 10⁴ and 1 solar masses (vertical axis in log scale). The lateral shocks can be identified in the wide spectra leading to a very broad line width.</figcaption></figure>

is curved enough to collapse in on itself, it will form a pillar. At high curvature, the tip of the curved shell is stopped by the collapse of the matter ahead of it, which leads to the formation of a pillar (as shown in the right part of Fig. 13, at early times the velocity field of the shocked gas shows that the shell collapses on the center of the initial structure). The shell around the “hill” will collapse quickly on the hill, with velocities that are nearly perpendicular to the direction of the expansion. This is why the motion of the gas at the top of the hill is small. When the curvature is low, the shell collapse in on the hill but with a velocity that has a component in the direction of the expansion. Therefore after the collapse, the shocked gas maintains its velocity in the direction of the expansion and propagates with the rest of the shell.

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Left: column density in a 0.7 × 1.4 × 0.5 pc³ box around the pillar structure at *t* ≈ 710 ky. Right: mass-weighted histogram of the line-of-sight velocity in the same box (similar to the optically thin observational line spectra). Each spectrum is made on a square of 0.1 × 0.1 pc² drawn on the column density map. The spectra are drawn between −4 and 4 km s⁻¹ in 80 bins (horizontal axis) and for mass between 10⁴ and 1 solar masses (vertical axis in log scale). The lateral shocks have collided to form the pillar, the spectra now being peaked around 0 km s⁻¹.</figcaption></figure>

The matter is accelerated with the shell and the density increases because the gas surrounding it converges smoothly at the center (left part of Fig. 13). This situation was studied and compared in the context of RCW36 in the Vela C molecular cloud (see Minier et al., in prep.). The new observations from the Herschel space telescope show the presence of dense clumps around the HII region RCW36. Dedicated numerical simulations show that it is rather unlikely that these clumps pre-existed because they would have triggered the formation of pillars. These dense clumps are located at the edge of the HII region, hence are formed by the lateral flows in the shell caused by curvature perturbations that are not high enough to form pillars. It is exactly the same process that is at work in the present turbulent simulations. From Paper I, it can be estimated that the transition from clump to pillar formation is at a curvature radius of between 2.5 × 10⁻² pc and 5 × 10⁻² pc in the present situation.

Thompson et al. (2012) showed that the dense clumps around HII regions are preferentially located at the interface between the ionized gas and the cold gas. An example of such dense clumps at the border of the ionized gas can also be found in Tothill et al. (2002). We propose here a formation mechanism for these dense clumps in the shell that is based on curvature perturbations. However, if the turbulent ram pressure is far greater than the ionized gas pressure, these dense clumps will not be present exactly at the interface, the cold gas will have sufficient kinetic energy to disrupt the interface, and a mixing region between the dense shocked gas and the ionized gas will be present (see Fig. 4). Therefore, the clear correlation in the observations between the position of the dense clumps and the interface suggests that in most of these regions the ionized-gas pressure dominates the turbulent ram pressure of the surrounding gas.
<figure><img src="image_14.png" /></figure>

**Fig. 13.** Illustration of the effect of the shell curvature on the formation of clumps and pillars. Density cut and velocity field (for the unionized gas) are taken from the simulations of Paper I. Initial curvature radius: left $5 \times 10^{-2}$ pc and right $1 \times 10^{-2}$ pc. *Top*: 20 ky after the beginning of the simulation, the shell is formed on the curved surface. *Bottom*: snapshots after 550 ky. A 1.5-pc long pillar is formed at high curvature, whereas a clump of the size of the initial structure is formed at low curvature.

### 5.3. Globules

In contrast to the dense clumps in the shell, globules are bubbles of gas disconnected from the molecular cloud and surrounded by ionized gas. The major difference between the turbulent simulations presented in this paper and the non-turbulent set-ups studied in Paper I is the appearance of these globules in the highly turbulent case. As we have already discussed above, globules emerge because the ram pressure of turbulence is far higher than the pressure of the ionized gas. These bubbles of cold dense gas have enough kinetic energy to penetrate within the HII region. An interesting consequence is that the motion of the globules is imposed by the turbulence. These globules could have by chance a motion aligned with the direction of the expansion of the HII region but most of the time, they should have a random motion direction set by the turbulence, which is indeed the case in our simulation. The direction of expansion of the HII region is from top to bottom in Fig. 3. Therefore, when projected along the line of sight that is perpendicular, the velocity in the expansion direction is zero. The motion of a pillar is typically aligned with the direction of the expansion of the HII region, as can be seen in Fig. 12. The spectra are centered on zero. However, this is not the case for a globule. Figure 14 shows a typical globule in the simulation at Mach four. The line-of-sight velocity spectra are systematically shifted by a velocity of +4 km s⁻¹, which is on the order of ten times the sound speed of the cold gas. The shift is also highlighted in Fig. 15. The globule does not have a motion that is aligned with the direction of the expansion. In a region where globules and pillars are present, the velocity of a pillar is the expansion velocity of the front projected along the line of sight, whereas the velocity of a globule is the signature of the initial turbulence. The dense clumps that are observed in the shell also move with the expansion of the HII region. The same velocity shift between these clumps and the globules nearby can

<figure><img src="image_15.png" /></figure>

**Fig. 14.** Left: column density on a $0.3 \times 0.3 \times 0.25$ pc³ box around a globule in the Mach 4 simulation. Right: mass-weighted histogram of the line-of-sight velocity in the same box (similar to optically thin observational line spectra). Each spectrum is made on a square of $0.05 \times 0.05$ pc² drawn on the column density map. The spectra are drawn between -6 and 6 km s⁻¹ in 80 bins (horizontal axis) and for mass between $10^4$ and 1 solar masses (vertical axis in log scale). All the spectra are red-shifted by a velocity of +4 km s⁻¹.

<figure><img src="image_16.png" /></figure>

**Fig. 15.** Mass fraction of the gas in the pillar and the globule as a function of the velocity. The bulk line-of-sight velocity of the pillar is at 0 km s⁻¹, whereas the bulk velocity of the globule is at +4 km s⁻¹.

be thus expected. In this picture, pillars and clumps at the edge of the HII region are the result of a structure dominated by the ionization dynamic whereas globules are turbulent dominated. This is quite different from the radiation-driven implosion scenario, in which there is no reason for the globules to have a velocity that is different from those of the pillars.

This signature of the turbulent nature of a globule is supported by recent observations in the Cygnus X region. Thanks to Herschel Open Time and SOFIA observations (Schneider et al. 2012b), 138 µm and CO maps are available and contain several pillars and globules. A pillar and a globule in this region were identified to be sufficiently close (and the HII region sufficiently big) for us to be able to assume that they should have the same direction of motion if the motion is only imposed by the HII region expansion. However the [CII] spectra show a line-of-sight velocity difference of two kilometers per second, which is clearly consistent with a turbulent scenario for the formation of the globule against the radiation implosion of an isolated clump.

The difference of $2$ km s⁻¹ is therefore an indication of the level of initial turbulence. For a sonic velocity of $0.4$ km s⁻¹, this shift typically indicates a motion at Mach 5. In our simulation, the globule is at Mach 10, whereas the mean Mach number of the simulation is 4. Therefore the velocity of the globule is an indication of the turbulence level but is not a precise measurement.
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Probability density function of the gas for the three simulations. Left *t* = 0 ky, middle *t* = 500 ky, right *t* = 1 My. The first peak is the ionized gas and therefore not seen in the observations. At *t* = 500 ky, the Mach-1 simulation presents a double peak, the first one corresponding to the unperturbed turbulent gas, the second one to the shocked region. In the Mach-2 simulations, the peak is barely visible because the turbulence is higher and the unperturbed region still present in the box becomes small. In the Mach-4 simulation, the turbulence is too high, and the second peak is hidden in the peak of the turbulent cold region.</figcaption></figure>

5.4. *Probability density function*

We investigate the density structure of the gas using probabil-
ity density functions (PDFs). Numerical models (e.g. Kritsuk
et al. 2007; Audit & Hennebelle 2010; Ballesteros-Paredes
et al. 2011) and observations of atomic and molecular gas
(e.g. Kainulainen et al. 2009) show that PDFs are log-normal
at low densities and can have more complex shapes at higher
densities. We here plot the mass fraction against the density
(Fig. 16) from our simulations at three different times find-
ing that the PDF shape depends on both the amount of turbu-
lence and time. At low turbulence (Mach 1), the probability den-
sity function is double-peaked because of the ionization, sim-
ilarly to the bimodal PDFs caused by the thermal instability
(see Sánchez-Salcedo et al. 2002; Gazol et al. 2005; Audit &
Hennebelle 2005). The dense compressed gas forms a new peak
in the distribution. At high turbulence, the probability density
function is shifted toward higher densities but does not produce
the double peak. This directly indicates that a double peak is the
signature of a region dominated by ionization, whereas a single
peak is the signature that turbulence predominates. This signa-
ture was already observed in the Rosette nebula by Schneider
et al. (2012a). The PDFs at the edge of the HII regions display
the same double-peaked signature. This suggests that at the bor-
der of the bubble, the turbulence is low compared to the ionized-
gas pressure.

On the basis of the PDFs of the Mach 2 simulation, it is clear
that the double peak disappears when there is no more unper-
turbed gas in the region for which the PDFs were compiled.
Therefore a good way of studying the PDFs around a HII re-
gion is to analyze concentric PDFs centered on the central ion-
izing cluster. By increasing the radius of the region, we first cap-
ture the peak of the shocked region, then a double peak should
appear when the unperturbed region is included. The first peak
should then disappear at some point, when the shocked region
becomes insignificant compared to the unperturbed region. This
point could indicate the scale on which the central cluster has an
impact on the global PDF of the cloud and therefore on which
scale the initial mass function is affected by the ionization. If the
turbulence is high, it is possible that the shocked peak is hidden
in the large peak of the turbulent cloud. In that case, no double
peak will appear and this will be an indication that the turbulence
in the region dominates over the ionized-gas pressure.

6. Conclusions

Structures at the interface between HII regions and molecu-
lar clouds can be classified into three main categories: pillars,

globules, and dense clumps. Various scenarios have been inves-
tigated in the past to explain their formation: collect and col-
lapse, radiation-driven implosion, shadowing effects, and turbu-
lence. We have presented a new model to explain the formation
of clumps and pillars in Paper I based on the curvature of the
dense shell formed by the collect processes. High curvature leads
to pillars, low curvature to instabilities forming dense clumps
and dips in the shell. This model was investigated using non-
turbulent set-ups. In the present paper, we have shown that the
same mechanisms are at work in simulations with a turbulent
cloud. In particular, the same spectral observational signatures
can be identified in the formation stages of the pillars.

Furthermore, we have shown that:

- Because of turbulence, hot ionized gas can get into the shadow of cold dense gas, leading to a very long recombination time for the gas. This gas can only be treated with a non-equilibrium model for ionization, where up to 20% of the box can get into that state. The equilibrium assumption is valid in situations where the ram pressure of turbulence is less significant than the ionized-gas pressure, i.e. low turbulence levels.

- Globules form preferentially when the turbulence in the cold gas predominates over the ionized-gas pressure. Bubbles of dense gas have sufficient kinetic energy to penetrate into the HII regions forming the globules. A signature of this scenario is the line-of-sight velocity spectrum of the globules, which is either blue-shifted or red-shifted relative to the spectrum of a pillar or a clump at the edge of the HII region. This can be directly observed in regions where pillars, clumps and globules are present.

- Probability density functions are double peaked when the turbulent ram pressure is low compared to the ionized-gas pressure. This result could be used in observations to establish the relative importance of the turbulence and the ionization.

An important step would be to analyze the statistics of the for-
mation of dense structures at the edge of HII regions. It will
be of great importance to assess the potential negative or pos-
itive effects of radiative feedback on the star-formation rates.
Furthermore, comparisons with observations will also help us
to determine the likelihood of the curvature scenario for the for-
mation of dense clumps and pillars and the turbulent scenario for
the formation of globules.

*Acknowledgements.* This work was granted access to the HPC resources of [CCRT/CINES/IDRIS] under the allocations c2011042204 and c2012042204 made by GENCI (Grand Équipement National de Calcul Intensif).
References

Arthur, S. J., Henney, W. J., Mellema, G., De Colle, F., & Vázquez-Semadeni, E. 2011, MNRAS, 414, 1747

Audit, E., & Hennebelle, P. 2005, A&A, 433, 1

Audit, E., & Hennebelle, P. 2010, A&A, 511, A76

Ballesteros-Paredes, J., Vázquez-Semadeni, E., Gazol, A., et al. 2011, MNRAS, 416, 1436

Bertoldi, F. 1989, ApJ, 346, 735

Dale, J. E., & Bonnell, I. 2011, MNRAS, 414, 321

Dale, J. E., Bonnell, I. A., Clarke, C. J., & Bate, M. R. 2005, MNRAS, 358, 291

Deharveng, L., Schuller, F., Anderson, L. D., et al. 2010, A&A, 523, A6

Elmegreen, B. G., & Lada, C. J. 1977, ApJ, 214, 725

Gazol, A., Vázquez-Semadeni, E., & Kim, J. 2005, ApJ, 630, 911

Gerin, M., Goicoechea, J. R., Pety, J., & Hily-Blant, P. 2009, A&A, 494, 977

González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429

Gritschneder, M., Burkert, A., Naab, T., & Walch, S. 2010, ApJ, 723, 971

Hester, J. J., Scowen, P. A., Sankrit, R., et al. 1996, AJ, 111, 2349

Kainulainen, J., Beuther, H., Henning, T., & Plume, R. 2009, A&A, 508, L35

Kritsuk, A. G., Norman, M. L., Padoan, P., & Wagner, R. 2007, ApJ, 665, 416

Lefloch, B., & Lazareff, B. 1994, A&A, 289, 559

Mackey, J., & Lim, A. J. 2010, MNRAS, 403, 714

Mellema, G., Arthur, S. J., Henney, W. J., Iliev, I. T., & Shapiro, P. R. 2006, ApJ, 647, 397

Peters, T., Mac Low, M.-M., Banerjee, R., Klessen, R. S., & Dullemond, C. P. 2010, ApJ, 719, 831

Sánchez-Salcedo, F. J., Vázquez-Semadeni, E., & Gazol, A. 2002, ApJ, 577, 768

Schmidt, W., Hillebrandt, W., & Niemeyer, J. 2006, Computers & Fluids, 35, 353

Schmidt, W., Federrath, C., Hupp, M., Kern, S., & Niemeyer, J. C. 2009, A&A, 494, 127

Schneider, N., Motte, F., Bontemps, S., et al. 2010, A&A, 518, L83

Schneider, N., Csengeri, T., Hennemann, M., et al. 2012a, A&A, 540, L11

Schneider, N., Güsten, R., Tremblin, P., et al. 2012b, A&A, 542, L18

Thompson, M. A., Urquhart, J. S., Moore, T. J. T., & Morgan, L. K. 2012, MNRAS, 421, 408

Tothill, N. F. H., White, G. J., Matthews, H. E., et al. 2002, ApJ, 580, 285

Tremblin, P., Audit, E., Minier, V., & Schneider, N. 2012, A&A, 538, A31

Williams, R. J. R. 1999, MNRAS, 310, 789

Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278

Zavagno, A., Russeil, D., Motte, F., et al. 2010, A&A, 518, L81