The first ALMA survey of protoplanetary discs at 3 mm: demographics of grain growth in the Lupus region

We present the first ALMA survey of protoplanetary discs at 3 mm, targeting 36 young stellar objects in the Lupus star-forming region with deep observations (sensitivity 20-50 microJy/beam) at ~0.35"resolution (~50 au). Building on previous ALMA surveys at 0.89 and 1.3 mm that observed the complete sample of Class II discs in Lupus at a comparable resolution, we aim to assess the level of grain growth in the relatively young Lupus region. We measure 3 mm integrated fluxes, from which we derive disc-averaged 1-3 mm spectral indices. We find that the mean spectral index of the observed Lupus discs is $\alpha_\mathrm{1-3 mm}=2.23\pm0.06$, in all cases $\alpha_\mathrm{1-3 mm}<3.0$, with a tendency for larger spectral indices in the brightest discs and in transition discs. Furthermore, we find that the distribution of spectral indices in Lupus discs is statistically indistinguishable from that of the Taurus and Ophiuchus star-forming regions. Assuming the emission is optically thin, the low values $\alpha_\mathrm{1-3 mm}\leq 2.5$ measured for most discs can be interpreted with the presence of grains larger than 1 mm. The observations of the faint discs in the sample can be explained without invoking the presence of large grains, namely through a mixture of optically thin and optically thick emission from small grains. However, the bright (and typically large) discs do inescapably require the presence of millimeter-sized grains in order to have realistic masses. Based on a disc mass argument, our results challenge previous claims that the presence of optically thick sub-structures may be a universal explanation for the empirical millimeter size-luminosity correlation observed at 0.89 mm.


INTRODUCTION
Protoplanetary discs are the birth place of planets. In the coreaccretion scenario (Safronov 1972;Wetherill 1980), the first step to form a terrestrial planet is the growth of the typically micronsized interstellar medium (ISM) dust grains to millimeter and centimeter-sized pebbles. Although this initial growth is favoured in the dense protoplanetary disc mid-planes, evidence of growth at early stages has been obtained also in very young (Class 0 and I) ★ Contact e-mail: mtazzari@ast.cam.ac.uk systems (Miotello et al. 2014;Agurto-Gangas et al. 2019;Galametz et al. 2019). Planet formation then proceeds with the assembly of kilometer-sized planetesimals, which eventually form rocky planets and the cores of giant planets. Thanks to their sensitivity to the thermal emission of dust grains, observations at sub-millimeter and millimeter wavelengths are a key probe of the first stages of the planet formation process , and references therein).
In a smooth protoplanetary disc, millimeter and centimeter sized grains are expected to undergo fast inward radial drift due to the aerodynamic drag exerted on them by the gas in sub-Keplerian motion (Weidenschilling 1977). Since radial drift occurs very quickly compared to the disc dynamical timescale, discs should be thoroughly depleted of large grains within the first 1 Myr of their life (Brauer et al. 2008).
Observations, however, are in contrast with this scenario. As known from seminal sub-millimeter and millimeter photometric observations (e.g., Beckwith et al. 1990), and recently confirmed by ALMA spatially resolved observations (e.g., ALMA Partnership et al. 2015;Ansdell et al. 2016;Isella et al. 2016;Andrews et al. 2018b), a large fraction of protoplanetary discs are bright and relatively extended (between 30 and 100 au in radius; see Hendler et al. 2020) when observed at 1 mm (a wavelength most sensitive to the thermal emission of millimeter sized grains). The measured disc fluxes and extents at 1 mm suggest that radial drift has to be effectively slowed down or even halted in most discs (Birnstiel et al. 2010;Pinilla et al. 2012a). Thanks to their ability to trap (or, at least, decelerate the inward drift of) dust grains, local maxima in the gaseous component of a protoplanetary disc are a ready solution to the radial drift problem. This scenario is supported by the recent high angular resolution ALMA observations (e.g., Isella et al. 2016;Fedele et al. 2018;Clarke et al. 2018;Andrews et al. 2018b;Long et al. 2018Long et al. , 2019 which showed a suggestive recurrence in the discs millimetre continuum emission of orderly and axisymmetric structures at small spatial scales. Whether the rings observed in the continuum emission are effective dust traps is a topic of active research (Dullemond et al. 2018). The origin of local maxima can be linked to a variety of physical mechanisms, such as the interaction with an embedded forming planet (e.g., Pinilla et al. 2012b;Clarke et al. 2018), zonal flows (Johansen & Klahr 2005;Flock et al. 2015), and the presence of vortices (Barge & Sommeria 1995;Klahr & Henning 1997).
The optical properties of dust grains can be used to probe their spatial distribution in discs: as grains grow to sizes close to the observing wavelength, their opacity changes significantly, leaving a signature of their presence in the disc spectral energy distribution (SED). Specifically, at millimeter and centimeter wavelengths, the spectral index of the emission of optically thin dust can be directly linked, for a given temperature, to the spectral dependence of the dust opacity coefficient, which in turn depends on the maximum grain size of the emitting dust Draine 2006). Extensive observational studies investigating the level of grain growth in discs have been conducted in the Taurus-Auriga (Andrews & Williams 2005;Rodmann et al. 2006;Ricci et al. 2010b) and Ophiuchus (Andrews & Williams 2007;Ricci et al. 2012) starforming regions. Other studies, with more limited sample sizes have also targeted southern star-forming regions such as Lupus and Chamaeleon (Lommen et al. 2007;Ubach et al. 2012) and the distant Orion Nebula Cluster (Mann & Williams 2010;Ricci et al. 2011). In the majority of discs these studies found relatively small millimeter spectral indices ( 0.89−3.1 mm ∼ 2.5), which were interpreted (by means of simplified disc modelling) in terms of emission from millimeter-sized grains. It is noteworthy that all these surveys aimed at measuring the spectral index of the dust opacity absorption coefficient are by no means complete in any star-forming region, and in all cases are spatially unresolved.
Evidence of enhanced grain growth in the inner disc region in line with the expectations for radial drift was obtained by Pérez et al. (2012Pérez et al. ( , 2015 and Tazzari et al. (2016), and more recently by Tripathi et al. (2018), Huang et al. (2018Huang et al. ( , 2020, Carrasco-González et al. (2019), and Long et al. (2020) for a dozen of discs in the Taurus-Auriga and Ophiuchus region for which spatially resolved observations in a wide wavelength range (0.88 mm to 1 cm) were available. Although the self-consistent modelling of disc structure and dust emission implemented in these studies constitutes a refinement over previous works and the evidence for grain size variations is robustly demonstrated, the statistical relevance of the results is limited by the small sample of discs with homogeneous observations in the 0.89-3 mm wavelength range.
In recent years, ALMA has been used to perform extensive surveys of Class II discs in several star-forming regions such as Lupus (Ansdell et al. 2016(Ansdell et al. , 2018, Chamaeleon I (Pascucci et al. 2016), Ori (Ansdell et al. 2017), Upper Scorpius (Barenfeld et al. 2016), Ophiuchus (Cox et al. 2017;Cieza et al. 2019), Taurus (Akeson & Jensen 2014;Akeson et al. 2019;Long et al. 2018), IC 348 (Ruíz-Rodríguez et al. 2018, the Orion Nebula Cluster (Eisner et al. 2018), and NGC 2024(van Terwisga et al. 2020. Although these studies provided a new wealth of information on the discs mass, size, gas-to-dust mass ratio, and spatial brightness distribution, they only probed the disc emission at 0.89 or 1.3 mm. Extensive observations at longer wavelength are needed to assess the level of grain growth in these regions. Here we present the first ALMA survey of protoplanetary discs at 3 mm, which targeted a sample of 36 objects (38% of all Lupus Class II discs), covering the brightest 55% of the Class II discs previously detected at 0.89 mm Ansdell et al. (2016). The moderate angular resolution (∼ 0.35 ) and high sensitivity (20-50 Jy beam −1 ) allowed us to detect all the discs at 3 mm and to resolve the largest ones.
In this paper we focus on the spatially-integrated analysis of fluxes and spectral indices, discussing the implications on the grain properties and on empirical relations such as the millimeter continuum size-luminosity relation (Tripathi et al. 2017). A forthcoming paper (Tazzari et al. 2021) will perform a homogeneous analysis of the multi-wavelength observations that are available for these Lupus discs at 0.89, 1.3, and 3 mm, with a particular focus on the disc sizes.
The paper is organised as follows. In Section 2 we introduce the source sample and in Section 3 we describe observational setup and calibration details. In Section 4 we present the measured 3 mm fluxes and the inferred 0.89−3.1 mm spectral indices, with a comparison with other regions. In Section 5 we discuss the results, presenting the implications for the level of grain growth in Lupus, and discussing the new insights on the interpretation of the millimeter size-luminosity relation in light of this new 3 mm data. In Section 6 we draw our conclusions.

SAMPLE
In this study we present 3 mm observations of 36 disc-bearing young stellar objects (YSOs) of the Lupus star-forming region. The objects belong to the Lupus I to IV clouds and are classified as Class II sources (Merín et al. 2008) or have flat infrared (IR) excess measured between the 2MASS Ks (2.2 m) and Spitzer MIPS-1 (24 m) bands (Evans et al. 2009). A near-complete census of the 0.89 mm brightness of Lupus Class II discs has recently been carried out with ALMA by Ansdell et al. (2016). To build the initial sample from which we select the sources for this study, we complement the Ansdell et al. 2016 sample with the 0.89 mm ALMA observations at comparable sensitivity and resolution by MacGregor et al. (2017) for Sz 75/GQ Lup, by Cleeves et al. (2016) for Sz 82/IM Lup, and by Canovas et al. (2016) for Sz 91. From the resulting sample of 92 Class II discs (65 detected at 0.89 mm), we select those that can be detected with a signal-to-noise ratio larger than 10 when observed at 3 mm with a 0.45 resolution and a nominal sensitivity between 30 and 70 Jy beam −1 (approximately 10 and 2 minutes on source, respectively). To extrapolate the predicted 3 mm brightness, we assumed a conservative spectral index of 3.0, compatible with the spectral indices found by Lommen et al. (2007) and Ubach et al. (2012) in their Lupus samples. The requirement of a 3 mm detection with signal-to-noise ratio larger than 10 ensured that the absolute uncertainty on the 0.89-1.3 mm integrated spectral index is on average smaller than 0.12, given the sensitivity of the 0.89 mm observations. This selection resulted in total 36 discs, which have been detected at 0.89 mm with a peak brightness larger than 14 mJy beam −1 at a resolution of ∼ 0.3 : 33 discs from Ansdell et al. (2016) , with spectral types from M5.5 up to K0. The sample is complete around the Solar-mass range between 0.8 and 1.2 . The stellar masses adopted in this paper are reported in Table 1. They were determined from optical-UV spectroscopic measurements by Alcalá et al. (2017) and updated by Alcalá et al. (2019) to account for the distances drawn from Gaia DR 2 (Gaia Collaboration et al. 2018) (see Table A.1 in Alcalá et al. 2019). For all but three sources (Sz 68, RY Lup, and J16083070-3828268) the stellar masses were derived using pre-main-sequence evolutionary tracks by Baraffe et al. (2015). For Sz 68, RY Lup, and J16083070-3828268 we use stellar masses derived using tracks by Siess et al. (2000). Figure 1 summarises the properties of the 3 mm survey sample in terms of 0.89 mm integrated flux as a function of stellar mass. The figure shows all the targets of the 0.89 mm survey (Ansdell et al. 2016) except for J15450634-3417378, J16070854-3914075, and J16011549-4152351 for which a stellar mass is not available.
The figure clearly shows that the 3 mm survey is 97% complete above the 0.89 mm flux median (14.4 mJy). Overall, the 3 mm survey sample covers the brightest 55% of the Class II discs detected at 0.89 mm by Ansdell et al. (2016) and is therefore biased towards the most massive discs.
Although the spectral setup was optimised to provide optimal continuum sensitivity, it covered the HNC(1-0) spectral line (90.663 GHz) at low spectral resolution, allowing for a serendipitous detection should the line have been bright enough. Detailed analysis of the HNC measurement from this dataset is provided in Long et al. (submitted).
The array configuration used forty-three 12m antennas with baselines of 17-3150 m, corresponding to 5.3-1075 k . To optimise the overall time needed for the observations, the correlator was set to integrate for 2 minutes per source on the eighteen brightest targets, 10 minutes per source on the five faintest ones (listed above), and between 4 and 6 minutes for the targets with intermediate brightness, achieving respectively an rms noise of 50, 17, and 30 Jy beam −1 . Data calibration and imaging were performed using CASA 4.7.2 (McMullin et al. 2007). The data was calibrated using the pipeline by ESO and included flux, gain and bandpass calibrations. Flux calibration used observations of J1427-4206, bandpass calibration used observations of J1517-2422, and gain calibration used observations of J1604-4228, J1534-3526, or J1610-3958. We estimated an absolute flux calibration error of 10% based on variations in the gain calibrators.
We successfully performed self-calibration on the four brightest sources in the sample (Sz 68, Sz 82, Sz 83, Sz 98). In all cases we performed two phase-only self-calibration steps: in the first step we sought solutions across the whole scan length, and in the second step over 60 seconds. We did not find any appreciable improvement in the noise and signal-to-noise properties for phase-only self-calibration steps with smaller time intervals or amplitude self-calibration. In all cases we improved rms noise by 20% and a signal-to-noise ratio by 60%.

Continuum emission at 3 mm
Figures 2 and 3 present the continuum images synthesised from the calibrated visibilities. To produce the images we apply the CASA tclean task to the full dataset before any channel averaging. Imaging was performed using the Briggs weighting scheme (robust parameter 0.5), which yields an optimal combination of resolution and sensitivity. We used the multi-frequency synthesis deconvolver with multiple spatial scales set to a point source, the size of the synthesised beam, and 3 times the synthesised beam in order to improve image fidelity. Given the large fractional bandwidth of the data (ratio of effective bandwidth to central average frequency, ≈ 8%), we performed a multi-term clean (nterms= 2), which better accounts for spectral index variations across the image. The average synthesized beam size is 0.39 × 0.29 , corresponding to approximately 60 × 46 au at the average distance of Lupus discs (160 pc; Manara et al. 2018).
When observed at 3 mm, most discs appear regularly shaped. In a few cases, they exhibit an unusual shape owing to a peculiar morphology or viewing geometry. Discs around J16070854-3914075 and J16083070-3838268 are seen almost edge-on, thus appearing particularly elongated. In addition, they are also transition discs (van der Marel et al. 2018). RY Lup is genuinely side-lobed (see also the 1.3 mm observations in Ansdell et al. 2018) and is observed at a high inclination. Finally, transition discs around Sz 84, Sz 100, Sz 111, Sz 118, and Sz 123 also exhibit some uneven structure. MY Lup is the only transition disc in the sample that does not show a clear cavity, however this is likely due to the high inclination (larger than 70 • ; Tazzari et al. 2017) and the limited resolution of these observations.
In Table 2 we report the integrated flux measured at 3 mm and the synthesized beam size and position angle. We detect all the sources except Sz 91, which was observed with an erroneous (too large) rms noise requirement due to a typo in the observational setup. From Canovas et al. (2016) observations at 2.7 mm we expect Sz 91 to be spatially resolved at the resolution of these new 3 mm observations (∼ 0.35 ), thus it is not possible to set a simple upper limit on its 3 mm flux from this data.
The reported integrated fluxes are measured from continuum images synthesized with natural weighting (in order to maximise signal-to-noise) using circular aperture photometry. The aperture radius for each source is determined by a curve-of-growth method in which a circular aperture of increasing radius is applied until the enclosed flux becomes constant at a 3 level. We find that typically the aperture radius at which the enclosed flux plateaus tightly encloses the 3 contours, indicating that the contribution to the flux from the outer disc regions emitting just below the sensitivity of the observations is minimal. The uncertainty on the total flux is measured as the standard deviation of the flux measured within 10 apertures of same size in a region of the field of view far from the source and with no other emission. The uncertainty on the total flux reported in Table 2 does not include the systematic 10% absolute flux uncertainty. We note that Ansdell et al. (2016) used a different method to measure the integrated fluxes: they did a two-dimensional Gaussian parametric fit of the visibilities with the uvmodelfit CASA task. As a check, we have used the same method on our 3 mm observations and found integrated fluxes that agree extremely well with those measured with the circular aperture method.

The 1-3 mm spectral indices
We now compute the disc-integrated spectral indices between 0.89 and 3 mm, deferring to a forthcoming paper (Tazzari et al. 2021) the study of their spatial variations. The integrated spectral indices 0.89−3.1 mm have been computed using the 0.89 mm fluxes measured by Ansdell et al. (2016) and assuming a linear spectral slope between 0.89 and 3.1 mm: ∝ 0.89−3.1 mm . The uncertainties on the spectral indices are propagated from the flux uncertainties at the two wavelengths and include a 10% flux calibration uncertainty at both wavelengths. In a few cases we obtain values 0.89−3.1 mm < 2.0 that are however still compatible at a 2 level with 0.89−3.1 mm = 2. In only one case (Sz 74) 0.89−3.1 mm < 2 at a 3 level. A measured value of 0.89−3.1 mm < 2 could be due to the contamination from non-thermal emission processes, due to the departure of thermal emission from the Rayleigh-Jeans regime, or due to extremely high optical depth. We can exclude the first scenario since the non-thermal contribution needed at 3 mm to obtain such a low spectral index would be much larger than the typical measured upper limits (Ubach et al. 2012). We can also exclude the second scenario as these low spectral indices below 2 are obtained for discs that are very small, in which the temperature is likely to be high enough to ensure the radiation is emitted in Rayleigh Jeans regime. It is more likely that in the case of Sz 74, a close binary, the low spectral index reflects the highly optically thick contribution from the primary disc, which is expected to be truncated due to tidal interaction.
Some of the targets in our sample were observed at millimetre wavelengths by Lommen et al. (2007) and Ubach et al. (2012) using the Australia Telescope Compact Array (ATCA) at lower resolution and more limited sensitivity. Both these studies targeted Sz 68, Sz 71, and Sz 98; Ubach et al. (2012) targeted also Sz 65, Sz 66, Sz 75, RY Lup, Sz 111 and Lommen et al. (2007) targeted also Sz 82 and Sz 83. These studies obtained observations at 1-3 resolution and a sensitivity of 0.5-2 mJy beam −1 , measuring 3.3 mm fluxes that are compatible with those presented in this paper except for Sz 82, Sz 83, RY Lup, and Sz 98, for which we obtain larger 3 mm fluxes and thus smaller 0.89−3.1 mm spectral indices. Lommen et al. (2007) and Ubach et al. (2012), however, do not report the total uncertainty on their measured integrated 3 mm flux: an inspection of their deprojected visibility profiles reveals that the data is indeed very noisy, and for Sz 82, RY Lup, and Sz 98 the short baseline visibilities appear compatible with the integrated flux that we measure from our observations. Figure 4 presents the histogram of the spectral index measurements, which clearly highlights that resolved transition discs (TDs) in Lupus defined as in Ansdell et al. (2018) (see also van der Marel et al. 2018) tend to have a higher spectral index compared to the bulk of the disc population. The same trend was observed in the Lupus discs survey by Ansdell et al. (2018) in the 0.89-1.3 mm wavelength range and by Pinilla et al. (2014) in a more heterogeneous set of observations at 0.89 and 3 mm of bright discs in the Taurus, Ophiuchus, and Orion star-forming regions. Pinilla et al. (2014) reported average spectral indices of TD = 2.7 ± 0.1 for transition discs and PPD = 2.2 ± 0.1 for protoplanetary discs (PPD) with no known cavities. In this study we find similar values of TD = 2.5 ± 0.1 and PPD = 2.14 ± 0.06, respectively. For completeness we note that Ansdell et al. (2018) considered a complete survey of the Lupus Class II discs, out of which this study selected the bright end of the population. The sample used by Pinilla et al. (2014) is also biased towards the brightest discs but encompasses discs from different regions and is much more dishomogeneous in terms of observational setups. Evaluating the universality of this trend would require a detailed physical modelling of these sources and a careful consideration of the observational biases, which is beyond the scope of this paper.
There are three binaries (Sz 65+Sz 66, Sz 74, and V856 Sco/J16083427-3906181) and one triple system (Sz 68/HT Lup) in our sample. They have very different separations: 6.4 (960 au) for Sz 65+Sz 66, 1.45 (225 au) for V856 Sco (Ansdell et al. 2018), and 0.3 (45 au) for Sz 74 (Ansdell et al. 2018). We note that Sz 65+Sz 66 and V856 Sco have both components individually detected at all ALMA Bands, while for Sz 74 the angular resolution gives us only tentative evidence of the detection of the secondary component (see Fig. 3). Sz 68 is a triple system with A+B components (not resolved individually) at 0.14 (20 au) separation (Kurtovic et al. 2018), which in turn are distant 2.8 (∼ 450 au) from the C component (detected at 11 ). For Sz 74 and Sz 68 Table 2 reports the integrated flux and the spectral index of their A+B components.

Comparison to other regions
In Figure 5 we compare the spectral indices of the Lupus discs with those of other nearby star-forming regions like Taurus and Ophiuchus as a function of their integrated 1 mm flux (scaled at a distance of 140 pc). The spectral indices in Taurus and Ophiuchus are taken from a compilation by Ricci et al. (2010a,b) who fitted the spectral energy distribution in the 0.89-3.0 mm range for a sample of Class II sources with known stellar properties, without envelope contamination, and with no evidence of companion at 10-400 au scale (for more details, we refer to their papers). The star-forming regions chosen for this comparison are all relatively young, with mean ages between 1 and 3 Myr Alcalá et al. 2017). Compared to the spectral index ISM 3.7 typical of the optically thin emission of interstellar medium grains, all these regions exhibit much lower spectral indices, the bulk of the measurements lying between 1.8 and 2.5. Moreover, the similarity of the spectral indices distribution among the different regions, as highlighted by the cumulative fractions plotted in the right panel of Fig. 5, is particularly striking. Two population Anderson-Darling tests applied to the spectral index distribution of these three regions confirm that they are indistinguishable in a statistically significant way ( -value> 0.25 for the null hypothesis of spectral index measurements being drawn from the same parent distribution). It is worth noting that the lack of discs with low sub-millimeter (or millimeter) flux and high spectral indices could in principle reflect an observational bias induced by the sensitivity limitations at the longest observing wavelength, 3 mm in this case. To evaluate whether we are affected by this issue, in Figure 5 we plot the sensitivity cut imposed by the combined sensitivity of the observations presented in this paper at 3 mm and those by Ansdell et al. (2016) at 1.3 mm. Even for the smallest 1.3 mm fluxes, the sensitivity cut is comfortably far from the maximum spectral index inferred from our measurements and we can safely assume that the lack of discs in such region is genuine. Provided that the sample of Lupus discs imaged at 3 mm is complete for discs brighter than the 0.89 mm flux median (14.4 mJy), we can exclude that other Lupus discs can populate such region of the plot. Although it may seem that Lupus discs are characterised by a weak positive correlation between 0.89−3.1 mm and 1 mm , this is mostly driven by the three TDs with large spectral index and large millimeter flux. Excluding TDs from the sample no statistically significant correlation can be found. Given the low number statistics (3 discs) driving the possible positive correlation, we do not explore this further.

DISCUSSION
Spectral indices at sub-millimeter and millimeter wavelengths provide us with valuable information on the optical properties of the population of large grains residing in the discs midplane. In Sect. 5.1 we discuss the constraints that our new 3 mm measurements set on the average grain properties in the Lupus discs assuming that most of the observed emission is optically thin. In Sect. 5.3 we discuss the millimeter continuum size-luminosity relation (Tripathi et al. 2017;Andrews et al. 2018a) in light of new constraints posed by the 3 mm observations presented here.

Implications for grain growth
Thanks to their sensitivity to the thermal emission of dust grains harboured in the dense and cold disc midplane, sub-millimeter and millimeter observations constitute a powerful probe of the early phase of grain growth from sub-micron to mm sizes (Testi et al. 2014, , and references therein). At sub-millimeter and millimeter wavelengths the dust emission is mostly optically thin and the slope mm of the (sub-)mm SED (namely, the spectral index) mm = d log d log ( 1) can be related in first approximation to the dust opacity power-law slope , being the opacity ∝ , as where the further assumption that the radiation is emitted in Rayleigh-Jeans regime was made. If we consider a power-law grain size distribution ( ) ∝ − for min ≤ ≤ max ( being the radius of the emitting grain), it is possible to show (Miyake & Nakagawa 1993) that the dust power-law index depends strongly on max (provided that min < 1 m). Typical values of 1−3mm ∼ 3.7 are found for the relatively small ISM dust grains ( max ≈ 0.25 m), while Natta & Testi (2004) showed that dust grains with sizes ≤ 1 ( 1−3mm ≤ 3) can be safely interpreted as evidence of large grains ( max ≥ 1mm) for a wide range of dust composition, porosity and size distribution (see also Draine 2006).
Using the new 3 mm data presented here, we find that the mean spectral index in the Lupus discs is 0.89−3.1 mm = 2.23 ± 0.06, which corresponds to a nominal average dust opacity = 0.23±0.06 according to Eq. (2). By comparing this result with theoretical dust opacity models based on Mie theory and a range of compositions (e.g., Draine 2006;Birnstiel et al. 2018), we find that ∼ 0.2 − 0.5 requires a grain populations dominated by large grains ( ≤ 3) with maximum grain size at least larger than 1 mm.
Global dust evolution models (Brauer et al. 2007(Brauer et al. , 2008 predict very short lifetimes for large grains: as soon as they grow past mm sizes at large distances from the star, they are expected to undergo rapid inward migration due to the loss of angular momentum to the gaseous component of the disc (Weidenschilling 1977). Observationally, the loss of large grains is expected to make discs evolve very quickly towards large values 0.89−3.1 mm > 3.5 (Birnstiel et al. 2010). However, the observational evidence of low values of 0.89−3.1 mm 2.23 (hence, 0.23 and max > 1 mm) gathered not only in the Lupus region, but also in the coeval (1-3 Myr old) Taurus and Ophiuchus regions (cf. Fig. 5), can be interpreted in terms of a high occurrence of dust retention mechanisms. Pressure maxima created by strong gas inhomogeinities (Pinilla et al. 2012a) or by a forming planet embedded in the disc (Whipple 1972;Brauer et al. 2008;Pinilla et al. 2012b), streaming instability (Youdin & Goodman 2005), and dust accumulation in vortices (Barge & Sommeria 1995;Klahr & Henning 1997;Lyra et al. 2009;Barge et al. 2017) are all effects that can potentially slow down or even completely halt the drift of large grains. In this work we focused on the spatially-integrated analysis of these new 3 mm observations at moderate resolution, which do not allow us to infer which of these mechanisms are shaping the Lupus discs. However, 7 discs in our sample have been targeted by ALMA at 1.3 mm and extreme angular resolution by the DSHARP (Andrews et al. 2018b) survey (see note in Table 1): while two discs (Sz 114 and MY Lup) exhibit a rather smooth surface brightness, five of them (Sz 68, Sz 71, Sz 82, Sz 83, and Sz 129) appear as highly structured, lending support to the scenario in which the large grains are mostly concentrated in long-lived disc structures.
A further argument supporting the presence of large grains is that even TDs have relatively small 0.89−3.1 mm 2.5 values. In these cavity-bearing discs, the spectral index is determined by the emission of the outer disc, i.e. mostly free from the contamination of a possibly optically thick inner disc. Even considering the caveat that TDs might still have very small inner discs and (or) very narrow outer rings, the fact that they exhibit 0.89−3.1 mm 2.5 suggests that optical depth is genuinely small and the inferred face-value 0.5 can be robustly interpreted with the presence of large grains.

Caveats
A natural caveat of interpreting the low spectral index values in terms of mm-sized grains is that an increased optical depth or any deviations from Rayleigh-Jeans regime would naturally alter the simple relation in Eq. (2) between 0.89−3.1 mm and , with both effects tending to bias 0.89−3.1 mm towards lower values (e.g., see the discussion in Huang et al. 2018). Tazzari et al. (2017) analysed the 0.89 mm ALMA observations of most of the discs in this sample with a two-layer disc model (Chiang & Goldreich 1997;Dullemond & Dominik 2004) finding that most discs should be emitting in the Rayleigh-Jeans regime. However, we cannot exclude that for large discs, the emission from the outer and colder parts of the disc with temperatures approaching 7 K could contribute a non negligible amount of emission. We note that the largest disc in our sample, Sz 82/IM Lup, has a relatively large spectral index 0.89−3.1 mm = 2.67 ± 0.16, suggesting that the contribution of non-Rayleigh-Jeans emission to its observed flux is small.
It is worth highlighting that in this analysis we have considered only the absorption component of dust opacity. Although at low optical depths this is essentially accurate, if the optical depth is high and large grains are present, dust scattering is expected to be an important contributor to the total millimeter opacity (Liu 2019; Zhu The disc-integrated analysis that we performed in this study does not allow us to robustly quantify the optical depth in these discs and we cannot therefore rule out this alternative explanation for the low spectral indices measured in these Lupus discs. It is worth mentioning that if optical depth will turn out to be small, then the above results based on absorption remain valid. A forthcoming spatially-resolved study of the multiwavelength observations gathered for the Lupus discs will allow us to put firmer constraints on the radial variations of the optical depth, on its dependency on frequency, and therefore on the dust properties in the Lupus discs (Tazzari et al. 2021). A third caveat is represented by the observational sample (Sect. 2), which includes only the brightest Lupus discs. Being the sub-millimeter and millimeter flux a rough proxy for the mass of the emitting dust, the sample is clearly biased towards the most massive discs, which can also be preferentially characterised by enhanced optical depth at small radii w.r.t. the bulk disc population. Extending the present 3 mm survey to the fainter disc population is the necessary next step to quantify this potential bias.

Disc dust masses
At millimeter wavelengths the disc emission is typically optically thin, allowing us to infer the disc dust mass from the integrated flux given some assumptions on the opacity and the temperature of the emitting dust (Beckwith et al. 1990). Despite admittedly simplified, this approach has been extensively used to infer disc masses from sub-millimeter and millimeter photometric surveys (see, e.g., Andrews & Williams 2005;Andrews et al. 2013), and more recently from ALMA continuum surveys (Ansdell et al. 2016(Ansdell et al. , 2017Pascucci et al. 2016;Barenfeld et al. 2016;Long et al. 2018;Cieza et al. 2019).
We compute disc dust masses ( dust ) from the measured 3 mm integrated fluxes as Hildebrand (1983): where is the distance inferred from Gaia DR2 measurements (Bailer-Jones et al. 2018), is the dust opacity and dust is the dust temperature. To ease the comparison with other measurements in literature (especially with the previous Lupus surveys at 0.89 and 1.3 mm), we adopt a power-law opacity = 3.37 ( /337 GHz) cm 2 g −1 , with measured for each disc as = 0.89−3.1 mm − 2 and a normalisation such that the opacity at 0.89 mm (337 GHz) matches that used by Ansdell et al. (2016). We further assume that the dust is isothermal, with dust = 20 K. The latter assumption is justified since the bulk of the mm emission is expected to originate in the nearly isothermal outer disc regions.
The dust masses that we infer from the 3 mm fluxes are reported in Table 2. In order to reduce the biases in comparing these masses with those inferred from previous 0.89 mm (Ansdell et al. 2016) and 1.3 mm (Ansdell et al. 2018) observations of the same Lupus discs, we re-compute the latter ones using Eq. (3) and up-to-date distances computed using Gaia DR2 measurements (Bailer-Jones et al. 2018). In general, the masses inferred from the 3 mm fluxes are smaller than those inferred from the 0.89 and 1.3 mm ones. Median dust values for this sample of Lupus discs inferred at 3, 1.3, and 0.89 mm are 22 ± 2, 24 ± 2, and 27 ± 2 Earth masses, respectively. Under the assumption of perfectly optically thin emission, this result would imply that the grains contributing to the 3 mm emission are on average ∼ 20% less abundant than the (roughly 3 times) smaller grains emitting at 0.89 mm. We shall note, however, that these dust masses are highly sensitive to the dust opacity spectral index (i.e., on 0.89−3.1 mm ): if some of the emission is optically thick (e.g., emission arising from narrow rings unresolved in the current observations), the value to be used in Eq. 3 would be larger, which would result in a smaller opacity at 3 mm and therefore in larger dust masses. Multi-wavelength observations at matching high angular resolution are needed in order to resolve the disc structure and break this degeneracy and improve the estimates on the dust mass.

Millimeter continuum size-luminosity relation: new insights on its interpretation
In Section 5.1 we have shown that, under the assumption of optically thin emission, the new 3 mm observations of the brightest 35 Lupus discs can be readily explained in terms of fluxes, spectral indices, and reasonable dust masses with the presence of large millimetresized grains. An alternative explanation that does not require the presence of large grains has been extensively discussed (e.g., Ricci et al. 2012) and posits that the disc emission has a non-negligible optically thick contribution. Lending support to this latter scenario is the correlation between the discs continuum size (the radius enclosing 68% of the disc flux) and their luminosity mm ∝ 2 eff (Tripathi et al. 2017;Andrews et al. 2018a) found in a sample of sub-arcsecond resolution observations at 0.89 mm of 50 nearby protoplanetary discs in the Taurus-Auriga and Ophiuchus star-forming region. The presence of such a relation has recently been confirmed also using ALMA observations (Tazzari et al. 2017;Andrews et al. 2018a;Hendler et al. 2020;Sanchis et al. 2020). Simple dust evolution models (Birnstiel et al. 2012) are able to reproduce the observed trend for reasonable initial conditions by assuming that grain sizes are set by radial drift (see also Rosotti et al. 2019). However, Tripathi et al. (2017) suggested that the observed relation can be alternatively explained if the discs are characterised by optically thick emission in narrow, spatially unresolved, regions, with filling factors (ratio of area covered by optically thick region to total disc area) of a few tens of percent. This would imply that the low spectral indices 0.89−3.1 mm are determined by these optically thick regions rather than by the presence of genuinely large grains.
Since Tripathi et al. (2017) has only considered observations at 0.89 mm, we aim to re-evaluate this latter conclusion in light of the new 3 mm observations presented here. We focus our analysis on the IQ Tau disc which is discussed in detail in Tripathi et al. 2017 and it is claimed that its emission can be explained with a 10-30% optically thick emission within 80 au. For the purpose of this simple modelling, let us assume that the over-density (hence, the optically thick emission) is concentrated at the core of the disc rather than being distributed in small-scale structures throughout the disc extent: we write the dust surface density as Σ dust ( ) = Σ 0 for ≤ core and Σ dust = 0 elsewhere. Although this may look less realistic than a myriad of localised rings emitting highly optically thick emission, the quantity of interest here is the total amount of optically thick contribution (and the mass it may be hidden in it) rather than its spatial distribution. The total disc mass is: where is the optically thick filling factor (0 ≤ ≤ 1), is the gas-to-dust mass ratio, and we assumed that in core . If we assume the same opacity 0.89 mm = 3.5 cm 2 g −1 used by Tripathi et al. (2017), a gas-to-dust mass ratio = 100, and we take IQ Tau's effective radius as the core radius ( core = eff = 80 au), the disc mass would be tot ∼ 6 × 10 −3 , which we derived assuming a small filling factor = 10% and the minimum surface density (Σ 0 = 0.285 g cm −2 ) required to have an 0.89 mm optically thick emission ( = Σ dust = 1). Considering that IQ Tau has an integrated flux of 170 mJy, the median of the Tripathi et al. (2017) sample, we conclude that half of the discs in that sample (which roughly includes the 50 brightest discs in the Taurus-Auriga and Ophiuchus regions) would have larger optically thick regions and, therefore, larger disc mass than IQ Tau. For reference, for a disc with an integrated flux at 1 mm of 1 Jy, the size-luminosity relation predicts an effective size of 160 au which would imply a total disc mass of tot ∼ 0.025 for = 10%. These disc masses are not unrealistically large at face value. However, they only represent the lower limit implied by the Tripathi et al. (2017) argument: larger filling factors > 10% or larger optical depths > 1 would both require more massive discs.
Crucially, the additional information that we now have from our 3 mm survey in terms of spectral indices suggest that the disc masses should be much larger than the nominal values that we have just derived, making the universal explanation of the sizeluminosity relation in terms of optically thick substructures less straightforward. Indeed, the disc masses that we derived according to the arguments in Tripathi et al. (2017) are the smallest compatible with their results, obtained assuming that the disc emission is marginally optically thick at 0.89 mm ( = 1). Assuming the disc is made of small grains ( max 1 mm, ∼ 2), if its core is marginally thick at 0.89 mm, it would naturally become optically thin at longer wavelengths ( ∝ ) and the disc would inescapably exhibit a large spectral index 0.89−3.1 mm 2. In order to reproduce not only the disc integrated flux at 1 mm, but also the typically observed 0.89−3.1 mm ≤ 2.5, the core needs to be optically thick also at 1.3 and 3 mm, thus implying an optical depth at 0.89 mm, and thus disc masses, 5-10 times larger ( 0.89 mm / 3 mm = (3/0.89) ≈ 7 for = 1.7) than suggested by the Tripathi et al. (2017) argument. We conclude that although the presence of optically thick regions can explain the size-luminosity relation on brightness grounds, it is much less viable in terms of required disc mass if the grains are small (i.e., with ∼ 2). Indeed, while the 0.89−3.1 mm values measured in the fainter fraction of discs can be well explained without the need of large grains, the low 0.89−3.1 mm values measured in the bright and large discs require the genuine emission of large grains in order to yield realistic (i.e., not leading to gravitational instability) disc masses.
We emphasise that, although this argument does not disfavour the presence of optically thick substructures in general, it does show that they are not compatible with the observations (spatiallyintegrated mm fluxes and spectral indices) if the dust component is mostly made of small grains ( ∼ 2). The presence of widespread optically thick substructures made of large mm-sized grains (with large mm albedo, see Zhu et al. 2019) still remains a viable scenario that will have to be investigated with multi-wavelength observations at higher resolution.

CONCLUSIONS
In this paper we have presented the first ALMA survey of protoplanetary discs at 3 mm, targeting the 36 brightest Class II protoplanetary discs in the Lupus star-forming region. The main results can be summarised as follows: (1) We obtained 3 mm ALMA observations at ∼ 0.35 resolution at a sensitivity between 20 and 50 Jy, which detected 35 out of 36 discs (Sz 91 was not detected due to an erroneous observational setup).
(2) By combining the new 3 mm observations with previous ALMA observations at 0.89 mm and a similar angular resolution, we find that all Lupus discs have spectral indices 0.89−3.1 mm < 3.0, with a tendency of larger values for transition discs. The mean spectral index for the entire Lupus sample is 0.89−3.1 mm 2.23 ± 0.06, while for the transition discs is TD 2.5 ± 0.1.
(3) Under the assumption of optically thin and Rayleigh-Jeans emission, the low spectral indices can be interpreted as evidence of large grains, with a median dust opacity power-law index 0.23 ± 0.06, which require grains with max > 1 mm for a range of reasonable dust composition and porosity.
(4) We find that the distribution of spectral indices in Lupus is statistically indistinguishable from that of the Taurus and Ophiuchus star-forming regions, suggesting that dust retention mechanisms are common in discs from their early stages of evolution. (5) The mean disc dust mass that we obtain from the 3 mm fluxes is 22 ± 2 ⊕ , smaller by ∼20% than those obtained from previous 0.89 mm observations, but in line with those obtained from 1.3 mm observations. (6) We revisit the claim that the millimeter continuum sizeluminosity relation can be explained by the widespread presence of localised substructures that emit optically thick radiation (Tripathi et al. 2017). In light of the new 3 mm measurements presented in this study we argue that the disc masses implied by such scenario would be possible, albeit very high if the grains were small.
In this paper we have focused on the spatially-integrated fluxes and spectral indices. A forthcoming paper will present the spatiallyresolved study of the multi-wavelength Lupus disc observations (Tazzari et al. 2021).