# Limits on Millimeter Continuum Emission from Circumplanetary Material in the DSHARP Disks

Sean M. Andrews, William Elder, Shangjia Zhang, Jane Huang, Myriam Benisty, Nicolás T. Kurtovic, David J. Wilner, Zhaohuan Zhu, John M. Carpenter, Laura M. Pérez, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Sean M. Andrews, William Elder, Shangjia Zhang, Jane Huang, Myriam Benisty, et al. Limits on Millimeter Continuum Emission from Circumplanetary Material in the DSHARP Disks. The Astrophysical Journal, 2021, 916, 10.3847/1538-4357/ac00b9. insu-03636633

HAL Id: insu-03636633

https://insu.hal.science/insu-03636633v1

Submitted on 6 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Limits on Millimeter Continuum Emission from Circumplanetary Material in the DSHARP Disks

Sean M. Andrews¹, William Elder¹, Shangjia Zhang², Jane Huang³,<fnref n="11" /> and Myriam Benisty⁴,⁵
David J. Wilner¹, Zhaohuan Zhu², John M. Carpenter⁷, Laura M. Pérez⁸, Richard Teague¹, Andrea Isella⁹
and Luca Ricci¹⁰

¹Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

²Department of Physics and Astronomy, University of Nevada, Las Vegas, 4505 S. Maryland Parkway, Las Vegas, NV 89154, USA

³Department of Astronomy, University of Michigan, 323 West Hall, 1085 South University Avenue, Ann Arbor, MI 48109, USA

⁴Unidad Mixta Internacional Franco-Chilena de Astronomía, CNRS/INSU UMI 3386, Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile

⁵Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁶Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

⁷Joint ALMA Observatory, Avenida Alonso de Córdova 3107, Vitacura, Santiago, Chile

⁸Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile
⁹Department of Physics and Astronomy, Rice University, 6100 Main Street, Houston, TX 77005, USA

¹⁰Department of Physics and Astronomy, California State University Northridge, 18111 Nordhoff Street, Northridge, CA 91130, USA

*Received 2021 April 5; revised 2021 May 10; accepted 2021 May 11; published 2021 July 26*

## Abstract

We present a detailed analysis for a subset of the high-resolution (~35 mas, or 5 au) ALMA observations from the Disk Substructures at High Angular Resolution Project (DSHARP) to search for faint 1.3 mm continuum emission associated with dusty circumplanetary material located within the narrow annuli of depleted emission (gaps) in circumstellar disks. This search used the Jennings et al. frank modeling methodology to mitigate contamination from the local disk emission and then deployed a suite of injection-recovery experiments to statistically characterize point-like circumplanetary disks in residual images. While there are a few putative candidates in this sample, they have only marginal local signal-to-noise ratios and would require deeper measurements to confirm. Associating a 50% recovery fraction with an upper limit, we find that these data are sensitive to circumplanetary disks with flux densities $\gtrsim 50-70 \ \mu\text{Jy}$ in most cases. There are a few examples where those limits are inflated ($\gtrsim 110 \ \mu\text{Jy}$) owing to lingering nonaxisymmetric structures in their host circumstellar disks, most notably for a newly identified faint spiral in the HD 143006 disk. For standard assumptions, this analysis suggests that these data should be sensitive to circumplanetary disks with dust masses $\gtrsim 0.001-0.2 \ M_{\oplus}$. While those bounds are comparable to some theoretical expectations for young giant planets, we discuss how plausible system properties (e.g., relatively low host planet masses or the efficient radial drift of solids) could require much deeper observations to achieve robust detections.

*Unified Astronomy Thesaurus concepts: Planetary-disk interactions (2204); Protoplanetary disks (1300); Planet formation (1241); Exoplanet formation (492)*

*Supporting material: figure set*

## 1. Introduction

In just the past few years, the observational landscape of planet formation research has expanded dramatically. New measurements at very high spatial resolution (~few au) have revealed that protoplanetary disks are richly *substructured* (see Andrews 2020 for a review). Observations with the Atacama Large Millimeter/submillimeter Array (ALMA) demonstrate that the (sub)millimeter continuum emission from small particles in these disks is frequently distributed in bright rings, separated by darker gaps (Huang et al. 2018b; Long et al. 2018; van der Marel et al. 2019; Cieza et al. 2021). Hydrodynamics simulations show that interactions between planets and their birth environments can create these kinds of perturbations (Dong et al. 2015; Kanagawa et al. 2015; Bae et al. 2017). Given the locations, widths, and depths of these gaps, such simulations suggest that planets with masses $M_p \approx 0.1-1 \ M_{Jup}$ orbiting their (roughly solar-mass) host stars at distances $a_p \approx 10-100$ au may be common at ages of only ~1–3 Myr

(e.g., Jin et al. 2016; Clarke et al. 2018; Zhang et al. 2018; Lodato et al. 2019).

This work continues to be refined and expanded upon with more (sub)millimeter/centimeter continuum data, as well as resolved measurements of scattered starlight (Avenhaus et al. 2018; Garufi et al. 2018) and spectral line emission, where the latter reveals substructures in both the intensities (Isella et al. 2016; Guzmán et al. 2018; Law et al. 2021) and kinematics (see Disk Dynamics Collaboration et al. 2020). Collectively, these constraints on the initial architectures of planetary systems can be compared to the properties of the (older) exoplanet population to inform models of planetary migration and add context to analyses of direct imaging (e.g., Bowler 2016; Nielsen et al. 2019; Vigan et al. 2021) and microlensing surveys (e.g., Gaudi 2012; Penny et al. 2019).

In addition to these more general insights into the evolution of planetary systems, we are also seeing detailed case studies that quantitatively link disk substructures to planetary perturbers. The directly imaged planets in the cleared disk cavity around the young star PDS 70 are especially exciting examples (Keppler et al. 2018, 2019; Haffert et al. 2019). Using hydrodynamics

<footnote n="11">¹¹NASA Hubble Fellowship Program Sagan Fellow.</footnote>
simulations of the planetary accretion process as a guide (e.g., Ayliffe & Bate 2009; Tanigawa et al. 2012; Szulágyi et al. 2014), emission models predict that such planets may be easiest to discover in the mid-infrared (~5–20 µm), where the dust in a circumplanetary disk (CPD)<fnref n="12" /> should outshine the planetary photosphere (Eisner 2015; Zhu 2015; Szulágyi et al. 2019). Direct imaging of that emission may be common in the near future (i.e., with JWST and ~20–40 m ground-based telescopes) and will provide access to the thermal structures of those CPDs. A helpful complement is available at (sub) millimeter/centimeter wavelengths, where the optical depths are low enough to give some insight on the CPD (dust) masses (e.g., Szulágyi et al. 2018; Zhu et al. 2018; Wang et al. 2021).

Deep imaging with ALMA can reach continuum sensitivities comparable to expectations for CPD masses. While few dedicated (sub)millimeter/centimeter CPD searches have been attempted, those available provide stringent upper limits (Isella et al. 2014; Pérez et al. 2019; Pineda et al. 2019). Recent ALMA observations have identified millimeter continuum emission from the CPD associated with the PDS 70 c planet (Isella et al. 2019; Benisty et al. 2021). This kind of information offers unique insights on how much material is available to form planetary satellites and provides important boundary conditions for understanding the planetary accretion process.

In this article, we describe an attempt to quantify the constraints on any millimeter continuum emission from CPDs associated with the disk gaps identified by the Disk Substructures at High Angular Resolution Project (DSHARP; Andrews et al. 2018). Section 2 presents the selection criteria and the data used in this effort. Section 3 develops a methodology to mitigate confusion from local circumstellar disk emission. Section 4 describes the techniques used to quantify the sensitivity to CPD emission, and Section 5 presents the results. Finally, Section 6 discusses the outcomes of this analysis in the context of simple CPD models and the future prospects for further work. Section 7 summarizes the principal results.

## 2. Data

### 2.1. Selection Criteria

The full DSHARP sample includes 20 targets observed with ALMA at 240 GHz (1.25 mm) to a roughly uniform depth and angular resolution (Andrews et al. 2018). Our focus in this article was a search for CPDs in the primarily symmetric examples of annular gaps in the millimeter continuum emission distributions of these disks. Therefore, we excluded targets where the emission has a dominant nonaxisymmetric morphology (i.e., the cases with global spiral patterns; see Huang et al. 2018c; Kurtovic et al. 2018). We further limited the sample to include only gaps that are at least marginally resolved over the full azimuthal range (Huang et al. 2018b). The latter criterion excluded cases with high inclination angles (e.g., the DoAr 25 or HD 142666 disks), presumably because their emission surfaces are vertically elevated.

These selection priorities were entirely practical, designed only to facilitate the CPD search methodology discussed below. They do not imply that excluded cases are any less likely to host CPDs. The resulting sample includes nine disk

<table><caption>**Table 1**
Sample and Data Summary</caption><thead><tr><th>Disk<br>(1)</th><th>d<br>(pc)<br>(2)</th><th>Gaps<br>(3)</th><th>θ<sub>bm</sub><br>(mas)<br>(4)</th><th>PA<sub>bm</sub><br>(deg)<br>(5)</th><th>rms<br>(μJy beam<sup>-1</sup>)<br>(6)</th></tr></thead><tbody><tr><td>SR 4</td><td>135</td><td>D11</td><td>35 × 34</td><td>5</td><td>18</td></tr><tr><td>RU Lup</td><td>158</td><td>D29</td><td>26 × 25</td><td>145</td><td>14</td></tr><tr><td>Elias 20</td><td>138</td><td>D25</td><td>31 × 23</td><td>76</td><td>10</td></tr><tr><td>Sz 129</td><td>160</td><td>D41, D64</td><td>44 × 31</td><td>87</td><td>12</td></tr><tr><td>HD 143006</td><td>167</td><td>D22, D51</td><td>48 × 45</td><td>51</td><td>10</td></tr><tr><td>GW Lup</td><td>155</td><td>D74</td><td>45 × 43</td><td>1</td><td>10</td></tr><tr><td>Elias 24</td><td>139</td><td>D57</td><td>37 × 34</td><td>88</td><td>12</td></tr><tr><td>HD 163296</td><td>101</td><td>D48, D86</td><td>48 × 38</td><td>79</td><td>19</td></tr><tr><td>AS 209</td><td>121</td><td>D61, D97<sup>a</sup></td><td>38 × 37</td><td>106</td><td>10</td></tr></tbody></table>

**Notes.** Column (1): target name. Column (2): distance based on the Gaia EDR3 parallax (Gaia Collaboration et al. 2016, 2021). Column (3): Huang et al. (2018b) gap designation (for reference, “D” refers to a “dark” feature, and the associated number corresponds to the radius of the gap center in au). Columns (4) and (5): FWHM dimensions and position angle of the synthesized ALMA beam. Column (6): rms noise in the synthesized image, measured in an annulus 1.2 R<sub>out</sub> (see Table 2 and associated discussion) to 4"/25 from disk center.

a Following the interpretations of Guzmán et al. (2018) and Zhang et al. (2018), we chose to designate the wide outer gap in the AS 209 disk as “D97,” effectively combining the D90 and D105 gap designations of Huang et al. (2018b).

targets but focuses the CPD searches on 13 individual gaps. Table 1 lists some basic properties of the data associated with the search.

### 2.2. Data Processing

For each disk, we retrieved the publicly available self-calibrated (pseudo)continuum visibilities from the measurement sets in the DSHARP data repository<fnref n="13" /> (see Andrews et al. 2018 for calibration details). Those visibilities were spectrally averaged to one channel per spectral window and time-averaged into 30 s intervals to reduce the data volume. This averaging is modest enough to avoid bandwidth- or time-smearing effects. All of the post-processing of these data was conducted with the CASA v5.7 package (McMullin et al. 2007).

Throughout the analysis described below, imaging associated with these visibility data was performed following the DSHARP script recommendations, but with three minor modifications. First, we used 2× larger pixel sizes (6 mas, or ~0.2× the typical FWHM of the point-spread function (PSF)). Second, we adopted deeper CLEAN thresholds, equivalent to twice the rms values listed in Table 1. And third, a final step to the imaging process was included to treat an intrinsic deficiency in the CASA/tclean algorithms associated with CLEAN beam restoration for data with complicated PSFs due to the combination of disparate antenna configurations. The issue and its solution are described by Jorsater & van Moorsel (1995) and Czekala et al. (2021). To summarize, an additional correction step rescales the CLEAN residuals by the ratio $\epsilon$ of the Gaussian CLEAN beam to the PSF before adding them to the convolved CLEAN model, to ensure that the appropriate units (Jy per CLEAN beam) are propagated. We found a typical rescaling factor $\epsilon \approx 0.65$ (with extremes of 0.4 and 0.8). The updated rms values are listed in Table 1, measured as described by Andrews et al. (2018).

<footnote n="12">¹² The morphology of the circumplanetary material is uncertain (Szulágyi et al. 2016), but for simplicity we refer to it as a “disk.”</footnote>

<footnote n="13">¹³ https://almascience.org/almadata/lp/DSHARP</footnote>
### 3. Mitigating Contamination from Circumstellar Disk Emission

Contrast with the local circumstellar disk emission is the most formidable challenge for identifying faint CPDs embedded in narrow gaps. The ambient material emits a comparatively bright continuum, and PSF convolution smears that emission into the gap and complicates the CPD search. Moreover, CLEAN artifacts produced by the sparse distributions of ALMA antennas with long baselines might mimic CPD emission (e.g., manifested as “spokes” in gaps; see Andrews et al. 2016, 2018). Ideally, these risks can be mitigated by searching for CPDs in *residual* products, where the emission from the circumstellar disk has been removed. The following sections describe our approach to this task and its outcomes.

#### 3.1. Modeling Procedure

In this specific context, the detailed morphology of the circumstellar disk emission is irrelevant so long as there is some model available that does a good job of accounting for (and thereby removing) it. This means that a sophisticated statistical inference of model parameters is not desirable. That is helpful, since developing an appropriate parameterization and comparing it to the data can be technically challenging and incredibly time-consuming. Instead, we adopted the approach introduced by Jennings et al. (2020, 2021) and implemented in the associated *frank* software package. Assuming an intrinsically one-dimensional (radial) emission distribution, *frank* uses a Gaussian process to rapidly optimize a model of the (deprojected, real) visibilities. The fundamental assumption is that this axisymmetric modeling excludes CPD emission by default: because it should be both faint and azimuthally localized, the one-dimensional nature of the model averages out any CPD contribution at that radius. We confirmed a posteriori that this is indeed the case (see Section 4).

To construct such a model, we first explored the effects of the *frank* hyperparameters, fixing the disk geometry parameters (see below) to the values measured by Huang et al. (2018b). The $R_{max}$ hyperparameter (beyond which *frank* assumes that there is zero emission) was set to twice the radius of the outer extent of the observed emission, $R_{out}$, measured in the synthesized image and defined by the contour that reaches twice the rms noise (see Table 1). The $N$ and $p_0$ hyperparameters, which define the radial gridding and regularize the emission power spectrum, were kept at their default values ($300$ and $10^{-15}$ Jy$^2$; see Jennings et al. 2020). Combinations over a range of the hyperparameters $\alpha$ and $w_{smooth}$, which act like prior weights to regularize variations in the model brightness distribution, were explored interactively. We favored combinations that damp high-frequency oscillations in the brightness profile and set $\alpha = 1.3$ and $w_{smooth} = 0.1$ for all disks except AS 209, where instead $w_{smooth} = 0.001$. In any case, the outcomes for this study are not affected by these selections.

Having set the hyperparameters, we then refined estimates for the four geometric parameters that describe the emission morphology projection into the sky plane: inclination ($i$, where $i = 0^\circ$ for face-on and $90^\circ$ for edge-on), position angle of the major axis (PA, measured E of N), and the R.A. and decl. offsets from the observed phase center ($\Delta x$ and $\Delta y$, where positive values are to the E and N, respectively). We relied on a visual selection of these parameters, examining the images synthesized from the residual visibilities for models that spanned small grids around the Huang et al. (2018b) estimates

<table><caption>**Table 2**

Circumstellar Disk Emission Geometries</caption><thead><tr><th>Disk<br>(1)</th><th>Δx<br>(mas)<br>(2)</th><th>Δy<br>(mas)<br>(3)</th><th>i<br>(deg)<br>(4)</th><th>PA<br>(deg)<br>(5)</th><th>R<sub>out</sub><br>(arcsec)<br>(6)</th></tr></thead><tbody><tr><td>SR 4</td><td>-60</td><td>-509</td><td>22</td><td>18</td><td>0.25</td></tr><tr><td>RU Lup</td><td>-17</td><td>+86</td><td>19</td><td>121</td><td>0.42</td></tr><tr><td>Elias 20</td><td>-52</td><td>-490</td><td>54</td><td>153</td><td>0.48</td></tr><tr><td>Sz 129</td><td>+5</td><td>+6</td><td>32</td><td>153</td><td>0.48</td></tr><tr><td>HD 143006</td><td>-6</td><td>+23</td><td>16</td><td>167</td><td>0.53</td></tr><tr><td>GW Lup</td><td>+0.5</td><td>-0.5</td><td>39</td><td>37</td><td>0.75</td></tr><tr><td>Elias 24</td><td>+107</td><td>-383</td><td>30</td><td>45</td><td>1.05</td></tr><tr><td>HD 163296</td><td>-3.5</td><td>+4.0</td><td>47</td><td>133</td><td>1.08</td></tr><tr><td>AS 209</td><td>+2</td><td>-3</td><td>35</td><td>86</td><td>1.20</td></tr></tbody></table>

**Note.** Column (1): target name. Columns (2) and (3): R.A. and decl. offsets from the phase center. Column (4): inclination angle. Column (5): position angle. Column (6): outer (radial) boundary of the observed emission (the deprojected contour radius for twice the rms noise in Table 1). Formal uncertainties on the geometric parameters were not measured, but approximate estimates based on the residual map appearances suggest a precision of ~1–2 mas in the offsets and ~1° in the projection angles (up to perhaps ~2°–3° for the more face-on geometries). The estimated $R_{out}$ values have a ~10% uncertainty.

(using 0.5 mas and 0°.5 steps for the offsets and projection angles, respectively). Unfortunately, automatic optimization approaches to this task performed poorly in the general case of lingering asymmetries, but such obstacles can be overcome with some (admittedly tedious) experience. Some insights on the relevant issues are shared in Appendix A. Table 2 lists the adopted geometric parameters.

Finally, the visibility data were modeled using *frank* with the adopted hyperparameters and disk geometry parameters fixed. For two of the nine targets, the disks around Elias 20 and GW Lup, we followed the “point-source correction” procedure outlined by Jennings et al. (2021). This involved removing central point-like contributions of 0.61 and 0.48 mJy (as measured for deprojected baselines longer than 5.7 and 7.0 M$\lambda$), respectively, performing the *frank* modeling, and then adding those contributions back into the best-fit models. These corrections improved the quality of the model brightness profiles by damping oscillatory artifacts, but since we are focused on the widest gaps, they had no real effect on the interpretations outlined in Section 4. For two other targets, the disks around HD 163296 and HD 143006, we needed an additional step to account for their bright, azimuthally localized asymmetries (Isella et al. 2018; Pérez et al. 2018). The emission from these features was first removed based on an excised portion of the CLEAN model, following the procedure described in Appendix B. The results were *revised* visibility data sets for “symmetric” emission distributions, which were then modeled with *frank* as outlined above.

#### 3.2. Modeling Results

This approach generally performed well. To illustrate the quality of the disk emission removal, Figures 1 and 2 compare the images synthesized from the data, model, and residual visibilities for each disk on the same spatial and intensity scales. A direct comparison of the data and model visibilities, along with the inferred model brightness profiles, is presented in Appendix C.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Comparisons of the data, models, and residuals for the smaller circumstellar disks (see Figure 2 for the larger disks) in the image plane. The model and residual images were constructed from the *frank* model and residual visibilities in the same way as for the data. The CLEAN beam dimensions are marked in the lower left corner of each panel. The data, model, and residual images of a given disk are shown on the same brightness temperature scale (with an asinh stretch), assuming the Rayleigh-Jeans approximation. The HD 143006 images were constructed from the *revised* visibilities (see the text and Appendix B). Figure 3 shows a more detailed examination of the residual images.</figcaption></figure>

We used the model brightness profiles derived with *frank* to estimate the locations, widths, and depths of the gaps of interest, presuming a Gaussian morphology superposed on a background power-law profile. The parameters of this model—particularly the means ($r_{gap}$), standard deviations ($\sigma_{gap}$), and depletion factors (i.e., amplitudes; $\delta_{gap}$) of those Gaussian gaps—are cataloged in Table 3 (see Appendix C for details). These are meant only as rough guidance: many of the gaps are not Gaussian, and the backgrounds are not all described well with power laws. The gap centers are accurate, but the widths and depths are crude estimates.

Close inspections of the residual images reveal some lingering structure that was imperfectly captured by the *frank* modeling. Figure 3 shows a more instructive view of the images synthesized from the residual visibilities to illustrate those features. These images are on a common S/N scale, normalized by the map rms listed in Table 1, and are annotated to reference the locations of the gaps of interest and $R_{out}$. The gap boundaries

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Additional image-plane comparisons of the data, models, and residuals for the larger disks in this sample, as in Figure 1 (note the 2x larger spatial extent in each panel). The HD 163296 images were constructed from the *revised* visibilities (see the text and Appendix B).</figcaption></figure>

<table><caption>**Table 3**
Estimated Gap Properties</caption><thead><tr><th>Disk<br>(1)</th><th>Gap<br>(2)</th><th>r<sub>gap</sub><br>(mas)<br>(3)</th><th>σ<sub>gap</sub><br>(mas)<br>(4)</th><th>δ<sub>gap</sub><br>(5)</th><th>T<sub>0</sub><br>(K)<br>(6)</th><th>q<br>(7)</th></tr></thead><tbody><tr><td>SR 4</td><td>D11</td><td>79</td><td>10</td><td>30</td><td>17</td><td>0.70</td></tr><tr><td>RU Lup</td><td>D29</td><td>179</td><td>10</td><td>2</td><td>25</td><td>0.80</td></tr><tr><td>Elias 20</td><td>D25</td><td>181</td><td>11</td><td>1.8</td><td>21</td><td>1.25</td></tr><tr><td>Sz 129</td><td>D41</td><td>243</td><td>20</td><td>1.3</td><td>12</td><td>0.90</td></tr><tr><td>Sz 129</td><td>D64</td><td>376</td><td>20</td><td>1.8</td><td>170</td><td>3.30</td></tr><tr><td>HD 143006</td><td>D22</td><td>140</td><td>35</td><td>20</td><td>7</td><td>0.50</td></tr><tr><td>HD 143006</td><td>D51</td><td>315</td><td>26</td><td>5</td><td>7</td><td>0.50</td></tr><tr><td>GW Lup</td><td>D74</td><td>485</td><td>14</td><td>10</td><td>10</td><td>1.25</td></tr><tr><td>Elias 24</td><td>D57</td><td>420</td><td>35</td><td>40</td><td>23</td><td>0.75</td></tr><tr><td>HD 163296</td><td>D48</td><td>490</td><td>40</td><td>50</td><td>39</td><td>0.70</td></tr><tr><td>HD 163296</td><td>D86</td><td>850</td><td>41</td><td>12</td><td>48</td><td>0.90</td></tr><tr><td>AS 209</td><td>D61</td><td>510</td><td>32</td><td>25</td><td>15</td><td>0.70</td></tr><tr><td>AS 209</td><td>D97</td><td>800</td><td>62</td><td>30</td><td>15</td><td>0.70</td></tr></tbody></table>

**Note.** Column (1): target name. Column (2): gap designation from Huang et al. (2018b). Column (3): gap center. Column (4): gap width (Gaussian standard deviation). Column (5): gap depth (multiplicative depletion factor). Column (6): normalization of the local background power-law profile (brightness temperature at 0"/1). Column (7): background profile power-law index. See Appendix C for the model and parameter definitions.

are intended to illustrate what is seen in the (beam-convolved) data images; they correspond to $r_{gap} \pm (\sigma_{gap} + \sigma_{bm})$, where $\sigma_{bm} = \langle \theta_{bm} \rangle / 2\sqrt{2\ln 2}$ represents the (geometric mean) width of the synthesized beam (see Table 1).
<figure><img src="image_6.png" /><figcaption>**Figure 3.** A more detailed view of the images synthesized from the residual visibilities. These are the same images as in the rightmost panels in Figures 1 and 2, but shown here together with the same linear stretch on a diverging color scale for each target. The scale represents the residual S/N, where each map has been normalized by the rms noise determined in an annulus spanning 1.2 R<sub>out</sub> to 4.25 from the disk center (with values marked in the lower right corner of each panel in µJy beam⁻¹ units; see Table 1). These representations better illustrate the low-level asymmetries that are not accounted for in the frank modeling. Solid black ellipses (and gray for the four cases with multiple gaps) mark the gap “boundaries,” defined as described in the text. The dashed gray ellipse denotes the outer boundary of the millimeter continuum emission (R<sub>out</sub> in Table 2). The purple dashed boxes mark regions 0.3 on a side centered on the peak residual in each gap; a gallery of zoomed-in views of these boxes is shown in Figure 7. The synthesized beam dimensions are shown as filled black ellipses in the lower left corner of each panel.</figcaption></figure>

Figure 3 demonstrates that the model residuals generally have low S/N (white is |S/N| ≤ 3 in these images), particularly within the gap regions of interest here. Most cases show lingering residuals at ~4–7× the rms noise near the disk centers (within ~0.1′). In some sense, this is expected for models with limited fidelity. We found that the *maximum* imaged residuals from this methodology were typically only (±) a few percent of the observed emission levels in the data inside 0.1′ (increasing to ~10% in the fainter outer regions). However, the emission peaks in these regions have S/N ≈ 200–400, so a few percent deviation will produce a residual S/N of order ~10. These inner disk features might be caused by radiative transfer effects (i.e., a vertically flared

emitting surface; see below and Appendix A) or genuine asymmetries, but in most cases it would require deeper observations (ideally at higher resolution) to know for sure.

### 3.3. Notable Residuals

The residual maps for three cases merit specific mention (see also Jennings et al. 2021). The SR 4 disk exhibits some of the strongest residuals (~8× the rms noise). This is the most compact disk in the sample, with the innermost gap of interest here. The scale of these features suggests that they are analogs of the inner disk residuals noted above for other targets, although it is interesting that they primarily reside at the outer
edge of the gap. Perhaps that is a clue that the residuals generally found close to the disk centers are associated with nonaxisymmetric substructures that are not easily recognized because they suffer from poor resolution and/or PSF convolution artifacts. The HD 163296 disk also exhibits strong and structured residuals, both in the inner ~0.′′1 and associated with the bright ring (B67) separating the D48 and D86 gaps. The morphological pattern of these residuals suggests a ring surface that is inclined and/or elevated with respect to the global mean values (e.g., Appendix A; see also Isella et al. 2018). Doi & Kataoka (2021) reached similar conclusions from a detailed modeling analysis of the same data. A comparable, albeit much weaker, pattern is visible well outside the D57 gap for the Elias 24 disk, perhaps associated with an analogous deviation in the projected emission surface height for the B77 bright ring.

The residuals for the HD 143006 disk are especially complex in Figure 3. A close examination of the residual map shows a noisy, one-armed spiral pattern (the ± residual patterns in Figure 3 are artifacts of the axisymmetric frank modeling: the negative [blue] residuals are the mirror image of the spiral feature). This is perhaps illustrated more clearly in Figure 4, where the polar deprojection (in the disk-plane polar coordinate system) of the image is also included. A (visually tuned, not fitted) Archimedean spiral<fnref n="14" /> morphology approximates this residual feature fairly well, with $r_{\text{spiral}} \approx 0.170 + 0.067\varphi$ (in units of arcseconds), where the azimuth $\varphi$ (in radians) follows the convention outlined by Huang et al. (2018b). This spiral apparently traverses through all of the rings and gaps in the disk, spanning ~420° in azimuth; it is accompanied by lingering clumpy residuals that are spatially coincident with the southeastern quadrants of the B41 and B65 rings. Pérez et al. (2018) noted the innermost portion of the feature as a “bridge” between the B8 and B41 rings (at PA ≈ 300°). Their residual map has a more muted peak at the innermost extension of the arm, in part because they modeled the B8 ring with a different *i* (i.e., a warp). For the orientation favored by interpretations of the scattered light shadows (Benisty et al. 2018) and gas kinematics (Pérez et al. 2018), with the east side of the disk nearer to the observer, the spiral is trailing. Given the low S/N of this spiral pattern, it is unclear what to make of this feature. We offer some brief speculation in Section 6, but deeper observations that can facilitate a more quantitative analysis are needed before drawing any conclusions.

## 4. Assessing Sensitivity to CPD Emission

After the bright (and presumed symmetric) emission from the circumstellar disk was modeled and removed, the residual images presented in Figure 3 could be used to search for faint CPD emission in the gaps of interest. In our analysis, we assumed that the regions most likely to exhibit CPD emission were the annular zones bounded by $r_{\text{gap}} \pm \sigma_{\text{gap}}$ in a disk-frame coordinate system. Moreover, we implicitly assumed that the optical depths in these gap regions were low enough that they did not obscure any potential CPD emission (which would then be missing after the frank model subtraction). This section presents the methodology adopted for a statistical assessment of our sensitivity to CPD emission in those search zones. An analysis of the outcomes from those assessments and evaluations of the residual peaks are presented in Section 5.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** The residual image (top) and its deprojection into disk-frame polar coordinates (bottom), following the formalism of Huang et al. (2018b), on a linear intensity scale to highlight the low-level spiral. The bright ring emission peaks and the region around the pronounced azimuthal asymmetry are marked in the sky-plane (polar-deprojected) image with white ellipses (vertical lines) and an arc (rectangle), respectively. For reference, the straight dashed lines in each panel mark zero azimuth ($\varphi = 0°$; with increasing $\varphi$ clockwise in the sky plane). The spiral model described in the text is overlaid as a dashed yellow curve, to guide the eye.</figcaption></figure>

A simplistic upper limit based only on the rms residual scatter in a gap annulus is generally an inappropriate estimator for the sensitivity to CPD emission, since the distribution of pixel intensities is often too correlated (i.e., the gap annuli are not covered by a sufficient number of independent resolution elements to use standard Gaussian statistics) and can be significantly skewed by lingering nonaxisymmetric residuals. Instead, we opted for a statistical approach that assessed the ability to faithfully recover injected (mock) CPD signals.

For this purpose, we assumed a point-source model for the mock CPD emission. Theoretical models suggest that CPDs are truncated at radii of ~0.1–0.5 $R_{\text{H}}$, where the Hill radius $R_{\text{H}} \approx a_p (M_p/3M_*)^{1/3}$ (e.g., Quillen & Trilling 1998; Martin & Lubow 2011). Presuming giant planets ($M_p \approx M_{\text{Jup}}$) near the gap centers ($a_p \approx r_{\text{gap}}$), any CPDs in this sample are indeed expected to have diameters smaller than the ALMA resolution.

<footnote n="14">¹⁴ We were unable to find a logarithmic spiral model that performed well across the full azimuthal range of the feature.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** An image-plane illustration of one iteration in the mock CPD injection-recovery analysis for the Elias 24 D57 gap. The three panels on the left show the data with an injected CPD signal, the axisymmetric frank model, and the corresponding residuals, as in Figure 2 (but with scale and stretch adjusted for clarity). The rightmost panel shows the residual S/N map, as in Figure 3. The search annulus corresponding to $r_{gap} \pm \sigma_{gap}$ is marked with black ellipses. This image has a more compact spatial range, as marked by the dashed square in the third panel from the left. This residual S/N map more clearly highlights the $F_{cpd} \approx 100 \mu Jy$ mock CPD to the northwest (circled in purple here, and white in the other panels). Note the low-level “spoke” artifacts crossing the gap in both the data and model panels, affirming the concerns that motivated the modeling efforts in Section 3. The mock CPD cannot be recovered in the original data image without that modeling and removal procedure.</figcaption></figure>

A mock CPD model has three parameters: a location ($r_{cpd}$, $\phi_{cpd}$; polar coordinates in the disk frame) and a flux density ($F_{cpd}$). The mock CPD visibilities for a given set of parameters were computed from the Fourier transform of an offset point source, as described in Appendix D.

Each iteration of the statistical assessment framework we followed has five basic steps:

1. Assign mock CPD parameters. $F_{cpd}$ was selected from a grid spanning 10–250 $\mu Jy$ at 10 $\mu Jy$ intervals. Locations were randomly drawn from (independent) uniform distributions such that $r_{cpd} \in (r_{gap} \pm \sigma_{gap})$ and $\phi_{cpd} \in (\pm 180^{\circ})$.

2. Compute mock CPD visibilities (see Appendix D) and add them to the observed ALMA visibilities.

3. Model these composite data and derive a set of residual visibilities as described in Section 3.

4. Image the residual visibilities (see Section 2).

5. Measure the peak in the search zone of the residual image and compare it to the mock CPD inputs.

To build up the mock catalog for each gap, these steps were iterated 500× per $F_{cpd}$ value. Figure 5 illustrates a typical outcome for this procedure.

While step 1 is self-explanatory, and steps 2 and 3 were described above, the details of steps 4 and 5 merit some discussion. Readers may (reasonably) wonder why we analyzed residuals in the image plane, rather than performing the search in the Fourier domain. Indeed, testing the recovery process for mock CPD injections onto a pure Gaussian noise distribution found that the latter option performed equally well, with a reduced computational cost. However, a visibility-based forward-modeling approach to the recovery became problematic in reality, where asymmetric residuals severely bias the outcomes (even at modest S/N). Attempts to guide the search to ignore those residuals (with increasingly sophisticated priors) ended up becoming complicated efforts to “model the noise,” which encouraged us to instead opt for the simpler image-based alternative.

Ideally, forward-modeling the visibilities has the advantage of a well-defined metric for a successful recovery (in terms of a posterior probability). But in the image-based case, “success” is subjective. We used an astrometric criterion, requiring that the peak residual lies within a (sky-projected) distance $\Delta R$ of the injected CPD location. The adopted distance tolerance was

$$ \Delta R = \max(2\sigma_{\text{ast}}, 12 \text{ mas}), \quad \text{where} \\ \sigma_{\text{ast}} = 70'' \left[ \left( \frac{\nu}{\text{GHz}} \right) \left( \frac{B_{\text{max}}}{\text{km}} \right) \left( \frac{I_{\text{peak}}}{\text{RMS}_{\text{gap}}} \right) \right]^{-1}. \quad (1) $$

Here, $\sigma_{\text{ast}}$ is the formal astrometric uncertainty,<fnref n="15" /> 12 mas corresponds to 2 image pixels, $\nu$ is the frequency, $B_{\text{max}}$ is the longest baseline, and $I_{\text{peak}}/\text{RMS}_{\text{gap}}$ is the recovered peak S/N (see below for details).

The residual imaging aspect of the recovery procedure (step 4 above) was especially time-consuming. It amounted to running the CASA/tclean algorithm some 12,500 times for each gap of interest. Extensive experimentation demonstrated that some helpful speed gains can be made by imaging with a mask that covered only around the appropriate search annulus ($r_{gap} \pm 1.5\sigma_{gap}$), without any negative effects on the recoveries. This required us to modify the adopted scales parameter in each case, starting from “point-like” (scales=0) and incrementing by 5-pixel intervals (~1 FWHM of the PSF) until exceeding the mask annulus width ($3\sigma_{gap}$).

## 5. Results

The key results of the mock CPD injection-recovery exercise described in Section 4 are summarized in Figure 6. These profiles show the fraction of mock CPDs that were recovered (using the criterion described in Equation (1)) as a function of the injected flux density.

As a guide to interpret these profiles, consider an oversimplified approximation where the non-CPD residuals are pure thermal noise (i.e., the frank modeling is perfect, and there is no phase noise or imaging artifacts). In that scenario, realizations of mock CPD injections result in point-like residuals with peak intensities consistent with random draws from a Gaussian distribution that has a mean equivalent to the injected flux density, $F_{cpd}$, and a standard deviation that reflects

<footnote n="15">¹⁵ The definition here is based on the ALMA Technical Handbook (https://almascience.nrao.edu/documents-and-tools/cycle7/alma-technical-handbook/view), although alternatives (e.g., Reid et al. 1988) give similar results.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 6. The fraction of injected mock CPDs that are successfully recovered as a function of their flux densities. Each point corresponds to 500 mock injections; error bars reflect only the Poisson uncertainties. In disks with multiple gaps under consideration, the recovery profiles for the innermost gaps are red (as labeled in the lower right corner of each panel). The 50% recovery fractions and their corresponding flux densities are marked with dotted lines for each gap of interest.</figcaption></figure>

the “local” noise, RMS<sub>gap</sub> (defined as the rms of pixel values within the search annulus of the residual map).<sup>16</sup> If a realization generates a mock CPD signal brighter than the *actual* peak residual in the search annulus, the recovery is successful. Therefore, the 50% recovery fraction occurs where the actual peak residual is ≈F<sub>cpd</sub>.

Despite the assumptions being technically invalid, the above approximation is a reasonable reflection of the actual recovery outcomes. Table 4 lists the locations and intensities for the actual peak residuals in both the sky-frame and disk-frame coordinate systems, along with the dispersions in the search annuli (rms<sub>gap</sub>) and some CPD flux densities that correspond to representative recovery fractions (interpolated from the profiles in Figure 6). The measured rms<sub>gap</sub> values are generally ~1.5–2× higher than found for a large, empty area (the rms in Table 1), primarily due to lingering nonaxisymmetric residuals from the disk modeling. However, especially in the cases where there are few independent resolution elements sampling the search annulus (i.e., for narrow gaps and/or gaps located at small radii), the pixel distributions are not Gaussian (and strongly covariant) and therefore also tend to artificially inflate rms<sub>gap</sub>.

We estimated the “false-positive” probability of a recovery that is actually a thermal noise peak (i.e., drawn from a Gaussian with mean zero and dispersion rms<sub>gap</sub>) and not the injected CPD signal as a function of F<sub>cpd</sub>. There are two factors to consider: the probability of the peak falling within ΔR (i.e., the areal ratio of the recovery region and the gap search annulus), and the probability of the peak being brighter than F<sub>cpd</sub>. The false-positive recovery fraction can be approximated as the product of these factors,

$$ \approx \left[ \frac{\max\left(\frac{\pi \langle \Delta R \rangle^2}{\Omega_{bm}}, 1\right)}{4\pi \sigma_{gap} r_{gap} \Omega_{bm}^{-1}} \right] \left[ \frac{1}{2} \operatorname{erfc}\left(\frac{F_{cpd}}{2 \text{ RMS}_{gap}}\right) \right], \quad (2) $$

where erfc is the complementary error function, ⟨ΔR⟩ is the average distance threshold from Equation (1) for each discrete flux density value, and Ω<sub>bm</sub> is the solid angle of the synthesized beam.<sup>17</sup> This kind of false positive from random noise peaks is unimportant in our analysis, typically accounting for ≪1% of the recoveries in the lowest flux density bins (F<sub>cpd</sub> ≲ RMS<sub>gap</sub>).

<footnote n="16">16 A note on nomenclature and units: for the point-like residuals we assume here, the flux density (in μJy) is equivalent to the peak intensity (in μJy beam<sup>-1</sup>) since the solid angle containing the signal is by default equivalent to the beam area.</footnote>

<footnote n="17">17 If each pixel were independent, we could write the first factor in Equation (2) more simply as (ΔR)<sup>2</sup>/4σ<sub>gap</sub>r<sub>gap</sub>. However, the version in the text reflects the spatial covariances in the residual images imparted by PSF convolution. The differences between the two approximations are negligible for the cases presented here.</footnote>
<table><caption>**Table 4**

Peak Residuals and Mock CPD Recoveries</caption><thead><tr><th rowspan="2">Disk<br>(1)</th><th rowspan="2">Gap<br>(2)</th><th rowspan="2">Δρ<sub>peak</sub><br>(mas)<br>(3)</th><th rowspan="2">PA<sub>peak</sub><br>(deg)<br>(4)</th><th rowspan="2">t<sub>peak</sub><br>(mas)<br>(5)</th><th rowspan="2">φ<sub>peak</sub><br>(deg)<br>(6)</th><th rowspan="2">I<sub>peak</sub><br>(μJy<br>bm)<br>(7)</th><th rowspan="2">rms<sub>gap</sub><br>(μJy<br>bm)<br>(8)</th><th colspan="5">F<sub>cpd</sub> (μJy) for a Given Recovery Fraction</th></tr><tr><th>0.5<br>(9)</th><th>0.6<br>(10)</th><th>0.7<br>(11)</th><th>0.8<br>(12)</th><th>0.9<br>(13)</th><th>1.0<br>(14)</th></tr></thead><tbody><tr><td>SR 4</td><td>D11</td><td>84</td><td>4</td><td>84</td><td>+105</td><td>127</td><td>35</td><td>118</td><td>127</td><td>138</td><td>150</td><td>178</td><td>230</td></tr><tr><td>RU Lup</td><td>D29</td><td>172</td><td>55</td><td>180</td><td>+157</td><td>54</td><td>19</td><td>55</td><td>60</td><td>68</td><td>75</td><td>86</td><td>130</td></tr><tr><td>Elias 20</td><td>D25</td><td>120</td><td>280</td><td>178</td><td>-24</td><td>55</td><td>18</td><td>60</td><td>65</td><td>72</td><td>81</td><td>88</td><td>120</td></tr><tr><td>Sz 129</td><td>D41</td><td>192</td><td>242</td><td>226</td><td>+1</td><td>44</td><td>16</td><td>44</td><td>49</td><td>56</td><td>63</td><td>72</td><td>160</td></tr><tr><td>Sz 129</td><td>D64</td><td>370</td><td>300</td><td>390</td><td>-53</td><td>51</td><td>15</td><td>53</td><td>56</td><td>60</td><td>68</td><td>82</td><td>150</td></tr><tr><td>HD 143006</td><td>D22</td><td>103</td><td>331</td><td>103</td><td>-73</td><td>70</td><td>25</td><td>71</td><td>79</td><td>88</td><td>97</td><td>111</td><td>220</td></tr><tr><td>HD 143006</td><td>D51</td><td>326</td><td>97</td><td>338</td><td>+161</td><td>78</td><td>23</td><td>85</td><td>92</td><td>98</td><td>107</td><td>120</td><td>185</td></tr><tr><td>GW Lup</td><td>D74</td><td>446</td><td>5</td><td>486</td><td>+129</td><td>46</td><td>16</td><td>55</td><td>61</td><td>68</td><td>77</td><td>91</td><td>205</td></tr><tr><td>Elias 24</td><td>D57</td><td>341</td><td>137</td><td>394</td><td>-2</td><td>55</td><td>13</td><td>62</td><td>68</td><td>75</td><td>82</td><td>91</td><td>175</td></tr><tr><td>HD 163296</td><td>D48</td><td>488</td><td>157</td><td>532</td><td>+57</td><td>118</td><td>36</td><td>124</td><td>135</td><td>147</td><td>162</td><td>188</td><td>>250</td></tr><tr><td>HD 163296</td><td>D86</td><td>591</td><td>226</td><td>861</td><td>-2</td><td>104</td><td>31</td><td>112</td><td>124</td><td>133</td><td>143</td><td>156</td><td>235</td></tr><tr><td>AS 209</td><td>D61</td><td>428</td><td>347</td><td>520</td><td>-173</td><td>52</td><td>16</td><td>56</td><td>60</td><td>68</td><td>77</td><td>86</td><td>175</td></tr><tr><td>AS 209</td><td>D97</td><td>689</td><td>330</td><td>814</td><td>-159</td><td>65</td><td>15</td><td>70</td><td>77</td><td>83</td><td>89</td><td>97</td><td>155</td></tr></tbody></table>

**Note.** Column (1): target name. Column (2): gap designation from Huang et al. (2018b). Columns (3) and (4): residual peak location in sky-frame coordinates (projected radial distance from the disk center and position angle E of N). Columns (5) and (6): residual peak location in disk-frame polar coordinates (following the azimuth convention of Huang et al. 2018b). Column (7): peak residual brightness. Column (8): standard deviation of pixel values in the gap. Columns (9)–(14): minimum mock CPD flux density that is recovered in 50%, 60%, 70%, 80%, 90%, and 100% (respectively) of the injection trials (based on spline interpolations of the profiles in Figure 6).

The only exception is for the narrow gap D11 in the SR 4 disk, where the false-positive rate is consistent with the recovery rate (~few percent) at $F_{cpd} \le 40$ μJy.

The underlying causes of the unsuccessful recoveries in this adopted methodology are simply the actual residual peaks in the data. These peaks are shown in the zoomed-in residual S/N maps of Figure 7 and have their properties cataloged in Table 4. These peaks can be attributed to various origins, including noise, features associated with the circumstellar disk, or genuine CPDs. As we noted in Section 3, the gaps in the SR 4, HD 143006, and HD 163296 disks exhibit residual peaks that seem linked to the local nonaxisymmetric behaviors of their circumstellar disks. Despite these being perhaps the most compelling targets for planet–disk interactions, differentiating between these measured features and any point-like residuals that could be associated with CPDs will remain a challenge without deeper observations and improved disk emission models.

The remaining cases suggest marginal residual peaks, often with a local S/N ≈ 3 in their search annuli. For the narrow gaps in the RU Lup, Elias 20, Sz 129, and GW Lup disks, peaks with such “significance” are not necessarily expected (assuming that the noise distribution is Gaussian) for search annuli with areas comparable to only ~30–60 synthesized beam solid angles. However, a close examination of the residuals in adjacent zones suggests that these are indeed consistent with local noise peaks. The peaks in the Elias 24 D57 and AS 209 D97 gaps are slightly higher S/N cases that merit some follow-up. We speculate that the peak in AS 209 D97 is probably associated with the very faint ring identified near the gap center (Guzmán et al. 2018; Huang et al. 2018b). We do not find any obvious residual signal near the location of the Elias 24 planet candidate identified by Jorquera et al. (2021), nor in the vicinities of the CO kinematic perturbations in the HD 143006, HD 163296, GW Lup, or Sz 129 disks identified by Pinte et al. (2020).

These residual peaks are considered candidate CPDs with marginal confidence until they can be pursued with more

sensitive measurements. Since their peak intensities are similar to the $F_{cpd}$ values that are recovered for 50% of the injection–recovery experiments (Table 4), we adopted the latter as a homogenized approximation for the “upper limits” on CPD flux densities.

# 6. Discussion

## 6.1. Limits on CPD Masses

These derived upper limits on the millimeter continuum flux densities can be translated into analogous constraints on CPD (dust) masses. We largely followed the approach of Isella et al. (2014, 2019), although comparable results were found based on the methods of Zhu et al. (2018). The simplified model for the CPD emission is described in Appendix E. This model permits an estimate of $M_{cpd}$, defined as the mass corresponding to the flux density upper limit (the $F_{cpd}$ with a 50% recovery probability using the analysis described in Section 4; see Column (9) in Table 4). We can alternatively think of $M_{cpd}$ as the “minimum detectable” mass, given the DSHARP data and our adopted assumptions and search methodology. These mass limits were calculated as a function of the two key unknowns: the planet mass, $M_p$, and the CPD radius (in units of the Hill radius), $R_{cpd}/R_H$.

The various other parameters of the CPD model were either fixed (see Appendix E) or crudely explored around representative boundary values. A more detailed discussion of the expected uncertainties associated with the “fixed” choices is available in Appendix E. The CPD model densities are defined by the key input ($M_p$, $R_{cpd}$) and output ($M_{cpd}$) parameters. The temperatures are based on approximations for irradiation heating by the planet and the host star (and local disk), as well as viscous heating from accretion. The irradiation contributions assumed planetary evolution models and an approximation for the local disk temperatures (e.g., Chiang & Goldreich 1997). We considered two cases for the viscous heating. In the first, accretion is not a significant contributor to the CPD thermal structure. We assumed a constant accretion
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Zoomed-in views of the purple dashed squares marked in the residual maps of Figure 3, showing the regions immediately surrounding the actual residual peaks in each gap. The annotations are as in Figure 3; purple dashed circles mark the peaks. The residual S/N scale is different here; the rms<sub>gap</sub> values in Table 4 are used to normalize the residual images, with a compressed S/N scale to emphasize the local background variations. Contours are drawn at 40, 60, 80, and 100 µJy beam<sup>-1</sup> (and their negative complements as dotted contours) to mark an absolute emission scale.</figcaption></figure>

rate for the planet “lifetime” ($t_p$), $\dot{M}_p \approx M_p/t_p$, with $t_p = 1$ Myr (the mean age of the host stars; Andrews et al. 2018). This is roughly consistent with the $\dot{M}_p$ limits derived from $H\alpha$ nondetections in recent direct imaging CPD searches (e.g., Cugno et al. 2019; Hashimoto et al. 2020; Zurlo et al. 2020). In the second case, we scaled up that rate $100 \times$ to simulate an active accretion phase with a thermal contribution more comparable to (or larger than) irradiation.

Moreover, we also explored two cases for the properties of the emitting dust grains in these model CPDs, as characterized by their absorption opacities ($\kappa_\nu$) and albedos ($\omega_\nu$). In one, we assumed that scattering is negligible ($\omega_\nu \approx 0$) and $\kappa_\nu = 2.4 \text{ cm}^2 \text{ g}^{-1}$ at the DSHARP data frequency (240 GHz). This is more or less the standard assumption used so far in the literature. The value for $\kappa_\nu$ is consistent with the classic Beckwith et al. (1990) opacity prescription. In another case, we explored the effects of strong scattering ($\omega_\nu = 0.9$) for the same $\kappa_\nu$, more appropriate when most of the mass is

concentrated in particles with sizes comparable to the observing wavelength (Birnstiel et al. 2018; Zhu et al. 2019).

The four components of Figure Set 8 illustrate how the CPD mass upper limits for each gap vary as a function of $M_p$ and $R_{\text{cpd}}/R_H$, corresponding to each of the four distinct sets of assumptions outlined above. Each panel also shows a magenta dotted curve that marks the (radial) resolution of the data (HWHM), and a black hatched region that denotes where the CPD size would be larger than twice the Gaussian standard deviation that describes the gap width ($2\sigma_{\text{gap}}$; see Table 3). Models that lie above the magenta curve violate our assumption of point-like CPD emission, and models within the black hatched region are “unphysical” in the sense that the CPD should not extend beyond the gap boundaries.<fnref n="18" /> The derived $M_{\text{cpd}}$ limits are in the range of 0.001–0.2 $M_\oplus$ (roughly 10% of a lunar mass to a Mars mass). There is only a modest dependence on $R_{\text{cpd}}/R_H$ (see Isella et al. 2019), but we see clear decreases in the limits for larger $M_p$ and $\dot{M}_p$. The latter behavior is caused by the associated heating: these CPD models are hotter, so less mass is required to produce a given flux.

Adjustments to the $M_{\text{cpd}}$ limits for different dust albedos are subtler (and still incomplete; e.g., scattering should also impact the CPD temperature structure, but this is not considered). Generally, slightly lower $M_{\text{cpd}}$ limits are found for higher $\omega_\nu$. The cause of that behavior is complicated, but related to the variation of the dust optical depths. In the optically thin limit, albedo has negligible effects. And in the purely optically thick case, more scattering (higher $\omega_\nu$) means less emission (e.g., Zhu et al. 2019). For the same physical parameters (temperatures and densities), a higher albedo means higher optical depths. Therefore, in the general case that spans from high to low optical depths with increasing radius, the transition from thick to thin occurs at larger radii for higher albedo. For our assumptions, this usually means higher flux densities for higher $\omega_\nu$, and therefore correspondingly lower $M_{\text{cpd}}$ limits.

For low $M_p$ and $R_{\text{cpd}}$ (the empty white regions in the lower left corners of these plots), the flux upper limits cannot be reproduced for any CPD mass with our adopted model assumptions. These swaths are larger for models that make disks colder or scattering more prevalent. In these regions of parameter space, the measured flux upper limits are higher than for a scenario where the entire CPD is optically thick: the data are not sensitive enough to find CPDs with these assumed properties.

To help contextualize these CPD mass limits, Figure Set 8 highlights two sectors of parameter space: the red hatched regions mark the $M_p$ ranges excluded (at 50% probability) by the direct imaging measurements of Jorquera et al. (2021) (for the AS 209 D97, Elias 24 D57, GW Lup D74, and Sz 129 D64 gaps), and the dashed red outlines correspond to the $M_p$ ranges predicted by the method of Zhang et al. (2018). The latter sectors are based on comparisons of the relative gap widths ($\Delta$; see Equation (21) in Zhang et al. 2018) inferred from the frank model brightness profiles (see Appendix C) and a suite of hydrodynamics simulations (both convolved with the same Gaussian kernel of width $0.06 r_{\text{gap}}$), assuming a dust population with a maximum particle size of 1 mm (commensurate with the

<footnote n="18">¹⁸ As we noted in Section 3.2, a Gaussian is not always the best representation of the gap emission distribution in the frank models. We hatched the region starting at twice the derived width (gap standard deviation) as a more conservative measure of the gap boundary to ensure that this ambiguity is not misleading.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The minimum detectable CPD dust masses as a function of the host planet mass and CPD radius (in Hill sphere radius units) for each target gap of interest, in this case presuming no scattering ($\omega_{\nu} \approx 0$) and modest accretion rates ($M_{\dot{p}} \approx M_{p}/1$ Myr). The regions above the dotted magenta curves correspond to CPDs larger than the resolution element. The black hatched regions correspond to CPDs larger than the gap size ($R_{\text{cpd}} > 2\sigma_{\text{gap}}$). The former technically violates our assumptions; the latter is unphysical. For the empty white regions, no amount of mass is sufficient to produce the corresponding millimeter flux density, given our assumptions (i.e., the optical depths are saturated). The red hatched regions correspond to the properties ruled out by the direct imaging searches by Jorquera et al. (2021). The red dashed rectangles mark the $M_{p}$ regions inferred with the simulations and methodology of Zhang et al. (2018) (see Table 5).</figcaption></figure>

(The complete figure set (4 images) is available.)

adopted $\kappa_{\nu}$) and a viscosity coefficient $\alpha = 0.001$. The corresponding $M_{p}$ (along with $\Delta$ and $a_{p}$) are compiled in Table 5, with uncertainties estimated as described by Zhang et al. (2018). For both the Jorquera et al. (2021) constraints and the Zhang et al. (2018) based masses, we highlighted the regions bound by $R_{\text{cpd}}/R_{\text{H}} \approx 0.1–0.5$ as representative of theoretical expectations (see Martin & Lubow 2011).

Generally, the mass limits derived here are similar to theoretical expectations for the CPDs predicted in simulations of gravitationally unstable disks (Stamatellos & Herczeg 2015) or the formation of the Galilean satellites (Canup & Ward 2002). For the SR 4, HD 143006, Elias 24, HD 163296, and AS 209 (D97) disk gaps, the Zhang et al. (2018) based predictions imply $M_{\text{cpd}}/M_{p} \lesssim 1\%–10\%$ (in some cases
<figure><img src="image_12.png" /><figcaption>**Figure 8.** (Continued.)</figcaption></figure>

≤0.01% or up to ~60%), presuming a global gas-to-dust ratio of 100. In the other cases, the flux limits are above the expectations for fully optically thick CPDs; limits on $M_{\text{cpd}}$ are unavailable if the planet masses in Table 5 are appropriate.

We should make an explicit reminder that these data and our analysis are only sensitive to CPD dust masses, and our opacity assumptions are specifically relevant for particles with approximately millimeter sizes. However, smooth CPD density (or, rather, pressure) structures like those assumed in our simple model are expected to facilitate the rapid inward radial drift of such particles (e.g., Shibaike et al. 2017). If there is no pressure modulation to “trap” those particles (Drążkowska & Szulágyi 2018; Batygin & Morbidelli 2020; Szulágyi et al. 2021), a CPD

could have a considerable total mass (in gas and small particles) but produce very little millimeter continuum emission. In that sense, ALMA observations might only be sensitive to the small grains still dynamically coupled to the gas, or to CPDs that have their own substructures.

## 6.2. *Empirical Context for CPD Limits*

Inferences of CPD mass constraints like those above should be considered with due caution: they still require many assumptions about physical properties where our knowledge is very limited. As an alternative, we can examine a more empirical context for these CPD constraints by comparing them with other samples in
<figure><img src="image_13.png" /><figcaption>**Figure 8.** (Continued.)</figcaption></figure>

the millimeter continuum luminosity–“host” mass domain, as illustrated in Figure 9. The upper limits derived in this article are comparable to those for planetary-mass companions (PMCs; Bowler et al. 2015; MacGregor et al. 2017; Ricci et al. 2017; Pineda et al. 2019; Wu et al. 2020), the isolated planet OTS 44 (Bayo et al. 2017), and the planet PDS 70 c (Isella et al. 2019), although often for much lower “host” masses (Zhang et al. 2018). While these limits rule out optically thick CPDs around the more massive putative planets (e.g., the gray dotted curves in Figure 9), they would lie well above (≥10–100×) the luminosities expected if the correlation identified for circumstellar disks (see the discussion by Andrews 2020) is extrapolated to lower host masses (black dashed line in Figure 9). Of course, since we do not yet

understand the physical origin of that correlation, there is not necessarily a good reason to expect that it would extend indefinitely into the planetary-mass regime.

Perhaps a more appropriate comparison to the PMC disk searches would be an exploration of putative CPDs in the DSHARP sample located at larger distances than were considered here, outside the continuum emission boundary ($a_p \gtrsim R_{\text{out}}$). In one sense, searches in these regions—where there are no nonaxisymmetric residuals from the host disk present—are much simpler. Using the formalism described in Section 4, we have confirmed that upper limits (50% recovery fractions) for point-like features are what we would naively expect, ~3× the map rms noise levels (listed in Table 1).
<figure><img src="image_14.png" /><figcaption>Figure 8. (Continued.)</figcaption></figure>

Those limits would be appropriate in the limiting case where efficient radial drift makes the distribution of emitting solids especially compact. But in the more general case, the real challenge is that theoretical models (e.g., Martin & Lubow 2011) suggest that CPDs at these distances will be spatially resolved with the DSHARP data: we would need to substantially modify the adopted search methodology to properly quantify flux limits for such CPDs.

One general assumption that we make is that CPDs will exhibit brighter millimeter emission if they have more massive planetary hosts, because there is more irradiation (and viscous) heating and (presumably) more mass. The dynamical perturbations such massive planets induce on the structures of the

circumstellar disks in which they are embedded are expected to be more pronounced than the narrow gaps observed in the DSHARP sample (i.e., see the *M*<sub>*p*</sub> values in Table 5). The large, cleared cavities of “transition” disks are nominally better search targets for these brighter CPDs around more massive planets (e.g., Salyk et al. 2009; Zhu et al. 2011). Clearly, the PDS 70 system offers the most compelling support for such a strategy (Isella et al. 2019; Benisty et al. 2021). However, it might also be the case that CPDs in transition disk cavities are relatively starved of the approximately millimeter-sized particles that emit most efficiently in the ALMA bands because the supply flow of those particles is throttled at the high-amplitude pressure maxima their host planets induce at the cavity edge
<table><caption>**Table 5**
Estimated Planet Properties</caption><thead><tr><th>Disk<br>(1)</th><th>Gap<br>(2)</th><th>a<sub>p</sub> (au)<br>(3)</th><th>log(M<sub>p</sub>/M<sub>lup</sub>)<br>(4)</th><th>Δ<br>(5)</th></tr></thead><tbody><tr><td>SR 4</td><td>D11</td><td>10.7</td><td>0.20<sup>+0.14</sup><sub>-0.17</sub></td><td>0.44</td></tr><tr><td>RU Lup</td><td>D29</td><td>28.3</td><td>-1.50<sup>+0.14</sup><sub>-0.17</sub></td><td>0.14</td></tr><tr><td>Elias 20</td><td>D25</td><td>25.0</td><td>-1.89<sup>+0.14</sup><sub>-0.17</sub></td><td>0.12</td></tr><tr><td>Sz 129</td><td>D41</td><td>38.9</td><td>-2.17<sup>+0.22</sup><sub>-0.16</sub></td><td>0.16</td></tr><tr><td>Sz 129</td><td>D64</td><td>60.2</td><td>...</td><td>...</td></tr><tr><td>HD 143006</td><td>D22</td><td>23.4</td><td>0.84<sup>+0.21</sup><sub>-0.50</sub></td><td>0.70</td></tr><tr><td>HD 143006</td><td>D51</td><td>52.6</td><td>-0.55<sup>+0.21</sup><sub>-0.50</sub></td><td>0.22</td></tr><tr><td>GW Lup</td><td>D74</td><td>75.2</td><td>-1.68<sup>+0.21</sup><sub>-0.50</sub></td><td>0.14</td></tr><tr><td>Elias 24</td><td>D57</td><td>58.4</td><td>-0.48<sup>+0.21</sup><sub>-0.50</sub></td><td>0.30</td></tr><tr><td>HD 163296</td><td>D48</td><td>49.5</td><td>0.02<sup>+0.21</sup><sub>-0.50</sub></td><td>0.33</td></tr><tr><td>HD 163296</td><td>D86</td><td>85.9</td><td>-0.94<sup>+0.21</sup><sub>-0.50</sub></td><td>0.15</td></tr><tr><td>AS 209</td><td>D61</td><td>61.7</td><td>-0.92<sup>+0.21</sup><sub>-0.50</sub></td><td>0.21</td></tr><tr><td>AS 209</td><td>D97</td><td>96.8</td><td>-0.06<sup>+0.21</sup><sub>-0.50</sub></td><td>0.42</td></tr></tbody></table>

**Notes.** Column (1): target name. Column (2): gap designation from Huang et al. (2018b). Column (3): estimated orbital semimajor axis, based on the derived distance (see Table 1) and $r_{gap}$ (see Table 3; note that these are not perfect matches for the designations in Column (2) owing to both the distance adjustments based on Gaia EDR3 and the modified approach for measuring gap centers). Column (4): estimated masses and uncertainties, based on the Zhang et al. (2018) approach and assuming a viscosity coefficient $\alpha = 0.001$ and a maximum particle size of 1 mm (see text). Column (5): relative gap width, as defined by Zhang et al. (2018) (their Equation (21)).

<sup>a</sup> The Sz 129 D64 gap is too narrow to reliably measure a planet mass using the Zhang et al. (2018) method.

(e.g., Pinilla et al. 2012, 2015). Contrary to our simple expectations, this might imply that CPDs in transition disk cavities will only emit weak millimeter continuum radiation despite their higher temperatures and (potentially) higher gas masses.

## 6.3. Future Directions

In any case, clear and meaningful constraints on the millimeter emission from CPDs for most of the DSHARP sample will require more sensitive measurements. This will likely mean a shift in strategy to target ALMA observations at higher frequencies, where the CPD is brighter (e.g., as advocated by Szulágyi et al. 2018). That comes with practical obstacles (more noise and less time available in favorable observing conditions), but also physical challenges. The local circumstellar disk emission will also be brighter at these frequencies and could be considerably more complex. At higher frequencies, the continuum emission traces smaller particles that are better coupled to the gas, which often results in emission gaps that are narrower and shallower (e.g., Tsukagoshi et al. 2016; Huang et al. 2018a, 2020; Carrasco-González et al. 2019; Long et al. 2020; Macías et al. 2021). Emission from such particles in the gaps would have higher optical depths and could contribute to circumstellar extinction of CPD signals. The search methodology advocated here would be beneficial with such data, although development of a more sophisticated recovery algorithm (rather than the simplistic peak/outlier identification that we have adopted) would likely be necessary.

It is especially striking that the disks that are predicted to host the most massive planets exhibit the most prominent

<figure><img src="image_15.png" /></figure>

**Figure 9.** The 1.3 mm continuum luminosity (flux normalized to a common distance of 150 pc) for disk material as a function of the stellar or planetary host mass; adapted from Andrews (2020). Upper limits are marked with downward-directed arrows. When 1.3 mm data were unavailable, the 0.8–0.9 mm measurements were scaled by the mean continuum spectrum $\nu^{2.2}$ (Andrews 2020; Tazzari et al. 2021). The black dashed line is the mean correlation derived from nearby circumstellar disks (in gray; Andrews et al. 2013; Ansdell et al. 2016; Pascucci et al. 2016; Barenfeld et al. 2017; Akeson et al. 2019; Cieza et al. 2019; Williams et al. 2019). The data points for PMCs are taken from MacGregor et al. (2017), Pineda et al. (2019), and Wu et al. (2020); the 2M1207 data from Ricci et al. (2017); the OTS 44 data from Bayo et al. (2017); and the PDS 70 c data from Isella et al. (2019). As a point of reference, the dotted gray curves show fully optically thick CPD models for the median properties of the DSHARP targets, $R_{cpd}/R_H = 0.3$, and the default (lower curve) and high (upper curve) accretion rate scenarios described above (assuming no scattering).

nonaxisymmetric residuals (HD 143006, SR 4, and HD 163296, according to Table 5; also Elias 24 at lower levels, which is made more intriguing by the tentative companion found by Jorquera et al. 2021). There is compelling indirect evidence for a perturber in the HD 143006 D22 gap (e.g., Ballabio et al. 2021), including a warp identified through extreme scattered light shadows (Benisty et al. 2018), the CO gas kinematics and continuum morphology at smaller radii (Pérez et al. 2018), and the faint $m = 1$ spiral identified here. As we noted above, these more massive companions nominally exhibit brighter millimeter continuum emission from their (warmer) CPDs, but their robust detection is unfortunately much more difficult owing to these asymmetries. More than raw sensitivity, this sort of nonaxisymmetric confusion limit could end up being the key obstacle in “blind” searches for CPDs like the effort presented here. It is possible that detectable CPD emission is present in these systems with the data used here, but we are incapable of differentiating it from asymmetric residuals that are not captured in the modeling.

In future work that pushes to improved sensitivity, it will be important to develop more flexible models of the circumstellar disk emission to help track down these faint CPDs. In some scenarios, particularly when the search region is confined by prior information, direct visibility modeling will outperform other approaches. Perhaps the biggest advance in hunting for millimeter continuum CPD emission will be forthcoming direct imaging planet detections that enable such targeted searches
and more cohesive explorations of the connections between planets, CPDs, and the disk substructures they create.

## 7. Summary

We used the high-resolution ALMA 1.25 mm continuum observations from the DSHARP survey to search for faint emission from circumplanetary material in the narrow gaps of circumstellar disks. Our key findings from this effort are summarized as follows:

1. We developed a prescription to mitigate contamination from the local circumstellar disk material, using the frank software package (Jennings et al. 2020) and a methodology for statistically quantifying the sensitivity to point-like CPD emission using injection-recovery experiments.

2. We found a few examples of pronounced asymmetric residuals in the target disks. The most interesting case is a faint, one-armed spiral that traverses across all of the axisymmetric substructures in the HD 143006 disk, possibly driven by a companion in the innermost (D22) gap.

3. There are a few peak residuals in these gaps that are marginal CPD candidates; deeper observations (preferably at better resolution) would be required to establish confidence that they are not merely local noise peaks. Upper limits on any CPD flux densities are 50–70 µJy in most cases, rising to 110 µJy in the few targets with lingering nonaxisymmetric features within their gaps.

4. If the gaps in these DSHARP disks are opened by giant planets with masses comparable to Jupiter, these constraints correspond to CPD (dust) mass upper limits of ~0.001–0.2 $M_⊕$. Alternatively, if the planet masses are much lower (as is the prediction for some targets based on the hydrodynamics simulations of Zhang et al. 2018), then considerably deeper observations may be required in future (sub)millimeter continuum CPD searches. Hopefully, those will be guided by direct imaging detections of the young planet hosts.

We are especially grateful to Kees Dullemond and Xue-Ning Bai for their helpful comments on the draft manuscript, Ryan Loomis for his technical advice and kindly providing access to CASA scripts, and Holly Thomas for her assistance facilitating computational efforts with the resources of the Smithsonian Radio Telescope Data Center. S.A. and W.E. acknowledge support from the National Aeronautics and Space Administration under grant No. 17-XRP17_2-0012 issued through the Exoplanets Research Program. The analysis in this article uses the following ALMA data sets:

1. ADS/JAO.ALMA #2016.1.00484.L

2. ADS/JAO.ALMA #2013.1.00226.S

3. ADS/JAO.ALMA #2013.1.00366.S

4. ADS/JAO.ALMA #2013.1.00498.S

5. ADS/JAO.ALMA #2013.1.00601.S

6. ADS/JAO.ALMA #2015.1.00486.S

7. ADS/JAO.ALMA #2015.1.00964.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 (JAO) is operated by ESO, AUI/NRAO, and NAOJ. This work used data from the European Space

Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Facility: ALMA.

Software: CASA (McMullin et al. 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), astropy (Astropy Collaboration et al. 2018), frank (Jennings et al. 2020), vis_sample (https://github.com/AstroChem/vis_sample).

## Appendix A
### Notes on Asymmetric Residuals

The modeling procedure outlined in Section 3 generally performs well in removing axisymmetric circumstellar disk emission that might contaminate the search for CPD emission. Understanding the patterned morphologies of lingering residuals can be useful for revising the (fixed) geometric parameters of the disk and interpreting other potential suboptimal assumptions made in the modeling. Using simple empirical models, we offer some brief guidance here to help identify the origins of such “artificial” residual features in future efforts (see also a similar discussion by Jennings et al. 2021).

We generated a series of synthetic data sets for this task. In each case, we assumed a fixed radial brightness temperature profile (in the Rayleigh–Jeans limit)

$$T_b(r) = 15 \left( \frac{r}{70 \text{ mas}} \right)^{-0.5} \text{K} \quad (A1)$$

out to a radius of 0.65, beyond which it decreases like $r^{-5}$. Two narrow gaps were imposed on this profile, from 100–120 mas and 470–540 mas, where the base profile above was multiplied by 0.01 and 0.05, respectively. This emission distribution in disk-frame polar coordinates $(r, \varphi)$ was used to make a synthetic image on a sky-frame coordinate system,

$$\begin{aligned} x_s &= r \sin \varphi \sin \text{PA} + (r \cos \varphi \cos i - z_s \sin i) \cos \text{PA} \\ y_s &= r \sin \varphi \cos \text{PA} - (r \cos \varphi \cos i - z_s \sin i) \sin \text{PA}, \quad (A2) \end{aligned}$$

where $z_s$ is the vertical height of the emission surface. In all of the analysis in the main text, and unless otherwise specified here, we have assumed $z_s = 0$. The center of the geometry specified in Equation (A2) could be offset from the image center by $(\Delta x, \Delta y)$, with positive values corresponding to shifts to the E and N, respectively. A synthetic visibility data set was generated from the Fourier transform of each image, sampled at the same spatial frequencies as the observations for the GW Lup disk, using the vis_sample software package.<fnref n="19" />

These synthetic visibilities were then modeled with frank following the same procedure outlined in Section 3. We assumed in this modeling the default geometric parameters $(\Delta x, \Delta y) = (0 \text{ mas}, 0 \text{ mas})$, $i = 35^{\circ}$, PA = 110$^{\circ}$, and $z_s = 0 \text{ mas}$ (the latter a requirement for frank). However, the synthetic visibilities were generated for a sequence of deviations from these default values. Figure 10 shows the imaged residual visibilities (as in Figure 3) to illustrate the

<footnote n="19">¹⁹ https://github.com/AstroChem/vis_sample</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 10.** Residual S/N images constructed from *frank* modeling with the default geometric parameters in the text, where the data were generated with the deviations from the default parameters noted in the upper left corner of each panel. All annotations are the same as in Figure 3. The “noise” in this case is set to the rms value for the GW Lup disk (Table 1).</figcaption></figure>

cases where the offsets, inclination, and/or PA have been slightly misassigned in the modeling with respect to the data. These residual images demonstrate that even small offsets produce strong ± residuals around the disk center, inclination mismatches are seen most prominently near “edges” (of the outer gap in this case) as quadrupolar ± residuals along the major and minor axes, and an incorrect PA assignment shows a similar behavior but rotated off-axis. Combinations of these features are somewhat more difficult to disentangle, although the asymmetric residual morphologies generated when an offset is misspecified are clear. When optimizing the geometry, a decision on the disk center can usually be made first and fixed before the projection angles are explored in more detail.

<figure><img src="image_17.png" /><figcaption>**Figure 11.** Residual S/N images as in Figure 10, but in this case associated with elevated emission surfaces.</figcaption></figure>

Figure 11 shows the analogous residual behavior when considering elevated emission surfaces, characterized with constant aspect ratio $z_s/r$, as might be expected in cases where continuum optical depths are high. Indeed, it was behavior like the top panels in some initial modeling exploration that led us to exclude some of the more highly inclined disk targets in the DSHARP sample—where such effects are most prevalent—from the analysis presented here. Modeling an elevated emission surface with a model that presumes an intrinsically flat morphology results in a pronounced, symmetric residual pattern along the minor axis. This behavior can be misinterpreted with other methodologies as a spatial offset in the minor-axis direction, but a careful examination of the residuals near the disk center can help distinguish the difference (see the middle right panel). Mixing surface and other geometric effects can create complicated residual patterns (see, e.g., the HD 163296 residuals).

# Appendix B

## Model Treatment of Confined Azimuthal Asymmetries

The HD 143006 and HD 163296 disks have pronounced, but spatially confined, “arc”-like azimuthal features in their emission distributions that present a challenge for the standard axisymmetric modeling methodology outlined in Section 3 (Isella et al. 2018; Pérez et al. 2018). Specifically, these asymmetries are sufficiently bright that ignoring them leads *frank* to derive axisymmetric models that overpredict the emission at comparable radii, resulting in a pronounced negative residual at azimuths that lie well away from the asymmetry. A demonstration of that effect is clear in the residual S/N images shown in the top right panels of Figure 12.
<figure><img src="image_18.png" /><figcaption>**Figure 12.** Graphical overviews of the workflow for removing the confined azimuthal asymmetries (the “arc” features confined to the regions outlined in each image) for the HD 143006 (top) and HD 163296 (bottom) disks. The upper left panels show the effects of modeling the emission without treating the asymmetries, in both the images (as in Figures 1 and 2) and visibilities (as in Figure 14). The bottom left panels illustrate how the asymmetry models are constructed. The bottom right panels demonstrate the results of modeling the revised data, where the asymmetries have been removed. The top right panels directly compare the residual S/N maps (as in Figure 3) for the uncorrected (with asymmetry) and revised (asymmetry removed) cases.</figcaption></figure>
We designed a workaround that first revises the data visibilities by removing a simple model for the asymmetry before performing the frank modeling described in Section 3. This asymmetry model was constructed first by spatially isolating the feature in the CLEAN model image, setting the CLEAN components outside a specific area associated with the asymmetry to zero. The area of interest was selected manually by specifying radial and azimuthal boundaries in the disk plane. For HD 143006, this corresponded to the annular arc spanning radii from 0.′37 to 0.′60 and azimuths from 90° to 142°, following the Huang et al. (2018b) azimuth convention (where 90° coincides with the major axis and azimuths increase clockwise on the sky). For HD 163296, the radial and azimuthal boundaries were 0.′48–0.′60 and 50°–150°, respectively. Next, the mean radial profile in the CLEAN model constructed from outside that region was subtracted from the asymmetry model image, leaving only the asymmetric contribution. These steps are illustrated in the bottom left parts of Figure 12.

Next, the Fourier transform of the asymmetry model image was sampled at the observed spatial frequencies and then subtracted from the original data visibilities. Finally, the resulting revised data visibilities were then modeled as described in Section 3. The entire process was iterated to settle on appropriate geometric parameters before finalizing the modeling outcomes.

The full workflows for treating these special cases are illustrated together in Figure 12. A direct comparison of the residual S/N images in the top right sections of these composite figures demonstrates a notable improvement in the fit quality achieved by first removing the confined azimuthal asymmetries. However, it is interesting to see that these two cases with the most pronounced nonaxisymmetric features still end up exhibiting considerable low-level asymmetric structure that is not captured by the modeling (see especially Figure 4).

Presumably that is a distinctive, albeit subtle, clue to the origins of the substructures in these cases.

## Appendix C: Model Radial Profiles

Figure 13 shows the model brightness profiles derived with frank following the procedure outlined in Section 3. The gray bands mark the search zones for CPD emission, as described in Section 4, and the dotted lines denote the gap centers identified by Huang et al. (2018b) (using a different approach). There is good agreement with those results, though the frank modeling hints at some annular features that are unresolved in the nominal DSHARP images (e.g., the gap at ~30 mas in the SR 4 disk first identified by Jennings et al. 2021).

These profiles were used to crudely measure the centers, widths, and depths of the gaps of interest, as noted in Section 3. These estimates were made by visual comparison to a simplistic model of a background (local) power-law profile with a Gaussian depletion,

$$ T_b = T_0 \frac{(r/0.′1)^{-q}}{(1 + \Gamma)}, \quad \text{where} \\ \Gamma = (\delta_{\text{gap}} - 1)\exp\left[-\frac{(r - r_{\text{gap}})^2}{2\sigma_{\text{gap}}^2}\right]. \quad (\text{C1}) $$

The adopted power-law ($T_0$, $q$) and gap ($r_{gap}$, $\sigma_{gap}$, $\delta_{gap}$) parameters are compiled in Table 3. In the context of Equation (C1), the depletion parameter $\delta_{gap}$ is a multiplicative scale (amplitude); e.g., $\delta_{gap} = 10$ implies relative depletion by an order of magnitude at $r_{gap}$.

Finally, to demonstrate the fit quality for the circumstellar disk emission in the native (Fourier) domain, Figure 14 compares the deprojected, azimuthally averaged visibilities with the corresponding frank models.
<figure><img src="image_19.png" /><figcaption>**Figure 13.** The brightness temperature radial profiles for the circumstellar disk emission derived in Section 3 using the `frank` package (assuming the Rayleigh-Jeans approximation). Gray bands mark the gap annuli ($r_{gap} \pm \sigma_{gap}$); dotted lines denote the gap centers measured in the image plane by Huang et al. (2018b). As a reminder, the gap feature at ~0.8 in the AS 209 disk (D97) was treated as two gaps separated by a faint ring by Huang et al. (2018b); see the note in Table 3.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 14. The shifted, deprojected (see Table 2 for geometric parameters), and azimuthally averaged (into 15 kλ wide bins) visibilities as a function of baseline length for the data (gray) and the optimized frank model (red). Note that by definition the imaginary components of the models are zero for all spatial frequencies.</figcaption></figure>

## Appendix D

### The CPD Model in the Fourier Domain

To quantify the sensitivity to CPD emission in the DSHARP data, we adopted the approach described in Section 4 that characterizes the recovery rate of injected CPD signals. In this framework, the mock CPD signal is computed in the Fourier domain and added to the observed (complex) “data” visibilities, $V_d(u, v)$, where $u, v$ are the Fourier spatial frequency coordinates (in wavelength units). Those composite visibilities, $V_d + V_{cpd}$, are modeled with frank, and the imaged residual visibilities are then searched for remnant CPD emission and compared with the known input parameters.

The CPD emission model is an offset point source with flux density $F_{cpd}$ at ($r_{cpd}$, $\varphi_{cpd}$) in the disk frame,

$$I_{cpd} = F_{cpd} \delta(r - r_{cpd}, \varphi - \varphi_{cpd}), \quad (D1)$$

where $\delta$ is the Dirac $\delta$-function. The Fourier transform of a point source at the origin is a constant (DC offset), in this case just $F_{cpd}$. But the offset position introduces an analytic oscillatory behavior. The sky-frame coordinates of the mock CPD ($x_s, y_s$) are given by Equation (A2) for $z_s = 0$. For a given set of geometric parameters (Table 2), the CPD visibilities can be expressed as

$$V_{cpd}(u, v) = F_{cpd} e^{2\pi i[u(x_s + \Delta x) + v(y_s + \Delta y)]}, \quad (D2)$$
where sky-projected terms are expressed in radians and in this case i is the imaginary unit (not inclination).

## Appendix E
CPD Mass Constraints

The framework we adopted to convert the derived flux constraints into limits on CPD masses follows closely the approach outlined by Isella et al. (2014, 2019). For convenience, we review the details here.

Physical models for the CPD continuum emission require assignments of the temperature and density structure. We treat the CPD as geometrically flat and axisymmetric. The CPD thermal structure can be approximated with the contributions of three mechanisms—irradiation by the planet, irradiation by the star and local disk (around the gap), and accretion heating—such that

$$T^4 = T_{\text{irr, p}}^4 + T_{\text{irr, *}}^4 + T_{\text{acc.}}^4 \quad (\text{E1})$$

There are a lot of parameters and assumptions hidden in Equation (E1); each will be clarified below.

The heating contribution of irradiation by the “host” planet can be approximated as

$$T_{\text{irr, p}}^4 = \frac{0.1 L_p}{4\pi\sigma r^2}, \quad (\text{E2})$$

where $L_p$ denotes the planet luminosity, the factor 0.1 corresponds to the fraction of energy absorbed by a CPD with a vertical pressure scale height ~10% of the radius, $\sigma$ is the Stefan-Boltzmann constant, and $r$ is the radial coordinate in the CPD frame. Since no planets associated with the disk gaps of interest here have been detected, we relied on theoretical models to estimate $L_p$ for any given $M_p$ (a free parameter). The curve in the top panel of Figure 15 shows the adopted relationship, derived from an interpolation of planetary evolution models (at a fixed age of 1 Myr) calculated by Linder et al. (2019) (solar-metallicity models from the petitCODE atmosphere modeling code) and Spiegel & Burrows (2012) (the “hot start” models shown in their Figure 5).

For models in “low” accretion states (see below), this planetary irradiation heating dominates at least the inner part of the CPD. The model planet luminosities decrease by a factor of 2–3 from 1 to 3 Myr, which amounts to only 20%–30% in $T_{\text{irr, p}}$. Since most of the emission comes from larger radii where other heating terms contribute, the corresponding decrease in the continuum emission (and thereby CPD mass estimates) is only ~5%–15%. For more massive planets, if we instead adopted the “cold start” models from Spiegel & Burrows (2012), the decrease in $L_p$ is more like a factor of 100. While that decreases $T_{\text{irr, p}}$ by a factor of ~3, the emission decrease is much more muted (~30%–40%) because stellar/disk irradiation heating then dominates (particularly at larger $r$, where most of the emission originates).

The irradiation heating from the central star and the local (circumstellar) disk is roughly constant in $r$, with

$$T_{\text{irr, *}}^4 = \frac{\phi L_*}{8\pi\sigma r_{\text{gap}}^2}, \quad (\text{E3})$$

where $L_*$ is the luminosity of the central star, $r_{\text{gap}}$ is the planet/CPD location within the host disk (see Table 3), and $\phi$ is the flaring angle of the host disk. We adopted the $L_*$ values

<figure><img src="image_21.png" /><figcaption>Figure 15. The adopted variations of planet luminosity and radius as a function of mass used in CPD model calculations (blue dotted curves), compared with predictions from detailed planetary evolution models at ages of 1 and 3 Myr (large and small symbols, respectively). The orange circles show the solar-metallicity predictions from Linder et al. (2019), using the petitCODE atmosphere modeling code and assuming cloud-free atmospheres. The red crosses and magenta plus signs denote the “hot start” and “cold start” models from Spiegel & Burrows (2012, their Figure 5). The adopted curves are derived from spline interpolations of the Linder et al. (2019) 1 Myr models and the Spiegel & Burrows (2012) 1 Myr “hot start” models at each $M_p$.</figcaption></figure>

cataloged by Andrews et al. (2018) and the $\phi$ values used by Huang et al. (2018b) (0.05 for the RU Lup disk, 0.02 otherwise). If we generously allow that both $\phi$ and $L_*$ are uncertain by a factor of 2, that corresponds to a ~40% ambiguity in $T_{\text{irr, *}}$. For lower $M_p$ where $T_{\text{irr, *}}$ dominates, this propagates almost directly into the CPD continuum luminosity; at larger $M_p$, the effect is considerably smaller (5%–10%).

The viscous heating term can be approximated as

$$T_{\text{acc}}^4 = \frac{3GM_p \dot{M}_p}{8\pi\sigma r^3} \left( 1 - \sqrt{\frac{R_p}{r}} \right) \quad (\text{E4})$$

(see D’Alessio et al. 1998), where $\dot{M}_p$ is the planet accretion rate, $R_p$ is the planet radius (presumed equivalent to the inner edge of the CPD), and $G$ is the gravitational constant. The adopted $R_p$ also depends on $M_p$ and an approximation of planetary evolution models, as shown in the bottom panel of Figure 15. We made the approximation that $\dot{M}_p = \eta M_p / t_p$, where $t_p$ is a characteristic timescale (assumed to be 1 Myr) and $\eta$ is an efficiency factor. The “low” (default) accretion state has $\eta = 1$, and a “high” state has $\eta = 100$. In the low state, irradiation heating dominates for all $M_p$; swapping $t_p$ to 3 Myr makes essentially no difference in the CPD fluxes. But in the high state, viscous heating dominates in many cases (at least for $M_p \gtrsim 0.1 M_{\text{Jup}}$), and the output
continuum fluxes are ~2× higher. There is considerable ambiguity associated with this heating term, since we do not yet understand the details of the CPD accretion process. We consider the low- and high-state cases reasonable boundary conditions on the associated uncertainties (a factor of ~two in the continuum emission levels).

The CPD density structure was described with a radial power-law for the surface densities,

$$ \Sigma = \Sigma_0 \left( \frac{r}{R_{\text{cpd}}} \right)^{-\gamma} \qquad (\text{E5}) $$

defined for $R_p \le r \le R_{\text{cpd}}$, where

$$ \Sigma_0 = (2 - \gamma) \frac{M_{\text{cpd}}}{2\pi R_{\text{cpd}}^{\gamma}} (R_{\text{cpd}}^{2-\gamma} - R_p^{2-\gamma})^{-1}. \qquad (\text{E6}) $$

Here, $M_{\text{cpd}}$ and $R_{\text{cpd}}$ are the CPD mass and outer radius, respectively. We set $\gamma = 0.75$, comparable to the analytical models of Canup & Ward (2002). The flux differences associated with different $\gamma$ (e.g., from 0 to 1.5) are relatively small (~20%) but depend in detail on the optical depths and heating terms for a given model.

The continuum flux from this CPD model is

$$ F_{\text{cpd}} = \frac{2\pi\mu}{d^2} \int B_{\nu}(T) [1 - e^{-\tau_{\nu}/\mu} + \omega_{\nu}\mathcal{F}] r dr \qquad (\text{E7}) $$

(Miyake & Nakagawa 1993; Sierra et al. 2019; Zhu et al. 2019), where $\mu = \cos i$ is the direction cosine (with $i$ assumed equivalent to the inclination in Table 2), $d$ is the distance to Earth (Table 1), and $B_{\nu}(T)$ is the Planck function. The dust grain emission properties are characterized with an absorption opacity $\kappa_{\nu}$ and albedo $\omega_{\nu}$. Then, the optical depth is $\tau_{\nu} = \kappa_{\nu}(1 - \omega_{\nu})^{-1}\Sigma$. The scattering correction term in Equation (E7) is

$$ \mathcal{F} = \frac{\left[ \frac{1 - e^{-\sqrt{3} \epsilon_{\nu} \tau_{\nu} - \nu \epsilon_{\nu} / \mu}}{\sqrt{3} \epsilon_{\nu} \mu + 1} + \frac{e^{-\nu \epsilon_{\nu} / \mu} - e^{-\sqrt{3} \epsilon_{\nu} \tau_{\nu}}}{\sqrt{3} \epsilon_{\nu} \mu - 1} \right]}{(\epsilon_{\nu} - 1)e^{-\sqrt{3} \epsilon_{\nu} \tau_{\nu}} - (\epsilon_{\nu} + 1)}, \qquad (\text{E8}) $$

where $\epsilon_{\nu} = \sqrt{1 - \omega_{\nu}}$. The main text explores the effect of albedo on the CPD flux (and mass) estimates for extreme boundary conditions. Of course, the ambiguities associated with the optical properties of the particles that emit the millimeter continuum studied here are relevant and are expected to dominate the $M_{\text{cpd}}$ uncertainties; more details on this general problem are discussed elsewhere (e.g., Birnstiel et al. 2018; Andrews 2020).

## ORCID iDs

Sean M. Andrews https://orcid.org/0000-0003-2253-2270

Shangjia Zhang https://orcid.org/0000-0002-8537-9114

Jane Huang https://orcid.org/0000-0001-6947-6072

Myriam Benisty https://orcid.org/0000-0002-7695-7605

Nicolás T. Kurtovic https://orcid.org/0000-0002-2358-4796

David J. Wilner https://orcid.org/0000-0003-1526-7587

Zhaohuan Zhu https://orcid.org/0000-0003-3616-6822

John M. Carpenter https://orcid.org/0000-0003-2251-0602

Laura M. Pérez https://orcid.org/0000-0002-1199-9564

Richard Teague https://orcid.org/0000-0003-1534-5186

Andrea Isella https://orcid.org/0000-0001-8061-2207

Luca Ricci https://orcid.org/0000-0001-8123-2943

## References

Akeson, R. L., Jensen, E. L. N., Carpenter, J., et al. 2019, ApJ, 872, 158

Andrews, S. M. 2020, ARA&A, 58, 483

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41

Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129

Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJL, 820, L40

Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44

Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 393, 49

Bae, J., Zhu, Z., & Hartmann, L. 2017, ApJ, 850, 201

Ballabio, G., Nealon, R., Alexander, R. D., et al. 2021, MNRAS, 504, 888

Barenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2017, ApJ, 827, 142

Batygin, K., & Morbidelli, A. 2020, ApJ, 894, 143

Bayo, A., Joergens, V., Liu, Y., et al. 2017, ApJ, 841, L11

Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924

Benisty, M., Bae, J., Facchini, S., et al. 2021, ApJL, submitted

Benisty, M., Juhász, A., Facchini, S., et al. 2018, A&A, 619, A171

Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJL, 869, L45

Bowler, B. P. 2016, PASP, 128, 102001

Bowler, B. P., Andrews, S. M., Kraus, A. L., et al. 2015, ApJL, 805, L17

Canup, R. M., & Ward, W. R. 2002, AJ, 124, 3404

Carrasco-González, C., Sierra, A., Flock, M., et al. 2019, ApJ, 883, 71

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

Cieza, L. A., González-Ruilova, C., Hales, A. S., et al. 2021, MNRAS, 501, 2934

Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698

Clarke, C. J., Tazzari, M., Juhasz, A., et al. 2018, ApJL, 866, L6

Cugno, G., Quanz, S. P., Hunziker, S., et al. 2019, A&A, 622, A156

Czekala, I., Loomis, R. A., Teague, R., et al. 2021, ApJS, submitted

D'Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411

Disk Dynamics Collaboration, Armitage, P. J., Bae, J., et al. 2020, arXiv:2009.04345

Doi, K., & Kataoka, A. 2021, ApJ, 912, 164

Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJL, 809, L5

Drăzkowska, J., & Szulágyi, J. 2018, ApJ, 866, 142

Eisner, J. A. 2015, ApJL, 803, L4

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Garufi, A., Benisty, M., Pinilla, P., et al. 2018, A&A, 620, A94

Gaudi, B. S. 2012, ARA&A, 50, 411

Guzmán, V. V., Huang, J., Andrews, S. M., et al. 2018, ApJL, 869, L48

Haffert, S. N., Bohn, A. J., de Boer, J., et al. 2019, NatAs, 3, 749

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hashimoto, J., Aoyama, Y., Konishi, M., et al. 2020, AJ, 159, 222

Huang, J., Andrews, S. M., Cleeves, L. I., et al. 2018a, ApJ, 852, 122

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018b, ApJL, 869, L42

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2020, ApJ, 891, 48

Huang, J., Andrews, S. M., Pérez, L. M., et al. 2018c, ApJL, 869, L43

Hunter, J. D. 2007, Computing in Science Engineering, 9, 90

Isella, A., Benisty, M., Teague, R., et al. 2019, ApJL, 879, L25

Isella, A., Chandler, C. J., Carpenter, J. M., Pérez, L. M., & Ricci, L. 2014, ApJ, 788, 129

Isella, A., Guidi, G., Testi, L., et al. 2016, PhRvL, 117, 251101

Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJL, 869, L49

Jennings, J., Booth, R. A., Tazzari, M., Clarke, C. J., & Rosotti, G. P. 2021, arXiv:2103.02392

Jennings, J., Booth, R. A., Tazzari, M., Rosotti, G. P., & Clarke, C. J. 2020, MNRAS, 495, 3209

Jin, S., Li, S., Isella, A., Li, H., & Ji, J. 2016, ApJ, 818, 76

Jorquera, S., Pérez, L. M., Chauvin, G., et al. 2021, AJ, 161, 146

Jorsater, S., & van Moorsel, G. A. 1995, AJ, 110, 2037

Kanagawa, K. D., Tanaka, H., Muto, T., Tanigawa, T., & Takeuchi, T. 2015, MNRAS, 448, 994

Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44

Keppler, M., Teague, R., Bae, J., et al. 2019, A&A, 625, A118

Kurtovic, N. T., Pérez, L. M., Benisty, M., et al. 2018, ApJL, 869, L44

Law, C. J., Loomis, R. A., Teague, R., et al. 2021, ApJS, submitted

Linder, E. F., Mordasini, C., Mollière, P., et al. 2019, A&A, 623, A85

Lodato, G., Dipierro, G., Ragusa, E., et al. 2019, MNRAS, 486, 453

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17
Long, F., Pinilla, P., Herczeg, G. J., et al. 2020, ApJ, 898, 36

MacGregor, M. A., Wilner, D. J., Czekala, I., et al. 2017, ApJ, 835, 17

Macías, E., Guerra-Alvarado, O., Carrasco-González, C., et al. 2021, A&A, 648, A33

Martin, R. G., & Lubow, S. H. 2011, MNRAS, 413, 1447

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Miyake, K., & Nakagawa, Y. 1993, Icar, 106, 20

Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13

Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125

Penny, M. T., Gaudi, B. S., Kerins, E., et al. 2019, ApJS, 241, 3

Pérez, L. M., Benisty, M., Andrews, S. M., et al. 2018, ApJL, 869, L50

Pérez, S., Marino, S., Casassus, S., et al. 2019, MNRAS, 488, 1005

Pineda, J. E., Szulágyi, J., Quanz, S. P., et al. 2019, ApJ, 871, 48

Pinilla, P., Benisty, M., & Birnstiel, T. 2012, A&A, 545, A81

Pinilla, P., de Juan Ovelar, M., Ataiee, S., et al. 2015, A&A, 573, A9

Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJL, 890, L9

Quillen, A. C., & Trilling, D. E. 1998, ApJ, 508, 707

Reid, M. J., Schneps, M. H., Moran, J. M., et al. 1988, ApJ, 330, 809

Ricci, L., Cazzoletti, P., Czekala, I., et al. 2017, AJ, 154, 24

Salyk, C., Blake, G. A., Boogert, A. C. A., & Brown, J. M. 2009, ApJ, 699, 330

Shibaike, Y., Okuzumi, S., Sasaki, T., & Ida, S. 2017, ApJ, 846, 81

Sierra, A., Lizano, S., Macías, E., et al. 2019, ApJ, 876, 7

Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174

Stamatellos, D., & Herczeg, G. J. 2015, MNRAS, 449, 3432

Szulágyi, J., Binkert, F., & Surville, C. 2021, arXiv:2103.12128

Szulágyi, J., Dullemond, C. P., Pohl, A., & Quanz, S. P. 2019, MNRAS, 487, 1248

Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853

Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65

Szulágyi, J., Plas, G. V. D., Meyer, M. R., et al. 2018, MNRAS, 473, 3573

Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47

Tazzari, M., Clarke, C. J., Testi, L., et al. 2021, MNRAS, in press

Tsukagoshi, T., Nomura, H., Muto, T., et al. 2016, ApJL, 829, L35

van der Marel, N., Dong, R., di Francesco, J., Williams, J. P., & Tobin, J. 2019, ApJ, 872, 112

Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Wang, J. J., Vigan, A., Lacour, S., et al. 2021, AJ, 161, 148

Williams, J. P., Cieza, L., Hales, A., et al. 2019, ApJL, 875, L9

Wu, Y.-L., Bowler, B. P., Sheehan, P. D., et al. 2020, AJ, 159, 229

Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJL, 869, L47

Zhu, Z. 2015, ApJ, 799, 16

Zhu, Z., Andrews, S. M., & Isella, A. 2018, MNRAS, 479, 1850

Zhu, Z., Nelson, R. P., Hartmann, L., Espaillat, C., & Calvet, N. 2011, ApJ, 729, 47

Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJL, 877, L18

Zurlo, A., Cugno, G., Montesinos, M., et al. 2020, A&A, 633, A119