Disk Wind Feedback from High-mass Protostars. IV. Shock-ionized Jets

Massive protostars launch accretion-powered, magnetically collimated outflows, which play crucial roles in the dynamics and diagnostics of the star formation process. Here we calculate the shock heating and resulting free–free radio emission in numerical models of outflows of massive star formation within the framework of the Turbulent Core Accretion model. We postprocess 3D magnetohydrodynamic simulation snapshots of a protostellar disk wind interacting with an infalling core envelope, and calculate shock temperatures, ionization fractions, and radio free–free emission. We find heating up to ∼107 K and near-complete ionization in shocks at the interface between the outflow cavity and infalling envelope. However, line-of-sight averaged ionization fractions peak around ∼10%, in agreement with values reported from observations of massive protostar G35.20-0.74N. By calculating radio-continuum fluxes and spectra, we compare our models with observed samples of massive protostars. We find our fiducial models produce radio luminosities similar to those seen from low- and intermediate-mass protostars that are thought to be powered by shock ionization. Comparing to more massive protostars, we find our model radio luminosities are ∼10–100 times less luminous. We discuss how this apparent discrepancy either reflects aspects of our modeling related to the treatment of cooling of the post-shock gas or a dominant contribution in the observed systems from photoionization. Finally, our models exhibit 10 yr radio flux variability of ∼5%, especially in the inner 1000 au region, comparable to observed levels in some hypercompact H ii regions.


INTRODUCTION
Massive stars, i.e., those with mass ≳ 8 M ⊙ , have a profound influence on the evolution of the universe via their radiative, mechanical and chemical feedback.However, the formation mechanism of massive stars remains under active debate, with two main theoretical paradigms being considered: Core Accretion and Competitive Accretion.Core Accretion is a scaled-up version of the standard model for low-mass star formation, in which self-gravity drives the formation of a concentrated core from which matter then accretes to a central protostar (Shu et al. 1987).The Turbulent Core Accretion (TCA) model (McKee & Tan 2003) extends this model by setting the boundary conditions of the initial pre-stellar core to pressure and virial equilibrium with a surrounding protocluster clump environment.Al-ternatively, Competitive Accretion (e.g., Bonnell et al. 2001;Wang et al. 2010;Grudić et al. 2022) proposes that massive stars form as part of a protocluster that globally funnels material to its central regions.Here protostars compete for the gas, with accretion proceeding in a chaotic manner without the presence of a massive, coherent core.Numerical simulations of competitive accretion generally predict that massive stars form relatively slowly compared to the TCA model, since the mass needs to be accumulated from larger scales.
While massive star formation models remain under debate, observations of their outflows can help constrain these models.One can simulate the outflow from a massive protostar under a given accretion scenario, and compare its predicted properties to observed systems.Understanding the launching and propagation of these outflows is crucial not only for testing formation models, but also to understanding how massive protostars impact their surrounding core and clump environments (e.g., Arce et al. 2007).
In this paper series we simulate disk-wind (Blandford & Payne 1982) feedback from a massive protostar forming from a 60 M ⊙ core in the context of the Turbulent Core Accretion (TCA) model of massive star formation (McKee & Tan 2003).In Paper I (Staff et al. 2019), MHD simulations at a given evolutionary stage were presented and basic properties, such as outflow cavity geometry and star formation efficiency were investigated.In Paper II (Staff et al. 2023), a full evolutionary sequence via a single simulation that followed growth of the protostar from low to high masses was presented.In Paper III (Xu et al. 2023, sub.) predictions for CO line emission from these simulations have been calculated.
Here, in Paper IV we present a model for the expected ionization due to shocks in the simulation and predictions for the associated radio continuum emission.The structure of the paper is as follows.In §2 we describe details of the simulations, calculations of shock parameters, and methods for predicting emission.In §3 we present the simulation shock-modeling results and compare them to observations.In §4 we discuss the caveats and limitations of the modeling and present the main conclusions of our study.

MHD Simulations
The post-processing analysis that we present in this paper uses snapshots from a 3D, ideal magnetohydrodynamic (MHD) simulation of a protostellar outflow interacting with its surrounding natal core infall envelope.The simulation domain includes one hemisphere of the protostellar core, from 100 au above the accretion disk mid-plane, up to a height of 25, 000 au (see Staff et al. 2023, for details).The outflow is injected into the simulation box at the lower z boundary.Mass can accrete from the envelope by flowing out of the lower z boundary, and this accreted mass results in the star and accompanying accretion disk growing (Staff et al. 2023).The rate at which mass flows through the z boundary is constrained so that the growth rate of the star matches the expectations of the TCA model, i.e., specifically the results of Zhang et al. (2014).Outflow boundary conditions allowing mass to flow out are used at all the other boundaries.
As the star grows, the mass and momentum input flux rates of the outflow also change with time, again following Zhang et al. (2014).Initially, the envelope has a mass of 60 M ⊙ and a radius of ∼ 12, 000 au.The envelope is initialized with a power-law density profile of the form ρ ∝ r −3/2 (McKee & Tan 2003), and is unstable to gravitational collapse.The simulation starts with a 1M ⊙ protostar at the center of the envelope.The envelope begins to collapse under the force of gravity, and its accretion and the propagation of the outflow are simulated for about 100,000 years, by which point the star has grown to just over 25 M ⊙ .
The simulation was run using ZEUS-MP (Norman 2000), with an isothermal equation of state and a fixed sound speed of 0.9 km s −1 .The simulation used a logarithmically stretched grid in all three dimensions, consisting of 168 × 280 × 280 cells.The cells are smallest near the outflow axis and near the lower z boundary, where the minimum cell size is ∼12 au.The initial core is threaded by a "Blandford-Payne" like poloidal magnetic field (Blandford & Payne 1982) plus a constant vertical component added to it to ensure a core flux of ∼ 1 mG.

Calculating Shock Temperatures, Ionization
Fractions and Free-Free Radio Emission A general overview of our method is the following.In order to model shock-ionization and the resulting emissions, we first calculate the shock velocities between cells and then use these to calculate post-shock temperatures (which is required because the MHD simulation is isothermal).Then, these temperatures are used to calculate the ionization fractions expected in the postshock gas.Two cases for the volume filling factor of the post-shock gas have been considered.Case A assumes the post-shock gas fills the entire cell into which the gas is flowing (averaging over all incoming flows, i.e., up to six).Case B, which is our fiducial model, estimates a cooling time in the post-shock gas and then assumes the filled volume is only a fraction of the cell equal to the ratio of cooling time to flow crossing time.Then, given the temperatures and ionization fractions, we evaluate the free-free emissivity in the radio regime and calculate the resulting radio spectrum from the computational domain.The details of these methods are explained in Appendix A.

Shock Structures
Snapshots from the simulation at various protostellar masses, ranging from 1.4 to 24 M ⊙ (corresponding to times from 4,000 to 94,000 yr) were selected for postprocessing to calculate shock heated temperatures, ionization fractions and radio free-free emission.Examples of the shock structures in a slice through the center of these snapshots (y = 0) are displayed in Fig. 1 for  4, 2, 4, 8, 12, 16, and 24 M⊙ at 4,000, 9,000, 21,000, 39,000, 54,000, 68,000, and 94,000 years, respectively.both an inner 4,000 au scale and the global 25,000 au scale.The density and z-velocity magnitude produced from the simulation are displayed in the first and second columns, respectively.The density map reveals a central low-density cavity, i.e., with n H ∼ 10 0 − 10 2 cm −3 , while the surrounding regions in the infall envelope reach densities of n H ∼ 10 7 − 10 9 cm −3 .This cavity increases in opening angle as it evolves, especially in the later stages from 8 to 24 M ⊙ .The z-velocities exceed 1000 km/s in the outflow cavity, but become much lower in the surrounding regions where the outflow is interacting with the infall envelope.
The third column shows the mass flux weighted postshock temperatures.Temperatures are highest, up to ∼ 10 7 K, on the boundary between the low-density, highvelocity outflow and the surrounding higher-density, low-velocity envelope.The fourth column shows these high-temperature regions are nearly fully ionized.Fi-nally in the fifth column the emissivities, j ν,ff at 5.3 GHz are shown.

Ionization Fractions
Ionization fraction in the outflow is one metric by which the results of our model can be compared to observations.The outflowing material was defined as including any cell with v z > v min , for v min = 1, 10, and 100 km s −1 .We averaged the ionization fractions over the line-of-sight (y) direction considering three different ways of doing the averaging: (1) mass-weighted, (2) volume-weighted; (3) emission-weighted.The resulting ionization fraction maps are shown in Fig. 2 and 3 for the 4,000 au and 25,000 au scales, respectively.The overall ionization fractions averaged by column density, area, and intensity weighting over the mass, volume, and emission-weighted maps, respectively, are summarized in Table 1.See Figs.B1 and B2 and Table B1 for the equivalent Case A results.
Higher velocity cutoffs lead to higher average ionization fractions because they exclude slower-moving, less ionized gas, particularly on the outskirts of the outflow.Weighting by radio free-free emissivity provides the highest and most high-contrast ionization fractions because, with ionization fraction and emissivity both dependent on shock temperatures and emissivity depending on χ 2 H+ , high emissivity correlates to high ionization fraction and vice versa.
For our fiducial case, Case B, which considers cooling, only the ionization from the portion of each cell flooded by the shock before cooling takes place is accounted for, leading to a reduction in the estimated ionization fractions, except for the emission-weighted metric.
Figure 4 shows the Case B profiles of ionization fraction along the outflow axis, averaging over strips of width ∆z = 1000 au (see Fig. B3 for Case A).Here we more clearly see that emission-weighting generally results in the highest ionization fractions, where regions above ∼4000 au are almost entirely ionized.Dips are present near 2000-4000 au, however, corresponding to a lack of shock-emission (see the emissivity column of Fig. 1).Higher velocity cutoffs correspond to higher ionization fractions, most noticeably in the mass-weighted case, attributable to the fact that the regions of slowermoving gas (1-10 km/s) on the outskirts of the outflow are denser than the fast-moving jet cavity.
Figure 4 also shows observed ionization fractions in four observed knots in the outflowing jet of the massive protostar G35.2-0.74N(Fedriani et al. 2019).The reported values are ∼ 0.1, with a slight decline as one goes from ∼ 3, 000 to 20,000 au.We note that these observations are based on measurements from regions that are emitting in the NIR via [FeII] emission lines, with this emission being used to derive the total number density of H nuclei.The electron density is derived independently from both the [FeII] emission and cm freefree emission from the knots, with each method giving similar results.The [FeII] line emission shows line of sight velocity shifts of ∼ 100 km s −1 .Thus we consider that the most relevant comparison between our results and the G35.2-0.74Ndata is for mass-averaged results for the case with outflow velocity ≥ 100 km s −1 .For this case, we see that our model ionization fractions can begin to match the observed data near 20,000 au when m * ≳ 12 M ⊙ .To match the inner knots requires somewhat higher protostellar masses.However, these results are sensitive to our choice of cooling model, i.e., Case B compared to Case A (see Fig. B3).The protostellar mass of G35.2-0.74N is quite uncertain, but from spectral en-ergy distribution (SED) fitting, Fedriani et al. (2023) estimate m * = 19 +9 −6 M ⊙ .Overall, we conclude that our simulation results for ionization fraction can match the observed values of G35.2-0.74N as long as m * ≳ 12 M ⊙ , although the profile along the jet and its detailed structure, i.e., presence of knots, show differences and/or are quite sensitive to modeling choices for the treatment of cooling and choice of velocity threshold to define the outflow material.

Intensity Maps
The free-free intensities at frequencies of 0.01, 0.05, 0.1, 0.5, 1, 5.3, 23, 43, 100, and 230 GHz were calculated for both Case A (no cooling) and Case B (with cooling).Examples are shown for the m * = 12 M ⊙ (54,000 yr) snapshot in Fig. 5 (see Fig. B4 for the equivalent Case A results).Although the morphology appears similar at each frequency, it is dimmest at 0.01 GHz, brightens significantly up to 0.5 GHz, then maintains similar brightness up to 230 GHz.A quantitative analysis of the radio spectral energy distribution is provided below in §3.3.2.
We also examine the time evolution of these intensity maps as the protostar grows in mass.Fig. 6 illustrates the evolution of the 5.3 GHz and 230 GHz intensity maps, along with mass-weighted, v cutoff = 100 km/s ionization fraction projections and ionized mass-weighted temperature projections.Fig. B5 shows the equivalent Case A results.The side-by-side intensity and projection maps show that the regions of highest intensity and ionization correspond to average ionized-gas temperatures of ∼ 10 7 K.The maps also show a gap in both 5.3 GHz and 230 GHz emissions corresponding to the low-shock activity region of 1000 au ≲ z ≲ 5000 au for snapshots up to 39,000 years, identifiable in Fig. 1 as having temperatures ≲ 10 4 K and little to no emissivity.
Profiles of the 5.3 and 230 GHz radio emission as a function of height were calculated by integrating the flux over 1,000 au wide regions in the z-direction, as shown in Fig. 7 (see Fig. B6 for Case A).These profiles clearly show the low radio flux region that is present around several thousand au in the early evolutionary stages (m * ≲ 8 M ⊙ ).This stage is associated with a relatively collimated outflow cavity, so that there is limited emission from shocks of the outflowing gas with the infall envelope at these locations.However, at later stages, when the opening angle becomes wider, then there are stronger shocks present at these boundaries.The averages are calculated over all cells in the y-column that exceed cutoff velocities of vz ≥ 1 km/s (columns 1, 4, 7 ), 10 km/s (columns 2, 5, 8 ), and 100 km/s (columns 3, 6, 9 ), plotted over −2000 au < x < 2000 au, 0 au < z < 4000 au.From top to bottom, the simulation grows in mass, reaching 1.4, 2, 4, 8, 12, 16, and 24 M⊙ at 4,000, 9,000, 21,000, 39,000, 54,000, 68,000, and 94,000 yr, respectively.Contours are overlaid displaying where the maximum z-velocity along the projection meets the cutoffs of 1 km/s (white), 10 km/s (pink ), and 100 km/s (purple).

Integrated Fluxes and Spectra
Integrated flux values at 5.3 GHz and 230 GHz were calculated for each evolutionary stage, i.e., m * = 1.4,2, 4, 8, 12, 16, 24 M ⊙ including flux within scales of r = 500, 1000, 2000, 4000, 8000, 16000, and 25000 au (doubled to account for the bipolar nature of the outflow).For this analysis, at each stage an average was made of 11 consecutive snapshots, each 10 years apart, covering the 100 yr period after each mass was achieved in the simulation.To assess the short timescale variability, we also measured the standard deviation of the values of log S ν (with these results discussed below in §3.3.3).The values of the fluxes and their standard deviations are presented in Table 2 for Case B (see Table B2 for Case A).
The evolution of the integrated fluxes as a function of increasing mass and bolometric luminosity is shown in Fig. 8.We see that the fluxes generally increase as m * grows, although there is a local minimum at m * ∼ 4M ⊙ , which is related to a stage of protostellar evolution when the star is relatively swollen so that its outflows are relatively slow (Zhang et al. 2014).At higher protostellar masses, the radio fluxes are seen to reach near constant values, although one should note that this is a regime where the photoionization process, which is not included in our modeling, will begin to play a more important role (Tanaka et al. 2016).bol with γ = 0.349, 0.563, 0.256, respectively.The semi-analytic photoionization model of Tanaka et al. (2016) for the fiducial TCA model, i.e., the same used in our numerical simulations, is shown by the red solid line.We see that in the luminosity range L bol ≲ 2 × 10 4 L ⊙ the contribution from photoionization is generally expected to be insignificant compared to shock ionization, with a minor exception near L bol ∼ 10 3 L ⊙ .This limited con-    Figure 9 also shows observational data from several samples of protostars.First, a sample of lowmass, lower-luminosity (≲ 10 3 L ⊙ ) protostars have been studied by Anglada (1995) and Anglada et al. (2015), including a power law fit (S ν d 2 /[mJy kpc 2 ]) = 8 × 10 3 (L bol /L ⊙ ) 0.6 .These sources are expected to be powered by shock ionization.We note that an extrapolation of this power law to higher bolometric luminosities predicts radio luminosities that are factors of several to ten times larger than our models of the entire 25,000 au Note-5.3GHz and 230 GHz fluxes are integrated over −r/2 < x < r/2, 0 < z < r for each scale r = 500, 1000, 2000, 4000, 8000, 16000, and 25000 au, multiplied by 2 to account for the counter jet.The exception is that for 25000 au, the integration is carried out over the whole simulation domain, which extends over −15000 au < x < 15000 au, 0 au < z < 25000 au.
domain.However, it is uncertain whether this extrapolation of this purely empirical relation is valid.The next data set is a sample of eight relatively massive protostars from the SOFIA Massive (SOMA) Star Formation Survey (De Buizer et al. 2017).These have had their 5.3 GHz fluxes measured on various scales (Rosero et al. 2019).The values shown in Fig. 9 are evaluated from the "Intermediate" scale, which is designed to capture any resolved radio jet (when there is no clear evidence for a resolved jet, the "Inner" scale is used).We see most of these protostellar sources have radio luminosities about 10 times greater than our Case B 25,000 au scale model.IRAS 20126+4104 sits close to this model, while G35.20-0.74Nand G45.47+0.05 are about 10 3 to 10 4 times more luminous than the shockionization model.A sample of ultracompact (UC) HII regions (Kurtz et al. 1994) have similarly elevated radio luminosities.
For G45.47+0.05 and many of the UC HII regions, which have high bolometric luminosities, ≳ 10 5 L ⊙ , photoionization is expected to be dominant.For the lower luminosity sources, it is generally unclear whether their radio emission is dominated by shock ionization or photoionization.However, in the case of G35.20-0.74N,Fedriani et al. (2019) argued that the modest ionization fractions of ∼ 0.1 indicate a preference for the shock ionization mechanism, at least in the jet knots.In this case, we see that our fiducial Case B models underpredict the required radio flux by a factor of ∼ 10 3 .We show in Appendix B that the Case A model, i.e., without a cooling limit, boosts the radio luminosities from the inner 1,000 and 4,000 au regions to ∼ 1 to 10 mJy kpc 2 levels, while the 25,000 au scale can reach the 100 mJy kpc 2 level that is relevant to G35.20-0.74.
The full spectrum of a protostar's radio emission provides additional information with which to constrain the models.The overall radio free-free spectra of the shockionized protostellar outflows are shown in Fig. 10, for both the inner (≤ 1, 000 au) and global (≤ 25, 000 au) scales.The spectra of the inner region show the expected optically thick behavior at low frequencies, i.e., S ν ∝ ν 2 , and then turnover at frequencies ∼ 0.1 GHz to the optically thin regime where S ν ∝ ν −0.1 .The spectra of the global region stay closer to the optically thin condition down to much lower frequencies, i.e., with flat or only slightly positive spectral indices down to 0.01 GHz.
We note that the radio spectral indices of G35.20-0.74N at the inner and intermediate scales are α inner = 0.7 ± 0.1 and α inter = −0.2± 0.1, respectively (Rosero et al. 2019).Such values are consistent with the general trends predicted by the shock ionization models shown in Fig. 10, where the inner scale shows a higher degree of optical depth.The radio flux variability was evaluated for a sample of 11 snapshots of the simulation separated by 10year periods, i.e., over a duration of 100 yr, at each fiducial evolutionary stage of the simulation, i.e., m * = 1.4,2, 4, 8, 12, 16, 24M ⊙ .These flux variations are shown in Fig. 11 (see Fig. B10 for Case A).

Flux Variability
At both 5.3 and 230 GHz, the average 10-year variability is σ(logS ν /[mJy]) ∼ 0.001 to 0.02, i.e., variations in flux of ∼ 0.2 to 5%.As expected, variations tend to be larger on smaller scales.On detailed examination of the radio intensity images, we find most variability is due to motion of particular hot spots of emission in the outflowing gas.Note, that a characteristic flow crossing timescale is t outflow ∼ 1000 au/1000 km s −1 ∼ 5 yr.
Figure 11 also shows the variability level of nine example observed hypercompact HII regions in W49A, which have measured 8.3 GHz flux variations over ∼ 20 yr timescales, i.e., from 1994 to 2015 (De Pree et al. 2018).Most of the sources show variations at the ∼ few percent level, which is quite comparable to the variability seen in our simulated outflows.However, the greatest change in flux identified by De Pree et al. (2018) was that of source W49A/G2, which decreased by 40% from 1994 to 2015.
We note that our model outflow has a smoothly varying outflow injection rate, while real sources may suffer some short timescale variation in accretion/outflow activity.We consider that such accretion/outflow variability is likely to help generate greater levels of radio flux variability, which is a subject worthy of future study.

DISCUSSION AND CONCLUSIONS
Shock ionization is expected to be the dominant source of ionized gas and associated radio emission in earlystage massive protostars.We have thus developed a model for shock-ionization of the disk wind driven outflow during an evolutionary sequence of massive star formation.In particular, the model has been applied as a post-processing step to the MHD simulations of Staff et al. (2023) of a protostar forming from an initial 60 M ⊙ core embedded in a Σ cl = 1 g cm −2 environment, i.e., the fiducial TCA model of massive star formation.The method involves estimating the post-shock conditions caused by converging flows between discrete cells in the simulation and then assessing the fraction of the cell that is filled based on the ratio of cooling to flow crossing time.Collisional ionization occurs in the shock heated gas.The free-free emission from this plasma was then calculated to predict fluxes and spectra in the radio regime.
The presented model is relatively simple and there are several caveats and limitations of the work.The first caveat is that the model involves post-processing, rather than a fully self-consistent MHD simulation with the full heating, cooling, and ionization processes.However, such numerical simulations would require very high resolution to be able to resolve the shock structures, highlighting the advantage and need of our subgrid model for the shocks.Nevertheless, there is still the question of the effect of numerical resolution on our results.To investigate this we have carried out "low" resolution simulations that have cells twice the size of those in our fiducial runs.We present their results for the radio spectra from the models in Appendix C. The overall result of increasing the resolution by a factor of two is typically to change the radio flux by within a factor of 2, although moderately larger differences can arise on certain scales and at certain frequencies (see Appendix C).
Another caveat is that the model is sensitive to how the cooling is treated for the post-shock gas.The simplest model with no cooling leads to much greater amounts of ionized gas and thus brighter radio emission, especially at higher frequencies where the structures are optically thin.The simplicity of our coolinglimited model, i.e., that sets the fraction of cells filled  with post-shock gas to be ∼ t cool /t flow , is such that one naturally expects systematic uncertainties in the quantitative model results at the level of a factor of a few.
As already mentioned, the model does not include any contribution to ionization from photoionization.Thus the model predictions should be viewed as lower limits, which we expect to be superseded when protostellar masses become large.Evaluating the contribution from photoionization as an additional postprocessing step is deferred to a future paper in this series.
Finally, when comparing to observed sources one must remember that the model has been applied only to a single evolutionary track in TCA model parameter space of inital core mass and clump environment mass surface density.Other parts of this parameter space involve protostars accreting at different rates and with different stellar radii, leading to variations in the outflow speeds, densities and thus shock conditions.Exploring a wider range of this parameter space is also deferred to a future study.
The main conclusions of our study are the following: 1. Shocks in a disk wind driven outflow, especially at its interface with the ambient infall envelope, produce localized regions that are heated to high temperatures, up to ∼ 10 5 to 10 7 K, with the maximum temperature reached increasing with increasing protostellar mass.
2. The resulting ionization fractions of the outflowing gas quickly reach levels of ∼ 0.1 when considering volume-weighted averages.However, massweighted averages only reach such levels when m * ≳ 16 M ⊙ and when selecting outflowing gas via a high velocity threshold of ∼ 100km s −1 .Yet, the ionization fraction of the material that contributes to the radio free-free emission (i.e., emission weighted) is close to fully ionized.The massweighted ionization fraction level of ∼ 0.1 is similar to that inferred in jet knots of the massive protostar G35.20-0.74Nby Fedriani et al. (2019).
3. In our cooling-limited model, the radio free-free emission from ∼ 1 to 300 GHz from the shockionized gas is in the range of ∼ 0.01 to 1 mJy, i.e., generally increasing with evolutionary stage from 1.4 to 24 M ⊙ , and increasing with the scale of the region considered from 1,000 to 25,000 au away from the protostar.At the inner scale of ∼ 1, 000au, these models become optically thick at frequencies below ∼ 0.1 GHz, while on the largest  scale, they remain mostly optically thin, i.e., with flat radio spectral indices.Models without the cooling limit have much more ionized gas and radio luminosities that are greater by factors of ∼ 10 to 1000, depending on the scale.
4. In all cases the emission from shock ionized gas is expected to dominate over that from photoionization for protostellar luminosities ≲ 10 4 L ⊙ .However, the precise transition from shock to photoionization dominance depends on details of the shock ionization model, i.e., treatment of cooling, and the protostellar formation conditions, i.e., accretion history as influenced by clump mass surface density.
5. A comparison of our models with observed protostars shows that our fiducial cooling limited case has radio luminosities similar to that seen in low and intermediate mass sources by Anglada (1995), but about 10 to 100 times lower than the massive protostars studied in the SOMA Survey (Rosero et al. 2019).Possible explanations for this include a more limited role for cooling in limiting the extent of post-shock ionization and/or a significant contribution from photoionization.More detailed comparisons, e.g., of full radio spectra and spectral index maps of jet structures, are needed to distinguish among these possibilities.
6. Our models produce 10-year variability in 5.3 GHz emission at the level of up to ∼ 5% on scales ∼ 1, 000 au to 25,000 au.Such short timescale variability is associated with changing gas conditions in the fast moving outflow.Similar levels of variability have been seen to be typical in a sample of HC HII regions in W49A by De Pree et al. (2018), which is strong evidence in favor of their ionized gas being associated with protostellar outflows.The higher level of variability, up to ∼ 40%, in a couple of sources may indicate a role for short timescale accretion/outflow variability, which has not yet been included in our modeling.
The authors thank Rubén Fedriani for discussions regarding his ionization fraction observations and Viviana Rosero for discussions regarding the SOMA-Radio observations.This research was initiated in the Virginia Initiative on Cosmic Origins (VICO) summer undergraduate student research program at the University of Virginia (UVA), Charlottesville, VA, USA, including a collaborative visit to the Chalmers Initiative on Cosmic Origins (CICO), Chalmers University of Technology, Gothenburg, Sweden.This work was also supported by the National Science Foundation (NSF) collaborative grant "Protostellar Jets Across the Mass Spectrum" (AST-1910675).The authors acknowledge the use of the Rivanna computing cluster, under the management of Research Computing at UVA. JCT acknowledges support from NSF grants AST-1910675 and AST-2206450 and ERC Advanced Grant 788829 (MSTAR).JES acknowledges support from NSF grant AST-1910675.JPR acknowledges support from VICO, the NSF under grant nos.AST-1910106 and AST-1910675, and NASA via  We refer to coordinates (z, x, y) indexed by (i, j, k) where z is the outflow direction.The methods are general for emission in any principle axial direction, but we focus specifically on adopting the y-direction as our line of sight direction and map x-z plane projections.
Shock velocities are calculated by taking, at each face of a cell, the velocity difference between the cell of interest and its neighbors, modulo only using converging velocity components.This yields the inward velocity to the current cell, in the reference frame of the gas in the cell.If this value is net positive in the inward direction (i.e., converging), it is set as the shock velocity at that face, otherwise, we take a shock velocity of zero.While these velocity differences will only be true shock velocities if they are larger than the sound speed, applying the shock calculation also to values below the sound speed has a negligible effect on the shock temperature.This process was replicated for every cell of the MHD simulation snapshot.
Given the shock velocities, a post-shock temperature is calculated for the shock at each face, according to the standard Rankine-Hugoniot shock jump conditions for the case that the Mach number, M, is ≫ 1: where γ = 5 3 is the ratio of specific heats, µ is the mass per particle and v s is the shock velocity.
The ionization fraction at each face from collisional ionization is calculated via (e.g., Draine 2011): where B = (157, 800K)k is the ionization energy for Hydrogen, n(H) is the number density of hydrogen, n(H + ) is the number density of hydrogen ions, ⟨σv⟩ rr and ⟨σv⟩ ci are the radiative recombination and collisional ionization rates, respectively.A temperature floor of 300 K for ionization is set, such that for cell temperatures below this value the ionization is taken to be zero since it can conservatively be assumed that any cell with a temperature below this value has zero ionization.With temperatures and ionization fractions calculated for the converging flow of each cell face, a single average value is determined for each cell by weighting each face's contribution by the associated mass flux.

A.1. Cooling Effects
The shocks in our simulation are modeled adiabatically, i.e., using an adiabatic equation of state to relate the pressure and temperature.For our initial approximation (Case A), the resulting temperatures and ionization fractions are assumed to fill the entire cell.This approximation is reasonable as long as the shocked gas fills the cell before significant cooling takes place.
However, for Case B, which is our fiducial case, we consider the effects of cooling.The cooling times can be calculated via (e.g., eq.35.34 of Draine 2011): where the radiative cooling function at low temperatures (T < 10 5 K) is approximated as Λ(T ) ≈ (3.98×10 −30 erg cm 3 s −1 )(T /K) 1.6 n H n e , (A5) and at high temperatures (T > 10 5 K) as: (A6) Accuracy of cooling times is more important for the high-temperature regions, as these are where most freefree emission occurs.
The post-shock cooling times are then compared to the flow times, i.e., the time it takes the shocked gas to fill the cell, given as a mass flux weighted average of the converging shock's timescale from each face: where v s /4 is the post-shock velocity and 1 2 ∆s is the distance from the face to the center of the cell.
Maps of the cooling time t cool , the flow time t flow , and their ratio t cool /t flow are shown in Fig. A1 for a slice from the snapshot corresponding to m * = 8 M ⊙ .It shows that there are significant regions with cooling times 10 times shorter than the corresponding flow times, thus motivating Case B in preference to Case A.
For Case B we modify the 1-D radiative transfer integration calculation in the following way.Every cell along the column of integration with a cooling time less than its flow time, has its integration depth ∆y scaled by that cell's ratio of t cool /t flow , as explained in Section A.2.In effect, this calculates emission only from the volume of the cell flooded by shocked gas before significant cooling has taken place.Results comparing Cases A and B are presented in Appendix B.

A.2. Calculating Free-Free Radio Emission
With the aforementioned mass flux weighted temperatures T and ionization fractions χ H+ , we then calculate the free-free radio emission due to shocks in the simulated outflow.The free-free emission coefficients j ν , absorption coefficients κ ν , and optical depths in each axial direction τ ν,z , τ ν,x , τ ν,y are (e.g., Draine 2011):  where g ff = 5.96(T /10 4 K) 0.15 (ν/GHz) −0.1 is the Gaunt factor and we assume n e = n p = n H+ = (χ H+ ρ)/µ H , where µ H = 1.4mH = 2.34 × 10 −24 g is the mass per H nucleus, i.e., assuming n He = 0.1n H . Note, our analysis ignores ionization of He, i.e., it assumes He is not ionized to a significant level.
For Case A, after determining the local emission and absorption coefficients for each cell, and assuming a uniform temperature within each cell, we calculate specific intensity, I ν , via integrating the 1D radiative transfer equation: In this work, we exclusively adopt the y-direction in the simulation snapshots as the line of sight.This analysis then yields maps of the emission from the simulation domain in the x-z plane.
For Case B, i.e., accounting for cooling effects, the radiative transfer equation is modified by scaling the integration depth through each cell by the ratio t cool /t flow if this ratio is less than 1.This serves to include emission only from the portion of each cell flooded by shocked gas before significant cooling takes place.Integrated fluxes were calculated for different-sized regions prescribed by −r/2 < x < r/2 and 0 < z < r, for r = 1000, 2000, 4000, 8000, and 16000 au, as well the full simulation region, which we refer to as 25000 au because it extends to 25000 au in the z direction.Although, the entire simulation region extends past 25000/2 au to ±15000 au in the x and y directions, by this point the x and y regions are large enough to cover the entire jet regardless, so it is the z range that is of most importance in characterizing how much of the jet is included.These fluxes were calculated as a surface integral of the intensity over the solid angle region Ω = A/d 2 , assuming a distance of d = 1 kpc, then doubled to account for the opposite hemisphere of the expected bipolar outflow that is not included in the MHD simulation.
In doing so for each snapshot and multiple frequencies, we obtain spectra and fluxes for a given frequency at a given time of the simulation.

A.3. Single Cell Examples
To demonstrate our methods on a simulation cell-bycell basis, we consider a partially ionized cell from the 39,000-year (8M ⊙ ) snapshot.This particular cell has indices (i, j, k) = (156,141,139) in the 3D snapshot, and is marked with a fuchsia × in Fig. A1.The relevant simulation data for cell (156,141,139) and each of its neighboring cells are given in columns 1 through 9 of the top section of Table A1.Only the neighboring cells' velocities that flow in the direction to/from the test cell are listed, as the other velocities have no impact on the test cell of interest.The velocity difference, shock velocity, shock temperature, and ionization fraction for each face are then given in columns 10-13 of Table A1, calculated according to Eqs. (A1)-(A3).Averaging the face temperatures T f and face ionization fractions χ H+,f by the mass-fluxes ρv f yields flux-weighted average values of T = 19, 300 K and χ H+ = 0.60 for the cell.
Since this temperature is less than 10 5 K, we use Eqs.(A5) and (A4) with n H = ρ/µ H to find a cooling time of t cool = 1025 yr.For the flow time, we average the timescales from all sides with non-zero shock velocities.In this example, those are cells i − 1, i + 1, and j − 1.This cell has ∆z = 7.24 × 10 10 km and ∆x = ∆y = 1.84 × 10 9 km.Eq.A7 gives time scales of 210 yr, 137 yr, and 614 yr, respectively.The flux-weighted average of these values yields t flow = 169 yr.Thus, we find a ratio t cool /t flow = 1025/169 ∼ 6.1, indicating that the cooling time is ∼ 6.1 times longer than the flow time, so it is reasonable in this case to assume the shock variables flood the cell before significant cooling occurs.
Cell (i, j, k) = (7, 122, 139) serves as a counterexample, where the flow timescale exceeds the cooling timescale.The cell is marked with a green × in Fig. A1 and the relevant data for each face are found in the middle section of Table A1.Averaging the T f and χ H+,f values by the mass flux yields T = 29, 300 K and χ H+ = 0.999 for the cell.Since this temperature is less than 10 5 K, we use Eqs.(A5) and (A4) to find t cool = 0.024 yr.To calculate the flow time, we consider the sides with non-zero shock velocity, i.e., j − 1, and j +1, with cell width ∆s = ∆x = 2.77×10 9 km.Eq. (A7 gives time scales of 4.88 yr and 4.87 yr, respectively, with a mass flux weighted average of t flow ∼ 4.9 yr.The ratio of t cool /t flow = 0.005 indicates that significant cooling will occur before the shocked gas floods the cell. Finally, we consider an especially hot, mostly (99%) ionized test cell with indices (i, j, k) = (165,191,139).The cell is marked in yellow in Fig. A1 and the relevant data for each face are found in the lower section of Table A1.Averaging the T f and χ H+,f values by the mass flux yields T = 4, 565, 000 K and χ H+ = 0.99 for the cell.Since this temperature is greater than 10 5 K, we use Eqs.(A4) and (A6) to find the cooling timescale, t cool = 4676 yr.The sides with non-zero shock velocity are now only i − 1, and i + 1, with cell height ∆s = ∆z = 7.24 × 10 10 km.This results in flow time scales of 15.39 yr and 12.35 yr, respectively, and a flux-weighted average of t flow = 12.75 yr.Thus, we find a ratio t cool /t flow = 4676/12.75∼ 367, indicating that the cooling time is ∼ 367 times the flow time, so it is again reasonable in this case to assume the shock floods the cell before significant cooling occurs.

B. RESULTS FOR CASE A (NO COOLING)
In Figures B1 to B10 we present results for Case A, i.e., without consideration of the cooling.

C. EFFECT OF VARYING RESOLUTION
The fiducial models in this paper use simulations with a numerical resolution of 168 × 280 × 280 cells logarithmically spaced over 25, 000 au × 30, 000 au × 30, 000 au.We conduct the same shock-modeling on a low resolution (84×140×140) version of this simulation, where all other parameters are kept the same, and present a comparison of the radio free-free emission spectra in Fig. C1.,191,139) Note-Each test cell of interest (i, j, k), is surrounded by neighboring cells i − 1 and i + 1 directly above and below it, respectively, in the z-direction, j − 1 and j + 1 left of and right of it in the x-direction, and k − 1 and k + 1 in front and behind it in the y-direction.For each neighboring cell, only data that contributes to the test cell of interest's flux-weighted temperature and ionization fraction is included, and all non-contributing entries are left blank.On the last line of each segment, the test cell's flux-weighted average temperature and ionization fractions are listed.
These results indicate varying effects of resolution.For Case B in both the inner (1000 au) and outer (25,000 au) scale regions, there are relatively modest variations in flux between the low and fiducial resolutions, i.e., most fluxes are the same within a factor of 2, although some larger variations can arise, especially at lower frequencies when differing levels of optical depth can amplify differences.For Case A, larger differences can appear, some of which are also amplified by the effects of optical depth.Some differences in spectral slope at high frequencies indicate the presence of individual cells with local conditions that lead to being optically thick even at 100 GHz.Such cells tend to be less common in the fiducial resolution simulations.

Figure 1 .
Figure 1.Time evolution of shock structures in slices through the center of the simulation in the x − z plane, over the ranges (left panel ) −2000 au < x < 2000 au, 0 au < z < 4000 au and (right panel ) −12500 au < x < 12500 au, 0 au < z < 25000 au.From left to right, the columns display number density, velocity in the z-direction, shock-heated temperature, ionization fraction, and free-free emissivity.Rows from top to bottom show the evolution of the simulation with protostellar masses of 1.4, 2, 4, 8, 12,  16, and 24 M⊙ at 4,000, 9,000, 21,000, 39,000, 54,000, 68,000, and 94,000 years, respectively.

Figure 9
Figure 9 presents the 5.3 GHz radio luminosity (measured by the metric S ν d 2 ) versus bolometric luminosity diagram.Our shock-ionization (Case B) models integrated on scales of 1000, 4000, and 25000 au are shown, along with power law fits S ν d 2 ∝ L γbol with γ = 0.349, 0.563, 0.256, respectively.The semi-analytic photoionization model ofTanaka et al. (2016) for the fiducial TCA model, i.e., the same used in our numerical simulations, is shown by the red solid line.We see that in the luminosity range L bol ≲ 2 × 10 4 L ⊙ the contribution from photoionization is generally expected to be insignificant compared to shock ionization, with a minor exception near L bol ∼ 10 3 L ⊙ .This limited con-

Figure 6 .
Figure 6.Time evolution of shock-emission variables projected along the y-direction for Case B (with cooling).Left set: Displayed domain is for inner 4,000 au scale.Right set: Displayed domain is for global 25,000 au scale.Within each set, from left to right, the columns show mass-weighted average ionization fraction of the jet with a velocity cutoff of 100 km s −1 , ionized-mass-weighted temperature, 5.3 GHz intensity, and 230 GHz intensity.tribution of photoionization is due to the typically low photospheric temperatures and relatively low luminosities of the TCA protostellar models for m * ≲ 10 M ⊙ , i.e., before Kelvin-Helmholtz contraction towards the zero age main sequence (ZAMS) (Fig.9also shows a reference model of photoionization from ZAMS stars).We also note that the predictions of the photoionization contribution at early stages are very uncertain since in this regime the accretion luminosity dominates L bol and the effective temperature of its emission depends on assumptions about the scale of the region where it is liberated: theTanaka et al. (2016) model assumes the protostar's internal and boundary layer accretion lumi-

Figure 7 .
Figure 7. Free-free emission fluxes for Case B (with cooling) at 5.3 GHz (top panel) and 230 GHz (bottom panel) as a function of height; fluxes are calculated by integrating the intensity Iν along the y-direction over 1000 au wide regions covering −25000 au < x < 25000 au, h − 500 au < z < h + 500 au for each height, h.Different colors represent the simulation at increasing protostellar masses: 1.4, 2, 4, 8, 12, 16, and 24 M⊙ are shown in pink, red, orange, yellow, green, blue, and purple, respectively.

Figure 9 .
Figure9.Radio continuum luminosity metric (Sν d 2 ) at 5.3 GHz versus bolometric luminosity.The solid, dotted and dashed lines correspond to the evolution in fluxes from 1.4 to 25 M⊙ on scales with r = 1000, 4000, 25000 au where the integration is carried out over −r/2 < x < r/2, 0 < z < r and then doubled to account for the counter jet (see Table2).Power law fits are shown by the thin lines of the corresponding line style.The red dashed line shows radio emission from an optically thin HII region photoionized by zero-age main sequence stars(Thompson 1984), while the red solid line is the equivalent model for the fiducial TCA model protostar (initial core mass Mc = 60 M⊙; mass surface density of clump environment Σ cl = 1 g cm −2 )(Tanaka et al. 2016).Yellow circles show low-mass protostar(Anglada 1995) with a power law fit of 8 × 10 3 (L bol ) 0.6 (black dashed line)(Anglada et al. 2015).Large circles are eight protostars from the SOFIA Massive Star Formation Survey(De Buizer et al. 2017;Rosero et al. 2019) (see legend).Crosses are ultracompact HII regions(Kurtz et al. 1994).

Figure 10 .
Figure10.Radio free-free emission spectra of the shockionized gas (Case B) in the 1,000 au region (top) and the full domain (bottom) for the different protostellar evolutionary stages (as labeled).Dotted line segments show the spectral indices of optically thick conditions (Sν ∝ ν 2 ) that are expected at low-frequencies and of optically thin condition (Sν ∝ ν −0.1 ) that are expected at high frequencies.

Figure 11 .
Figure 11.Flux variability, i.e., standard deviation of log Sν over 10-year intervals at 5.3 GHz (top panel ) and 230 GHz (bottom panel ) for 11 values spaced over 100 years, as a function of mass, for Case B. The black dotted lines in the top panel show the 20-year standard deviations of the 8.3 GHz flux from nine sources in massive star-forming region W49A; these regions are labeled from highest to lowest variability: G2, C1, G4, B, F, C, G5, G3, and D (De Pree et al. 2018).

Figure A1 .
Figure A1.Cooling time t cool (left), flow time t flow (center ), and their ratio t cool /t flow (right) for a central slice of the 39,000 yr, 8 M⊙ snapshot are presented for the region −2000 < x < 2000, 0 < z < 4000 in the top row and for −12500 < x < 12500, 0 < z < 25000 in the bottom row.Diverging color scales in the ratio maps show most of the region to have larger cooling times, with the exception of several strongly-emitting regions at the edges of the jet.The × symbols mark several test cells discussed in §A.3.

Figure B3 .
Figure B3.As Fig. 4, but now for Case A (no cooling).

Figure B4 .
Figure B4.As Fig. 5, but now for Case A (no cooling).

Figure B7 .
Figure B7.As Fig. 8, but now for Case A (no cooling).

Figure B8 .
Figure B8.As Fig.9, but now also showing the results of Case A (no cooling) (dark blue solid, dotted and dashed lines for 1000 au, 4000 au, and 25000 au scales, respectively).

Figure B9 .
Figure B9.As Fig. 10, but now for Case A (no cooling).

Figure B10 .
Figure B10.As Fig. 11, but now for Case A (no cooling).

Figure C1 .
Figure C1.Free-free emission spectra for low (dashed lines) and fiducial (solid lines) resolutions, for Case A (without cooling) in the top row and Case B (with cooling) in the bottom row.The left column represents flux from the inner 1000 au×1000 au and the right column represents flux from the entire domain.

Table 2 .
Radio Free-Free Fluxes at 5.3 GHz and 230 GHz for Case B, i.e., with Cooling the Astrophysics Theory Program under grant no.80NSSC20K0533.

Table A1 .
Test Cell Simulation Data