Radiation-MHD simulations of H II regions and PDRs

# Radiation-magnetohydrodynamic simulations of H  Ii regions and their associated PDRs in turbulent molecular clouds

S. J. Arthur, W. J. Henney, G. Mellema, F. De Colle and E. Vázquez-Semadeni
Centro de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Campus Morelia, Apdo. Postal 3-72, 58090 Morelia, Michoacán, México
Department of Astronomy & Oskar Klein Centre, AlbaNova, Stockholm University, SE 10691 Stockholm, Sweden
Department of Astronomy and Astrophysics, University of California Santa Cruz, Santa Cruz, CA 95064, USA
###### Abstract

We present the results of radiation-magnetohydrodynamic simulations of the formation and expansion of H  II regions and their surrounding photodissociation regions (PDRs) in turbulent, magnetised, molecular clouds on scales of up to 4 parsecs. We include the effects of ionising and non-ionising ultraviolet radiation and x rays from population synthesis models of young star clusters. For all our simulations we find that the H  II region expansion reduces the disordered component of the magnetic field, imposing a large-scale order on the field around its border, with the field in the neutral gas tending to lie along the ionisation front, while the field in the ionised gas tends to be perpendicular to the front. The highest pressure compressed neutral and molecular gas is driven towards approximate equipartition between thermal, magnetic, and turbulent energy densities, whereas lower pressure neutral/molecular gas bifurcates into, on the one hand, quiescent, magnetically dominated regions, and, on the other hand, turbulent, demagnetised regions. The ionised gas shows approximate equipartition between thermal and turbulent energy densities, but with magnetic energy densities that are 1 to 3 orders of magnitude lower. A high velocity dispersion ( km s) is maintained in the ionised gas throughout our simulations, despite the mean expansion velocity being significantly lower. The magnetic field does not significantly brake the large-scale H  II region expansion on the length and timescales accessible to our simulations, but it does tend to suppress the smallest-scale fragmentation and radiation-driven implosion of neutral/molecular gas that forms globules and pillars at the edge of the H  II region. However, the relative luminosity of ionising and non-ionising radiation has a much larger influence than the presence or absence of the magnetic field. When the star cluster radiation field is relatively soft (as in the case of a lower mass cluster, containing an earliest spectral type of B0.5), then fragmentation is less vigorous and a thick, relatively smooth PDR forms. Accompanying movies are available at http://youtube.com/user/divBequals0

###### keywords:
H II regions – ISM: kinematics and dynamics – magnetohydrodynamics – photon-dominated region (PDR) – Stars: formation
\setkeys

Ginwidth=

## 1 Introduction

H  II regions, that is, regions of photoionised gas, are among the most arresting astronomical objects at optical wavelengths. The basic theory behind their formation and expansion has been known for some time (Strömgren, 1939; Kahn, 1954; Spitzer, 1978). More recently, attention has turned to explaining the irregular structures, filaments, globules and clumps seen within and around these regions. These could be due to underlying density inhomogeneities in the molecular cloud into which the H  II region is expanding, or could be due to instabilities at the ionisation front itself (Garcia-Segura & Franco, 1996; Mellema et al., 2006a; Williams et al., 2001; Gritschneder et al., 2010; Mackey & Lim, 2010a). The effect of the stellar ionising radiation on these structures can lead to radiation driven implosion, and the subsequent compression could result in the formation of new stars (Bertoldi, 1989; Esquivel & Raga, 2007; Motoyama et al., 2007).

Magnetic fields are well known to pervade molecular clouds and are thought to play an important role in regulating the star-formation process by providing support to collapsing structures in molecular clouds (Mouschovias & Ciolek, 1999), although recently their relative importance in this process has become less clear (McKee & Ostriker, 2007). However, theoretical work on the interplay between ionising radiation and magnetic fields in and around H  II regions has only been possible recently, with the advent of radiation-magnetohydrodynamic codes. Krumholz et al. (2007) performed the first simulations of the expansion of an H  II region in a uniform, magnetised medium. At early times, the thermal pressure of the photoionised gas dominates over the magnetic field, which is pushed out of the expanding H  II region. At late times (several Myr), the photoionised region becomes elongated along the field lines and the magnetic field filled back in. Henney et al. (2009) and Mackey & Lim (2010b) have studied the effects of ionising radiation on magnetised globules, finding that the photoevaporation process is significantly altered only when the magnetic pressure exceeds the thermal pressure by more than a factor of ten in the initial globule. Peters et al. (2010) have recently studied the combined effects of magnetic fields and photoionisation on the formation of a massive star on scales  pc, finding that rotational winding of the field lines produces a magnetised ‘bubble’, whose magnetic pressure acts to confine the nascent H  II region during its ultracompact phase.

In this paper, we study the expansion of H  II regions in non-uniform (turbulent) magnetised molecular clouds. Previously, we studied H  II region expansion in a turbulent medium without a magnetic field and our results were strikingly similar to observed H  II regions (Mellema et al., 2006a). In the present work, we study the effect of the expanding photoionised gas on structures in the surrounding magnetised medium and also the effect that the turbulent magnetic field has on the H  II region itself. At least at the relatively early times studied in our simulations, most of the interesting effects due to the magnetic field will occur within the neutral, photon dominated region around the H  II region, and so we use a careful treatment of the heating and cooling processes in this gas, first discussed in Henney et al. (2009), in order to have a realistic portrayal of the dynamical processes here.

The wealth of new observational data at infrared wavelengths enables comparisons to be made of our simulations not only with observations of the photoionised gas, as done previously in Mellema et al. (2006a) but also with the PDR and molecular gas affected by the H  II region. In particular, the Spitzer GLIMPSE surveys of bubbles in the Galactic plane (Churchwell et al., 2006, 2007) highlight polycyclic aromatic hydrocarbon (PAH) emission at the outer edges of H  II regions and in PDRs. These surveys show that the thickness of the PAH-emitting shell varies between 0.2 and 0.4 of the outer radius and that shell thickness increases approximately linearly with radius. Moreover, the bubbles are generally non-spherical. Herschel and APEX observations reveal thermal emission from cold dust in the molecular gas outside the PDR (e.g. Anderson et al., 2010; Deharveng et al., 2010). These observations are useful in identifying condensations that could be collapsing to form new stars.

Studies of magnetic fields in and around molecular clouds are important for understanding the star-formation process. There is no consensus as to the importance of magnetic fields in the formation of clouds, cores, and ultimately stars. If magnetic fields are dynamically important, then star formation is regulated by ambipolar diffusion (e.g., Mouschovias & Ciolek, 1999). Alternatively, if magnetic fields are unimportant then other processes such as turbulence and stellar feedback provide support to molecular clouds and determine the star-formation efficiency (e.g., Ballesteros-Paredes et al., 1999; Elmegreen, 2000; Matzner, 2002; Krumholz, Matzner, & McKee, 2006; Tasker & Tan, 2009; Tasker, 2011). It is not known whether magnetic fields are strong enough to dynamically support molecular clouds and more evidence is needed. It is possible that reality lies somewhere between these extremes. Recent Zeeman observations of molecular clouds (Crutcher, Hakobian, & Troland, 2009; Crutcher et al., 2010) do not agree with predictions of idealized ambipolar diffusion models (strong magnetic field) but more observations are needed.

There are three main techniques for tracing the magnetic field in diffuse HI and molecular clouds: polarisation from aligned dust grains, linear polarisation of spectral lines, and Zeeman splitting of spectral lines (Heiles & Crutcher, 2005). Polarisation studies yield , the strength of the field projected onto the plane of the sky, whereas Zeeman splitting gives , the strength of the field parallel to the line of sight. In particular, the H  I Zeeman effect has been used to study the atomic gas (assumed to be in the photodissociated region) in the star-forming regions M17 (Brogan & Troland, 2001) and the Ophiuchus dark cloud complex (Goodman & Heiles, 1994), while the OH Zeeman effect is used to study the line-of-sight magnetic field in the molecular gas. The results of such studies suggest that there is rough equipartition of magnetic and turbulent (dynamical) energy densities in the neutral gas. In the Ophiuchus dark cloud region, the uniform field in the H  I gas is estimated to be 10.2 G (Goodman & Heiles, 1994), whereas the magnitude of the line-of-sight magnetic field strength in the H  I gas associated with M17 is found to be G (at resolution, or even twice as high at resolution) (Brogan & Troland, 2001), suggesting that there is small-scale structure in the magnetic field. Recent linear polarisation studies at submillimetre wavelengths of the magnetic field around the ultracompact H  II region G5.89–0.39 (Tang et al., 2009) provide evidence for compression of the B field in the surrounding molecular cloud and general disturbance of the magnetic field caused by the expansion of the H  II region.

For photoionised gas, Faraday rotation measures of line-of-sight extragalactic radio sources are used to calculate the magnetic field strengths in foreground H  II regions (Heiles et al., 1981). The Zeeman effect cannot be used in the ionised gas because the thermal linewidths are too large. Heiles et al. (1981) use double extragalactic radio sources to sample two positions in the H  II regions S117 and S264 and find identical rotation measures in each case of order 1 or G. In these regions, the thermal energy density dominates the magnetic energy density. In a third H  II region, S119, for which only a single rotation measure is available, the line-of-sight magnetic field strength is around 20 G and the ratio of magnetic to thermal energy density is around 0.4, hence the magnetic field could have had an effect on the expansion of this H  II region.

The remainder of this paper is organised as follows. In § 2 we describe our numerical algorithm and present two test problems: the classical Spitzer law for the expansion of a non-magnetised H  II region in a uniform medium, and the expansion of an H  II region in a uniform magnetised medium (Krumholz et al., 2007; Mackey & Lim, 2010b). In the second case, we draw attention to some interesting properties of the solution that have not been described previously.

In § 3 we describe the expansion of H  II regions in both magnetised and non-magnetised turbulent media, taking as initial conditions for the ambient medium the results of an MHD turbulence calculation (Vázquez-Semadeni et al., 2005), and considering both a strong (approximately O9 spectral type) ionising source and a weak (B0.5) one. In § 4 we make qualitative comparisons between simulated optical and long-wavelength images and observations. We also produce synthetic maps of the projected line-of-sight and plane-of-sky components of the magnetic field, which can be qualitatively compared with observations of the magnetic field in star-forming regions. We discuss the globules and filaments at the periphery of our simulated H  II regions in the context of recent numerical work on photoionised magnetised globules (Henney et al., 2009; Mackey & Lim, 2010b). Finally, in § 5 we outline our conclusions.

## 2 Numerical Model

### 2.1 Numerical Method

The numerical method used in this paper is the same as that described by Henney et al. (2009). The Phab-C code combines the ideal magnetohydrodynamic code described by De Colle & Raga (2006) with the C-Ray method (Conservative Causal Ray) for radiative transfer developed by Mellema et al. (2006b), and adds a new treatment of the heating and cooling in the neutral gas (Henney et al., 2009). The MHD code is a second-order upwind scheme that integrates the ideal MHD equations using a Godunov method with a Riemann solver, similar to that described by Falle et al. (1998) but with the constrained transport (CT) method (e.g., Tóth, 2000) used to maintain the divergence of the magnetic field close to zero. An artificial (Lapidus) viscosity is also included (Colella & Woodward, 1984) in order to broaden the discontinuities and reduce oscillations. The C-Ray method for radiative transfer is explicitly photon conserving (Mellema et al., 2006b) and the algorithm allows for timesteps much longer than the characteristic ionisation timescales or the timescale for the ionisation front to cross a numerical grid cell, by means of an analytical relaxation solution for the ionisation rate equations. This results in a highly efficient radiative transfer method.

The MHD and radiative transfer codes are coupled via operator splitting, the inclusion of equations for the advection of ionised and neutral density in the MHD code, and the energy equation, which includes both radiative cooling and heating due to the absorption of stellar radiation by gas or dust: ionising extreme ultraviolet radiation in the photoionised region, non-ionising far ultraviolet radiation in the neutral gas and x-rays in the dense molecular gas. A detailed discussion of the heating and cooling contributions is given by Henney et al. (2009). In particular, the treatment of heating in the neutral PDR region, taking into account FUV/optical dust heating and x-ray heating, is a substantial improvement over standard treatments in the literature. The heating term for the neutral gas is calibrated using Cloudy (Ferland et al., 1998) and is tailored to the ionisation source. Moreover, x-rays and FUV radiation will be enhanced due to the presence of associated low-mass stars in the star-formation environment, and we use the results of Fatuzzo & Adams (2008) to estimate the distributions of FUV and EUV fluxes appropriate to stellar clusters of membership sizes relevant to the two cases we study: an O9 and a B0.5 star. The energy equation is solved iteratively using substepping and the timestep for the whole calculation is determined solely by the magnetohydrodynamic timestep.

The individual components of the Phab-C code have been tested against standard problems (De Colle, 2005; De Colle & Raga, 2004; Iliev et al., 2006) and the whole code was applied to the problem of the photoionisation of magnetised globules Henney et al. (2009), work which has recently been independently corroborated by Mackey & Lim (2010b). In this paper, we have implemented an entropy fix (Balsara & Spicer, 1999) to correctly deal with the thermal and dynamical pressures in regions of low density where the magnetic pressure dominates. We have also improved the treatment of the boundary conditions because they can become important, particularly towards the end of the simulations. For our O star simulations we use outflow boundary conditions, since thermal pressure is expected to dominate throughout the evolution and the H  II region will rapidly grow to the size of the computational box. For the B star simulations, the expansion timescale of the H  II region is much longer and is comparable to the crossing time of the computational box corresponding to the initial turbulent velocity dispersion of the molecular gas. We therefore use periodic boundary conditions in this case, which ensures that the total mass in the computational box remains constant.

Parallelisation using Open-MP was implemented where possible, enabling the current simulations to be run in about two weeks on an Intel 8-core server with 32 GB of RAM. Preliminary work was performed on the Kan Balam supercomputer of the Universidad Nacional Autónoma de México.

### 2.2 Test Problems

#### 2.2.1 Hydrodynamic H  Ii Region Expansion in a Uniform Medium

Org/

The radiation part of the code has been extensively tested and documented (Mellema et al., 2006b; Iliev et al., 2006). In order to test the radiation (magneto-)hydrodynamics combined code we first consider the purely hydrodynamic evolution of a photoionised region around an ionising source in a uniform medium. The analytical solution to this problem is the familiar Spitzer (1978) expansion law

 RSpitzer=R0(1+74citR0)4/7 , (1)

where is the ionisation front radius, is the initial Strömgren radius, is the sound speed in the ionised gas and is time, is the ionising photon rate, is the ambient density and is the case B recombination rate, (Osterbrock & Ferland, 2006), where is the temperature in the ionised gas. For an ambient density  cm, stellar ionising photon rate  s and stellar temperature  K the results are shown in Fig. 1, together with the corresponding analytic solution. These parameters were chosen because they represent the mean values of the turbulent medium B star case described below in § 3. The lower panel shows radius against time, while the upper panel shows the relative error, defined as

 Error=Rion−RSpitzerRSpitzer . (2)

From this figure, we see that the largest discrepancies, of about 6%, occur at the beginning of the expansion, when a shock wave begins to propagate outwards ahead of the ionisation front and a corresponding rarefaction wave travels back towards the ionising source, which is not accounted for in the analytical solution. At later times, the error is never more than 1%. Our code, therefore, compares extremely favorably with other radiation-hydrodynamics codes on this problem (see, e.g. Arthur & Hoare, 2006; Krumholz et al., 2007; Iliev et al., 2009). There is no evidence for overcooling in the ionisation front in our simulations, as was found to be a problem by Krumholz et al. (2007). In the analytic Spitzer solution we assume an H  II region temperature  K, which gives the best fit to our numerical result. This is about  K hotter than the volume-averaged ionised gas temperature in our simulation, but is roughly equal to the ionised gas temperature just inside the ionisation front.

In Fig. 2 we plot the gas radial velocity from the numerical simulation and the corresponding analytical solution. The numerical solution closely follows the analytical solution, though with oscillations due to acoustic phenomena, such as the initial rarefaction wave, which propagates back into the H  II region when the ionization front begins to expand outwards.

#### 2.2.2 Mhd H  Ii Region Expansion

There is no corresponding analytical solution for the expansion of an H  II region into a magnetised medium. Therefore, in order to test the MHD part of the numerical code, we carry out a simulation of the expansion of an H  II region in a uniform magnetised medium using the same parameters as those of Krumholz et al. (2007). This problem was also used to test the radiation MHD code of Mackey & Lim (2010b). The simulation is performed in a computational box with constant ambient density and temperature of  cm and  K, respectively, a constant ambient magnetic field directed at to the -axis in the - plane111The original simulation performed by Krumholz et al. (2007) has the magnetic field parallel to the -axis. of strength G and an ionising source producing  ionising photons s. The magnetic field is placed at an angle to the grid lines in order to be able to distinguish between physical and computational effects. We run the simulation in a box and in a box, the former to study the early evolution  Myr) and the latter to highlight the late-time evolution. We note that the ionising source is very weak (later than spectral type B0.5, Vacca et al., 1996) and that the ambient density is lower than that typically found in star-forming regions.

Krumholz et al. (2007) found a critical radius and time for when magnetic pressure and tension start to become important for the expansion of the H  II region in a uniform magnetised medium, which corresponds to when the magnetic pressure in the neutral gas becomes of the same order as the thermal pressure in the ionised gas, i.e., . Here, , are the densities in the neutral and ionised gas, is the Alfvén speed, and is the sound speed in the ionised gas. Thus, magnetic effects become significant when the H  II region has expanded to a radius of roughly

 Rm≡(civA)4/3R0 , (3)

where is the initial Strömgren radius in a non-magnetised medium of the same uniform density. The corresponding time is

 tm≡47(civA)7/3t0 , (4)

where is the sound crossing time of the initial Strömgren sphere, and both of these equations implicitly assume that the classical Spitzer expansion law (Spitzer, 1978) holds until the magnetic field becomes important.

For the parameters in this test problem, the initial Strömgren radius is  pc, and the Alfvén speed in the ambient magnetised medium is  km s. Our simulations give a temperature in the photoionised gas of  K and hence a sound speed of  km s (compare the lower assumed ionised sound speed  km s in the Krumholz et al. 2007 models). Thus, the critical radius and time are  pc and  Myr, respectively, for our models.

In Figs. 3 and 4 we show snapshots of the magnetic field and the photoionised gas distribution in the central -, - and - planes of the computational box after 2 and 6 Myr of evolution. Â The coloured discs shown between the panels of the top row are keys to the interpretation of the magnetic field images. The hue indicates the projected field orientation in the plane of each cut (red for vertical, cyan for horizontal, etc.), while the brightness indicates the magnetic field strength, as is shown in the left-hand disc. The colour saturation diminishes as the out-of-plane component of the field increases, as shown in the right-hand disk. This non-standard representation allows the magnetic field to be sampled at every pixel. The lower panels of the figures show the ionisation fraction, temperature and density of the gas in the same central planes. Ionisation fraction is represented by colour: red indicates fully ionised gas while blue/green is partially ionised. The gas temperature varies beween  K (red) in the ionised gas to a few hundred Kelvin in the neutral (PDR) gas. The gas density is represented by the intensity in these figures, with the densest, cold molecular gas appearing white.

At early times, , the H  II region is essentially spherical and expands at approximately the same speed (half the sound speed) in all directions because the ionised gas thermal pressure is much higher than the ambient magnetic pressure. Expansion across the magnetic field lines compresses the B-field. The expanding H  II region and PDR are preceded by a fast-mode shock, which travels furthest perpendicular to the magnetic field direction. The fast-mode shock bends and compresses the field lines and so the PDR and H  II region propagate in a modified magnetic field. The slow-mode shock is strongest along the field lines and most evident at the end caps, where it produces a dense shell. It compresses the gas and demagnetises it (reduces the field). The slow-mode shock is almost isothermal and the dense shell behind it is pushed by the pressure in the photoionised region. At early times the magnetic field is reduced in the photoionised gas by the expansion of the H  II region. At later times, , the magnetic tension pulls the field back into the H  II region at the equator.

In Fig. 3 we show the situation after 2 Myr, when the magnetic field has started to pull back into the H  II region. The H  II region is elongated along the magnetic field direction but the photon dominated region (warm neutral gas) is approximately spherical. The magnetic field retains its original direction except around the edge of the H  II region. This is particularly evident at the poles of the H  II region. At the time shown in the figure, the magnetic field appears concentrated towards the centre, with regions of lower magnetic field around the equator or waist of the structure. Neutral gas is pulled into the H  II region along with the magnetic field, becomes ionised and flows away from the ionisation front into the H  II region. It then flows along the field lines and is channelled out towards the poles, where the ionisation front transforms into a recombination front. We note that at this time in the simulation, we do see instabilities form in the direction perpendicular to the magnetic field, as reported by Krumholz et al. (2007) and Mackey & Lim (2010b). These instabilities form at the equator where the ionisation front is parallel to the magnetic field direction because the slow-mode shock which precedes the ionisation front, is not permitted to cross the field lines. As a result, the ionisation front corrugates.

Our Fig. 4 shows the simulation after 6 Myr in a computational box. The H  II region has stopped expanding perpendicular to the magnetic field but the shocked neutral shell is still expanding along the field lines, leading to a very elongated structure. However, in the direction along the field lines, the flow has become one dimensional and so the density does not fall off as the H  II region expands. Thus, there is a limit as to how far the ionisation/recombination front can travel in this direction. There is an extended region of partially ionised gas beyond the recombination front. The warm, neutral gas (PDR) is still roughly spherical in shape but is pierced by the cigar-shaped H  II region along the field direction. The dense neutral knots formed as a result of the instability at the equator shadow parts of the PDR. The knots formed closest to the perpendicular direction remain there throughout the simulation, pinned by the magnetic field. The knots slightly off the equator experience the rocket effect due to the photoionisation of their outer skins. They move like beads on a wire along the magnetic field lines around the edge of the H  II region, a process which can be clearly observed in animations of this simulation. By the end of the calculation, the magnetic field has refilled the H  II region and become roughly uniform and returned completely to its original direction.

Finally, in Figs. 5 and 6, we examine the expansion of the magnetised H  II region. Fig. 5 shows the evolution of the radius of the photoionised region with time. The mean radius is consistently below the classical Spitzer expansion law, unlike the non-magnetic case (see Fig. 1). In the direction perpendicular to the magnetic field, the expansion slows with time, as shown by the evolution of the minimum radius. In the direction along the field lines, the expansion initially ( Myr) follows the classical Spitzer law but then accelerates as material is channelled along the field lines and out of the poles of the H  II region. This expansion abruptly stops just after 4 Myr, which corresponds to the time when the ionisation front can no longer expand along the field lines because the flow has become one dimensional.

In Fig. 6 the mean radial and root-mean-square (rms) velocities of the ionised gas are plotted together with the analytical value and also the mean velocity of the ionisation front. The discontinuity at 2 Myr is due to the change in resolution of the simulation. The calculated values are all above the analytical values and indicate that typical velocities in the photoionised gas are of the order 1–2 km s throughout the simulated lifetime of the H  II region. The mean radial velocity shows the same sort of ringing that we saw in the non-magnetised simulation (see Fig. 2). After about 3 Myr, the gas velocities become dominated by the one-dimensional flow being channelled along the field lines.

## 3 H  Ii Region Expansion in a Magnetised Turbulent Medium

In this section we describe the results of radiation MHD simulations of H  II region expansion in a more realistic medium where neither the initial density distribution nor magnetic field are uniform. We consider both an O star ionising flux and one more appropriate to a B star in order to study the effect of the magnetic field on the ionised gas and on the neutral gas of the PDR.

### 3.1 Initial Set Up

As initial conditions we take the results of a driven-turbulence, ideal magnetohydrodynamic simulation by Vázquez-Semadeni et al. (2005), specifically the simulation with initial isothermal Mach number , Jeans number and plasma beta . The original Vázquez-Semadeni et al. (2005) simulations are scale free and are characterised by the three nondimensional quantities (the rms sonic Mach number of the turbulent velocity dispersion ), (the Jeans number, giving the size of the computational box in units of the Jeans length ), and the plasma beta, (ratio of thermal to magnetic pressures).

We set the scaling by choosing the computational box size to be 4 pc, giving the ambient temperature as  K, the mean atomic number density as  cm and the mean magnetic field strength as G. We choose as our starting point the time in the evolution when there is one collapsing object. At this point, the mean values of the magnetic field the number density are unchanged, since flux and mass are conserved, while the mean plasma beta is now , which is consistent with the rms value of the magnetic field G. That is, the turbulence has enhanced the rms magnetic field by a factor , which leads to a change by a factor in the value of over the original, uniform conditions.

As in our previous paper (Mellema et al., 2006a), we place our ionising source in the centre of the collapsing object, remove of material (corresponding to the mass of the star formed) and take advantage of the periodic boundary conditions of the turbulence simulation to move the source to the centre of the grid. We also subtract the source velocity from the whole grid. We consider two ionising sources: one corresponds approximately to an O9 star (Vacca et al., 1996), having an effective temperature  K and an ionising photon rate of  photons s, while the other ionising source has  photons s, which corresponds to B0.5.

### 3.2 Results

In this section we present our results both for the O star and the B star and in each case for the magnetic and also for the purely hydrodynamic simulations. In this way we can assess the importance of the magnetic field and the ionising photon flux on the evolving structures in the expanding H  II regions and surrounding neutral gas. We first present emission-line images of the simulated H  II regions for both O and B stars in the magnetic and pure hydrodynamic cases when they have reached a similar size, which corresponds to 200,000 yrs of evolution for the O star case and  yrs of evolution in the B star case. This permits a direct comparison with observations at optical and longer wavelengths. We then examine the expansion and time evolution of the ionised, neutral and molecular components using global statistical properties. The importance of the magnetic field in the different components is then studied for the particular O and B star simulations presented earlier.

Although the ionised/neutral transition of hydrogen is treated explicitly in a detailed manner in our radiation-MHD code, the neutral/molecular transition is treated much more approximately. We simply assume that the molecular fraction is a predetermined function of the extinction from the central star to each point in our simulation: , which was determined by fitting to the Cloudy models described in the Appendix of Henney et al. (2009).

A rough estimate of the importance of the magnetic field in our turbulent media simulations can be obtained by calculating the critical radius and timescale, and (see Eqs. 3 and 4), for the initial value of the magnetic field and the mean density. Using G and  cm, and assuming a mean particle mass of 1.3 m, the initial representative Alfvén speed in the neutral gas is  km s. In the photoionised gas, the sound speed is  km s. For the strong ionising source (O star), the Strömgren radius in a uniform medium of density equal to the mean density would be  pc, giving a critical radius  pc and time  Myr for when the magnetic field starts to become important. If, instead of the rms magnetic field we use the mean magnetic field, G, then the corresponding values are  pc and  Myr. For the weak ionising source (B star), the numbers are  pc,  pc ( pc) and  Myr ( Myr) depending on whether the rms or mean magnetic field is considered.

Thus, for the strong ionising source, the critical diameter is larger than the dimension of the computational box, hence we do not expect that the magnetic field would have much effect on the expansion of the H  II region in this case. However, for the weak ionising source, it is not clear whether the rms or the mean value of the magnetic field is most important for determining the global evolution of the H  II region. If the rms field is the determining parameter then the fact that the critical diameter in this case is smaller than the size of the computational box leads us to expect that the magnetic field will affect the global expansion of the H  II region. On the other hand, if it turns out that the mean magnetic field is most important, then the effect will not be noticed within our computational size and timescales. Furthermore, given the non-uniform initial conditions, the timescale for the magnetic field to become important could vary depending on the local gas sound speed to Alfvén speed ratio.

#### 3.2.1 Morphologies and Images

We start our comparison by considering the morphological appearance of the H  II regions and the surrounding neutral and molecular material. In Figs. 7 and 8 we show two types of image data, one showing mostly the ionised gas, where the different colours represent emission from optical emission lines, the other shows synthetic, long-wavelength (infrared and radio) emission, emphasizing the neutral and molecular material.

The optical emission was calculated as described in our previous papers (Mellema et al., 2006a; Henney et al., 2009), assuming heavy element ionisation fractions that are fixed functions of the hydrogen ionisation fractions. For the B star models, these functions were recalibrated for the softer stellar ionising spectrum using the Cloudy plasma code (Ferland et al., 1998). For the O star models, we employ the “classical” HST red-green-blue colour scheme of [N  II ], H, [O  III ]. Since the [O  III ] emission from the B star models is predicted to be very weak, in this case we use the scheme [S  II ], [N  II ], H. In both schemes, the progression from red through green to blue corresponds to increasing degree of ionisation inside the H  II region. The emission from all the optical lines is negligible in the dense neutral/molecular zones, but the dust absorption associated with these dense regions is visible in the images. Scattering by dust is not included.

For the long-wavelength images the emission bands were chosen to give a more global view of the simulations than can be provided by the optical emission lines. The red band shows the total column density of neutral/molecular gas, which very crudely approximates far-infrared or sub-mm continuum emission from cool dust. The green band of the long-wavelength images shows a simple approximation to the mid-infrared emission from polycyclic aromatic hydrocarbons (PAH). The PAHs are assumed to reprocess a fixed fraction of the local far-ultraviolet radiation field. Detailed calculations (Bakes, Tielens, & Bauschlicher, 2001) show that this assumption is correct to within a factor of about two for most PAH emission features over a broad range of conditions. The local PAH density is assumed to be proportional to the neutral plus molecular hydrogen density, since there is strong evidence that PAHs are destroyed in ionised gas (Bregman et al., 1989; Lebouteiller et al., 2007). No attempt is made to discriminate between different PAH charge states. The blue band of the long-wavelength images shows the 6 cm radio continuum emission due to bremsstrahlung in the ionised gas. Apart from a slightly different temperature dependence this is very similar to the H emission, except that it does not suffer any dust absorption.

Starting with the O-star case, we see that the morphological appearance shows the typical characteristics of H  II regions in turbulent media found in earlier work (Mellema et al., 2006a; Dale et al., 2007; Mac Low et al., 2007), namely a fairly irregular structure with fingers or pillars, as well as bar-like features at the edge of the H  II region. Note how the most prominent finger in the bottom edge of the H  II region does not exactly point towards the source of ionising photons. In appearance the simulated H  II regions look strikingly like observed H  II regions.

Comparing the magnetic with the non-magnetic case, one sees an overall agreement in the structures; the presence of a magnetic field does not have a large impact on the global morphology of the H  II region. This is to be expected as for the O-star case the critical time (see Eq. 4) has not been reached by the time most of the volume has become ionised. However, the magnetic field does have an effect on small scale features in the H  II region. In the magnetic case the fingers and bars look smoother and broader, due to the magnetic pressure, which modifies the radiation-driven implosion of these structures (Henney et al., 2009; Ryutov et al., 2005). This behaviour is visible in both optical and long-wavelength images, but most clearly in the latter. In section 4 we discuss the behaviour of the globules in greater depth.

For the B-star case, the critical time does occur within the simulation. However, as is apparent from Fig. 8 this does not lead to a significant deformation of the H  II region. The reason is that magnetic field itself also has a turbulent structure, and there is no large-scale field which can impress its direction on the shape of the H  II region, as in the test case presented above (see § 2.2.2). Although the magnetic and non-magnetic results do differ, the overall impression is that they are still quite similar. In fact, the difference between the O-star and the B-star cases is much larger than between the magnetic and non-magnetic cases. For the B-star case, no pillars are found and the ionised region appears more spherical. The edge of the H  II region is often more fuzzy, although some sharper bar-like features are present. The long wavelength images show the reason for this: the H  II region is embedded in a thick PDR region (with an extent of of the radius of the H  II region) which erases much of the small-scale density fluctuations present in the original cold molecular medium. Close inspection reveals that the high-density regions already inside the PDR region photo-evaporate and thus lose their large density contrast. In the O-star case these are the features which develop into pillars. The fact that the temperature in the PDR is closer to 100 K than to 10 K in combination with the low ionising flux of the B-star should generally be less conducive to the formation of pillars, as shown by Gritschneder et al. (2010).

Careful examination of the images shows that otherwise the effect of the magnetic field is quite similar to that in the O-star case: small-scale structures are smoothed out in the magnetised case. However, the effect is quite marginal and does not give a clear observable diagnostic for the presence of magnetic fields or not.

#### 3.2.2 Global Properties of the H  Ii Region Evolution

The initial mean density in the computational box is  cm but the material is distributed inhomogeneously with the densest clumps having densities  cm. In Figs. 9 and 10 we can see how the mean densities in the ionised, neutral and molecular components evolve with time and the growth of the H  II regions around the O and B stars, respectively. In these figures, lines with symbols are for the MHD simulations and lines without symbols are the purely hydrodynamic case. We can see immediately that the magnetic field has a negligible effect on the densities of the different components in both the O and B star cases. As the evolution of the H  II region proceeds, the mean density of the molecular gas in the O star simulation increases slowly with time because the lower density molecular gas becomes ionised or incorporated into the neutral PDR and also because the dense clumps and filaments at the edge of the H  II region are compressed due to radiation-driven implosions, thereby increasing their density. By the end of the simulation, only these densest clumps remain, since the majority of the computational box has been ionised. In the B star case, the mean density of the molecular gas remains roughly constant throughout the simulation. This is because the H  II region does not break out into regions of lower density and remains confined within the molecular gas for the duration of the simulation. Also, there is less fragmentation in the B star simulation because density inhomogeneities in the neutral region of the PDR are smoothed out by photoevaporation flows due to FUV radiation before the ionisation front reaches them (see Fig. 8).

By the end of the respective simulations, the mean density in the ionised gas around the O star is  cm, while that around the B star is  cm even though the spatial extent is much smaller. This is because the latter is a much weaker ionising source. In both cases, the density of the neutral PDR gas tends to a roughly constant value of  cm.

In Figs. 11 and 12 we compare how the ionised, neutral and molecular gas component fractions vary throughout the evolution of the H  II regions. We calculate both volume and mass fractions for each component. The ionised gas in both the O and B star cases is thermally dominated and there is no discernible distinction between the MHD and purely hydrodynamic results. In the O star case, the ionised fraction comes to occupy 90% of the volume but only 10% of the mass of the computational box by the end of the simulation. In the B star case, the H  II region does not manage to globally break out of the computational box and so even at the end of the simulation the volume occupied by ionised gas is less than 30%. The mass fraction is in this case, since the density in the ionised gas is low. The neutral gas is most affected by the magnetic field for both O and B star simulations. Figs. 11 and 12 show that the MHD simulation has both higher mass and higher volume fractions of neutral gas than the purely hydrodynamic results. The neutral component is more important in the B star case than in the O star case, reaching of the mass and volume fractions after 500,000 yrs and remaining roughly constant thereafter. In the O star case, however, the mass fraction remains roughly constant for the second half of the simulation while the volume fraction decreases sharply. This difference in behaviour can be attributed to the greater fragmentation in the O star case (both MHD and HD) compared to the B star. The neutral material in the B star simulation is distributed in a quite broad, smooth, almost uniform density region around the H  II region, whereas in the O star simulation the neutral region around the photoionised gas is relatively thin but there is also neutral material inside the dense clumps and filaments that are formed by radiation-driven implosion, which have high density but low volume.

The expansion of an H  II region in a clumpy medium is not so easy to characterise. In Figs. 13 and 14 we plot a variety of different measures of the radius of the photoionised region as functions of time, together with the analytical solution obtained using Equation 1 for a medium of uniform density  cm. We calculate the mean radius222We define the mean radius as the radius of a sphere with the same volume as the real H  II region. and also the minimum and maximum radii of the ionisation front. For the H  II region around the O star, the mean radius closely follows the analytical solution for a long period of the simulation, possibly because the filling factor of the dense clumps and filaments which retard the ionisation front is small. The maximum radius for the ionisation front for the O star simulation quickly becomes equal to the distance to the edge of the computational box and therefore has no physical meaning after this point. Until it leaves the box, the maximum radius represents the direction in which the ionisation front is able to expand most quickly, i.e. a path of lowest density radially outward from the star.

For the B star simulation, the mean radius is always below the analytical solution. Since this is true for both the MHD and purely hydrodynamic simulations, it cannot be due to the magnetic field. The reason is that for the B star, the initial clump in which the star is embedded has a greater effect on the subsequent evolution of the H  II region than in the O star case. In fact, we see that the maximum radius of the ionisation front in the B star H  II region does follow the analytical solution. The maximum radius follows the expansion of ionisation front in the direction where it was first able to break out of the dense clump in which the star was embedded. We have already seen that the neutral medium around the H  II region in the B star case takes on a smooth and homogeneous appearance and that the mean molecular and neutral gas densities are roughly constant in time. Once the ionisation front has managed to break out of the clump, therefore, it follows the expansion law for a uniform medium in this mean density. In the majority of directions, however, the ionisation front takes much longer to break out of the initial clump and so the mean ionisation front radius is retarded compared to the analytical solution.

In both O and B star simulations, the minimum radius of the ionisation front indicates the presence of fingers and clumps of neutral material within the H  II region. For the O star simulation we have already mentioned how the magnetic field suppresses fragmentation, and this is reflected in the graph of minimum radius with time: fingers and globules formed in the purely hydrodynamic simulation are denser and survive longer closer to the star compared to the MHD simulation. In the B star case there is much less fragmentation and we see no fingers and clumps of neutral gas within the H  II region and the minimum radii of both MHD and hydrodynamic simulations seen in Fig. 14 bear this out.

In Figs. 15 and 16 we show the mean radial velocities for the three gas components. Again, we compare the results of both MHD (lines with symbols) and purely hydrodynamic simulations. Beginning with the ionised gas component, we see that for the O star the rms velocities reach a peak of about 8 km s, and then decreases very slowly, which is what we saw in our previous paper (Mellema et al., 2006a). The mean radial velocity, on the other hand, peaks at 8–10 km s then declines quite quickly. The rms velocities are due to interacting photoevaporated flows within the H  II region, which lead to gas flowing inwards as well as outwards (Henney, 2003). The mean radial velocity traces the general expansion of the photoionised gas. We see that for the non-uniform medium, the initial expansive motions are more rapid than for the simple analytical model but this is just because the gas first begins to expand into regions of lowest density (Shu et al., 2002). There is a small difference between the rms velocities in the MHD and hydrodynamic cases, with the latter being marginally higher. This could be due to the greater degree of fragmentation in the hydrodynamic simulations which, in turn, leads to stronger photoevaporated flows.

In the B star simulation the rms velocities are similar to the O star case just discussed, which suggests that photoevaporated flows are important in this H  II region, too. The pure hydrodynamic simulation shows slightly higher velocities than the MHD case, which indicates that the magnetic field plays some role in reducing the importance of photoevaporated flows within the ionised gas. The mean radial velocities show the same sort of ringing seen in the test problem of H  II region expansion in a uniform density medium (see Fig. 2). The period of the oscillations is of order  yrs, which is consistent with the sound crossing time of the photoionised region. Note that in both O and B star simulations, the actual mean radial velocities are considerably higher than the corresponding analytical values for expansion in a uniform medium with density equal to the mean density of the initial computational box ( cm).

The neutral gas velocities show more separation between hydrodynamic and MHD results, with the hydrodynamic velocities being consistently higher in both mean and rms velocities. Interestingly, the different simulation velocities appear to depart from each other at a specific point in time, which is different depending on whether the mean or rms velocity is considered. This is true for both O and B star simulations. In the O star case, after the initial acceleration, the neutral gas velocities reach 4–5 km s, while for the B star the velocities are 1–2 km s lower and show evidence for ringing as for the ionised gas.

The initial molecular gas velocities reflect the Mach 10 supersonic turbulence of the underlying density distribution. At early times, there is evidence of residual infalling motions, i.e. negative radial velocities, again from the initial conditions. These persist longer in the B star simulations than in the O star simulations. This is because the O star rapidly ionises out to the edge of the computational box, and the molecular material which remains is being accelerated radially away from the central star in clumps that are experiencing the rocket effect. In the B star case, most of the molecular material surrounding the H  II region remains undisturbed, and so the imprint of the initial conditions remains in this gas until the end of the simulation.

In summary, we find slight differences in the expansion of H  II regions between MHD and purely hydrodynamic simulations but the greatest differences are due to the different ionising and FUV fluxes of the central O or B star. Most of the differences between the MHD and hydrodynamic cases are due to the lesser amount of fragmentation in the former, indicating that magnetic fields provide some support to the neutral gas against radiation-driven implosions.

### 3.3 Magnetic Quantities

Even though we find that the presence of a magnetic field only has small effects on the growth of the H  II region, this does not mean that the field is uninteresting or unimportant. For one thing, the field is important in the still neutral medium, where it constitutes a significant fraction of the total energy.

In Fig. 17 we show the magnetic field and the distribution of the photoionised gas for the B star case in the same format as Figs. 3 and 4. We can see that the magnetic field looks very different. The initial conditions were the result of an MHD turbulence calculation and in the molecular gas (coloured grey/white in the lower panels of Fig. 17) we can see that the field appears to be randomly oriented.

In the warm, neutral (PDR) gas (coloured purple in the lower panels of the figure) there is strong evidence that the magnetic field is being aligned parallel to the ionisation front, particularly in the lower left region of the - and - plots. At some points in this image, the ionisation front has left the grid, and the periodic boundary condition in this simulation means that ionised material exits one side of the grid and reappears at the opposite boundary, where it will usually recombine because of the high column density through the molecular and neutral gas to the star at this position. At this point in the simulation, this is not affecting the global evolution of the H  II region and surrounding warm, neutral gas.

The photoionised region has a very low magnetic field intensity, showing that the field has been pushed aside by the expansion of the H  II region. Photoevaporated flows from density inhomogeneities in both the neutral and molecular gas drag the magnetic field with them into the photoionised region in long filamentary structures. However, there is little mass associated with these flows. In the - and - central planes there is a large region of partially ionised gas. This is also in the direction where the ionisation front has left the grid and could be due to a density or velocity gradient, which results in the ionisation front being unable to keep up with the gas, since it is no longer preceded by a neutral shock in this part of the computational grid.

In order to analyze the role of the magnetic field in our three components, ionised, neutral and molecular, we present colour plots showing the joint distribution of various quantities in these three components. Since the B-star and O-star results are quite comparable, we only show the B-star results.

We start by considering joint distribution of the magnetic field strength against the density. Fig. 18 shows this for the case of the B-star at  years. The distribution of the ionised material is shown in blue and displays a wide range of field strengths (0.1 to 10 G) across a narrow range of densities (10 to 20 cm; in the case of the O-star these values are some 5 times higher). The diagonal lines show the Alfvén speeds. If the sound speed in the gas is smaller than the flow will be magnetically dominated. As can be seen, the ionised material does not reach Alfvén speeds of 10 km s, the typical sound speed in the ionised medium. Consequently the ionised flow is never magnetically dominated, consistent with the results found in the previous sections.

Another way to represent the role of the magnetic field is to plot the joint distribution of magnetic and ram (turbulent) pressure, as is done in Fig. 19. In full equipartion, the turbulent, magnetic and thermal pressures should all be equal. Since the H  II region is expanding through its thermal pressure, it is the thermal pressure which is the overall dominating one. Fig. 19 shows that of the other two pressure components, turbulent pressure mostly dominates over magnetic pressure. In only a few places in the H  II region the magnetic pressure dominates.

In the neutral and molecular material (represented by green and red in the joint distribution figures) the situation is different. The joint distribution of density and magnetic field strength, Fig. 18, shows a wide range of densities ( cm) and a narrower range of magnetic field strengths (mostly around 10–30 G, although some areas, principally neutral ones, have magnetic fields as weak as 1 G). The molecular material generally has higher densities than the neutral material. Since the sound speed in the neutral and molecular regions is 1 km s or less, they are mostly magnetically rather than thermally dominated. In the distribution of pressures, Fig. 19, one sees that the neutral and molecular material clusters around equipartition of the magnetic and turbulent pressures. The molecular material has regions where the turbulent pressure dominates, and others where the magnetic pressure dominates. The neutral (PDR) material has generally somewhat lower pressures and is closer to equipartition.

Fig. 18 can be compared to observationally derived values for the magnetic field and densities in H  II regions and their surroundings. Results of Harvey–Smith et al. (in preparation) on large Galactic H  II regions show ranges in magnetic field strengths and typical densities for the ionised, neutral and molecular components very comparable to what we find in our simulation results.

## 4 Discussion

### 4.1 Comparison with observations: RCW 120

The H  II region RCW 120 and its surrounding medium have recently been extensively studied at near- to far-infrared wavelengths in the context of triggered star formation (Zavagno et al., 2007; Deharveng et al., 2009; Anderson et al., 2010; Martins et al., 2010; Zavagno et al., 2010). This small (3.5 pc diameter at 1.35 kpc distance, Zavagno et al., 2007) H  II region appears to have a single ionising source, identified from VLT-SINFONI near-IR spectroscopy and spectral line fitting to stellar atmosphere models as an O6–O8V/III star (Martins et al., 2010) with an effective temperature of  kK and ionising photon rate of . These observationally derived parameters are remarkably similar to those we have adopted for our O star simulations and so an opportunity for a direct comparison of RCW 120 with our results presents itself.

Morphologically, the observations and simulation look very similar. In particular, our Fig. 7, corresponding to the O star simulation after 200,000 yrs of evolution, is the same size as the observed bubble RCW 120. The upper panels of the figure show the optical emission (with extinction from dust), while the lower panels show synthetic 6 cm radio free-free emission from the ionised gas (blue), generic PAH emission from the PDR (green) and cold molecular gas column density (red). Figs. 1 and 2 of Deharveng et al. (2009) show the H emission from the ionised gas, the m emission from PAHs in the PDR, m emission from warm dust in the photoionised region, and m cold dust emission from the neutral and molecular gas. If the reader mentally rotates the simulated images so that the lower left corner moves to the top centre, then a certain correspondence between simulations and observations can be imagined. The ionised gas emission fills the central cavity. In the optical, the images appear criss-crossed by dust lanes. Some of these are foreground, others are the result of projection of filaments and ridges deeper inside the H  II region. The PAH emission comes from a thin shell around the periphery of the H  II region, with projection seeming to put some PAH emission in the interior of the H  II region, as can be seen in both the simulation and the observations. The simulated PAH emission presents the same sort of irregularities, filaments and clumps as are seen in the observations.

The initial total mass in the computational box is . After 200,000 yrs of evolution, the mass distribution, as seen in Fig. 11, is approximately 5% ionised gas, 55% neutral gas, and 40% molecular gas. Some of the molecular gas is undisturbed material at the edge of the computational box but quite a large proportion must be in the clumps and filaments in the H  II region and PDR, where it becomes compressed by radiation-driven implosion. The H  II region has influenced virtually all of the computational domain by this time, and the of neutral and molecular material is comparable to the 1100– derived by Zavagno et al. (2007) and Deharveng et al. (2009) from the m cold dust emission. This dust will be present in both the neutral and molecular gas. Hence, we can surmise that the average density in the vicinity of RCW 120 is similar to that in our computational box, that is  cm, though, of course, the density distribution is far from homogeneous. From this, we could go so far as to postulate that the age of RCW 120 must be similar to the time depicted in our simulation, that is, of order 200,000 yrs. Martins et al. (2010) were unable to assign an age to RCW 120 except to say that it must be younger than 5 Myr.

Although radiation-driven implosion does enhance the density of photoionised globules at the periphery of the simulated H  II region, the fact that we do not include self gravity in our simulations means that we are unable to model triggered star formation in the neutral shell around the H  II region. However, it is quite likely that radiation-driven implosion could trigger gravitational collapse in clumps that are already on the point of forming stars.

There are other potentially important physical processes that we do not include in our models. Dust grains are an important component of the interstellar medium in star-forming regions. They absorb ionising photons and both observations (e.g., Wood & Churchwell, 1989) and theory (e.g., Arthur et al., 2004) suggest that more than 50% of all the ionizing photons are absorbed by dust within the H  II regions in the initial stages of evolution when the ionised density is high. Radiation pressure on dust grains can form a central cavity in the H  II region (Mathews, 1967; Inoue, 2002), which can be as large as 30% of the ionisation front radius. Such cavities have been known for a long time from optical observations of H  II regions (e.g., Menon, 1962) and are also seen in more recent infrared observations (e.g., Watson et al., 2008). Krumholz & Matzner (2009) have studied the effect of radiation pressure on the dynamics of H  II regions and conclude that it is only an important effect for massive star clusters and not for H  II regions around individual massive stars. Stellar winds are also to be expected from stars whose effective temperatures are greater than 25,000 K and these winds will provide an alternative mechanism for evacuating a central cavity in the dust and gas distribution in the H  II region. The mechanical luminosity of the stellar wind is converted into thermal pressure by shock waves (Dyson & Williams, 1997) and this could affect the dynamics of the H  II region.

### 4.2 Predicted maps of projected magnetic field

Figs. 20 to 23 show visualisations of the projected integrated B-field (line-of-sight and plane-of-sky components) from our simulations. These visualisations represent generic idealised versions of maps that can be obtained from various observational techniques. The line-of-sight field component can be determined from the Faraday rotation measure (Heiles et al., 1981) for ionised gas or from Zeeman spectroscopy (Brogan & Troland, 2001) for neutral/molecular gas. The plane-of-sky field components can be determined from observations of polarized emission or absorption (Kusakabe et al., 2008) from aligned spinning dust grains. In all of these techniques, the determination of is not straightforward and relies on many auxiliary assumptions. In particular, it is very difficult to determine the absolute magnitude of the plane-of-sky field components except via statistical techniques such as the Chandrasekhar-Fermi method (Chandrasekhar & Fermi, 1953; Heitsch et al., 2001; Ostriker et al., 2001; Falceta-Gonçalves, Lazarian, & Kowal, 2008). Furthermore, many of the techniques do not uniformly sample all regions along the line of sight, but are biased towards regions with particular physical conditions. We therefore caution against direct comparison of observational results with our maps in any but a qualitative sense. The maps are nonetheless useful in showing an overview of the magnetic field geometry in our simulations. In the following discussion, we describe locations on the projected map using the points of the compass, assuming that north is up (positive ) and east is left (negative ).

The most striking aspect of Figs. 20 and 21 is the large-scale order that is apparent in the projected magnetic field. This is particularly visible in the neutral-weighted maps (top-right panels), where it can be seen that the field is frequently oriented parallel to the large-scale ionisation front, forming a ring around the H  II region. The effect is particularly strong along the north and south borders of the region because of the net positive field along the -axis (see § 3.1), but it is also seen to a lesser extent along the east and west borders, despite the fact that the mean -component of the field is zero. This is because of the fast-mode MHD shock that is driven into the surrounding gas by the expanding H  II region and which compresses both the gas and the field, tending to bend the field lines so that they are more closely parallel to the shock than they were in the undisturbed medium. A large-scale pattern is much harder to discern in the molecular-weighted maps (bottom-right panels), partly because the molecular column density is much less smoothly distributed, being concentrated in globules and filaments.

In most dense filaments, the field direction is parallel to the long axis of the filament, as can be seen particularly clearly in Fig. 22, which shows a detail view of the dense photoevaporating globule found to the south in Fig. 20 and which is fed by multiple neutral/molecular filaments. In the molecular gas at the head of the globule, the magnetic field is bent into a hairpin shape. The B-field in the ionised gas at the head of the globule tends to be oriented perpendicular to the ionisation front, as is the case with almost all the dense globules visible in Fig. 20. The same is true along much of the bar-like feature to the west of the globule. On the other hand, in a few regions, such as along the filament that extends south from the globule, the B-field in the ionised gas lies along the ionisation front.

Fig. 23 shows a detail view of the same dense southern globule, but from the B star simulation shown in Fig. 21. The molecular gas shows a similar field pattern to that seen in the O star case: the field goes up one of the feeding filaments and down the other, with a hairpin bend in between. Although this can be seen clearly in the molecular gas, the neutral and ionised gas show very different patterns, with magnetic field vectors that are generally perpendicular to the filament and continuous with the large-scale field pattern in the region. This is because the total ionised and neutral columns are dominated by diffuse material along the line of sight, rather than by material associated with the globule and filament.

### 4.3 Effects of magnetic field on globule formation and evolution

It is interesting to compare the properties of the globules generated in our turbulent simulations with the results of previous detailed studies of the photoionisation of isolated dense globules (Williams, 2007b; Henney et al., 2009; Mackey & Lim, 2010b). The principal findings of the earlier studies are that a sufficiently strong magnetic field ( in the initial neutral globule) will produce important qualitative changes in the photoevaporation process. Depending on the initial field orientation with respect to the direction of the ionising photons, either extreme flattening of the globule may occur (perpendicular orientations) or the radiative implosion of the globule may be prevented (parallel orientations). For all orientations the ionised photoevaporation flow cannot freely escape from the globule, leading to recombination at late times. On the other hand, a weaker field ( in the initial globule) leads only to a moderate flattening of the globule and does not prevent the free escape of the ionised photoevaporation flow.

In our turbulent simulations, the mean initial value of is (§ 3.1), which is intermediate between the two cases discussed above and might lead one to suspect that magnetic effects on the evolution of globules should be substantial. However, this is not the case. A careful examination of the three-dimensional globule morphologies in our simulations shows no evidence of magnetically induced flattening. Although many globules do show asymmetries, this is true equally of our non-magnetic simulations and is presumably due to their irregular initial shape and internal turbulent motions. The only difference in the globule properties between our non-magnetic and magnetic simulations is that the very smallest globules do not seem to form in the magnetic case.

In order to explain this apparent discrepancy with earlier work, it is necessary to examine in more detail the distribution of magnetic field in the initial conditions of our turbulent simulations. It turns out that the dense filaments of molecular gas (from which the globules will ultimately form) tend to be much less magnetically dominated than the more diffuse gas, typically showing . Furthermore, these filaments already show a magnetic geometry similar to that described above (§ 4.2), with the field running along the long axis of the filament and with the field changing sign as one moves across the short axis. This is very different from the initial conditions assumed in the earlier globule photoevaporation studies, which were a uniform magnetic field that threaded a spherical globule (e.g., Fig. 1 of Henney et al., 2009). As a result, the compressed globule heads in our turbulent simulations show approximately equal thermal and magnetic pressures () and the magnetic effects are very modest. One caveat to this result is that numerical diffusion due to our limited spatial resolution may be producing non-physical magnetic reconnection in the globule head when oppositely directed field lines are forced together. Higher resolution studies of globule implosion with realistic initial field configurations are required in order to clarify this and these should also include self gravity, which may be important during the phase of maximum compression (Esquivel & Raga, 2007).

## 5 Conclusions

We have performed radiation-magnetohydrodynamic simulations of the formation and expansion of H  II regions and PDRs around an O star and a B star in a turbulent magnetised molecular cloud on scales of up to 4 parsec. Our principal conclusions are as follows:

1. The expansion of the H  II region is little affected by the presence of the magnetic field, since the thermal pressure of the ionised gas dominates the dynamics on the timescales of our simulations (§ 3.2.2).

2. The O star simulations produce greater fragmentation and denser clumps and filaments around the periphery of the H  II region than the B star case (§ 3.2.1, Figs. 7 and 8).

3. For B stars the non-ionising far ultraviolet radiation plays an important role in determining the morphology of the region. H  II regions around such stars are surrounded by a thick, relatively smooth shells of neutral material (PDR), of the bubble radius (e.g., lower panels of Fig. 17). In the O star simulations, the PDR is thinner and more irregular in shape (e.g., Fig. 7).

4. The resemblance at optical and longer wavelengths of our simulations to observed bubbles is striking (§ 4.1 and Deharveng et al., 2009, 2010). Our  pc radius bubbles are a typical size compared to bubbles at known distances in the GLIMPSE surveys (Churchwell et al., 2006, 2007), and so comparisons are meaningful.

5. The expanding H  II region and PDR tend to erase pre-existing small-scale disordered structure in the magnetic field, producing a large-scale ordered field in the neutral shell, with orientation approximately parallel to the ionisation front (top panels of Fig. 17 and top-right panels of Figs. 20 and 21).

6. Dense evaporating globules, pillars, and elephant trunk structures tend to be fed by two or more neutral/molecular filaments, with magnetic fields running along their length (right panels of Fig. 22 and bottom-right panel of Fig. 23). The field geometry in the neutral and molecular gas at the bright head of the structure tends to be of a hairpin shape.

7. The weak magnetic field in the ionised gas also shows an ordered structure. On the largest scales and at the latest times, it tends to align (albeit weakly in the O star case), with the mean field direction of the simulation’s initial conditions (Fig. 21, top-left panel). On smaller scales, there is a general (but not universal) tendency for the field in ionised gas to be oriented perpendicular to the local ionisation front. This tendency is more pronounced in the O star simulation (Fig. 20, top-left panel) and in globule evaporation flows for both simulations (top-left panels of Figs. 22 and 23).

8. The highest pressure compressed neutral/molecular gas shows approximate equipartition between thermal, turbulent, and magnetic energy density, whereas lower pressure gas (either neutral or molecular) tends to separate into, on the one hand, magnetically dominated, quiescent regions, and, on the other hand, demagnetised, highly turbulent regions (§ 3.3, Fig. 19). The lower pressure gas also separates into low-, magnetically dominated regions (which are largely molecular) and high-, thermally dominated regions (which are largely neutral). The ionised gas, on the other hand, always shows approximate equipartition between thermal and turbulent energies, but with the magnetic energy being lower by 1 to 3 orders of magnitude.

9. Velocity dispersions in the ionised gas of   are maintained for the entire duration of all our simulations (§ 3.2.2, Fig. 15 and 16). This is 5 to 10 times higher than the value that would be predicted by expansion in a uniform medium. At early times ( yr for the O star, or  yr for the B star), this dispersion is mainly due to radial champagne-flow expansion as the H  II region escapes from its natal clump. At later times, the net radial expansion of ionised gas subsides, but the velocity dispersion is maintained by inwardly-directed photoevaporation flows from globules and pillars.

## Acknowledgments

SJA and WJH acknowledge financial support from DGAPA PAPIIT projects IN112006, IN100309 and IN110108. This work was supported in part by Swedish Research Council grant 2009-4088. Some of the numerical calculations in this paper were performed on the Kan Balam supercomputer maintained and operated by DGSCA, UNAM. This work has made extensive use of NASA’s Astrophysics Abstract Data Service and the astro-ph archive.

## References

• Anderson et al. (2010) Anderson L. D., et al., 2010, A&A, 518, L99
• Arthur & Hoare (2006) Arthur S. J., Hoare M. G., 2006, ApJS, 165, 283
• Arthur et al. (2004) Arthur S. J., Kurtz S. E., Franco J., Albarrán M. Y., 2004, ApJ, 608, 282
• Bakes, Tielens, & Bauschlicher (2001) Bakes E. L. O., Tielens A. G. G. M., Bauschlicher C. W., Jr., 2001, ApJ, 556, 501
• Ballesteros-Paredes et al. (1999) Ballesteros-Paredes, J., Vázquez-Semadeni, E., & Scalo, J. 1999, ApJ, 515, 286
• Balsara & Spicer (1999) Balsara D. S., Spicer D., 1999, Journal of Computational Physics, 148, 133
• Bertoldi (1989) Bertoldi F., 1989, ApJ, 346, 735
• Bregman et al. (1989) Bregman J. D., Allamandola L. J., Witteborn F. C., Tielens A. G. G. M., Geballe T. R., 1989, ApJ, 344, 791
• Brogan & Troland (2001) Brogan C. L., Troland T. H., 2001, ApJ, 560, 821
• Carlqvist et al. (2003) Carlqvist P., Gahm G. F., Kristen H., 2003, A&A, 403, 399
• Chandrasekhar (1960) Chandrasekhar S., 1960, Radiative transfer. New York: Dover
• Chandrasekhar & Fermi (1953) Chandrasekhar S., Fermi E., 1953, ApJ, 118, 113
• Churchwell et al. (2006) Churchwell E., et al., 2006, ApJ, 649, 759
• Churchwell et al. (2007) Churchwell E., et al., 2007, ApJ, 670, 428
• Colella & Woodward (1984) Colella P., Woodward P. R., 1984, Journal of Computational Physics, 54, 174
• Crutcher, Hakobian, & Troland (2009) Crutcher R. M., Hakobian N., Troland T. H., 2009, ApJ, 692, 844
• Crutcher et al. (2010) Crutcher R. M., Wandelt B., Heiles C., Falgarone E., Troland T. H., 2010, ApJ, 725, 466
• Dale et al. (2007) Dale J. E., Clark P. C., Bonnell I. A., 2007, MNRAS, 377, 535
• De Colle (2005) De Colle F., 2005, PhD thesis, UNAM
• De Colle & Raga (2004) De Colle F., Raga A., 2004, Ap&SS, 293, 173
• De Colle & Raga (2006) De Colle F., Raga A. C. 2006, A&A, 449, 1061
• Deharveng et al. (2009) Deharveng L., Zavagno A., Schuller F., Caplan J., Pomarès M., De Breuck C., 2009, A&A, 496, 177
• Deharveng et al. (2010) Deharveng L., et al., 2010, A&A, 523, A6
• Dyson & Williams (1997) Dyson J. E., Williams D. A., 1997, The Physics of the Interstellar Medium, 2nd Edition. (Bristol: Institute of Physics Publishing)
• Elmegreen (2000) Elmegreen B. G., 2000, ApJ, 530, 277
• Esquivel & Raga (2007) Esquivel A., Raga A. C., 2007, MNRAS, 377, 383
• Falceta-Gonçalves, Lazarian, & Kowal (2008) Falceta-Gonçalves D., Lazarian A., Kowal G., 2008, ApJ, 679, 537
• Falle et al. (1998) Falle S. A. E. G., Komissarov S. S., Joarder P., 1998, MNRAS, 297, 265
• Fatuzzo & Adams (2008) Fatuzzo M., Adams F. C., 2008, ApJ, 675, 1361
• Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761
• Garcia-Segura & Franco (1996) Garcia-Segura G., Franco J., 1996, ApJ, 469, 171
• Goodman & Heiles (1994) Goodman A. A., Heiles C. 1994, ApJ, 424, 208
• Gritschneder et al. (2010) Gritschneder M., Burkert A., Naab T., Walch S., 2010, ApJ, 723, 971
• Heiles et al. (1981) Heiles C., Chu Y.-H., Troland T. H., 1981, ApJ, 247, L77
• Heiles & Crutcher (2005) Heiles C., Crutcher R., 2005, Cosmic Magnetic Fields, 664, 137
• Heitsch et al. (2001) Heitsch, F., Zweibel, E. G., Mac Low, M.-M., Li, P., & Norman, M. L. 2001, ApJ, 561, 800
• Henney (2003) Henney W. J., 2003, Revista Mexicana de Astronomia y Astrofisica Conference Series, 15, 175
• Henney et al. (2009) Henney W. J., Arthur S. J., De Colle F., Mellema G., 2009, MNRAS, 398, 157
• Iliev et al. (2006) Iliev I. T., et al., 2006, MNRAS, 371, 1057
• Iliev et al. (2009) Iliev I. T., et al., 2009, MNRAS, 400, 1283
• Inoue (2002) Inoue A. K., 2002, ApJ, 570, 688
• Kahn (1954) Kahn F. D., 1954, Bull. Astron. Inst. Netherlands, 12, 187
• Kusakabe et al. (2008) Kusakabe N., et al., 2008, AJ, 136, 621
• Krumholz et al. (2007) Krumholz M. R., Stone J. M., Gardiner T. A., 2007, ApJ, 671, 518
• Krumholz & Matzner (2009) Krumholz M. R., Matzner C. D., 2009, ApJ, 703, 1352
• Li et al. (2004) Li Y., Mac Low M.-M., Abel T., 2004, ApJ, 610, 339
• Mac Low et al. (2007) Mac Low M.-M., Toraskar J., Oishi J. S., Abel T., 2007, ApJ, 668, 980
• Mackey & Lim (2010a) Mackey J., Lim A. J., 2010a, MNRAS, 403, 714
• Mackey & Lim (2010b) Mackey J., Lim A. J., 2010b, arXiv, arXiv:1012.1500
• Martins et al. (2010) Martins F., Pomarès M., Deharveng L., Zavagno A., Bouret J. C., 2010, A&A, 510, A32
• Mathews (1967) Mathews W. G., 1967, ApJ, 147, 965
• Matzner (2002) Matzner C. D., 2002, ApJ, 566, 302
• Krumholz, Matzner, & McKee (2006) Krumholz M. R., Matzner C. D., McKee C. F., 2006, ApJ, 653, 361
• McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565
• Lebouteiller et al. (2007) Lebouteiller V., Brandl B., Bernard-Salas J., Devost D., Houck J. R., 2007, ApJ, 665, 390
• Mellema et al. (2006a) Mellema G., Arthur S. J., Henney W. J., Iliev I. T., Shapiro P. R., 2006a, ApJ, 647, 397
• Mellema et al. (2006b) Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006b, New Astronomy, 11, 374
• Mellema & Lundqvist (2002) Mellema G., Lundqvist P., 2002, A&A, 394, 901
• Menon (1962) Menon T. K., 1962, ApJ, 135, 394
• Motoyama et al. (2007) Motoyama K., Umemoto T., Shang H., 2007, A&A, 467, 657
• Mouschovias & Ciolek (1999) Mouschovias T. C., Ciolek G. E., 1999, NATO ASIC Proc. 540: The Origin of Stars and Planetary Systems, eds. C. J. Lada, N. D. Kylafis. (Dordrecht: Kluwer),305
• Osterbrock & Ferland (2006) Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei, 2nd. ed. by D.E. Osterbrock and G.J. Ferland, Sausalito, CA: University Science Books, 2006
• Ostriker et al. (2001) Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980
• Peters et al. (2010) Peters, T., Banerjee, R., Klessen, R. S., & Mac Low, M.-M. 2010, arXiv:1010.5905
• Redman et al. (1998) Redman M. P., Williams R. J. R., Dyson J. E., Hartquist T. W., Fernandez B. R., 1998, A&A, 331, 1099
• Ryutov et al. (2005) Ryutov D. D., Kane J. O., Mizuta A., Pound M. W., Remington B. A., 2005, Ap&SS, 298, 183
• Shu et al. (2002) Shu F. H., Lizano S., Galli D., Cantó J., Laughlin G., 2002, ApJ, 580, 969
• Spitzer (1978) Spitzer L., 1978, New York Wiley-Interscience, 1978
• Strömgren (1939) Strömgren B., 1939, ApJ, 89, 526
• Tang et al. (2009) Tang Y.-H., Ho P. T. P., Girart J. M., Rao R., Koch P., Lai S.-P., 2009, ApJ, 695, 1399
• Tasker & Tan (2009) Tasker E. J., Tan J. C., 2009, ApJ, 700, 358
• Tóth (2000) Tóth G., 2000, Journal of Computational Physics, 161, 605
• Vacca et al. (1996) Vacca W. D., Garmany C. D., Shull J. M., 1996, ApJ, 460, 914
• Vázquez-Semadeni et al. (2005) Vázquez-Semadeni E., Kim J., Shadmehri M., Ballesteros-Paredes J., 2005, ApJ, 618, 344
• Watson et al. (2008) Watson C., et al., 2008, ApJ, 681, 1341
• Williams (2007a) Williams R. J. R., 2007a, Diffuse Matter from Star Forming Regions to Active Galaxies, 129
• Williams (2007b) Williams R. J. R., 2007b, Ap&SS, 307, 179
• Williams & Dyson (2001) Williams R. J. R., Dyson J. E., 2001, MNRAS, 325, 293
• Williams et al. (2000) Williams R. J. R., Dyson J. E., Hartquist T. W., 2000, MNRAS, 314, 315
• Williams et al. (2001) Williams R. J. R., Ward-Thompson D., Whitworth A. P., 2001, MNRAS, 327, 788
• Wood & Churchwell (1989) Wood D. O. S., Churchwell E., 1989, ApJS, 69, 831
• Zavagno et al. (2007) Zavagno A., Pomarès M., Deharveng L., Hosokawa T., Russeil D., Caplan J., 2007, A&A, 472, 835
• Zavagno et al. (2010) Zavagno A., et al., 2010, A&A, 518, L81

## Appendix A Integration of plane-of-sky field components

Polarisation-based methods for measuring the plane-of-sky components of the magnetic field cannot distinguish the sense of the field and are degenerate between a position angle and a position angle . It is therefore not sufficient to simply integrate and along the line of sight with appropriate weighting. Neither is it sufficient to integrate and , since this will always give a result in the first quadrant (). In order to calculate our projected magnetic field maps (§ 4.2), we therefore adopt a “Stokes parameter” approach (Chandrasekhar, 1960), whereby the plane-of-sky components and are first transformed to ()-space as , , where and . Then the line-of-sight integration is performed as , where is the relevant density (ionised, neutral, or molecular). Finally, the components of the projected field are transformed back to ()-space: , , where and .

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters