Multi-wavelength continuum sizes of protoplanetary discs: scaling relations and implications for grain growth and radial drift

We analyse spatially resolved ALMA observations at 0.9, 1.3, and 3.1 mm for the 26 brightest protoplanetary discs in the Lupus star-forming region. We characterise the discs multi-wavelength brightness profiles by fitting the interferometric visibilities in a homogeneous way, obtaining effective disc sizes at the three wavelengths, spectral index profiles and optical depth estimates. We report three fundamental discoveries: first, the millimeter continuum size - luminosity relation already observed at 0.9 mm is also present at 1.3 mm with an identical slope, and at 3.1 mm with a steeper slope, confirming that emission at longer wavelengths becomes increasingly optically thin. Second, when observed at 3.1 mm the discs appear to be only 9% smaller than when observed at 0.9 mm, in tension with models of dust evolution which predict a starker difference. Third, by forward modelling the sample of measurements with a simple parametric disc model, we find that the presence of large grains ($a_\mathrm{max}>1 $mm) throughout the discs is the most favoured explanation for all discs as it reproduces simultaneously their spectral indices, optical depth, luminosity, and radial extent in the 0.9-1.3 mm wavelength range. We also find that the observations can be alternatively interpreted with the discs being dominated by optically thick, unresolved, substructures made of mm-sized grains with a high scattering albedo.


INTRODUCTION
Multi-wavelength observations of protoplanetary discs at millimeter wavelengths have long been used to place constraints on the typical sizes of dust grains in discs (Testi et al. 2014, and references therein). This is important on several grounds. Firstly, if the grain size distribution is observationally calibrated then it improves estimates of the disc opacity; provided the emission is optically thin, it improves estimates of the mass of solid materials in discs, a quantity of obvious interest when assessing the potential of discs to form rocky planets. Secondly, grain size determines the strength of coupling between the dust and gas phases; for typical gas densities and radii in the disc, grains of size around a mm are subject to strong radial drift towards the star (Weidenschilling 1977). Thus measuring the grain size allows one to assess how tightly the gas ★ Contact e-mail: mtazzari@ast.cam.ac.uk and dust dynamics are coupled and thus to determine whether the dust distribution is likely to be a good indicator of the underlying gas profile.
Nevertheless, the information that can be gathered from unresolved multi-wavelength observations is limited for a number of reasons. Not only is there the obvious problem of trying to derive mean characteristics from emission arising from a range of spatial locations, but, without knowing the actual spatial extent of the emission, it is impossible to assess its optical depth. Ideally one would want to obtain well resolved radial profiles of disc emission at multiple mm wavelengths and hence derive local radial profiles of spectral index. This has however been performed only on few bright discs (AS 209 by Pérez et al. 2012; DoAr 25 and CY Tau by Pérez et al. 2015; FT Tau and DR Tau by Tazzari et al. 2016; UZ Tau E by Tripathi et al. 2018; TW Hya by Huang et al. 2018 andMacias et al. 2021; HL Tau by Carrasco-González et al. 2019; GM Aur by Huang et al. 2020; DS Tau, GO Tau, and DL Tau by Long et al. 2020), which makes it hard to derive statistical properties of disc populations.
An intermediate situation with regard to sample size and resolution is attained in the case of our ALMA surveys of the Lupus protoplanetary discs, for which we have gathered observations at intermediate resolution (0.25-0.35 ) across the 0.9-3.1 mm wavelength range for 36 objects. In Tazzari et al. (2021) we presented the 3 mm survey (with particular focus on the spatially-integrated fluxes and spectral indices), which complements the 0.9 mm (Ansdell et al. 2016) and 1.3 mm (Ansdell et al. 2018) surveys. Lupus is so far the only star forming region that has been targeted with ALMA survey-like observations at three wavelengths: the homogeneity of the data and the statistical relevance of the targeted sample make it the ideal benchmark to test our understanding of dust evolution and disc structure.
With an average distance from us of 160 pc (Manara et al. 2018), Lupus discs are partially resolved with typically 4-6 resolution elements across the disc diameter. While it is possible to derive spectral index profiles from such data, this typically extends over a limited radial range and lacks information on small scale structure. Notwithstanding, a quantity that can be robustly derived from even moderate resolution data is the disc radius (practically, the radius enclosing 68% of the flux, 68 ) at different wavelengths. The ratio of such radii at different wavelengths turns out to be a rather constraining quantity. For example, in the case of a disc where the maximum grain size decreases rapidly with radius, the less efficient radiation of small grains at long wavelengths would lead to the disc size being smaller as the wavelength of observation increases. On the other hand, if the entire disc is optically thick, or if the maximum grain size is large enough at all radii so that the disc emits efficiently at all the wavebands considered, the disc size would be expected to depend more weakly on radius.
Another quantity that can be assessed with data that is moderately resolved, as in the case of Lupus discs, is how close the disc is to being optically thick. Here we define the optically thick fraction as the ratio of the observed flux within 68 to the flux of a disc of size 68 where the emission was completely optically thick (having made an estimate of the likely temperature profile of the disc according to the stellar luminosity). We find that, in conjunction, the spatially averaged spectral indices (Tazzari et al. 2021), the ratio of disc sizes in the three ALMA wavebands, and the optically thick fraction provides information that is quite constraining of the required disc properties, even in cases where there is very little direct information about the radial profiles of disc emission.
In this paper we analyse spatially resolved ALMA observations at 0.9, 1.3, and 3.1 mm for the 26 brightest protoplanetary discs in the Lupus star-forming region. We characterise the discs multi-wavelength brightness profiles by fitting the interferometric visibilities in a homogeneous way, obtaining effective disc sizes at the three wavelengths, spectral index profiles and optical depth estimates.
The paper is structured as follows. In Section 2 we present the sample of young stellar objects that we consider. Section 3 describes how the visibility modelling is performed and Section 4 presents a worked example for one disc. Section 5 presents the results of the modelling in terms of demographic properties (disc radii, optical depth, continuum size-luminosity scaling relation). In Section 6 we employ simple toy models of the disc emission to interpret the sample of measurements in terms of disc structure and dust properties, discussing the implications for radial drift and grain growth. Finally, in Section 7 we draw our conclusions. Appendix A presents the linear regression of the millimeter size-integrated flux relation. Appendices B and C show detailed properties of some toy models representative of different regimes. Appendix D discusses the effects of scattering on the inferred dust properties. Appendix E reports the detailed multi-wavelength fit results for all the discs.

SAMPLE AND OBSERVATIONS
Protoplanetary discs in the Lupus star forming region have been recently targeted with extensive ALMA surveys at moderate sensitivity and resolution (0.25 -0.35 ) that provided the first systematic view on the structure of these discs. Here we assemble all the data into a multi-wavelength dataset collecting 0.9 mm observations (ALMA Band 7) from Ansdell et al. (2016), 1.3 mm observations (ALMA Band 6) from Ansdell et al. (2018), and 3.1 mm observations (ALMA Band 3) from Tazzari et al. (2020) for all the discs in common to these surveys. The targets originally selected by these ALMA surveys are sources typically classified as Class II discs (Merín et al. 2008) or with a flat infrared excess measured between the 2MASS Ks (2.2 m) and Spitzer MIPS-1 (24 m) bands (Evans et al. 2009). Please refer to the survey papers for full details on the sample selection. Note that IM Lup (Sz 82) was not targeted at 0.9 mm by Ansdell et al. (2016): we thus use the 0.9 mm observations at similar sensitivity and resolution by Cleeves et al. (2016).
The combined sample of these three surveys results in 35 sources targeted at 0.9, 1.3, and 3.1 mm. This constitutes the initial sample for the present study. Since Sz 91 has not been detected at 3.1 mm (Tazzari et al. 2021) we remove it from the sample. Moreover, we want to consider only discs that orbit around single stars or wide orbit binary companions, in order to avoid biases in the disc properties induced by the presence of the companion (e.g., tidal truncation effects; see Manara et al. 2019;Akeson et al. 2019). Since the interferometric visibilities contain contributions from all the sources in the field of view, analysing discs in multiple systems requires fitting all their components simultaneously (e.g., Manara et al. 2019), which is outside the scope of this paper. We thus exclude two close binaries (Sz 74, and V856 Sco/Lupus III 53) and a triple system (Sz 68/HT Lup) with separations less than 2 . The resulting sample is thus made of 31 sources: their coordinates and properties are summarised in Table 1. We note that we will remove 5 of these sources from the sample for different reasons (see Section 3.3), leaving 26 sources for the full multi-wavelength analysis.

Modelling the disc brightness profiles
To characterise the surface brightness profile of all discs we analyse the ALMA observations at all three observing wavelengths (0.89, 1.3, and 3.1 mm) in a homogeneous way. We fit the complex visibilities with a parametric brightness profile under the assumption that the disc emission can be represented well by an axisymmetric profile. At the resolution of these observations this assumption is typically well satisfied a posteriori, as the fit residuals confirm in all but a few cases (cf. detailed fit results in Appendix E). Inferring the brightness profile by fitting the visibilities (as opposed to extracting it from the corresponding synthesized image) allows for a more robust understanding of the uncertainties and for a spatial resolution typically better than the nominal synthesized beam.
Past studies adopted different parametrisations to fit protoplanetary discs brightness profiles: a sharply truncated power law (e.g.,  Alcalá et al. (2017) from spectroscopic observations and corrected for Gaia DR2 distances by Alcalá et al. (2019).
(1) Not found in Gaia DR2. For this source, we assume the average distance of the Lupus cloud complex (160 pc, Manara et al. 2018).
(3) Sources that will be removed from the multi-wavelength analysis for different reasons (see Sect. 3.3). A machine-readable version of this table is available online (see the Data Availability statement). Isella et al. 2009;Guilloteau et al. 2011), the self-similar solution to the viscous evolution equation by Lynden-Bell & Pringle 1974 (e.g., Andrews et al. 2009;Isella et al. 2009;Guilloteau et al. 2011;Tazzari et al. 2017), and more recently the profile introduced by Lauer et al. (1995) (e.g., Tripathi et al. 2017Andrews et al. 2018a;Hendler et al. 2020). These three parametrisations have an increasing degree of flexibility, with the Nuker profile being also the most general one (see Tripathi et al. 2017 for a discussion of its properties). The additional flexibility comes with an increased number of parameters (5, compared to 3 for the former two profiles) that allow for a decoupled description of the inner and outer disc region. The studies that employed the Nuker profile analysed disc observations with resolution and sensitivity comparable or identical to those of the observations we analyse here and they all show that at least one of the parameters ( , which controls the sharpness of the transition between inner and outer disc) is hardly constrained by them. In this study we adopt a modified version of the self-similar profile that has 4 parameters and, although lacks Nuker's capability to reproduce sharp broken power-law profiles, is still very flexible and retains the advantageous decoupling between inner and outer disc description: where ( ) is the scaling radius, 1 is the power-law slope in the inner disc, and 2 controls the slope in the exponentially tapered outer disc. Note that if 2 = 2 − 1 , then Eq. (1) corresponds to the canonical self-similar solution and in the case of 1 = 0, 2 = 2 it reduces to the Gaussian. This profile has already been used successfully by Long et al. (2019); Manara et al. (2019). The normalisation of ( ) is determined so that the integrated flux (mJy) matches the observations. We define the cumulative flux where the inclination and brightness ( ) have been inferred with the MCMC visibility modelling. By definition the integrated disc flux is = (∞). The functional form in Eq. (1) has 4 free parameters: , , angular offsets (ΔR.A.,ΔDec.) from the phase center. We thus have an 8-dimensional parameter space = ( , , 1 , 2 , , , ΔR.A., ΔDec.) , described by 4 parameters for the brightness profile plus 4 for the system geometry. We compare a model to the observations assuming a Gaussian likelihood where obs are the observed complex visibilities, mod ( ) are the model visibilities, is the total number of visibility points and is the weight 1 of the −th visibility. We denote j = ( , ) where ( , ) is the Fourier-plane coordinate of the −th visibility point.
For each model we compute the visibilities mod using the GALARIO Python package 2 (Tazzari et al. 2018), which first computes the 2D image of the disc for a given brightness ( ) and then samples its Fourier transform at the observed ( , ) locations.
We explore the parameter space with a Bayesian approach using an affine-invariant Markov chain Monte Carlo (MCMC) ensemble sampler implemented in the Python package v2.2.1 (Foreman-Mackey et al. 2013). We thus obtain samples of the posterior probability of the model parameters given the observations: where ( ) is the prior probability on the parameters and is a normalisation constant that can be neglected for the purposes of this study. Since the parameters are independent, the priors can be factored as ( ) = ( ). We adopt uniform priors ( , ) between and for all parameters, (log /mJy) = (−3, 5), For all fits we use 48 walkers (i.e., 6 chains for each free parameter) that we initially draw from the same uniform distributions that we use as priors. We then evolve the walkers for 10 5 steps, with convergence typically being achieved in the first 15-20 thousand steps (burn-in phase). To obtain independent samples of the posterior, we remove the burn-in steps from the entire 10 5 -steps long chain and we further thin the remaining samples by selecting only steps separated by one autocorrelation time (typically in the range 100-150 steps for all discs).
To prepare the observational data for the fits, we use the CASA split command to average the visibilities in spectrum (reducing the data to 1 channel per spectral window), and in time (over 30 seconds). By applying the CASA statwt command we also make sure that the absolute scale of the visibility weights is correct.

Measuring the disc size
The disc size is not quantified directly by the modelling presented above. Historically, studies that followed a similar approach of fitting the interferometric visibilities adopted the scaling radius associated 1 Theoretical visibility weights are computed by the CASA software package assuming Gaussian uncertainties (Wrobel & Walker 1999). 2 Code available at https://github.com/mtazzari/galario. with their parametrisation as the disc size. In the case of the functional form that we choose here (Eq. 1) this would translate into taking as the disc size. As pointed out by Tripathi et al. (2017), the choice of a punctual quantity such as a scaling radius is problematic for several reasons. A more robust way of measuring the disc size is through the definition of an effective radius eff enclosing a fraction of the total disc emission (such that ( eff ) = ), which accounts for the limited sensitivity and resolution of the data and is much less sensitive to the degeneracies between the parameters. We note that with this definition of disc size, the particular functional form chosen to fit the brightness profile becomes less important: any parametrisation that reproduces the data well will suggest the same effective size (Andrews et al. 2018a).
To ease the comparison with other studies, and given its immediate connection to a Gaussian standard deviation in the case of poorly resolved observations, throughout this study we use = 0.68: we will refer to the effective angular radius with 68 and with 68 = 68 × to its linear counterpart ( being the Gaia DR2 distance, see Table 1). Similarly to Tripathi et al. (2017), we report that choosing in the range between 0.5 and 0.8 has a negligible impact on the general conclusions of this work.

Modified modelling for special cases
There are 7 discs (J15450887-3417333, Sz 69, Sz 72, Sz 73, Sz 90, Sz 108B, Sz 110, J16085324-3914401, J16124373-3815031) for which the MCMC does not converge at one or more wavelengths. They are the faintest objects in the sample across the three wavelengths, with <= 3 mJy at 3.1mm. They appear unresolved in the CLEAN-synthesized images (Briggs weighting, robust 0.5) but the deprojected visibilities show that they are resolved, albeit with a larger uncertainty given the noisier data. The MCMC posteriors for the modified self-similar profile in Eq. (1) are highly degenerate and tend to prefer very steep profiles ( 1 < −2, 2 > 5), with extremely large scaling radii 5 ). Although it is possible that the underlying brightness of these sources is genuinely steeper than other discs, we found that equally good, and sometimes even better, solutions were found by adopting a Gaussian profile centered on the origin: where the normalisation (such that the integrated flux matches the observations) and the width ( ) are the only free parameters in the MCMC. In these cases we adopt uniform priors (log( /mJy) = (−3, 5) and ( ) = (0, 5 ).
There are 4 discs (J15450634-3417378, Sz 73, Sz 110, Sz 113) with very noisy observations for which neither the modified self-similar nor the Gaussian profiles converge at one or more wavelengths. We thus exclude them from the analysis.
In our sample there is one binary system that has both the components detected at the three wavelengths: the wide binary Sz 65+Sz 66 (separation 6.4 ). Sz 65+Sz 66 were targeted at all three wavelengths with two distinct scheduling blocks, each centered on one of the components. Given their large separation, Sz 65 and Sz 66 fall in the respective primary beam of both the scheduling blocks. For Sz 65 we fit the visibility data from the scheduling block centered on it, after the Sz 66 contribution to the visibilities was removed with CASA by subtracting the Fourier transform of its CLEAN components. Since Sz 66 is considerably fainter than Sz 65, the removal of Sz 65 emission from its scheduling block proved difficult. For this reason we do not report fit results for Sz 66. Out of the sample of 31 discs initially selected for the multi-wavelength analysis (Sect. 2), we are therefore left with 26 of them.

RESULTS
We fit the ALMA observations of 26 discs at 0.89, 1.3, and 3.1 mm using the modified self-similar parametrisation in Eq. (1) for 19 of them, and the Gaussian parametrisation in Eq. (6) for 7 of them. Table 2 summarises the statistics of the posterior distributions inferred for the free parameters of the modified self-similar brightness profile. Table 3 presents analogous results for the Gaussian fits. The tables provide also the posteriors of derived quantities such as the integrated flux , the effective (angular) radius 68 , the effective linear radius 68 , and the optically thick fraction F as defined in Sect. 5.2.
As an example, here we present the results obtained for Sz 83, and we collate the results for all the other discs in Appendix E. For each disc, the fits at 0.9, 1.3, and 3.1 mm are performed independently. For each wavelength, once the MCMC has converged, we build the posterior distribution of the brightness profile by drawing 200 random samples from the posterior of the free parameters. The top panel in Figure 1 shows the median brightness profiles that we obtain for Sz 83. To ease the comparison among the profiles at multiple wavelengths, we normalise them to 1.3mm ( 68 ), i.e. the brightness inferred at 1.3 mm at the location of the 1.3 mm effective radius. The location of 68 at each wavelength is marked as an empty circle on the brightness profiles, which are plotted for ≤ 95 . The bottom panel presents the radial profile of the spectral index between 0.9 and 3.1 mm and between 1.3 and 3.1 mm, derived from the median profiles shown in the top panel.
We note that these ALMA interferometric observations are affected by an absolute flux calibration uncertainty of 5% at Band 3 (Tazzari et al. 2021), 10% at Band 6 (Ansdell et al. 2018), and 10% for Band 7 (Ansdell et al. 2016). This reflects in a systematic uncertainty on the spectral index profiles of about ±0.09 for 0.9-3.1 mm and ±0.13 for 1.3-3.1 mm. The shaded areas around the median spectral index profiles visualise this systematic uncertainty. Figure 2 shows the comparison between the observed visibilities and those computed for the bestfit model (i.e., the median brightness). For each wavelength, the visibilities are normalised to the disc integrated flux , binned in 30k intervals, and deprojected according to the inferred , values. This plot provides a useful benchmark on the quality of the fits as it compares directly the fitted data with the inferred model. In the case of Sz 83, as for most of the discs in the sample, we obtain an excellent agreement between the model and the data. In many cases our results highlight a striking similarity between the profiles at different wavelengths.
As a further check on the goodness of the fits, in Figure 3 we present the synthesized images of the observed, bestfit model, and residual visibilities at 0.9 mm, 1.3 mm, and 3.1 mm. The images have been produced with the tclean CASA command using Briggs weighting and robust parameter 0.5. In the case of Sz 83, the excellent agreement of the visibility profiles manifests in negligible (<3 ) residuals at all the wavelengths.

Disc radii: dependence on frequency
For each disc, we derive the effective radius at 0.9, 1.3, and 3.1 mm from the inferred brightness posterior. Figure   Note The parameter values are the medians of their posterior distributions, and the uncertainties represent the 68% confidence interval. The posterior does not include systematic flux calibration uncertainty. Note that the parameters to the right of the vertical bar ( 68 , 68 , F) are not free parameters in the fit, but they are rather inferred from the joint posterior on { , , 1 , 2 }. Detailed fit results are presented in Appendix E1. A machine-readable version of this table is available online (see the Data Availability statement). wavelength effective disc radii as a function of frequency. To better visualise changes across wavelengths, the radii are normalised to their values at 0.9 mm. The most striking finding is that the majority of discs have a very similar radius across the 0.9-3.1 mm wavelength range. Out of 26 discs, 21 of them have a 3.1 mm radius that differs less than 20% from their 0.9 mm radius. Notably, the large Sz 82/IM Lup disc is the one with the largest difference: its 3.1 mm radius is 70% the 0.9 mm radius. Table 4 summarises the statistics on the multi-wavelength radii. Compared to the radius measured at 0.9 mm, the mean 1.3 mm radius is 4% smaller and the mean 3.1 mm radius ∼9% smaller, with both measurements being compatible with the 0.9 mm one within 1 . The mean of the size-frequency slopes for the whole sample is: namely compatible with a flat size-frequency relation. We note that there is a group of discs with a measured 1.3 mm radius that appears larger than that measured at 0.9 mm by more than 10%: in most cases they correspond to the Gaussian fits (represented in Figure 4 with smaller filled circles), which we employed for the discs with noisier observations. This group of discs has a flat sizefrequency relation: The 1.3 mm radii measurements for these discs are compatible with the 0.9 mm radii within 1 and are likely due to the fainter nature of their emission. We highlight that for many discs with high signal-to-noise observations (which typically have been fitted with the modified self-similar profile), the radius is essentially constant across wavelengths. This can be seen even before modelling, just by comparing the visibility profiles, which almost perfectly overlap in many cases. The mean slope only for the discs that have been fitted with the modified self-similar brightness profile ( Figure 4. Disc effective radius as a function of observing frequency, normalised to the value at 0.9 mm. The discs fitted with modified self-similar profile (Table 2) are shown as large filled circles. The discs fitted with Gaussian profile (Table 3) are shown as small filled circles. Measurements for the same disc are connected with a narrow gray line. The solid line represents the mean slope for the discs fitted with modified self-similar profile (Table 2), while the dotted line the mean slope for the Gaussian fits (Table 3).
The dashed line represents the radial drift prediction by Rosotti et al. (2019b) discussed in Sect. 6.2. brightest discs in the sample) is: which is steeper than the mean for the whole sample and but still compatible with a flat size-frequency relation.

Constraints on the optical depth
To quantify how much of the disc emission can be attributed to optically thick regions it is useful to introduce a new disc-averaged quantity, the optically thick fraction F , defined as the ratio between the integrated luminosity enclosed within ( being the fraction defined in Sect. 3.2) and the luminosity that would be emitted by a completely optically thick disc with a size : where mm is inferred from the visibility fits and we set = 0.68. Although a larger value (e.g., = 0.95) would presumably yield a closer approximation of the total contribution of the optically thick emission, it would be a problematic choice for the following reasons: first, since the 95% flux-containing radius ( 95 ) is much more sensitive to the slope of the outer disc than 68 is, a constraint on F evaluated at 95 would be inevitably much more uncertain. Second, due to the different signal-to-noise level of different disc observations, we would not be able to achieve a constraint on F at 95 for all the discs of the sample, and for different discs we would have constraints at different values, making the set of measurements inhomogeneous. For the dust temperature d we use the empirical parameterisation by Andrews et al. (2013): where the actual values = 0.5, 0 = 30 K 0 = 10 au were recently calibrated by Andrews et al. (2018b) using ALMA and SMA observations of discs in the Lupus and Taurus region. The stellar luminosities used for the Lupus sources are in Table 1. We ensure that the dust temperature does not reach unrealistically low values below the threshold of floor = 7 K induced by the typical interstellar radiation field in low mas star forming regions by using an effective dust temperature equal to 4 = 4 d + 4 floor . Note that F should not be regarded as the average optical depth of the disc because, by construction, it lies between zero and one. Nevertheless it is a measure of the dominance or otherwise of optically thick emission in the integrated flux.
We highlight that F is not directly measurable from the observations, as it requires knowledge of the size of the disc (which we obtained through the visibility fits) and the assumption of a dust temperature profile. Compared to a simple measurement of the integrated flux ( ), F is intrinsically more model-dependent. However, it is a useful observational quantity that can be determined robustly from spatially resolved observations and with reasonable assumptions on the dust temperature and, advantageously compared to , leverages the information on the spatial distribution of the disc brightness. The lack of spatial resolution that affected submm/mm observations until recent years made it rarely possible to characterise discs through F . Here, aiming to take full advantage of the resolving power of these ALMA observations, we will use F to gain insight into the structure of discs. Figure 5 summarises the distribution of optical depth fractions that we derive at the three wavelengths. The left panel shows, for each wavelength, F as a function of disc effective size. As expected, the contribution of optically thick emission to the total integrated flux decreases significantly at longer wavelengths, with median F values decreasing from 0.71 ± 0.06 at 0.9 mm, to 0.49 ± 0.05 at 1.3 mm, to 0.34 ± 0.06 at 3.1 mm. Moreover, at 3.1 mm we notice that there is a marked correlation for which the largest discs are also those with lowest F . The right panel shows, for each disc, F as a function of integrated flux at the three wavelengths, with both quantities normalised to their values at 0.9 mm. To ease the interpretation of the plot, measurements belonging to the same disc are connected with a thin grey line. The plot shows that there is a drop in optical depth at longer wavelengths.

Millimeter continuum size-luminosity relation
A correlation between the millimeter continuum disc sizes ( 68 ) and their flux at 0.9 mm was found by Tripathi et al. (2017) using SMA observations of a sample of bright discs in the Taurus region, for the discs in the Lupus region, was interpreted as a constant surface brightness (i.e., ∝ 2 ) with an average optically thick fraction of about 0.3. Here we revisit the size-luminosity correlation in the context of the multi-wavelength observations that we obtained at 0.9, 1.3 and 3.1 mm, looking for the presence of the same scaling relation at 1.3 and 3.1 mm. Figure 6 presents the Lupus discs radii ( 68 ) against their millimeter luminosity (re-scaled at the common distance of 150 pc) measured at 0.9, 1.3, and 3.1 mm. To quantitatively characterise the properties of the size-luminosity scaling relation, in the simple assumption of a linear correlation in the logarithmic space (i.e., a power-law correlation in the linear space), we parametrise the relation as log 68 au = A + B log cos 150 pc 2 + where A is the intercept (normalisation), B is the slope (powerlaw index) of the relation, and is a Gaussian scatter term along the y axis ( 68 , in this case). The implications of the cos( ) term are discussed at the end of this Section. We perform a Bayesian linear regression with a mixture of 2 Gaussian generative models following the method by Kelly (2007) and using the implementation in the linmix 3 Python package. Table 5 summarises the properties of the posterior inferred for the A, B parameters, as well as for (the standard deviation of the scatter term) and for the correlation coefficient . The corner plots of the MCMC used for the linear regression show that the posteriors of these parameters are singlepeaked and well behaved.
The results of the linear regressions performed at 0.9, 1.3, and 3.1 mm can be summarised as follows: (i) at 0.9 mm: we confirm the presence of a very tight correlation (correlation coefficient = 0.91 ± 0.04). Although we do not 3 Code available at https://github.com/jmeyers314/linmix. recover exactly the same slope reported by Andrews et al. (2018a), our regression is compatible within 2 with their results. The somewhat steeper slope that we find is likely due to the fact that we miss some of the faintest targets in the sample analysed by Andrews et al. (2018a), for which they typically infer uncertain disc sizes that are compatible with very large values: overall, these faint (and possibly large) discs are likely to have a flattening effect on the size-luminosity correlation. (ii) at 1.3 mm: we discover a correlation that is essentially identical to that at 0.9 mm, except for a larger normalisation. (iii) at 3.1 mm: we discover a new correlation, which is steeper than the ones at shorter wavelengths (slope is 0.83 as opposed to 0.63) and has a significantly larger normalisation. Due to the larger uncertainties of the fainter 3.1 mm emission, the slope inferred at 3.1 mm is still formally compatible within 2 with the slopes at shorter wavelengths.
An immediate interpretation for the increase in the intercept from 0.9 to 3.1 mm is that overall the emission becomes more optically thin at longer wavelengths: for a given disc size, discs emit fainter and optically thinner emission at longer wavelength. This is in line with the results discussed in the previous section and shown in Figure 5. The difference in slope, instead, suggests that there is a systematic effect for which the emission of large and small discs behaves differently across wavelengths, with most of the difference occurring between 1.3 and 3.1 mm. Using models of grain growth and drift, Rosotti et al. (2019a) found that a disc size -millimeter luminosity relation is to be expected if radial drift is the dominant mechanism setting the max- imum grain size. They predict that such relation should have the same slope if observed at different wavelength: although this is what we find for the 0.9-1.3 mm wavelength range, between 1.3 and 3.1 mm we observe a steeper slope. In the same scenario, they predict that the intercept of the size-luminosity relation should scale as 2 ( being the observing wavelength). We find that there is a broad agreement with this result, with the observed intercepts being within 20% of the expected values.
To further investigate the presence of a systematic effect across wavelengths, in the left panel of Figure 7 we plot the disc luminosity at the three wavelengths (normalised to the 0.9 mm value) as a function of the disc size measured at the same wavelength. It is evident that there is a systematic trend in which the luminosity at 3.1 mm consistently decreases by larger amounts as the disc size increases. Another way to look at the same data is shown in the right panel of the same Figure, which reports the 0.9-1.3 mm spectral index as a function of the disc size measured at 3.1 mm. Despite two outliers around 100 au, it is clear that the spectral indices tend to be larger for larger discs.
We shall note that, compared to Andrews et al. (2018a), here we re-scale the integrated flux by cos( ): effectively, we are considering the face-on luminosity of the discs (i.e. an intrinsic disc property) rather than their integrated flux (which inevitably includes the effects of the viewing geometry). It is worth noting that theoretically only the optically thick part of the emission scales with cos( ), while the radiation emerging from the optically thin parts is independent of the viewing angle. For completeness, we also perform the linear regression without the cos( ) correction (namely, as in Andrews et al. 2018a): the results are reported in Appendix A. Interestingly enough, we find that the relations between disc size ( 68 ) and face-on luminosity (presented in this section, results in Table 5) are much tighter (scatter is halved) than the relations between their size and integrated flux (results in Table A1). The tightening of the size-luminosity relation compared to the size-flux one suggests that the discs have a substantial optically thick contribution to their emission, broadly in line with the rather high optically thick fractions that we observe (Sect. 5.2). We defer the discussion of more detailed implications of this effect to dedicated future investigations.

INTERPRETATION OF THE DATA
The most striking result from our observational dataset is the fact that there is very little variation in disc size as measured at wavelengths in the range 0.9 − 3.1mm (see Figure 4). We first consider whether this result can be explained in terms of disc truncation and then discuss whether it is compatible with the expectations of grain growth models and radial drift in smoothly structured discs. Having found that neither of these scenarios are compatible with the data we undertake a more detailed exploration of the combination of parameters that are consistent with the observed properties of our sample.

Compatibility with truncated discs models
An obvious scenario that could produce a truncation in the disc brightness across multiple wavelengths would be a genuine drop in dust surface density, where 68 would reflect a physical disc outer edge. Our sample partially overlaps that of the high resolution DSHARP survey (7 common sources, see note 3 in Table 1) and so we can check whether such a steep decline in emissivity is indeed characterising our sources. We note that our finding that 68 is constant with wavelength is common to discs that are seen to be highly structured in DSHARP images (e.g., IM Lup) and also to those showing modest substructures and a smooth decline in surface brightness in the outer disc (e.g. MY Lup, Sz 114). We thus conclude that, although a surface density truncation could explain the constancy of 68 in some discs (namely, those that exhibit a sharp continuum outer edge), some further explanation is needed to account for the constancy of 68 in the whole sample.

Compatibility with radial drift models in smooth discs
Another way to produce an 68 that is roughly constant with wavelength is the case in which the surface brightness profiles at the three wavelengths are close to being scaled versions of each other, as appears to be the case, for example, in the case of Sz 83 (Figure 2). This is equivalent to saying that there is little radial variation in the spectral index, as is illustrated in the lower panel of Figure 2. This is however surprising in the context of models for grain growth and evolution which predict a rising spectral index with radius, due both to optical depth effects and the evolution of the grain size distribution towards a larger maximum grain size in the inner disc. For example, Rosotti et al. (2019b) found that, in such models, 68 at a given wavelength is closely related to the location in the disc at which the maximum grain size corresponds to an opacity resonance at that wavelength (i.e. where max ∼ /2 ), since the opacity drops steeply for grain populations with smaller values of max . This cliff in disc opacity (of amplitude 8 − 10 for compact grains) should imprint itself on the surface brightness profile regardless of the disc optical depth (unless the disc is so dense that it is optically thick even outside the opacity cliff, a situation that is never encountered in the case of smooth disc profiles with realistic parameters). In the case that the maximum grain size decreases with radius, the location of the opacity cliff is expected to move inwards at longer observing wavelengths, where the grain size corresponding to the opacity feature is larger. Thus the clear implication is that discs are expected to decline in size as the wavelength of observations changes from 0.9 to 3.1 mm.
In Figure 4 we plot the size-wavelength relation predicted in the radial drift regime, computed as the the average of the radial drift dominated models presented in Rosotti et al. (2019b): the expectation is indeed far from the observed properties of the sample, predicting much smaller radii at 1.3 and 3.1 mm than observed. For reference, with respect to 68 measured at 0.9 mm, this would correspond to a 17% smaller radius at 1.3 mm and a 40% smaller radius at 3.1 mm.
Our Lupus dataset indicates that the grain populations probed by 0.9-3.1 mm observations are typically co-located. In Figure 8 we depict the Lupus data (grey lines) extrapolated to longer wavelengths (dashed grey lines) and show that even at wavelengths of 1 cm, the predicted change in size is still relatively modest (less than a factor of 2).
In Figure 8 we also show the results of previous multiwavelength measurements of disc sizes from the literature (Pérez et al. 2012;Tripathi et al. 2018;Tazzari et al. 2016), which employed observations from the Submillimeter Array (SMA), the Combined Array for Research in Millimeter-wave Astronomy (CARMA) and the Karl G. Jansky Very Large Array (VLA), covering a wavelength range from ∼ 0.88 mm to 1 cm. It is striking that in the 0.88-3.1 mm wavelength range, all these literature measurements are consistent within the uncertainties with the distribution of size-frequency relations inferred for the Lupus sample.  Figure 8. Comparison of the size-frequency slopes for Lupus discs (grey lines) and literature measurements (coloured lines). Thick grey lines represent the slopes inferred for for Lupus discs in this study, with extrapolation to longer wavelengths as dashed grey lines.
Lupus slopes. Note that so far only the brightest (and thus largest) discs have been probed at long ( > 7 mm) wavelengths.
The homogeneous analysis that we present here for a sample of Lupus discs sheds new light on this issue by revealing that there is a sizeable fraction of the disc population with a significantly flat size-frequency relation, in which the bright objects probed so far at long wavelengths lie at the steeper end of the distribution.

An exploration of viable disc parameters
Given the difficulty of explaining the near wavelength-independence of 68 in this dataset we now resort to performing a suite of simply parametrised models in order to explore what emissivity properties could in principle explain our results. Given the limited resolution of our dataset, it is not fruitful to undertake detailed modeling of radial profiles of spectral indices of individual sources (cf. Carrasco-González et al. 2019). Instead we want to consider a set of global properties of each disc that can we can reliably extract from our observations and then consider what sorts of models are compatible with the ensemble. In order to assess the requirements that apply to discs with a variety of brightnesses and radial sizes we consider three dimensionless numbers: (i) the ratio between 68 at 3.1 mm and 68 at 0.9 mm; (ii) 0.9−3.1 mm , the disc-averaged spectral index; (iii) F 1.3 mm , the optically thick fraction at 1.3 mm.
As in many previous studies (Pérez et al. 2012;Tazzari et al. 2016;Tripathi et al. 2018;Huang et al. 2020) we can regard the spectral index as containing information both about the emissivity properties of the dust and/or the prevalence of optically thick emission, while the optically thick fraction primarily provides information on the latter issue. Since neither of these quantities can be simply mapped onto a unique disc emissivity profile, we conduct some simple forward modelling in order to constrain the types of parameters that match the values of the above mentioned three quantities in the Lupus dataset. Figure 9 presents the data in the plane of optically thick fraction at 1.3 mm versus spectral index between 0.9 and 3.1 mm and shows that it tends to be clustered at moderate values of the optically thick fraction and low values of spectral index, with the bulk of the population residing in the range: 0.2 ≤ F 1.3 mm ≤ 0.6 and 2.4 ≤ 0.9−3.1 mm ≤ 3.0 .

(14)
We now run a series of simple models with different optical depth radial profiles and search for those models that best reproduce these observational constraints. Note that in the following two sections our simple models neglect the role of scattering of mm emission. Inclusion of scattering would not change the predictions of the models that we find to be viable in the next two Sections since the modest observed optically thick fractions for most sources does not require the dominance of optically thick emission (where scattering modifies the spectral index). Nevertheless, an alternative explanation of the available data is one where the emission is dominated by optically thick emission with low areal filling factor where the scattering albedo is sufficiently high (requiring grains larger than 1 mm). We discuss the potential effect of scattering in more detail in Appendix D.

Smooth radial profiles
We first run models with the optical depth described by a modified self similar spatial profile, and a power-law spectral dependency: where ( ) increases linearly from 0 in the inner disc to an asymptotic value out at a radius , namely: This implementation gives us the flexibility to realise a simple scenario in which (and therefore the maximum grain size) is constant throughout the disc ( = 0, = out ), as well as the scenario in which grains are larger in the disc interior (within ). For the dust temperature we use the parameterisation in (11). The brightness of the disc at a given observing frequency is computed as Optically thick fraction at 1.3 mm as a function of 0.9-1.3 mm spectral index for a grid of models with constant throughout the disc ( = 0). Dashed gray lines connect models with same , from = 0 (left-most) to = 4 (right-most) in steps of 0.5. The coloured map shows the ratio of disc size at 3.1 and 0.9 mm, with labelled contours ranging from 1.1 to 0.8.
In these toy models we assume a face on disc ( = 0), a Solar luminosity star ( ★ = 1 ), 1 = −1, 2 = 1, and = 30 au (note that the model trajectories in the F 1.3 ,mm − 0.9−3.1 mm plane are only very weakly dependent on since the models are scalefree apart from the weak dependence on introduced by the temperature parametrisation. For a given value of , we compute a grid of models where we vary the optical depth normalisation 0 = 10 −3 , 10 −2 , 10 −1 , 1, 10, 100, the outer dust spectral index value out from 0 to 4 in steps of 0.5. Note that we set 0 = 345 GHz, corresponding to the observing wavelength of 0.869 mm. Figure 9 presents the model results for the case of = 0, i.e. is spatially constant. The dashed lines in Figure 9 each represent a sequence of models with increasing optical depth normalisation at fixed out . Each dashed line intersects the −axis at the spectral index for optically thin emission with = out . The fact that this intersection occurs at a value of 0.9−3.1 mm that is somewhat less than 2 + is a consequence of the fact that the outer parts of the disc are not entirely in the Rayleigh-Jeans regime. Figure 9 also depicts contours of the ratio of 68 at 3.1 to 0.88 mm, which is also denoted by the colour scale and where the pale shadings correspond to a ratio near unity. Detailed properties for some models representative of different regimes in the F 1.3 mm − 0.9−3.1 mm plane are presented in Appendix B. Figure 9 immediately demonstrates that the models that pass through the region of the F 1.3 mm − 0.9−3.1 mm plane occupied by the data also automatically satisfy the requirement of having similar 68 values in the two wavebands. Most of these successful models have input out values in the range 0.5 − 1, with a few sources being compatible with up to 1.5. We have experimented with a variety of monotonically declining surface density profiles and this conclusion remains robust (see below). Such a range of acceptable values is unsurprising given that the observed 0.9−3.1 mm values are rather low and yet the optically thick fraction is insufficiently high for these low 0.9−3.1 mm values to be explicable purely in terms of high optical depth. It is also unsurprising that a constant model with only moderate optical depth should yield a wavelength independent disc radius since the radial profile of spectral index is in this case rather flat and hence the emissivity profiles in the two wavebands are very similar.

Radially increasing
We now relax the assumption of radially constant input . Each panel in Figure 10 has a fixed value of and the various dashed lines in each panel correspond to different values of out .
It is immediately obvious that large values of (as in the lower two panels) fail to reproduce the data. Firstly, they have too much emission from low material at small radius to be able to replicate the larger values in the sample. However, what is more restrictive is the set of constraints imposed by the requirement that the radii at 3.1mm and 0.88mm should be nearly equal. As can be seen particularly in the lower-left hand panel, there are models that pass through the data in the F 1.3 mm − 0.9−3.1 mm plane but where the radius at 3.1 mm is substantially smaller (factor two) than that at 0.88 mm. These models fail because they predict too much variation in over a region of the disc that is still contributing significantly to the total flux. This translates into a radially variable spectral index and the resulting difference in brightness profiles in the two bands then leads to different corresponding 68 values. To test the robustness of this result we have computed the grid of models for a range of 1 , 2 values, reproducing discs with significantly flatter interior ( 1 = −0.3) and (or) steeper outer edge ( 2 = 3): although the actual value of that best fits the data changes slightly from case to case, the results that we have just presented do not change.
It is notable that this modeling rules out a scenario which has been proposed to explain spectral index data, where the inner disc is optically thick while the outer disc consists of a region of optically thin emission with small grains (for which = 2): see, e.g., Ricci et al. (2012). This combination can reproduce intermediate values of the spectral index as well as the observed values of the optical depth fraction, as seen, for example, in the lower left panel of Figure  10 where contours with out = 2 pass through the region occupied by the observational data. However, it can be seen that this region of the plot is shaded dark red, indicating that the predicted R68 values at 0.9 and 3.1 mm are very different. This can be readily understood in that the flux profiles at each wavelength are shaped by the radius at which the emission makes the transition from being optically thick to thin. In the case of high , the significantly lower opacity at 3.1mm drives this transition to smaller radius and hence results in a steeply declining disc size as a function of observing wavelength.
We therefore conclude from this simple exercise that the way to reproduce the typical spectral indices, optically thick fractions and multi-wavelength 68 ratios of the Lupus discs is to invoke a dust distribution where the value of the opacity index is in the range of 0.5 − 1 for most sources, at least over a substantial fraction of the disc, although some lower material at smaller radius is also allowed. We will discuss the implications of this result for the properties of grains in discs in Section 6.5 below but now turn to the question of how these conclusions would be modified in the case of discs with significant substructure.

Structured radial profiles: the case for small grains
We have shown that a class of models that satisfies the observational constraints listed above is one in which is in the range 0.5 − 1 over much of the disc, translating into a spectral index profile that plateaus with values in the range ∼ 2.5 − 3. If this is interpreted in terms of grain properties we will see in Section 6.5 that this corresponds to large grains ( max > 1mm).
However, before fixing upon this conclusion, we will now examine the alternative possibility that there is no significant grain growth to scales > 100 m and that the required spectral index is obtained via mixing the emission from such small dust grains with optically thick substructures. We have shown above that this does not work if the material with high spectral index is placed at large radius because this predicts a steeply decreasing disc size as a function of wavelength. However in this section we consider the case where optically thick and thin material are co-located at all radii. From the point of view of our analysis, it does not matter whether these substructures are large scale coherent structures that would be potentially resolvable with long baseline observations or whether they represent small scale dust condensations that would always remain at the sub-beam level.
Our purpose therefore is to try and construct a model that circumvents the conclusion that grain growth to mm scales in protoplanetary discs is inescapable. In this model we posit a background of small grains with = 2 and a population of optically thick substructures. The details of the model are given in Appendix C. We prescribe the radial variation of the area filling factor of optically thick substructures and of the optical depth of the medium and generate large suites of models that alter the balance between emission in the two components and the optical depth normalisation for the background. Given the degree of flexibility in the models, it is unsurprising that we find some parameters that work; successful models however occupy a narrow niche of parameter space. In practice we find that the only models that work are those where (i) there is a roughly equal balance between the emission from optically thick substructures and background emission and (ii) this balance shows little variation across the disc. The rough balance between the two emission components is set by the requirement that the integrated spectral index has an intermediate value between the values of 2 and 4 expected for optically thick emission and for small grains in the optically thin limit, respectively. The lack of radial variation in the balance between these components is set by the need to reproduce the similar 68 values at 0.9 and 3 mm. Figure 11 depicts a model in the small niche of parameter space that is able to reproduce typical properties of discs in Lupus. The top panel depicts a quarter of the disc, where the brightness of the background of small grains is represented with the colour scale, and the optically thick substructures as black circles: the area covering factor of optically thick structures is ∼ 0.04 at all radii. We emphasise that it makes no difference to the optical properties whether these optically thick structures are positioned randomly (as shown) or arranged in rings or spirals. The three lower panels represent the background optical depth, the total brightness (including emission from the optically thick structures), and the emerging spectral index profile, respectively.
We conclude from this that if there is no grain growth to max > 100 m, and if the observed spectral indices are explained in terms of optically thick substructures, this can be made to work only under contrived conditions regarding the balance of the emission components. Whereas we of course could not rule out such an interpretation for a particular source, the fact that all the sources would require this particular combination of parameters makes us disfavour this possibility.
We emphasise, however, that our analysis is not disfavouring the possibility of optically thick substructures in general but that this cannot readily be made to work if the distributed dust component is composed of very small grains (with = 2). We found that when we modeled a mixture of substructures and a background composed of large grains ( = 1), a wide variety of realisations were broadly compatible with the observed system properties.

Summary of constraints derived from the data
The data is rather tightly clustered around spectral index ∼ 2.4 − 3, with a moderate optically thick fraction (most sources in the range 0.2 − 0.6) and a ratio of radii at 3.1 mm to 0.88 mm that is close to unity. This combination of parameters can be easily realised by smooth disc profiles where the opacity index, , is spatially constant and with a value in the range 0.5 − 1 in most cases. For constant models, the main constraint is provided by optically thick fraction and spectral index data. Models that match this data automatically satisfy the requirement of having very similar radii in the two wavebands.
It is also possible to find models where is a smoothly varying function of radius. For example there is a wide range of profiles with low in the interior and large at large radius which can accommodate the observed values of spectral index and optical depth function. Crucially, however, most such solutions predict that the disc is significantly smaller at 3.1 mm than at 0.88 mm (see the red regions in Fig. 10). This is because if is increasing to high values over regions of the disc producing a significant fraction of the flux, the longer wavelength emission is being down-weighted at large radius, resulting in small 68 values. We find that the only profiles that are also consistent with the wavelength independence of 68 are those where the asymptotic value of (namely, out ) is in the range 0.5 − 1 and where only the innermost regions (10 − 15% of 68 ) can have lower values. The main constraint on the size of interior regions with low is set by spectral index data, since a large inner region of low would lead to predicted spectral indices that are too low. We have also explored whether the above constraints on the permissible profiles of necessarily constrain the microphysical dust properties (see Sect. 6.5 below).
Furthermore, we have investigated whether a mixture of optically thick regions and optically thin regions with much higher values (∼ 2) could be consistent with the data. We conclude that this can match the data only if the ratio of optically thin to optically thick emission does not vary with radius and is of order unity. In the absence of a reason for expecting such a universal distribution for the mixture of optically thin and thick material, we disfavour this possibility.
Finally, a further possibility is that the disc emission is dominated by optically thick substructures but that these consist of grains where the albedo is high and increasing towards longer wavelengths. Such substructures can produce emission with spectral index similar to that observed without the need to add an optically thin background and automatically predict that the disc radius is insensitive to wavelength.
It is worth emphasising two further points. First of all our conclusion that the data favours ∼ 0.5 − 1 or else optically thick emission from grains with high albedo that is an increasing function of wavelength refers to the bulk of the emission: we cannot rule out the presence of localised (and spatially unresolved) optically thin zones with higher (thus, higher 0.9−3.1 mm ) provided that they are not a major component of the flux. Secondly, while we conclude that the data cannot be readily accommodated by a mixture of high ∼ 2 material and optically thick zones, the data is readily fit by range of models combining optically thick emission with distributed low ∼ 1 material or by purely optically thick emission with suitable scattering properties.

Implications for dust properties
We have argued above that in the case of smoothly structured discs, = 0.5 − 1 over the bulk of the disc in most sources and that we cannot readily reconcile the data with a mixture of emission from regions with high combined with optically thick substructures. In this case, the value of can be immediately linked to the emission properties of the dust and hence the grain size distribution. Figure 12 depicts the theoretical value of the dust opacity spectral index between 0.9 and 3.1 mm as a function of maximum grain size max for a variety of assumptions about the grain size distribution and grain porosity.
The opacity curves were obtained using the dsharp_opac Python package 4 (Birnstiel et al. 2018), which implements a Mietheory dust opacity model, with appropriate mixing rules for composite materials. We compute the opacity for compact grains (with = 30 au, 1 = −1, 2 = 1) for the background and optically thick structures depicted as black dots. A dashed white line highlights the scale radius . The three panels at the bottom show the optical depth of the background ( ), the total brightness (including optically thick structures) normalised at the inner radius ( / 0 ), and the emerging spectral index profile. Radii enclosing 50, 68, and 95% of the emission are shown as thick circles in the brightness plot. a composition that is labelled default in the DSHARP analysis, see Sect. 2 in Birnstiel et al. 2018) and for porous grains (with same fractional abundances of compact grains, but with 80% porosity), and for different values of = 2.5, 3, 3.5, with being the slope of the power-law grain size distribution ( ) ∝ − for min ≤ ≤ max ( being the dust particle radius). The curves shown in Figure 12 refer to the absorption opacity, as is relevant for interpreting spectral indices of optically thin emission. The range of values that are consistent with the bulk of our Lupus data are denoted by the blue band. Figure 12 demonstrate that only compact grains with a relatively top-heavy grain size distribution (grain size power index = 2.5 or 3) are compatible with the bulk of data and then, only if the maximum grain size is in excess of 1 mm. Such a maxi- Lupus discs compact, q = 2.5 compact, q = 3 compact, q = 3.5 porous, q = 2.5 porous, q = 3 porous, q = 3.5 mum grain size is above the grain size corresponding to the opacity resonance feature at all the wavelengths studied. The fact that the predicted opacity curves are rather flat for larger maximum grain sizes shows that the observed value is compatible with a wide range of maximum grain sizes, consistent with the apparent ubiquity of this value for the Lupus sample. While top heavy ( = 2.5 − 3) distributions of compact grains with max larger than a few mm have a predicted spectral index which is nearly independent of max , the corresponding opacity value is steeply dependent on max in this range. Thus changing max from the mm to the cm range reduces the opacity by around an order of magnitude. Thus although our measurements have placed a lower limit on max , accurate estimation of disc mass requires further observations that can distinguish between max values in the mm and cm range. In Sect. 6.6 we discuss how future observations at longer wavelength will be able to distinguish these scenarios.
Models for dust growth in discs Birnstiel et al. (e.g., 2012) predict that grains rapidly grow to scales substantially more than a mm but that they are then subject to radial drift which drives the dust to gas ratios down to very low values (e.g. 10 −4 by an age of 3 Myr). Given typical dust masses estimates in protoplanetary discs, such low dust to gas ratios would imply unacceptably large gas masses so it is widely believed that radial drift is somehow inhibited, a popular option being the existence of dust traps in local pressure maxima (Pinilla et al. 2012), possibly associated with the presence of protoplanets. Our demonstration that emission from discs in Lupus is dominated by grains larger than 1mm (which would otherwise be subject to rapid radial drift) is strong, albeit indirect, evidence that dust trapping must be effective in these discs. The resolution of these observations does not allow the detection of sub-structure but some of the objects have been targeted by the DSHARP survey (Andrews et al. 2018b) which finds evidence of annular structures of various strengths in most of the sources targeted.

Predictions at longer wavelengths
The Q and Ka-Bands of the VLA offer the possibility to observe protoplanetary discs in the wavelength range between 7 mm and 1 cm (e.g., with spectral windows centered at 42 and 30 GHz), which will also be accessible in future via ALMA Band 1 (planned frequency coverage between 35 and 50 GHz). In future it will be possible to assemble information on disc radii, spectral indices and optical depth fractions for large samples of discs at wavelengths around three times longer than in the present study; to date, however, this information is available for relatively few discs (see Figure 8). By comparing such new measurements to the disc radii, masses, and spectral indices of discs in different regions or different evolutionary stages (e.g., the Class 0/I objects in the Orion and Perseus clouds; Segura-Cox et al. 2018;Tobin et al. 2020;Sheehan et al. 2020) it will be possible to trace the evolution of solids across different stages of disc evolution.
While future datasets should be subject to a modeling exercise similar to that conducted here, it is already possible to make some broad statements about how such measurements could be used to distinguish between the various possibilities that are compatible with the present dataset. At ALMA wavelengths, we have found two types of models that can readily fit the data i.e. smooth models where the emission is predominantly optically thin (see Section 6.3.1) and the grain size is in excess of ∼ 1mm and models where the bulk of the emission derives from optically thick material with a moderate area filling factor see Section 6.4; in the latter case, it is also necessary that the optically thick emission results from grains which are similarly large since such grains have a high scattering albedo which increases towards longer wavelengths and can explain the observed spectral indices.
In either scenario, if the emission is dominated by very large grains (i.e. pebbles on a cm scale or larger), both the absorption opacity and scattering opacity can be described as a single power law over wavelengths ranging from 1 mm to 1 cm (Carrasco-González et al. 2019). In this case we would expect little change in the spectral index in this wavelength range and that the 68th percentile flux radius would vary little as a function of wavelength (as in the grey lines in Figure 8). The smooth variation of optical properties of large grains over this wavelength range would also imply that the optical depth fraction would be slightly lower at 1 cm than at mm wavelengths, continuing the trend seen in Figure 5.
If instead the grain size were towards the lower end of the range allowed by our present modeling (i.e., around a few mm, so in excess of the size corresponding to the opacity resonance at 1 − 3 mm but below the resonant value at cm wavelength), the predictions for the structured models at a wavelength of 1 cm would be somewhat modified, since the spectral index declines towards longer wavelengths for optically thick emission when scattering is included (Zhu et al. 2019). However in the case of smoothly structured, largely optically thin models with maximum grain size of a few mm, the spectral index would be expected to rise significantly between mm and cm wavelengths, reflecting the abrupt reduction in opacity in the case that the grains are significantly smaller than the wavelength of emission. This rise in spectral index would produce a more marked decrease in the optical depth fraction and in the disc size as a function of wavelength.
We therefore conclude that observations at longer wavelengths have the capacity to further constrain the grain size distribution (discriminating between mm and cm scale grains) though detailed modeling would be required to firm up these expectations.

CONCLUSIONS
In this paper we have presented the analysis of multi-wavelength ALMA observations at 0.88, 1.3, and 3.1 mm of 26 protoplanetary discs in the Lupus star forming region. The observations have an angular resolution between 0.25 and 0.35 and a comparable sensitivity.
We have derived the multi-wavelength radial brightness profiles of these discs by fitting the interferometric visibilities with parametrised brightness models. At each wavelength, we derived the disc effective size ( 68 , enclosing 68% of the disc emission), and the spectral index radial profiles. The homogeneity of the observations (in terms of sensitivity and resolution) and of the analysis enabled us to characterise and compare the properties of all discs with a minimum relative bias across the wavelengths.
We emphasise that the fact that the discs are spatially resolved at multiple wavelengths (hence we have information on their size) allowed us to break degeneracies in interpreting the spectral index information alone. By forward modelling these observations with simple toy models for the disc emission enables us to present the strongest evidence to date that substantial grain growth (to scales > 1 mm) is required in a large sample of discs, irrespective of their fluxes and radii.
The main results can be summarised as follows: (1) millimeter continuum size-luminosity relation: we confirm the relation at 0.9 mm (Tripathi et al. 2017;Andrews et al. 2018a) and we discover that such relation is present with the same slope when observed at 1.3 mm, and with a steeper slope when observed at 3.1 mm, suggesting that large discs are preferentially characterised by a larger 0.9−3.1 mm spectral index.
(2) millimeter continuum size-frequency relation: in the 0.88-3.1 mm wavelength range this relation is flat (i.e., 68 is wavelength-independent), indicating that grains emitting at 0.9-1.3mm ( max ∼ 0.17 mm) are essentially co-located. The sizefrequency relations found for the Lupus discs are compatible with literature measurement, which typically map the mmbrightest objects and lie at the steeper end of the distribution of slopes. This is confirmed by the tentative evidence that the steepness of the millimeter continuum size-frequency relation correlates with the 3 mm disc luminosity. (3) except for the peculiar case of IM Lup, our analysis indicates that most of the Lupus discs require large grains ( max > 1 mm) at large radii (50-100 au), implying that radial drift has to be significantly halted. (4) using the optically thick fraction F to estimate the amount of emission that can be ascribed to optically thick regions, we prove that Lupus discs are systematically optically thinner at longer wavelengths. Lupus discs are clustered around optically thick fractions 0.2 ≤ F 1.3mm ≤ 0.6, spectral indices 2.4 ≤ 0.9−3.1 mm ≤ 3.0, and multi-wavelength size ratios 68,3.1mm / 68,0.9mm (0.92 ± 0.10). (5) by modelling observations with simple models of smooth discs we conclude that a ready way to reproduce the Lupus measurements of F 1.3mm , 0.9−3.1 mm , and size ratios is to invoke a dust distribution where the opacity index is in the range of 0.5 − 1 over a substantial fraction of the disc, although some lower ∼ 0 material at smaller radius is also allowed. (6) we also model observations in an alternative scenario in which discs are populated by small grains ( max < 100 m) and a large number of optically thick substructures: we find that this model can work only under contrived conditions regarding the balance between the emission from the optically thin and thick regions.
Although we could not rule out such an interpretation for a particular source, the fact that all the sources would require this particular combination of parameters makes us disfavour this possibility. (7) a further possible scenario is that the emission is instead entirely dominated by optically thick substructures with a small area filling factor. In this case it is necessary, in order to explain the fact that the spectral indices are significantly larger than 2, that these substructures are composed of large grains with a high scattering albedo. The constraints on required grain size in this scenario are similar to those for the smoothly structured models described above. (8) in terms of grain growth, the observations of the bulk of the Lupus sample can be explained with (0.75 ± 0.25), which can be produced only by compact grains with a top-heavy grain size distribution ( = 2.5 − 3) and max > 1 mm. Weidenschilling S. J., 1977, MNRAS, 180, 57 Wrobel J. M., Walker R. C., 1999, in Taylor

APPENDIX A: MILLIMETER CONTINUUM SIZE-INTEGRATED FLUX RELATION:
Here we present the results of the linear regression between the disc size ( 68 ) and the integrated flux ( )), as opposed to Sect. 5.3 where we tested the correlation against the disc face-on luminosity.
The linear regression is now parametrised as: namely, without the cos( ) re-scaling term in (13). Except from the slightly different re-scaling distance (150 versus 140 pc), this is the same parametrisation used in Andrews et al. (2018a). By using the same Bayesian linear regression described in Sect. 5.3 we obtain the results reported in Table A1. The 0.9 and 1.3 mm relations are very similar to those found for the face-on luminosity (Sect. 5.3), while the 3 mm slope is significantly flatter in this case. Most notably, in these latter linear regressions we find a significantly larger scatter. Figure A1 displays the three correlations with the same colour and line conventions used in Figure 6.

APPENDIX B: DETAILED TOY MODEL PROPERTIES: SMOOTH
Here we present the detailed properties of the toy model with smooth structure used in Sect. 6.3.1. To document the behaviour of the model, in Figure B1 we show three models (a, b, and c), that are representative of different regimes. The three models have different input ,0 and , which allows us to reproduce different disc-integrated values of F 1.3 mm and 0.9−3.1 mm . All the models have radially uniform and assume = 30 au, 1 = −1, 2 = 1. For each model we present detailed properties: brightness profiles, optical depth and effective dust temperature, observed spectral index and input dust opacity spectral index. The 50%, 68%, 95% flux enclosing radii are highlighted as filled circles in the brightness profiles.
Model (a) has ,0 = 0.07 and = 0.7: the low opacity spectral index produces slowly varying 0.9−3.1 mm ( ) profile, resulting in a disc radius that is essentially constant across 0.9-3.1 mm wavelength range. Model (b) has higher = 1.3 and larger variations in 0.9−3.1 mm ( ): its 3.1 mm radius is ∼90% the 0.9 mm radius (as typically observed for the Lupus discs). Model (c) has a very large = 2, which makes 0.9−3.1 mm ( ) reach large values ∼ 3.5: this produces a strong reduction of the disc size from 0.9 to 3.1 mm.
The location of these three models in the same F 1.3 mm − 0.9−3.1 mm plane as in Figure 9 are highlighted in Figure B2.

APPENDIX C: DETAILED TOY MODEL PROPERTIES: STRUCTURED
In structured models we posit a background of small grains (< 100 m) with optical depth profile given by (15) with = 2. We then introduce the quantity ( ) = ( ) (where is a scaling factor that can in principle be a function of radius) such that the fraction of the disc area at a given location that is occupied by optically thick substructure is given by 1 − exp[− ( )]. The total emission at that location is then the sum of that from the optically thick substructures and the background of small grains (which occupy an area filling factor of exp[− ( )]). In regions of the disc where the background of the disc is optically thin, the composite emission has radially constant 0.9−3.1 mm if is independent of radius; the value of 0.9−3.1 mm is controlled by the scaling factor which is the model parameter that controls the relative dominance of emission from the substructures and from the background of small grains.
With this prescription we can generate large suites of models that alter the balance between emission in the two components and the optical depth normalisation for the background; in order to find a model that predicts similar radii at 0.9 and 3.1 mm we consider the case of radially constant since these models have little radial variation of 0.9−3.1 mm . Given the degree of flexibility in the models, it is unsurprising that we find some parameters that work; successful models however occupy a narrow niche of parameter space. Not only is it necessary to invoke radially constant (as noted above) but the value of has to be of order unity, since otherwise the emission would be dominated by one of the components and the resulting spectral index would be either too high (close to 4 if is too low) or too low (close to 2 if is too high).

APPENDIX D: THE ROLE OF SCATTERING
The toy-model analysis we presented in Section 6 relied on considering thermal emission without scattering and a prescribed frequency dependence for the optical depth associated with absorption (Eq. 15). Inclusion of scattering opacity has no effect on the emission properties at low optical depth. However, in optically thick regions, scattering reduces the emission below its black body value: the increased path length of photons undergoing multiple scattering means that the effective surface of last emission moves higher in the disc atmosphere, implying that radiation derives from a smaller column of emitting material (Carrasco-González et al. 2019;Zhu et al. 2019). The effect on the spectral index from optically thick regions depends on the wavelength dependence of the albedo which can have either sign. Consequently, the canonical value of = 2 for optically thick emission in the Rayleigh Jeans limit can be modified by scattering to lie in the range ∼ 1.6 − 2.5 (Zhu et al. 2019).
These considerations do not significantly modify the conclusions from the modelling Sections 6.3.1,6.3.2. In the case of the smooth radial profiles considered in Section 6.3.1, the models that match the rather low optically thick fraction values of much of the data contain only a minor contribution from optically thick emission and would therefore not be affected by inclusion of scattering. The models considered in Section 6.3.2 contained only small grains (< 100 m) for which the albedo is very low and thus would not be affected by scattering even within the optically thick substructures.
Consideration of scattering however permits another interpretation of the distribution of the observational data in the plane of spectral index versus optically thick fraction. If the albedo rises with wavelength in the range 0.88 − 3.1mm, then this suppresses the flux at 3 mm in regions of high optical depth, i.e. it increases the spectral index in optically thick regions to > 2. For sufficiently large grains, 0.9−3.1 mm can attain a value of 2.3 − 2.5 (Zhu et al. 2019), which is typical of the observed values in our sample. In this case, the data would be consistent with emission deriving almost exclusively from optically thick structures composed of large grains, with the observed optically thick fraction requiring an area filling factor of optically thick emission of tens of percent. Such a scenario would automatically satisfy the requirement of having the same radii at different mm wavebands, since emission in each of these bands would derive from the same structures. Carrasco-González et al.
(2019) successfully model the resolved multi-wavelength profiles of HL Tau with such a scenario (where in this case the optical depth at 1.3 mm is at least modestly larger than unity area over a large radial range).

APPENDIX E: DETAILED FIT RESULTS
In this Appendix we report the detailed fit results for the 26 discs that we considered for the multi-wavelength visibility modelling (Sect. 3). In Sect. E1 we present the fits performed with the self-similar brightness profile, and in Sect. E2 the fits performed with the Gaussian profile.
(This appendix is available online as supplementary material)

E1 Fits with modified self-similar profile E2 Fits with Gaussian profile
This paper has been typeset from a T E X/L A T E X file prepared by the author.