DETAILED NUMERICAL SIMULATIONS ON THE FORMATION OF PILLARS AROUND H ii REGIONS

, , , and

Published 2010 October 18 © 2010. The American Astronomical Society. All rights reserved.
, , Citation Matthias Gritschneder et al 2010 ApJ 723 971 DOI 10.1088/0004-637X/723/2/971

0004-637X/723/2/971

ABSTRACT

We study the structural evolution of turbulent molecular clouds under the influence of ionizing radiation emitted from a nearby massive star by performing a high-resolution parameter study with the iVINE code. The temperature is taken to be 10 K or 100 K, the mean number density is either 100 cm−3 or 300 cm−3. Furthermore, the turbulence is varied between Mach 1.5 and Mach 12.5, the main driving scale of the turbulence is varied between 1 pc and 8 pc. We vary the ionizing flux by an order of magnitude, corresponding to allowing between 0.5% and 5% of the mass in the domain to be ionized immediately. In our simulations, the ionizing radiation enhances the initial turbulent density distribution and thus leads to the formation of pillar-like structures observed adjacent to H ii regions in a natural way. Gravitational collapse occurs regularly at the tips of the structures. We find a clear correlation between the initial state of the turbulent cold cloud and the final morphology and physical properties of the structures formed. The most favorable regime for the formation of pillars is Mach 4–10. Structures and therefore stars only form if the initial density contrast between the high-density unionized gas and the gas that is going to be ionized is lower than the temperature contrast between the hot and the cold gas. The density of the resulting pillars is determined by a pressure equilibrium between the hot and the cold gas. A thorough analysis of the simulations shows that the complex kinematical and geometrical structure of the formed elongated filaments reflects that of observed pillars to an impressive level of detail. In addition, we find that the observed line-of-sight velocities allow for a distinct determination of different formation mechanisms. Comparing the current simulations to previous results and recent observations, we conclude that, e.g., the pillars of creation in M16 formed by the mechanism proposed here and not by the radiation driven implosion of pre-existing clumps.

Export citation and abstract BibTeX RIS

1. INTRODUCTION

Stars are known to form in turbulent, cold molecular clouds. When massive stars ignite, their UV radiation ionizes and heats the surrounding gas, leading to an expanding H ii bubble. As soon as the H ii region breaks through the surface of the molecular cloud, a low-density, optically thin hole is formed, which reveals the otherwise obscured interior. At the interface between the H ii region and the molecular gas, peculiar structures, often called pillars, are found. The most famous examples are the "pillars of creation" in the Eagle Nebula (M16; Hester et al. 1996). There is also wide-spread evidence for star formation at the tips of the pillars (e.g., Sugitani et al. 2002, 2007; Thompson et al. 2002; Oliveira 2008). Since the launch of the Spitzer Space Telescope a wealth of highly resolved observations of the peculiar, pillar- or trunk-like structures observed around the hot, ionized H ii regions around massive stars, and the star formation in these trunks, has become available, e.g., in the Orion clouds (Stanke et al. 2002; Lee & Chen 2007; Bowler et al. 2009), the Carina nebula (Smith et al. 2000), the Elephant Trunk Nebula (Reach et al. 2004), the Trifid Nebula (Lefloch et al. 2002), the Rosette Nebula (Schneider et al. 2010), M16 (Andersen et al. 2004), M17 (Jiang et al. 2002), 30 Dor (Walborn et al. 2002), and the Small Magellanic Cloud (Gouliermis et al. 2007). In addition, several recent observations of bright rimmed clouds (BRCs; Urquhart et al. 2009; Chauhan et al. 2009; Morgan et al. 2009, 2010) have been carried out. An interesting aspect is the surprisingly spherical shape of many observed nebulae, especially in RCW 120, "the perfect bubble" (Deharveng et al. 2009; Zavagno et al. 2010). Other regions, such as, e.g., RCW 79 (Zavagno et al. 2006), RCW 82 (Pomarès et al. 2009), RCW 108 (Comerón & Schneider 2007), and Sh 104 (Deharveng et al. 2003; Rodón et al. 2010), share this morphology.

In general, the pillars point like fingers toward the ionizing source and show a common head-to-tail structure. Most of the mass is concentrated in the head which has a bright rim facing the young stars (e.g., Gahm et al. 2006). Thin, elongated pillars connect the head with the main body of the molecular cloud. They have typical widths of 0.1–0.7 pc and are 1–4 pc long (Gahm et al. 2006; Schuller et al. 2006). The observations show that the pillars are not smooth, but show small-scale structure, filaments, and clumps (Pound 1998). Some filaments run diagonal across the pillars, suggesting a complex twist into a helical structure (Carlqvist et al. 2003). This is also supported by spectroscopic measurements of the line-of-sight (LOS) gas velocity: the pillars show a bulk motion away from the ionizing stellar sources with a superimposed complex shear flow that could be interpreted as corkscrew rotation (Gahm et al. 2006). Occasionally, close to the tip of the head small spherical gas clumps are observed to break off and float into the hot H ii region. These so-called evaporating gaseous globules (EGGs) have been found with Hubble Space Telescope (HST), e.g., in the Eagle Nebula (McCaughrean & Andersen 2002). If stars with surrounding gas disks happen to form in these clumps, they transform into evaporating proto-planetary disks, so-called proplyds. More recent observations (Gahm et al. 2007) have revealed a wealth of even smaller sized globules or globulettes, which are in general not bright rimmed and show no detectable sign of star formation. Direct signatures of star formation are found in the head of the pillars, e.g., through jets from obscured proto-stars piercing through the surface into the H ii region (e.g., in Eta Carina; Smith et al. 2000). Whether these jets are preferentially aligned perpendicular to the trunk is a matter of current debate (Raga et al. 2010; Smith et al. 2010).

On the theoretical side, early models of pillar formation suggested that they form by Rayleigh–Taylor instabilities when the expanding hot, low-density H ii region radially accelerates the cold, dense gas (Frieman 1954). This has been ruled out by the observations of the complex flows inside the pillars (Pound 1998).

Another possibility is the collect and collapse model. Here, the radiation sweeps up a large shell, which then fragments to form stars and pillars (Elmegreen & Lada 1977; Klein et al. 1980; Sandford et al. 1982). However, the timescales (>5 Myr) and masses (>1000 M) involved (Elmegreen et al. 1995; Wünsch & Palous 2001) are much larger than in M16. Therefore, this is a more likely scenario for supernova-driven shells.

A third scenario is the radiation-driven implosion (RDI, e.g., Bertoldi 1989) of pre-existing dense cores. This has been studied in great detail with numerical simulations (e.g., Lefloch & Lazareff 1994; Williams et al. 2001). More recently, Kessel-Deynet & Burkert (2003) presented three-dimensional RDI simulations with a smoothed particle hydrodynamics (SPH) code and were able to show that an otherwise gravitationally marginally stable sphere can be driven into collapse by ionizing radiation. In Gritschneder et al. (2009a, hereafter G09a), we showed that marginally stable density enhancements get triggered into forming stars in cases with high as well as low ionizing flux. Miao et al. (2009) further analyzed this RDI scenario with an SPH-based radiative transfer scheme. They show that there is an evolutionary sequence, depending on the initial size of the cloud, as suggested by Lefloch & Lazareff (1994). Bisbas et al. (2009) studied the implosion of a single clump with a new ray-tracing scheme, based on the HEALPix algorithm. A new implementation in the adaptive-mesh-refinement (AMR) code FLASH using the hybrid characteristics ray tracing was achieved by Peters et al. (2009, 2010). Very recently, Mackey & Lim (2010) were able to reproduce the density structure of the main pillar in M16 by setting up several pre-existing dense cores in a triangular way. They investigate the effects of different cooling recipes in a grid code but they do not include self-gravity.

Recently, the focus has moved toward the ionization of the turbulent interstellar medium (ISM). Mellema et al. (2006) reproduced the observed morphologies of H ii regions by ionizing a turbulent medium using a grid code without the inclusion of gravity. Dale et al. (2007) used an SPH code to compare the gravitational collapse of a molecular cloud with and without ionization. They found slightly enhanced star formation in the simulation with ionization. The inclusion of ionization in a grid code in combination with a magnetic field was discussed by Krumholz et al. (2007). A homogenous magnetic field leads to a non-spheric H ii region, as the gas is held back by the magnetic field lines and an oval-shaped bubble develops. However, all these simulations focused on the evolution of the entire H ii region and therefore lack the resolution for a detailed kinematical and structural analysis of the pillars. For a comprehensive quantitative study, we investigated the radiative ionization of a turbulent molecular cloud from a nearby star cluster with so far unprecedented resolution by zooming into a subregion of the cloud. First, we focused on the evolution of the power spectrum in the interaction zone of a turbulent cloud affected by ionization (Gritschneder et al. 2009b, hereafter G09b). Recently, Lora et al. (2009) further investigated the ionization of a turbulent cloud by the combination of a two-temperature equation of state with gravitational forces and transfer of ionizing radiation. They produce pillar-like structures including cores. The angular momentum of these cores is preferentially aligned perpendicular to the direction of the ionizing field.

A different approach was presented by Nayakshin et al. (2009). They focused on the momentum transfer in regimes of high radiation pressure by combining SPH with a Monte Carlo approach. Another main motivation of developing software able to treat ionizing radiation has been investigations of the re-ionization of the early universe (see, e.g., Iliev et al. 2006, 2009; Pawlik & Schaye 2008; Altay et al. 2008, and references therein).

In this work, we investigate for the first time the effect of different initial conditions on the ionization of turbulent molecular clouds. We vary the initial temperature, the mean number density, the level of turbulence as well as the turbulent scale and the ionizing flux. The structure of this paper is as follows. In Section 2, we briefly review the concept of ionizing radiation, followed by a short summary of the iVINE code. After that we present the set of initial conditions for the parameter study. In Section 3, the outcome of the different simulations is discussed in detail. Section 4 is dedicated to the stability and shape of the structures in dependence on the initial conditions. A close comparison to the observed masses, morphologies, and LOS velocities is done in Section 5. We draw the conclusions in Section 6.

2. BASIC APPROACH AND INITIAL CONDITIONS

2.1. Ionizing Radiation

As soon as an O star is born, it ionizes its surroundings with its UV radiation. This leads to an ionized, hot H ii region (Tion ≈ 104 K). In the beginning, the "R-type" ionization front travels with a speed vR which is larger than the sound speed of the hot gas aion. This phase ends as soon as ionization is balanced by recombination in the H ii region. The ionized volume VS, the so-called Strömgren sphere (Strömgren 1939), is then given by

Equation (1)

Here, JLy is the flux of ionizing photons of the source, which is assumed to be constant and monochromatic, αB is the recombination coefficient, ρ0 is the density of the pre-existing, homogeneous gas, and mp the proton mass. At a larger distance from the star, the radiation can be assumed to impinge onto a surface in a plane-parallel way and thus Equation (1) simplifies to

Equation (2)

where xS is the thickness of the ionized region and FLy is the infalling ionizing flux per unit time and unit area.

After a hydrodynamical timescale, the hot gas reacts to its increased temperature and pressure. The pressure of an ideal one-atomic gas is given by

Equation (3)

where ρ is the density, kB is the Boltzmann constant, μ is the mean molecular weight, and cs is the isothermal sound speed. Now the evolution is characterized by an isothermal shock followed by a weaker, "D-type" ionization front. The front velocity is now vD < aion. For a full analysis see, e.g., Shu (1991). As the hot gas expands, its density is reduced. At the same time, the cold surrounding gas is compressed. Under the assumption that the homogeneously ionized region consumes all UV photons of the source it follows from Equation (2) that the density of the hot gas for a constant flux and a constant temperature Thot at any given time is

Equation (4)

To calculate the front position x(t), we follow the approach of Dopita & Sutherland (2003). Under the assumption of a thin shock which is traveling at a speed vs, the ram pressure in the hot, ionized gas has to be equal to the ram pressure in the cold gas:

Equation (5)

where

Equation (6)

The pressure on the ionized side of the shock is mainly given by the thermal pressure of the hot gas:

Equation (7)

where we have already introduced a constant fitting factor f to account for the approximations made (e.g., the one leading to Equation (4)). Combining Equations (6) and (7) and using Equations (2) and (4) yield

Equation (8)

With the initial condition x(t0) = xs, we can solve by integration and obtain

Equation (9)

Using Equation (4), it follows that

Equation (10)

for a plane-parallel infall of a constant flux onto a homogeneous medium.

2.2. Numerical Method and First Tests

In order to investigate the effect of different initial conditions and levels of UV radiation on the formation of pillars, we conduct a parameter study. All simulations were performed with iVINE (G09a), an implementation of ionizing radiation in the tree-SPH code VINE (Wetzstein et al. 2009; Nelson et al. 2009). As the mean free path of the atoms and electrons in the cold and the hot phase is of order λ ≈ 1012 cm and λ ≈ 1014 cm, respectively (see, e.g., Shu 1991), which is much smaller than the length scales involved in our simulations, the gas can be treated by the means of fluid dynamics. The evolution of a turbulent molecular cloud under the influence of ionization spans several orders of magnitude in density. We therefore chose to solve the hydrodynamic equations with the method SPH, a Lagrangian approach with adaptive resolution. To prevent artificial fragmentation (Bate & Burkert 1997), the Jeans length has always to be resolved with at least 50 particles. This leads to a resolution limit which we ensure to be small enough by using a sufficient amount of total particles (see Section 3.7).

In iVINE, the ionizing radiation is assumed to impinge plane parallel onto the simulated volume from the negative x-direction. From the surface of infall, the radiation is propagated along the x-direction by a ray-shooting algorithm. Along these rays, the ionization degree ηi is calculated for each particle i. According to the ionization degree, the pressure Pi of the particle is calculated by a linear interpolation between the temperature Thot of the hot, ionized and the temperature Tcold of the cold, unionized gas. Here, we assume both gas components to be isothermal, since for the density range considered in our simulations heating and cooling should balance each other to approximate isothermality (see, e.g., Scalo et al. 1998). Following Equation (3), the new pressure in our simulation is given as

Equation (11)

where ρi is the SPH density of the particle i and μion = 0.5 and μnion = 1.0 are the mean molecular weights of the ionized and the unionized gas in the case of pure hydrogen, respectively.

As a first test, we verify Equation (9) and fit a value for f by ionizing a slab of atomic hydrogen with a constant homogeneous density of ρcold = 300 mp cm−3 and temperature of Tcold = 10 K. We perform three different runs, corresponding to a low flux (FLy = 1.66 × 109 γ cm−2 s−1), an intermediate flux (FLy = 5 × 109 γ cm−2 s−1), and a high flux (FLy = 1.5 × 1010 γ cm−2 s−1). This corresponds to the ionization penetrating immediately into the first 0.55%, 1.67%, and 5% of the region, respectively (see Equation (2)). At the same time, this is equal to placing the simulation volume further away or closer to the source, e.g., the O star. The simulations are conducted with the same accuracy and setup as in the parameter study given below (see Section 2.3). Figure 1 shows the resulting evolution of the front. As one can clearly see, the approximations leading to Equation (2) (f = 1, the dotted lines in Figure 1) do not produce satisfactory results. Instead, assuming

Equation (12)

(i.e., $f=\sqrt{\frac{5}{4}}$, the dashed lines in Figure 1) perfectly matches the simulations during the entire simulated time of tsim = 500 kyr. Thus, we keep this assumption for this work.

Figure 1.

Figure 1. Front position vs. time for the three test simulations with a different flux impinging on a homogeneous medium. Green, blue, and red line: simulations with a low, intermediate, and high flux, respectively. Black lines: solution according to Equation (9), dotted f = 1, dashed $f=\sqrt{\frac{5}{4}}$.

Standard image High-resolution image

2.3. Initial Conditions

To produce different turbulent initial conditions, we use the same approach as in G09b. We set up 2 × 106 particles to resemble a homogeneous medium with ρcold = 300 mp cm−3 in a volume of (4 pc)3. Then, we add a supersonic velocity field (Mach 10 in most cases) with a steep power law E(k) ∝ k−2 on the modes k = 1...4. Before switching on the ionizing source, each setup decays freely under the influence of isothermal hydrodynamics and periodic boundary conditions until the desired initial Mach number is reached (after t ≈ 0.8–1.0 Myr, depending on the specific simulation). This self-consistent evolution of the turbulence leads to a combination of solenoidal and compressible modes, which seems to be the case in turbulent molecular clouds.7 The individual particle time steps in iVINE are determined as in G09b by using an accuracy parameter of τacc = 1.0 and a Courant–Friedrichs–Lewy (CFL) tolerance parameter of τCFL = 0.3. An additional time step criterion based on the maximum allowed change of the smoothing length with an accuracy parameter of τh = 0.15 is also employed.

After the desired Mach number is reached, the ionization is turned on. Now, the boundaries are still periodic in the y- and z-directions. In the negative x-direction, the boundary is reflecting to represent conservation of flux toward the star, whereas in the positive x-direction the gas is allowed to stream away freely. To test this approach, we perform one simulation with open boundaries in both x-directions. Gravitational forces are calculated without boundaries. This is valid as the free-fall time of the whole simulated area is tff ≈ 3 Myr, which is much longer than the simulation time of tfinal = 0.5 Myr. For the tree-based calculation of gravitational forces, we use a multi-pole acceptance criterion (MAC; Springel et al. 2001) with a tree accuracy parameter of θ = 5 × 10−4. The correct treatment of the ionization and the resulting acceleration of the particles are guaranteed by the modified CFL-condition discussed in G09a. The recombination of the hot gas is modeled assuming αB = 2.59 × 10−13 and the cross-section for the ionizing photons is set to σ = 3.52 × 10−18 cm2. The initial properties of the different simulations are listed in Table 1. The simulations were performed on an SGI Altix 3700 Bx2 supercomputer, the calculation of each setup took approximately 100 wall clock hours on 16 CPUs.

Table 1. Listing of the Different Initial Conditions

Simulation M $\bar{\rho }$ lbox FLy Mach k Tnion log10(Npart)
  (M) ( mp cm−3) ( pc) (γ cm−2 s−1)     (K)  
Low resolution 474 300 4 5 × 109 5 1–4 10 5.4
Open boundaries 474 300 4 5 × 109 5 1–4 10 6.3
Mach 1.5 474 300 4 5 × 109 1.51 1–4 10 6.3
Mach 4 474 300 4 5 × 109 4 1–4 10 6.3
Mach 5 (G09b) 474 300 4 5 × 109 5 1–4 10 6.3
Mach 7 474 300 4 5 × 109 7 1–4 10 6.3
Mach 12.5a 474 300 4 5 × 109 12.5 1–4 10 6.3
Mach 5 warm 474 300 4 5 × 109 5 1–4 100 6.3
M5 warm high flux 474 300 4    1.5 × 1010 5 1–4 100 6.3
Low flux 474 300 4 1.7 × 109 5 1–4 10 6.3
High flux 474 300 4    1.5 × 1010 5 1–4 10 6.3
Low density 158 100 4 1.7 × 109 5 1–4 10 6.3
Smaller kmax 474 300 4 5 × 109 5 4–8 10 6.3
Small boxb 119 300 2 5 × 109 5 1–4 10 6.6
Big box 3795 300 8 5 × 109 5 1–4 10 6.3

Notes. Given are initial mass, average density, and size of the simulation. In addition, the impinging flux, turbulent Mach number, the largest driving mode of the turbulence, and the temperature are listed. Mach 5 (G09b) is the standard case as presented in G09b. aIn the case Mach 12.5 we start with an initial velocity field of Mach 20 which is then allowed to decay freely as prescribed in Section 2.3. bTo allow for comparable results the box is 2pc in the y- and z-direction but 4pc in the x-direction. Therefore, the particle number is increased.

Download table as:  ASCIITypeset image

3. RESULTS OF THE PARAMETER STUDY

3.1. General Properties

The time evolution of the density for turbulent regions with Mach numbers of 1.5, 5, 7, and 12.5, respectively, is shown in Figure 2. In all of the simulations forming structures, the same effect as described in G09b takes place. The ionizing radiation penetrates deeper into the turbulent cloud along low-density channels. As the ionized gas reacts to its increase in pressure, it starts to compress the adjacent, unionized, higher density regions, thereby widening the channels of low density and thus allowing the ionization to penetrate even further. We call this process "radiative round-up." At t = 250 kyr, the pre-existing high-density structures have been enhanced by the outside compression and pillars start to become visible. After t = 500 kyr, the pillars have achieved characteristic shapes which match the observations remarkably well. Figure 2 (row 3) and Figure 3 show the density projected along the z-axis at this stage for all simulations of the parameter study.

Figure 2.

Figure 2. Time evolution of four different initial conditions. We show the density projected along the z-axis of four different simulations at subsequent stages, from left to right column by increasing Mach number, as indicated. Color-coded is the surface density, each figure is 4 pc × 4 pc. The significant structures only form above a certain level of turbulence (M ⩾ 2) and get less stable with increasing Mach number.

Standard image High-resolution image
Figure 3.

Figure 3. Final stage of the different simulations. Color-coded is the surface density of the simulations given in Table 2 as indicated in the top left corner, each figure is 4 pc × 4 pc. Note that in panels 6, 8, and 9 an earlier stage is depicted since the evolution of the simulation is more rapid. Since all panels have the same physical size panel 12 shows only the top left quarter of that simulation which includes the most significant structure.

Standard image High-resolution image

For a quantitative discussion, we investigate the most prominent pillar structure in each simulation in more detail. To define the tip, we take the particle closest to the source of radiation8 above a threshold density of ρtresh = 104mp cm−3. We then take its surrounding, the region spanning 1 pc in the x-direction and 0.3 pc in the negative and positive y- and z-directions. The cold, unionized gas in this region is defined as the pillar.9 This definition allows us to extract the important quantities of the most prominent structure by the same algorithm for all simulations. The characteristic values, which allow for a comparison with the observations as well as a deduction of the underlying physics, are given in Table 2 for the defined pillar at tfinal = 500 kyr. Note that due to the adaptive nature of SPH we have a very high resolution in the prominent structures. The pillar in Mach 5 (G09b), e.g., contains 5.6 × 104 particles, which corresponds roughly to a spatial resolution of 178 × 17.8 × 17.8 on the 1 pc × 0.1 pc × 0.1 pc of this pillar. This resolution is up to now unprecedented and allows for a detailed comparison of the kinematics of this pillars with the observations (see Section 5).

Table 2. Results of the Parameter Study at t = 500 kyr

Simulation M $\bar{\rho }_{\rm p}$ $\bar{\Sigma }$ log10[N(H2) σ $\bar{v}_{\rm x}$ d $\bar{\rho }_{\rm ion}$ ΔP
  (M) (104mp cm−3) (10−3gcm−2) (cm−2) (km s−1) (km s−1) ( pc) ( mp cm−3)  
Low resolution 14.26 5.66 1.69 20.5 0.9 ± 0.5 4.8 ± 0.5 0.11 42.5 1.50
Open boundaries 11.0 5.06 1.44 20.5 1.2 ± 0.7 5.4 ± 1.1 0.10 28.6 1.13
Mach 1.5 6.2 3.46 1.18 20.4 1.6 ± 1.7 4.5 ± 1.5 0.10 41.8 2.42
Mach 4 12.6 3.62 1.32 20.4 1.3 ± 0.6 3.9 ± 0.9 0.13 39.8 2.20
Mach 5 (G09b) 250 kyr 17.0 3.60 1.43 20.5 1.8 ± 1.7 4.0 ± 1.7 0.16 60.6 3.46
Mach 5 (G09b) 12.6 4.56 1.52 20.5 1.1 ± 0.7 4.8 ± 0.9 0.12 43.4 1.90
Mach 7 14.0 5.36 1.71 20.6 1.1 ± 0.5 4.0 ± 0.9 0.11 39.2 1.46
Mach 12.5 8.1 5.66 2.56 20.7 1.3 ± 1.9 4.0 ± 1.7 0.09 43.5 1.54
Mach 5 warm 10.2 0.37 0.23 19.7 2.0 ± 0.8 3.3 ± 1.7 ... 43.6 2.36
M5 warm high flux 250 kyr 7.1 2.59 1.0 20.3 3.5 ± 3.6 7.0 ± 3.7 0.11 116 0.90
M5 warm high flux 20.5 2.38 0.84 20.2 2.3 ± 1.7 8.0 ± 2.1 0.22 92.0 0.77
Low flux 12.4 3.75 1.38 20.5 0.9 ± 0.5 3.3 ± 0.7 0.13 26.5 1.41
High flux 250 kyr 10.9 6.73 2.80 20.8 2.4 ± 3.3 6.3 ± 3.3 0.09 108 3.22
High flux 10.9 7.23 2.72 20.8 2.3 ± 1.3 8.2 ± 2.2 0.09 73.5 2.02
Low density 250 kyr 3.1 4.46 1.82 20.6 3.2 ± 3.8 7.0 ± 4.0 0.06 38.7 1.74
Low density 5.79 2.80 0.99 20.3 2.0 ± 0.9 6.1 ± 1.9 0.10 19.0 1.36
Smaller kmax 2.78 3.44 1.16 20.4 1.2 ± 0.7 6.2 ± 0.9 0.06 42.2 2.45
Small box 3.39 2.76 1.00 20.3 1.6 ± 0.8 3.4 ± 1.5 0.08 38.2 2.77
Big box 51.1 4.76 1.94 20.6 1.0 ± 0.5 2.9 ± 0.7 0.24 42.5 1.78

Notes. Listed are the mass, mean density, mean surface density, corresponding column depth, velocity dispersion, and the x-velocity away from the source of the most prominent structure. Then, the mean diameter (see Equation (13)) of the pillar and the mean density of the hot gas are given. Finally, the pressure difference ΔP (see Equation (14)) is listed. For the fiducial simulation as well as the more rapidly evolving simulations, these quantities are given at an earlier stage (t = 250 kyr) as well.

Download table as:  ASCIITypeset image

We calculate the diameter of a pillar via

Equation (13)

with x = 1 pc as the length of the pillar. In addition, the mean density10 of the hot gas is listed in Table 2. It is notable that all simulations with the same impinging flux share a similar density of the pillar as well as of the hot gas.

As both gas phases are treated isothermal (cf. Equation (3)), the pressure difference at tfinal = 500 kyr is

Equation (14)

ΔPfinal is very close to unity for all simulations.11 Therefore, we derive as a first result that the pillars are in thermal equilibrium with the hot surrounding gas.

3.2. Resolution and Boundary Conditions

The first two simulations do not address different physical properties but rather numerical details. Simulation low resolution was performed with exactly the same setup as Mach 5 (G09b), but with eight times less particles. This leads to a two times lower spatial resolution. Nevertheless, the morphology is comparable to the high-resolution case (Figure 3, panel 2). The only noticeable difference is that the second largest structure in this case has already merged with the third structure. Furthermore, tiny structures are less frequent. The physical properties (see Table 2) are similar as well. The structures in the low-resolution case tend to be a bit more massive, which can be expected. Altogether, the morphology and the global physical properties are comparable and thus, we conclude that Mach 5 (G09b) is reasonably converged.

In the other test case, we investigate the boundary condition in the negative x-direction. In open boundaries, this boundary is not reflecting. Instead, the gas is allowed to stream away freely. This leads subsequently to a lower density in the ionized region. As a consequence, open boundaries (Figure 3, panel 3) look similar to high flux (panel 8, see Section 3.5), as the radiation is able to penetrate further into the computational domain. Nevertheless, the formation of pillars still takes place and is not strongly affected. Even density and mass assembly of the most prominent structure are alike. Therefore, the choice of the boundary condition does not influence the overall scenario significantly. As it is more realistic to assume hot gas is already present in the region between the ionizing source and the simulated part of the molecular cloud, we keep the reflecting boundary condition in all other simulations. These reflections can be interpreted as flux conservation at the left border of the simulation: as much gas streams from the area toward the source into the region as is streaming outward.

3.3. Turbulent Mach Number

A main purpose of this study is to disentangle the effects of the initial turbulent density distribution on the formation of the pillars. Therefore, the level of turbulence is changed. We take different turbulent setups: one has been evolved from a very high Mach number (Mach 20) to Mach 12.5. The other four represent different stages of the decay starting from our fiducial turbulent setup (G09b) at Mach 10. They are taken at Mach 7, Mach 5 (G09b), Mach 4, and Mach 1.5, respectively. When non-driven turbulence decays, most power is lost on the large-scale modes. This can be seen in Figure 2. In Mach 5 (G09b) (Column 2) and Mach 7 (Column 3), the surface density is clearly dominated by the large modes, which form the prominent fingers. In contrast, in Mach 1.5 (Column 1) no significant pillars evolve, since the initial density distribution is already too smooth and the dominant mode has decayed too far. This trend can already be seen in Mach 4 (Figure 3, panel 4), where the structures are less distinct. Mach 12.5 (Column 4) is a much more violent case. Since there is a lot of power on the largest density scale, structures are evolving. However, these are already being torn apart at the same time, as discussed in Section 4.

Overall, the evolution is mainly dominated by the pressure differences between hot and cold gas. Compared to the increase in the pressure due to the ionization (3 orders of magnitude), the differences from varying the Mach number are small. However, a small trend is visible in the average density of the assembled structure (see Table 2). The higher the Mach number, the higher the density of the formed structure. That can be directly related to the density of the initial turbulent filament, which is denser at a higher Mach number. This effect is also visible from the first row of Figure 2, where the absolute value in the area of highest density is increasing from left to right.

3.4. Temperature and Pressure

The most striking difference can be seen between Mach 5 (G09b) (Figure 3, panel 1) and Mach 5 warm (panel 5). Both initial conditions are self-similar, both were set up with the same initial random seed for the turbulence, and they are relaxed until their velocities resemble Mach 5 at their respective temperature. Consequently, at the time the ionization is switched on, their density distribution is identical. Since the impinging flux is the same, the radiation ionizes the same regions in both cases. Nevertheless, in Mach 5 warm no evolution of any filamentary structure is visible. This leads to the conclusion that the pressure balance between the hot, ionized and the cold, unionized gas plays a crucial role in the formation of structures.

Taking Equation (3) and the straightforward assumption that only regions with a pressure lower than the pressure of the hot gas can be compressed gives

Equation (15)

and thus

Equation (16)

If we assume ρion = 100 mp cm−3 in the beginning, as the ionization mainly penetrates the lower density regions, this equation yields ρnion,10 K ⩽ 3.6 × 105mp cm−3 and ρnion,100 K ⩽ 3.6 × 104mp cm−3. The maximum density ρmax = 8.8 × 104mp cm−3 in both simulations lies in between these thresholds. Thus, in Mach 5 (G09b), the pressure of the ionized gas is high enough to compress even the densest structures, whereas in Mach 5 warm several regions are able to resist the compression. Therefore, the ionized "valleys" are not expanding significantly in the tangential direction, the density is not lowered as much and the ionization is not able to reach much further. At the same time, the unionized "hills" are less compressed, but since they are closer to the front they are accelerated more strongly in the x-direction and the initial differences in the front position are leveled out.

In general, the formation of pillars critically depends on whether the density contrast between the dense regions (ρhigh), which cannot be ionized, and the less dense regions (ρlow), which can be ionized, is lower than the temperature ratio between ionized and unionized gas. By defining the density contrast in the initial conditions as Δρinit = ρhighlow and taking into account Equation (15), the critical criterion for the formation of pillars can be written as

Equation (17)

Since stars will only form in compressed regions, e.g., pillars, this gives an estimate if a region will undergo triggered star formation.

To test this condition, we used the same warm initial conditions but increased the ionizing flux in M5 warm high flux. The flux is now at the level of high flux (see Section 3.5). This is equivalent to increasing the lowest density which can still be ionized, i.e., to increasing ρlow. As a result, Δρinit is decreased and thus structure formation is possible again in the warm case according to Equation (17). In fact, it can be seen in Figure 3 (panel 6) that pillar formation is indeed triggered again.

Unfortunately, there is no straightforward way to define the density contrast in a turbulent medium, especially since the impinging flux plays a major role in defining "high" and "low" density as seen in M5 warm high flux. However, Equation (17) already shows that increasing the mean density $\bar{\rho }$ while keeping the temperature constant will not help to hinder the formation of pillars. This is supported by the fact that pillars form in Mach 5 (G09b) as well as in low density (see Section 3.5).

3.5. Initial Flux and Density

In the next test, we vary the impinging photon flux. As in the simulations performed in Section 2.2, the flux is able to ionize immediately 0.55%, 1.67%, and 5% of a medium at a constant density of ρ = 300 mp cm−3 in low flux, Mach 5 (G09b) (intermediate flux), and high flux, respectively. The evolution of the density in the hot component is shown in Figure 4. Although the medium is highly turbulent Equation (10) still gives a very good estimate of this density, especially after an initial phase. Only the case with the high flux differs from the analytical solution. This can be understood as Equation (4) depends on both the penetration depth and the density of the ionized gas. A higher flux can ionize a larger fraction of the computational volume straight away. Equivalently, the evolution of high flux would follow Equation (4) more closely if the mean gas density would be higher, thus resulting in a shorter penetration length. As long as the penetration length is relatively small (<5%), the turbulent case is still comparable with the case of a constant flux and Equation (10) still represents a valid description of the evolution of a turbulent H ii region.

Figure 4.

Figure 4. Change of the mean density in the hot gas during the simulated time. Solid lines: green low flux, blue Mach 5 (G09b), red high flux, and yellow low density. Dotted lines: the comparison simulations of Section 2.2 at the respective flux. Dashed black lines: the analytical solution according to Equation (4) with $f=\sqrt{\frac{5}{4}}$.

Standard image High-resolution image

From a morphological point of view, the pillars in the simulations with a higher flux (i.e., when the computational volume is located closer to the ionizing source) are smaller than in the case of a lower flux (Figure 3, panel 8, panel 1, and panel 7 in decreasing flux order). In addition, they gain more momentum away from the source and move faster away from the source as the photoevaporation is stronger. At the same time, the density of the hot gas is higher, leading to denser, more compressed structures with a smaller diameter. Due to a higher photoevaporation rate, their average masses are as well lower (see Table 2).

Changing the initial flux is expected to have a similar effect as changing the mean density, as ρ ∝ x2 (cf. Equation (2)). In low density, we reduced the mean density by a factor of 3. At the same time, we reduced the flux by a factor of 3 to avoid an extremely high level of ionization degree. In total, this corresponds to the same penetration length as in high flux. Thus, we expect a similar morphology to evolve, but the densities should be lower. In Figure 3 (panel 9), this can be clearly seen. The morphology is similar to high flux, the front is at a similar position. Again, the density (Figure 4) in the hot gas evolves similarly to the expectation for a homogeneous medium. The mass assembled in the most prominent structure (Table 2) is lower and the density of the structure fits the findings of pressure equilibrium (see Section 3.4).

Combining these findings with the results of Section 3.1 allows us to make an interesting prediction. As the density of the hot gas behaves similarly to the case of a homogeneous medium and as the structures are in approximate pressure equilibrium with the surrounding hot gas, we can predict the density of the structures from the initial mean density of the medium, the flux of the source, and the time since the ignition of the source or the position of the ionization front. The density of the forming structures is thus given as

Equation (18)

where we used Equation (10). xs depends on the initial density and the impinging flux (see Equation (2)). As we expect the assumption of pressure equilibrium to hold in the case of a point source, the three-dimensional expression taking into account geometrical dilution is given by

Equation (19)

where RS is the Strömgren radius (for a detailed derivation of the three-dimensional front position see, e.g., Shu 1991).

3.6. Turbulent Scale

To study the effect of the largest scales of the initial turbulence (the turbulent input scales), we compare Mach 5 (G09b) with smaller kmax, a run in which we populate modes k = 4...8, instead of k = 1...4 as usual. The resulting surface density in the first 2 pc facing the star is shown in Figure 5. Already in the initial conditions (left column), a clear difference can be seen. Whereas the power on the larger k-modes leads to large, distinct structures (top panel), power on the smaller modes already shows a much more diversified density distribution (lower panel). After tfinal = 500 kyr (right column), the ionization leads to an enhancement of the pre-existing structure. The densest filaments survive, while the other material is swept away by the ionization. In Mach 5 (G09b) (top panel), this leads to an excavation of the few, but bigger, structures and thus to the creation of few, but distinct, pillars. On the other hand, in smaller kmax (bottom panel) more structures, but of smaller scales, survive, which leads to several, but more diffuse, structures.

Figure 5.

Figure 5. Projected surface density along the x-axis. The projected slice is always 2 pc thick. Left: t = 0 kyr and the slice starts at the surface facing the O star, right: 500 kyr and the slice is adjusted to encompass the sub-structures. Top row: only modes k = 1–4 are populated initially; bottom row: only modes k = 4–8 are populated initially.

Standard image High-resolution image

Together with Section 3.3 this shows that only a strong enough turbulent driving on a large enough driving scale leads to the formation of coherent structures as seen in observations. As has been shown in Section 3.2, this is not an effect of the resolution. The turbulence is well enough resolved to allow for small enough modes to produce fuzzy structure in Mach 5 (G09b), but the evolution under the influence of UV radiation is dominated by the larger modes.

Another possibility to change the input scale of the turbulence is to simply increase or decrease the size of the simulation domain. In big box, the box size is doubled to 8 pc. Since the particle number is kept constant, this leads to a factor of 2 lower spatial resolution. So, the resolution in the part of the domain shown in Figure 3 (panel 12) is comparable to the low-resolution case low resolution (panel 2). In small box, the box size is halved to 2 pc. This corresponds to doubling the spatial resolution. The domain has a smaller extent in the x-direction as well, which we compensate by taking two times the evolved turbulent box in the x-direction. This is valid, since the initial conditions were evolved with periodic boundary conditions. The particle number is thus 4 × 106, twice as high as in most other cases. Small box is the situation in between Mach 5 (G09b) and smaller kmax, as the largest mode is 2 pc, which corresponds to a k = 2 mode in a 4 pc domain.

As the density distribution and therefore the density contrast are self-similar in all three cases, we expect that the same regions of the initial conditions will form the dominant structures. Only the size of the encompassed region and therefore the size and mass of the pillars should change, if the process is indeed scale free. In Figure 3, all three simulations (panel 1, panel 11, and panel 12) show a clear sign of the largest k = 1 mode. The size of the structures formed is linearly dependent on the initial box length or size of the largest k-mode. The values of Table 2 confirm the importance of the largest mode. In the assembled structures, the estimated diameters are roughly a factor of 2 different and the masses vary by a factor of 4. This is as expected since the regions initially encompassed by the radiation should differ by a factor of 2 in the y- and the z-direction.

Taking these results on the turbulent scale into account, we conclude that the mass and size of the pillars are directly dependent on the input scale of the turbulence, e.g., the size of the driving process or the size of the pre-existing molecular cloud. On average, the most prominent structures in our simulations with an intermediate flux are $d_{\rm pil}\approx \frac{1}{40}x_{\rm turb}$, where xturb is the largest turbulent input mode.12

3.7. Star Formation

In several simulations, triggered dense regions form cores and are driven into gravitational collapse. Since star formation is not the main goal of this study we do not replace them by sink particles. Instead, we remove the particles forming a core from the simulation to avoid a considerable slowdown of the calculation. Following G09a, we define a core as all gas with a density above ρcrit = 107mp cm−3 in the region around the density peak and remove the particles representing this core. The core formation is not a numerical effect, since the resolution limit as given by Bate & Burkert (1997) is ρnum = 3 × 108mp cm−3 in the lowest resolution case. We give the simulation, mass,13 formation time, average speed away from the source, and positions in Table 3. If we assume the cores to be decoupled from the rest of the cold gas then their position at the end of the simulation can be estimated by these positions and velocities. All of them are still close to the prominent structures, some are traveling further inside the structures, and some are lagging behind and would by now be slightly outside of the pillar, closer to the source.

Table 3. Listing of the Proto-stellar Cores Forming in the Different Simulations

Simulation tform M vx x y z
  t(kyr) (M) (km s−1) ( pc) ( pc) ( pc)
Big box 305.3 0.86 4.37 ± 1.4 −2.83 −1.84 1.26
Big box 353.3 0.90 4.16 ± 0.7 −2.60 −1.86 1.34
Big box 403.7 0.78 4.53 ± 1.1 −2.58 0.29 0.41
Big box 469.3 0.83 4.99 ± 0.7 −2.06 −1.88 1.48
High flux 429.5 0.57 11.27 ± 1.0 1.86 −1.05 0.64
Mach 5 (G09b) 493.3 0.72 3.92 ± 0.8 0.19 0.48 −0.79

Note. Given are mass, formation time, average speed away from the source, and their formation position.

Download table as:  ASCIITypeset image

As can be expected, in big box, where the compressed structures are most massive, star formation is most frequent and happens earliest. There is an age spread present as well—the earlier a core forms the closer to the source it will be. The other simulations where cores form during the first 500 kyr after the ignition of the O star are high flux and Mach 5 (G09b). The higher flux leads to a higher compression and thus an earlier formation of a core in high flux. Due to a higher photoevaporation rate, this core is also less massive and moves faster compared to the core in Mach 5 (G09b).

Altogether, triggered star formation is very likely in this scenario. The cores form at the center of the structures, but since their velocities differ from the velocity of their parental structure they can be decoupled and either wander further into the trunk or lag behind. This depends on the specific environment and does show no correlation with, e.g., the initial flux. At the time they become observable, there might be no clear correlation to their birthplace anymore.

4. PILLARS, CAPS, AND GLOBULES

Another interesting feature can be found when looking at different projections of the simulation Mach 7 (Figure 6). By comparing both surface densities, it becomes clear that the density enhancements inside the pillars are caps, which are exposed directly to the UV radiation in the other projection. Thus, the small steps and wiggles seen in the observations can be explained directly by the sweeping up of smaller caps, which are not shadowed by the leading tips of the pillars.

Figure 6.

Figure 6. Projected surface density of the pillars in the simulation with Mach 7. (a) x–y projection, (b) x–z projection. The surface density is color-coded, each figure is 1 pc × 1 pc. The density enhancements inside the pillars at the left-hand side directly correspond to caps, which are not shadowed by the leading tip on the right-hand side and vice versa.

Standard image High-resolution image

As a last point, it is interesting to examine the stability of the forming structures. In Table 4, we show the correlation between the formation of globules or pillars and the initial turbulent velocity perpendicular to the plane of irradiation of the particles forming these structures. Here, we determine the mass of a globule by integrating over its material with a density higher than nc = 104 cm−3. A tip is defined similarly by taking all the material above nc in a sphere of Rc = 0.025 pc around the local density maximum. Table 4 lists the most important features seen in the final stage of Figure 2. In the physical properties of the tips and globules, a striking correlation can be seen. The tips have lower initial and final tangential velocities, the globules much higher ones. Thus, a pillar can only be assembled when the leading blob is moving slow enough perpendicular to the irradiation that a sub-structure can survive in its shadow. Therefore, an environment with lower Mach number favors the formation of stable, massive structures, as in the simulation with Mach 5 (G09b). Mach 7 represents the intermediate case where tips and globules form. The case of Mach 12.5 corresponds to a rather violent scenario. Here, globule 1 has already decoupled from the rest of the gas. In addition, even the most prominent pillar (Tip 1) has a high velocity and represents a transient stage, which soon decouples from its stem and gets disrupted during the next 100 kyr. In contrast, Tip 2 is an exceptional case—due to its very low tangential velocity it allows for a stable pillar even in this violent environment.14

Table 4. Properties of the Globules and Tips in the Four Simulations

Simulation Mach 5 Mach 7
  Tip 1 Tip 2 Tip 3 Tip 1 Tip 2 Tip 3
M (M) 0.62 1.87 0.35 0.79 0.45 0.48
vT0(km s−1) 0.60 1.08 0.54 0.98 0.56 1.04
vT500(km s−1) 0.10 0.47 0.25 1.24 0.32 0.45
Simulation Mach 7 Mach 12.5
  Glob. 1 Glob. 2 Glob. 3 Glob. 1 Tip 1 Tip 2
M (M) 0.32 0.15 0.20 0.12 1.12 0.54
vT0(km s−1) 2.44 2.31 2.20 3.98 1.36 0.55
vT500(km s−1) 3.00 2.85 2.05 2.19 2.42 0.78

Notes. The structures are chosen from the final stage at t = 500 kyr, the numbers are given from top to bottom. Listed are the mass, the initial tangential velocity (vT0) of the particles that are going to form the structure later, as well as their final tangential velocity (vT500).

Download table as:  ASCIITypeset image

Furthermore, it is worth noting that the tangential velocity is not affected strongly by the UV radiation. This is reasonable, since the surviving structures experience tangential compression from all sides, so the net effect of the ionization balances out. In general, our simulations are able to explain the formation of the observed low mass globules or globulettes. The similarity between the structures forming in Mach 7 and the structures observed in the Rosette Nebula (Gahm et al. 2007, their Figure 3) is especially striking.

5. COMPARISON TO OBSERVATIONS

5.1. General Properties

At first, we compare the density of the hot gas to observations. Lefloch et al. (2002) estimate the electron density of the H ii region in the Trifid Nebula15 from O iii as ne = 50 cm−3. In a fully ionized region, this corresponds directly to the density given in Table 2, since ne = nH = ρ/mP. Thus, the observed value is very similar to the density found in all simulations with an intermediate flux at t = 500 kyr. The average density of the pillars is around 104mp cm−3, depending on the individual simulation. This is in very good agreement with recent results from Herschel by Schneider et al. (2010), where they find a typical average density of 1.1 × 104mp cm−3 in the cold structures in the Rosette Nebula.16 In the Trifid Nebula, Lefloch et al. (2002) estimate N(CS) ≈ 1.8 × 1013 cm−2 for the column density of the dense core, which corresponds to log10[N(H2)/cm−2] ≈ 22.5 with their conversion factor. To compute the H2 column density, we apply a conversion factor of χ = N(H2)/Σ = 0.35, using a hydrogen abundance of X = 0.7 and assuming that all hydrogen is molecular at these densities. In all cases where structures form the column densities are around log10[N(H2) cm−2] ≈ 20.5 (Table 2). As this is the averaged surface density of the entire structure in our simulations, it is 2 orders of magnitude lower than the observed values for the dense cores. In the tips of the pillars (see, e.g., Figure 7), the peak surface density is log10max/(g cm−2)] = −1.2 which leads to a column density log10[N(H2)/cm−2] = 22.12, which is in good agreement with observations. Thompson et al. (2002) estimate log10[N(H2)/cm−2] ≈ 21.3 for the most prominent of the pillars in M16. In RCW 120 (e.g., Deharveng et al. 2009) condensation 4 seems to be a good candidate for triggered star formation. The peak surface density is log10[N(H2)/cm−2] = 22.15–22.49. In addition, Urquhart et al. (2009) investigate a sample of 60 BRCs and find the column densities in cases of triggered star formation (that is in the cases with photodissociation regions (PDRs)) to be log10[N(H2)/cm−2] = 20.9–22.8.

Figure 7.

Figure 7. Most prominent pillar of the simulation with Mach 5 in its reference frame. Color-coded is the surface density, each figure is 0.4 pc × 0.8 pc. (a) x–y projection, (b) xz projection. The superimposed grid denotes the bins along which the LOS velocities in Figure 9 are taken. In order to match the observational beam-size Gahm et al. (2006), each bin is 0.04 pc × 0.04 pc.

Standard image High-resolution image
Figure 8.

Figure 8. For a better comparison, the DQ trunk in NGC 7822 as observed by Gahm et al. (2006) is shown. Depicted is a 0.8 pc × 0.8 pc subset of their Figure 1, roughly to scale with our Figure 7.

Standard image High-resolution image

As all these observations match our simulations, we conclude that the evolution of the density of the hot gas is in good approximation given by the estimate in the case of a homogeneous medium (Equation (10)). Furthermore, since the densities of the compressed structures are reproduced as well, we can assume that the structures are indeed in pressure equilibrium with the hot, surrounding gas. This provides the opportunity to determine the density of the hot gas and the compressed structures directly from the initial mean density, the flux from the source and the time since the ignition of the source or the position of the ionization front, respectively.

5.2. Velocity Field of a Singular Pillar

Another property to compare between simulations and observations is the details of the velocity distribution in the pillars. Since turbulence is a highly complex process it is hard to precisely predict the outcome of simulations. Therefore, we do not set up a simulation to match some specific observation but instead start with realistic initial conditions and then look for an observed counterpart of the outcome of our simulation. In the following, we analyze the velocity structure of a typical pillar in Mach 5 simulation which has a similar morphology and mass as the only trunk with well-observed LOS velocities, the Dancing Queen (DQ) trunk in NGC 7822 (Gahm et al. 2006; see Figure 8). The authors give the diameter as dobs ≈ 0.12–0.15 pc, the total estimated mass from 12CO is Mobs ≈ 9.2 M. This is similar to the simulated pillar (dsim ≈ 0.12 pc, Msim ≈ 12.6 M, see Table 2). If we subdivide the pillar into a head and a body, the mass splits up into Mhead,s ≈ 7.2 M and Mbody,s ≈ 5.3 M (compared to Mhead,o ≈ 5.7 M and Mbody,o ≈ 3.5 M in the DQ trunk).

To enable a more detailed comparison, we impose a grid across the head in the y–x plane (see Figure 7) to resemble the beams along which the LOS velocity was taken. We divide the LOS velocity, ranging from vz = −4 km s−1 to vz = 4 km s−1 into 80 equally sized bins. In each of the velocity bins the mass is integrated. Figure 9 shows the profiles obtained in that way. As we do not take any radiative transfer, temperature dependencies, or chemistry into account in our profiles, they are not as smooth and symmetric as the observed HCO+ profiles. Nevertheless, the similarities are striking.

Figure 9.

Figure 9. LOS velocities along the most prominent pillar in the simulation with Mach 5. The LOS velocities are taken in bins of the size 0.04 pc × 0.04 pc (see Figure 7) to match the observational resolution. The mean velocity in each bin is depicted by the dotted lines. In Figure 10, this velocities are plotted for Columns 2, 3, and 4. Note that in addition to the overall pattern, the line width is in very good agreement with the observations (Gahm et al. 2006).

Standard image High-resolution image

In our simulation, the standard mean deviation for all lines with a peak higher than 5 × 10−3M is on average 0.38 km s−1, which corresponds to an FWHM of 0.94 km s−1. This is very similar to the observed FWHM of about 1.2 km s−1. Thus, the irregular motions in the simulations and observations correspond very well. Our simulations are missing the offset of v0 ≈ −16.4 km s−1, the speed at which the observed DQ is moving with respect to us. In addition, the profile is reminiscent of a rotational pattern, as the peak is shifting from left to right, as well as a gradient along the x-direction, since the peak is shifting from top to bottom. This so-called "cork-screw" pattern (see Figure 10) has often been attributed to magnetic fields. However, our pure hydrodynamical simulations reproduce the pattern, indicating that magnetic fields might play a minor role. Instead, the pattern is produced by confining the motion of the turbulent cold gas inside the pillars by the hot gas surrounding it.

Figure 10.

Figure 10. LOS velocity along the pillar. We show the LOS velocity of the most prominent pillar in the simulation with Mach 5 as function of position along three longitudinal cuts. The x-axis is parallel to the major axis. Stars, diamonds, and triangles correspond to cuts parallel to the x-direction through the center (diamond), the left (star), and right (triangle) sides. For more details, see Figures 7 and 9. The location of the head is at x = 0 pc. Top panel (a): projection along the y-axis. Bottom (b): projection along the z-axis. A velocity gradient along the pillar is clearly visible. The gradient matches the observations of the DQ trunk very well (Gahm et al. 2006). The fact that the velocities for the different cuts rarely cross, especially in the bottom plot, is a signature of rigid rotation. This in combination with the overall gradient produces the so-called corkscrew pattern, as observed.

Standard image High-resolution image

As elongated, pillar-like structures can also be the result of the ionization of pre-existing clumps (RDI model, see, e.g., Mackey & Lim 2010), we applied the same analysis to our previous simulation of this scenario (G09a). We took the case of a low ionizing flux impinging onto a marginally stable Bonnor–Ebert Sphere. The low flux case was taken to allow for more moderate velocities. This scenario results in an elongated feature depicted in Figure 11 at t = 600 kyr. The corresponding LOS profiles are given in Figure 12. Here, the profiles differ significantly from the observations as well as from Figure 9. First of all, there is a double peak. This shows that the structure is produced by the collision of two shock fronts, one moving away from us and one moving toward us. These two shock fronts, which are encompassing the original density enhancements, can be seen directly colliding in the perpendicular projection (Figure 11(a)).

Figure 11.

Figure 11. Main structure in the RDI model (G09a). Color-coded is the surface density, each figure is 0.4 pc × 0.8 pc. (a) x–y projection, (b) x–z projection. The superimposed grid denotes the bins along which the LOS velocities in Figure 12 are taken.

Standard image High-resolution image
Figure 12.

Figure 12. LOS velocities along the main structure in the RDI model (G09a). The LOS velocities are taken in bins of the size 0.04 pc × 0.04 pc (see Figure 11) to match the observational resolution. Neither the overall pattern nor the line width is in agreement with the observations (Gahm et al. 2006).

Standard image High-resolution image

Even if the two peaks become indistinguishable, the lines are much broader. Furthermore, there is no detectable density gradient or rotational pattern visible. In addition, the veil seen between the two smaller pillars in M16 (Hester et al. 1996) poses a real challenge. Even the sophisticated simulations of Mackey & Lim (2010) do not reproduce it. In our simulations (see, e.g., Figure 6) veils arise naturally due to the turbulent motions. This is a strong indication that the pillars in M16 or in NGC 7822 are produced by the interplay of pre-existing turbulent structures and ionizing radiation.

6. CONCLUSIONS AND DISCUSSION

With iVINE, a tree-SPH code including ionization, we perform a parameter study on the formation of pillar-like structures and triggered star formation in H ii regions. First, we show that our simulations are converged and the choice of boundary conditions does not affect the outcome significantly (Section 3.2). After that, we show that ionizing radiation imposed on a turbulent molecular cloud can result in the formation of pillar-like structures which resemble observed pillars in size, mass, density as well as velocity structure. The rotational pattern observed in several pillars especially indicates that these were formed by the ionization of a turbulent cloud and not by the RDI of pre-existing clumps.

We thus conclude the following.

  • 1.  
    Varying the turbulent Mach number between 1.5 and 12.5 changes the morphology. If the Mach number is too low, there are no pre-existing structures that can be enhanced, if it is too high, the structures disperse quickly. The most favorable regime for the formation of pillars is Mach 4–10. However, the physical quantities such as the mass assembled and the column density are only weekly dependent on the Mach number (see Section 3.3).
  • 2.  
    The formation of pillar-like structures in the case with a favorable Mach number critically depends on the initial density contrast. Structures and therefore stars only form if the density contrast is lower than the temperature contrast between the hot and the cold gas: $\Delta \rho _\mathrm{init}\le \frac{2T_\mathrm{ion}}{T_\mathrm{nion}}$ (see Section 3.4).
  • 3.  
    The evolution of the ionized mass, density, and the front position in a turbulent medium under the influence of different initial fluxes is remarkably similar to the evolution in a homogeneous case. Thus, the size and density inside an H ii region solely depend on the initial flux, density, and the time since the ignition of a star or the distance from the star and globally follows the simplified analytical prescription (see Section 3.5).
  • 4.  
    The density of the resulting pillars is determined by a pressure equilibrium between the hot and the cold gas. Therefore, the expected density of the structures can be calculated as well (see Section 3.1).
  • 5.  
    The size of the evolving structures critically depends on the driving modes of the turbulence. Smaller driving modes lead to smaller structures. In our simulations, the relation is roughly dpilxdriving/40 (see Section 3.6).
  • 6.  
    Core and star formation are likely to occur. The higher the mass in the structures and the higher the initial flux, the earlier cores form (see Section 3.7).

Combining (3) and (4) allows us directly to determine the density of the forming structures as function of the initial mean density of the medium, the flux of the source, and the time since the ignition of the source or the position of the ionization front (see Equation (18)).

One has to keep in mind that our approach is simplified. First of all, no scattering is taken into account. Once an electron recombines, the emitted photon is assumed to be absorbed in the direct neighborhood (on the spot approximation). Thus, the reheating of shadowed regions by the adjacent hot gas is not taken into account. How much this affects the formation of pillars is the topic of ongoing research (B. Ercolano & M. Gritschneder 2010, in preparation). In addition, we focus on atomic hydrogen only, which makes it impossible to follow the precise temperature evolution as well as the PDRs. On the other hand, our simulations indicate that the pillars are in pressure equilibrium with the hot gas. Therefore, the PDRs might be transition regions comparable to a thin shock layer which is not resolved. Furthermore, we do not take magnetic fields into account. These might have implications on the global shapes of the H ii region (see, e.g., Krumholz et al. 2007). Nevertheless, we are able to reproduce the corkscrew morphologies in the pillars which were sometimes attributed to magnetic fields (see Section 5.2). Another aspect we neglect is stellar winds. Although there is clear observational evidence the winds of a massive star interact with the surrounding ISM (e.g., Westmoquette et al. 2010), it is up to now still unclear how effectively they affect its surroundings. From our simulations, we would estimate that stellar winds are of minor importance, maybe mainly enhancing the shock front as soon as a lower density in the hot gas allows for the effective driving of a stellar wind.

Altogether, our simulations are able to reproduce even the detailed fine structure of the pillars that have been observed with high resolution. In addition, we find that the observed LOS profiles allow for a clear distinction between the RDI scenario and the "radiative round-up" presented here. Current observations are in favor of our "radiative round-up" mechanism. Besides, our simulations give a deeper insight on the tight correlation between the parental molecular clouds size, density and turbulence, and the structures excavated by the ionizing radiation. The ionization acts as a magnifying glass, revealing the condition of the molecular cloud previous to the ignition of the massive star.

We thank the referee for valuable comments which helped to improve the manuscript. This research was funded by the DFG cluster of excellence "Origin and Structure of the Universe." M.G. acknowledges additional funding by the China National Postdoc Fund Grant No. 20100470108 and the National Science Foundation of China Grant No. 11003001. All simulations were performed on an SGI Altix 3700 Bx2 supercomputer that was partly funded by the DFG cluster of excellence "Origin and Structure of the Universe."

Footnotes

  • For a more detailed investigation of the role of compressional and solenoidal modes, see Federrath et al. (2010).

  • In low density, Mach 7, and small box, the second tip was taken to produce comparable results.

  • Since we only take the unionized gas inside the surrounding the actual volume of the pillar changes from simulation to simulation.

  • 10 

    Note that we always give the real density ρ of the hot gas, not the number density n to avoid the factor of μ = 0.5 when comparing the low-density, unionized gas to the ionized gas.

  • 11 

    In fact, the value is always slightly above one, but this can be attributed to the complete neglection of the turbulent motion in the cold gas in Equation (14).

  • 12 

    In fact the value for low density from Table 2 does not match precisely, but from Figure 3 the factor of ≈4 between Mach 5 (G09b) and low density can be seen.

  • 13 

    The masses do not differ significantly, as we do not follow the further accretion process and at the moment of formation the cores are still similar.

  • 14 

    Both statements have been verified by continuing this simulation until t = 600 kyr.

  • 15 

    The exciting source of the Trifid Nebula is HD 164492A an O7V star (Lynds et al. 1985), so the UV flux is comparable to our simulations.

  • 16 

    The structures are in the "Extended Ridge" which is roughly 10 pc away from the several O stars (O4–O9) of NGC 2244. Thus, the situation there is as well comparable to our simulations.

Please wait… references are loading.
10.1088/0004-637X/723/2/971