IX. Multiband, multiscale dust study of L1527 IRS. Evidence for variations in dust properties within the envelope of a class 0/I young stellar object

Context. Early dust grain growth in protostellar envelopes infalling on young disks has been suggested in recent studies, supporting the hypothesis that dust particles start to agglomerate already during the class 0/I phase of young stellar objects. If this early evolution were confirmed, it would impact the usually assumed initial conditions of planet formation, where only particles with sizes ≲ 0 . 25 µ m are usually considered for protostellar envelopes. Aims. We aim to determine the maximum grain size of the dust population in the envelope of the class 0/I protostar L1527 IRS, located in the Taurus star-forming region (140 pc). Methods. We use Atacama Large millimeter/submillimeter Array and Atacama Compact Array archival data and present new observations, in an effort to both enhance the signal-to-noise ratio of the faint extended continuum emission and properly account for the compact emission from the inner disk. Using observations performed in four wavelength bands and extending the spatial range of previous studies, we aim to place tight constraints on the spectral ( α ) and dust emissivity ( β ) indices in the envelope of L1527 IRS. Results. We find a rather flat α ∼ 3.0 profile in the range 50–2000 au. Accounting for the envelope temperature profile, we derived values for the dust emissivity index, 0 . 9 < β < 1 . 6, and reveal a tentative, positive outward gradient. This could be interpreted as a distribution of mainly interstellar medium like grains at 2000 au, gradually progressing to (sub)millimeter-sized dust grains in the inner envelope, where at R = 300 au, β = 1 . 1 ± 0 . 1. Our study supports a variation of the dust properties in the envelope of L1527 IRS. We discuss how this can be the result of in situ grain growth, dust differential collapse from the parent core, or upward transport of disk large grains.


Introduction
Circumstellar disks orbiting class II young stellar objects (YSOs) are also commonly referred to as protoplanetary disks even though a large fraction of them already display structures thought to be shaped by embedded planets (Andrews et al. 2018).Rings, gaps, and spirals have also been observed in younger (<1Myr) disks (ALMA Partnership et al. 2015, Sheehan & Eisner 2018, Segura-Cox et al. 2020, Nakatani et al. 2020).Several mechanisms have been proposed to explain the formation of structures during the early stages of the disk lifetime, including gravitational instabilities (Takahashi & Inutsuka 2014), disk winds (Johansen et al. 2009, Takahashi & Muto 2018), and the evolution of dust (Okuzumi et al. 2016).Along with these possibilities, the early formation of large planetesimals that gravitationally interact with the disk remains a viable explanation for the observations.
The early formation of planetary embryos is also suggested by the dust mass budget of evolved disks when compared to their younger progenitors (Testi et al. 2014, Ansdell et al. 2016).Manara et al. (2018) estimated that the solid (dust) mass observed in evolved class II disks of the Ophiuchus (Sanchis et al. 2020) and Lupus star-forming-region is 1 or 2 orders of magnitude lower than the mass of the expected exoplanet population.Conversely, the class 0/Is disks of Perseus contain roughly the same mass in solids as the known exoplanets (Williams et al. 2019, Tychoniec et al. 2020.One possible explanation is that the missing mass in the class II disks is the early conversion of class 0/I solids into planetesimals (e.g.Testi et al. 2022, Bernabò et al. 2022, Xu & Armitage 2023).In a separate study, Mulders et al. (2021) suggest that the solid masses in class II disks and in exoplanets might be of the same order of magnitude.Nevertheless, an unrealistic planet formation efficiency of 100% would be required to place the beginning of planet formation during the class II stage.
Direct observational evidence for early-formed, kilometersized planetesimals is out of reach for any radio interferometer and thus this problem remains open.Furthermore, observational constraints on the initial properties of the dust population are needed for simulations to reconstruct the pathways that lead to planetesimal formation, starting from submicronsized interstellar medium (ISM) dust particle interactions in a core accretion paradigm (e.g., Safronov & Zvjagina 1969, Goldreich & Ward 1973, Ormel et al. 2009, Birnstiel et al. 2016, Drazkowska et al. 2022).Even more importantly in this context, in an effort to investigate planetesimal formation at early stages, Table 1.Main properties of our target, L1527 IRS.The age of the object has been roughly estimated by [1] (Tobin et al. 2012) based on the mass loss rate of the source.The mass of the protostar has been estimated by means of kinematical analysis by [1] and [2] (Aso et al. 2017).The bolometric luminosity was reported by [3] (Karska et al. 2018) and might suffer up to a factor two uncertainty due to the high inclination of the disk (e.g., Whitney et al. 2003).The envelope mass was derived in [4] (Motte & André 2001) based on the IRAM 30 m telescope and MPIfR bolometer arrays ' 1.3  < 3 • 10 5 [1]   M * [M ⊙ ] 0.2 [1] -0.45 [1,2]   L bol [L ⊙ ] 1.6 [3]   M env 4200au [M ⊙ ] 0.8 [4]   Cridland et al. (2022) find that differential gas and dust replenishment of a protoplanetary disk from its surrounding envelope sets favorable conditions for planetesimal formation.However, grains that are moderately coupled to the gas (>30 µm) are required to meet the streaming instability conditions that lead to the birth of these planetary embryos (Youdin & Goodman 2005).Moreover, constraining grain sizes during the protostellar collapse phase has proven to be important for understanding the role of magnetic breaking during the early stages of disk formation.Since small grains are the main charge carriers, their distribution is paramount for the coupling of the infalling material with magnetic fields.In turn, such a coupling seems to regulate disk masses and sizes (Zhao et al. 2016, Lebreuilly et al. 2020).
An additional reason to investigate the properties of dust grains in these extended envelopes is to constrain the role that dust grains play as sites -and catalysts -of molecular reactions in diffuse environments.The most commonly accepted pathways to produce complex organic molecules all rely on the presence of icy mantles on dust grains to capture the building blocks of these complex molecules and facilitate reactions among them (e.g., Tielens & Hagen 1982).
In the recent past, several works have attempted to measure the dust emissivity index β at the large scales of class 0/I protostellar envelopes to constrain whether dust growth might be significant at these very early stages of star and planet forma-tion.Many have found surprisingly low β values (≤ 1) and interpreted this result as evidence for early dust growth (Kwon et al. 2009, Shirley et al. 2011, Chiang et al. 2012, Miotello et al. 2014, Le Gouellec et al. 2019, Valdivia et al. 2019).Following similar methods, other works have not found hints of such growth (Agurto-Gangas et al. 2019).And, in their sample study of ten CALYPSO (Maury et al. 2019) class 0/I sources, Galametz et al. (2019) found examples of both relatively low and large β values.
The hypothesis of early grain growth in the envelopes of class 0/Is has been considered to be challenging from a theoretical perspective: growth to millimeter-sized particles seems to require environments characterized by higher density and/or longer timescales than the average class 0/I envelope (e.g., Ormel et al. 2009, Guillet et al. 2020, Lebreuilly et al. 2023, Bate 2022and Silsbee et al. 2022).Using both analytical models and numerical simulations, these authors find that dust particles cannot grow larger than ∼ 2µm in collapsing envelopes.It will be crucial for next simulations to incorporate generally disregarded effects, like the dust back-reaction on the turbulence through gasdust friction and dust-magnetic-field interaction.Early growth is not the only possible explanation for the mentioned observations.Lebreuilly et al. (2020) proposed a scenario in which the differential collapse of dust grains of different mass through the prestellar core leads to a stratification of the dust sizes since larger grains collapse faster.Finally, another possibility could be an uplifting of grown disk grains to the inner envelope by the protostellar outflows and/or jets (Wong et al. 2016, Tsukamoto et al. 2021).
It is thus imperative to characterize dust properties and evolution in the very early stage of star and planet formation.
In this paper we focus on a single source in Taurus, L1527 IRS, (140 pc, e.g., Torres et al. 2007, Zucker et al. 2019), hereafter L1527, to benchmark a dust continuum study across multiple frequencies and physical scales.Due to its distance and brightness, L1527 is one of the best studied class 0/I YSOs.The presence of an edge-on circumstellar disk was first suggested by Bontemps et al. (1996) and Ohashi et al. (1997) based on observations of the orientation of the bipolar outflow.It has been later confirmed by means of Spitzer and Gemini-North scattered light observations as an edge-on dark lane obscuring the central protostar (Tobin et al. 2008 andTobin et al. 2010, respectively) as well as high-resolution Submillimeter Array (SMA) and Combined Array for Research in Millimeter Astronomy (CARMA) interferometric observations of its continuum emission (Tobin et al. 2013).The Atacama Large Millimeter/submillimeter Array (ALMA) C 18 O kinematic detection of the disk was reported by Ohashi et al. (2014) and Aso et al. (2017).Furthermore, this edge-on, warped disk presents asymmetries at ∼ 20 au that are consistent with spiral structures, the physical origin of which has been proposed to be gravitational instability (Sakai et al. 2019, Nakatani et al. 2020, Ohashi et al. 2022, Sheehan et al. 2022).Observations of L1527 also display a bipolar jets and outflows extending to 20000 au perpendicular to the disk plane, and carving a cavity in the collapsing envelope (Bontemps et al. 1996, Ohashi et al. 1997, Hogerheijde et al. 1998, Tobin et al. 2010, Oya et al. 2015).Being this object in a late stage of accretion, its molecular jet and outflow are not very energetic (Podio et al. 2021).Table 1 summarizes the main properties of this young source.
We aim to study the dust continuum emission of L1527 using the Atacama Large millimeter/submillimeter Array (ALMA) and the Atacama Compact Array (ACA) to probe the outer (10 3 au) and inner (10 2 au) envelope, as well as the disk (10 au).We consider observations of our target in four different ALMA Table 2. ALMA observations.For projects with more than one execution block, we only report the total integration time of the project.The rescaling factors are relative to the dataset whose factor is exactly 1.000.Bands (3,4,6,7).Section 2 presents the ALMA and ACA observations that we have used throughout the analysis and the data reduction process.In Section 3, we justify the uv-plane geometrical modeling used for the source.In Section 4, we analyze the spectral and dust emissivity indices of L1527 at different scales.We discuss our results in Section 5, and wrap up our conclusions in Section 6.

ALMA observations
We here summarize the details of the observations used throughout this work.This source has been extensively studied with both the extended and compact arrays of ALMA.Archival data include high-resolution observations designed to detect disk substructure as well as lower-resolution, high-sensitivity observations of the extended envelope.We collect every available ALMA and ACA dataset suitable for our study of L1527, aiming to combine them in the deepest continuum analysis of the envelope so far.A summary of the datasets used in this work can be found in what follows and in Tables 2 and 3.

ALMA FAUST B3 and B6
Here, we present new B6 and B3 data from the ALMA Large Program Fifty AU Study of the chemistry in the disk/envelope system of Solar-like protostars (FAUST, PI: Yamamoto, S.).The main goal of FAUST is to investigate the gas chemical composition of the environments surrounding young, Sun-like protostars including their embedding envelopes, their outflows, and disks (see Codella et al. 2021 for further details).
The observations aimed at L1527 were centered at right ascension α(2000) = 04h39m53.878sand declination δ(2000) = +26 • 03 ′ 09.56 ′′ .B3 observations were run on December 14, 2018 and August 25, 2019 with the 12m array in ALMA configurations C43-3 and C43-6, respectively.The baselines of these observations thus ranged 15-2500 meters.ALMA observed L1527 in the 93-95 and 104-108 GHz ranges.The source was observed for approximately 35 minutes in B3, for a continuum sensitivity of 0.025mJy.B6 observations were performed between October 20th, 2018 and December 15, 2019 with both the 7m and 12m arrays.The baselines of the 7m array ranged from 8.9 to 48.9 meters.The 12m array observations were taken in configurations C43-2 and C43-5, with baselines ranging from to 15-314 meters and 15-1400 meters, respectively.ALMA observed Table 3. ACA observations.The total integration times refer to the sum of the integration times within the entire project, in the cases with more than one execution block.The rescaling factors are relative to the 12m array dataset observed with the 12-m array, whose rescaling factor is 1.000 (cfr Table2).
in the 244-247 and 258-262 GHz ranges with a spectral resolution of approximately 62.5 MHz.The source was observed for approximately 75 minutes in B6, for a continuum sensitivity of 0.026mJy.
The FAUST data were calibrated using a modified version of ALMA pipeline version 42866, using the Common Astronomical Software Applications (CASA, McMullin et al. 2007), version 5.6.1-8.This included a correction for errors introduced by the per-channel normalization of data by the ALMA correlator 1 .Line-free Local Standard of Rest Kinematic (LSRK) frequency ranges were identified by visual inspection and averaged per spectral window, and initial continuum images were produced for each separate ALMA configuration.These were then used as initial models for subsequent per-configuration phase-only selfcalibration, followed by amplitude and phase self-calibration.Great care was taken to ensure that the models were as complete as possible to avoid changing the overall flux density scale of the data when doing amplitude self-calibration.L1527 is sufficiently strong that per-integration phase-only self-calibration was possible, while for amplitude self-calibration per-scan selfcal was used.The per-configuration datasets were then aligned across configurations in both phase and amplitude, again using a self-calibration technique.Corrections to the amplitude scale of up to 10% were found to be needed for some datasets.Improvements in the dynamic range (peak/RMS away from emission) of more than an order of magnitude for the final images were achieved using this technique for setup 1 and 2. The improvement for setup 3 was ∼ 35%.Finally, we imaged the visibilities with the CASA tclean function, using the hogbom deconvolver, a "briggs" weighting scheme with robust parameter set to 0.5.The 1 https://help.almascience.org/kb/articles/what-errors-could-originatefrom-the-correlator-spectral-normalization-and-tsys-calibrationresulting beam in B3 is 0 ′′ .44x0 ′′ .27 wide with PA=-29 • ; while the synthesized beam in B6 is 0 ′′ .42x0 ′′ .28, PA=19 • .We present the new continuum maps in Fig. 1 and Fig. 2.

ALMA Band 3
We gathered ALMA Band 3 (B3) observations of L1527 spanning from 2017 to 2020.Given the different scientific aims of the projects, the data collected with the 12-m array has a resolution between 0 ′′ .09 (project 2017.1.00509.S, PI: Sakai, N.) and 2 ′′ .4 (project 2015.1.00261.S, PI: Ceccarelli, C.).The sensitivity of the 12-m array observations that have been used in this work ranges from 0.03 to 0.06 mJy and the total time on source is ∼210 minutes.The frequency of the side bands ranges from 85 to 115 GHz.The 7-m array has been pointed at L1527 in B3 for a total ∼102 minutes (projects: 2016.1.01541.S and 2018.1.00799.S; PIs: Harsono, D. and Pineda, J.).The observations have resolutions 12 ′′ .0 -13 ′′ . 2 and have been carried out in a range of frequencies from 92 to 105 GHz with sensitivities of 0.3-0.6 mJy.

ALMA Band 4
L1527 has been observed in Band 4 (B4) throughout several months in 2016 and 2018.The data collected with the 12-m array has a resolution between 0 ′′ . 1 (project 2016.1.01203.S, PI: Oya, Y.) and 2 ′′ . 1 (project 2016.1.01541.S, PI: Harsono, D.).The sensitivity of the 12-m array data that have been used in this work ranges from 0.03 to 0.3 mJy and the total time on source is ∼82 minutes.The frequency of the side bands ranges from 138 to 154 GHz.The 7-m array observed L1527 in B4 for a total ∼107 min-

ALMA Band 7
L1527 has been observed in Band 7 (B7) from 2014 and 2018.The data collected with the 12-m array has a resolution between 0 ′′ .075 (project: 2016.A.00011.S, PI: Sakai, N.) and 0 ′′ .5 (project 2011.0.00604.S; PI: Sakai, N.).The sensitivity of the 12-m array observations that have been used in this work ranges from 0.06 to 0.2 mJy and the total time on source is ∼155 minutes.The frequency of the side bands ranges from 218 to 263 GHz.The 7-m array observed L1527 in B6 for a total ∼46 minutes during project 2016.2.00117.S (PI: Yoshida, K.).The observations have resolutions 3 ′′ . 1 and have been carried out in a range of frequencies from 352 to 363 GHz with sensitivity of 1.44 mJy.

Data reduction and self-calibration
Since we work with data taken at several frequencies, we avoid frequency smearing by splitting datasets with a relatively large frequency offset.This is especially important in B3.Thus, we split its datasets into two groups.This way, after calibration, we work with five final frequencies: 88, 100, 141, 249, and 348 GHz.These frequencies roughly correspond to 3.4, 3.0, 2.1, 1.2 and 0.88 mm wavelengths.
The first round of calibrations was performed using the standard CASA pipeline methods provided by the ALMA Regional Centers (ARC).Considering the data were acquired over a period of several years, different versions of CASA were used for these calibrations.The version used for each dataset is listed in Table 2.This first round includes system temperature, phase, amplitude, and bandpass calibration, along with corrections to account for atmospheric water vapor.
We carried out further data reduction and calibration steps using CASA version 6.2.1.7.First, the ALMA cubes for each execution block were inspected and additional flagging was applied when necessary, that is to mask out spectral lines.We are only interested in the continuum emission, thus, to speed up further operations, we channel-averaged the spectral windows of every execution block.We used a common width to level out the S/N among spectral windows and we were careful to avoid bandwidth smearing effects by limiting the averaging based on the maximum baselines and spectral window bandwidth (Bridle & Schwab 1999).Before combining the datasets for the analysis, we performed several additional operations.First, we ran the statwt CASA routine in order to fix the visibility weights.This is especially important for datasets from old ALMA Cycles that have been reduced with CASA 4.2.1 or earlier versions 2 .We then imaged the averaged-channel visibilities with the tclean task and fitted a Gaussian with imfit to pinpoint the maximum of the emission.Next, we used the fixvis function to shift the phase center to the position of the peak of L1527 at the epoch of every observation.Since the target coordinates were slighlty different across datasets, we also set the source sky position (in the metadata) to a common center across execution blocks using fixplanets.Once the coordinates of each dataset were phase-centered Fig. 3. L1527 ALMA 88 GHz (purple), 100 GHz (orange), 141 GHz (green), 249 GHz (blue), 348GHz (red) visibility amplitudes up to 1000 kλ.We also compare the 94 and 231 GHz observations (white points) performed with the Plateau de Bure Interferometer (PdBI) and anlayzed in Galametz et al. (2019).The slope of our observed Spectral Energy Distribution (black line) has been obtained by fitting a line to the fluxes at all wavelengths in each uv-distance bin.While 1.5 < α f it < 2.5, most of the contribution to the flux -and to α -at any uv-distance is due to the disk (see Section 4).and aligned following this procedure, we rescaled the fluxes following the procedure of Andrews et al. (2018): for every frequency, selected a range of uv-distances where the uv coverage of two datasets overlaps, binned the visibilities in that range, and compared them to obtain a rescaling factor where the visibilities overlapped.The rescaling factors were in a 1-15% range for all execution blocks.The rescaling correction was then computed using the task gencal and the rescaling itself was achieved by applying the correction factors with applycal.We report the rescaling factors in Tables 2 and 3.
Based on this model, we then computed the phase corrections on each execution block using gaincal and applied them with applycal.This step aimed to correct phase errors between executions and between spectral windows.We repeated this procedure three times.In the first round, we set the solution interval of gaincal to inf, combined the scans within each execution block and set gaintype = "G", that is the gains were determined for each polarization and spectral window.In the second and third rounds, we combined the spectral windows of each block, shortened the solution intervals to 60 and 10 seconds, respectively, and set gaintype = "T", to obtain one solution for both polarizations.The self-calibration yielded peak S/N improvements in the range 10-500 % when evaluated on the images obtained with the combined datasets.We did not find any appreciable improvement in the noise and signal-to-noise properties for phaseonly self-calibration steps with even smaller time intervals.At this point, we combined the corrected datasets using the concat CASA task.
Since we analyze extended emission, we estimated the effect of the antennas primary beam (PB) correction on the recovered flux.In B7, the ALMA 7m and 12m antennas PB full width at half maximum is 33 ′′ .0 and 19 ′′ .0, respectively.This results in a loss of about 10% of flux at 20 kλ (1000 au scales), for the 12m antennas.This only marginally affects the recovered spectral index we are interested in, increasing it by a factor of <5% even at the largest scales.

L1527 as seen in the visibility space
The final binned visibility amplitudes obtained with our data reduction and self-calibration are illustrated in Fig. 3.We also show that the associated spectral index lies in the range 1.5 < α 0.88−3.4mm< 2.4.This low α is not surprising as the disk emission contribution dominates at every uv-distance (see Section 3), thus the total spectral index is strongly influenced by the optically thick disk (cfr.Section 4).
Among other class 0/Is of the CALYPSO survey, the continuum emission of the envelope infalling on L1527 has been studied in Galametz et al. (2019), hereafter G19, who used PdBI observations at 1.3 and 3.2 millimeter.The data they used spans a range in uv-distances from 20 to 200 kλ at both frequencies.Fig. 3 shows a comparison between the ALMA observations we used and G19's PdBI fluxes (white points).
We cut the ALMA datasets at the shortest common baseline, set by B7 at 9.8 kλ, hence a Maximum Recoverable Scale (MRS) of roughly 26 ′′ .0, or 3500 au at the distance of the source.We aim to more robustly determine the spectral and dust-emissivity indices of the source by extending the frequency range, that dominates the error on α (see Section 4).Furthermore, a denser sampling of frequencies is critical to understand the optical depth of the source across wavelengths, which is vital to constrain where the spectral index is a good proxy for dust grain sizes.On top of enriching the studied frequency range, we used a combination of ALMA and ACA data to study the dust continuum emission with high S/N at scales 2 times larger than G19, toward both shorter and longer baselines.In fact, although our primary attentions are devoted to the envelope, long-baseline visibilities are of great importance if we want to properly model the compact emission.G19 subtracted a constant amplitude (F 200kλ ) from each visibility at the two frequencies, assuming an unresolved disk that would contribute a constant flux at baselines shorter than 200 kλ (∼1 ′′ .0, or 140 au).
Considering the wider uv-coverage of the data we used (9.8-1000 kλ), we can more precisely model and subtract the compact emission by fitting a Gaussian component within a model that separately accounts for envelope and disk (see Section 3).This procedure also makes a significant difference when computing the errors on α and β, which in turn depend on the relative errors of the amplitudes.Since the PdBI amplitudes in G19 were roughly constant and only slightly above F 200kλ , the relative error on the fluxes after the subtraction of F 200kλ became large already at a uv-distance of ∼ 100 kλ.Given the high S/N of our combined datasets, the relative flux uncertainty remains relatively contained after disk subtraction (see Fig. 6).In the upcoming Section, we lay out the methods we followed to model the full visibilities and subtract the disk component to obtain an envelope-only spectral index.

Disentangling disk and envelope emission
To investigate dust properties in the envelope of L1527, we worked in the visibility (u, v) plane.This approach has the main advantage of keeping the analysis clear of the artifacts and nonlinearity from which image reconstruction algorithms suffer.
The visibility amplitudes are the combined flux of disk and envelope at all uv distances.Since we are mainly interested in studying the spectral index of the envelope, we have to disentangle the compact emission flux from the extended one.Hence, we need to fit the data with a model appropriate for a compact disk that, at the same time, self-consistently accounts for excess emission at the shortest baselines (largest physical scales).
Brightness profiles of protostellar envelopes of diverse class 0/I YSOs have been found to be well fit by broken power laws with exponents in the (-2.5,-1) range (Plummer 1911, Motte & André 2001, Shirley et al. 2002, Maury et al. 2019).Thus, we consider a geometrical model in which a Gaussian, that will trace the compact emission, is summed to a power law (Plummer profile): The Gaussian is defined by a peak f ′ 0 = rI 0 , a width σ, an inclination inc and position angle PA.The envelope power law has peak flux f ′′ 0 = (1 − r)I 0 and three free parameters: R i , R out , p.These are such that the envelope emission is constant within a radius R i , it follows Eq. 1 between R i and R out and it is null beyond R out .The peak fluxes of the extended and compact emission are modulated by the free parameter 0 < r < 1.We have performed the fitting with galario, a library that exploits the power of modern graphical processing units (GPUs) to accelerate the analysis of observations from radio interferometers (Tazzari et al. 2018).In our framework, galario would compute a model image given our total profile, then it would Fouriertransform it into synthetic visibilities and sample them at the uv-points covered by the antenna configurations with which the observations were performed.Then, it would run a minimum-χ 2 fit between the data and the model visibilities.The 10 fit parameters were the flux amplitude I 0 and the r factor; the width σ, the inclination (inc) and position angle (PA) of the Gaussian; the inner and outer envelope radii, R in and R out , and the power law exponent p. Finally two parameters, dRA and dDec, fit the offset of the peak from the phase center.For each fit, 80 walkers were set to run for 3,000 burn-in steps and 7,000 more steps after the burn-in stage.The best fit results are summarized in Table C.We show, as an example, the single best-fit model for the B3 88 GHz data, decomposed in its power law and Gaussian components (Fig. 4).We report the rest of the fits in  C) in Appendix C. In particular, we find that the compact emission has a deprojected width of 0 ′′ .15 < σ f it < 0 ′′ .28 on the sky, depending on the Band.If we define the radius of the disk to be the 2σ contour of the Gaussian component, then R disk ∼ 75 au in B7, consistent with what was found in the kinematical analysis of Aso et al. (2017).In addition, the Gaussian component contributes to a minimum of 75% in B4, up to 98% in B3 (see Appendix C for the detail).The derived disk fluxes are listed in Table 4.This fitting procedure allows us to disentangle disk and envelope emission and study the spectral energy distribution of the disk and envelope separately in the next section.

Spectral and dust emissivity indices
In this Section, we first derive the spectral indices of the observed emission at both disk and envelope scales, then we discuss to what extent these can be used to evaluate dust grain properties.4).The protoplanetary disk of L1527 is optically thick up to 3.4 mm (α B7−B3 ∼ 2), likely due to its edge-on nature.

The circumstellar disk
To study possible early dust grain growth in the envelope of L1527, we made use of long baseline data to properly model the compact emission to be subtracted from the short baseline amplitudes.In doing so, we obtained estimates of the flux of the compact (0.15" < σ f it < 0.27") emission of L1527 in four ALMA bands.We here show the spectral energy distribution of the disk we obtained by combining and modeling all suitable ALMA archival data for L1527 (Figure 5, Tab. 4).To offer a more comprehensive view to the reader, we extended the range of wavelengths by including literature Very Large Array (VLA) measurements from Melis et al. (2011), up to 6 cm.
Fitting the SED in the 0.88 to 7-millimeter range using our measurements along with the values of Melis et al. (2011), we find an optically thick disk with α = 2.1 ± 0.1.This is not surprising: the disk orbiting L1527 is very inclined (i ∼ 80 • ), almost edge-on to the sky: this geometric factor contributes to increase its optical depth along our line of sight.Furthermore, Ohashi et al. (2022) used high-resolution, multiple frequency observations of the disk and found that the brightness temperature toward the mid-plane of the disk obtained by converting the flux at 0.88 mm (42 K) is much lower than that found at 2 and 3 mm (60 K and 90 K, respectively).This progressive increase in brightness temperature with wavelength means that the disk is at least partially optically thick at 0.88 and 2 mm.
The high optical depth of the circumstellar disk of L1527 prevents us from concluding anything about its dust grain properties based on the spectral index alone.Finally, the 3.5 and 6 cm fluxes can be used to estimate the free-free emission that affects the measurements at the shorter studied wavelengths.The slope at these longer wavelengths is α 35−60 = 0.33 ± 0.17, as reported in Fig. 5 and is compatible with what is expected for free-free emission (Panagia & Felli 1975).However, it is worth noting that Melis et al. (2011) showed how the 3.5 cm flux of L1527 is highly variable in time.Indeed, they listed eight measurements (references therein) at this wavelength taken between 1996 and 2010 where the 3.5 cm flux shows significant variations between 0.55 and 0.81 mJy, with errors of 1-5%.In Fig. 5, we show Melis et al. (2011) points taken at roughly simultaneous epochs (2010 July 30th).The contribution of the free-free emission is lower than 10% in B3 and lower than 1% in B7.

The envelope
We investigate changes in the spectral index throughout the spatial scales of the envelope in the visibility plane.Starting from the common shortest baseline of ∼ 9 kλ, we log-uniformly bin the visibilities across the available uv-distances and then measure the spectral index in each bin.Since we are interested in studying the spectral index of the envelope, we have subtracted the galario best-fit Gaussian component of our model (see Section 3) from the data.This way, the remaining flux is the contribution of the extended emission only (Fig. 6).
Here, the amplitude and related error bars account for both statistical and calibration uncertainties.We set the latter to 10% for B7 and 5% for bands 3, 4 and 6, following the prescriptions of the ALMA Handbook and on-sky analysis (respectively; Remjian et al. 2019;Francis et al. 2020).While the error on α is dominated by flux calibration errors on the observed amplitudes (see Fig. 3), the relative amplitude error becomes important after the subtraction of the disk model.Not only is the relative error overall larger, it also visibly increases with uv-distance.This effect is clear in Fig. 6.Using almost every archival dataset for this source, we made an effort to maximize the S/N of the faint envelope of L1527 in what is its deepest continuum analysis.
To make the best use of our multiwavelength data, we fitted a line to the fluxes along all wavelengths in each uv-distance bin and we more robustly determined the spectral index.Since we combined a number of datasets for each band, the frequency of each flux point in each bin is taken to be the weighted average of the frequencies in that bin.The uncertainty on the spectral index obtained by fitting all bands is the error on the slope obtained with the weighted linear regression.Figure 6 shows that the fitted spectral index of the envelope of L1527 appears roughly flat in the studied range that spans from ∼ 2000 down to 50 au, α ∼ 3.For completeness, we show the spectral indices computed between adjacent bands in Appendix B. These are flat as well, although they show systematic differences in their average value, which clearly shows the need for a multifrequency (>2) approach.
Finally, we measured the spectral index in different angle bins on the sky, after removal of the compact component.We divided the uv-plane into three different regions: the envelope (−20 • < PA < 20 • ), the outflow (70 • < PA < 110 • ) and the cavity walls (the remaining zones).We do not observe any significant difference among the three spectral indices at scales larger than 200 au (Fig. 7).

Dust emissivity index
The simple link between the exponent of the dust opacity power law and the spectral index, β = α + 2, only holds for optically thin emission when the RJ approximation is valid.Thus, before interpreting the values of the envelope's spectral index in terms of dust properties, we check whether these necessary conditions are met at the wavelengths we probe.First, we evaluate the optical thickness of the emission.The specific intensity I λ , or flux per unit solid angle that we receive from the source (F λ /Ω) can be generally expressed as an absorbed black body: where the optical depth τ λ modulates the difference between the observed flux and the optically thick black body emission.If τ λ ≪ 1, then F λ /Ω ∼ τ λ B(T dust ), while if τ λ ≫ 1, then the observed emission tends to a black body spectrum.
To evaluate the optical thickness, we consider a typical dust temperature prescription for a centrally illuminated protostellar envelope (Shu 1977, Butner et al. 1990, ?, Motte & André 2001): Considering this relation and the range of scales that we probe, the temperature gets as low as 10 K in the outermost radii of the envelope (∼ 2000 au).We find that, even at the shortest of our wavelengths, the envelope emission is optically thin at all scales as F 0.88mm /Ω < B(T dust ) by a factor of 10 (Fig. 8).The optically thin regime is naturally satisfied at longer wavelengths.Secondly, we check if the RJ approximation is valid.Based on the same temperature profile for the dust, in the cold outskirts of the envelope (∼ 20 K at 500 au), hν/k B T ∼ 0.8 at our mean B7 frequency of 348 GHz.This violates the RJ condition hν/k B T ≪ 1.Thus, while the envelope of L1527 is in the optically thin emission regime, it does not satisfy the RJ approximation and β α − 2.
Thus, using the same profile, we can obtain a temperaturecorrected β profile as: where A is a constant (similarly to what done in G19, but for multiple frequencies).Figure 9 shows how the dust emissivity index is nearly α − 2 for the smaller spatial scales, where T(R) is large enough, while it starts to show diskrepancies after 200 au (where T∼ 25 K).We also investigated the effects on β of changing the temperature profile's radial power law exponent.We found that exponents in the interval [-0.3,-0.6]do not affect our conclusions.We exclude steeper profiles as they lead to unreasonably low values of the temperature at the studied scales.
The error on the dust emissivity index has been considered to be the same one that affects the spectral index since we are assuming an exact temperature profile, thus ∆α = ∆β.
In Figure 9, we report the measured β as a function of uvdistance, and we overplot a linear fit to the points.We have run the fit by means of the Bayesian method described in Kelly (2007) and implemented in linmix3 .The priors were uniformly distributed.We plot some single chain results in Fig. 9, along with the best fit model.We find evidence for an outward increasing β, rising from 0.8 to 1.6 between 50 and 2000 au radial scales:  (5) Fig. 7. Spectral index of L1527 as measured along directions of the outflow (violet) and envelope (green), cavity walls (orange).The total α is also plotted (gray).The way the uv-coverage has been sampled is shown in the upper inset (B7 as an example).The scarce uv-coverage of the ACA B7 observations cause some bins to be undefined at the short baselines, hence the gaps.In all cases, α is only computed starting from 20 kλ.
In turn, Eq. 5 yields: These measurements represent tentative evidence for dust optical properties variations throughout the envelope of L1527.Furthermore, low β at scales lower than ∼ 300 au are consistent with the presence of submillimeter sized grains in the inner envelope.

Literature comparison
G19 studied the envelope continuum emission of a few class 0/I sources, including L1527.Here, we briefly compare our findings to theirs.First, the wide uv-distance range of the datasets we used allowed us to model the disk emission spatial profile, while only a constant amplitude was subtracted in the amplitudes in G19, who assumed an unresolved compact disk.This spatial modeling is important to properly subtract the compact emission from the total amplitudes and work out the rest of the analysis on the flux of the envelope alone.G19 reported a slope of β of 0.18 ± 0.39.Noticeably, this flat β profile of G19 is consistent with what we observe across the same region they studied, 70-700 au radii.Extending the studied scale by a factor 2, we are able to explore the outer envelope and determine the positive outward gradient of Eq. 5.At the outer probed radii (∼ 2000 au), the derived values for β are consistent with typical ISM measurements.Secondly, G19 reported a fitted β 500au = 1.41 ± 0.16 while we find β 500au = 1.15 ± 0.08 based on Eq.4.These improvements have been obtained thanks to higher sensitivity data as well as the extension of the studied frequency ranges.

Discussion
We here discuss our results in further detail, how they compare with existing theoretical and observational literature, and the caveats relevant to our analysis.is consistently lower by at least an order of magnitude than a black body spectrum of emitting dust with a radial temperature profile T ∝ r −0.4 (green).The optical depth (violet) is much smaller than 1 at all envelope scales, thus we can consider the envelope emission at 0.88 mm as optically thin.

Dust growth in envelopes is challenging
The low β value we measured at scales of a few hundred au from L1527 is a tentative evidence for grown dust grains in the inner envelope of the class 0/I YSOs L1527.If solids indeed grew in the envelope, this would represent a shift in the usually assumed initial conditions for planet formation, which only include < 0.25µm grains at such scales.However, a challenge to our measurements rises when trying to link these low β with the in situ formation of relatively large dust grains in the envelope environment.Indeed, theoretical studies of grain agglomeration have not been able to reproduce such growth in envelope-like environments.The numerical simulations of Ormel et al. (2009), Lebreuilly et al. (2023) and Bate (2022) have shown that it is very unlikely that solid dust grains could grow larger than about 1 µm at the typical densities and timescales of envelopes such as the one infalling on L1527 and its young disk.These works included a thorough treatment of grain-grain collision energetics, accounting for relative grain velocities due to turbulence, Brownian motions, ambipolar diffusion and hydro-dynamical drift.Additionally, Silsbee et al. (2022) proposed a simpler, analytical model to place strict upper bounds on the maximum grain size that can be reached in extended protostellar envelopes.They considered a coagulation model for the growth of spherical, yet fractal grains whose relative motions are driven by turbulence and the fragmentation of which happens for velocities above a quite generous threshold of v rel = 10 m/s.Considering that the optical properties of dust grains depended on the product of their radius (a) with a filling factor (or porosity) ϕ, they find that it is not possible to grow grains with an optical radius a opt = aϕ of 1 millimeter, in typical protostellar envelope conditions.They find that, in the 10 5 yr lifetime of a typical class 0 source whose volume density is approximately 10 7 cm −3 , grains grow up to only about 2.5 µm.

Grain size, composition and temperature effects on β
While these theoretical works strongly disfavor the in situ growth of (sub)millimeter dust grains in the envelopes of young objects, the observational evidence of the low dust emissivity indices found therein (Kwon et al. 2009, Miotello et al. 2014, Galametz et al. 2019) remains largely uncontested.
First, an objection to the link between low β and large grains might arise from the dependency of the dust-emissivity index on other dust properties.In fact, the dust emissivity index is sensitive to changes in dust composition, porosity and the presence (or absence) of an icy mantle (e.g., Beckwith & Sargent 1991, Testi et al. 2014).Although variations in β are evident and wellcharacterized across different grain properties, β < 1 remains a strong evidence for the presence of (sub)millimeter grains in the studied distribution ( Natta & Testi 2004, Draine 2006, Natta et al. 2007, Testi et al. 2014and Köhler et al. 2015, Ysard et al. 2019).Thus, although further laboratory studies will be necessary to more precisely interpret β in terms of both grain sizes and compositions, dust particles significantly larger than what expected for the typical ISM grains (0.2-2 µm) remain a robust explanation for our observations.Secondly, if the compact emission contribution is not properly subtracted from the visibilities, a bias can be introduced when determining the slope of the SED.Previous studies (e.g., Miotello et al. 2014 and G19) have subtracted the compact emission assuming unresolved disks that would only contribute an offset at all uv-distances.They determined this value by fitting a constant to the visibilities amplitudes after some arbitrary long baseline and subtracted it from the visibilities at shorter uvdistances.Here, we have tried to more properly model and subtract the disk emission as explained in Section 3.
Finally, perhaps the temperature correction discussed in Section 4.3 could be revisited.If the emission is optically thin, Equation 4 reduces to: where the second term is two in case the RJ approximation is satisfied.The temperature power law we consider (cfr.Eq.3) yields temperatures of about 25 K at 200 au, where we find β ≲ 1.
Here, the correction term of Eq.7 amounts to 1.85.If the true temperature of the medium at 300 au were lower than what is predicted with our simple power law, then the real β would be higher.However, to produce β ∼ 1.7, the temperature of the envelope at 300 au would need to be lower than 5 K, which is far below a reasonable value for a protostellar envelope.

Alternatives: differential collapse or outflow transport
Although simulations predict so far that large grains cannot grow in situ, observations strongly hint to their presence in the inner few hundreds au of these protostellar envelopes.Different mechanisms that can justify their presence might be at play.First, Lebreuilly et al. (2020) investigated how the distribution of dust grains of a prestellar core evolves during the early phases of the protostellar collapse and how this evolution depends on the initial conditions of the cloud and the dust distribution.They found a significant decoupling between gas and dust for grains of a few 10 µm.Moreover, observations of scattered light from molecular cloud cores in the 3-5 µm range ("coreshine") have provided evidences for the presence of grains larger than usually expected ISM ones: up to 10 µm compared to about 0.25 µm (Steinacker et al. 2010, Steinacker et al. 2014, Steinacker et al. 2015).These larger grains would tend to settle more efficiently in the first-core, leading to a radial stratification of dust properties that could reproduce a gradient of β as we observe here, but hardly explain the low β values at the inner envelope scales.Further refinement of this kind of model will require deeper knowledge of prestellar core dust properties.
A second scenario to explain large grains in the envelopes of class 0/I YSO was studied by Wong et al. (2016), who proposed that such grains might be transported to the envelope after the growth has happened in the disk in the very early stages of the system.Their model suggested that a typical protostellar outflow (T∼ 10 K, v ∼ 10 km/s) could lift grains as large as 1 mm in the first 10 4 yr of the protostellar lifetime if the mass loss rate of the protostar was high enough (∼ 10 −6 Ṁ⊙ /yr).Among others, Brauer et al. (2008), Kawasaki et al. (2022), Bate (2022), Lebreuilly et al. (2023) all find that growth up to cm-sized pebbles is extremely fast (a few ∼ 10 3 yr) at the high density of the disk environment, thus producing the large grains that the outflow would entrain in the first place.Based on three-dimensional magneto-hydro-dynamical simulations, Lebreuilly et al. (2020) and Tsukamoto et al. (2021) found that the large dust grains grown in the inner region of a disk can be entrained by an outflow up to the envelope scales (up to 100 µm according to the first work, and even >1 mm for the second).These grains then decouple from the gas and are ejected from the outflow into the envelope itself, enriching its dust population before falling back to the disk.Chen et al. (2016) found evidence for low β correlated with protostellar outflow locations in Perseus.More recently, G19 observed a correlation between β and the envelope mass of the CA-LYPSO sources, which might be in turn caused by a correlation between envelope mass and CO outflows fluxes (Bontemps et al. 1996).This two relationships together might support the scenario of disk dust transport into class 0/I protostellar envelopes.
While a more thorough treatment of this link for L1527 is outside of the scope of this work, we note that L1527 indeed hosts a large scale outflow structure in the east-west direction, perpendicular to the edge-on disk.Hogerheijde et al. (1998) observed the source with the James Clerk Maxwell Telescope and detected in 12 CO (J=3-2) line emission an outflow extending over about 20000 au.Moreover, Tobin et al. (2010) imaged this source with the Gemini North telescope (in its L' band) and with the Infrared Array Camera (IRAC) on the Spitzer Space Telescope (2.15-8.0µm), detecting the outflow cavity out to roughly 10000 au in both cases.Future investigations should consider whether the L1527 outflow can transport (sub)millimeter grains from the inner disk and lift them in the envelope.Such research is vital to our understanding the origin and consequences of relatively large grains at envelope scales, and may even shed light on this evolutionary stage of disk-based planet formation models.

Caveats
Although we aim to present a more robust way of measuring the spectral index of the faint extended emission from the envelope of class 0/I YSOs, our analysis is not free of caveats.Additional work remains to further improve the reliability of our claims.
First, while we present a study of the spectral-and dust emissivity indices as a function of uv-distance and physical source scales, we have to keep in mind that a particular baseline does not only probe an annular region at some distance from the center, rather it probes any physical structure within a scale given by θ ∼ λ/B.It is for this reason that we have to model and subtract the disk from the visibility amplitudes before evaluating the spectral index for the envelope and, for the same reason, the α and β we compute for the envelope at different scales represent a flux-weighted average over the relevant envelope spatial scales.If the inner envelope is much brighter than the outer one, as expected, the flux at even the shorter baselines will be dominated by the inner envelope.Thus, the spectral index value at these scales would still be mainly describing the inner region.
Second, in Section 4 we have calculated the dust emissivity index β of the envelope taking into account diskrepancies from the RJ approximation, hence accounting for a temperature profile, T ∝ R −0.4 .This power-law was derived for a dusty envelope illuminated by a central protostar (e.g., ?).While the values we obtain using this proportionality are reasonable, as they do not reach values lower than typical cores temperatures at typical cores scales (10 K at 2000 au;e.g., Ferrière 2001), a more thorough treatment that accounts for the temperature structure and its associated uncertainties can be achieved by post-processing 3D envelope models with radiative transfer in order to.

Conclusions
We aim to characterize the maximum grain size of the dust distribution in the envelope of a class 0/I YSO: L1527 IRS, or L1527 for short.Given its vicinity (140 pc), this source has been extensively studied over decades with different types of data and methodologies to investigate its star, extended envelope, and circumstellar disk both in the dust continuum and line emission.In this context, we exploit the richness of data from the Atacama Large (sub)millimeter Array to greatly enhance the S/N of the extended emission in the continuum and determine the peak of the dust grain population in the envelope that surrounds and infalls onto the young circumstellar disk of L1527.
We find that: • The spectral index of the compact disk, i.e. within a radius of ∼ 75 au, is α 0.88−3.4mm= 2, consistent with the expected high optical depth of its edge-on geometry.No dust grain properties information can be extracted in such a condition without thorough radiative transfer modelling, outside of the scope of this work.• After subtraction of the inner disk emission and correction for the RJ approximation, we find tentative evidence for a positive outward gradient of the dust emissivity index of the envelope of L1527 (β(R) ∼ 0.37 log(R au ), which can be interpreted as variations of dust properties from ISM-like particles (β ∼ 1.6) at 2000 au down to grown grains below 300 au (β ≤1).• The implications of physical and chemical properties of dust grains on the values of β as well as the possible impact of the temperature profile used to calculate the correction for departure from the RJ approximation in computing the dust emissivity index (cfr Section 5.1) is carefully discussed.None of these possibilities alone are able to justify β values as low as 1 at hundreds of au, thus strengthening the hypothesis of submillimeter grains at these scales.While in situ formation seems disfavored by most theoretical studies, differential core dust collapse or outflow dust transport from the disk could explain our observations.
• We argue that a multiscale (10 1 − 10 3 au), multifrequency (n>2), and high-sensitivity study is necessary to tightly constrain the spectral-and dust emissivity index profiles of faint class 0/I YSO envelopes.Using ALMA B3, B4, B6 and B7 effectively halves the error bars of measurements made in previous studies that made use of two wavelengths only.
ALMA has provided us with the necessary capabilities in terms of resolution, recoverable scales, sensitivity, and frequency ranges that are critical to the study of the continuum emission of faint, extended envelopes such as the one infalling onto the young protostar L1527.It is now the time to exploit this infrastructure to conduct sample studies of these objects with state-ofthe-art data sets to finally pinpoint the initial conditions of grain growth and planet formation in both space and time.

Fig. 1 .
Fig. 1.Zoom in to the inner 6 ′′ .0 of the ALMA B3 continuum of L1527 obtained with the 12m array setups of the FAUST Large Program.The color map shows only flux densities higher than ∼5σ, with the white contours highlighting the [5, 15, 30]σ levels.The ALMA synthesized beam is shown in black in the lower left corner.

Fig. 2 .
Fig.2.ALMA B6 continuum of map of L1527 obtained imaging the 7m array setup of the FAUST Large Program.The inner 4 ′′ .0 inset has been imaged using ALMA B6 FAUST 12-m array setup data.The color map for the ACA image shows only flux densities higher than 5σ, and the white[5,10, 50, 150]σ level contours.The color map of the 12marray image shows only flux densities higher than 20σ, and the white[20,60,100]σ level contours.The synthesized beams are shown in black in the lower left corner in both cases.We note that the colorbar is different for each map.

Fig. 4 .
Fig. 4. Final Plummer+Gaussian best fit (orange) is overplotted on the real part of the visibilities for the B3 (88 GHz) data (upper panel, black points).The fit on the imaginary is in the lower panel.The Plummer only (violet line) and Gaussian only (green line) components of the total model are also shown.In the case of the 88 GHz emission, the disk contributes up to 98% of the total emission.The wiggles in the model are due to its sampling on the uv points of the observations.

Fig. 5 .
Fig. 5. Spectral energy distribution of the circumstellar disk around L1527.Our galario best-fit total flux of the compact Gaussian component for each band, in green, along with longer wavelength VLA measurements from Melis et al. (2011) (see Table4).The protoplanetary disk of L1527 is optically thick up to 3.4 mm (α B7−B3 ∼ 2), likely due to its edge-on nature.

Fig. 6 .
Fig. 6.L1527 ALMA B3 (orange), 4 (green), 6 (blue), 7 (red) visibilities after removal of the fitted compact Gaussian component.The spectral index of the envelope emission (gray) has been obtained by fitting a line to the fluxes at all wavelengths in each uv-distance bin.

Fig. 8 .
Fig.8.Flux density of the envelope emission in ALMA B7 (orange) is consistently lower by at least an order of magnitude than a black body spectrum of emitting dust with a radial temperature profile T ∝ r −0.4 (green).The optical depth (violet) is much smaller than 1 at all envelope scales, thus we can consider the envelope emission at 0.88 mm as optically thin.

Fig. 9 .
Fig. 9. Dust emissivity index (purple dots) of the envelope of L1527 as a function of radial scales.The solid dots indicate the β computed fitting all available bands.The dashed orange line is the best model for the profile of β along scales, while the light orange lines are a subsample of the results of some chains of the fitting procedure.The black line shows the approximation in which β = α f it − 2, where α is the slope of the SED considering all available bands.

Fig. C. 1 .
Fig. C.1.The galario fit of the B3 (100 GHz) real and imaginary part of the visibilities.The best model (orange) is composed of a compact Gaussian model (green) and an outer power law (violet).The wiggles in the plotted model are due to its sampling on the uv points of the observations.