Locating Hidden Exoplanets in ALMA Data Using Machine Learning

We run 1000 3D smoothed particle hydrodynamics (SPH) simulations using the PHANTOM (Price et al. 2018) code. System properties, such as stellar mass, disk properties, observational resolution, etc., were LHC-sampled from Table 1. A significant percentage (25%) of these simulations were withheld for testing. Each simulation uses 10⁶ SPH particles and evolves for 100 orbits. A given planet's accretion radius is set to one-fourth of its Hill radius. The surface density profile is set to Σ ∝ R^−p , and the sound speed is set to c_s ∝ R^−q .

The simulation results are used to create velocity channel maps in only ¹³CO to keep the size of the parameter space manageable, although ¹²CO and C¹⁸O are also reasonable choices. The channel maps are created using the MCFOST radiative transfer code (Pinte et al. 2006, 2009). For each simulation, multiple points between 10 and 100 orbits are chosen randomly. Since one-third of the simulations have no planets, six snapshots are chosen from each of those simulations while only two snapshots are chosen for simulations with planets. This eliminates the class imbalance. From each snapshot, a velocity channel cube is made using either the ¹³CO J 2 → 1 (220.4 GHz) or 3 → 2 (330.6 GHz) transition, which is randomly sampled. These transitions are used because they were also used to identify the planets in HD 163296 and HD 97048 (Pinte et al. 2018, 2019). The final result is a velocity cube, where each velocity channel is an image with 600 × 600 pixels at a resolution of 1 au pixel⁻¹. This velocity cube is used as the model input. Channel maps are calculated using 10⁸ photon packets. The dust-to-gas ratio is 1:100, and dust is composed of carbon and silicates (Draine & Lee 1984). Stellar radii and temperatures are calculated with isochrones given by Siess et al. (2000) using the LHC-sampled masses and ages. A spherical beam is used.

The velocity channel maps are convolved spatially and spectrally, and noise is added by LHC sampling from Table 1. rms noise is added independently for each pixel, which more faithfully replicates observational effects such as hot pixels. This method results in an average rms of ∼10% for each image.

Except for the specific values, this process largely follows that done in Terry et al. (2022). Figures 1 and 2 show the results of convolving selected velocity channels for a disk with and without a planet, respectively. Figure 1 shows that, while the perturbation is strongest in the velocity channel that covers the planet (rightmost column), the perturbation is not necessarily completely localized and can be visible in distant regions of the disk. This is an expected behavior. The deviation is strongest at the location of the planet, but the spiral wake leaves further imprints as it disperses throughout the disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By itself, a single velocity channel is poorly suited for our method. Only a few channels display kinks, and channels with kinks can neighbor the exact channel containing the planet. We therefore stack all C velocity channels into a single image of dimension (H × W × C), where H and W are the height and width of the image (600 × 600 pixels for these purposes).

The number of velocity channels in an observation varies but mostly falls in the range ∼40–80 channels that cover the disk (Andrews et al. 2018; Öberg et al. 2021). We therefore use C = 47, 61, and 75 channels for each observation. The number of channels must be odd to cover the disk symmetrically and still retain the systemic (Δv = 0) channel. The input velocity of the simulated channels ranges from ∣Δv∣ = ∣v_systemic − v∣ ≤ 2.5 km s⁻¹ to fully cover the disk.

The input for all models is a (600 × 600 × C) image, and their output is a two-component vector with softmax activation (Goodfellow et al. 2016). The softmax function maps the output vector such that the sum of the components equals one. Each component can be interpreted as the predicted probability that the input belongs to the corresponding class (i.e., "contains at least one planet" class or "does not contain a planet" class). An Adam optimizer (Kingma & Ba 2014) is used to add momentum, which remembers previous gradients, and gradient scaling to a simple gradient descent algorithm. This helps escape local minima and improves the speed and range of convergence. The cross-entropy loss (Hinton et al. 1995) penalizes confident predictions that are incorrect. Training stops when there is no improvement in the validation loss after eight complete epochs, i.e., eight sweeps through the training data.

Multiple architectures are used to compare performance on two tasks: predicting if the system hosts at least one planet and determining its azimuthal location.

We use two different models based on PyTorch (Paszke et al. 2019) implementations: EfficientNetV2 (Tan & Le 2021) and RegNet (Xu et al. 2022). Neither model uses the default hyperparameters or pretrained weights.

EfficientNetV2 is an updated version of EfficientNet (Tan & Le 2019) that maximizes performance and minimizes the number of trainable parameters and training time. RegNet, on the other hand, seeks to maximize performance by adding variants of convolutional recurrent neural networks to ResNet (He et al. 2016). This can extract greater spatial information at the cost of significantly more trainable parameters.

We perform Bayesian hyperparameter tuning using WANDB (Biewald 2020) to find separate sets of hyperparameters that minimize the validation loss for C = 47, 61, and 75. The default parameters for these sweeps are based on default versions "EfficientNetV2 S" and "RegNetY 16GF."

Each architecture is evaluated for accuracy and area under the receiver operating characteristic curve (AUC; Hanley & McNeil 1982). Accuracy in this context is defined as whether a model correctly predicted that a system with planet(s) has at least one planet and similarly for systems with no planets. We choose AUC because it gives a robust measure of the performance of the model across different decision thresholds, i.e., the softmax threshold required to claim the presence of a planet. Substantial confidence is required to claim a planet detection, but setting the decision threshold too high can lead to false negatives. False positives may be more acceptable because further analysis can overrule them, whereas a false negative may lead to no further analysis and result in the planet not being detected. The optimum threshold for this boundary is open for discussion by the community. We show multiple thresholds (50%, 75%, 90%, and 95%) in Figures 3 and 4 to demonstrate how the threshold for claiming the presence of at least one planet affects the results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While some protoplanetary disks exhibit gaps in the continuum that are highly suggestive of a planet, this does not help one find the azimuthal location. Planet-induced non-Keplerian kinks tend to be strongest near the planet, so our model must determine if and where the kink is in our data in order to give information on the azimuthal location.

One simple approach to this is to look at the activations inside the model itself. The activation strengths highlight regions that the model finds particularly informative for classification. EfficientNetV2 uses SiLU activations (Elfwing et al. 2018), and RegNet uses ReLU (Glorot et al. 2011). We use mean-subtracted activations to achieve this. For a given pixel with an activation value x, the mean-subtracted activation $x^{\prime}$ is defined as $x^{\prime} =| x-\langle x\rangle |$, where 〈x〉 is the mean activation over the entire image. Activations occur most strongly at the location of the strongest deviation from Keplerianity, i.e., the kink.

We train six models in total: an EfficientNetV2 and a RegNet, each for 47, 61, and 75 velocity channels. Using a typical decision threshold of 0.5, each model, other than the 47-channel RegNet, has an accuracy of at least 93%. When this decision threshold is increased to 0.95, four models have an accuracy of at least 92%. Five models have an AUC of at least 0.98. This confirms that our models have learned to correctly detect the presence of planets using kinematic signatures.

Table 2 shows the number of parameters as well as relevant performance measures for each final model. All EfficientNetV2 models converge to a similar number of trainable parameters, whereas the number of parameters for the RegNets increases sharply with the number of channels. RegNet is a larger and more complicated architecture, so this is expected. The performance of the RegNet models varies widely compared to the EfficientNetV2 models. The accuracy of EfficientNetV2 models varies by at most 8%, whereas RegNet varies by as much as 31%. By all performance metrics, while the average EfficientNetV2 outperforms the average RegNet, the 61-channel RegNet outperforms all other models.

Figure 3 shows the confusion matrix, which describes the distribution of correct and incorrect predictions, for all models applied to the withheld test set. Panel (a) shows the results using a 50% decision threshold, i.e., a planet is predicted if the confidence is over 50%. Panels (b), (c), and (d) show the same for cutoffs of 75%, 90%, and 95%. Figure 4 gives several relevant metrics for model performance and all ROC curves. Claiming a planet is present should come with high confidence, so a 50% decision threshold is not necessarily what should be used in practice. It is encouraging that, for all other models than the 47-channel RegNet, there are no qualitative changes in the results when a higher decision threshold is used; the main effect is simply increasing the number of false negatives and decreasing the number of false positives.

If there are multiple planets in a disk, then multiple kinks may be present. It is therefore prudent to check that the quality of the predictions does not depend on the number of planets in the disk. We check the possibility of this dependence by calculating the accuracy for each model for disks with a given number of planets. Figure 5 gives the results of these calculations. The p-values from Pearson correlation tests shown in Figure 6 demonstrate that only the 47-channel RegNet using a decision threshold of 50% has a statistically significant relationship between accuracy and the number of planets (p < 0.05). While the 47-channel RegNet using a decision threshold of 95% does not have a statistically significant relationship, it is clear that the performance is worse when there are planets, i.e., there are many false negatives. This is in line with this model's relatively low AUC and Figure 3. The performance of all other models does not differ significantly with the number of planets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Many disks do not show obvious signs of planets, either in continuum or line emission. Inspecting the activations of the disks helps overcome this. Figure 7 shows the result from a disk that is difficult to classify by eye. The models predict the presence of an exoplanet with high confidence, activating on a subset of channels that are deemed important by the model. This gives information on the azimuthal location, and the radial location may be inferred by either the location with the strongest kink or coincident gaps in continuum images.

Zoom In Zoom Out Reset image size

Panels (c) and (d) in Figure 7 show the raw and convolved velocity channels, respectively. If one were simply inspecting the velocity channels by eye, it would be trivially easy to miss the kink even in the unconvolved data. The convolved, noisy data show essentially no sign of the planets. Despite lacking any clear signatures that can be observed visually, the models correctly classify the disks and locate the planets. This demonstrates that our models are able to accurately find planets that the by-eye method may fail to find.

Figure 8 shows a typical activation for a disk with no planet. The entire disk is activated equally, indicating that there is no localized region that exhibits behaviors indicative of a planet.

Zoom In Zoom Out Reset image size

Next, we demonstrate the effectiveness of our model by applying it to real telescope observations. Kinematic detections of planets have been demonstrated for many disks (Pinte et al. 2020), but here we focus on the HD 163296 and HD 97048 disks as a proof of concept (Pinte et al. 2018, 2019). Using the same data as Pinte et al. (2018, 2019), we demonstrate that our models replicate the prediction and estimated location of the forming exoplanet. Figure 9(a) shows the HD 97048 in continuum, and Figure 9(b) shows the velocity channel with the planet's signature circled.

Zoom In Zoom Out Reset image size

Table 3 shows the softmax values for the observations. For HD 97048, all RegNet models correctly predict a planet with >95% confidence. Two EfficientNetV2 models correctly predict a planet with ≥85% confidence. Figure 9 shows that the velocity channels that contain the planet are strongly activated. The activation in Figure 9(d) is strong on the planet itself rather than simply the entire velocity channel.

The results for HD 163296 are not as confident, but the two best-performing models predict the presence of at least one planet with over 99% confidence.

Both classes of the models perform better when there are more input channels, a trend not seen for HD 97048. A possible explanation is the fact that there are about twice as many velocity channels that cover HD 163296 when compared to HD 97048. While HD 97048 has roughly 75 channels covering the disk, HD 163296 has roughly 150. All models applied to HD 163296 must therefore combine significantly more velocity channels, which may decrease their ability to extract information from the images. The activation structure tends to be less informative, perhaps for the similar reasons, so it is omitted. Despite this limitation, the high confidence of some of the models should encourage a human observer to seriously scrutinize the kinematic structure.

Our work comes with several limitations. All simulations were run with SPH, so we did not demonstrate that our methods would also be applicable to data from grid-based codes. However, there is typically a qualitative agreement between results from SPH- and grid-based methods (De Val-Borro et al. 2006; Agertz et al. 2007). A good direction for future work would be to train on simulation data from a variety of hydrodynamics codes and check that the same results can be obtained.

Any differences in performance of the models between the codes, e.g., if the models consistently performed better on SPH data, would be difficult to interpret, and the cause could be hard to identify. It is possible that the method of domain adaptation (Ben-David et al. 2010), which has already found success in astronomy (Vilalta et al. 2019; Alexander et al. 2021; Ćiprijanović et al. 2022), could be used to encourage the models to overcome any differences between data sets. This is an avenue that is ripe for exploration in future work.

In order to keep the parameter space manageable and demonstrate a proof of concept in this work, we therefore trained on data produced by only one hydrodynamics code (PHANTOM).

The absence of true image segmentation (Minaee et al. 2022) or object detection (Zhao et al. 2019) during training is another limitation that requires addressing in future work. While we present results in the form of activation structures that suggest the location(s) of the planet(s), these activations neither always perfectly coincide with the planets nor capture all planets. User analysis is required in order to pinpoint the radial (and sometimes azimuthal) location(s) of the planets even with the information given by the activation structures. The main utility of the presented models is therefore to give a user an indication as to whether a planet is present and motivate and direct the analysis of the velocity channels. Segmentation and object detection would address these shortcomings by precisely pinpointing the location(s) of the planet(s) without any necessary further analysis.

The assertion of a vertically isothermal disk is a simplifying assumption that is commonly used in simulations of protoplanetary disks (see, e.g., Dipierro et al. 2015; Pinte et al. 2019). It has been demonstrated that the dispersion relation and morphology of the planet's signature can change significantly if the disk is vertically stratified (Ogilvie & Lubow 2002; Juhász & Rosotti 2018).

However, the presence of any disturbance is more important than its specific morphology. There does exist the possibility that vertical stratification dampens the signatures to an extent that the models' performances appreciably decrease, but our results demonstrate that even extremely subtle signals can be found by the models. Regardless, a more realistic treatment of the thermal structure may be useful in future works.