Using Bayesian Deep Learning to Infer Planet Mass from Gaps in Protoplanetary Disks

Observed substructures, like concentric rings and gaps, from near-infrared images of protoplanetary disks (PPDs) in dust emission (ALMA Partnership et al. 2015; Andrews et al. 2016; Huang et al. 2018a) are often interpreted as the signatures of an embedded planet. However, owing to limitations in the current planet search techniques (Fischer et al. 2014; Lee 2018), the direct detection of planets in disks is still a challenge.

Only recently have there been measurements using gas kinematics (Pinte et al. 2018; Teague et al. 2018; Pinte et al. 2019, 2020) that have found planetary signatures, such as "kinks" due to velocity perturbations, inside the observed dust gaps. Additionally, the direct detection using H α emission of an accreting planet (Keppler et al. 2018; Wagner et al. 2018; Haffert et al. 2019) inside the cavity of the transition disk PDS-70 has provided compelling evidence in favor of planetary gaps.

Thus, characterizing the large-scale planet-induced distortions observed in PPDs offers a unique opportunity to probe unseen young planets during the epoch of their formation. This is achieved by identifying and comparing the features with theoretical models of planetary gaps generated using customized simulations (Crida et al. 2006; Paardekooper & Papaloizou 2009; Duffell & Macfadyen 2013; Fung et al. 2014; Duffell 2015; Kanagawa et al. 2015a; Zhang et al. 2018; Ilee et al. 2020) or empirical relations (Kanagawa et al. 2016; Lodato et al. 2019). However, owing to the increasing sample size of the observed disks and the complexity of disk–planet simulations, one needs to run hundreds of customized simulations to analyze each observed system. This is too expensive and inefficient to be of practical use as a detection tool. The alternative is to use deep-learning models, trained with simulation data (generated only once), to detect/characterize exoplanets from the observed data in a more efficient and accurate way.

In Auddy & Lin (2020) and Auddy et al. (2021), hereafter Papers I and II, respectively, we successfully designed deep-learning models using a multilayer perceptron (MLP) and/or a Convolutional Neural Network (CNN) to predict planet mass from the planetary gaps in PPDs. Both models (DPNNet-1.0 and DPNNet-2.0) are end-to-end pipelines that are capable of predicting the mass of an exoplanet from an observed PPD harboring a gap-inducing planet. Recently, Zhang et al. (2022) have also developed a CNN-based model to directly infer the planet mass from radio dust continuum images. However, one of the main challenges with traditional deep-learning models is quantifying the uncertainties associated with their predictions. This has direct implications for how a deep-learning model performs when trained on limited and/or noisy data. Furthermore, the inability of the traditional deep-learning model to answer "I do not know" or "I am not confident" is concerning, particularly when treating out-of-training distribution points (i.e., input data points that are outside the parameter space explored in the training data set). Thus, it is crucial to understand the limitations of a deep-learning model by quantifying the errors associated with its predictions and identifying the source of the uncertainty.

Recently, some of these shortcomings have been addressed by introducing probabilistic methods, such as Bayesian techniques, to deep-learning architectures (Kendall & Gal 2017; Jospin et al. 2020; Wilson & Izmailov 2020). A Bayesian Neural Network (BNN) improves on traditional deep-learning models by providing uncertainty estimates associated with the predictions of these models. BNNs are truly unique in their ability to distinguish between the uncertainty that is caused by the limitations of the model and the uncertainty that is inherent to the input data. More so, these architectures can even respond to pathological out-of-training cases with a large error bar.

For a classical neural network, one would have deterministic point estimates from a single optimized network as the output of the network (e.g., see Paper I). A BNN considers an aggregation of predictions (also known as marginalization) from multiple independent predictor networks and performs a model average to obtain the output as a distribution. A BNN is a type of stochastic neural network, in which the weights are given as a distribution. It enables us to understand the uncertainty that is associated with the network architecture. This is particularly powerful in constraining the error when performing a regression task to estimate parameters (such as planet mass, in our case) using deep-learning models. Thus, any prediction comes with a confidence interval, meeting the requirement of a scientific experiment, and one that can be compared with other methods of parameter estimation.

In this paper, we develop and introduce "DPNNet-Bayesian," a standard MLP-based model designed using Bayesian formalism to predict planet mass from disk gaps as well as to quantify the uncertainty associated with the prediction. Our network is trained with synthetic data from explicit planet–disk simulations using the fargo3d hydrodynamics code, with a newly implemented fixed grain size module and improved initial conditions.

The paper is organized as follows. In Section 2, we give an overview of Bayesian deep learning, its implementation, and its advantages over traditional deep-learning approaches. In Section 3, we describe the disk–planet hydrosimulations as well as the parameter space considered for the current study. We introduce the BBN architecture implemented in DPNNet-Bayesian in Section 4, and discuss in detail the steps used to preprocess the data. In Section 5, we give the predictions of, along with the error estimates from, our trained DPNNet-Bayesian, based on simulated data. In Section 6, as a proof of concept, we deploy DPNNet-Bayesian to observed PPDs around HL Tau and AS 209 to infer the planet masses and the associated errors. Finally, in Section 7, we discuss the results and draw conclusions.

A neural network architecture consists of a few thousands (or even millions) of nodes/neurons that are densely interconnected. For this architecture, let

• 1.  
  θ represent all the internal model parameters, and
• 2.  
  D represent the training data, which consists of a specified set of observed features D_x (e.g., gap size) and their corresponding output labels D_y (e.g., planet mass).

Deep learning is a technique for optimizing the model parameters θ based on a training data set D (Goodfellow et al. 2016; Shrestha & Mahmood 2019). This is done by computing the minimum "cost" using some method of backpropagation. Usually, the cost function is defined in terms of the log likelihood of the training data set, and thus the training process can be defined as a maximum-likelihood estimation. This approach might suffer from overconfidence, as it is deterministic in nature. An alternative to this is to implement the use of stochastic components in the neural network. This can account for the uncertainties in the determination of the network parameters θ, which finally shows up as an estimated error in the predictions made by a deep-learning architecture (Jospin et al. 2020). For these types of neural nets, the stochastic elements can be introduced in terms of stochastic model weights or activation. Thus, instead of a single deterministic output from the traditional deep-learning model, for the stochastic networks a marginalization of the output is considered, based on some statistical principle (Wilson & Izmailov 2020).

In order to design a predictive model based on BNN, one starts with a suitable choice of deep-learning architecture. The next step is to choose a systematic prescription, in order to introduce the stochastic elements into the network architecture. This can be modeled in terms of some possible prior distribution of the model parameters, p(θ).

From Bayes' theorem, we get

$$\begin{eqnarray}&&p(\theta | D)\propto p({D}_{y}| {D}_{x},\theta )p(\theta ),\end{eqnarray} \tag{ 1 }$$

where D_y denotes the training labels and D_x are the training features, as mentioned earlier. The Bayesian posterior distribution, p(θ∣D), depicts the model parameters as a probabilistic distribution, rather than some deterministic value.

The Bayesian approach allows us to compute the model averages by essentially integrating over the models, parameterized by the model parameters θ. The distribution of the possible values of y (e.g., planet mass) given some input x (e.g., dust gap width) is given by

$$\begin{eqnarray}&&p(y| x,D)={\int }_{\theta }p(y| x,{\theta }^{{\prime} })p({\theta }^{{\prime} }| D)d{\theta }^{{\prime} }.\end{eqnarray} \tag{ 2 }$$

In Equation (2), the variance of $p(y| x,{\theta }^{{\prime} })$ captures the uncertainty due to the measurement noise of the input data, while the variance of $p({\theta }^{{\prime} }| D)$ corresponds to the uncertainty associated with the neural network architecture.

The classical training would correspond to p(θ∣D) being replaced by a constant value of θ, but in the Bayesian case this would typically be a probability distribution with some finite variance (Wilson & Izmailov 2020; Abdar et al. 2021). Also, computing p(θ∣D) from the training data is one of the most challenging and difficult aspects of the BNN architectures, due to the large size of the parameter space. This is further discussed in the following sections.

2.1. The Variational Inference Approach to Bayesian Deep Learning

Due to the large dimensionality of the sampling space, it can be quite challenging to merge the concepts of a neural network with the Bayesian approach. As the posterior distribution is extremely complex, directly estimating p(θ∣D) to perform the integral in Equation (2) is not always possible. Different methods, involving various approximations, have been implemented to do this. One of the most popular and direct methods for sampling the posterior distribution, p(θ∣D), is using a Markov Chain Monte Carlo (MCMC) method (Bardenet et al. 2017). However, MCMC often has issues with scalability, owing to the large model sizes, and it can be computationally expensive as well.

Instead of sampling from the posterior directly, the method of variational inference (Blei et al. 2017; Mohan et al. 2022) can be used as an alternative. A distribution (referred to as the variational distribution), parameterized by an arbitrary set of parameters ϕ, is used to approximate the posterior. The parameter ϕ is learned in a way that the distribution q_ϕ (θ) becomes close to the exact posterior p(θ∣D). The closeness between the posterior and q_ϕ (θ) is quantified in terms of the Kullback–Leibler (KL) divergence (Zhang et al. 2019; Jospin et al. 2020). More precisely, the formalism of the evidence lower bound (ELBO) is used for this purpose:

$$\begin{eqnarray}&&\mathrm{ELBO}={\int }_{\theta }{q}_{\phi }({\theta }^{{\prime} })\mathrm{log}\left(\displaystyle \frac{p\left({\theta }^{{\prime} },D\right)}{{q}_{\phi }\left({\theta }^{{\prime} }\right)}\right)d{\theta }^{{\prime} }.\end{eqnarray} \tag{ 3 }$$

Minimizing the ELBO is equivalent to minimizing the KL divergence, as $\mathrm{ELBO}=\mathrm{log}(p(D))-{D}_{\mathrm{KL}}({q}_{\phi }| | p)$, where log (p(D)) depends on the prior and D_KL is the KL divergence. Such a variational approach has been used in a variety of other astronomical applications, such as in Roberts et al. (2013), Perreault Levasseur et al. (2017), Leung & Bovy (2019), and Spindler et al. (2021).

2.2. Features and Advantages of BNN

Once the ELBO function is optimized, and the variational distribution is thus learned/trained, we can use it to predict the labels of new (unseen) observations. The predictive posterior distribution p(y^*∣x^*, D) is given by integrating out the model/variational parameters (trained using the data set D):

$$\begin{eqnarray}&&\begin{array}{l}p({y}^{* }| {x}^{* },D)={\displaystyle \int }_{\theta }p\left({y}^{* }| {x}^{* },{\theta }^{{\prime} }\right)p\left({\theta }^{{\prime} }| D\right)d{\theta }^{{\prime} },\end{array}\end{eqnarray} \tag{ 4 }$$

where (x*, y*) represent the features and labels, respectively, and the superscript "∗" denotes unseen data. In order to obtain this predictive distribution, we use Monte Carlo sampling. First, we sample from the distribution of model parameters N times:

$$\begin{eqnarray}&&{\theta }^{j}\sim p(\theta | D)\,j\in \{1,2,...,N\}.\end{eqnarray} \tag{ 5 }$$

Then x_i is sampled n times from the distribution of the input features, whose width is determined by the errors associated with the feature itself, to estimate the predictive distribution $p({y}_{i}^{* }| {x}_{i}^{* },{\theta }^{j})$. For N × n iterations, we get

$$\begin{eqnarray}&&\begin{array}{l}p\left({y}^{* }| {x}^{* },D\right)\approx \displaystyle \sum _{i=1}^{n}\displaystyle \sum _{j=1}^{N}p\left({y}_{i}^{* }| {x}_{i}^{* },{\theta }^{j}\right)p\left({\theta }^{j}| D\right).\end{array}\end{eqnarray} \tag{ 6 }$$

The error/uncertainty in the predictions of the BNN can be broadly classified into two categories: epistemic and aleatoric uncertainty (Abdar et al. 2021; Hüllermeier & Waegeman 2021). "Aleatoric" uncertainty (AL) represents the uncertainty inherent in the data that is given as the input for the model (e.g., the measurement noise). In contrast, "Epistemic" uncertainty (EP) represents the limitations of the machine-learning (ML)–based model that usually stem from a lack of information about the system. Taking a larger data set for the training of the model (spanning a larger parameter space) or improving the ML architecture can typically reduce the EP. On the other hand, as the AL is intrinsic to the measurement noise of the input data, it cannot be reduced by solely improving the ML model. Both these uncertainties propagate to the output to give the net predictive uncertainty (Abdar et al. 2021), as shown in Figure 1.

Zoom In Zoom Out Reset image size

We can measure the EP, , from the variance of the posterior distribution of the model parameters, p(θ^j ∣D), given as:

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{\mathrm{EP}}^{2}\sim \mathrm{Variance}\ \mathrm{of}\ p({\theta }^{{\prime} }| D).\end{array}\end{eqnarray} \tag{ 7 }$$

σ_EP quantifies the intrinsic uncertainty associated with the model. It enables us to understand the correlation between the model's performance, the quantity/quality of the data, and the parameters of the ML architecture. Thus, monitoring σ_EP can give insights into improving the overall ML model.

The AL, is obtained from the variance of $p({y}_{i}^{* }| {x}_{i}^{* },{\theta }^{j})$, given as:

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{\mathrm{AL}}^{2}\sim \mathrm{Variance}\ \mathrm{of}\ p(y| x,\theta ).\end{array}\end{eqnarray} \tag{ 8 }$$

σ_AL quantifies the statistical uncertainty/error associated with the various measurements that go as input features into the neural net at the training stage, as well as when the model is deployed for real-world application. Knowing σ_AL, one can monitor the errors in the predictions of the model as the measurement noise of the input features changes. The advantage of knowing both σ_EP and σ_AL is that when testing the model, one can trace the source of the uncertainty and improve the model/data set accordingly.

We model the structural variations of a dusty PPD due to an embedded gap-opening planet using the fargo3d hydrodynamic simulation code (Benítez-Llambay & Masset 2016). We update the public version of the code by parameterizing the dust–gas interaction with a constant particle size, instead of a fixed Stokes number. This includes a generalization of the steady-state profiles derived by Benítez-Llambay et al. (2019), which are presented in Appendix. We generate the data set with the updated simulations to train the BBN model. The simulation setup and the parameter space are discussed in the next section.

3.1. PPD Setup

We model a simple 2D razor-thin disk with an embedded planet of mass M_P located at R = R₀. The planet is on a Keplerian orbit around a central star of mass M_* and it is not migrating. The unit of time is the planet's orbital period P₀ = 2π/Ω_K0, where Ω_K0 is the Keplerian frequency. The subscript "0" denotes evaluation at the planet's location R₀.

The disk's initial gas surface density profile is given as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{g}}}(R)={{\rm{\Sigma }}}_{{\rm{g}}0}{\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{-{{\sigma }}_{0}},\end{eqnarray} \tag{ 9 }$$

where σ₀ is the exponent of the surface density profile and Σ_g0 = 10⁻⁴ is a constant in the code unit, which is arbitrary for non-self-gravitating simulations without orbital migration. For reference, however, this surface density scale corresponds to Toomre parameters ${ \mathcal Q }\geqslant 80$ for all of our disk models, which are sufficient for the disk to be gravitationally stable (Toomre 1964). The disk is flared with a constant flaring index F, such that the aspect ratio varies as

$$\begin{eqnarray}&&h=\displaystyle \frac{H}{R}={h}_{0}{\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{F},\end{eqnarray} \tag{ 10 }$$

where H is the disk scale height and h₀ is the aspect ratio at R = R₀. The disk is locally isothermal, with a sound-speed profile of c_s (R) = h₀ R Ω _K. To mimic the disk's turbulence, we include viscous forces in the gas disk using the standard α-prescription (Shakura & Sunyaev 1973), where the kinematic viscosity ν = α c_s ²/Ω _K and α is set by Equation (A21).

We consider a single population of dust particles modeled as a pressureless fluid (Jacquet et al. 2011). The dust particle size is parameterized using the Stokes number S ≡ t_s Ω _K, where t_s is the stopping time characterizing the strength of the dust–gas coupling (Weidenschilling 1977). The initial dust surface density is Σ_d = Σ_g, where is the dust-to-gas ratio set by Equation (A19). The dust's backreaction onto the gas is included in the simulations.

3.2. Simulation Setup and Parameter Space

In this section, we introduce the Bayesian implementation of the DPNNet model.

4.1. Network Architecture

4.2. Data Acquisition and Preprocessing

For the training and testing of the model, we produce the new data set following the same procedure as used in Paper I. We generate a uniform random distribution of the parameters, each centered within the sampling interval, using the Latin hypercube sampling method (McKay et al. 1979; Iman et al. 1981). We run a total of 1800 simulations with a range of input parameters, as described in Section 3.2. From each simulation, we obtain the dust and gas surface density maps of the PPD with the embedded planet. Figure 2 shows sample images of the dust surface density distributions and the corresponding radial profiles of the azimuthally averaged surface densities from some of the simulations. These improved disk–planet simulations (with fixed grain sizes and initial equilibrium disk setups) result in a much more physical model for dust–gas interaction compared to Paper I and Paper II.

Zoom In Zoom Out Reset image size

Following Paper I and Kanagawa et al. (2016), we define the dust gap width as

$$\begin{eqnarray}&&w=\displaystyle \frac{{R}_{{\rm{d}},\mathrm{out}}-{R}_{{\rm{d}},\mathrm{in}}}{{R}_{0}},\end{eqnarray} \tag{ 11 }$$

where R_d,in and R_d,out are the inner and outer edge of the gap where the dust surface density Σ_d equals a predefined threshold fraction of Σ_d0. For consistency with Paper I, we adopt a threshold fraction of half the initial surface. We acknowledge that this definition of w_d has its limitations, as it fails for shallow gaps (for low-mass planets) when the gap depth is less than the threshold function. However, this does not limit the model's application, as DPNNet-Bayesian can be retrained for any threshold function or other definition of the gap width.

In the current paper, we only consider simulations with either one dust gap or two dust gaps. We remove models with >2 dust gaps, as the number of samples is not adequate to effectively train the model. To build our feature space, we simply define two dust gap widths for each simulation, w_d1 and w_d2, but set w_d2 = 0 if only a single dust gap is formed.

The initial data set consists of 1800 simulations. We preprocess the data by eliminating runs that do not open up detectable axisymmetric gaps or that have more than two gaps. After the initial filtering, we have data from 1098 simulations. Each simulation is characterized by its input parameters (M_P, α₀, h₀, σ₀, F, S_t, and ₀) and the associated gap width(s) in dust (w_d1, w_d2) and gas (w_g). In this paper, we will not use the gas gap width as an input parameter for training the model, since this is not directly measured from the observations.

Once the gap widths have been measured, we assign each simulation a set of feature variables (w_d1, w_d2, ₀, α₀, S_t, h₀, σ₀, and F) that will go as the inputs for our model, and we label the planet mass, M_P, as the target variable. The data for each feature variable are then normalized using the standard (Z-score) scaling, by subtracting the mean and scaling it by the standard deviation.

4.3. Training DPNNet-Bayesian

Once the model has been trained, it is ready to be deployed to predict the planet mass from observed disk features. We test the network's performance on the test data set that the model has not been exposed to. This enables us to quantify the network's accuracy in predicting the planet mass from unseen data. The test data set comprises all the feature variables along with the target variable (the true value of the planet mass). The feature variables (w_d1, w_d2, ₀, α₀, S_t, h₀, ${\sigma }_{0}$ , and F) are normalized with the same scaling that the network was trained on. The trained DPNNet-Bayesian network takes the feature variables as inputs and predicts the corresponding planet mass. Unlike the standard MLP model, for the same set of input feature variables, DPNNet-Bayesian gives the predicted output as a distribution.

Figure 3 shows the correlation between the simulated and predicted planet mass in Earth-mass ( M_⊕) units for all the test samples. For each set of input variables, we conduct the experiment (i.e., deploy the trained network) 200 times, to obtain the marginalization of the predicted output. Thus, we obtain a distribution for the predicted planet mass (as demonstrated in Figure 4) capturing the uncertainties associated with the parameters of the trained model. The error bars give the standard deviation associated with each predicted value. The predicted values lie along the black line, indicating a strong correlation with the true planet mass. This is reflected by the r2 score of 0.82 computed over the average prediction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.1. Error Analysis Using BNN

As a first test of our newly improved DPNNet-Bayesian model, we deploy it to estimate the planet masses from the observed dust gaps in the PPDs around HL Tau and AS 209. The results are summarized in Table 1 and discussed below.

Table 1. Inferred Planet Masses from the Gaps in HL Tau and AS 209

				Observed Features					DPNNet-			Other Models
									Bayesian
Name	M*	Rgap	wd1	h0	${{\alpha }}_{0}$	St	${{\epsilon }}_{0}$	σ0	Mp0	Mp1	Mp2	Mp3	Mp4	Mp5	Mp6
	(M⊙)	(au)							(MJ)	(MJ)	(MJ)	(MJ)	(MJ)	(MJ)	(MJ)
	1.00	10	0.81	0.05	10−3	0.005	0.01	1.0	0.27 ± 0.02	1.40	0.20	0.35	0.20	⋯	⋯
HL Tau	1.00	30	0.23	0.07	10−3	0.005	0.01	1.0	0.14 ± 0.01	0.20	0.20	0.17	0.27	⋯	⋯
	1.00	80	0.29	0.10	10−3	0.005	0.01	1.0	0.30±0.02	0.50	0.20	0.26	0.55	⋯	⋯
AS 209	0.83	9	0.42	0.04	10−4	0.016	0.01	1.0	0.06 ± 0.01	⋯	⋯	⋯	⋯	0.37	⋯
	0.83	99	0.31	0.08	10−4	0.016	0.02	1.0	0.12 ± 0.01	⋯	⋯	⋯	⋯	0.18	0.21

Note. M_P0 is the mass (in the Jupiter-mass unit M_J) predicted by the DPNNet-Bayesian model using the input features w_d1, h₀,α₀ , S_t,ε₀ , and σ₀, as given in columns 4–9, and setting w_d2 = 0. M_P1 is the mass predicted by Kanagawa et al. (2016), implementing the empirical relation. M_P2, M_P3, and M_P4 are the masses inferred using customized simulations of HL Tau by Dong et al. (2015), Jin et al. (2016), and Dipierro et al. (2015), respectively. M_P5 is the predicted planet mass of AS 209 by Zhang et al. (2018). M_P6 is the predicted planet mass of AS 209 by Fedele et al. (2018), using specialized simulations.

Download table as:  ASCIITypeset image

6.1. HL Tau

HL Tau is a well-studied system, which shows distinct axisymmetric—possibly planet-induced—gaps in the dust emission (ALMA Partnership et al. 2015). Several studies involving hydrodynamic simulations (Dipierro et al. 2015; Dong et al. 2015; Jin et al. 2016) have suggested that each of the gaps are induced by a single planet. With the assumption that the observed gaps are all planet-induced, we use the properties of the observed gaps as input features for our BNN model to predict the planet mass. DPNNet-Bayesian predicts a distribution for the planet masses for each of the gaps. For the three identified gaps at R = 10 au, 30 au, and 80 au, we obtain the gap widths (w_d1 = 0.81, 0.23, and 0.29) and adopt the disk aspect ratios of h₀ = 0.05, 0.07, and 0.10, respectively, from Kanagawa et al. (2015b, 2016). HL Tau is mostly laminar, with viscosity α₀ = 10⁻³ (Dipierro et al. 2015) and a central star mass of M_* = 1 M_⊙. We set a canonical value of the dust-to-gas ratio ₀ = 10⁻² (ALMA Partnership et al. 2015; Dong et al. 2015) and adopt a mean Stokes number S_t = 0.005 (Dipierro et al. 2015) for all three gaps. Since the uncertainty associated with the input parameters is not well constrained across the literature, we take a conservative approximation and consider a 1%–5% variation in the gap widths only. The top panel of Figure 5 shows the distribution of the predicted masses in the three identified gaps for the given disk parameters and gap widths. The mean values along with the standard deviations of the predicted planet masses at 10 au, 30 au, and 80 au are 86.0 ± 5.5 M_⊕, 43.8 ± 3.3 M_⊕, and 92.2 ± 5.1 M_⊕, respectively. The standard deviation captures the total uncertainty, which includes the uncertainty associated with the BNN architecture (σ_EP is estimated as 4.0 M_⊕, 3.3 M_⊕, and 5.0 M_⊕, respectively, for the three gaps) along with the AL. For the assumed variation of 1%–5% in the gap widths, σ_AL is comparable to σ_EP. However, σ_AL can vary, depending on the noise in the input data. These estimates are mostly consistent with the inferred planet masses from direct numerical simulations like those of Dong et al. (2015), Jin et al. (2016), Dipierro et al. (2015), and DPNNet-1.0 in Paper I.

Zoom In Zoom Out Reset image size

6.2. AS 209

In this paper, we have introduced DPNNet-Bayesian, a BNN-based model that characterizes the masses of unseen exoplanets from observed PPDs. Compared to its predecessor, DPNNet-1.0, the Bayesian approach has the following advantages. First, it quantifies the errors associated with the predictions of the network and traces the source of the uncertainty. Second, it gives a framework to understand various regularization techniques that are used in traditional deep-learning approaches (Jospin et al. 2020). Finally, compared to traditional deep-learning models, where the weights are initialized by some implicit prior, the Bayesian framework accounts for this prior distribution explicitly.

For a BNN, one trains a distribution/ensemble of deep-learning models and then, based on the Bayesian formalism, performs model averaging/marginalization of the output. It is to be noted that the BNN approach is fundamentally different from "deep ensembles," where different initializations are used to train a set of networks, with the outputs being combined/optimized. More so for the BNN approach, a weighted averaging is achieved by taking into account the posteriors of the network parameters as weights. This is different from deep ensembles, where the averaging is directly performed with the output of the network.

Furthermore, the error/uncertainty associated with the output of a BNN model is different from the standard metrics, like the rms error (RMSE) and/or the mean absolute error (MAE) that are used in traditional deep-learning models. Both the MAE and the RMSE are measures of the error (∣M_P,simulation − M_P,predicted∣) when the model is deployed on a validation/testing data set for which the "true" target values are known. However, they do not necessarily quantify the errors accurately (or provide a confidence interval) when applied to an unknown data set for which the target values are not known. This can be problematic, as there is no way of knowing how accurate the prediction is, so one needs to trust the validation/testing error as an indication of the accuracy of the prediction. Since a BBN model works on the principle of marginalization, rather than optimization, the output is obtained by sampling from a posterior distribution, thus providing the necessary confidence interval.

The performance of any ML-based model depends on the quality and/or quantity of the training data set. The Bayesian approach helps us to understand the correlation between the quality/quantity of the data set and the model performance by monitoring the change in model uncertainty. The knowledge of the EP and the AL enables the model to learn from a smaller data set, by regulating the overfitting or underfitting. Thus, the use of BNN architecture plays a pivotal role in improving the training process, as well as determining the optimal deep-learning architecture.

We conclude by highlighting the key results:

• 1.  
  DPNNet-Bayesian estimates the planet mass from the gap width observed in the dust emission. It takes as its input dust properties (abundance, Stokes numbers) and gas disk properties (aspect ratio, surface density profile, viscosity), along with dust gap widths from observations, and predicts the embedded planet mass.
• 2.  
  For the training of DPNNet-Bayesian, we generated a data set from disk–planet simulations using the fargo3d hydrodynamics code, with a newly implemented fixed grain size module and improved initial conditions. The trained model was tested on the simulated data set in order to illustrate its applicability. The predicted mass closely correlated with the "true" value of the simulated planet mass.
• 3.  
  Compared to DPNNet-1.0 (Paper I), the DPNNet-Bayesian architecture additionally estimates the error associated with the prediction of the model. It further distinguishes between the uncertainty pertaining to the deep-learning architecture (EP) and the errors associated with the input variables due to measurement noise (AL).
• 4.  
  We deploy DPNNet-Bayesian to dust gaps observed in the PPDs around HL Tau and AS 209. Our network predicts masses of 86.0 ± 5.5 M_⊕, 43.8 ± 3.3 M_⊕, and 92.2 ± 5.1 M_⊕ for the gaps at 10 au, 30 au, and 80 au, respectively, at HL Tau. Similarly, for the AS 209 disk, we find planet masses of 20.6 ± 2.2 M_⊕ and 36.3 ± 2.5 M_⊕ at the 9 au and 99 au gaps, respectively. These estimates are in close agreement with the results from other studies based on specialized simulations (Dong et al. 2015; Zhang et al. 2018).

We thank the anonymous referee for the constructive comments. S.A. and J.B.S. acknowledge support from NASA under Emerging Worlds through grant 80NSSC20K0702. D.C. acknowledges support from NASA under Emerging Worlds through grant 80NSSC21K0037. M.K.L. is supported by the Ministry of Science and Technology of Taiwan (grants 107-2112-M-001-043-MY3, 110-2112-M-001-034-, and 110-2124-M-002-012-) and an Academia Sinica Career Development Award (AS-CDA-110-M06). Numerical simulations were performed on the TIARA cluster at ASIAA, as well as the TWCC cluster at the National Center for High-performance Computing in Taiwan.

The current public version of the fargo3d code has a multifluid setup, with pressureless dust species interacting with the gas by means of a drag force parameterized by a constant Stokes number. We update this version of the code to incorporate a steady background solution for the radial drift of dust of fixed particle sizes instead of a fixed Stokes number. Following the generalized steady-state solutions from Benítez-Llambay et al. (2019) for a vertically integrated disk with an isothermal equation of state, with pressure P = c_s ²Σ, we rewrite the exact background solutions for the gas and dust (i.e., in the limit of zero coupling) as

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{g}={\beta }(R){{\boldsymbol{v}}}_{K},\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{d}={{\boldsymbol{v}}}_{K},\end{eqnarray} \tag{ A2 }$$

respectively, where ${{\boldsymbol{v}}}_{K}=R{{\rm{\Omega }}}_{K}(r)\hat{{\boldsymbol{\phi }}}$ is the Keplerian velocity, $\hat{{\boldsymbol{\phi }}}$ is the unit vector in the azimuth, and

$$\begin{eqnarray}&&{\beta }=\sqrt{1+2{\eta }(R)},{\eta }=\displaystyle \frac{{h}^{2}}{2}\displaystyle \frac{d{\rm{l}}{\rm{o}}{\rm{g}}P}{d{\rm{l}}{\rm{o}}{\rm{g}}R}.\end{eqnarray} \tag{ A3 }$$

When the dust and gas are coupled, we write their velocities as v _d,g → v _d,g + δ v _d,g . Considering the steady-state axisymmetric Navier–Stokes equations for perturbed velocities (see Equations (63)–(66) in Benítez-Llambay et al. 2019), the gas and dust velocity perturbations in the radial direction are given as

$$\begin{eqnarray}&&{\delta }{v}_{{\rm{g}}{\rm{R}}}(R)=-2{\beta }{{ \mathcal Q }}_{N}{\rm{\Psi }}\left({\beta }-1\right){v}_{{\rm{K}}},\end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray}&&{\delta }{v}_{dR}(R)=\displaystyle \frac{2{S}_{t}}{1+{S}_{t}^{2}}\left({\beta }-1\right){v}_{{\rm{K}}}+\displaystyle \frac{{\delta }{v}_{{\rm{g}}{\rm{R}}}+2{S}_{t}{\delta }{v}_{g{\phi }}}{1+{S}_{t}^{2}},\end{eqnarray} \tag{ A5 }$$

and the azimuthal counterpart is

$$\begin{eqnarray}&&{\delta }{v}_{g{\phi }}(R)=-\left[\left({{ \mathcal S }}_{N}+2{\xi }\right){{ \mathcal S }}_{N}+{{ \mathcal Q }}_{N}^{2}{\rm{\Psi }}\left({\beta }-1\right){v}_{{\rm{K}}}\right],\end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray}&&{\delta }{v}_{d{\phi }}=\displaystyle \frac{1}{1+{S}_{t}^{2}}\left({\beta }-1\right){v}_{{\rm{K}}}+\displaystyle \frac{2{\delta }{v}_{g{\phi }}-{S}_{t}{\delta }{v}_{{\rm{g}}{\rm{R}}}}{2\left(1+{S}_{t}^{2}\right)},\end{eqnarray} \tag{ A7 }$$

where

$$\begin{eqnarray}&&{{ \mathcal S }}_{N}\equiv \displaystyle \frac{\epsilon }{1+{S}_{t}^{2}},{{ \mathcal Q }}_{N}\equiv \displaystyle \frac{\epsilon {S}_{t}}{1+{S}_{t}^{2}},\end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray}&&\begin{array}{c}{\rm{\Psi }}=\left[\left({{ \mathcal S }}_{N}+{\beta }\right)\left({{ \mathcal S }}_{N}+2{\xi }\right)+{{ \mathcal Q }}_{N}^{2}\right]\\ \times \,{\rm{a}}{\rm{n}}{\rm{d}}\,{\xi }\equiv {\beta }\left(\displaystyle \frac{1}{2}+\displaystyle \frac{d{\rm{l}}{\rm{o}}{\rm{g}}P}{d{\rm{l}}{\rm{o}}{\rm{g}}R}\right).\end{array}\end{eqnarray} \tag{ A9 }$$

Expanding Equations (A4) to (A7) using Equations (A8) and (A9), we get

$$\begin{eqnarray}&&{\delta }{v}_{{\rm{g}}{\rm{R}}}(R)=\displaystyle \frac{-2{\beta }\left({\beta }-1\right){v}_{{\rm{K}}}{S}_{t}{\epsilon }}{{{\epsilon }}^{2}+{\epsilon }{\beta }\left(2{\gamma }+3\right)+2{{\beta }}^{2}\left({\gamma }+1\right)+2{{\beta }}^{2}{St}^{2}\left({\gamma }+1\right))},\end{eqnarray} \tag{ A10 }$$

$$\begin{eqnarray}&&{\delta }{v}_{g{\phi }}(R)=-\displaystyle \frac{\left({\beta }-1\right){\epsilon }{v}_{{\rm{K}}}\left(2{\beta }{\gamma }+{\epsilon }+2{\beta }\right)}{{{\epsilon }}^{2}\,+\,2{{\beta }}^{2}{St}^{2}\left({\gamma }+1\right)+{\beta }{\epsilon }\left(2{\gamma }+3\right)+2{{\beta }}^{2}\left({\gamma }+1\right)},\end{eqnarray} \tag{ A11 }$$

$$\begin{eqnarray}&&{\delta }{v}_{dR}\,=\,\displaystyle \frac{4{{\beta }}^{2}\left({\beta }+1\right){S}_{t}{v}_{{\rm{k}}}\left({\gamma }+1\right)}{{{\epsilon }}^{2}+{\epsilon }{\beta }\left(2{\gamma }+3\right)+2{{\beta }}^{2}\left({\gamma }+1\right)+2{{\beta }}^{2}{S}_{t}^{2}\left({\gamma }+1\right)},\end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray}&&\delta {v}_{d\phi }=\displaystyle \frac{{v}_{{\rm{k}}}\left(\beta -1\right)\left(\epsilon \beta +2{\beta }^{2}\left(\gamma +1\right)\right)}{{\epsilon }^{2}+\epsilon \beta \left(2\gamma +3\right)+2{\beta }^{2}\left(\gamma +1\right)+2{\beta }^{2}{S}_{{\rm{t}}}^{2}\left(\gamma +1\right)}.\end{eqnarray} \tag{ A13 }$$

The above solutions, Equations (A10)–(A13), are solutions to the steady-state drift solution only if they satisfy the following continuity equations:

$$\begin{eqnarray}&&{{\rm{\partial }}}_{R}\left(R{{\rm{\Sigma }}}_{g}{\delta }{v}_{{\rm{g}}{\rm{R}}}\right)={{\rm{\partial }}}_{R}\left(R{{\rm{\Sigma }}}_{d}{\delta }{v}_{d}\right)=0.\end{eqnarray} \tag{ A14 }$$

The particle size is parameterized by the Stokes number, S_t = t_s Ω_K, where the stopping time characterizing the strength of the dust–gas coupling is given as

$$\begin{eqnarray}&&{t}_{s}={\rho }_{{\rm{p}}}s/\left({\rho }_{{\rm{g}}}{c}_{{\rm{s}}}\right)=\displaystyle \frac{\pi s{\rho }_{{\rm{p}}}}{2{{\rm{\Sigma }}}_{{\rm{g}}}}.\end{eqnarray} \tag{ A15 }$$

Here, ρ_p is the density of the dust particles, s is the radius of the dust particles, and Σ_g is the gas surface density. For a fixed particle size, s = constant while the Stokes number S_t ∝ Σ_g ⁻¹ is a function of disk radius, since Σ_g ∝ ${R}^{-{{\sigma }}_{0}}$ , unlike the default setup, where S_t is set to constant.

To maintain a constant mass flux, i.e., ${{\rm{\partial }}}_{R}\left(R{{\rm{\Sigma }}}_{g}{\delta }{v}_{{\rm{g}}{\rm{R}}}\right)\,={{\rm{\partial }}}_{R}\left(R{{\rm{\Sigma }}}_{d}{\delta }{v}_{d}\right)=0.$, we set

$$\begin{eqnarray}&&R{{\rm{\Sigma }}}_{g}{\delta }{v}_{{\rm{g}}{\rm{R}}}={\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{t}}{\rm{a}}{\rm{n}}{\rm{t}}.\end{eqnarray} \tag{ A16 }$$

Here, we only consider the gas component, but same result can be arrived at for dust as well. Furthermore, since $R{{\rm{\Sigma }}}_{{\rm{g}}}\propto {\left(\tfrac{R}{{R}_{0}}\right)}^{1-{{\sigma }}_{0}}$, ${\delta }{v}_{{\rm{g}}{\rm{R}}}$ must scale as

$$\begin{eqnarray}&&{\delta }{v}_{{\rm{g}}{\rm{R}}}\propto {\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{{{\sigma }}_{0}-1},\end{eqnarray} \tag{ A17 }$$

such that the net mass flux is constant. Thus, from Equations (A10) and (A17) we can establish

$$\begin{eqnarray}&&{\delta }{v}_{{\rm{g}}{\rm{R}}}(R)={\delta }{v}_{{\rm{g}}{\rm{R}}0}{\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{{{\sigma }}_{0}-1}=\displaystyle \frac{a{\epsilon }}{{{\epsilon }}^{2}+b{\epsilon }+c},\end{eqnarray} \tag{ A18 }$$

where a = −2β (β − 1)v_K S_t , b = β (2γ + 1), c = 2β² (γ + 1) + 2β² St²(γ + 1), and ${\delta }{v}_{{\rm{g}}{\rm{R}}0}=\tfrac{{a}_{0}{{\epsilon }}_{0}}{{{\epsilon }}_{0}^{2}+{b}_{0}{{\epsilon }}_{0}+{c}_{0}}$, with the subscript "0" referring to the values at $R={R}_{0}$. Rearranging Equation (A18), we solve for the dust-to-gas ratio and arrive at

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{-\left({{bd}}_{0}-a\right)\pm \sqrt{({\left({{bd}}_{0}-a\right)}^{2}-4{d}_{0}^{2}c}}{2{d}_{0}},\end{eqnarray} \tag{ A19 }$$

where ${d}_{0}=\tfrac{{a}_{0}{{\epsilon }}_{0}}{{{\epsilon }}_{0}^{2}+{b}_{0}{{\epsilon }}_{0}+{c}_{0}}{\left(\tfrac{R}{{R}_{0}}\right)}^{{{\sigma }}_{0}-1}$ and the negative root gives a physical solution. Thus, in order to satisfy the condition for constant mass flux for a fixed particle size, we adjust the dust-to-gas ratio accordingly.

Furthermore, the viscous velocity of the gas V^2D _vis without dust grains is given as (Lynden-Bell & Pringle 1974)

$$\begin{eqnarray}&&{V}_{\mathrm{vis}}^{2{\rm{D}}}=-\displaystyle \frac{3\nu }{R}\displaystyle \frac{d\,\mathrm{ln}\left(\nu {{\rm{\Sigma }}}_{{\rm{g}}}{R}^{1/2}\right)}{d\,\mathrm{ln}R},\end{eqnarray} \tag{ A20 }$$

where ν = αc_s h is the kinematic viscosity. For the steady-state background solution, i.e., V^2D _vis = 0, we set ν Σ_g R^1/2 = constant in Equation (A20), which is equivalent to

$$\begin{eqnarray}&&\alpha ={\alpha }_{0}{(R/{R}_{0})}^{{\sigma }_{0}-2F-1},\end{eqnarray} \tag{ A21 }$$

where the subscript "0" refers to the values at $R={R}_{0}$ .

We initialize the disk–planet fargo3d simulation with the updated dust-to-gas ratio in Equation (A19) and fix the particle size to a constant value using the fixed particle module. Figure 6 shows the azimuthally averaged initial dust and gas radial surface densities and the final surface density profiles after the system has evolved 1000P₀. As demonstrated, the steady-state profile is maintained with little evolution, in the absence of an embedded planet.

Zoom In Zoom Out Reset image size