© G. Guidi et al. 2022

1 Introduction

The evolution of protoplanetary disks during their lifetime of a few million years (Hernández et al. 2007) is intrinsically connected with the birth of planets from their solid and gaseous material. The multiplicity of processes at play make a global description of disks extremely challenging, although the increasing number of disk observations taken in the recent years – from radio to optical and infrared frequencies – are providing new and precious information that is helping to unravel the puzzle.

The emerging picture is that all the disks observed at high resolution present, even at early stages (e.g., 1 Myr for HL Tau, ALMA Partnership 2015), a variety of small-scale structures (cavities, rings, spirals), indicating that local processes play a major role in the disk evolution (e.g. Andrews et al. 2018a; Long et al. 2018). This opposes the previous view of disks with "smooth" and monotonically decreasing surface density profiles, which rapidly disperse under the action of global mechanisms such as viscous evolution (Lynden-Bell & Pringle 1974), consistent with low to medium resolution (≥0.″75) observations in the submillimeter/millimeter available until ~ 10 years ago. Substantial effort in the theoretical interpretation of this plenitude of new observations is being carried out, but the origin of the observed rings and spiral structures is still under debate. The most popular scenarios invoke planets perturbing the disk dust and gas distribution (e.g., Zhang et al. 2018), photoevaporation (Ercolano et al. 2017), and gravitational instabilities (Pérez et al. 2016), but also several magnetized, non-self-gravitating disk instabilities, such as vortices via Rossby waves (Huang et al. 2018c), vertical shear instabilities (Pfeil & Klahr 2020; Blanco et al. 2021), Magnetohydrodynamics (MHD) winds (Riols et al. 2020), warps and their induced instabilities (Deng et al. 2020). Chemical effects, such as rapid dust growth around condensation fronts of the most abundant molecular species (Ros & Johansen 2013) or a pileup resulting from the sintering effect (Okuzumi et al. 2016), have been proposed as responsible for the observed ringed structures as well. A correlation between the position of dust substructures and the location of the snowlines of major volatiles has been observed in some protoplanetary disks, supporting this scenario (e.g., Zhang et al. 2015; Carrasco-González et al. 2019; Sierra et al. 2021).

At the current stage, the question of the origin of the substructures observed in protoplanetary disks is still open. With the notable exception of PDS70, where two planets and corresponding circumplanetary disks were imaged within the disk cavity (Müller et al. 2018; Isella et al. 2019; Benisty et al. 2021), no conclusive explanation can be linked to the observed morphologies, even within the single sources. The comparison with simulations is complicated by the complex radiative processes at play, which make the inference of the gas and dust density structures from the emitted radiation nontrivial. Dust in particular, despite being the most accessible observational tracer in disks, is not so straightforward to account for in global hydrodynamic simulations as it is subject to a variety of processes such as radial drift, coagulation, and fragmentation (e.g., Birnstiel et al. 2010). The presence of rings and gaps in disks implies that dust could accumulate in certain regions which could promote further dust growth. It has been shown that grain growth has strong influences on understanding the lifetime and appearance of rings and vortices in 2D disks (Drazkowska et al. 2019; Laune et al. 2020; Li et al. 2020), and it can affect the optical depth and spectral index of multiple ringed structures (Li et al. 2019).

In this paper we focus on the bright disk around the Herbig Ae star HD 163296, and we use multiwavelength observations in the millimeter and centimeter range to reveal the radial variation in the dust properties with a nominal spatial resolution of about 5–8 au. At a distance of 101.2 ± 1.2 pc (Bailer-Jones et al. 2018), this disk is an ideal candidate for studying planet formation in progress: a ringed structure in the dust emission of this disk was revealed by ALMA (Atacama Large Millimeter/submillimeter Array) (Isella et al. 2016; Zhang et al. 2016), and the presence of multiple planets perturbing the dust and gas dynamics was proposed to explain the observations (Isella et al. 2016). Several additional studies estimated 0.3–1 M_Jup planets at separations of about 50, 80, and 140 au (Liu et al. 2018; Teague et al. 2018) and a ~2 M_Jup planet at 260 au (Pinte et al. 2018). Direct imaging of the giant planets was attempted at infrared wavelengths (Guidi et al. 2018; Mesa et al. 2019), but without any robust detection.

In this work, we take advantage of new high resolution observations from ALMA and VLA (Very Large Array), in addition to archival DSHARP data at 1.3 mm (Andrews et al. 2018a), covering a wide spectral range (from 0.8 mm to 3 cm) to reconstruct the dust distribution in the midplane of the HD 163296 disk. In Sect. 2, we describe in detail the observations and the calibration procedures, we show the continuum intensity maps and corresponding brightness temperature profiles. Section 3 describes the methods we use for the analysis: the extraction of the non-dust component and the radial profiles, and the models setups for fitting the spectral energy distribution. In Sect. 4, we present our results: the spectral index maps and the dust properties we derive from the different models we apply to our multiwavelength analysis. In Sect. 5, we discuss the results in terms of grain size distribution and dust mass, we compare them to previous studies and discuss the implications for the origin of the ringed structure in HD 163296. Finally, in Sect. 6 we summarize our main findings and draw our conclusions.

2 Observations and data reduction

2.1 ALMA Observations

HD 163296 (also known as MWC 275) was observed with ALMA at Band 3 on September 7 and 13, 2016 with an array configuration of 42 antennas covering baselines from 15 m to 2.5 km, resulting in an angular resolution of ~0.3″ and a maximum recoverable scale of 40″. The correlator was set up with four spectral windows in dual polarization mode: three SPWs for the continuum detection were set to cover the total bandwidth of 1.875 GHz, and were centered at 91.142, 103.006 and 104.694 GHz. One SPW was set to have a higher spectral resolution of 61 kHz (~0.2 km s⁻¹ after Hanning smoothing) and was centered at 93.165 GHz to include the N₂H⁺(1–0) emission line. J1924-2914, J1733-1304 and J1751-1950 were observed to calibrate for bandpass, flux and phase, respectively. J1753-1934 was additionally observed every 20 minutes as a check source for the complex gains calibration. Data were calibrated using the ALMA pipeline in the CASA software package (version 5.4.1), then self-calibration was performed on the lower side band (LSB) and upper side band (USB) separately, as the wavelength separation between the two resulted in an appreciable difference in flux (average flux difference of 35% across the uv-space). Phase only self-calibration was carried out for each dataset, amplitude calibration was avoided in order to preserve the flux for multiband analysis. The final S/N (signal-to-noise) was improved by 60%.

Additional observations in Band 3 at a higher angular resolution were taken by ALMA on 16, 19 and 21 September 2017 (PI Isella). The array configuration consisted of 45, 44 and 41 antennas respectively, covering baselines from 41 m to 12 km. The same calibrators listed for the more compact configuration were used for these observations. These were combined with the shorter baselines observations, again separating LSB from USB. Because of the different measurements of the flux calibrators, a flux-scaling was performed before combining the datasets: the reference flux was chosen by looking at the measurements of the correspondent ALMA flux calibrators, selecting the most recent in relation to the date of the observations.

Observations at Band 4 were taken by ALMA in two complementary configurations: on November 13 and 23, 2017, HD 613296 was observed in the more extended configuration (C43–8), covering baselines from about 100 m to 12 km, using 43 and 50 antennas, respectively. The corresponding angular resolution was of 0.″06, with maximum recoverable scale of 1.″7. The total on-source time for HD 163296 was 90 min; J1924-2914 was used for flux and bandpass calibration, J1753-1843 for the phase calibration. The correlator was set with three spectral windows for the continuum centered at 133.987, 135.925 and 145.976 GHz, with a channel width of 976 kHz (~2 km s⁻¹). One spectral window was centered around the CS (3–2) line at 146.963 GHz, with a spectral resolution of 30.5 kHz (~60 m s⁻¹). In order to recover the large-scale emission that is filtered out by the extended configuration, HD 163296 was also observed in the more compact configuration C43–4 on January 17, 2018. Covering baselines from 15 m to 1.4 km and making use of 44 antennas, the target was observed for 15 min with final angular resolution of ~0.″5. Phase self-calibration was performed first on the short baseline dataset (configuration C43–4), then the dataset were combined, flux scaled as for the 3 mm datasets, and self-calibrated together.

Observations at Band 7 were carried out on August 15, 2017 in the C40-7 configuration, using 42 antennas with a maximum recoverable scale of 0.″935 and longest baseline of 3.6km. The on-source time on HD 163296 was 28 min, J1733-1304 was used as flux calibrator, J1751-1950 as phase calibrator and J1924-2914 as bandpass calibrator. The correlator was set with four spectral windows with 1.875 GHz bandwidth, centered at 330.588 GHz, 329.331 GHz, 342.883 GHz and 341.000 GHz. In order to recover the flux at the shorter spatial frequencies, this dataset was combined with lower angular resolution data (project 2013.1.00053).

Finally, we used the calibrated Band 6 observation of HD 163296 made available by the DSHARP collaboration (Andrews et al. 2018a), with three spectral windows for the continuum centered at 232.6, 245.0, and 246.9 GHz. The high resolution DSHARP observations were combined with two additional datasets from previous ALMA programs, for a total baseline coverage spanning from 20 m to 5.8 km (configurations C34-4, C34-7, C40-8), see also Table 2 and 3 in Andrews et al. (2018a). The final beam size is 0.″048 × 0.″038 and maximum recoverable scale 13'″.

2.2 VLA observations

HD163296 was observed with the VLA at Ka band in the A and B configurations during semesters 16B and 17B, respectively. The bandwidth covered the range between 29.04 and 36.96 GHz (corresponding to wavelengths of 8.1-10.3 mm). 3C286 was used for absolute flux and bandpass calibration, J1820-2528 for pointing and J1755-2232 for complex gain (amplitude and phase) calibration. In the more compact configuration the maximum baseline was 11 km and the total time on source was 2 h 40 min. In the extended configuration the baselines ranged over up to 36.6 km, with a total time on the science target of 7.5 h. The datasets from the two configurations were first calibrated using the VLA pipeline script in CASA 5.6.1, epoch 17B was then self-calibrated with one round of phase calibration improving the signal-to-noise by ~30%. We combined the two configurations after scaling the flux of the A configuration, that resulted about 10% lower than the B configuration (this was expected as the disk emission is likely larger than the maximum recoverable scale of ~ 1.″6 in the extended array, so that spatial filtering can occur).

Additional observations were carried out in the X band (frequencies from 8 to 10 GHz, corresponding to ~3 cm) in the A configuration with maximum baseline 36.6 km, on October 1, 2016 and January 19, 2017 with a total of 1 h 50 min on the science target. The flux and bandpass calibrator was 3C286, the complex gain calibrator was J1820-2528 and After calibration using the CASA pipeline, we applied one round of self-phase calibration improving the S/N of about 20%.

2.3 Bandwidth and time smearing

For each of the ALMA and VLA datasets, a partial averaging in frequency and time was applied after the calibration. In order to avoid significant effects of chromatic aberration (also called bandwidth smearing, that is a radial smearing of the intensity distribution before it is convolved with the beam), we calculated the frequency bins corresponding to a reduction on the peak response of 1% at the edge of the primary beam for a point source, as derived in Mangum 2016¹ (public NRAO documentation). We then applied the corresponding frequency averaging in each dataset with the task mstransform and using the closest integer number of channel bins that corresponded to the derived ∆v.

Similarly, an excessive time binning can produce smearing of the intensity, but in the azimuthal direction. We estimated the time bins to apply in order to have a reduction of the intensity peak up to 1% anywhere in the image, writing the reduction of the peak response as $$R_{\text{a}}\simeq 1-\frac{1}{3}{\left(\frac{0.83{\omega }_{\text{e}}{\tau }_{\text{a}}}{{\theta }_{\text{b}}}\right)}^2\left(l_1^2+m_1^{2}si{n}^2{\delta }_0\right)$$

as derived in Sect. 6.4 of Thompson et al. (2001). We chose a time bin τ_a to have a peak reduction of 1% (R_a = 0.99), using the synthesized beamwidth θ_b, angular velocity of Earth ω_e, image plane coordinates (l₁,m₁) and declination of the source δ₀.

2.4 Systematic errors

The final uncertainty of the flux density measurements is mostly given by the accuracy of flux calibrators measurements: as they are in most cases quasar type objects, their intensity is intrinsically variable and this makes the evaluation of a flux density more difficult. Measurements of such calibrators are taken periodically by ALMA, but they are observed more frequently in certain bands (e.g. Band 3 and Band 6), while for other frequencies often an extrapolation of the flux from another wavelength is necessary (it is important to note that not only the intensity but also the spectral slope of a quasar can vary with time). Therefore the resulting calibration uncertainty will vary between ALMA Bands, and will depend on how close in time the calibrator measurement that is taken as reference was performed. Often a conservative calibration error of 10% is used for ALMA observations, but the analysis of ALMA calibrator measurements spanning several years has shown that a value of 5% is consistent with flux differences measured for short time spans in Band 3 and Band 6 (Bonato et al. 2018). More recently, Farren et al. (2021) compared ALMA and Planck observations and suggested that the calibration accuracy for ALMA Band 3 should be taken as 2.4% rather than the nominal 5%.

As it was pointed out that ALMA flux calibration accuracy could be poorer than the nominal value when the calibrator

catalog is not up-to-date (Francis et al. 2020), we check the individual calibrator measurements that were used as a reference during the calibration of our ALMA datasets. We find that the all the flux calibrations rely on recent measurements (taken within a few days of the science observations) of the corresponding calibrators, with deviations lower than than 3–4% (see Appendix A). Therefore we adopt the nominal flux uncertainties for the ALMA bands, corresponding to 5% at Band 3 and Band 6, and to 10% at Band 4 and Band 7.

For the VLA observations at the Ka and the X band, 3C286 was used as flux calibrator, known to have flux densities and spectral index constant in time. Based on the indications by NRAO², the single-epoch absolute accuracy is 10% for the Ka band and 5% for the X band.

Fig. 1 Continuum maps of HD 163296 from ALMA observations (top row) and VLA (bottom row). The synthesized beams of the final maps are drawn as white ellipse at the bottom left corner of each panel and are listed in Table 1. The images are available in electronic form at the CDS as FITS files.

2.5 Continuum maps

In Fig. 1, we show the continuum maps of HD 163296 analyzed in this work, with wavelengths from 0.9 mm to 3.4 cm. All the observations come from new datasets at high resolution, with the exception of the image at 1.3 mm taken from the DSHARP survey (Andrews et al. 2018a) and presented in Isella et al. (2018). The datasets at wavelength >3 mm are spanning a larger bandwidth in terms of Δv/v, so they have been split into high and low frequencies as the fluxes can have a non-negligible variation within the baseband (an average difference across the spatial frequencies of 35% is measure for ALMA Band 3 between Upper and Lower Side Band, 13% between high and low frequencies of VLA Ka band, 10% between high and low frequencies in VLA X band). The images were then produced with the CASA task tclean in multifrequency mode with nterms = 2. The properties derived from each image are listed in Table 1. The rms

noise was derived as the standard deviation of the flux in the signal-free regions of the images. The integrated fluxes where derived through aperture photometry inside ellipses of 46° inclination and 133° position angle with steps corresponding to ~1 beam FWHM. The final flux correspond to the aperture at which successive variations of the flux remain within a 3 σ level, where σ is calculated as the standard deviation of the flux inside signal-free regions of the same area as the aperture.

3 Methods

The main scope of this work is to constrain the dust properties in the HD 163296 disk using the high-resolution observations of the continuum emission presented in Sect. 2. We describe in this section the methodology used for the determination of the non-dust contribution, the extraction of the flux profiles, and the modeling setup for the spectral analysis.

3.1 Contamination from ionized gas emission

At radio frequencies, the emission from young stellar objects contains not only the thermal continuum from dust grains in the disk/envelope, but includes additional contributions that are often associated with free-free emission from ionized gas in the close surrounding of the star (e.g. Pascucci et al. 2012).

The exact origin and therefore spatial extent of such emission in Herbig and T-Tauri stars is still not clear: free-free electrons can be generated in different environments, such as an envelope produced by either mass loss or mass accretion, ionized wind from the disk's atmosphere (i.e., photoevaporation) or collimated jets/ouflows from the star. An additional nonthermal process that produces a radio continuum is gyrosynchrotron emission from flares in the corona of magnetically active stars (Güdel et al. 2002). The spectral slope at cm-wavelengths can help in identifying the responsible mechanism: emission from an optically thick, expanding shell at constant velocity is characterized by a slope of 0.6 (Panagia & Felli 1975). Free-free from disk atmospheres is expected to have a spectral index α between −0.1 and 2 for the optically thin and optically thick case respectively (Ubach et al. 2017). Nonthermal emission from corona flares have been observed to have a slope of around 2 from solar and stellar studies (Giidel et al. 2002). Surveys at cm-wavelengths in nearby star-forming regions seem to indicate a correlation of the spectral index α with the evolutionary stage, consistent with free-free emission dominating in early type YSOs (Class 0 to Class II, getting more optically thin toward the Class II stage), and gyrosynchrotron emission dominating in Class Ills (Dzib et al. 2013; Liu et al. 2014).

This contamination in the HD 163296 system was estimated in previous works, and the slope of the free-free emission was found to be ~0.6, using integrated fluxes at long wavelengths (Natta et al. 2004) and ~−0.2 using VLA resolved observations (Guidi et al. 2016). Both studies indicated that the contribution is negligible at wavelengths shorter than ~7 mm.

In Fig. 2, we plot the Spectral energy distribution in the mm/cm range for HD 163296: the change of slope of the SED when approaching cm wavelengths is clearly visible, and hints to the transition between different mechanisms responsible for the emission.

To get an estimate of the dust contribution at the long wavelengths, we can fit only the ALMA integrated fluxes (λ ≤ 3.3 mm) with a single slope, obtaining an integrated spectral index α = 2.6 ± 0.1 (blue dotted curve in Fig. 2). We note that such functional form seems to overpredict the fluxes at λ ~ 9 mm, that we know contain a fraction of non-dust emission. We demonstrate later that a single α slope is not a good description of the dust emission, as the spectral index tends to get artificially lowered at shorter wavelengths by the high percentage of optically thick emission. For this reason, and given the higher resolution of the datasets presented in this work, we get an estimate of the nonthermal contribution based on the resolved data themselves, instead of extrapolating the dust flux from the integrated SED.

We analyze the datasets at wavelengths as short as 3 mm, and we fit for a constant emission in the visibility plane using the highest spatial frequencies (smaller spatial scales), in the assumption that the free-free is a point-like source. The fits for each dataset are shown is Appendix B, and the estimated free-free contributions are displayed in Table 2.

As a further step, we employ a different method to get upper limits on the non-dust emission and check the consistency with the point-source approach. We produce images from the VLA datasets at the highest possible resolution, and extract the flux of the deconvolved models within a radius of 5 au. These measurements will contain both the entire free-free contamination and the dust emission, and are displayed as green triangles in Fig. 2.

The procedure and the motivation of the choice of 5 au is presented in detail in Appendix B, along with the model images.

Based on our estimates in Table 2, we derive a free-free spectrum as a power-law as a function of frequency, that results y^0.11±0.17 (dashed line in Fig. 2). The fit is done with a least-squares method (using the python routine scipy. optimize. curve_f it) and the error on α is defined as the square root of the variance as estimated from the fit. We can observe that the value of 0.11 for the spectral slope is consistent with optically thin free-free emission from a stellar or disk wind (Pascucci et al. 2014, and references therein), while it seems to rule out the hypothesis of gyrosynchrotron emission for which an α_cm around 2 is predicted. This would be also in agreement

with what reported for Class II systems in nearby star forming regions (Dzib et al. 2013).

We note that the assumption of a constant spectrum across cm wavelengths can be incorrect: models of free-free emission from disk winds ionized by X-rays or EUV predict for example a variation of the spectral index with frequency in the 1 to 100 GHZ range (Owen et al. 2013). Indeed, the intra-band spectral index that we measure at the peak of our VLA multifrequency images, where the emission is likely dominated by the free-free (Appendix B) is different for 33 GHz (α = 0.92 ± 0.07) and 10 GHz (α = 0.24 ± 0.03), where the low uncertainties are due to the fact that they consist of a relative measurement within the same frequency band and therefore not affected by absolute amplitude calibration errors.

While it is not possible at the moment to link such spectral indexes to a specific underlying mechanisms of the free-free emission, they can provide some useful reference for future studies of photoevaporation/accretion mechanisms in this system, subject that is beyond the scope of this paper.

Fig. 2 Spectral energy distribution of HD 163296. The blue circles represent the integrated fluxes from the data analyzed in this work, while empty squares are values from the literature (Isella et al. 2007; Natta et al. 2004; Guidi et al. 2016). The blue dotted curve corresponds to a curve with spectral index of 2.6. Red stars are the contamination from a central compact emission, determined by the asymptotic values of the visibilities as described in Sect. 3.1, and that are best fitted by a power law with index 0.11. Green triangles are upper limits for the free-free emission (see Sect. 3.1), and correspond to the emission within a radius of 5 au.

Fig. 3 Radial profiles of the continuum emission derived with frank for each dataset, the shaded regions represent the statistical error computed by the frank fit. The star symbol in the plots relative to the VLA datasets indicates that these have been first corrected for the free-free/centrally peaked contamination, as described in Appendix C.

3.2 Extraction of radial profiles

The more complex the sky brightness distribution, the less trivial is the image reconstruction from interferometric measurements through aperture synthesis, in our case by means of the CLEAN algorithms in CASA. Besides requiring some initial assumptions on the clean components model (e.g. multiscale cleaning), the image reconstruction requires a deconvolution, which is a nonlinear operation that can introduce similar nonlinear alteration on the flux distribution. In order to mitigate the uncertainties of aperture synthesis and bypass the deconvolution operation, we additionally analyze the observations in the visibility (or Fourier) plane. It has been recently showed that the resolution of interferometric datasets can be “pushed” beyond the limit obtained with synthesized imaging, by analysing the data in the visibility space (see e.g. Tazzari et al. 2018; Jennings et al. 2021). By fitting the visibilities of the DSHARP observations at 1.3 mm with an axisymmetric model, Isella et al. (2018) characterized the continuum emission in HD 163296 with 5 Gaussian rings and highlighted two minor asymmetric features: one at about 4 au and a crescent in the south-east direction at about 55 au separation. This latter is partially visible also in our images at 0.9, 2 and 3 mm displayed in Fig. 1.

We perform a similar analysis in this work, but we use the python package frank (Jennings et al. 2020), that allows the recovery of brightness profiles without assuming any functional form for the intensity as a function of the radius. The main assumption is always that the brightness is axisymmetric, therefore only the real part of the visibilities is modeled and the Fourier transform can be simplified as an Hankel transform. We refer to Jennings et al. (2020) for an exhaustive description of the metodology and convergence criteria of this tool, and to Appendix C for the details of our modeling. The extracted brightness profile for each dataset is displayed in Fig. 3.

The VLA datasets have been corrected for the strong central emission likely due to free electrons (see Sect. 3.1), to avoid artificial oscillations in the brightness profiles (see Appendix C).

In Table 3, we show the position of the peaks in the radial brightness profiles at each wavelength ranging from 0 to 1725, as found with the python routine scipy.signal.find_peaks. The most pronounced ringed-structure is shown by the ALMA Band 4 dataset (2.1 mm) where we can clearly identify 6 pronounced peaks, i.e. one more than what found by previous studies (Isella et al. 2018). This likely results from this datasets having the best combination of angular resolution and intrinsic width of the peaked emission. In fact, the DSHARP Band 6 dataset that has the best resolution shows less pronounced peaks in the inner disk with respect to Band4. This can be either an optical depth effect, or could reflect the different radial distribution of smaller dust particles. We use the radii identified for Band 4 as reference for the ring positions in the analysis carried out in this paper.

Worth noticing is also the difference in the relative peak intensity of the outer rings: while the intensity at the 67 au ring is much higher than the one at the 100 au for the shortest wavelengths, this difference tends to disappear as the wavelength increases. In the next section, we show that this is due to the higher optical depth in the outer ring that is also increasing with frequency.

For the purpose of the multiwavelength analysis described in Sect. 3.3, we need to compare the emission at the same resolution at all wavelengths. To make sure this is verified, we perform a second round of fits with frank, where we truncate the visibility distribution for the ALMA tables at shorter spatial frequencies, in order to obtain the same accuracy B₈₀ = (2100 ± 100) kλ from 0.9 mm up to 9 mm (see Appendix C for the details). The 30 mm dataset has a considerably lower resolution (see Table C.1), and since degrading the spatial resolution at the shorter wavelengths to match the ~0.″3 of 30 mm would result in a consistent loss of information, we do not include the 30 mm observation in the multiwavelength analysis described in the next section.

3.3 Spectral analysis setup

Resolving the disk structure at multiple wavelengths can help us constraining the dust properties as a function of the radius. We use the brightness profiles extracted with frank at a matched resolution (see Sect. 3.2) to fit the Spectral energy distribution with three different analytic models for the dust emission. We start by using a simple parametric model where we describe the optical depth as a power law in function of the frequency. We then employ a physical model that includes only the absorption opacity for the dust, and finally we introduce the contribution from scattering.

In the parametric model, we assume that the intensity emitted from dust is regulated by an opacity with a power-law dependency from the frequency. Under the assumption of LTE (Local Thermodynamic Equilibrium), the dust thermal emission can be written as $$I_v\left(r\right)=B_v\left(T\left(r\right)\right)\times \left({\text{1e}}^{-{\tau }_v\left(r\right)/\mu }\right),$$(1)

with the optical depth given by τ_v = τ₀(v/v₀)^β, and the inclination parameter µ = cos(i). At a given radius, the intensity in Eq. (1) depends on three parameters: the midplane temperature T, the optical depth τ₀ (at a reference frequency v₀), and the spectral index β.

We use the python package UltraNest (Buchner 2021) a Monte Carlo Nested Sampling tool based on the MLFriend method (Buchner 2014), that computes both the posterior probability of the three free parameters and the marginal likelihood of the model. The convergence and termination criteria relative to this tool are described in detail in Buchner (2021), and references therein. We fit Eq.(1) independently at each separation, i.e. without assuming any trend as a function of radius, with bins of 2 au starting from 0 up to 120 au. We use a flat prior for all parameters, corresponding to log₁₀τ₀: [−4, 3], β: [0, 5] and T: [3, T_up], where T_up is an upper limit variable with radius (see Appendix D for the details). The likelihood function in this Bayesian framework is defined at each radius as the sum of the χ² values at the different wavelengths, calculated weighting the squared residuals by 1/σ² as ${\displaystyle {\sum }_{v_i}{\left(F_{\text{obs},\text{r,}{v}_i}-F_{\mathrm{mod},\text{r,}{v}_i}\right)}^2/{\sigma }_{\text{r,}{v}_i}^2}$, where σ is the error on the fluxes given by the statistical errors estimated in Sect. 3.2. As the measurements are affected by a systematic error on the order of 5-10% related to the flux calibration, we perform 30 independent fits where we scale the flux densities at all radii and at each wavelength by a random offset, generated from a normal distribution with standard deviation corresponding to the flux calibration uncertainty (10% for ALMA Band 7, Band 4 and VLA Ka band, and 5% for ALMA Band 3 and Band 6). The resulting posterior probabilities at each radius are obtained by merging the posteriors of the 30 fits.

To obtain dust physical properties such as the surface density and grain size, we need to employ a physical model that relates these quantities to the observed intensity. A necessary step in this direction involves the computation of the dust opacity as a function of the grain size.

The main caveat is that the exact constituents of the dust grains are not known, and different properties (especially composition and porosity) determine dramatic differences in terms of opacity and albedo of the dust (e.g. Min et al. 2016). An overview of the effects of different compositions assumed for grains in protoplanetary disks on their optical properties is given in Birnstiel et al. (2018). The result is that each analysis can lead to very different conclusions in terms of the quantities of interest (opacity, mass, grain size), depending on the initial assumption of dust composition/porosity. Furthermore, the dust in protoplanetary disks does not consist of a single-size population, but rather an ensemble of grains of different sizes. This is usually described with a power-law distribution dn(a)/da ∝ a^−q, with a representing the particle size and n the number of particles with a size a, between a minimum size a_min and a maximum a_max. The opacity of such an ensemble is mostly sensitive to the maximum grain size and it will depend on the size distribution power-law q. This latter is estimated by theoretical and experimental studies of collisional and coagulation processes in dust grains (Testi et al. 2014), resulting in typical values of q ~ 2–4. Often a value of q = 3.5 is assumed, following studies characterizing the interstellar dust (see e.g. Draine 2006; Testi et al. 2014, and references therein). Ultimately, the value of q for protoplanetary disks is not known as it is expected to vary with time evolution (see e.g. Testi et al. 2014) and location within the disk.

With these uncertainties in mind, we consider a set of different grain populations for a total of 8 initial models. We compute a SED fit for each model and we compare the evidence to assess which one is more representative of the dust population at each radius. We compute the opacities from the DIANA project (Min et al. 2016), through the DIANA Opacity Tool Fortran package. The authors define a standard grain composition for protoplanetary disks as a mixture of amorphous silicates (Dorschner et al. 1995) and amorphous carbonaceous materials (Zubko et al. 1996) in a volume fraction of 60% and 15%, respectively, with a remaining 25% of vacuum. In the opacity calculations, grains are modeled as distributions of hollow spheres (Min et al. 2005), overcoming the assumption of spherical grains used in the Mie scattering theory (Mie 1908). We generate a series of opacities for a range of wavelengths covering our observations, varying the maximum grain size from 10⁻⁴ to 10⁴ cm and keeping a minimum size of 0.05 µm. This procedure is repeated for four different slopes q of the size distribution, with values of 2.5, 3, 3.5, 4. The absorption and scattering opacities as function of the maximum grain size at 1.3 mm are plotted in Fig. 4, upper panel, and labelled as “standard”. We compute a further set of opacities enhancing the grain porosity to 80%, while keeping the same relative fraction of silicate and carbons as in the standard composition (Fig. 4, lower panel labelled “porous”). This results in eight different dust populations: two compositions (compact and porous) with four size distribution each.

The optical depth in Eq. (1) can then be written explicitly as function of the opacity and surface density: $$I_v\left(r\right)=B_v\left(T_{mid}\left(r\right)\right)\left[1-\exp (-{\kappa }_{v,abs}\left(a_{\max }\left(r\right)\right)\times {\Sigma }_{\text{d}}\left(r\right)/\mu \right]$$(2)

where K_abs is the absorption opacity and Σ_d the dust surface density. For a given dust composition and size distribution with spectrum q, we have that κ_v,abs depends only on the maximum grain size a_max (see Fig. 4). At a given radius r we still assume that the temperature is the same for the emission at all wavelengths. As for the parametric fit, we use the Monte Carlo nested sampling algorithm with the UltraNest software to fit Eq. (2) independently at radial bins of 2 au from the star, and using a flat prior of [3, T_up] for T, [−4, 3] for log₁₀ a_max and [0.0001, 10] for Σ_d.

Finally, to include the possible effects of dust self-scattering, we use the analytic expression for the emergent intensity given in Zhu et al. (2019), that is valid in the assumption of isotropic scattering from an isothermal slab: $$\begin{array}{l}I\left(r\right)=B_v\left(T\left(r\right)\right)\times \left(1-\exp \left(-{\tau }_{\text{d}}\left(r\right)/\mu \right)\right)\hfill \\ \left(1-{\omega }_v\left(r\right)\frac{\exp \left(-\sqrt{3\left(1-{\omega }_v\left(r\right)\right)}{\tau }_{\text{d}}\left(r\right)\right)+\exp \left(\sqrt{3\left(1-{\omega }_v\left(r\right)\right)}\left(\tau \left(r\right)-{\tau }_{\text{d}}\left(r\right)\right)\right)}{\exp \left(-\sqrt{3\left(1-{\omega }_v\left(r\right)\right)}{\tau }_{\text{d}}\left(r\right)\right)\times \left(1-{\omega }_v\left(r\right)\right)+\left(\sqrt{1-{\omega }_v\left(r\right)}+1\right)}\right).\hfill \end{array}$$(3)

The deprojected optical depth in this description is calculated as τ = 2µτ_d/(3τ_d + 1), from the total optical depth τ_d that includes both absorption and scattering, given by τ_d = Σ_d(κ_abs + κ_sca,eff) (Zhu et al. 2019). Here Σ_d represents the dust surface density, κ_abs is the absorption opacity, and the effective scattering opacity is obtained with a correction by the asymmetry parameter g as κ_sca,eff = (1 − g)κ_sca. Finally, the albedo ω is the ratio between the scattering and the total opacity κ_sca,eff/(κ_abs + κ_sca,eff).

We can observe that for a given dust composition and size distribution with spectrum q, we have that κ_abs, κ_sca and g at a certain wavelength depend only on the maximum grain size a_max. Therefore, as in the non-scattering case described by Eqs. (2), (3) depends only on T(r), a_max(r) and Σ_d(r). We recall here that the underlying assumption is that at a given radius r the temperature at the emitting layer is the same at all wavelengths. Following the same procedure as for the first two models, we find the best-fit parameters and relative uncertainties with ultranest, and we perform 30 additional fits at each radius to account for the systematic calibration offset.

Fig. 4 Dust opacity at 1.3 mm as a function of the maximum grain size. Upper panel: absorption (solid lines) and scattering (dashed lines) opacities for the standard dust composition and different size distribution q. Lower panel: same as the upper panel but for grains with a porosity of 80%.

3.4 Model comparison with Bayes factor K

We take advantage of our Bayesian framework to compare the performances of the different analytical models we use to fit our observations of HD 163296. The Monte Carlo nested sampling routine we apply provides not only the posterior probabilities but also the Bayesian evidence (or marginal likelihood) Z. This corresponds to the normalization factor in the Bayes equation, or the integral over the whole parameter space of the likelihood times the prior density Z = ∫ L(D|θ)π(θ)dθ. It follows that we can compute the Bayesian factor K between different models, defined as the radio of the marginal likelihoods, and assess which model is more compatible with our datasets. To interpret the Bayes factors, we refer to the scale proposed originally by Jeffreys (1939),

that associates different ranges of the K factor to a strength of a change in evidence. Defining K₁₂ as Z₁/Z₂, where Z₁ and Z₂ are the marginal likelihood of model 1 and 2, respectively, we use the following scale:

Fig. 5 Flux density spectral index. Left: map of the spectral index computed from the ALMA continuum images at 0.9 mm, 1.3 mm, 2.1 mm, 2.8 mm and 3.2mm with a matching beam of 0.″095 × 0.″065 and position angle of 61.24°. Right: radial profile of the spectral index in the northwest and southeast sides (blue and red curve, respectively), error bars are shown for each bin and calculated as described in the text. Overplotted with dashed lines is the corresponding flux density profiles in Band 4 (right y axis). The grey shaded region on the left denotes the angular resolution of the map as half of the average beam FWHM.

4 Results

4.1 Spectral index

The measure of the flux spectral index of the dust emission at millimeter wavelengths has been the primary tool for deriving information on the grain properties, through its relation with the dust opacity spectral index β (where κ_v ∝ v^β). Multiple surveys of protoplanetary disks allowed in the past to estimate integrated values of β - in the Rayleigh-Jeans and optically thin assumptions – that resulted systematically lower than the values measured for the insterstellar medium β_ISM ≃ 1.7 (e.g. Beckwith et al. 1990; Testi et al. 2003; Natta et al. 2004; Ricci et al. 2010). This trend was confirmed by more recent ALMA surveys targeting disks in nearby star forming regions (e.g. Ansdell et al. 2018; Tazzari et al. 2021).

A common explanation for these measurements (typically β < 1) was found in the growth of solids, since other mechanisms such as chemical composition of the grains and porosity, are expected to have only moderate effects on the total opacity, within certain limits of the grain size distribution. Specifically, when such distribution follows a power law of the form dn/da ∝ a^−q, a size distribution with q = 3.5 will have β ≲ 1 for a_max ≳ 3λ (Draine 2006).

In the ALMA era, more accurate measurement of the spectral index are possible: by spatially resolving the continuum emission from the disk, the radial variation of the spectral index can be computed. Studies employing medium-resolutions ALMA observations showed for different disks a monotonically increasing spectral index from small to larger radii: this was interpreted as a signature of larger grains in the inner regions, as expected from the differential action of radial drift (e.g. Pérez et al. 2012, 2015; Tazzari et al. 2016; Guidi et al. 2016). More recently and thanks to higher spatial resolution observations (≤20 au) it has been possible to measure the spectral index variations on smaller scales and show that it deviates from a pure monotonic behavior in several disks, such as HLTau (Liu et al. 2017; Carrasco-González et al. 2019), TWHya (Macías et al. 2021), HD 169142 (Macas et al. 2019), GMAur (Huang et al. 2020) HD 163296 (Dent et al. 2019; Ohashi & Kataoka 2019) and a few disks in the Taurus association (Long et al. 2020). These improved measurements are showing that a low α in the millimeter range does not necessarily coincide with the presence of large (millimeter) grains, but it is in some cases artificially lowered by the high optical depth of the continuum emission.

With the extended and improved set of ALMA observations of HD 163296 we present in this work, we can build a detailed map of the spectral index between 0.9 mm and 3.2 mm and measure the spectral index with higher resolution and smaller uncertainties compared to previous works.

After producing images with a matching beam of ~0.″08 and centered on the same pixels, corresponding to the lowest resolution available (in this case the 3 mm observations), we can compute the spectral index α of the flux (assuming F_v ∝ v^α), with the least-squares method for each pixel. The resulting α map and corresponding radial profile are shown in Fig. 5: the profiles are computed inside segments within a 45° (deprojected) angle centered on the disk major axis (PA =133°) with bins taken every half-beam (4 au) and treating the NE and SW side separately. For each bin we computed the weighted mean of the values, with weights given by the inverse squared error of each pixel σ_α, derived analytically from the linear least-squares regression. The uncertainties on the flux densities at each pixel were taken as ${\sigma }_{\text{F}}=\sqrt{rm{s}^2+{\left(\text{Δ}{F}_{\text{cal}}\right)}^2}$, i.e. the sum in quadrature of the rms of the image and the flux calibration error, corresponding to 5% or 10% depending on the wavelength (see Sect. 2.4). The error on

each bin is then computed as the standard error of the weighted mean, accounting for the fact that the pixels are not independent.

The spectral index map appears overall azimuthally symmetrical and the profiles do not show strong deviations between the south-east and north-west sides; some appreciable differences are observed inside the gaps – where we are dominated by the noise – and at the location of the crescent at 55 au on the south-east side (red curve in Fig. 5): here the spectral index is consistent with what is measured in the adjacent ring at 67 au (α = 2.8 ± 0.1). Interestingly, the outermost ring at 100 au shows a value of α = 2.4 ± 0.1, lower with respect to the ring at 67 au. This is in agreement with previous recent measurements of the millimeter spectral index at lower resolution (Dent et al. 2019; Sierra et al. 2021). We compare in Fig. 6 the spectral index obtained from the images with the one computed from the radial profiles (see Sect. 3.2). While there is generally a very good agreement between the two profiles, we note how the fluxes extracted with frank allow us to reveal with more detail the variation of the α index in the inner region of the disk (inside ~40 au).

Finally, we note that α reaches values lower than 2 in the innermost regions (r ≲ 10 au), i.e. below the black-body limit for an optically thick emission. This has been already observed in HD 163296 using a smaller set of observations (Dent et al. 2019), and in other sources (e.g. TWHya Huang et al. 2018a). This feature can be due to multiple effects, such as the self-scattering reducing the emission at shorter wavelength in the optically thick inner regions, or the free-free emission increasing the contribution of longer wavelength in the central flux (e.g. the contamination calculated in Sect. 3.1 accounts for ~20% of the flux at 3 mm in the central beam, and reduces to about 1% at 0.9 mm).

Fig. 6 Flux spectral index between 0.9 mm and 3.2mm computed from the radial profiles extracted with frank (Sect. 3.2, black curve), with shaded areas representing the errors derived from the linear least-squares regression (with the uncertainties for the single fluxes defined as the sum in quadrature of the statistical and calibration errors). The profile is masked at the locations where the flux S/N < 3 for at least one wavelength. Overplotted is the azimuthal average of the map in Fig. 5 (green curve), calculated with bins of 1 beam (~0.″08) and with errorbars computed as for the profiles in Fig. 5, right panel.

4.2 Parametric model

We show in Fig. 7, the best-fit parameters of our parametric model described in Sect. 3.3 as a function of the radius;

the results are plotted starting from 8 au as for lower separations some of the parameters resulted highly unconstrained (see Appendix D). The parameters estimates from our nested sampling method that includes the statistical error are shown with a dashed red line, and we overplot the total normalized posterior obtained merging the 30 fits with a random calibration offset, as described in Sect. 3.3.

The results indicate that the 1.3 mm emission is moderately optically thick: while the average value in the inner disk inside 35 au is τ_1.3mm ≃ 1.3, we find ${\tau }_{1.3\mathrm{m}\text{m}}={0.8}_{-0.2}^{+0.2}$ at the 66au ring, and ${\tau }_{1.3\text{\hspace{0.17em}mm}}={1.5}_{-0.3}^{+0.3}$ in the outer ring at 100 au, where the uncertainties are given as the 16th and 84th percentile of the total posterior distribution. The temperature is well-constrained across the disk except in the dust gaps, where on the contrary the optical depth is lower and the temperature uncertainties are higher. Nevertheless, the temperature profile obtained with this approach deviates from a smooth decreasing power-law (the functional form that is generally assumed for the radial dependence of T when analysing medium-resolution observations), and in particular it shows enhanced values in correspondence of the dust gaps. This temperature increase in the dust gaps of HD 163296 was already pointed out by van der Marel et al. (2018) and Rab et al. (2020), and is typically explained with the higher penetration of the scattered light from the disk surface into the midplane, because of the dust depletion. This is also consistent with the effect of a planet on the disk temperature structure as shown by hydrodynamic simulations (e.g. Isella & Turner 2018). From our modeling, an increase in temperature is appreciable in the large dust gaps at about 50 and 85 and 115 au, with temperature peaks corresponding to 2.5, 2.7 and 3.7 times the values at the closest inner rings (32, 66 and 100 au, respectively). A tentative increase in T is also observed in the innermost gap at ~10 au.

Beyond 30 au, the β index seems consistent or larger than the ISM value, and no difference in β is measured between the ring at 66 au (β = 1.9 ± 0.3) and the ring at 100 au (β = 2.0 ± 0.3), i.e. no indication of larger grains at 100 au, as the spectral index shown in Sect. 4.1 might have suggested. The lower α can be explained by the higher optical depth of the outer ring as mentioned above.

Fig. 7 Best fit parameters for the power law model: temperature, optical depth, and opacity spectral index β as function of the radius. The dashed red curve represents the estimates obtained as the mean values of the posterior distributions from the Monte Carlo nested sampling fit. The color map is the normalized probability after merging the 30 posteriors obtained introducing a random offset in the fluxes according to their flux calibration accuracy. The dashed-dotted curve in the upper panel shows the upper limit of the Temperature prior for each radius used in the fit (see Appendix D). The vertical dashed lines correspond to the dust ring found in Sect. 3.2, while the white dashed line in the bottom panel is drawn for β = 1.7, corresponding to ISM dust grains.

Fig. 8 Dispersion of the best-fit temperature, maximum grain size and surface density in the non-scattering model, calculated as the mean value at each radius of the four different size distributions (the single best-fit models are show in Appendix D), considering standard and porous grains separately. The shaded regions correspond to the standard deviation of the four best-fit values at each radius.

4.3 Physical models

We fit Eq. (2) at each radius for all the 8 dust models, and we show in Fig. 8 the best-fit values averaged on the four size distributions (q ranging from 2.5 to 4) for both the standard and porous dust grains. In Appendix D we show the results for all the single models. We observe that while the temperature profile is consistent between the two sets of dust, the maximum grain size for porous grains is on average a factor of 3–6 higher (with a ratio increasing with the size distribution spectral index q), and the surface density a factor of 5 higher (roughly consistent across the size distributions). Similarly, in the model that includes scattering (Eq. (3)), assuming porous grains we predict a a_max between 2 and 4 times larger (again increasing with the q spectral index) and a surface density that is on average ≈5 times larger (see Fig. 9).

To determine the most probable model between the two different dust compositions (standard and porous), we look at the Bayesian evidence of the Monte Carlo nested sampling fits and identify the models with the highest evidence. We illustrate the results for the non-scattering and scattering case in Fig. 10, showing the Bayes factors between the standard and porous models for each value of q. Referring to the scale reported in Sect. 3.4, we find that overall standard grains seem to better reproduce our data, with K > 10 (strong evidence for standard grains) in 61% of the cases (cases corresponding to the number of q values times the number of separations), while for porous grains we have K < 1/10 only in 16% of the cases. Moreover, standard grains seem to have a decisive superior evidence (K > 100) at 57% of the radii, while this goes down to 5% for porous grains. These values are calculated over all the four size distributions, but they present similar values when looking at the single q models. In the scattering case we find similar results,

with a strong evidence for standard grains 55% of the times (decisive 50% of the times), and 15% for porous grains (decisive for 3%). Interestingly, porous grains result a better model only in correspondence of the gaps in the dust distribution (although for the few radii in the inner disk where the rings are not resolved, this is less clear).

Since we determined that the models with standard grains are preferred to the ones with porous grains, we focus on the former to infer the final estimates of our best-fit parameters for the HD 163296 disk. The further step consists in estimating the best size distribution as function of the radius, following the same procedure based on the K factor. In Fig. 11, we plot the Bayes factor K relative to the q model with the larger evidence: this shows simultaneously the preferred size distribution at each radius and the strength of the evidence ratios. We show the corresponding final best-fit parameters, relative to the best size distribution at each radius, in Figs. 12 and 13 for the nonscattering and scattering model, respectively. We draw with a red curve the best-fit values for the Temperature, maximum grain size and surface density obtained from the Monte Carlo fit with statistical error only, and the total posterior distribution obtained with the 30 additional fits to account for the calibration error as a color map.

The best-fit temperature is consistent with the one derived from the simple power-law model (Fig. 7), with hints of higher values of T in the dust gaps at about 50, 85 and 115 au. In addition, in both the nonscattering and scattering model we see a steep increase in the temperature in the innermost gap, with T ≃ 167 K and 180 K respectively. If we compare our midplane temperature profile with the one derived by Dullemond et al. (2020) using CO emission lines (gray curve in the upper panels of Figs. 12 and 13), we note that while at the 66 au ring there is only a small difference between the two studies, at the 100 au ring we get a significantly lower T (≃12K) in both the nonscattering and scattering case. Since this is below its condensation temperature, the CO would be frozen-out onto dust grains at this location in the midplane. This invalidates the assumption made by Dullemond et al. (2020), that would find a higher temperature as their measured CO emission is coming from higher vertical layers at this separation.

The maximum grain size radial profile is interestingly flat outside ~40 au, roughly consistent with a constant a_max ≃ 200 µm from 40 to 120 au with no significant difference between

gaps and rings. While the grain size is well constrained in the outer rings (at r > 40 au), in the inner disk (in particular between 15 and 30 au) a degeneracy in a_max results in higher uncertainties on this parameters, with values ranging from 10⁻² to 10² cm. The surface density drops, i.e. indicates dust depletion, at the same locations where the temperature increases (see previous paragraph).

Comparing our result with the work by Isella et al. (2016), that performed radiative transfer modeling using a smooth surface density profile with rectangular gaps, we find that outside 20 au the two surface densities are consistent within a factor of 2. A significant exception is the ring at 100 au, where we find a higher Σ_dust by a factor of 7. The optical depth at the five wavelengths from our scattering model is shown in Fig. 14. The errors are estimated by taking the 2.5/16th and 97.5/84th percentiles of the distribution of τ obtained drawing 1000 random samples from a Gaussian distribution of a_max and Σ_d (with σ corresponding to the standard deviation of their total posterior distribution at each radius). In the outer rings at 66 and 100 au, the emission is still moderately optically thick at a wavelength of 1.3 mm, with τ of order unity. It is worth noticing that the optical depth in the outer ring at 100 au results larger than the one at the 66 au ring at all wavelengths, with τ1₀₀ ~ 1.7τ₆₆. At the higher frequencies

this is even comparable with the optical depth of the inner rings, with ${\tau }_{0.9\mathrm{m}\text{m}}={4.6}_{-1.1}^{+1.6}$ and ${\tau }_{1.3\mathrm{m}\text{m}}={17}_{-0.3}^{+0.4}$ at 100 au, with uncertainties corresponding to ~ 1σ (16th and 84th percentiles) of the posterior distribution and showed as shaded areas in Fig. 14.

In Fig. 15, we plot the albedo at the different wavelengths in correspondence of the rings and at the innermost gap at 10 au: while in the outer disk the albedo decreases with wavelengths, in the inner disk (r ~ 10 au) it shows the opposite trend. Across the intermediate rings (at 14, 24 and 32 au) the trend appears similar to the innermost part of the disk, but the uncertainties on the albedo are too high to draw a robust conclusion. This “spectral inversion” of the albedo from the outer to the inner disk is related to the transition from small grains in the outer rings (ω_{1 mm} ~ 0.3) to large grains inside (ω_{1 mm} ~ 0.7), with the peak of the albedo falling at wavelengths λ ~ 2πa_max.

Fig. 9 Same as Fig. 8, but for the model that includes scattering (the single best-fit models are show in Appendix D).

Fig. 10 Bayes factor K calculated between the standard and the porous models as Zstandard/Zporous, for the nonscattering model (upper panel) and scattering model (lower panel). The different markers correspond to the different size distribution coefficients q. The points are color-coded according to their values: K > 1 (larger evidence for the standard composition) are drawn in cyan and K < 1 (larger evidence for porous composition) in magenta. Empty markers correspond to 1 < K < 3 or 1/3 < K < 1, full markers to K > 3 or K < 1/3. The horizontal dashed line is drawn at K = 10, so that all points above this line correspond to a strong evidence in favor of standard grains. The points larger than 100 are drawn at the location of 100, since for K ≥ 100 the interpretation in terms of evidence strength does not change (see Sect. 3.4).

Fig. 11 Bayes factor K computed between the model with the highest evidence (identified by the color of the marker) and the remaining 3 size distributions at each radius, for standard grains in the non-scattering (top panel) and scattering (bottom panel) case. The full color indicates that the Bayes factors K at that specific radius are all >10, i.e. the is a strong evidence for that size distribution compared to the other three discrete values.

5 Discussion

5.1 Grain size

Previous works that analyzed mm-observations of HD 163296 at low and moderate spatial resolution indicated the presence of ~millimeter/centimeter sized grains (Beckwith et al. 1990; Natta et al. 2004; Isella et al. 2007; Guilloteau et al. 2011; Guidi et al. 2016) based on the measured millimeter flux spectral index. However, in the recent years the robustness of a straightforward interpretation of this parameter as a direct proxy for the grain size has been disputed: observational evidence suggest that the assumption of optically thin emission at millimeter wavelengths is no longer justified for protoplanetary disks, and that dust self-scattering can therefore play a significant role in regulating the emitted intensity (e.g. Kataoka et al. 2016; Liu 2019; Zhu et al. 2019; Ueda et al. 2020). Accounting for this effect, several multiwavelength studies have been carried out for other bright disks (Huang et al. 2020; Maclas et al. 2021; Sierra et al. 2021), and indicate that in fact the optical depth and albedo can assume high values especially in the inner disk, but also at large separation from the central star (e.g. HLTau Carrasco-González et al. 2019).

Similarly to the studies reported above, in our multiwavelength analysis of the continuum emission we include self-scattering from an isothermal slab within an analytical formulation of radiative transfer. Using this physical model we could successfully fit the observed profiles of HD 163296 disk deriving the dust temperature, maximum grain size and surface density at separations larger than 8 au. At smaller separations, the surface density could not be constrained, likely because of the high optical depth that results in the Planck term in Eq. (3) dominating the emitted intensity. The fact that at shorter radii we could not constrain the surface density and maximum grain size could indicate that the underlying assumption of emission from an isothermal surface at all wavelengths does not hold in the very inner disk. The best-fit model indicates that the outer rings (at 66 and 100 au) are composed by grains on the order of 200 µm, with no significant variations in grain size between these rings and the adjacent gaps. In the inner disk (r ≲ 40 au) the a_max is less well constrained: the total posterior shows a degeneracy in the grain size, that can take values from ~200 µm to ~ 100 cm, and the best estimates in this region are more dependent from the size distribution spectral index q (see Fig. D.2). Despite these localized degeneracy, we note that a_max overall shows larger values in the inner disk (e.g. >millimeter-size in the rings at ~16 and 24 au, see Fig. 13). We note that this degeneracy is not present in the outer disk where we spatially resolve the rings at all wavelengths.

In the outer rings at 66 and 100 au, the solutions are strongly dominated by the steeper size distributions (q = 4) in both the non-scattering and the scattering model, which corresponds to the bulk of the dust mass contained in the small grains. On the contrary, a prevalence of flatter distribution with q = 2.5/3 is found at small separations (<40 au), but only for the scattering model. Although the explored grid of values is very limited and therefore we cannot determine the size distribution slope with high accuracy, we note that the steeper size distributions in the outer rings of HD 163296 are consistent with what found in the outer disks of HD 169142 (Macias et al. 2019) and TW Hya disk (Macias et al. 2021), by fitting multiwavelength observations with q as a free parameter.

In a recent study, Sierra et al. (2021) analyzed ALMA observations of HD 163296 at medium resolution (~ 19 au) to constrain the dust properties across this disk. Using a physical model that includes scattering (as in Eq. (3) in this work) and assuming a power-law size distribution with q = 2.5, they find two possible families of solution for the a_max parameter in the inner disk (r ≤ 40 au): one with grains on the order of 100 µm and one with millimeter grains, similarly to what found in this work. At larger separation they find a maximum grain size on the order of millimeter, with a local increase at the 100 au ring. These latter values differ by about one order of magnitude with our findings of ~200 µm grains outside 40 au. This is not entirely surprising, as we showed in this work how the initial assumptions, such as the dust composition and size distribution, can affect the derived parameters in these analyses. We note that these two studies differ for the choice of dust opacities and the fixed temperature profile that Sierra et al. (2021) assume while fitting only for the maximum grain size and surface density. To understand how each of these factors affects the final results, we performed some tests varying the dust opacities and introducing a fixed temperature profile. We describe the procedure and results in Appendix E, where we find that both the chosen composition and the fixed temperature profile significantly affect the final estimates of maximum grain size.

At this stage we still lack information on the typical dust composition in protoplanetary disks, therefore it is important to test a wider range of dust properties when fitting disk observations. In this direction, we note that Zormpas et al. (2022) recently found that dust opacities including such amorphous carbons from Zubko et al. (1996, the same used in this study) can better reproduce the size-luminosity relation observed in nearby star forming regions (Tripathi et al. 2017), with respect to the DSHARP opacities.

Another important measure of the dust properties in disks comes from polarization studies: ALMA polarimetric observations of HD 163296 indicated that the grain size across this disk is smaller than 100–150 µm, when interpreting the polarized emission in terms of dust self-scattering (Dent et al. 2019; Lin et al. 2019; Ohashi & Kataoka 2019). While our results are in agreement with dust polarization measurements in the outer disk (r ≥ 40 au), we find evidence for larger grain sizes in the inner disk. However, it was recently pointed out that a mix of dust scattering and magnetic alignment could be responsible for the detected polarized signal in the disk of HL Tau (Mori & Kataoka 2021). If this is the case for HD 163296, this would loosen the constraint on the maximum grain size of 100 µm in the inner disk and mitigate the discrepancy with our results.

The bimodal distribution of the maximum grain size between the inner and outer disk could have multiple explanations. As a consequence of the aerodynamical friction with the surrounding gas moving at sub-Keplerian velocities, dust grains lose angular momentum and drift toward smaller radii, with larger particles drifting inward more efficiently than smaller particles (Weidenschilling 1977). As a result we expect to find larger particles in the inner disk and smaller outside. If pressure traps are present in the disk, they could stop or slow down radial drift and retain some large particles within localized structures. Because of the poor constrains on a_max in the inner disk, we cannot tell whether this is the case for the inner rings (r ≤ 40 au). On the opposite, we find a more robust evidence of no differential trapping (no change in grain size between rings and gaps) in the outer rings at 66 and 100 au. This could mean that the timescales of radial drift are shorter than the ones relative to the mechanism that created the rings, or the rings could have form recently and not have enough time to trap the particles. Another reason could lie in the sticking properties of dust grains, hindering their growth outside the water snowline: new laboratory measurements are showing that water ice-coated grains have a lower sticking force than previously reported and especially at the low temperatures of protoplanetary disks (e.g. Gundlach et al. 2018), whereas dry grains (bare silicate or refractory carbonaceous) have an increase in sticking force at high temperatures (Kimura et al. 2015, e.g.), leading to a “sweet spot” for grain growth around 1200-1400 K (Bogdan et al. 2020; Pillich et al. 2021). These predictions are also consistent with a recent multiwavelength study of FU Ori (Liu et al. 2021), showing mm-sized grains in the hot inner disk and grain size ≤200 µm at larger separations.

Finally, another scenario involves the replacement of the original dust particles by a second-generation of dust: N-body simulations presented in Turrini et al. (2019) show that within an evolved disk – where dust has already evolved into larger bodies and planetesimals – the formation of giant planets in a disk can generate dynamical perturbations and create a highly collisional environment in its surroundings. Depending on the mass of the planets and on the original size distribution of the planetesimals, these collisions would generate a populations of rejuvenated dust (with sizes from 100 µm to centimeters) that could account for a large fraction of the dust that is measured in evolved disk (50–100% of the dust mass in the case of HD 163296 as estimated in Turrini et al. 2019). The same scenario is proposed to explain the trend of dust mass as function of age observed in nearby star forming regions in Testi et al. (2022): after an initial decrease with ~1/t the solid mass increases again at 2–3 Myr, which could be a sign of early planet formation and production of reprocessed dust. Recently, Doi & Kataoka (2021) claimed that the ring at 68 au exhibits an increased dust-to-gas scale height h_d/h_g with respect to the inner disk and to the outer ring at 100 au. A large scale height of smaller (micron-sized) dust in the two outer rings (corresponding to aspect ratios h/r ~ 0.25-0.3) was suggested as well in Guidi et al. (2018) to interpret the scattered light emission from HD 163296 in the thermal infrared. This hints to the presence of dynamically excited dust at this location, that we recall is predicted by planet-disk interaction simulations (Turrini et al. 2019; Bi et al. 2021; Binkert et al. 2021). This could also reconcile the dust temperature at 66 au derived in this work being closer to the gas temperature measured by Dullemond et al. (2020) compared to the 100 au ring (see Fig. 13). We note that the DSHARP observations of HD 163296 at 1.3 mm were carefully analyzed to search for localized emission from circumplanetary material in the main gaps at 48 and 86au in Andrews et al. (2021), who did not find any detection for this or other DSHARP disks included in the study. One of the possible causes could be the aforementioned higher scale height of the dust rings in inclined disk such as HD 163296, that would increase the extinction of the circumplanetary disk emission along the line of sight (e.g. see Fig. 4 in Guidi et al. (2018) in relation to the small grains).

Fig. 12 Temperature, maximum grain size and surface density as a function of the radius from the SED fitting with a model without scattering. The dashed red curve represents the estimates from the Monte Carlo fit including the statistical error only. The color map shows the normalized posterior distribution obtained merging 30 additional fits after introducing a random offset in the fluxes according to their flux calibration accuracy. The resonance in the opacity as function of amax (see Fig. 4) can result in a degeneracy of this parameter, e.g. in the inner disk inside ~40 au. In the top and bottom panels we overplot the temperature and surface density from previous studies.

Fig. 13 Same as Fig. 12, but for the model that includes scattering.

Fig. 14 Optical depth at the five different wavelength, resulting from the best fit model shown as a dashed curve in Fig. 13. The uncertainties (shaded regions) correspond to 1 and 2σ (16th/84th and 2.5th/97.5th percentiles, respectively) of the distribution at each radius (see main text).

Fig. 15 Albedo at the different wavelengths from our best-fit physical model, shown at specific locations in the disk. The error bars correspond to the 16th and 84th percentile of the set models, the vertical dashed line are drawn at the position of the flux peaks at 2 mm.

5.2 Dust mass

We calculate the dust mass from our derived surface density profile as $M_{\text{d}}={\displaystyle {\int }_{{\mathrm{r}}_{\mathrm{i}}}^{{\mathrm{r}}_{\mathrm{f}}}{\displaystyle \sum \left(r\right)2\pi r\text{d}r}}$. With r_i = 8 au and r_f = 120 au this is equal to 8.0 × 10⁻⁴ M_⊙ or 265 M_⊕ for the scattering model, and 1 × 10⁻⁴ M_⊙ 337 M_⊕ for the nonscattering model. We note that since we are pushing the spatial resolution of our observations to characterize in detail the inner structure of the HD 163296 disk, we have typically a low signal-to-noise in the outer disk, that was not included in our SED fitting. Therefore we are not sensitive to the dust mass contained at radii > 120 au, that shoul

anyway represent only a few percent of the total dust mass (see next paragraph).

We can compare our result with previous studies considering the same separation range: integrating the surface density from Isella et al. (2016) between 8 and 120 au, we find a dust mass of 0.47 × 10⁻³ M_⊙, about 1.7 times smaller than what found in this work. We note that the dust surface density in Isella et al. (2016) was derived using RADMC3D (Dullemond et al. 2012) without including scattering opacity on a single wavelength ALMA dataset at 1 mm. The difference can therefore be due to the fact that we fit for the opacity at each radius and we include dust self-scattering, although with a simplified analytical description. We note that integrating the surface density function found in Isella et al. (2016) up to radii of 400 au, we find that the remaining mass is only a small fraction of the total, specifically M_dust [8–120] au ≃ 97% M_dust [0.5–400] au.

We can estimate the mass of the spatially resolved rings at 66 and 100 au, where we define the radial limits for each ring as the closest minima to the ring peaks, calculated with scipy.signal.argrelmin on the surface density profile obtained in Sect. 4.3. The values are displayed in Table 4: we note that the mass of the outer ring at 100 au is larger than the one at 66 au and accounts for about a third of the total mass. This is opposite to what found in previous studies, for example, Dullemond et al. (2018) found 56.0 M_⊕ for the 66 au ring (consistent with our result), but only 43.6 M_⊕ for the 100 au ring, using only the DSHARP dataset at 1.3 mm. The discrepancy at 100 au remains if we compute the ring masses using the same radius intervals as in the cited paper, as we find 57 M_⊕ in the limits 52–82 au and 84 M_⊕ between 94-104 au.

In the context of disk demographic studies, the dust mass is typically calculated using a simple scaling of the flux density with a fixed dust opacity, as originally described by Hildebrand (1983): $$M_{\text{d}}=\frac{F_v{d}^2}{{\kappa }_v{B}_v\left(T\right)}.$$(4)

If we use this relation to derive the dust mass in HD 163296 from the integrated fluxes from the ALMA observations (Table 1), using T = 21K as the median disk temperature for the Black Body term, and an opacity of κ_v = 10 cm² g⁻¹ (v/1000GHz)^β with β = 1, we obtain dust masses on the order of 60–80% of the one derived in our multiwavelength analysis. Specifically, we obtain M_d,0.9mm = 63% M_d,mwle, M_d,1.3mm = 73% M_d,mwle, M_d,2.1mm = 64% M_d,mwle and M_d,3.0mm = 75% M_d,mwle.

We must stress here that the value for the dust mass we obtain from our spectral fit are heavily dependent from the choice of the dust composition. Within the different models we explored, we showed that a difference of a factor of 5 in surface density – that corresponds to the same factor in terms of dust mass – is already present varying only the porosity of the grains. Our dust mass estimates from the porous grains models are listed in Appendix D and indicate a total dust mass of 4.26 × 10⁻³ M_⊙ or 1417 M_⊕, i.e. about 5 times larger than our reference model with standard grains. Changing the dust constituents would also cause a difference: for example, the amorphous carbons from the DIANA project that we employ in this study have absorption opacities that are a factor of a few larger than the standard opacities used in the DSHARP modeling (Birnstiel et al. 2018). We showed in Appendix E how for the case of HD 163286 this translates in surface densities (i.e. dust mass) larger by a factor of ~3.

In relation to demographic studies that found a linear size-luminosity relationship for protoplanetary disks at the frequency of 340 GHz (Tripathi et al. 2017; Andrews et al. 2018b), we recall here that one of the scenarios proposed to explain such scaling relation invoked the presence of sub-structures in the dust that would generate optically thick emission with filling factors of a few tens percent (see Ricci et al. 2012). In the case of HD 163296 this is verified, as we measure a filling factor flux_thick/flux_tot of 82% at Band 7 in the range 8–115 au, with the upper limit in separation taken as the effective radius R_eff given in Tripathi et al. (2017).

5.3 Origin of the ringed structure

Revealing the grain size and distribution across a structured-disk can provide useful information on the mechanisms that generate such substructures. Theoretical studies show that various forms of instabilities, including the presence of embedded planets, generate pressure traps in the gas distribution of different shape and intensity; these, in turn, cause different levels of segregations of dust grains depending on their size (e.g. Pinilla et al. 2012; Nazari et al. 2019). So far, several studies have proposed the presence of planets in the disk of HD 163296, based on the observed depletion in the dust and gas surface density (Isella et al. 2016; Liu et al. 2018) and signatures in the gas kinematics (Pinte et al. 2018; Teague et al. 2018).

Sierra et al. (2021) found evidence of dust trapping in a pressure bump at the 100 au ring (hint by the higher grain size at this ring). We do not find such an indication in this study, but we showed how this result is highly dependent on the assumptions in the dust composition (see Appendix E). Even if we do not find evidence of such a differential trapping of dust grains in the rings, we note that we see an increase in the surface density at the 100 au ring, with respect to the inner 66 au ring. This is generally predicted by disk-planet interaction simulations, where the presence of a planet results in an accumulation of dust in the ring external to the planet, and a starving of material in the internal ring (e.g. Zhu et al. 2012; Rosotti et al. 2016; Binkert et al. 2021).

Finally, as discussed in Sect. 5.2 the relatively small grain size in the outer rings could indicate the effect of collisions excited by an embedded planet. In this case, we note that the size distribution of the particles could be significantly different than the power-law assumed in this study, and this can in turn affect the derived maximum size of the distribution.

Fig. 16 Model comparison. Top panel: Bayes factors computed between the evidence of the fits with the physical model without scattering (Eq. (2)) and the parametric model (Eq. (1)). Bottom panel: Bayes factor between the scattering (Eq. (3)) and nonscattering model.

5.4 Model comparison

In this study, we fit observations of HD 163296 at five wavelengths using a parametric model and a physical model (with and without scattering). To determine the relative strength of each model in predicting our datasets, we use again the Bayes factor K (see Sect. 3.4). In the top panel of Fig. 16 we confirm the expected outcome that a physical model is a better description of our submillimeter data with respect to a simple power-law. In particular, there is decisive evidence (K > 100) in favor of the physical model at 62% of the radii, a strong evidence (K > 10) at 64% or the radii and a moderate evidence (K > 3) at 72% of the radii. On the opposite, strong evidence goes down to only 4% and moderate evidence at 6% of the radii for the power-law model. This indicates that a simple prescription of the optical depth as a power-law function of the frequency is less-likely to predict the observed emitted intensity (see also Fig. D.3 in the Appendix), compared to the analytic physical model described in Eq. (2).

When comparing the physical models with and without scattering, we do not find a strong change in evidence in favor of the former (Fig. 16, bottom panel). We also find that the relative evidence between scattering and nonscattering model varies significantly among the single size distribution models. We can only conclude that a nonscattering model can equally reproduce the observations, generally predicting a lower temperature with respect of the scattering model.

6 Conclusions

We analyzed new ALMA and VLA observations to study the dust properties in the rings of HD 163296. We employ parametric and physical descriptions of the flux density to reproduce the Spectral energy distribution as function of the radius, and compare the performances of the different models in a Bayesian framework. We summarize here our main results:

• We estimate a non-dust contribution from the inner disk (r ≲ 5 au) that accounts for about 5% of the total flux at 9 mm and 40% of the total flux at 3 cm. Based on its mm-cm spectrum, its origin is consistent with free-free emission from a disk wind.

• We extract the radial brightness profiles at all wavelengths with the python tool frank fitting the original visibilities. We recover a total of six bright peaks in the flux distribution of HD 163296 within a separation of 120 au, i.e. one more ring than what derived from the convolved ALMA images at the highest available resolution (e.g. Huang et al. 2018b).

• We fit a simple opacity power-law prescription to the extracted radial profiles at wavelengths from 0.88 mm to 9 mm and independently at each radius. We find that the dust temperature increases in the gaps at average separations of 10, 50 and 83 au, and the opacity spectral index β is consistent with the ISM values (micron-sized grains) at most radii, while it decreases to values below 1.7 in the inner disk (r ≲ 30 au)

• We follow the same procedure using a physical model with an analytical expression for radiative transfer with and without dust scattering. We consider a grid of eight different dust populations, where we vary the size distribution slope and the composition. The best-fit models indicate the presence of 200 µm sized grains in the outer disk at r ≳ 40 au. The grain size is less constrained in the inner disk, but there are local indications of larger grains (≳millimeter) at separation smaller than ~30 au.

• We find that our observations are generally better described by compact grains (porosity of 25%) than by porous grains (80% porosity). We also observe that a steeper size distribution (with q = 3.5/4) better describes the outer rings at 66 and 100 au.

• By comparing the evidence of the parametric and physical model, we confirm that the latter is better at predicting our datasets with respect to the simple power-law.

• Finally, we note that different choices of dust opacities change the estimates of the dust mass of a factor of a few (a factor of 5 just varying the dust porosity), and locally affect the derived grain size in a less systematic but significant manner.

Our results confirm that the presence of optically thick structures can artificially lower the millimeter spectral index, and therefore in general effects of dust self-scattering should not be neglected when interpreting the continuum millimeter emission in terms of dust properties. This factors, combined with the highly non homogeneous spatial distribution of dust grains, accentuates the necessity of both high spatial resolution and spectral information, as they are crucial for resolving the small-scale structures and removing the degeneracies between parameters such as dust density and grain size. It is not always possible to fully break such degeneracy, despite the large spectral information available (wavelengths from 0.9 to 9 mm), especially where the substructures are not resolved.

We stress that a major caveat is still represented by the unknown composition of the dust grains that can significantly affect the final dust parameters and make the comparison with theoretical prediction more challenging.

Regarding the origins of the observed dust rings in the HD 163296 disk, the most accredited scenario invokes the presence of three to four giant planets perturbing the dust and gas structure. Our results are in general agreement with this scenario: the higher values of the dust temperature we derive in the optically thin gaps at 10, 50, 85 and 115 au is expected if a planet is present in the gap as resulting from hydrodynamical calculations (Isella & Turner 2018). A higher dust scale height in the 66 au ring reported independently by Doi & Kataoka (2021) would explain the temperature values we derive at this location being closer to the gas temperature, and would be consistent with vertical stirring of the dust caused by planets (e.g. Binkert et al. 2021).

We do not detect a significant size difference between the grains in the two outer rings and the adjacent gaps, as one could expect in presence of planets, as the pressure trap generated by the planets would result in larger grains being trapped at the pressure maxima in the gap edges (e.g. Pinilla et al. 2012; Rosotti et al. 2016; Nazari et al. 2019). However, the small 200 µm grains found in the rings at 66 and 100 au could belong to a second-generation dust population resulting from collisions of large km-sized bodies, that could account for a large fraction of the dust mass in the rings (Turrini et al. 2019). Furthermore, we derive a higher surface density in the 100 au ring with respect to the 66 au ring, which is interestingly consistent with some theoretical predictions of disk with embedded planets and could be related to the protoplanet predicted at a radius of 80-90 au (e.g. Isella et al. 2016; Teague et al. 2018; Izquierdo et al. 2022). In this scenario the dust grains responsible for the (sub)millimeter emission could have a more complex three-dimensional structure rather than being distributed in a thin midplane. A future 3D Monte Carlo radiative transfer modeling (including full scattering) of multiwavelength observations will help in reconstructing the full structure of the solids in this disk.

Acknowledgements

This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation. G.G. acknowledges the financial support of the SNSF. A.I. acknowledges support from the National Science Foundation under grant No. AST-1715719 and from NASA under grant No. 80NSSC18K0828. H.B.L. is supported by the Ministry of Science and Technology (MoST) of Taiwan (Grant Nos. 108-2112-M-001-002-MY3). G.R. acknowledges support from the Netherlands Organisation for Scientific Research (NWO, program number 016.Veni.192.233) and from an STFC Ernest Rutherford Fellowship (grant number ST/T003855/1). I.dG.M. is partially supported by MCIU-AEI (Spain) grant AYA2017-84390-C2-R (co-funded by FEDER). H.L. gratefully acknowledges the support by the LANL/CSES program and the NASA/ATP program. M.T. has been supported by the UK Science and Technology research Council (STFC) via the consolidated grant ST/S000623/1, and by the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 823823 (RISE DUSTBUSTERS project). This paper makes use of the following ALMA data: ADS/JAO.ALMA #2015.1.00725.S, ADS/JAO.ALMA #2016.1.01086.S, ADS/JAO.ALMA #2016.1.00484.L, ADS/JAO.ALMA #2017.1.01682.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Appendix A ALMA flux calibration errors

To get an estimate of the calibration error that affects our ALMA datasets, we compared the amplitudes of the flux calibrators that were set in the ALMA pipeline or calibration script with the measurements of the same calibrators taken before and after the science observations from the ALMA Calibrator Source Catalogue³. We note that, as described in Sect. 2, when combining multiple epochs and/or antenna configurations we scale the amplitudes to match the execution block that has been calibrated using the most recent calibrator measurements. Therefore for each ALMA band we plot only the fluxes from the epoch we used as reference, instead of plotting all the single scheduling blocks.

At Band 3 we find that frequent measurements of the flux calibrator J1733-1304 were taken around the date of our HD 163296 observations, as we show in Figure A.1, upper panel. The quasar shows a flux increasing with time within the selected interval (1.5 months before and after the science observations), so we can get an estimate of its flux in function of time by performing a linear fit on the data points. The error bars on the calibrator fluxes in the plotted time range are on the order of 1–3%, and the deviation of the flux used during ALMA calibration (diamond markers in Figure A.1) lie within this range (1% for the lower side band and 3% for the upper side band).

At Band 4 the amplitude calibration relied on measurements of the flux calibrator J1924-2914 at other frequencies: therefore we plot the observed data points at the two sidebands of Band 3 were close to the Band 4 observation date. We scale the flux used during calibration to the Band 3 frequencies using the spectral index specified in the corresponding CASA calibration script in the setjy task. We note that the measured fluxes present a rapid temporal variability within the selected time frame, so that a linear fitting is no longer adequate to interpolate the data points. If we estimate the calibrator fluxes at the time of our science observations by a simple linear regression between the precedent and following measurements (Figure A.1, middle panel), we find a deviation on the order of 1% and 2.6% in the upper and lower side-band, respectively.

Band 6 calibrated data were taken from the DSHARP program⁴: amplitude scaling to match the most recent catalog entries was performed (as described in Andrews et al. 2018a) and the final dataset results from the combination of several datasets from multiple epochs and configurations, so that we can reasonably assume that the calibration error was partially attenuated from the nominal 5% (Bonato et al. 2018).

Finally, we find several measurements of the J1733-1304 calibrator at Band 7, so we can compare directly the flux used during calibration with the measurements before and after (Figure A.1, lower panel). We find a deviation <1% between the flux used during calibration and the estimated flux at that date, with the closest measurements taken only one day before the science observations and reported with an error of 3.7%.

With this exercise we do not aim to estimate accurately the calibration error at the different bands, but only to verify that there are no major issues affecting the absolute flux calibration, mostly out-of-date reference values for the calibrators. We note that using an updated catalog is fundamental to attenuate the calibration offset uncertainty when employing a quasar as a flux calibrator, because of their intrinsic temporal variability. We verified that this is the case for our datasets, so that the flux calibration accuracy for the single-epoch observations can be taken

Fig. A.1 Calibrators measurements. Upper panel: ALMA measurements of the J1733-1304 calibrator are plotted as blue and red circles, with corresponding error bars, at the two frequencies of Band 3, respectively. The diamond markers represent the amplitudes set for the flux calibrator in the CASA calibration script. Dashed lines are the linear regression between the plotted catalogue measurements Middle panel: same as for the upper panel, but with the values of the flux calibrator for Band 4 (J 1924-2914), scaled at the frequencies of Band 3 with the spectral index used during calibration. Lower panel: measurements of the J1733-1304 calibrator used for the observations at Band 7. The time range on the x-axis is three months in all the panels.

indicatively as the nominal value, that is ~5% for Bands 3 to 6, and ~10% for Band 7.

Appendix B Free-free emission estimate

Fig. B.1 Visibilities in the VLA X Band: the green solid line represents the best fit for a gassian + a constant in the deprojected real visibilities.

We inspected our ALMA and VLA datasets to detect and estimate a possible contamination from free-free electrons or other nonthermal processes. Despite the contribution below 7 mm should be negligible (Natta et al. 2004; Guidi et al. 2016), we included in the analysis our ALMA dataset at the longest wavelengh (≥3 mm), taking advantage of the moderately high resolution to check for possible point-like sources in the center.

If we assume that such emission comes from a very compact region close to the star (a Delta Dirac-like function in the image plane), we can describe it as a constant in the Fourier plane. i.e. we can fit an horizontal line in the deprojected visibility profile. We start analyzing the deprojected visibilities at Band 3: to choose the lower end of the spatial frequency range to include in the fit, we first fit a line of the form $y{\left(\rho \right)}_{\rho >{\rho }_0}=a·{\rho }_{\rho >{\rho }_0}+b$ where ρ is the uv-distance defined as $\sqrt{u^2+v^2}$. By varying the starting point ρ₀, we observe that at 3000 kλ the slope starts to approach zero, i.e. the visibilities reach an asymptotic value. This latter is similar for both Band 3 sidebands and corresponds to 0.52 ± 0.20, where the error is calculated as the sum in quadrature of the covariance of the parameter and the systematic calibration error.

We used the same method and starting point (3000 kλ) for the VLA data at 8.6 and 9.7 mm, and we find values of 0.17 ± 0.03 and 0.16 ± 0.20, respectively.

The VLA data at 3 cm cannot reach the high resolution of the observations at shorter wavelengths, and are sampled up to ~1000 kλ only. As the visibility profile appears as a combination of a compact source plus a gaussian-like distribution (Figure B.1), we fit a function composed by gaussian plus a constant over the whole range of visibilities. That would give us an estimate of the dust emission (gaussian) and a compact central contribution (constant). This results in a compact source of intensity (0.19 ± 0.02) for the low frequency side of the X band (33 mm), and (0.20 ± 0.02) for the high frequencies (27 mm).

As an additional step, we inspect the images to try and obtain further information on the spatial extent of a central contaminating emission. We produce images with uniform weighting (robust = −2 and nterms = 2 in the tclean task in CASA), to achieve the maximum resolution from aperture synthesis. The model images in Figure B.2 show the clean components resulting at the end of the cleaning iterative process, for the three wavelengths we considered: 3.2 mm, 9.1 mm and 30 mm. In the VLA observations (central and right panels in Fig. B.2, a bright central source is clearly dominating the emission, while this is not the case for ALMA Band 3 (left panel). In particular at 10 GHz (X band), the free-free is expected to dominate the total flux (see. e.g. Pascucci et al. 2012), and since the observations are marginally resolved we can obtain an estimate of the extent of this central bright emission. By deconvolving the image with the task imfit in CASA, we obtain a deconvolved size with FWHM ≃ 0.″1, corresponding to 10 au. Therefore we can conclude that the free-free is dominating the 3 cm flux within a radius of about 5 au. Similarly, the deconvolved size at 33 GHz (Ka band) results ~0.″09 or 9 au. The flux in the deconvolved model images inside this area (within a 5 au radius) can be translated in an upper limit for the free-free (as it will contain both dust emission and free-free). This results 2.97 ± 0.30 mJy at 3 mm, 0.55 ± 0.06 mJy at 9 mm and 0.21 ± 0.01 mJy at 3 cm, where the uncertainties are dominated by the calibration error. In the case of 3 mm, we expect a putative free-free emission to be very small compared to the total flux, and the dust is likely to dominate the emission in the central regions.

Finally, we subtract the fourier components corresponding to this central model emission from the original VLA visibilities, to obtain “corrected” visibilities that are used to extract the dust brightness profiles (see Sect. 3.2). Given the uncertainties in the free-free estimation, the inner regions will be in any case masked for our multiwavelength analysis, but the removal of the central point-source is necessary for a good convergence of the forward modeling with the frank python package.

Fig. B.2 Model images obtained with uniform cleaning for the three long wavelength datasets. The field of view s 0.″2 for all images, while the pixel size depends on the dataset's resolution and is chosen according to the Nyquist sampling theorem, in order to have 2-3 pixel across one beam.

Appendix C Brightness profiles extraction

We use the python package frank (Jennings et al. 2020) to derive the brightness profiles of HD 163296 from the ALMA and VLA observations, in the approximation that the emission is azimutally symmetric. This method takes advantage of the radial-only dependence of the intensity I(r), which allows one to simplify the relation between this latter and the visibility function V(ρ) (where ρ is the uv-distance or spatial frequency) with a Hankel transform instead of a Fourier transform (Pearson 1999). This consistently speeds up the inversion procedure from the visibility to the real plane, without the need of a discrete Fourier transform, and each fit is normally completed within minutes.

The frank method follows a statistical approach and taking as an input the visibility data points, it retrieves a final brightness profiles with an associated confidence interval, that represent the statistical error resulting from the initial uncertainties on the visibility data. These latters are provided in the input tables as their square inverse value or weight = 1/σ², as extracted from the ALMA MS tables: such weights are associated with the RMS noise of each visibility and are first initialized with a value of 2∆v∆t (for the most recent CASA versions), and subsequently modified within the CASA software during the calibration process, when they are multiplied by the antenna-based gain factors calculated at each calibration step, and scaled for the system Temperature. Additional parameters entering the fit are related to the geometry of the system: for consistency we fix the inclination and position angle of the disk to the same values for each wavelength (46° and 133°, respectively), while we use the internal routine in the frank package to determine the offset in Right Ascension and Declination of each dataset, providing an initial guess obtained by comparing the positions of the phase center and the brightness peak in the synthesized images (see Figure 1) derived with a gaussian fit to the central regions of the maps with the CASA tool imfit. Finally, the outer radius is set to 2.2 arcseconds for the ALMA datasets and lower values for the VLA datasets (175 for 9 mm and 172 for 3 cm), and the number of radial bins is set to 300. We note that these last two parameters have very little influence on the final fit, as described in Jennings et al. (2020). For each dataset, we generate a grid of

models varying the hyperparameters α and w_smooth, which corresponds to a change in the prior in this statistical framework. Using values of α = [1.05, 1.1, 1.2, 1.3] and w_smooth = [0.0001, 0.001, 0.01, 0.1] we obtained 16 models for each dataset. The best-fits are chosen as the best compromise between reaching high spatial resolution (i.e. fitting to the long uv-distances) and limit the inclusion of noisy data points in fit, that generally correspond to the longest baselines and produce artificial oscillations in the extracted brightness profiles.

A particular care is needed with the VLA datasets: as described in Sect. 3.1 these observations present a strong compact emission in the center, that in the Fourier space translates into a nonzero asymptotic value at the highest spatial frequencies.

This impacts the frank solutions introducing strong oscillations in the intensity profiles, at scales corresponding to the uv-distance where the power spectrum (or the frank prior) drops to zero. The reason of this behavior is explained in (Jennings et al. 2021), and is related to the frank power spectrum that needs to converge to zero by construction. For this reason, we correct the input visibilities for this compact emission, which we already discussed it is not associated with dust (see Sect. 3.1). We therefore subtract from the VLA 9 mm dataset the constant value previously

Fig. C.1 Observations and residual images for all wavelengths fit with frank, displayed with the same color scale.

obtained by fitting the real visibilities at the high uv-distances (see Appendix B). This approach is suggested by the authors of the package in Jennings et al. (2021), and seems to be working fairly well for this dataset, allowing us to unveil the 9 mm emission in the 100 au ring, already hinted by the images produces with CASA tclean, but too noisy to be appreciated in the azimuthally averaged profiles of the cleaned images. We note that the central emission (≲10 au) remains affected by the uncertainties on the free-free estimate and subtraction, so that we do not consider it a reliable representation of the dust emission. At 3 cm, even after such a point-source correction, we still get significant oscillations and negative values in the final solutions. We find the best result by first subtracting from the visibilities the central emission obtained by the clean model with robust = −2 described in Appendix B. This is done by transposing the clean model of the central source into the Fourier place with the task ft in CASA, and then subtracting it from the data. We note that this means subtracting the dust central emission as well, but since the uncertainty on the free-free contribution is anyway limiting our capability to give a precise estimate of the dust central flux, this does not represent an issue, as our goal is trying to unveil the emission at larger radii.

The final models are displayed in Figure 3, the corresponding values of the hyperparameters are listed in Table C.1, along with a metric proposed in Jennings et al. (2021) to characterize the spatial resolution accuracy of the fits: we measure B₈₀ as the shortest baseline beyond which the difference between the model and the observations has values ≥20% for at least 200kλ consecutively, and it is meant to represent the baseline where the fits starts to depart appreciably from the data. We note that we choose to fit the whole frequency band even at long wavelengths (larger than 3 mm), as the fit obtained separating the datasets in higher and lower frequencies (as for the synthesized images in Figure 1) results in a much lower quality solutions (higher oscillations at large radii and larger RMS error). In Figure C.1 we show the residual maps from the modeling: these are produced from the frank residual visibilities and using the same parameters in CASA tclean as for the observations images (displayed next to the residual maps).

As explained in Sect. 3.2, for the purpose of the spectral modeling a second round of fits is performed so that all the datasets are sampled at the same spatial scales. This way we obtained five spectral profiles with same accuracy B₈₀ = (2100 + 100) kλ from 0.9 mm up to 9 mm. Such a resolution could not be achieved for he 30 mm dataset, that is sampled up to lower spatial scales (as reported in Table C.1), and was therefore not included in the SED modeling. However, we show the prediction of our models at this wavelength in Figure D.5.

Appendix D SED fitting

In Sect. 4.2 we show the results of fitting a power-law opacity (eq. 1) to the brightness profiles from 0.8 mm to 9 mm. In all the fits the uniform prior probabilities are set for all the three

Fig. D.1 Results of the Monte Carlo fits for the nonscattering models. Left panel: best-fit parameters for the standard composition and the four different size distributions for the nonscattering model. The shaded regions represent the 16th and 84th percentile of the posterior distributions at each radius, from the fits perfomed including only the statistical error in the likelihood evaluation. Right panel: same as the left panel but for porous grains.

Fig. D.2 Same as in Figure D.1 for the scattering model.

parameters, in the intervals log₁₀τ = [−4, 3], β = [0, 5], and T (K) = [3, T_up]. Similarly, the best-fit parameters of the physical model T, a_max and Σ_d (Sect. 4.3) are obtained using a uniform prior probabilities for T (K) = [3, T_up], log₁₀a_max (cm) = [−4, 3] and Σ_d (g/cm²) = [0, 10]. The temperature prior upper limit varies with the radius, and corresponds to the temperature of the surface layer in the optically thin approximation, computed as in Dullemond et al. (2001), Equation 7. For the effective temperature and stellar radius we use the values computed by Setterholm et al. (2018) (R_* = 1.6 R_⊙ and T_*,eff = 9250 K). The Planck mean opacities for protoplanetary disks are taken from Semenov et al. (2003) and computed with the fortran script opacity. f made available by MPIA Heidelberg⁵. These variable upper limits are useful especially where the Temperature is poorly constrained such as in the dust gaps, as they allow us to obtain physically meaningful values for the temperature that could otherwise undertake indefinitely high values. At each radius, we include the fluxes above 3 times the rms value at the single wavelengths, calculated as the standard deviation in the signal-free region of the disk. If this criteria is satisfied for less than 3 wavelengths we do not perform the fit at that location.

At radii smaller than 8 au, some of the fitted parameters result highly unconstrained: in particular the posterior distributions of the surface density (in the physical model) and the β index (in the parametric model) tend to accumulate on the upper and lower edge of the prior limits, respectively. This can be expected in the case of a highly optically thick emission, when the surface density or the opacity index are very loosely constrained because of the dominance of the Planck term in equation 1 and 2. In fact, even enhancing the upper limit of the prior to very high values (10⁵ g/cm²), we still observe the same behavior of the posterior distribution. Another explanation could be the failure of the simplified radiative transfer formulation in describing the emission in the very inner regions. In Figure D.1 and D.2 we show the

Fig. D.3 Absorption opacity at three different location (inner disk at 20 au, gap at 80 au and ring at 100 au) from the nonscattering model plotted as circles. The dashed lines represent the power-law from the parametric fit, with κabs = κ0(v/v0)β where we used κ0 = κ1.3mm from the physical model.

best fits for all 8 dust populations introduced in Sect. 3.3 for the nonscattering and scattering model, respectively.

In Figure D.3 we show the absorption opacities derived from the physical model compared to a power-law prediction. Given the lower evidence of the parametric model shown in Sect. 5.4, this indicates that the dust opacities are not well approximated with a power-law.

In Sect. 4.3 we show how the standard grains result to be a better model for describing our observations, and therefore we choose these latter as our reference model. The comparison between the predicted and observed brightness profiles is displayed in Figure D.4, while we show in Figure D.5 the predicted flux at 3 cm. The scattering albedo at each wavelength for the scattering model are shown in Figure D.6.

Finally, we show in Table D.1 the masses we derived with the porous-grains model, as a comparison with the reference values displayed in Table 4. While the total mass is higher by a factor of four, the mass in the 100 au ring results 40% of the total, similarly to what derived for the standard grains in Sect. 5.2.

Fig. D.4 Flux densities at the five different wavelengths predicted by the best-fit physical models using standard grains in the nonscattering and scattering case. The red and black dashed-dotted curve shows the models predictions at the radii where the flux at each wavelength was included in the fitting procedure (only fluxes larger than 3 times the rms error were considered at each radius). The grey curve shows the scattering model prediction at all radii.

Fig. D.5 Flux densities at 30 mm predicted by the best-fit models in the nonscattering and scattering case, shown as red and black dashed-dotted lines, respectively. The observed brightness profiles at VLA Band X are plotted with the corresponding 1, 2 and 3 σ errors (computed as the sum in quadrature of the statistical error and the flux calibration error) as shaded regions.

Fig. D.6 Scattering albedo as function of the radius for the five considered wavelengths. The shaded areas are drawn between the 16th and 84th percentile of the model distribution, as described in Sect. 3.3.

Appendix E Comparison with Sierra et al. (2021)

In a recent study Sierra et al. (2021) found the presence of millimeter-grains in the disk of HD 163296 outside 40 au and a local increase in grain size at ~ 100 au. Since in this work we do not find millimeter grains or signature of dust trapping in the outer rings, we run a series of test to determine what are the main factor responsible for this discrepancy. First, we note that there are some difference between these two studies: in terms of datasets, this work relies on a larger wavelength coverage (4 ALMA and one VLA band, compared to 2 ALMA bands (1 and 3 mm) in Sierra et al. (2021)), and a higher spatial resolution of a factor of about 2.5. They also assume a midplane temperature profile from Zhang et al. (2021), obtained by modeling the CO emission lines and continuum emission at 1 mm. Finally, the assumptions on the dust composition are different: our absorption opacities are on average a factor 2–3 higher than the DSHARP opacities used in Sierra et al. (2021), while the scattering opacities are comparable (this is related to the presence of amorphous carbon from Zubko et al. (1996) in the DIANA opacities, that produce higher absorption coefficient partially suppressing the role of dust self-scattering).

To test whether the choice of dust opacity plays the major role in the final estimates of a_max, we performed a quick check by taking the flux profiles at 1, 2.3 and 2.6 mm published in Sierra et al. (2021) and performing the same procedure that the author described but using the standard DIANA opacities described in this work - instead of the DSHARP - and with a size distribution q = 2.5. We fix the temperature profile as in Sierra et al. (2021), and approximate the statistical error to 5% at each radius, summing in quadrature an additional 5% for the calibration error (and giving a double weight to the Band 3 data). Using a two-dimensional grid with 200 bins for both log a_max [−4, 3] and log Σ_d [−3,1], we calculate the likelihood as L = exp(−0.5 · χ²), with ${\chi }_r^2={\displaystyle {\sum }_{v_i}{w}_i{\left(F_{obs,r\text{,}{v}_i}-F_{\mathrm{mod},r\text{,}{v}_i}\right)}^2/{\sigma }_{\text{r,}{v}_i}^2}$ and the model fluxes computed with eq. 3. We derive the estimates for the two parameters as the values at the peak of the marginalized posterior distribution. We focus on the outer disk (r ≳ 40 au) and find a well defined posterior distribution for a_max and Σ_d at all radii, as described in Sierra et al. (2021). The maximum grain size found using DIANA opacities are systematically lower than the one found using DSHARP opacities, by a factor of 3 on average (see Figure E.1, top panel). The DIANA surface densities are lower by a similar factor (Figure E.1, bottom panel), which is expected because of the higher values of absorption opacities mentioned above. As a check, we plot the final values from Sierra et al. (2021) in the scattering case, and confirm that we get the same results (the small difference with our best-fit values in orange is likely due to our approximation of the flux statistical error). We note that the maximum grain size calculated with the DIANA opacities does not seem to have a local increase at the 100 au ring, opposite to what found with the DSHARP grains and in Sierra et al. (2021).

Estabilished that the grain composition and dust opacities can significantly affect the grain size and surface density values, we need to consider that signal dilution - where the disk structures are not resolved - and the assumption of a fixed dust temperature can play a additional role. To get an approximate idea of these effects, we run our nested sampling routine fitting the datasets used in this work, first using the DSHARP opacities, and then adding the fixed temperature profile as in Sierra et al. (2021). We note that this is just a simple test that relies on the data statistical error only (we explained how we additionally computed the effect of calibration error with 30 additional fits in Sect. 3.3), so that it does not give a clear indication of the absolute uncertainties, but it can be useful for a relative comparison. We show the results in Figure E.2: we find that the spatial resolution and the assumed temperature profile have both an impact on the final parameters. We note that fixing the temperature profile (that at 100 au corresponds to a ~50% increase compared to out bestfit values) can induce a significant variation in a_max resulting in a larger grain size in this ring.

Fig. E.1 Parameters estimates using two different dust opacities. Top panel: maximum grain size in the outer disk of HD 163296 computed from lower resolution data (see main text), using the DIANA and DSHARP opacities, with a grain size distribution with q = 2.5. We over-plot with empty black circles the values from Sierra et al. (2021). Bottom panel: same as for the top panel, but showing the surface density estimates.

Fig. E.2 Comparison of the dust parameters in the outer disk estimated with a Monte Carlo nested sampling tool in the scattering case, using statistical error only and for a size distribution with q = 2.5. Empty markers represent the fixed Temperature n the upper panel, and the values of the best-fot parameters from Sierra et al. (2021) in the middle and bottom panel.

References

All Tables

All Figures

	Fig. 1 Continuum maps of HD 163296 from ALMA observations (top row) and VLA (bottom row). The synthesized beams of the final maps are drawn as white ellipse at the bottom left corner of each panel and are listed in Table 1. The images are available in electronic form at the CDS as FITS files.
In the text

Fig. 2 Spectral energy distribution of HD 163296. The blue circles represent the integrated fluxes from the data analyzed in this work, while empty squares are values from the literature (Isella et al. 2007; Natta et al. 2004; Guidi et al. 2016). The blue dotted curve corresponds to a curve with spectral index of 2.6. Red stars are the contamination from a central compact emission, determined by the asymptotic values of the visibilities as described in Sect. 3.1, and that are best fitted by a power law with index 0.11. Green triangles are upper limits for the free-free emission (see Sect. 3.1), and correspond to the emission within a radius of 5 au.

In the text

	Fig. 3 Radial profiles of the continuum emission derived with frank for each dataset, the shaded regions represent the statistical error computed by the frank fit. The star symbol in the plots relative to the VLA datasets indicates that these have been first corrected for the free-free/centrally peaked contamination, as described in Appendix C.
In the text

	Fig. 4 Dust opacity at 1.3 mm as a function of the maximum grain size. Upper panel: absorption (solid lines) and scattering (dashed lines) opacities for the standard dust composition and different size distribution q. Lower panel: same as the upper panel but for grains with a porosity of 80%.
In the text

Fig. 5 Flux density spectral index. Left: map of the spectral index computed from the ALMA continuum images at 0.9 mm, 1.3 mm, 2.1 mm, 2.8 mm and 3.2mm with a matching beam of 0.″095 × 0.″065 and position angle of 61.24°. Right: radial profile of the spectral index in the northwest and southeast sides (blue and red curve, respectively), error bars are shown for each bin and calculated as described in the text. Overplotted with dashed lines is the corresponding flux density profiles in Band 4 (right y axis). The grey shaded region on the left denotes the angular resolution of the map as half of the average beam FWHM.

In the text

Fig. 6 Flux spectral index between 0.9 mm and 3.2mm computed from the radial profiles extracted with frank (Sect. 3.2, black curve), with shaded areas representing the errors derived from the linear least-squares regression (with the uncertainties for the single fluxes defined as the sum in quadrature of the statistical and calibration errors). The profile is masked at the locations where the flux S/N < 3 for at least one wavelength. Overplotted is the azimuthal average of the map in Fig. 5 (green curve), calculated with bins of 1 beam (~0.″08) and with errorbars computed as for the profiles in Fig. 5, right panel.

In the text

Fig. 7 Best fit parameters for the power law model: temperature, optical depth, and opacity spectral index β as function of the radius. The dashed red curve represents the estimates obtained as the mean values of the posterior distributions from the Monte Carlo nested sampling fit. The color map is the normalized probability after merging the 30 posteriors obtained introducing a random offset in the fluxes according to their flux calibration accuracy. The dashed-dotted curve in the upper panel shows the upper limit of the Temperature prior for each radius used in the fit (see Appendix D). The vertical dashed lines correspond to the dust ring found in Sect. 3.2, while the white dashed line in the bottom panel is drawn for β = 1.7, corresponding to ISM dust grains.

In the text

Fig. 8 Dispersion of the best-fit temperature, maximum grain size and surface density in the non-scattering model, calculated as the mean value at each radius of the four different size distributions (the single best-fit models are show in Appendix D), considering standard and porous grains separately. The shaded regions correspond to the standard deviation of the four best-fit values at each radius.

In the text

	Fig. 9 Same as Fig. 8, but for the model that includes scattering (the single best-fit models are show in Appendix D).
In the text

Fig. 10 Bayes factor K calculated between the standard and the porous models as Zstandard/Zporous, for the nonscattering model (upper panel) and scattering model (lower panel). The different markers correspond to the different size distribution coefficients q. The points are color-coded according to their values: K > 1 (larger evidence for the standard composition) are drawn in cyan and K < 1 (larger evidence for porous composition) in magenta. Empty markers correspond to 1 < K < 3 or 1/3 < K < 1, full markers to K > 3 or K < 1/3. The horizontal dashed line is drawn at K = 10, so that all points above this line correspond to a strong evidence in favor of standard grains. The points larger than 100 are drawn at the location of 100, since for K ≥ 100 the interpretation in terms of evidence strength does not change (see Sect. 3.4).

In the text

Fig. 11 Bayes factor K computed between the model with the highest evidence (identified by the color of the marker) and the remaining 3 size distributions at each radius, for standard grains in the non-scattering (top panel) and scattering (bottom panel) case. The full color indicates that the Bayes factors K at that specific radius are all >10, i.e. the is a strong evidence for that size distribution compared to the other three discrete values.

In the text

Fig. 12 Temperature, maximum grain size and surface density as a function of the radius from the SED fitting with a model without scattering. The dashed red curve represents the estimates from the Monte Carlo fit including the statistical error only. The color map shows the normalized posterior distribution obtained merging 30 additional fits after introducing a random offset in the fluxes according to their flux calibration accuracy. The resonance in the opacity as function of amax (see Fig. 4) can result in a degeneracy of this parameter, e.g. in the inner disk inside ~40 au. In the top and bottom panels we overplot the temperature and surface density from previous studies.

In the text

	Fig. 13 Same as Fig. 12, but for the model that includes scattering.
In the text

	Fig. 14 Optical depth at the five different wavelength, resulting from the best fit model shown as a dashed curve in Fig. 13. The uncertainties (shaded regions) correspond to 1 and 2σ (16th/84th and 2.5th/97.5th percentiles, respectively) of the distribution at each radius (see main text).
In the text

	Fig. 15 Albedo at the different wavelengths from our best-fit physical model, shown at specific locations in the disk. The error bars correspond to the 16th and 84th percentile of the set models, the vertical dashed line are drawn at the position of the flux peaks at 2 mm.
In the text

	Fig. 16 Model comparison. Top panel: Bayes factors computed between the evidence of the fits with the physical model without scattering (Eq. (2)) and the parametric model (Eq. (1)). Bottom panel: Bayes factor between the scattering (Eq. (3)) and nonscattering model.
In the text

Fig. A.1 Calibrators measurements. Upper panel: ALMA measurements of the J1733-1304 calibrator are plotted as blue and red circles, with corresponding error bars, at the two frequencies of Band 3, respectively. The diamond markers represent the amplitudes set for the flux calibrator in the CASA calibration script. Dashed lines are the linear regression between the plotted catalogue measurements Middle panel: same as for the upper panel, but with the values of the flux calibrator for Band 4 (J 1924-2914), scaled at the frequencies of Band 3 with the spectral index used during calibration. Lower panel: measurements of the J1733-1304 calibrator used for the observations at Band 7. The time range on the x-axis is three months in all the panels.

In the text

	Fig. B.1 Visibilities in the VLA X Band: the green solid line represents the best fit for a gassian + a constant in the deprojected real visibilities.
In the text

	Fig. B.2 Model images obtained with uniform cleaning for the three long wavelength datasets. The field of view s 0.″2 for all images, while the pixel size depends on the dataset's resolution and is chosen according to the Nyquist sampling theorem, in order to have 2-3 pixel across one beam.
In the text

	Fig. C.1 Observations and residual images for all wavelengths fit with frank, displayed with the same color scale.
In the text

Fig. D.1 Results of the Monte Carlo fits for the nonscattering models. Left panel: best-fit parameters for the standard composition and the four different size distributions for the nonscattering model. The shaded regions represent the 16th and 84th percentile of the posterior distributions at each radius, from the fits perfomed including only the statistical error in the likelihood evaluation. Right panel: same as the left panel but for porous grains.

In the text

	Fig. D.2 Same as in Figure D.1 for the scattering model.
In the text

	Fig. D.3 Absorption opacity at three different location (inner disk at 20 au, gap at 80 au and ring at 100 au) from the nonscattering model plotted as circles. The dashed lines represent the power-law from the parametric fit, with κabs = κ0(v/v0)β where we used κ0 = κ1.3mm from the physical model.
In the text

Fig. D.4 Flux densities at the five different wavelengths predicted by the best-fit physical models using standard grains in the nonscattering and scattering case. The red and black dashed-dotted curve shows the models predictions at the radii where the flux at each wavelength was included in the fitting procedure (only fluxes larger than 3 times the rms error were considered at each radius). The grey curve shows the scattering model prediction at all radii.

In the text

	Fig. D.5 Flux densities at 30 mm predicted by the best-fit models in the nonscattering and scattering case, shown as red and black dashed-dotted lines, respectively. The observed brightness profiles at VLA Band X are plotted with the corresponding 1, 2 and 3 σ errors (computed as the sum in quadrature of the statistical error and the flux calibration error) as shaded regions.
In the text

	Fig. D.6 Scattering albedo as function of the radius for the five considered wavelengths. The shaded areas are drawn between the 16th and 84th percentile of the model distribution, as described in Sect. 3.3.
In the text

Fig. E.1 Parameters estimates using two different dust opacities. Top panel: maximum grain size in the outer disk of HD 163296 computed from lower resolution data (see main text), using the DIANA and DSHARP opacities, with a grain size distribution with q = 2.5. We over-plot with empty black circles the values from Sierra et al. (2021). Bottom panel: same as for the top panel, but showing the surface density estimates.

In the text

	Fig. E.2 Comparison of the dust parameters in the outer disk estimated with a Monte Carlo nested sampling tool in the scattering case, using statistical error only and for a size distribution with q = 2.5. Empty markers represent the fixed Temperature n the upper panel, and the values of the best-fot parameters from Sierra et al. (2021) in the middle and bottom panel.
In the text