## Abstract

The effects of initially uniform magnetic fields on the formation and evolution of dense pillars and cometary globules at the boundaries of H ii regions are investigated using 3D radiation-magnetohydrodynamics simulations. It is shown, in agreement with previous work, that a strong initial magnetic field is required to significantly alter the non-magnetized dynamics because the energy input from photoionization is so large that it remains the dominant driver of the dynamics in most situations. Additionally, it is found that for weak and medium field strengths an initially perpendicular field is swept into alignment with the pillar during its dynamical evolution, matching magnetic field observations of the ‘Pillars of Creation’ in M16 and also some cometary globules. A strong perpendicular magnetic field remains in its initial configuration and also confines the photoevaporation flow into a bar-shaped dense ionized ribbon which partially shields the ionization front and would be readily observable in recombination lines. A simple analytic model is presented to explain the properties of this bright linear structure. These results show that magnetic field strengths in star-forming regions can in principle be significantly constrained by the morphology of structures which form at the borders of H ii regions.

## 1 INTRODUCTION

The study of elephant trunks, pillars and globules, found at the borders of H ii regions around massive stars, has received significant attention in recent years, both from an observational perspective and in theoretical and computational models. The famous ‘Pillars of Creation’ in the Eagle Nebula M16 were observed at optical wavelengths with *HST* (Hester et al. 1996), in the IR (Indebetouw et al. 2007; Sugitani et al. 2007) and in sub-mm/radio (e.g. Pound 1998; White et al. 1999), showing that these are dynamic structures with ongoing star formation which may or may not have been triggered by the radiation which has shaped the pillars. Smith et al. (2010b) provide strong evidence that pillars in the Carina nebula are significant sites of sequential star formation propagating away from the older star clusters in this region, building on previous observations of synchronized star formation around the periphery of the nebula by Smith & Brooks (2007). On smaller scales, studies of T-Tauri star ages in the Orion nebula (Lee et al. 2005; Lee & Chen 2007) show decreasing stellar ages moving away from massive stars and towards bright-rimmed clouds at the H ii region/molecular cloud interface, again strongly suggesting at least sequential and possibly triggered star formation. The clear relationship between pillars/globules and second-generation star formation around OB associations, and the question of the extent to which this star formation is triggered, provides strong motivation to understand the formation and evolution of these structures.

In a previous paper (Mackey & Lim 2010,, hereafter ML10) we investigated the formation and evolution of dense pillars of gas and dust – elephant trunks – on the boundaries of H ii regions using 3D hydrodynamical simulations including photoionizing radiative transfer (R-HD). It was found that shadowing of ionizing radiation by an inhomogeneous density field naturally forms elephant trunks without the assistance of self-gravity, or of ionization front and cooling instabilities. A combination of radiation-driven implosion (RDI; Bertoldi 1989) and acceleration due to the rocket effect (Oort & Spitzer 1955) produce elongated structures: RDI compresses neutral gas until pressure equilibrium with ionized gas is achieved; the rocket effect accelerates gas away from the radiation source producing dynamic elongated structures with lifetimes of a few hundred kyr (depending on clump masses/densities). Strong neutral gas cooling was found to enhance this formation mechanism, producing denser and longer-lived pillar-like structures compared to models with weak cooling.

Models such as these for the formation of bright-rimmed clouds, globules and pillars have been considered for many years (e.g. Pottasch 1958; Marsh 1970; Bodenheimer, Tenorio-Tagle & Yorke 1979; Sandford, Whitaker & Klein 1982); much of this work is summarized by Yorke (1986). The RDI of a photoionized clump and its subsequent acceleration and evolution was calculated analytically (Bertoldi 1989; Bertoldi & McKee 1990) and subsequently numerically by Lefloch & Lazareff (1994). Williams, Ward-Thompson & Whitworth (2001) considered a range of axisymmetric models showing that multiple scenarios could form long-lived pillar-like structures. It has been shown (Kessel-Deynet & Burkert 2003; Miao et al. 2006; Pound et al. 2007; Bisbas et al. 2010) that RDI of single clumps can generate cometary globules and trigger gravitational collapse, but it is more difficult to form pillars like those in M16 because gas must accumulate to a high density in the shadowed tail region. Lim & Mellema (2003) showed how photoionization of multiple clumps which partially shadow each other leads to dense gas accumulating in shadowed tail regions; it was suggested by Pound et al. (2007) that multiple clumps are required to quantitatively match observations of the pillars in M16. Recent models (Mellema et al. 2006a; Gritschneder et al. 2009, 2010; Lora, Raga & Esquivel 2009; Raga et al. 2009; ML10) showed that elephant trunks can form quite naturally from the photoionization of a clumpy density field under a range of initial conditions. Gritschneder et al. (2010) extended our previous work (ML10) which used static initial conditions by considering the photoionization of dynamic density fields generated by isothermal decaying turbulence, measuring the formation of pillar-like structures as a function of turbulent Mach number.

Filamentary structure in an apparently helical geometry was found in the shadowed tail regions of a number of elephant trunks (Carlqvist, Kristen & Gahm 1998; Carlqvist, Gahm & Kristen 2002, 2003), and it was suggested that this could arise due to twisting of magnetic field lines. Velocity profiles for some trunks were shown to be consistent with solid-body rotation (Gahm et al. 2006), again consistent with a magnetic origin of the observed structure, although the actual magnetic field orientation and strength has not yet been measured for these structures. The magnetic field in cometary globules has been measured by optical polarimetry (Sridharan, Bhatt & Rajagopal 1996; Bhatt 1999; Bhatt, Maheswar & Manoj 2004), showing in two cases a field orientation along the head–tail axis of the globule, but in one case (Bhatt 1999) a perpendicular field was found. Sugitani et al. (2007) used near-IR polarimetry of background stars to measure the magnetic field in M16, finding an ordered large-scale field in the H ii region, but within the pillars the field is aligned with the pillar axes and significantly misaligned with the ambient field by θ∼ 30°–40°. They suggest this should constrain the magnetic field strength since it has not been strong enough to resist reorientation during the formation and evolution of the pillars, a suggestion we explore in more detail in this work.

Theoretical calculations of the effects of magnetic fields on the expansion of H ii regions were first considered by Lasker (1966b). Using analytic calculations Redman et al. (1998) studied ionization fronts with a plane-parallel magnetic field in the plane of the ionization front. They found that the velocity separation between R-type and D-type solutions decreases with increasing field strength, and that D-critical ionization fronts also advance into the neutral gas more rapidly for increasing field strength. In Williams, Dyson & Hartquist (2000) jump conditions for ionization fronts with oblique magnetic fields were presented, together with 1D numerical models showing how an ionization front could progress from fast-R-type through fast-D-type, slow-R-type and finally to slow-D-type. These extra modes are allowed because fast and slow shocks detach from the ionization front at different times and propagate into the neutral medium. This work was extended by Williams & Dyson (2001), who calculated the internal structures of stationary 1D ionization fronts. They found that oblique fields could produce significant transverse velocities and regions where the transverse field component is significantly modified.

Williams (2007) used 2D slab-symmetric radiation-magnetohydrodynamics (R-MHD) simulations to study the photoevaporation flows from magnetized globules. In these models clumps with an initial density of *n*_{H} = 2 × 10^{5} cm^{−3} were placed in an ambient medium 10× less dense. For simplicity a uniform magnetic field with various strengths and orientations was used. Plane-parallel radiation was assumed, with a thermal model where the temperature relaxed to a value between 100 and 10 000 K according to its ionization fraction. It was found that for a weak field the ionized gas pressure dominates the dynamics and the field was swept into a configuration where it was parallel to the column of neutral gas behind the dense clump. For a sufficiently strong field, however, the field determined the dynamics and made the flow almost one-dimensional along field lines.

The first 3D R-MHD calculation including non-equilibrium photoionizing radiative transfer was performed by Krumholz, Stone & Gardiner (2007). They used the Athena code with a new ray-tracing module to simulate the expansion of an H ii region around a point source in a uniform magnetic field. The overall expansion is now axisymmetric rather than spherically symmetric (cf. Lasker 1966a) with a dense shell forming in directions parallel to the field, and a possibly numerical instability developing for expansion perpendicular to the field.

Using similar methods, Henney et al. (2009) model the photoionization of a dense clump of gas in 3D with an initially uniform magnetic field. They use the photon-conserving C^{2}-ray ray-tracing algorithm (Mellema et al. 2006b) which allows large time-steps to be taken without loss of accuracy. They found that the evolution of a photoionized globule can be significantly altered by the presence of a strong field, and in some cases a recombining shell forms at the termination shock of the photoevaporation flow when it is confined by a transverse field. The implosion phase is strongly asymmetric since the clump compresses much more readily along field lines than across them. These authors introduce a detailed thermal model to model the dynamics as realistically as possible for conditions in the Orion nebula, finding that X-ray heating keeps the neutral gas temperature at *T*≳ 50 K, but the cooling in neutral gas is somewhat stronger than that considered in ML10.

In this work we add magnetic fields of various strengths and orientations to some of the models considered in ML10 to study their effects on the dynamics of the dense neutral gas. The simulation code is described in Section 2: Section 2.1 describes the MHD dynamics algorithm; Section 2.2 reviews the microphysical processes included; Section 2.3 presents a brief comparison of results obtained with our code to the test problem described by Krumholz et al. (2007). The 3D simulations presented here are introduced in Section 3 and our results are presented in Section 4. Section 4.1 shows the evolution of the projected gas density and magnetic field orientation; Section 4.2 shows the emission from ionized gas; Section 4.3 evaluates boundary effects in strong-field simulations. A calculation explaining the presence of a bright ionized linear ridge/ribbon in the strong-field simulations is presented in Section 5. These results are compared to observations of H ii regions and their magnetic fields in Section 6, where we also discuss our work in the context of other recent computational research. Our conclusions are presented in Section 7, and some technical details and tests for the simulations are given in the appendices.

## 2 NUMERICAL METHODS

### 2.1 MHD algorithms and tests

These calculations are performed on a uniform finite-volume grid, with the hydrodynamics and ray-tracing/microphysics modules as previously described in ML10. The simulations presented here use the same code, but with an MHD module coupled to the ray-tracing. The integration scheme was described by Falle, Komissarov & Joarder (1998) and is second-order accurate in time and space, using the symmetric van Albada slope limiter to ensure monotonicity near shocks (van Albada, van Leer & Roberts 1982). To calculate the MHD fluxes a number of Riemann solvers were tried; we use the Roe solver in conserved variables described by Cargo & Gallice (1997) as it was the most robust for these calculations (with corrections to typos obtained by comparison with Stone et al. 2008). Magnetic field units are used in the code such that factors of 4π do not appear in the dynamical equations (see e.g. Falle et al. 1998); cgs units can be obtained by scaling the field strength by . To avoid confusion we will always quote field strength in cgs units (usually micro-Gauss, μG) rather than the code units. We use notation such that the total energy per unit volume is , being the sum of the kinetic, internal (thermal) and magnetic energies. The gas pressure is *p*_{g}, magnetic pressure is *p*_{m}=*B*^{2}/8π and the ideal gas adiabatic index γ is set to 5/3. Multidimensional calculations use the artificial diffusion prescription of Falle et al. (1998) to fix the ‘carbuncle problem’ and to prevent rarefaction shocks developing; typically a viscosity parameter η_{v} = 0.1 is used.

The code uses cell-centred values for vector quantities; the divergence constraint in the magnetic field is addressed using the ‘mixed-GLM’ divergence cleaning method of Dedner et al. (2002). This algorithm introduces an extra scalar field which couples to ∇ · ** B**, advecting and diffusing divergence errors so that they do not build up near shocks or other discontinuities. While it does allow non-zero divergence of

*(which can introduce small force errors), the method has some advantages. (1) It is fully conservative in all of the physical variables; only the evolution equation for the unphysical scalar field, ψ, contains a source term. (2) The total energy and magnetic field updates are calculated in the same algorithmic step, ensuring they are consistent. An energy correction is still required, but this can be calculated from the fluxes generated by divergence cleaning. This means the internal energy (as a fraction of the total energy) is treated quite accurately, making the code robust for a wide range of field strengths. (3) It involves relatively little computational and memory overhead.*

**B**A number of 1D and multidimensional adiabatic MHD code tests have been performed to validate the code and for comparison to other astrophysical MHD codes. Results from a number of these test problems are shown in Appendix A and can also be found in more detail at http://www.astro.uni-bonn.de/~jmackey/jmac/ together with HD and photoionization test problems (some of which were presented in ML10).

Three modifications have been made to the mixed-GLM divergence cleaning method. (1) The parameter *c _{r}* is set to

*c*= 4Δ

_{r}*x*(where the grid cell diameter is Δ

*x*) instead of the recommended value of

*c*= 0.18 (see Dedner et al. 2002, equation 45). A typo in their paper (A. Dedner, private communication) meant

_{r}*c*should have been defined as

_{r}*c*=Δ

_{r}*x*/0.18, similar to the value used here. (2) The algorithm induces extra magnetic field transport across cell boundaries; the energy associated with this is accounted for by an extra term in the energy flux:

*F*[

*X*] represents the flux of conserved variable

*X*across a cell boundary,

*B*is the component of

_{x}*normal to the cell boundary and*

**B***B*

_{x,m}is the value of

*B*in the interface state calculated using the mixed-GLM algorithm (Dedner et al. 2002, equation 42). This correction removed gas pressure dips which appeared ahead of oblique shocks in axisymmetric models of magnetized jets. (3) Motivated by numerical difficulties encountered for subsonic flow across boundaries in magnetically dominated media, the boundary condition for the scalar field ψ in the mixed-GLM method is modified for zero-gradient boundaries. This is described in Appendix B showing the improvement obtained with the new boundary condition.

_{x}### 2.2 Coupling MHD to photoionization

We briefly review the microphysics algorithm here; the reader is referred to ML10 for further details. The dynamics and microphysics are solved separately in turn by operator splitting: in a given time-step the adiabatic MHD equations are solved and updated first (with the ion fraction advected with the gas flow), and then the internal energy per unit mass (*e*_{int}) and ion fraction of hydrogen (*x*) are updated by an adaptive substepping integration to a user-specified relative accuracy. In the absence of photoionization this uses a fifth-order Runge–Kutta integration. Ray-tracing uses the short characteristics method with the photon-conserving discretization of the photoionization rate given by Mellema et al. (2006b). When a cell is strongly photoionized explicit integration is very unstable and so the ‘constant electron density’ approximation of Schmidt-Voigt & Koeppen (1987) is used (see also Mellema et al. 2006b) to analytically integrate the ion fraction when *x* > 0.95. The time-averaged fraction of photons which pass through the cell is calculated simultaneously to obtain the appropriate cell optical depth for calculating rays to more distant grid cells. The algorithm used in ML10 was designed to work in an identical way for R-MHD as for R-HD. As such, the non-dynamical test problems presented in ML10 produce identical results when run with a non-zero magnetic field.

We continue to use the on-the-spot approximation for diffuse radiation; this remains a potential limitation of our results which should be checked with future work. Recent calculations including diffuse radiation (Raga et al. 2009; Williams & Henney 2009) have shown that its role in the dynamics of H ii regions may not be as significant as was suggested by e.g. Ritzerveld (2005). We plan to test this in dynamical simulations using an algorithm such as that described by Kuiper et al. (2010); it should be possible to extend their method to treat ionizing photons.

Radiative cooling processes are modelled using the C2 cooling function from ML10 – in photoionized gas the coolants considered are collisionally excited emission from oxygen and nitrogen (Osterbrock 1989) and recombining hydrogen (Hummer 1994), and neutral gas cooling is assumed to follow Newton's Law of exponential cooling with a cooling time *t*_{c} = 10 kyr. Photoionization assumes a monochromatic ionizing source with 5.0 eV of thermal energy added to the gas per ionization. Collisional ionization is calculated using the rate of Voronov (1997).

### 2.3 H ii region expansion with a magnetic field

Krumholz et al. (2007, hereafter KSG07) added a photoionization module to the Athena code to calculate the 3D R-MHD expansion of an H ii region into a uniform medium with an initially uniform magnetic field. We have modelled this problem both in axisymmetry and in 3D to compare results from our code to those of KSG07. Our 3D model used a smaller domain of diameter 14 pc resolved by 160^{3} cells (12 per cent larger cell diameters than KSG07). The only other difference in our simulation was to add a low level of noise to the otherwise uniform initial density and gas pressure fields. An initial ambient temperature of *T* = 100 K was used rather than the 11 K used by KSG07 and we used the C2 cooling function described above, but tests showed this does not affect the results significantly.

Gas density and velocity field are shown for a slice through the *x**y* plane at *t*≃ 1.58 Myr in Fig. 1. Comparison to fig. 18 of KSG07 shows that very similar results are found here, notably the dense shell which forms for expansion in the directions and its absence in perpendicular directions. We also find the instabilities which were noted by KSG07 in the perpendicular direction; we suggest they are real (i.e. not numerical artefacts) and related to the R-HD ionization front instabilities studied by previous authors (e.g. Garcia-Segura & Franco 1996; Williams 2002; Whalen & Norman 2008) on the basis of tests which showed the instabilities are only present when neutral gas can cool strongly. The velocity arrows show that the magnetic field significantly alters the ionized gas dynamics: at late times the strong outflow along field lines transforms the ionization front to a recombination front in the directions. We plan to study these features in more detail in future work.

## 3 3D R-MHD SIMULATIONS

To focus specifically on the effects of magnetic fields on the photoionization process, the dense clump configurations in models 17 and 18 from ML10 were used rather than the random clump distributions in models 1–16. The two clump configurations as well as the simulation domain and radiation source properties are described in Table 1. Configuration 1 is basically the same as model 17 in ML10; configuration 2 is changed slightly – the two front clumps are slightly larger and further apart, and an extra clump is added further from the source. The radiation source located at the origin has the same properties as before: an ionizing photon luminosity *L*_{γ} = 2 × 10^{50} s^{−1} with monochromatic photon energy *h*ν_{0}− 13.6 eV = 5.0 eV. The simulation domain used here is twice as large in every dimension as was used for models 17 and 18 in ML10; this is because boundary effects were found to be much more significant with strong magnetic fields (cf. Henney et al. 2009). The simulation domain contains 384 × 256^{2} cells, giving a physical resolution of Δ*x*≃ 0.012 pc per cell. Results of 3D R-MHD simulations with various domain sizes are described in Appendix C and demonstrate that the domain used here is sufficient to prevent significant boundary effects for all but the strongest fields modelled.

Config. | Object | x | y | z | δ_{0} | r_{0} | M |

1,2 | X_{min} | 1.5 | −1.5 | −1.5 | n/a | n/a | n/a |

1,2 | X_{max} | 6.0 | 1.5 | 1.5 | n/a | n/a | n/a |

1,2 | Source | 0 | 0 | 0 | n/a | n/a | n/a |

1 | Cl 1 | 2.30 | 0 | 0 | 500 | 0.09 | 28.4 |

1 | Cl 2 | 2.75 | 0 | 0.12 | 500 | 0.09 | 28.4 |

1 | Cl 3 | 3.20 | 0 | −0.12 | 500 | 0.09 | 28.4 |

2 | Cl 1 | 2.30 | 0 | −0.174 | 250 | 0.12 | 33.7 |

2 | Cl 2 | 2.30 | 0 | 0.174 | 250 | 0.12 | 33.7 |

2 | Cl 3 | 2.73 | 0 | 0 | 500 | 0.09 | 28.4 |

2 | Cl 4 | 4.20 | 0.12 | 0.12 | 400 | 0.12 | 53.9 |

Config. | Object | x | y | z | δ_{0} | r_{0} | M |

1,2 | X_{min} | 1.5 | −1.5 | −1.5 | n/a | n/a | n/a |

1,2 | X_{max} | 6.0 | 1.5 | 1.5 | n/a | n/a | n/a |

1,2 | Source | 0 | 0 | 0 | n/a | n/a | n/a |

1 | Cl 1 | 2.30 | 0 | 0 | 500 | 0.09 | 28.4 |

1 | Cl 2 | 2.75 | 0 | 0.12 | 500 | 0.09 | 28.4 |

1 | Cl 3 | 3.20 | 0 | −0.12 | 500 | 0.09 | 28.4 |

2 | Cl 1 | 2.30 | 0 | −0.174 | 250 | 0.12 | 33.7 |

2 | Cl 2 | 2.30 | 0 | 0.174 | 250 | 0.12 | 33.7 |

2 | Cl 3 | 2.73 | 0 | 0 | 500 | 0.09 | 28.4 |

2 | Cl 4 | 4.20 | 0.12 | 0.12 | 400 | 0.12 | 53.9 |

These models were run with zero-gradient boundary conditions, which perfectly model supersonic outflow and damp reflected waves in subsonic outflow. There is nothing to stop inflow, however, if that is what the solution tends towards, and some of the R-MHD simulations do develop strong inflows at later times. An alternative boundary condition was imposed in the 2D slab-symmetric R-MHD simulations of Williams (2007): when the velocity was outward the usual zero-gradient boundary was applied, but if the velocity changed sign then inflow was suppressed. This ‘only-outflow’ boundary condition was also tested in our simulations to study how significantly inflows affected the solution, and in particular its observable properties. This was initially tested using axisymmetric simulations with a magnetic field parallel to the radiation propagation vector; the results showed that for a weak or medium field up to *B*∼ 50 μG the inflow made little difference. Its only effect was to confine the photoevaporation flow from the clump to a smaller volume by the inflow's ram pressure. The higher mean density of ionized gas (and associated stronger recombination) had no discernible dynamical effect. With a strong field of *B*≃ 160 μG, however, the inflow could overrun the photoevaporation flow and the ionization front, transforming it into a recombination front and hence strongly affecting the solution. For this reason the strong-field simulation R8 described below was run with both a zero-gradient (as R8) and an only-outflow (as R8a) boundary condition in the direction back towards the radiation source.

### 3.1 Description of simulations

The simulations are listed in Table 2. R1 and R11 are R-HD simulations of clump configurations 1 and 2, respectively, while the other simulations are R-MHD calculations of the same clump configurations with varying field strengths and orientations. The field strengths were approximately *B*≃ 18, 53 and 160 μG and are referred to as weak, medium and strong, respectively. R2–R10 are R-MHD simulations of configuration 1: R2–R4 are weak-field models, R5–R7 have a medium field and R8–R10 a strong field. R12–R15 are R-MHD simulations of configuration 2: R11–R14 have a weak field and R15 a medium field.

Model | Cl. | BCs | B _{x} | B _{y} | B_{z} | t_{sim} |

R1 | 1 | 1 | 0 | 0 | 0 | 800 |

R2 | 1 | 2 | 1.77 | 1.77 | 17.7 | 525 |

R3 | 1 | 1 | 3.1 | 13.4 | 11.2 | 300 |

R4 | 1 | 2 | 3.1 | 13.4 | 11.2 | 300 |

R5 | 1 | 1 | 0 | 0 | 53.2 | 800 |

R6 | 1 | 1 | 53.2 | 0 | 0 | 800 |

R7 | 1 | 1 | 9.23 | 40.1 | 33.7 | 450 |

R8 | 1 | 1 | 14.2 | 7.1 | 158.9 | 400 |

R8A | 1 | 3 | 14.2 | 7.1 | 158.9 | 400 |

R9 | 1 | 1 | 159.5 | 7.1 | 10.6 | 800 |

R10 | 1 | 1 | 27.7 | 120.3 | 101.0 | 280 |

R11 | 2 | 3 | 0 | 0 | 0 | 700 |

R12 | 2 | 3 | 1.77 | 1.77 | 17.7 | 375 |

R13 | 2 | 3 | 1.77 | 17.7 | 1.77 | 700 |

R14 | 2 | 3 | 17.7 | 1.77 | 1.77 | 700 |

R15 | 2 | 3 | 5.3 | 5.3 | 53.2 | 525 |

Model | Cl. | BCs | B _{x} | B _{y} | B_{z} | t_{sim} |

R1 | 1 | 1 | 0 | 0 | 0 | 800 |

R2 | 1 | 2 | 1.77 | 1.77 | 17.7 | 525 |

R3 | 1 | 1 | 3.1 | 13.4 | 11.2 | 300 |

R4 | 1 | 2 | 3.1 | 13.4 | 11.2 | 300 |

R5 | 1 | 1 | 0 | 0 | 53.2 | 800 |

R6 | 1 | 1 | 53.2 | 0 | 0 | 800 |

R7 | 1 | 1 | 9.23 | 40.1 | 33.7 | 450 |

R8 | 1 | 1 | 14.2 | 7.1 | 158.9 | 400 |

R8A | 1 | 3 | 14.2 | 7.1 | 158.9 | 400 |

R9 | 1 | 1 | 159.5 | 7.1 | 10.6 | 800 |

R10 | 1 | 1 | 27.7 | 120.3 | 101.0 | 280 |

R11 | 2 | 3 | 0 | 0 | 0 | 700 |

R12 | 2 | 3 | 1.77 | 1.77 | 17.7 | 375 |

R13 | 2 | 3 | 1.77 | 17.7 | 1.77 | 700 |

R14 | 2 | 3 | 17.7 | 1.77 | 1.77 | 700 |

R15 | 2 | 3 | 5.3 | 5.3 | 53.2 | 525 |

For the field orientation, ‘parallel’ denotes parallel to the radiation propagation vector along the centre of the simulation (), and ‘perpendicular’ denotes a field in the *y**z* plane. Field orientations were chosen to be either parallel or perpendicular, or oriented 80° from the *x*-axis and 50° from the *y*-axis in the *y**z* plane. The simulations with the latter field configuration (R3, R4, R7, R10) produced very similar results to the perpendicular models (R2, R5 and R8) for configuration 1 so they were not run for as long as the other models and were not run at all for configuration 2.

Whether a field is weak or strong is largely determined by its dynamical importance, set by the plasma parameter β≡ 8π*p*_{g}/*B*^{2}, the ratio of thermal to magnetic pressure. Since ionized gas in H ii regions is largely isothermal, as is dense molecular gas, the thermal pressure is approximately proportional to density. Hence a field of 50 μG could be strong or weak depending on whether the gas is ionized or neutral, and on its density. The plasma parameter is shown in Table 3 for a range of gas pressures encountered in the photoionization simulations: (1) the initial conditions have a constant pressure of *p*_{g} = 1.38 × 10^{−11} dyne cm^{−2} (or *p*_{g}/*k*_{B} = 10^{5} cm^{−3} K); (2) gas at the background density of *n*_{H} = 200 cm^{−3} and the ionized gas temperature of *T*≃ 8000 K has *p*_{g} = 4.42 × 10^{−10} dyne cm^{−2}; and (3) the peak pressure typically encountered in the simulations is *p*_{g}∼ 10^{−8} dyne cm^{−2}. For weak-field simulations the gas pressure clearly dominates and for the strong field the situation is reversed; for the medium field case the initial conditions are magnetically dominated but, once ionized, the gas pressure is larger. Thus it is expected that the weak-field results will largely follow the R-HD results, the strong-field simulations should show very different behaviour and the medium field models will lie somewhere in between. The gas pressures in these simulations are comparable to those estimated for the pillars and their environment in M16 (see the discussion in ML10), but it should be borne in mind that other massive star-forming regions can have significantly higher (or lower) gas pressures.

|| (μG) B | p_{g} = 1.38 × 10^{−11}, | 4.42 × 10^{−10} | 10^{−8} |

18 | 1.1 | 34 | 800 |

53 | 0.12 | 4.0 | 89 |

160 | 0.014 | 0.43 | 10 |

|| (μG) B | p_{g} = 1.38 × 10^{−11}, | 4.42 × 10^{−10} | 10^{−8} |

18 | 1.1 | 34 | 800 |

53 | 0.12 | 4.0 | 89 |

160 | 0.014 | 0.43 | 10 |

## 4 RESULTS

Three main observable consequences of the magnetic field are expected (cf. Williams 2007; Henney et al. 2009): (1) for weak magnetic fields the field orientation will be changed by the dynamics of the photoionization process; (2) when the field is sufficiently strong the density structure of the neutral gas will be significantly altered because gas can only move along field lines; RDI produces a sheet rather than an axisymmetric compression; (3) strong magnetic fields will confine the photoevaporation flow changing its observable properties. We first plot projections through the simulation domain showing column density of neutral gas and the projected magnetic field. Following this we show emission maps in recombination radiation. R2, R5 and R8 proved most useful for showing the effects of increasing field strength so we focus most of our analysis on these three simulations.

### 4.1 Projected density and magnetic field

The column density is calculated as in ML10 by integrating the neutral gas number density along the line of sight (LOS). Note that because we do not consider molecules explicitly, the column densities should be divided by 2 to give N(H_{2}). The projected magnetic field is more difficult to calculate since it must be a weighted integral, for example to calculate the polarisation of background starlight induced by aligned dust grains (cf. observations of M16 by Sugitani et al., 2007). For our integration we assume a constant gas-to-dust ratio and weight the integral by gas density. To allow for the possibility that grain alignment may be less effective at high densities (e.g. Goodman et al., 1995) we limit the density weighting to a maximum density, *n*_{max} = 2.5 × 10^{4} cm^{−3}; this is rather ad-hoc but it has a very limited effect on the projected field orientation as long as *n*_{max}≳ 10^{4} cm^{−3}. The projected field may be calculated by integrating the Stokes *Q* and *U* parameters along the LOS and subsequently recovering the field orientation using trigonometric relations (see Arthur et al., 2010). With the approximations described above the projected *Q* and *U* parameters of the magnetic field are given by

**B**

_{p}=[

*B*

_{1},

*B*

_{2}]. Technically this integration has units of G.cm

^{−2}but the normalisation was not considered because we only use the projected orientation of the field in this work and not its magnitude. In principle the relationship between polarisation and magnetic field is a very complex function of temperature, density, local radiation field and dust properties (e.g. Spitzer, 1998) but the microphysical processes included in our simulations are not sufficient to model this in any detail.

Results from the perpendicular field models R2, R5 and R8 are shown in Figs 2 and 3 for simulation times 100, 200, 300 and 400 kyr. The LOS is the simulation *y*-axis so the initial field in all three simulations is almost vertical and perpendicular to the LOS. At 100 kyr the leading dense clump is just past the point of maximum compression due to RDI. The partially shadowed clumps are being compressed asymmetrically and are not yet at maximum compression. At this stage the structure of the dense gas is very similar in all three models, although it should be noted that the low-density tail of the strong-field model (R8) is actually a thin sheet seen edge-on, whereas it is closer to axisymmetric in R2. Already the magnetic field has been significantly altered in R2 due to RDI; this is less apparent in R5 while the field in R8 is almost unaffected. By 200 kyr the clumps have re-expanded somewhat (limited because much of the compressional heating was radiated away) but the rocket effect has taken hold and the first clump has almost fully merged with the second. At this stage the magnetic field reorientation in R2 is more pronounced and in R5 is beginning to become significant. There are also differences in the morphology among the three models, quite slight between R2 and R5 but significant for R8. After 300 kyr the differences have become even more pronounced although the three clumps have basically merged to one in all three simulations. By 400 kyr R2 and R5 have become cometary globules rather than pillar-like structures. R5 is also fragmenting – the small protrusion at the top of the dense region eventually detaches from the main clump and is rapidly accelerated. R8 at this stage looks completely different from the other models. Low-density gas has recombined behind the very broad ionization front formed by gas expansion along the field lines. The field remains almost unchanged from its initial value even at this late stage.

The field evolution is shown in a more quantitative way in Fig. 4, where the ratio of the mean (volume averaged) parallel field 〈|*B _{x}*|〉 to perpendicular field is plotted as a function of time for both the full simulation domain and for only those cells with gas density

*n*

_{H}≥ 5000 cm

^{−3}. This ratio changes relatively little in the full domain because ionized gas takes up the overwhelming majority of the simulation volume. For the dense gas, however, a very strong effect is seen in the evolution as the field strength increases. The weak field in model R2 is swept into alignment with the pillar by the dynamics of RDI and the rocket effect, but a sufficiently strong field (here ≃160 μG in R8) is unaffected.

In the parallel field models R6 and R9, there was very little field evolution because the rocket effect acts in the same direction as the initial field orientation. The results obtained are comparable to the S00L model in Henney et al. (2009), except that for R9 a strong inflow from the boundary closest to the source overran the ionization front after *t*∼ 200 kyr, leading to a stable 1D recombination front. This would not be expected to form if an only-outflow boundary is imposed, as in Williams (2007), although we note that for 1D photoionization where ionized gas is confined to stay in the LOS back towards the source, a recombination front is a more stable and natural solution than an ionization front.

### 4.2 Emission maps in recombination lines

The gas emission in recombination radiation (e.g. Hα) was also calculated from the simulation outputs as in ML10, using the emissivity formula from Henney, Arthur & García-Díaz (2005) and dust opacity as in Mellema et al. (2006a) to attenuate radiation along a LOS. Emission maps are shown at times 100, 250 and 400 kyr for models R2, R5 and R8 in Figs 5–7. Here, as in the column density maps, there are small differences between R2 and R5, but R8 shows very different emission. This can be understood from Table 3 which shows that only R8 is magnetically dominated in ionized gas. The initial field is along the LOS for the left-hand projections (LOS is simulation *z*-axis) and vertical for the right-hand projections (LOS is simulation *y*-axis).

The field is so weak in R2 that the evolution is very similar to the purely hydrodynamical evolution in R1 (cf. ML10) except that the low-density shadowed tail region is not quite axisymmetric. A strong photoevaporation flow from the front clump expands spherically with velocity *v*≃ 30–35 km s^{−1} until the ram pressure equals the total pressure of the ambient medium at which point a standing shock is established. For the full evolution of this model the dense gas at the ionization front is by far the brightest structure in recombination emission. There are small differences in model R5, notably that the photoevaporation flow is confined to a smaller volume by the higher ambient total pressure, but the evolution is largely the same.

In R8 the total pressure of the ambient medium is magnetically dominated and much higher than in R2 or R5. This confines the photoevaporation flow much more effectively so that the region of spherical expansion is significantly smaller. Instead shocked photoionized gas is diverted along field lines into a bar-shaped linear structure with a density 2.5–4× the background density. In Henney et al. (2009) this structure is termed a *dense ribbon* to avoid confusion with photon-dominated regions such as the Orion Bar (see e.g. O'Dell 2001). We will avoid the term ‘bar’ for the same reason, instead referring to the structure as a ribbon or ridge (cf. Bohigas et al. 2004). This dense ionized gas is quite spatially extended and so has a comparable recombination radiation intensity to the bright rim at the ionization front. Because of the high column density in the ribbon there is significant recombination in its shadow at later times. This evolution is similar to that seen in the strong-field models of Henney et al. (2009).

The left-hand plots for R8 (with * B* along the LOS) show wave-like features which appear to be partially developed Kelvin–Helmholz instabilities. These move back towards the mid-plane (

*y*= 0) due to magnetic tension and have collided by

*t*= 400 kyr (Fig. 7); similar but more pronounced evolution was found by Henney et al. (2009) in their simulations.

In simulation R8 at 400 kyr the dense ridge/ribbon has become optically thick, enhancing recombination in its shadow and dramatically weakening the photoevaporation flow. The result of this is that the standoff distance between the ridge and the original ionization front of the dense clump has decreased almost to zero. This was found to be a short-lived phase in the models of Henney et al. (2009) and we expect the same to be true here because a strong photoevaporation flow is required to generate the overdense ionized gas.

### 4.3 Boundary effects for simulation R8

It was found that a reasonably strong inflow of gas developed at the boundary *x* = 1.5 pc in R8 at late times, with velocity *v _{x}* = 1–3.5 km s

^{−1}. The ram pressure from this inflow could raise the density in the shocked region and hence affect the resulting evolution, so the simulation was repeated with this boundary set to only-outflow, thereby preventing the inflow developing. This simulation, denoted by R8a, is shown in Fig. 8 for times 100, 250 and 400 kyr, again in recombination line emission. Comparison to Figs 5 and 6 shows that the evolution is almost identical up to

*t*= 250 kyr except that the bright ridge is slightly farther from the ionization front because now there is no ram-pressure confinement. At

*t*= 400 kyr the differences are more significant (see Fig. 7): the ionized ridge is broader and less well defined in R8a than in R8, and while the dense neutral gas has a similar structure there are small differences. The general features remain the same, however, and the confined photoevaporation flow still shields the ionization front significantly and outshines it in recombination radiation.

A similar but weaker inflow was also found in simulation R5. The only effect this has on the results is that the ionized ambient medium is denser than in R2, as seen by comparing the mean emission between the two simulations in Figs 5–7.

## 5 IONIZED RIDGE FORMATION

The presence of a bright ionized ridge/ribbon of dense gas in simulation R8 which is absent in R1, R2 and R5 is easily understood quantitatively by evaluating the jump conditions for the isothermal termination/reverse shock of the photoevaporation flow (the temperature is always in the range of 7500–8500 K in the ionized gas). If we assume that the system has reached an equilibrium, then the total pressure in the upstream (I), shocked (II) and downstream (III) regions will be equal. To avoid dealing with the spherical expansion of the photoevaporation flow, we consider region I to consist only of the conditions immediately upstream from the shock (referred to below as subscript ‘1’). If we further assume that the only significant magnetic field is transverse to the velocity of expansion, we can consider only its pressure effects and define a speed *b*^{2}≡*B*^{2}/(8πρ) =*p*_{m}/ρ so that the total pressure in a region is

We wish to calculate the density in region II as a function only of the upstream velocity *v*_{1} and the ambient medium properties (*b*_{3}, *v*_{3}). Only when the ratio ρ_{2}/ρ_{3} > 1 will a bright ribbon or ridge be observable. Thus the equations to be solved are (cf. de Hoffmann & Teller 1950)

*k*

_{1}≡ρ

_{1}/ρ

_{3}and

*k*

_{2}≡ρ

_{2}/ρ

_{3}. The isothermal sound speed is

*c*

^{2}≡

*p*

_{g}/ρ=

*kT*/(μ

*m*

_{p}), where μ is the mean mass per particle in units of the proton mass

*m*

_{p}and is μ = 0.5 in ionized gas in our models. For

*T*≃ 8300 K this gives

*c*= 11.7 km s

^{−1}. In all simulations the magnetic field is negligible in region I (

*b*

_{1}≪

*c*) and smaller than the other terms in region II (

*b*

_{2}<

*c*). The flow velocity in the reverse shock reference frame is

*v*

_{1}≃ 33–35 km s

^{−1}in all simulations giving a Mach number

*M*≡

*v*

_{1}/

*c*≃ 3.

The hydrodynamic jump conditions for the isothermal reverse shock are

When*b*

_{1}≃ 0 but

*b*

_{2}> 0, the density jump is given by the solution to the quadratic equation This equation has real solutions only for some values of

*b*

_{2}, leading to the constraint For

*M*= 3 this gives

*b*

^{2}

_{2}≤ 16

*c*

^{2}/9 but in the simulations it is actually significantly less than this limiting value.

Using these equations it is easy to show that in the hydrodynamic limit *k*_{2}=*M*^{2}/(1 +*M*^{2}) ≃ 0.9 and *k*_{1} = 1/(1 +*M*^{2}) ≃ 10. No bright overdensity is formed, and the photoevaporation flow shocks at a density about 10× below the ambient density. This is indeed seen in the R-HD simulation R1 and also the weak perpendicular field simulation R2.

With *b*_{1}=*v*_{3} = 0 and *b*_{3}, *b*_{2} > 0 we obtain

*b*

_{2}given by equation (7). With

*M*= 3 this gives where β

_{3}≡

*p*

_{g}/

*p*

_{m}=

*c*

^{2}/

*b*

^{2}

_{3}. The overdensity in the post-shock region is therefore determined by the value of (1 + 1/β) in the ionized ambient medium, showing that the ambient medium must be magnetically dominated to give an overdense ridge. For R5 β

_{3}= 4.0 and for R8 β

_{3}= 0.43, leading to maximum overdensities

*k*

_{2}of 1.1 and 3.0, respectively. Actual post-shock densities in R5 correspond very well to this value, and in R8/R8a the ribbon density is

*n*

_{H}≃ 550–700 cm

^{−3}depending on the time. The ambient density is

*n*

_{H}= 200 cm

^{−3}so the overdensity of

*k*

_{2}= 2.5–3.5 is again close to the predicted value.

Similarly, it is easy to show that *k*_{1} = (1 + 1/β_{3})/(1 +*M*^{2}) giving *k*_{1} = 0.125 and 0.33 for R5 and R8. According to this equation the photoevaporation flow should shock for simulations R1, R2, R5 and R8 at *n*_{H} = 20, 21, 25, 67 cm^{−3}, respectively. In the simulations at *t* = 150 kyr the actual densities are *n*_{H}≃ 25, 22, 31, 75 cm^{−3}, whereas at *t* = 200 kyr they are *n*_{H}≃ 18, 19, 30, 100 cm^{−3}, in both cases comparable to the values predicted.

This simple model has obvious limitations, in particular we have not specified the form of the boundary between regions II and III – it is not a shock, but if the magnetic field strength changes significantly then a contact discontinuity forms to equalize the total pressure. Also, we have not considered any parallel field component which, if present, would alter the pressure balance. Despite this the model quantitatively captures the properties of the overdense ridge. Addition of an opposing ram pressure in the ambient medium can clearly have a similar effect on a downstream magnetic field in that it increases the density at which the photoevaporation flow termination shock forms (and hence also the post-shock density ρ_{2}). The differences between R8 and R8a (Figs 6–8) show that even an inflow of *v*_{3}≃ 3 km s^{−1} can have an observable effect. Further work will be required to investigate observational differences (possibly in the shape and velocity of gas within the ridge and/or whether it is sheet-like or bar-like) which could be used to distinguish between a magnetically confined bright ridge and one which is ram-pressure confined.

## 6 DISCUSSION

This work confirms the suggestion by Sugitani et al. (2007) that their observations can in principle constrain the field strength in M16 when combined with detailed simulations. On the evidence of these simulations (and those of Henney et al. 2009) it appears that the magnetic field in the region around the pillars in M16 cannot be as large as the 160 μG field in model R8/R8a, and is likely |** B**| ≲ 50 μG. This can be deduced by comparing observations of the field orientation in the H ii region and Hα emission around the pillars with the emission maps presented here and in Henney et al. (2009). The recombination line intensity in the images from Hester et al. (1996) clearly decreases with distance from the ionization front and was shown to be consistent with a spherically expanding photoevaporation flow, which matches only our weak and medium field strength models. Additionally, no prominent ribbon/ridge or sheet-like features are present near the pillars.

Interpreting the polarization measurements within the pillars in M16 as further evidence for a weak field implicitly assumes that the pillars formed via a mechanism similar to the one considered here. Because the rocket effect always accelerates gas away from the radiation source and into the shadowed region, we believe that this alignment of weak fields with the pillar axis is a generic feature of all such models, but the initial conditions for pillar formation are not static (cf. Gritschneder et al. 2010) and it is uncertain how well correlated the initial magnetic field orientation would be between dense and rarefied gas.

Compared to the simulations of Henney et al. (2009), the RDI of clumps in our models is weaker leading to less flattened structures. The main reason is that here the clumps are considerably more centrally concentrated with higher initial densities, making them less susceptible to compression. Additionally, the clump in their model is much closer to the ionizing source and subject to higher incident radiation flux (only partially offset by the higher source luminosity in our model), increasing the effect of RDI. Thermal physics models and differing geometric effects due to source proximity may also have an effect.

The weak-field model W80L in Henney et al. (2009) should be roughly comparable to our model R5. Fig. 4 shows that the ratio in dense neutral gas remains at ∼0.3 for 300 kyr in R5, similar to the ratio in neutral gas shown in Henney et al. (2009, fig. 12) for W80L. The subsequent increase in this ratio in R5 may represent a later evolutionary phase but, because of the simulation differences already mentioned, it is difficult to quantify this.

Our results are otherwise in good agreement with those of Henney et al. (2009) considering the differences in the simulation initial conditions. We also find that the photoionization of dense clumps is very different in strongly magnetized media than in the non-magnetized case, with the dynamics becoming more planar than axisymmetric. Photoevaporated ionized gas tends to form ribbon-like structures which can (transiently) have significant optical depths and lead to recombination behind them. Clump compression along field lines creates flattened sheet-like structures and an ionization front which is far from axisymmetric.

Despite the simplistic initial conditions used here and in ML10, qualified support for these models comes from the 3D R-HD simulations of Gritschneder et al. (2010). They used an isothermal self-gravitating decaying turbulence model as the initial conditions, varying the turbulent Mach number at which the radiation source is switched on. In agreement with the trends we reported they also found that to form pillars the density field was required to include both large low-density regions and reasonably massive dense regions of sufficient overdensity to trap the ionization front. Because of the dynamic initial conditions in their simulations they were also able to show the dependence on initial gas motions, finding that high Mach number turbulence did not produce pillar-like structures as successfully as models with Mach number *M*∼ 4–10 because structures formed and dispersed too rapidly. When the Mach number was too low a sufficient density contrast to form pillars was not obtained. Gritschneder et al. (2010) also found that ionized gas pressure is the dominant driver of the dynamics and gravity appears to play a smaller role, an explicit assumption in our work based on previous calculations by Esquivel & Raga (2007).

In our opinion Gritschneder et al. (2010) overstate the differences between our results and theirs – in both models a combination of RDI and the rocket effect reinforces pre-existing inhomogeneities in the ISM and generates elongated structure along the radiation propagation direction. We do not see our results as being in conflict with those of Gritschneder et al. (2010); rather that they have confirmed many of our conclusions and also extended our work by studying pillar formation in a more realistic initial density field. These results suggest that the initially static models considered in ML10 and in this work can capture the essential physics of the pillar formation and evolution, once the initial conditions of dense regions surrounded by a less dense ambient medium have been set up. As noted in ML10, our models can also explain the LOS velocities seen in the pillars in M16 (Pound 1998; White et al. 1999), although we do not see the helical morphology and apparent rotation found in some elephant trunks (Gahm et al. 2006). This may require an initially dynamic density field since such structures seem to occur quite readily in the simulations of Gritschneder et al. (2010).

Our simulations (and those of Henney et al. 2009) suggest that ionization fronts with a strong perpendicular magnetic field should have clear observational consequences in the morphology of both the recombination emission and the dense gas. Henney et al. (2009) note that while the Orion Bar has a linear shape similar to the ribbon-like structures seen here, it is associated with the main ionization/dissociation front itself rather than being an overdensity within the H ii region. The extended photoionization front and photon-dominated region seen edge-on in M17 was modelled by Pellegrini et al. (2007) as being supported by magnetic pressure (observations by Brogan & Troland 2001 suggest an LOS magnetic field *B* > 100 μG). Unfortunately, much of this H ii region is heavily obscured at optical wavelengths and so straightforward inspection of Hα images is not particularly revealing. In the extensive Hα survey of the Carina nebula by Smith, Bally & Walborn (2010a), only the elephant trunk in ‘Pos. 23’ of fig. 1 shows evidence of a bright ridge in front of the trunk, suggesting that such features are not common. A clearer example is seen in NGC 6357 (Bohigas et al. 2004) where a prominent bright ridge in front of an elephant trunk is the brightest Hα structure in the field of view. A similar structure is found in front of a massive pillar in NGC 3603 (Brandner et al. 2000). In both of these observations, however, the adjacent massive star cluster is expected to drive outflowing gas and so it is unclear if the bright ridges are ram-pressure confined or magnetically confined.

The larger simulation domains used here (compared to ML10) allow us to follow the evolution of the dense gas for significantly longer, up to 700 kyr in some simulations. We find that in the absence of dense gas further from the star, a pillar-like structure will be flattened into a cometary globule with a dense head and low-density tail. To study this in more detail, clump configuration 2 was set up with an extra clump a further parsec from the radiation source. It was found that the pillar-like stage of the evolution lasted ∼600 kyr in this model, until all of the dense gas reached the position of the furthest clump. Subsequently a cometary morphology developed again. This experiment shows that while the pillar's lifetime can be extended somewhat, it cannot survive indefinitely unless there is a long dense filament pointing away from the radiation source. This supports the general picture of elephant trunks shown in fig. 15 of Smith et al. (2010b) as structures which form and disperse over about 10^{6} yr. With our non-gravitating simulations we cannot model the star formation which also occurs, but the simulations by Gritschneder et al. (2010) and Bisbas et al. (2010) show that the compression induced by RDI produces gravitationally unstable fragments which would likely form stars. While it may seem unlikely that the three pillars in M16 should all have originated from elongated overdensities pointing back towards the ionizing stars, recent observations from the Herschel observatory (e.g. Molinari et al. 2010) have shown that molecular clouds appear to have a distinctly filamentary structure. Simulations of MHD turbulence generated by colliding flows (Banerjee et al. 2009) show similar structures, and simulations of small-scale non-ideal isothermal MHD turbulence (Downes & O'Sullivan 2009) also show significantly more linear structure compared to purely hydrodynamical models.

## 7 CONCLUSIONS

We have performed a series of R-MHD simulations of the photoionization of dense clumps of gas and their evolution from pillar-like to cometary globule-like structures. Our results for the emissivity of ionized gas agree very well with those of Henney et al. (2009), showing that a dense, ionized, bar-shaped region standing off from the ionization front is a generic feature of strongly magnetized photoionization in a clumpy medium for a perpendicular field orientation. This ridge can be as bright as, or even brighter than, the photoionization front when observed in recombination radiation (e.g. Hα) and its presence or absence can be used as a diagnostic of the strength of any large-scale magnetic field which may be present. Bright ridges or ribbons are observed in some H ii regions (e.g. Brandner et al. 2000; Bohigas et al. 2004; Smith et al. 2010a), although they are not common. An overdense ridge could also be produced by ram-pressure confinement, and more detailed modelling is required to find observational signatures which could distinguish these different confinement mechanisms.

Comparing to observations of M16 (Hester et al. 1996) there is no such ribbon or ridge, suggesting that the ambient field measured by Sugitani et al. (2007) is not dominant in the ionized gas. This conclusion is strengthened when we consider the magnetic field orientation observed in M16 (Sugitani et al. 2007, fig. 9). The results presented here show that both RDI and acceleration of clumps by the rocket effect tend to align the magnetic field in dense neutral gas with the radiation propagation direction. In our models a field configuration similar to the observed one is clearly seen when the initial field strength is |** B**| ≃ 20 μG; the simulation with |

**| ≃ 50 μG is consistent with observations, and the simulation with |**

*B***| ≃ 160 μG is not consistent. Our simulations thus suggest an ambient field strength of |**

*B***| ≲ 50 μG around the M16 pillars.**

*B*The morphology of the structures which develop due to RDI and the rocket effect is also affected by a strong magnetic field, partly due to shielding by the dense ionized ridge and partly by the effect of the field within the pillar or globule. An inspection of Hα images of elephant trunks and globules in the literature (e.g. Hester et al. 1996; Smith et al. 2010a) suggests that the features seen in the strong-field simulations are not common, although the uniform initial field configurations considered here are certainly somewhat artificial. Additionally, many H ii regions have significantly higher gas pressure than that modelled in our simulations, in which case the magnetic field must also be correspondingly stronger to dominate the dynamics.

Comparing these results with our earlier simulations in ML10, the larger simulation domains used here show that the pillar-like structures which form will ultimately evolve to cometary structures in the absence of dense gas further from the star. The lifetimes of pillars in our models are *t*≲ 500 kyr, although this depends significantly on the initial mass and concentration (and presumably velocity; cf. Gritschneder et al. 2010) of the dense gas clumps.

Finally we emphasize, in agreement with previous authors (Krumholz et al. 2007; Williams 2007; Henney et al. 2009), that a strong magnetic field has a very significant influence on the dynamics of the photoionization process, and many of these effects should be easily observable. Given the difficulty of measuring the full 3D magnetic field in the ISM, a comparison to the detailed numerical simulations such as these offers an indirect means to constrain the field strength and orientation in and around H ii regions.

JM's work has been partly funded by the Irish Research Council for Science, Engineering and Technology; also by a grant from the Dublin Institute for Advanced Studies, and by Science Foundation Ireland. AJL's work was funded by a Schrödinger Fellowship from the Dublin Institute for Advanced Studies. JM acknowledges support from an Argelander Fellowship during the writing of this paper. Figures were generated using the VisIt visualization tool. The authors wish to acknowledge the SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities and support. We thank the referee for useful suggestions and for pointing out an error in an earlier draft.

## REFERENCES

*et al.*,

*et al.*,

*et al.*,

## Appendices

### APPENDIX A: TESTS OF ADIABATIC MHD ALGORITHMS

Results from the following MHD test problems can be found at http://www.astro.uni-bonn.de/~jmackey/jmac/ together with the HD test problems. These adiabatic calculations are run using dimensionless units. The shock-tube tests of Brio & Wu (1988) and Falle et al. (1998) were performed in 1D and on a 2D grid with the shock-normal at an angle θ = tan^{−1} 0.5 to the grid axes, ensuring that the problem is genuinely multidimensional. 2D Results are shown for the Brio & Wu (1988) test in Fig. A1 and are comparable to those obtained by Falle et al. (1998). The propagation of a polarized Alfvén wave was also calculated in 2D, using the initial conditions from Tóth (2000) and set up on a grid as in Stone et al. (2008). Running this test at different resolutions shows the degree of numerical diffusion in the algorithm and also the rate of convergence of the solution. Results are shown in Fig. A2; comparison of the first panel to results from Stone et al. (2008) shows that the third-order Athena algorithm is significantly more accurate for this test at a given resolution, as expected. The second plot in this figure shows the L1 error:

*N*cells, where φ

^{0}

_{i}is the reference state for each cell and φ

_{i}is the approximate numerical solution. This verifies that the code is indeed second-order accurate for continuously differentiable density fields.

Results for the Orszag–Tang vortex problem (Orszag & Tang 1979) were also very similar to those from the Athena code, with slightly more diffusion for smooth waves because of the lower order of accuracy.

The advection of a magnetic field loop (e.g. Stone et al. 2008) across a periodic domain is a challenging test for grid-based codes. When run with a non-zero velocity in the third dimension (*v*_{z}≠ 0) this test shows a weakness of the mixed-GLM divergence cleaning method – small divergence errors lead to small force errors and some of the initial magnetic field leaks unphysically into the direction (∼5 per cent). Results are shown in the first three panels of Fig. A3 for the evolution of the magnetic pressure. The left-hand panel shows the initial conditions; the centre panel shows the state at *t* = 2 for a model with no advection (to show the numerical diffusion) and the right-hand panel shows the state at *t* = 2 for a model advected twice across the domain. The time decay of magnetic pressure for the advected models is shown in the right-most panel of Fig. A3. The constrained transport method employed in the Athena code, together with its third-order accuracy, leads to less diffusion and a more symmetric solution to this problem than is obtained with our code (cf. Stone et al. 2008), but our results are comparable to those of other second-order astrophysics MHD codes. Tests also showed that for this test the van Albada slope-limiter (van Albada et al. 1982) produced a noticeably superior solution than the commonly used MinMod limiter, probably because the latter is significantly more diffusive.

The expansion of a multidimensional adiabatic blast wave provides a good test of a code's robustness in modelling 3D shocks and rarefactions. We ran the same 3D problem as described in Stone & Gardiner (2009, section 6.5), with the gas pressure in a sphere at the centre of the domain set to 100× the ambient pressure, and with a uniform magnetic field at 45° to the grid axes. Very similar results were found for the case where the initial plasma beta parameter β = 0.2. Because of the robustness of the Dedner et al. (2002) algorithm, and presumably also because the integration algorithm used is slightly more diffusive than the Athena algorithm, our code can also simulate this test problem with a field 10× stronger (β = 0.002), shown in Fig. A4. This extra robustness is important for the strong-field simulations run here.

### APPENDIX B: BOUNDARY CONDITION FOR MIXED-GLM DIVERGENCE CLEANING ALGORITHM

Dedner et al. (2002) suggest that for zero-gradient boundary conditions (BCs), using zero gradient for the unphysical scalar field ψ is sufficient. We have found, however, that more care is needed for magnetically dominated regions (plasma parameter β≡ 8π*p*_{g}/*B*^{2}≪ 1) when gas is flowing subsonically near the boundary. In what follows is the normal vector to the boundary, the computational domain is ‘left’ of the boundary and off the domain is ‘right’. Note that ψ is not a physical variable so its gradient is not constrained. In fact if the field is constant across the boundary then the formal solution to the 1D Riemann problem is that the flux, *F*[*B _{x}*] = 0 (as it is for all 1D Riemann problems). So ideally a boundary condition for ψ should be selected which ensures

*F*[

*B*] = 0. From Dedner et al. (2002, equations 41 and 42) it is easy to calculate this: the zero-gradient condition states that

_{x}*B*−

^{r}_{x}*B*= 0 so if additionally ψ

^{l}_{x}^{r}=−ψ

^{l}then

*F*[

*B*] = 0 is obtained.

_{x}A sequence of figures from a 2D slab-symmetric test simulation (this time using cgs units) are shown in Fig. B1 to demonstrate the effectiveness of the new boundary condition. The simulation domains used have sizes [0.2, 0.1] pc (large) and [0.1, 0.1] pc (small) and the radiation source is at infinity in the direction. A group of randomly located dense clumps are placed between *x* = 0.02 and 0.08 pc on a uniform background density of *n*_{H} = 100 cm^{−3}. The initial magnetic field is G, and the initially constant gas pressure is set to give β = 0.017 (*T* = 1500 K in the lowest-density gas). Zero-gradient BCs are enforced in the directions, and periodic in .

Fig. B1 shows the gas pressure at time *t* = 12 kyr using cooling model C2. The problem with the zero-gradient BC for ψ is very apparent, as is the dramatic improvement on switching to the ψ^{r}=−ψ^{l} BC. The black regions are hiding cells which obtained negative pressures and were set to an artificial pressure floor value in order to allow the simulation to continue. Regarding the other boundaries, the *x* = 0.2 pc boundary only has waves moving perpendicular to it and so there is negligible net flow across the boundary. The *x* = 0 pc boundary has strong supersonic outflows, so the details of the BC for ψ have little effect on it.

### APPENDIX C: BOUNDARY EFFECTS AS A FUNCTION OF DOMAIN SIZE

To study boundary effects on the results presented here, test simulations using clump configuration 1 were performed with increasingly larger domains, both in 2D and 3D with a medium perpendicular magnetic field (model R5 in Table 2). The same physical resolution was used for each simulation; boundary positions for five 3D simulations are listed in Table C1. Various global quantities were measured within the volume common to all simulations. Neutral mass, mean neutral gas velocity and position of the first neutral cell on the domain are plotted as a function of time in Fig. C1. The small and medium simulations are clearly significantly affected by the boundaries whereas the large model is essentially identical to the X-long and X-wide models, indicating that the results have converged at least in this limited subdomain of the simulation. The position of the leading edge of the ionization front is unaffected by the position of the boundaries, but the neutral gas mass and velocity are significantly affected. It was found that flow of gas through the distant boundary was impeded in these R-MHD simulations, something which did not happen in R-HD simulations.

Name | x_{min} | x_{max} | y_{min} /z_{min} | y_{max} /z_{max} | cells |

Small | 2.0 | 4.0 | −0.75 | 0.75 | 128 × 96^{2} |

Medium | 1.5 | 4.5 | −1.125 | 1.125 | 192 × 144^{2} |

Large | 1.5 | 5.5 | −1.5 | 1.5 | 256 × 192^{2} |

X-Long | 1.5 | 6.5 | −1.5 | 1.5 | 320 × 192^{2} |

X-Wide | 1.5 | 5.5 | −2.0 | 2.0 | 256^{3} |

Name | x_{min} | x_{max} | y_{min} /z_{min} | y_{max} /z_{max} | cells |

Small | 2.0 | 4.0 | −0.75 | 0.75 | 128 × 96^{2} |

Medium | 1.5 | 4.5 | −1.125 | 1.125 | 192 × 144^{2} |

Large | 1.5 | 5.5 | −1.5 | 1.5 | 256 × 192^{2} |

X-Long | 1.5 | 6.5 | −1.5 | 1.5 | 320 × 192^{2} |

X-Wide | 1.5 | 5.5 | −2.0 | 2.0 | 256^{3} |