© The Authors 2024

1 Introduction

Over the past ten years, the Atacama Large Millime-ter/submillimeter Array (ALMA) has made significant contributions to our understanding of protoplanetary disks by revealing substructures within them. These studies (e.g., ALMA Partnership 2015; Isella et al. 2016; Long et al. 2018) have shown the frequent appearance of substructure in proto-planetary disks with different morphologies (van der Marel et al. 2013; Pérez et al. 2016; Huang et al. 2018). These substructures play a significant role in protoplanetary disks: they allow dust grains to accumulate and grow, facilitating the formation of planets by providing regions of higher dust density and trapping, and they slow down the radial drift of dust grains, which would otherwise hinder the planet formation process (Pinilla et al. 2012a).

Most ALMA studies of protoplanetary disks are performed at wavelengths of about 1 mm, where the dust emission is relatively bright and high angular resolutions of about 30 mas can be achieved. In some disks, similar studies have also been performed at longer wavelengths, where the emission is known to be optically thinner, up to 3 mm with ALMA (e.g., Tazzari et al. 2021; Cazzoletti et al. 2018; Long et al. 2020), but also up to 1 cm with the Very Large Array (VLA; e.g., Carrasco-González et al. 2016; Hashimoto et al. 2023). These observations at long wavelengths allow us to study the properties of the dust grains, but at the expense of angular resolution. Observations of protoplanetary disks with ALMA at shorter wavelengths provide the highest possible angular resolution with ALMA (~10 mas), but they are usually avoided because of the higher atmospheric opacity and turbulence, which lead to necessarily longer integration and calibration times. Moreover, dust emission is expected to be optically thick. However, complementary information can also be extracted from studies at these wavelengths. For instance, optically thick emission imposes better constraints on the temperature and on the albedo of the dust particles, which improves the results from multiwavelength modeling of disks. Finally, observations at short wavelengths are the only way to explore the dust content in the gaps, whose emission is usually too optically thin to be detected at long wavelengths.

Recent multiwavelength analyses of Class II protoplane-tary disks (Carrasco-González et al. 2019; Macías et al. 2021; Sierra et al. 2021; Guidi et al. 2022) have measured dust properties within these substructures. The dust sizes determined from these observations are about 1 mm. However, there is a discrepancy with measurements obtained from polarization observations (Kataoka et al. 2017), which suggest dust sizes on the order of a few hundred micrometers. To reconcile these contradicting results, Kataoka et al. 2017 and Ohashi et al. 2020 have proposed the existence of two dust populations. The larger (millimeter) particles would experience more efficient vertical settling, become decoupled from the gas, and would gather toward the midplane of the disk. On the other hand, smaller particles would remain mixed with the gas in higher layers of the disk (Barrière-Fouchet et al. 2005). More recently, Zhang et al. 2023 proposed that large (larger than 1 mm) and porous (porosities higher than 70%) particles could reproduce both mul-tiwavelength and polarization results. A very recent polarization multiwavelength study strongly supported the presence of porous particles in HL Tau (Lin et al. 2024)

Pinte et al. (2016); Villenave et al. (2020, 2022); Doi & Kataoka (2021); and Pizzati et al. (2023) have focused on understanding the vertical structure of protoplanetary disks. Although limited in number, these studies suggest that Class II disks exhibit significantly low dust scale heights, particularly at the longer wavelengths observed with ALMA. This is a consequence of the change in the distribution of dust grains in the vertical direction, which in turn is due to settling (Dullemond & Dominik 2004; Pinilla et al. 2021). Most constraints presented in these papers remain upper limits (except for HD 163296), with suggestions such as H(100 au) < 1 au. Additionally, Pizzati et al. 2023 discovered values <4 au for certain disks from the Disk Substructures at High Angular Resolution Project (DSHARP), particularly those with less favorable orientations. These findings were predominantly derived from the presence and depth of gaps and rings in the outer regions of the disks (outer 100 au), where the scale height in the inner regions remains mostly unknown.

In this work, we present high-resolution Band 9 observations of HL Tau, a young stellar object in the Taurus-Auriga molecular cloud (~ 147 pc, Galli et al. 2018) with an estimated age younger than one million years (Liu et al. 2017). HL Tau is considered a Class I–II disk (Furlan et al. 2008), and it was the first proto-planetary disk in which dark and bright rings were discovered (ALMA Partnership 2015). It was intensively observed in the last years with ALMA and VLA, and therefore, the currently available high-quality data cover a wide range of wavelengths between 0.9 mm and 1 cm (Carrasco-González et al. 2016, 2019). Several previous studies also modeled its continuum emission to infer the dust properties (e.g., Jin et al. 2016; Pinte et al. 2016; Liu et al. 2017). These studies clearly established that the dark rings are gaps in the dust distribution, that is, regions with a density lower than in the adjacent bright rings. The origin of the substructures is still debated, but the most promising suggestions are planet-disk interactions (e.g., Jin et al. 2016; Dipierro et al. 2016). Our new sensitive and high-resolution observations at 0.45 mm are combined with previous multiwavelength observations to analyze the dust properties in the disk. The addition of shorter-wavelength observations allowed us to improve previous modeling of the dust properties in the disk and also to obtain new insight into the dust temperature very close to the protostar and the vertical structure of the disk.

2 Observations

We analyzed archival observations obtained with the Atacama Large Millimeter/submillimeter Array (ALMA) at a wavelength of 0.45 mm. One of the archival data sets that we used was part of the test observations for the band-to-band calibration during the High-Frequency Long Baseline Campaign (Asaki et al. 2020, project code: 2011.0.00005.E). The observations were taken during two execution blocks on November 3, 2017. The total observing time was 65 and 100 min for the first and second execution, respectively, which combined, gives a 45-min on-source time for HL Tau. The final data has eight spectral windows, each of which covers a bandwidth of 2 GHz. We also used additional data at a lower resolution (project code: 2017.1.01178.S, PI: E. Humphreys) to complete the uv coverage at shorter baselines. These data have eight spectral windows covering frequencies from 646–662 GHz with a bandwidth of 1.875 GHz each. The observations were taken from October 2, 2018, to October 18, 2018, with a total observing time on-source of 95 min.

We calibrated the data by using the calibration pipeline and scripts provided by the ALMA staff. We used version 5.6.1 of Common Astronomy Software Applications (CASA; McMullin et al. 2007) for self-calibration and cleaning in both data sets. We identified and flagged the spectral lines that were found in both data sets. We also averaged channels for the low- and high-resolution data; the final self-calibrated data sets contain spectral windows of eight channels of 200 MHz each.

We performed phase-only and amplitude self-calibration in both data sets. First, self-calibration was performed on the short-baseline data. During the self-calibration process, we created the images using the task CLEAN in CASA. The MTMFS decon-volver (Rau & Cornwell 2011) was used with scales of 0, 1, 3, and 5 times the beam size. We also used a Briggs weighting with a robust parameter of 0.5. For the low-resolution data set, the self-calibration process resulted in a clear improvement in the image quality. We applied six iterations of phase-only self-calibration in which we decreased the time-solution interval to 18 s. Furthermore, we were also able to apply two rounds of amplitude and phase self-calibration (calmode=ap) until no improvement in the signal-to-noise ratio was found. By the end of the whole self-calibration process, we achieved an increase in the S/N by a factor of ~17A We combined our short-baseline data together with the long-baseline data after aligning and correcting for proper motions and astrometric errors between our data sets. This was done by applying a Gaussian fit using CASA in the high-resolution image and using the tasks FIXVIS and FIXPLANETS to apply shifts to the low-resolution data. After this, self-calibration was attempted on the combined dataset. We started with phase-only self-calibration, but due to the low S/N at the long baselines, no improvement was found.

Last, self-calibration separately on the long-baseline data was also attempted by using the model obtained from the self-calibrated low-resolution data, but no improvement was found either. Since not enough solutions were found with any method to improve the quality of the combined data set, the self-calibration process was stopped.

With our final data set at hand, we recalculated the weights using the task STATWT in CASA. For the final images, we used the TCLEAN task in CASA with the MTMFS deconvolver (Rau & Cornwell 2011) using nterms = 1 and scales at 0, 1, 3, and 5 times the beam size. We used a pixel size that was ten times smaller than the beam size for all of the images. We explored different weightings: natural, uniform, and Briggs weighting, with different robust values. Finally, two images were used for this study. The first image was made by setting the parameter robust to 0 and using baselines shorter than 12 kλ in length. This is a good compromise between resolution and sensitivity and was used to study the large-scale emission of the disk. This image has an rms noise of 270 µJy beam⁻¹ and a beam size of ~23 mas. A second image was obtained by setting the parameter robust to −0.5, which resulted in a beam size of ~12 mas and an rms noise of 380 µJy beam⁻¹. The higher angular resolution of this image allowed us a detailed study of the innermost part of the disk (≲20 au). Both images are shown in Fig. 1.

In order to study the radial properties of the disk, we also obtained radial profiles of the emission and convolved the images to a round beam (robust 0, 24 mas, and robust −0.5, 13 mas). We used the known inclination and position angles of the disk for this, 46.72° and 138.02°, respectively (ALMA Partnership 2015). However, due to the inclination of the disk, there is a loss in angular resolution in the NE–SW direction. Since the S/N of the image is very high, we were able to avoid this loss in resolution by simply averaging the emission azimuthally within ±0.2 rad of the major axis of the disk. The radial profiles of the brightness temperature obtained from the intensity profiles of both images are shown in Fig. 2.

Fig. 1 HL Tau Band 9 images. Left panel: image using a robust parameter of 0 with a beam of 0.0226″ × 0.021″. Top right panel: zoom into the inner part of the HL Tau disk of the robust 0 images showing the asymmetry with gray contours at 11 times the rms and 33 times the rms. Bottom right panel: inner HL Tau disk from the robust −0.5 image with a beam of 0.0116″ × 0.010″.

3 Results

Figures 1 and 2 show that our images at 0.45 mm recover emission from the whole disk with a radius of ~ 1″ or ~ 140 au, similar to previous ALMA images at longer wavelengths (see, e.g., Carrasco-González et al. 2019; ALMA Partnership 2015). Our 0.45 mm image also shows some substructures, but as expected, substructures are detected with a lower contrast between dark and bright rings because of the higher optical depth at this wavelength. This suggests that the gaps contain small particles traced by the 0.45 mm image. These particles are more strongly coupled with the gas and might fill the gaps more than the larger particles that are observed at longer wavelengths (Pinilla et al. 2012b, de Juan Ovelar et al. 2013). Furthermore, we observed a slight deviation in the positions of the gaps and rings at 0.45 mm when compared to the radial profiles at other wavelengths (see Fig. 2), which suggests that there might be slight variations in the depth and width of the disk substructures at different wavelengths.

In the top right of Fig. 1, we show a zoomed-in view of the inner 50 au. An asymmetry is visible in the 0.45 mm data that is not seen at longer wavelengths: the NE part of the first ring is more intense than the SW part of the ring.

In the bottom right panel of Fig. 1, we show the image with the highest angular resolution, which allowed us to resolve the most internal part of the disk (≲20 au) for the first time. This image shows the first ring and a very deep gap.

This first gap, centered at 13 au, shows some emission in the low-resolution 0.45 mm, but appears very dark in the image with the highest angular resolution (see Figs. 1 and 2).

These observations provide a clearer view of the gap structure. They reveal that it is likely deeper and more depleted of dust particles than was initially observed with a lower resolution at 0.87 mm and 1.3 mm. Moreover, the very deep gap in the 0.45 mm image strongly suggests that it is indeed devoid of small dust particles as well.

The inner panel of Fig. 2 shows that the brightness temperature increases within 2.5 au at the highest angular resolution. This might be attributed to substructures within the inner disk that become resolved for the first time with the higher angular resolution of our 0.45 mm data.

In the following section, we first discuss the updated modeling of the radial dust properties using all the currently available data of HL Tau. We then describe the origin of the new features found in the 0.45 mm images, the asymmetry, and the substructure within the inner disk.

Fig. 2 Comparison of the HL Tau Band 6–7 (ALMA Partnership 2015) and the HL Tau Band 9 radial profiles of the brightness temperature. We plot the Band 6–7 ALMA image with the lowest resolution (33 mas, yellow), the Band 9 image with 22 mas (black), and the Band 9 image with 13 mas (green). The substructures are shallower in the radial brightness temperature for Band 9 than in the Band 6–7 image, and the emission from the Band 9 image extends up to 140 au. In the inset, we show the first 20 au of the HL Tau disk. We plot the brightness temperatures of the three images, and assuming that Tdisk = TB because of the optically thick regime, we display the position (in temperature) of the water snowline.

4 Discussion

4.1 Modeling of the dust properties

Several studies have modeled the dust properties in HL Tau using the available data (e.g., Pinte et al. 2016; Liu et al. 2017; Carrasco-González et al. 2019). We used a new modeling procedure that included the new 0.45 mm image we report here, which provides valuable constraints on the dust temperature and enhances the comprehensiveness of our multiwavelength analysis. Thus, this new analysis includes images at 0.45, 0.87, 1.3, 2.1, and 7.0 mm. We created a new 0.45 mm image with a robust parameter of 1.0 using baselines of <12 kλ, as previously described for the robust 0.0 image. Then, we convolved it with a 50 mas Gaussian beam and extracted radial profiles in wedges along the major axis within radial beams of approximately 0.01 arcsec. We also extracted radial profiles at the other wavelengths by applying the same procedure to the images presented in Carrasco-González et al. (2019). Subsequently, each point in the radial profiles was modeled independently. We note that the points in the model are not independent from the neighboring points due to the effects of the beam. Consequently, the contrast of the model parameters between rings and gaps is slightly smoothed in our results. We followed an approach similar to that of Carrasco-González et al. (2019), that is, at each radius, we fit the spectral energy distribution (SED) in the millimeter wavelength range using a 1D slab model given by $$I_v=B_{\nu }(T)\left[\left(1-\exp \left({\tau }_v/\mu \right)\right)+{\omega }_{\nu }F\left({\tau }_{\nu },{\omega }_v\right)\right],$$(1)

where $$\begin{array}{l}F\left({\tau }_v,{\omega }_v\right)=\frac{1}{\exp \left(-\sqrt{3}{ϵ}_v{\tau }_v\right)\left({ϵ}_v-1\right)-\left({ϵ}_v+1\right)}\\ \times \left[\frac{1-\exp \left(-\left(\sqrt{3}{ϵ}_v+1/\mu \right){\tau }_v\right)}{\sqrt{3}{ϵ}_v\mu +1}+\frac{\exp \left(-{\tau }_v/\mu \right)-\exp \left(-\sqrt{3}{ϵ}_{\nu }{\tau }_{\nu }\right)}{\sqrt{3}{ϵ}_{\nu }\mu -1}\right],\end{array}$$(2)

where ${\omega }_v=\frac{{\sigma }_v}{{\kappa }_v+{\sigma }_v}$ is the albedo, which is defined by the scat-tering coefficient (σ_v) and the absorption coefficient (κ_v), τ_v = Σ_dustχ_v, where χ_v = κ_v + σ_v, µ = cos(i) considers the inclination effects by correcting the optical depth by this inclination (i), and ${ϵ}_v=\sqrt{1-{\omega }_v}$. This model assumes that most of the dust content is settled in the midplane (see Sierra et al. 2019 and Carrasco-González et al. 2019 for details). While this assumption holds well for the outer HL Tau disk, it is not met within the initial ~32 au. A particle size distribution that follows a power-law distribution was also assumed, specifically, n(a) ∝ a^−p. We set the value of the exponent p to 3.5. This value is commonly used for the interstellar medium (ISM) (Mathis et al. 1977).

The dust opacity values were obtained assuming the DSHARP dust composition and following the effective medium approximation (Birnstiel et al. 2018). As usual in the study of protoplanetary disks, we first considered compact solid spherical particles. However, we also performed an analysis by considering spherical porous dust particles with a porosity of 90%. This is motivated by very recent results that strongly suggested that the HL Tau disk contains porous dust particles (Zhang et al. 2023 and Lin et al. 2024).

The opacity tables for the porous particles were obtained using the same procedure as in Birnstiel et al. (2018): the grains were assumed to be composed of small monomers within a void matrix. The optical properties of the monomers were first computed using the Bruggeman rule. This same rule was also used for the compact grains. The monomers were then mixed with a void matrix using the Maxwell-Garnett rule. It is well established that the fluxes, and consequently, dust masses significantly depend on the assumed dust composition, in particular, on the quantity of carbon. For instance, when the opacities presented by Ricci et al. (2010) and Stadler et al. (2022) are used, the fluxes they inferred were higher, which facilitated an explanation of the spectral indices of typical disks. This shows that the composition of the grains plays a crucial role, and the total dust mass derived in the subsequent section should be approached with caution.

We modeled the SED of the disk at each radius independently, and we fit the dust parameters (T_dust, a_max, and Σ_dust). We assumed a flux calibration uncertainty of 10% at 0.45 mm, 0.9 mm, 1.3 mm, and 7 mm, and 5% at 2.1 mm. We computed the posterior probability distribution of our model at each radius using the Markov chain Monte Carlo (MCMC) implementation within emcee (Foreman-Mackey et al. 2013). We assumed uniform priors in the dust surface density and maximum grain size. Following Macías et al. 2021, we used a conservative (i.e., wide) prior for the dust temperature based on the expected radial temperature profile for a passively irradiated disk, T_d = (ϕL_⋆/8πr²σ_SB)^0.25, assuming L_star = 6 ± 5 L_⊙ (i.e., a very conservative uncertainty on the stellar luminosity) and a range in ϕ (i.e., flaring angle) between 0.005 and 0.3. This prior is wide enough to have no effect on the posterior of the temperature at most radii. It only helps to avoid unphysically high dust temperatures at a few radii within the wide gap between 60 and 90 au, where the emission is optically thin and the dust temperature and dust surface density are highly degenerate.

Figures 3 and 4 show the dust parameters for the HL Tau disk obtained from the fitting with the solid and porous particles, respectively. Appendix A shows the radial intensity profiles of the model at each wavelength and their uncertainties, obtained from the median, 16th, and 86th percentiles of 500 random SEDs computed from the MCMC chains. As anticipated, the dust temperature is very well constrained for the compact and porous particles. We also note that the surface density is well constrained in both cases. Moreover, our findings reveal that the dust particles in the HL Tau disk are relatively large. Their sizes range from millimeters to centimeters. These sizes resemble well with models of dust evolution with traps located at different distances (Pinilla et al. 2015).

Our models find a sudden, and likely unphysical, increase in the particle size at the 13 au gap. However, Figs. A.1 and A.2 show that our models can’t reproduce the observed emission at these radii. The reason probably is our limited angular resolution, paired with the fact that the 13 au gap is highly depleted of large dust particles. The SED at high frequencies is therefore dominated by the small grains in the gap, while at longer wavelengths, the emission mostly comes from the beam smearing of the emission at the 20 au ring. Observations with a higher angular resolution at long wavelengths, which are only feasible with the upcoming ng VLA, are therefore necessary to accurately constrain the dust properties at these radii.

A distinct gap at 70 au, wide and almost without large particles, is evident in Figs. 3 and 4, consistent with prior observations by Carrasco-González et al. 2019. Interestingly, an increase in the dust surface density is then present at the 84 au and 96 au rings, with grain sizes that still remain at a few hundred microns. We note that an infalling streamer has been found to interact with the HL Tau disk at these radii (Garufi et al. 2022). We speculate that the small dust particles of this late-infalling material might be responsible for the creation of this dust ring, which would explain why the particles at these radii have not yet grown to millimeter/centimeter sizes. More detailed modeling of the interaction between the disk and infalling streamers is required to confirm this scenario.

Furthermore, clear distinctions emerge between the compact and porous particle solutions. In general, the model using compact particles exhibits lower temperatures and dust surface densities across the entire range of radii. The particle sizes predicted by the porous dust composition appear to be slightly smaller, although we note that they have a significantly wide uncertainty range.

Our new 0.45 mm data have allowed us to better constrain the radial dust temperature profile. We can therefore compare the position of the rings and gaps in HL Tau with the expected positions of the snowlines of some of the most relevant volatiles in protoplanetary disks more robustly. We plot in Figs. 3 and 4 the expected radii of the snowlines of H₂O, NH₃, CO₂, H₂S, and CO. Except for the H₂O snowline, the uncertainties in the sublimation temperatures of these volatiles now dominate the uncertainty in their positions. These large uncertainties still make it hard to extract any robust conclusion, despite our improved temperature profile. The abundance of substructures in HL Tau means that these snowlines will always fall close to a ring or gap. This makes it hard to confidently reject the hypothesis that snowlines cause at least a few of the ring substructures in HL Tau, as previously proposed (Zhang et al. 2015; Carrasco-González et al. 2019). However, many other substructures appear to be far from the expected snowlines, as in other protoplanetary disks (Bae et al. 2023). This mostly applies to the outer substructures, between 65 and 98 au, but also some between 36 and 45 au. We therefore conclude that mechanisms other than snowlines (e.g., MHD effects or planet-disk interactions) must cause most, if not all, the ring substructures in HL Tau.

We infer total dust masses of ${2.1}_{-0.84}^{+1.88}\text{ and }{6.3}_{-1.04}^{+1.57}$ M_J from the compact and porous particles, respectively. The dust masses when porous particles are assumed are thus higher by a factor of ~3, as previously reported by Zhang et al. (2023). Our measurement of the dust mass with compact grains is twice higher than the mass obtained in Carrasco-González et al. (2019) assuming compact particles with the same composition as ours (1.04 M_J). The primary difference between their analysis and ours is that we included Band 9 in our sample. These additional data have allowed us to more robustly constrain the dust temperature, which in turn affects the rest of the modeling. We obtained lower dust temperatures at all radii. The surface density is comparable at the inner radius, but it declines more rapidly than in their study. Additionally, our analysis reveals more pronounced substructures in the surface density. Notably, we observe dust particle sizes that are larger by approximately one order of magnitude within 60 AU, but beyond this point, the particle sizes decrease more significantly compared to the findings in Carrasco-González et al. (2019). These larger particles inside 60 au have lower absorption opacities, which explains that our dust mass is higher than that in Carrasco-González et al. (2019). Additionally, Booth & Ilee (2020) estimated a lower limit for the total gas mass of the disk of 209.5 M_J. This agrees well with our results of the compact dust particles when we consider a gas-to-dust ratio of 100. However, the elevated mass associated with porous dust particles from our analysis might imply a reduced gas-to-dust ratio of ~33 when compared to the total disk mass of 209.5 M_J. This suggests that the total gas mass of the disk might indeed still be higher than 209.5 M_J.

We note that our assumption of a 1D slab model might not be accurate in the inner radius <32 AU, where the emission at 0.45 mm is so optically thick that it might trace layers of the disk above the midplane (see Sect. 4.2). If all wavelengths effectively probe adjacent and close layers in the disk, indicating very little variation in temperature and dust size distribution between them, then the assumption could still be valid. The midplane of protoplanetary disks is generally vertically isothermal. It is, therefore, most likely that, if our observations trace significantly different layers of the disk, this only affects the maximum grain sizes traced at each wavelength. Sierra & Lizano (2020) explored these effects in more detail and found that the high optical depths in the innermost regions of disks could result in an incorrectly measured maximum grain size of a few hundred microns. While we cannot discard the possibility that our modeling is affected by these effects, the fact that we find maximum grain sizes of ~ 5 mm in the inner region instead of the 200–300 microns found by Sierra & Lizano (2020) suggests that our modeling is not severely affected by the high optical depths.

Fig. 3 Dust parameters at each radius using compact dust particle opacities. The positions of the most common snowlines are plotted in color. The vertical dotted and dashed lines are the positions of the gaps and rings, respectively, as seen from the surface density parameter. Top panel: dust temperature at each radius. Middle panel: dust surface density at each radius. Bottom panel: maximum particle size at each radius. The particle sizes are about 1 mm−1 cm.

4.2 Azimuthally asymmetric emission in the inner disk

As mentioned in Sect. 3, the first ring, centered at 21 au, shows an azimuthal asymmetry at 0.45 mm that becomes fainter at longer wavelengths. This might be contingent upon the resolution and sensitivity of the observations. This asymmetry appears as an increase in brightness at the northeastern side of the ring, covering almost 180 degrees and centered on the minor axis of the disk (see Fig. 1). The position of this asymmetry strongly suggests that it is caused by the combination of the disk inclination and its high optical depth at the observed wavelength (see Fig. 5). Specifically, in the NE part of the disk, we directly observe the emission from the internal wall of the first ring, which is directly irradiated by the central star, and this direct irradiation increases its temperature. On the other hand, in the SW part of the disk there is no direct line of sight toward the wall. The wall emission is therefore blocked by the optically thick material in the line of sight (see Fig. 5).

The HL Tau outflow is blueshifted in the NE and redshifted in the SW (Lumbreras & Zapata 2014). We can therefore infer that the NE corresponds to the far side of the disk and the SW corresponds to the close side of the disk. This explains the position of the asymmetry in that side of the minor axis of the disk.

The asymmetry is in this case fundamentally similar to the asymmetries found in Class 0/I objects (e.g., Guerra-Alvarado et al. 2024; Lin et al. 2023; Ohashi et al. 2023), where the dust particles did not have time to settle or are affected by significant turbulence that prevents the grains from falling onto the midplane (larger dust scale height), so that they remain optically thick even at longer wavelengths with ALMA. On the other hand, this asymmetry is fundamentally different from the large-scale asymmetries observed in other disks, where the effects become more pronounced at longer wavelengths (e.g., IRS48, van der Marel et al. 2015; and HD 142527, Casassus et al. 2015) which contradicts the observed behavior in HL Tau. These asymmetries are commonly associated with dust traps in which dust can grow to centimeter sizes.

In order to confirm that the asymmetry is caused by a geometric effect, we used a toy model to reproduce the observed asymmetry. We used RADMC-3D (Dullemond et al. 2012) to perform radiative transfer models at 0.45 mm while varying parameters such as the disk dust mass and scale height. We then compared the model emission with the observed data to gain insight into the physical properties and vertical structure of the disk. In this specific radiative modeling, a generic protoplane-tary disk model was used within RADMC-3D. The disk dust masses and scale height were adjusted to explore different scenarios and to reproduce the observed asymmetry. The optool package (Dominik et al. 2021) was employed to compute the DSHARP dust particle opacities using a_min as 0.050 µm and a_max as 3 mm. For this toy model, we were only interested in the inner regions of the disk. We therefore only considered emission up to a radius of 32 au (the location of the second gap). We assumed that all the grains are at the same height, and we set a small inner radius R_in = 1 au. A gap between 8 au and 17 au was included in the model to capture the effects of the asymmetry at the outer edge of the gap.

The generic protoplanetary disk model uses a density distribution, $$\rho (r,z)=\frac{\text{Σ}(r)}{H_{\text{p}}\sqrt{2\pi }}\exp \left(-\frac{z^2}{2H_{\text{p}}^2}\right),$$(3)

where r is the distance to the star from the disk, H_p (r) is the dust scale height of the disk, and Σ(r), is the dust surface density, defined as $$\text{Σ}(r)={\text{Σ}}_0{\left(\frac{r}{r_{\mathrm{out}}}\right)}^{-1}.$$(4)

The dust scale height in the radiative modeling is defined as a power-law-like dependence, given by $$H_{\text{p}}(r)=H_{100}{\left(\frac{r}{100\text{\hspace{0.17em}AU}}\right)}^{1+\Psi },$$(5)

where Ψ is the value of the flaring index set by default to 0.14, and H₁₀₀, is the value of the scale height at a distance of 100 AU from the central star.

Finally, a total stellar mass of 1.7 M_⊙ (Pinte et al. 2016), a stellar luminosity of 6 L_⊙, and an effective temperature of 4395 K (White & Hillenbrand 2004) were assumed. In addition, the inclination angle for HL Tau was set to 46.72°, and the distance was 147 pc.

Several models were run. We varied the dust scale height and the dust disk mass within 32 au $\left(M_{\text{Disk }}^{<32\text{au}}\right)$ based on previous studies (Pinte et al. 2016; Liu et al. 2017). We show these parameters and their combination in Table 1.

The model images were then convolved using a Gaussian beam with the same beam as for the 0.45 mm image. The final convolved model images that seem to reproduce the emission of the observation better are shown in Fig. 6. We also extracted the deprojected emission from each pixel in a ring ranging from 0.12″ to 0.15″ and around 360° to capture the increased emission in the asymmetry. We show the normalized azimuthal profiles from the observation and all models in Fig. 6. We note that our goal with this toy models was to show that this type of asymmetry can be obtained at 0.45 mm, and that the dust mass and scale height can be constrained from these images. We therefore did not try to reproduce the exact flux of our observations because that would require a much more detailed modeling in which the disk parameters would need to be varied more strongly. The azimuthal profiles for 1, 3, and 5 Jupiter masses (inner disk mass + ring mass) are also shown in Fig. 6.

Our modeling shows that an asymmetry such as that observed in HL Tau can be produced by the combined effects of disk orientation, high optical depth, and larger scale heights at 0.45 mm. Figure 6 shows that the dust mass plays a significant role in reproducing the observed asymmetry at 0.45 mm. In models in which the scale height is not large enough, the emission no longer produces an asymmetry.

If the asymmetry arises from the scenario we outlined before, we find that a disk dust mass (within 32 au) of approximately 1 M_J, combined with a scale height H/R of 0.08, a disk dust mass of 3 M_J with a scale height H/R of 0.05, or a disk dust mass of 5 M_J with a scale height H/R slightly smaller than 0.05 reproduces the observed asymmetry in the observations best. In other words, the higher the dust mass, the smaller the dust scale height required to reproduce an asymmetry. Because of this, an independent estimate of the dust mass within 32 au would allow for a better constraint on the scale height of the disk. We employed the mass obtained from our multiwavelength analysis in Sect. 4.1, acknowledging its known limitations, for an initial approximation of the scale height in this specific region of the disk. From this analysis, we obtained a disk dust mass inside the first 32 au of ${0.86}_{-0.368}^{+1}$ M_J. If our proposed scenario is the cause of the asymmetry, our derived mass would then imply a scale height ~≥0.08 H/R in that specific region of the disk. As Fig. 6 shows, to determine the scale height more effectively, we need a more restrictive constraint on the disk dust mass in subsequent studies, together with a model that reproduces the flux and the actual structure of the HL Tau disk.

We note that the scale height necessary to reproduce the observed asymmetry is significantly larger than the estimates from radiative transfer modeling of the 0.9–3 mm emission, H/R < 0.01 (Pinte et al. 2016). One possible explanation for this inconsistency is that in Pinte et al. (2016), the focus was primarily on gaps and rings in the outer disk, rather than on those within 30 au as constrained in this paper. This result then suggests that the HL Tau disk may exhibit a greater vertical extent (thickness) in the inner regions while maintaining a flatter profile in the outer regions. Radial variations in scale height have indeed been identified for the two rings of HD 163296, with the 67 au ring being very thick (>4 au) and the 100 au ring being very flat (<1 au) (Doi & Kataoka 2021; Liu et al. 2022).

While there are limited constraints on the radial variation in the scale height, this outcome implies that turbulence in disks might be higher in the inner than in the outer regions. This aligns well with the high accretion rates on some disks that cannot be explained by their current low turbulence values. In this context, the magnetically driven accretion disk models suggested by Delage et al. (2022) could explain the heightened turbulence parameter and increased scale height, leading to the observed asymmetry beyond the dead zone (at 23 AU) and within the MRI-active layer in the HL Tau disk. Additionally, this is consistent with what is expected in the case of vertical shear instability with dust coagulation, where the transition between the thick and thin disk occurs at around 100 au (Pfeil et al. 2023). Finally, it also remains inconclusive whether the streamer observed by Garufi et al. (2022) in HL Tau could have an effect at smaller radii through the flow of material into the inner disk. The findings in Jiang et al. (2024) suggest a potential influence of this streamer on the turbulence at the outer radius (70 au), raising questions about a possible impact on smaller scales. Further analysis is essential to discern the factors contributing to the increased turbulence, if any, in the HL Tau inner disk and to unravel the specifics of this phenomenon.

Fig. 4 Same as Fig. 5, but using a 90% porosity for the dust particles and for the dust opacities. The dust temperature is lower and the surface density is higher around the gaps, and the maximum particle sizes are similar, but at the 13 au gap, the particle size increases strongly.

Fig. 5 Scenario proposed for the origin of the asymmetry seen in HL Tau at 0.45 mm.

Fig. 6 RADMC-3D models of the inner HL Tau disk compared with the Band 9 observation. Top panel: convolved images of the nine different models with RADMC-3D that reproduce the observation better, corresponding to masses of 1 MJ, 3 MJ, and 5 MJ, paired with scale height values of 0.08, 0.05, and 0.03 Hp. The disk mass is defined as the dust mass within 32 au. Bottom panels: normalized azimuthal profiles considering the emission between 0.12 and 0.15 arcseconds in the HL Tau Band 9 observation and all the RADMC-3D models. The lowest panels depict the difference, measured in sigmas, between the observation and the respective models.

4.3 Substructure in the inner disk: Possible sign of the traffic-jam effect

As mentioned, our image with the highest angular resolution at 0.45 mm revealed an increase in the brightness temperature within ~2.5 au (see the inset in Fig. 2). There are several possible explanations for this.

First, we discard the hypothesis that this increase is due to free-free contamination. It is well known that HL Tau drives a powerful jet whose emission has been previously modeled using high angular resolution VLA data at long wavelengths (Carrasco-González et al. 2019). Thus, we know that in the 2343 GHz frequency range, the expected free-free emission from the radio jet is given by $$\left[\frac{F_v}{\mu {\text{Jy beam }}^{-1}}\right]=175\times {\left[\frac{v}{32.5\text{\hspace{0.17em}GHz}}\right]}^{0.7},$$(6)

where v is the frequency of the image at which we wish to calculate the free-free emission. In this frequency range, the free-free emission was found to be partially optically thin, which is reflected in the 0.7 spectral index. Because the optical depth of the free-free emission is proportional to v^−2.1 (e.g., Reynolds 1986), we expect at the frequency of our 0.45 mm image, that is, ~650 GHz, that free-free emission from the jet is optically thin. Thus, there must be a turning point between 50 and 650 GHz. A conservative assumption is that emission from the jet becomes totally optically thin at 100 GHz. Then, using Eq. (6) to estimate the free-free emission at 100 GHz, and extrapolating this with a spectral index of −0.1 to higher frequencies, we obtain that an estimate of the free-free contamination at 650 GHz would be ~300 µJy, which is ~2% of the peak emission at the center of the disk in the 0.45 mm image. The actual contribution at this frequency is probably lower since it is very likely that the free-free emission would become optically thin at a frequency lower than 100 GHz. Another potential explanation for a significant free-free emission at the central part of the disk is the possibility that it becomes highly optically thick due to photoionization. However, this scenario faces challenges because HL Tau is a K-type star, and these stars typically do not emit sufficient ionizing photons to fully ionize a very dense gas. All of these factors lead us to the conclusion that the elevated emission in the inner region of the disk is unlikely to be primarily caused by free-free contributions.

An interesting explanation for the increased dust emission in the inner part of the disk is the possibility of a dust particle pile-up in that region. Our unprecedented angular resolution allows us to probe the dust emissions down a few au, inside the estimated water snowline at ~5 au when we assume T_disk = T_B (Zhang et al. 2015, Fig. 2). Inside the water snowline, theoretical models predict the appearance of a traffic-jam effect (e.g., Banzatti et al. 2015; Pinilla et al. 2017), as dry silicates are more sensitive to fragmentation than ice-rich grains (Dominik & Tielens 1997; Supulver et al. 1997; Wada et al. 2013; Gundlach & Blum 2014). This leads to lower particle sizes and to a reduced drift speed. The traffic-jam effect offers a favored pathway for the formation of dry planetesimals in the inner disk, but this was recently questioned in the light of new experiments to determine the stickiness of ice-rich grains (e.g., Gundlach et al. 2018; Musiolik & Wurm 2019) and in light of observations that showed the resilience of icy pebbles to sublimation (Houge et al. 2024). Interestingly, Facchini et al. (2024) recently detected water emission at 321 GHz from the inner disk of HL Tau, but their spatial resolution was not enough to robustly measure the position of the water snowline.

Our results could be consistent with the traffic-jam effect because a pile-up of solids like this would increase the dust density and optical depth inside the water snowline (e.g., Banzatti et al. 2015), causing observations to probe higher layers of the disk with more elevated dust temperatures, thereby manifesting a higher brightness temperature, as shown in Fig. 2 (see Drazkowska & Alibert 2017; Ros et al. 2019). Moreover, an increase in the population of smaller particles will also result in a decrease in the albedo, which in turn results in an increase in the emission.

A second possibility is that the increase in the brightness temperature is attributed to viscous heating in the inner disk. Overall, the heating mechanisms in protoplanetary disks significantly influence the temperature profile of these. The two dominant heating processes in these disks are stellar radiation and viscous heating (D’Alessio et al. 2005). Throughout most of the disk, stellar radiation contributes most to the thermal structure of the disk. However, in the innermost regions of the disk close to the midplane, where dust surface densities reach notably high values, it is highly likely that viscous heating becomes the primary contributor to the temperature profile of the disk (D’Alessio et al. 2001). More recently, Chambers (2009) found that the midplane temperature may be affected by viscous heating within the radial range of 1–40 au in protoplanetary disks. This phenomenon also influences the position of the snowline, which varies depending on the age of the disk. The position where viscous heating takes place, however, depends upon the characteristics of the star and on many other parameters that are hard to constrain.

It is currently challenging to distinguish between the traffic-jam and the viscous heating scenarios. We strongly encourage future works to focus on this aspect, however, given the implications it could have concerning planet formation in the inner disk. For example, observations with a high angular resolution at other wavelengths could help us to constrain the dust emission spectral index and particle size, which would help us to understand the underlying mechanisms that drive the observed changes in brightness temperature.

5 Summary and conclusions

We reported new sensitive ALMA images at 0.45 mm of the HL Tau disk with an unprecedentedly high angular resolution of 12 mas or 1.7 au. At this wavelength, the detected emission is known to be highly optically thick. This is evidenced by the shallower substructures that we observed. The extended nature of the disk suggests that we trace smaller dust particles that are more closely coupled to the gas.

By combining the Band 9 image with the images at other wavelengths, we updated the model of the dust properties in the disk. The 0.45 mm image has allowed us to better constrain the dust temperature of the disk, which provided better constraints on the dust density and particle size. The particle sizes are still predicted to be larger than millimeter sizes until ~60 au, but there is a significant decrease (to ~200 µm) between 60 and 100 au. These smaller particles coincide with the region in which an accretion shock of infalling material has been found, traced by SO and SO₂. The final derived disk dust masses are 2.1 M_J for compact particles and 6.3 M_J for porous dust particles. These dust masses translate into gas-to-dust ratios ≤100.

Our 0.45 mm data clearly show an asymmetry in the emission from the first ring. The emission is considerably more intense in the NE part of the disk than in the SW part. We interpreted this as a combination of very optically thick emission, a moderately inclined disk, and a large dust scale height. In the NE part of the disk, we directly view the emission from the internal wall of the first ring, which is illuminated by the radiation from the central star. This interpretation is supported by our radiative transfer model, which demonstrates the dependence of this asymmetry on the dust mass of the inner regions of the disk and the dust scale height. We propose that this type of asymmetry in other disks can be used to constrain the dust scale height if an independent estimation of the dust mass of the disk has been made. Based on our proposed scenario and the multiwavelength analysis model, we implied the possibility of a dust scale height H/R ∽≥ 0.08 for the HL Tau disk at least in the inner 32 au. This approximation appears to be in strong contrast with the constraints identified by Pinte et al. (2016), who determined H/R < 0.01 in the outer disk. We note that this discrepancy might stem from higher turbulence in the inner regions of HL Tau compared to its outer region, leading to radial variations in the vertical structure of the HL Tau disk. It may be necessary to incorporate the vertical structure in the multiwave-length analysis to determine the disk mass and to employ more quantitative modeling approaches, instead of the qualitative ones in RADMC-3D, to obtain a more realistic scale height of the disk.

We have resolved the innermost part of the HL Tau disk and showed that there is some previously unseen substructure. It is revealed as a steep increase in the brightness temperature very close to the central star, at radii smaller than 5 au. Interestingly, the radius at which this increase starts appears to coincide with the position of the water snowline. Therefore, we discussed the possibility of a pile-up of small particles due to a decrease in the drift velocity of the small grains after they cross the water snowline and are fragmented more efficiently. We also discussed a different possibility in which we resolve the region in which the heating of the disk changes from being dominated by irradiation to being dominated by viscous heating. Previous theoretical models predicted that this change should take place very close to the protostar. To further confirm either of these scenarios, additional modeling and observations at other wavelengths are necessary for the HL Tau disk.

Acknowledgements

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00005.E and ADS/JAO.ALMA#2017.1.01178.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. C.C.-G. acknowledges support from UNAM DGAPA-PAPIIT grants IG101321 and IG101224 and from CONACyT Ciencia de Frontera project ID 86372. The research of M.V. was supported by an appointment to the NASA Postdoctoral Program at the NASA Jet Propulsion Laboratory, administered by Oak Ridge Associated Universities, under contract with NASA.

Appendix A Markov chain Monte Carlo model versus observations: VLA fitting

Fig. A.1 Comparison of the brightness temperatures of the model (dash lines) with observations (solid lines) at each radius and wavelength. Compact dust particles were used, and the VLA emission is plotted and fit within the MCMC. The vertical lines are the positions of the gaps and bright rings taken from the dust parameters, and the size of the beam is plotted as a Gaussian in the top right corner for the five wavelengths (50 mas).

Fig. A.2 Same as Figure B.1, but with the porous particle model instead.

References

All Tables

All Figures

	Fig. 1 HL Tau Band 9 images. Left panel: image using a robust parameter of 0 with a beam of 0.0226″ × 0.021″. Top right panel: zoom into the inner part of the HL Tau disk of the robust 0 images showing the asymmetry with gray contours at 11 times the rms and 33 times the rms. Bottom right panel: inner HL Tau disk from the robust −0.5 image with a beam of 0.0116″ × 0.010″.
In the text

Fig. 2 Comparison of the HL Tau Band 6–7 (ALMA Partnership 2015) and the HL Tau Band 9 radial profiles of the brightness temperature. We plot the Band 6–7 ALMA image with the lowest resolution (33 mas, yellow), the Band 9 image with 22 mas (black), and the Band 9 image with 13 mas (green). The substructures are shallower in the radial brightness temperature for Band 9 than in the Band 6–7 image, and the emission from the Band 9 image extends up to 140 au. In the inset, we show the first 20 au of the HL Tau disk. We plot the brightness temperatures of the three images, and assuming that Tdisk = TB because of the optically thick regime, we display the position (in temperature) of the water snowline.

In the text

Fig. 3 Dust parameters at each radius using compact dust particle opacities. The positions of the most common snowlines are plotted in color. The vertical dotted and dashed lines are the positions of the gaps and rings, respectively, as seen from the surface density parameter. Top panel: dust temperature at each radius. Middle panel: dust surface density at each radius. Bottom panel: maximum particle size at each radius. The particle sizes are about 1 mm−1 cm.

In the text

	Fig. 4 Same as Fig. 5, but using a 90% porosity for the dust particles and for the dust opacities. The dust temperature is lower and the surface density is higher around the gaps, and the maximum particle sizes are similar, but at the 13 au gap, the particle size increases strongly.
In the text

	Fig. 5 Scenario proposed for the origin of the asymmetry seen in HL Tau at 0.45 mm.
In the text

Fig. 6 RADMC-3D models of the inner HL Tau disk compared with the Band 9 observation. Top panel: convolved images of the nine different models with RADMC-3D that reproduce the observation better, corresponding to masses of 1 MJ, 3 MJ, and 5 MJ, paired with scale height values of 0.08, 0.05, and 0.03 Hp. The disk mass is defined as the dust mass within 32 au. Bottom panels: normalized azimuthal profiles considering the emission between 0.12 and 0.15 arcseconds in the HL Tau Band 9 observation and all the RADMC-3D models. The lowest panels depict the difference, measured in sigmas, between the observation and the respective models.

In the text

Fig. A.1 Comparison of the brightness temperatures of the model (dash lines) with observations (solid lines) at each radius and wavelength. Compact dust particles were used, and the VLA emission is plotted and fit within the MCMC. The vertical lines are the positions of the gaps and bright rings taken from the dust parameters, and the size of the beam is plotted as a Gaussian in the top right corner for the five wavelengths (50 mas).

In the text

	Fig. A.2 Same as Figure B.1, but with the porous particle model instead.
In the text