Formation of Dust Clumps with Sub-Jupiter Mass and Cold Shadowed Region in Gravitationally Unstable Disk around Class 0/I Protostar in L1527 IRS

The L1527 protostellar disk was observed with ALMA Band 4 on 2021 September 15 and 21, 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₂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 minutes 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 15,238 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 ∼15 (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.

The L1527 protostellar disk was also observed with ALMA Band 7 on 2021 September 26, 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 minutes 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 14,361 m. The maximum recoverable scale (θ_MRS) was ∼16 (corresponding to ∼220 au). The data reduction was performed using a pipeline script provided by ALMA.

We performed the self-calibration of the continuum observations in CASA. The entire imaging process was carried out with the CASA tclean task. To construct CLEAN maps, we adopted Briggs's weighting, with robust parameters 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.

Figure 1 shows the dust continuum images toward the L1527 disk obtained at submillimeter 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 that it is elongated in the north–south direction. Because the disk has an almost edge-on 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 decl. (r_proj ² = Δx² + Δy²). The radial direction is mostly consistent with the Δx direction.

Zoom In Zoom Out Reset image size

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 Bands 7, 4, and 3, and the Q-band emission, are ∼40, ∼60, ∼70, and ∼50 K, respectively. The highest temperature (∼70 K) is found for the Band 3 emission, which suggests that the 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 rms noise levels are shown in Table 1. The total flux for the Band 7 emission (430 ± 11 mJy) is slightly lower than that of the 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 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 the 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's 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.

Zoom In Zoom Out Reset image size

The most prominent substructures were identified in the VLA Q-band observations. Three clumps (clump-N, clump-C, and clump-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 (the 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 the absolute position error. We found that the ALMA images have no counterparts 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 shorter wavelengths suggest that the shorter-wavelength emission is affected by the higher optical thickness in regions at the 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's 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 coincides 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.

Zoom In Zoom Out Reset image size

The spectral index, α, was calculated as

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{\mathrm{ln}{I}_{1}-\mathrm{ln}{I}_{2}}{\mathrm{ln}{\nu }_{1}-\mathrm{ln}{\nu }_{2}},\end{eqnarray} \tag{ 1 }$$

where I and ν are the intensity and frequency at each band, respectively. The intensity (I_ν ) of the dust continuum emission was calculated using the radiative-transfer equation:

$$\begin{eqnarray}&&{I}_{\nu }=[{B}_{\nu }(T)-{B}_{\nu }({T}_{\mathrm{bg}})][1-\exp (-{\tau }_{\nu })],\end{eqnarray} \tag{ 2 }$$

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 the 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 multifrequency synthesis for the combination of ALMA Band 3 and the 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 over the entire disk indicates that the low spectral index is due to the optical depth effect, rather than to 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 the emission at shorter wavelengths traces a colder region. The optically thick emission allows only the foreground annulus of the edge-on 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 Figures 1 and 2. This effect was previously reported by Galván-Madrid et al. (2018).

Zoom In Zoom Out Reset image size

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 (ISM; 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.

3.3.1. 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 confirm 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 nondetection 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 of the ISM, 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.

Zoom In Zoom Out Reset image size

3.3.2. Vertical Dependence

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 over 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°.

Zoom In Zoom Out Reset image size

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 the 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 side-lobe effect. The point-spread function for our observations shows that the side-lobe levels are ∼0.2. Because the peak emission was detected with a signal-to-noise 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 portion from the midplane was caused by a slight inclination of the disk. We discuss disk inclination in the following section (Section 4.1).

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 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) have 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).

Zoom In Zoom Out Reset image size

Figure 8 shows the profiles of the brightness temperatures in Bands 4 and 3 along the disk's vertical direction (the 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 the 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 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 (Lin et al. 2021), 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).

Zoom In Zoom Out Reset image size

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 (the near side) due to the edge-on geometry, and it may be emitted from the vicinity of the protostar toward the western side (the 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 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 western 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.

Zoom In Zoom Out Reset image size

The multiwavelength 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 the disk temperature (see Equation (2)). Figure 10 shows the brightness temperatures of 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 calculated from the intensity F_d (in units of mJy beam⁻¹) without the Rayleigh–Jeans approximation.

Zoom In Zoom Out Reset image size

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 the Band 7 emission extending beyond the disk component, as discussed in Section 3.1, these results suggest that the Band 7 temperature is contaminated by the temperatures 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, then 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 at the outer part of the disk, the brightness temperature for Band 3 is higher than the temperatures 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 a 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 (the 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 from the nondetection of the N₂D⁺ (J = 3 − 2) emission and the detection of the optically thick emission of the ¹³CO (J = 2 − 1) and C¹⁸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 power-law 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 mainly determined by 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 temperatures of these heating mechanisms are suggested to have power-law profiles. In 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 viscous heating, the disk temperature is suggested to be T ∝ r^−9/10 (Baumann & Bitsch 2020).

The slope of ${r}_{\mathrm{proj}}^{-2.0}$ in the L1527 disk is much steeper than the slopes for the irradiation model (r^−3/7) (Chiang & Goldreich 1997; Chiang et al. 2001), the optically thin disk model (r^−1/2), and the 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 resolved-out effects due to the ALMA extended configurations. We discuss these effects 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 those 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 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 the disk being locally flared in the southern VLA clump, the sudden drop in temperature is caused by the shadowing of the VLA clumps.

Zoom In Zoom Out Reset image size

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 pileup at an inner rim due to a dead zone. They showed that the temperature drops in the shadow region, which is consistent with our results. 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\propto {r}_{\mathrm{proj}}^{-2}$, 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 for the Class 0/I protostars in the VLA 1623−2417 region (Murillo et al. 2015). They found a cold ringlike structure at the disk–envelope interface in the DCO⁺ (J = 3−2) emission and suggested that the temperature drop was 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.

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 a 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}_{\mathrm{proj}}\propto {r}_{\mathrm{proj}}^{1.0}$, which is similar to that for previous observations (r^1.2−1.3; Tobin et al. 2013; Aso et al. 2017).

Zoom In Zoom Out Reset image size

It has been suggested that the 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 calculated as

$$\begin{eqnarray}&&{c}_{s}/{{\rm{\Omega }}}_{{\rm{K}}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{K}}}=\sqrt{\displaystyle \frac{{{GM}}_{\star }}{{r}^{3}}}=2.0\times {10}^{-7}{\left(\displaystyle \frac{r}{1\,\mathrm{au}}\right)}^{-3/2}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{1/2}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 4 }$$

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⁻²) 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 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 }\sim 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 (the 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⁻² 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 (Lin et al. 2021). 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 the Band 3 and 4 emission. The linewidth of the Gaussian fit is larger, due to the 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 dust scale 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 was 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 the 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. 2008, 2010; Flores-Rivera et al. 2021). These outflows are launched perpendicular to the disk elongation. Recent three-dimensional 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 radius. Figure 13 shows the observed aspect ratio for the disk as a function of 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. A recent simulation also shows that the whole outer disk can be shielded because the aspect ratio (H/r) decreases outwards in a gravitationally unstable Class 0/I disk (Xu & Kunz 2021).

Zoom In Zoom Out Reset image size

It is important to investigate the grain size in the disks, to reveal the first steps 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 wavelength similar to their size (e.g., Draine 2006; Ricci et al. 2010; Kataoka et al. 2014; Testi et al. 2014; Birnstiel et al. 2018). The dust opacity index, β, provides information on grain size.

The dust opacity index was calculated as

$$\begin{eqnarray}&&\beta =\displaystyle \frac{\mathrm{ln}{\tau }_{1}-\mathrm{ln}{\tau }_{2}}{\mathrm{ln}{\nu }_{1}-\mathrm{ln}{\nu }_{2}},\end{eqnarray} \tag{ 5 }$$

where τ and ν are the optical depth and frequency at each band, respectively. To obtain the dust opacity index, the optical depth needs 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 that it becomes optically thin (τ_7mm ∼ 0.4−0.6) at r_proj ≳ 40 au.

Zoom In Zoom Out Reset image size

The derived opacity index based on these optical depths is $\beta \sim {1.67}_{-0.7}^{+1.2}$. 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 }\lesssim 100$ μm, by assuming that the dust grains are simple spherical particles without porosity (Birnstiel et al. 2018). Here, we note that the grain size estimation strongly depends on the dust models. For example, dust grains with porosity keep β ∼ 1.67 up to ${a}_{\max }\lesssim 1\,\mathrm{cm}$ (Kataoka et al. 2014; Birnstiel et al. 2018). However, it is difficult for the grain growth to occur over such a short timescale 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 (Kobayashi & Tanaka 2021; Ohashi et al. 2021). Ohashi et al. (2021) showed that the critical radius (R_c) where grain growth proceeds can be calculated as

$$\begin{eqnarray}&&\begin{array}{l}{R}_{{\rm{c}}}=22\ {\left(\displaystyle \frac{{M}_{\star }}{0.45\ {M}_{\odot }}\right)}^{1/3}{\left(\displaystyle \frac{{\zeta }_{{\rm{d}}}}{0.01}\right)}^{2/3}\\ \,\times \,{\left(\displaystyle \frac{{t}_{\mathrm{disk}}}{0.037\ \mathrm{Myr}}\right)}^{2/3}\ \mathrm{au},\end{array}\end{eqnarray} \tag{ 6 }$$

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.

The origin of the VLA clumps has been discussed in Nakatani et al. (2020) and 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) discuss the snowline ring, because the disk temperature of T ∼ 60 K at the VLA clumps coincides with the CO₂ snowline, while Ohashi et al. (2021) suggest that the radius of r_proj ∼ 22 au at the VLA clumps coincides with the growth-front ring. However, neither the snowline 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

$$\begin{eqnarray}&&Q=\displaystyle \frac{{c}_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}}}{\pi G{\rm{\Sigma }}},\end{eqnarray} \tag{ 7 }$$

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

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{dust}}=\displaystyle \frac{{\tau }_{7\mathrm{mm}}}{{\kappa }_{7\mathrm{mm}}},\end{eqnarray} \tag{ 8 }$$

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 shows 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 01. Because the Q-band emission has contaminations from 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). The subtraction of the free–free emission in the Q-band data is then calculated as

$$\begin{eqnarray}&&{I}_{7\mathrm{mm},\mathrm{dust}}={I}_{7\mathrm{mm}}-{I}_{1.3\ \mathrm{cm}}{\left(\displaystyle \frac{{\nu }_{7\mathrm{mm}}}{{\nu }_{1.3\mathrm{cm}}}\right)}^{-0.1},\end{eqnarray} \tag{ 9 }$$

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).

Zoom In Zoom Out Reset image size

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 ISM in the outer region. However, the grain size in the VLA clumps may grow, due to the 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⁻¹ cm² for a dust size distribution of q = 3.5 and κ_7mm ∼ 0.16 g⁻¹ cm² for a distribution of q = 2.5.

Zoom In Zoom Out Reset image size

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 assumed to be governed by the model shown in the upper panel of Figure 16. Based on these parameters and Equation (8), we estimate the dust surface density as shown in the lower panel of Figure 16, which constrains the lower limit of the surface density to Σ_dust ≳ 14 g cm⁻² 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⁻². 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 is found to be Q ≲ 0.25 for q = 3.5 and Q ≲ 0.36 for q = 2.5. We show the surface densities where the Q values are Q = 0.1 and Q = 1, respectively, in the lower panel of Figure 16. 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 to estimate the Q value here is that for an edge-on disk. It is difficult to estimate the surface density for this disk when it is 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 densities 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; Veronesi et al. 2021; Paneque-Carreno et al. 2021). We investigated whether spiral arms can be observed as clumps in an edge-on disk, 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 of the surface density obtained by integrating the density along the disk's vertical direction. The right panel shows the column density along the line of sight when we observe the model disk at an inclination of 5° and a position angle of 5°. Both the surface density and the column density are normalized by their maximum values. The column density has three clumps along the disk, as in the brightness distribution for the VLA Q band. Thus, the spiral arm model is applicable to the inner part of the disk.

Zoom In Zoom Out Reset image size

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 the 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 could be clumps formed by gravitational instability.

If the VLA substructures (clump-N and clump-S) are real clumps, the estimation of these masses is important for planet formation or disk evolution. We estimate the dust mass of the VLA clumps (clump-N and clump-S) using

$$\begin{eqnarray}&&{M}_{\mathrm{dust}}=\displaystyle \frac{{F}_{7\mathrm{mm}}}{{\kappa }_{7\mathrm{mm}}{B}_{\nu }({T}_{\mathrm{dust}})},\end{eqnarray} \tag{ 10 }$$

where F_7mm is the flux density, κ_7mm is the dust opacity, and B_ν is the Planck blackbody function. We use the maximum value κ_7mm ∼ 0.11 g⁻¹ cm² for the dust opacity, to constrain the lower limit of the dust mass. The derived dust masses are 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 obtained 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 circumplanetary 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 indicate that planet formation may already be occurring in the Class I disk-forming stage (Alves et al. 2020; Segura-Cox et al. 2020). Our results may support such early planet formation in Class 0/I protostellar disks. The L1527 disk could be massive enough to enhance the gravitational instability, which may 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 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 YSO that shows an accretion outburst with an accretion rate that increases from $\dot{{M}_{\star }}\sim ({10}^{-7}-{10}^{-8})$ M_⊙ yr⁻¹ to (10⁻⁵ − 10⁻⁴) M_⊙ yr⁻¹ over a period from several decades to 100 yr (Audard et al. 2014). The masses of the VLA clumps are sufficient to increase the accretion rate, should 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.

To 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's gravitational instability. These VLA clumps may be the origin of gas giant planets or accretion bursts, by accreting onto the central star.

In the above discussion, we have shown 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 have used the standard opacity model developed by Birnstiel et al. (2018). There are some opacity models that have larger dust opacities, such as those of Jager et al. (1998), Draine (2003), Zubko et al. (1996), and some of the carbonaceous dust models in Birnstiel et al. (2018). If we apply these opacity models, it is only possible for the disk to be gravitationally stable with Q > 2.0 when the dust grains are larger. Although we have 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 from the VLA clumps, a dust pileup mechanism is needed to form the disk substructures, such as a dead zone (Flock et al. 2015). Further observations at higher spatial resolution and greater sensitivity are needed to confirm the origin of the substructures.