The Dynamics of Ultracompact HII Regions
AAccepted for publication in MNRAS
The Dynamics of Ultracompact HII Regions
Nathaniel Roth Steven W. Stahler Eric Keto [email protected] ABSTRACT
Many ultracompact HII regions exhibit a cometary morphology in radio continuumemission. In such regions, a young massive star is probably ablating, through its ul-traviolet radiation, the molecular cloud clump that spawned it. On one side of thestar, the radiation drives an ionization front that stalls in dense molecular gas. On theother side, ionized gas streams outward into the more rarefied environment. This windis underpressured with respect to the neutral gas. The difference in pressure draws inmore cloud material, feeding the wind until the densest molecular gas is dissipated.Recent, time-dependent simulations of massive stars turning on within molecular gasshow the system evolving in a direction similar to that just described. Here, we explorea semi-analytic model in which the wind is axisymmetric and has already achieved asteady state. Adoption of this simplified picture allows us to study the dependence ofboth the wind and its bounding ionization front on the stellar luminosity, the peakmolecular density, and the displacement of the star from the center of the clump.For typical parameter values, the wind accelerates transonically to a speed of about15 km s − , and transports mass outward at a rate of 10 − M (cid:12) yr − . Stellar radiationpressure acts to steepen the density gradient of the wind. Subject headings:
ISM: HII regions, clouds, jets and outflows — stars: formation, early-type Dept. of Physics, U. of California, Berkeley, CA 94720 Dept. of Astronomy, U. of California, Berkeley, CA 94720 Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138 a r X i v : . [ a s t r o - ph . GA ] N ov
1. Introduction1.1. Observational Background
An ultracompact HII (UCHII) region is one of the earliest signposts for the presence of a young,massive star (for reviews, see Churchwell 2002 and Hoare et al. 2007). While the star itself is stilltoo embedded in its parent molecular cloud to be detected optically, it heats up surrounding dust.A small region, some 10 cm in extent, glows brightly in the far infrared. Ionization of ambient gasalso creates free-free emission in the radio continuum, with electron densities in excess of 10 cm − and thus emission measures of 10 pc cm − or more. It is through this radio emission, relativelyminor in the overall energy budget, that UCHII regions have been classified morphologically.In their pioneering radio interferometric survey, Wood & Churchwell (1989) found that 30 per-cent of the spatially resolved regions have a cometary shape (see also Kurtz et al. 1994, Walsh et al.1998). One sees a bright arc of emission filled in by an extended, lower-intensity lobe that fadesaway from the arc. OH masers may be found along the bright rim. Other UCHII regions exhibitonly the emission arc; presumably the interior lobe in these cases is undetectably faint. Still otherclasses identified by Wood & Churchwell (1989) include: spherical, core-halo (a bright peak sur-rounded by a fainter envelope), shell (a ring of emission), and irregular (multiple emission peaks).All told, the cometary morphology is the most common one found, and needs to be explained byany viable theoretical model.Spectral line studies have been used to probe the kinematics of these regions. Observations ofradio recombination lines (e.g. Afflerbach et al. 1996; Kim & Koo 2001), and infrared fine-structurelines (e.g. Zhu et al. 2005) reveal large line-widths, indicative of supersonic flow. In cases wherethe flow can be spatially resolved, one also sees a velocity gradient. This gradient is steepest in the“head-to-tail” direction (Garay et al. 1994; Garay & Lizano 1999).Very often, the observed peak in radio continuum or OH maser emission, either of whicheffectively locates the star, does not coincide with the peak in molecular lines or submillimetercontinuum emission, which trace the densest molecular gas and dust (e.g. Mueller et al. 2002;Thompson et al. 2006). This gas is located within infrared dark clouds, currently believed to be thebirth sites of all massive stars (Beuther et al. 2007). The clumps within these clouds have typicalsizes of 1 pc, number densities of 10 cm − and masses of about 10 M (cid:12) ; some qualify as hot cores(Hofner et al. 2000; Hoare et al. 2007). The offset of the peak radio emission of the UHCII regionfrom the center of this clump is typically a few arcseconds, corresponding to approximately 0.1 pcfor a distance of 1 kpc. Both the cometary morphology and the acceleration of ionized gas arelikely related to this physical displacement, as was first emphasised by Kim & Koo (2001). 3 – The foregoing observations, taken together, show convincingly that the cometary structuresrepresent ionized gas accelerating away from the densest portion of the nearby cloud material.Theorists have long considered photoevaporating flows created by a massive star illuminating oneface of a cloud (Kahn 1954; Oort & Spitzer 1955). The most well-studied classical HII region of all,the Orion Nebula, is a prime example of the resulting “blister,” here formed by the massive star θ C on the surface of the Orion A molecular cloud (Zuckerman 1973). A tenuous, hemisphericalbody of ionized gas surrounds θ C and the other Trapezium stars, and is flowing away from thebackground cloud.When massive stars are still deeply embedded in the densest portion of their parent cloud, itis not obvious how such photoevaporating flows can be maintained. Wood & Churchwell (1989)pointed out that the small size of UCHII regions suggests a brief dynamical lifetime. If they undergopressure-driven expansion at ∼
10 km s − , they will expand to a size greater than 0.1 pc in roughly10 yr, or one percent of the lifetime of the host star. In reality, some 10 percent of O stars areassociated with UCHII regions, suggesting that the lifetime of these regions is larger by an orderof magnitude. Confinement by thermal pressure alone would result in an emission measure evenhigher than is observed (Xie et al. 1996). In the context of cometary structures, this venerable“lifetime problem” raises a fundamental question. What reservoir of matter can feed the ionizedflows over a period of 10 yr?Hollenbach et al. (1994) suggested that the star may be photoevaporating its own accretiondisk. When ultraviolet radiation from a massive star impinges on the disk, gas streams off at thesound speed, at least in that outer region where this speed exceeds the local escape velocity. Lugoet al. (2004) have analyzed this launching process in more detail. The outer accretion disk radiusof approximately 10 cm is much smaller than the 10 cm size of UCHII regions. Thus, while themodel views a disk as the ultimate source of matter for the ionized flow, it does not address theflow’s cometary morphology.One possibility is that the prominent arc represents the shock interface between a high-velocitystellar wind and the parent cloud. Massive stars indeed generate, through radiative acceleration,winds with terminal velocities of about 1000 km s − . If the star itself moves through the cloudat supersonic speed, e.g. 20 km s − , then the curved bowshock has the right form (van Burenet al. 1990; Mac Low et al. 1991). Moreover, a star that is moving toward the observer createsa “core-halo” structure, also commonly seen. The relatively fast stellar velocity, however, impliesthat the whole interaction would last for a few times 10 years if the star indeed travels througha molecular clump 1 pc in size. Consequently, this bowshock model may fall prey to the lifetimeproblem. Another concern is that massive stars, with the exception of runaway objects, do notmove at such high speeds relative to parent molecular gas. For instance, the Trapezium stars inthe Orion Nebula Cluster have a velocity dispersion of a few km s − (F˝ur´esz et al. 2008).Suppose instead that the embedded star is not moving with respect to the densest gas, but is 4 –displaced from it, as is suggested by the observations. Then one side of its expanding HII regioneventually erupts into the surrounding low-density medium. Such a “champagne flow” modelwas first explored by Tenorio-Tagle (1979), Bodenheimer et al. (1979), and Whitworth (1979) toexplain asymmetric, classical HII regions. The high internal pressure of the HII region acceleratesthe ionized gas to supersonic speed away from the ambient cloud, whose density was taken to be10 cm − . Meanwhile, the ionization front steadily advances at several km s − into this cloud,generating volumes of ionized gas several parsecs in diameter.Another model, that of the mass-loaded wind (Dyson 1968), invokes champagne-type dynamicsin combination with a stellar wind. The idea is that the wind, perhaps in combination with thestellar radiation field, ablates the cloud, entraining the gas within it. Pressure gradients may thenaccelerate the gas in a champagne flow. In some versions of the model, the new mass originates indense globules that are continually being ionized (Williams et al. 1996; Lizano et al. 1996; Redmanet al. 1998). While successful in explaining the morphologies and lifetimes of UCHII regions, themodel does not address in detail how the dense globules enter the ionized flow. Thus, the mass-loading prescription adopted was somewhat ad hoc.Keto (2002a) modeled the growth of an HII region inside a Bondi accretion flow. If theionization front is located inside the Bondi radius r B ≡ GM ∗ / a I , where M ∗ is the stellar massand a I the ionized sound speed, then gas simply crosses the ionization front and continues towardthe star as an ionized accretion flow. Until its size reaches r B , the HII region can expand only asthe ionizing flux from the star increases along with the stellar mass. The HII region then exists insteady state, continuously fed by the molecular accretion. Thus, the lifetime of the HII region istied to the accretion time scale of the star.Observations of the UCHII region around the cluster of massive stars G10.6-0.4 indicate thatthis kind of molecular and ionized accretion could be occurring (Keto 2002b). Subsequent ob-servations of the same cluster (Keto & Wood 2006) also show an asymmetric bipolar outflow ofionized gas. These observations motivated Keto (2007) to develop a model in which inflow andoutflow occur simultaneously in a rotationally flattened geometry. The ionization creates an HIIregion elongated perpendicular to the accretion flow. Where the ionization extends beyond r B , theHII region can expand hydrodynamically as a pressure-driven Parker wind (see also McKee & Tan(2008)). Along the equatorial plane, dense molecular gas continues to flow into the HII region.One of the motivations of the present paper is to detail exactly how an ionized outflow may besupplied by inflow of molecular gas, although here we consider a mechanism to draw in this gasthat is separate from the gravitational attraction of the star.Another model combining ionization and gravitational accretion was that of Mac Low et al.(2007), who suggested that an UCHII region may form as a gravitationally collapsing substructurewithin a larger, expanding HII region. More detailed simulations by Peters et al. (2010a,b) foundthat an UCHII region embedded in an accretion flow rapidly changes morphologies through all theobserved types, and could be sustained by the addition of infalling gas from the parent cloud. The 5 –collapsing ionized gas in their simulations creates bipolar molecular outflows, as are observed toaccompany some UCHII regions (Beuther & Shepherd 2005). Peters et al. (2011) also claimed thatmagnetic pressure might play a significant role in confining UCHII regions. On the other hand,Arthur et al. (2011) performed their own radiation-MHD simulation of an HII region expandinginto a turbulent cloud. They found that magnetic pressure plays only a minor role in confinement.Other recent simulations have combined various elements in the original models to account forthe proliferating observational results now becoming available. Henney et al. (2005) simulated achampagne flow in which the ionization front is stalled as it climbs a density gradient into a neutralcloud supported by turbulent pressure. Comeron (1997) explored the coupled roles of champagneflows and stellar winds. Finally, Arthur & Hoare (2006) revised the bowshock idea by introducinga finite but modest stellar velocity with respect to the cloud.In the present paper we pursue a different goal. Accepting that every cometary UCHII regionis an ionized flow, we elucidate the basic issue of how such a flow may continue to draw in gas fromthe neutral cloud. We agree with previous researchers that the momentum of the ionized flow isnot supplied by the star, but by the thermal pressure gradient of the gas itself. A key result of ourown study will be to demonstrate that the expansion of this flow causes the ionized gas to become underpressured with respect to the neutral cloud. This pressure differential draws neutral materialfrom the cloud into the champagne flow in a self-sustaining manner.We assume at the outset that the flow is quasi-steady, as was found in the late stages of thesimulations of Henney et al. (2005) and Arthur & Hoare (2006). This basic simplification allows usto rapidly explore parameter space, as we detail below. Thus, we can assess how well our ratherminimal set of physical assumptions explains the basic characteristics of cometary UCHII regions.Our quasi-one dimensional model does not allow us to include evacuation by a stellar wind, but wedo account for the effect of radiation pressure and show that it is appreciable.In Section 2 below, we introduce our steady-state model for cometary UCHII regions. Section 3develops the equations governing the density and velocity of the ionized flow. In Section 4, we recastthese equations in nondimensional form and then outline our solution strategy. Section 5 presentsour numerical results, and Section 6 compares them to observations. Finally, Section 7 indicatesfruitful directions for future investigations.
2. Steady-state Model2.1. Physical Picture
We idealize the molecular cloud as a planar slab, in which the cloud density peaks at themidplane (see Figure 1). The choice of planar geometry is made for computational convenienceand is probably not a realistic representation of the clouds of interest. However, our qualitativeresults are not sensitive to the adopted cloud density profile, as we verify explicitly in Section 5 6 –below. We envision a massive star embedded within this cloud, but offset from the midplane, inaccordance with the observations mentioned previously. This offset creates an asymmetric, ionizedregion of relatively low density. Toward the cloud midplane, this ionized gas resembles a classical HIIregion. Its boundary, a D-type ionization front, advances up the density gradient. In the oppositedirection, the front breaks free and gas streams away supersonically. The pressure gradient withinthe ionized gas creates an accelerating flow away from the cloud. Apart from the higher-densitycloud environment, our model is identical in spirit to the champagne flows proposed in the past.Once the velocity of the advancing ionization front falls significantly below the ionized soundspeed, the flow becomes steady-state. We adopt this steady-state assumption in our model, andassume that the structure of the background cloud is evolving over a long timescale compared tothe ionized sound crossing time. In the simulations of Arthur & Hoare (2006), a steady flow isapproached some 10 years after the massive star turns on. By this point, the ionization front hasvirtually come to rest in the frame of the star. We will later determine more precisely the speed ofthe ionization front in our model and show that it is consistent with the steady-state assumption.The shock that originally formed ahead of the ionization front, as it transitioned from R-type toD-type, has by now advanced deep into the cloud and died away.The HII region facing the cloud midplane expands, albeit slowly. While the ionized gas isoptically thick to ultraviolet radiation from the star, some does leak through and strikes the cavitywall within the neutral cloud. Additional gas is thus dissociated and ionized, and streams offthe wall to join the outward flow. While the injection speed at the wall is subsonic, the thermalpressure gradient within the ionized gas accelerates it to supersonic velocity. Both this flow andthe advancing front erode the cloud, whose structure gradually evolves in response. We will always consider systems that possess azimuthal symmetry. We establish a sphericalcoordinate system whose origin is at the star. The polar direction, θ = 0, coincides with the centralaxis of the flow depicted in Figure 1. Let N ∗ denote the total number of ionizing photons per timegenerated by the star. If we further assume ionization balance within the volume of the ionizedcavity, and make use of the on-the-spot approximation, then the ionizing radiation extends out tothe cavity wall r f ( θ ), given implicitly by N ∗ π = α B (cid:90) r f dr r [ n I ( r, θ ) / + r f F ∗ , wall ( θ ) . (1)Here, α B is the case B recombination coefficient (Osterbrock 1989 Chapter 2). We have pulled thisfactor out of the integral since it depends only on the ionized gas temperature, which we assumeto be spatially uniform. The term n I is the number density of ionized gas, including both protonsand electrons. For simplicity, we have assumed that the composition of the gas is pure hydrogen,so the number of free protons is equal to the number of free electrons and both of these number 7 –densities are equal to n I /
2. Finally, F ∗ , wall ( θ ) is that portion of the star’s ionizing photon fluxwhich escapes the HII region and strikes the wall of the cavity. This remnant flux is critical formaintaining the flow.But how important is the second term of equation (1) in a quantitative sense? Each photonstriking the front ionizes a hydrogen atom, itself already dissociated from a hydrogen molecule,thereby creating a proton and an electron. Let f be the number flux of protons and electronsinjected into the flow at each position along the front. Thus, F ∗ , wall = f / θ , and wemay write the foregoing equation as N ∗ π = α B (cid:90) r f dr r ( n I / + r f f . (2)The flux f is of order n I a I , where n I is now a typical value of the ionized gas number density.If r f in equation (2) now represents the size scale of the flow, then the ratio of the second to thefirst righthand term in this equation is of order a I α B n I r f ∼ (cid:18) a I a cl (cid:19) a cl α B n cl r f , (3)where a cl and n cl are the effective sound speed and number density within the cloud (see below).Here we have assumed that the ionized material is at least in rough pressure balance with the cloud,so that n I a I ≈ n cl a . For the dense clumps that harbor young massive stars, n cl ∼ cm − and a cl ∼ − (Garay & Lizano 1999). Further using r f ∼ cm, a I ∼
10 km s − , and α B ∼ − cm s − , we find that the ratio of terms is of order 10 − . In practice, therefore, weneglect the second term entirely. We trace the ionization front by finding that function r f ( θ )which obeys the approximate, but quantitatively accurate relation, N ∗ π = α B (cid:90) r f dr r ( n I / . (4) Even within our simplifying assumption of axisymmetry, it would be a daunting task to trace afully two-dimensional flow. While we accurately follow the ionization front bounding the flow in twodimensions, we further assume that the density and velocity are only functions of z . In this quasione-dimensional model, we are effectively averaging the ionized density n I and velocity u laterally ateach z -value. This simplification is innocuous in regions where the lateral change in these quantitiesis relatively small. It is more problematic when the ionized gas becomes overpressured with respect In the terminology of Henney (2001), our photoevaporation flow is recombination-dominated . If the secondrighthand term in equation (2) were relatively large, the flow would be advection-dominated . According to Henney(2001), knots in planetary nebulae fall into this category. a cl . In our model, it is only the mass of the cloud that affectsthe flow gravitationally. That is, we neglect both the self-gravity of the ionized gas, and the pull ofthe massive star. The latter force is negligible outside the Bondi radius R B ≡ G M ∗ / a I , where M ∗ is the stellar mass. For representative values M ∗ = 20 M (cid:12) and a I = 10 km s − , R B = 90 AU,much less than the size of UCHII regions.Our cloud is an isothermal self-gravitating slab whose midplane is located at z = − H ∗ (seeFig. 1). Solving the equations of hydrostatic equilibrium along with Poisson’s equation yields theneutral cloud density n cl and the gravitational potential Φ cl n cl = n sech (cid:18) z + H ∗ H cl (cid:19) , (5)Φ cl = − a ln sech (cid:18) z + H ∗ H cl (cid:19) . (6)Here, n is the midplane number density of hydrogen molecules, each of mass 2 m H , and H cl is thescale thickness of the slab: H cl ≡ (cid:20) a π G (2 n m H ) (cid:21) / . (7)A representative value for n is 10 cm − . Combining this with a turbulent velocity dispersionof a cl = 2 km s − yields a Jeans mass on the order of 10 M (cid:12) , far less than the 10 M (cid:12) clumpsin giant molecular clouds spawning massive stars and their surrounding clusters. We stress againthat our cloud is a relatively small fragment containing the newborn star. In our view, it is thehigh density of this gas that determines the morphology of the UCHII region, and the dispersal ofthe clump that sets the characteristic lifetime. Quantitative modeling of HII regions, and UCHII regions in particular, has generally neglectedthe dynamical effect of radiation pressure from the massive star (see, however, Krumholz & Matzner2009). For the very dense environments we are now considering, the radiative force has substantialinfluence on the ionized flow, as we shall demonstrate through explicit calculation.Suppose the star emits photons with mean energy (cid:15).
Those traveling in the direction θ withrespect to the central axis carry momentum ( (cid:15)/c ) cos θ in the z -direction. If we assume that the gas 9 –is in ionization equilibrium, then the number of photons absorbed per unit volume of gas equalsthe corresponding volumetric rate of recombination, ( n I / α B . The assumption of ionizationequilibrium is justified because the typical recombination time, t rec ∼ ( n I α B ) − ∼ s, is muchless than the flow time t flow ∼ r f /a I ∼ s.The radiative force per volume is the product of the photon momentum and the volumetricrate of ionization, or ( (cid:15)/c ) cos θ ( n I / α B . To obtain f rad , the radiative force per unit mass of gas,we divide this expression by the ionized mass density, m H n I /
2. We thus find f rad = (cid:15) (cid:104) cos θ (cid:105) n I α B m H c . (8)This expression is consistent with that given in Krumholz & Matzner (2009, Section 2) if we addthe condition of ionization balance. Note, finally, that (cid:15) in equation (8) is implicitly a function of N ∗ , a fact that we shall use later.In accordance with our quasi one-dimensional treatment, we have laterally averaged cos θ atfixed z . Explicitly, this average, weighted by the cross-sectional area is (cid:104) cos θ (cid:105) = 2 η z R (cid:32)(cid:114) R z − (cid:33) . (9)Here, R ( z ) is the cylindrical radius from the central axis to the ionization front at each z , and η = +1 or − z , respectively. At z = 0, the level of the star, we set (cid:104) cos θ (cid:105) = 0. At the base of the flow, (cid:104) cos θ (cid:105) = −
1. We define the distance from this point to thestar as H b ≡ r f ( π ), and show it in Figure 1.
3. Flow equations3.1. Mass and Momentum Conservation
To obtain laterally averaged dynamical equations, consider a control volume of height ∆ z spanning the flow, as pictured in Figure 2. This volume has cylindrical radius R and R + ∆ R atits lower and upper surfaces, respectively. Let u I ( z ) be the average flow speed in the z -direction.Similarly, let n inj and v inj represent the number density and speed, respectively, of ionized gas beinginjected into the flow just inside the ionization front. Then the requirement of mass conservationis ddt ( π R ∆ z n I ) = π R n I u I − π ( R +∆ R ) ( n I +∆ n I )( u I +∆ u I )+2 π R (cid:112) ∆ R + ∆ z n inj v inj . (10)Here, the first two terms of the righthand side are the rate of mass advection through the bottomand top layers of the control volume, respectively. The final righthand term is the rate of mass 10 –injection. We assume that the injected flow direction is normal to the ionization front. Henceforth,we will drop the subscript I when referring to the density and velocity of the ionized flow.Under our steady-state assumption, the lefthand side of equation (10) vanishes. Dividingthrough by ∆ z , and taking the limit ∆ z →
0, leads to ddz ( R n u ) = 2 R (cid:115) (cid:18) dRdz (cid:19) n inj v inj . (11)We will later derive an expression for v inj from the jump conditions across the ionization front andshow that this velocity is subsonic, i.e., v inj < a I at any z . Although our derived v inj is properlymeasured in the rest frame of the ionization front, we will show in Section 3.2 that the front ismoving very slowly compared to a I . Thus, v inj is also, to a high degree of accuracy, the speed ofthe injected gas in the rest frame of the star.Another crucial assumption we will make is that n inj = n , i.e., the injected density is thesame as the laterally averaged flow value at any height. This condition seems to hold at leastapproximately in the simulations of Arthur & Hoare (2006), and is physically plausible when oneconsiders that lateral density gradients will tend to be smoothed out if the velocity components inthat direction are subsonic. After making this assumption, we are left with ddz ( R n u ) = 2 R (cid:115) (cid:18) dRdz (cid:19) n v inj . (12)Requiring conservation of z -momentum for the control volume leads to ddt ( π R ∆ z n u ) = π R n u − π ( R + ∆ R ) ( n + ∆ n )( u + ∆ u ) + 2 π R n v ∆ R + π R n a I − π ( R + ∆ R ) ( n + ∆ n ) a I + 2 π R ∆ R + ∆ z n a I ∆ R − π R ∆ z n d Φ cl dz + π R ∆ z n f rad , (13)where Φ cl is the gravitational potential from equation (6) and f rad is the radiative force per massfrom equation (8). Included are terms representing both static pressure and the advection ofmomentum through the top, bottom, and sides of the control volume. For the advective terms,we have again replaced n inj by n . The geometric factor (cid:112) dR/dz ) present in equation (12)disappears because its inverse is used when projecting the injected momentum into the z -direction.We apply the steady-state condition and divide equation (13) through by πR n ∆ z . After takingthe ∆ z → u dudz = − a I n dndz − d Φ cl dz + f rad + 2 R dRdz v − R (cid:115) (cid:18) dRdz (cid:19) v inj u . (14)This equation resembles the standard Euler momentum equation in one dimension, but has twoadditional terms. The first accounts for injection of z -momentum via ram pressure. The secondrepresents the inertial effect of mass loading. 11 – In the rest frame of the ionization front, upstream molecular gas approaches at speed v (cid:48) cl andleaves downstream as ionized gas, at speed v (cid:48) inj . We assume that the intermediate photodissociationregion, consisting of neutral hydrogen atoms, is geometrically thin. Conservation of mass andmomentum across the ionization front is expressed in the jump conditions(1 / n v (cid:48) inj = 2 n cl v (cid:48) cl (15a)(1 / n ( a I + v (cid:48) inj2 ) = 2 n cl ( a + v (cid:48) cl2 ) . (15b)The factors of 1 / m H , while each particle of ionized gas has a mean mass of m H /
2. We solve equation(15a) for v (cid:48) cl v (cid:48) cl = 14 nn cl v (cid:48) inj , (16)and use this result in equation (15b) to derive an expression for v (cid:48) inj : v (cid:48) inj = a I (cid:115) n /β − n n cl n n cl − n , (17)where β ≡ a I /a (cid:29)
1. As long as the ionized gas is underpressured with respect to the neutralmolecular gas, β n < n cl , and the quantity inside the square root in equation (17) is positive,guaranteeing a solution for v (cid:48) inj . Using this inequality, equation (16) tells us that v (cid:48) cl < v (cid:48) inj /β . Inall our solutions, v (cid:48) inj is subsonic. Thus, the ionization front moves relatively slowly into the cloud,and we may set the lab frame injection velocity v inj equal to v (cid:48) inj . Since we know n cl at all z ,equation (17) gives us v inj as a function of z and n , which we may use in the mass and momentumequations (12) and (14). The mass and momentum conservation equations can be combined to solve separately for thederivatives of the velocity and density. These decoupled equations are dudz = (cid:18) a I − u (cid:19) − u R dRdz ( a I + v ) + 2 v inj R (cid:115) (cid:18) dRdz (cid:19) ( a I + u ) + u d Φ cl dz − u f rad (18) dndz = (cid:18) a I − u (cid:19) n R dRdz ( u + v ) − n R (cid:115) (cid:18) dRdz (cid:19) u v inj − n d Φ cl dz + n f rad . (19) Roshi et al. (2005, Section 4.2) estimate the PDR thickness in G35.20-1.74 to be of order 10 − pc.
12 –The righthand sides of both equations have denominators that vanish when u = a I . Thus, wemust take special care when integrating through the sonic point, as is also true in steady-state windsand accretion flows of an isothermal gas. If we were to ignore the terms relating to mass injection,radiation pressure, and self-gravity, then the sonic transition would occur when dR/dz = 0, as in ade Laval nozzle. In our case, the extra terms cause the sonic transition to occur at other locations.
4. Nondimensionalization and solution strategy4.1. Characteristic scales
A represenative ionizing photon emission rate is 10 s − , which corresponds to an O7.5 star(Vacca et al. 1996). We denote this emission rate by N . We define a nondimensional emissionrate normalized to that value: ˜ N ∗ = N ∗ N . (20)To find characteristic density and length scales for the flow, consider first H cl , the scale heightof the neutral cloud. According to equation (7), this quantity depends on both the effective soundspeed a cl , which we fix at 2 km/s, and on the midplane cloud density n , which will be a freeparameter. A second length of importance is the Str¨omgren radius of fully ionized gas of uniformparticle number density n , given by R S ≡ (cid:18) N ∗ πα B ( n/ (cid:19) / . (21)Notice the appearance of n/ R S ∼ H cl . If R S (cid:28) H cl , the HII regionwould be trapped within the cloud and unable to generate the observed flow. If, on the otherhand, R S (cid:29) H cl , the ionized region would be free to expand in all directions, again contrary toobservation. For the purpose of defining a characteristic ionized density scale, we first set N ∗ = N in equation (21). We then set n = βn/ n thatsatisfies the relation H cl = R S . We label this density n , and find that it can be expressed as n = 9 π β (cid:18) G m H a (cid:19) (cid:18) N α B (cid:19) = 1 . × cm − . (22)We insert this value of n into equation (21) and set the resulting R S equal to the characteristiclength scale Z . We find for this length Z = 13 π β (cid:18) G m H a (cid:19) − (cid:18) N α B (cid:19) − = 0 .
18 pc . (23) 13 –It is encouraging that our values for n and Z match typical observations for UCHII regions(Churchwell 2002). We first normalize all lengths by Z :˜ r ≡ r/Z (24a)˜ z ≡ z/Z , (24b)and all densities by n : ˜ n ≡ n /n (25a)˜ n ≡ n/n . (25b)Since the ionized flow is transonic, we normalize all velocities to the ionized sound speed a I , whichwe fix at 10 km/s: ˜ u = u/a I . (26)We introduce a nondimensional expression for the force per mass due to radiation pressure:˜ f rad ≡ f rad (cid:18) a I Z (cid:19) − . (27)This nondimensional quantity is the radiative force relative to that from thermal pressure. Then,using equation (8) for f rad , we have˜ f rad = (cid:104) cos θ (cid:105) (cid:18) (cid:15)/cm H a I (cid:19) (cid:26) Z /a I [( n/ α B ] − (cid:27) . (28)The second factor on the right is the ratio of the momentum of an ionizing photon to the thermalmomentum of a gas particle. The third factor is the ratio of the sound crossing time in the flow tothe local photon recombination time.Since we are fixing the ionized sound speed, the only quantities that vary with z in equation(28) are n and (cid:104) cos θ (cid:105) . Thus, we are motivated to write˜ f rad = γ ˜ n (cid:104) cos θ (cid:105) , (29) Strictly speaking, (cid:15) also varies with position. Higher-energy photons have a lower photoionization cross section,and thus travel farther from the star before they are absorbed, leading to a gradual hardening of the radiation withincreasing distance. In our wavelength-independent analysis, we ignore this effect.
14 –where γ is also nondimensional. After setting the righthand sides of equations (28) and (29) equalto each other, and making use of equations (22), (23) and (25b), we find that γ can be expressedas γ ≡ (cid:18) L c (cid:19) (cid:18) a G (cid:19) − (cid:18) (cid:15)(cid:15) (cid:19) . (30)Here L ≡ N (cid:15) is the luminosity (in erg s − ) of a star of spectral type O7.5, and (cid:15) is themean energy of ionizing photons emitted from such a star. We show in Appendix A that (cid:15) , andhence γ , is a weak function of ˜ N ∗ . The function γ ( ˜ N ∗ ) is plotted in Figure 3. In summary thethree free parameters that we vary between calculations are ˜ n , ˜ N ∗ , and the star’s displacementfrom the midplane, ζ ≡ H ∗ /Z .We now summarize our nondimensional equations. After dropping the tilde notation for therest of this section, the decoupled equations of motion are dudz = (cid:18) − u (cid:19) − u R dRdz (1 + v ) + 2 v inj R (cid:115) (cid:18) dRdz (cid:19) (1 + u ) + u d Φ cl dz − u n γ (cid:104) cos θ (cid:105) (31) dndz = (cid:18) − u (cid:19) n R dRdz ( u + v ) − n R (cid:115) (cid:18) dRdz (cid:19) u v inj − n d Φ cl dz + n γ (cid:104) cos θ (cid:105) . (32)The expression for (cid:104) cos θ (cid:105) is still given by equation (9) if we use the appropriate nondimensionallengths. From equation (17), the injection velocity is v inj = (cid:115) n /β − n n cl n n cl − n , (33)where β is fixed at (10 / = 25.The nondimensional cloud density and gravitational potential are n cl = n sech (cid:20) ( z + ζ ) (cid:114) n β (cid:21) , (34)Φ cl = − β ln sech (cid:20) ( z + ζ ) (cid:114) n β (cid:21) . (35)Finally, the simplified ray tracing equation (4) becomes13 N ∗ = (cid:90) r f ( θ )0 n ( z ) r dr . (36) The shape of the ionization front can be obtained through equation (36), but only after thedensity n ( z ) is established. Since we do not know this density a priori, we begin with a guessed 15 –function. We then trace out the ionization front, and thus establish R ( z ). We next calculate both n ( z ) and u ( z ) by integrating the coupled equations (31) and (32). This procedure yields a newdensity distribution n ( z ) which we use to retrace the locus of the ionization front, again usingequation (36). The process is repeated until convergence is reached.To integrate equations (31) and (32), we must specify values of n and u at the base of theflow, where z = R = 0. These two initial values are not independent. We show in Appendix B that u ( − H b ) = v inj ( − H b ), where v inj ( − H b ) is the injection speed at the base, as found from equation(33). Since we know n cl ( − H b ) from equation (34), v inj ( − H b ) is solely a function of n ( − H b ), andthus u ( − H b ) is too. In practice, therefore, we need only guess n ( − H b ). Also derived in AppendixB are expressions for du/dz and dn/dz at the base, where the righthand sides of equations (31)and (32) have divergent terms.What is the correct value of n ( − H b )? The righthand sides of both equations (31) and (32) haveprefactors that diverge when u = 1. Thus, crossing the sonic point requires special care. For anarbitrarily guessed n ( − H b ), u ( z ) either diverges upward or declines toward zero as the sonic pointis approached. This behavior is a generic feature of wind problems, and a bifurcation procedureis often employed to pinpoint the physical flow. We use the method of “shooting and splitting”(Firnett 1974). Here, we repeatedly guess n ( − H b ), and successive densities farther downstream,until the velocity profile is established to within a preset tolerance. Specifically, iterations continueuntil the range of this accurate profile include u -values sufficiently close to 1, typically 0.98 or 0.99.To jump over the sonic point, we use the current values of du/dz and dn/dz to perform asingle first-order Euler integration step, typically with a z -increment of 0.01 – 0.03, depending onthe values of the derivatives. Once we are downstream from the sonic point, we revert to directintegration of equations (31) and (32).A key feature of the flows we generate is that the density near the base is low enough that theionized gas is underpressured with respect to the neutral cloud. The pressure drop causes neutralgas to be drawn into the flow and thereby replenish it. Driven by a combination of thermal andradiative forces, the flow accelerates. Thus, its density falls, but the decline is mitigated by thecontinual influx of fresh gas. On the other hand, the cloud density always falls sharply (see equation34). Eventually, n cl ( z ) reaches the value βn ( z ) /
4, at which point the ionized and neutral gas haveequal pressures. According to equation (33), no more neutral gas is drawn into the flow beyondthis point. In reality, the flow diverges laterally and its density also falls steeply. We do not followthis spreading process, but end each calculation at the point where the pressures cross over. 16 –
5. Results5.1. Fiducial model characteristics
For our fiducial model, we set the three nondimensional parameters to: ˜ n = β/
4; ˜ N ∗ = ζ = 1.The value of ˜ n is chosen so that the cloud midplane is in pressure balance when the flow density˜ n is unity. Dimensionally, the midplane density is 8 . × cm − , the star’s photon emission rateis N ∗ = 1 × s − , and the star is displaced from the midplane by H ∗ = 0 .
18 pc. Figure 4 showsthe converged shape of the ionization front for this case. The rapid flaring of the base essentiallyreproduces the typical cometary shapes observed. Of course, a more precise comparison is betweenthe predicted and observed emission measures. Here, too, the qualitative agreement is good, as wewill later demonstrate.The lowest dashed line in Figure 4 represents the cloud midplane. Note that the base of theionization front lies slightly below it. We have also displayed, as the middle dashed line, the sonictransition. In this particular model, the flow speed reaches a I close to the z -position of the star.This near match does not hold throughout most of parameter space. For example, lowering thestellar emission rate N ∗ moves the sonic transition farther above the star. Finally, the uppermostdashed line marks the height where the internal pressure of the flow overtakes that of the parentmolecular cloud. In this case, the crossover point is about the same distance from the star as thebase of the flow. As explained previously, we end our calculation at the crossover point, and do notattempt to track the complex dynamics of the flow as it continues to spread laterally and depositits momentum into the lower-density surrounding gas. In any event, there are fewer observationalconstraints on this more diffuse flow.Figure 5 shows the run of the density and velocity of the ionized gas with z -position. Thevelocity displays a smooth sonic transition, and a nearly constant acceleration throughout. Thevelocity at the base is very close to zero, but this is not the case for other parameter choices, aswe will describe in the following section. As the velocity rises, the density falls, creating a pressuregradient that works to accelerate the flow.To analyze more quantitatively the flow dynamics, it is helpful to gauge the relative contri-butions of the various terms contributing to the overall momentum balance. The nondimensionalversion of equation (14) is u dudz = − dndz − d Φ cl dz + f rad + 2 R dRdz v − R (cid:115) (cid:18) dRdz (cid:19) v inj u . (37)Recall that the fourth and fifth righthand terms represent, respectively, the ram pressure from theinjected gas and the retarding effect of mass loading.Figure 6 plots all five righthand terms as a function of z -position in the flow. The injected rampressure and mass loading dominate at first, and nearly balance one another. Soon, however, thethermal pressure gradient takes over and remains dominant thereafter. The radiation force at first 17 –retards the flow, and is a minor contributor to its acceleration above the star’s position. Finally, itis noteworthy that the gravitational force is negligible for most of the flow. The top panel of Figure 7 shows how the shape of the converged ionization front changes withthe stellar photon emission rate N ∗ , while the cloud midplane density n and stellar displacementfrom the midplane ζ are held fixed at their fiducial values. In this case the variation is predictable.A higher photon emission rate allows the ionizing photons to penetrate farther into the neutralcloud, in a manner similar to how the Str¨omgren radius of a spherical HII region increases withemission rate. A larger luminosity also leads to more flaring of the ionization front at large valuesof z .We did not consider systems with photon emission rate N ∗ < . n and ζ , we could not achieve transonic solutionswith N ∗ (cid:38)
1. As can be seen from the plot of velocity, u is already quite close to 0 for N ∗ = 1.Attempting to raise the photon emission rate any higher without changing the other parameterscreates such a large ionized density near the base that the flow is no longer underpressured withrespect to the cloud. It is possible to create flows with higher values of N ∗ if, for example, ζ isincreased at the same time.Generally, we find that the density profile of the ionized flow mimics that of the neutral cloud.Figure 8 shows how the density and velocity profiles vary when the stellar luminosity is changed.The flow structure remains strikingly similar as the luminosity is varied, changing primarily inspatial extent, because a higher luminosity allows ionizing radiation to penetrate farther into thecloud. The fact that the density and velocity profiles do not show other significant variations withluminosity reinforces the conclusion that it is the density structure of the neutral cloud that setsthe spatial variation in the ionized flow.The second panel of Figure 7 shows how the shape of the converged ionization front changeswith the stellar displacement ζ . Here, the effect is similar to that of increasing N ∗ . With higher ζ ,the ionizing photons penetrate a larger distance into the cloud, in this case because they encountera lower density when the star is displaced farther from the cloud midplane. Note also the flaringat large z , which increases sharply for higher ζ .Figure 9 shows how the density and velocity profiles vary when ζ is changed. Again, theionized gas density tracks that of the neutral cloud. For larger offsets, the flow begins in a lessdense portion of the cloud, and the density of the ionized flow is smaller. Since the ionizing stellarphotons penetrate farther into the cloud when ζ is larger, the flow begins at more negative valuesof z . The velocity profiles shift along the z -axis as ζ is varied, but the acceleration remains roughlythe same. 18 –Attempting to lower ζ below a value of unity, while keeping n and N ∗ fixed at their fiducialvalues, results in the same problem as attempting to increase N ∗ on its own. Solutions with lowervalues of ζ can be achieved only if n is simultaneously increased or N ∗ is decreased. Increasing ζ alone also leads to a large amount of flaring of the ionization front. For ζ = 2, the flaring at thelocation of pressure crossover is such that R/H b = 5 .
3. As we shall see in section 6.1, this aspectratio is large compared to observations.The third panel of Figure 7 shows how the shape of the ionization front changes when n isvaried on its own. In this case, the variation is more complex. As can be seen from equations (5)and (7), changing n changes not only the overall magnitude of the density profile, but also thesteepness of its falloff. When n is increased from β/ β/
2, the rise in the midplane densitymeans that radiation cannot penetrate as far, and the base of the ionized flow moves closer to thestar. However, as n is increased further to 3 β/
4, the steeper falloff of the density comes into play,and the distance between the star and the base of the flow remains nearly constant.We may effectively factor out the increasing cloud scale height if, instead of increasing n onits own, we simultaneously decrease ζ so that H ∗ /H cl remains equal to unity. This result of thisexercise is shown in the bottom panel of Figure 7. As n increases and ζ decreases, the shape of theflow remains strikingly similar, changing primarily in spatial extent. In all cases, H b / H cl remainsclose to unity, ranging from 1 .
21 when n = β/ .
94 when n = 3 β/ n . For a denserneutral cloud, the ionized flow also begins at a larger density. At the same time, the density ofboth the neutral and ionized gas drop more rapidly for a larger n , and the z -extent of the flowdiminishes. Finally, since the ionized gas is primarily accelerated via pressure gradients, a steeperdecline in density also corresponds to a more rapid acceleration.We may again factor out the effect of scale height by varying n and ζ simultaneously in themanner described previously. The resulting density and velocity profiles are shown in Figure 11.In this case the density at the base of the ionized flow again increases with rising n , but nowin a more systematic manner. In fact, this base density remains almost exactly equal to 4 n /β ,reflecting the fact that the flows are beginning at pressures very close to the neutral cloud pressure.This near-pressure equality also accounts for the fact that the velocities begin near zero in all ofthese flows. The velocity profiles are remarkable for the fact that they share an anchor point nearthe stellar position at z = 0. Their z -length scales correlate tightly with the changing cloud scaleheight. 19 – We next consider the total rate at which the neutral cloud adds mass to the ionized flow.Dimensionally, this rate, from the base to any height z , is˙ M ( z ) ≡ π (cid:90) z − H b ρ v inj R (cid:115) (cid:18) dRdz (cid:19) dz . (38)Figure 12 displays ˙ M ( z ) in our fiducial model. The cross-sectional area of the flow starts atzero and monotonically increases with z . The mass injection rate also climbs from zero, but levelsoff at the pressure crossover point, where no new mass is being added. We have also plotted theinjection speed, v inj ( z ). This starts out relatively small, peaks at about half the ionized soundspeed at z ≈ − .
5, and then eventually falls to zero. As the figure also shows, d ˙ M /dz attains itsmaximum close to where v inj peaks.We stop our calculation at the pressure crossover location z f , where neutral material ceases tobe drawn into the ionized flow. In steady state, the rate at which ionized mass crosses z f shouldequal the total rate of mass injection up to this point. That is, the mass outflow rate in ionizedgas should be ˙ M out = ˙ M ( z f ) (39)We may also calculate ˙ M out using the local dimensional relation˙ M out = π R f ρ f u f . (40)We find that we obtain the same mass outflow rate using these two methods to within a precisionof 0.1 percent, providing a useful check on mass conservation.Dimensionally, the mass outflow rate is˙ M out = 2 π (cid:16) m H (cid:17) n Z a I I f (41)= a Ga I I f (42)= 3 . × − M (cid:12) yr − I f , (43)where the nondimensional quantity I f is I f ≡ (cid:90) z f n v inj R (cid:115) (cid:18) dRdz (cid:19) dz . (44)The dependence of ˙ M on our three parameters N ∗ , n , and ζ is entirely contained within I f . 20 –We also define a momentum outflow rate ˙ p out ≡ ˙ M out u f . This quantity may be written˙ p out = 2 π (cid:16) m H (cid:17) n Z a I u f I f (45)= (cid:18) a G (cid:19) (cid:18) u f a I (cid:19) I f (46)= 2 . × dyne (cid:18) u f a I (cid:19) I f . (47)Finally, the dimensional kinetic energy outflow rate, ˙ E out ≡ / M out u f , is˙ E out = (cid:18) (cid:19) π (cid:16) m H (cid:17) n Z a I u f I f (48)= (cid:18) a a I G (cid:19) (cid:18) u f a I (cid:19) I f (49)= 1 . × erg s − (cid:18) u f a I (cid:19) I f . (50)Table B displays the values of ˙ M out , ˙ p out and ˙ E out for a variety of parameter choices. The massloss rate is far more sensitive to changing parameters than is u f , so the trends in momentum andenergy transport can be explained almost entirely by the trends in ˙ M out . Moreover, as equation(40) demonstrates, ˙ M out is principally affected by changes in the cross-sectional area and densityat z f .As before, we first examine the effect of varying the photon emission rate. With increasing N ∗ , R f increases as well, driving up ˙ M out . On the other hand, spreading of the flow results in aslight decline of the density n f . This latter effect weakens the sensitivity of ˙ M out to N ∗ . As thedimensional N ∗ rises from 10 to 10 s − , ˙ M out scales as N ∗ . .Similar trends appear when we increase the stellar offset ζ , leaving other parameters fixed (seesecond level of Table 1). With larger ζ , the ionization flow flares out. At the same time, the flowbegins in a less dense portion of the cloud, so that n f falls. Nevertheless, the net effect is for ˙ M out to increase. In the range 1 < ζ <
2, we find that ˙ M out scales as ζ . .Increasing n on its own leads to complicated behavior similar to that we encountered previ-ously. As n rises from β/ β/
2, ˙ M out first falls because the ionization front shrinks in size (seethird panel of Figure 7). However, as n further increases from β/ β , the cloud scale heightcontinues to shrink. The rapidly declining cloud density at any fixed z causes the ionized outflowto broaden, and ˙ M out rises.Finally, we may vary n and ζ simultaneously, so as to keep H ∗ /H cloud = 1. The fourth panelof Figure 7 shows that the ionization front retains its shape but shrinks in scale. As Figure 11shows, increasing n under these circumstances raises the ionized density at any z . Consequently,the stellar radiation cannot penetrate as far. In particular, R f becomes smaller at the pressure 21 –crossover point. Table 1 verifies that this shrinking of the cross-sectional area causes ˙ M out todecline. Quantitatively, we find ˙ M out ∝ n − . .The fact that the mass outflow rate generally decreases with increasing cloud density is signif-icant, and calls for a more basic physical explanation. The lateral size R of the outflow is aboutthat of a pressure-bounded HII region. From equation (21), we have R ∝ n − / I , where n I is arepresentative ionized density in the flow. From equation (40), the mass loss rate ˙ M out is propor-tional to R n I . Together these two relations imply ˙ M out ∝ n − / I . Finally, if n I is proportional tothe peak neutral density n , as follows from the condition of pressure equilibrium, then we have˙ M out ∝ n − / , which is close to our numerical result.This argument also helps to explain why the mass loss rates that we obtain are smaller thanthose of past champagne flow calculations. For example, Bodenheimer et al. (1979) found that,for a star with N ∗ = 7 . × s − embedded in a slablike molecular cloud of uniform numberdensity 10 cm − , the mass loss rate is 2 × − M (cid:12) yr − . For the much denser molecular cloud inour fiducial model, extrapolation of the Bodenheimer et al. (1979) rate using ˙ M out ∝ n − / yields˙ M out = 4 × − M (cid:12) yr − , quite close to our calculated result. It should be kept in mind that thiscomparison is only meant as a consistency check, given the many differences in detail of the twocalculations.We may combine our scaling results for the mass outflow rate in order to show explicitly itsdependence on N ∗ and n . We have˙ M out = 3 . × − M (cid:12) yr − (cid:18) N ∗ s − (cid:19) . (cid:16) n cm − (cid:17) − . , (51)More precisely, the two power-law indices are 0 . ± .
05 and − . ± . H ∗ varies with n so that the former equals onecloud scale height. We also caution that all these numerical results are based on a slab modelfor the parent cloud. Bearing these caveats in mind, equation (51) may prove useful for futuremodeling of massive star formation regions.The ionized wind in an UCHII region represents a large increase in mass loss over what thedriving star would achieve on its own. A bare, main-sequence O7.5 star with N ∗ = 1 × s − radiatively accelerates its own atmosphere to create a wind with 10 − to 10 − M (cid:12) yr − (Mokiemet al. 2007; Puls et al. 2008), two to three orders of magnitude below our fiducial value. However,one cannot completely ignore the dynamical effect of the stellar wind in an UCHII region, as wewill discuss shortly. Even younger stars of the same luminosity drive bipolar outflows that havefar greater mass loss rates than we compute, typically 10 − to 10 − M (cid:12) yr − (Churchwell 2002).While the exact mechanism behind these molecular outflows is uncertain, they might result fromthe entrainment of cloud gas by massive jets emanating from an accreting protostar (Beuther &Shepherd 2005). 22 –Over the inferred UCHII lifetime of 10 yr, a star driving an outflow at the rate of 4 × − M (cid:12) yr − disperses 40 M (cid:12) of cloud mass. This figure is tiny compared with 10 M (cid:12) , the typical mass of ahigh-density clump within an infrared dark cloud (Hofner et al. 2000; Hoare et al. 2007). Thus theUCHII region, which arises in the densest part of the cloud structure, represents an early stage inits clearing. Presumably, a compact HII region, of typical size 0.1 to 0.3 pc, appears as the starbegins to clear less dense material. We expect the outflow region to broaden and the mass loss rateto increase during this longer epoch.Our O7.5 star drives a wind with an associated momentum output of about L ∗ /c = 3 . × dyne(Vacca et al. 1996). This figure is remarkably close to the ˙ p out of 4 × dyne in our fiducial model.The numerical agreement is fortuitous, since the ionized outflow represents material drawn in fromthe external environment and accelerated by thermal pressure. The correspondence between thesetwo rates reflects the numerical coincidence that a cl /G ∼ L ∗ /c for an embedded O star. In anyevent, a more complete treatment of the flow in an UCHII region would also account for the addi-tional forcing from the stellar wind (see, e.g., Arthur & Hoare 2006 for a simulation that includesthis effect). Note finally that the ionizing photons themselves impart momentum to the flow. Theresulting force is f rad , which we introduced in Section 2.4, and whose quantitative effect we discussbelow.The kinetic energy transport rate for our fiducial model is only 4 × − times L ∗ , the bolometricluminosity of the star. Given that the flow is transonic, the thermal energy carried in the outflow isof comparable magnitude. The vast bulk of the stellar energy is lost to radiation that escapes fromthe HII region during recombination, through line emission from ionized metals, and continuum,free-free emission (Osterbrock 1989). For a time span of 10 yr, the total energy ejected in ourfiducial model is a few times 10 erg, a figure comparable to the turbulent kinetic energy in amolecular clump of mass 10 M (cid:12) and internal velocity dispersion of 2 km s − . This match broadlysupports the contention of Matzner (2002) that HII regions provide the ultimate energy source ofturbulence in molecular clouds large enough to spawn massive stars. One of the novel features in this analysis of UCHII regions is our inclusion of the momentumdeposition by ionizing photons. We derived f rad , the radiative force per mass, in Section 2.4 (seeeq. (8)), and expressed it nondimensionally in equation (29). How significant is this term in theoverall flow dynamics?To assess the role of f rad , we recalculated our fiducial model after artificially removing theforce from the equations of motion, (31) and (32). Figure 13 shows the result of this exercise. Thedashed curve in the top panel is the altered density profile, n ( z ), while the analogous curve in thebottom panel is the altered velocity, u ( z ).Near the base, the stellar flux is directed oppositely to the ionized gas velocity. Thus, the 23 –radiation force decelerates the flow in this region. The bottom panel of the figure shows how thestarting velocities are lower when radiation pressure is included. As a result of this force, gas pilesup near the base, and therefore develops a larger thermal pressure gradient. The enhanced gradientdrives the flow outward in spite of the retarding stellar photons. The density pileup is evident inthe solid curve within the top panel of Figure 13.The effect of radiation pressure diminishes at higher values of z . Recall that f rad is proportionalto the flow density n . Gas of higher density has a greater volumetric recombination rate, and thusabsorbs more ionizing photons per time to maintain ionization balance. Since the density falls withdistance z , so does the magnitude of the force.A second contributing factor is the changing direction of photons emanating from the star.Near the base, the incoming flux is nearly all in the negative z -direction, so that (cid:104) cos θ (cid:105) is close to −
1. This geometric term vanishes at the level of the star ( z = 0), and then climbs back up toward+1. However, the rapidly falling density overwhelms the latter effect, and the magnitude of theforce still declines (see the dotted curve in Figure 6).Finally, the radiation force does not significantly alter the shape of the ionization front. Inour fiducial model, the base of the flow is farther from the star when f rad is omitted. However, thefractional change in this distance from the case with the force included is only 0.05. The outflow topology sketched in Figure 1 does not exist for arbitrary values of our parameters.If the star is situated too close to the densest portion of the cloud, or if the ionizing luminosityis too weak, then a transonic flow, steadily drawing in neutral material, cannot develop. At asufficiently low luminosity or high neutral density, the HII region first undergoes pressure-drivenexpansion, but then remains trapped, i.e., density-bounded, on all sides.Alternatively, with a much higher luminosity and/or lower cloud density, an outflow may erupton both sides of our model planar slab. Such a bipolar morphology is harder to achieve, in partbecause of the values of cloud density and stellar luminosity required, and in part because realclouds do not have a slab geometry, i.e., they are not infinite in lateral extent. Increasing thestellar luminosity in an originally monopolar flow usually just widens the ionization front, withoutcreating another outflow lobe. It is therefore not surprising that out of the hundreds of UCHIIregions that have been identified, only a handful have a bipolar morphology (Garay & Lizano 1999;Churchwell 2002).Given their relative scarcity, we forgo a parameter study of this type of UCHII region, andfocus instead on the simplest example. Here we place the star exactly at the cloud midplane, ζ = 0.We set the cloud midplane density to n = β/
4, as in our fiducial, monopolar model. However, wefind that N ∗ = 1 does not lead to a bipolar flow. To be safe, we have raised N ∗ to 3. 24 –In order to generate this solution, we began with the fact that u (0) = 0, as demanded bysymmetry. We guessed n (0), the flow density at the midplane, and then used these two startingvalues as the initial conditions for integrating the coupled first-order equations (31) and (32). Asbefore, we use equation (33) to obtain the injection velocity. Since R does not vanish at the base,the righthand sides of (31) and (32) are well-behaved at the start, and the integration is relativelystraightforward. Using the method of shooting and splitting, we refine our guess for the startingionized density until we approach the sonic point. We jump over this point in just the mannerdescribed in Section 4.3.Figure 14 shows the symmetric ionization front for this model. The two sets of horizontallines show the location of the sonic transitions and the pressure crossover points, respectively. Wenotice immediately how close these points are to one another (compare Fig. 4). This same featureis apparent in Figure 15, which displays the density and velocity profiles. The velocity, which isan odd function of the height z , reaches unity just before the curve ends at the pressure crossover.The ionized density n ( z ) is symmetric about z = 0 and has a shape similar to that of the neutralcloud, peaking at the midplane.In this outflow, neutral gas with a pressure exceeding that of the ionized gas is drawn inlaterally through the ionization front. The ionized gas flows away in a symmetric manner from itsregion of maximal density at the midplane. Radiation pressure never acts to decelerate the flow,as in the monopolar case. Rather, it contributes to the acceleration, consequently reducing theionized density gradient. Nevertheless, the thermal pressure gradient remains the strongest drivingforce.As mentioned previously, attempts to lower N ∗ closer to 1 resulted in failure to obtain atransonic solution. Specifically, the pressure crossover point was reached before the sonic transition.We found this result for any guess of the starting ionized density below 4 n /β , i.e. for ionizedpressures at the midplane less than the corresponding cloud pressure. Our inability to find transonicsolutions in this regime suggests that the true steady-state solution is not an outflow, but a trappedHII region in which the interior and cloud pressures match. Up to this point we have taken our molecular cloud to be a planar slab, a geometry that isconvenient mathematically, but not truly representative of the clouds found in nature. There existno detailed studies of infrared dark cloud morphologies. In the realm of low-mass star formation,recent observations of Gould Belt clouds reveal that a tangled network of filaments creates thedense cores that in turn collapse to form stars (Andr´e et al. 2010). Thus, a more realistic modelfor our background cloud might be a cylinder in force balance between self-gravity and turbulentpressure. If the interior velocity dispersion is again spatially uniform, as we have assumed, then thedensity falls with a power law with distance from the central axis, as opposed to the much steeper, 25 –exponential falloff we have employed until now.Our model is flexible enough that we can explore various functional forms for the cloud densityprofile. However, our cloud is still a one-dimensional slab, so that the density ρ is a function of z ,the distance from the midplane. We construct a pseudo-cylindrical cloud, which has the densityprofile of a self-gravitating, isothermal cylinder (Ostriker 1964), but with the cylindrical radius R replaced by z . In this model, the dimensional number density and gravitational potential are n cl = n (cid:34) (cid:18) z + H ∗ H cl (cid:19) (cid:35) − , (52)Φ cl = − a ln (cid:34) (cid:18) z + H ∗ H cl (cid:19) (cid:35) . (53)Here H ∗ is again the distance between the star and the cloud midplane, and H cl is the scale heightin equation (7) .Choosing N ∗ = 1 and ζ = 1, we could not achieve a transonic flow, if we also used the fiducial n = β/
4. In retrospect, this result could have been anticipated, since the pseudo-cylindrical modelhas a smaller column density, as measured from the midplane, than the slab model with the same n . We therefore utilized n = β/ n = πβ/
2, which has thesame column density as the pseudo-cylinder. The two density profiles, both non-dimensional, aredisplayed together in Figure 16.Figure 17 compares the ionization front shapes for the two cloud models. Because of itsrelatively high n -value, the ionization front in this slab model is more flared than in the fiducialone (recall Fig. 7). The degree of flaring in the pseudo-cylindrical case is much less, a consequenceof the gentler falloff in the cloud density profile.Finally, Figure 18 compares the density and velocity profiles within the flows themselves. Theflow density within the pseudo-cylinder starts out lower and falls off more slowly. Since the thermalpressure gradient is reduced, so is the acceleration of the ionized gas. As seen in the lower panel ofthe figure, the velocity begins at a more subsonic value and thereafter climbs less steeply.
6. Comparison to Observations6.1. Emission measure maps
One way to compare our model with observations is to generate synthetic contour maps of theradio continuum emission measure. The latter is the integral of n with respect to distance along Our choice of formula for H cl is a factor of √
26 –each line of sight that penetrates the outflow. The resulting maps may also be compared to thosegenerated by other theoretical models, such as the ones of Redman et al. (1998; Fig. 3) and Arthur& Hoare (2006; Fig. 7),Figure 19 displays a series of emission measure maps at different viewing angles, all using ourfiducial model. Here, the inclination angle θ is that between the flow’s central ( z -) axis and the lineof sight. The star, as always, lies at the origin, and the spatial coordinates are non-dimensional.For a flow oriented in the plane of the sky ( θ = π/ z , while the n decreases. In practice, the peak emission is located about midway between theflow base and the star. The latter two points are separated by a physical distance of about Z / . × cm − pc. When N ∗ is loweredby a factor of 10, as in Table 1, the peak emission measure drops to 1 . × cm − pc. Suppose, onthe other hand, we fix N ∗ at its fiducial value while increasing n and concurrently decreasing ζ inthe manner described below equation (51). Then we find that the peak emission measure reaches1 . × cm − pc for an n of 3 β/ .
75. Of the 15 cometaries in Wood & Churchwell (1989)with estimated peak emission measures, the figure varies from 2 × to 3 × cm − pc. Of the 12cometaries in Kurtz et al. (1994) with estimated peak emission measures, they vary from 1 × to6 × cm − pc. Thus our model, including reasonable parameter variations, yields peak emissionmeasures that fall within the middle range of observed values.As the z -axis tips toward the line of sight, the emission becomes fainter. The reason is thatthe total emission measure along any line of sight is increasingly weighted by more rarefied gasat higher values of z . Thus, as seen in Figure 19, the contours representing higher values of theemission measure successively disappear as θ decreases.For θ close to zero, the shape of the remaining contours is nearly spherical. However, our syn-thetic maps fail to reproduce the shell or core-halo morphologies identified by Wood & Churchwell(1989). The latter have pronounced limb brightening. As previously noted by Zhu et al. (2005),champagne-flow models have difficulty explaining this feature.Returning to generic cometary regions, our model fits their morphology quite well. Figure 20compares our fiducial flow (with θ = π/
2) to an emission map taken from the survey of Wood &Churchwell (1989). This specific UCHII region has an estimated distance of 16.1 kpc. Accordingly,our synthetic map is displayed in angular coordinates. We could, in principle, make more precisefits to both this region and many others, by tuning our three free parameters and the viewing angle.Here, we do not attempt such a detailed, comprehensive matching. 27 –
As explained at the outset of Section 2.3, we calculate a laterally-averaged ionized gas velocityat each z -location within our model HII region. For this reason, we cannot attempt to matchdetailed observations of UCHII velocities along particular lines of sight. Nevertheless, we do re-produce several broad characteristics of some UCHII regions. Our model yields a transonic flowthat proceeds along the cometary axis, with steep acceleration from head to tail. A similar patternis seen in the objects observed by Garay et al. (1994). For our fiducial cometary UCHII region,the velocity gradient we compute is 8 km s − pc − , about a factor of 2 lower than the gradientmeasured in G32.80 ± − pc − , about one third the gradientmeasured in the compact bipolar HII region K3-50A (Depree et al. 1994).Like all champagne-flow models, ours predicts that the ionized and molecular gas should beclosest in velocity near the head of the cometary region, and that the two velocities should growfarther apart down the tail of the ionized region. This is not the case for G13.87 ± .
28, for whichthe ionized velocity closely approaches the molecular cloud velocity in the tail, a characteristic ofbowshock models (Garay et al. 1994). However, in G32.80 + .19, the central line emission of theionized gas in the tail shows accelerates to a speed at least 6 km s − greater than that of themolecular cloud (Garay et al. 1994). This systematic increase is more consistent with a champagneflow.For completeness, we note that there exist a number of cometary UCHII regions for whichneither a champagne flow nor bowshock accurately accounts for the detailed velocity structure.One well-studied example is G29.96-02 (Mart´ın-Hern´andez et al. 2003). Objects such as these havemotivated numerical simulations in which a champagne flow is combined with a stellar wind andmotion of the star itself through a molecular cloud (Arthur & Hoare 2006).In summary, a pure champagne flow model such as ours is inadequate to explain all the observedfeatures of cometary UCHII regions. We are encouraged, however, that the main characteristics arewell produced. In addition, the simplicity of our quasi-one dimensional model in comparison withmultidimensional numerical simulations, allows a much broader exploration of parameter space thatwill aid in a general understanding of the phenomenon. 28 –
7. Summary and Discussion
We have explored a quasi-one dimensional, steady-state wind model for UCHII regions. Thispicture offers a natural explanation for the cometary morphology that is frequently observed. Weare also able to match, at least broadly, such basic observational quantities as the diameter of theregion, its peak emission measure, and the ionized gas velocity. Additionally, we have employed themodel to see how the density and velocity structure of the region respond to changes in the stellarluminosity, molecular cloud density, and displacement of the exciting star from the peak moleculardensity. Checking these trends against observations is a task for the future.In our view, the very young, massive star creating the ionization is still buried deep withinwithin a large molecular cloud. Ionizing photons from the star steadily erode the cloud, expellinggas at a rate of order 10 − M (cid:12) yr − . Consequently, over the inferred UCHII lifetimes of order10 yr, the star drives out only a relatively small amount of cloud gas, representing the denseclump in which it was born. As time goes on, we expect the flow to grow in size, with the massexpulsion rate concurrently rising.The ionized gas just inside the ionization front is underpressured with respect to the surround-ing neutral material. Indeed, it is just this pressure difference that draws in more neutral materialand sustains the wind. Roshi et al. (2005) observe just such a difference in the UCHII regionG35.20-1.74. In our model, this pressure difference is not due to the generation of a shock by theexpanding HII region as in the classical model for D-type fronts (Spitzer 1978, Chapter 12). Sucha shock probably did form earlier in the evolution of the UCHII region, but has died away by thetime a steady flow is established.Our semi-analytic model necessarily entails simplifying assumptions. One is our neglect of thestellar wind, which could partially excavate the ionized region. In over half the objects studied inZhu et al. (2008), the ionized gas appears to be skirting around a central cavity, possibly indicatingthe influence of the wind. Note also that a wind would produce an additional shock near the front,elevating the pressure difference.We have also assumed a steady-state flow in which we neglect the relatively slow advance of theionization front into the molecular cloud. Recall that equation (16) gives v (cid:48) cl , the incoming velocityof cloud gas in the frame of the ionization front. Using a typical injection speed of v (cid:48) inj = (1 / a I ,along with n = n and n cl = ( β/ n , we find that v (cid:48) inj = a I / β = 0 . − , which is indeedsmall compared to the effective cloud sound speed of 2 km s − . However, over the 10 yr lifetime ofthe UCHII region, the front advances by 0.02 pc, an appreciable fraction of the region’s diameter.A more complete calculation would account for the motion of this front, through retention of theflux leakage term in equation (2).In the future, both our type of semi-analytic model and multidimensional numerical simulationscan play valuable roles. Extending the present calculation, one could follow the spread of ionizationto see how the HII region grows and disperses the bulk of the molecular cloud. Simulations can 29 –eventually provide a more realistic picture of the parent cloud. In particular, we see a need toexplore the role of turbulence, which we have treated simplistically as providing an enhanced,isotropic pressure. As Arthur et al. (2011) have shown, the advance of an HII region into aturbulent cloud is quite different from the classical account. Within our model, the wind itself maynot be a laminar flow, as we have assumed. Inclusion of these effects may be necessary to explainthe more extreme members of the irregular UCHII class identified by Wood & Churchwell (1989). Acknowledgements
Throughout the work, NR was supported by the Department of Energy Office of ScienceGraduate Fellowship Program (DOE SCGF), made possible in part by the American Recoveryand Reinvestment Act of 2009, administered by ORISE-ORAU under contract no. DE-AC05-06OR23100. SS was partially supported by NSF Grant 0908573.
REFERENCES
Afflerbach, A., Churchwell, E., Acord, J. M., et al. 1996, ApJS, 106, 423Andr´e, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102Arthur, S. J., Henney, W. J., Mellema, G., de Colle, F., & V´azquez-Semadeni, E. 2011, MNRAS,414, 1747Arthur, S. J., & Hoare, M. G. 2006, ApJS, 165, 283Beuther, H., Churchwell, E. B., McKee, C. F., & Tan, J. C. 2007, Protostars and Planets V, 165Beuther, H., & Shepherd, D. 2005, in Cores to Clusters: Star Formation with Next GenerationTelescopes, ed. M. S. N. Kumar, M. Tafalla, & P. Caselli, 105–119Bodenheimer, P., Tenorio-Tagle, G., & Yorke, H. W. 1979, ApJ, 233, 85Churchwell, E. 2002, ARA&A, 40, 27Comeron, F. 1997, A&A, 326, 1195Depree, C. G., Goss, W. M., Palmer, P., & Rubin, R. H. 1994, ApJ, 428, 670Dyson, J. E. 1968, Ap&SS, 1, 388F˝ur´esz, G., Hartmann, L. W., Megeath, S. T., Szentgyorgyi, A. H., & Hamden, E. T. 2008, ApJ,676, 1109Garay, G., & Lizano, S. 1999, PASP, 111, 1049 30 –Garay, G., Lizano, S., & Gomez, Y. 1994, ApJ, 429, 268Hansen, C. J., Kawaler, S. D., & Trimble, V. 2004, Stellar interiors : physical principles, structure,and evolutionHenney, W. J. 2001, in Revista Mexicana de Astronomia y Astrofisica Conference Series, Vol. 10,Revista Mexicana de Astronomia y Astrofisica Conference Series, ed. J. Cant´o & L. F.Rodr´ıguez, 57–60Henney, W. J., Arthur, S. J., Williams, R. J. R., & Ferland, G. J. 2005, ApJ, 621, 328Hoare, M. G., Kurtz, S. E., Lizano, S., Keto, E., & Hofner, P. 2007, Protostars and Planets V, 181Hofner, P., Wyrowski, F., Walmsley, C. M., & Churchwell, E. 2000, ApJ, 536, 393Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654Kahn, F. D. 1954, Bull. Astron. Inst. Netherlands, 12, 187Keto, E. 2002a, ApJ, 568, 754—. 2002b, ApJ, 580, 980—. 2007, ApJ, 666, 976Keto, E., & Wood, K. 2006, ApJ, 637, 850Kim, K.-T., & Koo, B.-C. 2001, ApJ, 549, 979Krumholz, M. R., & Matzner, C. D. 2009, ApJ, 703, 1352Kurtz, S., Churchwell, E., & Wood, D. O. S. 1994, ApJS, 91, 659Lizano, S., Canto, J., Garay, G., & Hollenbach, D. 1996, ApJ, 468, 739Lugo, J., Lizano, S., & Garay, G. 2004, ApJ, 614, 807Mac Low, M.-M., Toraskar, J., Oishi, J. S., & Abel, T. 2007, ApJ, 668, 980Mac Low, M.-M., van Buren, D., Wood, D. O. S., & Churchwell, E. 1991, ApJ, 369, 395Mart´ın-Hern´andez, N. L., Bik, A., Kaper, L., Tielens, A. G. G. M., & Hanson, M. M. 2003, A&A,405, 175Matzner, C. D. 2002, ApJ, 566, 302McKee, C. F., & Tan, J. C. 2008, ApJ, 681, 771Mokiem, M. R., de Koter, A., Evans, C. J., et al. 2007, A&A, 465, 1003 31 –Mueller, K. E., Shirley, Y. L., Evans, II, N. J., & Jacobson, H. R. 2002, ApJS, 143, 469Oort, J. H., & Spitzer, Jr., L. 1955, ApJ, 121, 6Osterbrock, D. E. 1989, Astrophysics of gaseous nebulae and active galactic nucleiOstriker, J. 1964, ApJ, 140, 1056Peters, T., Banerjee, R., Klessen, R. S., & Mac Low, M.-M. 2011, ApJ, 729, 72Peters, T., Banerjee, R., Klessen, R. S., et al. 2010a, ApJ, 711, 1017Peters, T., Mac Low, M.-M., Banerjee, R., Klessen, R. S., & Dullemond, C. P. 2010b, ApJ, 719,831Puls, J., Vink, J. S., & Najarro, F. 2008, A&A Rev., 16, 209Redman, M. P., Williams, R. J. R., & Dyson, J. E. 1998, MNRAS, 298, 33Roshi, D. A., Goss, W. M., Anantharamaiah, K. R., & Jeyakumar, S. 2005, ApJ, 626, 253Spitzer, L. 1978, Physical processes in the interstellar mediumTenorio-Tagle, G. 1979, A&A, 71, 59Thompson, M. A., Hatchell, J., Walsh, A. J., MacDonald, G. H., & Millar, T. J. 2006, A&A, 453,1003Vacca, W. D., Garmany, C. D., & Shull, J. M. 1996, ApJ, 460, 914van Buren, D., Mac Low, M.-M., Wood, D. O. S., & Churchwell, E. 1990, ApJ, 353, 570Walsh, A. J., Burton, M. G., Hyland, A. R., & Robinson, G. 1998, MNRAS, 301, 640Whitworth, A. 1979, MNRAS, 186, 59Williams, R. J. R., Dyson, J. E., & Redman, M. P. 1996, MNRAS, 280, 667Wood, D. O. S., & Churchwell, E. 1989, ApJS, 69, 831Xie, T., Mundy, L. G., Vogel, S. N., & Hofner, P. 1996, ApJ, 473, L131Zhu, Q.-F., Lacy, J. H., Jaffe, D. T., Richter, M. J., & Greathouse, T. K. 2005, ApJ, 631, 381—. 2008, ApJS, 177, 584Zuckerman, B. 1973, ApJ, 183, 863
This preprint was prepared with the AAS L A TEX macros v5.2.
32 –
A. Relation between (cid:15) and N ∗ Consider a massive star of bolometric luminosity L ∗ , radius R ∗ , and temperature T ∗ . Assumingthe star emits like a blackbody, the mean energy of its ionizing photons is (cid:15) = (cid:82) ∞ ν crit B ν dν (cid:82) ∞ ν crit B ν / ( hν ) dν = k T ∗ (cid:82) ∞ x crit x / ( e x − dx (cid:82) ∞ x crit x / ( e x − dx , (A1)where x ≡ hν/ ( kT ∗ ). Here, ν crit denotes the Lyman continuum frequency, and x crit ≡ hν crit / ( kT ∗ ).We also may write N ∗ = 4 π R ∗ (cid:90) ∞ ν crit B ν hν dν , (A2)and so N ∗ ∝ R ∗ T ∗ (cid:90) ∞ x crit x e x − dx . (A3)We next use the standard scaling relations R ∗ ∝ M . ∗ , L ∗ ∝ M . ∗ for massive main sequence stars(e.g. Hansen et al. 2004) as well as the blackbody relation, L ∗ ∝ R ∗ T ∗ . We find N ∗ N = (cid:18) T ∗ . × K (cid:19) (cid:90) ∞ x crit x e x − dx . (A4)Here, we have implicitly assumed that the bulk of the star’s luminosity is in ionizing photons. Wehave used results of Vacca et al. (1996) to set the proportionality constant in the last relation. Forany T ∗ , equation (A4) then gives us N ∗ , while equation (A1) gives (cid:15) . We thus obtain the functionaldependence of (cid:15) on N ∗ . Finally, equation (30) in the text gives γ as a function of ˜ N ∗ , as plotted inFigure 3.As a check on the accuracy of our prescription, Figure 21 shows N ∗ as a function of T ∗ , bothfrom our analysis and from the more detailed calculations of Vacca et al. (1996). Our approximationis most accurate at the highest luminosities, but tends to overestimate the ionizing photon emissionrate for the lowest luminosities we consider. B. Boundary conditions at base of flow
We expand the derivative on the lefthand side of equation (12), and factor out dR/dz fromthe radical. We thus find2
R dRdz n u + R dndz u + R n dudz = 2 R dRdz (cid:115) (cid:18) dzdR (cid:19) n v inj . (B1) 33 –At the base, both R and dz/dR tend to zero. Thus, the second and third terms on the lefthandside of equation (B1) vanish, as does the square root on the righthand side. We therefore have u = v inj , (B2)exactly at the base. Recall that v inj is a function of n . We choose n at the base such that we obtaina transonic solution via the method of shooting and splitting.We also need to know the values of dn/dz and du/dz at the base. We begin by rewriting thedecoupled, nondimensional equations of motion as(1 − u ) dudz − u d Φ cl dz + uf rad = 2 (cid:113) dz/dR ) (1 + u ) v inj − v ) uR dz/dR (B3)(1 − u ) dndz + n d Φ cl dz − nf rad = 2( u + v ) n − n (cid:113) dz/dR ) u v inj R dz/dR , (B4)where f rad is given by equation (8). In the limits R → dz/dR →
0, the righthand sides ofboth equations can be evaluated using L’Hˆopital’s rule. In order to compute quantities such as d z/dR , we use the fact that, to lowest order in R , the cavity wall is parabolic, i.e. z ≈ A R − H b where A is a constant. At each timestep, we perform a numerical fit to determine A . The finalresults for the derivatives at the base are dndz = (cid:20) − u (cid:21) (cid:20) − n d Φ cl dz + n f rad − A n u (cid:21) (B5)and dudz = (cid:20) − u (cid:21) (cid:20) u d Φ cl dz − u f rad + A u (1 + u ) + 12 dv inj dz (1 − u ) (cid:21) . (B6)We obtain the term dv inj /dz in equation (B6) by differentiating equation (33): dv inj dz = 2( β n − n n cl + 16 n )( n dn cl /dz − n cl dn/dz ) β n ( n − n cl ) v inj . (B7) 34 –Table 1: Mass, momentum and energy transport rates N ∗ ζ n ˙ M out ˙ p out ˙ E out ( M (cid:12) yr − ) (dyne) (erg s − )0.1 1.0 β/ . × − . × . × β/ . × − . × . × β/ . × − . × . × β/ . × − . × . × β/ . × − . × . × β/ . × − . × . × β/ . × − . × . × β . × − . × . × (cid:112) / β/ . × − . × . × (cid:112) / β/ . × − . × . × (cid:112) / β/ . × − . × . ×
35 –Fig. 1.— Our flow schematic. A massive star is situated to one side of an initially neutral cloud.It creates, as in the original champagne model, an asymmetric HII region of relatively low density.Toward the midplane, the ionization front is stalled by rising cloud density. In the opposite direc-tion, the front has broken free. The high pressure of ionized gas creates an accelerating flow awayfrom the densest gas. At the base of the ionized flow, the pressure has been relieved sufficientlysuch that the neutral cloud gas is slightly overpressured with respect to the ionized gas. 36 –Fig. 2.— Control volume diagram for a segment of the ionized flow. The trapezoidal sectiondisplayed here represents a meridional slice of the ionised flow - one should imagine rotating thisdiagram about the central vertical symmetry axis in order to generate the represented volume. 37 –Fig. 3.— The variation of the radiative force coefficient γ defined in equation (30) with ˜ N ∗ , thestar’s ionizing photon emission rate, normalized to 10 s − . 38 –Fig. 4.— Converged ionization front shape for our fiducial UCHII region, with the star’s positionindicated. In this and succeeding figures, all physical variables are displayed nondimensionally.The lowest dashed line represent the z -position of the cloud midplane. The middle dashed linerepresents the z -position of the sonic point. Finally, the upper dashed line is the endpoint of oursolution, where the pressures of the ionized and neutral gas cross over. 39 –Fig. 5.— Ionized gas density and velocity as a function of z -position for our fiducial model. 40 –Fig. 6.— The magnitude of each of the terms in the momentum equation (37) as a function of z -position. The thermal pressure gradient, along with the ram pressure from injection of ionizedgas, act to accelerate the flow, while mass loading and the gravity of the parent cloud decelerate it.Radiation pressure mostly acts to decelerate the flow, but for z > β is set to 25. 42 –Fig. 8.— Density and velocity profiles for various ionizing photon emission rates. The stellardisplacement ζ and the midplane cloud density n are held fixed at their fiducial values of 1 and β/
4, respectively. 43 –Fig. 9.— Density and velocity profiles for various values of the stellar displacement ζ . The photonemission rate N ∗ and the midplane density n are held fixed at their fiducial values of 1 and β/ n . Thenondimensional ionizing photon emission rate N ∗ and stellar offset from the midplane ζ are heldfixed at their fiducial values of 1. 45 –Fig. 11.— Density and velocity profiles for various flow models. The midplane density n andstellar displacement ζ are varied simultaneously so as to keep the ratio H ∗ / H cl = 1. 46 –Fig. 12.— Nondimensional mass injection rate and related quantities as functions of z . The massinjection rate ˙ M ( z ) has been normalized to ˙ M = 3 . × − M (cid:12) yr − (see equation 43). Itsderivative d ˙ M /dz is normalized to ˙ M /Z = 2 × − M (cid:12) yr − pc − . Finally, the injection speed v inj has been normalized to the ionized sound speed, assumed to be 10 km s − . 47 –Fig. 13.— Density and velocity profiles both for the fiducial model ( solid curves ) and with radiationpressure omitted ( dashed curves ). 48 –Fig. 14.— The ionization front for a symmetric, bipolar transonic outflow. Here N = 3, n = β/ ζ = 0. The two horizontal dashed lines closer to the star mark the sonic transitions, while theouter pair indicates the pressure crossover. 49 –Fig. 15.— Density and velocity profiles for the bipolar outflow shown in Figure 14. 50 –Fig. 16.— Cloud density profiles for a slab and pseudo-cylindrical model of equal column density. 51 –Fig. 17.— Ionization fronts for the slab and pseudo-cylindrical cloud models shown in Figure 16.The horizontal dashed line corresponds to the cloud midplane for both models. 52 –Fig. 18.— Density and velocity profiles for the slab and pseudo-cylindrical cloud models shown inFigure 16 53 –Fig. 19.— Emission measure contours for our fiducial UCHII model viewed at various angles θ with respect to the flow’s central z -axis. The peak emission measure is 4 . × cm − pc, and thecontours in the top panel represent 0.95, 0.85, ... 0.15 times that value. A grey scale has beenapplied consistently to all the maps to indicate which emission measures correspond to the plottedcontours. 54 –Fig. 20.— Contours of radio continuum emission. The left panel uses our fiducial theoreticalmodel, with inclination angle θ = π/ . × pc cm − , whilethe inferred peak emission measure for G12.21-0.10 is 5 . × pc cm −6