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

TL;DR

This paper models the early phase of ionized region breakout in massive star formation, showing how the region rapidly expands and brightens, with variability driven by accretion bursts detectable over decades.

Contribution

It introduces a detailed simulation of the breakout phase using 3D MHD and radiative transfer, highlighting the rapid expansion and observational signatures of the ionized region.

Findings

Ionized region appears at 10Msun and confined within 10-100AU.

Breakout occurs in less than 10,000 years, significantly increasing radio brightness.

Radio spectral index shifts from 2 to 0.6 during breakout.

Abstract

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.

Equations

$$\frac{1}{r^2}\frac{d}{dr}\left(r^2{F}_{\mathrm{EUV}}\right)=-(1-x_{\mathrm{II}})n_{\mathrm{H}}{\sigma }_{\mathrm{H}}{F}_{\mathrm{EUV}},$$

$$F_{\mathrm{EUV}}=\frac{S_{\ast ,\mathrm{acc}}{e}^{-{\text{uptau}}_{\mathrm{EUV}}}}{4\pi r^2},$$

$${\text{uptau}}_{\mathrm{EUV}}(r)={\int }_{r_{\ast }}^r(1-x_{\mathrm{II}})n_{\mathrm{H}}{\sigma }_{\mathrm{H}}d{r}^{\prime },$$

$$(1-x_{\mathrm{II}})n_{\mathrm{H}}{\sigma }_{\mathrm{H}}{F}_{\mathrm{EUV}}={\alpha }_{\mathrm{B}}(T)x_{\mathrm{II}}^2{n}_{\mathrm{H}}^2,$$

$$\frac{d{I}_{\nu ,\mathrm{ff}}}{ds}=-{\kappa }_{\nu ,\mathrm{ff}}\left(I_{\nu ,\mathrm{ff}}-B_{\nu }\right)$$

$$F_{\nu ,\mathrm{ff}}={\int }_{{\Omega }_s}{I}_{\nu ,\mathrm{ff}}({\theta }_{\mathrm{view}})d\Omega .$$

Full text

1. 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 (Tan et al., 2014). 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\rm\>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 ($\leq 0.1\rm\>pc$) and, especially, hyper-compact ($\leq 0.01\rm\>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 star-forming 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 observational signatures of the Turbulent Core model of McKee & Tan (2003). In the same context, Zhang et al. (2014) and (Zhang & Tan, 2017) have developed a suite of radiative transfer models from infrared (IR) to sub-millimeter (sub-mm) wavelength (see also Zhang & Tan, 2011; Zhang et al., 2013b; Zhang & Tan, 2017, 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}\rm\>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 $\Delta m_{*}\geq 4\>M_{\odot}$ 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 protostellar-mass resolution of our gird model to $\Delta m_{*}=1\>M_{\odot}$ 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 magneto-hydrodynamics (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}\rm\>yr$ of the photoionized region and its free-free emission, but also their potential short-timescale variability over $10$–$100\rm\>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 magnetocentrifugally-driven 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.

2. Methods

2.1. Protostellar models

We calculate protostellar evolution based on the semi-analytic method of Tanaka et al. (2017), which self-consistently included mass accretion regulated by multiple feedback processes (see also the series of papers by Zhang & Tan, 2011; Zhang et al., 2013b, 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 $\Sigma_{\rm cl}$, and the core’s initial rotational to gravitational energy ratio $\beta_{c}$. In this paper, we set fiducial core mass and surface density to $M_{c}=60\>M_{\odot}$ and $\Sigma_{\rm cl}=1\>{\rm g\>cm^{-2}}$ (McKee & Tan, 2003). The core radius is then determined to be $0.057{\rm\>pc}$. Based on observations of both low-mass and high-mass cores, $\beta_{c}$ is found to have a typical value of $\sim 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 $\beta_{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 $\dot{m}_{*}$ is shown in panel (a) of Figure 1. The typical value of $\dot{m}_{*}$ is $\sim 10^{-4}\>M_{\odot}\>{\rm 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_{\odot}$, i.e., the star formation efficiency from the initial core is $24M_{\odot}/60M_{\odot}=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_{\odot}$ when $m_{*}=5M_{\odot}$, due to the relatively high mass-accretion rate. Then the protostar undergoes Kelvin-Helmholtz (KH) contraction, reaching the zero-age main sequence (ZAMS) phase at around $20\>M_{\odot}$. The accretion luminosity, $L_{\rm acc}$, dominates before KH contraction sets in, while the internal stellar luminosity, $L_{*}$, becomes dominant later. The total luminosity, $L_{\rm*,acc}$, increases almost monotonically. The effective temperature, $T_{\rm*,acc}$, has a local minimum when $m_{*}=5\>M_{\odot}$ because of the large stellar radius. By the combination of variations of $L_{\rm*,acc}$ and $T_{\rm*,acc}$, the ionizing-photon rate, $S_{\rm*,acc}$, drops down to $10^{33}\>{\rm s^{-1}}$ at around $5M_{\odot}$. It then rises dramatically to over $10^{48}\>{\rm 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_{\odot}$ 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_{\odot}$ is $6,800{\rm\>yr}$ in this evolutionary calculation.

2.2. Outflow density profile

The outflow density profile is needed to calculate the photoionized structure. In Paper I, we adopted the semi-analytic 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_{\rm 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 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.

2.2.1 Models for fiducial accretion rate

Staff et al. (2015) studied disk winds from a disk using the ZEUSMP code (Norman, 2000). This is a three-dimensional 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 self-similar 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 mass-accretion 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/\Omega_{\rm K*}\simeq 0.18\left(m_{*}/10M_{\odot}\right)^{-1/2}\left(r_{*}/10R_{\odot}\right)^{3/2}$. The density unit is set to match the simulated disk-wind rate to the analytic value of $f_{\rm w}\dot{m}_{*}$.

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_{\rm p}\propto\varpi^{-1.25}$, where $B_{\rm p}$ is the initial poloidal magnetic field strength at the disk surface and $\varpi$ is a cylindrical radius. Panel (f) of Figure 1 shows a density slice for the outflow from our $m_{*}=12\>M_{\odot}$ protostar model. Mirror symmetry at the disk mid-plane is assumed. The density profile is quite close to being axisymmetric, but with smaller non-axisymmetric 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 $100{\rm AU}$ (panel (f) in Fig. 1), while, in reality, the outflow is expected to be larger than the core scale of $0.1\>{\rm\>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 $\sim 100{\rm\>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.

2.2.2 Models for variable accretion bursts

In this paper, we also address the potential short-timescale 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}{\rm\>yr}$. However, some numerical simulations of massive star formation predict that the accretion rate can be highly time-variable 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 ($\sim$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{\rm\>AU}$, the outflow density would quickly respond to the accretion variation since the crossing timescale of outflow is only $\lesssim 10{\rm\>yr}$. Here we have used the outflow velocity of $\sim 250$–$600{\rm\>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_{\odot}$ 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}\simeq 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.

2.3. 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

[TABLE]

where $F_{\rm EUV}$ is the radial number flux of ionizing photons, $r$ is the distance from the center of the star, $n_{\rm H}$ is the number density of hydrogen nuclei, $x_{\rm II}$ is the ionization fraction of H, and $\sigma_{\rm H}$ is the mean photoionization cross-section of the hydrogen atom. The transfer equation (1) can be solved as,

[TABLE]

where $\uptau_{\rm EUV}$ is the optical depth for the ionizing radiation. The ionization fraction $x_{\rm II}$ is obtained from the balance between photoionization and recombination,

[TABLE]

where $\alpha_{\rm 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{\rm\>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, $\sigma_{\rm 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\rm\>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.

2.4. 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,

[TABLE]

$I_{\nu,{\rm ff}}$ is the free-free intensity, $\kappa_{\nu,{\rm ff}}(T,n_{\rm H})$ is the opacity of free-free emission, and $B_{\nu}(T)$ is the Planck function. The observed flux density is calculated by integrating over source solid angle of $\Omega_{s}$,

[TABLE]

Here the far field limit is assumed. See also §2.3.1 of Paper I for more details of the radiative transfer of free-free emission.

3. Results

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

3.1. Earliest break-out of Outflow-Confined H II Regions

The evolution from $m_{*}=10M_{\odot}$ to $12M_{\odot}$ 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_{\rm II}^{2}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_{\rm II}\ll 1$) or neutral ($x_{\rm II}\ll 1$). The photoionized region, where the temperature is about $10,000{\rm\>K}$, is confined on $100{\rm\>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{\rm\>AU}$ when $m_{*}=10\>M_{\odot}$. Then, the photoionized region elongates near the outflow axis at $11\>M_{\odot}$, since the density is lower in this direction than that near the equatorial plane. The photoionized region starts to expand widely at $12\>M_{\odot}$. 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 $10\rm\>GHz$. The viewing angle from the outflow axis is assumed to be $60\arcdeg$. The spiky features in the temperature slices are now less apparent because of their overlap along the line of sight. 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\rm\>kpc$ and a viewing angle of $60\>\arcdeg$ 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_{\odot}$ 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_{\odot}$ model, the free-free emission evolves in the following manner. At the moment when the photoionized region is formed (M10), the free-free flux is lower than the dust flux at most frequencies. Then, the ionized region expands to $100\rm\>AU$ scales in about $10^{4}\rm\>yr$, and the radio free-free flux starts to dominate the dust flux at around $10\rm\>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{\rm\>AU}$ scale, most of the outflow cavity is ionized, and the free-free flux density becomes as high as $100\rm\>mJy$ at $10\rm\>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 $10\>M_{\odot}$ to that of $16\>M_{\odot}$, the flux at $10\rm\>GHz$ 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\rm\>AU$-region is very high. On the other hand, the shallower index of $\sim 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 $\sim 2$ to $\sim 0.6$ as the protostellar mass increases from $m_{*}=10$ to $16M_{\odot}$. The main difference in free-free evolution of the simulated and analytic outflow models is the degree of increase from $10$ to $12M_{\odot}$ 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 $12\>M_{\odot}$ at $\sim 10{\rm\>GHz}$. On the other hand, in the case of the analytic outflow model, the free-free flux is weaker than $\sim{\rm\mu Jy}$ at $10M_{\odot}$ and then increases to the $\sim{\rm mJy}$ level at $12\>M_{\odot}$. 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 $\nu<10{\rm\>GHz}$ in the case of $m_{*}=12M_{\odot}$: the index is close to 2 in the analytic model since the photoionized region is confined within a $100\rm\>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.

3.2. 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 continuum images of M11L, M11, and M11H at $10{\rm\>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_{\rm b}\simeq 10^{4}\rm\>K$. On the other hand, in the accretion-burst model M11H, the photoionized region is confined to within $\sim 10\rm\>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\rm\>GHz$ than model M11H, because it has a larger photoionized region. On the other hand, the spectral index at $10\rm\>GHz$ of M11L is $\sim 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 time-variability 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.

4. Discussions and Conclusions

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

4.1. 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 $11\>M_{\odot}$ (§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 $\sim 10^{-4}$–$10^{-3}\>M_{\odot}\>{\rm yr}^{-1}$, the mass of the onset of photoionization is $\sim 10$–$20M_{\odot}$ (see Paper I). The photoionized region is first confined by the dense MHD outflow wall at $<100{\rm\>AU}$ scales. At this stage, the typical radio luminosity is $\sim{\rm\>mJy\>kpc^{2}}$, with a spectral index of $\sim 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}\rm\>yr$ after its formation. At this stage, the free-free luminosity becomes $\sim 100\rm\>mJy\>kpc^{2}$, with a spectral index of $\sim 0.6$ at around $10\rm\>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\arcsec$ at $10{\rm\>GHz}$ or $0.043\arcsec$ at $45{\rm\>GHz}$ and typical distances to massive protostars of several $\rm kpc$, it is usually not possible to observationally resolve the photoionized region when it is confined within $\lesssim 100\rm\>AU$ scales. However, this becomes possible for sources closer than $1\rm\>kpc$, like Source I in the Orion KL region, which is $0.4\rm\>kpc$ away (Menten et al., 2007). In Figure 5, we show the image of M11 at $45\rm\>GHz$, convolving by a beam of $0.05\>\arcsec$ FWHM diameter, and assuming a distance of $0.5\rm\>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 $\lesssim 2,000\rm\>K$, which is lower than the ionized gas temperature of $\sim 10,000\rm\>K$, even though it is almost completely optically thick. Such a model may be considered as an alternative explanation for the low brightness temperature, $\sim 1,500\rm\>K$, of Source I in Orion KL, rather than being due to $\rm H^{-}$ free-free emission from a collisionally ionized disk (Reid et al., 2007; Hirota et al., 2015). However, we note that, while the spectral index of Source I is $1.86\pm 0.26$ which is consistent with our model, its radio luminosity of $0.2\rm\>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\rm\>GHz$ (Plambeck & Wright, 2016), we need to expand the range of protostellar source model parameters, which we defer to a future paper.

4.2. Short-time variabilities induced by accretion variation

Apart from the break-out on $\sim 10^{3}\rm\>yr$ timescales, the photoionized structure can also change on shorter timescales due to accretion rate variability. Accretion bursts on $\sim 10$–$100\rm\>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\rm\>yr$. This is because the timescale of protostellar evolution is described by the averaged accretion rate over the local KH time, $t_{\rm KH}=Gm_{*}^{2}/(r_{*}L_{\rm*acc})$. The KH timescale of the protostar considered in this study is about $5,000\rm\>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.

We thank Viviana Rosero, Kazuhito Motogi, and Taishi Nakamoto for fruitful discussions and comments. J.C.T. acknowledges support from NSF grants AST 1212089 and AST 1411527. J.E.S acknowledges support from NASA grant NNX15AP95A.