Outflow-Confined HII regions. II. The Early Break-Out Phase

In this series of papers, we model the formation and evolution of the photoionized region and its observational signatures during massive star formation. Here we focus on the early break out of the photoionized region into the outflow cavity. Using results of 3-D magnetohydrodynamic-outflow simulations and protostellar evolution calculations, we perform post-processing radiative-transfer. The photoionized region first appears at a protostellar mass of 10Msun in our fiducial model, and is confined to within 10-100AU by the dense inner outflow, similar to some observed very small hypercompact HII regions. Since the ionizing luminosity of the massive protostar increases dramatically as Kelvin-Helmholz (KH) contraction proceeds, the photoionized region breaks out to the entire outflow region in<10,000yr. Accordingly, the radio free-free emission brightens significantly in this stage. In our fiducial model, the radio luminosity at 10 GHz changes from 0.1 mJy kpc2 at m=11Msun to 100 mJy kpc2 at 16Msun, while the infrared luminosity increases by less than a factor of two. The radio spectral index also changes in the break-out phase from the optically thick value of 2 to the partially optically thin value of 0.6. Additionally, we demonstrate that short-timescale variation in free-free flux would be induced by an accretion burst. The outflow density is enhanced in the accretion burst phase, which leads to a smaller ionized region and weaker free-free emission. The radio luminosity may decrease by one order of magnitude during such bursts, while the infrared luminosity is much less affected, since internal protostellar luminosity dominates over accretion luminosity after KH contraction starts. Such variability may be observable on timescales as short 10-100 yr, if accretion bursts are driven by disk instabilities.


INTRODUCTION
Massive stars are the main sources of energetic and chemical feedback in most galaxies via their radiation, wind, and supernovae. However, the formation of massive stars is not understood well compared to low-mass star formation . One widely discussed class of models is based on the Core Accretion scenario, e.g., the Turbulent Core model of McKee & Tan (2003). These models are scaled-up versions of the standard scenario of low-mass star formation (Shu et al. 1987): massive stars are formed from massive prestellar cores that are supported largely by nonthermal pressure forces, i.e., due to turbulence and magnetic fields. However, there is also significant differences compared to low-mass star formation, such as the degree of photoionization around the protostar. Massive stars are expected to start nuclear burning while still accreting. Their luminosity and photospheric temperature become high and they then radiate vast amounts of Lyman continnum photons with energies above 13.6 eV, which ionize H to create H II regions. These protostellar H II regions emit radio free-free continuum, providing observational clues of the massive star formation process.
The radio emission from massive protostars has been studied for a long time (e.g., Mezger & Henderson 1967;Wood & Churchwell 1989;Kurtz et al. 1994;Hoare et al. 2007). Ultra-compact (≤ 0.1 pc) and, especially, hypercompact (≤ 0.01 pc) H II regions may include sources that are in the protostellar phase. Elongated radio continuum sources, refereed to as radio jets, have been seen in a number of sources (e.g., Guzmán et al. 2012). Recently, Rosero et al. (2016) performed a high-sensitivity radio continuum survey with the Karl G. Jansky Very Large Array (VLA) toward early stages of massive starforming molecular clumps. They showed that the radio detection rates increase from cold infrared dark clouds (IRDCs) to molecular cores, which suggests the radio emission is a useful tracer of the evolutionary stages of massive star formation. Theoretical models are essential to interpret these observational signatures of evolution.
Here we present the second of a series of papers on modeling the photoionization process and its observa-tional signatures of the Turbulent Core model of Mc-Kee & Tan (2003). In the same context, Zhang et al. (2014) and  have developed a suite of radiative transfer models from infrared (IR) to submillimeter (sub-mm) wavelength (see also Zhang & Tan 2011;Zhang et al. 2013b;, as their paper series). In our first paper (Tanaka et al. 2016, hereafter Paper I), we calculated the evolutionary sequence of photoionized structures using the infall-and-outflow structures developed by Zhang et al. (2014). We showed that at early stages the photoionized region is tightly confined by the dense wall of the inner outflow. Later the photoionized region breaks out to fill the entire outflow cavity in a timescale of 10 3 -10 4 yr. The predicted free-free flux and the width of Hydrogen recombination lines from these models are consistent with the trend of observed radio jets. However, our model grid in Paper I has the protostellar-mass resolution of ∆m * ≥ 4 M ⊙ which is too coarse to follow the break-out phase of the photoionized region in detail. The break-out phase is important as the first signature of the evolution of the photoionized region, i.e., the first observational signature that a massive star is actually forming.
In this second paper, we thus increase the protostellarmass resolution of our gird model to ∆m * = 1 M ⊙ around the break-out phase. We also upgrade the outflow density structure from the analytic model by Zhang et al. (2014) to the simulated model based on the high-resolution three-dimensional magnetohydrodynamics (MHD) code by Staff et al. (2015). An additional improvement from Paper I is that we discuss not only the longterm evolution over 10 4 yr of the photoionized region and its free-free emission, but also their potential short-timescale variability over 10-100 yr. Recent observations of massive protostars have reported short-timescale variability observed in radio free-free continuum emission (Carrasco-González et al. 2015) and in thermal continuum emission from dust (Caratti o Garatti et al. 2016;Hunter et al. 2017). Also theoretical simulations have suggested accretion bursts may occur via disk instabilities in massive star formation (Krumholz et al. 2009;Kuiper et al. 2011;Meyer et al. 2017;Matsushita et al. 2017). The photoionized structure is sensitive to the outflow density profile, which is expected to scale with the mass accretion rate, e.g., in magnetocentrifugallydriven disk wind models. Therefore, we investigate the potential of the short-timescale variability of the photoionized structure and its free-free emission as a response to variable accretion rates.

Protostellar models
We calculate protostellar evolution based on the semianalytic method of Tanaka et al. (2017), which selfconsistently included mass accretion regulated by multiple feedback processes (see also the series of papers by Zhang & Tan 2011;Zhang et al. 2013bZhang et al. , 2014. This model assumes a prestellar core collapses to form one massive star. The initial cloud core is assumed to be spherical and in quasi-virial equilibrium, supported by turbulence and/or magnetic fields (McKee & Tan 2003). The core properties are determined by three main parameters; core mass M c , mass surface density of the ambient clump Σ cl , and the core's initial rotational to gravitational energy ratio β c . In this paper, we set fiducial core mass and surface density to M c = 60 M ⊙ and Σ cl = 1 g cm −2 (McKee & Tan 2003). The core radius is then determined to be 0.057 pc. Based on observations of both low-mass and high-mass cores, β c is found to have a typical value of ∼ 0.02 (e.g., Goodman et al. 1993;Li et al. 2012;Palau et al. 2013). Therefore, we fix the fiducial value of 0.02 for β c .
The cloud core undergoes the inside-out collapse which is treated by a self-similar solution (McLaughlin & Pudritz 1997). The accretion rate onto the central starṁ * is shown in panel (a) of Figure 1. The typical value ofṁ * is ∼ 10 −4 M ⊙ yr −1 . The accretion is regulated by multiple feedback processes, i.e., MHD disk winds, radiation pressure, photoevaporation, and stellar winds (Tanaka et al. 2017). Eventually, the main accretion phase finishes at a stellar mass of 24M ⊙ , i.e., the star formation efficiency from the initial core is 24M ⊙ /60M ⊙ = 0.4.
The protostellar evolution is obtained based on the model of Hosokawa & Omukai (2009) and Hosokawa et al. (2010). This model solves for the internal structure of the protostar, assumed to be spherical, including effects such as hydrogen and deuterium burning, convection, and radiative transfer. The evolution of the protostar depends on the accretion history, which is self-consistently calculated as explained above. The evolution of the accretion rate and protostellar properties are shown in panels (b) -(e) in Figure 1. The stellar radius, r * , first increases to be as large as 45 R ⊙ when m * = 5M ⊙ , due to the relatively high mass-accretion rate. Then the protostar undergoes Kelvin-Helmholtz (KH) contraction, reaching the zeroage main sequence (ZAMS) phase at around 20 M ⊙ . The accretion luminosity, L acc , dominates before KH contraction sets in, while the internal stellar luminosity, L * , becomes dominant later. The total luminosity, L * ,acc , increases almost monotonically. The effective temperature, T * ,acc , has a local minimum when m * = 5 M ⊙ because of the large stellar radius. By the combination of variations of L * ,acc and T * ,acc , the ionizing-photon rate, S * ,acc , drops down to 10 33 s −1 at around 5M ⊙ . It then rises dramatically to over 10 48 s −1 during KH contraction. This leads to the formation and expansion of a photoionized region in the KH contraction phase (Paper I). To study the earliest phase of photoionization, the protostellar models in the KH phase at stellar masses of m * = 10, 11, and 12 M ⊙ are adopted for the following photoionization calculations, which are named as models M10, M11, and M12, respectively. The time duration from m * = 10 to 12M ⊙ is 6, 800 yr in this evolutionary calculation.

Outflow density profile
The outflow density profile is needed to calculate the photoionized structure. In Paper I, we adopted the semianalytic disk-wind model, which is modified from the magnetocentrifugally-driven outflow model of Blandford & Payne (1982) (Zhang et al. 2013b). The mass-loading rate of the wind is proportional to the mass accretion rate with a ratio f w = 0.1 which is a typical value for disk winds (Königl & Pudritz 2000). Though our analytic disk-wind model appears to have some agreement with observed sources (Zhang et al. 2013a, Paper I), it (Left panels a -e:) Evolution of the accretion rate (ṁ * ), protostellar radius (r * ), total luminosity including that from accretion (L * ,acc), effective surface temperature (T * ,acc), and total H-ionizing photon production rate (S * ,acc), as functions of protostellar mass. The properties of the zero-age main sequence (ZAMS) are shown by dashed lines (Table 1 of Davies et al. 2011). In panel (c), the accretion and intrinsic luminosities (Lacc and L * ) are also shown by thin green and orange lines, respectively. The vertical dotted lines indicate the stages at m * = 10, 11, and 12 M ⊙ , which are inputs for post-processing photoionization calculations. The protostellar evolution is calculated by the method of Tanaka et al. (2017) (see text). (Right panel f:) Outflow density structure (in a thin slice) at the end of the simulation by Staff et al. (2015). Here the density structure at a protostellar mass of 12 M ⊙ is shown, however the structure is also scalable to other cases (see §2.2). Mirror symmetry is assumed about the disk mid-plane, which is shown by the white dotted line.
is highly simplified and idealized, e.g., assuming an axisymmetric smooth structure. In this study, therefore, we apply the density profile from the three-dimensional MHD simulation of disk winds by Staff et al. (2015) for a post-processing photoionization calculation. Staff et al. (2015) studied disk winds from a disk using the ZEUSMP code (Norman 2000). This is a threedimensional ideal MHD simulation with a Keplerian accretion disk as a fixed boundary condition, along with a magnetic field that is initially poloidal and in a selfsimilar configuration. The wind launched from the disk surface is magneto-centrifugally accelerated and thus creates the outflow. Although the original goal of Staff et al. (2015) was the simulation of jets from low-mass protostars, this simulation is also fully scalable to the outflows around massive protostars. This is because they solved the dimensionless form of the equations, and thus the results can be re-scaled based on following fundamental protostar properties (more details please see Ouyed & Pudritz 1997). The unit length (and also the highest resolution) is set by the the innermost radius of the accretion disk. The innermost radius may be larger than the stellar radius if the stellar magnetosphere is strong enough to truncate the disk (Ghosh & Lamb 1978;Shu et al. 2000). However, in the massaccretion phase which we are interested in, the accretion disk may continue to the surface of the star, because the kinematic and gaseous pressures are higher than the magnetic pressure due to the high accretion rate. Therefore, we chose the disk innermost-radius as the stellar radius, r * . The time unit is the inverse of the Keplerian angular velocity at the length unit r * , i.e., 1/Ω K * ≃ 0.18 (m * /10M ⊙ ) −1/2 (r * /10R ⊙ ) 3/2 . The density unit is set to match the simulated disk-wind rate to the analytic value of f wṁ * .

Models for fiducial accretion rate
While Staff et al. (2015) have investigated several initial magnetic field configurations, we use the result of the "BP" configuration, which is consistent to the model of Blandford & Payne (1982), i.e., B p ∝ ̟ −1.25 , where B p is the initial poloidal magnetic field strength at the disk surface and ̟ is a cylindrical radius. Panel (f) of Figure 1 shows a density slice for the outflow from our m * = 12 M ⊙ protostar model. Mirror symmetry at the disk mid-plane is assumed. The density profile is quite close to being axisymmetric, but with smaller nonaxisymmetric structures, which could not be treated in the previous semi-analytic model we used in Paper I. Note that the geometric thickness of the disk is assumed to be zero in this simulation, so the model applies to the region of outflow cavity away from the disk region.
We note a caveat on the usage of this simulated outflow. As a result of its high spatial resolution, the simulated time is limited to be shorter than the expected accretion time. Therefore, the size of the simulated outflow is only about 100AU (panel (f) in Fig. 1), while, in reality, the outflow is expected to be larger than the core scale of 0.1 pc. On the other hand, thanks to this high resolution, we are able to model structures down to the vicinity of the protostellar radius where the stellar radiation is mainly processed. Thus, in this paper, we focus only on the earliest break-out phase of the photoionized region while it is confined on scales of ∼ 100 AU or smaller. In this earliest phase, the rest of the outer outflow would be essentially neutral and would not contribute significantly to the free-free emission due to photoionization.

Models for variable accretion bursts
In this paper, we also address the potential shorttimescale variability induced by accretion bursts. The accretion rate in our evolutionary model described in §2.1 changes smoothly on a timescale of 10 4 -10 5 yr. However, some numerical simulations of massive star formation predict that the accretion rate can be highly timevariable due to disk self-gravity instabilities (Krumholz et al. 2009;Kuiper et al. 2011;Meyer et al. 2017;Matsushita et al. 2017). In the episodic accretion scenario, the protostars spend most of their time in a low accretion-rate, quiescent phase, with short-timescale (∼tens of years) accretion bursts in which accretion rates are one or two orders of magnitude higher (Vorobyov & Basu 2015).
The outflow is powered by accretion and thus its density is enhanced in the accretion burst phase (e.g., Tomida et al. 2017;Matsushita et al. 2017). Especially, on the small scales of 100 AU, the outflow density would quickly respond to the accretion variation since the crossing timescale of outflow is only 10 yr. Here we have used the outflow velocity of ∼ 250-600 km s −1 from the MHD simulation of Staff et al. (2015).
To investigate the variation of the photoionized structure induced by episodic accretion, we mimic the burst/quiescent phases by enhancing/reducing the outflow-density by a factor of 3.16 from the standard value for the 11M ⊙ model, which we refer to as M11H and M11L, respectively. The density variation between models M11H and M11L is one order of magnitude, i.e., 3.16 2 ≃ 10. Note that the accretion bursts may in general induce density changes that are larger or smaller than this factor. In such cases, the resulting variation in radio flux would be greater or smaller than those shown in §3.2, respectively. We use the same protostellar model as the standard accretion case, although the episodic accretion may also affect the protostellar properties. We discuss this point further in §4.

Photoionizing calculation
We calculate the photoionized structure using the protostellar properties from the stellar evolution calculation ( §2.1) and the outflow density profile from the MHD simulation ( §2.2). In Paper I, under the assumption of axisymmetry, the full transfer of ionizing photons was solved, which allows the accurate treatment of both the direct and diffuse radiation fields. As the density profile from the 3-D simulation is used in this study, it is computationally much more expensive to perform the full radiative transfer calculation. Thus we simplify the method of the photoionizing calculation. Tanaka et al. (2013) showed that the direct stellar radiation is more important than the diffuse radiation from recombination to the ground state. In this study, we use the on-the-spot approximation to treat the diffuse radiation field, i.e., the diffuse radiation is assumed to be absorbed at the point of emission. Then, the photoionized radiation field can be obtained by where F EUV is the radial number flux of ionizing photons, r is the distance from the center of the star, n H is the number density of hydrogen nuclei, x II is the ionization fraction of H, and σ H is the mean photoionization crosssection of the hydrogen atom. The transfer equation (1) can be solved as, where τ EUV is the optical depth for the ionizing radiation. The ionization fraction x II is obtained from the balance between photoionization and recombination, where α B (T ) is the recombination coefficient to all states except the ground state (so-called case B) which depends on the local gas temperature T . The temperature of photoionized gas is usually regulated around 10, 000 K. We set the temperature as a function of density and the stellar spectrum, based on parameterized spherical calculations using CLOUDY (Ferland et al. 2013) (see §2.2.3 of Paper I). The mean values of the hydrogen cross-section, σ H , is evaluated from the stellar spectra for each model. Solving equations (2) -(4) over the entire solid angle, the photoionized structure is obtained. In this paper, the dust scattering and absorption of ionizing photons are ignored which was included in Paper I. We expect this approximation to be valid, since here we focus on the innermost region of < 100 AU, where the dust grains have been destroyed, i.e., if most of the volume is filled by outflows launched from the region of the disk inside the dust destruction front.

Free-free emission
We calculate the free-free radio continuum emission from the obtained photoionized outflow structures. The method for the free-free transfer calculation is the same as Paper I, but now applied to the 3-D calculation of this study. The equation for the free-free radiation transfer is, I ν,ff is the free-free intensity, κ ν,ff (T, n H ) is the opacity of free-free emission, and B ν (T ) is the Planck function. The observed flux density is calculated by integrating over source solid angle of Ω s , Here the far field limit is assumed. See also §2.3.1 of Paper I for more details of the radiative transfer of freefree emission.

RESULTS
We here describe the results of photoionizing and freefree emission calculations. First, we show the earliest break-out evolution of the photoionized region which occurs over a time of about 10 4 yr. Then, we explain the variation of the photoionized structure due to accretion bursts, which could occur in tens of years.

Earliest break-out of Outflow-Confined H II Regions
The evolution from m * = 10M ⊙ to 12M ⊙ is shown in Figure 2 (models M10, M11, and M12). In the top panels, the outflow densities of different models are similar, since they are simply scaled based on the protostellar parameters as described §2.2. In the middle panels, the displayed temperature structure is weighted by the square of the ionization fraction, i.e., x 2 II T , because we are only interested in the photoionized region and do not solve for the temperature of the neutral region in this study. Recall that, in the photoionization calculation, most of the region is either fully-ionized (1 − x II ≪ 1) or neutral (x II ≪ 1). The photoionized region, where the temperature is about 10, 000 K, is confined on 100 AU scales by the high-density outflow, i.e., an "outflow-confined" H II region (Tan & McKee 2003). Especially, the size of the photoionized region is only about 10 AU when m * = 10 M ⊙ . Then, the photoionized region elongates near the outflow axis at 11 M ⊙ , since the density is lower in this direction than that near the equatorial plane. The photoionized region starts to expand widely at 12 M ⊙ . The shape of the photoionized region is very spiky. This is because of the shadows cast by the high-density clumpy structures.
The bottom panels of Figure 2 show the free-free emission maps of models M10, M11, and M12 at 10GHz. The viewing angle from the outflow axis is assumed to be 60 • . The spiky features in the temperature slices are now less apparent because of their overlap along the line of sight. Figure 2. The earliest evolution of the outflow-confined H II region for models M10, M11, and M12 (from left to right), shown by density slices (top), ionized-fraction-weighted temperature slices (middle), and the resolved free-free emission images at 10 GHz (bottom). The scale-length is shown in each panel. In the bottom panels, the intensity is shown by the brightness temperature and a viewing angle of 60 • is assumed.
However, the bipolar shape and the elongated structure are still visible, especially in model M11. Figure 3(a) shows the SEDs for models M10, M11, and M12. A distance of 1 kpc and a viewing angle of 60 • are assumed. The dust thermal emission is evaluated based on the model of Zhang et al. (2014), which typically has a spectrum index of 3.3 at radio frequencies. For comparison, the SED of a later stage model with m * = 16 M ⊙ from Paper I is also shown in Figure 3(a) (model A16 in Paper I). Combining our new results of M10, M11, and M12 and the previous 16M ⊙ model, the free-free emission evolves in the following manner. At the moment when the photoionized region is formed (M10), the freefree flux is lower than the dust flux at most frequencies.
Then, the ionized region expands to 100 AU scales in about 10 4 yr, and the radio free-free flux starts to domi- Figure 3. (a) SEDs for models M10, M11, and M12, together with the case of m * = 16M ⊙ from Paper I (orange, blue, green, and gray, respectively). (b) For comparison, SEDs for m * = 10, 11, 12, and 16 M ⊙ based on the analytic outflow density model calculated as in Paper I (orange, blue, green, and gray, respectively). In both panels, the thick dashed lines represent free-free emission, while the thick solid lines show the total emission including dust emission. The black thin-dashed lines are plotted for illustrative purposes, to show spectral indices of 0.6 and 2. Note that the dust thermal emission is obtained from , and a distance of 1 kpc and the viewing angle of 60 • are assumed in all models. The free-free emission for the 10M ⊙ model in panel (b) is not seen because it is weaker than 10 −3 mJy. nate the dust flux at around 10 GHz or lower frequency. (M11 and M12). At this stage, the free-free flux has a spectral index of about 2. Later, the photoionized region breaks out beyond the 1000 AU scale, most of the outflow cavity is ionized, and the free-free flux density becomes as high as 100 mJy at 10 GHz with a spectral index of 0.6 (models of Paper I).
While the infrared emission at the peak of the SED increases by a factor of two from the case with a protostellar mass of 10M ⊙ to that of 16M ⊙ , the flux at 10GHz jumps by three orders of magnitude because of the sharp rise of the ionizing photon rate due to the KH contraction. The spectral index of near 2 in the earlier stage indicates the ionized region is mostly optically thick at these wavelengths. This is because the density of the inner 100 AU-region is very high. On the other hand, the shallower index of ∼ 0.6 in the later stage represents partially optically thin conditions, since the outer region has lower density (for more details about the partially optically thin case, see §4.1 of Paper I and a classic model by Reynolds 1986).
In Figure 3(b), for comparison, we present the SEDs calculated using the analytic outflow density structures from Paper I. The basic evolution of free-free emission is the same as that in the simulated outflow case: the flux rises up by orders of magnitude and the spectral index changes from ∼ 2 to ∼ 0.6 as the protostellar mass increases from m * = 10 to 16M ⊙ . The main difference in free-free evolution of the simulated and analytic outflow models is the degree of increase from 10 to 12M ⊙ at radio wavelength. In the case of the simulated outflow model, the free-free flux changes about two orders of magnitude from m * = 10 to 12M ⊙ at ∼ 10 GHz. On the other hand, in the case of the analytic outflow model, the free-free flux is weaker than ∼ µJy at 10M ⊙ and then increases to the ∼ mJy level at 12M ⊙ . This is because the analytic outflow has a steeper density gradient near the outflow axis than the simulated outflow and thus the photoionized region breaks out more quickly. The difference of break-out timing also causes differences in the spectral indices at ν < 10 GHz in the case of m * = 12M ⊙ : the index is close to 2 in the analytic model since the photoionized region is confined within a 100 AU scale, where it is optically thick at this frequency, while the index of the simulated model is close to 0.6 because the photoionized region is extended to include optically thinner outer regions. In this way, the details of photoionized-region break-out are sensitive to the outflow structure, while the general trends of the evolution are not.

Variations Induced by Accretion Bursts
Next, we examine the potential variability induced by accretion bursts. Model M11L represents the quiescent phase which has lower outflow density by a factor of 0.316 than the standard case of M11, and model M11H represents the burst phase having the higher density by a factor of 3.16 than the standard case.
The left panels in Figure 4 show the radio free-free . The intensity is shown as brightness temperature for emission at a frequency of 10 GHz. The linear scale is same in all panels, with sources assumed to be at 1 kpc distance and a viewing angle of 60 • from the outflow axis. (Right:) SEDs for models M11L, M11, and M11H (orange, blue, green, and gray, respectively), again for a distance of 1 kpc and viewing angle of 60 • to the outflow axis. The dashed lines show free-free emission, and the solid lines represent the total emission including dust emission. The thermal dust emission is obtained from . Note that the radio flux at frequencies lower than 10 GHz decreases by an order of magnitude from a higher value in the low accretion phase to a lower value during the accretion burst, while the peak of the thermal dust emission changes by a much smaller factor.
continuum images of M11L, M11, and M11H at 10 GHz. In the quiescent model M11L, the photoionized region extends to larger distances than the standard model M11 because the outflow density is lower. However, due to the low density, the outer part of the photoionized region becomes optically thin and the free-free emission from there is fainter than the optically thick value of T b ≃ 10 4 K. On the other hand, in the accretion-burst model M11H, the photoionized region is confined to within ∼ 10 AU by the high density outflow and is optically thick at this frequency. In this way, the quiescent phase has a larger size and partially optically-thin emission, while the burst phase shows smaller and optically-thick emission. These features are also seen in the SEDs shown in the right panel of Figure 4. Model M11L has about one order of magnitude higher free-free flux at 10 GHz than model M11H, because it has a larger photoionized region. On the other hand, the spectral index at 10 GHz of M11L is ∼ 1, indicating it is partially optically thick, and is shallower than M11H's index of 2. The free-free flux from outflow-confined H II regions is thus sensitive to the timevariability of the outflow density. Also note that since the internal luminosity from the protostar dominates over accretion luminosity after KH contraction starts, the flux from thermal dust emission does not change as significantly due to such accretion bursts.

DISCUSSIONS AND CONCLUSIONS
Here we discuss the observational predictions based on our results of the earliest break-out phase of outflowconfined H II regions and their time-variability due to variable accretion.

Break-out phase of photoionized regions
The photoionized region starts to form as KH contraction proceeds. In the model of this paper, the onset of photoevaporation is at a stellar mass of about 11M ⊙ ( §3). However, we note that the onset mass depends on the accretion history because a higher accretion rate leads to a higher mass at which KH contraction starts, and also causes the density of outflow to be higher (Hosokawa & Omukai 2009;Zhang et al. 2014). For typical accretion rates of ∼ 10 −4 -10 −3 M ⊙ yr −1 , the mass of the onset of photoionization is ∼ 10-20M ⊙ (see Paper I). The photoionized region is first confined by the dense MHD outflow wall at < 100 AU scales. At this stage, the typical radio luminosity is ∼ mJy kpc 2 , with a spectral index of ∼ 2, i.e., optically thick conditions. Ilee et al. (2016) reported a candidate of this type of hypercompact H II region around the young massive star G11.92−0.61 MM1. The radio flux increases as the photoionized region expands, and almost the entire outflow is ionized about 10 4 yr after its formation. At this stage, the free-free luminosity becomes ∼ 100 mJy kpc 2 , with a spectral index of ∼ 0.6 at around 10 GHz which is similar to observed radio jets around massive protostars, e.g., the central radio source of IRAS16547−4247 (Rodríguez et al. 2008). See Paper I for more discussion about this later stage.
Considering the highest resolution of the Very Large Array (VLA) of 0.2 ′′ at 10 GHz or 0.043 ′′ at 45 GHz and typical distances to massive protostars of several kpc, it is usually not possible to observationally resolve the photoionized region when it is confined within 100 AU scales. However, this becomes possible for sources closer than 1 kpc, like Source I in the Orion KL region, which is 0.4 kpc away ). In Figure 5, we show the image of M11 at 45 GHz, convolving by a beam of 0.05 ′′ FWHM diameter, and assuming a distance of 0.5 kpc. While elongation along the outflow axis can be resolved, in the perpendicular direction the source is unresolved since the photoionized region has a very long, thin shape. This causes the peak brightness temperature to be 2, 000 K, which is lower than the ionized gas temperature of ∼ 10, 000 K, even though it is almost completely optically thick. Such a model may be considered as an alternative explanation for the low brightness temperature, ∼ 1, 500 K, of Source I in Orion KL, rather than being due to H − free-free emission from a collisionally ionized disk Hirota et al. 2015). However, we note that, while the spectral index of Source I is 1.86±0.26 which is consistent with our model, its radio luminosity of 0.2 mJy kpc 2 is smaller than our model prediction (Forbrich et al. 2016). This low-radio luminosity with the optically-thick spectral-index might be explained by the outflow-confined H II model with a lower accretion rate. However, we have studied only one initial condition in this paper. For more detailed discussion of this photoionized outflow hypothesis to explain Source I including non-detection of hydrogen recombination lines at the frequencies of 353.6 and 662.4 GHz (Plambeck & Wright 2016), we need to expand the range of protostellar source model parameters, which we defer to a future paper.

Short-time variabilities induced by accretion
variation Apart from the break-out on ∼ 10 3 yr timescales, the photoionized structure can also change on shorter timescales due to accretion rate variability. Accretion bursts on ∼ 10-100 yr timesscales induced by disk instability may be expected during massive star formation (e.g., Krumholz et al. 2009;Kuiper et al. 2011;Meyer et al. 2017;Matsushita et al. 2017). Since the outflow is powered by accretion, the outflow density is higher (lower) in the accretion burst (quiescent) phase (see simulations by Tomida et al. 2017;Matsushita et al. 2017). Due to this density variation, the photoionized region is smaller in the burst phase than that in the quiescent phase ( §3.2). However, there are as yet no clear observational examples that fully support this theoretical expectation. Carrasco-González et al. (2015) reported that, around the massive protostar W75N(B)-VLA 2, the radio continuum emission changed from a compact to an elongated morphology within a period of 18 years. Although the timescale is consistent with those that may result from accretion and outflow variability, other observational aspects are not particularly supportive of this being an example of photoionized-structure variation induced by accretion/outflow fluctuations. In particular, water masers, which are expected to trace shocked regions, also change their morphology systematically with the radio continuum, potentially indicating that this radio flux is induced from shocked ionization rather than photoionization. However, following the recently reported dust-reprocessed luminosity bursts around massive protostars (Caratti o Garatti et al. 2016;Hunter et al. 2017), which are expected to be due to accretion bursts, we predict that the radio flux variations induced by the accretion/outflow fluctuations affecting inner photoionized regions will eventually be observed.
Note that the free-free flux from the photoionized region is smaller in the accretion burst phase, which is the opposite behavior to that of the dust luminosity. Additionally our model shows that, when the accretion rate changes by one order of magnitude, the free-free emission also varies by one order of magnitude, while the peak of the thermal dust emission varies much less (Fig. 4). This is because the total luminosity is dominated by the internal protostellar luminosity after KH contraction starts (see Fig. 1). The accretion luminosity is of greater relative importance in the earlier protostellar phase. Therefore, massive protostars associated with the dust-emitted luminosity bursts (Caratti o Garatti et al. 2016;Hunter et al. 2017) may be more likely to be in the pre-KH contraction phase. We thus suggest that radio free-free observations will be a better method for tracing accretion variability in the later, more evolved stages.
Finally, we note the caveat that we used the same protostellar parameters for the different accretion rate models in §3.2, although protostellar evolution is expected to be related to accretion history. Smith et al. (2012) showed that, in the case of primordial star formation, protostellar evolution under the influence of a variable accretion rate is consistent with the evolution calculated with accretion rates averaged over the last 100 yr. This is because the timescale of protostellar evolution is described by the averaged accretion rate over the local KH time, t KH = Gm 2 * /(r * L * acc ). The KH timescale of the protostar considered in this study is about 5, 000 yr, which is longer than the typical duration of accretion bursts, and thus we expect the accretion burst does not alter these protostellar properties significantly.