# Leveraging space-based data from the nearest Solar-type star to better
understand stellar activity signatures in radial velocity data.

Tamar Ervin Department of Physics and Astronomy, University of California, Los
Angeles, Los Angeles, CA 90024, USA  Jet Propulsion Laboratory, California
Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA  Samuel
Halverson Jet Propulsion Laboratory, California Institute of Technology, 4800
Oak Grove Drive, Pasadena, CA 91109, USA  Abigail Burrows Department of
Physics and Astronomy, Dartmouth College, Hanover, NH 03755, USA  Jet
Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive, Pasadena, CA 91109, USA  Neil Murphy  Jet Propulsion Laboratory,
California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109,
USA  Arpita Roy Space Telescope Science Institute, 3700 San Martin Dr.,
Baltimore, MD 21218, USA  Department of Physics and Astronomy, Johns Hopkins
University, 3400 N. Charles St., Baltimore, MD 21218, USA  Raphaelle D.
Haywood Astrophysics Group, University of Exeter, Exeter EX4 4QL, UK  Federica
Rescigno Astrophysics Group, University of Exeter, Exeter EX4 4QL, UK  Chad F.
Bender Steward Observatory, The University of Arizona, 933 N. Cherry Ave.,
Tucson, AZ 85721, USA  Andrea S.J. Lin Department of Astronomy & Astrophysics,
525 Davey Laboratory, The Pennsylvania State University, University Park, PA,
16802, USA  Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory,
The Pennsylvania State University, University Park, PA, 16802, USA  Jennifer
Burt Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak
Grove Drive, Pasadena, CA 91109, USA  Suvrath Mahadevan Department of
Astronomy & Astrophysics, 525 Davey Laboratory, The Pennsylvania State
University, University Park, PA, 16802, USA  Center for Exoplanets and
Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University,
University Park, PA, 16802, USA

###### Abstract

Stellar variability is a key obstacle in reaching the sensitivity required to
recover Earth-like exoplanetary signals using the radial velocity (RV)
detection method. To explore activity signatures in Sun-like stars, we present
SolAster, a publicly-distributed analysis pipeline
111https://tamarervin.github.io/SolAster/ that allows for comparison of space-
based measurements with ground-based disk-integrated RVs. Using high spatial
resolution Dopplergrams, magnetograms, and continuum filtergrams from the
Helioseismic and Magnetic Imager (HMI) aboard the _Solar Dynamics Observatory_
(_SDO_), we estimate ‘Sun-as-a-star’ disk-integrated RVs due to rotationally
modulated flux imbalances and convective blueshift suppression, as well as
other observables such as unsigned magnetic flux. Comparing these measurements
with ground-based RVs from the NEID instrument, which observes the Sun daily
using an automated solar telescope, we find a strong relationship between
magnetic activity indicators and RV variation, supporting efforts to examine
unsigned magnetic flux as a proxy for stellar activity in slowly rotating
stars. Detrending against measured unsigned magnetic flux allows us to improve
the NEID RV measurements by $\sim$20% ($\sim$50 $\mathrm{cm\ s^{-1}}$ in a
quadrature sum), yielding an RMS scatter of $\sim$60 $\mathrm{cm\ s^{-1}}$
over five months. We also explore correlations between individual and averaged
spectral line shapes in the NEID spectra and _SDO_ -derived magnetic activity
indicators, motivating future studies of these observables. Finally, applying
SolAster to archival planetary transits of Venus and Mercury, we demonstrate
the ability to recover small amplitude ($<$50 $\mathrm{cm\ s^{-1}}$) RV
variations in the _SDO_ data by directly measuring the Rossiter–McLaughlin
(RM) signals.

††software:  astropy (Astropy Collaboration et al., 2018), barycorrpy (Kanodia
& Wright, 2018), matplotlib (Hunter, 2007), numpy (van der Walt et al., 2011),
pandas (McKinney, 2010), scipy (Virtanen et al., 2020), corner (Foreman-
Mackey, 2016), SunPy (SunPy Community et al., 2020)

##  1 Introduction

The field of exoplanet science has drastically grown in popularity and fervor
since the first confirmed exoplanet discovery around a Sun-like star (Mayor &
Queloz, 1995). This has led to significant advancements in both
instrumentation and data analysis techniques as we push towards the detection
of a Earth-like planets (Fischer et al., 2016). The radial velocity (RV)
technique has been a cornerstone of exoplanet science since the first
exoplanet discovery using Doppler velocimetry, and is credited with over 1000
additional planet discoveries (Hatzes, 2016; Fischer et al., 2016).

The RV technique searches for periodic Doppler shifts in the host star’s
spectra (Hatzes, 2016). These periodic variations are driven by the presence
of a planetary companion whose motion shifts the spectroscopic signature of
its host star. As the field strives towards the detection of smaller,
terrestrial-mass planets, improvements in RV measurement precision are
required (Fischer et al., 2016) to push beyond the current $\sim$1 $\mathrm{m\
s^{-1}}$ measurement floor. Recent advancements in RV instrumentation,
culminating in the delivery of a new generation of Doppler measurement
facilities such as NEID (Schwab et al., 2016), ESPRESSO (Pepe et al., 2021),
EXPRES (Jurgenson et al., 2016), HARPS3 (Thompson et al., 2016), and KPF
(Gibson et al., 2018) aim to push down to the $\sim$30 $\mathrm{cm\ s^{-1}}$
range. Detecting an Earth-like planet orbiting a Sun-like star requires
additional improvement down to the 10 $\mathrm{cm\ s^{-1}}$ level (Wright,
2018; Fischer et al., 2016; Hatzes, 2016).

The next challenge to improving detection sensitivity lies largely in the
removal of stellar variability, which leads to noise that can often dominate
measured RV variability (Saar & Donahue, 1997) at the $\mathrm{m\ s^{-1}}$
level (Crass et al., 2021). The signal from stellar activity can often mask or
even masquerade as planetary signals (Robertson et al., 2015; Wright, 2018).
Stellar activity signals are due to a combination of (super) granulation
(Meunier et al., 2015; Dumusque et al., 2011), oscillations (Palle et al.,
1995), meridional circulation (Meunier & Lagrange, 2020), magnetic activity,
and photospheric motion (Meunier et al., 2010a; Haywood et al., 2016; Crass et
al., 2021). These phenomena are often periodic, aligning with the stellar
rotation period and subsequent harmonics, and can consequently be mistaken for
planetary signals (Boisse et al., 2011). The lack of temporal stability across
the stellar surface, coupled with inhomogeneous stellar intensity and
differential rotation of the star, makes it difficult to robustly disentangle
stellar signals from planetary ones when studying disk-integrated spectra.

To improve our understanding of stellar activity and its effects in Sun-like
stars, we turn to our closest Solar-type star – the Sun. The sheer amount of
available solar data products, combined with established abilities to leverage
high-cadence images of the solar surface to produce maps of solar velocity,
intensity, and magnetic field strength, make the Sun the perfect candidate for
studying activity-induced temporal variability (Pesnell et al., 2012; Scherrer
et al., 2012). Using the Sun as a test case allows us to more cleanly separate
the various component velocities and analyze relationships between measured
disk-integrated RV variations and calculated solar observables.

NASA’s Helioseismic and Magnetic Imager (HMI), an instrument aboard the _Solar
Dynamics Observatory_ (_SDO_), was built as the successor to the Michelson
Doppler Imager (MDI) to study the solar surface magnetic field. Launched in
2010, it continuously observes the Sun in the spectral region of the Fe I 6173
Å line, providing four high resolution data products: line-of-sight and vector
magnetograms, continuum filtergrams, and Dopplergrams (Pesnell et al., 2012;
Scherrer et al., 2012). These measurements of the magnetic field variability,
intensity continuum, and velocity profile across the solar disk allow us to
study the Sun’s temporal variability and the effect these variations have on
solar RV’s (Pesnell et al., 2012; Scherrer et al., 2012; Haywood et al.,
2016).

In this study, we develop a _SDO_ /HMI data analysis pipeline to compliment
future extreme-precision RV (EPRV) studies of the Sun. Originally developed by
Fligge et al. (2000), this technique has been adapted by Meunier et al.
(2010b), Haywood et al. (2016), and Milbourne et al. (2019) to extract disk-
averaged quantities from spatially resolved solar observations. Our publicly
available Python pipeline, SolAster, uses data products from _SDO_ /HMI to
better characterize a suite of solar magnetic activity parameters, and
performs a simple decorrelation analysis on disk-integrated solar RV
measurements (now available from a number of RV facilities). There are two
primary activity effects that strongly impact the measured RV (on timescales
of days to months): the velocity variation due to the traversing motion of
sunspots and faculae across the rotating solar surface, and the variation due
to the suppression of the convective blueshift by active regions (Aigrain et
al., 2012). When linearly combined, these velocity components can be used to
generate an independent estimate of the disk-integrated solar RV (Haywood et
al., 2016). The individual velocity components serve as a strong proxy for
surface magnetic activity, providing a independent window into the stellar
surface that can aid in interpreting ground-based RV measurements (Haywood et
al., 2016; Milbourne et al., 2019; Haywood et al., 2020). Additionally, from
the _SDO_ /HMI data we calculate an array of magnetic observables that can be
used to gauge the effects of the size and intensities of active regions on
measured RV variations.

The paper is organized as follows. In Section 2, we describe our SDO/HMI data
processing pipeline and analysis products. We outline the data correction
process and methodology for classifying different magnetically active regions.
In Section 2.7 we discuss the calculation of the full ‘Sun-as-a-star’ RVs,
outlining how each of the velocity components are independently calculated. We
also describe the calculation of solar magnetic observables (Sections 3.3 and
3.4) and compare these results with our space-based measurements from HMI and
ground-based RV measurements from the NEID instrument (4). Finally, we apply
these calculation techniques to archival planetary transits in Section 5 to
highlight the precision of the reconstructed RVs delivered by the pipeline,
and demonstrate that magnetic variability can affect precision RV measurements
at the 10’s of $\mathrm{cm\ s^{-1}}$ level over multi-hour timescales.

##  2 SolAster \- An _SDO_ /HMI analysis pipeline

The plethora of data available from the Helioseismic and Magnetic Imager
aboard _SDO_ allows us to calculate space-based, ‘Sun-as-a-star’ radial
velocity estimates that can be directly compared to ground-based measurements.
Here we describe the underlying data products and techniques used to calculate
various solar observables using SolAster.

Before computing the RVs from the _SDO_ /HMI data, there are a number of data
preparation steps required. The suite of _SDO_ /HMI images used in this study
were: wide-band continuum filtergrams (intensity grams), line-of-sight
longitudinal magnetic field measurements (magnetograms), and maps of solar
surface velocity (Dopplergrams) (Pesnell et al., 2012; Scherrer et al., 2012)
(see Figure 1 for example images). These three data products provide the
necessary intensity, magnetic field strength, and velocity information for
active regions to be detected, tracked, and accurately integrated into the
full RV model.

HMI data products are publicly available and can be queried from the data
archive using Sunpy, a community based Python package for solar data analysis
(SunPy Community et al., 2020). In addition to providing archive querying
capabilities, Sunpy includes user-friendly methods for accessing and
visualizing solar data.

Using Sunpy, SolAster calculates a combination of velocities and magnetic
observables using the _SDO_ /HMI intensity, velocity, and magnetic field data.
Photometric and convective velocity components are independently calculated
and then linearly combined to generate ‘Sun-as-a-star’ RVs. Additionally, we
calculate both unsigned magnetic flux and filling factor, which can be used to
study the correlation between disk-integrated radial velocity and measures of
magnetic activity. These magnetic observables are calculated contemporaneously
to the RV calculations and include unsigned flux and filling factor
measurements specific to all relevant active regions (plage, intranetwork, and
sunspots). We also look at flux due to convective regions, and area cuts to
study the differing effects between large and small active regions on the
modeled RVs.

###  2.1 Coordinate transformations

Before calculating the three-dimensional Heliocentric velocity, the data must
be transformed from the Helioprojective Cartesian Frame into the Heliographic
Carrington frame, then corrected for line-of-sight projections and relative
positioning of the spacecraft. This transformation is based off the
description in Thompson (2006) and is necessary to ensure the images are
centered on the solar surface and independent of the Carrington rotation
cycle, the 25.38 day solar sidereal rotation period (Carrington, 1859). To
transform coordinate systems, we rotate the image grid from Cartesian pixel
coordinates to Heliographic Carrington coordinates by building a rotation
matrix calculated from the reference coordinates listed in each image’s FITS
header (Ulrich & Boyden, 2006). Each HMI image’s relative pixel locations in
the Heliographic Carrington frame are specified by
$(w_{ij},\,n_{ij},\,r_{ij})$ denoting the direction westward, northward, and
radially outward from disk center, respectively. This coordinate system fixes
the image onto the solar surface and allows for a determination of the
relative position of the spacecraft with respect to the Sun using only the
radial coordinate. Additionally, we constructed an array of
$\mu\;(\cos\theta)$ values for each pixel in each image, which determines the
position of the pixel relative to disk center. Flux values for pixels with
$\mu$ values below 0.3 were set to zero in all images as the limb-brightening
model is often unreliable far from disk center (Haywood et al., 2016), and
projection issues can cause non-physical fluctuations in measured values.

###  2.2 Spacecraft Velocity Correction

To isolate the solar velocity component due to magnetic activity alone, we
first corrected the Dopplergrams for the relative motion of the spacecraft.

We corrected the Dopplergrams for the motion of the spacecraft relative to the
Sun by building a pixel-wise mask of the relative spacecraft velocity. The w,
n, and r components of the relative spacecraft velocity are read in from the
FITS headers with a quoted precision of $0.01$ $\mathrm{m\ s^{-1}}$ (Hoeksema
et al., 2018). We then calculated the position of the spacecraft relative to
each pixel $ij$, which combined with the velocity components from the FITS
header, determined the required velocity correction. After the coordinate
transformation, the spacecraft is located at position
$(0,\,0,\,r_{\mathrm{sc}})$ where $r_{\mathrm{sc}}$ is the radial position of
the spacecraft relative to disk-center, and can be determined by dividing the
net distance to the Sun by the solar radius (both these values are found in
the FITS header of all _SDO_ /HMI images with keywords dsun_obs and rsun_ref).

$\begin{split}\delta w_{ij}=w_{ij}-0\\\ \delta n_{ij}=n_{ij}-0\\\ \delta
r_{ij}=r_{ij}-r_{\mathrm{sc}}\end{split}$ (1)

The $(w,\,n,\,r)$ components due to the relative motion of the spacecraft are
found in the FITS header of the Dopplergram files (obs_vw, obs_vn, obs_vr
respectively). We then project these components such that each pixel $ij$ has
a Heliocentric velocity magnitude of (Haywood et al., 2016):

$\begin{multlined}v_{\mathrm{sc},\,ij}=-\frac{\delta
w_{ij}\,v_{\mathrm{sc},\,w_{ij}}+\delta
n_{ij}\,v_{\mathrm{sc},\,n_{ij}}+\delta
r_{ij}\,v_{\mathrm{sc},\,r_{ij}}}{d_{ij}}\end{multlined}v_{\mathrm{sc},\,ij}=-\frac{\delta
w_{ij}\,v_{\mathrm{sc},\,w_{ij}}+\delta
n_{ij}\,v_{\mathrm{sc},\,n_{ij}}+\delta
r_{ij}\,v_{\mathrm{sc},\,r_{ij}}}{d_{ij}}$ (2)

where $d_{ij}$ is the distance between the spacecraft and pixel $ij$.

$d_{ij}=\sqrt{\delta w^{2}_{ij}+\delta n^{2}_{ij}+\delta r^{2}_{ij}}$ (3)

###  2.3 Solar Rotational Velocity Correction

Next, we turn our attention to the differential rotation of the solar disk
which must be accounted for when correcting the measured Doppler maps.
Differential rotation is the result of turbulent motion and convective
activity due to temperature gradients permeating outwards from the stellar
core (Schou et al., 1998). This produces a latitude-based rotation profile,
where the rate of surface rotation is maximized at the equator
($\phi=0^{\circ}$) and is inversely proportional to latitude (Schröter, 1985).
The angular velocity due to rotation in the photospheric layer ranges from
14.1-14.4 deg day-1 at the equator to 10.07 deg day-1 at the poles (Snodgrass,
1984). The sidereal rotation period for a Carrington rotation is 25.38 days,
which is the rotation rate at a latitude of 26∘, where sunspots are most often
found. This rotation period is accounted for in our coordinate transformation
to the Heliographic Carrington frame (Thompson, 2006).

Snodgrass & Ulrich (1990) used full-disk Magnetograms and Dopplergrams from
the Mount Wilson Observatory to track magnetic features on the solar surface
over time in order to build a model of the solar differential rotation
profile. They determined three constants
$\alpha_{1},\,\alpha_{2},\,\alpha_{3}$, all of which are in units of deg
day-1. The parameterization of the differential rotational profile is as
follows:

$\omega(\phi)=\alpha_{1}-\alpha_{2}\,\sin^{2}(\phi)-\alpha_{3}\,\sin^{4}(\phi).$
(4)

Using this parameterization with coefficients 14.713, 2.293, and 1.787 as
$\alpha_{1},\alpha_{2},\alpha_{3}$ respectively, we calculate our differential
rotation profile and project this onto the solar disk to build a map of solar
rotational velocity at each latitude. We then project this differential
rotation profile into the Heliographic Carrington frame to determine the
rotational velocity component at each pixel. Finally, we calculate the full
rotational velocity array based on the methods of Haywood et al. (2016) and
Milbourne et al. (2019):

$v_{\mathrm{rot}}=-\frac{\delta w_{ij}\,v_{\mathrm{rot},\,w_{ij}}+\delta
n_{ij}\,v_{\mathrm{rot},\,n_{ij}}+\delta
r_{ij}\,v_{\mathrm{rot},\,r_{ij}}}{d_{ij}}.$ (5)

###  2.4 Foreshortening Correction

The line-of-sight magnetograms measure the longitudinal surface magnetic
field. Foreshortening causes a decrease in observed spatial resolution
relative to the distance from disk center due to the geometric projection and
must be accounted for when estimating the true magnetic flux (Zhao et al.,
2016). This measured magnetic field is less than the true radial solar
magnetic field by a factor of $\mu=\cos{(\theta)}$, where $\theta$ is the
center-to-limb angle (Zhao et al., 2016). To calculate the true field
strength, we divide the observed field ($B_{\mathrm{obs}}$) by $\mu$ and
recover the full radial field:

$B_{r,\,ij}=B_{\mathrm{obs},\,ij}/\mu_{ij}.$ (6)

Additionally, we set all pixels with magnetic field strengths below the noise
threshold ($\sigma_{B_{\mathrm{obs},\,ij}}$) to zero to account for instrument
noise as described in Yeo et al. (2013). Although the noise does increase as a
function of angle from disk center ($\mu$) we take a constant minimum noise
threshold of 8G based on Yeo et al. (2013). Therefore, pixels with
longitudinal magnetic field strengths ($B_{\mathrm{obs},\,ij}$) below 8G are
set to 0 for both $B_{\mathrm{obs},\,ij}$ and $B_{r,\,ij}$. This ensures our
magnetic measurements are not contaminated by instrument noise, which would
otherwise propagate through many aspects of the analysis pipeline.

###  2.5 Limb Darkening Correction

Similar to the effect foreshortening has on the HMI magnetograms, the
continuum images are also affected by limb-darkening. We correct for this by
using a static fifth-order circularly symmetric polynomial brightness function
($L_{ij}$), with scaling coefficients determined through empirical methods by
Allen (1973). This polynomial produces a pixel-wise array of correction values
and the base intensity image is divided by these correction factors.

$I_{\mathrm{flat},\,ij}=\frac{I_{ij}}{L_{ij}}$ (7)

The flattened intensity image can now be used to classify bright and dark
regions (Figure 1), which are separated via thresholding.

###  2.6 Region Identification

The underlying assumption in our space-based RV calculation is that magnetic
activity is the primary driver of bulk RV variability in the Sun. We identify
magnetically active regions, and differentiate between regions of bright
faculae and dark sunspots, to distinguish the impact of different types of
magnetic activity on RVs.

Active regions are detected using a thresholding identification scheme
described in Yeo et al. (2013). Regions above the threshold are marked as
‘active’ and regions below the threshold are stored as ‘quiet-Sun’ pixels.
Additionally, we remove pixels near the solar limb ($\mu<0.1$) and ignore
pixels with $\mu$ values below 0.3 since the limb-darkening model is often
flawed near the limb, as was done in Haywood et al. (2016) and Milbourne et
al. (2019). We apply the same magnetic threshold described in Yeo et al.
(2013), where pixels three times the noise cutoff in unsigned radial magnetic
field strength (8G) are considered active:

$|B_{r,\,\mathrm{thresh}}|=3\sigma_{B_{\mathrm{obs},\,ij}}/\mu_{ij}$ (8)

where $\sigma_{B_{\mathrm{obs},ij}}$ is the magnetic noise level of 8G from
Yeo et al. (2013). We set any isolated active pixels, i.e. those with no
identified neighboring active pixels, to 0 (‘quiet-Sun’) as these can often be
misidentified as sunspots and may instead be instrumental artifacts.

Once active regions are identified (Fig 1), we then apply intensity
thresholding to differentiate between faculae and sunspot regions. Similar to
the magnetic thresholding previously described, we base our intensity
thresholding on values determined by Yeo et al. (2013) and used by Haywood et
al. (2016) and Milbourne et al. (2019). Pixels with flattened intensity values
above the threshold are denoted as faculae and those below the threshold are
sunspots:

$I_{\mathrm{thresh}}=0.89\;I_{\mathrm{quiet}}.$ (9)

The intensity threshold is based on $I_{\mathrm{quiet}}$, the mean flattened
pixel intensity of quiet-Sun pixels, and is calculated by summing the
flattened intensity of quiet-Sun pixels with a binary weighting array based on
magnetic thresholding:

$I_{\mathrm{quiet}}=\frac{\sum_{ij}I_{\mathrm{flat},\,ij}\,W_{ij}}{\sum_{ij}W_{ij}},$
(10)

where $W_{ij}$ is set to 1 for quiet-Sun pixels
($|B_{r,\,ij}|<|B_{r,\,\mathrm{thresh},\,ij}|$) and 0 for active pixels.

Figure 1: Example of corrected _SDO_ /HMI data products during the highest
activity day (largest magnetic filling factor) in our studied time frame (June
30, 2021). _Top left:_ Flattened continuum intensity after limb-darkening
correction (used to identify bright and dark regions) _Top right:_ Doppler
velocity map after removal of spacecraft and solar rotational velocity.
_Bottom left:_ Unsigned magnetic flux density ($|B_{r}|$). In this image,
3.25% of the Sun is covered by faculae/plage and 0.026% by sunspots, as
identified by our pipeline. _Bottom right:_ Thresholded image showing
magnetically active regions identified by our algorithms. Faculae/plage
regions are red and sunspot regions are blue.

###  2.7 Radial Velocity Calculation

Following Haywood et al. (2016), we parameterize the full, disk-integrated
solar radial velocity as a linear combination of contributions from the quiet-
Sun and active regions. Active regions produce RV variations through two
primary mechanisms: photometric effect and convective effect. Meanwhile, the
quiet-Sun RVs are primarily driven by granulation.

###  2.8 Photometric Contribution

The photometric velocity traces the rotational Doppler imbalance caused by
bright faculae and dark sunspots. The presence of bright and dark active
regions leads to an inhomogeneity across the solar disk, altering the Doppler
balance between redshifted and blueshifted hemispheres. This leads to RV
shifts of up to several percent depending on spot size and stellar activity
levels (Saar & Donahue, 1997). Sunspots are generally the dominant source of
variability in the photometric RV signal, although plage regions also
contribute to the signal (Lagrange et al., 2010). The photometric effect due
to these two factors is accounted for in$\Delta\hat{v}_{\mathrm{phot}}$, which
we calculate based on the methodology outlined in Haywood et al. (2016) and
Milbourne et al. (2019):

$\Delta\hat{v}_{\mathrm{phot}}=\frac{\sum_{ij}v_{\mathrm{rot},\,ij}\,(I_{ij}-\hat{K}\,L_{ij})\,W_{ij}}{\sum_{ij}I_{ij}},$
(11)

where $\hat{K}$ is a scaling factor based on the limb darkening correction
polynomial:

$\hat{K}=\frac{\sum_{ij}I_{ij}\,L_{ij}\,W_{ij}}{\sum_{ij}L^{2}_{ij}\,W_{ij}}.$
(12)

$I_{ij}-\hat{K}\,L_{ij}$ is the intensity map corrected for limb-darkening, as
seen in the top left panel of Figure 1. An example of the $W_{ij}$ weighting
array can be seen in the bottom right panel of Figure 1, where $W_{ij}=0$ for
quiet-Sun pixels and $W_{ij}=1$ for active pixels.

The velocity perturbations from faculae and sunspots are approximately anti-
correlated due to their opposing flux signs. When calculating the photometric
velocity component, we find that this velocity perturbation is almost entirely
driven by sunspots (see Figure 6), corroborating the results of Meunier et al.
(2010a). This is likely due to the Sun’s geometric configuration and the ratio
of bright/spot regions at the time, meaning this may not be the case for other
stars with different filling factors and distributions of bright and dark
regions.

###  2.9 Convective Contribution

Active magnetic regions have different velocity amplitudes and surface areas
distributed between upward and downward flows of solar granulation (Dravins,
1990). In the photosphere, active magnetic regions inhibit granular convective
motions of the quiet-Sun, and these convective motions manifest as wavelength
shifts of photospheric lines (Dravins et al., 1981). While the photometric
velocity variation is driven by sunspots, the convective velocity variation is
driven by larger brighter faculae regions and thus these drive the overall RV
signal.

_SDO_ /HMI images can resolve these granules, allowing us to calculate the
velocity contribution specifically due to suppression of the convective
blueshift. In convective cells, dark outward flowing plasma at the cell’s
center and downward flowing bright plasma on the cells edge leads to overall
convective blueshifts on the order of 0.5 $\mathrm{km\ s^{-1}}$ (Dravins et
al., 1981). This correlation weakens across the solar disk as we see primarily
horizontal velocity flows on the solar limb (Dravins, 1990).

The effect of the suppression of the convective blueshift varies with line
depth (Gray, 2009) and thus we expect to see different temporal convective
velocity shifts across different wavelengths. For this reason, we do not
expect perfect correlation between NEID observations and the space-based
convective velocity, as the ground-based RV measurements utilize thousands of
spectral features while _SDO_ /HMI observes velocities only in the
magnetically sensitive 6173.3 Å Fe I line. We use linear regression to scale
the _SDO_ /HMI derived convective velocity to account for this difference (see
Section 2.10).

The convective velocity is then calculated by taking the disk averaged Doppler
velocity, $\hat{v}$, and subtracting from it the disk averaged quiet-Sun
velocity, $\hat{v}_{\mathrm{quiet}}$. We subtract the quiet-Sun velocity
because _SDO_ /HMI Dopplergrams are not well calibrated nor stable over long
timescales (Haywood et al., 2020).

$\Delta\hat{v}_{\mathrm{conv}}=\hat{v}-\hat{v}_{\mathrm{quiet}}.$ (13)

We calculate the disk-averaged and quiet-Sun velocities following the
methodology of Haywood et al. (2016) and Milbourne et al. (2019):

$\hat{v}=\frac{\sum_{ij}(v_{ij}-v_{sc,\,ij}-v_{rot,\,ij})\,I_{ij}}{\sum_{ij}I_{ij}}$
(14)

where $v_{sc,\,ij}$ is the relative spacecraft velocity and $v_{rot,\,ij}$ is
the solar rotational velocity. An example of the corrected spacecraft velocity
($v_{ij}-v_{sc,\,ij}-v_{rot,\,ij}$) can be seen in the top right panel of
Figure 1. The quiet-Sun velocity is calculated by using the corrected Doppler
velocity and weighting by the intensity of quiet-Sun pixels:

$\hat{v}_{\mathrm{quiet}}=\frac{\sum_{ij}(v_{ij}-v_{sc,\,ij}-v_{rot,\,ij})\,I_{ij}\,W_{ij}}{\sum_{ij}I_{ij}\,W_{ij}}$
(15)

where $W_{ij}$ is the magnetic weighting array.

###  2.10 RV Reconstruction from Velocity Features

To estimate the full disk-integrated _SDO_ /HMI space-based radial velocities,
we follow the methodology outlined in Milbourne et al. (2019) and Haywood et
al. (2020), adapted from Haywood et al. (2016). We build a model radial
velocity variation $\Delta RV_{\mathrm{model}}$ assuming a linear combination
of $\Delta\hat{v}_{\mathrm{conv}}$ and $\Delta\hat{v}_{\mathrm{phot}}$:

$\Delta
RV_{\mathrm{model}}=A\,\Delta\hat{v}_{\mathrm{phot}}+B\,\Delta\hat{v}_{\mathrm{conv}}+RV_{0}$
(16)

where $A$ and $B$ are independent scaling factors, and $RV_{0}$ is the
relative RV offset parameter. Similar to Milbourne et al. (2019), these
coefficients are determined by linear least-squares optimization using the
ground-based RV measurements, assuming the two RV components are orthogonal
and do not have correlated noise. These scaling factors account for the
systematic differences between observations taken using _SDO_ /HMI in one line
($\mathrm{\lambda}=6173.3\,\mathrm{\AA}$), and ground-based spectra using
thousands of lines.

###  2.11 Validation of SolAster using previously published _SDO_ /HMI
measurements

As our methodology for calculating the velocity components and full model RVs
is based on the methods outlined in Haywood et al. (2016) and Milbourne et al.
(2019), we analyzed the same solar data used in these studies to verify our
performance. We use our pipeline to calculate RVs for the time frame in
Milbourne et al. (2019) to compare our model with published reference values.
When comparing our derived _SDO_ /HMI measurements with the equivalent
measurements in Milbourne et al. (2019) we find excellent agreement. The
Spearman correlation coefficients for $\Delta\hat{v}_{\mathrm{conv}}$,
$\Delta\hat{v}_{\mathrm{phot}}$, and $\Delta RV_{\mathrm{model}}$ are 0.97,
0.93, and 0.97 accordingly (Figure 2). Our model RVs and unsigned flux both
show strong correlation with the HARPS-N RV measurements with Spearman
correlation coefficients of 0.80 and 0.74 respectively. These strong
correlations show the reproducibility of the Haywood et al. (2016) and
Milbourne et al. (2019) results using our analysis pipeline.

Figure 2: Comparison between our pipeline-derived _SDO_ /HMI derived
velocities and observables (red) with reference calculations from Milbourne et
al. (2019) (blue). Correlation coefficients are shown in the upper left of
each plot. _First panel:_ Filling factor estimate, defined as the percentage
of magnetically active pixels on the solar surface. The filling factor, $f$,
for this period reflects the period of higher solar activity during which this
data was taken. _Second/Third panel:_ Comparison of unscaled velocity
components $\Delta\hat{v}_{\mathrm{conv}}$ and
$\Delta\hat{v}_{\mathrm{phot}}$. _Fourth panel:_ Total $\Delta
RV_{\mathrm{model}}$ from Milbourne et al. (2019) compared with our results,
showing excellent agreement.

##  3 Comparison of _SDO_ /HMI derived observables with ground-based
measurements

Our ultimate goal in calculating these solar observables is to gain insight
into the physical mechanisms driving measured RV variability seen in ground-
based Doppler measurements. We use data from the recently-commissioned NEID
instrument (Schwab et al., 2016), which has a dedicated solar feed that
delivers disk-integrated sunlight to the RV spectrometer (Lin et al. 2021
submitted.). NEID records high signal-to-noise (SNR$\sim$600) spectra every
$\sim$90 seconds throughout the day. NEID spectra are reduced using the
standard NEID pipeline, which delivers both integrated RVs and cross-
correlation functions (CCFs) for each frame recorded throughout the day
111https://neid.ipac.caltech.edu/search_solar.php. We filtered for days with
low cloud coverage, using data from the pyheliometer atop the NEID Solar
telescope (Lin et al. 2021 submitted.), and good instrumental drift
correction.

Using our _SDO_ analysis pipeline, we then computed the component RVs
($\Delta\hat{v}_{\mathrm{phot}}$ and $\Delta\hat{v}_{\mathrm{conv}}$) during
periods when NEID spectra were being collected. The independent amplitudes of
$\Delta\hat{v}_{\mathrm{phot}}$ and $\Delta\hat{v}_{\mathrm{conv}}$, are
significantly lower for the period of NEID data collection (December 2020 -
May 2021) than the values shown in Haywood et al. (2016) and Milbourne et al.
(2019), reflecting the low level of magnetic activity over the period for
which NEID has been observing the Sun (Table 1).

We find that the resultant ‘Sun-as-a-star’ _SDO_ /HMI computed RVs are largely
dominated by the convective velocity signal, likely due to the low level of
surface features (spots, plages, etc.) on the solar surface during the period
of analysis. This is further highlighted when comparing the modeled RVs and
convective velocity component with unsigned flux and filling factor, which are
all strongly correlated (Figure 3).

Velocity component | Velocity RMS | Haywood et al. (2016) | Milbourne et al. (2019)
---|---|---|---
$\Delta\hat{v}_{\mathrm{phot}}$ | .06 m s-1 |  $0.17$ m s-1 |  $0.21$ m s-1
$\Delta\hat{v}_{\mathrm{conv}}$ | 0.62 m s-1 |  $1.30$ m s-1 |  $0.88$ m s-1
Ground-based RVs | 0.8 m s-1 (NEID) | 3.12 $\mathrm{m\ s^{-1}}$ (HARPS) | 1.64 $\mathrm{m\ s^{-1}}$ (HARPS-N)

Table 1: RMS amplitudes for full time series of NEID data collection,
(December 2020 - May 2021) compared with archival results from Haywood et al.
(2016) and Milbourne et al. (2019). The $\Delta\hat{v}_{\mathrm{phot}}$ and
$\Delta\hat{v}_{\mathrm{conv}}$ values are derived from _SDO_ /HMI products.
The lower levels of solar activity during the period of NEID data collection
are readily apparent when comparing to the HARPS-N data from Haywood et al.
(2016) and Milbourne et al. (2019), which covered a more active span of the
Sun.

###  3.1 Comparison of SDO Model RVs with NEID Ground Based RVs

We then proceeded to compare the results of the _SDO_ /HMI measurements to
NEID ground-based RVs collected during instrument commissioning between
December 2020 and May 2021 (the nominal NEID commissioning period). Using
these data, we refit for the linear coefficients in equation 16 for the
convective and photometric velocity components. Table 2 lists our derived
scaling factors using the NEID RVs for calibration in comparison to the
scaling factors derived in Haywood et al. (2016) and Milbourne et al. (2019).
Both scaling factors characterize the impact of a measurement taken using one
line (_SDO_ /HMI data) to ground-based RVs measured across many lines, and we
expect these factors to change over time due to fluctuations in spot coverage
and activity level (Milbourne et al., 2019). Scaling factor A traces the
contribution of rotating active regions (primarily spots) to the bulk RVs. Our
time frame of interest had consistently low spot coverage, leading to a very
low amplitude for the photometric component, and thus large uncertainty in the
determination of scaling factor A (see Table 2). Scaling factor B is a measure
of the systematic difference in the convective blueshift due to varying
spectral line formation depths.

As established by Meunier et al. (2010a), Haywood et al. (2016), and Milbourne
et al. (2019) we expect the suppression of the convective blueshift to
dominate the overall RV, which is consistent with our measurements during this
time period (see Table 1). Our ground-based RVs from NEID do not show
measurable correlation with the photometric velocity signals calculated from
the _SDO_ /HMI images (Figure 3), which is consistent with the current phase
of solar activity (minimum) and in agreement with Lagrange et al. (2010). The
lack of sunspots during the NEID commissioning period leads to the low
variability in the photometric velocity, which further minimizes the
contribution of the photometric component to the overall model RVs ($\Delta
RV_{\mathrm{model}}$), as seen by the low value of scaling factor A in Table
1.

Figure 3: Comparison between _SDO_ /HMI derived observables and NEID RV
observations. _Top panel:_ Filling factor computed from _SDO_ /HMI
magnetograms. _Second panel:_ Unsigned magnetic flux computed from _SDO_ /HMI
magnetograms. _Third/Fourth panel:_ Unscaled basis velocity components
$\Delta\hat{v}_{\mathrm{conv}}$ and $\Delta\hat{v}_{\mathrm{phot}}$. _Fifth
panel:_ Total $\Delta RV_{\mathrm{model}}$ computed for the NEID observational
time span (in light lavender) with the ground-based NEID RVs filtered for bad
weather and cloud cover overlaid (in dark red). The average photon-limited
error on the daily binned $RV_{\mathrm{NEID}}$ measurements is $\sim$5
$\mathrm{cm\ s^{-1}}$. _Sixth panel:_ Residuals between $\Delta
RV_{\mathrm{model}}$ and $\Delta RV_{\mathrm{NEID}}$ from panel five. The
residual rms scatter is 78 $\mathrm{cm\ s^{-1}}$.

Parameter | This work | Haywood et al. (2016) | Milbourne et al. (2019)
---|---|---|---
A (photometric) | $1.10\pm 1.07$ | $2.45\pm 2.02$ | $2.24\pm 0.60$
B (convective) | $1.42\pm 0.11$ | $1.85\pm 0.27$ | $0.93\pm 0.11$
RV0 ($\gamma$) |  $-646.08\pm 0.38$ m s-1 |  $99.80\pm 0.28$ m s-1 |  $102.51\pm 0.06$ m s-1

Table 2: Scaling Factors for _SDO_ /HMI derived RVs using NEID data compared
with previous results from Haywood et al. (2016) and Milbourne et al. (2019)
using HARPS-N. In the period studied here, the overall activity level and
amplitude of the velocity components was lower along with the rms amplitude of
the RVs (see Table 1). We see a $\sim$20% difference between the $A$ and $B$
parameter values between the different time periods. This discrepancy is
likely due to the significantly lower level of activity in our recent data
from 2020-2021.

We are able to improve our results when restricting our analysis to days with
the most reliable ground-based observations. These dates have both excellent
observing weather (very low or no cloud coverage) and verified wavelength
solutions from the NEID laser frequency comb. For these dates we see very
strong correlation between our ground based RV measurements and both the model
RVs and unsigned magnetic flux. We find that unsigned magnetic flux is a
strong proxy for RV variation supporting the conclusion of Haywood et al.
(2020).

###  3.2 Active Region Area Dependence

In addition to the comparison of $\Delta RV_{\mathrm{model}}$ with the
convective velocity and ground-based measurements, we examine the area
dependence of active regions on the suppression of the convective blueshift.
To differentiate between large and small convective regions, we build an area
plot as a function of latitude, similar to Milbourne et al. (2019) and use 20
$\mu$Hem as our area threshold. We find network regions below the cutoff
across the solar disk, while larger plage and spot regions are found near the
solar equator at latitudes $\pm 30^{\circ}$. Milbourne et al. (2019)
implemented the same area cutoff and found similar latitude cuts of large
active regions, $0.75<\sin\Phi\leq 1.0$, where $\Phi$ is active region co-
latitude.

We then compute the $\Delta\hat{v}_{\mathrm{conv}}$ for the full time series
using the contributions from the small network regions, and large plage/spot
regions separately to compare with the $RV_{\mathrm{NEID}}$ from NEID. In
support of Milbourne et al. (2019) we find that large plage and spot regions
drive the RV variability more so than smaller network regions (Figure 6).
Additionally, we find that the power spectral density (PSD) contributions due
to small active regions do not contribute on timescales relevant to this study
(rotation timescales) as seen in Figure 4, similar to the results found in
Milbourne et al. (2019), while large active regions (plage/spots) show
rotational modulation as seen in Figure 4. While the convective velocity
component is calculated with the corrected Doppler map (rotational velocity
removed) there is still periodic structure on rotation timescales (see Figures
2 & 3). Large active magnetic regions (plage, spots) are expected to drive
variability on shorter timescales such as the period in this study, while the
network is expected to show more impact on longer timescales such as a Solar
Cycle (Milbourne et al., 2019).

Figure 4: Periodogram showing the contributions due to the small active
regions (green) and large active regions (blue) in comparison with model (red)
and ground-based RV measurements (black). The small active regions do not
contribute to power on rotation timescales, supporting results of Milbourne et
al. (2019), while large active regions (spots/plage) show measurable power at
periodicities near $P_{\mathrm{rot}}$.

###  3.3 Unsigned Magnetic Flux

Unsigned magnetic flux, $|\hat{B}_{\mathrm{obs}}|$, has been shown to be a
valuable proxy for stellar activity in the Sun (Haywood et al., 2020). Our
pipeline calculates unsigned magnetic flux using the methodology posed in
Haywood et al. (2016) and Haywood et al. (2020), which yields a high signal-
to-noise, independent activity metric to guide our analysis of ground-based
RVs. Figure 3 shows the calculated $|\hat{B}_{\mathrm{obs}}|$ time series
during the full NEID commissioning period. $|\hat{B}_{\mathrm{obs}}|$ is
determined by performing an intensity weighted sum of each pixel in the
magnetogram:

$|\hat{B}_{\mathrm{obs}}|=\frac{\sum_{ij}|B_{\mathrm{obs},\,ij}|I_{ij}}{\sum_{ij}I_{ij}}$
(17)

For our five month span (December 2021 - May 2021) of ‘best’ weather dates,
defined as days with no measurable cloud cover and minimal extinction at Kitt
Peak, we find a Spearman correlation coefficient of 0.43 between
$|\hat{B}_{\mathrm{obs}}|$ and measured NEID RVs and a correlation of 0.29
between measured NEID RVs and model RVs. We see strong correlation between
$|\hat{B}_{\mathrm{obs}}|$ and both the space-based RVs and
$\Delta\hat{v}_{\mathrm{conv}}$ ($\sim$ 0.90) throughout the entire time span,
emphasizing the dependence of magnetic activity on the suppression of the
convective blueshift. This is expected based on previous works such as Meunier
et al. (2010a), Meunier et al. (2010b), Haywood et al. (2016), and Haywood et
al. (2020).

These observations were taken during the solar minima of Solar Cycle 25 and
thus there was very little to no sunspot activity on the solar surface. The
convective velocity is primarily driven by large magnetic structures,
specifically plage regions. We see strong correlations between the plage
filling factor, convective velocity, and unsigned flux. Additionally, we find
the unsigned flux due to active regions strongly correlates with the _SDO_
/HMI model RVs.

Performing a linear regression of NEID RVs against measured unsigned magnetic
flux (see Figure 5), the NEID RV measurement RMS decreases from $\sim$80
$\mathrm{cm\ s^{-1}}$ to $\sim$60 $\mathrm{cm\ s^{-1}}$ over the five month
span (Figure 5) – an improvement of $\sim$50 $\mathrm{cm\ s^{-1}}$ in a
quadrature sum sense. This simplistic activity decorrelation implies there is
much promise for using unsigned magnetic flux to reduce RV variability, and we
await additional observations to ensure this correlation remains consistent as
the Sun enters a more active phase of the magnetic cycle.

Figure 5: Ground-based RV measurements before and after correction for
correlation with unsigned magnetic flux. For these ‘good’ weather days, where
there was low measurable cloud cover at WIYN, the correlation between
$RV_{\mathrm{NEID}}$ and $|\hat{B}_{\mathrm{obs}}|$ is 0.43, while the
correlation between $\Delta RV_{\mathrm{model}}$ and $RV_{\mathrm{NEID}}$ is
0.29, and the average daily binned error for $RV_{\mathrm{NEID}}$ is $\sim$5
$\mathrm{cm\ s^{-1}}$. We linearly fit the correlation between
$|\hat{B}_{\mathrm{obs}}|$ and $RV_{\mathrm{NEID}}$ and subtract this from
$RV_{\mathrm{NEID}}$ to reduce the scatter due to the unsigned magnetic flux.
The scatter reduces from $\sim$80 $\mathrm{cm\ s^{-1}}$ to $\sim$65
$\mathrm{cm\ s^{-1}}$ after this correction.

###  3.4 Filling Factor

The magnetic filling factor, defined as the fraction of the observed solar
disk that is active (corrected for foreshortening), is the second metric we
studied and compared to the NEID RVs:

$f=\frac{\sum_{ij}W_{ij}}{N_{\mathrm{pix}}}$ (18)

where $W_{ij}$ is the magnetic weighting array and $W_{ij}=1$ for active
pixels and 0 for pixels identified as inactive.

The Spearman correlation coefficient between $|\hat{B}_{\mathrm{obs}}|$ and
$f$ is 0.93 and we also find strong correlation between $f$ and
$\Delta\hat{v}_{\mathrm{conv}}$ (0.97), and $\Delta RV_{\mathrm{model}}$
(0.96) for the full five month time span. We compare the unsigned flux due to
active regions with the three filling factors calculated (see Figure 6). The
dominant feature driving the active region flux is the condensed faculae
regions known as plage. Plage are bright regions on the solar surface
typically found near sunspots that make up the majority of polarity in solar
active regions (Buehler et al., 2019). MHD simulations show that the amount of
thick flux tubes in plage regions is small, while flux tubes in the faculae
network expand more quickly than those in plage regions, and that the
continuum intensity of bright regions strongly correlates with magnetic field
strength (Röhrbein et al., 2011; Danilovic et al., 2013).

Figure 6: Comparison of faculae in condensed plage areas, faculae in the
diffuse network, and sunspot filling factors (%) with $\Delta
RV_{\mathrm{model}}$ (m s-1) and $|\hat{B}_{\mathrm{obs},\,\mathrm{active}}|$,
the unsigned flux due to active regions. We see the strongest correlation
between $f_{\mathrm{plage}}$ and $|\hat{B}_{\mathrm{obs},\,\mathrm{active}}|$
(Spearman coefficient of 0.97) and note the low percentage of sunspots on the
surface as this time span falls during the minimum of Solar Cycle 25. The
maximum sunspot percentage on the solar surface in this time span is 0.026%
while faculae cover between 1.5 - 3% of the solar disk.

##  4 Comparing _SDO_ measurements with NEID CCF metrics

In addition to comparing our space-based observables to integrated ground-
based RVs, we also compared them to RVs calculated using different spectral
line masks. These filtered masks were used to try and isolate the most and
least activity-sensitive features in the NEID solar spectrum. To test this, we
built several physically-motivated line masks to compare the effects of
magnetic activity on different groups of spectral lines using our _SDO_
-derived solar activity proxies.

We derived RVs using these tailored masks by calculating cross-correlation
functions (CCF) in the same manner as the standard NEID pipeline. These masks
were selected based on filtering for line species and line depth, and in all
cases our starting line list was the standard ESPRESSO G2 line mask
111https://www.eso.org/sci/software/pipelines/espresso/espresso-pipe-
recipes.html. Once computed, the CCFs were fit with Gaussian functions to
determine the averaged spectrum velocity for that particular mask.

###  4.1 Exploration of line depths and the relationship with magnetic
observables

The shifts in spectral line profile shape are primarily driven by convection
(Dravins et al., 1981) associated with velocity and intensity variations in
stellar granules (Hathaway et al., 2000). As the plasma from active regions
(sunspots, faculae, plage) interacts with the solar magnetic field, we see an
inhibition of this convective blueshift (see section 2.9) and this effect
varies with line depth (Gray, 2009). Deeper lines are less blue-shifted than
shallow lines (Gray, 2009) and thus we expect to see a stronger inhibition of
the convective blueshift from the mask built with shallow lines.

We study the effects of deep lines versus shallow lines to determine how these
correlate with both the bulk NEID RVs, as calculated by the NEID pipeline
using the ESPRESSO G2 mask, and unsigned magnetic flux values computed by
SolAster. We first applied a binary split between ‘deep’ lines and ‘shallow’
lines (deep lines having binary mask weights $\geq$ 0.5 and shallow having
weights $<$ 0.5) to see if the correlation between the NEID RVs and the
unsigned magnetic flux changed measurably (see Figure 7). The RVs calculated
using the mask of deep lines show strong correlation with NEID RVs (Spearman
correlation coefficient of 0.92), but show weaker correlation with unsigned
magnetic flux (correlation coefficient of 0.33). For shallow lines, we
similarly found a strong correlation with the NEID RVs (Spearman coefficient
of 0.98) and weak correlation with unsigned flux (correlation coefficient of
0.28). The rms ampltiude for the RVs calculated using the depth cut line mask
is $\sim$90 $\mathrm{cm\ s^{-1}}$ for the deep lines and $\sim$85 $\mathrm{cm\
s^{-1}}$ for the shallow lines. We also look at the residual scatter between
the integrated ground-based RVs (RVNEID) and the RVs from our curated masks,
finding rms scatter of $\sim$22 $\mathrm{cm\ s^{-1}}$ for deep lines and
$\sim$10 $\mathrm{cm\ s^{-1}}$ for shallow lines. We find that the RVs
calculated using either deep or shallow lines show similar correlations with
NEID RVs and unsigned flux, along with similar RMS scatter (Figure 7).

###  4.2 Derivation of ground-based RVs using physically-motivated masks

In addition to studying the variability as a function of line depth, we also
constructed masks based on previously published studies that isolate active
lines. To attempt to enhance the activity signature in the NEID spectra, we
adapted the line list from Wise et al. (2018) of activity-sensitive lines and
recomputed the NEID RVs. This list contains only those lines that showed
significant correlations between their line depths and the chromospheric Ca II
H&K index for active stars (Wise et al., 2018). We then compared these RVs
with our _SDO_ /HMI pipeline measurements, as well as the ‘standard’ NEID
pipeline RVs. We found the RVs of these select $\sim$20 lines (those found in
the Wise list and ESPRESSO G2 mask) correlate well with the NEID pipeline RVs,
and do not show a strong correlation with the _SDO_ -derived unsigned magnetic
flux. This may be unsurprising, given that Wise et al. (2018) built this list
of activity-sensitive lines through observations of stars that are more active
than the Sun’s current activity level and show significantly higher
rotationally-modulated behavior. As such, the observed relationship between
activity and line variability may not hold nor be measurable during this time
of low solar activity.

Figure 7 shows the correlation between the RVs derived using different line
masks and the _SDO_ -derived unsigned magnetic flux. The left panel shows the
correlations for the activity-sensitive line list (Wise et al., 2018) while
the second column displays the correlations for the mask built using only Fe I
lines within the ESPRESSO G2 mask. We filtered specifically for Fe I lines to
compare against the SDO-derived RVs, which use a single Fe line to compute the
majority of the observables used in SolAster. This Fe I line list contains
significantly more features than the mask constructed using features from Wise
et al. (2018), yet we find similar scatter in the RV time series in both
masks, implying neither list is significantly more activity sensitive. We do
find that the integrated RVs from the Fe I lines show mildly stronger
correlation with both the bulk NEID RVs and unsigned magnetic flux in
comparison to the activity-sensitive lines from Wise et al. (2018).

Figure 7: Comparison of NEID RVs derived using different cross-correlation
line masks (columns) compared to the _SDO_ unsigned magnetic flux (top row)
and the bulk NEID RVs (bottom row). Spearman correlation coefficients are
shown in the upper left in each panel along with fitted slope values and
uncertainties on both measurements. _First column:_ Correlations for RVs
calculated using activity-sensitive lines from Wise et al. (2018). There is
moderately strong correlation with the NEID pipeline RVs, but only a weak
correlation with unsigned magnetic flux. _Second column:_ Correlation for RVs
derived using only Fe I lines in the ESPRESSO G2 mask, showing stronger
correlation with both the bulk NEID RVs and the unsigned magnetic flux. _Third
column:_ Correlation for RVs derived using deep lines in the ESPRESSO G2 mask,
showing very strong correlation with the bulk RVs. _Fourth column:_
Correlation for RVs derived using shallow lines in the ESPRESSO G2 mask,
showing the strongest correlation with the bulk RVs. All error bars report
only photon noise, excluding any instrument systematics or additional stellar
jitter (5-7 $\mathrm{cm\ s^{-1}}$).

###  4.3 Comparing measurables from NEID Cross Correlation Functions

Beyond integrated RV measurements, we also studied an assortment of CCF
parameters to determine which metrics best trace magnetic field variability.
The primary CCF metrics we compared are: fitted amplitude, full width at half
maximum (FWHM), skew, and integrated area below the line profile, in addition
to calculating the RV shift via Gaussian fit of our CCF. Comparing the
variation of these metrics with magnetic observables over time provides
insight into the effects driving the CCF shape changes we observe.

Figure 8 shows the full suite of computed CCF measurements compared to the
SolAster data products. We looked at CCF metrics using the full G2 ESPRESSO
mask which covers wavelengths from 3700 - 7900 Å to calculate RV variations
and the list of metrics outlined previously. We find that certain CCF metrics
serve as better proxies for magnetic activity due to their strong correlation
with unsigned magnetic flux. The SDO data allow us to isolate days with higher
magnetic activity, namely days with larger magnetic filling factors, and
compare these results to days of ‘quiet-Sun’.

Figure 8: Corner plot showing correlation between space-based magnetic
observables and ground-based CCF metrics derived using the full G2 ESPRESSO
mask. The time frame for the dataset shown is December 2020 - May 2021. Days
of high activity (measurable spots) are shown as dark pink dots and days of
quiet-Sun are shown in light lavender. Correlation coefficients are listed in
the legends in each sub-plot for the full dataset (top value), high activity
days (middle), and quiet-Sun days (bottom), respectively. When comparing the
space and ground-based data, we see the strongest correlation between unsigned
magnetic flux, $|\hat{B}_{\mathrm{obs}}|$, and the CCF amplitude and
integrated area measurements.

We find the strongest proxies for magnetic activity to be integrated area and
amplitude, each showing moderate correlation coefficients ($\sim$0.5, see
Figure 8) when compared to $|\hat{B}_{\mathrm{obs}}|$. This broadly supports
the results of Costes et al. (2021), specifically in regards to the CCF FWHM
and integrated area measurements, though we note the base level of solar
activity is different between these two studies. The moderately strong
negative correlation between integrated area and unsigned magnetic flux is
also supported by the work of Collier Cameron et al. (2019) and Costes et al.
(2021), where integrated area was found to be a strong tracer for the
evolution of the magnetic network due to the variation in CCF area showing
little to no rotational modulation. The lack of rotational modulation of the
integrated area metric implies that its variation is likely driven by
axisymmetrically distributed structures over the solar surface, specifically
from the circulation of dispersed magnetic flux elements from regions that
were once active (Collier Cameron et al., 2019). For days with sunspots, the
correlation between photometric velocity and unsigned flux increases while the
correlation between convective velocity and unsigned flux decreases. This
shows the effect of the rotationally modulated spots on the photometric
velocity component as discussed in section 2.8 and supports the correlation
between spot factor and jumps in unsigned active magnetic flux as seen in
Figure 6.

Figure 9: Comparison of the two CCF metrics (amplitude and integrated area)
that correlate most strongly with unsigned magnetic flux as a function of line
depth (deep having CCF mask weights $\geq$ 0.5, and shallow having weights $<$
0.5). We show a comparison of the percent change of these metrics against the
unsigned magnetic flux. The correlation coefficients are listed in the legend
for the full dataset (top value), high activity days (middle), and quiet-Sun
days (bottom), respectively. The dark red dots show days of high activity,
while the light lavender points are quiet-Sun days.

To explore the activity behavior as a function of spectral line depth, we
created several CCF masks using different depth cuts. As previously described
(section 4.1), we use the CCF mask line ‘weights’ listed in the ESPRESSO mask
as proxies for relative depth. We then recompute the CCF and resultant RVs for
all lines within each depth cut. We find that across all depth masks, we see a
strong correlation between the derived RVs and bulk NEID pipeline RV
measurements (indicating the majority of lines are shifting in a similar
manner). For CCF metrics derived with the shallow line mask, we see stronger
correlation between unsigned flux and both amplitude and integrated area in
comparison to the metrics calculated with the deep line mask as shown in
Figure 9. In addition, our observed correlations broadly supports previous
results from Meunier et al. (2017), and Reiners et al. (2016) of stronger
correlation between the CCF RVs from only the shallow lines and the calculated
convective velocity component.

##  5 Application of SolAster to archival planetary transits

Since the launch of _SDO_ in 2010, three planetary solar transits have been
imaged by _SDO_ instruments – Mercury in 2016 and 2019, and Venus in 2012.
Here we detail our attempts at recovering the corresponding
Rossiter–McLaughlin (RM) signals due to the transits of Venus and Mercury. The
measurement of these transits serve as proof of concept of using spatially
resolved disk images to calculate disk-integrated RVs at precisions currently
unattainable from the ground. We apply our SolAster pipeline methods to these
time frames in the attempt to 1) recover the effect of the planetary transit
on the RV measurements and 2) empirically estimate the precision floor of our
constructed ‘Sun-as-a-star’ RVs. Beyond testing our pipeline, these RV
measurements also showcase the magnitude of the effect that stellar activity
has on our ability to detect small RV signals over short timescales.

For the all three transits, we used _SDO_ /HMI data products at a two minute
cadence and computed the unsigned magnetic flux, convective and photometric
velocities. We then reconstruct the overall model RV variation based on the
parameterization outlined in Section 2.7.

###  5.1 Mercury Transits

Since the start of the _SDO_ observation period, there have been two Mercury
transits on May 9, 2016 and on November 11, 2019. Using SolAster we attempt to
recover the RM signal of these planetary transits to understand the noise
floor (both astrophysical and instrumental) of RV measurements due to the
extremely low amplitude expected from these transits.

####  5.1.1 Mercury 2016 Transit

The RV signal induced during the Mercury 2016 transit was expected to be on
the $\sim$5 $\mathrm{cm\ s^{-1}}$ level, and our constructed RV model from
_SDO_ /HMI (using the weighting factors in Table 2) shows significantly higher
variability. Through observations of the Sun, we have established that solar
RV variations are driven by bright faculae (in regions of concentrated plage).
From long-term surveys, such as the Mt Wilson HK project (Baliunas et al.,
1988), we know the surfaces of old, slowly rotating Sun-like stars are
faculae-dominated (Radick et al., 1983; Lockwood et al., 1984). These surveys
monitored the optical photometric variations and the Ca II H&K over decades
for FGK stars noting that brightness increases as a function of activity, just
as observed for the Sun throughout its magnetic cycle. Therefore, these stars
are faculae rather than spot dominated and we expect RV variations to be
driven by the suppression of the convective blueshift, as seen for the Sun.
This is generally consistent with our observed $\Delta\hat{v}_{\mathrm{conv}}$
in Figure 10, which has significantly more scatter
than$\Delta\hat{v}_{\mathrm{phot}}$. This result serves as a glaring example
of the ability of stellar activity to degrade RV sensitivity to planetary
signals, even over short timescales.

Figure 10: Calculated RVs for the eight hour Mercury transit on May 9, 2016.
_Top panel:_ Unsigned magnetic flux over the period of the transit. _Second
panel:_ Photometric velocity from the 2016 transit. _Third panel:_ Convective
velocity from the 2016 transit. _Fourth panel:_ Model RV variation constructed
from _SDO_ /HMI data during the planet’s transit. For all panels, the red
diamonds show the 5.25 minute binned RVs (averaging out p-mode oscillations)
while the background grey signal shows the 2 minute cadence of our
calculations. Even with the low-level of underlying surface activity, the
integrated RV model is dominated by $\Delta\hat{v}_{\mathrm{conv}}$.

####  5.1.2 Mercury 2019 Transit

Identical to the construction of the model RVs for the Mercury 2016 transit,
we also look at the November 11, 2019 transit of Mercury. This transit period
was clear of sunspots and occurred while the Sun moved out of the absolute
solar minimum of Solar Cycle 24 which occurred in October 2019, just prior to
the transit. These low activity conditions provided an ideal background for
the recovery attempt of this $\sim$5 $\mathrm{cm\ s^{-1}}$ signal. We show our
recovery attempt along with the component velocities in Figure 11. This
transit occurred at an even lower level of solar activity than the 2016
Mercury transit, but the overall RV variability is still largely dominated by
the convective velocity component.

Figure 11: Calculated RVs and magnetic observables for the Mercury transit on
November 11, 2019. _Top panel:_ Unsigned magnetic flux over the period of the
transit showing the low level of activity during this time period. _Second
panel:_ Photometric velocity component from the Mercury 2019 transit. _Third
panel:_ Convective velocity component from the Mercury 2019 transit. _Fourth
panel:_ Model RV variation constructed from _SDO_ /HMI data during the
planet’s transit using the SolAster pipeline as outlined in Section 2.7. For
all panels, the red diamonds show the 5.25 minute binned RVs (averaging out
p-mode oscillations) while the background grey signal shows the 2 minute
cadence of our calculations.

###  5.2 Venus 2012 Transit

The Venus transit occurred during a period of high activity in the Solar
Cycle, however this specific day had relatively low solar activity with a
small sunspot and facular filling factor which remained consistent throughout
the transit period. The transit occurred from 22:09 UTC on June 5, 2012 until
04:49 UTC on June 6, 2012. Using the SolAster pipeline, we calculate the RV
components and reconstruct the model RV variation as outlined in Section 2.10.
We show the overall velocity signal is largely dominated by the RM signal,
rather than the convective velocity component (see Figure 12), allowing for a
clean recovery of the RM waveform.

Figure 12: Calculated RVs for the Venus transit from June 5 to June 6, 2012.
_Top panel:_ Unsigned magnetic flux over the period of the transit. _Second
panel:_ Photometric velocity component from the Venus 2012 transit. _Third
panel:_ Convective velocity component from the Venus 2012 transit. _Fourth
panel:_ Model RV variation constructed from _SDO_ /HMI data during the
planet’s transit using the SolAster pipeline as outlined in Section 2.7. For
all panels, the red diamonds show the 5.25 minute binned RVs (averaging out
p-mode oscillations) while the background grey signal shows the 2 minute
cadence of our calculations.

##  6 Discussion

Our Python based, publicly available _SDO_ /HMI analysis pipeline (SolAster)
allows us to calculate both magnetic observables and ‘Sun-as-a-star’ RV
variations using space-based data for comparison with ground-based
measurements. Moving forward, these data products will aid in studies aimed at
deriving new stellar activity indicators explorations in ground-based spectra.
By looking at correlations between the space-based data and ground-based
measurements from RV facilities such as NEID, we hope to find stronger proxies
for stellar activity in Sun-like stars, which would help to improve planet
detection sensitivity in future RV surveys. Leveraging the _SDO_ /HMI and NEID
data, there are a variety of paths we aim to explore in future studies,
including:

  * •

Comparing line-by-line (Dumusque, 2018) and integrated CCF metrics with space-
based observables, and more closely explore the metrics that show the most
promise as activity proxies (CCF integrated area and amplitude being two
examples that show promise based on our preliminary study).

  * •

Exploring the wavelength dependence of correlations between NEID measurements
and _SDO_ -pipeline calculations both in RVs and CCFs. This could allow us
determine whether the line shape variability is driven by Zeeman effects
(Reiners et al., 2013). If this is the case, we would expect to see a stronger
correlation between unsigned magnetic flux and the CCF depth/integrated area
in the redder portions of the spectrum. If instead the variation is driven by
the inhibition of the convective blueshift, then we expect to see stronger
correlation in the blue lines (Reiners et al., 2013).

  * •

Detrending NEID RVs against the unsigned magnetic flux
($|\hat{B}_{\mathrm{obs}}|$) using more complex parameterizations: including
the FF’ (Aigrain et al., 2012) technique and Gaussian Processes (Haywood et
al., 2014) and others. While our preliminary linear detrending (Figure 5) did
improve the scatter for a subset of the NEID data, a more robust use of the
_SDO_ data, looking at the solar rotation period and life time of sunspots, to
detrend the RVs could further decrease the activity signal.

  * •

Revisiting our analysis on times with heightened solar activity. The time span
used in this study (December 2020 - May 2021) was during a period of low solar
activity. While there were times with sunspots and slightly increased activity
levels, the overall activity level was very low, both for the Sun itself and
in comparison with other Solar-type stars. Studying periods of higher solar
activity (which we are entering now) that have complimentary ground-based
spectra (e.g. 2013-2015, with HARPS-N (Dumusque et al., 2021)) would be useful
for quantifying how the level of solar activity affects the correlations
between the ground-based metrics and space-based observables. This would also
allow us to fine tune our scaling factors (Table 2) and derive more precise
model RVs. This would also allow us to better understand how the amplitudes
and relative contributions of the two velocity components,
$\Delta\hat{v}_{\mathrm{conv}}$ and $\Delta\hat{v}_{\mathrm{phot}}$, vary as a
function of time and activity level.

##  7 Conclusion

Using _SDO_ /HMI data products and ground-based spectra, we studied the ‘Sun-
as-a-star’ to estimate RV variations due to forms of solar activity. In doing
so, we developed a standalone solar data analysis package, SolAster, using the
methods outlined by Haywood et al. (2016) and Milbourne et al. (2019) to
compute model RVs, and validated our results against these previously
published datasets. We also calculated additional magnetic observables to
search for any correlations between these measurements and ground-based RVs
from NEID. We found that while the RV variation is driven by the combination
of convective blueshift suppression and the rotational Doppler imbalance due
to sunspots and plage, the dominant component is due to the convective
blueshift suppression, confirming previous results by Meunier et al. (2010b)
and others. We found the RMS scatter of both $\Delta\hat{v}_{\mathrm{conv}}$
and $\Delta\hat{v}_{\mathrm{phot}}$ are lower than the results of Haywood et
al. (2016) and Milbourne et al. (2019) due to our study taking place during a
period of minimal solar activity.

We found that plage regions are the dominant driver of the observed magnetic
flux and convective blueshift suppression, as seen in Figure 6, and the
overall RV variation is dominated by the convective blueshift suppression. The
photometric contribution, primarily affected by bright active regions and
sunspots, is quite minimal and strongly dependent on sunspot filling factor.
These conclusions support previous work on this topic by Meunier et al.
(2010b), Haywood et al. (2016), Milbourne et al. (2019), and Haywood et al.
(2020).

Filtering the NEID data for only days with optimum observing conditions, we
find a strong correlation between the ground-based RVs, model RVs, and
unsigned magnetic flux. There is a strong correlation between
$B_{\mathrm{obs}}$ and the NEID RVs, which we removed via a linear
decorrelation to improve the RMS scatter down to the $\sim$60 $\mathrm{cm\
s^{-1}}$ level.

To better understand which surface features drive active region magnetic flux,
we compared a variety of activity indicators and proxies for solar activity.
We found that large facular regions known as plage are the dominant component
driving the temporal flux variation. Additionally, we found a strong
correlation between these magnetic observables and RV variations, aligning
with current work in this field (Haywood et al., 2016, 2020). This provides
additional evidence for the exploration of unsigned magnetic flux as a proxy
for stellar activity, especially in Sun-like stars.

We also investigated correlations between spectral line shape and magnetic
activity indicators. We built a variety of physically motivated line masks to
better quantify the affect of magnetic activity on different spectral line
parameters. We found that RVs calculated using masks built from either deep
lines or shallow lines show similar correlation with both ground-based NEID RV
pipeline measurements and unsigned magnetic flux in comparison to RVs from
shallow lines.

Using SolAster, we are able to recover the planetary RM signal for the Venus
2012 transit at high SNR (see Figure 12). However, the planetary RM signals in
both the 2016 and 2019 Mercury transits are dwarfed by activity signals, which
dominate the RV noise floor over during transit. While the 2019 Mercury
transit took place during a period of low solar activity, we are still unable
to measure the RM signal due to colluding noise from solar activity at levels
$\sim$2-4$\times$ the RM signal. We used these transit events to both
empirically gauge the noise floor of our integrated RV measurements, as both
signals are expected to be of low amplitude ($<$ 50 $\mathrm{cm\ s^{-1}}$),
and study activity signals at short ($<$8 hour) timescales, finding that even
at low activity levels we were unable to recover the RM signal due to the
Mercury transit highlighting the incredibly small amplitude of the signal.

While a simple linear detrending with unsigned flux enabled a measurable
improvement in a subset of NEID RVs down to the $\sim$60 $\mathrm{cm\ s^{-1}}$
level, there are still additional avenues of exploration required in order to
reach the $\sim$10 $\mathrm{cm\ s^{-1}}$ sensitivity needed to detect Earth-
like planets orbiting Sun-like stars. We are able to recover a planetary
transit (RM signal) with an amplitude on the 10’s of $\mathrm{cm\ s^{-1}}$
level, showing the promise of using this method to improve our understanding
of stellar activity and its effect on RV measurements at a variety of
timescales. With more ground-based observations from NEID and the upcoming
increase in solar activity, we will be able to continue improving our modeling
methods and more finely hone our understanding of how specific types of
stellar activity affect ground-based high resolution spectra.

##  8 Acknowledgements

The research was carried out at the Jet Propulsion Laboratory, California
Institute of Technology, under a contract with the National Aeronautics and
Space Administration (80NM0018D0004). This research used version 3.1 (SunPy
Community et al., 2020) of the SunPy open source software package.

Based on observations at Kitt Peak National Observatory, National Optical
Astronomy Observatory, which is operated by the Association of Universities
for Research in Astronomy (AURA) under a cooperative agreement with the
National Science Foundation. These results are based on observations obtained
with NEID on the WIYN 3.5m Telescope. WIYN is a joint facility of the
University of Wisconsin–Madison, Indiana University, NSF’s NOIRLab, the
Pennsylvania State University, Purdue University, University of California,
Irvine, and the University of Missouri. The authors are honored to be
permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a
mountain with particular significance to the Tohono O’odham.

## References

  * Aigrain et al. (2012) Aigrain, S., Pont, F., & Zucker, S. 2012, MNRAS, 419, 3147, doi: 10.1111/j.1365-2966.2011.19960.x
  * Allen (1973) Allen, C. W. 1973, Astrophysical quantities
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Baliunas et al. (1988) Baliunas, S. L., Bennett, W., Briggs, J., et al. 1988, in Bulletin of the American Astronomical Society, Vol. 20, 994
  * Boisse et al. (2011) Boisse, I., Bouchy, F., Hébrard, G., et al. 2011, A&A, 528, A4, doi: 10.1051/0004-6361/201014354
  * Buehler et al. (2019) Buehler, D., Lagg, A., van Noort, M., & Solanki, S. K. 2019, A&A, 630, A86, doi: 10.1051/0004-6361/201833585
  * Carrington (1859) Carrington, R. C. 1859, MNRAS, 19, 81, doi: 10.1093/mnras/19.3.81
  * Collier Cameron et al. (2019) Collier Cameron, A., Mortier, A., Phillips, D., et al. 2019, MNRAS, 487, 1082, doi: 10.1093/mnras/stz1215
  * Costes et al. (2021) Costes, J. C., Watson, C. A., de Mooij, E., et al. 2021, Monthly Notices of the Royal Astronomical Society, 505, 830–850, doi: 10.1093/mnras/stab1183
  * Crass et al. (2021) Crass, J., Gaudi, B. S., Leifer, S., et al. 2021, arXiv e-prints, arXiv:2107.14291.  https://arxiv.org/abs/2107.14291
  * Danilovic et al. (2013) Danilovic, S., Röhrbein, D., Cameron, R. H., & Schüssler, M. 2013, A&A, 550, A118, doi: 10.1051/0004-6361/201219726
  * Dravins (1990) Dravins, D. 1990, Mem. Soc. Astron. Italiana, 61, 513
  * Dravins et al. (1981) Dravins, D., Lindegren, L., & Nordlund, A. 1981, A&A, 96, 345
  * Dumusque (2018) Dumusque, X. 2018, A&A, 620, A47, doi: 10.1051/0004-6361/201833795
  * Dumusque et al. (2011) Dumusque, X., Udry, S., Lovis, C., Santos, N. C., & Monteiro, M. J. P. F. G. 2011, A&A, 525, A140, doi: 10.1051/0004-6361/201014097
  * Dumusque et al. (2021) Dumusque, X., Cretignier, M., Sosnowska, D., et al. 2021, A&A, 648, A103, doi: 10.1051/0004-6361/202039350
  * Fischer et al. (2016) Fischer, D. A., Anglada-Escude, G., Arriagada, P., et al. 2016, PASP, 128, 066001, doi: 10.1088/1538-3873/128/964/066001
  * Fligge et al. (2000) Fligge, M., Solanki, S. K., & Unruh, Y. C. 2000, A&A, 353, 380
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Gibson et al. (2018) Gibson, S. R., Howard, A. W., Roy, A., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, ed. C. J. Evans, L. Simard, & H. Takami, 107025X, doi: 10.1117/12.2311565
  * Gray (2009) Gray, D. F. 2009, ApJ, 697, 1032, doi: 10.1088/0004-637X/697/2/1032
  * Hathaway et al. (2000) Hathaway, D. H., Beck, J. G., Bogart, R. S., et al. 2000, Sol. Phys., 193, 299, doi: 10.1023/A:1005200809766
  * Hatzes (2016) Hatzes, A. P. 2016, The Radial Velocity Method for the Detection of Exoplanets, ed. V. Bozza, L. Mancini, & A. Sozzetti, Vol. 428, 3, doi: 10.1007/978-3-319-27458-4_1
  * Haywood et al. (2014) Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517, doi: 10.1093/mnras/stu1320
  * Haywood et al. (2016) Haywood, R. D., Collier Cameron, A., Unruh, Y. C., et al. 2016, Monthly Notices of the Royal Astronomical Society, 457, 3637, doi: 10.1093/mnras/stw187
  * Haywood et al. (2020) Haywood, R. D., Milbourne, T. W., Saar, S. H., et al. 2020, arXiv e-prints, arXiv:2005.13386.  https://arxiv.org/abs/2005.13386
  * Hoeksema et al. (2018) Hoeksema, J. T., Baldner, C. S., Bush, R. I., Schou, J., & Scherrer, P. H. 2018, Sol. Phys., 293, 45, doi: 10.1007/s11207-018-1259-8
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jurgenson et al. (2016) Jurgenson, C., Fischer, D., McCracken, T., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, 99086T, doi: 10.1117/12.2233002
  * Kanodia & Wright (2018) Kanodia, S., & Wright, J. 2018, Research Notes of the American Astronomical Society, 2, 4, doi: 10.3847/2515-5172/aaa4b7
  * Lagrange et al. (2010) Lagrange, A. M., Desort, M., & Meunier, N. 2010, A&A, 512, A38, doi: 10.1051/0004-6361/200913071
  * Lockwood et al. (1984) Lockwood, G. W., Thompson, D. T., Radick, R. R., et al. 1984, PASP, 96, 714, doi: 10.1086/131410
  * Mayor & Queloz (1995) Mayor, M., & Queloz, D. 1995, Nature, 378, 355, doi: 10.1038/378355a0
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. S. van der Walt & J. Millman, 51 – 56
  * Meunier et al. (2010a) Meunier, N., Desort, M., & Lagrange, A. M. 2010a, A&A, 512, A39, doi: 10.1051/0004-6361/200913551
  * Meunier & Lagrange (2020) Meunier, N., & Lagrange, A. M. 2020, A&A, 638, A54, doi: 10.1051/0004-6361/201937354
  * Meunier et al. (2015) Meunier, N., Lagrange, A. M., Borgniet, S., & Rieutord, M. 2015, A&A, 583, A118, doi: 10.1051/0004-6361/201525721
  * Meunier et al. (2010b) Meunier, N., Lagrange, A. M., & Desort, M. 2010b, A&A, 519, A66, doi: 10.1051/0004-6361/201014199
  * Meunier et al. (2017) Meunier, N., Mignon, L., & Lagrange, A. M. 2017, A&A, 607, A124, doi: 10.1051/0004-6361/201731017
  * Milbourne et al. (2019) Milbourne, T. W., Haywood, R. D., Phillips, D. F., et al. 2019, ApJ, 874, 107, doi: 10.3847/1538-4357/ab064a
  * Palle et al. (1995) Palle, P. L., Jimenez, A., Perez Hernandez, F., et al. 1995, ApJ, 441, 952, doi: 10.1086/175414
  * Pepe et al. (2021) Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96, doi: 10.1051/0004-6361/202038306
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3, doi: 10.1007/s11207-011-9841-3
  * Radick et al. (1983) Radick, R. R., Mihalas, D., Lockwood, G. W., et al. 1983, PASP, 95, 621, doi: 10.1086/131229
  * Reiners et al. (2016) Reiners, A., Mrotzek, N., Lemke, U., Hinrichs, J., & Reinsch, K. 2016, Astronomy & Astrophysics, 587, A65
  * Reiners et al. (2013) Reiners, A., Shulyak, D., Anglada-Escudé, G., et al. 2013, A&A, 552, A103, doi: 10.1051/0004-6361/201220437
  * Robertson et al. (2015) Robertson, P., Roy, A., & Mahadevan, S. 2015, ApJ, 805, L22, doi: 10.1088/2041-8205/805/2/L22
  * Röhrbein et al. (2011) Röhrbein, D., Cameron, R., & Schüssler, M. 2011, A&A, 532, A140, doi: 10.1051/0004-6361/201117090
  * Saar & Donahue (1997) Saar, S. H., & Donahue, R. A. 1997, ApJ, 485, 319, doi: 10.1086/304392
  * Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, The Helioseismic and Magnetic Imager (HMI) Investigation for the Solar Dynamics Observatory (SDO), ed. P. Chamberlin, W. D. Pesnell, & B. Thompson (New York, NY: Springer US), 207–227, doi: 10.1007/978-1-4614-3673-7_10
  * Schou et al. (1998) Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390, doi: 10.1086/306146
  * Schröter (1985) Schröter, E. H. 1985, Solar Physics, 100, 141, doi: 10.1007/BF00158426
  * Schwab et al. (2016) Schwab, C., Rakich, A., Gong, Q., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, 99087H, doi: 10.1117/12.2234411
  * Snodgrass (1984) Snodgrass, H. B. 1984, Solar Physics, 94, 13, doi: 10.1007/BF00154804
  * Snodgrass & Ulrich (1990) Snodgrass, H. B., & Ulrich, R. K. 1990, ApJ, 351, 309, doi: 10.1086/168467
  * SunPy Community et al. (2020) SunPy Community, Barnes, W. T., Bobra, M. G., et al. 2020, ApJ, 890, 68, doi: 10.3847/1538-4357/ab4f7a
  * Thompson et al. (2016) Thompson, S. J., Queloz, D., Baraffe, I., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, 99086F, doi: 10.1117/12.2232111
  * Thompson (2006) Thompson, W. T. 2006, A&A, 449, 791, doi: 10.1051/0004-6361:20054262
  * Ulrich & Boyden (2006) Ulrich, R. K., & Boyden, J. E. 2006, Sol. Phys., 235, 17, doi: 10.1007/s11207-006-0041-5
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wise et al. (2018) Wise, A. W., Dodson-Robinson, S. E., Bevenour, K., & Provini, A. 2018, The Astronomical Journal, 156, 180, doi: 10.3847/1538-3881/aadd94
  * Wright (2018) Wright, J. T. 2018, Radial Velocities as an Exoplanet Discovery Method, ed. H. J. Deeg & J. A. Belmonte (Cham: Springer International Publishing), 619–631, doi: 10.1007/978-3-319-55333-7_4
  * Yeo et al. (2013) Yeo, K. L., Solanki, S. K., & Krivova, N. A. 2013, A&A, 550, A95, doi: 10.1051/0004-6361/201220682
  * Zhao et al. (2016) Zhao, J., Stejko, A., & Chen, R. 2016, Solar Physics, 291, 731, doi: 10.1007/s11207-016-0864-7