A VLA View of the Flared, Asymmetric Disk around the Class 0 Protostar L1527 IRS

We present high-resolution Karl G. Jansky Very Large Array ( VLA ) observations of the protostar L1527 IRS at 7 mm, 1.3 cm, and 2 cm wavelengths. We detect the edge-on dust disk at all three wavelengths and ﬁ nd that it is asymmetric, with the southern side of the disk brighter than the northern side. We con ﬁ rm this asymmetry through analytic modeling and also ﬁ nd that the disk is ﬂ ared at 7 mm. We test the data against models including gap features in the intensity pro ﬁ le, and though we cannot rule such models out, they do not provide a statistically signi ﬁ cant improvement in the quality of ﬁ t to the data. From these ﬁ ts, we can, however, place constraints on allowed properties of any gaps that could be present in the true, underlying intensity pro ﬁ le. The physical nature of the asymmetry is dif ﬁ cult to associate with physical features owing to the edge-on nature of the disk, but it could be related to spiral arms or asymmetries seen in other imaging of more face-on disks.


INTRODUCTION
Protostellar disks are a natural consequence of conservation of angular momentum during the star formation crocess when the natal cloud collapses to form a young star (e.g.Terebey et al. 1984;Ulrich 1976).Such * NSF Astronomy & Astrophysics Postdoctoral Fellow disks are thought to be established early in this process (e.g.Tobin et al. 2012;Eisner 2012;Sheehan & Eisner 2017a;Tobin et al. 2020;Garufi et al. 2021) and are important for setting the stage for planet formation (e.g.Sheehan & Eisner 2017a;Sheehan 2018;Tychoniec et al. 2020;Segura-Cox et al. 2020).Recent high resolution imaging of protoplanetary disks, the more evolved siblings of protostellar disks, has uncovered a diversity arXiv:2206.13548v1[astro-ph.SR] 27 Jun 2022 of "substructures", features that deviate from an otherwise smooth, monotonically decreasing intensity profile.These substructures typically come in the form of narrow rings and gaps (e.g.ALMA Partnership et al. 2015;Isella et al. 2016;Andrews et al. 2016;Long et al. 2018;Huang et al. 2018a;Cieza et al. 2021;Sierra et al. 2021) that are typically associated with the depletion of dusty material in certain regions of disks, but also as large scale asymmetries (e.g.van der Marel et al. 2013;Casassus et al. 2013;Boehler et al. 2018;Cazzoletti et al. 2018;Dong et al. 2018) and spiral arms (e.g.Pérez et al. 2016;Huang et al. 2018bHuang et al. , 2021)).
Disk substructures are frequently tied to the presence of planets hiding in the disks, and carving the material (e.g.Dodson-Robinson & Salyk 2011;Kley & Nelson 2011;Zhu et al. 2012;Dong et al. 2015;Zhang et al. 2018), leading many to use these features as signposts of young planets.Few planets have, however, as of yet been found hiding within substructures directly (e.g.Keppler et al. 2018), though the presence of some have been inferred through kinematic features (e.g.Teague et al. 2018;Pinte et al. 2018;Isella et al. 2019;Pinte et al. 2020).As such, a number of other mechanisms have therefore been proposed to explain their presence (e.g.Cuzzi & Zahnle 2004;Zhang et al. 2015;Flock et al. 2015;Okuzumi et al. 2016;Suriano et al. 2018;Takahashi & Muto 2018;Ohashi et al. 2021).Regardless of whether they are carved by planets or other mechanisms, substructures likely represent current over-densities of dust that may be conducive to further planet formation within.
The ubiquity of substructures in more-evolved disks (e.g.Huang et al. 2018a;Andrews et al. 2018;Long et al. 2018) begs the question of when substructures first arise in disks.As a result of the connection between substructures and planets or planet formation, a few corollaries to this question are when planet formation begins, and at what time might planets or planetary embryos be hiding in disks.Some evidence exists that such features may be present in Class I protostellar disks (Sheehan & Eisner 2017b, 2018;Sheehan et al. 2020;Segura-Cox et al. 2020;de Valon et al. 2020), but whether substructure can be present in the earliest phase of starformation, the Class 0 phase (e.g.Andre et al. 1993) is as of yet uncertain.Though spiral arms have been identified in such young disks (e.g.Tobin et al. 2016;Lee et al. 2020;Takakuwa et al. 2020), those sources are multiple systems, and the features are likely the result of companion formation in gravitationally unstable disks, gravitational interactions with a circum-multiple disk, or accretion from the envelope onto the disk, rather than the direct impact of ongoing planet formation.
In this paper, we present new, high resolution and high sensitivity observations of L1527 IRS taken with the National Science Foundation's Karl G. Jansky Very Large Array (VLA) at 7 mm, 1.3 cm, and 2 cm, as a part of the Early Planet Formation in Embedded Disks with the VLA (eDisk@VLA) program, a companion to the forthcoming Atacama Large Millimeter Array (ALMA) Large Program of the same name (Ohashi et al., in prep).The goal of the eDisk Large Program is to conduct a systematic search for substructures in young ( 1 Myr-old; e.g.Evans et al. 2009;Dunham et al. 2015) disks to search for evidence of the early onset of planet formation.The companion VLA program is meant to image many of those same sources at long wavelengths in order to help constrain dust grain properties and better characterize substructures at wavelengths where optical depth should be less significant.
L1527 IRS is a well-known edge-on (e.g.Ohashi et al. 1997;Tobin et al. 2008), single (Nakatani et al. 2020, c.f. Loinard et al. (2002)) Class 0 protostar in the Taurus star-forming region and has been extensively studied across a range of wavelengths.Though formally classified as a Class 0 protostar (e.g.Kristensen et al. 2012), it is difficult to tie this classification to an evolutionary stage for the system due to the edge-on geometry of the disk (e.g.Crapsi et al. 2008).Though ages of protostars are notoriously difficult to measure, a more physically motivated system for classifying the evolutionary stage of a system, such as considering how the envelope mass compares with the protostellar mass (e.g.Robitaille et al. 2006), suggests that L1527 IRS is indeed quite young, even if it is not among the youngest of the Class 0 systems (see e.g.Tobin et al. 2013, and references therein).
Early observations with the VLA suggested that the disk is likely optically thick out to at least 1 mm, but that there could be significant optically thick material even at wavelengths of ∼ 1 cm (Melis et al. 2011).It was the first Class 0 protostellar source identified to have a Keplerian-rotation-supported disk (Tobin et al. 2012;Ohashi et al. 2014).Aso et al. (2017) estimated a protostellar mass of 0.45 M and a disk radius of 74 au from further ALMA kinematic observations.The disk is also warm, with midplane temperatures exceeding 20 K out to 75 au (van't Hoff et al. 2018) based on the presence of CO emission at large radii, and molecular line observations also show interesting features near the disk-envelope transition (e.g.Sakai et al. 2014a,b).Finally, and most relevant to the work presented here, it was recently suggested to have substructures in the form of three clumps spread across the disk (Nakatani et al. 2020) from earlier observations with the VLA at 7 mm.The structure of this work is as follows: we describe the observations and data reduction in Section 2. In Section 3, we present our observations and perform careful modeling using analytic intensity profiles to characterize the structure of the disk at each wavelength included in our observations.Finally, we discuss the implications of the disk structure found as a result of our modeling in Section 4, and summarize our results in Section 5. We further apply the same modeling framework to the observations of L1527 IRS at 7 mm reported in Nakatani et al. (2020) to test how our results compare with that work in Appendix A.

OBSERVATIONS & DATA REDUCTION
L1527 IRS was observed by the VLA in four epochs between 7 January 2021 and 14 January 2021 (Program 20B-322) in the A-configuration, with baselines ranging from 680 m -36.4 km.The pointing center for the observations was set at α(J2000) = 04 h 39 m 53.9 s δ(J2000) = 26 • 03 09.6 based on van't Hoff et al. (2018).The observations were taken using the Q-band (44 GHz, 7 mm), K-band (22 GHz, 1.3 cm), and Ku-band (15 GHz, 2 cm) receivers.The Q-band and K-band observations were taken during the same epochs using band switching and were observed three times on 7, 8, and 12 January 2021, while the Ku-band observations were taken during a single epoch on 14 January 2021.The Q-band and K-band receivers were configured in wideband continuum mode, with four 2 GHz wide basebands centered at 41,43,45,47 GHz and 19,21,23,25 GHz, respectively.The Ku-band receivers were also configured in wideband continuum mode but with only three 2 GHz wide basebands centered at 13, 15, 17 GHz.
The science target was observed along with the quasars 3c147 as the flux calibrator and 3c84 as the bandpass calibrator for all bands.For complex gain calibration, the quasar J0440+2728 was used as the calibrator for Q and K-band, and J0403+2600 was the phase calibrator in Ku-band.The total time on-source for L1527 IRS was 190 minutes in Q-band, 50 minutes in K-band, and 18 minutes in Ku-band.The data were reduced, including flux, bandpass and phase calibration along with automatic flagging for radio-frequency interference (RFI) using the VLA pipeline in the CASA software package (McMullin et al. 2007), version 6.1.2.If additional flagging was found to be necessary after a pipeline run, the necessary flags were applied to the data and the pipeline was re-run.
For all bands, we image the data using the tclean routine within CASA using multi-frequency synthesis mode and a robust parameter of 0.5.We did repeat the imaging using a range of robust parameters and found that 0.5 produced the best balance between resolution and sensitivity, but also note that the bulk of the remainder of our analysis is done in the uv-plane and is therefore unaffected by this choice.The resulting Q-band continuum image has a beam size of 0. 045 × 0. 043 with a position angle of 51.6 • and an RMS of 8.8 µJy beam −1 .The K-band continuum image has a beam size of 0. 087 × 0. 084 with a position angle of -58.0 • and an RMS of 4.6 µJy beam −1 .Finally, the Ku-band continuum image has a beam size of 0. 123 × 0. 120 with a position angle of -89.5 • and an RMS of 5 µJy beam −1 .We show these images in Figure 1.We note that there is a small systematic spatial offset of ∼ 0. 035 between our Q and K-band data with respect to the Ku-band data.This is due to using a different phase calibrator at Ku-band data relative to Q and K-bands.This offset does not, however, affect our analysis as we consider each wavelength separately.
As much of our analysis will be done in the uv-plane, we check the uncertainties on the visibilities by comparing σ vis = 1/ Σ(1/σ 2 i ), where 1/σ 2 i = w i is the weight for an individual integration as given by the VLA pipeline, with the RMS of a naturally weighted image, and scale the weights on the visibilities until they match.Though it is not necessarily the case that σ vis must match the RMS of a naturally weighted image exactly, as systematic uncertainties and imaging artefacts could affect this comparison, they should be in reasonable agreement.This results in scaling the weights on the visibilities, w i , by a factor of 0.125.We note that the correction factor is small, which also increases the uncertainties, and is, therefore, a conservative estimate of the uncertainties on the data.
Finally, there is an additional uncertainty on the overall flux calibration of the data, typically on the order of 10% of the flux.We do not, however, consider this uncertainty in our analysis as it only affects the flux scale of separate observations relative to each other, and does not have an impact on the uncertainties within an observation, i.e. on source structure.Any flux density uncertainties reported within this work, however, are purely statistical, so one should include an additional 10% uncertainty on reported fluxes when used beyond this work.

ANALYSIS & RESULTS
In agreement with previous observations (Ohashi et al. 1997;Loinard et al. 2002;Tobin et al. 2008Tobin et al. , 2010Tobin et al. , 2012Tobin et al. , 2013;;Sakai et al. 2014b;Aso et al. 2017;Nakatani et al. 2020), our new data (shown in Figure 1) of L1527 IRS show an edge-on disk elongated in the North-South direction.This can be seen clearly at 7 mm, but the To the right of each image we show the intensity profile from a one-dimensional pixel slice in the North-South direction through the center of the disk, along with a shaded region representing the 3σ confidence interval.The edge-on disk is clearly visible at 7 mm but also can be made out at longer wavelengths, as well.An asymmetry in the disk, with the southern side brighter than the northern side, can also be made out in the all three images, though is most prominent at 7 mm and 1.3 cm.
North-South disk can also be seen peeking out from the bright central point source at 1.3 and 2 cm.The extent of the disk at 7 mm is consistent with previous 7 mm observations by Loinard et al. (2002) and Nakatani et al. (2020), while it is significantly smaller than the disk extent at shorter millimeter wavelengths (e.g.Sakai et al. 2014b;Aso et al. 2017;Nakatani et al. 2020).This discrepancy is likely due to the radial drift (e.g.Weidenschilling 1977) and/or preferential growth (e.g.Birnstiel et al. 2010) of the larger dust grains traced by these long wavelength observations.The bright central emission seen at all three wavelengths is likely a combination of dust and free-free emission from the protostellar jet, with the fraction increas-ing with wavelength.Indeed, a weak East-West protrusion can be made out in both the 1.3 and 2 cm observations.Moreover, the spectral index of the central peak in each of our images, as calculated from a simple two-dimensional Gaussian fit limited to the central region of the emission at each wavelength, is consistently < 1.5 and is equal to ∼ 0.8 between our longest wavelength images, indeed suggesting some amount of freefree emission mixed with the dust thermal emission.
We find that in both the 7 mm and 1.3 cm images, the disk appears to be asymmetric with the southern half of the disk brighter than the northern half.At 7mm the difference in intensity is approximately 3σ when measured 0. 1 north (204 µJy beam −1 ) and south (232 µJy beam −1 ) of the emission peak, though the emission is resolved so this may be more significant when integrated over larger regions.At 1.3 cm the asymmetry is more pronounced, with an intensity of 190 µJy beam −1 measured 0. 1 south of the peak and 98 µJy beam −1 to the north, a difference of ∼ 20σ.This asymmetry may also be present in the 2 cm observations, particularly visible on the southern side of the disk, but is more difficult to pick out by-eye.In contrast with Nakatani et al. (2020), we find that the disk has no apparent large, clumpy structures; as our new observations have ∼ 10× higher sensitivity than the images previously presented, such features should have been present at the > 10σ level.Instead, we find that the northern and southern clumps previously identified correspond to the northern and southern shoulders of the disk in our imaging.For interested readers, we compare our results with their work in greater detail in Appendix A. On the southern half of the disk at 7 mm, around ∼ 0.1 south of the central peak, there is a weak feature that might deviate from the smoothly decreasing emission.From the images, however, it is not immediately obvious whether this feature is real, and if so, what it might represent.If we consider the difference between the lowest and highest values found near the feature in the image plane, then the significance is < 3σ, and so quite tentative.
To better characterize the disk emission at all three wavelengths and determine whether such features seen in these images are real in a statistically rigorous way, we fit analytic models to our observational data.We describe the model, fitting procedure, and results below.

Analytic Models
To characterize the brightness distribution of L1527 IRS in each image, we fit analytic models to the visibility data to try and reconstruct the image.The goal of this modeling is not to necessarily provide a fully physically motivated model, such as a radiative transfer model (e.g.Sheehan & Eisner 2017a), but rather to produce a simple model that can reproduce the features seen in the image such that we can test whether those features are statistically significant in a more quantitative way.We provide further interpretation of these model features in Section 4. We start with a simple model but add additional components motivated by the features seen in each image.
To find the best-fit set of parameters for each model, we use the dynesty package (Speagle 2020) to sample the posterior distributions using Nested Sampling (Skilling & John 2004;Skilling 2006).Nested Sampling is a method for estimating the Bayesian Evidence for a model.Bayesian Evidence, or the marginal or integrated likelihood, is given by where P (D|Θ, M ) is the likelihood of the data given the parameters Θ in the model M , P (Θ|M ) is the prior for the parameters in the model, and Ω Θ represents the entire parameter space (for further details see Speagle 2020).The Bayesian Evidence is an important tool in model selection, as it provides a quantitative way to compare models through the Bayes Factor, or the ratio of Bayesian Evidences, i.e.Bayes Factor = Z M1 /Z M2 .Nested sampling calculates the Bayesian Evidence by sampling randomly from increasingly small nested shells of constant likelihood, and integrating the prior over these nested shells. 1 By using Nested Sampling, we can simultaneously sample the posterior for our models, but also calculate the Bayesian Evidence to allow us to quantitatively compare the quality of fit of our different model choices.
Our base model is a rectangle, motivated by an inspection of the data in Figure 1, in which the 7 mm data appears to be broadly rectangular in outline.A simple rectangular model would be that of a two-dimensional top hat, but it can have significant ringing due to the sharp edges, so we introduce an exponential taper in both directions to smooth the profile and minimize such effects.The analytic prescription for this rectangle is therefore given by, where the x coordinate is defined to be along the major axis of the rectangle, and the y coordinate is defined to be along the minor axis of the rectangle.In our base model, we fix γ x = γ y = 4 as this prescription provides a sharper truncation and a more rectangular appearance than a more traditional two-dimensional Gaussian function (γ x = γ y = 2), which would appear more like an oval.We do, however, allow both to vary, independently, in subsequent fits to allow for the possibility of smoother or sharper truncation.We also vary the centroid of the rectangle (x 0 , y 0 ), and the width of the rectangle along the major (x w ) and minor (y w ) axes.Finally, though the emission is close to due North-South, we also allow the position angle (p.a.) of the rectangle to shift to match the slight offset.
To reproduce the North-South emission profile, which is sharply peaked at the center and not uniform across the extent of the disk, we add a power-law brightness profile along the major axis.We initially use a single power-law, where γ is the power-law slope and is allowed to vary.To smooth between a point-source-like central component, where I ν −→ ∞ as x −→ 0, or a broader central peak, we add a small constant value to x, such that x = x + δ.
We initially fix δ = 0.1 dx, where dx is the pixel size in the model image, but also allow it to vary in later fits.
After an initial round of fitting, we also explored using a smoothly broken power-law profile, ) to better fit the central peak.Here, γ in and γ out are the inner and outer power-law slopes and x b is the point where the transition from inner to outer slope occurs.∆ controls the transition smoothness, with small values indicating a sharp transition and large values a slow, smooth transition.
To test whether the disk is asymmetric, we allow the model to have differing power-law indices for x > 0 (south, towards the excess emission) and x < 0 (north).For the power-law model, this means that instead of γ, we have γ + and γ − parameters.For the broken-power law model, because the central peak appears to be more or less symmetric, we fix γ in,+ = γ in,− , but allow the outer power-law indices to vary independently as γ out,+ and γ out,− .
We also note that the disk appears to be flared in the 7 mm image.To model this flaring, we allow the width of the rectangle in the y-direction to vary as a function of position in the x-direction, Here, A controls how much wider the disk is at x = x w compared with at x = 0, with the width scaling linearly along the major axis.Though our observations presented in Figure 1 ostensibly show a disk with no clear gaps in the intensity profile, with the exception of the tentative feature on the southern side of the disk, visibility data can encode information at smaller spatial scales than are recovered by CLEANed images (e.g.Jennings et al. 2022).Moreover, clumps that appear gap-like were previously reported in 7 mm imaging of L1527 by Nakatani et al. (2020).Therefore, to search for substructures, we add a gap to the model.To prevent ringing in the model from sharply truncating the intensity in the gap, we use a smooth gap model parameterized as a Gaussian subtracted from the gap-free model.
(6) Here x gap represents the center of the gap, ∆ gap the multiplicative factor by which the intensity is reduced at the center of the gap, i.e., and w gap the width of the gap.We did also consider other prescriptions for the gap, such as a simple one in which the density is reduced by ∆ gap within |x −x gap | < w gap /2, and found consistent results regardless of our exact choice of how to represent such a feature.We consider models with just a single gap on one half of the disk, motivated by the feature seen in the southern half of the disk in the 7 mm image, but also a model in which the gap feature is symmetric across the center of the disk.Rather than fit for the peak intensity, or intensity normalized to a specific location in the disk, we instead integrate the emission from the entire model over all space and re-scale the model image to a given total flux, F ν .We then, in our model fitting, use F ν as a free parameter to report the total flux in each image.
Finally, in our initial fits we found that the shortest baseline data could typically not be fit well with the models that reproduce the disk emission alone.This is likely due to the presence of a low surface-brightness envelope around a young protostar that cannot be otherwise seen below the noise in the image.To account for this emission in our model, we also add a simple, large scale envelope component represented by a symmetric, two dimensional Gaussian with total flux F ν,env and width σ env .Though this is a relatively simple parameterization of what might otherwise be a more complex structure (e.g.Ulrich 1976), the main intention is to provide a reasonable approximation of the large scale data so the "disk components" of the model need not try to fit such emission and can focus on the disk.
In the list above, we report only broken power-law models (e.g.Equation 4; or in the case of Model 1, Equation 2 with no power-law) because these consistently fit the observations better than a single power-law model, with typical Bayes Factors of ∼ 6 − 45 in favor of the models with the broken power-law intensity profile.We do, however, note that the choice does not impact any of our conclusions.We also note that for simplicity, throughout much of the remainder of the text we will refer to these models by their number in the above list rather than the full name describing the model, though we will also mention the pertinent features of the model along with the number as well.
We provide a summary of all of the analytic model parameters, including a short description of what they represent, in Table 1.We note that a number of parameters could span orders of magnitude in value, and we fit the base-10 logarithm of the parameter rather than its linear value.For most parameters, we assume a simple uniform prior that limits the value to a reasonable range, but do put further restrictions on some.In particular, we require that the break in the power-law slope occurs "inside" the rectangle, i.e. x b < x w , and that the gap falls between the break and the edge, when present.We also require that γ in > 0 to represent the steep increase in brightness at the center of the disk, but only require that γ out > −1 to allow for a brightness profile that is flat or even increasing with radius over some portion of the disk.We also require that for the asymmetric models γ + > γ − so that the posterior is not bi-modal with the brighter half of the disk occurring at either positive or negative x values.When considering size scales in the model, in particular x w , y w , and x b , we use 0. 005 as a lower limit on those sizes as features on such scales should be well below the resolution of our observations  and therefore difficult to distinguish.For w gap we use a smaller limit of 0. 001 as the Gaussian gap is actually wider than this by a factor of a few.We include the priors for each parameter in Table 1.
We fit these models to the data directly in the twodimensional visibility plane, where the errors are best calibrated.To do so, we generate models in the image plane with 1024 2 pixels that are smaller than the pixels of the observed images in Figure 1 by a factor of 4 (i.e.dx = 0.25 dx image ).We then use the galario package (Tazzari et al. 2018) to Fourier transform the model image into the visibility plane and sample at the baselines of the observations.

Analytic Modeling Results
The results of our model fitting are presented in Table 2. Though we hesitate to define a "best" model, as our analysis does not always provide a singular best model, for the purposes of presenting a reasonable amount of information we select a "reference" model from among our model fits.To select our reference model, we choose the simplest model for which the Bayes Factor compared with the previous model indicates strong evidence for the new model (Bayes Factor > 2.5 following the Jeffrey's scale) with 3σ significance.We list the parameters for that representative dataset in Table 2, calculated as the peak of the marginalized posterior for each individual parameter with the uncertainties representing the range around this peak that contains 68% of the posterior samples.We also show an image of that model compared with the observations in Figure 2 using the maximum-likelihood parameters from the fit for that model.We also list the Bayes Factors calculated relative to this representative model in Table 2.
We find that at 7 mm, Model 4, which includes a broken power-law profile with an asymmetry and flaring, provides the last significant increase in Bayesian Evidence.Models that include an asymmetry (Models 3+) provide a substantial improvement in the quality of fit as compared with models that do not have an asymmetry in the disk brightness (Models 1/2).This improvement demonstrates that the asymmetry that could be made out visually is indeed a statistically significant feature in the data.Including the parameters that control flaring  A comparison of our observations of L1527 IRS with our reference model at each wavelength using the maximum likelihood parameters from the respective fit, with the 7 mm data shown in the top row, the 1.3 cm data in the center row, and the 2 cm data on the bottom.The left column shows the one-dimensional azimuthally averaged visibilities compared with the reference model curve.Though the visibilities are shown averaged radially for ease of viewing, all fits were done to the full two dimensional data.In the middle three columns we show the images of our data, the reference model, and the residuals.The model and residual images and visibilities were generated by Fourier transforming a model image, sampling at the same baselines as the data in the uv-plane using GALARIO, subtracting these synthetic visibilities from the data in the case of the residuals, and re-imaging with a CLEAN implementation built into the pdspy package.Finally, in the rightmost column we show an intensity profile from a one-dimensional slice through the center of both the data and model images.We also show curves for other models with equivalent Bayesian evidence to our reference model in the leftmost and rightmost panels.
of the disk (Model 4) also provides a strong increase in the Bayesian Evidence as compared with models that do not include flaring (Models 1 -3), indicating that the flaring that can be made out in Figure 1 is also real.
We further find that adding either a single gap or a symmetric gap (Models 5/6) produces Bayes Factors that are consistent with zero indicating that we cannot determine whether one of Models 4, 5 or 6 is a better representation of the data than the others.In other words, our observations are perfectly consistent with a gap-free intensity profile, but the presence of gaps in the intensity profile is also allowed.We can, however, use the posterior distributions inferred from fits of mod-els with gaps to characterize the properties of gaps that would be allowed by our observations.In Figure 3, we show the posterior distribution for the width and depth of the gap (w gap and ∆ gap ), marginalized over all other parameters, for the best fit single-gap model (Model 5).To better explore, visually, the range of possible solutions, we also randomly sample 20 models from the posterior distribution and plot the intensity profile for those models along the major axis in Figure 3.We find that solutions with a range of gap widths could be consistent with our observations.Wide gaps are required to be quite shallow, otherwise they would have been detectable by our observations.In fact, for the most shal- (right) 20 realizations of intensity profiles for the L1527 IRS disk at 7 mm drawn from the posterior of the fit, and zoomed in on the gap to see the structure better.We note that the intensity profile is normalized at the location of the gap, before the gap is added, for easier comparison of the gap shape.The solution is quite degenerate, with both narrow and deep or shallow and wide gaps allowed.In the latter case, the feature may be more of a "shoulder", with the emission dropping and flattening out rather than falling to a true local minimum, instead of an actual "gap" in the emission.
low "gaps", the feature is barely even present and may be more of a shoulder that drops and flattens out rather than descending to a local minimum, similar to what has been seen in other high resolution imaging of older disks (e.g.Huang et al. 2018b).Narrow gaps, however, particularly those with widths ∼ 0. 01, below the resolution of our observations, could potentially be quite deep.
At 1.3 cm and 2 cm we find that Model 3, which includes an asymmetry but no further additional model components, provides the last significant improvement to the quality of fit to both sets of observations.Adding additional components or parameters does not result in a statistically significant improvement in the Bayesian evidence.We note that, strictly speaking, adding flaring to the disk (Models 4 -6) improves the quality of the model fit to the 2 cm data substantially.On inspecting those models, however, we found that the flaring of the disk was actually being used to fit the East-West jet feature that can be seen in the 2 cm image.We therefore discard these models but also note that adding the gap feature within them did not increase the Bayesian evidence either, so there does not appear to be evidence for the presence of a gap in the 2 cm observations.This is somewhat unsurprising, however, given the lower resolution and sensitivity to dust emission of these data.
As such, we can confidently say that the disk is indeed asymmetric at all three wavelengths of our imaging.
To summarize, we find that disk-like emission is detected at all three wavelengths, and furthermore, models for which the disk is asymmetric provide statistically significant increases in the evidence for those models.This demonstrates, in a statistically rigorous way, that the disk of L1527 IRS is indeed asymmetric with a brighter southern side, as can be seen by-eye in the images in Figure 1.We find no conclusive evidence, however, that gaps are present in the data.That said, we cannot rule such structures out, but can provide constraints on the sorts of features that might still be consistent with our observations.
The edge-on nature of the disk makes it difficult to disentangle the underlying physical nature of this feature, however.One would typically assume that thermal dust emission at 7 mm is optically thin, and therefore traces dust surface density.If that was the case then this asymmetry would seem to indicate some enhancement of dust surface density on the southern side of the disk.Such a density enhancement could be related to a number of physical mechanisms that have been seen in more face-on images of protoplanetary disks as well as other protostellar disks, such as vortices or one-armed spirals (e.g.van der Marel et al. 2013van der Marel et al. , 2016;;Dong et al. 2018;Boehler et al. 2018;Cazzoletti et al. 2018;Sheehan et al. 2020), that produce pressure bumps and thereby cause dust grains to pile up (e.g.Barge & Sommeria 1995;Birnstiel et al. 2013;Meheut et al. 2012).
That said, with the disk edge-on it is not entirely clear that the disk should be optically thin, even at such long wavelengths, as we may be looking through the entire column of the disk.We do note, though, that the disk asymmetry becomes more pronounced at 1.3 cm than it is at 7 mm, with the brightness a factor of ∼ 1.93× brighter 0. 1 to the southern side than to the northern side at 1.3 cm but only a factor of ∼ 1.13× brighter at 7 mm.As dust is more optically thin at longer wavelengths, this would be consistent with a scenario where the disk is at least partially optically thin at longer wavelengths and we therefore see deeper into the density enhancement there.If the disk is optically thick, the asymmetry could still be related to a dust density enhancement, though it would need to be far enough out to be in the optically thin region.
To further estimate the optical depth of the 7 mm emission, we compare the brightness temperature of the observations with the estimated temperature profile.We use two separate estimated temperature profiles to do so.First we consider the temperature profile expected for a ∼ 2 L protostar (T = (L * /4πσR 2 ) 0.25 ).We find that the brightness temperature of the disk around ∼ 0. 1 = 14 au of ∼ 60 − 75 K falls below the expected temperature of ∼ 125 K, indicating an optical depth of ∼ 0.6 − 0.9.On the other hand, if we consider the temperature profile from from van't Hoff et al. (2018, T = 33 (R / 38.5 au) −0.35 K), which is based on measurements of the temperature of L1527 IRS's disk using CO isotopologues, we find an expected temperature of ∼ 50 K, suggesting an optical depth ∼ 1.We note that the temperature profile from van't Hoff et al. ( 2018) is based on measurements at larger disk radii, where heat-ing from the envelope (e.g.Agurto-Gangas et al. 2019) is important.This profile may, however, be too shallow when extrapolated inwards to smaller radii where direct heating from the protostar becomes increasingly important.As such, this latter value is likely an upper limit on the optical depth.Collectively, it seems likely that the disk is perhaps partially optically thin at 7 mm, and increasingly optically thin at longer wavelengths, consistent with the appearance of the asymmetry across images at these wavelengths.
Another interesting result is that the well-resolved 7 mm emission appears to trace the north-south oriented disk well, indicating that it is dominated by dust thermal emission rather than free-free emission.If this is true, it would imply that a significant amount of relatively large grains already exists in this deeply embedded Class 0 disk since mm/cm-sized grains (sometimes referred to as "astrophysical pebbles") are generally thought to be optimal for producing dust emission at VLA Bands.
We estimate the amount of dust present in the disk using the standard assumption of optically thin dust emission such that the dust mass can be calculated from the millimeter flux, (e.g.Hildebrand 1983): We use a distance of 140 pc (e.g.Torres et al. 2007;Zucker et al. 2019) and a temperature of 51 K typical of protostellar disks (e.g.Tobin et al. 2020) based on a suite of radiative transfer models that find that average protostellar disk temperatures follow T = 43 K(L/L ) 0.25 and a ∼ 2 L protostar (Kristensen et al. 2012).For the flux of the disk, we use the disk flux estimated from the rectangle model (Model 1), of 3.3 mJy.The rectangle model without the power-law component to match the central peak that is likely dominated by free-free emission should provide the best estimate of the disk-only flux from our modeling.The dust opacity is the largest source of uncertainty, as it depends significantly on the dust grain size distribution.We assume that where we adopt a 230 GHz opacity of 2.3 cm 2 g −1 following Andrews et al. (2013), and where small, micron-sized grains typical of the interstellar medium have β ≈ 1.5−2 while grains grown to sizes similar to the observed wavelength have β ≈ 0 (e.g.Hartmann & Lee 2008).We find that, depending on the value of β, L1527 IRS has between 15 and 411 M ⊕ for β = 0 and β = 2, respectively.This is lower than the dust mass found by Nakatani et al. (2020), of ∼ 866 M ⊕ , likely due to the difference in 7 mm dust opacity; for β = 0 our dust opacity is 0.08 cm 2 g −1 while theirs is 0.02 cm 2 g −1 , though the different ways that we treat the temperatures may also play a role.Despite the differences in masses, our recovered fluxes are in good agreement, with Nakatani et al.
(2020) finding a total flux of 3.7 mJy to our 3.6 mJy when including the central peak.The mass found by Tobin et al. (2013) from more careful modeling of 870 µm and 3.4 mm observations, of ∼ 25 M ⊕ is in good agreement with the lower end of our range.
It is interesting to note that the 7 mm emission is vertically extended, with a best fit width of 14 au assuming a distance of ∼ 140 pc (e.g.Torres et al. 2007;Zucker et al. 2019), implying that the grains responsible for its emission have yet to settle to the disk midplane.We would naively expect that the large mm/cm-sized dust grains that are primarily probed by our 7 mm observations should settle to the midplane on a timescale relatively short compared with the age of L1527 IRS.We estimate that the scale height of the gas, calculated as h = c s /Ω at ∼ 25 au with c s using both temperature profiles described previously, is ∼ 2.3 − 3.5 au, depending on temperature profile, or a width of ∼ 5 − 7 au, consistent with the dust extending vertically up to a few scale heights in the disk.The timescale for settling is given by (see, e.g., the Armitage et al. 2015, review article on"Physical Processes in Protoplanetary Disks", Section 7.2): where ρ is the local gas density, ρ m the dust material density, v th the gas thermal speed, Ω K disk rotation angular frequency, R the local radius, and M * the central stellar mass.
The most uncertain quantity is the gas density at R ∼ 25 au near the outer edge of the 7 mm disk.It is constrained by the Toomre parameter It would be difficult for the gas density to go well above the fiducial value of 10 −12 g cm −3 , which corresponds to a Toomre Q value that is already close to unity.As such, any dust settling time we calculate above is likely an upper limit on the true timescale for dust to settle.
For the fiducial values of the gas density and other quantities, the dust settling time is about 6, 000 years for mm-sized grains (and 10 times shorter for cm-sized grains), which is significantly shorter than the time scale for the Class 0 (∼ 1.6 × 10 5 years) and I (∼ 5.4×10 5 years) stages of star formation (e.g.Evans et al. 2009;Dunham et al. 2015).One would therefore expect such large grains to have settled to the midplane unless they are continuously stirred up by some kind of "turbulent" flows in the disk meridional plane.The fact that the large grains emitting at 7 mm do not appear to be settled may indicate the presence of a significant level of disk turbulence, which would be consistent with active accretion that is needed to transport the fast envelope infall through the disk to the central protostar, through mechanisms such as the magneto-rotational instability (e.g.Balbus & Hawley 1991).
Finally, we also note that Bae & Zhu (2018) found that the number of spiral arms driven in a disk by a planet is determined, in part, by the disk's aspect ratio (h/r, where h is the scale height and r is the radius within the disk).They found that higher h/r typically led to only a single spiral in the outer disk.If the relatively large vertical extent of the large grains in this disk, as evidenced by the East-West extent, is indicative of similarly high values of h/r for the gas, this would then be consistent with the presence of a single spiral arm in the disk driving the asymmetry.
It is, of course, quite speculative to assume a planetary origin for the asymmetry, and indeed there are other plausible mechanisms by which such asymmetries might be formed.Single spiral arms are to be expected for massive disks with large aspect ratios where self-gravity is dominant, leading to local gravitational instabilities within the disk (e.g.Kratter & Lodato 2016).Infall from envelope to disk may also drive spiral arms in the disk (e.g.Tomida et al. 2017), or alternatively could incite Rossby wave instabilities that form vortices within the disk (e.g.Bae et al. 2015).Asymmetric "streamers" of infalling material that have recently been found towards some protostars (e.g.Alves et al. 2020;Pineda et al. 2020;Thieme et al. 2022) could also preferentially deposit material asymmetrically into the disk.Interestingly, the ALMA Band 4 observations of L1527 IRS shown in Nakatani et al. (2020) show what appears to be a slight asymmetry on the northwestern side of the disk.Based on its location high in the disk and towards the outskirts, this feature could also be associated with infall, if real.Regardless of the origin, however, such over-densities of material could potentially serve as sites with conditions favorable for planet formation.

CONCLUSIONS
In this work, we have presented new, high sensitivity VLA observations at 7 mm, 1.3 cm, and 2 cm of the Class 0 protostar L1527 IRS.Our observations show an edge-on protostellar disk visible at all three wavelengths, along with a central point-like feature that increases in brightness relative to the disk at longer wavelengths.This central point source is likely a combination of both dust emission as well as free-free emission associated with a jet, with the contribution from the jet increasing at longer wavelengths, and indeed an East-West protrusion can be seen perpendicular to the disk at 1.3 cm and 2 cm.With our new, order of magnitude higher sensitivity observations, we do not find evidence of the clumps reported by Nakatani et al. (2020), instead finding that the disk is gap-free to the limit of our sensitivity and resolution.We do find, however, that the disk is asymmetric at all three wavelengths, with the southern half of the disk appearing brighter than the northern half.
To confirm these features visible in the image-plane, we conduct careful modeling of the observations in the uv-plane.We find that models that include an asymmetry in the disk provide statistically significant improvements to the Bayesian evidence in favor of those models compared with models that do not include an asymmetry at all three wavelengths, indicating that the asymmetry is indeed a real feature of the observations.We also find that at 7 mm, models that include flaring of the disk provide statistically significant increases in evidence in favor of those models.Models that include a gap feature have equivalent Bayesian evidence to gapless models.As such, we cannot rule out that such features might be present in the intensity profile, though our modeling does provide constraints on the properties of such putative features.
The origin of the asymmetry is unclear, and particularly difficult to interpret due to the edge-on appearance of the disk, but could be associated with spiral arms or vortices in the disk, or other features that might produce an asymmetric density enhancement.The large vertical extent of the disk is consistent with simulations of both planet-disk interactions and also gravitationally unstable disks.Such a spiral could, in turn, produce an asymmetry like the one seen here.Infall from envelope to disk could viably produce such a feature, as well.The large vertical extent inferred for the large dust grains in the disk also suggests a significant level of disk turbulence, consistent with active accretion through the disk at early times.
Regardless of the origin of the asymmetry, its presence provides an interesting look at the conditions in a particularly young disk around a single protostar, where such observations are sorely lacking.Further observations like these will be critical for understanding the onset of planet formation in the youngest disks.2020), the "clumpy" features that they find in their image are reproduced qualitatively here.(right) Robust = 0.5 weighted image of L1527 IRS at 1.3 cm using the observational data presented in Nakatani et al. (2020), but centered at the same position as the 7 mm to demonstrate the offset between Q-band and K-band.In both images we show the intensity profile from a one-dimensional pixel slice along the North-South direction through the center of the disk to the right of the image in blue, along with a shaded region representing the 3σ confidence interval.The same one-dimensional slice for our new observations, scaled by a factor of the ratio between old and new beam sizes to account for difference in beam area, is shown in orange for a direct comparison with our newer data.We also show contours at intervals of 5σ (7 mm) and 10σ (1.3 cm) for our new observations with significant transparency so as to not obscure the underlying image.
Though our modeling cannot rule out that gap features are present on either the north or south side of the disk (models with gaps (Models 5/6) have comparable Bayesian Evidence to gapless models (Model 4) for our 7 mm observations) the restrictions that are placed on such features are qualitatively quite different from the large clumps proposed by Nakatani et al. (2020).As the sensitivity of our new image is a factor of 10× higher than those presented previously, and the previously reported clumps were ∼ 2σ peaks above the background of the disk, we should have detected such features with high significance in our data.Such features would likely be relatively wide and deep, i.e. down and to the right in Figure 3, which we can confidently rule out.
Instead, we believe that the substructures previously reported are the result of noise peaks or troughs in a noisy image combined with poorer coverage of the uv-plane.We note that the flux difference between the gap and the clump is ∼ 2σ when measured at the highest intensity value within the gap.The total L1527 IRS disk has an area of ∼ 15 beams, so assuming Gaussian noise statistics we would expect 0.025 × 15 =∼ 0.4 noise peaks (or troughs) with a > 2σ significance within the disk.We would also expect numerous ∼ 1σ peaks or troughs that could work together to create the appearance of clumps.
The clumps seen by Nakatani et al. (2020) were likely further emphasized by the subtraction of the compact central free-free emission seen in their K-band observations when imaged with robust=-2 weighting.Doing so may have enhanced the appearance of the clumps in two ways: first there is a systematic spatial offset of 0. 05 between the Q-band emission and K-band emission in the Nakatani et al. (2020) data.As the astrometry of interferometric images is tied to the position of the gain calibrator, which was different for the Q-and K-band observations presented in Nakatani et al. (2020), systematic offsets in the positions of said calibrators could lead to systematic astrometric offsets between two images of the same source but made with different calibrators.As there is no such offset between our new Q-and K-band images that employed the same phase calibrator between them (and also different from either calibrator used in Nakatani et al. (2020)), the offset is likely due to the different phase calibrators used between the Qand K-band observations from Nakatani et al. (2020).The subtraction of the K-band source performed by Nakatani et al. (2020) was done without correcting for this offset and over-emphasized clump C.Moreover, both our new K-band image (see Figure 1) and also the previous data (see Figure 4) show extended emission from the disk.Thus, at the   2020) at 7 mm with the reference model using the maximum likelihood parameters from the respective fit.The left column shows the one-dimensional azimuthally averaged visibilities compared with the reference model curve.Though the visibilities are shown averaged radially for ease of viewing, all fits were done to the full two dimensional data from Nakatani et al. (2020).In the middle three columns we show the images we made from those data, the reference model, and the residuals.The model and residual images and visibilities were generated by Fourier transforming a model image, sampling at the same baselines as the data in the uv-plane using GALARIO, subtracting these synthetic visibilities from the data in the case of the residuals, and re-imaging with a CLEAN implementation built into the pdspy package.Finally, in the rightmost column we show an intensity profile from a one-dimensional slice through the center of both the data and model images.We also show curves for other models with equivalent Bayesian evidence to our reference model in the leftmost and rightmost panels.
same time, the subtraction of the K-band image also likely subtracted both dust and free-free emission and not only free-free emission.
Finally, to test for the presence of gaps in the intensity profile quantitatively, we repeat our modeling analysis for the previous observations.Following the same conventions regarding the reference model for this fit, we list the results of this modeling in Table 3 and we show a comparison of the reference model with the observations in Figure 5.We find that Model 3, which includes an asymmetry but no further components, provides the last significant strong increase in Bayesian evidence.This suggests that with these older data it would have been possible to confidently identify the disk as having a North-South asymmetry.While the addition of further components does increase the Bayesian evidence in some cases, these increases are typically insignificant when compared with the previous model.The most significant increase as compared with the reference model is Model 5, which includes flaring and a gap, with Bayes factor of 1.56 ± 0.55 that does not provide strong evidence in favor of that model, particularly not when compared with other more extensive models (e.g.Model 4).We cannot, however, rule out models with either flaring, a gap, or both, either.Indeed, the posteriors from gapped-disk models find constraints on the gap features that might be present that are similar to our own observations, though somewhat less constraining due to the lower sensitivity of the observations.

Figure 1 .
Figure 1.Robust=0.5 weighted images of our L1527 IRS VLA observations at 7 mm (top left), 1.3 cm (top right), and 2 cm (bottom left).To the right of each image we show the intensity profile from a one-dimensional pixel slice in the North-South direction through the center of the disk, along with a shaded region representing the 3σ confidence interval.The edge-on disk is clearly visible at 7 mm but also can be made out at longer wavelengths, as well.An asymmetry in the disk, with the southern side brighter than the northern side, can also be made out in the all three images, though is most prominent at 7 mm and 1.3 cm.

Figure 2 .
Figure2.A comparison of our observations of L1527 IRS with our reference model at each wavelength using the maximum likelihood parameters from the respective fit, with the 7 mm data shown in the top row, the 1.3 cm data in the center row, and the 2 cm data on the bottom.The left column shows the one-dimensional azimuthally averaged visibilities compared with the reference model curve.Though the visibilities are shown averaged radially for ease of viewing, all fits were done to the full two dimensional data.In the middle three columns we show the images of our data, the reference model, and the residuals.The model and residual images and visibilities were generated by Fourier transforming a model image, sampling at the same baselines as the data in the uv-plane using GALARIO, subtracting these synthetic visibilities from the data in the case of the residuals, and re-imaging with a CLEAN implementation built into the pdspy package.Finally, in the rightmost column we show an intensity profile from a one-dimensional slice through the center of both the data and model images.We also show curves for other models with equivalent Bayesian evidence to our reference model in the leftmost and rightmost panels.

Figure 3 .
Figure 3. (left) The two-dimensional, marginalized posterior for log wgap and log ∆gap parameters of our best-fit singlegapped model (Model 7) for the 7 mm data.(right)20 realizations of intensity profiles for the L1527 IRS disk at 7 mm drawn from the posterior of the fit, and zoomed in on the gap to see the structure better.We note that the intensity profile is normalized at the location of the gap, before the gap is added, for easier comparison of the gap shape.The solution is quite degenerate, with both narrow and deep or shallow and wide gaps allowed.In the latter case, the feature may be more of a "shoulder", with the emission dropping and flattening out rather than falling to a true local minimum, instead of an actual "gap" in the emission.
Software:CASA (McMullin et al. 2007), dynesty(Speagle 2020), matplotlib(Hunter 2007), corner.py(Foreman-Mackey 2016) , galario(Tazzari et al. 2018) ACKNOWLEDGMENTSWe would like to thank the anonymous referee whose feedback helped to focus and clarify the manuscript.P.D.S is supported by a National Science Foundation Astronomy & Astrophysics Postdoctoral Fellowship under Award No. 2001830.J.J.T. acknowledges funding from NSF grant AST-1814762.Z.Y.L. is supported in part by NASA 80NSSC18K1095 and NSF AST-1910106.M.L.R.H. acknowledges support from the Michigan Society of Fellows.J.K.J. and S.G. acknowledge support from the Independent Research Fund Denmark (grant No. 0135-00123B).L.W.L. acknowledges support from NSF AST-2108794.N.O.acknowledges support from the Ministry of Science and Technology (MOST) of Taiwan (MOST 109-2112-M-001-051 and MOST 110-2112-M-001-031).S.T. is supported by JSPS KAKENHI grant Nos.21H04495 and 21H00048.J.P.W. acknowledges support from NSF grant AST-2107841.C.W.L. is supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF-2019R1A2C1010851). K.T. is supported by JSPS KAKENHI Grant Numbers JP16H05998 and JP21H04487.The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Figure 4 .
Figure 4. (left) Robust = -1.0weighted image of L1527 IRS at 7 mm using the observational data presented in Nakatani et al.(2020).Though not identical to the image presented inNakatani et al. (2020), the "clumpy" features that they find in their image are reproduced qualitatively here.(right) Robust = 0.5 weighted image of L1527 IRS at 1.3 cm using the observational data presented inNakatani et al. (2020), but centered at the same position as the 7 mm to demonstrate the offset between Q-band and K-band.In both images we show the intensity profile from a one-dimensional pixel slice along the North-South direction through the center of the disk to the right of the image in blue, along with a shaded region representing the 3σ confidence interval.The same one-dimensional slice for our new observations, scaled by a factor of the ratio between old and new beam sizes to account for difference in beam area, is shown in orange for a direct comparison with our newer data.We also show contours at intervals of 5σ (7 mm) and 10σ (1.3 cm) for our new observations with significant transparency so as to not obscure the underlying image.

Figure 5 .
Figure 5.A comparison of the of L1527 IRS from Nakatani et al. (2020) at 7 mm with the reference model using the maximum likelihood parameters from the respective fit.The left column shows the one-dimensional azimuthally averaged visibilities compared with the reference model curve.Though the visibilities are shown averaged radially for ease of viewing, all fits were done to the full two dimensional data fromNakatani et al. (2020).In the middle three columns we show the images we made from those data, the reference model, and the residuals.The model and residual images and visibilities were generated by Fourier transforming a model image, sampling at the same baselines as the data in the uv-plane using GALARIO, subtracting these synthetic visibilities from the data in the case of the residuals, and re-imaging with a CLEAN implementation built into the pdspy package.Finally, in the rightmost column we show an intensity profile from a one-dimensional slice through the center of both the data and model images.We also show curves for other models with equivalent Bayesian evidence to our reference model in the leftmost and rightmost panels.

Table 1 .
Summary of Analytic Model Parameters and Priors Center of the model in the East-West direction, with positive x0 to the East x0,guess − 0.3 < x0 < x0,guess + 0.3 y0 Center of the model in the North-South direction, with positive y0 to the North y0,guess − 0.3 < y0 < y0,guess + 0.3

Table 2 .
Best-fit Analytic Model Parameters

Table 3 .
Best-fit Analytic Model Parameters for Archival Data