Modeling CN Zeeman Effect Observations of the Envelopes of a Low-Mass Protostellar Disk and a Massive Protostar

We use the POLARIS radiative transfer code to produce simulated circular polarization Zeeman emission maps of the CN $J = 1 - 0$ molecular line transition for two types of protostellar envelope magnetohydrodynamic simulations. Our first model is a low mass disk envelope system (box length $L = 200\text{ au}$), and our second model is the envelope of a massive protostar ($L = 10^4\text{ au}$) with a protostellar wind and a CN enhanced outflow shell. We compute the velocity-integrated Stokes $I$ and $V$, as well as the implied $V/I$ polarization percentage, for each detector pixel location in our simulated emission maps. Our results show that both types of protostellar environment are in principle accessible with current circular polarization instruments, with each containing swaths of envelope area that yield percentage polarizations that exceed the 1.8\% nominal sensitivity limit for circular polarization experiments with the Atacama Large Millimeter/submillimeter Array (ALMA). In both systems, high polarization ($\gtrsim$1.8\%) pixels tend to lie at an intermediate distance away from the central star and where the line-center opacity of the CN emission is moderately optically thin ($\tau_{LC} \sim 0.1-1$). Furthermore, our computed $V/I$ values scale roughly with the density weighted mean line-of-sight magnetic field strength, indicating that Zeeman observations can effectively diagnose the strength of envelope-scale magnetic fields. We also find that pixels with large $V/I$ are preferentially co-located where the absolute value of the velocity-integrated $V$ is also large, suggesting that locations with favorable percentage polarization are also favorable in terms of raw signal.


INTRODUCTION
Magnetic fields are pervasive across many scales in star formation environments and are expected to play an important role in regulating gas flows and structure formation at a variety of stages of the star formation process.On the molecular cloud scale (≳1 pc), magnetic fields both restrict gas flows and provide opposition to gravitational collapse (Mestel & Spitzer 1956;Mouschovias & Spitzer 1976).Cloud-scale fields are well studied with far-infrared/sub-mm polarimetry experiments such as Planck, BLASTPol, and SOFIA (Planck Collaboration Int.XIX 2015; Galitzki et al. 2014;Fissel et al. 2016;Chuss et al. 2019).These experiments leverage the well-established "radiative torques" (RAT) grain alignment theory, wherein spinning, effectively oblate dust grains preferentially align with their short axes parallel to the local magnetic field direction (Lazarian & Hoang 2007), yielding linear polarization that is perpendicular to the magnetic field (Davis & Greenstein 1951).
While dust grain alignment with the magnetic field is understood to be the dominant source of linearly polarized far-IR emission at larger scales, at smaller scales (∼100s-1000 au), however, the interpretation of linear polarization as a magnetic field tracer is more ★ E-mail: mazzei@virginia.edutenuous because other sources of polarization become important.In disks particularly, self-scattering of thermal dust emission (Kataoka et al. 2015;Yang et al. 2016), gas flow alignment (Kataoka et al. 2019), and "k-RAT" radiation field alignment (Kataoka et al. 2017;Tazaki et al. 2017) may all produce linear polarization signatures.Nonetheless, observing magnetic fields on these smaller scales remains of great interest.In older sources (Class II), magnetic fields are predicted to launch jets and winds along the disk axis (Blandford & Payne 1982), generate magneto-rotational instability (MRI; Balbus & Hawley 1991), and produce flows that contribute to the formation of rings and gaps (Suriano et al. 2017) or planetesimals (via zonal flows; Johansen et al. 2009).In younger (Class 0/I) protostars, magnetic fields may also hinder disk growth on 100+ au scales through magnetic braking, especially when the magnetic field aligns with the rotation axis of the surrounding envelope (Mellon & Li 2008;Hennebelle & Fromang 2008).
Another method for accessing magnetic fields in small-scale sources is to observe molecular line transitions in species that are sensitive to the Zeeman Effect (e.g., CN, OH, HI).In the presence of a magnetic field, the energy levels of these lines are split into higher and lower energy circularly polarized components, and the degree of splitting is proportional to the line-of-sight magnetic field strength (Crutcher et al. 1993).Historically, Zeeman observations have mainly been performed using single-dished telescopes to probe core-scale (∼0.1 pc) line-of-sight magnetic field strengths (see, e.g., Heiles & Troland 2004;Falgarone et al. 2008;Troland & Crutcher 2008;Crutcher et al. 2010).However, in recent years with the advent of a circular polarization observing mode on the Atacama Large Millimeter/submillimeter Array (ALMA), higher resolution (≲1 ′′ ) Zeeman observations are now possible.Disk-scale observations have proven challenging, with programs to-date yielding upper limit constraints only (Vlemmings et al. 2019;Harrison et al. 2021).Though disk-scale magnetic field strengths of up to ∼3 mG or more are predicted, the line-of-sight strength is likely suppressed in many disks due to cancellation within the torodial -field component (Mazzei et al. 2020;Lankhaar & Teague 2023).
Given the complicated situation inside disks, an alternative approach for searching for magnetic field signatures in protostellar environments is to probe the inner envelope regions beyond the edge of the disk.The advantages of the envelope are two-fold.First, the inner envelope-scale magnetic field is expected to be almost as strong as in a disk, while still remaining relatively unaffected by disk-scale tangling and toroidal wind-up.Linear dust polarization studies of these environments have revealed polarization percentages up to ∼5%, with relatively uniform geometry in some cases (Cox et al. 2018).Second, some envelope sources have been observed to have bright emission of Zeeman-sensitive molecules such as CN (see, e.g., Tychoniec et al. 2021).
In this work, we conduct 3D radiative transfer simulations to produce simulated emission maps of circularly polarized emission of the CN  = 1 − 0 transition for two different simulated protostellar envelope sources.Our first test case is a protostellar disk around a low-mass protostar, and the second the envelope around a massive protostar.For each simulation, we calculate the Stokes  and Stokes  obtained from many beams across the envelopes of our sources.We report the percentage polarization (/) and compare the implied line-of-sight magnetic field strength ( LOS ) from the Zeeman fitting with the actual mean (density-weighted)  LOS .There is substantial variability in the results for different beams due to non-uniform local magnetic field structure, but in some cases we find percentage polarization values ≳2%.
Given the nominal ALMA sensitivity limit of 1.8%, this suggests that Zeeman experiments with current instruments can be a useful way to study magnetic fields in (proto)stellar envelopes.However, program success depends critically on source selection and beam placement.It should be noted that our comparisons in this work are based strictly on radiative transfer simulations; we do not use the ALMA simulator to produce simulated observable products.Direct observational comparison of simulations to ALMA program data should involve use of the CASA simulator (McMullin et al. 2007).This is outside the scope of our work here.Nonetheless, from our results we are able to identify some general criteria for setting up observations to maximize detectability.Particularly, we find that regions some distance away from the central star at intermediate line-center optical depth (  ∼ 0.1 − 1) tend to be favorable.Furthermore, our results suggest that continued improvement of circular polarization instruments will be extremely fruitful; we predict that a factor of ten improvement in sensitivity (i.e., to a 0.18% limit) will make sightlines across nearly the entire inner envelope detectable for sources with typical magnetic field strengths (i.e., ≳ 1 mG, comparable to those found in our simulations).
This paper is organized as follows.In Section 2 we introduce the two model types we consider in this work, then in Section 3 we discuss our methods for making our simulated Zeeman emission maps.Section 4 presents the results of our simulations, including maps of the observables (Sec.4.1), computations of percentage polarization values through the envelopes of our sources (Sec.4.2), and comparisons of our pixel-by-pixel percentage polarizations to other associated observables as well as the underlying magnetic field information from our MHD simulations (Sec.4.3).We then provide some discussion in Section 5 on some additional criteria that can affect whether a given sightline or location is favorable for Zeeman experiments, including some auxiliary observational factors of which observers should also be aware (Sec. 5.4).Finally, the main conclusions of this work are summarized in Section 6.

NUMERICAL SIMULATIONS
Both sets of simulations studied in this work were performed using 3D grid-based MHD codes.

Low-Mass Protostellar Disk Envelope Simulation
Our first model type is an envelope and disk system around a lowmass protostar (henceforth known as lmde).The turbulent, non-ideal MHD model we used was the reference model of Tu et al. (2023).Here we highlight a few salient features of the model and refer the reader to Tu et al. (2023) for details.We use the ATHENA++ code (Stone et al. 2008(Stone et al. , 2020)), and include adaptive mesh refinement (AMR) and full multigrid (FMG) self-gravity (Tomida & Stone 2023).
The physical setup is similar to that used in Lam et al. (2019).We initialize a pseudo-Bonner-Ebert sphere prescribed by We set the central density  0 = 4.6 × 10 −17 g cm −3 , and characteristic radius   = 670 au.The sphere is embedded in a diffuse low-density gas ( amb = 4.56 × 10 −20 g cm −3 ).The total box size of the simulation domain is 10,000 au per side length, and the total gas mass is 0.56  ⊙ .The simulation is initialized with  −2 power spectrum turbulence (Kolmogorov 1941;Gong & Ostriker 2011) with an rms Mach number of unity, solid body rotation rate  = (  ,   ,   ) = (0, 0, 6.16 × 10 −13 s −1 ), and uniform magnetic field  = (  ,   ,   ) = (0, 0, 2.2 × 10 −4 G).The rotation rate corresponds to a ratio of rotational to gravitational energies of 0.03, which is in the range of the values inferred observationally by Goodman et al. (1993) and Caselli et al. (2002) through velocity gradients across dense cores.However, what fraction of the measured gradient is contributed by rotation remains uncertain.The adopted field strength corresponds to a dimensionless mass-to-flux ratio of 2.6, close to the median value for dense cores inferred by Troland & Crutcher (2008) through OH Zeeman measurements after geometric corrections.We note that although the magnetic field is initially uniform, it is significantly distorted by the imposed turbulence, gravitational collapse, and rotation at later times, particularly in the inner envelope surrounding the disk [see, e.g., Fig. 4a of Tu et al. (2023)].
We perform our polarization analysis on a simulation frame at  = 3×10 4 yr.By this time in the simulation, a stable disk has formed.
The disk remains stable for a while after this as well, but slowly decreases in mass at later times (Tu et al. 2023).Our chosen frame corresponds roughly to the time of maximum disk mass.Column density plots of the lmde model for each the −, −, and −linesof-sight are shown in Figure 1, along with 2-dimensional histograms of magnetic field strength versus gas density and distance from the (central) sink particle1 .
In Figure 2, we provide a 3-dimensional view of the magnetic field lines in the envelope.Whereas the magnetic field inside the disk is dominated by the toroidal component because of rotational wrapping [not shown here, but see discussion in Tu et al. (2023)], there are regions in the envelope where the field lines are more uniform and well-behaved.

Massive Star-forming Envelope (MSEnv) Simulation
The massive star formation simulation (henceforth referred to as model MSEnv) considered here was conducted using the Radiation-Magnetohydrodynamic (RMHD) code Orion2 (Li et al. 2021) with 5 levels of AMR (highest resolution Δ ≈ 18 au).The setup is similar to those described in Cunningham et al. (2011) but now with magnetic fields.The simulation was initialized as a massive cloud core of mass  core = 103 M ⊙ with a power-law density profile () ∝  −3/2 up to the radius of the core  core = 0.18 pc, which gives an average column density of 2.0 g cm −2 .We then add a turbulent velocity field to the core with rms Mach number M = 2.43, which corresponds to a virial parameter  vir = 1.67 for the cloud core.This value was chosen for the balance between gravity and turbulent motions during the protostellar system evolution; see McKee & Tan (2003) for more detailed discussions.The initial magnetic field strength of the core is chosen so that the dimensionless mass-to-flux ratio is 3.0 inside the core, slightly larger than those observed in high-mass star-forming clumps (∼ 1.5 − 2; see e.g., Crutcher 2012;Pillai et al. 2016;Motte et al. 2018) to ensure gravitational collapse.The initial gas temperature inside the core is  gas = 35 K, which is also set to be the temperature floor of the simulation box to avoid numerical rarefaction.A frequency-integrated flux limited diffusion (FLD) algorithm is adopted to approximate the radiation transport (see Cunningham et al. 2011 for more details).
In Orion2, the formation and evolution of protostars are handled by the star particle model described in Offner et al. (2009), which includes launching protostellar outflows.A tracer particle routine is implemented to record the properties of the ejected gas (Offner et al. 2009;Cunningham et al. 2011), and we utilized this function to trace the effective region of protostellar outflows in our simulated emission maps (see Sec. 3).
We investigate the frame when the central star is about 29.8 M ⊙ in mass, focusing on a  box = 10 4 au region around the massive protostar.For easier computation we reproject the AMR data to 256 3 unigrid array in python using yt (Turk et al. 2011).In Figure 3 we present column density plots and 2-dimensional histograms of magnetic field and density data for model MSEnv.
in a complex way on several factors, including the initial core mass and field strength, level of turbulence, and magnetic diffusivities.For example, Mignon-Risse et al. ( 2021) simulated magnetized disk formation in more massive cores with turbulence.They found milli-Gauss (mG) magnetic fields up to 1000 au, with the field strength reaching 100 mG or more on the inner 100 au scale (see their Fig.13).Their magnetic field is stronger than ours, making it more detectable in principle.We postpone a detailed exploration of parameter space to a future investigation.

SIMULATED ZEEMAN EMISSION MAPS
To characterize the observable magnetic field strength in each of the physical environments discussed above in Section 2, we produce simulated circular polarization emission maps of the CN  = 1 − 0 molecular line at 113 GHz.This line comprises a suite of nine hyperfine components 2 that are non-blended and stackable (Mazzei et al. 2020).For the majority of this work we choose to only simulate one representative sub-transition, the 113.144GHz component with relative intensity RI = 8 and Zeeman factor   = 2.18.In Section 5.4.2 we consider the potential signal boost that can be gained from stacking.
We perform Zeeman-splitting line emission calculations using the POLARIS raditaive transfer code with the ZRAD extension (Brauer et al. 2017), which incorporates data from the Leiden Atomic and Molecular DAtabase (LAMDA; Schöier et al. 2005) and the JPL spectral line catalog (Pickett et al. 1998).A Faddeeva function solver 3 is used to compute the final line shape, and included in the calculations are considerations for natural, collisional, and Doppler broading, as well as the magneto-optic effect (Larsson et al. 2014).
Our POLARIS line radiative transfer calculations are computed on 128 3 and 256 3 fixed resolution grids for models lmde and MSEnv.These translate to resolutions of 1.56 au and 39.1 au, respectively.For each cell, local gas density, gas velocity components, and magnetic field components are supplied from our MHD simulations.We also must specify a CN abundance value for each cell.For our lmde model we choose constant CN abundance  CN = 10 −9 .This value ensures a moderately optically thin envelope beyond the disk edge at  ∼30 au, and is similar to estimates for CN abundances in disk modeling (Cazzoletti et al. 2018).In the massive star (MSEnv) model, we also set  CN = 10 −9 , but only within a "shell" region that has cells containing tracer particles with values in a given range.This is meant to simulate a CN enhancement at the edge of a cavity blown out by a protostellar wind that is more directly exposed to the UV radiation from the central massive protostar system.Such regions are observed in some cases to be enhanced due to far-ultraviolet dissociation chemistry (Arulanantham et al. 2020).We set the range of tracer particle that define the shell by hand, choosing values that empirically produce a reasonable result.Slice and projection plots of our shell are shown in Figure 4.For the remainder of the massive star envelope (i.e., all cells not within the shell), we set a lower ambient CN abundance  CN,amb = 10 −12 .
Each line emission simulation converts the position-positionposition (PPP) data cube into a position-position-velocity (PPV) data cube.Values for the Stokes , Stokes , and optical depth  are recorded in each pixel for 169 velocity channels in velocity range [−3 km/s, 3 km/s] with respect to the rest frame of source.For all of our simulated emission, we assume local thermodynamic equilibrium (LTE) and set  = 10 K.A total number of 10 5 unpolarized background photons are initialized in each run, which are ray-traced to a 256 2 detector.For simulations of our lmde model, we place the detector at a distance  = 150 pc from the box center, and for simulations of our MSEnv model we use  = 1 kpc.These values are based on typical distances to nearby regions of low-mass star formation (e.g., the Taurus molecular cloud) and massive star formation, respectively.
For a locally uniform magnetic field that is weak enough such that any Zeeman splitting is unresolved (Crutcher et al. 1993), the following relationship applies: where  LOS is the line-of-sight component of the magnetic field.
Using this relationship, we may compute the line-of-sight magnetic field strength required (under uniform magnetic field conditions) to achieve the nominal ALMA percentage polarization limit of 1.8%.
For this, we ran a sample POLARIS simulation with a uniform magnetic field (oriented along the line-of-sight) in a box with uniform gas density.We then calculated / from the Stokes  profile and found its maximum value, (/) max , as well as the maximum value of the Stokes  profile,  max .Under uniform conditions, the required line-of-sight magnetic field for 1.8% is then which yields a value of  LOS ≈ 3.1 mG in our case.Throughout this work, we refer to this 3.1 mG estimate (in addition to the 1.8% observational limit) as guidance for assessing potentially detectable magnetic field conditions in our simulations.We also sometimes refer to a hypothetical factor-of-ten improved limit at 0.18% and 0.31 mG, respectively.This exercise is undertaken to evaluate the potentially improved utility of conducting Zeeman observations with a next generation circular polarization instrument of the future.
In Table 1, we report the fraction of cells in each of our grids with magnetic field component values above 3.1 mG and 0.31 mG.

RESULTS
Presented in this section are our circular polarization results for both the lmde and MSEnv models.We include integrated emission maps (Section 4.1) and computations of the local polarization percentage in locations across different cuts of the simulated observer-space (Section 4.2).We also use 2-dimensional histograms to relate our polarization data to local magnetic field information (Section 4.3).

Maps
In Figure 5 and Figure 6 we show maps (for each of the −, −, and −sightlines) of the observable data obtained from our POLARIS simulations of model lmde and model MSEnv, respectively.Included Table 1.Percentage of 3D cells on our re-projected fixed grids that have magnetic field component values above 3.1 mG and 0.31 mG.In addition to the full domain, we also perform computations that only include cells within a given radius; the chosen radii are the same we use for cuts in our polarization analysis in Section 4.2.For model MSEnv, we also report values when only the cells within our CN enhanced shell are considered.are the velocity-integrated Stokes  and Stokes , the CN optical depth at line-center (  ), and a derived percentage polarization quantity that we refer to as the "Median /".To calculate this value, we bin down our data from 169 channels to 13 channels, resulting in a velocity resolution of 0.46 km/s that is comparable to that which is typical of Zeeman observations with ALMA.For each pixel in the observer space, we then calculate the ratio of the Stokes  to the Stokes  in each of the 13 re-binned channels.The "Median /" for a given pixel is then the 7th highest of these 13 values.Functionally, this parameter produces similar percentage polarization values to comparing the peaks of the Stokes  versus Stokes , while also imposing the requirement that the majority of channels must have detectable emission (which is a soft criterion for being able to reasonably estimate the local line-of-sight magnetic field strength, e.g., by applying equation 3).Generally, the maps for both models show that the percentage polarization is maximized in at an intermediate distance from the central protostar.In model MSEnv, for all three viewing angles there is a low polarization "hole" near the star (within  ≲ 500 au) where the line-center optical depth is high (  > 1).Moreover, the locations of high polarization are spatially correlated with the location of the enhanced CN outflow shell defined by our tracer particle constraint.
Our results from model lmde also reveal that there is polarization sub-structure when we focus in on the  ≲ 100 au scale near the disk.In this case, however, rather than a roughly spherically symmetrical low polarization "hole," the low polarization region is defined by the shape of the (optically thick) disk.For the edge-on views (i.e., the − and −sightlines), the percentage polarization is low along the disk midplane, with pixels that exceed 1.8% polarization fraction mainly occurring between 5 − 75 au above or below the midplane.In the  line-of-sight, there is also clearly a region of high polarization that spatially overlaps with the location of a disk wind ( ≳ 10 au along the disk axis).Looking face-on at the disk ( line-of-sight), there is a circular low polarization region near disk center, again corresponding to where   > 1.From this view, pixels in excess of 1.8% tend to occur in a radial band  ∼ [30 au, 50 au] surrounding the edge of the optically thick disk.
The polarization percentage is also, unsurprisingly, low in the outskirts of the box domains, in the areas far beyond the CN enhanced shell or disk edge, for models MSEnv and lmde respectively, where the magnetic field tends to be weaker in general.

Percentage Polarization Statistics
As a corollary to Table 1, where we listed the fraction of 3D cells with magnetic field strengths above 3.1 mG and 0.31 mG, in Table 2 we calculate the fraction of pixels from our simulated emission maps that have percentage polarization ("Median /") values above 1.8% and 0.18%.For both simulation types we perform these computations for the full observer-space, as well as for cuts that restrict the domain to more favorable sub-regions (from a percentage polarization perspective).The locations of these cuts are also noted in Table 2.
For the low mass disk model lmde, we find that across the whole observer space the fraction of cells with fractional polarization above 1.8% is between ∼20 − 40% depending on viewing geometry.Interestingly, the  line-of-sight in this case has a higher fractional polarization value largely due to a swath of high polarization pixels along a disk wind.When Cut 1 is applied (i.e., the highest column density cells corresponding to the high  disk are removed from consideration), there is effectively no change in the fraction of cells with percentage polarization above 1.8%.This is sensible, since the disk only accounts for a small area of the observer space.Visually, as we see in Figure 5, it is clear however that the inner disk has quite low polarization.When we add in the restriction to limit the probed space to  < 75 au (Cut 2), the percentage polarization increases by factors of ∼2, 1.5, and 2 for the −, −, and −lines-of-sight, respectively.This calculation aligns well with what is suggested in the maps; that the regions of high polarization tend to lie spatially in an intermediate range just beyond the disk edge, where the magnetic field is relatively strong and the optical depth is not too high.
The results for the massive star case are qualitatively similar, in that the fraction of pixels above 1.8% polarization is highest when only an annulus is considered.For the whole domain, the percentage of high polarization (>1.8%) pixels is about ∼15% for all three Cartesian sightlines.This value increases by about a factor of 2 up to ∼30-35% when only pixels within 2500 au of the central star are included (Cut 1).We also define a cut (Cut 2) that additionally removes the polarization "hole" in the  line-of-sight view.The location of this mask is set manually to be centered on the lowest polarization pixel in that region, with a radius of 900 au.Removing the hole from consideration increases the fraction of pixels with percentage polarization above 1.8% from 35% to 41% for the  line-of-sight.
Finally, it is also worth noting that for both models and all cuts, the fraction of cells with percentage polarization above 0.18% always exceeds 90%.This suggests that a hypothetical factor of ten improvement in sensitivity to the fractional circular polarization would result in the vast majority of the envelopes becoming accessible with Zeeman observations.

2-Dimensional Histograms
Here we relate our percentage polarization / results to the underlying magnetic field information as well as other observable quantities of interest (i.e., Stokes  and CN   ).Two-dimensional histogram plots for models lmde and MSEnv are presented in Figure 7 and Figure 8, respectively.
Each plot has significant scatter, suggesting that high polarization can in principle be obtained under a wide variety of envelope conditions.However, there are also some trends apparent in our data that show some characteristics are more favorable than others.As expected, Median / scales positively with the density weighted average line-of-sight magnetic field strength (first column of Figures 7  and 8).Moreover, for both simulations the associated mean trend line tends to pass through a polarization percentage of 1.8% at roughly 3.1 mG.This suggests that this equivalent value for a uniform magnetic field we calculated in Section 3 is a useful guide for interpreting envelope Zeeman signatures.That is to say, it is reasonable to expect that robust detections with ALMA should require average magnetic fields strengths ≳ 3 mG.Note, however, that this value applies only in envelope regions like those we simulate here, where the magnetic Table 2. Computations of the fraction of pixels in each of our simulated emission maps with percentage polarization values above 1.8% and 0.18%.We report values for the entire box domain for each of the Cartesian lines-of-sight, as well as for some selected cuts within the observer-space.field is mostly uniform or only moderately tangled.In regions with more complex (i.e., tangled or wound-up) magnetic field geometry, larger average field strengths will be required to achieve the same polarization percentage.It is also the case that the maximum magnetic field strength (along the line-of-sight) scales positively with median / (second column of Figures 7 and 8).For pixels with percentage polarization greater than 1.8%, this maximum value always exceeds 3.1 mG and can reach as high as ∼10-30 mG, which is consistent with the expectation that the Zeeman-inferred field strength is a lower limit to the maximum field strength found along the line-of-sight.Another notable trend is that pixels with strong velocity-integrated | | signals tend to also have large /.This is especially true at large | |.This result is convenient for observers, because it means that locations with high percentage polarization will also be the most likely to have a detectable Stokes  signal.
For both model types, there is also a turnover in percentage polarization at an intermediate radius.The peaks occur at approximately ∼30 au for model lmde and ∼1500 au for model MSEnv.As we touched on in Section 4.1, visual inspection of the maps of the observables reveal that these radii correspond approximately to the edge of the optically thick disk and the polarization "hole" for the two model types, respectively.
Assessing percentage polarization trends with line-center optical depth is less clear.Overall, there is a wide range of   values that can produce potentially detectable percentage polarization levels in both envelopes.In model lmde there tends to be a moderate anti-correlation between optical depth and /, especially in the opacity regime where the majority of the pixels tend to be located (log [  ] ∼ [−2, 0]).For the face-on case, there is also a local maximum (in the mean curve) around   ∼ 0.1, which lends some credence to the idea the intermediate optical depth may sometimes be favorable.For model MSEnv there tends to be low polarization at very small   ∼ −4.This is due to low-density ambient region far outside the CN enhanced shell.Within and near the shell (where log [  ] ∼ [−2, 1]) the percentage polarization tends to be close to 1.8% with very little discernible trend, except for a downturn at the highest opacities (log [  ] ≳ 0) for the  and  lines-of-sight.This modest anti-correlation corresponds to the polarization "hole," which is evidently not as visible in the −view.
One further caveat that we should note about optical depth is that our simulations are tailored toward the optically thin scenario; we chose our CN abundances such that we would generally have   < 1 throughout the envelopes of our modeled sources.This is in part because we expect modeled optically thin emission to be more directly comparable to observations than optically thick emission.In practice, optically thick regions are subject to additional line effects that make drawing comparisons to modeled emission more tenuous, such as continuum over subtraction and absorption by resolved-out, cold outer envelope material.

DISCUSSION
Generally, the results of our simulations suggest that observable / signals (with e.g., ALMA) are available in both massive and low-mass protostellar envelopes.Therefore, probing magnetic fields in these environment with the Zeeman effect is a tractable goal.Our results also show, however, that there is significant variation in percentage polarization across the envelopes; that is to say, some sight-lines are more favorable than others.
In this discussion, we investigate more deeply some factors that can influence whether a given location is a good candidate to have high percentage polarization: optical depth (Section 5.1), CN enhancement on a wind-blown shell (Section 5.2), and disk inclination (Section 5.3).We then highlight in Section 5.4 some additional technical factors that observers should be mindful could also affect the percentage polarization results.Finally, for our lmde simulation we provide example lines (Stokes  and Stokes  profiles) from a selected envelope location to give a sense of typical line morphology and linewidth in this model (Section 5.5).We also compare the magnitude of the integrated flux found in this simulation to TMC-1, an example Class I source that is a good candidate for inner envelope Zeeman observations.

Optical Depth
From our 2-dimensional histogram analysis in Section 4.3 we concluded that global trends in / versus line-center optical depth are not especially clear for either model type.Here, we re-visit the topic with an alternative approach.In Figure 9 we plot the radial profiles of two quantities, the fraction of cells with polarization percentage above 1.8% and the average   , within thin annuli of radial size [ − Δ/2,  + Δ/2].For model lmde we set Δ = 6 au and for model MSEnv we set Δ = 300 au.
Our plots for the model lmde show different trends depending on line-of-sight.For the edge-on disk views (− and − lines of sight), the percentage of pixels above 1.8% polarization percentage is relatively low at small radius ( ≲ 15 au) then peaks around  = 20 − 50 au.In this radial range the average optical depth hovers around   ∼ 0.2 − 0.4.Visual inspection of the maps (top two rows of Fig. 5) shows us that these results are due to high polarization swaths above and below the disk midplane.The face-on view of model lmde, by contrast, shows that the annuli with the highest fraction of high polarization cells also have the highest optical depth.While it is certainly the case that some pixels in the optically thick disk have high polarization, taken at face value this result is perhaps somewhat misleading.Looking at Figure 5 again (bottom row), we can see that most of the pixels with median / > 1.8% are either in the transition region near the disk edge (where  is rapidly decreasing) or further away.Furthermore, in the face-on view of the disk there are opposing radial asymmetries in both the optical depth map and the / map in the range  ∼ 30 − 50 au.Particularly, there is a lobe of high   at positive , whereas the high / extension is at negative  where the CN optical depth is actually lower, at around   ∼ 0.1 − 1.
For model MSEnv, all three lines-of-sight show similar trends.Though at small  there is some difference in where the high polarization regions lie as a function of radius (this is due to slight differences in the location of the central polarization "hole" for each projection), within  ≲ 2500 au it is generally the case that peaks of the fraction of pixels with / > 1.8% are accompanied by dips in the average line center optical depth (to   ≲ 1).Outside  = 2500au, both quantities gradually taper off toward the outskirts of the envelope.

CN Enhancement
In our simulations of model MSEnv we included a non-uniform abundance prescription, with 3D cells lying within a prescribed shell region (see Fig. 4) being supplied with a CN abudance a factor of 1000 larger ( shell = 10 −9 ) than the surrounding ambient material ( shell = 10 −12 ).The results of our radiative transfer simulations for this model therefore give us some opportunity to comment on how For each pixel these quantities include, from left-to-right, the density weighted average line-of-sight magnetic field strength, the maximum magnetic field strength in any 3D cell along that sightline, the absolute value of the velocity-integrated Stokes  signal, the line-center optical depth, and the distance (in observer space) from the central sink particle.the locations of CN enhanced regions, particularly those which may be driven by UV irradiation on the cavity wall from a protostellar wind, can inform observational choices when designing a Zeeman experiment.
In our plot of the observable data from the massive star simulation (Figure 6), we annotate reference lines on the Stokes  maps (white dashes) and median / maps (black dashes).These lines are co-located and drawn in by hand strictly for visual reference.From inspection of the maps it is clear that the regions of high polarization are well-correlated with the CN enhanced shell.Particularly high polarization occurs near the shell edges surrounding the central polarization "'hole" or in its lower-intensity (as measured by Stokes  value) fringes.For example, the blown out "tip" of the shell located at (, , ) ∼ (−3000 au, −1000 au, 3000 au) contains many high polarization pixels.These results suggest that interfaces between outflow structures and ambient gas (e.g., due to outflow-swept shells) are prime areas to target, due both to favorable CN abundances and relatively low to intermediate optical depth compared to (low polarization) central regions.

Disk Inclination
So far in this work we have only considered simulations viewed from the Cartesian lines-of-sight (with respect to the frame of our MHD simulations).While this may be sufficient to obtain a good understanding of the MSEnv model, which is more-or-less spherically symmetric in its interior region, our results for the lmde model show significant contrast in polarization map morphology between the face-on and edge-on disk views.
Figure 10 presents percentage polarization results for a series of intermediate disk inclination views of our low mass protostar simulation.The general picture is similar for all viewing orientations; there is a central (highest intensity in Stokes  and highest optical depth) region that is low in polarization, and it is surrounded by a relatively high polarization envelope.The shapes of the regions with large / vary as a function of inclination.As the angle is adjusted from  = 90 • (edge-on) to  = 0 • (face-on), the high polarization swaths gradually progress from appearing shaped like blocks above and below the disk midplane into a ring near the disk edge.From this experiment we can see that intermediate inclination low mass protostellar envelopes are potential targets for Zeeman experiments, as they also contain many high polarization pixels.They just perhaps have more complex morphology, with the high polarization regions appearing as some combination of "block-like" and "ring-like."

Additional Observational Considerations
In Sections 5.4.1 and 5.4.2 below, we briefly assess the impact beam convolution and stacking of sub-transitions can have on our observational results by considering a few case studies.

Beam Convolution
We test the effect of applying 0.5 ′′ FWHM beam convolution to the  line-of-sight view of our MSEnv model.For this experiment we bin our 169 channels of Stokes  and  data to 13 channels (corresponding to an observed velocity resolution of 0.46 km/s) as usual, then at that stage apply a Gaussian filter to each of those 13 channels.We then use those data to compute velocity-integrated  and , line-center CN optical depth, and the median /.
A side-by-side comparison of our high resolution maps (i.e., with no beam convolution) and our  FWHM = 0.5 ′′ maps are presented in Figure 11.We also include histograms of each observable to compare their values in aggregate.In addition to the obvious morphological changes from Gaussian smoothing, the main effect is in the / map.Clearly, the values of the pixels with the highest percentage polarization are significantly reduced.This is an expected result, as beam convolution effectively averages adjacent pixels.Therefore, any high polarization regions that are smaller than the beam size will have some contribution from low polarization pixels post-convolution.This effect is captured in the histograms as well.The overall Stokes  and  distributions are similar in shape, but the convoluted / distribution is notably narrower than the unconvolved version.The percentage of pixels above 1.8% fractional polarization is ∼15% for the unconvolved map and ∼3% for the convolved map.Meanwhile, the percentage of pixels above 0.18% fractional polarization is ∼94% for the unconvolved map and ∼99% for the convolved map.It should also be noted that the two / distributions peak at roughly the same value, / ∼ 0.8%.
These results suggest that either high sensitivity (/ limit < 1.8%) or high resolution ( FWHM < 0.5 ′′ ), or a combination of both, would be of great advantage for Zeeman experiments in this type of environment.

Stacking Sub-Transitions
The CN  = 1 − 0 molecular line comprises nine velocity-resolved hyperfine sub-transitions.For the preceding sections of this work we only considered one representative transition, the 113.144GHz component.Here, we perform simulations of all seven4 of the subcomponents that are available by default in POLARIS.We run these computations on the face-on ( line-of-sight) view of our low mass disk envelope simulation, which serves as a good test case because it contains a wide variety of optical depth conditions.
We find that each of the seven components have nearly identical Stokes , , and  morphology, so the main effect of stacking is a boosted gain in signal.In Figure 12 we compare velocity-integrated  and , and median / maps obtained from the stacked transitions versus those from the representative transition.We also compute the ratio between the stacked value and representative transition value for each of these quantities.Stacking produces a (per pixel) boost in the  signal by a factor of ∼7-10, with the distribution strongly peaked at roughly factor of 8.The  signal is boosted by roughly the same factor on average, but the distribution is a bit broader; some pixels reach as high as a factor of ∼25 brighter.Finally, the median / is about the same for the stacked data as it is for the single subtransition.Interestingly, the regions with the lowest values for this quantity [(/) stacked /(/) Rep. ∼ 0.8] tend to be located within the optically thick disk, and the regions where stacking increases / tend to be in the envelope.

Example lines in model lmde and comparison of integrated Stokes 𝐼 flux to TMC-1
From the perspective of an observer, it is also useful to get a sense of what the Stokes  and  line profiles look like in a typical envelope beam.In Figure 13, we provide the velocity-integrated brightness temperature (in K km s −1 units) for each pixel from the face-on view of model lmde, as well as Stokes  and  profiles taken from an example location in the envelope.The line shapes observed in this location are generally representative of what we see across the envelope, with roughly Gaussian Stokes  emission (the particular shape of the line, however, will of course depend on the velocity structure of the chosen location, and the viewing angle).Notably, the Stokes  morphology matches well with the / fit, which is consistent with the expectation of a fairly uniform, well-behaved line-of-sight magnetic field in the envelope (see e.g., Fig. 2).Incidentally, the percentage polarization we obtain at this particular location is ∼1%.
Using the velocity-integrated surface brightness in each pixel (i.e., the left panel of Figure 13), we can also compare the fluxes obtained in our modeling with real observed fluxes.Particularly, TMC-1 is a Class I source with known bright CN emission (Tychoniec et al. 2021).It has been observed to have CN  = 2 − 1 ( 0 = 226 GHz) emission with velocity-integrated surface brightness values between ∼50-150 mJy beam −1 km s −1 in its inner envelope.Since the line transition with which TMC-1 was observed is different from that we used in our modeling (CN  = 1 − 0, with  0 = 113 GHz), we convert to velocity-integrated brightness temperature for this comparison.Noting that the observation of TMC-1 in Tychoniec et al. ( 2021) used a beam size of ≈0.16 arcsec 2 , we calculate the velocity-integrated specific intensity and then use the Rayleigh-Jeans equation to compute the corresponding velocity-integrated brightness temperature.The result we obtain is that 50-150 mJy beam −1 km s −1 corresponds to approximately 8-25 K km s −1 .These values are roughly comparable to those found in our simulated disk-envelope system, which has velocity-integrated brightness temperatures of ∼8-20 K km s −1 in the inner envelope.The modeling we've conducted throughout this work suggests that envelope locations with CN surface brightnesses in this general range are able to yield detectable emission with fractional polarization of order ∼2%.Sources with known bright CN emission (such as TMC-1) should be prioritized when considering targets for envelope Zeeman experiments.

SUMMARY
In this work, we produced simulated Zeeman emission maps of the CN  = 1 − 0 molecular line transition from MHD simulations of low-mass star formation and massive star formation regions.We calculated the / percentage polarization throughout the envelopes of each environment and then placed these results in context by comparing them with the local magnetic field data supplied from the 3-dimensional inputs.We also compared our results with nominal instrumental limits to assess the current and future feasibility of using Zeeman observations of protostellar envelopes to broadly diagnose magnetic field character during the early embedded phase of star formation.Our principal conclusions are summarized below: (i) In 3D, both models contain cells with magnetic field strengths that are in principle sufficiently large to produce circularly polarized emission that may be detectable with current instruments (Table 1).Intrinsically, roughly ∼45-60% of cells (depending on viewing orientation) in our low mass disk envelope system and ∼25-45% of cells in our massive star envelope system have local line-of-sight magnetic field strengths that exceed 3.1 mG, the line-of-sight magnetic field strength we estimate to in principle be needed to reach a polarization percentage of 1.8%, the nominal ALMA limit.Furthermore, if we ignore the low density diffuse outskirts (by placing cuts at  < 75 au and  < 2500 au for models lmde and MSEnv, respectively), these percentages increase to ∼60-70%.(ii) Our simulated emission maps yield pixel-by-pixel polarization re-sults that vary significantly across the observer space.Each simulation contains some regions with very low polarization and others that well exceed 1.8% (see, e.g., the final column of Figures 5 and  6).Broadly, the low polarization regions tend to be near the edges of the simulation box or in the central, highest intensity and optical depth sight lines.This leaves an intermediate range (in both radius We generally find that the amount of signal boost gained from stacking has some modest pixel-by-pixel variation, with both distributions centered on a factor of ∼8. Figure 13.Left: Velocity-integrated brightness temperature for each pixel in the face-on view of our lmde model, expressed in K km s −1 units.Also annotated on the plot is a translucent white circle, which represents an example beam location.Middle and Right: Stokes  and Stokes  profiles calculated from the annotated example beam, to demonstrate the line morphology and linewidth for a typical envelope location.Also included on the  profile plot is a fit curve using /.For this particular location, we obtain a polarization percentage of ∼1%. and optical depth) where high polarization cells are preferentially located.For our low mass disk envelope model this favorable area corresponds to the regions just beyond the optically thick disk, and for our massive protostar envelope model it mainly corresponds to the regions just surrounding the central low-polarization "hole."Though it is difficult to identify a clearly optimal optical depth that maximizes /, overall it appears that   ∼ 0.1 − 1 tends to produce the best conditions for a protostellar envelope Zeeman experiment.
(iii) For both model types, there are significant portions of the respec-tive envelopes that produce percentage polarizations that in principle are accessible with current instruments (i.e., above the nominal ALMA limit, see Table 2).For our low mass disk envelope model, we found the percentage of pixels above 1.8% fractional polarization to be between ∼20 − 40% (depending on line-of-sight).This increases to ∼40 − 60% when the low intensity ambient material outside  = 75 au is excluded.Meanwhile, for the massive protostellar envelope model the percentage of pixels above 1.8% fractional polarization is ∼15% overall, and ∼30 − 40% when we consider only pixels within  = 2500 au.Furthermore, both models have > 90% pixels above 0.18% polarization, meaning that a factor of ten improvement in sensitivity to the magnetic field strength beyond the current ALMA limit would in principle make essentially the entirety of these envelopes accessible to Zeeman experiments.One caveat that should be noted, however, is that in this work we study emission maps produced directly from radiative transfer software.This type of simulation is useful for offering theoretical guidance, but a true direct observational comparison would require producing synthetic observations [e.g., using the CASA observing tool (McMullin et al. 2007)], which is outside the scope of this work.(iv) We find that the percentage polarization in a given pixel is positively correlated with the density weighted average line-of-sight magnetic field strength (first column of Figures 7 and 8).For each of our simulations the associated mean trend line passes through a polarization percentage of 1.8% at roughly 3.1 mG, the reference value we predicted using Eq. 3. (v) We find that percentage polarization is positively correlated with the magnitude of the Stokes  signal, especially at large | | (third column of Figures 7 and 8).This suggests that the regions that are most favorable in terms of polarization percentage are also the most favorable in terms of raw Zeeman signal.(vi) Some regions of high polarization in our massive protostar model are spatially correlated with the location of the shell we included to simulate a cavity from a protostellar wind (see first and last columns of Figure 6).Such sites are predicted to be locations of CN enhancement due to far-ultraviolet dissociation chemistry, and may therefore be good targets for Zeeman experiments with CN. (vii) Stacking the seven main sub-components of the CN  = 1 − 0 transition in our low mass disk envelope simulation yielded an average signal boost of a factor of about 8, with some pixel-to-pixel variation (between about ∼6-10 for the Stokes  and ∼7-25 for the Stokes ).Stacking affects the percentage polarization results only modestly, decreasing / in some parts of the optically thick disk and increasing / in parts of the envelope by up to ∼20% (Figure 12).(viii) Convolution of the emission from our massive protostellar envelope with a  FWHM = 0.5 ′′ beam resulted in a narrower distribution of observed / values, reducing the fraction of pixels with percentage polarizations above 1.8% (Figure 11).Both the unconvolved and convolved maps have an average / ∼ 0.8%.These results indicate that high-resolution observations offer significant advantages to detection prospects for Zeeman experiments in protostellar envelope environments.

Figure 1 .
Figure 1.Top row: Column density plots for each of the Cartesian axis views of our lmde model.Middle row: 2-dimensional histogram plots of local magnetic field component strengths (in the final column,  = √︃  2  +  2  +  2  ) versus distance from the central sink particle.Annotated in red are lines corresponding to 3.1 mG and 0.31 mG line-of-sight magnetic field strength (Section 3 discusses why these values are highlighted).Bottom row: 2-dimensional histogram plots of local magnetic field component strengths versus local density.

Figure 2 .
Figure 2. 3-dimensional view of magnetic field lines in the inner protostellar envelope region of our lmde model.Also included are 3D contours of the dimensionless plasma- to highlight the flattened dense demagnetized structures that dominate the dynamics of the inner envelope [see Tu et al. (2023) for detailed discussion].

Figure 3 .
Figure 3. Top row: Column density plots for each of the Cartesian axis views of our MSEnv model.2nd row: 2-dimensional histogram plots of local magnetic field component strengths (in the final column,  = √︃  2  +  2  +  2  ) versus distance from the central sink particle.Annotated in red are lines corresponding to 3.1 mG and 0.31 mG line-of-sight magnetic field strength (Section 3 discusses why these values are highlighted).3rd row: 2-dimensional histogram plots of local magnetic field component strengths versus local density.Bottom row: 2-dimensional histogram plots of local magnetic field component versus density, only including cells that are given a factor of 1000 enhancement in CN abundance in our model (from  amb = 10 −12 to  shell = 10 −9 ).

Figure 4 .
Figure 4. Top row: Midplane slice plots of the wind tracer value for the MSEnv model, as viewed from the −, −, and − lines-of-sight (left-to-right).Black contours are placed at wind tracer values of 0.001 and 0.2 cm 3 g −1 , the limits that define our CN enhanced shell.Bottom row: Projection plots of the line-of-sight integrated wind tracer value.

Figure 5 .
Figure 5. Velocity-integrated Stokes  and Stokes , line center optical depth   , and derived polarization quantity "Median /" for each Cartesian view of the lmde model.The circle annotated on the last column corresponds to  = 75 au.

Figure 6 .
Figure6.Same as Figure5, but now for the MSEnv model.The circle annotated on the last column corresponds to  = 2500 au.We also draw additional reference lines on the Stokes  panels (in white) and on the Median / panels (in black).These sets of lines are co-located with each other and are drawn in by hand to guide visual inspection of the polarization in the CN enhanced shell and its nearby surroundings.

Figure 7 .
Figure 7. Model lmde 2-dimensional histograms (for each Cartesian line-of-sight) comparing the median / quantity (vertical axis) in each pixel with several other quantities derived from the same pixel.For each pixel these quantities include, from left-to-right, the density weighted average line-of-sight magnetic field strength, the maximum magnetic field strength in any 3D cell along that sightline, the absolute value of the velocity-integrated Stokes  signal, the line-center optical depth, and the distance (in observer space) from the central sink particle.

Figure 8 .
Figure 8. Same as Figure 7, but now for model MSEnv.

Figure 9 .
Figure9.Percentage of pixels with median polarization percentage above 1.8% (blue curve) and average (mean) optical depth (pink curve) calculated in thin radial shells, computed for both of our models as observed along each Cartesian line-of-sight.A horizontal line corresponding to a CN line-center optical depth of unity is also annotated on each plot.

Figure 10 .
Figure 10.Line-center optical depth and median / calculated for several intermediate inclination views of model lmde.Over-plotted on each panel in cyan are contours corresponding to 5 × 10 −4 and 10 −3 Jy/px, which trace out the brightest parts of the central disk.For all viewing angles, these brightest parts of the disk tend to have low polarization and are immediately surrounded by high polarization regions.For views that are close to edge-on (∼60 • -90 • ) the high polarization regions tend to be above and below the disk midplane, and for views close to face-on (∼0-30 • ) the high polarization regions roughly form a ring around the optically thick disk.)

Figure 11 .
Figure 11.Left: Comparison of velocity-integrated Stokes  and , line-center optical depth, and median polarization percentage between pixel resolution and  FWHM = 0.5 ′′ beam convolved cases, for model MSEnv as viewed from the  line-of-sight.The cyan contours correspond to Stokes  values of 5 × 10 −4 and 10 −3 Jy/px.Right: Histogram plots comparing the distribution of each observable for the unconvolved versus convolved cases.The main effect of applying beam convolution is that the / distribution narrows, resulting in fewer pixels having percentage polarization values above 1.8%.The distribution remained peaked at roughly the same value, however (/ ∼ 0.8%).

Figure 12 .
Figure 12.Left: Comparison of the velocity-integrated  and  signals and median polarization percentage calculated from the single 113.144GHz representative transition versus those calculated after sub-transition stacking.Over-plotted contours correspond to   = 1.Center: Pixel-by-pixel ratios of the observables for the stack versus the representative transition.Right: Histogram plots of the same ratios.We generally find that the amount of signal boost gained from stacking has some modest pixel-by-pixel variation, with both distributions centered on a factor of ∼8.
tot        tot