© M. Villenave et al. 2020

1 Introduction

The process of planet formation requires small submicron-sized particles to grow up to large pebbles and boulders that will eventually aggregate to form planetesimals and planets. Given the short lifetimes of disks, such efficient growth has to occur in the highest density regions of protoplanetary disks, that is to say, the inner regions and/or the disk midplane. Gas drag, the interaction of dust (in Keplerian rotation) with the gas orbiting around the central star at slower (sub-Keplerian) velocities, is at the origin of the vertical settling to the midplane and of the inward radial drift of large (e.g., millimeter-sized) dust grains (Weidenschilling 1977; Barrière-Fouchet et al. 2005). Unlike the larger grains, micron-sized particles are well coupled to the gas and are located in similar regions, well-mixed with the gas. The characteristic timescale of radial drift is predicted to be about a hundred times longer than that of vertical settling (Laibe et al. 2014). However, the strength of these effects is not yet well constrained and depends in particularon the disk viscosity and/or turbulence, on the gas-to-dust ratio, and on the initial grain size distribution (Fromang & Papaloizou 2006; Dullemond & Dominik 2004; Mulders & Dominik 2012). The comparison of observations at widely different wavelengths, for example optical-near-infrared (NIR) and (sub)millimeter, can help to constrain the settling intensity and radial drift of dust grains. Moreover, by performing multi-wavelength observations in the millimeter, one can achieve spectral index measurements and, assuming optically thin emission, interpret the findings in terms of grain growth (e.g., Guilloteau et al. 2011; Pérez et al. 2012).

Up to now, most studies of protoplanetary disks have focused on radial features, such as rings, gaps, and spirals, which naturally led to observations of low inclination systems (≤70°), where they are more readily visible. Some of these studies showed that the gas is often more radially extended than the millimeter-sized dust (Ansdell et al. 2018; Facchini et al. 2019), which is likely a combined effect of optical depth and dust radial drift (Facchini et al. 2017). However, though it is important to constrain radial drift, dedicated comparisonsof the radial distribution of different dust grains sizes remain sparse (see e.g., Pinilla et al. 2015; Long et al. 2018).

For relatively face-on disks, it is difficult to estimate the difference in vertical extent between gas and dust grains. Detailed modeling of ring and gap features is required (e.g., Pinte et al. 2016). On the other hand, edge-on disks offer a unique perspective, as they allow us to directly observe their vertical structure. Previous studies of edge-on disks at different wavelengths indicate that the grains are stratified, with larger dust closer to the disk midplane (Glauser et al. 2008; Duchêne et al. 2003, 2010; Villenave et al. 2019), as predicted by models. However, even for these very inclined systems, the vertical extent of the midplane remains poorly constrained in early studies, limited by the resolution of the observations.

In this work, we present a survey of edge-on disks observed with the Atacama Large Millimeter Array (ALMA) at high angular resolution (~0.1′′). The sample was selected based on Hubble Space Telescope (HST) images and most of the targets were observed at multiple millimeter bands. After describing the sample and the data reduction in Sect. 2, we present the fluxes, brightness temperatures, spectral indices, and the radial and vertical extents of all disks in Sect. 3. In Sect. 4, we compare our results with a toy model. We focus on the vertical and radial extent of the disks, and investigate optical depth in the disks by studying the brightness temperature and spectral indices of the sources. Finally, we summarize our conclusions in Sect. 5.

2 Observations and data reduction

2.1 Target selection

In this study, we selected a sample of twelve highly inclined disks. All sources were identified as candidates from their spectral energy distribution (SED, see Stapelfeldt et al. 2014, and in prep.) and confirmed as edge-on disks (EOD) according to their opticalor NIR HST scattered-light images. At these wavelengths, edge-on disks are inclined enough so that the star is not in direct view of the observer. The ALMA observations targeted 8 sources located in the Taurus star-forming region, 3 in Chamaeleon I, and 1 in Ophiuchus. Most of the observations presented in this work were part of our cycle 4 and 5 survey of edge-on disks (Project 2016.1.00460.S, PI: Ménard), but we also included complementary observations of Tau 042021, HH 30, and Oph 163131 from previous programs (Projects 2013.1.01175.S, 2016.1.01505.S, and 2016.1.00771.S, PIs: Dougados, Louvet, and Duchêne).

We report the coordinates, spectral types, and stellar masses of the sources in Table 1. As the spectral features associated with the central source are still visible for edge-on disks through the scattered light (Appenzeller et al. 2005), the spectral types could be determined from spectroscopy (Luhman 2007; Luhman et al. 2010, Flores et al., in prep.). However, the edge-on configuration does not allow a direct estimate of the stellar luminosity. Thus, we report dynamical stellar masses estimated from CO emission.

At optical and NIR wavelengths, edge-on disk systems highlight extended nebulosities on both sides of a dark lane, tracing the disk midplane. Because of severe extinction in the dark lane, the central star is usually undetected at optical-NIR wavelengths, also resulting in fainter systems for a given spectral type and distance. In a few cases, the brightness distribution of the nebulosities has also been observed to vary significantly (Stapelfeldt et al. 1999). These facts combine to render parallax measurements difficult and the distances determined by Gaia can be plagued by large errors. For our targets, we decided to adopt the average distance of the parent star-forming regions instead; 140 pc for the sources in Taurus and Ophiuchus (Kenyon et al. 2008; Ortiz-León et al. 2018; Cánovas et al. 2019) and 192 pc for those in Chamaeleon I (Dzib et al. 2018).

Four of the observed disks are part of multiple systems: HK Tau B, HV Tau C, Haro 6-5B, and HH 48 NE (Stapelfeldt et al. 1998, 2003; Krist et al. 1998; Haisch et al. 2004) with apparent companion separations larger than 2.4′′. While HH 30 has been suggested to be a binary on the basis of jet wiggles and a disk central hole in lower resolution and signal-to-noise ratio millimeter continuum maps (Guilloteau et al. 2008), higher resolution ALMA observations do not confirm the central hole (Louvet et al. 2018). Neither do they exclude the possibility of spectroscopic binary. The other targets in the sample are not known to be in multiple systems.

2.2 ALMA edge-on survey

Our ALMA cycle 4 and 5 observations (Project 2016.1.00460.S, PI: Ménard) were divided into two groups, targeting respectively seven sources located in Taurus and three in Chameleon I. The observations were performed in band 7 (0.89 mm) and in band 4 (2.06 mm). We present the different setups in the following section.

2.2.1 Band 7 survey

The band 7 observations of the Taurus sources (Tau 042021, IRAS 04302, HK Tau B, HV Tau C, IRAS 04200, Haro 6-5B, and IRAS 04158) were performed with both a compact and an extended array configuration. For the Chamaeleon sources (ESO-Hα 569, ESO-Hα 574, and HH 48 NE), only the compact configuration was observed. The observational setup is summarized in Table 2. The spectral setup was divided into three continuum spectral windows, with dual polarization, 2 GHz bandwidth spectral windows centered at 344.5, 334.0, and 332.0 GHz, and one spectral window set to observe the ¹²CO J = 3–2 transition at 345.796 GHz. In this paper, we focus on the continuum data which has a geometric mean frequency of 336.8 GHz (0.89 mm). The reduction and analysis of the CO data will be presented in a separate study. Because the observationswere performed over two cycles, different versions of CASA have been used for the calibration.

We calibrated the raw data of the compact array executions using the pipeline from CASA package version 4.7. The raw data of the extended configuration observations of the Taurus sources was manually calibrated, using CASA version 5.1. Whenever possible we used the supplied ALMA phase calibrator in the dataset (for Tau 042021, Haro 6-5B, IRAS 04158). However, for four targets (namely HK Tau B, HV Tau C, IRAS 04200, and IRAS 04302), the spectral window setup used for the supplied phase calibrator was different from that of the science target. These data could not be reduced using the standard pipeline. For HV Tau C and IRAS 04200, it was possible to use the check source as a phase calibrator, and for HK Tau B, the more distant bandpass calibrator could be used. In these cases, the calibrator was only observed once before each of the science targets (rather than bracketing it in time); this increased the phase calibration uncertainties for these objects. IRAS 04302 did not have a usable phase calibrator. However, this source is bright and extended, and self-calibration using the average of all spectral windows could be performed without an initial phase reference. Consequently for this target, there was no absolute astrometric solution in the extended configuration data; also the achievable angular resolution was somewhat worse than the other sources of the sample.

2.2.2 Band 4 survey

The band 4 observations of the edge-on survey were only performed for the sources located in Taurus (Tau 042021, IRAS 04302, HK Tau B, HV Tau C, IRAS 04200, Haro 6-5B, and IRAS 04158). The continuum spectral windows were centered on 138, 140, 150, and 152 GHz, with a geometric mean frequency of 145.0 GHz (2.06 mm). The relevant parameters of the observations are reported in Table 2. The raw data were pipeline calibrated using the CASA package, version 5.1.

2.3 Archival and literature data

2.3.1 ALMA observations of HH 30

We include multi-wavelength, band 4, 6, and 7, observations of HH 30 in the present study. The observational setup and data reduction of the band 6 observations are presented in Louvet et al. (2018, Project 2013.1.01175.S, PI: Dougados). We also use new band 4 and band 7 cycle 4 observations (Project 2016.1.01505.S, PI: Louvet), and present the data reduction in the following paragraphs.

The band 7 observations of HH 30 were performed with only one array configuration, in four executions between October 14 and October 21, 2016 (see Table 2). The dual polarization spectral setup included two 2 GHz bandwidth spectral windows for the continuum emission, centered at 331.6 and 344.8 GHz, and two other spectral windows set to observe the ¹³CO and C¹⁸ O J = 3–2 emission. Here we present only the continuum observations. The observations were calibrated by the ALMA pipeline using CASA 4.7.

The band 4 observations were performed with two configurations: a compact configuration and an extended configuration (see Table 2). The ALMA correlator was configured to record dual polarization with four separate continuum spectral windows of 2 GHz each centered at 138, 140, 150, and 152 GHz. The observations were calibrated by the pipeline using CASA 4.7.

2.3.2 ALMA observations of Tau 042021 and Oph 163131

We also include band 6 cycle 4 observations of Tau 042021 and Oph 163131 (Project 2016.1.00771.S, PI: Duchêne). Although the spectral setup targeted emission lines of three CO isotopologues, we focus here on the continuum data. An analysis of the line emission of Oph 163131 will be presented in a separate study (Flores et al., in prep.).

For both sources, the two continuum spectral windows were centered on 216.5 and 232.3 GHz leading to a geometric mean frequency of 224.4 GHz (1.34 mm). The detailed observational setup is presented in Table 2. Data from both configurations were reduced using the pipeline from CASA package version 4.7.

2.4 ALMA imaging

We constructed the images from the calibrated visibilities with a Briggs robust weighting parameter of 0.5, except for the band 6 of HH 30 for which we used the reduction presented in Louvet et al. (2018, robust of 0.56). When the data were taken with several array configurations (see Table 2), the images were generated by combining the visibilities from the compact and extended configurations. The only exception is IRAS 04158 for which we show the band 7 compact configuration observations only. We applied self-calibration on the sources with the highest signal-to-noise ratio (S/N) (on at least one array configuration). In band 7 this corresponds to Tau 042021, IRAS 04302, HK Tau B, and Haro 6-5B, in band 6 to Tau 042021 and Oph 163131, and to only HH 30 in band 4. For the Taurus survey, this leads to typical beam sizes of about 0.11′′ × 0.07′′ in band 7 and about 0.11′′ × 0.04′′ in band 4. For the other sources, the beam sizes range from 0.23 × 0.13′′ to 0.48 × 0.30′′. The beam sizes of each observations are reported in Table A.1, and we present the images in Fig. 1.

Additionally, we also recompute the maps to get a unique angular resolution for each source observed both in band 4 and band 7 (except for IRAS 04158, for which the disk is not detected in band 4). To do so, we first re-imaged each source using the same uvrange and applying a uvtaper to limit the effect of flux filtering and of different uvcoverage. Then, to obtain exactly the same angular resolution between both bands, we used the imsmooth CASA task. The restored beam sizes are reported in the last column of Table A.1.

2.5 Astrometric accuracy and map registration

ALMA maps are subject to astrometric uncertainties due to limited signal-to-noise on the phase calibrator and errors in the phase referenced observations related to weather or antenna position errors. To compute accurate spectral index maps (see Sect. 3.3), the images at the different wavelengths have to be accurately aligned. In this section, we discuss the registration of our images.

For each source and band, we performed simple Gaussian fits in the image plane to estimate the centroid position. The offsets between band 4 and band 7 ranged from 2 mas (milliarcseconds) for Tau 042021 and IRAS 04200, up to 60–90 mas for HH 30, HV Tau C, and IRAS 04302. Most are larger than expected from source proper motion between the different executions and from the astrometric accuracy of ALMA (10% of the synthesized beam, see ALMA Technical Handbook). However, for the three disks with the largest offsets (HH 30, HV Tau C, and IRAS 04302), different phase calibrators were used between the observations. In addition, the standard phase calibration could not be performed for the band 7 data of IRAS 04302 and HV Tau C (see above), and the weather conditions during the band 7 extended configuration observations were relatively poor, with phase rms of ~ 64° on the longest baselines. This will also deteriorate the positional accuracy. Overall, the astrometric accuracy of our observations is not sufficient to confirm any significant motion or shift of the sources between the two bands. Because models of edge-on disks also show that offsets between bands should remain minimal compared to the beam size, we registered our images to a common center in both bands. We aligned the band 4 and band 7 observations using the fixplanets task in CASA¹.

3 Results

3.1 Continuum emission and brightness temperatures

All millimeter-wavelength continuum images are presented in Fig. 1. The majority of the disks in our sample show an elongated emission shape, with large axis ratios and in several cases roughly constant surface brightness along the major axis, confirming that they are highly inclined.

Two disks in the sample, however, present a different shape. Haro 6-5B and IRAS 04158 show the presence of a ring and central emission peak. We evaluate the position of the rings in Sect. 3.2.1. As in previous CARMA observations (Sheehan & Eisner 2017), the band 7 image of IRAS 04158 reveals a highly asymmetric ring (brighter on the western side) with a compact source toward the center. This ring is not detected in our band 4 observations, mostly due to flux dilution intosmall beams. However, the band 4 observations clearly resolve the emission at the center of the ring into two-point sources. IRAS 04158 is a binary system. A detailed analysis of this source will be presented in a dedicated study (Ragusa et al., in prep.).

The total continuum flux density of each source is measured by integrating the signal, down to the 3σ noise level, withinelliptical apertures tailored to each source. The flux densities are reported in Table 3. We use 10% error values throughout Table 3, which correspond to the typical flux calibration errors of ALMA (see ALMA Technical Handbook²). We note that although the phase calibration method used for the extended configuration observations in band 7 was non-standard (see Sect. 2.2.1), the flux calibration followed the usual procedure and, accordingly, the flux calibration uncertainty should be nominal. We checked that the integrated fluxes recovered using the compact configuration band 7 observations or both compact and extended configurations jointly are consistent within 10%.

To further investigate the disk properties, we calculate the brightness temperature, T_B, maps of each source. To ease the comparison between bands, we estimate the brightness temperatures from the band 7 and band 4 maps computed at the same angular resolution (see Sect. 2.4). We do not use the Rayleigh-Jeans approximation. We show the brightness temperature profiles measured along the major axis in Fig. B.1 and report the peak brightness temperatures in Table 4.

We note that, even for the most resolved sources (i.e., least impacted by beam dilution), the inferred brightness temperatures are lower or comparable to those required to be in the Rayleigh–Jeans regime, respectively T > 16.2 K for band 7 and 6.8 K for band 4. This strengthens our choice not to use this approximation. Additionally, we find that the band 4 brightness temperatures are systematically lower than the band 7 ones.

Fig. 1 Images of all sources included in this study normalized to their peak intensity (reported in the bottom right corner of each image). Each column corresponds to a different band and each line displays two sources, separated by the vertical dashed line. We show the beam size (ellipse) and a 0.5′′ scale (dark line) in the bottom left corner of each panel.

3.2 Surface brightness profiles

3.2.1 Radial extent

To characterize the radial extent of the disks, we present cuts along the major axis of each disk in Fig. 2. The cuts are normalized to their maximum intensity. For Tau 042021, IRAS 04302 and HV Tau C, for which the band 4 emission is very noisy, we estimate the normalizing factor as the amplitude of the Gaussian best fitting the curves. All sources are well resolved along their major axis.

Morphology

For several sources, the brightness profile is flat along the major axis direction and drops steeply at the edges. This is particularly true for HH 30, HK Tau B, HV Tau C, and HH 48 NE. The lack of a central brightness peak further supports the idea that these disks are optically thick and highly inclined, so that we are viewing only the outer optically-thick edge. Conversely, the disks of Tau 042021, IRAS 04302, IRAS 04200, ESO-Hα 569, and ESO-Hα 574 show more centrally peaked emission without a clear plateau, suggesting that they are optically thinner and/or viewed with a lower inclination, less edge-on. We note that the radial brightness profile of HK Tau B shows hints of rising at the edges. While this is a marginal (3σ) feature, this may indicate the presence of a ring (or transition disk).

Finally, three disks (Haro 6-5B, IRAS 04158, and Oph 163131) show symmetric shoulders or clear evidence of ring like features. The main ring of Haro 6-5B peaks at ~0.29′′ (41 au) in both bands. Furthermore, a shoulder seen in band 7 is associated with a peak in the higher-resolution band 4 cut; this would correspond to a fainter ring located at ~0.77′′ (108 au). The brightness asymmetry of the outer ring in IRAS 04158 is very clear in the major-axis cut: its western side is about 3 times brighter than the eastern side. We also note that the ring is slightly off-centered compared to the central binary. The western side of the disk peaks at ~1.71′′ (239 au), while the eastern side peaks at ~2.28′′ (319 au) from the center. We fit Gaussians on each side of the disk in the radial profile and find that the eastern ringis about 1.7 times wider than the western ring (with full width half maximum, FWHM, of 1.98′′ and 1.18′′ respectively). We also find that Oph 163131 displays a relatively flat profile in the inner 0.5′′, but has a sharp central peak. It has symmetric shoulders at ~0.70′′ radius, which suggests that this disk contains two rings and is viewed slightly away from edge-on. Further modeling of this source, focusing on dust emission, will be presented in a separate analysis (Wolff et al., in prep.).

Size estimations

To determine the sizes, we first normalized the major axis cuts presented in Fig. 2, as mentioned in the beginning of this section. Then, for each source, we estimate the relevant 3σ noise level in the image with the worst S/N, either band 4 or band 7. This 3σ level is converted into a fraction of the peak, in percent. The sizes in both bands are then measured at the same level, which means at the same fraction of the peak. The size of the disks along their major axis are reported in Table 5. The errors correspond to a tenth of the beam size along the major axis direction. We verified that using the FWHM of the cut profiles instead of the 3σ levels yields similar conclusions. This is also the case for the cuts obtained from the maps computed with a similar uv- coverage.

Along with the major axis size, we also estimate the apparent sharpness of the disk edges, as measured in the image plane. To do so, we measure the fractional range of radius over which the millimeter emission decreases from 80% to 20% of the peak emission, called Δr∕r. We define Δr by Δr = |r_80% − r_20%| and the normalization radius by r = |r_80% + r_20%|∕2. We chose this flux range because the radial profiles can usually be well approximated by straight lines in this interval. We report the values of Δr∕r in the last column of Table 5. The values are the mean of the estimations in band 4 and band 7, and they are deconvolved by the beam size. Sharp outer edges have small Δr and hence small Δr∕r. For example, a typical beam with a FWHM of 0.1′′ would have Δr∕r ≃ 0.1 when calculated in a small disk similar to HK Tau B (r = 0.5′′), or Δr∕r ≃ 0.04 in a disk as large as IRAS 04302 (r =1.2′′). All the disks in our sample have Δr∕r values larger than 0.2, which corresponds to shallower outer edges than the typical beam.

We find two different families of objects when comparing the radial extent of the disks in band 4 and band 7. Four sources (Tau 042021, IRAS 04302, IRAS 04200, and Haro 6-5B) show a band 7 size more extended than the band 4 by more than 2 beams. When compared with the same angular resolution and uv-coverage, these sources have band 7 major axis size on average 12% larger than that of the band 4. As the band 7 shorter wavelength traces smaller particles than those probed by band 4, the smaller sizes in band 4 suggest that the larger particles have drifted inward relative to the smaller ones. The outer edges of these four sources are well resolved and they have an average apparent sharpness of Δr∕r ~ 0.8, much shallower than the beam.

For the three remaining sources (HH 30, HK Tau B, HV Tau C), no difference in radial extent is found between band 7 and band 4. This might suggest the presence of dust traps at the outer edges of these disks, which can slow radial drift and lead particles to stop at particular radial locations (Powell et al. 2019; Long et al. 2020). Including HH 48 NE, these four sources have the sharpest edges, with Δr∕r between 0.2and 0.5. The edges of these disks are only marginally shallower than the typical beam. We note that 3 out of these four systems are known binaries and dynamical interactions may also lead to sharp outer edges.

Fig. 2 Normalized major axis intensity profiles. Data in the three bands are represented by the green (band 4), orange (band 6) and blue (band 7) lines. The light shading corresponds to the normalized rms in each band. The beam sizes in the direction of the cut are shown in the left part of each plot as horizontal lines. We indicate the shoulders of Haro 6-5B and Oph 163131 by an arrow. We smoothed the cuts through Tau 042021 by convolving them by a 1-D Gaussian of the beam width, to reduce the noise and make the plot clearer.

Fig. 3 Normalized averaged minor axis intensity profiles. For the most inclined sources (resp. least inclined: IRAS 04200, Haro 6-5B, and IRAS 04158), these were obtained by averaging minor axis intensity cuts over the whole major axis extent of the disk (resp. over the central 0.3′′). The light shading corresponds to the normalized rms in each band. The beam sizes in the direction of the minor axis are shown as dotted Gaussian of the corresponding color.

3.2.2 Disk extent perpendicular to the midplane

For the most inclined systems, the brightness maps shown in Fig. 1 have very elongated, linear shapes rather than elliptical ones, so we generate the minor axis profiles by taking the mean of the cuts at all distances along the major axis. In that case, the size of the minor axis is dominated by the vertical extent of the disk perpendicular to the midplane. For the less inclined sources where a clear ellipticity is visible in the image (namely IRAS 04200,Haro 6-5B, and IRAS 04158), we generate the minor axis profiles by averaging over a restricted range, only ± 0.15′′ around the center of the disk. In that case, the minor axis is dominated by the projection of the disk radius. We show the averaged brightness profile along the minor axis for all sources in Fig. 3. Dashed lines trace the Gaussian beam profiles along the cut direction.

We find that six out of the twelve disks of our sample are well-resolved along the minor axis, having a width more than twice the beam width. They are Tau 042021, IRAS 04200, Haro 6-5B (in both band 7 and band 4), IRAS 04158 (in band 7), and IRAS 04302,HK Tau B (in band 4). Additionally, the minor axis profiles of Haro 6-5B and IRAS 04158 reveal the presence of rings, and we see a clear asymmetry in the band 7 cut of Haro 6-5B.

We measure the minor axis sizes by fitting Gaussians to the generated profiles, Gaussians being good first approximations. For Haro 6-5B and IRAS 04158 for which the cuts show ring features, we estimate the FWHM directly through the cuts, without fitting a Gaussian. We report the resulting FWHM in Table 6, where the errors correspond to a tenth of the beam size projected in the direction of the cut³.

In order to extract the intrinsic vertical extent of the disks, we also deconvolve the minor axis sizes by the beam, assuming that both profiles are Gaussian. These values are presented in the last columns of Table 6. For the least resolved disks (i.e., for Oph 163131, ESO-Hα 569, ESO-Hα 574, and HH 48 NE), we only report upper limits. We estimate the uncertainties on the deconvolved minor axis sizes by propagating the errors, assuming that the error on the beam size is 10% of the beam major axis size.

We also estimate the disk inclinations from the measured axis ratio in all bands, and report them in Table 5. For consistency, we use the FWHM of the major axis sizes (cut at 50% of the peak flux) to estimate the inclination, as opposed to the size at 3σ reported in Table 5. For the sources that are not resolved in the vertical direction (Oph 163131, ESO-Hα 569, ESO-Hα 574, HH 48 NE) and those with the smallest axis ratio (Tau 042021, IRAS 04302, HK Tau B, HV Tau C), we present these values as lower limits. Indeed, for the most inclined systems (i.e. those with the smallest axis ratio), the minor axis size might not be dominated by the inclination but by the actual vertical thickness of the disks, which leads to lower apparent inclinations based on the axis ratio only. Except IRAS 04200, Haro 6-5B, and IRAS 04158, all resolved disks have an inclination larger than 75°. We note that Tau 042021 is the only highly inclined disk resolved along its minor axis in both band 7 and band 4. For this disk, the band 7 appears about 1.5 times more extended vertically than the band 4.

3.3 Estimation of spectral indices

The millimeter spectral index, defined as F_ν ∝ ν^α_mm, can be usedto study optical depth and grain growth in a disk (see e.g., review by Williams & Cieza 2011). Indeed, assuming scattering is negligible, the millimeter intensity can be expressed as I_ν = B_ν(T)(1 − e^−τ_ν), in which B_ν(T) is the Planck function and τ_ν the dust optical depth (which is proportional to the dust absorption coefficient, κ_ν ∝ ν^β). In the Rayleigh-Jeans regime and for optically thin emission, we expect α ≈2 + β ≥ 2. Small ISM-like grains have a β parameter of 1.5-2 (e.g., Li & Draine 2001), while grains of millimeter or centimeter sizes are expected to have a β parameter closer to 0 (e.g., Pavlyuchenkov et al. 2019). Low values of spectral indices (α ≤ 3) are usually interpreted either in terms of the emission being optically thick or that the dust grain size distribution has grown significantly to reach millimeter/centimeter for the maximum sizes (e.g., Testi et al. 2014).

Using the continuum fluxes from our survey and (sub)millimeter fluxes from the literature, we estimate the global millimeter spectral index α_mm for each source. We use a least-squares fit of all photometric points between 800 μm and 3.3 mm. We find α_mm ≤ 3 for all disks,as can be seen in Table 7.

Global spectral indices do not take into account the spatial distribution of the emission. Thus, for sources with multiple millimeter images, we computed spectral index maps using the band 4 and band 7 observations. We generated the spectral index maps pixel-by-pixel by applying the CASA task immath on the band 7 and band 4 maps computed to a unique resolution (see Sect. 2.4). Finally, to lower the noise level in the spectral index map, we applied a filter to keep only the pixels with emission well above the noise level (5σ) in both the band 4 and the band 7 images. The final maps and cuts along the major and minor axis are displayed in Fig. 4.

We find that the spectral index increases with radius for most sources, albeit with larger uncertainty at larger radii due to the lower S/N. Similarly, the spectral index also rises along the minor axis direction for two disks, Tau 042021 and IRAS 04200. Previous studies identified similar increases in the radial direction in several disks seen at lower inclinations (e.g., Pinte et al. 2016; Dent et al. 2019). The increases were attributed to changes in the dust size distribution and/or to lower optical depths at large radii. Very inclined systems on the other hand appear optically thicker than low inclination ones for the same mass. This is due to projection effects, because the line-of-sight crosses the disk over a longer distance. This suggests that, for the most inclined disks of our sample, spectral index variations are dominated by changes of optical depth insidethe disk (i.e., opacity effects) rather than by grain growth (see also our radiative transfer model, Appendix C). However, at the outer edges (both radially and vertically), where the disks become optically thinner, spectral index variations are enhanced by changes in the dust size distribution (see e.g., Tau 042021, IRAS 04302, IRAS 04200,and Haro 6-5B, where we found large differences in major axis sizes between band 7 and band 4). We also note that HH 30, IRAS 04302, HK Tau B, and HV Tau C show the opposite behavior of Tau 042021 and IRAS 04200along the minor axis, their spectral indices decrease. However, none of them is well resolved along the minor axis at the resolution of the restoring beam, so variations can be more affected by beam dilution and are less reliable.

Finally, we point out that we obtain spectral index values lower than 2 in the innermost regions of two disks: namely Tau 042021 and IRAS 04200 (see also our radiative transfer model, Appendix C). Such low values have also been reported in the inner regions of other disks (e.g., Huang et al. 2018a; Dent et al. 2019), and have often been interpreted as flux calibration errors because in the Rayleigh–Jeans regime α should not be smaller than 2. However, considering a 10% flux calibration error (yellow shaded regions in Fig. 4), the lowest α values measured in Tau 042021 and IRAS 04200 cannot be reconciled with α = 2. Recent studies showed that low dust temperature (outside the Rayleigh Jeans regime) or dust scattering in optically thick regions can reduce significantly the emission of a disk. In both cases, the spectral index can be even lower than 2 (e.g., Liu 2019; Zhu et al. 2019). In the context of highly inclined disks such as Tau 042021 and IRAS 04200, which are optically thick, both explanations are equally valid to explain low spectral indices observed. Modeling is needed to determine which one is dominant.

To summarize, because of the high inclination of our sources, we interpret the observed variation of spectral indices along the major axis as being dominated by optical depth changes in the disks. Additionally, we find that the low spectral index values measured in Tau 042021 and IRAS 04200 can either be related to low dust temperatures or to dust scattering in optically thick regions.

Fig. 4 Top row: spectral index maps, applying a filter to keep only the pixels with more emission than 5σ in both bands. The beam size is shown in the bottom left corner, along with a dark line representing a 0.5′′ scale. Middle row: spectral index cuts along the major axis. Bottom row: spectral index profiles along the minor axis, averaged as done for Fig. 3 (see text Sect. 3.2.2). For all cuts, yellow errors correspond to a flux calibration error of 10% in both bands, while orange errors are estimated from the signal-to-noise in each band. The x-axis corresponds to the offset to the center of the disk in arcseconds.

3.4 Comparison with optical and NIR images

We present overlays of optical and NIR images with our band 7 observations in Fig. 5. For most disks, we use HST opticalimages but prefer (space-based or ground-based) NIR images in a handful of cases to reduce confusion with extended nebulosity(see references in Fig. 5). All scattered light images show the same characteristic features, with two bright reflection nebulae separated by a dark lane tracing the obscuration of direct starlight by the edge-on disk. As opposed to the scattered light images, the millimeter data appear as very flat disks. All sources are found much less extended vertically in the millimeter than in scattered light, the result of a combination of opacity effects and vertical settling. Most of them also appear less extended radially, which can be linked to dust radial drift or opacity effects.

We estimate the scattered light major axis sizes by following the spine of each nebula (see method in Appendix D) and report the inferred radial sizes in Table 5. We could not estimate the scattered light sizes for two disks of the sample: for IRAS 04302 because we do not see the disk but the envelope; and for IRAS 04200 because the bright point source in the northern nebulae prevented the method to converge. Also, we indicate that our scattered light radial sizes might be under-estimated because lower illumination or lower sensitivity in the outer regions might reduce the apparent optical-NIR size (see for example Muro-Arena et al. 2018). A complete analysis will require the use of tracers of the gas distribution, which we postpone to a future paper.

Despite this, we find that most sources are larger radially in scattered light than in the millimeter, albeit with a few exceptions, the most obvious case being HV Tau C. ESO-Hα 574, and HH 48 NE, although formally more compact in scattered as estimated with our algorithm, have very similar sizes and will require deeper millimeter data for confirmation.

The ratios of scattered light over thermal continuum band 7 sizes are between 0.7 and 2.0. This is in general consistent with predictions from radial drift theory: the objects with clear evidence for radial drift being those in which scattered light disk (small grains) is significantly larger than the millimeter continuum disk (large grains). While we do not quantify explicitly the vertical extent of the scattered light images, Fig. 5 clearly shows that all disks are significantly more extended vertically in scattered light than at millimeter wavelengths, which indicates vertical settling has occurred in each disk.

Fig. 5 Overlay of scattered light (colors, in logarithmic scales) and ALMA band 7 continuum images (5, 10, 20, and 50σ contours) for all sources in this study, except for Oph 163131 for which we show the band 6 image. The scattered light images are plotted between 3σ (except for IRAS 04302 and IRAS 04158, where we use respectively 60 and 1σ to increase the contrast) and their maximum brightness (except for IRAS 04158, for which we take a lower value to increase the contrast). The ellipse and horizontal line in the bottom left corner indicate the beam size of the ALMA image and a 0.5′′ scale. References. Scattered light images, 2.2 μm: HV Tau C: Duchêne et al. (2010); 1.9 μm: IRAS 04302: Padgett et al. (1999); 1.6 μm: Haro 6-5B: Padgett et al. (1999); 0.8 μm: Tau 042021: Duchêne et al. (2014), HH 30: Watson & Stapelfeldt (2004), IRAS 04158: Glauser et al. (2008), ESO-Hα 569: Wolff et al. (2017), HH 48 NE: Stapelfeldt et al. (2014); 0.6 μm: IRAS 04200: Stapelfeldt et al. (in prep.), Oph 163131: Stapelfeldt et al. (2014), ESO-Hα 574: Stapelfeldt et al. (2014); 0.4 μm: HK Tau B: Duchêne et al. (in prep.).

Fig. 6 Radiative transfer models computed at 0.89 mm (band 7), for different inclinations. Top and third row: high and low mass settling models (h1mm = 0.70 au at r = 100 au), second and bottom row: high and low mass No Settling models (h1mm = hgas = 10 au at r = 100 au). The beam size and a 0.5′′ scale bar are shown in the bottom left corner of each panel.

4 Discussion

In this section, we use the brightness profiles of the highly inclined disks of our sample to discuss critically the amplitude of vertical dust settling, radial drift, and the effects of the enhanced optical depth due to projection effects. To this effect, we constructed several toy disk models and performed radiative transfer with the MCFOST code (Pinte et al. 2006, 2009) to produce synthetic images for comparison with the data. The toy models include a disk without rings or gaps and have an outer radius of 140 au (1′′). We assume that the surface density follows a truncated power law. The synthetic images are computed with and without vertical dust settling.

For the model without settling (NS), the dust is well mixed with the gas and assumed to have a scale height of 10 au at a radius of 100 au. For the model with dust settling (S), we assume that, as a function of size, dust follows the vertical density profileprescribed by Fromang & Nelson (2009). Following the results from Pinte et al. (2016) for HL Tau, a very flat disk when observed with ALMA, we set the vertical distribution of the millimeter dust (expressed in terms of a “scale height”) to be h_1mm = 0.7 au at 100 au. This corresponds to a disk viscosity coefficient of α = 3 × 10⁻⁴.

The toy model is also calculated for two different dust masses in order to probe the effect of optical depth. We consider intermediate to high mass disk models based on the lower limits on the dust mass derived for our sample of highly inclined disks (Sect. 4.3.1). We use M_dust = 1 × 10⁻³ M_⊙ for high mass disk models (HM), and M_dust = 5 × 10⁻⁵ M_⊙ for low mass models (LM). In total four different sets of images are calculated. We use a distance of 140 pc.

We computed the models at 0.89 mm (resp. 2.06 mm) for inclinations between 75° and 90°, and produced synthetic images using the CASA simulator with the same uv-coverage as our band 7 (resp. band 4) data. The synthetic band 7 images (0.89 mm) of each model are presented in Fig. 6. The band 4 images are qualitatively similar as those presented in Fig. 6 but slightly less extended in the vertical direction.

We also generated band 7 and band 4 model images with a unique angular resolution using a uv-taper. Using these maps, we computed the brightness temperatures and spectral index maps of the models, which are presented in Appendices B and C.

Fig. 7 Normalized averaged minor axis profiles for the high and low mass settling models (solid lines) and No Settling models (dashed lines) computed at 0.89 mm (band 7) and 2.06 mm (band 4), at 80° (left panel) and 90° (right panel).

4.1 Constraints from the surface brightness profiles

4.1.1 Vertical extent

For the high and low mass models presented in Fig. 6, we see that the disks appear less extended in the minor axis direction when settling is included. For the high mass model without settling (model NS-HM), we find that as the inclination approaches 90°, the high optical depth in the midplane results in a clear low intensity lane, separating the two sides of the disk. On the other hand, at lower inclinations the bottom (far) side of the disk is about 5–6 times fainter than the top (near) side. This is visible directly in Fig. 6 and highlighted in a different way by showing cuts along the minor axis in Fig. 7 for inclinations of 80° (fainter back side) and 90° (dark lane). As the angular resolution is similar and the signal-to-noise is greater than 14 for all targets included in this study, such an asymmetry would be detectable easily in our observations. However, neither the vertical asymmetry nor the dark lane in the midplane are present in our data. This configuration (NS-HM) can be ruled out. Interestingly, the features can hardly be seen in the high mass model which includes settling, nor in any of the low mass models. If the disks included in our sample are as massive as the high mass model, this would indicate that they have a small millimeter dust scale height, closer to 1 au than to 10 au at a radius of 100 au. This is similar to the results of Pinte et al. (2016) for the disk of HL Tau.

To go further, the deconvolved minor axis sizes of the models can be compared with the average size obtained for the most inclined disks of our sample (i.e., with i > 80°). For the data, the mean deconvolved minor axis size is about 0.18′′ in band 7 (0.17′′ in band 4, see Table 6). Except for the high mass model without settling (model NS-HM), which is more than twice as thick as the observations (deconvolved minor axis size ~ 0.4′′ at 80° and 90°), all models studied in this section are compatible with the vertical extent measured in the data (in band 7, typical deconvolved sizes of ~ 0.2′′ for S-HM, NS-LM at 80° and 90°, and for S-LM at 80°, and ~ 0.1′′ for S-LM at 90°). Thus, the vertical extent alone is not sufficient to distinguish whether vertical settling is required to explain the observations or not.

Fig. 8 Normalized major axis intensity profiles for the high mass (left) and low mass (right) models at 0.89 mm, for different inclinations between 75° and 90°.

4.1.2 Radial brightness profile

We now investigate the effect of inclination on the observed radial brightness profiles. To do so, we produce major axis cuts of our radiative transfer models, as in Sect. 3.2.1, and present them in Fig. 8. We also estimate the apparent sharpness between 20% and 80% of the peak flux (Δr∕r) for our models as in Sect. 3.2.1, and report them in Table 8.

Along the major axis, the effect of inclination on the shape of the brightness profile is very clear both in the images (Fig. 6) and in the cuts (Fig. 8). For all models, at the lowest inclinations,the images and cuts show a steep increase in intensity at the center, and low apparent edge sharpness (Δr∕r > 1 for i < 80°). On the contrary, for the fully edge-on configuration (i = 90°), the major axis brightness profile of the high mass models is flat at all radii and drops steeply (Δr∕r ~ 0.2 at 90° for both high mass models S-HM and NS-HM). Between these extreme cases, the cuts show a less extended plateau and shallower outer edges than for the 90° case, relatedto lower optical depth than in the edge-on case. We note that in the low mass models, which are optically thinner, the flat plateau along the major axis direction is never reached and the apparent disk sharpness is always larger than Δr∕r > 0.6. Nevertheless, we find that, independently of the dust mass assumed for the models, the apparent sharpness of the disk outer edge increases with increasing inclination.

Three disks in our sample show edges as sharp as 0.3. They are HH 30, HK Tau B, and HH 48 NE. The comparison with models suggests that these disks are the most inclined of the sample. Additionally, these disks present flat radial brightness profiles that the low mass models are unable to reproduce. This indicates that they are optically thicker than the low mass model, likely because they are more massive. In Sect. 4.1.1, we have shown that a more massive disk, while leading to a flat profile along the major axis direction also presents a larger vertical extent. Our high mass model without settling is inconsistent with the observations (see Fig. 6, NS-HM model). This implies thatvertical settling is needed to explain both the flat radial profile and small apparent vertical extent of these three disks. For these systems the millimeter scale height would be closer to 1 au than to 10 au at 100 au, therefore increasing significantly the concentration of dust mass in the disk midplane. Although we cannot confirm vertical settling for the other sources, we believe that their vertical structure is likely similar and governed by settling, in particular because of the large vertical size difference between the scattered light and thermal emission images.

Vertical settling models show that the turbulence generated self-consistently by ideal MHD or vertical shear instabilities are likely too strong to lead to millimeter scale height as small as 1 au at r = 100 au for gas scale heights of 10 au at 100 au (Flock et al. 2017, 2020). On the other hand, non-ideal effects such as ohmic resistivity or ambipolar diffusion lead to lower turbulence and thus to very thin millimeter grain layers (Riols & Lesur 2018). Those mechanisms may well be dominant in the sample of disks analyzed here. Detailed modeling of each individual object is needed to obtain a quantitative estimate of the millimeter dust and gas scale height of the disks.

4.2 Comparison with a radial drift model

Similarly to previous multi-wavelength studies (e.g., Pérez et al. 2012; Tripathi et al. 2018; Powell et al. 2019), our observations show that 4 disks have major axis sizes that decrease with wavelength (namely Tau 042021, IRAS 04302, IRAS 04200, Haro 6-5B, see Table 5 and the radial variation of their spectral index in Fig. 4). For these disks, the average size difference between band 7 and band 4 observations is about 12%. Estimating the major axis sizes from our radiative transfer models described in Sect. 4 (which do not include radial drift), we find that opacity effects alone predict a difference of only a few percent between bands, which is not sufficient to reproduce the observations. In this section, we compare the measured radial differences with an analytical radial drift model, presented by Birnstiel & Andrews (2014). Similarly to other theoretical models, they showed that inward dust migration of single size particles spontaneously produces a sharp edge in the dust density distribution (Birnstiel & Andrews 2014; Facchini et al. 2017; Powell et al. 2019). Birnstiel & Andrews (2014) computed an analytical formula to infer the position of the disk outer edge in a disk with radial drift only (see their Eq. (B9)). The vertical extent of the grains is not considered in their model. They assume a smooth tapered-edge gas surface density profile and parametrize the turbulence following Shakura & Sunyaev (1973).

Assuming that the band 7 and band 4 emission only originate from grains of the optimal size (a ≈ λ∕2π), the model by Birnstiel & Andrews (2014) predicts a size difference in surface density between band 7 and band 4 of about 25% (after 0.1 Myr). While the predicted effect is marginally too strong, the absolute disk sizes are more problematic. When grains of 0.1 μm detectable in scattered light are expected to be found up to 135 au after 0.1 Myr, the model predicts that grains emitting most at 0.89 mm should have drifted to 42 au. This corresponds to a micron/millimeter disk radius ratio greater than 3, which is more than 1.5 times larger than the largest ratio measured in our data. Part of these differences might be explained because the observations are not probing the surface density of the disk, because several grain sizes (that might have drifted to different radii) have to be considered rather than a unique one, or because our disks are older than 0.1 Myr so they might be affected by viscous spreading as well. Besides, as discussed previously (Sect. 3.4), we note the scattered light images might not always trace the whole disks, since they require illumination (and sufficient optical depth) to trace the disk all the way to the edge. This can lead to apparent sizes in scattered light that are smaller than the real radial extent of small grains. However, the small expected sizes at millimeter wavelengths suggest that this radial drift model is too efficient to reproduce the observations (see also Brauer et al. 2007). The existence of pressure fluctuations in the disks (e.g., rings and gaps tracing pressure bumps and dust traps) rather than smooth power law surface density profiles are expected to slow down the drift efficiently (see, e.g., Gonzalez et al. 2017; Pinilla et al. 2012) and help reconcile models with observations. We speculate that the disks included in our study may have complex radial structures to slow down the radial migration of dust.

4.3 Effect of inclination on global values

Most Class II disk surveys in close-by star-forming regions estimate dust masses directly from the measured integrated fluxes assuming optically thin disks, irrespective of the disk inclinations (e.g., Andrews et al. 2013; Ansdell et al. 2016; van der Plas et al. 2016). However, because of the projection effect, the optical depth along the line-of-sight increases with inclination and will affect the observed flux. In this section, we investigate the impact of inclination on integrated fluxes, dust mass estimations, and integrated spectral indices. Finally, we also discuss the measured brightness temperatures obtained for our disks in the context of optical depth.

4.3.1 Flux density and derived masses

We present the variation of the integrated band 7 and band 4 fluxes in our high mass radiative transfer models (with and without settling) as a function of inclination in Fig. 9. The effects of dust scattering are fully included in the radiative transfer calculations. Overall, the emitted flux density of the disk becomes attenuated by up to an order of magnitude with increasing inclination. This is due to a combination of the increasing optical depth, a lower average dust temperature seen by the observer for high inclinations, and to geometrical effects (reduced emitting surface with increasing inclination). In this model, from ~60° to 90°, the flux density and thus the derived disk mass would appear significantly smaller than for the same disk viewed at intermediate inclinations (by up to ~10 times for the high mass settled model). The amplitude of the attenuation depends on the parameters of the model. We note that in Fig. 9, the model without settling is brighter than the corresponding settled model for all inclinations. This is due to temperature differences between the settled and no settling models, to differences in optical depths and to geometrical effects. Similar behaviors are seen in the low mass models which are not represented.

This raises the question of the reliability of disk masses derived purely from millimeter fluxes, in particular for surveys where the disks are not well resolved and inclinations cannot be estimated. We focus here on the highly inclined disks of our sample. For a direct comparison with previous studies, we assume optically thin and isothermal dust emission at sub-millimeter wavelengths. In that case the flux (F_ν) is directly related to the dust mass following: M_dust = F_νd²∕κ_νB_ν(T_dust), where d is the disk distance to the Sun, B_ν(T_dust) the Planck function evaluated at a mean representative dust temperature T_dust, and κ_ν the dust grain opacity. We use typical values of T_dust = 20 K and κ_ν = 3.4 cm² g⁻¹, as in Ansdell et al. (2016, and references therein) to estimate the dust mass of the disks in the sample, using the band 7 fluxes. For the data, the results are reported in Table 9. As noted above, these masses are almost certainly underestimated and so we quote them as lower limits.

By applying the formula directly we find a mean dust mass of about 25 M_⊕ for the 6 edge-on disks more inclined than 75°. For comparison, at an inclination of 90°, the estimated dust mass of both the high mass and low mass settled models are more than 3 times lower than the real dust mass. Without full modeling, it is difficult to find a reliable correction factor for individual sources to compensate the attenuation due to inclination. The disk masses of our highly inclined sample are probably a few times larger, that is up to ~75 M_⊕ when applying the same factor of 3. This is on the high end of the dust mass distributions of Taurus or Lupus star-forming regions (mean of 15 M_⊕, Ansdell et al. 2016).

This result suggests that our sample of edge-on disks is biased toward more massive disks. This is a direct consequence of the important attenuation caused by the high inclination, making the starlight and the disk emission fainter. This may also explain, at least partly, why the number of known edge-on disks with resolved images remains sparse, even today. A comprehensive study of the biases affecting the edge-on disk population and their detection will be presented in Angelo et al. (in prep.).

Fig. 9 Top: variation of the band 7 and band 4 integrated fluxes of the high mass radiative transfer models with settling (thick lines) and without settling (thin lines) as a function of the inclination. Bottom: variation of the integrated spectral index with inclination for both high and low mass radiative transfer models, with or without including settling.

4.3.2 Integrated spectral indices

The median spectral index of the edge-on disks in our study, ${\overline{\alpha }}_{\text{mm}}=2.5\pm 0.3$, is similar to previous measurements of inter- mediate-inclination disks, for example 2.3 ± 0.1 in the Lupus and Taurus star-forming regions (Ribas et al. 2017; Ansdell et al. 2018). At first glance, this is surprising as edge-on disks are expected to appear more optically thick, but Sierra & Lizano (2020) showed that the spectral index is only mildly dependent on inclination based on an extensive study.

To further understand the integrated α_mm for edge-on disks, we computed them for our radiative transfer models. The expected values at different inclinations are shown in the bottom panel of Fig. 9. We find that the variation of spectral index with inclination is relatively small (1.9 ≲ α_mm ≲ 2.2 for the high mass settling model, and 2.1 ≲ α_mm ≲ 2.5 for the low mass settling model for inclinations lower than 87°), in agreement with Sierra & Lizano (2020). This is because the models are (at least partially) optically thick even at low inclination and because large grains are included (maximum grain size is 3 mm). The small difference in integrated spectral index observed between our sample of edge-on disks and disks at lower inclinations can be explained if most disks are partially optically thick at millimeter wavelengths and/or if grains have grown to millimeter/centimeter sizes. Grain growth is known to have occurred in Class II disks, leading to low spectral indices (Ricci et al. 2010; Testi et al. 2014). Similarly, recent imaging campaigns have revealed that ring structures are very common in protoplanetary disks (Huang et al. 2018b) and generally associated with optically thick regions with large grains (e.g., Dent et al. 2019). So far, the edge-on disks in our sample appear similar to the disks observed on other surveys, except from their viewing angle.

Interestingly, for inclinations greater than 70° for the high mass or 87° for the low mass model, we find a clear difference in the integrated spectral index between the models with and without settling. While at these inclinations the spectral index decreases in the models without settling (see also Galván-Madrid et al. 2018), we find that the integrated spectral index increases for the settled models, reaching about 2.2 at 90° for the high mass model (2.8 for the low mass model). This can be explained by the contribution of the optically thick midplane that decreases for increasing inclination (as it appears colder) while the contribution of the upper layers (rich in small grains with higher β values) increases, leading to higher values of the integrated spectral indices for the settled model (see, e.g., Fig. C.1). Said differently, for very large inclinations several line-of-sights do not cross the disk midplane containing large grains but only the high-altitude layers above and below it. Since these layers are optically thinner and contain only small grains, the spectral index increases. However, further studies are required to investigate this effect (e.g., Sierra & Lizano 2020) and its applicability to specific objects.

4.3.3 Brightness temperatures

Assuming that scattering is negligible and that the dust temperature is high enough (e.g., in the Rayleigh-Jeans limit) and constant over the emitting region, the brightness temperature of dust emission at frequency ν can be expressed as T_B = η_c T_p (1 − e^−τ_ν), where τ_ν is the optical depth of the medium, T_p the mean dust physical temperature (see e.g., Wilson et al. 2009), and η_c the fractionof beam covered by the source. In the isothermal optically thick limit (τ_ν ≫ 1), for a source fillingthe beam (but smaller than the largest angular scale of the interferometric observations) the brightness temperature corresponds to T_p − 2.7 K, because the cosmic microwave background (CMB) is resolved out by the interferometer. For compact sources, beam dilution would reduce the observed brightness temperatures below T_p. Scattering is also known to decrease dust emission from very optically thick regions, which would also effectively lead to lower observed brightness temperatures (Zhu et al. 2019). The brightness temperatures were estimated in Sect. 3.1. Major axis cuts are presented in Fig. B.1 and we report the peak brightness temperatures in Table 4.

We estimate a mean peak brightness temperature of 10.6 K for the three best resolved sources in band 7 (Tau 042021, IRAS 04200,and Haro 6-5B). These values are particularly low compared to previous estimates on other disks around stars of similar spectral types. As an example, Andrews et al. (2018) derived a mean brightness temperature peak of 66.5 K for the 17 disksof the DSHARP sample around K & M stars observed in band 6. These disks have a mean inclination of 42°, while all disks of our study are more inclined than 62°. The brightness temperatures we derive for our disks are also much lower than the traditional T_dust = 20 K assumed in flux-to-mass conversions (see e.g., Sect. 4.3.1).

Although the measured brightness temperatures are integrated over some vertical extent because of the beamsize, and therefore include a vertical temperature gradient, for the well-resolved disks T_B does provide a reasonable estimate of the temperature of the outer midplane where the line-of-sight optical depth reaches unity. This is the case for Tau 042021 for example. The low brightness temperatures measured in these optically-thick edge-on systems is likely to reflect the midplane temperature in the cold outer radii of the disks. This in agreement with the idea that the disks appear optically thick because of projection effect.

Fig. 10 Left panels: ALMA data of HK Tau B. In all the figure, top row is band 7 and bottom row is band 4. Middle-left panels: millimeter images of the best model C of Stapelfeldt et al. (1998), assuming well mixed grains. They were computed for the same uv-coverage as in the data. Middle-right panels: major axis cuts of the model and the data. Right panels: averaged minor axis cuts of the data and the models, performed as in Sect. 3.2.2.

4.4 Individual targets

4.4.1 HK Tau B: Comparison with a published model

HK Tau B is the only Class II disk of the sample which is resolved by ALMA in the minor axis direction (in band 4) and for which the scattered light image has been modeled previously. Stapelfeldt et al. (1998) estimated the scale height of micron-sized grains to be on the order of 3.9 au at 50 au (8.3 au at 100 au, assuming a flaring exponent of 1.1 as in Stapelfeldt et al. 1998). Aiming for a quantitative comparison of scale height between millimeter and micron-sized grains, we computed their best model C at millimeter wavelengths (at 0.89 mm and 2.06 mm) using the radiative transfer code MCFOST. The parameters of the model are reported in Table 1 of Stapelfeldt et al. (1998) and the 0.8 μm image of the model can be found in their Fig. 4. Following the results of Duchêne et al. (2003), we adopt revised opacities in our model, leading to a revised disk dust mass of 5.9 × 10⁻⁵ M_⊙ to match the observed millimeter fluxes. We use silicate dust grains with sizes from 0.01 μm to 15 cm and assume a number density described with a power law of the grain size dn(a) ∝ a^−3.5da. This yields opacities of 4.8 and 2.6 cm² g⁻¹ (per dust mass) respectively at 0.97 and 2.06 mm. Assuming a gas-to-dust mass ratio of 100, this is comparable to standard assumptions (e.g., Beckwith et al. 1990, 0.048 and 0.026 cm² g⁻¹ of gas and dust, respectively). For a direct comparison, we assume that grains of all sizes are fully mixed and produce synthetic images using the CASA simulator with the same uv coverage as in our observations. We show the results in Fig. 10.

We find that the model does not reproduce well the millimeter images of HK Tau B. Along the major axis, the model is too peaked at the center (less flat-topped) and more extended in the radial direction than the actual millimeter data. As discussed in Sect. 4.1.2, the more peaked profile suggests that the observed millimeter disk is more optically thick than what is predicted by the model. While the best model of Stapelfeldt et al. (1998) contains a full disk at an inclination of 85°, we suggest that a higher inclination, closer to 90° (see Sect. 4.1.2), or the presence of a large but unresolved inner cavity would reproduce better the major axis brightness profile.

Along the minor axis, the band 4 cut indicates that the model is also more extended vertically than the data (see Fig. 10). As previously proposed by Duchêne et al. (2003) this hints for vertical segregation of dust grains in this disk, millimeter grains being located in a vertically thinner layer than micron-sized grains. However the low signal-to-noise of the band 4 data and the low angular resolution in band 7 prevent us from a strong conclusion.

4.4.2 Vertical settling in Tau 042021

Tau 042021 is the only edge-on disk clearly resolved vertically in both band 7 and band 4. For this disk, the band 7 appears about 1.5 times more extended vertically than the band 4 (see Table 6). We note that Tau 042021 is only marginally resolved in band 6, so the size at this wavelength is uncertain. For our high mass radiative transfer model without settling (i.e., including only optical depth effects), we find a band 7 to band 4 minor axis ratio of about 1.1 at i = 90° (respectively 1.2 for the low mass settled model), significantly smaller than the observed value. Although this model is not unique, it indicates that opacityeffects are not sufficient to produce the large difference in minor axis size observed, which suggests that grain-size-dependent vertical settling is occurring in this disk as well. In this section, we compare the measured minor axis ratio with predictions from several vertical settling models. We assume that the disk is perfectly edge-on so that variations along the minor axis are dominated by differences in vertical extent between bands rather than projections of the disk radius.

From the minor axis ratio measured in Tau 042021, one can estimate the scaling of the minor axis size (S_d) with grain size (a) assuming S_d ∝ a^−m. Assuming that most of the emission comes from grains of the optimal size (a ≈ λ∕2π), we obtain: m = −log(S_d,b7∕S_d,b4)∕log(a_b7∕a_b4) ~ 0.5.

If we additionally assume that S_d is directly proportional to the dust “scale height” (h_d), an exponent of m = 0.5 has been predicted for large grains in the context of a 1-D diffusion theory (Dubrulle et al. 1995). Performing numerical simulations including non-ideal MHD effects such as ambipolar diffusion, Riols & Lesur (2018) also estimated a relationship of dust scale height with grain size with an exponent of 0.5, valid for large grains (t ≳ 10⁻²). For comparison, using the standard disk model presented in Riols & Lesur (2018, Eq. (13)), at 100 au, grains of the optimal size emitting in band 7 would have Stokes numbers larger than 2.5 × 10⁻².

On the other hand, Fromang & Nelson (2009) estimated that in the case of ideal MHD, the dust scale height varies as: h_d ∝ a^−0.2. Settling obtained with ideal MHD is expected to be less efficient than other models previously discussed. However, we note that this expression was estimated for small grains (St ≤ 10⁻²) and might not apply for the grain sizes that we are probing at millimeter wavelengths.

To summarize, we find that the large minor axis difference observed in Tau 042021 at different millimeter wavelengths has to be associated with strong settling. To compare our observations with settling models, we make the assumption that the measured minor axis extent is proportional to the real dust scale height. While ideal MHD simulations predict a settling less efficient than observed, a simple 1-D diffusion theory (in the case of strong settling, Dubrulle et al. 1995) or numerical simulations including non-ideal MHD effects (Riols & Lesur 2018) provide for now the best consistency with the observations.

Fig. 11 Minor axis size as a function of radius for IRAS 04302 (band 4) and Tau 042021 (band 7). The errors correspond to the standard deviation of the curves.

4.4.3 Radial variation of vertical extent in IRAS 04302

The high-resolution band 4 disk image of IRAS 04302 shows evidence of flaring: the extent perpendicular to the disk midplane at large radii is broader than at the center (see Fig. 1). This is in contrast with the other disks in our survey, which are flatter or unresolved. Fig. 11 compares the minor axis sizes at different radii in IRAS 04302 with the other well-resolved edge-on disk with large signal-to-noise: Tau 042021. In the latter case, the vertical extent appears quasi constant as a function of radius with only a slight increase in the outer regions, whereas it is ~50% larger in IRAS 04302 at a projected radius of 100 au (0.7′′) than at the center. This might be related to differences in optical depth or settling between the disks. Additionally, IRAS 04302 is the only Class I object in our sample (Gräfe et al. 2013) and the scattered light image does not show the flared upper surface of a disk, like in the other objects (Class II, see Fig. 5), but instead traces an envelope. Objects like this one, and HH 212 a Class 0 (Lee et al. 2017), open the possibility to do comparative studies of disk evolution in the early phases.

5 Summary and conclusions

We presented high angular resolution ALMA band 7 and band 4 observations of 12 highly-inclined disks previously selected from the shape of their scattered light images. All disks are well resolved along the major axis and 6 are also resolved in the direction perpendicular to the disk midplane in at least one millimetric band. Several disks show flat surface brightness profiles along their major axis with a steep drop off at their outer edge, indicating their large inclination and significant optical depth.

Haro 6-5B and IRAS 04158 are the least inclined disks of the sample (less than 75°) and both show a distinct ring and an isolated emission peak at the center. At the highest angular resolution, the point source at the center of IRAS 04158 is a binary source.

The analysis of global quantities such as integrated fluxes, spectral indices, and brightness temperatures shows that the highly inclined disks of our sample have (at least partly) optically thick emission. Because of the low brightness temperatures and small beam sizes, we conclude that the emission originates from the outer radii of the disks, with a peak brightness temperature below 10 K for half the sources in band 7.

We also found that the median spectral index in our disk sample is similar to that of disks seen at lower inclinations.This can be explained if disks at intermediate inclination are already partly optically thick (which implies significant scattering even at millimeter wavelengths) and/or if grains have grown to millimeter/centimeter sizes in all cases.

All disks were observed at several wavelengths with similar angular resolution, from the optical to the millimeter range. This enables a comparison of the radial extent of different grain populations in the disks (i.e., grain sizes). We assumed that the small dust responsible for the scattered light in the optical is tracing closely the gas distribution. Most disks have larger radial sizes in the optical-NIR than at millimeter wavelengths indicative of dust radial drift, the larger particles having drifted inward. Three of the disks have the same radial extent in both millimetric bands; theseones also have the sharpest apparent outer edges (estimated between 20% and 80% of the peak flux): Δr∕r ~ 0.3. Four sources have band 7 emission which is radially more extended than band 4, by about 12% on average. However, current radial drift models predict larger differences – both between optical and millimeter, and between band 7 and band 4 – than we actually observe. This suggests that other mechanisms such as pressure traps are likely present in these disks to slow down or halt the radial drift.

The peculiar viewing angle of the disks presented in this survey allows us to obtain more direct information on their vertical structures. First of all, the direct comparison of the ALMA observations with scattered light data shows that these disks have larger vertical sizes in the optical-NIR than at millimeter wavelengths, indicative of the different optical depths and of vertical dust settling. To further estimate the vertical distribution of millimeter grains (parametrized as a “scale height”), we compared the shape of ALMA observations with four radiative transfer toy models of different mass, which include or not vertical settling. We computed the models at high inclinations, with the same angular resolution as our data, and considered two different values for the scale height of millimeter grains. We find that at least three disks of our survey require that the millimeter dust “scale heights” is low, of order of a few au at r = 100 au: these disks are vertically thin at millimeter wavelengths. This is much thinner than the gas traced by the small dust, which has a typical scale height of 10 au (Burrows et al. 1996; Stapelfeldt et al. 1998).

On a case by case basis, for HK Tau B, a more detailed comparison of the ALMA images with a published scattered light model (re-computed at millimeter wavelengths) also suggests differences in vertical extent between millimeter and micron-sized grains, as previously suggested by Duchêne et al. (2003). Also, for Tau 042021, the only edge-on disk well resolved in the two millimetric bands, we find that the band 7 emission is about 1.5 times more extended vertically than the band 4. Assuming that the measured vertical extent is directly proportional to the dust scale height, this ratio is expected for relatively large dust grains in the simple 1-D diffusion theory or numerical simulations including non-ideal MHD effects (Dubrulle et al. 1995; Riols & Lesur 2018), further supporting the idea that strong vertical dust settling has taken place, leading to an increase in dust concentration in the disk midplane.

Finally, we find evidence of a more flared structure in IRAS 04302, suggesting that the millimeter grains in this Class I source are less settled. The millimeter dust in this disk may be in transition between the vertically unsettled structures seen in some Class 0 objects, and the flatter dust found in our Class II disks.

In forthcoming studies, we will present the CO gas distribution measured in these disks. We will also produce more detailed case by case models of selected targets to quantify the dust and gas density distribution profiles.

Acknowledgements

The authors thank the referee for the constructive comments which improved significantly the paper. The authors also thank Laura Perez and Eric Villard for useful discussions. M.V., F.M., M.B., G.vdP., C.P. acknowledge funding from ANR of France under contract number ANR-16-CE31-0013. This project was financially supported in 2019 and 2020 by the CNRS as part of its programme 80|PRIME. GD and ZT acknowledge support from NASA under grant 80NSSC18K0442 and NSF grant AST-1 518 332. KRS acknowledges funding from the Space Telescope Science Institute under GO grant 12514. CP acknowledges funding from the Australian Research Council via FT170100040 and DP180104235. F.L. acknowledges the support of the Fondecyt program N°3170360. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00460.S, ADS/JAO.ALMA#2016.1.01505.S, ADS/JAO.ALMA#2016.1.00771.S and ADS/JAO.ALMA#2013.1.01175.S. A.L.M.A. 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 ofKorea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ.

Appendix A Beam sizes

We present the beam sizes obtained for the different observations in Table A.1. The last column describes the beam size of the maps used to compute the brightness temperatures and the spectral index maps.

Appendix B Brightness temperature cuts

We present brightness temperature cuts along the major axis for all disks in Fig. B.1. Solid lines represent the band 7 and band 4 cuts computed with the same beam, for which we reported the peak temperature in Table 4. We also show the brightness temperatures calculated with the resolution of the images presented in Fig. 1 in light colors. We note that, as the beam sizes of the band 6 observations of Tau 042021 and HH 30 are significantly larger than the restored beam (see Table A.1), it is not possible to compare directly the band 6 brightness temperatures to band 4 and band 7, the beam dilution will be different. Moreover, a direct comparison of brightness temperature between sources is difficult as they would have different levels of beam dilution.

Fig. B.1 Observed brightness temperatures as a function of the radial distance to the central star. Solid lines profiles (band 7 and band 4 only) were computed with the same angular resolution (restored beam in Table A.1), while the light lines have the original angular resolution. We report the beam sizes in the direction of the cut as horizontal lines in the left part of each plot. The restored beam corresponds to the black line. Error bars correspond to 10% of the brightness temperature.

Fig. B.2 Brightness temperatures as a function of the radial distance for the radiative transfer models presented in Sect. 4, computed at different inclinations. Top and third row: high and low mass settling models, second and bottom row: high and low mass No Settling models. We report the beam size in the direction of the cut in the left part of each panel.

In Fig. B.1, we find that most sources do not show strong differences in brightness temperatures between the different angular resolution studied. Only the band 4 observations of Haro 6-5B and IRAS 04302 show important variations of the brightness temperature with the beam size. The temperatures obtained with the restored beam are smaller than those with the original angular resolution. For both objects, this likely indicates that the disk (or central point source for Haro 6-5B) is significantly less resolved with the restored beam than at the original angular resolution. This leads to a lower brightness temperaturedue to beam dilution.

In Fig. B.2, we present brightness temperature cuts at several inclinations computed for our radiative transfer models withand without including vertical settling. We also report the peak brightness temperatures obtained for our models in Table B.1. In each model, we find that as the inclination increases, the peak temperature decreases. This indicates that the disk becomes optically thicker with inclination. At the highest inclinations millimeter emission likely originates in the outer radii of the disk. We note that this is presumably also the case for several disks of our survey in which low brightness temperatures are measured even if they are well resolved (for example for Tau 042021).

Additionally, we find that, in the low mass models, the band 4 brightness temperatures are lower than the band 7 brightness temperatures. This is not the case for the high mass model where the band 4 brightness temperature can sometimes be higher than the band 7 temperature. In the observations, we detect large variations in brightness temperatures between bands (Fig. B.1, Table 4). While this suggest that the disks are optically thinner than the high mass models, further studies are needed to explain the differences in brightness temperature measured in the observations.

Appendix C Spectral index maps of the high mass settled model

In order to interpret the spectral index maps shown in Fig. 4, we computed the spectral index map of our high mass settled radiative transfer model presented in Sect. 4 (the most representative of our observations, see Sect. 4.1), for inclinations between 75° and 90°. We show the maps on the top row of Fig. C.1, and display the major axis cuts in the disk midplane in the middle row and the averaged minor axis profiles in the bottom row. As in Sect. 3.2.2, we obtained the minor axis profiles over the full major axis for the disks that do not show ellipticity (i.e., more inclined than 80°), and over ± 0.15′′ for the least inclined models.

Fig. C.1 Top row: spectral index maps of the high mass settled model, applying a filter to keep only the pixels with more emission than 5σ in both bands. The beam size is shown in the bottom left corner, along with a dark line representing a 0.5′′ scale. Middle row: spectral index cuts along the major axis. Bottom row: spectral index profiles along the minor axis, averaged along the major axis. For all cuts, yellow errors correspond to a flux calibration error of 10% in both bands, while orange errors are estimated from the signal-to-noise in each band. The x-axis corresponds to the offset to the center of the disk in arcseconds.

The spectral index maps of the models show several features also seen in the observations. First, for most inclinations, it reaches values lower than two in the midplane (as for Tau 042021 and IRAS 04200), which can indicate low temperatures or important scattering in the disk midplane. Second, at the highest inclinations, there is a slight increase in α_mm along the major axis, which is related to optical depth effects only. This variation is smaller than for several sources of our sample (e.g., Tau 042021, IRAS 04302, IRAS 04200 and Haro 6-5B) which suggests that these disks are affected by radial drift. Finally, for inclinations larger than 85°, the spectral index increases along the minor axis, similarly to what is measured in Tau 042021, which is likely related to optical depth effect and vertical settling.

Appendix D Estimation of scattered light sizes

Here we describe a method to estimate the radial extent of edge-on disks as observed in scattered light. To measure the extents of the scattered light images of the disks, we first find a set of points to define the spine of each nebula. To this end, we extract cuts along the minor axis of the disk. When this cut contains two clear local maxima, we compute the centroid in a small region surrounding the peaks. If the two nebula are partially blended, or hard to disentangle, we perform a two-Gaussian fit, forcing the same FWHM of the two nebulae. The center of each Gaussian then defines the location of the spines. The process is initiated at the disk’s axis of symmetry and we proceed outward on each side until the surface brightness drops below 1− 7% of the peak surface brightness, depending on the signal-to-noise of the image and adjusted to match the visually detectable edge of the disk. The result is two sets of points that define the spines of the bright and counter nebulae (blue and red points in Fig. D.1). Outliers, due to bad pixels or substructure in the disk, are excluded based on deviations from a running median. We then fit second-order polynomial functions to each nebula, which provides adequate morphological information while minimizing sensitivity to noise and low-level departures from symmetry. The radial extent (diameter) of the disk is then defined as the maximum distance between identified spine points along the major axis. We present an example of fit in Fig. D.1.

Fig. D.1 Example of scattered light size estimate in Tau 042021. We represent the two spines of the nebula by the blue and red curves. The adopted radial size corresponds to the horizontal line, the maximum distance between spine points along the major axis.

References

All Tables

All Figures

	Fig. 1 Images of all sources included in this study normalized to their peak intensity (reported in the bottom right corner of each image). Each column corresponds to a different band and each line displays two sources, separated by the vertical dashed line. We show the beam size (ellipse) and a 0.5′′ scale (dark line) in the bottom left corner of each panel.
In the text

Fig. 2 Normalized major axis intensity profiles. Data in the three bands are represented by the green (band 4), orange (band 6) and blue (band 7) lines. The light shading corresponds to the normalized rms in each band. The beam sizes in the direction of the cut are shown in the left part of each plot as horizontal lines. We indicate the shoulders of Haro 6-5B and Oph 163131 by an arrow. We smoothed the cuts through Tau 042021 by convolving them by a 1-D Gaussian of the beam width, to reduce the noise and make the plot clearer.

In the text

Fig. 3 Normalized averaged minor axis intensity profiles. For the most inclined sources (resp. least inclined: IRAS 04200, Haro 6-5B, and IRAS 04158), these were obtained by averaging minor axis intensity cuts over the whole major axis extent of the disk (resp. over the central 0.3′′). The light shading corresponds to the normalized rms in each band. The beam sizes in the direction of the minor axis are shown as dotted Gaussian of the corresponding color.

In the text

Fig. 4 Top row: spectral index maps, applying a filter to keep only the pixels with more emission than 5σ in both bands. The beam size is shown in the bottom left corner, along with a dark line representing a 0.5′′ scale. Middle row: spectral index cuts along the major axis. Bottom row: spectral index profiles along the minor axis, averaged as done for Fig. 3 (see text Sect. 3.2.2). For all cuts, yellow errors correspond to a flux calibration error of 10% in both bands, while orange errors are estimated from the signal-to-noise in each band. The x-axis corresponds to the offset to the center of the disk in arcseconds.

In the text

Fig. 5 Overlay of scattered light (colors, in logarithmic scales) and ALMA band 7 continuum images (5, 10, 20, and 50σ contours) for all sources in this study, except for Oph 163131 for which we show the band 6 image. The scattered light images are plotted between 3σ (except for IRAS 04302 and IRAS 04158, where we use respectively 60 and 1σ to increase the contrast) and their maximum brightness (except for IRAS 04158, for which we take a lower value to increase the contrast). The ellipse and horizontal line in the bottom left corner indicate the beam size of the ALMA image and a 0.5′′ scale. References. Scattered light images, 2.2 μm: HV Tau C: Duchêne et al. (2010); 1.9 μm: IRAS 04302: Padgett et al. (1999); 1.6 μm: Haro 6-5B: Padgett et al. (1999); 0.8 μm: Tau 042021: Duchêne et al. (2014), HH 30: Watson & Stapelfeldt (2004), IRAS 04158: Glauser et al. (2008), ESO-Hα 569: Wolff et al. (2017), HH 48 NE: Stapelfeldt et al. (2014); 0.6 μm: IRAS 04200: Stapelfeldt et al. (in prep.), Oph 163131: Stapelfeldt et al. (2014), ESO-Hα 574: Stapelfeldt et al. (2014); 0.4 μm: HK Tau B: Duchêne et al. (in prep.).

In the text

	Fig. 6 Radiative transfer models computed at 0.89 mm (band 7), for different inclinations. Top and third row: high and low mass settling models (h1mm = 0.70 au at r = 100 au), second and bottom row: high and low mass No Settling models (h1mm = hgas = 10 au at r = 100 au). The beam size and a 0.5′′ scale bar are shown in the bottom left corner of each panel.
In the text

	Fig. 7 Normalized averaged minor axis profiles for the high and low mass settling models (solid lines) and No Settling models (dashed lines) computed at 0.89 mm (band 7) and 2.06 mm (band 4), at 80° (left panel) and 90° (right panel).
In the text

	Fig. 8 Normalized major axis intensity profiles for the high mass (left) and low mass (right) models at 0.89 mm, for different inclinations between 75° and 90°.
In the text

	Fig. 9 Top: variation of the band 7 and band 4 integrated fluxes of the high mass radiative transfer models with settling (thick lines) and without settling (thin lines) as a function of the inclination. Bottom: variation of the integrated spectral index with inclination for both high and low mass radiative transfer models, with or without including settling.
In the text

Fig. 10 Left panels: ALMA data of HK Tau B. In all the figure, top row is band 7 and bottom row is band 4. Middle-left panels: millimeter images of the best model C of Stapelfeldt et al. (1998), assuming well mixed grains. They were computed for the same uv-coverage as in the data. Middle-right panels: major axis cuts of the model and the data. Right panels: averaged minor axis cuts of the data and the models, performed as in Sect. 3.2.2.

In the text

	Fig. 11 Minor axis size as a function of radius for IRAS 04302 (band 4) and Tau 042021 (band 7). The errors correspond to the standard deviation of the curves.
In the text

Fig. B.1 Observed brightness temperatures as a function of the radial distance to the central star. Solid lines profiles (band 7 and band 4 only) were computed with the same angular resolution (restored beam in Table A.1), while the light lines have the original angular resolution. We report the beam sizes in the direction of the cut as horizontal lines in the left part of each plot. The restored beam corresponds to the black line. Error bars correspond to 10% of the brightness temperature.

In the text

	Fig. B.2 Brightness temperatures as a function of the radial distance for the radiative transfer models presented in Sect. 4, computed at different inclinations. Top and third row: high and low mass settling models, second and bottom row: high and low mass No Settling models. We report the beam size in the direction of the cut in the left part of each panel.
In the text

Fig. C.1 Top row: spectral index maps of the high mass settled model, applying a filter to keep only the pixels with more emission than 5σ in both bands. The beam size is shown in the bottom left corner, along with a dark line representing a 0.5′′ scale. Middle row: spectral index cuts along the major axis. Bottom row: spectral index profiles along the minor axis, averaged along the major axis. For all cuts, yellow errors correspond to a flux calibration error of 10% in both bands, while orange errors are estimated from the signal-to-noise in each band. The x-axis corresponds to the offset to the center of the disk in arcseconds.

In the text

	Fig. D.1 Example of scattered light size estimate in Tau 042021. We represent the two spines of the nebula by the blue and red curves. The adopted radial size corresponds to the horizontal line, the maximum distance between spine points along the major axis.
In the text