A Highly Settled Disk around Oph163131

We present the continuum observations of Oph 163131 inFigure 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 (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.

Zoom In Zoom Out Reset image size

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 export_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: α_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 < α_frank < 1.4 and $-2\lt {\mathrm{log}}_{10}({w}_{\mathrm{smooth}})\lt -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.

Zoom In Zoom Out Reset image size

Starting from the outside, we find that ring 2 peaks at 073 ± 002, gap 1 is lowest at 063 ± 002, and ring 1 peaks at 051 ± 002. Inward of ∼05, 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 ∼018 (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.

We display the ¹²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 ¹²CO emission, using the CASA imstat function on our ALMA moment 0 map. We measure the flux over a rectangle of 09 along the minor axis direction and 44 along the major axis, and obtain ${F}_{{}^{12}\mathrm{CO}}=6.1\pm 0.6$ Jy km s⁻¹. We report a 10% uncertainty due to flux calibration errors.

Zoom In Zoom Out Reset image size

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, 14), 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.

Zoom In Zoom Out Reset image size

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.

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 25.

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σ_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 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 ¹²CO emission. We find that, up to about 150 au, there is a clear increase of the ¹²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 ¹²CO emission surface as a function of radius with a linear regression and obtain z/r ∼ 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.

Zoom In Zoom Out Reset image size

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 ¹²CO moment map, which is the integration of the ¹²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 ¹²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.

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:

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)\propto {r}^{p},\mathrm{for}\ {R}_{\mathrm{in}}\lt r\lt {R}_{\mathrm{out}}.\end{eqnarray} \tag{ 1 }$$

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

$$\begin{eqnarray}&&H{(r)={H}_{100\,\mathrm{au}}(r/100\,\mathrm{au})}^{\beta },\mathrm{for}\ {R}_{\mathrm{in}}\lt r\lt {R}_{\mathrm{out}},\end{eqnarray} \tag{ 2 }$$

where β is the flaring exponent. For simplicity, we use astronomical silicates (similar to those shown in Figure 3 of Draine & Lee 1984) ¹⁴ with a power-law distribution of the grain size following n(a)da ∝ a^−3.5 da.

The main free parameters of each region are thus R_in, R_out, H_{100 au}, and M_dust. In all our modeling, we fixed a flaring exponent of β = 1.1 for large grains, and β = 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_⊙, stellar radius of 1.7 R_⊙, 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_{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.

Table 1. Parameters for Our Radiative Transfer Model

Inclination	(°)	84
PA	(°)	49
		Central point source
${a}_{\min }-{a}_{\max }$	(μm)	10–1000
Rin − Rout	(au)	0.07–3.5
Mdust	(M⊙)	4 × 10−6
Hld,100 au, β, p	(au, −, −)	0.5, 1.1, −1
		Region 1
${a}_{\min }-{a}_{\max }$	(μm)	10–1000
Rin − Rout	(au)	3.5–13
Mdust	(M⊙)	3 × 10−7
Hld,100 au, β, p	(au, −, −)	0.5, 1.1, −1
		Region 2
${a}_{\min }-{a}_{\max }$	(μm)	10–1000
Rin − Rout	(au)	13–60
Mdust	(M⊙)	1.3 × 10−5
Hld,100 au, β, p	(au, −, −)	0.5, 1.1, +1
		Ring 1
${a}_{\min }-{a}_{\max }$	(μm)	10–1000
Rin − Rout	(au)	60–87
Mdust	(M⊙)	4 × 10−5
Hld,100 au, β, p	(au, −, −)	0.5, 1.1, −1
		Ring 2
${a}_{\min }-{a}_{\max }$	(μm)	10–1000
Rin − Rout	(au)	98–150
Mdust	(M⊙)	4 × 10−5
Hld,100 au, β, p	(au, −, −)	0.5, 1.1, −6
		Small grains
${a}_{\min }-{a}_{\max }$	(μm)	0.01–10
Rin − Rout	(au)	0.07–170
Mdust	(M⊙)	2 × 10−5
Hsd,100 au, β, p	(au, −, −)	9.7, 1.2, −1

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

Download table as:  ASCIITypeset image

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_{ld,100 au} = 0.5 au at 100 au. Then, in Section 5.1, we explore and discuss this scale height assumption in more detail.

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 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 ± 1° 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° 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 (∼05). 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 009 and 045), 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 (009 and 045), region 2 needs to have an increasing surface density with radius, with p = +1. Such a 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 × 10⁻⁴ M_⊙, which is in agreement with the results from MCMC modeling by Wolff et al. (2021) constraining M_dust > 3 × 10⁻⁵ M_⊙ 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 055 (∼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 055 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} ≤ 0.5 au or less.

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_{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.

In the previous section, we have constrained the vertical extent of the millimeter dust particles to be H_{ld,100 au} ≤ 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_{g,100 au} = H_{sd,100 au} = 9.7 ± 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 α/St, where α 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 ≪ 1), and under the assumption of z ≪ H_g , the dust scale height can be written as follows (e.g., Youdin & Lithwick 2007; Dullemond et al. 2018):

$$\begin{eqnarray}&&{H}_{d}={H}_{g}{\left(1+\displaystyle \frac{{StSc}}{\alpha }\right)}^{-1/2},\end{eqnarray} \tag{ 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_{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 α/St of ${[\alpha /{{St}}_{\mathrm{ld}}]}_{100\,\mathrm{au}}\,\lt 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 α/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 [α/St]_{100 au,HD163296} < 1 × 10⁻² for the outer ring of HD 163296 (and [α/St]_{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 [α/St]_{100 au,HD163296} ≥ 4 × 10⁻² for the second ring of HD 163296, and [α/St]_{120 au,AS209} ≥ 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.

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_ld,out ∼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_{ld,100 au} ∼ 20 (or at least H_{g,100 au}/H_{ld,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.

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:

$$\begin{eqnarray}&&\rho (r,z,a)\propto {\rm{\Sigma }}(r)\exp \left[-\displaystyle \frac{{\rm{\Omega }}{\tau }_{S}(a)}{\tilde{D}}\left({e}^{\tfrac{{z}^{2}}{2{H}_{g}^{2}}}-1\right)-\displaystyle \frac{{z}^{2}}{2{H}_{g}^{2}}\right],\end{eqnarray} \tag{ 4 }$$

with ${\rm{\Omega }}=\sqrt{\tfrac{{{GM}}_{\star }}{{r}^{3}}}$, $\tilde{D}=\tfrac{{H}_{g}\alpha }{{S}_{c}{\rm{\Omega }}}$, and ${\tau }_{S}(a)=\tfrac{{\rho }_{d}a}{{\rho }_{g}{c}_{S}}$, and where c_S = H_g Ω 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), ρ_d is the dust material density, and ρ_g is the gas density in the midplane. We note that the Stokes number described in Section 5.2 corresponds to St = Ωτ_S (a). The degree of settling is set by varying the α parameter. We also note that, when z ≪ 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 α = 10⁻³ − 10⁻² 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 α = 10⁻⁴–10⁻⁶ (Bai & Stone 2013; Flock et al. 2020). In this context, we produced models for different settling strengths, obtained for respective α parameters of 10⁻⁵, 10⁻⁴, and 10⁻³. For this modeling, we assume that the grain size integrated over the whole disk follows a simplified power-law distribution of n(a)da ∝ 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 β = 1.1, but the large grains will have a smaller scale height than the gas because of the settling prescription considered (H_{ld,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 × 10⁻⁵ M_⊙. 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⁻⁵, 10⁻⁴, and 10⁻³ 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 α = 10⁻⁵ 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 α = 10⁻⁵ provides the best match with the observations.

Zoom In Zoom Out Reset image size

Thus, we find that the millimeter observations of Oph 163131 are consistent with a turbulent viscosity coefficient of α ≤ 10⁻⁵, at least in the outer regions of the disk (r ∼ 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 × 10⁻³ 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.

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_{ld,100 au}/H_{g,100 au} ∼ 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 (≳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 (α ≲ 10⁻⁵) 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.

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_⊙. 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.

Zoom In Zoom Out Reset image size

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⁻² 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 ∼ 0.05, consistent with H_{ld,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.