# A Highly Settled Disk around Oph163131

M. Villenave, K. R. Stapelfeldt, G. Duchêne, F. Ménard, M. Lambrechts, A.

Sierra, C. Flores, W. R. F. Dent, S. Wolff, Á. Ribas, et al.

► To cite this version:

M. Villenave, K. R. Stapelfeldt, G. Duchêne, F. Ménard, M. Lambrechts, et al. A Highly Settled Disk around Oph163131. The Astrophysical Journal, 2022, 930, 10.3847/1538-4357/ac5fae. insu-03705352

HAL Id: insu-03705352

https://insu.hal.science/insu-03705352

Submitted on 27 Jun 2022

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_1.png" /></figure>
OPEN ACCESS

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

CrossMark

# A Highly Settled Disk around Oph163131

M. Villenave¹, K. R. Stapelfeldt¹, G. Duchêne²,³, F. Ménard³, M. Lambrechts⁴,⁵, A. Sierra⁶, C. Flores⁷, W. R. F. Dent⁸, S. Wolff⁹, Á. Ribas¹⁰, M. Benisty³,¹¹, N. Cuello³, and C. Pinte³,¹²

¹ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA; marion.f.villenave@jpl.nasa.gov

² Astronomy Department, University of California, Berkeley, CA 94720, USA

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

⁴ Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, SE-22100 Lund, Sweden

⁵ Centre for Star and Planet Formation, Globe Institute, University of Copenhagen, Øster Voldgade 5-7, DK-1350 Copenhagen, Denmark

⁶ Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile

⁷ Institute of Astronomy, University of Hawaii, 640N Aohoku Place, Hilo, HI 96720, USA

⁸ Joint ALMA Observatory, Alonso de Córdova 3107, Vitacura 763-0355, Santiago, Chile

⁹ Steward Observatory, University of Arizona, Tucson, AZ 85721, USA

¹⁰ European Southern Observatory, Alonso de Córdova 3107, Vitacura, Casilla 19001, Santiago 19, Chile

¹¹ Unidad Mixta Internacional Franco-Chilena de Astronomía, CNRS, UMI 386 Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile

¹² School of Physics and Astronomy, Monash University, Clayton Vic 3800, Australia

Received 2021 November 16; revised 2022 March 10; accepted 2022 March 11; published 2022 April 28

## Abstract

High dust density in the midplane of protoplanetary disks is favorable for efficient grain growth and can allow fast formation of planetesimals and planets, before disks dissipate. Vertical settling and dust trapping in pressure maxima are two mechanisms allowing dust to concentrate in geometrically thin and high-density regions. In this work, we aim to study these mechanisms in the highly inclined protoplanetary disk SSTC2D J163131.2-242627 (Oph 163131, *i* ~ 84°). We present new high-angular-resolution continuum and ¹²CO ALMA observations of Oph 163131. The gas emission appears significantly more extended in the vertical and radial direction compared to the dust emission, consistent with vertical settling and possibly radial drift. In addition, the new continuum observations reveal two clear rings. The outer ring, located at ~100 au, is well-resolved in the observations, allowing us to put stringent constraints on the vertical extent of millimeter dust particles. We model the disk using radiative transfer and find that the scale height of millimeter-sized grains is 0.5 au or less at 100 au from the central star. This value is about one order of magnitude smaller than the scale height of smaller micron-sized dust grains constrained by previous modeling, which implies that efficient settling of the large grains is occurring in the disk. When adopting a parametric dust settling prescription, we find that the observations are consistent with a turbulent viscosity coefficient of about α ≲ 10⁻⁵ at 100 au. Finally, we find that the thin dust scale height measured in Oph 163131 is favorable for planetary growth by pebble accretion: a 10 M_E planet may grow within less than 10 Myr, even in orbits exceeding 50 au.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Planet formation (1241); Radiative transfer (1335); Dust continuum emission (412)*

*Supporting material: data behind figure*

## 1. Introduction

In recent years, a number of studies have presented direct or indirect detections of planets embedded in protoplanetary disks (e.g., Keppler et al. 2018; Pinte et al. 2018; Teague et al. 2018), indicating that planet formation is a fast process. In the core-accretion paradigm, micron-sized particles that are inherited from the interstellar medium need to grow quickly to form a planetary core, which can attract the surrounding gas before the disk has dissipated. The streaming instability (Youdin & Goodman 2005) is one of the currently favored mechanisms allowing a sufficiently rapid growth, from pebbles to planetesimals, in the disk phase. For this instability to develop, the dust-to-gas ratio needs to be high in the disk midplane, of the order of unity, and large pebble-sized particles need to be present. After this planetesimal formation stage, planetary embryos can continue

to form the cores of giant planets by accreting remaining inward-drifting pebbles (Lambrechts & Johansen 2012).

Various studies have shown that substructures, in particular rings and gaps, are ubiquitous in protoplanetary disks in the Class II phase (Huang et al. 2018; Long et al. 2018), and already present in some young Class I disks (ALMA Partnership et al. 2015; Segura-Cox et al. 2020). In many cases, they act as dust traps and grain growth is associated with rings (e.g., Carrasco-González et al. 2019; Macías et al. 2021; Sierra et al. 2021). On the other hand, the efficiency of dust vertical settling remains largely unconstrained. Similarly to radial drift (Weidenschilling 1977), this mechanism is a balance between the stellar gravity and the interaction of dust with gas. Large grains (e.g., millimeter sizes) are expected to settle most efficiently toward the midplane, while smaller grains (e.g., micron sizes) are predicted to remain well-coupled to the gas and are co-located with it. At the same time, because of radial drift, large grains are also predicted to drift toward the star.

Currently, the vertical scale height of the gas in protoplanetary disks around T Tauri stars has been constrained in a

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
number of studies to be about 10 au at a radius of 100 au (e.g., Burrows et al. 1996; Watson et al. 2007; Wolff et al. 2017). This was mainly done by modeling scattered-light images, which probe small dust grains assumed to be well-coupled to the gas. Some other studies also estimated near-infrared scattering surface extent (e.g., Ginski et al. 2016; Avenhaus et al. 2018) or the height of gas emission layers (e.g., Pinte et al. 2018; Law et al. 2021; Rich et al. 2021), using less model-dependent techniques. However, these latter techniques do not probe the pressure scale height of the disk; instead, they probe the scattering or emission surfaces, which may be several times higher than the physical scale height.

On the other hand, only a few studies have estimated the scale height of millimeter dust particles, which are most affected by vertical settling. This is in part because the majority of observed protoplanetary disks are moderate-inclination systems in which it is difficult to constrain the thin vertical extent of millimeter grains. For these systems, detailed modeling of gaps and rings is needed (e.g., Pinte et al. 2016; Doi & Kataoka 2021). On the other hand, disks seen edge-on offer the most favorable orientation to study their vertical extent. Villenave et al. (2020) presented the first high-angular-resolution (~0.′′1) millimeter survey of edge-on disks. By comparing the vertical and radial extent of three systems with radiative transfer models, they showed that the scale height of the millimeter dust particles is on the order of a few au at 100 au, significantly smaller than the typical gas scale height. These results indicate that efficient vertical settling is occurring in protoplanetary disks. Further studies are needed to increase the number of disks with known gas and millimeter dust vertical extent in order to better understand the efficiency of vertical settling.

In this paper, we focus on SSTC2D J163131.2-242627 (hereafter Oph 163131), a highly inclined protoplanetary disk located in the Ophiuchus star-forming region, at a distance of about 147 ± 3 pc (Ortiz-León et al. 2017, 2018). Oph 163131 was included in the survey of Villenave et al. (2020), and has recently been studied by two companion papers (Flores et al. 2021; Wolff et al. 2021), which presented ALMA observations at an angular resolution of ~0.′′2, as well as scattered-light HST and Keck images. Wolff et al. (2021) used radiative transfer to model both the 1.3 mm ALMA image and the 0.8 µm HST images with an extensive MCMC framework. They found that some degree of vertical settling is needed to reproduce both observations. On the other hand, Flores et al. (2021) analyzed the ALMA $^{12}$CO and $^{13}$CO maps to characterize the 2D temperature structure of the disk. They obtained a dynamical mass estimate of 1.2 ± 0.2 $M_\odot$ from the ALMA observations, and characterized the spectral type of the source to be K4, using optical and near-infrared spectroscopy.

In this work, we present new high-angular-resolution observations of Oph 163131 at 1.3 mm (resolution of ~0.′′02 or 3 au). The new images reveal a number of rings that were not detected with previous millimeter observations. We derive a detailed radiative transfer model of the new millimeter observations to characterize the physical structure of the disk and add constraints on vertical settling. In Section 2, we present the observations and data reduction. In Section 3, we present the main features seen in the images. In Section 4, we describe the disk model and present the results. We discuss the implication of our results in Section 5. Finally, our conclusions are presented in Section 6.

## 2. Observation and Data Reduction

We present new cycle 6 observations of Oph 163131 in band 6 (1.3 mm, Project: 2018.1.00958.S; PI: Villenave). The spectral setup was divided in three continuum spectral windows of rest frequency 229.0, 243.5, and 246.0 GHz, and a fourth spectral window including the $^{12}$CO $J=2-1$ transition at 230.538 GHz. The line spectral window has a native velocity resolution of 0.64 km s⁻¹. The data were obtained on June 8, 2019, with baselines ranging from 80 m to 16 km. The total observing time on source was about 2 hr and 30 minutes. The raw data were calibrated using the CASA pipeline version 5.4.0 (McMullin et al. 2007), and the rest of the data processing was done with CASA 5.6.1.

To produce the final images, we combined our cycle 6 observations with lower-angular-resolution observations from cycle 4 (Project 2016.1.00771.S; PI: Duchêne) that were published previously (Villenave et al. 2020; Flores et al. 2021; Wolff et al. 2021). We note that a shift of about 0.′′05 (~20% of the cycle 4 beam) was present between the cycle 4 and cycle 6 observations, which is roughly consistent with the astrometric accuracy of ALMA (see ALMA technical handbook<fnref n="13" />). Thus, before producing the combined image, we aligned both observations using the fixplanet CASA task (with the option fixuvw as True). To maximize the dynamic range of the final image, we performed phase self-calibration on the cycle 4 continuum observations as mentioned in Flores et al. (2021) and Villenave et al. (2020). Due to the limited signal-to-noise per beam, no self-calibration could be performed on the cycle 6 observations, nor on the combined cycle 4 and cycle 6 data. We then produced the continuum image using the CASA tclean task on the combined data set, with a Briggs weighting (robustness parameter of +0.5) and using the multiscale deconvolver, with scales of 0, 1, 5, and 10 times the beam FWHM. The resulting continuum beam size is 0.′′024 × 0.′′020 (~3.5 × 3 au) with a major axis at PA = −81°, and the resulting continuum rms is 9.3 µJy beam⁻¹.

After applying the continuum self-calibration solutions to all spectral windows, we subtracted the continuum emission using the uvcontsub task in CASA. Then, we derived the emission line maps from the calibrated visibilities using the tclean function. We use 0.7 km s⁻¹ velocity resolution and a Briggs weighting (robustness parameter of +0.5) to create the line images. Additionally, to increase the signal-to-noise of the $^{12}$CO observations, we applied a UV taper while generating the images. The combined resulting $^{12}$CO beam is 0.′′081 × 0.′′072 with a major axis at PA = 89°, and the average rms of the line is 0.8 mJy beam⁻¹ per channel. Finally, we generate moment 0 and 1 maps including only pixels above three times the rms.

## 3. Results

### 3.1. Continuum Emission

We present the continuum observations of Oph 163131 in Figure 1. The high angular resolution of the image reveals several rings even though the disk is highly inclined. We highlight the structures in the right panel of Figure 1. From outside in, we see two rings (ring 2 and ring 1) separated by a clear gap (gap 1), some emission inside the second ring (that we call region 2 in Figure 1), and an inner central emission

<footnote n="13">¹³ https://almascience.eso.org/documents-and-tools/cycle8/alma-technical-handbook</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Left: continuum image of Oph 163131. The beam size is indicated by an ellipse in the bottom left corner of the plot. North is up and east is left. Right: labeled zoomed-in continuum image of Oph 163131, scaled so that the rings are more visible (see Section 4.2 and Figure 8 for details on the identification of region 1). The central region is saturated (black), and pixels with less than $2\sigma$ appear in white. The beam size is indicated by an ellipse in the bottom left corner of the plot. The continuum image is available as the data behind the figure in FITS format.

*(The data used to create this figure are available.)*</figcaption></figure>

(region 1 and central point source; see Section 4.2). We note that the existence of the outer gap was hinted by a shoulder detected in previous observations (Villenave et al. 2020), but it was not resolved as with the current image.

To characterize the substructures, we fit the deprojected visibilities using the `frank` package (Jennings et al. 2020). We exported the visibilities using a modified version of `expor-t_uvtable` from Tazzari (2017). In our version, instead of using the average wavelength from all data, each baseline is normalized using the wavelength associated with each spectral window and channel. The visibilities from cycles 6 and 4 are independently obtained, shifted, and deprojected using an inclination and position angles of 84° and 49°, respectively (see Section 4). Then, they are concatenated. The radial profile obtained from `frank` depends on two hyperparameters: $\alpha_{frank}$ and $w_{smooth}$. The former controls the maximum baseline that is used in the fitting process (long baselines where the signal-to-noise ratio is small are not included), and the latter controls how much the visibility model is allowed to fit all the visibilities structures from the data (see Jennings et al. 2020 for more details). We ran 1000 fits using different combinations of these two hyperparameters varying between $1.2 < \alpha_{frank} < 1.4$ and $-2 < \log_{10}(w_{smooth}) < -1$. We chose these ranges to avoid artificial high-frequency structure, and we observe no significant difference between the profiles. In addition, we chose not to allow the resulting radial profiles from `frank` to have negative values.

The average resulting radial profile and the average real part of the 1000 visibility fits are presented in Figure 2. In the top panel, we also include the cut along the major axis obtained from the image (Figure 1). We find that the structures obtained from the visibility fit and the major axis cut are in very good agreement. They share the same location, and the relative brightness between ring 1 and ring 2 is similar in both cases. Gap 1 appears deeper in the major axis cut than in the visibility fit, likely because the visibility fit is averaged over all angles

and the gaps are filled along the minor axis due to projection effects.

Starting from the outside, we find that ring 2 peaks at $0^{\prime\prime}73 \pm 0^{\prime\prime}02$, gap 1 is lowest at $0^{\prime\prime}63 \pm 0^{\prime\prime}02$, and ring 1 peaks at $0^{\prime\prime}51 \pm 0^{\prime\prime}02$. Inward of $\sim 0^{\prime\prime}5$, there is a flat region with increasing surface brightness with radius, that we name region 2. In the `frank` profile, we also detect a small peak at $\sim 0^{\prime\prime}18$ (inside region 2) that is not visible in the major axis cut. Finally, we detect some emission from the inner regions of the disk. In the major axis cut, the central emission can be described by a central marginally resolved Gaussian plus some extended emission (possibly a ring; see Section 4.2 and Figure 8). The two components are not detectable in the visibility fit, likely because this is an azimuthally averaged profile and the inner ring is not resolved in the minor axis direction, even in the visibility domain. In addition, we find that the central flux of the `frank` fit is greater than that of the major axis cut. This is because, contrary to the major axis cut, the visibility fit is an azimuthal average, which is differently affected by the central beam smearing.

In Section 4, we present a radiative transfer modeling of the continuum millimeter emission. By reproducing the features characterized with the visibility fit, we aim to constrain the physical structure of the disk and of the different grain populations. In particular, we use the presence and shape of gap 1 to constrain the vertical extent of millimeter-sized particles in the disk. As both the `frank` profile and major axis cut show similar substructures, we use the major axis profile for the rest of the analysis.

## 3.2. CO Emission

We display the $^{12}$CO moment 0 and 1 maps in Figure 3. Those represent the integrated intensity and the velocity map, respectively. Before discussing the shape of the maps, we first estimate the integrated flux of the $^{12}$CO emission, using the CASA `imstat` function on our ALMA moment 0 map. We measure the flux over a rectangle of $0^{\prime\prime}9$ along the minor axis
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Top: mean radial intensity profile of millimeter-continuum emission obtained from the visibilities with *frank* (in blue), and major axis cut obtained from the image (in orange). The major axis cut is normalized to its peak, and the *frank* profile is normalized so that the normalized intensity of ring 2 coincides in both profiles. The uncertainties on the visibility fit correspond to the typical 10% flux calibration error that dominates over the fit's uncertainty, while the uncertainty on the major axis cut corresponds to ±3σ. Bottom: observed averaged deprojected visibilities (black) and model obtained with *frank* (red).</figcaption></figure>

direction and 4′′4 along the major axis, and obtain $F_{12\text{CO}} = 6.1 \pm 0.6$ Jy km s⁻¹. We report a 10% uncertainty due to flux calibration errors.

From the overlay of the ¹²CO moment 1 map and continuum emission in Figure 3, we see that the ¹²CO emission is significantly more extended than the millimeter-continuum emission, both in the radial and vertical directions (see also Flores et al. 2021). This is consistent with large grains being affected by vertical settling (Barrière-Fouchet et al. 2005) and possibly radial drift (Weidenschilling 1977).

Thanks to the increased angular resolution, the ¹²CO moment 0 map (see left panel of Figure 3) shows significantly more details than previously detected in Flores et al. (2021), though the same overall features are present. The bottom southeast side appears brighter than the northwest side of the disk, which indicates that the southeast side is closest to us (see also Wolff et al. 2021). In addition, the disk emission can be described by two distinct regions. First, the inner 150 au (~1′′) of the disk describe an X shape, typical of very inclined systems. With the new observations, we clearly resolve the cold

midplane with little CO emission, which separates the two bright sides of the disk. This inner region is co-located with the continuum millimeter emission. On the other hand, at larger radii ($R > 200$ au, 1′′4), the disk appears flat, with a nearly constant brightness profile with radius. We find that the transition from the inner flaring X-region to the outer flat region starts at a radius of ~150 au, which roughly corresponds to the radius where the millimeter-continuum emission disk ends and the disk's scattered light stops (see Section 3.3 and Figure 4). This change in the dust structure seems to also correspond to a change in the disk gas structure.

Using lower-angular-resolution observations, Flores et al. (2021) found that this outer region is associated with a uniform temperature both radially and vertically. In Appendix A, we present an updated temperature map of the disk, resulting from the Tomographically Reconstructed Distribution method (TRD; Dutrey et al. 2017; Flores et al. 2021) applied to the new high-angular-resolution observations. We also recover the isothermal outer region. We note that Flores et al. (2021) suggested that it was due to external UV illumination through a mostly optically thin region, and refer the reader to their analysis for a more detailed interpretation.

### 3.3. Comparison with Scattered-light Image

In this section, we compare the high-angular-resolution millimeter-continuum and ¹²CO observations with an HST scattered-light image of Oph 163131 modeled in Wolff et al. (2021). The scattered-light image is expected to trace the scattering surface of small micron-sized particles, well-mixed with the gas. We present an overlay of the HST image at 0.6 µm, the ¹²CO moment 0 emission, and millimeter-continuum contours in Figure 4. As first shown in Stapelfeldt et al. (2014), the scattered-light image of Oph 163131 presents the characteristic features of edge-on disks, with two parallel nebulosities separated by a dark lane. Villenave et al. (2020) estimated that the disk size (diameter) in scattered light is 2′′5.

We first compare the apparent radial and vertical extent of the disk in the different tracers. In the radial direction, we find that the millimeter dust emission is less extended than both the scattered-light and ¹²CO emission. The outer radius of the disk millimeter dust emission (at ~150 au) appears to coincide with the $15\sigma_{\text{gas}}$ contours that also mark the limit of the inner region of ¹²CO emission (characterized by the bright X pattern split by a cold midplane; see Section 3.2). Further away, the scattered light is fainter and comes from a more diffuse region that seems to extend as far out radially as the gas emission (out to approximately 400 au). The change of intensity in scattered light suggests that small grains are either not illuminated by the central star, possibly due to a decrease of the height of their scattering surface, or that they are depleted in the outer regions of the disk (e.g., Muro-Arena et al. 2018).

Both the scattered-light and the ¹²CO emission appear significantly more extended vertically than the millimeter-continuum image, which is consistent with vertical settling occurring in the disk. The disk in scattered light seems to be as extended vertically as the ¹²CO emission, but with an additional fainter halo extending to significantly higher levels, which is due to PSF convolution (see, for example, the vertical halo also present in the PSF-convolved model of the HST image in Appendix B).

Interestingly, we also see significant differences between the morphologies seen in scattered-light and ¹²CO emission. In
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left: $^{12}$CO continuum-subtracted moment 0 map, Right: $^{12}$CO moment 1 map and continuum $5\sigma$ and $10\sigma$ contours in white. We illustrate the beam sizes with ellipses in the bottom left corner of each panel ($^{12}$CO beam in black—both panels, and continuum beam in white within the CO beam—right panel only). The cube used to make these figures is available in Zenodo at 10.5281/zenodo.6371758.</figcaption></figure>

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

<figure><img src="image_8.png" /><figcaption>**Figure 4.** HST 0.6 µm scattered-light image (colors), $^{12}$CO moment 0 ($5\sigma_{gas}$, $15\sigma_{gas}$, and $30\sigma_{gas}$ contours, in blue), and $5\sigma$ and $10\sigma$ contours of the continuum map (white). The beam sizes are shown by ellipses in the bottom left corner ($^{12}$CO beam in white, and continuum beam in black, inside the CO beam).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Position of the spines of the 0.6 µm image (green lines), and $^{12}$CO moment 0 maps (blue lines), overlaid on top of either the moment 0 map (top panel) or the scattered-light image (bottom panel).</figcaption></figure>

Figure 5, we applied the method described in Appendix D of Villenave et al. (2020) to show the spine of each nebula (peak intensity) as seen in scattered-light and $^{12}$CO emission. We find that, up to about 150 au, there is a clear increase of the $^{12}$CO apparent height with radius (blue lines). When deprojected, Flores et al. (2021) found that this increase is roughly linear. We fitted the vertical extent of the $^{12}$CO emission surface as a function of radius with a linear regression and obtain $z/r \sim 0.3$, for $r < 150$ au. The variation of the CO scattering surface with radius in Oph 163131 is relatively small compared to similar estimates in other protoplanetary disks (e.g., Pinte et al. 2018; Flores et al. 2021; Law et al. 2021). On the other hand, the opposite is seen for the scattered-light nebulae (green lines), which are extremely flat, in the sense that their apparent height does not vary significantly with projected distance. This difference in behavior is particularly clear on the least bright

side of the disk, which is at the bottom of Figure 5, and is likely due to optical depth effects.

Differences between the height of the dust scattering surface and CO emission surfaces have already been identified in other protoplanetary disks (e.g., Rich et al. 2021), and can be due to different optical depth in the different tracers. In the case of our study, we are comparing the integrated scattered-light intensity to the $^{12}$CO moment map, which is the integration of the $^{12}$CO channel maps. In the channel maps, higher velocities detected at small projected distances allow us to probe the warmer (brighter) disk interior. On the other hand, scattered light is optically thick and comes from the disk outer edge even at small projected distances, as indicated by its non-variation with projected distance. Thus, at small projected distances, it is expected that the scattered-light surface, probing the disk outer edge, appears higher than the $^{12}$CO moment 0 peak of emission, which traces the inner regions. The agreement in
height toward the outer radii, however, indicates that both components are emitted at similar altitudes above the midplane.

## 4. Radiative Transfer Modeling

In the previous section, we presented high-angular-resolution observations of Oph 163131, which reveal a highly structured disk. Given the high inclination of the disk, the presence of rings provides constraints on the vertical extent of millimeter dust particles, and on the efficiency of vertical settling. We aim to use these additional morphological constraints to refine the millimeter-continuum radiative transfer modeling of the source presented in Wolff et al. (2021), and to constrain the physical structure of the millimeter grains in the disk of Oph 163131. Building on the existing model, and to obtain a more complete view of the disk, we also compute the spectral energy distribution (SED) and the 0.6 µm image.

### 4.1. Methodology

To model the disk of Oph 163131, we use the radiative transfer code MCFOST (Pinte et al. 2006, 2009). We assume that the disk structure is axisymmetric, and we model the disk using several regions to reproduce all substructures seen in the continuum observations. Given the complexity of the ALMA continuum image, we do not aim to find a unique model, but rather a well-fitting one.

For each region, we assume a power-law distribution for the surface density:

$$ \Sigma(r) \propto r^p, \text{ for } R_{\text{in}} < r < R_{\text{out}}. \quad (1) $$

We parameterize the vertical extent of the grains by the scale height, such that

$$ H(r) = H_{100\,\text{au}}(r/100\,\text{au})^{\beta}, \text{ for } R_{\text{in}} < r < R_{\text{out}}, \quad (2) $$

where $\beta$ is the flaring exponent. For simplicity, we use astronomical silicates (similar to those shown in Figure 3 of Draine & Lee 1984)<fnref n="14" /> with a power-law distribution of the grain size following $n(a)da \propto a^{-3.5} da$.

The main free parameters of each region are thus $R_{\text{in}}$, $R_{\text{out}}$, $H_{100\,\text{au}}$, and $M_{\text{dust}}$. In all our modeling, we fixed a flaring exponent of $\beta=1.1$ for large grains, and $\beta=1.2$ for small grains (as in Villenave et al. 2019, for example). In addition, we fixed the surface density exponent to $p=-1$ in all regions, except for region 2 and ring 2, where we needed to adjust this exponent (see Section 4.2). The inclination is a global disk parameter that we constrain with this modeling to be 84° (see Section 4.2, ring 2). We assume that all regions are coplanar with this inclination. Following Wolff et al. (2021) and Flores et al. (2021), we adopt a distance of 147 pc, a stellar mass of $1.2 M_{\odot}$, stellar radius of $1.7 R_{\odot}$, and a stellar effective temperature of 4500 K, and we assume 2 mag of interstellar extinction in the SED.

For each set of parameters, we compute the 1.3 mm continuum image, SED, and 0.6 µm scattered-light image. All maps are convolved by a representative PSF before being compared to the data. We use the HST PSF generated by the TinyTim software package (Krist et al. 2011) for the 0.6 µm image. For the ALMA data, we simulate the real interferometric response by producing synthetic images using the CASA simulator.

For each model, we first generate synthetic visibility files for all three observing times and configurations, using the CASA task simobserve. We then produce synthetic images of these visibilities using the tclean function with the same weighting parameters as for the observations.

Previous modelings of Oph 163131 by Wolff et al. (2021) have demonstrated that some degree of vertical settling is needed to reproduce both the scattered-light and millimeter images. In this section and Section 4.2, we mimic dust settling by considering two separated grain populations (e.g., Keppler et al. 2018; Villenave et al. 2019). One layer is composed of large grains (10–1000 µm) aimed to model the ALMA continuum map, and the other layer includes smaller dust particles (0.01–10 µm) necessary to reproduce the SED and scattered-light image. In addition, we note that we present a complementary model in Section 5.3, using a different (parametric) settling prescription and considering a continuous distribution of dust particles. This second approach allows us to test the level of turbulence in the disk.

For the model with two layers of dust grains (this section and Section 4.2), we build the layer of small grains based on the results of the comprehensive MCMC fitting presented in Wolff et al. (2021). We have used a model where the small grains extend from 0.07 to 170 au, they have a scale height of $H_{\text{sd,100 au}} = 9.7$ au, and their surface density and flaring exponents are -1 and 1.2, respectively. The parameters of the small-grain layer are summarized in Table 1. They were chosen from the range of allowed parameters in Wolff et al. (2021) and to provide a reasonable match of the SED and scattered-light image (see Appendix B). We note that, for all regions of our model, we assumed that the gas is co-located with the dust with a gas-to-dust ratio of 100.

For the new millimeter-continuum data, our strategy followed an iterative process, looking for a representative model. We adjusted the parameters to reproduce both the surface brightness cut along the major axis and minimize the residual map, by way of visual inspection. In Section 4.2, we focus on the radial structure and adopt a scale height for the large grains of $H_{\text{Id,100 au}} = 0.5$ au at 100 au. Then, in Section 5.1, we explore and discuss this scale height assumption in more detail.

### 4.2. ALMA Continuum Model

In the right panel of Figure 1, we presented a labeled version of the continuum image of Oph 163131. We highlighted the main features that we aim to reproduce with radiative transfer modeling, namely ring 2, ring 1, the partially depleted region 2, and the central region. To reproduce the ALMA continuum emission, we thus consider five regions in our model: (1) a central point source, representing the central peak; (2) a first ring, needed to reflect the elongated structure around the central point source, labeled region 1 in Figure 1; (3) an inner region, to reproduce region 2; (4) a second ring, for ring 1; and (5) an outer ring, reproducing ring 2.

In this section, we present the main characteristics of each region, starting from the outside. Our strategy was to first fix the inner and outer radii of each region, and then we followed an iterative process to estimate the best combination of dust masses in each region to reproduce the major axis cut and residual maps. We present the model in Figures 6, 7, and 8, and its parameters in Table 1.

<footnote n="14">¹⁴ The dust properties are evaluated from: ftp://astro.princeton.edu/draine/dust/diel/eps_suvSil.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Left: continuum image of Oph 163131, Middle: Synthetic observations of the model using the CASA simulator, Right: residual map. Deviations from $3\sigma$ and $5\sigma$ are indicated by black and gray contours. In the scale bar, $\sigma$ corresponds to the rms of the data image. The maximum of the scale corresponds to the peak signal-to-noise in the data (~48.7$\sigma$).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Top: cut along the major axis, with each map normalized to its maximum and the contours corresponding to $3\sigma$ levels. Bottom: residual cut. The blue region in the residual map corresponds to $3\sigma$, and the green to $5\sigma$. Green vertical lines are at $0^{\prime\prime}63$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Zoom-in on the cut along the major axis (Figure 7), to show the central region. We show the beam size by a dashed blue Gaussian.</figcaption></figure>

**Ring 2.** We reproduce the outer emission by implementing a broad ring between 98 and 150 au. The brightness profile does not drop steeply at the edges of the disk, which we model with

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Dust surface densities of the model presented in Table 1.</figcaption></figure>

a steeper power-law exponent of $p = -6$. Because this region is the best-resolved, we used it to constrain the inclination of the system. We find that an inclination of $84 \pm 1^{\circ}$ best-matches the major axis profile, the residual map, and the minor axis size of the disk. Both a higher and lower inclination lead to axis ratios that are not compatible with the data. We note that an inclination of $84^{\circ}$ is consistent with the results of Wolff et al. (2021). We used this inclination and the observed position angle of the disk in the visibility fit presented in Section 3.1.

**Ring 1.** The millimeter-continuum image and frank profile (Figures 1 and 2) clearly show a ring centered at 73 au (~$0^{\prime\prime}5$). Because of the high inclination of the system, the ring is less clear in the major axis cut, as it is affected by projection effects. We reproduce this ring with a disk region between 60 and 87 au. After including ring 1 in our model, it became clear both by looking at the model images and the major axis cut that the region inward of it is not devoid of dust. To reproduce this feature, we introduced a region of large grains in region 2.

**Region 2.** To reproduce the increase of surface brightness with distance to the star in region 2 (between $0^{\prime\prime}09$ and $0^{\prime\prime}45$), we included a disk region between 13 and 60 au. However, based on the shape of the major axis profile, we found that a surface density exponent of $p = -1$ (as in most other disk regions) does not provide a good match to this region. With a decreasing surface density with radius, the profile of the region would not be flat but instead would drop steeply after a few au. Our results indicate that, to reproduce the major axis profile between 13 and 60 au ($0^{\prime\prime}09$ and $0^{\prime\prime}45$), region 2 needs to have an increasing surface density with radius, with $p = +1$. Such a
<table><caption>**Table 1**
Parameters for Our Radiative Transfer Model</caption><thead><tr><th></th><th>(°)</th><th></th></tr></thead><tbody><tr><td>Inclination</td><td>(°)</td><td>84</td></tr><tr><td>PA</td><td>(°)</td><td>49</td></tr><tr><td colspan="3"><em>Central point source</em></td></tr><tr><td>a<sub>min</sub> - a<sub>max</sub></td><td>(μm)</td><td>10-1000</td></tr><tr><td>R<sub>in</sub> - R<sub>out</sub></td><td>(au)</td><td>0.07-3.5</td></tr><tr><td>M<sub>dust</sub></td><td>(M<sub>☉</sub>)</td><td>4 × 10<sup>-6</sup></td></tr><tr><td>H<sub>ld,100 au</sub>, β, p</td><td>(au, -, -)</td><td>0.5, 1.1, -1</td></tr><tr><td colspan="3"><em>Region 1</em></td></tr><tr><td>a<sub>min</sub> - a<sub>max</sub></td><td>(μm)</td><td>10-1000</td></tr><tr><td>R<sub>in</sub> - R<sub>out</sub></td><td>(au)</td><td>3.5-13</td></tr><tr><td>M<sub>dust</sub></td><td>(M<sub>☉</sub>)</td><td>3 × 10<sup>-7</sup></td></tr><tr><td>H<sub>ld,100 au</sub>, β, p</td><td>(au, -, -)</td><td>0.5, 1.1, -1</td></tr><tr><td colspan="3"><em>Region 2</em></td></tr><tr><td>a<sub>min</sub> - a<sub>max</sub></td><td>(μm)</td><td>10-1000</td></tr><tr><td>R<sub>in</sub> - R<sub>out</sub></td><td>(au)</td><td>13-60</td></tr><tr><td>M<sub>dust</sub></td><td>(M<sub>☉</sub>)</td><td>1.3 × 10<sup>-5</sup></td></tr><tr><td>H<sub>ld,100 au</sub>, β, p</td><td>(au, -, -)</td><td>0.5, 1.1, +1</td></tr><tr><td colspan="3"><em>Ring 1</em></td></tr><tr><td>a<sub>min</sub> - a<sub>max</sub></td><td>(μm)</td><td>10-1000</td></tr><tr><td>R<sub>in</sub> - R<sub>out</sub></td><td>(au)</td><td>60-87</td></tr><tr><td>M<sub>dust</sub></td><td>(M<sub>☉</sub>)</td><td>4 × 10<sup>-5</sup></td></tr><tr><td>H<sub>ld,100 au</sub>, β, p</td><td>(au, -, -)</td><td>0.5, 1.1, -1</td></tr><tr><td colspan="3"><em>Ring 2</em></td></tr><tr><td>a<sub>min</sub> - a<sub>max</sub></td><td>(μm)</td><td>10-1000</td></tr><tr><td>R<sub>in</sub> - R<sub>out</sub></td><td>(au)</td><td>98-150</td></tr><tr><td>M<sub>dust</sub></td><td>(M<sub>☉</sub>)</td><td>4 × 10<sup>-5</sup></td></tr><tr><td>H<sub>ld,100 au</sub>, β, p</td><td>(au, -, -)</td><td>0.5, 1.1, -6</td></tr><tr><td colspan="3"><em>Small grains</em></td></tr><tr><td>a<sub>min</sub> - a<sub>max</sub></td><td>(μm)</td><td>0.01-10</td></tr><tr><td>R<sub>in</sub> - R<sub>out</sub></td><td>(au)</td><td>0.07-170</td></tr><tr><td>M<sub>dust</sub></td><td>(M<sub>☉</sub>)</td><td>2 × 10<sup>-5</sup></td></tr><tr><td>H<sub>sd,100 au</sub>, β, p</td><td>(au, -, -)</td><td>9.7, 1.2, -1</td></tr></tbody></table>

**Note.** Each parameter was adjusted during the modeling, except for the grain size ($a_{min} - a_{max}$) and the flaring exponent of each region ($\beta$). The parameters $H_{100 au}$ and $\beta$ were defined in Equation (2), and the surface density exponent in Equation (1).

profile is not expected in protoplanetary disks with a smooth pressure gradient, and it might indicate the presence of more complex structure, such as additional unresolved rings. We note that self-induced dust traps at the position of ring 1 might be able to reproduce a profile of increasing surface density with radius at millimeter wavelengths, corresponding to region 2 (Vericel et al. 2021; Gonzalez et al. 2017).

**Central point source and region 1.** During our modeling, we found that the innermost region of Oph 163131 cannot be correctly modeled by an unique region extending from 0.07 to 13 au. This is because, as seen in the zoomed-in major axis cut on Figure 8, the central emission is not properly reproduced by an unresolved point source convolved with the corresponding beam size: it can be decomposed into two Gaussians with different widths. We reproduce the central region with a slightly resolved central point source, plus a second ring extending up to 13 au. We fixed the inner radius of the central point source to be a conservative estimate of the sublimation radius (see Wolff et al. 2021, their Section 3.2).

To summarize, we present the millimeter model images and model parameters in Table 1 and Figures 6, 7, and 8. The SED

and scattered-light model are shown in Appendix B, and we also display the surface density of the model in Figure 9 (see Appendix C for details). The model has a total dust mass of $1.2 \times 10^{-4} M_{\odot}$, which is in agreement with the results from MCMC modeling by Wolff et al. (2021) constraining $M_{dust} > 3 \times 10^{-5} M_{\odot}$ for the scattered-light model. However, we note that, with the current modeling, the large-grain layer is partially optically thick, which implies that the dust masses presented in Table 1 are formally lower limits.

# 5. Discussion

## 5.1. Vertical Extent of Millimeter-sized Particles

Despite the high inclination of Oph 163131, the new millimeter-continuum image presented in this work reveals interesting radial substructures. Those can be used to constrain the vertical extent of large dust grains in this disk. In particular, the fact that the outer gap is well-resolved readily indicates that the millimeter-emitting dust is constrained to a very thin midplane layer. In Section 4.2, we presented a model with an extremely thin large-grain scale height of $H_{ld,100 au} = 0.5$ au. Here, in order to constrain the vertical extent of millimeter dust particles, we follow the methodology of Pinte et al. (2016) on HL Tau and look for geometrical constraints on the vertical extent of the millimeter dust in the rings. We modify our model presented in Section 4 to test larger scale heights for the large-dust layers, with $H_{ld,100 au} = 0.5, 1, 2$, and $3$ au at $100$ au. We keep the same disk radial structure, dust mass, and small-grain layer as presented in Table 1, and only vary the value of $H_{ld,100 au}$ in the large-grain layers. In Figure 10, we compare the data with the model images, a zoom-in on the major axis cut at the position of the outer gap, and minor axis cuts across ring 1.

From the model images and major axis cut, one can clearly see that, when $H_{ld,100 au}$ increases, the outer gap gets filled in due to projection effects (see also Pinte et al. 2016). By observing the major axis profiles displayed in the top right panel of Figure 10, we find that the gap is clearly less deep for the model with $H_{ld,100 au} = 3$ au than in the data, which indicates that this model is too vertically thick to reproduce the observations.

In addition, we computed minor axis cuts at $0^{\prime\prime}55$ (~80 au) from the central point source along the major axis direction, to assess how well the ring/gap is resolved in the data and in the model. In the observations, the cut along the minor axis at $0^{\prime\prime}55$ reveals three peaks: one for each side of ring 2 and the central peak corresponding to ring 1. We also see a brightness asymmetry between the two peaks corresponding to the outer ring, which is mostly due to optical depth / geometrical effects. This asymmetry varies between 20% and 40% along the major axis. When compared to the models, we find that the models with a millimeter scale height of $H_{ld,100 au} = 3$ au and $H_{ld,100 au} = 2$ au do not show three peaks in the cut. Their minor axis profiles are smoother, with only one peak, which indicates that they are too thick vertically to reproduce the observations. On the other hand, the models with a scale height of $H_{ld,100 au} = 1$ au and $H_{ld,100 au} = 0.5$ au show the three peaks. The peaks are significantly clearer for the model with a large-grain scale height of 0.5 au, which provides the best match to the data. This indicates that (at least) the outer region of the disk is extremely thin vertically, with a vertical scale height on the order of $H_{ld,100 au} \le 0.5$ au or less.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Scale height comparison. Top and bottom left panels: both sides of the continuum image of Oph 163131, Middle panels: model images of Oph 163131, with $H_{\text{ld,100 au}}$ varying from 0.5 to 3 au. Top right panel: zoom-in on gap 1 (0.63 is marked by the green vertical line) along the major axis of the disk, for the data (averaged from both sides) and models. Bottom right panel: minor axis cuts of the data and the models, at 0.55 (80 au) from the central point source, marked by a black cross on the data and model images. The uncertainty on the cuts corresponds to $1\sigma$.</figcaption></figure>

Our results are in agreement with findings from Pinte et al. (2016) who modeled the gaps and rings of HL Tau, and also with the conclusions of Villenave et al. (2020) based on the comparison of the major and minor axis profiles of HK Tau B, HH 30, and HH 48 with some fiducial radiative transfer modeling. Both studies find a scale height of millimeter grains $H_{\text{ld,100 au}}$ of about or less than one au in these systems. Our analysis of Oph 163131 thus seems to generalize the finding that grains emitting at millimeter wavelengths are extremely thin vertically in the outer regions of protoplanetary disks to a larger number of disks, pointing to the possibility that this is the case in most protoplanetary disks. We discuss some implications of this result for the formation of wide-orbit planets in Section 5.4.

## 5.2. Implications on the Degree of Dust-Gas Coupling

In the previous section, we have constrained the vertical extent of the millimeter dust particles to be $H_{\text{ld,100 au}} \le 0.5$ au at 100 au. This value is about one order of magnitude smaller than the small grains and the gas scale height obtained by Wolff et al. (2021) using MCMC fitting of the scattered-light data (of $H_{\text{g,100 au}} = H_{\text{sd,100 au}} = 9.7 \pm 3.5$ au), and is indicative of efficient settling of millimeter grains in the disk. From the difference of scale height between the gas and millimeter dust grains, one can estimate the degree of coupling of this dust with gas, which is what we aim to do in this section. More specifically we characterize the ratio $\alpha/St$, where $\alpha$ represents the turbulence strength (Shakura & Sunyaev 1973), and $St$ is the Stokes number, which describes the aerodynamic coupling between dust and gas.

In the classical 1D prescription of settling (Dubrulle et al. 1995), the vertical transport of particles is related to turbulence through a diffusion process. Assuming a balance between settling and turbulence, for grains in the Epstein regime ($St \ll 1$), and under the assumption of $z \ll H_g$, the dust scale

height can be written as follows (e.g., Youdin & Lithwick 2007; Dullemond et al. 2018):

$$H_d = H_g \left( 1 + \frac{StSc}{\alpha} \right)^{-1/2}, \quad (3)$$

where $Sc$ is the Schmidt number of the turbulence, characterizing the ratio of the turbulent viscosity over turbulent diffusivity. In this section, we follow Dullemond et al. (2018) and Rosotti et al. (2020) and assume that the Schmidt number of the turbulence is equal to $Sc = 1$.

Using Equation (3) for $H_d = H_{\text{ld,100 au}} = 0.5$ au and $H_{g,100 au} = 6.2$ au (the lowest limit from Wolff et al. 2021), we obtain an upper limit for the ratio $\alpha/St$ of $[\alpha/St]_{\text{ld,100 au}} < 6 \times 10^{-3}$, suggesting a strong decoupling between the large grains in this model and the gas in the vertical direction.

This value is relatively small compared to previous estimates of the $\alpha/St$ ratio, which have been measured both in the vertical and in the radial directions. In this paragraph, we compare our estimates with results from Doi & Kataoka (2021), also studying the strength of the coupling in the vertical direction. Doi & Kataoka (2021) estimated $[\alpha/St]_{\text{ld,100 au,HD163296}} < 1 \times 10^{-2}$ for the outer ring of HD 163296 (and $[\alpha/St]_{\text{ld,67 au,HD163296}} > 2.4$ for the inner ring). Their relatively looser constraint on the coupling strength of the outer ring likely comes from the fact that HD 163296 is significantly less inclined than Oph 163131 (47° versus 85°), making the study of its vertical structure more difficult.

On the other hand, Dullemond et al. (2018) and Rosotti et al. (2020) used the dust and gas width of several ringed systems to constrain the coupling strength in the radial direction. They consistently obtained $[\alpha/St]_{\text{ld,100 au,HD163296}} \ge 4 \times 10^{-2}$ for the second ring of HD 163296, and $[\alpha/St]_{\text{ld,120 au,AS209}} \ge 0.13$ for AS 209, which is significantly higher than our estimate of the vertical coupling on Oph 163131 and that of Doi & Kataoka (2021) in HD 163296.
<figure><img src="image_15.png" /><figcaption>Figure 11. Data (left panel) and models using the settling prescription of Fromang & Nelson (2009), with respective α parameters of 10⁻⁵ (second left panel), 10⁻⁴ (second right panel), and 10⁻³ (right panel).</figcaption></figure>

In the radial direction, it is also clear that the rings in Oph 163131 are not particularly thin (and do not appear Gaussian). In particular, there is only a factor $R_{CO,out}/R_{Id,out} \sim 2.8$ in radial size between ¹²CO and the millimeter dust emission, while the difference in the vertical direction reaches about ten times that value: $H_{g,100\,au}/H_{Id,100\,au} \sim 20$ (or at least $H_{g,100\,au}/H_{Id,100\,au} > 12$ if we assume the lowest limit from Wolff et al. 2021). As previously suggested by Doi & Kataoka (2021), these apparent inconsistencies might indicate that the turbulence level is different in the radial and vertical direction, or that the ring formation mechanism is different from what was assumed by Dullemond et al. (2018) and Rosotti et al. (2020). The first possibility is also consistent with the results from Weber et al. (2022), who produced hydrodynamical simulations of V4046 Sgr and found that, to reproduce all their observations, the turbulence in the vertical direction must be reduced compared to its value in the radial direction.

### 5.3. Vertical Settling

The modeling presented in Sections 4.2 and 5.1 and used in Section 5.2 was based on the simplified assumption of two independent layers of dust particles, where the small grains (<10 µm) have a large scale height and the large grains (>10 µm) are concentrated into the midplane. We now aim to test a more realistic settling prescription (previously used by, e.g., Pinte et al. 2016 and Wolff et al. 2021), which allows to consider a continuous distribution of dust particles. We use the Fromang & Nelson (2009) prescription for settling (their Equation (19)), the implementation of which in mcfost is described in Pinte et al. (2016). We assume that the gas vertical profile remains Gaussian and that the diffusion is constant vertically. With this prescription, each dust grain size follows its individual vertical density profile, which depends on the gas scale height, local surface density, and the turbulent viscosity coefficient α (Shakura & Sunyaev 1973). This profile reproduces relatively well the vertical extent of dust grains obtained with global MHD simulations (Fromang & Nelson 2009). For completeness, we report the vertical density profile for a grain size a that we adopt in this section:

$$ \rho(r, z, a) \propto \Sigma(r) \exp \left[ -\frac{\Omega \tau_S(a)}{\tilde{D}} (e^{2H_g^2} - 1) - \frac{z^2}{2H_g^2} \right], \quad (4) $$

with $\Omega = \sqrt{GM_*/r^3}$, $\tilde{D} = H_g/\rho_c\Omega$, and $\tau_S(a) = \rho_d a/\rho_c c_S$, and where $c_S = H_g\Omega$ is the midplane sound speed, $H_g$ is the gas scale

height, $S_c$ is the Schmidt number that is fixed to 1.5 in mcfost (Pinte et al. 2016), $\rho_d$ is the dust material density, and $\rho_g$ is the gas density in the midplane. We note that the Stokes number described in Section 5.2 corresponds to $St = \Omega\tau_S(a)$. The degree of settling is set by varying the α parameter. We also note that, when $z \ll H_g$, the dust density defined in Equation (4) reduces to a simple Gaussian function with a scale height defined by Equation (3) (e.g., Riols & Lesur 2018).

The value of the turbulence parameter is usually predicted to range from $\alpha = 10^{-3} - 10^{-2}$ when driven by the MRI (magnetorotational instability, e.g., Balbus & Hawley 1991). On the other hand, in regions with low ionization and weak coupling between the gas and magnetic fields, recent studies have showed that hydrodynamic or non-ideal MHD effects may dominate, and could lead to turbulent parameters as low as $\alpha = 10^{-4}-10^{-6}$ (Bai & Stone 2013; Flock et al. 2020). In this context, we produced models for different settling strengths, obtained for respective α parameters of $10^{-5}$, $10^{-4}$, and $10^{-3}$. For this modeling, we assume that the grain size integrated over the whole disk follows a simplified power-law distribution of $n(a)da \propto a^{-3.5}da$, with minimum and maximum grain sizes of 0.01 and 1000 µm, respectively. We note that the power-law size distribution may vary as a function of vertical height, due to the presence of vertical settling in our model (Sierra & Lizano 2020). We distribute the grains over five radial regions coincident with the large-grain regions presented in Table 1. All five regions are normalized to a gas scale height of $H_{g,100\,au} = 9.7$ au at 100 au with a flaring of $\beta = 1.1$, but the large grains will have a smaller scale height than the gas because of the settling prescription considered ($H_{Id,100\,au} < H_{g,100\,au}$). We adjust the mass of each region to obtain a model that represents the millimeter image relatively well, and more specifically the levels of the surface brightness profile of the major axis. The total dust mass in these models varies between 7.3 and $9.9 \times 10^{-5} M_\odot$. As for Section 4.2, we produced synthetic images of the models using the CASA simulator.

We present the 1.3 mm model images obtained for an α parameter of $10^{-5}$, $10^{-4}$, and $10^{-3}$ in Figure 11. For comparison with Section 5.1, the respective scale heights of 1 mm-sized particles obtained in these models are $H_{1\,mm,100\,au} = 0.95$ au, 2.7 au, and 5.7 au. Similarly to the previous section, the effect of vertical settling is clearly visible in these models. Focusing on the outer regions of the disk (ring 1, gap 1, and ring 2), we find that the highest turbulence level leads to a high millimeter dust scale height, making the gap blurred due to projection effects.
On the other hand, the lowest tested turbulence of $\alpha = 10^{-5}$ efficiently settles the large grains to the midplane and reproduces well the shape of the millimeter emission. When performing a similar study for various other gas scale heights, between $H_{g,100\,au} = 6.2\,au$ and $13.2\,au$, compatible with the range of posterior values estimated by Wolff et al. (2021), we also found that $\alpha = 10^{-5}$ provides the best match with the observations.

Thus, we find that the millimeter observations of Oph 163131 are consistent with a turbulent viscosity coefficient of $\alpha \lesssim 10^{-5}$, at least in the outer regions of the disk ($r \sim 100\,au$). Together with results from turbulent line broadening in several protoplanetary disks (e.g., Flaherty et al. 2017, 2018, 2020), as well as those of Pinte et al. (2016) using a similar technique on HL Tau, our study suggests that the turbulence is low in the outer regions of protoplanetary disks. Additionally, we note that, when combining these results with the coupling estimate from Section 5.2, we find that the Stokes number of particles emitting at millimeter wavelengths must be greater than $St > 1.7 \times 10^{-3}$ at $100\,au$.

Finally, we note that the settling prescription implemented in this section depends on many parameters (e.g., grain size distribution, settling prescription, same turbulence in all the disk, constant grain opacity, ...) and that it might be necessary to add more complexity to the models in order for them to be able to simultaneously reproduce the millimeter, scattered-light images, and SED. Nevertheless, the high-resolution data sets available for the highly inclined disk in Oph 163131 allows us to obtain independent measures of scale height for multiple dust sizes, which is extremely valuable to test and improve current dust settling models.

## 5.4. Potential for Wide-orbit Planet Formation in Oph 163131

The exceptionally well-characterized outer disk of Oph 163131 is strongly indicative of dust growth and settling. A significant fraction of the dust mass appears to have grown to near-mm sizes and to have settled toward a dense disk midplane layer, with $H_{1d,100\,au}/H_{g,100\,au} \sim 0.05$ at $100\,au$. These two processes, dust growth and settling, are widely recognized as the key first steps for planetesimal formation and planetary growth (for example, see reviews by Youdin & Kenyon 2013; Johansen & Lambrechts 2017; Ormel et al. 2017). Nevertheless, the outer parts of protoplanetary disks ($\gtrsim 50\,au$) are a challenging environment for planet formation, due to inherently low solid surface densities and long orbital timescales. In this section, we address the potential for wide-orbit planet formation in Oph 163131, motivated by the ring-like features and dust depletion pattern that could be indicative of ongoing planet formation.

A common suggestion for the formation of wide-orbit giant planets is that such planets emerge through the direct fragmentation of young, massive, gravitationally unstable disks (Helled et al. 2014). However, the apparent quiescent nature of Oph 163131, with signs of extremely low vertical stirring ($\alpha \lesssim 10^{-5}$) and little-to-no accretion onto the star (Flores et al. 2021), rules out this scenario. Possibly, the formation of such gravitational-instability planets could have taken place earlier in the disk lifetime, but the early formation of such massive planets—exceeding Jupiter in mass (Kratter et al. 2010)—would have resulted in a deep inner gas and dust cavities, which are not seen. Instead, as we argue below, planet formation could have proceeded more gently through core accretion via pebble accretion.

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Time to grow a planetary embryo of $10^{-2} M_E$ to a fully grown core of $10 M_E$ at $r = 100\,au$, as given by the gray contours, depending on the pebble scale height and local pebble-to-gas surface density ratio, assuming pebbles $a = 0.5\,mm$ in size. The vertical blue dashed line gives the estimate for the pebble scale height in Oph 163131. The red point corresponds to a disk with solar dust-to-gas ratio, with a 75% mass fraction locked in pebbles. The vertical orange dashed line corresponds to the assumption of a massive gas disk, with a total mass of $0.1 M_\odot$. The gray area shows where the midplane pebble-to-gas ratios is equal to, or above, unity, which corresponds to the parameter region in line with planetesimal formation through the streaming instability. Pebble scale heights exceeding the one inferred here for Oph 163131 would suppress planet formation in wide orbits.</figcaption></figure>

In the core-accretion scenario (Mizuno 1980; Pollack et al. 1996), the formation of a giant planet is initiated by the formation of a solid icy/rocky core that accretes a H/He gaseous envelope of varying mass, resulting in either an ice giant—with an envelope mass smaller than the core mass—or a gas giant. However, the formation of a typical giant-planet core of about 10 Earth masses ($M_E$) cannot occur solely through the process of planetesimal accretion at orbits outside 10 au: formation times would exceed disk lifetimes even when assuming a planetesimal mass budget exceeding a solar dust-to-gas ratio by factor 10 (Rafikov 2004; Kobayashi et al. 2010). On the other hand, core growth can be accelerated by the direct accretion of pebbles, provided that, as is the case in Oph 163131, they are abundantly present and concentrated toward the disk midplane.

The pebble accretion efficiency is highly dependent on the vertical pebble scale height, because only when the pebble accretion radius exceeds the pebble scale height does the accretion reach maximal efficiency (in the so-called 2D Hill accretion branch; Lambrechts & Johansen 2012). Figure 12 illustrates the growth timescale as function of the pebble scale height and the pebble-to-gas ratio (following Lambrechts & Johansen 2012; Lambrechts et al. 2019). Pebbles are assumed to be spherical with a radius of $a = 0.5\,mm$. We take for the pebble surface density the values inferred in Section 4.2 (see Table 1 and Figure 9). However, we note that regions in the disk may be optically thick, which would result in higher dust-to-gas ratios than used here. We present two assumptions on the gas surface density. One where the dust-to-gas ratio is the nominal solar value of 0.01, with a 75%-mass fraction locked
in pebbles (consistent with our model from Section 4.2; see also Appendix D) and one where we assume a total disk mass of $0.1 M_\odot$. The latter is an upper limit, given that the gas disk is likely less massive: there is no strong gas accretion onto the star detected (Flores et al. 2021) and no evidence for asymmetries in the disk (potentially triggered by gravitational instabilities; see, e.g., Hall et al. 2019, 2020). We choose a location of $r=100$ au, corresponding to the inner edge of outer ring 2.

The disk around Oph 163131 appears to be conductive to core growth via pebble accretion within a range of reasonable values for the pebble-to-gas surface density, as shown in Figure 12. We find that a $10^{-2} M_E$ embryo can grow to up to $10 M_E$ at 60 au from the central star, within less than 10 Myr (the maximum lifetime of protoplanetary disks; see Ribas et al. 2014). Importantly, a pebble scale height of $H_{pebbles}/H_g \sim 0.05$, consistent with $H_{1d,100\,au}/H_{g,100\,au}$ (see Sections 4.2 and 5.1, and also Appendix C), strongly promotes core growth, while pebble scale heights that are a factor 10 larger would largely suppress core growth by pebble accretion. Furthermore, in Appendix D, we show that core-growth timescales do not strongly depend on the chosen location. We present results near the inner and outer edges of ring 1 at 61 and 85 au, respectively. In contrast, our results do depend on the chosen particle size. In Appendix D, we illustrate that particles need to grow beyond 0.1 mm in size for pebble accretion to drive core growth, which is consistent with what is generally inferred in protoplanetary disks (e.g., Carrasco-González et al. 2019; Macías et al. 2021; Tazzari et al. 2021).

A possible issue is that core growth via pebble accretion needs to be initiated by a sufficiently large seed mass that results from an episode of planetesimal formation. Planetesimals are believed to be the results of the gravitational collapse of pebble swarms that self-concentrate through the streaming instability in the pebble midplane (Youdin & Lithwick 2007; Johansen & Youdin 2007). The streaming instability requires dust growth to pebble sizes and a high degree of pebble settling (approaching midplane pebble-to-gas ratios near unity). Although planetesimal formation may be unlikely to be ongoing now (see the gray upper triangle in Figure 12 marking where the midplane pebble-to-gas ratios is larger than one), the combined lack of strong vertical pebble stirring and the possibility of local moderate increases in pebble concentrations, would be conductive to planetesimal formation (Johansen et al. 2009; Bai & Stone 2010; Carrera et al. 2021). When planetesimal formation is triggered in simulations of the streaming instability, the analysis of the typical planetesimal mass distribution argues for large planetesimals in the outer disk (Schäfer et al. 2017). Following the mass scaling of Liu et al. (2020), we find planetesimals in the exponential tail of the mass distribution could reach masses exceeding $0.01 M_E$ beyond 60 au (as assumed in Figure 12).

In summary, the low observed pebble scale height of the protoplanetary disk around Oph 163131 is conducive to planetary growth by pebble accretion, even in wide orbits, exceeding 50 au. If indeed planet formation was initiated in these outer regions, it could possibly be consistent with the depletion of pebbles in the inner disk ($<60$ au) and the ring-like features identified in this work. Thus, the Oph 163131-disk provides us with an exciting opportunity to study possibly ongoing planet formation at large orbital radii.

# 6. Summary and Conclusions

In this paper, we present new high-angular-resolution ALMA continuum and $^{12}$CO millimeter images of the highly inclined disk Oph 163131. The $^{12}$CO image can be described by two distinct regions: (1) an inner region ($R \lesssim 150$ au) where the disk describes a clear X shape, with a linear increase of the gas emission height with radius as of $z/r \sim 0.3$, and (2) a flat outer region (for $R > 200$ au), with uniform brightness and temperature. In contrast with the $^{12}$CO emission, the scattered-light surface brightness of HST $0.6 \ \mu$m data appears flat at all radii. This indicates that scattered light comes from the outer regions of the disk, due to the high optical depth of micron-sized particles in Oph 163131. We find that the heights of scattered-light and $^{12}$CO emission are similar at large radii, which suggests that both components are emitted at similar altitudes above the midplane. In addition, we find that the millimeter-continuum emission from larger grains is less extended in both the radial and vertical direction compared to the scattered-light and gas emission, which is in agreement with expectations from vertical settling and possibly radial drift.

Our millimeter-continuum observations of Oph 163131 reveal clear rings, which remained undetected with previous lower-angular-resolution observations. From the outside in, the disk shows two rings separated by a clear gap, some emission inside the first ring, and some bright central emission. We performed comprehensive radiative transfer modeling of the ALMA continuum image in order to constrain the physical structure of the source. In particular, we used the resolved outer ring, located at $\sim 100$ au, to add strong geometrical constraints on the vertical extent of millimeter-sized particles. Our modeling of Oph 163131 indicates that the grains emitting at 1.3 mm are extremely thin vertically, with a vertical scale height at 100 au of $H_{1d,100\,au} \lesssim 0.5$ au. Because the vertical extent of small dust particles (and indirectly the gas) has been constrained to be $\sim 10$ at 100 au (Wolff et al. 2021), this is clear evidence that vertical settling is occurring in the disk.

Using a classical 1D prescription of settling and our estimate of dust and gas scale height, we estimate the degree of coupling of dust and gas to be $[\alpha/St]_{100\,au} < 6 \times 10^{-3}$. This is value is particularly low compared to previous estimates in protoplanetary disks.

We also aimed at constraining the turbulence parameter $\alpha$ in Oph 163131. To do so, we produced three additional radiative transfer models assuming the settling model from Fromang & Nelson (2009), with $\alpha = 10^{-3}, 10^{-4}$, and $10^{-5}$, and a gas vertical extent of 9.7 au at 100 au. We find that the coefficient of turbulent viscosity needs to be extremely low in the disk to reproduce the observations, on the order of $\alpha \lesssim 10^{-5}$ in the outer regions of the disk.

Finally, we used our results to test the pebble-accretion scenario in the outer regions of the disk. The remarkably small pebble scale height of Oph 163131 is particularly favorable for pebble accretion: we show a 10 Earth mass planet can form in the outer disk, between approximately 60 and 100 au, within less than 10 Myr. If, on the other hand, the dust scale height had been 10 times larger than our observational constraint, then this would have largely suppressed core growth via pebble accretion. Thus, the extreme vertical settling measured in Oph 163131 could be the origin of the formation of wide-orbit planets. Further constraints on vertical settling in a larger
number of protoplanetary disks might provide more insights
into the dominant formation mechanism of wide-orbit planets.

We thank the anonymous referee for providing a detailed
revision of our work, which helped to improve the quality of
this study. This paper makes use of the following ALMA
data: ADS/JAO.ALMA#2018.1.00958.S and 2016.1.00771.
S. ALMA is a partnership of ESO (representing its member
states), NSF (USA) and NINS (Japan), together with NRC
(Canada), MOST and ASIAA (Taiwan), and KASI (Republic
of Korea), in cooperation with the Republic of Chile. The Joint
ALMA Observatory is operated by ESO, AUI/NRAO and
NAOJ. The National Radio Astronomy Observatory is a
facility of the National Science Foundation operated under
cooperative agreement by Associated Universities, Inc. The
research of M.V. was supported by an appointment to the
NASA Postdoctoral Program at the NASA Jet Propulsion
Laboratory, administered by Universities Space Research
Association under contract with NASA. This project has
received funding from the European Union's Horizon 2020
research and innovation program under the Marie Skłodowska-
Curie grant agreement No 210021. G.D. acknowledges support
from NASA grants NNX15AC89G and NNX15AD95G/
NExSS as well as 80NSSC18K0442. M.B. acknowledges
funding from the European Research Council (ERC) under the
European Union's Horizon 2020 research and innovation
program (grant PROTOPLANETS No. 101002188). C.P.
acknowledges funding from the Australian Research Council
via FT170100040, and DP180104235. K.R.S. acknowledges
support from the NASA Exoplanet Exploration Program
Office. A.S. acknowledges support from ANID/CONICYT
Programa de Astronomía Fondo ALMA-CONICYT 2018
31180052.

Software: CASA (McMullin et al. 2007), mcfost (Pinte
et al. 2006, 2009), frank (Jennings et al. 2020), Matplotlib
(Hunter 2007), Numpy (Harris et al. 2020), TinyTim (Krist
et al. 2011).

# Appendix A
## Tomographically Reconstructed Distribution

Flores et al. (2021) first obtained a 2D temperature map of
Oph 163131 using the Tomographically Reconstructed Dis-
tribution method (TRD; Dutrey et al. 2017). In Figure 13, we
present an updated temperature map of Oph 163131 obtained
with the same technique from our higher-angular-resolution
observations, after continuum subtraction. As in Section 4, we
assumed a distance of 147 pc, a mass of 1.2 M☉, a beam size of
0.′′081 × 0.′′072, position angle of 49°, and an inclination angle
of 84°. Each pixel is about 5 au wide.

Our new temperature map of Oph 163131 presents the same
global features as previously obtained by Flores et al. (2021)
from lower-angular-resolution observations, but with signifi-
cantly more details. In particular, the new map now clearly
resolves the cold midplane in the inner regions (R < 200 au) of
the disk. The new high-resolution data show no evidence of a
thin, cold midplane at radii beyond 200 au. This is comparable
with the outer radius of the millimeter dust, confirming that the
shielded, high optical depth region of the disk ends at a radius
of ~150–200 au.

The increased angular resolution also reduces the impact of
beam smearing in the inner regions of the disk, which explains
why peak temperatures almost twice as large are obtained with
this new map compared to the previously published map. On
the other hand, the reduced spectral resolution implies a larger
effect of Keplerian velocity shear at large radii compared to the
previous observations (see radial extent of the ellipses in
Figure 13). As a reminder, we note that the Keplerian shear
velocity is calculated as $dr = 2rdv/v_K(r)$, where $v_K(r)$ is the
Keplerian velocity of the disk, and $dv$ is the local line width. In
our data set, the local line width is dominated by the velocity
resolution of 0.7 km s⁻¹. In addition, we recover the isothermal
region in the outer 200 au of the disk. We refer to Flores et al.
(2021) for detailed analysis and interpretation of the full
temperature structure of this disk.

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Tomographically reconstructed distribution of the ¹²CO. The reconstructed temperature maps are shown in physical units of radius and altitude, both in au. The gray ellipses shown on the bottom left part of the plot illustrate how the Keplerian velocity shear affects the resolution as a function of radius.</figcaption></figure>
## Appendix B
### SED and Scattered-light Model Image

In Section 4, to mimic dust vertical settling, we have added a layer of small dust particles located higher up in the disk than the large grain dust (see Table 1). This allows us to obtain a

relatively good agreement with the 0.6 µm maps and SED, which are presented in Figures 14 and 15, respectively. We note that the photometry points used for the SED were obtained from Wolff et al. (2021).

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Left: HST 0.6 µm image, from Wolff et al. (2021). Middle: 0.6 µm model obtained after convolution by a representative PSF. The x- and y-ticks indicate 0.5σ levels of the data and model, respectively. Right: residual map.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Spectral energy distribution of Oph 163131 (blue circles) and our model prediction (red line). The model is corrected for interstellar extinction with an Av of 2 (see Section 4.1).</figcaption></figure>
## Appendix C
### Model Surface Density

In Figure 9, we presented the surface density of our model from Section 4.2. We estimated it directly from our mcfost grid of the dust density (obtained with the command `-disk_struct`) following

$$ \text{Surface density}(r) = \sum_z \left( \sum_a \rho(a, r, z) \times m(a) \times z_{\text{cell}}(r) \right), \quad (C1) $$

where $\rho(a, r, z)$ is the dust density of a particle of size $a$ at a radius $r$ and a height $z$, $m(a)$ is the mass of a particle of size $a$, and $z_{\text{cell}}(r)$ is the height of one mcfost cell at the radius $r$.

We note that the gas surface density is not displayed in this plot, as we simply assumed a gas-to-dust ratio of 100 in our modeling. In addition, the large dust is likely optically thick, which implies that the dust surface densities are formally lower limits.

## Appendix D
### Pebble Accretion Timescales: Parameter Dependency

In this section, we explore the dependency of the core-growth timescale on two model parameters: the orbital radius of the core and particle size. Results inside 100 au, for the inner and outer edges of ring 1, respectively, are shown in Figure 16. We have taken for the pebble surface densities the values shown in Figure 9 at 61 and 85 au, respectively. For the gas scale height, we have used values based on the small dust component radius scaling presented in Wolff et al. (2021). We furthermore have assumed that the pebble-to-gas scale height is constant with orbital radius, and set the value to the observed ratio of $H_{\text{ld,100 au}}/H_{\text{g,100 au}} = 0.05$. We do, however, caution that some disks, such as HD 163296, have strong indications of radial changes in the pebble scale height between rings located in a similar orbital range between 60 and 100 au (Doi & Kataoka 2021).

The weak orbital dependency in growth timescales, from approximately 60–100 au (see also Figure 12), is the result of the here-assumed fixed particle size and surface density profile: the decrease in Stokes numbers on closer orbits reduces the accretion efficiency, which is compensated for by the increase in the pebble surface density. We further note that, in the 60 au region, a planet with a mass exceeding ~10 $M_{\text{E}}$ can already accrete more than 50% of the mass flux of pebbles drifting

<figure><img src="image_20.png" /><figcaption>Figure 16. Dependency of the pebble accretion core-growth rates on the orbital radius. The left panel shows results at $r = 61$ au and the right panel at $r = 85$ au, corresponding to the inner and outer edges of ring 1.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 17. Dependency of the pebble accretion core-growth rates on the particle radius. The left panel shows results for pebbles $a = 0.1$ mm in size, and the right panel for larger pebbles exclusively $a = 1$ mm in size.</figcaption></figure>

inward. This would cause a corresponding dip in the pebble surface density, similar to the one seen in Figure 9, which could be further deepened as the planet reaches its pebble isolation mass (Lambrechts et al. 2014; Bitsch et al. 2018).

The model presented in Section 4.2 has a dust size distribution where 75% of the total mass is in dust particles larger than 0.1 mm. We therefore assume the pebble surface density to make up 75% of the total dust density. The exact mass fraction and upper limit of the mass are, however, not fully constrained. We show in Figure 17, our results when all pebbles have a size of 1 mm or 0.1 mm, covering the size range beyond particles can grow in observed protoplanetary discs (e.g., Carrasco-González et al. 2019; Macías et al. 2021; Tazzari et al. 2021). Core growth within disk lifetimes of approximately 10 Myr does require that dust particles grow beyond 0.1 mm in size, consistent with the here-presented model of the disk of Oph 163131.

In summary, we find that the timescale to grow an embryo to a $10 M_{\text{E}}$ core depends weakly on the orbital radius, but strongly on the maximal particle size for typical top-heavy particle size distributions.

## ORCID iDs

M. Villenave https://orcid.org/0000-0002-8962-448X

K. R. Stapelfeldt https://orcid.org/0000-0002-2805-7338

G. Duchêne https://orcid.org/0000-0002-5092-6464

F. Ménard https://orcid.org/0000-0002-1637-7393

M. Lambrechts https://orcid.org/0000-0001-9321-5198

A. Sierra https://orcid.org/0000-0002-5991-8073

C. Flores https://orcid.org/0000-0002-8591-472X

S. Wolff https://orcid.org/0000-0002-9977-8255

Á. Ribas https://orcid.org/0000-0003-3133-3580

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

N. Cuello https://orcid.org/0000-0003-3713-8073

C. Pinte https://orcid.org/0000-0001-5907-5179

## References

ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJL, 808, L3

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

Bai, X.-N., & Stone, J. M. 2010, ApJ, 722, 1437

Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Barrière-Fouchet, L., Gonzalez, J. F., Murray, J. R., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185

Bitsch, B., Morbidelli, A., Johansen, A., et al. 2018, A&A, 612, A30

Burrows, P. E., Stapelfeldt, K. R., Watson, A. M., et al. 1996, ApJ, 473, 437

Carrasco-González, C., Sierra, A., Flock, M., et al. 2019, ApJ, 883, 71
Carrera, D., Simon, J. B., Li, R., Kretke, K. A., & Klahr, H. 2021, AJ, 161, 96

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

Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icar, 114, 237

Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJL, 869, L46

Dutrey, A., Guilloteau, S., Piétu, V., et al. 2017, A&A, 607, A130

Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150

Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117

Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155

Flores, C., Duchêne, G., Wolff, S., et al. 2021, AJ, 161, 239

Fromang, S., & Nelson, R. P. 2009, A&A, 496, 597

Ginski, C., Stolker, T., Pinilla, P., et al. 2016, A&A, 595, A112

Gonzalez, J. F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984

Hall, C., Dong, R., Rice, K., et al. 2019, ApJ, 871, 228

Hall, C., Dong, R., Teague, R., et al. 2020, ApJ, 904, 148

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

Helled, R., Bodenheimer, P., Podolak, M., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 643

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

Hunter, J. D. 2007, CSE, 9, 90

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

Johansen, A., & Lambrechts, M. 2017, AREPS, 45, 359

Johansen, A., & Youdin, A. 2007, ApJ, 662, 627

Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJL, 704, L75

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

Kobayashi, H., Tanaka, H., Krivov, A. V., & Inaba, S. 2010, Icar, 209, 836

Kratter, K. M., Murray-Clay, R. A., & Youdin, A. N. 2010, ApJ, 710, 1375

Krist, J. E., Hook, R. N., & Stoehr, F. 2011, Proc. SPIE, 8127, 81270J

Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32

Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35

Lambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, A&A, 627, A83

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

Liu, B., Lambrechts, M., Johansen, A., Pascucci, I., & Henning, T. 2020, A&A, 638, A88

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

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

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. F. Hill & D. J. Bell (San Francisco, CA: ASP), 127

Mizuno, H. 1980, PThPh, 64, 544

Muro-Arena, G. A., Dominik, C., Waters, L. B. F. M., et al. 2018, A&A, 614, A24

Ormel, C. W., Pessah, M., Gressel, O., et al. 2017, in Astrophysics and Space Science Library, Formation, Evolution, and Dynamics of Young Solar Systems, ed. M. Pessah & O. Gressel, Vol. 445 (New York: Springer), 197

Ortiz-León, G. N., Loinard, L., Dzib, S. A., et al. 2018, ApJL, 869, L33

Ortiz-León, G. N., Loinard, L., Kounkel, M. A., et al. 2017, ApJ, 834, 141

Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJL, 860, L13

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icar, 124, 62

Rafikov, R. R. 2004, AJ, 128, 1348

Ribas, Á., Merín, B., Bouy, H., & Maud, L. T. 2014, A&A, 561, A54

Rich, E. A., Teague, R., Monnier, J. D., et al. 2021, ApJ, 913, 138

Riols, A., & Lesur, G. 2018, A&A, 617, A117

Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 495, 173

Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69

Segura-Cox, D. M., Schmiedeke, A., Pineda, J. E., et al. 2020, Natur, 586, 228

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Sierra, A., & Lizano, S. 2020, ApJ, 892, 136

Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, ApJS, 257, 14

Stapelfeldt, K. R., Duchêne, G., Perrin, M., et al. 2014, in IAU Symp. 299, Exploring the Formation and Evolution of Planetary, ed. M. Booth, B. C. Matthews, & J. R. Graham (Cambridge: Cambridge Univ. Press), 99

Tazzari, M. 2017, mtazzari/uvplot: v0.1.1, Zenodo, doi:10.5281/zenodo.1003113

Tazzari, M., Clarke, C. J., Testi, L., et al. 2021, MNRAS, 506, 2804

Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJ, 860, L12

Vericel, A., Gonzalez, J.-F., Price, D. J., Laibe, G., & Pinte, C. 2021, MNRAS, 507, 2318

Villenave, M., Benisty, M., Dent, W. R. F., et al. 2019, A&A, 624, A7

Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164

Watson, A. M., Stapelfeldt, K. R., Wood, K., & Ménard, F. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. Arizona Press), 523

Weber, P., Casassus, S., & Pérez, S. 2022, MNRAS, 510, 1612

Weidenschilling, S. J. 1977, MNRAS, 180, 57

Wolff, S. G., Duchêne, G., Stapelfeldt, K. R., et al. 2021, AJ, 161, 238

Wolff, S. G., Perrin, M. D., Stapelfeldt, K., et al. 2017, ApJ, 851, 56

Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459

Youdin, A. N., & Kenyon, S. J. 2013, in Planets, Stars and Stellar Systems, ed. T. D. Oswalt, L. M. French, & P. Kalas (Dordrecht: Springer), 1

Youdin, A. N., & Lithwick, Y. 2007, Icar, 192, 588