Formation of dust clumps with sub-Jupiter mass and cold shadowed region in gravitationally unstable disk around Class 0/I protostar in L1527 IRS

We have investigated the protostellar disk around a Class 0/I protostar, L1527 IRS, using multi-wavelength observations of the dust continuum emission at $\lambda=0.87$, 2.1, 3.3, and 6.8 mm obtained by the Atacama Large Millimeter/submillimeter Array (ALMA) and the Jansky Very Large Array (VLA). Our observations achieved a spatial resolution of $3-13$ au and revealed an edge-on disk structure with a size of $\sim80-100$ au. The emission at 0.87 and 2.1 mm is found to be optically thick within a projected disk radius of $ r_{\rm proj}\lesssim50$ au. The emission at 3.3 and 6.8 mm shows that the power-law index of the dust opacity ($\beta$) is $\beta\sim1.7$ around $ r_{\rm proj}\sim 50$ au, suggesting that grain growth has not yet begun. The dust temperature ($T_{\rm dust}$) shows a steep decrease with $T_{\rm dust}\propto r_{\rm proj}^{-2}$ outside of the VLA clumps previously identified at $r_{\rm proj}\sim20$ au. Furthermore, the disk is gravitationally unstable at $r_{\rm proj}\sim20$ au, as indicated by a Toomre {\it Q} parameter value of $Q\lesssim1.0$. These results suggest that the VLA clumps are formed via gravitational instability, which creates a shadow on the outside of the substructure, resulting in the sudden drop in temperature. The derived dust masses for the VLA clumps are $\gtrsim0.1$ $M_{\rm J}$. Thus, we suggest that Class 0/I disks can be massive enough to be gravitationally unstable, which might be the origin of gas-giant planets in a 20 au radius. Furthermore, the protostellar disks can be cold due to shadowing.


INTRODUCTION
The formation of protostellar disks is an important process not only for the physical processes of accretion to a protostar and removal of angular momentum, but also for the initial conditions of planet formation. Recent high-resolution observations of molecular lines have revealed the transition of gas motion from accretion to Keplerian rotation toward Class 0/I young stellar objects (YSOs), which allows us to investigate the disk formation process (Tobin et al. 2012;Murillo et al. 2013;Ohashi et al. 2014;Sakai et al. 2014a;Yen et al. 2014;Aso et al. 2015). In addition, recent observations using dust continuum emission have revealed the beginning of the formation of structures such as warped, spiral, and ring structures in Class 0/I protostellar disks (Sheehan & Eisner 2017Sakai et al. 2019;Sheehan et al. 2020;Segura-Cox et al. 2020;Nakatani et al. 2020).

satoshi.ohashi@riken.jp
It has been suggested that ring structures are related to planet formation in the early stage. It is thus important to clarify the formation mechanism. Various mechanisms for ring formation have been studied, including those that involve unseen planets (Goldreich & Tremaine 1980;Nelson et al. 2000;Zhu et al. 2012), snowlines of volatile gas species freezing-out onto dust grains (Zhang et al. 2015;Okuzumi et al. 2016), magneto-rotational instability (Flock et al. 2015), disk winds in unevenly ionized disks (Takahashi & Muto 2018), a growth front (Ohashi et al. 2021), secular gravitational instability (Takahashi & Inutsuka 2014), and coagulation instability (Tominaga et al. 2021). In addition to the ring formation mechanisms described above, gravitational instability has been theoretically suggested to be important in Class 0/I protostellar disks (e.g., Nakamoto & Nakagawa 1994;Vorobyov & Basu 2006Machida et al. 2011;Tsukamoto et al. 2017), and is a possible mechanism for gas-giant planet formation (e.g., Gammie 2001;Kratter & Lodato 2016;Vigan et al. 2017). The fragmentation of a gravitationally unstable disk and the associated formation of gas-giant planets could be the origin of the existence of planets in Class II protoplanetary disks, as implied by several observations (Zhang et al. 2018;Teague et al. 2018;Pinte et al. 2018). Therefore, understanding the grain growth and gravitational instability in protostellar disks will help us determine when and how planet formation begins. Note that it has been suggested that the disk masses of Class II protoplanetary disks are insufficient for planet formation (e.g., Manara et al. 2018), and that planet formation in Class 0/I disks may be important from the perspective of the mass budget problem.
Dust temperature is one of the most important parameters for understanding how and where planets begin to form because it is related to the chemical composition of planetary atmospheres. Dust temperature has been investigated in protoplanetary disks using various molecular lines (e.g., Kamp & Dullemond 2004;Notsu et al. 2020;Öberg et al. 2021;Liu et al. 2021). It has recently been suggested that the temperature in Class 0/I disks is hotter than that in Class II protoplanetary disks (van 't Hoff et al. 2018;Zamponi et al. 2021;Liu 2021).
To explore the physical conditions for the earliest embedded protostellar system, we study a representative Class 0 low-mass protostar, namely IRAS 04368+2557 in L1527 IRS. L1527 IRS is located in the Taurus molecular cloud at a distance of d = 137 pc (Torres et al. 2007), one of the closest star-forming regions. It is a young protostar with a bolometric luminosity of L bol = 2.75 L (Tobin et al. 2008) and a bolometric temperature of T bol = 44 K (Kristensen et al. 2012). In the recent study, they are updated to L bol = 1.6 L and T bol = 79 K, respectively by using the Herschel-PACS far-IR spectroscopy (Karska et al. 2018).
The existence of a large Keplerian disk with a size scale of 100 au was reported by Tobin et al. (2012) based on a kinematic analysis of the 13 CO (J = 2 − 1) line. Ohashi et al. (2014) identified a Keplerian disk associated with an infalling envelope based on ALMA observations. A dramatic change in chemical composition across the centrifugal barrier has also been discovered, which suggests that the disk is still growing (Sakai et al. 2014a,b). Aso et al. (2017) investigated the rotation and infalling kinematics around L1527 IRS with much higher sensitivity as well as higher angular resolution and estimated the disk radius and stellar mass to be ∼ 74 au and ∼ 0.45 M , respectively.
Recent high-resolution observations have revealed complicated structures of the disk. Sakai et al. (2019) identified different orbital planes between the inner and outer parts of the disk with a boundary at 40 − 60 au, suggesting a warped structure. In addition, Nakatani et al. (2020) identified several clumps inside the 20 au radius of the disk using VLA archival data. Although the origin of these substructures is still under debate, the VLA clumps may be related to planet formation in the protostellar disk.
Regarding the envelope scale ( 100 au), Ohashi et al. (1997) identified an edge-on flattened envelope associated with infalling motion, which was confirmed by Yen et al. (2013). Radio continuum jets have been identified by the VLA observations (Loinard et al. 2002;Reipurth et al. 2004). A bipolar outflow was identified in 12 CO molecular lines with single-dish telescopes (Tamura et al. 1996;Hogerheijde et al. 1998). Mid-infrared observations with the Spitzer Space Telescope and L band imaging with the Gemini North telescope show bright bipolar scattered light nebulae along the outflow axis on the 10 3 − 10 4 au scale (Tobin et al. 2008). The protostellar disk is highly embedded in the infalling-rotating envelope, proving that the disk around L1527 IRS is still growing.
In this paper, we analyze the multi-wavelength dust continuum emission from sub-millimeter to millimeter wavelengths toward the protostellar disk around L1527 to investigate the origin of the substructures and early planet formation in the disk-forming stage. We achieved a resolution of 0. 16 (corresponding to 22 au) at all observation wavelengths, allowing us to study the detailed structure of the disk. In particular, the resolution in the disk height direction is 5 au, which is sufficient to discuss the vertical structure of the young growing disk. The rest of this paper is organized as follows. The data reduction is described in Section 2, and the resulting continuum images are presented in Section 3. In Section 4, the physical conditions for the disk, such as inclination (Sec 4.1), low temperature due to the shadowing effect (Sec 4.2), dust scale height (Sec 4.3), grain growth (Sec 4.4), and gravitational instability required for the origin of the substructures, leading to planet formation (Sec 4.5), are discussed. A summary of the results and a discussion are given in Section 5.

OBSERVATIONAL DATA
We use the dust continuum data for L1527 taken with ALMA Band 7, Band 4, Band 3, and JVLA Q-Band over the period from 2011 to 2021. The detailed observations and data reduction for ALMA Band 7, Band 3, and JVLA Q-Band are described in Nakatani et al. (2020). In this paper, we present new observations of ALMA Band 4 and Band 7 with a higher spatial resolution. The details are summarized in Table 1.

ALMA Band 4 Observations
The L1527 protostellar disk was observed with ALMA Band 4 on September 15 and 21, 2021, with a configuration of C43-9/10. Four spectral windows with a bandwidth of 1.85 GHz each were used for the dust continuum and molecular lines (e.g., H 2 CO). We used all of the spectral windows for the dust continuum emission as a representative frequency for 146 GHz because no molecular lines were detected. The total bandwidth of the continuum map was 7.4 GHz. The observations were composed of three tracks, each with an integration time of ∼ 40 min on the source. J0237+2848 and J0510+1800 were observed as the passband and flux calibrator. J0433+2905 was observed as the phase calibrator. The antenna configuration covered baseline lengths of 178 to 15238 m. The maximum recoverable scale (θ MRS ) was estimated as ∼ 0.6λ/L min , where λ is the observation wavelength and L min is the minimum baseline. Because the maximum recoverable scale is ∼ 1. 5 (corresponding to ∼ 200 au), any structures that extend beyond this size will be resolved out. The data reduction was performed using a pipeline script provided by ALMA.

ALMA Band 7 Observations
The L1527 protostellar disk was also observed with ALMA Band 7 on September 26, 2021, with a configuration of C43-9/10. Four spectral windows with a bandwidth of 1.875 GHz each were used for the dust continuum and molecular lines. We used all of the line free channels for the dust continuum emission as a representative frequency for 345 GHz. The total bandwidth of the continuum map was 7.5 GHz. The observations were composed of one track with an integration time of ∼ 50 min on the source. J0510+1800 was observed as the passband and flux calibrator. J0438+3004 was observed as the phase calibrator. The antenna configuration covered baseline lengths of 70 to 14361 m. The maximum recoverable scale (θ MRS ) was ∼ 1. 6 (corresponding to ∼ 220 au). The data reduction was performed using a pipeline script provided by ALMA.
We performed self-calibration for the continuum observations in CASA. The entire imaging process was carried out with the CASA tclean task. For constructing CLEAN maps, we adopted Briggs' weighting with a robust parameter of 0.5 and −2. Images with higher sensitivity and a larger beam size were obtained for a robust parameter of 0.5 and images with lower sensitivity and a smaller beam size were obtained for a robust parameter of −2. The beam sizes and sensitivities are summarized in Table 1. 3. RESULTS Figure 1 shows dust continuum images toward the L1527 disk obtained at sub-millimeter and millimeter wavelengths of 0.87 to 6.8 mm. The disk structure is well resolved and similar to those reported in previous studies (e.g., Ohashi et al. 2014;Aso et al. 2017;Nakatani et al. 2020). We confirmed that the disk is close to being edge-on and is elongated in the northsouth direction. Because the disk has an almost edgeon geometry aligned from north to south, we assume that the ∆δ and ∆α directions correspond to the projected distances in the disk horizontal (∆x) and vertical (∆y) directions, respectively. The projected distance in the disk radial direction (r proj ) was measured based on the peak position in the Gaussian fit at each declination (r 2 proj = ∆x 2 + ∆y 2 ). The radial direction is mostly consistent with the ∆x direction.

Observed Disk Images
The dust emission in Band 7 extends to r proj 100 au, and that in Bands 4 and 3 extends to r proj ∼ 80 au. In contrast, the Q-Band emission is detected only within r proj ∼ 30 au. The smaller disk size in the VLA observations is attributed to the sensitivity limit due to weaker dust emission at longer wavelengths. The extended emission in Band 7 could be from an envelope component as well as disk emission.
The derived peak brightness temperatures (T B ) for Band 7, 4, and 3, and Q-Band emission are ∼ 40, ∼ 60, ∼ 70, and ∼ 50 K, respectively. The highest temperature (∼ 70 K) is found for Band 3 emission, which suggests that Band 7 and 4 emission is obscured in the edge-on geometry. In contrast, the Q-Band emission is affected by beam dilution (see Section 3.2). The detailed disk structure is discussed in Section 4.1. The total flux, beam size, and root-mean-squared (rms) noise level are shown in Table 1. The total flux for Band 7 emission (430 ± 11 mJy) is slightly lower than that in previous dust polarization observations (488±14 mJy) derived by Harris et al. (2018). Because the previous dust polarization observations were conducted with a compact configuration (C40-5), more of the envelope emission may have been recovered compared to that in our observations.
The high spatial resolution and sensitivity of the observations allow the vertical distributions to be determined. In particular, in the images for Bands 4 and 3, the disk appears to be wider toward the outer edge. This is consistent with the models developed by Tobin et al. (2013) and Aso et al. (2017), who suggested that the disk is flared with H ∝ r 1.2−1.3 , where H and r are the dust-scale height and radius of the disk, respectively. Note that the vertical structure with broadening around 40 − 60 au in radius was also reported by Sakai et al. (2017).
Here, we note that the Band 7 emission is detected over a wide spatial range along the disk vertical direction, extending to ±50 au. In contrast, the Band 4 and 3 emission is detected within ±20 au. These different spatial distributions suggest that the Band 7 emission traces not only the disk component but also envelope emission because the Band 7 emission would be sensitive to trace dust thermal emission, and the protostellar disk grows via infalling materials from the envelope. Figure 2 shows enlarged views of the disk images with higher spatial resolutions than those in Figure 1. The Band 7 image was produced using high-resolution data. The Band 4 and 3 and Q-Band images were produced by applying weighting with a robust parameter of −2 in the CASA tclean task. We found that about 10% of the total flux was resolved when the weighting was adjusted. Disk substructures within r proj ≤ 40 au were identified.

Disk Substructures in Inner 40-au Radius
The most prominent substructures were identified in the VLA Q-Band observations. Three clumps (clump-N, -C, and -S) were recognized around r proj ∼ 10 − 20 au along the disk elongation, which is consistent with Nakatani et al. (2020), who discussed these substructures as ring or spiral structures in the disk. We discuss the origin of these substructures in Section 4.5.
The VLA clump locations (cross symbols in the ALMA images) seem to be slightly offset from the disk midplane in the Band 3 image, which may be caused by the correction of the proper motion or absolute position error. We found that the ALMA images have no counterpart to the VLA clumps even though the spatial resolutions are comparable (the ALMA images have slightly larger beam sizes). This may have been caused by the high optical depth of the dust emission at the ALMA wavelengths. The maximum temperature in the VLA observations is T B ∼ 100 K; it decreases to T B ∼ 90 K in Band 3, 60 K in Band 4, and 50 K in Band 7. The lower maximum brightness temperatures at the disk center for the shorter wavelengths suggest that the shorter-wavelength emission is affected by the higher optical thickness in regions with outer part of the disk, where the temperature is lower.
In addition to the substructures in the radial direction, we resolved the vertical structure of the disk in all bands. The Band 7 image shows an asymmetric structure in the brightness temperature even though the VLA clumps are located at almost the disk midplane. The eastern side seems to be slightly warmer than the western side. We discuss this temperature asymmetry in terms of a flared disk model in Section 4.1. The Band 4 and 3 images show local enhancements in the vertical direction at the southern VLA clump. To confirm these substructures, Figure 3 plots the vertical distribution (∆y) of the normalized intensity of the Band 4 and 3 images at ∆x = −22 au, where the maximum local enhancements were found, and at ∆x = 22 au for comparison. The Band 4 and 3 intensity distributions both show skewed profiles at ∆x = −22 au. There is an additional component at ∆x ∼ −5 au shifted from the midplane. The location of ∆x = −22 au coincide with the southern VLA clump, which suggests that this local enhancement of the disk height is related to the VLA clumps, as discussed in Section 4.5.

Spectral Index Analysis
The spectral index, α, was calculated as where I and ν are respectively the intensity and frequency at each band. The intensity (I ν ) of the dust continuum emission was calculated using the radiative transfer equation.
where B ν is the Planck function, T bg = 2.7 K is the temperature of the cosmic background radiation, and τ is the optical depth of the dust continuum emission. We mapped the spectral index α by applying multifrequency synthesis (mtmfs) in the CASA tclean task with nterms = 2. The three α maps were produced using two band combinations, namely Band 7 − Band 4, Band 4 − Band 3, and Band 3 − Q-Band. Note that the spectral index for the Band 7 − Band 4 combination may have a first-order Taylor expansion uncertainty because nterms = 2 may produce a relatively lower spectral index if band separation is wide, as pointed out by Tsukagoshi et al. (2022). Therefore, the spectral index for the Band 4 − Band 3 combination is more reliable. We note that multi-frequency synthesis for the combination of ALMA Band 3 and VLA Q-Band is not available in CASA. This may be due to the different data setups for the ALMA and VLA observations. The α 7mm−3mm map obtained using Band 3 and Q-Band was constructed from the image plane produced with a natural weighting to recover the extended emission. The total flux of Q-Band increased to 4.8 mJy, indicating that the emission was significantly recovered by natural weighting (see also Table 1). Figure 4 shows spectral index maps for the three combinations. The α 2mm−0.9mm value is 2.0 in the inner disk (r proj 60 au), which indicates that the emission in Bands 7 and 4 is highly optically thick. Even in the outer part of the disk (r proj 60 au), the α 2mm−0.9mm value is lower than < 2.5, suggesting that the Band 7 emission is highly optically thick and the Band 4 emission may be moderately optically thick. Although larger dust grains also lead to a lower spectral index, the α 7mm−3mm value being larger than α 2mm−0.9mm and α 3mm−2mm in the entire disk indicates that the low spectral index is due to the optical depth effect rather than dust grain size. In particular, α 2mm−0.9mm and α 3mm−2mm are as low as ∼ 1.5 in the central region (r proj 40 au). A spectral index lower than 2.0 suggests that emission at shorter wavelengths traces a colder region. The optically thick emission allow only the foreground annulus of the edgeon disk to be observed. In contrast, emission at longer wavelengths can penetrate the inner radius of the disk because of its lower optical depth, which results in the spectral index being lower than 2.0. This is confirmed by the different brightness temperatures obtained at different wavelengths in Figure 1 and 2. This effect was previously reported by Galván-Madrid et al. (2018).
The spectral index α 7mm−3mm in the right panel of Figure 4 is ∼ 2.0 toward the central region, which may have been caused by the high optical depth. However, α 7mm−3mm becomes as large as 3.6 in the outer radius of r proj ∼ 50 au. The α 7mm−3mm value of 3.6 is consistent with that of the interstellar medium (ISMs) (e.g., Draine 2006), which suggests a small dust grain size in the outer part of the disk. The grain size distribution is discussed in Section 4.4.

Radial Dependence
To investigate the variation of the spectral index in the disk, we plot the spectral indices as a function of the horizontal distance (∆x) on the midplane in Figure  5. We confirmed that α 2mm−0.9mm and α 3mm−2mm are 2.0 in the disk for ∆x 60 au, indicating that the emission in Bands 7 and 4 is optically thick. In contrast, α 7mm−3mm becomes as high as 3.5 at ∆x ∼ 60 au. The lower limit of α 7mm−3mm > 3.4 was derived for ∆x ∼ 70 au from the non-detection of the VLA Q-Band. These large spectral index values suggest that the emission in Band 3 and Q-Band becomes optically thin in the outer part of the disk (∆x 60 au). Furthermore, α 7mm−3mm ∼ 3.6 is consistent with that for the ISMs, suggesting that grain growth has not yet begun, at least in the outer disk. We discuss the grain size distribution in more detail in Section 4.4. Figure 6 shows the vertical distributions (∆y) of the spectral index α 3mm−2mm at distances of ∆x = 0 au and ±60 au. We found that α 3mm−2mm increases with increasing ∆y in a radius of ∆x = 0 au, which indicates that the emission becomes optically thinner toward the disk surface. Note that we ignored the inclination effect because the disk has an almost edge-on geometry with an inclination of ∼ 5 • .  Profiles of spectral indices (α2mm−0.9mm, α3mm−2mm, and α7mm−3mm) versus horizontal distance in disk (∆x) at disk midplane. The error bars represent ±1σ, which was calculated using the 5% flux error and the rms values of the noise levels for the two-band data.

Vertical Dependence
The peak value of α 3mm−2mm is ∼ 3.4, found at ∆y ∼ ±20 au above the midplane, suggesting that the emission becomes optically thinner and dust grains are not significantly large ( mm). In contrast, α 3mm−2mm shows a constant value of ∼ 2.0 in the vertical direction of ∆y ± 30 au at a radius of ∆x = ±60 au. This suggests that the emission remains optically thick even at ∆y ∼ ±20 au above the midplane. The outer part of the disk remaining optically thick in the upper layers is consistent with a flared disk structure. Note that α 3mm−2mm decreases again in ∆y 20 au at a radius of ∆x = 0 au. This might be an artifact due to the Figure 6. Profiles of spectral index (α3mm−2mm) versus vertical distance in disk (∆y) at ∆x = 0 and ±60 au. The shaded regions represent ±1σ. The error bars were calculated using the 5% flux error and the rms values of the noise levels for the two-band data. The data are plotted only where the emission is more than 4σ.
side-lobe effect. The point spread function for our observations show that the side-lobe levels are ∼ 0.2. Because the peak emission was detected with a signal-tonoise ratio of > 500, the side lobe may produce artificial emission in the outer region.
The lowest value of α 3mm−2mm was ∼ 0.8, found at ∆y ∼ −10 au at a radius of ∆x = 0 au. If the disk is totally edge-on (inclination angle of 0 • ), α is expected to be a minimum at the midplane (Galván-Madrid et al. 2018). The offset of the lowest potion from the midplane was caused by a slight inclination of the disk. We discuss disk inclination in the following section (Section 4.1).

DISCUSSION
We have shown that the dust continuum emission from the protostellar disk of L1527 was detected in a wide range of observation wavelengths (0.87, 2.1, 3.3, and 6.8 mm). Based on these results, we discuss the disk geometry and physical conditions. Then, we discuss the origin of the substructures and the early planet formation scenario for a Class 0/I protostar.

Disk Inclination
The protostellar disk of L1527 is nearly edge-on. The inclination of the disk has been estimated to be ∼ 5 • based on a spectral energy distribution (SED) fitting (Tobin et al. 2013) and a kinematic analysis of the infalling-rotating motions of the envelope (Oya et al. 2015). In addition, Oya et al. (2015) suggested that the western side of the envelope faces the observer, indicating that the western side is the far side and the eastern side is the near side.
This disk geometry can be assessed using the asymmetric structure of the temperature along the vertical direction of the disk. If a flared disk is inclined, the near side of the disk will have a temperature gradient where the outer cold region is in front of the inner hotter region. The far side of the disk will have the opposite temperature gradient. For an optically thick emission, the near side should be fainter than the far side because it is obscured (Figure 7). Figure 8 shows the profiles of the brightness temperatures in Bands 4 and 3 along the disk vertical direction (east-west direction) across the protostar location of ∆x = 0 au. The enlarged profile of the Band 4 emission shows a central dip of ∼ 2 K (corresponding to ∼ 4σ). The peak temperature of the western side is slightly higher (∼ 1 K) than that of the eastern side. In addition, the Band 3 profile shows that the peak position is slightly shifted to the western side. These different temperature profiles are caused by the different optical depths of these bands. Due to higher optical depth toward the central region in Band 4, the emission is efficiently obscured. Furthermore, the asymmetric profile of the temperature in Band 3 indicates that the Band 3 emission is also optically thick (but less thick than that for the Band 4 emission). Similar asymmetric profiles have been found for other Class 0 protostellar disks, such as IRAS 16253−2429 (Hsieh et al. 2019), HH 212 , and OMC3-MMS6 (Liu 2021). The temperature of the western side being higher than that of the eastern side indicates that the western side is the far side and the eastern side is the near side of the disk. The disk geometry suggested by this profile analysis is consistent with that suggested by Oya et al. (2015).  The spectral index map could also be affected by the inclination. As shown in Figures 5 and 6, the lowest position of α 3mm−2mm is shifted by ∼ 10 au to the western side from the midplane. The Band 3 emission could be highly obscured toward the eastern side (near side) due to the edge-on geometry, and may be emitted from the vicinity of the protostar toward the western side (far side) (Figure 7). In this configuration, the western side shows a steep temperature gradient. Therefore, the lowest value of the spectral index is much lower than 2.0 on the western side. This scenario should be further assessed using radiative transfer calculations.
We plot the temperature profile of the Band 7 emission in Figure 9 in the same way as that in Figure 8. Interestingly, the temperature profile for Band 7 is opposite to that for Band 3. The eastern side of the disk is hotter than the wester side for Band 7, which is confirmed by the Band 7 image shown in Figure 2. If we apply the above discussion to the temperature profile for Band 7, the disk inclination needs to be opposite. Different inclinations might be the case for this disk because the disk has a warped structure with a boundary at 40 − 60 au, as suggested by Sakai et al. (2019). The warp angle is ∼ 5 • . If the disk is warped, it is also possible that the inclinations of the inner and outer parts of the disk are different . Because the Band 7 emission is highly optically thick, only the outer part of the disk can be observed, which leads to the different temperature profiles between the Band 7 and Band 3 emission. However, the Band 7 emission is highly optically thick and is contributed by the envelope components. The above discussion may not be applicable to the Band 7 profile. Further studies of gas kinematics are needed to assess this scenario.

Disk Temperature
The multi-wavelength observations of the dust continuum emission allow us to estimate the disk temperature at the midplane. In particular, the sum of the brightness temperature for an optically thick emission and the cosmic microwave background (T = 2.7 K) reflects disk temperature (see Equation (2)). Figure 10 shows the brightness temperatures for the ALMA Band 3, 4, and 7 emission shown in Figure 1 as a function of the projected radius (r proj ). The brightness temperature was calcu-lated from the intensity F d (in units of mJy beam −1 ) without the Rayleigh-Jeans approximation.
The brightness temperatures for Bands 3 and 4 are almost the same value for r proj 50 au, indicating optically thick emission. This is confirmed by the low spectral index (α 3mm−2mm ) in Figure 5. We also found that the brightness temperature for Band 7 is systematically higher (by a few degrees Kelvin) than that for Bands 3 and 4 in the outer region of the disk (r proj 30 au). Taking into account that the Band 7 emission extends beyond than the disk component, as discussed in Section 3.1, these results suggest that the Band 7 temperature is contaminated by those of other components. Envelope emission is a possible contributor because the disk is highly embedded in the parent envelope. If the Band 7 emission consists of the disk and envelope components, the Band 7 temperature is higher than the disk temperature by an amount due to envelope emission. It may also be possible that the Band 7 emission is contaminated by the emission around the centrifugal barrier, where the temperature is higher due to accretion shock, because the Band 7 emission traces the outer part of the disk. Unlike in the outer part of the disk, the brightness temperature for Band 3 is higher than those for Bands 4 and 7 in the inner region (r proj 20 au). The location where the emission becomes optically thick is closer to the central protostar at longer wavelengths, resulting in higher brightness temperature for longer-wavelength emission.
The brightness temperature for Band 4 seems to be more reliable for tracing the disk temperature in the region of 30 au r proj 50 au compared to that for Band 7 because α 3mm−2mm ∼ 2 in this region and the Band 7 emission are contaminated by the envelope component. The similar brightness temperatures for Band 3 and Band 4 suggest that the obscuring effect is insignificant. The sum of the brightness temperature and T bg = 2.7 K reflects the disk temperature at the midplane because the emission is optically thick. Note that the beam size in the direction normal to the midplane is sufficiently small to resolve the vertical distribution of the disk. In contrast, the disk temperature at r proj 50 au (outer part of the disk) cannot be determined using the brightness temperature for Band 4 because the emission would become marginally optically thick. However, the disk temperature should be 10 K in the outer part because the brightness temperatures for the Band 4 and 7 emission are T B < 10 K. On the other hand, van 't Hoff et al. (2018) suggested the warm disk scenario in L1527 by the non-detection of the N 2 D + (J = 3 − 2) emission and the detection of the optically thick emission of the 13 CO (J = 2 − 1) and C 18 O (J = 2 − 1) lines. The disk temperature is suggested to be above 20 K around a radius of 75 au. However, the spatial resolutions (∼ 25 − 130 au) of these molecular line observations may not be sufficient to resolve the disk. In addition, the CO gas would be emitted from the disk surface or envelope. Therefore, it may be the case that the midplane of the disk is cold (T ∼ 10 K), while the disk surface and the envelope are warm (T 20 K).
The radial distributions of the brightness temperature show a steep power-law profile in all bands. The powerlaw index of r −2.0 was found by fitting the brightness temperature for Band 4 in r proj of the ∼ 20 − 90 au region shown by the red lines in Figure 10. The disk temperature is determined by mainly two heating mechanisms: irradiation and accretion heating. Irradiation from the protostar directly heats the surface and determines the bulk disk temperature (Chiang & Goldreich 1997). Accretion heating is caused by viscous dissipation mediated by turbulence (Shakura & Sunyaev 1973). The disk temperature of these heating mechanisms are suggested to have power-law profiles. In the irradiation heating, the disk temperature shows T ∝ r −3/7 for optically thick and T ∝ r −1/2 for optically thin radiation (Chiang & Goldreich 1997;Chiang et al. 2001). In the viscous heating, the disk temperature is suggested to be T ∝ r −9/10 (Baumann & Bitsch 2020).
The slope of r −2.0 proj in the L1527 disk is much steeper than those for the irradiation model (r −3/7 ) (Chiang & Goldreich 1997;Chiang et al. 2001), optically thin disk model (r −1/2 ), and viscous heating model (r −9/10 ) (Baumann & Bitsch 2020). The irradiation temperature model is also plotted in Figure 10. It was calculated by adopting a luminosity of 2.75 L and a stellar mass of 0.45 M , respectively, as done by Baumann & Bitsch (2020). We note that the viscous heating model has a large uncertainty for the temperature profile because of various dependencies, such as the radial mass inflow rate and dust opacity. One may consider that the steep temperature profiles are caused by the resolved-out effects due to the ALMA extended configurations. We discuss this effect in Appendix A using simulated observations, based on which we conclude that the steep temperatures are not caused by resolved effects.
The temperature distributions in protoplanetary disks obtained using various molecular lines have recently been reported by the ALMA Large Program MAPS (Öberg et al. 2021). Using MAPS data, Law et al. (2021) showed that the temperature at a disk radius of 100 au is 20 − 30 K with a radial power-law index of q = −0.01 to −0.23. In contrast, the L1527 disk shows temperatures of 5−10 K with a radial power-law index of q = −2.0. Although these temperatures were estimated using different observations of the dust continuum emission and molecular lines, the temperature distribution for the L1527 Class 0/I protostellar disk is very different from those for the MAPS protoplanetary disks. The L1527 disk is colder and has a steeper power-law temperature profile than that for the protoplanetary disks in the outer region (r proj 20 au) of the disk.
To investigate the temperature distribution within r proj 20 au, we plot the brightness temperatures for Bands 3 and 4 as a function of the radius in Figure  Figure 11. Brightness temperatures for Band 3 and 4 emission versus radius based on higher-spatial-resolution data shown in Figure 2. The gray lines show the irradiation temperature model (T ∝ r −3/7 proj ) and the red lines show the power-law best-fit model (T ∝ r −2 proj ). The inner part of the disk is affected by beam dilution, as indicated by the gray regions. The error bars represent ±1σ, which was determined using the 5% flux error and the rms values of the noise levels.
11 based on the higher-resolution data obtained with a robust parameter of −2. As shown, the change in the power-law index around r proj ∼ 20 au is not caused by beam dilution. The temperature profile within r proj 20 au seems to be consistent with the irradiation model for Band 3; the slope of the temperature changes after r proj ∼ 20 au. Interestingly, the transition radii in the temperature slopes spatially coincide with the locations of the VLA clumps. Furthermore, the southern part of the outer disk region (r proj 20 au) is systematically colder than the northern part. Taking into account that the disk is locally flared in the south VLA clump, the sudden drop in temperature is caused by shadowing of the VLA clumps.
The shadowing effect in protoplanetary disks has been investigated using radiative transfer calculations (e.g., Dullemond et al. 2001;Dullemond & Dominik 2004;Ueda et al. 2019;Ohno & Ueda 2021;Okuzumi et al. 2022). These simulations assumed that a shadow is caused by a planet-induced gap or a dust pile-up at an inner rim due to a dead zone. They showed that the temperature drops in the shadow region, which is consistent with our result. However, for the simulated temperature, the power-law profile returns to normal with T ∝ r −0.5 in the outer radius. This is because the outer region is flared and the stellar radiation can heat the flared outer region in the simulations. The absorbed energy at the flared surface is transferred toward the midplane via diffusion. In contrast, the disk temperature for L1527 maintains a steep slope, T ∝ r −2 proj , to the outer edge of the disk, suggesting that the entire disk may be shadowed by the VLA clumps.
It is worth mentioning that the shadowing of the envelope by the embedded disk was reported in the Class 0/I protostars in the VLA 1623−2417 region (Murillo et al. 2015). They found a cold ring-like structure at the disk-envelope interface in the DCO + (J = 3 − 2) emission and suggested that the temperature drop is caused by the shadowing due to the disk. These results suggest that the shadowing effect may be a common possible structure in the disks as well as the envelopes.

Disk-Scale Height
The dust-scale height of the disk is an important parameter for constraining the shadowing effect and disk evolution. The almost edge-on disk geometry of the L1527 disk allows the vertical structure to be investigated. The projected dust-scale height (H proj ) was derived from the vertical intensity profile with the standard deviation (σ) of Gaussian fitting. Figure 12 shows the projected dust-scale heights derived from the Band 3 and 4 emission for the higher-sensitivity images. The increases in the dust-scale heights are confirmed in these plots. The power-law index is roughly estimated to be H proj ∝ r 1.0 proj , which is similar to that for previous observations (r 1.2−1.3 ) (Tobin et al. 2013;Aso et al. 2017).
It has been suggested that dust grains in protoplanetary disks settle onto the disk midplane due to grain growth and low turbulence (Dubrulle et al. 1995;Youdin & Lithwick 2007). The ratio of the dust-and gas-scale heights allows us to investigate the turbulence strength and grain growth (e.g., Pinte et al. 2016;Ohashi & Kataoka 2019;Villenave et al. 2020;Doi & Kataoka 2021). Therefore, we also derive the projected gas-scale heights for comparison. The gas-scale height was calcu-lated as c s /Ω K , where c s is the sound speed, Ω K is the Keplerian frequency, G is the gravitational constant, and M is the central stellar mass, respectively. We assume that the shadowed disk temperature (T ∝ r −2 ) is that shown by the red lines in Figure 10. To derive the projected gas-scale height, we performed radiative transfer calculations using RADMC-3D (Dullemond et al. 2012) assuming an inclination of 5 • and a position angle of 5 • . Based on the obtained image smoothing with the beam sizes of the Band 3 and 4 observations, we derived the gas-scale heights for the Band 3 and 4 images. The setup was the same as that in Ohashi & Kataoka (2019) but the grain size was simply assumed to be a max = 10 µm because α 7mm−3mm suggests that the dust grain size is not large. According to the spectral index and grain size relation (Birnstiel et al. 2018), the grain size is expected to have a range of a max ∼ 0.1 − 100 µm. The emission region was assumed to be optically thin (τ < 0.1) regardless of the radius to avoid the obscuring effect of the edge-on disk. Note that the gas-scale height is independent of the surface density as long as the emission region is optically thin. Figure 12 plots the projected gas-scale height (purple solid lines). We found that the gas-scale height is mostly consistent with the dust-scale height at r proj 50 au, which is consistent with the steep temperature gradient of T ∝ r −2 and suggests that the dust grains have not settled onto the midplane. No vertical dust settling is consistent with another protostellar disk around the Class 0 protostar HH 212 . The slight differences between the dust-and gas-scale heights in the inner region (r proj 25 au) may be caused by the high optical depth associated with Band 3 and 4 emission. The linewidth of the Gaussian fit is larger due to saturation of the emission profile toward the inner radius of the disk midplane. The substructures of the VLA clumps also contribute to the larger dust-scale heights because additional components are found in the vertical direction, as discussed in Section 3.2 and shown in Figures 2 and 3.
In the outer part of the disk (r proj 60 au), the dustscale height seems to be larger than the gas-scale height. Sakai et al. (2017) found a similar trend, with the CCH emission showing a flared disk structure. Their interpretation of this result is that the accretion gas stagnates around the centrifugal barrier and moves in vertical directions. This scenario is possible for the dust because the gas and dust motions are coupled in the disk-forming region.
The increase in dust-scale height may also be caused by inclination, temperature, and outflow contamination. For example, if the disk is slightly face-on at r ≥ 60 au, the projected scale height becomes larger. This scenario might be consistent with the warped structure because the transition of the inner and outer orbits has been suggested to be r proj ∼ 40 − 60 au (Sakai et al. 2019). It is reasonable that the warped structure may change the disk inclination. The temperature distribution may also increase the scale height of the gas and dust because the accretion shock caused by infalling materials from the envelope may raise the temperature around the centrifugal barrier (Sakai et al. 2014a). The outflow emission may also affect the dust distribution because energetic outflows have been observed toward L1527 (Tamura et al. 1996;Hogerheijde et al. 1998;Tobin et al. 2008Tobin et al. , 2010Flores-Rivera et al. 2021). These outflows are launched perpendicular to the disk elongation. Recent threedimensional magnetohydrodynamical simulations show that dust grains can also be launched by outflows, similar to gas motion (Tsukamoto et al. 2021). Therefore, it is possible that the dust-scale height is affected by the entrainment of outflows. To investigate these scenarios, detailed modeling and molecular line observations are needed to reveal the gas kinematics in the disk.
To shadow the outer disk (r proj 20 au), the aspect ratio needs to decrease with a radius. Figure 13 shows the observed aspect ratio for the disk as a function of the radius. In the outer disk (r proj 30 au), the aspect ratio seems to be constant (∼ 0.2). This is consistent with the increasing trend of the dust-scale height (H ∝ r 1.0 ). In contrast, the aspect ratio shows a decreasing trend in the inner region (r proj 30 au). Even though the spatial resolution may be insufficient to precisely measure the aspect ratio in the inner 10-au region, the positions of the VLA clumps (r proj ∼ 10 − 20 au) have larger aspect ratios than in the outer region. This is confirmed by the discovery of a substructure in the vertical direction, as discussed in Section 3.2 (e.g., Figures 2 and 3). Therefore, the shadowing by the VLA clumps is a possible mechanism.

Grain Size in Protostellar Disk
It is important to investigate the grain size in the disks to reveal the first step of planet formation. One way of measuring the dust grain size is to derive the frequency dependence of the thermal dust continuum emission, because larger grains efficiently emit thermal radiation at a Figure 12. Dust-scale height derived from Band 3 and 4 images versus projected radius (rproj). The higher-sensitivity images created with a robust parameter of 0.5 were used. The purple lines indicate the gas-scale height estimated using the images from RADMC-3D assuming a power-law best-fit temperature (T ∝ r −2 ). The green and blue circles represent the north and south directions, respectively. wavelength similar to their size (e.g., Draine 2006;Ricci et al. 2010;Testi et al. 2014;Kataoka et al. 2014;Birnstiel et al. 2018). The dust opacity index, β, provides information on grain size.
The dust opacity index was calculated as where τ and ν are respectively the optical depth and frequency at each band. To obtain the dust opacity index, the optical depth need to be measured at multiple wavelengths. Here, we used the emission of Band 3 and Q-Band to derive the optical depth because the emission at Bands 4 and 7 was thought to be highly optically thick at r proj 50 au. We assume that the disk temperature is the sum of the brightness temperature for Band 4 and the background temperature of 2.7 K (T B + T bg ). Data for all bands were smoothed to a ∼ 30-au resolution to detect the Q-Band emission. The spatial resolution is limited by the natural weighting of the VLA Q-Band data. Figure 14 shows radial plots of the brightness temperature, optical depth, and dust opacity index. In the central region (r proj 30 au), the optical depth for the Band 3 emission cannot be derived because the Band 3 temperature is higher than that of Band 4. The optical depths for Band 3 and Q-Band in the outer region were derived by using Equation (2), then the dust opacity indices were derived by using Equation (5). The derived Band 3 emission is moderately optically thick (τ 3mm ∼ 1−2), which is consistent with the low spectral index (α 3mm−2mm ∼ 2.0). We found that the Q-Band emission is also optically thick (τ 7mm ∼ 1.2) around r ∼ 30 au, and become optically thin (τ 7mm ∼ 0.4 − 0.6) at r proj 40 au.
The derived opacity index based on these optical depths is β ∼ 1.67 +1.2 −0.7 . Even though the uncertainty is large, the dust opacity index, β ∼ 1.67, is the same as that for the ISM, which constrains the grain size to a max 100 µm by assuming that the dust grains are simple spherical particles without porosity (Birnstiel et al. 2018). Here, we note that grain size estimation strongly depends on dust models. For example, dust grains with porosity keeps β ∼ 1.67 up to a max 1 cm (Kataoka et al. 2014;Birnstiel et al. 2018). However, it is difficult to proceed the grain growth in such a short time scale in the L1527 protostellar disk because no substructure induced by dust accumulation is found in the outer part of the disk, r proj 40 au. Therefore, we suggest that significant grain growth has not yet begun in the outer region (r proj 50 au) of the disk.
Recent coagulation simulations show that the dust grains become larger from the inner region to the outer region of a disk because the growth timescale is roughly proportional to the orbital period (Ohashi et al. 2021;Kobayashi & Tanaka 2021). Ohashi et al. (2021) showed that the critical radius (R c ) where grain growth proceeds can be calculated as where the protostellar mass M , gas-to-dust mass ratio ζ d , and disk age t disk are roughly the same values as those for L1527. This means that grain growth does not occur in the outer region (r 22 au), which is consistent with the large dust opacity index found around 50 au. Figure 14. Upper panel shows brightness temperatures for Band 4 and 3 and Q-Band emission versus radius. Middle panel shows radial plots of optical depths for Band 3 and Q-Band emission derived by assuming that dust temperature corresponds to Band 4. Lower panel shows dust opacity index, β, derived from optical depths for Band 3 and Q-Band emission. The error bars represent ±1σ, which was determined using the 5% flux error and the rms values of the noise levels.

Gravitational Instability and Substructure Formation
The origin of the VLA clumps was discussed in Nakatani et al. (2020); Ohashi et al. (2021). Due to the edge-on disk geometry, the substructures can be explained either by ring or spiral structures. Nakatani et al. (2020) discussed the snow line ring because the disk temperature of T ∼ 60 K at the VLA clumps coincides with the CO 2 snow line, while Ohashi et al. (2021) suggested that the radius of r proj ∼ 22 au at the VLA clumps coincides with the growth front ring. However, neither the snow line nor the growth front mechanism explain the shadowed region outside the ring. To shadow the outer region, the density needs to be enhanced in the VLA clumps.
The local enhancement in density seems to be consistent with the material spiral arms formed via gravitational instability (Tomida et al. 2017). Such spiral arms form in marginally unstable disks, where Toomre's Q parameter (Toomre 1964) is as low as 2.
The Q parameter is defined as where c s is the sound speed, Ω K is the Keplerian frequency, and Σ is the surface density of the disk. To measure the Q parameter, the surface density, Σ, needs to be obtained. It can be constrained by our observations. The dust surface density, Σ dust , was calculated using where τ 7mm is the optical depth of the 7-mm dust continuum emission and κ 7mm is the absorption opacity at a wavelength of 7 mm. Equation (8) indicates that the optical depth and the absorption opacity are needed for the surface density. First, we estimate the optical depth for the Q-Band emission by assuming that the disk temperature corresponds to the brightness temperature of the Band 3 emission. This assumption is reasonable for the VLA clump regions because the Band 3 and 4 emission is optically thick and show similar temperatures that follow the irradiation disk model. Figure 15 shows radial plots of the Band 3 and Q-Band emission by smoothing the beam size of 0. 1. Because the Q-Band emission has contaminations of the dust thermal emission and free-free emission, Figure 15 also plots the K-Band 1.3 cm emission to investigate the free-free contamination. The K-Band data are taken from Nakatani et al. (2020). Then, the subtraction of the free-free emission in the Q-Band data is calculated as where the adopted spectral index, −0.1, is typical for optically thin free-free emission (e.g., Anglada et al. 2018) and is the case for L1527 (Nakatani et al. 2020). The optical depth of the Q-Band dust thermal emission around the VLA clumps is derived to be τ 7mm ∼ 1.5 after the subtraction of the free-free emission. Even though the temperature should be measured more precisely in future observations, the similar brightness temperatures for Band 3 and Q-Band suggest that the 7 mm emission is optically thick at the clumps (τ 7mm > 1). Figure 15. Radial plots of brightness temperature of Band 3, Q-Band, and K-Band emission. The brightness temperature was derived using the higher-spatial-resolution images with a robust parameter of −2 shown in Figure 2, and smoothing to 0. 1 spatial resolution. The error bars represent ±1σ, which was calculated using the 5% flux error and the rms values of the noise levels.
Second, we consider the absorption opacity, κ 7mm , which depends on the grain size and dust components (e.g., Draine 2006). We found that the dust opacity index, β, is the same as that for the ISMs in the outer region. However, the grain size in the VLA clumps may grow due to local enhancement of the density. Because the grain size could not be determined, we applied the upper limit of the standard dust opacity model given in Birnstiel et al. (2018), which constrains the lower limit of the dust surface density. Figure 16 (a) shows the dust absorption opacity as a function of grain size at an observation wavelength of 7 mm (Birnstiel et al. 2018). As reproduced in the top panel of Figure 16, the model given by Birnstiel et al. (2018) predicts the maximum value κ 7mm ∼ 0.11 g −1 cm 2 for a dust size distribution of q = 3.5 and κ 7mm ∼ 0.16 g −1 cm 2 for that of q = 2.5.
Finally, we estimate the surface density and discuss Toomre's Q parameter. In the above discussion, the optical depth is constrained to be τ 7mm ∼ 1.5 in the VLA clump regions. In addition, the absorption opacity is Figure 16. Upper panel shows dust absorption opacity models for dust size distribution for q = 3.5 and q = 2.5 at observation wavelength of 7 mm (Birnstiel et al. 2018). Lower panel shows possible dust surface density depending on dust absorption opacity in VLA clumps. The dotted lines represent the dust surface densities where Q = 0.1 and Q = 1, respectively. The Q value is Q 0.25 for q = 3.5. assumed to be governed by the model shown in Figure  16 (a). Based on these parameters and Equation (8), we estimated the dust surface density as shown in Figure  16 (b), which constrains the lower limit of the surface density to Σ dust 14 g cm −2 for a dust grain size distribution of q = 3.5. Even if we assume a dust grain size distribution of q = 2.5, the lower limit is Σ dust 9.6 g cm −2 . By using the above dust surface densities and assuming a dust-to-gas mass ratio of 0.01, the Q value can be calculated using Equation (7). The Q value was found to be Q 0.25 for q = 3.5 and Q 0.36 for q = 2.5. We show surface densities where the Q values are Q = 0.1 and Q = 1, respectively, in Figure 16 (b). The upper limit of Q 0.25 indicates that the disk is gravitationally unstable at the VLA clumps.
It should be noted that the surface density used to estimate the Q value is that for a face-on disk. In contrast, the surface density used here to estimate the Q value is that for an edge-on disk. It is difficult to estimate the surface density for this disk when observed face-on. However, if the clumps have a spherical structure, the surface density will be mostly constant. Even if the clumps are ring structures, the difference in the surface density observed edge-on and face-on would be similar to the aspect ratio (H/r ∼ 0.25 − 0.35). In this case, the Q value is Q 1.0 for q = 3.5. Therefore, we conclude that the disk is gravitationally unstable in the VLA clump locations.
The gravitationally unstable condition suggests that the VLA clumps are formed by gravitational instability, such as that caused by spiral arms and fragmentation. Spiral arms are commonly observed in gravitationally unstable disks (Rice et al. 2003;Kratter & Lodato 2016). We investigated whether spiral arms can be observed as clumps in an edge-on disk as similar to the VLA clumps in the L1527 disk. We applied the toy model of spiral arms by Sakai et al. (2019) in which a geometrically thick disk contains two loosely wound spiral arms. The assumed density distribution is given in Appendix B. The left panel of Figure 17 displays the spiral pattern by the surface density obtained by integrating the density along the disk vertical direction. The right panel denotes the column density along the line of sight when we observe the model disk at inclination of 5 • and the position angle of 5 • . Both the surface density and column density are normalized by their maximum values. The column density has three clumps along the disk as in the brightness distribution at the VLA Q-Band. Thus, the spiral arm model is applicable to the inner part of the disk.
Fragmentation may also occur in the spiral arms and produce several clumps if the disk is sufficiently massive. Takahashi et al. (2016) investigated the conditions required for fragmentation in a disk using two-dimensional numerical simulations. They suggested that Toomre's Q parameter in the spiral arms is essential for fragmentation. The spiral arms fragment when Q < 0.6 in the spiral arms. If the VLA clumps are fragments in spiral arms rather than a ring structure, the derived Q value is Q 0.25, consistent with the fragment threshold of Q < 0.6. These results suggest that the substructure observed in the VLA can be the clumps formed by gravitational instability.
If the VLA substructures (Clump-N and S) are real clumps, the estimation of these masses is important for planet formation or disk evolution. We estimated the dust mass of the VLA clumps (Clump-N and S) using where F 7mm is the flux density, κ 7mm is the dust opacity, and B ν is the Planck blackbody function. We used the maximum value κ 7mm ∼ 0.11 g −1 cm 2 for the dust opacity to constrain the lower limit of the dust mass. The derived dust mass is M dust 0.12 M J and M dust 0.17 M J for the northern and southern clumps, respectively. Here, the thermal dust component at 7 mm is used by subtracting the free-free emission component (Nakatani et al. 2020). The existence of planets in Class II protoplanetary disks was suggested by a kinematic analysis of molecular lines (Pinte et al. 2018;Teague et al. 2018). In addition, a circum-planetary disk has been identified in the protoplanetary disk around PDS 70 (Keppler et al. 2018). These results indicate that planet formation starts in Class I disks or at a much earlier stage. Indeed, recent ALMA observations of Class I protostellar disks indicated that planet formation may already occur in the Class I disk-forming stage (Segura-Cox et al. 2020;Alves et al. 2020). Our results may support such early planet formation in Class 0/I protostellar disks. The L1527 disk can be massive enough to be enhance the gravitationaly instability, which might induce planet formation in future (e.g., Gammie 2001). Xu (2022) also suggested that the majority of disks are likely to be gravitationally unstable.
Rather than the planet formation, it may also be possible for the VLA clumps to accrete to the central star and cause an accretion burst such as in the case of FU Orionis, which is a low-mass young stellar object that shows an accretion outburst with an accretion rate that increases fromṀ ∼ (10 −7 − 10 −8 ) M yr −1 to (10 −5 −10 −4 ) M yr −1 over several decades to 100 years (Audard et al. 2014). The masses of the VLA clumps are sufficient to increase the accretion rate if these clumps accrete to the central star via migration (e.g., Lin & Papaloizou 1986;Ward 1997;Tanaka et al. 2002;Vorobyov & Basu 2015). Thus, the VLA clumps might be the origin of future accretion bursts.
As summarize the discussion of this subsection, the L1527 disk is found to be gravitationally unstable (Q 1.0). Therefore, the VLA clumps would be formed by the disk gravitational instability. These VLA clumps may be the origin of gas-giant planets or accretion burst by accreting onto the central star.

Caveat and Other Possibilities for Substructure Formation
In the above discussion, we showed that the disk is gravitationally unstable because the Q value is Q 1.0. However, the disk surface density for estimating the Q value strongly depends on the dust opacity (κ abs ). In this study, we used the standard opacity model developed by Birnstiel et al. (2018). There are some opacity models that have larger dust opacities, such as that in Jager et al. (1998);Draine (2003); Zubko et al. (1996) and some carbonaceous dust models in Birnstiel et al. (2018). If we apply these opacity models, it is possible for the disk to be gravitationally stable with Q > 2.0 only when the dust grains are larger. Although we suggested that gravitational instability is the likely origin of the substructure, other mechanisms may have created the shadowing effect at r ∼ 20 au. By taking into account the shadowing by the VLA clumps, a dust pile-up mechanism is needed to form the disk substructures such as the dead zone (Flock et al. 2015). Further observa-tions with higher spatial resolution and sensitivity are needed to confirm the origin of the substructures.

SUMMARY
We have analyzed multi-wavelength dust continuum emission toward the protostellar disk around the Class 0/I protostar L1527 IRS. The dust continuum data cover a wide wavelength range of 0.87, 2.1, 3.3, and 6.8 mm obtained by the ALMA and VLA observations, and have angular resolutions of 0. 025 − 0. 16 (corresponding to 3.4 − 20 au). The main results are listed below.
• The protostellar disk was shown to be an edgeon structure aligned in the north-south direction based on the ALMA observations. This structure is similar to that derived from previous observations (e.g., Aso et al. 2017;Sakai et al. 2017). In contrast, the VLA observations show that the 6.8mm emission is only detected within 30 au, which is attributed to the sensitivity limit.
• The peak temperature increases with increasing wavelength, which is explained by the optical depth effect; emission at shorter wavelengths is optically thicker and is obscured by the outer cold annuli. In contrast, the higher brightness temperature found for the Band 7 emission in the outer part of the disk (r proj 50 au) suggests the existence of an envelope component or accretion shock to the edge of the infant disk.
• Based on an analysis of the optical depths associated with the Band 3 and Q-Band emission, the dust opacity index, β, was determined to be β ∼ 1.7 at r proj ∼ 50 au, suggesting that significant dust grain growth has not yet begun. This is consistent with the growth front model developed by Ohashi et al. (2021) because the growth front model predicts that the critical radius (R c ) where the grain growth proceeds can be calculated as where the protostellar mass M , gas-to-dust mass ratio ζ d , and disk age t disk are roughly the same values as those for L1527. This means that grain growth does not occur in the outer region (r 22 au).
• The disk temperature seems to show a power-law profile with an index of T ∝ r −3/7 proj inside a 20au radius, and a steeper power-law profile with an index of T ∝ r −2 proj at r proj 20 au. The inner temperature can be explained by an irradiation model, while the steep drop in temperature toward the outer region is caused by a shadowing effect by the VLA clumps.
• The VLA clumps form via gravitational instability because the Q value is estimated to be Q 1.0. The derived dust mass for the VLA clumps is 0.1 M J . Thus, we suggest that Class 0/I disks can be sufficiently massive to be gravitationally unstable, which might be the origin of gas-giant planets in a 20-au radius.
We gratefully appreciate the comments from the anonymous referee that significantly improved this article. This paper makes use of the following ALMA data: In Section 4.2, we discuss the steep temperature of T ∝ r −2 due to the shadowing. However, one may consider that the steep temperatures are caused by the resolved-out effect of the ALMA extended configurations.
To investigate the resolved-out issue, we conduct the simulation observations in CASA simobserve at ALMA Band 4. The disk model is created by RADMC-3D (Dullemond et al. 2012). The temperature and dust surface density are assumed to be T = 224(r/AU) −3/7 and Σ dust = 10 g cm −2 , respectively, to mimic the observed optically thick emission. The disk size is assumed to be r out = 100 au. Note that the temperature profile assumes the irradiation model, same as Figures 10 and 11, but the dust surface density should be made with any optically thick model. We simply assume a constant dust surface density independent of radius. The dust size is simply assumed to be a max = 0.1 µm because the grain growth does not proceed yet discussed in Section 4.4. The inclination is set to be 15 • to mimic the large dust scale height.
The CASA simobserve is applied to the model image with a total integration time of 120 min with the most extended configuration in Cycle 8 (alma.cycle8.10.cfg), which is the same observing setup with our observations. A Precipitable Water Vapor (PWV) of 0.5 mm is used to input the noise. After the measurements sets (MS file) are created by simobserve, we map the simulation image in CASA tclean. Figure 18 (a) and (b) show the model image and the radial plot of the normalized intensity, respectively. Furthermore, Figure 18 (c) and (d) shows the same image and radial plot but for the simulation observations. The beam size of the simulation observations is 0. 049 × 0. 036 shown in the bottom-left corner of Figure 18 (c). Figure 18 (d) also plots the radial distribution of the model smoothed by the beam size for comparison. Both of the model and simulation observations show almost the same intensity distributions, although some intensity fluctuations and flatter power-law profile are found in the simulation observations due to the noise and beam convolution. Because intensity profile does not change by the simulation observations, we conclude that the observed steep temperature profile is not caused by the resolved-out issue.

B. SPIRAL MODEL
Our spiral model is essentially the same as that of Sakai et al. (2019), although the model parameters are slightly different. To describe the model, we use the Cartesian coordinates, (x, y, z), of which origin coincides with the protostar. The x-axis is parallel to the disk major axis on the sky, while the z-axis is perpendicular to the disk plane. The inclination and position angle of the disk are assumed to be i = 5 • and P.A. d = 5 • , respectively. Then the Cartesian coordinates are related to the observation by where D and s denote the distance to L1527 and coordinate along the line of sight, respectively. For later convenience, we define the spherical coordinates, The right panel denotes the column density along the line of sight, N (∆α, ∆δ) = ρ(x, y, z)ds, which mimics the Q-band image. The surface density is normalized by the maximum value on each panel. The contours denote log 10 Σ(x, y)/ max Σ and log 10 N (−∆α, ∆δ)/ max N in the left and right panels, respectively.