Observations of edge-on protoplanetary disks with ALMA I. Results from continuum data

We analyze a sample of 12 HST-selected edge-on protoplanetary disks for which the vertical extent of the emission layers can be constrained directly. We present ALMA high angular resolution continuum images (0.1arcsec) of these disks at two wavelengths, 0.89mm and 2.06mm (respectively ALMA bands 7 and 4), supplemented with archival band 6 data (1.33mm) where available. For most sources, the millimeter continuum emission is more compact than the scattered light, both in the vertical and radial directions. Six sources are resolved along their minor axis in at least one millimeter band, providing direct information on the vertical distribution of the millimeter grains. For the second largest disk of the sample, the significant difference in vertical extent between band 7 and band 4 suggests efficient size-selective vertical settling of large grains. Furthermore, the only Class I object in our sample shows evidence of flaring in the millimeter. Along the major axis, all disks are well resolved. Four of them are larger in band 7 than in band 4 in the radial direction, and three have a similar radial extent in all bands. For all disks, we also derive the millimeter brightness temperature and spectral index maps. We find that the disks are likely optically thick and that the dust emission reveals low brightness temperatures in most cases (<10K). The integrated spectral indices are similar to those of disks at lower inclination. The comparison of a generic radiative transfer disk model with our data shows that at least 3 disks are consistent with a small millimeter dust scale height, of a few au (measured at r=100au). This is in contrast with the more classical value of h_g\sim10au derived from scattered light images and from gas line measurements. These results confirm, by direct observations, that large (millimeter) grains are subject to significant vertical settling in protoplanetary disks.


Introduction
The process of planet formation requires small submicronsized 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 Article number, page 1 of 24 arXiv:2008.06518v1 [astro-ph.SR] 14 Aug 2020 A&A proofs: manuscript no. edge-on-survey 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 particular on 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-NIR (near-infrared) 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 comparisons of 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. 2003Duchêne et al. , 2010Villenave 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 Section 2, we present the fluxes, brightness temperatures, spectral indices, and the radial and vertical extents of all disks in Section 3. Section 4 compares 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 Section 5.

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 Stapelfeldt et al., in prep) and confirmed as edge-on disks (hereafter EOD) according to their optical or near-infrared 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. 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 edgeon 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 ). 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 .
Four of the observed disks are part of multiple systems: HK Tau B, HV Tau C, Haro 6-5B, and HH 48 NE 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.

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 References. band 7 (0.89 mm) and in band 4 (2.06 mm). We present the different setups in the following section.

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 12 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 observations were 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 aver-age 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.

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.

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. 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 13 CO and C 18 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 con- figuration (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.

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.

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 selfcalibration on the sources with the highest signal-to-noise ratio (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 reimaged 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 immath CASA task. The restored beam sizes are reported in the last column of Table A.1.

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 Section 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 Notebook). 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 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 Section 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 into small 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, within elliptical 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 2 ). We note that although the phase 1 We note that we set the option fixuvw to True when applying the fixplanets task, which is similar to using the fixvis task. 2 https://almascience.eso.org/documents-andtools/latest/documents-and-tools/cycle8/alma-technicalhandbook Notes. The uncertainties are limited by errors in flux calibration and can thus be estimated by 10% of the reported brightness temperature peak. The reported peak brightness temperatures were computed using the same beam in both bands (see Table A.1). We only report the peak band 7 brightness temperature of IRAS 04158 because the disk is not detected in band 4. calibration method used for the extended configuration observations in band 7 was non-standard (see Section 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 Section 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.

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  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 correspond 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. Table 5. Position angle, inclination, major axis sizes, and estimates of the mean deconvolved major axis sharpness of the millimeter data. Notes. For each source, the errors in the millimeter sizes correspond to one tenth of the beam width in the major axis direction. When possible, the ∆r/r ratios correspond to their averaged values between the band 7 and band 4. All are deconvolved by the beam. We also indicate an estimate of the inclination inferred from the millimetric mean axis ratio. These are lower limits when the disks are not resolved along their minor axis in all bands by more than 2 beams, or when the measured inclination is too high not to be influenced by the physical vertical extent of the disk. ( †) ∆r/r of the band 7 only.

Sources
References. (a) Louvet et al. (2018) 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 higherresolution 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 ring is 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 signal-to-noise ratio, 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.2 and 0.5. The edges of these disks are only marginally shallower than the typical beam. We note that 3 out of these 4 systems are known binaries and dynamical interactions may also lead to sharp outer edges.

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

Estimation of spectral indices
The millimeter spectral index, defined as F ν ∝ ν αmm , can be used to 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 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 Section 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 signal-to-noise ratio. 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 Article number, page 9 of 24 A&A proofs: manuscript no. edge-on-survey Notes. The uncertainties on the measured minor axis sizes correspond to one tenth of the beam size along the cut direction.
( * ) Disks with inclination lower than 75 • . Their minor axis size is likely dominated by the radial extent of the disk over its vertical extent, as opposed to more inclined disks. ( †) Resolved by more than 2 beams in the minor axis direction. These are the ones for which the deconvolved minor axis size should be the most reliable.
References. 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 inside the 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    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 can not 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.

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 optical images but prefer (space-based or ground-based) near-infrared (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.

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(Pinte et al. , 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 profile prescribed 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 −4 .
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 upper limits on the dust mass derived for our sample of highly inclined disks (Section 4.3.1). We use M dust = 1 · 10 −3 M for high mass disk models (HM), and M dust = 5·10 −5 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 syn-thetic 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 Appendix B and C.

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 di- rection 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.

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 Section 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 Section 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, related to 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 Section 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 that vertical 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 can not confirm vertical settling for the other Normalized intensity

Incl (deg) S -HM NS -HM S -LM NS -LM 75
1.7 ± 0.3 1.7 ± 0.2 1.5 ± 0.4 1.5 ± 0.3 80 1.7 ± 0.2 1.2 ± 0.1 1.6 ± 0.3 1.6 ± 0.3 85 1.5 ± 0.1 0.4 ± 0.1 1.6 ± 0.2 1.2 ± 0.1 87.5 0.9 ± 0.1 0.2 ± 0.1 1.6 ± 0.2 1.0 ± 0.1 90 0.2 ± 0.1 0.2 ± 0.1 0.6 ± 0.1 1.0 ± 0.1 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(Flock et al. , 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.

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 Section 4 (which do not include radial drift), we find that opacity effects alone predict a difference of only a few percents 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 equation 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 (Section 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.

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.

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: 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 2 g −1 , 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. Notes. Dust masses are estimated from the band 7 fluxes (Table 3), assuming optically thin emission, and therefore are lower limits (see text for details).
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).

Integrated spectral indices
The median spectral index of the edge-on disks in our study, α mm = 2.5 ± 0.3, is similar to previous measurements of intermediate-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 edgeon 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 settled 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 par-tially) 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.

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 fraction of beam covered by the source. In the isothermal optically thick limit (τ ν 1), for a source filling the 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 Section 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 disks of 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 65 • . 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., Section 4.3.1).
Although the measured brightness temperatures are integrated over some vertical extent because of the beam size, 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.

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 micronsized 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. 3. 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 −5 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.5 da. This yields opacities of 4.8 cm 2 /g and 2.6 cm 2 /g (per dust mass) respectively at 0.97 mm 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 cm 2 /g and 0.026 cm 2 /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 Section 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 Section 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.

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 opacity effects are not sufficient to produce the large difference in minor axis size observed, which suggests that grain-sizedependent 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 (St 10 −2 ). For comparison, using the standard disk model presented in Riols & Lesur (2018, equation 13), at 100 au, grains of the optimal size emitting in band 7 would have Stokes numbers larger than 2.5 · 10 −2 .
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 −2 ) 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,

Data B7
Model B7  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 Section 3.2.2. Dubrulle et al. 1995) or numerical simulations including non-ideal MHD effects (Riols & Lesur 2018) provide for now the best consistency with the observations.

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-tonoise: 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.

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 millimeter 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; these ones 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 millimeter 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 proportionnal to the dust scale height, this ratio is expected for relatively large dust grains in the simple 1-D diffusion theory or numerical sim-ulations 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.   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.