Outflow-Confined HII Regions

# Outflow-Confined H ii Regions. I. First Signposts of Massive Star Formation

Kei E. I. Tanaka Department of Astronomy, University of Florida, Gainesville, FL 32611, USA; ktanaka@ufl.edu Jonathan C. Tan Departments of Astronomy & Physics, University of Florida, Gainesville, FL 32611, USA; jt@astro.ufl.edu Yichen Zhang Departamento de Astronomía, Universidad de Chile, Santiago, Chile; yzhang@das.uchile.cl
###### Abstract

We present an evolutionary sequence of models of the photoionized disk-wind outflow around forming massive stars based on the Core Accretion model. The outflow is expected to be the first structure to be ionized by the protostar and can confine the expansion of the H ii region, especially in lateral directions in the plane of the accretion disk. The ionizing luminosity increases as Kelvin-Helmholz contraction proceeds, and the H ii region is formed when the stellar mass reaches depending on the initial cloud core properties. Although some part of outer disk surface remains neutral due to shielding by the inner disk and the disk wind, almost the whole of the outflow is ionized in after initial H ii region formation. Having calculated the extent and temperature structure of the H ii region within the immediate protostellar environment, we then make predictions for the strength of its free-free continuum and recombination line emission. The free-free radio emission from the ionized outflow has a flux density of for a source at a distance of 1 kpc with a spectral index , and the apparent size is typically 500 AU at 10 GHz. The line profile has a width of about . These properties of our model are consistent with observed radio winds and jets around forming massive protostars.

###### Subject headings:
stars: formation, evolution - accretion, accretion disks.

## 1. Introduction

The formation mechanism of massive stars remains much debated with Core Accretion, Competitive Accretion and even Protostellar Collisions still discussed (see Tan et al., 2014, for a recent review).

One of the key differences between high and low-mass protostars is that the former are expected to complete the later stages of their accretion while fusing hydrogen on the main sequence. Their photospheres should be hot and strong sources of Lyman continnum photons with energies  eV, which can photoionize H to create an H ii region. These protostellar H ii regions are important for at least two main reasons. First, by their long wavelength radio free-free and Hydrogen recombination line (HRL) emission they provide diagnostics of the massive star formation process. Second, the ionized gas reaches temperatures of  K and its thermal pressure could be an important feedback process that regulates massive star formation, either by disrupting the wider scale accretion reservoir (i.e., core envelope or Bondi-Hoyle competitive accretion region: e.g., Dale et al., 2005; Peters et al., 2010, 2011) or by photoevaporation of the accretion disk (Hollenbach et al., 1994; Yorke & Welz, 1996; Richling & Yorke, 1997; McKee & Tan, 2008; Tanaka et al., 2013).

Our goal in this study is to present model calculations of the structure of the earliest-stage protostellar H ii regions in the context of the Turbulent Core Model of McKee & Tan (2003) (hereafter, MT03). Early work in this direction was carried out by Tan & McKee (2003), who described these structures as “Outflow-Confined” H ii regions, since the magneto-centrifugally-launched protostellar outflow (either disk-wind, X-wind or both) will be the first barrier to the ionizing photons.

Keto (2007) studied the structure and evolution of H ii regions that are confined by the infalling flow onto the protostellar disk. However, their model did not include the effect of the MHD-driven outflow. If massive stars form in a scaled-up, but similar manner to low-mass stars, then we expect such outflows to always be present around massive protostars during their main accretion phase. Indeed, larger-scale collimated outflows are often observed to be driven from massive protostars (e.g., Beuther et al., 2002; Garay et al., 2007; López-Sepulcre et al., 2009)

A series of massive protostellar evolution calculations to predict continuum radiative transfer, especially infrared to sub-mm morphologies and spectral energy distributions (SEDs), has been carried out by Zhang & Tan (2011), Zhang et al. (2013), and Zhang et al. (2014) (hereafter, ZTH14). Here we will calculate the ionization structure and associated radio emission predicted by these protostellar evolutionary sequences.

Study of radio emission from massive protostars has a long history (see, e.g., Wood & Churchwell, 1989; Hoare et al., 2007; Tan et al., 2014). Ultra-compact () and, especially, hyper-compact () H ii regions may include sources that are in the protostellar phase. Elongated radio continuum sources, e.g., radio “jets”, have been seen in a number of sources (e.g., Guzmán et al., 2012). The eventual goal of our work will be to compare the developed models with such observations, however, detailed comparison with individual sources, requiring tailored protostellar models, is deferred to future works in this series.

This paper is organized as follows. In §2 we describe our methods: in particular, we briefly review the protostellar evolution models, describe calculation of the photoionized structure and the radiative transfer of free-free continuum and HRL emission. Next, in §3, we present our results: especially, we explore the evolution of the photoionized regions for various initial core conditions and predict their observational features. In §4 we discuss the implications of our results, especially summarizing the evolution of the outflow confined H ii regions, and making an initial comparison of our models with the general characteristics of observed HC/UC H ii regions and radio winds/jets. We conclude in §5.

## 2. Methods

### 2.1. Evolution of protostars, disks, envelopes and outflows

We calculate the photoionization structure and present synthetic observations in the context of the Turbulent Core Accretion Model of MT03. We use a series of models of protostellar evolution and surrounding gas structure constructed by ZTH14, which self-consistently include the evolution of the protostar and disk formed from a rotating parent core embedded in a particular clump environment, and mass accretion regulated by disk-wind outflows. We note that these models are highly idealized in that they assume axisymmetry about a single protostar, while real objects will be more irregular, especially if the initial core is highly turbulent, and may be forming in a cluster environment with other protostars in the vicinity. Even if there is a relatively ordered and isolated core that is collapsing to a central disk, this collapse may lead to binary or low-order multiple star formation. Still, by first presenting the simplest case of formation of an isolated massive protostar, we can then see how relevant such models are as approximations to real systems.

Here we briefly review the evolution of the protostar and its surrounding gas structures. The initial core is assumed to be spherical and in near virial equilibrium via support by turbulence and/or magnetic fields. The evolutionary track of a core is determined by three environmental initial conditions: core mass , mass surface density of ambient clump (which sets surface pressure of the core), and ratio of the core’s initial rotational to gravitational energy . We fix the rotational parameter as a typical value from observations of lower-mass cores (Goodman et al., 1993; Li et al., 2012; Palau et al., 2013). The core density structure is described by a power law in spherical radius, . Values of are inferred from observations of dense cores in Infrared Dark Clouds (IRDCs) (Butler & Tan, 2012; Butler et al., 2014). We adopt as a fiducial value, following MT03. Then, the radius of a core is . Since the core is in pressure equilibrium with the surrounding clump, it is more compact for higher .

Table 1 presents the parameter sets of for which we study the ionization structures. These are the same as studied by ZTH14. We refer to model group A as a fiducial group with and . We also investigate high- and low- cases, model groups Al and Ah, and and 240 cases, as model groups B and C, respectively.

The core undergoes inside-out collapse (Shu, 1977; McLaughlin & Pudritz, 1997) with the effect of rotation described with the solution by Ulrich (1976) (see also Tan & McKee, 2004). Given the high accretion rates associated with collapse of dense cores, the protostellar disk is expected to be massive and self-gravitating such that angular momentum may be transported by torques due to, e.g., spiral arm structures. Considering this efficient angular momentum transfer, the mass ratio between disk and protostar is set to be a constant (e.g., Kratter et al., 2008). The disk structure is described by a standard -viscosity disk model (Shakura & Sunyaev, 1973), but also including the effects of infall and outflow (conservations of mass and angular momentum). The temperatures at the midplane are evaluated by vertical energy balance, and the vertical structure of the disk is determined by this temperature. The disk surface is set by assuming continuity with the density structure of the outflow. More details about the disk model are described in Zhang et al. (2013).

The mass accretion rate onto the forming protostar is

 ˙m∗ = 9.2×10−4ϵ∗¯ϵ0.5∗(m∗Mc)0.5 (1) ×(Mc60M⊙)3/4(Σcl1 gcm−2)3/4M⊙yr−1,

where and are the instantaneous and averaged star formation efficiencies. The instantaneous star formation efficiency is the ratio of the accretion rate onto the star to the idealized spherical infall rate from the initial cloud core in the absence of feedback: a certain percentage of infalling material is swept up by the disk wind (see next paragraph) and a fraction, , of the accreted mass remains in the disk during the protostellar phase. The averaged star formation efficiency is the ratio of the stellar mass to the idealized collapsed mass of the initial cloud core at the same stage of infall in the absence of feedback. The accretion rate is high for high clump mass surface densities, since high results in a more compact core and thus shorter collapse time. Also, this formula indicates that in the no-feedback limit of , the accretion rate increases with time, i.e., as (see MT03).

The magnetically driven bipolar outflows sweep up part of the core and thus help set the star formation efficiency . The density and velocity distributions of the outflow are described by a semi-analytic disk wind solution that is modified from the centrifugally driven outflow model of Blandford & Payne (1982) (for details, see Zhang et al., 2013, and ZTH14). The mass loading rate of the wind is proportional to the mass accretion rate, with a ratio assumed as a typical value (Königl & Pudritz, 2000). The velocity is as high as near the rotation axis and then smoothly decreases at lower latitudes. This outflow velocity structure represents features of both a collimated fast jet and a wider-angle, slower wind. The opening angle of the outflow cavity is estimated following the method of Matzner & McKee (2000): if the outflow momentum is strong enough to accelerate the core material to its escape speed, the outflow extends in that direction. The outflow momentum, , is evaluated by integrating the momentum rate over time. The momentum injection rate is typically in our model, where is the Keplerian speed at the stellar radius. This process sets the star formation efficiency (and also ). For more details of the outflow model and the evaluation of the star formation efficiency, see Zhang et al. (2013) and ZTH14.

Figure 1 shows an example of the evolution of the density and velocity structures in the fiducial case. The gradual opening-up of the outflow can be seen at the scale of with typical outflow velocity of which is high near the axis. This leads to a decrease in the instantaneous star formation efficiency. At the end of accretion in this model, the stellar mass reaches and the final averaged star formation efficiency is . For the different models, the averaged final star formation efficiencies are similar, with a range of . While we expect the momentum flux from protostellar outflows to be the dominant source of feedback, we note that additional feedback processes due to ionization, (line-driven) stellar winds and radiation pressure are not yet included in our modeling, so the above efficiencies may be somewhat overestimated. The growth of the disk radius can be seen in the panels showing the scale of . Note, partly for the reason of being near the end of its accretion phase and having a relatively low accretion rate, the disk in the case appears relatively thin compared to earlier stages (ZTH14). On the scales, we also see the narrow, inner, on-axis cavity of the disk wind, which is set to have zero (i.e., negligible) density. This cavity is too thin to see in the larger scale panels, but it is resolved by the logarithmic grid spacing.

The evolution of the protostar depends on its accretion history. ZTH14 calculated the protostellar properties with a detailed stellar evolution model by Hosokawa & Omukai (2009) and Hosokawa et al. (2010), self-consistently adapted to the Turbulent Core Model accretion rates. When the flow reaches the stellar surface, the accretion energy, , is released there. Following ZTH14, we treat this accretion luminosity and the internal stellar luminosity as being emitted isotropically with a single effective temperature, , i.e., . We will refer to as the stellar effective temperature. ZTH14 also included the effects of the accretion energy released at the disk surfaces in the dust temperature calculation. However, in this study, we neglect the ionizing photons from the disk surface, since its surface temperature is not high enough to emit significant amounts of ionizing flux. If the protostar rotates at almost its breakup velocity, it distorts the structure of the star. This effect leads to nonuniform surface temperatures: the polar region is hotter than the equatorial region (Meynet et al., 2010). However, the degree of rotation of massive protostars is quite uncertain, and thus for simplicity we adopt a spherical model with a single surface temperature as a first approximation.

Figure 2 shows the evolution of protostellar properties with protostellar mass: the stellar radius, luminosity and effective temperature were calculated by ZTH14, while the evaluation of ionizing flux will be described in §2.2.1. In the fiducial case with and , the stellar radius reaches a maximum at when the Kelvin-Helmholz time becomes as short as the accretion time, after which the protostar is able to begin contraction towards the zero-age main sequence (ZAMS) (Hosokawa & Omukai, 2009). Due to this variation in radius, the effective temperature increases strongly for . Comparing to the fiducial case, the mass at which the stellar radius has a maximum is larger in the higher case and smaller in the lower case. This is because a higher causes higher accretion rates (eq. 1), which make the accretion timescale shorter and also lead to increased swelling of the protostar.

### 2.2. Photoionization

#### 2.2.1 Ionizing photon luminosity

The total stellar H-ionizing photon luminosity is evaluated as

 S∗,acc=∫∞ν0Lν,∗,acchνdν, (2)

where is the frequency, is the Lyman limit frequency, and is the luminosity per frequency interval (i.e., the stellar spectrum). The stellar atmosphere model of Castelli & Kurucz (2004) is adopted to estimate , which covers wide ranges of effective temperature and surface gravity. The variation of ionizing luminosity is sensitive to the effective temperature and varies by orders of magnitude (bottom panel of Figure 2). In the fiducial case of group A, the ionizing luminosity is as small as at , and quickly increases to or larger in the Kelvin-Helmholz contraction phase. Even at the same mass of , the ionizing luminosity of model Ah12 is five orders of magnitude smaller than that of model Al12. Therefore, the accretion history and the detailed stellar evolution calculation are essential to study the ionization structure. We note that the obtained ionizing photon luminosities are somewhat lower than those in Figure 2 of Tan et al. (2014) due to improved calculation of the effects of line blanketing.

#### 2.2.2 Transfer of ionizing photons

We calculate the extent of the ionized regions using our previously developed radiative transfer code (Tanaka et al., 2013), which allows a treatment of both direct and diffuse ionizing radiation fields. We consider only the ionization of Hydrogen, ignoring Helium and other elements. If the stellar effective temperature is high, , we would underestimate the electron density by about due to singly ionized Helium, which would then enhance free-free emission by about .

The gray approximation for the photoionization calculation is adopted, where the frequency of ionizing photons is represented by the mean value. The equation of radiative transfer along a ray is

 dIds=−xInHσHI+α1x2IIn2H⟨hν⟩4π −nHσd,aI+nHσd,s(J−I), (3)

where is the number density of hydrogen nuclei, is the irradiance intensity, is the mean intensity, is the mean energy of ionizing photons, is the mean photoionization cross-section of the Hydrogen atom, is the radiative recombination coefficient to the ground state, and are neutral (atomic) and ionized fractions of H, and are absorption and scattering cross sections of dust grains per H nucleon. The first term of the right-hand side represents the consumption of ionizing photons by neutral Hydrogen, the second term describes the re-emission of ionizing photons by recombinations directly to the ground state, and the third and fourth terms correspond to absorption and scattering by dust grains, respectively. The local ionization fraction of H, , is calculated from the balance between photoionization, recombination and advection,

 ∇⋅(nHxIIv)=4πxInHσHJ⟨hν⟩−αAx2IIn2H, (4)

where is the velocity of the disk wind gas, is the recombination coefficient for all levels (so-called case A). The advection term on the left-hand side of this equation can be especially important near the ionization boundary. The recombination rates, and , are functions of temperature of the ionized gas. The mean values of the photon energy and Hydrogen cross-section are evaluated from the stellar spectra for each model,

 ⟨hν⟩ = 1S∗,acc∫∞ν0Lν,∗,accdν, (5) σH = 1S∗,acc∫∞ν0σH,ν(Lν,∗,acchν)dν, (6)

where is the frequency dependent cross-section of the Hydrogen atom. We adopt the fitting formulas from Draine (2011) for the recombination rates and , and the Hydrogen cross-section . The absorption and scattering properties of dust in these H ii regions are quite uncertain. We adopt and equal to zero in regions where the temperature history has exceeded the dust destruction temperature (1,400 K) in the thermal radiation calculation by ZTH14 (i.e., there is a dust-free inner region of the disk wind), and elsewhere they are fixed at from a typical value of the diffuse interstellar medium (Weingartner & Draine, 2001). The effect of variation of dust grain properties, i.e., sizes, abundances, compositions, is deferred to a future paper.

We conduct axisymmetric two dimensional calculations to solve equations (3) and (4), using our developed ray-tracing radiative transfer code. In the previous work of Tanaka et al. (2013), the code solved the transfer equation by the method of short characteristics, which involves only the cells next to the point under consideration (Stone et al., 1992). Here, we update the code implementing the method of long characteristics. With this method, the transfer problem is solved along the rays from the domain boundary to each point, which leads to less numerical diffusion and treats the diffuse radiation more accurately than the short characteristics method. The computational domain is a cylindrical region with radial and vertical coordinates and , assuming equatorial plane symmetry. The ionizing photon source is a sphere located at with a radius of , which represents the accreting star. We denote the distance from the center of the star as . The spatial grids for radial and vertical directions are set to be spaced logarithmically since the inner region has finer structures. The entire computational domain is resolved with a grid. The radiation intensity is a function of location and also direction: where and are the polar and azimuthal angles, respectively. The radiation emitted from the vicinity of the star has very fine angular structure when it reaches the outer boundary, i.e., . To resolve these fine structures while still saving computational costs, the angular resolutions of and are configured to vary depending on location . In the direction, the radiation arriving from directions with , which comes from the vicinity of the star, is resolved with a high resolution of , while the radiation arriving from is resolved by . In the direction, based on the tangential plane method (Stone et al., 1992), the radiation emitted from near the axis and from the opposite direction, and , is resolved with high resolution of , while the radiation from is resolved by . We also conduct resolution tests for the fiducial case with grids of and . The obtained ionized structures are consistent at these resolutions and free-free fluxes change by less than at radio frequencies from .

In this first work we study the photoionization of the disk wind, and do not allow photoionization of gas in the disk or infall envelope components of the model. Such photoionization would enhance the mass loss rate that is being injected into the disk wind, i.e., by photoevaporation, and would thus change its density structure and also the star formation efficiency. This feedback of photoevaporation will be discussed in a future paper.

#### 2.2.3 Temperature of the ionized gas

The temperature of the ionized gas is important to evaluate recombination rates and to study free-free emission for observational diagnostics. The photoionized gas temperature is regulated to around by the balance of heating and cooling processes. The dominant heating process is photoionization and the main cooling processes are radiative cooling of collisionally excited metal lines, recombination of H, and free-free emission.

It is computationally expensive to take into account all heating and cooling processes at every step of the radiative transfer calculation, and thus here we introduce an approximate treatment. The heating rate by photoionization depends on the stellar spectrum. Also, heating and cooling rates are both mostly sensitive to the density. Therefore, we treat the photoionized gas temperature, for each stellar model, as a function of the local density, i.e., . We evaluate from the volume average temperature of static H ii regions with uniform densities in the range using version 13.03 of CLOUDY (Ferland et al., 2013). This approximate approach is appropriate since a uniform density H ii region for each stellar model is almost isothermal as shown below. It should be noted that the ionized gas temperature is also sensitive to the metallicity of the gas, and we assume solar metallicity in this study.

The top panel of Figure 3 shows the radial structure of an H ii region with uniform density of formed by the stellar radiation of the stellar model A24. It can be seen that the H ii region has a sharp edge at where the ionization degree drops from near unity to near zero. The temperature is almost isothermal in the ionized region. We evaluate the typical photoionized gas temperature at as being from the volume average of this temperature for this stellar model.

In this way, we obtain photoionized gas temperatures over the density range of for each stellar model. The bottom panel of Figure 3 shows the photoionized gas temperature as functions of density, , for the stellar models A12, A16, and A24. The temperatures of all models are around for orders of magnitude of density variation. The temperature is typically higher at higher density. This is because the cooling is suppressed if the density exceeds the critical density of some of the cooling lines. We adopt these averaged temperatures to set as the local temperature of the ionized gas in our radiative transfer calculations.

Since the adopted temperature is evaluated based on static H ii region calculations, adiabatic cooling has not been included in this study. The adiabatic cooling rate, , is proportional to . On the other hand, the radiative cooling rates in ionized regions are almost proportional to the square of density. Since the density gradient in the disk wind is as steep as to , adiabatic cooling may not be negligible in low density regions far from the star. We compare the adiabatic cooling rate to the cooling rate evaluated by the static calculations. In the fiducial model of A16, adiabatic cooling would be of total cooling rate when and . As we will see in §3.2.1, the free-free flux at frequencies lower than is mainly provided from . Therefore, we note that the SEDs obtained by our simplified modeling would be overestimated somewhat at .

### 2.3. Observables

#### 2.3.1 Free-free emission

Using the obtained structures of the photoionized outflow, we estimate the free-free radio continuum emission. The equation of transfer for free-free emission is

 dIν,ffds=−κν,ff(Iν,ff−Bν), (7)

where and are the intensity and opacity of free-free emission, and is the Planck function. The opacity due to free-free absorption is given by

 κν,ff=43(2π3)1/2e6m3/2ehcn2Hx2HII√kBT ×[1−exp(−hν/kBT)ν3]gff, (8)

where and are the charge and mass of the electron, is the Planck constant, is the speed of light, is the Boltzmann constant, and is the Gaunt factor for free-free emission. The observed flux density is obtained by integrating over source solid angle of ,

 Fν,ff=∫ΩsIν,ffcos(θ−θview)dΩ=∫ΩsIν,ffdΩ, (9)

where is the viewing inclination measured from the outflow axis. In the last formula, we assumed that the distance is much farther than the object size, i.e., . In this study, for quantitative results, we assume protostars are located at a distance of , and total fluxes are obtained by integrating over a radius out to from the central star (for all our models this region dominates the total flux from free-free emission). The free-free emission is typically observed at , while the thermal dust emission would usually be dominant at higher frequencies. However, for completeness, we calculate free-free emission over the frequency range of .

#### 2.3.2 Hydrogen recombination lines

We also evaluate the HRL emission from the photoionized outflow. Considering a H line, i.e., the line from the transition , the opacity is

 κν,HRL = nnλ2n8π(n+1n)2 (10) ×An+1,n(1−e−hνn/kBT)ϕν,

where is the normalized line profile (), is the population number in energy level , is the Einstein coefficient of spontaneous emission, and and are the energy and wavelength of the H transition, respectively. Here, for simplicity, we assume the case of local thermodynamic equilibrium (LTE). However we note that, in the particular case of the massive star MWC 349A, maser amplification of HRLs are observed and thus non-LTE effects can be important (for modeling of MWC 349A including non-LTE effects, see Báez-Rubio et al., 2013).

The transfer for the total intensity, , i.e., the sum of the free-free and HRL emission, is

 dIν,ff,HRLds = −(κν,ff+κν,HRL)(Iν,ff,HRL−Bν), (11)

and the total observed flux is

 Fν,ff,HRL = ∫ΩsIν,ff,HRLcos(θ−θview)dΩ (12) = ∫ΩsIν,ff,HRLdΩ.

Here, we assumed that the source distance is greater than its size as in equation (9). We evaluate the HRL flux as the excess from the continuum flux at , i.e., . Of course, observationally it is not possible to observe the free-free flux exactly at since the line exists there and we cannot divide the total flux into two components. However, the value of is still able to evaluated from the far tail of the line profile or from interpolating from continuum observations at other frequencies. Since the typical line width is and the typical free-free emission spectral index in this study is , the variation of the continuum flux in this width is less than . Thus, it is safe to assume the continuum is constant over the line width. In this paper we will illustrate these results with the H line at . The flux ratio is typically for the H line in our models.

We assume that locally the gas has a Maxwellian velocity distribution, ignoring the intrinsic width of the line relative to the thermal and bulk-motion broadening. The velocity dispersion of Hydrogen is , and the corresponding dispersion in frequency is . The local normalized line profile combining the bulk gas velocity along the line of sight, , is

 ϕν=1√2πσ2νexp[−{(1+v/c)ν−νn}22σ2ν]. (13)

The frequency observationally represents a radial velocity of .

### 2.4. Escape fraction of ionizing photons

As we will see in next section, the ionized region, after initial early expansion, becomes unconfined in the radial direction along the outflow axis. Thus, some fraction of the ionizing photons escapes from the protostellar core and can ionize the surrounding, ambient clump gas, which may be forming a stellar cluster. Although the main focus of this work is the ionized structure of disk winds and their observational properties, ambient gas ionization is also important for its effect on star cluster formation and production of additional radio emission. Therefore, we provide estimates of the ionization of escaping photons and radio emission from the ionized clump gas.

The escape fraction of ionizing photons depends on the polar angle , i.e., in response to the ionized disk wind structure. We calculate the escape fraction as

 fesc(θ)=4πr2Fr,out(θ)S∗,acc⟨hν⟩, (14)

where is the radial energy flux of ionizing photons at the outer boundary, i.e., or . We note that the escape fraction is insensitive to the choice of the radius at which escape is defined as long as it is significantly larger than the disk radius. Integrating over the boundary, we obtain the total escape fraction,

 ¯¯¯fesc=14π∫fesc(θ)dΩ. (15)

## 3. Results

### 3.1. Structure and evolution of the photoionized outflow

First we report results of the fiducial case, model group A. Figure 4 shows the evolution of the photoionized structure. The high temperature region with is the photoionized gas, i.e., the H ii region, and elsewhere is neutral gas (with temperatures calculated from the ZTH14 radiative transfer models; note, we expect there would be some localized changes to the ZTH14 results for the neutral gas due to the presence of the H ii region and its associated photodissociation region (PDR), but we ignore those here—the PDR region will be considered in a future paper). Note that the white regions seen along the rotation axes in the -scale panels are disk wind inner cavities, where the density is negligible (formally zero in these models since we do not include a line-driven wind from the stellar surface) and thus where ionizing photons can travel freely. The ionized gas temperature is slightly higher in the inner region of the disk wind due to the higher density (see Fig. 3). When , the ionizing photon luminosity is too low to ionize the disk wind. As the ionizing luminosity increases with mass, the photoionized region is formed. When , the photoionized region is confined to intermediate zenith angles of . This is because densities are higher near the outflow axis at and near its launching zone for . Since the radial density gradient in the outflow is typically steeper than , the photoionized region is not confined in the radial direction and breaks out. The opening angle of the photoionized region increases as the protostellar mass grows. The ionized region eventually reaches the flared disk surfaces, however some parts of the disk wind remain neutral due to the shielding by the inner disk and disk wind even though the diffuse radiation is considered. Therefore, photoevaporation from the disk surface is suppressed, while photoevaporation from the envelope may still be efficient.

The ionizing radiation has two components: one is the direct stellar radiation and the other is the diffuse radiation from recombinations to the ground state and dust scattering. Figure 5 shows the fraction of diffuse radiation, , at each evolutionary stage. The direct radiation dominates over most of the outflow. However, in some locations, due to the direct radiation being shielded by the inner disk at , the diffuse radiation becomes more important, i.e., near the ionized-neutral boundary, e.g., as seen in model A24. In this way, as suggested by Hollenbach et al. (1994), the diffuse radiation is important to determine the structure of the H ii region and rates of photoevaporation from the neutral components. For these reasons, accurate treatment of the radiative transfer including the diffuse component is necessary.

To illustrate the structure near the rotation axis, Figure 6 shows the density and ionization fraction of model A16 with a logarithmic scale. The disk wind inner cavity, with zero density inside its innermost streamline, is seen as a very collimated structure with at (the black region in left panel and the white region in right panel). The density is high along the inner edge of the outflow. Due to the shielding by this high density inner wall, the gas in the region at high latitude with and remains neutral. This shielding is the same mechanism seen in model A12 in which the neutral region is at . The strip of ionized gas entering this neutral region is caused by high velocity advection near the axis (eq. 4). These thin structures near the axis are sensitive to our adoption of an idealized axisymmetric model. In reality, the outflow winds are not expected to be perfectly axisymmetric. For example, the jet in IRAS 16547-4247 appears to undergo precession (Rodríguez et al., 2008). However, the idealized structures in our model illustrate the wind shielding effect and also demonstrate the accuracy of our numerical calculation even at such fine angular scales.

Next we show results with different initial cores with lower and higher . The evolution of density and photoionization (as revealed by gas temperature, ) structures in groups Al and Ah are shown in Figures 7 and 8. In group Al, the photoionized region breaks out when , as in the fiducial case. On the other hand in group Ah, the photoionized region is tightly confined to a thin layer at the inner wind wall even at and only expands later. This is mostly because higher leads to higher accretion rates (eq. 1), which result in a higher mass at which the ionizing luminosity shows its dramatic increase due to Kelvin-Helmholz contraction (Fig. 2). We discuss the breakout of outflow-confined HII regions in §4.1.

Figures 9 an 10 show the evolution of density and photoionization (as revealed by gas temperature, ) structures in higher models of groups B and C. In both cases, the photoionized region breaks out when . This is similar to the high case, as the higher-mass core leads to delayed (in ) contraction via higher accretion rates. If we ignore feedback effects, the accretion rate is proportional to at the same stellar mass (eq. 1). In this way, the mass at which the photoionized region emerges varies from depending on its parent cloud core.

### 3.2. Free-free continuum emission

We conduct the radiative transfer calculation of free-free emission from the obtained photoionized structures. The total flux, spectrum, and size of the free-free emission of our model are similar to those of the radio winds and jets associating with massive star forming regions as we will see in §4.2. However, our ionized region model has a limited cylindrical volume (§2.2.2) and the contribution from external ionized regions is not included. Thus, it should be kept in mind that the free-free flux would be enhanced by the emission from the ionized clump if it is on the line of sight, i.e., .

#### 3.2.1 Spectral Energy Distributions

For all SED calculations, a distance of and an aperture diameter of (20,000 AU) are assumed. Figure 11 shows full SEDs including thermal dust and free-free emission for model A16 at viewing inclinations of and of the line of sight to the rotation/outflow axis. The flux from thermal dust emission was calculated by ZTH14. The free-free emission dominates over the thermal dust emission at radio frequencies of , while the dust flux dominates at higher frequencies. We note that, at wavelengths shorter than , we do not count the free-free emission in the total flux since we do not treat the dust extinction in the free-free radiative transfer calculation (eq. 7).

The infrared (IR) dust flux (at relatively short wavelengths) depends strongly on viewing angle. This is because the stellar radiation is absorbed by the dense disk and core envelope along directions close to the equatorial plane, and the re-emitted IR radiation, still suffering from high optical depths, escapes along the polar, low-density outflow cavity. Therefore, the observed flux is higher if the line of sight is passing through the outflow cavity (ZTH14). This process is known as the “flashlight effect,” which is also important for helping to overcome radiation pressure feedback (Nakano et al., 1995; Yorke & Bodenheimer, 1999; Krumholz et al., 2005, 2009; Kuiper et al., 2010, 2011).

On the other hand, the radio flux is almost independent of viewing angle since it typically experiences much smaller optical depths than the IR radiation. We can approximate the spectrum of free-free emission as being a power law in frequency, , with spectral index . From , the free-free flux from model A16 is fitted as . The spectral index of the free-free emission decreases from at 0.1 GHz to at 100 GHz. An index smaller than two indicates that the photoionized disk wind is only partially optically thick with respect to its free-free emission. The overall spectral index (including thermal dust emission) has a minimum of at , and quickly increases at higher frequencies approaching the dust spectrum value of . Figure 11 also indicates that the free-free flux at  GHz is mostly emitted from  AU scales, while that at  GHz is mostly probing  AU scales. Thus, the higher frequencies trace the structure of the inner wind.

Figure 12 shows SEDs for all models. The features of the free-free emission are similar to those discussed for model A16: the free-free fluxes at radio frequencies are almost independent of viewing angle, and are with the spectral indices of at . These common features are caused by the fact that the outflows of all models have similar structures, i.e., the amount of gas and the density gradient. If the photoionized region is small, the free-free emission is not observable, as in, e.g., model Al08. Once the photoionized region breaks out, the free-free emission becomes bright, as in, e.g., model Al12. Since the flux of free-free emission depends on the density structure of the outflow, the mass of the outflow that is emitting free-free flux (within a given distance from the protostar) may be estimated from its free-free flux.

With our focus in this paper on radio wavelengths, we note again that we have not accounted for the free-free emission for the total flux at wavelengths less than . However, in some cases the free-free flux might make a significant contribution in this range: e.g., at wavelengths shorter than of model C16 viewed at . The potential importance of free-free emission at shorter wavelengths will be explored in a future paper.

#### 3.2.2 Images

Next, we present images of the free-free emission. We show both resolved images and those after convolving with finite beam sizes of various telescopes. Detailed comparisons of models with individual sources are deferred to future papers in this series.

The top three panels of Figure 13 show the resolved free-free emission images of model A16 at the inclination of at three different frequencies, 1.5, 10, and 870 GHz. These frequencies correspond to the L and X bands of the Karl G. Jansky Very Large Array (J-VLA) and band 10 of the Atacama Large Millimeter/submillimeter Array (ALMA). Box sizes are , , and , respectively. The apparent size of the free-free images is larger at lower frequencies since free-free absorption is stronger (eq. 8). The effective radius at which half of the total flux is emitted is , with for . Thus these would typically be categorized as HC H ii regions (see also §4.2). At 1.5 GHz and 10 GHz, the apparent shape is similar to an hour-glass shape with a jet, because the disk wind has a self-similar structure at the scale of . The surface brightness is high in the central wind region and at the jet axis where the density is high, and the maximum brightness temperature of indicates optically thick conditions. At the higher frequency of 870 GHz, the shape is now more rounded in the inner region. This is because of the accretion disk with radius of . The collimated on-axis jet is not seen, since it becomes optically thin at this frequency. For the total free-free flux in Figure 11, the innermost wind region dominates the on-axis jet since its area is bigger. However, the bright linear feature is still very interesting, since it is potentially related to observed radio jets (Gibb et al., 2003; Rodríguez et al., 2008; Guzmán et al., 2012, 2014).

The bottom three panels of Figure 13 show the free-free images convolved with the finite beam sizes of the telescopes at these frequencies. The full-width at half-maximum (FWHM) beam sizes are 1.3, 0.2, and 0.01  for 1.5, 10, and 870 GHz, respectively, i.e., the highest resolutions in the L and X bands of the J-VLA and in band 10 of ALMA. A distance of 1 kpc is assumed, and the point-spread function is approximated as a Gaussian function. The jet feature with high brightness temperature along the axis is diluted at 1.5 GHz and 8.4 GHz, since it is narrower than the beam size. However, the linear feature is still observable at the level of or at both frequencies. The sensitivity with root-mean-square value of can be attained with an observing time of 1 hour and a bandwidth of 500 MHz by the J-VLA. The inner brightest region with is smaller than the beam size, and thus the highest brightness temperature decreases to at 870 GHz. The dust emission cannot have a brightness temperature above the dust sublimation temperature of . Therefore, even though the total flux is dominated by dust emission at the frequencies (Fig. 11), the free-free emission may be observable even at this high frequency. Since most massive protostars are at somewhat larger distances, this emission may need to be observed in more massive, denser sources.

Figure 14 shows resolved and convolved images at 10 GHz as the protostar evolves, viewed at inclination angles of and for all models. From the resolved images (top panels of Fig. 14), it can be seen that the outflow opening angle gradually increases, especially at the inclination of , once the photoionized region is formed. When the line of sight passes through the outflow cavity (), like model A16 at inclination, then the image becomes elliptical rather than hourglass shaped. The outflow axis always exhibits a high brightness at this frequency. As stated above, the linear jet-like feature is observable at the level of in the convolved image of model A16. In the case of higher- or , the brightness temperature is higher since the outflow rate is higher. Thus the jet feature is observable even at in these models, such as models Ah24 and C32. The small spots in convolved images of models Ah16 and B12 are the sign of unresolved outflow-confined H ii regions.

### 3.3. Hydrogen recombination lines

We perform radiative transfer of the H line to obtain simulated observations that probe the kinematics of the photoionized wind. Here, a distance of and aperture diameter of (20,000 AU) are assumed.

Figure 15 shows profiles of H flux from model A16. The line profiles are very similar for the inclination angles of and : the FWHM are , the peak fluxes are , the ratios of line flux to the free-free continuum flux are and , respectively. Although we used a Gaussian line profile of at each position for the radiative transfer calculation (eq. 13), the obtained total spectrum has broader wings. This is because the disk wind has higher velocity components, , near the axis.

Figure 16 shows resolved and convolved maps of H of model A16 at three different velocities, , at inclinations of and . The map is convolved with a beam size of 0.1 ″ which is the highest resolution in band 3 of ALMA. For this ionized bipolar outflow, the blueshifted component appears on the north (upper) side and the redshifted component in south (lower) side. However, the maps are tilted slightly to northwest-southeast direction. This is due to the rotation of the outflow. Rotating outflows have been found around low-mass and intermediate-mass protostars (Klaassen et al., 2013; Codella et al., 2014), but not yet for high-mass objects. In resolved maps at an inclination of , the line of sight is outside the cavity and thus the rim of the outflow is brightest. On the other hand, in resolved maps at an inclination of , the line of sight passes through the cavity and the maps are elliptical illustrating the outflow cone. This difference becomes less clear in the convolved maps at . However, it is still visible in the maps at the level of . The sensitivity of could be attained with an observing time of 8 hours with a channel width of 6 MHz by ALMA. The effective spectral resolution is two times of the channel width which corresponds to 12 MHz or (the H40 HRL is at 99 GHz). We should mention that again a distance of 1 kpc is assumed here, and most massive protostars are typically at larger distances. Thus, this emission may need to be observed in more massive and denser sources.

### 3.4. Escape fraction of ionizing photons

So far we presented results for the ionization of disk winds. Next, we consider the ionizing photons that escape along the outflow cavities created by these winds. Since the obtained ionized regions are unconfined in the radial direction, some fraction of ionizing photons escape and are expected to ionize the ambient clump gas. This ionization feedback, in addition to that of the disk wind protostellar outflows, may affect star formation around the forming massive star and will also lead to additional radio emission.

The top panel of Figure 17 shows the escape fraction of ionizing photons to each polar angle, , along with the opening angle of the disk wind cavity, . In the case of model A12, the escape fraction is at its maximum value and almost zero for or . This is because the ionized region is restricted to (Fig. 4). As the ionized region expands in models A16 and A24, the angular range in which photons escape also expands and the maximum value of reaches about unity. However, the escape fraction is still in directions near the boundary of the core infall envelope (i.e., : see also Fig. 5) or near the outflow axis (: see also Fig. 6) due to strong attenuation in these higher density regions.

The bottom panel of Figure 17 shows the evolution of the outflow cavity opening angle, , and the total escape fraction of ionizing photons, , for all model groups. It can be seen that the total escape fraction increases as the opening angle grows in every case. However, we note that, due to Hydrogen ionization and dust absorption, the total escape fraction is lower than the fraction of solid angle extended by the disk wind cavity. The maximum value of in our models is with , so more than half of ionizing photons do not escape from the system (including the directions blocked by the disk and infalling core envelope).

The total ionizing photon rate escaping from the protostellar core is . Most of these photons will be able to ionize the ambient clump or larger-scale cloud gas, especially the gas in front of the head of the outflow bowshock. Our model does not explicitly include the structure of the outflow bowshock or the ambient clump gas. Thus, we do not study details of the geometry or kinematics of the ionized clump. Instead, as a simple first approximation, we use a Strömgren sphere analysis to make some simple estimations. In this case, the approximate size scale of the ionized region created by escaping photons is

 RSt,cl=⎛⎝3¯¯¯fescS∗,acc4πn2H,clαB⎞⎠1/3, (16)

where is the density of the ambient clump, and is the case B recombination rate (this estimate ignores absorption of photons by dust). The mean clump density is given by . The recombination rate, , is a function of temperature, which is calculated from the stellar spectrum and the clump density (§2.2.3). The free-free emission is evaluated by solving equation (7) assuming spherical geometry. The FWHM of the velocity distribution of the ionized clump gas, which will be used as part of the diagnostic process to classify H ii regions (§4), is estimated simply from thermal broadening, . This estimation gives the minimum value of the line width, since the clump gas could also be disturbed by outflows and additional turbulent motions that broaden the line profile.

We note some caveats of this Strömgren sphere approximation for ambient clump ionization. First, the geometry of the ionized region is likely to be affected by the outflow cavity, leading to deviations from a simple sphere. Second, the ambient clump density is expected to be inhomogeneous, which can lead to various effects, including enhanced effective densities at ionization fronts that can lead to smaller H ii region sizes (e.g., Tan & McKee, 2004). Thus, the free-free flux evaluated by the Strömgren sphere analysis is only expected to be a very crude approximation. If the ambient clump region is optically thin, its total free-free flux from the ionized region is just proportional to the total escaping photon rate, , regardless of its shape or density. However, if the ionized region is optically thick, the free-free flux is smaller than the optically thin limit by the factor of its mean optical depth. Therefore, the error in the free-free flux would be small if the ionized clump region is optically thin.

## 4. Discussion

### 4.1. Outflow-confined H ii regions and their break out

Here we summarize the evolution of ionized outflows. Note, the increase of ionizing flux in the Kelvin-Helmholz contraction phase is very rapid (Fig. 2) and our numerical results have only been carried out with relatively coarse sampling across the evolutionary sequence. However, we can still present a description of the following evolutionary sequence. At the first stage, when the ionizing flux is small, the MHD-driven outflow is the first barrier to the propagation of ionizing photons, except in the very narrow range of directions within the inner disk wind cavity. The H ii region is confined in a very thin layer of the inner outflow wall, e.g., model Ah16 in Figure 8. As the ionizing flux increases, the H ii region extends to greater distance from the protostar. At a critical ionizing flux, the H ii region breaks out at a middle latitude (see model A12 in Fig 4). At the lower latitude, the density is higher than that in the middle latitude since it is closer to the launching zone of the disk wind. On the other hand, the density gradient at the higher latitude is flatter than that at the middle latitude due to the collimated structure of the wind. Therefore, the H ii region breakout starts from . The critical ionizing flux for the radial breakout depends on the detailed structure of outflow, i.e., mass loss rate, opening angle and velocity structure. From a qualitative standpoint, the mass loss rate is roughly proportional to the mass accretion rate, and thus the critical ionizing flux for the breakout is higher for higher- cases. The breakout does not occur at even in the low- case of model Al08. The typical flux for breakout is in our numerical calculations, which is as high as the flux from the ZAMS at about . All models in this work reach at the final stage, however, it should be noted that there are almost as many stars between and as above . In the formation of such stars with , the H ii region is confined during the main accretion phase and breakout occurs only towards the end of accretion when the outflow rate also decreases.

After first breakout at , the photoionized region extends in polar angle. To ionize almost the entire outflow, the ionizing flux needs to increase by about one order of magnitude from the first breakout ionizing flux. This takes only since the ionizing flux increases rapidly in the KH contraction phase. Once the entire outflow is ionized, the expansion in polar angle of the H ii region slows down matching the increase of the outflow cavity opening angle, which has a timescale similar to the accretion time, i.e., .

Our model of photoionized disk winds may also help to explain the short timescale variability that is seen in some UC/HC H ii regions (Franco-Hernández & Rodríguez, 2004; Rodríguez et al., 2007; Galván-Madrid et al., 2008) and which has also been modeled in the radiative hydrodynamical simulations of Peters et al. (2010). Fluctuations in the accretion rate would also lead to changes in the outflow massloss rate, which would change the density and thus the extent of the photoionized region. Changes could happen on outflow crossing timescales, i.e., as short as  pc /  yr.

Our model makes prediction for the spectral and the size indices, i.e., and , of the free-free emission. The free-free emission spectral index at 10 GHz is about (see §3.2.1), which indicates that the outflow is partially optically thick. The apparent size of the observed free-free emission (see §3.2.2) is as small as at 10 GHz and smaller at higher frequencies, scaling with size as , with the index of , even though the ionized region is larger than the core scale of after breakout. These frequency dependences are related to the density structure of the ionized outflow. A simplest analytical structure is a spherical wind with constant velocity, i.e., where and are the wind mass loss rate and velocity (Panagia & Felli, 1975). Since the opacity for the free-free emission at radio frequencies is approximatly if the gas is isothermal, the optical depth rapidly decreases as . Therefore, the outer region becomes optically thin and has less contribution to the total free-free flux even though the ionized gas extends to larger scales. In this spherical model, the radius of the optically thick region scales as , and thus the total free-free flux, which is proportional to the flux from the optically thick region, is . Additionally, as shown by Reynolds (1986), even with a bipolar jet geometry, the set of indices will be similar to the simple spherical wind if the jet is conical (constant aperture angle), of constant velocity and ionization fraction, and is isothermal. Thus the departure from the indices of is related to acceleration, collimation, thermal structure, or variation of ionization fraction (Reynolds, 1986; Guzmán et al., 2014). Our wind model is more complex compared to the spherical wind or the conical jet, however the fundamental physics is the same and thus the values of its spectral and size indices are similar.

### 4.2. Comparison with observed radio sources around massive protostars

Here, we compare general properties of our model with observed radio sources around massive protostars and show that the features of our theoretical models are similar to those of observed radio winds and jets.

Figure 18 shows the properties of size, line width, radio and bolometric luminosities, and outflow momentum rate of our models and of observed radio sources. For our models, size and free-free luminosity are measured at and line width is evaluated from H assuming an inclination angle of . The size is evaluated by the diameter at which half of the total radio flux of the system is emitted. As we have seen, the inclination angle does not affect these properties significantly at radio frequencies. Note, the outflow momentum rate, , of our model shown in Figure 18 is obtained from the protostellar evolution calculation of ZTH14 rather than from an estimation via a synthetic observation.

In panel (a) of Figure 18, we show how Hoare et al. (2007) distinguished objects into three categories based on size and line width: UC H ii regions are and , HC H ii regions are and , and massive young stellar object (MYSO) wind sources are and . It should be noted that this classification of the radio sources is somewhat arbitrary. For example, Kurtz (2002, 2005) defined UC H ii regions as having diameters and HC H ii regions as with emission measures . Tan et al. (2014) adopts a simpler definition only based on the above sizes. When we refer to data from Hoare et al. (2007) and Hoare & Franco (2007), we follow Hoare et al.’s classification. The line width of UC and HC H ii regions are derived from radio HRLs, while those of wind sources are from IR HRLs. The wind sources (BN, W33A, S140 IRS 1, NGC 2024 IRS 2, GL 989 and GL 490 are plotted here) are the smallest and fastest objects. Comparing in those plots, our models reproduce the properties of wind sources quantitatively, with sizes of and line widths of . Two of our models have broader line widths than ( for A12 and for C16). These two models are at the beginning of H ii region expansion, and the ionized region is confined only near the axis where the wind velocity is highest. Thus, their line width is exceptionally broad.

The difference between UC/HC H ii regions and MYSO wind sources is more clearly seen in the ratio of radio and bolometric luminosities in panel (b) of Figure 18. UC/HC H ii regions are more radio-loud than MYSO wind sources by two orders of magnitude or more. Even though the average radio to bolometric luminosity ratio of our models are about one order of magnitude higher than those of MYSO wind sources, there is still overlap given their dispersions.

Anglada et al. (2015) found two empirical correlations related to the radio luminosity of observed radio jets: one is a correlation with bolometric luminosity; the other is with outflow momentum rate:

 Fνd2mJykpc2 = 8×10−3(LbolL⊙)0.6; (17) Fνd2mJykpc2 = 190(˙pwM⊙yr−1kms−1)0.9. (18)

Interestingly, those two correlations are consistent over a wide bolometric luminosity range of to (dotted lines in panels (c) and (d) of Fig. 18). MYSO wind sources from Hoare & Franco (2007) (black circles) also fall on the first correlation. The solid line in panel (b) is the expected radio luminosity in the optically thin limit using the ionizing photon rate of ZAMS models by Thompson (1984), which represents the maximum radio luminosity from photoionized gas. At , however, there are radio jets that exceed this maximum value (Anglada et al., 1992). It has been suggested that the single power law correlations of equations (17) and (18) indicates a common nature of (ionized) outflows from low-mass and high-mass protostars (Rodríguez et al., 1989b; Anglada, 1995). Since the photoionization properties differ greatly from low-mass to high-mass stars, this might imply that a different mechanism, i.e., shock ionization, is responsible for the observed radio emission along these correlations. However, we find that the radio luminosities of our model agree reasonably well with both the correlations of eqs. (17) and (18), even though they are photoionized structures. We predict that there should be a transition to photoionization-dominated winds for massive protostars, but that in this regime the radio/bolometric luminosities and the outflow momentum rate may not differ too much from this empirical relation. The systematic differences of our models are towards somewhat higher values of radio luminosities by factors of 1.5–13. This may be a signature of the initial stage of photoionization leading to enhanced radio fluxes. Note, our models tend to be more massive than the observed sources, so more data in this range is needed.

If the photoionized region is optically thin, the free-free emission is simply proportional to the ionizing photon flux rate (Hummer & Seaton, 1963). However, our model indicates much lower free-free flux than the optically thin limit (panel (c) in Fig. 18). First, this is because more than half of the emitted ionizing photons are absorbed at the innermost region of . This region is optically thick for the radio free-free emission, and thus the absorption there does not contribute to the total free-free flux. In the case of model A16, of photons are absorbed inside . Second, most of the rest of the photons, , escape from the system. Those escaping photons do not induce free-free flux from the core scale, but contribute to the flux from the ambient clump gas, which is discussed below (§4.3). The final reason is that the photoevaporation flow has so far not been included in our model. If the ionizing photon rate is high enough to ionize the disk and envelope surfaces, e.g., in the case of model A24, the outflow rate of ionized gas increases and the free-free flux should be enhanced. While the first and second reasons are physical, the last one is due to the incompleteness of our model. We will address photoevaporation and its feedback effects in a subsequent paper.

Our ionized outflow models typically have spectral indices with (Fig. 12), and size index with frequency of from free-free emission at 10 GHz (Fig. 13). As discussed above, those indices are related to the density gradient of the outflow. The indices obtained by our calculations are consistent with some observed radio jets, such as MWC 349 ( and , Tafoya et al., 2004) and W75N(B) ( and , Carrasco-González et al., 2015). G345.4938+01.4677, however, has the slightly steeper indices of and than our model, which may be interpreted as being due to a stronger degree of collimation and acceleration (Guzmán et al., 2014, Guzmán et al., in prep.).

We should notice that, even though our theoretical model has similar properties to those of the observed radio jets, there is still a lack of some physical processes in our model. The velocity of some jets is estimated to be as high as from the proper motion of associated lobes along the axes, e.g., HH 80-81 (Martí et al., 1998, 1995; Masqué et al., 2015) and Cep A HW 2 (Curiel et al., 2006). Similar to these sources, our wind models have high radio brightnesses and high velocities of near the axes (Figs. 1, 7, etc.). However, unlike our model, some of lobes are found to have the radio emission with negative spectral indices, such as Serpens (Rodríguez et al., 1989a), HH 80-81 (Martí et al., 1993), Ceph A HW 2 (Garay et al., 1996), W3(OH) (Wilner et al., 1999), and IRAS 16547-4247 (Rodríguez et al., 2005). These negative spectral indices can be interpreted as due to nonthermal synchrotron emission induced by strong shocks with the ambient gas (Araudo et al., 2007), and this has been confirmed by the detection of linearly polarized emission in HH 80-81 (Carrasco-González et al., 2010). Additionally, our simple axisymmetric model has not accounted for the possibility of outflow axis precession, which is apparent in IRAS16547-4247 (Rodríguez et al., 2008). Although the current number of sources with these features is relatively small, inclusion of such additional physics in our theoretical models, i.e., shock ionization and outflow axis precession, are likely to be needed for a complete understanding of jet properties in massive star formation. Such analysis is also deferred to a future paper.

As we have seen, the properties of our models are consistent with the so-called “MYSO wind sources” in Hoare et al. (2007) and Hoare & Franco (2007), and “radio jets” with high bolometric luminosity in Anglada et al. (2015). In general, the main difference of “winds” and “jets” is the shape, i.e., winds have relatively wide opening angles and jets are more collimated bipolar outflows. This difference could be related to evolutionary stage (see also Fig. 14). To clarify our view of the evolutionary sequence of massive star formation and its outflow launching mechanism, further study is needed utilizing multi-wavelength observations to be compared with theoretical models.

### 4.3. Ionization of ambient clump gas by escaping photons

Finally, we comment on the ionization of the ambient clump. Tens of percent of the total ionizing photon output escape along the outflow cavity (Fig. 17). These photons are likely to ionize the ambient clump which would lead to enhanced radio luminosity. We estimate features of these extended ionized regions, as described in §3.4. Figure 19 shows the comparison of our ionized clump models and observed sources (the observational data are the same as in Fig. 18). The Strömgren-based models of ionized ambient clumps typically have the size of and line widths of , and large radio luminosities of , which are larger in size, narrower in line width, and more radio loud than the ionized outflow models. These features are roughly consistent with observed UC/HC H ii regions, except the line widths are broader in the observed systems. However, this difference is not surprising as the theoretical models include only thermal broadening and so give minimum widths.

We conclude that some UC/HC H ii regions may be formed by the escaping photons from massive protostars with radio jets or winds. In this study, the ionized clump features are evaluated based on the simple Strömgren analysis (§3.4). For more detailed analysis, planned for future work, one needs to model more accurately the structure and kinematics of ambient clump gas, including its perturbation by the outflow, to better understand the connection of radio winds/jets and UC/HC H ii regions.

## 5. Summary

We have studied the structure, evolution and observational properties of photoionized disk wind outflows that are expected to be associated with massive protostars forming by Core Accretion. Specifically, we have followed the quantitative properties of the Turbulent Core Accretion Model of McKee & Tan (2003) and follow-up protostellar evolution, disk wind and IR continuum radiative transfer models of Zhang & Tan (2011) and Zhang et al. (2013, 2014).

We calculated the transfer of ionizing photons and obtained the photoionized structure. The disk wind starts to be ionized when the stellar mass reaches about , depending on the properties of the initial cloud core, due to the stellar temperature increase resulting from Kelvin-Helmholz contraction. At this moment, the photoionized region is confined to zenith angles of about , while the high density regions near the outflow axis and equatorial plane remain neutral. On the other hand, in the radial direction, the photoionized region is unconfined, since the density gradient of the disk wind is steep. The disk wind is almost entirely ionized in about , and the photoionized region extends as the outflow opening angle increases. The outer disk surface is not ionized even though the diffuse radiation has been considered. This is because of the shielding effects of the inner disk and the disk wind.

Using the obtained photoionized structures, we calculated radiative transfer of free-free continuum and Hydrogen recombination line emission. The free-free emission exceeds dust emission at the radio frequencies lower than once the ionized region is formed. While shorter wavelength infrared flux from dust emission can depend on inclination of viewing angle, the radio flux from free-free emission is almost independent of inclination and its typical flux density is at a distance of 1 kpc with a spectral index The apparent size depends on frequency as approximately . The profile of has a FWHM of with broad wings. These properties of our model are similar to those of observed radio winds and jets from massive protostars. Detailed comparison with individual sources is deferred to future papers in this series.

## Acknowledgments

The authors thank Christopher McKee, Andres Guzmán, Taishi Nakamoto, Hiroyuki Kurokawa, Kohei Inayoshi, and Shuo Kong for fruitful discussions and comments. The authors also thank the anonymous referee for comments, which were useful to improve the original manuscript. JCT acknowledges support from NSF grant AST 1411527.

## References

• Anglada et al. (1992) Anglada, G., Rodríguez, L. F., Canto, J., Estalella, R., & Torrelles, J. M. 1992, ApJ, 395, 494
• Anglada (1995) Anglada, G. 1995, Revista Mexicana de Astronomia y Astrofisica Conference Series, 1, 67
• Anglada et al. (2015) Anglada, G., Rodríguez, L. F., & Carrasco-Gonzalez, C. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 121
• Araudo et al. (2007) Araudo, A. T., Romero, G. E., Bosch-Ramon, V., & Paredes, J. M. 2007, A&A, 476, 1289
• Báez-Rubio et al. (2013) Báez-Rubio, A., Martín-Pintado, J., Thum, C., & Planesas, P. 2013, A&A, 553, AA45
• Beuther et al. (2002) Beuther, H., Schilke, P., Sridharan, T. K., et al. 2002, A&A, 383, 892
• Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
• Bunn et al. (1995) Bunn, J. C., Hoare, M. G., & Drew, J. E. 1995, MNRAS, 272, 346
• Butler & Tan (2012) Butler, M. J., & Tan, J. C. 2012, ApJ, 754, 5
• Butler et al. (2014) Butler, M. J., Tan, J. C., & Kainulainen, J. 2014, ApJ, 782, LL30
• Carrasco-González et al. (2010) Carrasco-González, C., Rodríguez, L. F., Anglada, G., et al. 2010, Science, 330, 1209
• Carrasco-González et al. (2015) Carrasco-González, C., Torrelles, J. M., Cantó, J., et al. 2015, Science, 348, 114
• Castelli & Kurucz (2004) Castelli, F., & Kurucz, R. L. 2004, IAU Symp. No 210, Modelling of Stellar Atmospheres, eds. N. Piskunov et al. 2003, poster A20, arXiv:astro-ph/0405087
• Codella et al. (2014) Codella, C., Cabrit, S., Gueth, F., et al. 2014, A&A, 568, L5
• Curiel et al. (2006) Curiel, S., Ho, P. T. P., Patel, N. A., et al. 2006, ApJ, 638, 878
• Dale et al. (2005) Dale, J. E., Bonnell, I. A., Clarke, C. J., & Bate, M. R. 2005, MNRAS, 358, 291
• Davies et al. (2011) Davies, B., Hoare, M. G., Lumsden, S. L., et al. 2011, MNRAS, 416, 972
• Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium by Bruce T. Draine. Princeton University Press, 2011. ISBN: 978-0-691-12214-4
• Dupree & Goldberg (1970) Dupree, A. K., & Goldberg, L. 1970, ARA&A, 8, 231
• Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, RMxAA, 49, 137
• Figer (2005) Figer, D. F. 2005, Nature, 434, 192
• Franco-Hernández & Rodríguez (2004) Franco-Hernández, R., & Rodríguez, L. F. 2004, ApJ, 604, L105
• Galván-Madrid et al. (2008) Galván-Madrid, R., Rodríguez, L. F., Ho, P. T. P., & Keto, E. 2008, ApJ, 674, L33
• Garay et al. (2007) Garay, G., Mardones, D., Bronfman, L., et al. 2007, A&A, 463, 217
• Garay et al. (1996) Garay, G., Ramirez, S., Rodriguez, L. F., Curiel, S., & Torrelles, J. M. 1996, ApJ, 459, 193
• Gibb et al. (2003) Gibb, A. G., Hoare, M. G., Little, L. T., & Wright, M. C. H. 2003, MNRAS, 339, 101
• Goodman et al. (1993) Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528
• Guzmán et al. (2012) Guzmán, A. E., Garay, G., Brooks, K. J., & Voronkov, M. A. 2012, ApJ, 753, 51
• Guzmán et al. (2014) Guzmán, A. E., Garay, G., Rodríguez, L. F., et al. 2014, ApJ, 796, 117
• Hoare et al. (2007) Hoare, M. G., Kurtz, S. E., Lizano, S., Keto, E., & Hofner, P. 2007, Protostars and Planets V, 181
• Hoare & Franco (2007) Hoare, M. G., & Franco, J. 2007, A Volume Honouring John Dyson, Edited by T.W. Hartquist, J. M. Pittard, and S. A. E. G. Falle. Series: Astrophysics and Space Science Proceedings. Springer Dordrecht, 2007, p.61, arXiv:0711.4912
• Hollenbach et al. (1994) Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654
• Hosokawa & Omukai (2009) Hosokawa, T., & Omukai, K. 2009, ApJ, 703, 1810
• Hosokawa et al. (2010) Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478
• Hosokawa et al. (2011) Hosokawa, T., Omukai, K., Yoshida, N., & Yorke, H. W. 2011, Science, 334, 1250
• Hummer & Seaton (1963) Hummer, D. G., & Seaton, M. J. 1963, MNRAS, 125, 437
• Keto (2007) Keto, E. 2007, ApJ, 666, 976
• Klaassen et al. (2013) Klaassen, P. D., Juhasz, A., Mathews, G. S., et al. 2013, A&A, 555, A73
• Königl & Pudritz (2000) Königl, A., & Pudritz, R. E. 2000, Protostars and Planets IV, 759
• Kratter et al. (2008) Kratter, K. M., Matzner, C. D., & Krumholz, M. R. 2008, ApJ, 681, 375
• Krumholz et al. (2009) Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science, 323, 754
• Krumholz et al. (2005) Krumholz, M. R., McKee, C. F., & Klein, R. I. 2005, ApJ, 618, L33
• Krumholz et al. (2007) Krumholz, M. R., Stone, J. M., & Gardiner, T. A. 2007, ApJ, 671, 518
• Kuiper et al. (2010) Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010, ApJ, 722, 1556
• Kuiper et al. (2011) Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2011, ApJ, 732, 20
• Kurtz (2002) Kurtz, S. 2002, Hot Star Workshop III: The Earliest Phases of Massive Star Birth, 267, 81
• Kurtz (2005) Kurtz, S. 2005, Massive Star Birth: A Crossroads of Astrophysics, 227, 111
• Lanz & Hubeny (2007) Lanz, T., & Hubeny, I. 2007, ApJS, 169, 83
• Li et al. (2012) Li, J., Wang, J., Gu, Q., Zhang, Z.-y., & Zheng, X. 2012, ApJ, 745, 47
• López-Sepulcre et al. (2009) López-Sepulcre, A., Codella, C., Cesaroni, R., Marcelino, N., & Walmsley, C. M. 2009, A&A, 499, 811
• Martí et al. (1993) Martí, J., Rodríguez, L. F., & Reipurth, B. 1993, ApJ, 416, 208
• Martí et al. (1995) Martí, J., Rodríguez, & Reipurth, B. 1995, ApJ, 449, 184
• Martí et al. (1998) Martí, J., Rodríguez, L. F., & Reipurth, B. 1998, ApJ, 502, 337
• Martins & Plez (2006) Martins, F., & Plez, B. 2006, A&A, 457, 637
• Masqué et al. (2015)