Evolution and Kinematics of Protostellar Envelopes in the Perseus Molecular Cloud

All protostars in this study are in the Perseus molecular cloud, which has the largest count of Class 0/I protostars of all clouds within 350 pc of Earth (Dunham et al., 2015). For this analysis, we take the nominal distance of Perseus from the Earth to be 300 pc for all of our calculations (Zucker et al., 2018). We use observations from the full MASSES survey, which combines the SMA’s subcompact and extended array configuration (Stephens et al., 2019). The observations predominately used in this study are of the C¹⁸O(2-1) spectral line, which traces the protostellar envelopes (Frimann et al., 2017). The C¹⁸O envelopes were studied alongside CO(2–1) (henceforth CO) outflows, which were also observed using the SMA through the MASSES survey (Stephens et al., 2017). Many protostars in this study are part of multiple systems, as detected and cataloged by the VLA Nascent Disk and Multiplicity (VANDAM) Perseus survey down to a resolution of $\sim$19 au (Tobin et al., 2016). Data on source bolometric temperature and luminosity ($T_{\text{bol}}$ and $L_{\text{bol}}$) and their associated errors were taken from the VANDAM survey (Tobin et al., 2016). Figure 2 shows a typical Class 0 object, Per-emb-5, in the C¹⁸O and CO lines.

Table 1 gives the identifying characteristics of the 54 protostars used in this study, collected from (Tobin et al., 2016) and (Stephens et al., 2018). Of these, 26 are Class 0, 11 are Class I, and seven are at the boundary of these two classifications (denoted as Class 0/I in Table 1 and the rest of this paper). The following analysis does not explicitly use class as a data point of interest or comparison, instead choosing to use $T_{\text{bol}}$ for the “age" of the system (see the following paragraphs for a discussion of $T_{\text{bol}}$). Class 0 and Class I protostars were deliberately chosen because more evolved Class II/III sources (i.e., T Tauri stars) have little or no envelope detectable in the C¹⁸O spectral line.

$T_{\text{bol}}$ is used as a proxy for the age of each protostar (Myers & Ladd, 1993). $T_{\text{bol}}$ is defined to be the temperature of a blackbody which has the same flux-weighted mean frequency as the spectral energy distribution (SED) of the source. While we report the errors in $T_{\text{bol}}$ from Tobin et al. (2016), we do not show or use these errors in our plots. These errors only propagate flux uncertainties, while the dominant error is expected to be due to the incomplete SED coverage (Enoch et al., 2009). These errors only suggest whether the $T_{\text{bol}}$ measurements are significant purely based on flux measurements, but the quoted uncertainties have no bearing on whether we trust the $T_{\text{bol}}$ measured toward one protostar more than the other. When considering incomplete SED coverage, the typical error in $T_{\text{bol}}$ is expected to be about 20% (Enoch et al., 2009).

The age separation between Class 0 and I protostars is typically taken to be at $T_{\text{bol}}$ = 70 K (Chen et al., 1995). However, the observed $T_{\text{bol}}$ can change considerably depending on the inclination of the source. Tobin et al. (2016) gave the approximate classification of each protostar, which we report in Table 1. The lifetime of each protostellar class is highly uncertain. The best estimates based on population synthesis is about 0.15 and 0.3 Myr for Class 0 and I protostars, respectively (Dunham et al., 2015). However, these ages could be considerably shorter if populations are not in steady state ($\sim$0.05 and 0.09 Myr, respectively; Kristensen & Dunham, 2018).

As mentioned, the C¹⁸O(2-1) data comes from the full MASSES data survey (Stephens et al., 2019). The data products are position-position-velocity cubes in FITS format. We specifically used the data cubes mapped using the Briggs’ robust parameter of 1. The pixel size for all maps is 0.4$\arcsec$, and the velocity channels have a width of 0.2 km/s. The angular resolution (measured as the geometric mean of the major and minor axes of the synthesized beam) ranged from observation to observation (Stephens et al., 2019), with resolutions ranging from 1.4$\arcsec$ to 4.0$\arcsec$. The median resolution was 2.6$\arcsec$, corresponding to $\sim$780 au,

We made moment 0 (integrated intensity) maps of each source using the moment() method from the Python package SpectralCube¹¹1See spectral-cube.readthedocs.io. This function calculates the moment identically to the immoments method in Common Astronomy Software Applications package (CASA; McMullin et al., 2007), integrating over velocity channels for each pixel. We then calculated the velocity and line width maps by fitting a 1D Gaussian to each pixel’s spectrum along the velocity axis. We only kept pixels with less than 50% line width errors. To extract useful data from the velocity and line width maps, we needed to isolate the C¹⁸O envelope, so we masked the emission using the astropy.stats.sigma_clip method with a 3-$\sigma$ cut above the RMS noise in the moment 0 map. The moment 0, velocity, and line width maps for the protostars studied in this paper are shown in Appendix A.

To measure the shape and flux of the C¹⁸O envelopes, we fit a 2-dimensional Gaussian model onto each unmasked moment 0 map using CASA. We used the fitting functionality of the CASA viewer() on a polygon region around the C¹⁸O envelope. From this, we obtained the position angle of the elongated axis of the envelope and the full-width half maximum (FWHM) major and minor axes of the Gaussian fit. We only consider fits in which the protostar is located within the FWHM of the C¹⁸O envelope. Some position angles have errors higher than 45^∘ (typically circular envelopes); while we include these fits in our tables and figures, they are not used in any of our statistical analyses. It may be noted that some envelopes have an X-like morphology (see for e.g., Per-emb-2 in Fig. 10) which could affect the fitting. The Gaussian fits, however, are only sensitive to the central emission, and thus the outer emission shape has little to no effect on the fit. This can be best shown by again referring to Fig. 10.

Before calculating the total integrated flux of each source, we removed small regions of emission that were not associated with the protostar’s envelope. We did this by using the remove_small_holes() method in the Python package skimage.morphology²²2scikit-image.org. The combination of the 3-$\sigma$ mask and hole removal we henceforth call the “clean mask." For the remainder of this study, all analysis was done using the “clean" 3-$\sigma$ masked versions of the moment 0, velocity, and line width maps, unless explicitly stated otherwise.

To calculate total integrated fluxes, we used the clean mask on the moment 0 maps and calculated the total integrated flux by summing the flux of each pixel. The peak integrated flux for each source was also taken from these clean mask moment 0 maps by choosing the pixel with the maximum pixel flux value. The errors on the peak integrated flux are taken to be the standard deviation of all background pixels (i.e., pixels that are masked in the 3-$\sigma$ maps), while the errors of the total integrated flux were calculated as the peak integrated flux error multiplied by the square root of the number of beams per envelope area.

From the total integrated flux maps, in principle it is straightforward to calculate the local thermodynamic equilibrium gas mass of the envelope by assuming an excitation temperature, an abundance ratio, and optically thin emission (i.e., a local thermodynamic equilibrium mass). However, the abundance of carbon monoxide in the gas phase in dense and cold circumstellar envelopes is poorly constrained, as it differs significantly from the standard values used in the interstellar medium due to freeze-out onto dust grains (e.g., Yıldız et al., 2013). Indeed, there may not be an accurate canonical abundance ratio for envelopes, as episodic accretion may be ubiquitous and will change the fraction of carbon monoxide gas that will be depleted (e.g., Frimann et al., 2017). Moreover, C¹⁸O may be optically thick, and we do not have an accurate way to estimate the C¹⁸O optical depth with MASSES data. Because of the considerable uncertainties in both the abundances and optical depths, we do not estimate envelope gas masses from C¹⁸O data. Envelope masses from the 1.3 mm continuum are given in Stephens et al. (2018); note that these masses should be updated based on the revised distance estimate of Perseus of $\sim$300 pc, where Stephens et al. (2018) uses a distance of 235 pc.

We follow the same method as Goodman et al. (1993) to calculate the velocity gradients for the C¹⁸O envelopes. We used the clean mask velocity maps discussed above, and applied the curve_fit() method from the Python package scipy.optimize to fit the velocity map with a simple 2-dimensional linear model:

where $\Delta\alpha$ is the displacement in right ascension and $\Delta\delta$ the displacement in declination from the center of the envelope in the image (Goodman et al., 1993). After fitting $a$ and $b$ for each source, we calculated the gradient $\mathcal{G}$ and gradient direction (measured east of north) $\theta_{\mathcal{G}}$:

The units for $a$ and $b$ are the same as that for the gradient, which for our study is km s^-1 pc^-1.

We used the multiplicity of the protostellar systems in order to classify the sources in our sample. We defined three classifications: single, close binary, and medium binary. A close binary is a system with two continuum sources detected by the VANDAM survey within 200 au of each other, a medium binary is a system with two sources greater than 200 au but closer than 3000 au from each other, and any protostars with no other protostars within 3000 au are considered single sources. Close binaries are approximate disk-scale binaries, while medium binaries share (or previously shared) a common envelope.

In §1, we presented the empirical model for envelope evolution proposed by Arce & Sargent (2006). At early stages, the model predicts gas to be entrained along the outflows; nevertheless, it is possible the density of C¹⁸O is too low to accurately trace this effect in some of these Class 0 sources. As the outflows widen and push the gas away from the outflow axis, we would expect Class I protostars to have C¹⁸O envelopes that are elongated perpendicular to the outflow direction. We tested for this evolutionary trend by calculating the axial ratio and angular differences between a few key position angles, namely the C¹⁸O elongation, the outflow axis direction, and the gradient direction.

As mentioned in Section 2.2, we fit 2D Gaussian profiles to the integrated intensity (moment 0) maps and found the major and minor axes of the model for each source. Figure 3 shows the distribution of axial ratios for the sources in relation to $T_{\text{bol}}$, where the error bars are from the elliptical Gaussian fits. (For Fig. 3 and all subsequent figures, the number of sources of each class plotted is mentioned in the caption, where ‘C0’ stands for Class 0, ‘CI’ for Class I, and ‘C0/I’ for Class 0/I.) No significant trends are found in axial ratio, although it should be noted that there are no high-$T_{\text{bol}}$ protostars with near-circular envelopes (axial ratio > 0.85). The shapes of the envelopes can be visually seen in Appendix A, where the 2D Gaussian fits are overlaid on the moment 0 maps for each source. For several protostars, the C¹⁸O moment 0 maps have envelope morphologies that are indicative of gas entrained along the outflow, especially those with an x-like morphology centered on the outflow axis (e.g., Per2, Per20, and Per35).

The total integrated flux over the solid angle of the source is vital information for understanding the relative amount of C¹⁸O in the protostellar system. In §2.2 we discussed the method used to calculate the total integrated flux of the envelopes. Figure 4 shows the total integrated flux with respect to $T_{\text{bol}}$.

We detect a C¹⁸O envelope for every Class 0 (as well as for every intermediate Class 0/I) source in our survey, but only detect envelopes in 13 of 23 Class I sources, for a $57\pm 10$% detection rate. This may indicate that some older protostars have already lost the majority of their natal envelopes. Note that these statistics exclude MASSES targets that were deemed to not be actual protostars (Per-emb-4, Per-emb-39, Per-emb-43, Per-emb-45, Per-emb-59, Per-emb-60, and Per-emb-61) and candidate first cores.

In Section 2.2, we discussed how we calculated the velocity gradients of the envelopes. Analyzing the velocity gradients of C¹⁸O envelopes may give a good picture of how the dense gas is moving and how the outflows are affecting the flow of material in the system. (Note, the velocity maps and gradient fit are shown in Appendix A). Table 2 gives the velocity gradients, gradient directions, and axial ratios for sources with “good" velocity gradient fits. Here “good" fits are defined as those fits that both had a velocity map with a visually clear gradient across the envelope and had errors less than 50% for both coefficients $a$ and $b$ in Eqn. 1. The velocity gradient values range from 2 km/s/pc up to 50 km/s/pc.

In Figure 5, we plot the gradient with respect to $T_{\text{bol}}$, where the error bars come from the least squares fits of the velocity maps. No obvious relationship is found. Though not shown, we also looked at the relationship between the velocity gradient and the total integrated flux and the peak integrated flux. We did not find any significant trends, and there was no relationship between the gradient and flux and the multiplicity of the system.

As shown, we did not find a significant relationship between integrated envelope C¹⁸O flux and $T_{\text{bol}}$. Also, we did not find any relationship between the axial ratio of the C¹⁸O envelope and $T_{\text{bol}}$ (Figure 3). This implies that regardless of mass, at young ages C¹⁸O envelopes have a wide variety of shapes, ranging from nearly circular to having an axial ratio of over 2:1. The lack of a relationship appears to be independent of multiplicity.

Figure 6 plots against $T_{\text{bol}}$ the angular difference between (a) the gradient direction ($\theta_{\mathcal{G}}$) and the outflow axis ($\theta_{\text{OF}}$), (b) the gradient direction ($\theta_{\mathcal{G}}$) and the C¹⁸O envelope elongation position angle ($\theta_{\text{Envelope}}$), and (c) the outflow axis ($\theta_{\text{OF}}$) and the C¹⁸O envelope elongation ($\theta_{\text{Envelope}}$) (see Table 1 for values). The plots are color-coded based on the multiplicity of the system. Figure 6(a) shows a possible trend of decreasing angle difference $|\theta_{\mathcal{G}}-\theta_{\text{OF}}|$ with increasing $T_{\text{bol}}$ in Class 0 sources, but this disappears with the Class I sample. As a test of this trend, we plotted box plots of the angle differences for sources in multiple $T_{\text{bol}}$ ranges (10-30K, 31-50K, 51-70K, and $>$70K) on Figure 6(a). While this hints at a weak trend for sources with $T_{\text{bol}}$$<$70K, there is not enough evidence to support any hypothesis of such a trend. For the $|\theta_{\mathcal{G}}-\theta_{\text{Envelope}}|$ vs $T_{\text{bol}}$ panel (b), the early Class I stages all have low absolute differences, which could indicate aligning of the gradient and envelope as protostars evolve into Class I stage. This trend disappears (or becomes harder to define) in the later Class I stage. Nevertheless, we do not have a large enough statistical sample to further investigate this possible trend.

A possible trend between $|\theta_{\text{Envelope}}-\theta_{\text{OF}}|$ and $T_{\text{bol}}$ is evident in Figure 6(c), with an increasing angular difference between the C¹⁸O elongation and the outflow axis direction with increasing $T_{\text{bol}}$, during the Class 0 phase ($T_{\text{bol}}$ $<$ 70 K), which flattens to an approximately constant value of $|\theta_{\text{Envelope}}-\theta_{\text{OF}}|$ for larger $T_{\text{bol}}$ (Class I phase). To highlight this trend, we plotted the same boxplot analysis on 6(c) as we did in panel (a) using the same $T_{\text{bol}}$ ranges. With increasing $T_{\text{bol}}$, we find that $|\theta_{\text{Env}}-\theta_{\text{OF}}|$ increases towards 90 degrees, supporting the concept of flattening envelopes as the protostars evolve. The data shown in Figure 6(c) are plotted again in Figure 7 with a linear scale in $T_{\text{bol}}$ for clarity of the subsequent model fitting.

To quantify the sharp change in envelope orientation, we need to fit a function that has a location where the curve turns over and asymptotes. One function that can be used is a tanh function:

The turnover point of the fit can be estimated as the point with maximum curvature, $x_{c}$, i.e. where the function $f(x)=\tanh x$ maximizes the value of

For $f(x)=\tanh x$, the point of maximum curvature is at $x_{c}=0.919$, which equates to a turnover temperature of $52\pm 15$ K. This estimate requires a small correction, which must be computed numerically, so that $x_{c}$ is independent of changes in axis scaling. For the tanh function, this correction reduces $x_{c}$ by about 3% (Christopoulos, 2014), which is negligible compared to our fit errors. We define $T_{\text{flex}}$ as the $T_{\text{bol}}$ corresponding to $x_{c}$, i.e.,

Using a tanh function is not unique, but we simply want a method that has a knee transition, as it is clear that the data make a sharp transition from low to high values at a $T_{\text{bol}}$ between the typical value for class 0 and class I.

We found a best fit of

when both SVS13C and Per-emb-36 were removed from the fit data, for a value of $T_{\text{flex}}=53\pm 20$ K. The fit without both of these points has the lowest reduced $\chi^{2}$ goodness-of-fit, although this is mostly due to the removal of the “outlier" sources in the calculation of the test statistic. Here we chose to plot $T_{\text{bol}}$ on a linear scale to highlight the asymptotic behavior of the angle differences of the Class I sources. We will discuss the implication of this fit in §4.1.

Because our data is skewed towards younger sources with lower $T_{\text{bol}}$, the above tanh fit may be similarly biased by the majority Class 0 population. To compensate for this, we used a bootstrapping method to iteratively select random samples of N Class 0 sources with replacement, where N is the total number of Class I sources. (For this purpose, any sources classified as Class 0/I in Table 1 were classified by their $T_{\text{bol}}$ value in relation to the separation value of $T_{\text{bol}}$ = 70 K.) We then fit the combined 2N Class 0 and I sources with the tanh model and calculated $T_{\text{flex}}$ for each fit. We found a $T_{\text{flex}}$ inter-quartile range of $[30,72]$ K and a median of 48 K, which overlaps considerably with the error-bars for $T_{\text{flex}}$ calculated from the simple tanh fit described in the preceding paragraph.

We also color-code each target by multiplicity in Figures 6. The relations for each panel in Fig. 6 do not seem to be significantly different based on the multiplicity properties of each protostar.

In this subsection, we analyze MASSES velocity gradient data in the context of gradient analysis from other surveys of larger scale molecular line data. Specifically, we analyze velocity gradients and specific angular momenta with respect to the effective radii.

We calculated the effective radius $R_{\text{eff}}$, defined to be the radius of a circle with area equal to the projected area of the C¹⁸O envelope. For this study, we took the boundary of an envelope to be at the 3-$\sigma$ contour of its moment 0 map. After tabulating this area $A$, we then calculated the radius by deconvolving the beam:

where $A_{\text{bm}}$ is the area of the synthesized beam. Table 2 lists the effective radii for each envelope with a good gradient fit (see §3.1 for definition of a “good" fit). The radius values range from 0.002 to 0.011 pc, giving us about an order of magnitude range of radii.

The findings from the plot of $|\theta_{\text{Envelope}}-\theta_{\text{OF}}|$ versus $T_{\text{bol}}$ in Figure 7 support the evolutionary process proposed in Arce & Sargent (2006). As the outflow opening angle increases, the amount of C¹⁸O along the outflow axis direction decreases, as most of the denser gas is pushed out and away by the outflows, creating cavities in the parent core. This forms a flattened envelope structure that is mostly perpendicular to the outflow axis, and hence we would expect the angular difference between the envelope elongation and the outflow axis direction to be close to 90° at later ages.

We should note that for low $T_{\text{bol}}$ (i.e., younger sources), the distribution of angular differences between the envelope elongation and outflow axis direction is fairly random. This is somewhat expected since the shape of the envelope at early stages can be highly affected by the initial structure.

Via our tanh fit to the data in Figure 7, we estimated the parameter $T_{\text{flex}}$ (where the fit “turns over" and asymptotes) to be 53 $\pm$ 20 K. This bolometric temperature is similar to the typical evolutionary separation of Class 0 and I sources of 70 K given in Chen et al. (1995). Further, the Monte Carlo tanh fitting method resulted in a $T_{\text{flex}}$ quartile range of [30,72] K and median of 48 K, which again contains the conventional 70 K value for $T_{\text{bol}}$ class separation. We should note, however, that the range of uncertainty for $T_{\text{flex}}$ contains much of the low-$T_{\text{bol}}$ sources in our sample, which indicates that the tanh model fitting is not robust enough to conclude on its own a definitive estimate of a separation $T_{\text{bol}}$ value for Class 0 and Class I sources. Nevertheless, coupled with the visible trend in increased $\Delta\theta_{\text{Env - OF}}$ with $T_{\text{bol}}$ as seen in the boxplots in Fig. 6(c), the tanh fit indicates an increasingly perpendicular envelope in relation to the outflows with increasing $T_{\text{bol}}$. Further study via this method and better measurements of $T_{\text{bol}}$ for protostellar systems will potentially provide a more accurate estimate of a separation of Class 0 and Class I sources by $T_{\text{bol}}$.

For a source with an envelope that has an original spheroidal morphology with its major axis perpendicular to the outflow, it may be hard for the outflow to entrain enough gas along its axis to produce an envelope elongated along the outflow axis. As mentioned in Section 3.1, while C¹⁸O does appear to trace entrained gas for some outflows, due to its low abundance compared to ¹²CO along the outflow it may not be sensitive to the entrained gas for every outflow at the early stages. Thus, the orientation of the C¹⁸O envelope axis may be weakly associated with the outflow in Class 0 envelopes.

In Figure 3, we showed that although C¹⁸O envelopes are typically elongated, there is no trend between their axial ratio and $T_{\text{bol}}$. Together with the results above, this says that despite elongation during Class 0 stage (e.g., due to collapse along magnetic field lines), the direction of the envelope’s elongation can later evolve to be perpendicular to the molecular outflow.

It is important to note that our MASSES sample is inherently biased toward younger protostars, i.e. Class 0 and Class I sources that still have an intact, detectable C¹⁸O envelope. This places definitive limits on the timescales of protostellar evolution that we can study, namely any protostars that have shed their envelopes and moving into the Class II phase of evolution (Figure 1). Nevertheless, any study such as this one that is concerned with the dynamics, morphology, and evolution of protostellar envelopes will necessarily have to look at a sample that is younger than the general population of protostars, and this needs to be addressed. This bias does not affect the strength of the results in this paper, as we are not claiming any evolutionary trends beyond the confines of protostars with envelopes. As a test of our Class 0 and Class I samples, we performed a 2-sample Kolmogorov-Smirnov test and found a KS statistic of 0.472, which rejects the null hypothesis of the two samples being drawn from the same distribution at a value of $\alpha=0.15$.

Intuitively, we expect the velocity gradient of a system to increase with a decrease in radius, as material falls in towards the center and picks up velocity. In Figure 8, we plot the velocity gradient versus the effective radius for MASSES envelopes, and data presented in Goodman et al. (1993) and Chen et al. (2019a, b). The inverse relation between gradient and size is shown to many orders of magnitude of radii for different gaseous structures.

Compared to the relationship found in Goodman et al. (1993) and Chen et al. (2019b), we found a significantly higher slope, $|\mathcal{G}|\propto R_{\text{eff}}^{-0.72\pm 0.06}$. Goodman et al. (1993) found the gradient goes as $R_{\text{eff}}^{-0.4}$ while Chen et al. (2019b) (combining data from Goodman et al. 1993) found $R_{\text{eff}}^{-0.45}$ (if one does not correct the radius as we did above). The difference in these previous studies and our model fit is the availability of a larger sample size, resulting in a wider range of values for $R_{\text{eff}}$ for fitting.

The relationship also does not appear to be greatly affected by multiplicity. However, we do note that all the close binaries tend to lie above our best fit, indicating that close multiples may have higher velocity gradients than the simple trend indicates.

One question that might come naturally from this analysis is what this power law fit physically means. Should conservation of angular momentum hold, we would expect that as an envelope increases in radius, its velocity would decrease as $R^{-1}$. However, this assumes solid-body rotation, which is a very crude approximation for the turbulent, accreting envelope environment. Our result is not consistent with solid-body rotation, and we do not observe conservation of angular momentum over these size scales. Since our result of $R^{-0.72}$ is shallower than the $R^{-1}$ relation, this suggests a dissipation of angular momentum as material falls in towards the protostellar disk. There are many physical processes that dissipate angular momentum from the system, such as bipolar outflows and tension in the magnetic field lines. Further, gas is constantly being accreted onto the disk and being pulled into the envelope from the larger cores, and even the cores are accreting gas from the molecular cloud while the outflows are pushing gas out of the system. Thus, it is very difficult to understand a gaseous structure as its own entity without considering the processes at larger and smaller scales that affect it.

It is also important to note that although the data appear to clearly follow a power law fit, we do not claim that our fit (or any power law, for that matter) accurately describes the physics happening in the protostellar systems. In fact, we would expect different efficiencies in angular momentum conservation for different length scales as different processes become more dominant. For example, while bipolar outflows could potentially disperse gas traced by C¹⁸O and affect the calculated velocity gradients at the envelope scale, they probably do not grossly affect the molecular cloud ($>$10 pc scale).

Rather than measuring a single velocity gradient value for each protostar at an effective radius as we did in this paper, Gaudel et al. (2020) calculated velocity gradients at multiple radii from 50 to 5000 au and fit this relationship for each of their 12 protostellar envelopes. In their results, the relationships between the velocity gradient and radius varies from protostar to protostar, but the median velocity gradient is proportional to radius to the –0.85 power. This value is similar to the $R_{\text{eff}}^{-0.72\pm 0.06}$ relation we find. However, we note that the power-law index measured by Gaudel et al. (2020) varies markedly from protostar to protostar.

In contrast to the velocity gradient, we expect the specific angular momentum ($J/M$) to increase with increasing radius, recalling Equation 9: $J/M\propto|\mathcal{G}|R^{2}$. Noting the power law definition of $J/M$ in terms of $R_{\text{eff}}$, then, it is no surprise and in fact expected that we found a power law relationship in Figure 9. As mentioned before, several previous studies find power-law indices for $J/M\propto R^{q}$ to be between 1.5 and 1.8.

Our analysis, coupled with past results, allows us to study this relationship at a wider range of core sizes, from 1 pc down to $\sim$0.0025 pc (500 au). We find a power law fit of $J/M\propto R^{1.83\pm 0.05}$ (see Figure 9), which lies between the relations for solid-body rotation ($J/M\propto R^{2}$) and for turbulence ($J/M\propto R^{1.5}$) (Burkert & Bodenheimer, 2000). This is to be expected, as the envelope gas accrues the angular momentum of the infalling gas from the core, while processes such as bipolar outflows and infalling gas onto the disk provide turbulence in the system.

Past papers, including Belloche (2013), have posited that at scales smaller than $\sim$5000 au, the specific angular momentum will be constant. This turnover is necessary to couple the turbulent dynamics of the core and envelope with the flat inner disk. Indeed, Gaudel et al. (2020) seems to find this turnover at around $\sim$1000 au. Conversely, Pineda et al. (2019) did not find a break in the power law fits for any of their 3 sources’ radial $J/M$ profiles down to 1000 au. Our data, which contains 7 sources with $R_{\text{eff}}<1000$ au and 25 sources with $R_{\text{eff}}<2000$ au, cannot rule out such a flattening occurring in the region of $\sim$1000 au and consequently cannot confirm nor refute such a break in the data. As mentioned in Section 3.3.2, this leveling out of the specific angular momentum as seen in Figure 9 is statistically similar in our data to a simple power law relation.