Using Deep Neural Networks to compute the mass of forming planets

TL;DR

This paper introduces a deep neural network approach to efficiently predict planetary envelope masses and critical core masses during planet formation, significantly reducing computation time while maintaining high accuracy.

Contribution

The authors developed a large database of internal structure models and trained DNNs to accurately approximate planetary masses, outperforming other machine learning methods and analytical formulas.

Findings

DNN predictions match detailed models within a few percent.

DNN outperforms Random Forest, Gradient Boost, Support Vector Regression.

Analytical formulas can overestimate planetary masses, affecting formation predictions.

Abstract

Computing the mass of planetary envelopes and the critical mass beyond which planets accrete gas in a runaway fashion is important when studying planet formation, in particular for planets up to the Neptune mass range. This computation requires in principle solving a set of differential equations, the internal structure equations, for some boundary conditions (pressure, temperature in the protoplanetary disk where a planet forms, core mass and accretion rate of solids by the planet). Solving these equations in turn proves being time consuming and sometimes numerically unstable. We developed a method to approximate the result of integrating the internal structure equations for a variety of boundary conditions. We compute a set of planet internal structures for a very large number (millions) of boundary conditions, considering two opacities,(ISM and reduced). This database is then used to train Deep Neural Networks in order to predict the critical core mass as well as the mass of planetary envelopes as a function of the boundary conditions. We show that our neural networks provide a very good approximation (at the level of percents) of the result obtained by solving interior structure equations, but with a much smaller required computer time. The difference with the real solution is much smaller than the one obtained using some analytical formulas available in the literature which at best only provide the correct order of magnitude. We compare the results of the DNN with other popular machine learning methods (Random Forest, Gradient Boost, Support Vector Regression) and show that the DNN outperforms these methods by a factor of at least two. We show that some analytical formulas that can be found in various papers can severely overestimate the mass of planets, therefore predicting the formation of planets in the Jupiter-mass regime instead of the Neptune-mass regime.

Figures

Tables

Table 1. Table 1: Validation set accuracy for different machine learning methods to predict the critical mass. All DNN lines refer to fully connected Deep Neural Network, the numbers giving, respectively, the number of hidden layers, and the number of units per layer.

baseline DNN ($3\times 128$ )	1.67 %
linear regression	30.5 %
Random Forest (10 trees)	10.5 %
Random Forest (100 trees)	8.61 %
Random Forest (1000 trees)	8.18 %
Random Forest (10000 trees)	8.11 %
Support Vector Regression	12.5 %
Gradient Boosting	13.5 %
best Gradient Boosting after random search	6.82 %
DNN $2\times 128$	1.87 %
DNN $1\times 128$	13.1 %
DNN $3\times 64$	2.31 %
DNN $3\times 32$	3.49 %

Equations

$$t_{\mathrm{gas}}=10^b{\left(\frac{M_{\mathrm{planet}}}{{\mathrm{M}}_{\oplus }}\right)}^{-c}\left(\frac{\kappa }{1\mathrm{c}{\mathrm{m}}^2/g}\right),$$

$$M_{\mathrm{crit}}=10{\left(\frac{{\dot{M}}_{\mathrm{core}}}{10^{-6}{\mathrm{M}}_{\oplus }{\mathrm{yr}}^{-1}}\right)}^{0.2-0.3}{\left(\frac{\kappa }{1\mathrm{c}{\mathrm{m}}^2/g}\right)}^{0.2-0.3}{\mathrm{M}}_{\oplus },$$

$${\dot{M}}_{\mathrm{gas}}$$

$${\dot{M}}_{\mathrm{gas}}=\frac{M_{\mathrm{planet}}}{t_{\mathrm{KH}}},$$

$$\frac{dP}{dv}=$$

$$\frac{d{M}_{\text{r}}}{dv}=$$

$$\frac{dT}{dv}=$$

$$\nabla_{\text{ad}}=\bigg{(}\frac{\partial\text{ln}T}{\partial\text{\text{ln}}P}\bigg{)}_{\text{S}}$$

$${\nabla }_{\text{rad}}=\frac{3\kappa LP}{64\pi \sigma G{M}_{\text{r}}{T}^4},$$

$$R_{\mathrm{P}}=\frac{G{M}_{\mathrm{planet}}}{\left(C_{\mathrm{S}}^2+4G{M}_{\mathrm{planet}}/R_{\mathrm{H}}\right)}$$

$$\kappa =\max (f_{\mathrm{dust}}{\kappa }_{\mathrm{dust}},{\kappa }_{\mathrm{gas}})$$

$$ReLU(x)=\max (x,0)$$

$$C(w)=\frac{1}{N}\sum _{k=1}^N{\left(M_{\mathrm{critical},\mathrm{predicted},\mathrm{k}}-M_{\mathrm{critical},\mathrm{target},\mathrm{k}}\right)}^2$$

$$28{\left(\frac{a}{AU}\right)}^{-1/2}\mathrm{K}<T<2800{\left(\frac{a}{AU}\right)}^{-1/2}\mathrm{K}$$

$$3.65{\left(\frac{a}{AU}\right)}^{-3.25}{dyn/cm}^2<P<365{\left(\frac{a}{AU}\right)}^{-3.25}{dyn/cm}^2$$

Full text

*Aims. *The aim is to provide a way to approximate the result of integrating the internal structure equations for a variety of boundary conditions.

*Methods. *We compute a set of planet internal structures for a very large number (millions) of boundary conditions, considering two opacities, namely the ISM one and a reduced one. This database is then used to train Deep Neural Networks in order to predict the critical core mass as well as the mass of planetary envelopes as a function of the boundary conditions.

*Results. *We show that our neural networks provide a very good approximation (at the level of percents) of the result obtained by solving interior structure equations, but with a much smaller required computer time. The difference with the real solution is much smaller than the one obtained using some analytical formulas available in the literature which at best only provide the correct order of magnitude. We compare the results of the DNN with other popular machine learning methods (Random Forest, Gradient Boost, Support Vector Regression) and show that the DNN outperforms these methods by a factor of at least two.

*Conclusions. *We show that some analytical formulas that can be found in various papers can severely overestimate the mass of planets, therefore predicting the formation of planets in the Jupiter-mass regime instead of the Neptune-mass regime. The python tools that we provide allow to compute the critical mass and the mass of planetary envelopes in a variety of cases, without having to solve the internal structure equations. These tools can easily replace the aforementioned analytical formulas, while providing much more accurate results.

††offprints: Y. Alibert

1 Introduction

Understanding the formation of planets (from terrestrial to gas giants) requires the development of theoretical models whose outcome can be compared with observations (e.g, Mordasini et al. 2009; Alibert 2019). One primary outcome of planet formation models is the mass and internal structure (in particular mean density) of planets, as these two quantities can be obtained from radial velocity and transit observations of extrasolar planets.

In the core accretion model (e.g., Bodenheimer & Pollack 1986; Pollack et al. 1996; Hubickyj et al. 2005), the mass growth of planets results from two phenomena. The first is the accretion of solids (planetesimals or pebbles, e.g Alibert et al. 2013; Ormel 2017), which depends on the properties of the accreted solids, the planetary mass, and the disc thermodynamical properties (Fortier et al. 2013; Ormel & Klahr 2010). The second is the accretion of gas which occurs in two regimes (see Helled et al. 2014, for a review). In the first regime, when the planet is small, the planetary gas envelope extends until the protoplanetary disc, at a radius that is generally between the Bondi and Hill radius (Lissauer et al. 2009). The planet is said to be attached to the disc, and the accretion rate of gas depends on the cooling efficiency of the planet. Indeed, planets must cool in order to accrete gas (Lee & Chiang 2015), and the mass of the gas envelope in this regime must be computed by solving the equations for planetary interior structure (Ikoma et al. 2000; Venturini et al. 2015). When the mass of the planetary core is larger than the so-called critical mass, the radiative loss at the surface of the planet cannot be compensated for by the accretion energy of solids (Mizuno 1980), and the planet starts to contract, in a runaway fashion111This does not imply that the accretion rate of gas is large, but that it is accelerating. (Bodenheimer & Pollack 1986; Ikoma et al. 2000)

For sub-critical planets (when the core mass is smaller than the critical mass), the computation of the planetary envelope requires solving a set of differential equations (see Sect. 2). Although standard procedure, solving these internal structure equations can require a non-negligible amount of computer time, and lead in some cases (e.g. close to the critical core) to some numerical instability. For this reason, some authors have developed fitting formula which allows estimating the envelope mass as a function of different parameters (e.g. opacity, core mass, etc.), as well as some analytical approximation of the critical core mass. For example, Ikoma et al. (2000, hereafter I00), showed that the gas accretion rate is given by:

[TABLE]

where $b=8$ and $c=2.5$. $M_{\rm planet}$ is the planetary mass and $\kappa$ the envelope opacity. This formula was derived in the absence of solid accretion, but the authors showed that the accretion timescale in the presence of a constant solid accretion follows a similar law, the pre-factor being multiplied by six. Certainly, the accretion energy of solids provides some thermal support that slows down the accretion of gas.

Using the same calculations, Ida & Lin (2004, hereafter IL04) showed that the critical core mass can be approximated as:

[TABLE]

where $\dot{M}_{\rm core}$ is the accretion rate of solids.

More recently, based on the results from Piso & Youdin (2014), Bitsch et al. (2015, hereafter B15) used an accretion rate of gas given by:

[TABLE]

where $\kappa$ is the opacity in the envelope, $M_{\rm core}$ and $M_{\rm env}$ are the core and envelope mass respectively, $\rho_{c}$ is the core density, $T$ the temperature at the planet’s location in the disk, and $f$ is a numerical factor of the order of 0.2. This formula is in principle valid in the absence of any solid accretion, but following the argument of I00, this accretion rate can be used in the case of a constant sold accretion rate, by reducing factor $f$ from $0.2$ down to $\sim 0.033$.

For super-critical planets, the accretion rate of gas can be approximated using the Kelvin-Helmoltz timescale $t_{\rm KH}\sim t_{\rm gas}$ as:

[TABLE]

$t_{\rm KH}$ itself depending on the planetary mass (see Eq. 1, I00, IL04).

The use of these simplified formulas, although convenient from the numerical point of view, is questionable when computing the formation of low mass planets, in the Super-Earth up to the Neptune mass range. Indeed, as the envelope mass of these planets is small (smaller than the core mass), an error of a couple of Earth masses can notably change their radius and other properties like their habitability (Alibert 2014). This was for example demonstrated in the context of planet formation by pebble accretion by Brügger et al. (2018), who compared the resulting mass function obtained using Eq. 3 on one side, and computed solving internal structure equations on the other side. These authors showed that the resulting planetary mass was very different in the two cases. When using the approximate formula (Eq. 3), many planets were in the gas giant regime (Jupiter mass and beyond), whereas when solving properly the internal structure equations, the planet mass was in the Neptune regime.

In the present paper, we focus on these sub-critical planets, and we use deep learning techniques to compute the critical core mass as well as the envelope mass as a function of the relevant parameters. Deep learning, specifically Deep Neural Networks (DNNs - see Goodfellow et al. 2016) is a sub-branch of machine learning, which has recently been applied to different problems in planetary sciences, from exoplanet search by transit or direct imaging method (Pearson et al. 2018; Gomez Gonzalez et al. 2018), to dynamical analysis (Tamayo et al. 2016; Smirnov & Markov 2017; Lam & Kipping 2018), to spectrum analysis (Márquez-Neila et al. 2018). The advantage of using DNNs is two fold. First, the resulting critical core mass and envelope mass are close to the ‘real’ one (meaning computed by solving the internal structure differential equations). Second, the computer time required to compute these masses is orders of magnitude smaller.

The plan of the paper is as follows. We first present the models we use in Sect. 2, in particular regarding the treatment of the micro-physics and boundary conditions. We then present the architecture of the DNNs we use, as well as the data we use to train it (Sect. 3). Finally, in Sect. 4 we compare the results of the DNN with the ones obtained by solving the internal structure equations. All the tools developed to use the DNNs are provided as Jupyter notebooks at https://github.com/yalibert/DNN_internal_structure/ for an easy implementation in python codes.

2 Internal structure models

2.1 Planetary growth

In the core accretion model, a planet is formed first by the accumulation of heavy elements in a central core, followed by gas accretion, which builds an envelope on top of the core. The boundary between the protoplanet and the disc is in principle arbitrary, usually taken as a mean between the Bondi and the Hill radius (Lissauer et al. 2009). Between two time steps the envelope mass grows by two mechanisms:

(i) The increase in core mass allows to bind, in hydrostatic equilibrium, a more massive envelope around it.

(ii) The envelope radiates energy into the disk, cooling, and therefore, contracting. This allows for more gas to fill the radius of the protoplanet.

The standard procedure to compute the growth of a planet is to solve the internal structure equations (Sect.2.2) from the top of the envelope to the core-envelope boundary. The energy radiated away in each time step is calculated from energy conservation. This provides the updated luminosity $L$ required to solve the equations (see Venturini et al. 2016).

2.2 Equations

In order to generate the data used to train and validate our DNNs, we solve first the internal structure equations (Kippenhahn et al. 2013) for a large set of boundary conditions. For this, the protoplanet is modelled as composed by an inert core of mean density $\rho_{\rm core}$ = 5 g/cm3, surrounded by a gaseous envelope with a certain composition (see below). For the structure equations, we use $v=r^{3}/M_{\text{r}}$ as the independent variable, where $r$ is the radial coordinate of the planet, and $M_{\text{r}}$ the mass enclosed inside a sphere of radius $r$. This variable is convenient to stop the integration of the equations at the exact core density and avoid interpolations. Thus, the equations are written as (Venturini et al. 2015):

[TABLE]

where $P$ is the pressure, $\rho$ the density, $T$ the temperature, and $G$ the gravitational constant. Following the Schwarzschild criterion, the gradient in Eq.(5) is: $\nabla={\rm min}(\nabla_{\rm conv},\nabla_{\rm rad})$, where

[TABLE]

is the adiabatic gradient, and the radiative gradient is given by:

[TABLE]

$\sigma$ being the Stefan-Boltzmann constant, $\kappa$ the opacity and $L$ the luminosity.

The luminosity is assumed uniform throughout the envelope, and results from the energy released by the accretion of solids and the contraction of the envelope (see detailed explanation in Venturini et al. 2016).

To close the system of structure equations, an equation of state (EOS) must be provided. The EOS gives $\rho$ and $\nabla_{\text{ad}}$ as a function of $T$, $P$ and composition.

2.3 Boundary conditions

The outer boundary conditions to solve Eqs. 5 are set by the temperature ($T_{\rm out}$) and pressure ($P_{\rm out}$) at the planetary radius, defined as a combination of the Bondi and Hill radius (Lissauer et al. 2009):

[TABLE]

where $C_{\rm S}^{2}$ is the square of the sound velocity in the protoplanetary disk at the planet’s location $a$, and $R_{\rm H}=a_{\rm p}\left(\frac{M_{\rm planet}}{3M_{\odot}}\right)^{1/3}$.

The inner boundary conditions are given by the planetary core mass and core radius. When solving these differential equations in planet formation models (see e.g. Alibert et al. 2005), the total mass of the planet is iteratively adjusted until all the boundary conditions are fulfilled.

2.4 Microphysics

In this work we assume the envelope to be composed by a mixture of hydrogen and helium, with a helium mass fraction of $Y=0.28$. We do not consider the effect of pollution by accreted solids (Venturini et al. 2015; Venturini & Helled 2017). Such an effect will be included in a future work. The density and adiabatic gradient of the internal structure equations are determined by $T$, $P$ and composition, through the equation of state (EOS). Since we assume a H-He composition for the envelope, we adopt the EOS of Saumon et al. (1995). For temperatures and pressures below the limits of the Saumon et al. (1995) table, the EOS of ideal gas for the mixture of H2 and He is assumed, accounting for the corresponding degrees of freedom for the entropy and energy of the mixture.

The opacities are defined as:

[TABLE]

where the dust opacities are taken from Bell & Lin (1994) (hereafter BL94), and the gas opacities from Freedman et al. (2014). The factor $f_{\rm dust}$ accounts for the possibility of the dust opacities to have reduced values compared to the ISM. This can result due to grain growth and settling (Movshovitz et al. 2010; Mordasini 2014). In this work we present calculations for $f_{\rm dust}=1$ and $f_{\rm dust}=0.01$. Since the dust opacity values peak at the top of the envelope and the gas ones at the deeper, hotter regions, the definition in Eq.(9) is equivalent to considering $\kappa=f_{\rm dust}\,\kappa_{\rm dust}+\kappa_{\rm gas}$, as shown in Mordasini (2014).

The domain of validity of the Freedman gas opacities in pressure and temperature is, respectively, 1 to 300 bars and 75 K to 4000 K. At the core-envelope boundary, higher values of temperature and pressure typically exist, but in that domain of high temperature and pressure the envelopes are always convective. Therefore, even taking a constant extrapolation in these region does not affect the internal structure computation. An analytical fit to the Freedman gas opacities is given in Freedman et al. (2014). This is implemented to reduce computational time.

3 Deep Neural Network

3.1 Principle

In the last decades, machine learning has seen a tremendous development, thanks to the advancement of novel techniques, and to the availability of very large amount of data. In particular, neural networks have been more and more used during the last decade to achieve very good performances in tasks such as image classification and text translation. This method has also been applied to regression, which is linked to the prediction of a continuous variable (see Goodfellow et al. 2016), and is the case of interest in this paper. In regression, the task of a neural network is to use some input variables (like temperature, pressure at the planetary outer radius, etc.) to predict one or more variables (the mass of the gas envelope in our case). In the machine learning literature, the input variables are referred as features and the predicted variable as target variable.

A neural network is made of different layers of units, each unit performing a simple task: linear combination of its inputs, followed by a non-linearity (see below). Each layer of a neural network takes typically as input the outputs of the previous layer to produce a series of outputs that will be processed by the next layer. The first layer takes as input the features, and the last layer outputs a prediction of the target variable that is, hopefully, close to the actual target variable. A neural network is said to be deep (named DNN for Deep Neural Network) when it has at least one hidden layer (i.e., a layer that is not the first or the last one). A DNN is trained by providing a (large) number of data (features and target variable), and adjusting the coefficients (also named weights) of all linear combinations of all units in an iterative way, in order to achieve the best match between the predicted target value and the actual one. The matching between the predicted and actual values is quantified using a cost function (see below), which is a function of all the parameters of the DNN. The gradient of the cost function is then computed (as a function of the parameters of the DNN) using a method called backpropagation (Goodfellow et al. 2016), which consists basically in the application of the chain rule for derivation. Finally, gradient descent is used in order to find the optimum parameters. Note that the cost function is in general a non-convex function of the DNN parameters, and many local minima can exist. There is therefore no guarantee that the parameters found are the ones leading to the absolute best performances of the DNN.

We have considered two DNNs. The first one, DNN1, to predict the critical core mass (as a function of four parameters, see next section). The second one, DNN2, to predict the gas envelope mass, as a function of five parameters (the same four parameters plus the core mass). The architecture of the DNNs for the different values of $f_{\rm dust}$ is the same, but the resulting weights are of course different.

DNN1 has 3 hidden layers, each having 128 units. These numbers were found by a series of trial and error tests where we varied the number of layers and the number of units per layer in a grid search (see Sect.5). We also compare in Sect. 5 the results obtained using other machine learning methods. The DNN is fully connected, meaning that each unit of a layer is connected to each layer of the previous and next layer. For each unit, we chose to use the ReLU function for the non-linearity, given by:

[TABLE]

The cost function we minimise is the mean square error between the predicted and the actual critical mass:

[TABLE]

where $N$ is the total number of data points we use, $M_{\rm critical,target,k}$ is the critical core mass computed using the internal structure differential equations (see previous section), $M_{\rm critical,predicted,k}$ is the value predicted by the DNN, and $w$ are the weights.

DNN2 has five hidden layers, each having 128 units, the non-linearity is also the $ReLU$ function, and the cost function is similar, except that we compare the predicted envelope mass and the real one.

Both DNN were trained using the ADAM optimiser (Geron 2017) for 10000 training epochs. The learning rate is chosen equal to 0.0001 for the training. Larger learning rate could be chosen, but this would translate to noisy decrease of the cost function. The other parameters of the ADAM optimiser are taken as the default values in TensorFlow (see https://www.tensorflow.org/api_docs/python/tf/train/AdamOptimizer). We did not used mini-batch since we could fit the whole training dataset (of size 25000 and 800000 respectively for DNN1 and DNN2) in the memory at once. The sizes of the test and validation sets are both equal to 1000 for DNN1 and 300000 for DNN2. Finally, the bias for all units was initialised at 0, whereas the weights were randomly initialised using a gaussian distribution centered on 0, with a variance scaling with the square root of the number of input and output connections (Geron 2017).

3.2 Data

To generate the data, we selected, in a first step, $\sim$10000 points in a four dimension space, by drawing at random a semi-major axis $a$ (uniform in log between 0.1 and 30 AU), a pressure $P_{\rm out}$ (uniform in log between $10^{-2}$ and $1500$ dyn/cm2.), a temperature $T_{\rm out}$ (uniform between $30$ K and $1500$ K), and a luminosity $L$ (uniform in log between $10^{22}$ erg/s and $10^{28}$ erg/s). These luminosities encompass the energy per unit time that would be given by an accretion rate of solids of $\dot{M}_{Z}=10^{-10}\,\rm{M}_{\oplus}$/yr onto a core of 1 $\rm{M}_{\oplus}$, and $\dot{M}_{Z}=10^{-5}\,\rm{M}_{\oplus}$ /yr onto a core of 20 $\rm{M}_{\oplus}$.

For each set of $(a,L,P_{\rm out},T_{\rm out})$, we compute, by solving the differential equations presented in Sect. 2, the core mass for a given planetary mass. $P_{\rm out}$ and $T_{\rm out}$ are used as the outer boundary conditions (Sect. 2.3). For a planetary mass growing from $M_{\rm planet}$ = 1 $\rm{M}_{\oplus}$ onwards, the core mass first increases, then reaches a local maximum and then decreases. The critical core mass is simply the first local maximum core mass obtained in such a way (Mizuno 1980; Papaloizou & Terquem 1999).

For each of these $(a,L,P_{\rm out},T_{\rm out})$, we compute a full $M_{\rm planet}-M_{\rm core}$ curve, for a planetary mass ranging from 1 $\rm{M}_{\oplus}$ up to the corresponding critical mass. Each of these curves is sampled with a step of 0.01 $\rm{M}_{\oplus}$ in $M_{\rm planet}$, and the resulting $M_{\rm envelope}$ as a function of the five parameters $(a,L,P_{\rm out},T_{\rm out},M_{\rm core}$) is finally computed. In addition to the relation between the envelope mass and the five input parameters $(a,L,P_{\rm out},T_{\rm out},M_{\rm core})$, these curves provide us with the critical core mass as a function of the four first parameters $(a,L,P_{\rm out},T_{\rm out})$.

Importantly, since the critical core mass depends on the $(a,L,P_{\rm out},T_{\rm out})$ parameters (it is e.g., smaller for smaller $L$), the number of points on a $M_{\rm planet}-M_{\rm core}$ is not the same for each $(a,L,P_{\rm out},T_{\rm out})$ choice. Since the cost function $C(w)$ defined above for DNN2 gives the same importance to each data points, taking all the points of all the $M_{\rm planet}-M_{\rm core}$ curves would result in increased performances for large critical core masses than for small ones. Indeed, data points for which $M_{\rm critical}(a,L,P,_{\rm out}T_{\rm out})$ is large have more importance in the cost function, since they are more numerous. Such an approach could lead to a DNN predicting accurately the envelope masses for, say, large $L$ (large critical cores) and very poorly for low $L$ (small critical cores). To avoid this effect, and have a similar prediction accuracy for all parts of the $(a,L,P_{\rm out},T_{\rm out},M_{\rm core})$ space, we use the same number of points for each of the $M_{\rm planet}-M_{\rm core}$ curves (this technique is known as stratification). Finally, before implemented, all the data are pre-processed by removing the mean and scaling to unit variance each in each dimension separately.

We obtain $\sim$27000 data points to be used with DNN1 (for predicting the critical core) for each choice of the opacity parameter $f_{\rm dust}$, and $\sim 10$ millions data points to be used to train DNN2 (to predict the envelope mass), again for each choice of the opacity parameter $f_{\rm dust}$. Fig. 1 shows the distribution of $(a,L,P_{\rm out},T_{\rm out})$ for a subset of the datapoints.

These data points are split randomly in a training set, a validation set and a test set. The training set is used to adjust the weights of the DNNs, the validation set is used to estimate the generalisation performances of the DNNs, and avoid, in particular, overfitting. Once the optimum weights of the DNNs are computed (in practical when the cost function on the validation set is the smallest), we use the test set to present the results shown in the rest of the paper (the data points used to generate the different plots shown below are therefore never used in any part of the optimisation process of the DNNs). It is important to check that the distribution of the parameters is the same for the three sets (training, test and validation). This should be the case, as the sets were generated randomly, and this can be checked looking at the cumulative distribution of the input parameters for the three sets as shown in Fig. 2.

The $(a,L,P_{\rm out},T_{\rm out},M_{\rm core})$ parameters are drawn at random independently in our data generation procedure. It is clear, on the other hand, that some of these parameters are strongly correlated, for exemple the pressure and the temperature are a function of the semi-major axis in actual proto-planetary discs. This raises two questions: 1- what kind of internal structure could unrealistic input parameters (e.g. high temperature and low pressure) lead to, and 2- what is the effect of the these unrealistic data on the global optimisation process of the DNN. Indeed, the cost function as presented above takes into account all considered data with an equal importance. Taking into account unrealistic data could therefore degrade the performances of the whole DNN.

Regarding the first question, we have checked that our internal structure code converges even for unrealistic cases. As an example, Fig. 3 shows the internal structure (envelope mass as a function of the radius) for two cases. The planetary core mass is equal to 2 $\rm{M}_{\oplus}$, the semi-major axis equal to 0.5 AU, and the luminosity equal to $10^{25}$ erg/s, which corresponds to an accretion rate of $\sim 5.5\times 10^{-8}\rm{M}_{\oplus}/$yr . The temperature is equal to 400K and the pressure is equal to 10 dyn/cm2 in the first (unrealistic) case, and equal to 350 dyn/cm2, in the second (more realistic) case. These values are typical for a protoplanetary disc (see Venturini & Helled 2017). The envelope mass in the unrealistic case is equal to 0.261 $\rm{M}_{\oplus}$, to be compared to 0.437 $\rm{M}_{\oplus}$ in the realistic case.

We address the second question (the potential effect of unrealistic data on the DNN performances) in Sect. 5 below.

4 Results

4.1 Critical mass

The DNNs are trained using the Tensorflow library (https://www.tensorflow.org), the notebooks to use the DNN are available at https://github.com/yalibert/DNN_internal_structure/.

As mentioned above, the DNNs are trained using only a sub-sample of the generated data, the other part is split between a test set, and a validation set. We use the test set to compute the difference between the predicted critical core mass and the one derived using internal structure equations. The correlation between the two masses is shown in Fig. 4, our results are plotted in blue whereas the critical mass from IL04 is shown in red (Eq. 2). We emphasise the fact that, in the resolution of the internal structure equations, the envelope opacity is not constant but rather a function of the local temperature and pressure, and varies between $\sim\!\!1$ and $\sim\!\!10$ cm2/g in the radiative layers of the atmosphere ($T\lesssim 1000$ K). Hence, we use these two limiting values to compare our results with the ones of IL04, and plot the resulting critical mass as a vertical line in Fig. 4. As can be observed, the critical mass computed using Eq. 2 has a much larger variance compared to the results obtained using the DNN, and is general under-estimated. Using our DNN, the average difference between the predicted critical core mass and the one obtained by solving the internal structure equations is $\sim\!\!1.67\%$, with a majority ($99.3\%$) of the critical masses being predicted with an accuracy better than 5%, and a maximum error of $\sim\!\!14.2\%$ for critical masses larger than 1 $\rm{M}_{\oplus}$.

We perform the same comparison on the right panel of Fig. 4, now assuming that the opacity is reduced by a factor 100 compared to the one of BL94. When using the IL04 formula, the opacity ranges from 0.01 to 0.1 cm2/g. In this case, the critical core mass obtained using Eq. 2 is substantially biased towards smaller values.

4.2 Envelope mass

Similar to the critical mass prediction, we use the test set of envelope masses to compute the difference between the predicted envelope mass and the one derived using internal structure equations. The correlation between the two masses is shown in Fig. 5, where we also show the result estimated using Eqs. 1 (I00) and 3 (B15). In these two latter cases, we need to specify the timescale over which gas is accreted, as these two formulas only provide a gas accretion rate. For the simple estimates provided in Fig. 5, we have assumed that gas is accreted on a timescale that is uniformly distributed between 0 and 5 Myr. This timescale should be equal to the disk lifetime. We emphasis the fact that, in a real situation, the core mass of the planet could increase during this timescale, therefore modifying the accretion rate of gas. In the estimates provided in the figures, this effect is not included.

Comparing the three different cases, we see that the variance is much larger when using the B15 and I00 formulas. Some of this variance of course results from the assumption of the distribution of the time during which gas is accreted. We remark that, as for the case of the critical mass, it is not obvious which value of the envelope opacity should be used in order to obtain the most meaningful comparison between the DNN and the fitting formulas of I00 and B15.

Using our DNN, the masses are predicted with an accuracy better than 20% in more than $97\%$ of the cases.

4.3 Planet formation tracks

In order to compare in a more realistic way the result of the DNN and the one coming from the resolution of internal structure equations, we present in Fig. 6 planetary growth tracks for different cases. We consider different formation locations (5 and 10 AU), different accretion rates ($\dot{M}_{\rm core}=10^{-6}\rm{M}_{\oplus}/$ yr and $\dot{M}_{\rm core}=5\times 10^{-7}\rm{M}_{\oplus}/$), different temperatures (50 and 150 K) and two dust opacities. For each case, we also compare the evolution of the envelope mass with the one that would be obtained with the approach of I00 (Eq. 1) and B15 (Eq. 3). For these two latter cases, we consider in each case two possible values of the envelope opacity, namely $1$ and $10$ g/cm2 when comparing with the BL94 opacity, and $0.01$ and $0.1$ g/cm2 for the case where the BL94 opacity is reduced by a factor 100.

In all the cases, it is clear that the results of the DNNs are much closer to the results obtained by solving the internal structure equations, than by implementing the I00 and B15 formulas. It is also notable that the envelope mass is extremely over-estimated by Eq. 3. In all the cases presented here, the planet ends up as a Neptune mass planet when using the DNN (and therefore also when solving the internal structure equations), whereas it ends-up in the Jupiter mass range using Eq. 3. This result is compatible with the findings of Brügger et al. (2018). The results obtained by the formulas of I00 in the case of BL94 opacity are similar but still somewhat larger than the one obtained from the internal structure equations. In the case of the reduced opacity, the corresponding results (using Eq. 1) are rather different from what is obtained by solving internal structure equations. We finally note that the difference between the DNNs and the internal structure equations seems somewhat larger in the case of the low temperature. This could be solved by using DNNs with larger layers and/or number of units per layer. The difference is however much smaller than the one obtained using either I00 or B15.

5 Discussion

5.1 Overfitting

The very high accuracy that we obtain may seem surprising, and it is important to check that the network is not overfitting the training data. In order to check this, we plot in Fig. 7 the training and validation accuracy as a function of the training epoch for the two DNNs that we have considered. As can be seen, the validation accuracy decreases and starts to plateau after 8000-10000 epochs. On the contrary, the training accuracy still decrease, indicating that the network start to overfit the data at this point. For that reason, the results presented in this paper are the ones obtained after 10000 epochs. This number is high, but results from the choice of a small learning rate in the ADAM optimiser.

5.2 Effect of the dataset

As mentioned above, our input parameters are generated randomly in an independent way. In reality, one expects that, for exemple, the temperature and pressure are strongly correlated with the semi-major axis. As discussed previously, including such unrealistic data-points could degrade the performances of the DNN. In order to test this, we have considered two cases where we trained DNN1 over 10000 epochs, with the same initialisation of the weights. In the first case, we have considered all the data we computed (including the unrealistic ones), whereas in the second case, we have only considered data-points fulfilling the two following conditions:

[TABLE]

and

[TABLE]

These two conditions were chosen allowing for a one order of magnitude variation of the pressure and temperature below and above the values used in Venturini & Helled (2017).

We then compared the output of the two trained DNN on the same set of validation data-points (304 points). In order to have a fair comparison, all validation data-points fulfil the above conditions and are therefore ‘realistic’, and both DNNs were trained using the same number of training data-points (5738). We obtained an accuracy of 4.86 % in the case of the full dataset, and 4.38 % in the case of the reduced dataset. Using only ’realistic’ data-points indeed improves the performances of the DNN as expected, although the change is not massive. The accuracies reached in these two cases are not as good as in the previous section, since both DNNs are trained using a reduced number of points, which degrades their performances.

5.3 Comparison with other machine learning methods.

We compared our results with the ones obtained using other machine learning methods:

• •

linear regression

• •

random forest with 10,100, 1000 and 10000 trees

• •

gradient boosting (using the XGBoost library - see https://xgboost.readthedocs.io/en/latest/)

• •

Support Vector Regression

• •

different architectures of the DNN

For each of these methods, we used the logarithm of the input quantities (except the temperature), and the temperature, as features.

All these methods are described in the literature (e.g. Geron 2017). For each of these methods, a number a hyper-parameters can be tuned in order to improve the performances. For all the methods we considered, we kept the default parameters as defined in scikit-learn (see https://scikit-learn.org/stable/), except when explicitly mentioned. An exhaustive evaluation of the influence of all these hyper-parameters is out of the scope of this paper.

We have in addition, in the case of gradient boosting, performed a random grid search in order to infer if some hyper-parameters could lead to better results. We have ran 100 models, varying the boosting method (‘gbtree’, ‘gblinear’ or ‘dart’, see https://xgboost.readthedocs.io/en/latest/index.html), the maximum depth of the trees (from 1 to 8), the learning rate (0.05, 0.1 or 0.2), the number of estimators (30, 100 or 300). For each of these 100 cases, we have used a 5-fold cross-validation. The best model has been found to correspond to the dart boosting method, a learning rate of 0.05, a maximum depth of trees of 6, and 300 estimators. In this case, the error is of 6.82 %. This is in any case around 3 times worse than the best result we obtained with our baseline DNN.

In all the cases, we used the same data-points for training and testing (including some ’unrealistic’ datapoints, see above), and we compared the resulting accuracy to the one obtained with our baseline DNN1 model (3 hidden layers of 128 units each). The results are presented in Table 1. As can be seen in the table, the DNN outperforms all other methods we considered by at least a factor 2 in accuracy. Interestingly enough, the DNN with 3 and 2 hidden layers gives similar results, with however some slightly better results with 3 hidden layers.

It should be noted, however, that the computer time required to train the DNN is much larger than for the other methods. For inference (computing the target value given a set of input values), the required computer power using DNN is also larger than for other methods. It remains, however, orders of magnitude smaller than the time required to solve differential equations.

6 Conclusion

We trained Deep Neural Networks to compute the critical core mass and envelope masses of forming planets, for a variety of conditions (formation location, temperature and pressure in the disc, core mass, solid accretion rate). The resulting DNNs, which can be easily implemented with the tools we provide on github (https://github.com/yalibert/DNN_internal_structure/), give very similar results to the ones obtained by solving the internal structure equations, using a much reduced computer time. We also showed that some fitting formulas found in the literature in general over-estimate the resulting planetary envelope mass modestly (e.g. I00 in the case of large opacity) or very largely (e.g. using the B15 approach). Using the aforementioned formulas can therefore completely overestimate the resulting mass of a forming planet, a caveat that can be avoided using the Deep Neural Networks we provide here.

Acknowledgements.

This work was supported in part by the European Research Council under grant 239605. This work has been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge the financial support of the SNSF. The plots shown in this paper made use of the seaborn software package https://seaborn.pydata.org.