Identifying Exoplanets with Deep Learning. VI. Enhancing Neural Network Mitigation of Stellar Activity RV Signals with Additional Metrics

The RV data were obtained from the HARPS-N (High Accuracy RV Planet Searcher- North) optical spectrograph, which observes the Sun continuously with 5 minute integration times to mitigate short-term stellar activity p-mode oscillations (X. Dumusque et al. 2015; D. F. Phillips et al. 2016). HARPS-N is a temperature stabilized, cross-dispersed, R ∼ 115,000, echelle spectrograph, covering an optical wavelength range of 383–690 nm (R. Cosentino et al. 2012).

The solar data from the 5 minute exposures of the Sun, taken throughout the day using HARPS-N, are reduced using the HARPS-N Data Reduction Software (DRS-2.3.5; X. Dumusque et al. 2021). In brief, the DRS extracts a one-dimensional sky-background-subtracted spectrum for each echelle order and solves for the spectrograph’s wavelength solution. The data undergoes cross-correlation with a digital mask based on solar absorption lines. This results in a 49-element array, which we refer to as the CCF for each echelle order. After corrections for instrumental drift are applied, the CCFs from each echelle order are summed to produce a white light high signal-to-noise ratio (SNR) CCF, which is used as the input representation for the ML method employed in the analysis described in Section 4. The DRS then extracts the RVs by fitting the white light CCF with a Gaussian function. However, stellar variability including spots and faculae causes time-varying small shape deviations in the white light CCF from a perfect Gaussian, which result in the measured RVs including contributions from both Doppler shifts and stellar activity signals.

Unlike distant stars, the Sun is a resolved disk in the sky, so additional steps must be taken to correct for the fact that light from different regions of the solar disk pass along different lines of sight through Earth’s atmosphere. First, we account for the effect of differential atmospheric extinction following A. Collier Cameron et al. (2019). Then, we identify and remove observations where clouds obscured at least part of the solar disk following a multistep process. First, we calculate a quality factor, based on the SNR and airmass of each observation, using a mixture model as described in A. Collier Cameron et al. (2019). We only include data where this quality factor (which spans 0 and 1) is above 0.99. As this technique does not detect all low-quality observations, we additionally use information from the HARPS-N exposure meter to more sensitively identify observations contaminated by cloud coverage. We define the exposure meter quality factor E as the ratio between the maximum and mean count of the exposure meter. This ratio does not pile up at 1 due to a slight delay in opening the shutter but can be well modeled by a Gaussian function to identify outlying observations where E is large. After fitting the distribution of E, we remove all data where E is higher than 3σ from the mean value of the entire dataset. Finally, as a last cut, we also remove observations where the resulting RV is more than 5σ from the mean value. This rather strict approach of removing low-quality data may indeed remove good data too, but is purposefully strict to only have the highest quality of data. Although this process removes 35% of the spectra, the number of days in the dataset is reduced by only 20% since we take a daily average for each day. The dataset from HARPS-N after applying these quality cuts contains 1148 days of solar observations between 2015 July 29 and 2021 November 12.¹⁶

We use multiple data products from the HARPS-N solar observations in our analysis, which we describe in detail in Section 3.

In addition to the HARPS-N solar spectra, we used observations of the TSI of the Sun. The TSI is the intensity of solar radiation integrated over the solar surface, and over the entire spectrum of the Sun. The TSI can be thought of as an analogous data product to a light curve observed by Kepler or TESS and, therefore, contains information about stellar activity (e.g., S. Aigrain et al. 2012). For example, increases in the TSI can be caused by the presence of faculae, while decreases can indicate the presence of sunspots.

We downloaded TSI datasets from the Total Irradiance Monitor (TIM) instruments on the SORCE satellite (G. Kopp et al. 2005; G. Rottman 2005; G. Kopp 2020) and on the TSIS-1 sensor on the International Space Station (G. Kopp 2023). The TIM on SORCE covers the entire solar spectrum—from X-ray to far-infrared—and has a high absolute accuracy of 350 ppm, and relative accuracy of 10 ppm yr⁻¹. It takes measurements at a 50 s cadence, which are combined to produce a daily averaged TSI value. The dataset includes daily averaged TSI data from 2003 February 25 until 2020 February 25.¹⁷ On dates where the quality of data is low due to shielding of the Sun from the Earth, the TSI is recorded as 0.0. We removed these observations from the dataset.

The data from the TIM on the TSIS-1 sensor covers the same wavelength range, but with an improved absolute accuracy of 100 ppm, and relative accuracy of 10 ppm yr⁻¹ (G. Kopp 2023). The dataset includes daily averaged TSI data from 2018 January 11 until 2023 July 19.¹⁸

We cross-referenced the dates with measurements of TSI with HARPS-N Solar Telescope observations and found that all but 36 days with HARPS-N observations had a corresponding TSI measurement. We excluded those 36 days from our ML analysis to ensure that we had values for all input features included in our models.

To measure the disk-averaged, unsigned, unpolarized magnetic flux, we used data from the Helioseismic and Magnetic Imager (HMI; P. H. Scherrer et al. 2012) aboard NASA’s Solar Dynamics Observatory (SDO). The HMI observes the solar disk in the Fe I absorption line at 617.3 nm. It achieves a high image resolution of 1″ by observing the solar disk using two 4096 × 4096 pixel CCD cameras, and an image stabilization system, with a data rate of 55 Mbps. We used HMI’s line-of-sight magnetogram data, which are maps of the Sun’s photospheric magnetic field. The magnetograms have a precision of 10 G, a zero-point accuracy of 0.05 G, and are taken with a cadence of 45 s.

The line-of-sight magnetogram data were downloaded and processed using the open-source python package SolAster¹⁹ (T. Ervin et al. 2022), which we also used to calculate the unsigned magnetic flux from the magnetogram data, as described in Section 3.4. SolAster downloads the data using the Sunpy (The SunPy Commun et al. 2020) query software. For each date, three files are downloaded, a Dopplergram, a Magnetogram, and a Continuum Intensity. Each of these files correspond to the measurement at 12:00:45 UTC for that day. We downloaded daily observations for the full time period of the HARPS-N Solar observations. During the HARPS-N observational baseline, there were nine dates where the HMI data were incomplete and not all three files were available, which meant the unsigned magnetic flux could not be calculated for those dates. We therefore excluded these dates and removed them from all of the datasets used in our analysis.

ML models learn most effectively when data is pre-processed and scaled in a uniform format to ensure each feature is of the same scale and considered equally (e.g., D. Singh & B. Singh 2020). This pre-processing allows our ML models to capture and learn patterns in our data more readily. For the white light CCFs, we design the input representation such that the model becomes sensitive to shape changes and not translational shifts.²⁰ The measured RVs from the CCFs contain not only true velocity shifts from Doppler reflex motion but also contain stellar activity contributions and instrumental systematics. This is because the RVs from the DRS are computed by fitting a Gaussian to the white light CCFs. However, stellar activity in the form of spots and faculae, as well as instrumental systematics, can cause shape changes to CCFs that result in asymmetries. This means that the center measured by the Gaussian fit is not the true center (i.e., the true velocity shift), but instead the true velocity shift and some additional deviation. We want to train our model to predict this deviation between these center-of-the-CCF measurements and the true velocity shift. In this way, our model then aims to predict the contribution of CCF asymmetries to the RV measurements, and thereby models shape-driven stellar activity and instrumental noise, leaving us with shape-driven corrections for our RV time series.

The steps we take to design the input representation to predict differences between the Gaussian fit CCFs and the true velocity shifts are described in detail in Z. L. de Beurs et al. (2022, 2024), and summarized here:

1. 1.  
   Weighted average CCFs. We compute an SNR weighted average of the CCFs as described in detail in Section 2.1
2. 2.  
   Remove the RV signatures of solar system planets. Since the RVs and CCFs from the DRS contain both the Doppler reflex motion of the solar system planets and stellar variability, we transform both data types from the barycentric to the heliocentric reference frame using the JPL Horizons ephemeris (J. D. Giorgini et al. 1996). We tested multiple interpolation methods (linear, quadratic, cubic) and found that cubic was most optimal. After performing this transformation, the RVs and CCFs should contain only shifts due to stellar activity and instrumental systematics.
3. 3.  
   Center the CCFs to the mean. We fit a Gaussian to each CCF and shift all of the CCFs to a common center using the expression x_i = mean(μ) − μ_i, where x_i is the amount shifted, μ_i is the center of the Gaussian fit for the CCF spectrum i, and mean(μ) is the mean value of μ across the entire dataset. This step is taken because we know that planets cause translation shifts whereas stellar variability causes shape changes to the CCFs. By centering all of the CCFs, we remove the translational differences between the observations and help the model focus on learning shape change patterns instead. We note that centering the CCFs does not require knowing the planetary reflex motion a priori and can be done for stars with unknown planetary contributions by using either pipeline-provided RVs or fitting a Gaussian and shifting by μ.
4. 4.  
   CCF normalization. First, we normalize the CCFs by their continuum level. Then, to prepare the CCFs for the neural network, we subtract the median of the CCF time series and divide by the standard deviation. This normalization is performed across all input parameters as described below so that the scale of variations of each input parameter is approximately equal and speeds up training and optimization of the model.

In addition to the standard white light CCF that uses all 69 spectral orders to compute a weighted average CCF, we generate three chromatic CCFs that use only a subset of the 69 orders: (a) a blue CCF that is computed using the first 23 orders (387.4–459.0 nm), (b) a yellow CCF that uses the middle 23 orders (453.9–554.0 nm), and (c) a red CCF that uses the last 23 orders (547.9–690.9 nm). These chromatic CCFs are pre-processed in the same way as described above to ensure that they have the same normalized input representation. However, these chromatic CCFs each probe only a subset of spectral lines compared to the standard white light CCF and, thus, could encode activity information that is wavelength dependent. The presence of H-opacity, particularly in the blue wavelengths, allows us to probe deeper into the photosphere compared to red wavelengths. Therefore, variations in RVs across the chromatic CCFs could reveal differences in the hemisphere-averaged convective blueshift at different depths in the photosphere (K. Al Moulla et al. 2022).

In addition to the CCFs, we included the EWs of the Hα line and the sodium D lines in the solar spectra as features to our ML model. We computed the EWs of these lines in each HARPS-N spectrum using the python package specutils.analysis.equivalent_width (N. Earl et al. 2023), which uses as inputs the continuum normalized spectra, and the small wavelength regions centered around each absorption line summarized in Table 1. The code outputs EWs for each absorption line, calculated as the numerically integrated area between the observed spectrum and the normalized continuum level within the specified wavelength ranges.

We calculated the mean EW per day for each of the three absorption lines. We normalized the remaining EW values (subtracting the mean of the values and dividing by the standard deviation, similar to the CCF normalization), leaving us with three new features corresponding to the three absorption lines.

We have two TSI datasets from TSIS-1 and SORCE, which span different dates as described in Section 2.2. These two datasets, containing daily averaged TSI values, overlap from 2018 January 11 until 2020 February 25, and there is an offset between them, with TSIS-1 having higher TSI values compared to SORCE, as seen in Figure 1. We also noticed that the first 100 days of the TSIS-1 data had a significantly larger scatter compared to the SORCE data, which may be caused by initial instrument calibrations during the first period of operations for TSIS-1.

Zoom In Zoom Out Reset image size

To merge the datasets and remove the anomalous measurements at the beginning of the TSIS-1 time series, we took the steps described below. To start, we removed the first 110 days of the TSIS-1 data for our analysis, which results in a smaller overlapping region between the two datasets that spans 2018 May 1 to 2020 February 25. We then computed the mean of the difference between the two datasets within this region and added this value to all of the SORCE data, as seen in Figure 1. Then, for the overlapping region, we took the mean of the two daily averaged TSI values for those days. For dates where we only had values for one of the two datasets, those values were used. Lastly, we normalized all of the TSI data in the same way as all of the other input features (subtracting the mean and dividing by standard deviation) before feeding them into the neural network.

In addition to using the TSI values, we also computed the time derivative of the TSI, since this has also been shown to be important for predicting stellar activity signals, such as when using the FF’ method (S. Aigrain et al. 2012). In particular, sunspots, which move across the surface of the Sun, could potentially be identified through changes in TSI.

To find the change in the TSI, we first fit a spline to the un-normalized TSI data using keplerspline,²¹ as shown in Figure 2, where we used on average a 3.4 day break point spacing. Then, we used the central difference scheme to numerically differentiate the spline. The derivative values were normalized, again by subtracting the median value and dividing by the standard deviation.

Zoom In Zoom Out Reset image size

Another feature used in our model is the disk-averaged, unsigned, unpolarized, magnetic flux of the Sun, which we obtain using the SolAster package, as described in 2.3. The package identifies and removes regions without significant magnetic activity, where the unsigned magnetic field, B_obs, is less than a threshold chosen to be 8G. SolAster then calculates the disk-averaged, unpolarized, unsigned magnetic flux, , using

where α and β are the coordinates on each point on the solar disk, I_αβ is the uncorrected continuum intensity map, and where B_obs,αβ is the corrected observed magnetic field strength from the magnetogram (R. D. Haywood et al. 2016; T. Ervin et al. 2022).

After computing the unsigned magnetic flux using SolAster, we finally normalize this feature (again by subtracting the median and dividing by the standard deviation) such that it can be used as an NN input feature.

We also used the S-Index (the ratio of the flux in the core of the calcium II H and K lines to that in the continuum; e.g., A. H. Vaughan et al. 1978) as a feature. These values were produced by the HARPS-N DRS (see footnote 16). In the DRS, ghost contamination is automatically corrected (see X. Dumusque et al. 2021). Before inputting into the neural network, the S-Index required some pre-processing. The raw data contained a significant number of outliers, which we identified using the same quality factor we used to identify good CCFs (A. Collier Cameron et al. 2019). We only kept data where the quality flag is over 0.99. This cut removed 42% of the data points, but only 31% of dates, most of which correspond to dates when the CCF data were also low quality. This additional cut therefore removed only 90 additional dates from the dataset. We then normalized the S-index measurements in the same way as all other input features.

We only retained dates in our dataset for training and evaluation of the neural network when there were measurements for all input features we considered. After removing all dates with either missing or outlier values from each independent metric, the final dataset consisted of daily average observations on 978 individual dates between 2015 July 29 and 2021 September 08.

We use feedforward neural networks, where the network takes an input x, and attempts to approximate a function f^*, so that it can output the value y, where y = f(x;p), and p are the parameters it learns to produce the best approximation of f^*. In our case, f^* is the true stellar activity correction, and f is our neural network approximation. The network is composed of many simpler functions. For example, the network might be the chain

where f⁽¹⁾, f⁽²⁾, and f⁽³⁾ are three functions, and f⁽¹⁾ is the first layer, f⁽²⁾ is the second layer, and f⁽³⁾ is the final output layer. The output from one layer is the input to the next layer. The length of the chain, in this case three, is called the depth of the model. We want f(x) to match f^*(x) as closely as possible.

To train the neural network, we used labeled training data that consists of input and output pairs (x, y). The neural network determines how to model the other functions in between the input and output, which are called hidden layers. The width of the model is given by the dimensions of the hidden layers. We can think of each layer as being composed of many units, where each layer takes many inputs from the previous layer, and produces a single output.

CNNs are commonly used for pattern recognition in structured data where the proximity of the dimensions contains relevant information. CNNs learn to identify local patterns across the entire input space. There are three types of layers in a CNN: convolutional layers, pooling layers, and fully connected layers. With each layer i, the complexity of the CNN increases, a larger portion of the input can be identified, and the final layer produces an output based on the entire input.

4.2.1. Convolutional Layers

We can first consider the convolutional layers, where we apply a cross-correlation operation. For each layer, i, we apply a one-dimensional discrete CCF from a stack of T vectors for t = 1 to T of length n_i−1 called the input, to the stack of L output vectors for l = 1 to L, or feature map,

where the convolution kernel, , is a vector of length m_i of learned parameters during training, is a vector of length n_i of learned bias parameters, and ϕ is an activation function as described in Section 4.2.4. The * represents the convolution operator where for a vector x, feature map s, and kernel w,

The kernel size is less than the input vector size so that it can be applied several times across regions of a. It is usually small, around m_i = 3 or 5, so that it can effectively detect local feature changes of a. Typically, m_i is odd, as then the units from the previous layer can be centered around an output pixel, which would otherwise not be possible and cause distortion for an even valued m_i.

4.2.2. Pooling Layers

We can next consider pooling layers. Pooling layers compress the information output from the previous layer into a lower-dimensional space. We use a specific version of pooling called max pooling, where given the output of the previous layer as the input, the output of the pooling layer is a vector of the maximum values of small regions along the input vector. The stride length is how many units apart each region is, and the pooling width is the number of neurons summarized in each region, as can be seen in Figure 3. Pooling helps to keep the output invariant to small translations of the input, which is useful when we solely want to know if a characteristic is present, rather than its exact location.

Zoom In Zoom Out Reset image size

4.2.3. Fully Connected Layers

Finally, we consider the fully connected layers, where the last layer outputs the final prediction. Every unit in the fully connected layer takes the entire output from the previous layer as the input. The activation is defined by

where i is, again, the layer number in the CNN, c_i is a vector of length n_i of activations in layer i, W_i is an n_i × n_i−1 matrix of learned weights, d_i is a vector of length n_i of learned bias parameters, and ϕ is an activation function as described in Section 4.2.4.

4.2.4. Activation Functions

An activation function is a nonlinear mathematical operation that determines the output from a unit, introducing nonlinearity to enable the network to learn complex patterns and make better predictions. The choice of ϕ depends on the application. The most popular activation function, rectified linear unit (ReLu), is given by

For many activation functions, the derivative vanishes as x → ±∞, which can lead to slow convergence for gradient-based optimization algorithms. The derivative of the ReLu function on the other hand does not face this issue, and so tends to converge quicker and with better performance. We used ReLu activation functions in our CNN layers.

Neural networks are trained to minimize a loss function. We provide the neural network with a training set containing examples of inputs and true values, and the loss function is a measure of the difference between the prediction from the model, and the true value. The mean squared error (MSE) is a common loss function used for regression tasks, and is given by

where y₁, y₂, …, y_M are the true labels of all M examples in the training set, and are the predicted outputs given parameters p, where p is the vector of the parameters in the model. The values of p are learned during training. For the CNN, in the convolutional layers, these are the elements of all convolutional kernels and bias vectors b_n from Equation (3), and for the fully connected layers, these are the weight matrices W_n and bias vectors b_n from Equation (5), where n corresponds to the index of the final output layer.

We minimize the loss function through an optimization algorithm called gradient descent. To reach this minimum, we first start with a random set of parameters p, and then iteratively update them, so that we descend along the gradient with respect to the parameters. If we write the loss function as f(x), then this can be expressed as

where ∇_xf(x) is the gradient, and is a positive scale factor called the learning rate, which determines the step size between each iteration. The learning rate is a tunable hyperparameter as further described in Section 4.4. We tune this rate until we reach a suitable minimum value of the loss function.

An important class of minimization algorithms is called Stochastic Gradient Descent (SGD). To find the true gradient, the entire training set of size M must be used. However, this is computationally expensive, and so in SGD, a random subset of the data is used of size B_SGD, where B_SGD is called the ’batch size’, and 1 ≤ B_SGD << M. Then the minimization is performed using this approximate gradient. We used a constant batch size of B_SGD = 1024, as B_SGD does not affect the performance of the model if the other hyperparameters are well tuned, as demonstrated by C. J. Shallue et al. (2019). We used a variant of the SGD algorithm, SGD with momentum (B. Polyak 1964), and fixed the momentum parameter at 0.9.

In ML, it is the gold standard to split your data into three subsets: a training, validation, and testing set. This is done to ensure generalization, which is the ability of a model to perform well on new, previously unseen data (I. Goodfellow et al. 2016). An ML model is first trained using the training set, which commonly contains 80% of the data. The validation set contains data that the model has not seen yet and is used to evaluate the performance of the model on new inputs at every step of the learning process. For example, this can be used to optimize the model architecture. This process allows us to prevent under- and overfitting on the observations, which should be avoided since it can result in poor out-of-sample performance. Finally, after the model has been tuned using the training and validation sets, the final model’s performance can be evaluated using the testing set.

Since we are using a relatively small dataset (hundreds of examples compared to the thousands commonly used in ML), we use a k − fold cross-validation method on our training dataset, as this allows us to make out-of-sample performance for the entire dataset, rather than only on the relatively small validation and test sets. In our model, we split the data as 80% in training, 10% in validation, and 10% in testing. The steps for the k-fold cross-validation method performed on the training set are as follows.

1. 1.  
   The training data is randomly assigned into k subsets.
2. 2.  
   For each training data subset:

   • a.  
     take the current subset as the “hold-out” dataset,
   • b.  
     take the other k − 1 groups as the new training dataset,
   • c.  
     fit a model using the new training dataset, and evaluate it using the hold-out dataset,
   • d.  
     record the model performance.

Using this strategy, each of the k subsets is used in the hold-out group once. The value of k must be chosen so that there is not a high variance (lower variance implies model performance changes a lot based on training data, which can occur for higher k values) or a high bias (overestimation of the model performance occurs for high k). We used a 10-fold cross-validation method, as k = 10 has been shown empirically to result in test error estimates that do not give high bias or variance (G. James et al. 2023). We then optimized our model by first training using the 10-fold cross-validation, and using the 10% validation set to tune the hyperparameters as detailed in Section 4.6. Finally, we evaluated the optimized model performance on the 10% test set.

Compared to other optimization algorithms, neural network methods are especially capable of picking up on subtle patterns in datasets. However, they can also be especially prone to overfitting. To prevent this and ensure proper generalization, which is the ability of a model to perform well on new, previously unseen inputs, we implemented regularization methods. The generalization of a model is often described by the training error and test error. The error of predictions on the training set is called the training error, and we try to make this low for the best results to reduce underfitting. We also compute the generalization error or the test error, which is the expected value for the error on a new input, and minimization reduces overfitting. Methods involving changing the learning algorithm to reduce the generalization error but not the training error are called regularization methods. The method we used is weight decay regularization, which reduces the complexity of the model by limiting the values the parameters can take. If we write the parameters as a vector p, then on each iteration, it updates as

where is the change computed by the original algorithm at iteration i, and 1 − ε is a factor to reduce the parameter vector by. The value for ε was optimized during hyperparameter tuning.

We implemented the CNN model in TensorFlow, an open-source software library for ML (M. Abadi et al. 2015).

To minimize the loss function, we used SGD with momentum on the cross-validation set. To find the optimum hyperparameters for the model, we performed ∼300 random searches across the parameter space over the cross-validation set for the learning rate, kernel size, filters, convolutional, fully connected, and max-pooling layers, units, pool size, pool strides, weight decay, and the number of epochs as listed in Table 2. We did this for ∼10 different models, with each model containing the original white light CCF data, and one additional feature.

We evaluate a model’s performance using the rms error (RMSE) between the labels and predictions on the validation set. For a more detailed description of the RMSE, see Appendix A.1. To find the optimum hyperparameters, we first found the hyperparameter values that corresponded to the lowest RMSEs for each of the models, and then found the best value across all of the random searches and over all of the models, by visualizing how the RMSE varied using the TensorBoard²² software included in the TensorFlow library (M. Abadi et al. 2015). These values are listed as the “Optimized Values” in Table 2. These optimized values were then used for our best CNN model architecture.

We then trained several different models using the best set of hyperparameters. The input to the first model trained was only the white light CCF. Then, we included one additional feature in each model, which was added after the max-pooling layer, as illustrated in Figure 6. The additional features were the white light CCF bisector and FWHM, S-Index, the EW of the Hα lines and of the sodium D1 and D2 lines, the TSI and TSI derivative, the unsigned magnetic flux, and the chromatic CCFs. Then, for the best-performing features, we trained models including several of these features simultaneously. In the case of the blue, yellow, and red CCFs, adding one of these would result in 49 features since the CCFs are 49-dimensional arrays rather than one-dimensional arrays such as features like s-index and EW of Hα.

Although RMSE is an easily accessible metric in our Tensorflow implementation of the CNN and can be a helpful metric to assess model performance, RMSE is not resistant to outliers in the dataset, and our dataset does have some non-Gaussian characteristics and outliers²³ (see Figure 4). Therefore, to determine which stellar activity features work the best in the CNN in providing predictions of the RV, we calculate the metric σ_percentile instead. We do this using the 84th and 16th percentiles, which we describe in detail in Appendix A.2.

Zoom In Zoom Out Reset image size

We calculate σ_percentile for both the raw RVs (without any activity mitigation), which correspond to the RV labels used in the neural network, and the “corrected” RVs, RV_corrected, where we have removed the predicted activity signal from our neural networks. In particular, we define:

To determine which combination of features produces the best stellar variability model, we trained 100 different instances of each model architecture, where each instance is initialized with different random values for each parameter before optimization. We chose to train such a large number of models because any given architecture can perform better or worse depending on the input initialization, which will vary due to randomness in data ordering during training. By training many models and averaging the results, we can ensure that our model architecture comparisons are robust to these random variations.

For each of the 100 models trained for each architecture, we find the σ_percentile for the corrected RVs. Then, we calculate the median and standard error of the median (see Appendix A.3 for details) over the 100 models. By using these metrics for each model architecture, we can robustly and fairly compare each model’s performance.

The median and uncertainty on the σ_percentile scatter in the corrected RV time series for each model architecture are shown in Figure 5 and Table 3. The figure shows the results across the test set. It includes models using only one feature in addition to the white light CCFs, as well as models with multiple additional features. Table 3 focuses on models with multiple additional features. Additional figures and tables comparing results across the cross-validation, validation, and full datasets individually for both RMSE and σ_percentile, for all model architectures, can be found in the Appendices (Figure 15, Table 4, Figure 16, and Table 5).

Zoom In Zoom Out Reset image size

Overall, while the relative performance of the models varies slightly depending on which datasets and scatter metrics are used for evaluation, a few consistent patterns emerge. In particular, we see that the unsigned magnetic flux, TSI, TSI derivative, chromatic CCFs, S-Index, Hα EW, and the contrast and FWHM of the white light CCFs consistently reduce the σ_percentile scatter most compared to the models with no additional features added (only the white light CCFs are input). We also note that the bisector of the white light CCFs, and EW of the sodium D1 and D2 absorption lines do not seem to improve the model performance, suggesting that they may not contain any additional information that is not already captured by the white light CCFs.

Many of the architectures we considered performed similarly well, so it is not possible to truly identify which particular features are optimal. Nevertheless, for the sake of simplicity in understanding the results of our tests, we identify a nominal “best-performing” model and use its predictions in our analysis for the rest of this paper. We chose the model with the lowest σ_percentile value for the corrected RVs in the full dataset (shown in Table 5). This best model uses the unsigned magnetic flux and TSI derivative as additional features, and reduces the σ_percentile from 154.5 to 91.6 cm s⁻¹. The architecture of the final CNN using this model is shown in Figure 6. The results of this best model across the full dataset are summarized in Figure 7. For our best architecture for a model with no additional features (i.e., only the white light CCFs), we find that the corrected RV scatter is 102.6 cm s⁻¹. However, adding in the TSI derivative and unsigned magnetic flux improves the results significantly; we measure a corrected RV scatter of 91.6 cm s⁻¹. Similarly across the test set, for a model with no additional features, the corrected RV scatter reduces from 147.1 to 95.5 cm s⁻¹, and using the best model, this reduces to 93.3 cm s⁻¹. This highlights that these additional features further help the CNN to be able to model and predict the RV signals originating from stellar activity better. For the best model architecture, the raw RVs, CNN-predicted stellar activity contributions, and CNN stellar activity corrected RVs are plotted over time in Figure 8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To understand which RV signals are being modeled and removed by the best model, and whether these signals correspond to stellar activity, we examined the RVs in the Fourier domain. The raw and corrected RVs are shown in Figure 9 in the form of a Lomb–Scargle periodogram (N. R. Lomb 1976; J. D. Scargle 1982). To obtain the Lomb–Scargle periodogram, we used the LombScargle²⁴ function included in the astropy.timeseries package (J. T. VanderPlas & Ž. Ivezić 2015; Astropy Collaboration et al. 2013, 2018, 2022). We used the periodogram normalization outlined in M. Zechmeister & M. Kürster (2009).

Zoom In Zoom Out Reset image size

The peaks indicated in Figure 9 as P_rot and P_rot/2 correspond to the Sun’s mean synodic Carrington rotation period of ∼27 days and half this rotation period, respectively. The peaks at periods >200 days correspond to long-term magnetic cycles. The leftmost peaks correspond to aliases resulting from sampling of the RVs once a day.

A similar Lomb–Scargle periodogram was produced in Z. L. de Beurs et al. (2022) using only 3 yr of HARPS-N data, which we include here in Figure 9(a). For this subset of the data, the peaks in the raw RVs at P_rot and P_rot/2 appear much stronger compared to our full 6 yr of data in Figure 9(b). This is likely due to us using a dataset that is twice as long. For the 3 yr dataset, the spots that we see at the rotation period coming in and out of view at that cadence will be relatively stable. However, for the 6 yr dataset, with twice as much longitudinal data, it is more likely that there are spots that are not stable and could be out of phase, weakening the rotation signal. There is also more time at solar minimum within our larger dataset, and this is likely the primary cause of the weakened rotational signal strength.

The corrected RVs from the CNN model in the bottom panel no longer have the peaks that correspond to these stellar activity signals. These results demonstrate that the CNN can detect and remove the quasiperiodic variability arising from spots and plague networks using only the white light CCF, the unsigned magnetic flux, and the TSI derivative, and no information about the timing of the observations.

Our primary goal in developing our CNN stellar variability model is to improve our ability to detect planetary Doppler shifts in RV datasets. To demonstrate that our method preserves planetary signals while modeling stellar variability, we performed several end-to-end planet injection-recovery tests. In these tests, we injected simulated Keplerian planetary signals into the solar data at the CCF level, pass these data through our full pipeline, and then evaluate how well we can recover the orbital parameters of the simulated planets.

Using this approach, we performed four independent planet injection-recovery tests spanning a range of RV semiamplitudes. For simplicity, we assumed circular orbits with orbital periods of P = 400 days and with semiamplitudes of K = 0.10 m s⁻¹, 0.20 m s⁻¹, 0.30 m s⁻¹, and 0.40 m s⁻¹. We simulate the planetary RV signals using radvel’s kepler.rv_drive function (B. J. Fulton et al. 2018) and inject these planetary signals into our CCFs and our HARPS-N radial velocities. Although our best NN model includes unsigned magnetic flux data and TSI derivatives, we do not inject planetary signals into these datasets because we expect unsigned magnetic flux to be unaffected by the presence of planets and only 0.1%–10% of exoplanets transit (depending on orbital orientations). We do no want to limit our model to only be able to detect transiting exoplanets, so we intentionally do not inject planetary transits in our TSI measurements.

To prepare our CCFs and RVs for the neural network model and inject the planet signals into them, we follow the same steps as described in Section 3.1, with the modification that we inject the planetary signal between Step 2 and Step 3 and perform one extra step after Step 3. For the CCFs, the planet signal is injected by applying a translational Doppler shift to the CCFs, and for the RVs, the planetary signal is simply added to the HARPS-N RV values. As part of our normal pipeline, we then perform Step 3 described in Section 3.1, which is to shift all CCFs to a common center by fitting each CCF to a Gaussian and extracting the μ values from the Gaussian fit. Notably, this measures and removes translational shifts and is how we would approach any dataset with unknown planetary signal(s). After Step 3, we perform one extra step for the injection-recovery process, which is to subtract the μ value from our RV values so that the RV values used to train our neural network do not contain translational shift information and focus on stellar variability contributions. Next, we proceed with Step 4 of our pipeline and finally train our neural network on the CCFs, TSI values, unsigned magnetic flux values, and the RVs. After training, we apply the neural network model’s stellar activity correction to the RVs to get our corrected RV time series. We add our μ values back to the now stellar activity corrected RVs since those are our estimates of planetary reflex motion.

Next, we compute Lomb–Scargle periodograms on the corrected RVs with the planetary Doppler motion. We require that for a planet signal to be detected, its peak in the periodogram must exceed the false-alarm probability of 0.1%. With this criteria, we found that the signals corresponding to semiamplitudes of K = 0.30 m s⁻¹ and 0.40 m s⁻¹ on 400 days orbits produce significant peaks in the periodogram and the lower-amplitude signals are not detected. We performed a Markov Chain Monte Carlo (MCMC) orbital fit for both planets using the software edmcmc (A. Vanderburg 2021), which incorporates differential evolution to speed up the MCMC fitting process. We applied uniform priors to the period, semiamplitude, and phase that were broadly centered around the injected values. We fixed the eccentricity to zero. A comparison of the injected planetary orbit and MCMC fit are shown in the left and right panels of Figure 10 for K=0.30m s⁻¹ and K=0.40m s⁻¹, respectively. For both planets, the semiamplitudes estimated from the MCMC, m s⁻¹ and m s⁻¹ agree within error with the true values. However, as can be seen in the left panel of Figure 10, there is a mismatch between the injected and retrieved planet curves for the 0.30 m s⁻¹planet, which is due to the MCMC best-fit value of the period ( days) being overestimated. We therefore consider the K = 0.30m s⁻¹ case a weak detection. The lowest-amplitude signal we are able to confidently retrieve is the K = 0.40 m s⁻¹ at a 400 day orbital period, which corresponds to a 5.16 M_E planet around the Sun.

Zoom In Zoom Out Reset image size

To qualitatively understand how the magnetic flux and TSI derivative inform the neural network’s predictions of stellar variability, we plot in Figures 11 and 12 the time series for the values of the unsigned magnetic flux, TSI, TSI derivative, raw RVs, predicted RVs, and corrected RVs from our best model during periods of low and high magnetic activity on the Sun. We focus on a period of 100 days, centered around months with one of the highest and lowest number of sunspots in our dataset, respectively. The high activity period occurred on 2015 October, with a mean of ∼64 sunspots over the month. The low activity period occurred on 2019 December, with a mean of ∼1.5 sunspots ²⁵ (F. Clette & L. Lefèvre 2015).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the high-magnetic activity plot shown in Figure 11, we can see that the unsigned magnetic flux follows the variations in the raw RVs, with similar peaks at around the dates 5260, 5295, and 5330, which correspond to 2457620 Julian Date (JD), 2457295 JD, and 2457330 JD, respectively. This is expected due to an increase in magnetic flux during active periods. The TSI instead dips around these dates, indicating a decrease in brightness due to sunspots. The TSI derivative similarly has shape changes at the same positions as the raw RVs, with the most prominent at around the dates 5260 and 5330 in the figure, which correspond to 2457260 JD and 2457330 JD, respectively. We can see that the predicted RVs from our best model capture these RV variations well, with predicted peaks at the same dates. However, the amplitude of the peaks predicted by our neural network are often lower compared to the raw RVs, which is likely due to our dataset mostly containing periods of low solar activity. We suspect that the neural network struggles with predicting the amplitude of the largest activity signals because there are not many examples of such an active Sun in the training set. We anticipate that the network’s ability to predict these high-activity time periods will improve if we increase the number of high activity examples in the training set. The value of σ_percentile for the raw RVs in this period is 154.0 cm s⁻¹ and for corrected RVs is 107.8 cm s⁻¹, leading to a reduction of 30.0%.

For the low-magnetic activity time period shown in Figure 12, the predicted RVs follow the amplitude of the raw RVs well, where the value of σ_percentile for the raw RVs is 90.6 cm s⁻¹ and for the corrected RVs is 81.5 cm s⁻¹, leading to a reduction of scatter of 10.0%. The network is performing well, even when there are no strong activity signals in the dataset, and it knows to ignore the information from the features (the unsigned magnetic flux, TSI, etc.), leading to predicting very low activity signals. Interestingly, the network predicts a jump in the RV time series at Julian date 2458840. This is likely a correction for an instrumental systematic (see Section 6.2).

In addition, we visualize the entire dataset for each of the one-dimensional features in Figure 13. The raw RVs are shown in pink. We can see that some of the features that did not perform well, such as white light CCF bisector span and white light CCF contrast do not appear to correlate well with the variations in the raw RVs. Bisector span is known to be particularly sensitive to spots. However, over the course of the last 3 yr included in this analysis, the Sun was mostly in solar minimum and had few spots and plagues, so this may potentially explain the low correlation of the bisector span with the RVs. On the other hand, features that perform well such as unsigned magnetic flux, TSI, TSI derivative, and s-index have similar overall shape and patterns to the raw RVs. The predicted RVs from our best model are shown in blue, and capture the variations of the raw RVs well, especially at periods of high stellar variability.

Zoom In Zoom Out Reset image size

Although our main goal in developing this neural network was to correct astrophysical noise sources, we found that it was also able to correct instrumental artifacts. In particular, the network learned and corrected for the impact of cryostat warm-ups, which are events that occurred approximately twice a year between 2015 and 2021 and introduced small jumps in measured radial velocities. These warm-up dates are marked on the raw RVs in the top panel of Figure 14, and on the predicted RVs from the CNN in the second panel. During periods of higher magnetic activity (the first half of the dataset), the RV offsets introduced by the warm-ups are smaller than the activity variations, so it is not easy to notice their effect. However, during periods of lower activity (the second half of the dataset), there are apparent jumps in the measured RVs occurring at the warm-up dates. This is illustrated by the two examples in the bottom two panels of Figure 14. (A similar effect is also visible in Figure 12 at time 6840 days, or Julian Date 2458840). These two figures also show the neural network’s predicted RVs at the times of the cryostat warm-ups, which also show jumps, indicating that the network learned to predict instrumental artifacts in addition to the stellar activity signal.

Zoom In Zoom Out Reset image size

The fact that the neural network can predict both instrumental and activity signals is interesting and suggests avenues for future improvement. The neural network was designed primarily to model stellar variability, and we therefore did not include features that directly trace the instrument’s performance. As a result, we did not expect the neural network to mitigate instrumental-driven RV variations. Instead, it must have identified instrumental artifact events based on shape changes in the white light CCFs, since that is the only input into our best model that is derived from HARPS-N. In the future, we may be able to enhance the ability of the network to remove artifacts by inputting instrument diagnostics like measurements from temperature and pressure sensors to the network. This potentially makes neural networks an even more powerful tool for post-processing RV measurements, since identifying and removing instrumental artifacts will become even more critical as we try to find increasingly small planets with the next generation of extreme precision spectrographs.

In this work, we have demonstrated an improvement in RV scatter on the Sun to about 90 cm s⁻¹. If we could achieve similar improvements in RV measurements of other stars, it would considerably improve our ability to detect small planets with the RV method. However, the 90 cm s⁻¹ scatter in the corrected RVs is still far larger than the 10 cm s⁻¹ amplitude of an Earth-like planet around a Sun-like star.

We therefore wish to understand how to improve the precision further and why it is being limited. The corrected periodogram in Figure 9 shows that most of the quasiperiodic stellar signals at the stellar rotation period have been effectively removed by the CNN model predictions, although some residual still remains near 1/2 the rotation period. Moreover, in simulated data with only magnetic activity signals, Z. L. de Beurs et al. (2022) were able to use CNNs to achieve RV precision of a few centimeters per second using this method. This implies that we are able to remove most of the contribution of magnetic activity signals either at the stellar rotation period or in the long-term activity cycle, so much of the remaining RV scatter must be from some other source.

There are two possible factors dominating the remaining RV scatter in the HARPS-N data: (i) instrumental noise and (ii) supergranulation.²⁶ In terms of instrumental noise, σ_instrumental, HARPS-N requires frequent calibrations to ensure accuracy of the generated wavelength solutions, which limit the precision of the RVs. X. Dumusque et al. (2021) found that the wavelength solutions from the HARPS-N DRS we use for our data change by around 49 cm s⁻¹ on day-to-day timescales. While we showed that our neural network can predict and correct at least some instrumental artifacts, we do not expect that most of our input features to the CNN are able to trace the day-to-day variations in the instrument wavelength solution that dominate the scatter found by X. Dumusque et al. (2021).²⁷ The second factor likely dominating our remaining RV scatter is supergranulation, σ_supergran, which occurs on inactive regions of the solar surface and would not be traced well by many of the activity indicators in our model.²⁸ Additionally, supergranulation may manifest in the CCFs as translational shifts (N. Meunier & A. M. Lagrange 2020), to which our model is not sensitive due to some of our data pre-processing steps described in Section 3.1.

Supergranulation results in RV variations of a similar scale to instrumental noise, with estimates ranging from 68 cm s⁻¹ by K. Al Moulla et al. (2023) to 86 cm s⁻¹ by B. S. Lakeland et al. (2024). Likewise, we do not expect that most of the supergranulation signal can be predicted given our inputs; it is not magnetic, so most of our additional features do not apply, and it should produce only small changes in the shape of the white light CCF. The overall contribution to RV signals from these two factors is given by

which, for σ_instrumental = 49cm s⁻¹ and σ_supergran = 68–86cm s⁻¹, gives σ_tot = 84–99cm s⁻¹. This is similar to the 92 cm s⁻¹ residuals of our CNN model.

One avenue for future work is to apply the lessons learned from this investigation of solar data to observations of other stars. Already, the neural network white light CCF-based method of Z. L. de Beurs et al. (2022) has been modified to search for planetary signals around other stars (L. L. Zhao et al. 2022b; Z. L. de Beurs et al. 2024). Other stars have smaller datasets, so modifying our method requires a linear regression approach, where a white light CCF-based stellar activity model and Keplerian signals are simultaneously fitted. This approach has the advantage of not requiring the prior removal of stellar activity signals before performing a planet search and allowing us to immediately apply these methods to other stars without requiring the extensive training set necessary for a neural network. Thus far, these methods have been successfully applied to five systems and reduced their RV scatter (L. L. Zhao et al. 2022b; Z. L. de Beurs et al. 2024). In the future, these linear regression methods could incorporate some of the features that we found to contain complementary information to the white light CCFs in this work. For example, we could add the chromatic CCFs, a proxy for the unsigned magnetic flux for stars beyond the Sun (e.g., F. Lienhard et al. 2023), and Hα into this linear regression and test how much this helps our stellar activity models for stars spanning the HR diagram. However, since other stars have smaller datasets, we will need to investigate dimensionality reduction techniques if we wish to add all of these additional features.

In the future, we are also interested in looking for new input features that can be used to predict granulation and supergranulation for the Sun. Although testing the applicability of these methods developed for the Sun to other stars is essential, developing new methods by focusing on solar observations can also give insight into new ways of mitigating variability, which can later be applied on other targets. In this way, the Sun is a great test-bed for trying out new features to see if they can help improve RV precision, and identifying new correlations to stellar variability. Since we find that the current noise floor for solar observations appears to be dominated by supergranulation (B. S. Lakeland et al. 2024), we plan to look for new correlations with these phenomena in the future by using extremely high-resolution spectra at higher resolutions than HARPS-N’s (see, e.g., M. L. Palumbo et al. 2024), or other observables, such as SDO’s Dopplergrams or UV continuum observations. Although many of the types of observations SDO provides for the Sun may not be readily available for other stars, they can still help probe the underlying physics, find new correlations, and potentially inform the design of future extreme precision spectrographs to incorporate our knowledge of stellar variability.

The rms error (RMSE) is a standard metric used to quantify the performance of a model by assessing the differences between the values predicted by the model and the values actually observed. The RMSE is defined as

where M represents the total number of samples being input or predicted. RV_predicted,i denotes the predicted value for the ith sample from the model, which in our analysis is an NN, which utilizes input activity indicators for its predictions. RV_labels,i denotes the actual, observed values for the ith sample, which, for our NN models, are the raw RVs without any activity corrections applied.

A lower RMSE value indicates a model that closely predicts the observed outcomes, reflecting its higher accuracy and performance. Therefore, in our analysis, we conclude that the NN models with the lowest RMSE have the optimal model architecture or most effective input activity features.

When selecting the optimal model architecture and identifying the most effective stellar activity features for our NN in predicting RVs, we initially rely on the RMSE metric as described above. However, RMSE is not resistant to the influence of outliers. We therefore also wish to consider outlier resistant metrics to assess how effective our models are.

The σ percentile metric, σ_percentile, is a robust statistical measure used to assess the variability of a dataset while minimizing the influence of outliers. We define this as

where RV_84th% and RV_16th% are the 84th and 16th percentiles of the RVs, respectively. This metric measures the spread of the middle 68% of the data, providing a summary of data variability similar to the standard deviation in a normal distribution. However, unlike the standard deviation, σ_percentile is less sensitive to extreme values or outliers in the distribution, making it particularly useful for datasets with non-Gaussian characteristics.

When we determine the most effective stellar activity features, we find the σ_percentile of the corrected RVs 100 times for each model architecture. This is to ensure our results minimize the effect of random variations due to data ordering during training. We then calculate the median and standard error of the median of the 100 σ_percentile values for each architecture.

The standard error of the median is a statistical measure used to quantify the uncertainty of the median estimate of a dataset. It is especially useful in situations where data distributions are not symmetric or contain outliers. It can be approximated for a dataset of N samples by the formula:

where ς is the standard deviation of the dataset (J. Maindonald & W. J. Braun 2010).