EXPRES. II. Searching for Planets around Active Stars: A Case Study of HD 101501

Fairborn Observatory's Automatic Photoelectric Telescopes (APTs; Henry 1999) observed HD 101501 for 28 seasons (spanning 1993 April 18–2020 June 21). Observations have about 1 day typical cadence for 6–7 months each year, totalling 2673 data points. The data display a significant correlated structure arising from stellar activity in the target and have a standard deviation of 7.1 mmag in the V band. We focus on activity from spots, which are modulated by the rotational period, and remove long-term brightness variations from the light curve. The long-term trend of the light curve that cannot be simply attributed to rotational modulation is removed by smoothing the light curve over 100 days and subtracting the resultant trend. The light-curve standard deviation is reduced to 4.5 mmag, following detrending (Figure 1). The periodogram has complex structure around 17 days, close to the 16.18 day rotation period estimated by Donahue et al. (1996). The maximum power is at 17.51 days.

Zoom In Zoom Out Reset image size

TESS (Ricker et al. 2014) observed HD 101501 during Sector 22 (2020 February 18–2020 March 18; Figure 2). The 2 minute cadence, simple aperture photometry (SAP) light curve was used and the first five cotrending basis vectors were applied to remove instrumental signatures in order to preserve stellar astrophysics (following the process used in Roettenbacher & Vida 2018). Both the light curve and cotrending basis vectors were obtained through the Barbara A. Mikulski Archive for Space Telescopes (MAST).

Zoom In Zoom Out Reset image size

We analyze 76 RVs (over 33 distinct nights) of HD 101501 obtained between 2018–2020. Data were collected with EXPRES commissioned at the 4.3 m Lowell Discovery Telescope (Levine et al. 2012; observing program: The 100 Earths Survey). EXPRES achieves ∼30 cm s⁻¹ measurement precision for a pixel signal-to-noise ratio of 250 at 5500 Å. EXPRES has typical resolving power of R ∼ 137,500 and spans a wavelength range of 380–780 nm. More details regarding EXPRES may be found in recent studies investigating performance benchmarks and detailing the RV extractions (Blackman et al. 2020; Brewer et al. 2020; Petersburg et al. 2020). We also extract activity indicators (CCF FWHM, CCF BIS, and Hα equivalent width) for each exposure, derived from the spectrum and CCF. The 100 Earths Survey targets chromospherically quiet stars without close-in gas giants in order to search for terrestrial planets. HD 101501 represents one of several additional, active stars observed for purposes of investigating and mitigating activity signals. Physical attributes of HD 101501 are listed in Table 1, and summary statistics for the EXPRES exposures are in Table 2. For convenience, we subtract the mean of the RV data (about −5.55 km s⁻¹) corresponding to the systemic velocity.

Table 1.  Physical Properties of HD 101501

Property	Symbol	Units	Valuea
Visual Magnitude	V	mag.	5.31
Distance	d	pc	9.61
Effective Temperature	Teff	K	5502
Surface Gravity	log g	cm s−2	4.52
Metallicity	[Fe/H]	[Fe/H]⊙	−0.04
Age	tage	Gyr.	${3.5}_{-2.2}^{+2.8}$
Proj. Rotation Speed	$v\sin i$	km s−1	2.2
Luminosity	$\mathrm{log}$ L*	$\mathrm{log}$ L⊙	−0.21 ± 0.02
Radius	R*	R⊙	0.86 ± 0.02
Mass	M*	M⊙	0.90 ± 0.12

Note.

^aAll values and uncertainties are as reported by Brewer et al. (2016).

Download table as:  ASCIITypeset image

Table 2.  Summary of EXPRES Radial Velocity Measurements Used in This Study

Property	Symbol	Units	Value
Number of Exposures	Nexp	⋯	76
Number of Nights	Nnight	⋯	33
Time Baseline	⋯	days	796
Average RV	$\widehat{\mathrm{RV}}$	m s−1	−5554.43
RV rms	rms	m s−1	6.1
Median Measurement Uncertainty	${\tilde{\sigma }}_{n}$	m s−1	0.38

Download table as:  ASCIITypeset image

The stochastic nature of stellar activity makes it difficult to model analytically, and has motivated the use of GPs (Haywood et al. 2014; Rajpaul et al. 2015). A GP is a flexible model which assumes that data points are drawn from a multivariate normal distribution (Rasmussen & Williams 2006). Under a GP model, RV measurements y at times t have the joint distribution

$$\begin{eqnarray}&&p({\boldsymbol{y}};{\boldsymbol{t}})\sim { \mathcal N }(m({\boldsymbol{t}}),K({\boldsymbol{t}},{\boldsymbol{t}})+{\sigma }_{n}^{2}I),\end{eqnarray} \tag{ 1 }$$

where m is a mean function and K is a covariance matrix. The mean and covariance functions have hyperparameter vectors θ and ϕ, respectively. The white noise term involves measurement uncertainties σ_n. GPs can serve as predictive models for estimating values and uncertainties at times between measurements. The agreement between a GP model and observed data may be quantified by the logarithm of the marginal likelihood

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{ \mathcal L }(\theta ,\phi ) & = & -\displaystyle \frac{1}{2}{{\boldsymbol{r}}}^{\top }{\left(K+{\sigma }_{n}^{2}I\right)}^{-1}{\boldsymbol{r}}\\ & & -\ \displaystyle \frac{1}{2}\mathrm{log}| K+{\sigma }_{n}^{2}I| -\displaystyle \frac{n}{2}\mathrm{log}2\pi ,\end{array}\end{eqnarray} \tag{ 2 }$$

where the vector of residuals is

$$\begin{eqnarray}&&{\boldsymbol{r}}={\boldsymbol{y}}-m({\boldsymbol{t}}).\end{eqnarray} \tag{ 3 }$$

GPs have been used extensively in the literature for modeling correlated noise in RV data sets, and enabling detections of low-mass planets (e.g., Haywood et al. 2014; Grunblatt et al. 2015; Rajpaul et al. 2015; Cloutier et al. 2017; Benatti et al. 2020; Faria et al. 2020; Suárez Mascareño et al. 2020). In these cases, m as defined above takes the form of a Keplerian signal or sum of multiple Keplerians. The Keplerians increase $\mathrm{log}{ \mathcal L }(\theta ,\phi )$ by decreasing the residuals in r and changing the optimal hyperparameters. Some restrictions in the GP framework make it more appropriate for modeling stellar noise in RV data sets, and prevent it from absorbing planetary signals; namely, a quasiperiodic covariance function

$$\begin{eqnarray}&&{K}_{{ij}}={a}^{2}\exp \left[-\displaystyle \frac{| {t}_{i}-{t}_{j}{| }^{2}}{{\lambda }_{e}^{2}}-\displaystyle \frac{{\sin }^{2}(\pi | {t}_{i}-{t}_{j}| /{P}_{\mathrm{GP}})}{{\lambda }_{p}^{2}}\right],\end{eqnarray} \tag{ 4 }$$

which is the product of squared-exponential and sinusoidal covariance functions. The hyperparameters ϕ = {a², λ_e, λ_p, P_GP} correspond to the magnitude of covariance, a decay parameter for the overall GP evolution, a dimensionless smoothing parameter for the periodic component, and the period of oscillations, respectively. This choice is motivated by the underlying physics when P_GP equals P_rot, the rotation period of the star, and λ_e is related to the typical lifetimes of spots; however direct interpretations of hyperparameters beyond the rotation period are tenuous (Rajpaul et al. 2015). Often a jitter parameter s is added in quadrature with measurement uncertainties (Grunblatt et al. 2015). Foreman-Mackey et al. (2017) show the covariance function,

$$\begin{eqnarray}&&{K}_{{ij}}=\displaystyle \frac{B}{2+C}{e}^{-| {t}_{i}-{t}_{j}| /L}\left[\cos \displaystyle \frac{2\pi | {t}_{i}-{t}_{j}| }{{P}_{\mathrm{GP}}}+(1+C)\right],\end{eqnarray} \tag{ 5 }$$

behaves similar to Equation (4), and allows faster matrix inversion (Equation (2)). It has been used in recent RV studies (Robertson et al. 2020; Suárez Mascareño et al. 2020) and compared to other available kernels (Espinoza et al. 2020). Hyperparameters ϕ = {B, C, L, P_GP} correspond to the magnitude of covariance, weighting of the sinusoidal term, decay parameter, and period, respectively. Most studies which use GPs to model correlated RV noise adopt one of the above two covariance functions, and for our analyses we use the george implementation (Ambikasaran et al. 2015) for Equation (4) and the celerite implementation (Foreman-Mackey et al. 2017) for Equation (5). We briefly note a few differences between the kernels. The celerite GP is not mean-square differentiable (Rasmussen & Williams 2006) and is less smooth than the george GP. Also the covariance decreases faster on short timescales compared to the george GP for equal λ_e = L. RV variations due to the star are stochastic on many different timescales and are not necessarily a smooth or coherent process. However, the actual power spectrum of high-frequency variations will depend on the RV signatures of granulation and oscillations, which are not well understood. Given its smoothness, Equation (4) is a more attractive model for activity associated with faculae and spots, which themselves evolve on timescales comparable to the rotation period. We applied the GP framework with both kernels on the identical CoRoT-7 data set analyzed by Faria et al. (2016). The GP + two-planet model favored use of the george GP over the celerite GP with ${\rm{\Delta }}\mathrm{ln}{ \mathcal Z }\approx 8$ (this metric is described in the following subsection). In both cases we found that the sampler tended to converge to an alias of the 0.85 day period planet at ∼6 days, but that the shorter period planet is visible in the residuals and periodogram following subtraction of the 3.7 day planet and GP. The correct orbital period was retrieved after imposing a prior restricted to orbital periods of P < 5 days. Figure 3 shows our best fit. Moving forward we adopt the george GP for modeling RV noise.

Zoom In Zoom Out Reset image size

We couple the GP framework with the importance nested sampling algorithm MultiNest (Feroz & Hobson 2008; Feroz et al. 2009, 2019) via the PyMultiNest implementation (Buchner et al. 2014). The parameter space includes GP hyperparameters {a², λ_e, λ_p, P_GP}, jitter parameter s, systemic offset γ, and orbital parameters {K_s, ϕ₀, P, ω, e} for N planets, corresponding to semi-amplitude, phase of first epoch, orbital period, longitude of periastron, and eccentricity, respectively. We adopt ϕ₀ as a boundary condition instead of time of periastron (T_p) since there are no degeneracies in the [0, 2π] prior. MultiNest returns the log evidence of the model, $\mathrm{ln}{ \mathcal Z }$, which may be used for model comparison for N = {0, 1, 2, 3, ...} planets. Given data set d and a model ${ \mathcal M }$ parameterized by a vector θ, the Bayesian evidence is defined as

$$\begin{eqnarray}&&{ \mathcal Z }\equiv p({\boldsymbol{d}}| { \mathcal M })=\int p({\boldsymbol{d}}| {\boldsymbol{\theta }},{ \mathcal M })p({\boldsymbol{\theta }}| { \mathcal M })d{\boldsymbol{\theta }}\end{eqnarray} \tag{ 6 }$$

(Bayesian evidence in the context of RV analyses is discussed at length by Nelson et al. 2020). Assuming equal prior odds on given models ${{ \mathcal M }}_{0}$ and ${{ \mathcal M }}_{1}$ (e.g., zero planets and one planet, respectively), it is generally agreed that a difference in corresponding evidences $\mathrm{ln}{{ \mathcal Z }}_{1}-\mathrm{ln}{{ \mathcal Z }}_{0}\geqslant 5$ indicates that ${{ \mathcal M }}_{1}$ is strongly preferred over ${{ \mathcal M }}_{0}$ (Kass & Raftery 1995). Our model evidences and uncertainties correspond to the median and standard deviation of evidences from five runs of the sampler. Repeated runs are known to provide more reliable uncertainty estimates than single-run output (Nelson et al. 2020). Our choices of priors are listed in Table 3. The priors on orbital parameters represent our expectations for what could be detected in the data, given the sparse sampling and large fluctuations.

Table 3.  Priors on Parameters for the Light-curve-conditioned, One-planet RV Model

Parameter	Definition	Units	Distribution
GP
a	Amplitude	m s−1	${ \mathcal L }{ \mathcal U }(0.1,1000)$
λe	Decay Timescale	days	δ(18.74)
PGP	Periodic Timescale	days	δ(16.28)
λp	Smoothing Parameter	⋯	δ(0.76)
Global
s	Jitter	m s−1	${ \mathcal U }(0.01,5)$
γ	Systemic Offset	m s−1	${ \mathcal U }(-20,20)$
Orbital
Ks	Semi-amplitude	m s−1	${ \mathcal L }{ \mathcal U }(0.1,10)$
ϕ0	Phase of First Epoch	rad.	${ \mathcal U }(0,2\pi )$
P	Period	days	${ \mathcal L }{ \mathcal U }(0.5,20)$
ω	Longitude of Pericenter	rad.	${ \mathcal U }(0,2\pi )$
e	Eccentricity	⋯	${ \mathcal L }{ \mathcal U }(0,0.99)$

Note. We sample george GP covariance parameters (Equation (4)), global parameters including a jitter term and RV offset, and orbital parameters for a single planet. The table lists each parameter, its corresponding units, and prior distribution. For uniform (${ \mathcal U }$) and log-uniform (${ \mathcal L }{ \mathcal U }$) distributions, we specify upper and lower bounds. For this model, three of the four GP parameters are fixed (denoted by δ) to values inferred from photometry.

Download table as:  ASCIITypeset image

We perform a similar analysis as Haywood et al. (2014) by conditioning GP hyperparameters based on simultaneous photometry. We fit a george GP model to the most recent three seasons of photometry (2018–2020), including a jitter parameter and constant offset. Note, the matrix inversion becomes intractable for many data points, so we restrict this step to the timeframe overlapping with EXPRES RV observations, neglecting pre-2018 data. We use log-uniform priors spanning multiple orders of magnitude on all parameters with two exceptions: the constant offset is drawn from a uniform prior, and the periodic timescale is bounded by 12 days and 22 days. This latter constraint was chosen to prevent convergence on a harmonic or multiple of the rotation period. Afterwards, when we use a GP model to fit the RVs, the hyperparameters λ_e, λ_p, and P_GP are fixed to the maximum a posteriori (MAP) values from photometry fitting. The amplitude a is allowed to be different. Fixing model parameters decreases the sampling dimensionality and computation time. It also informs the model in case the RVs alone are insufficient to constrain the stellar rotation period. It is expected that GP fits to RV and photometry time series should have similar λ_e, λ_p, and P_GP, as demonstrated by Kosiarek & Crossfield (2020) with solar data (see their Figure 9).

In fitting a george GP model to 2018–2020 photometry, we obtain MAP hyperparameter values a = 3.40 mmag., λ_e = 18.74 days, P_GP = 16.28 days, and λ_p = 0.76. The george GP with MAP hyperparameters is plotted against the recent photometry data in Figure 4. The joint distributions between GP hyperparameters and their marginalized histograms are shown in Figure 5. All of the hyperparameters are constrained around well-defined peaks. Both λ_e and λ_p have an effect on whether the GP rapidly varies or gradually changes (Rasmussen & Williams 2006). Larger λ_e allows the GP to repeat itself more times before it loses coherence. Larger λ_p forces the repeating signal to be smoother, whereas smaller λ_p allows a more fine structure. While not strictly enforced by our priors, the evolutionary timescale converged to a value larger than the periodic timescale (λ_e > P_GP). This is a realistic constraint in regressing quasiperiodic GPs to photometry (Kosiarek & Crossfield 2020).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The zero-planet (GP-only) model is most favored at a log evidence of $\mathrm{ln}{ \mathcal Z }=-149.26\pm 0.15$, compared to the GP + one-planet model at $\mathrm{ln}{ \mathcal Z }=-149.72\pm 0.06$. We additionally restricted the Keplerian to a sinusoid (circular orbit) by fixing ω = 0 and e = 0. This model returned $\mathrm{ln}{ \mathcal Z }=-149.75\,\pm 0.18$. For the one-planet model, the MAP planet period is at 15.24 days, which is close to the adopted stellar rotation period, whereas the MAP period is at 6.68 days in the sinusoid model. The MAP eccentricity reaches 0.12 in the one-planet model, which also suggests that the Keplerian component is conforming to signals associated with stellar activity and rotation. Fit results for the three models are in Table 4.

Table 4.  Results of Our george GP Retrieval

Parameter	Units	GP-only		GP + One-planet		GP + Sinusoid
a	m s−1	${5.94}_{-1.43}^{+11.51}$	(5.06)	${5.38}_{-0.70}^{+6.91}$	(4.49)	${5.47}_{-0.80}^{+7.98}$	(5.04)
λe	days	18.74*	∗	∗	∗	∗	∗
PGP	days	16.28*	∗	∗	∗	∗	∗
λp	-	0.76*	∗	∗	∗	∗	∗
s	m s−1	${0.38}_{-0.20}^{+1.42}$	(0.27)	${0.30}_{-0.10}^{+0.81}$	(0.16)	${0.31}_{-0.11}^{+0.87}$	(0.31)
γ	m s−1	$-{0.41}_{-4.89}^{+5.27}$	(-0.41)	$-{0.42}_{-2.84}^{+2.86}$	(-0.70)	$-{0.44}_{-3.11}^{+3.22}$	(-0.60)
Ks	m s−1	-	-	${0.63}_{-0.43}^{+1.19}$	(4.02)	${0.59}_{-0.39}^{+1.18}$	(2.53)
ϕ0	rad.	-	-	${3.11}_{-1.87}^{+1.93}$	(3.59)	${3.03}_{-1.85}^{+2.03}$	(3.87)
P	days	-	-	${3.95}_{-2.64}^{+13.03}$	(15.24)	${4.20}_{-3.00}^{+13.63}$	(6.68)
ω	rad.	-	-	${3.06}_{-1.96}^{+2.01}$	(3.90)	0*	∗
e	-	-	-	${0.03}_{-0.03}^{+0.27}$	(0.12)	0*	∗
$\mathrm{ln}{ \mathcal Z }$	-	−149.26 ± 0.15		−149.72 ± 0.06		−149.75 ± 0.18
$\mathrm{ln}{{ \mathcal L }}_{\mathrm{MAP}}$	-	−139.80		−134.09		−135.81
rms	cm s−1	0.45		0.46		0.45

Note. The three models correspond to GP-only (no planet), GP + one-planet, and GP + sinusoid (one planet, restricted w and e). Columns contain the median of the marginalized distribution of each sampled parameter, and uncertainties correspond to 16th and 84th percentiles. Values in parentheses " () " denote the maximum a posteriori (MAP). Asterisks " * " denote fixed values, either from conditioning GP hyperparameters on the photometry data, or from restricting the Keplerian to a circular orbit. Missing values " - " denote that the parameter is not used in the model. The bottom rows contain the log-evidences returned by the nested sampler, the log-likelihood of the MAP vector, and the rms of the residuals. Uncertainties on model evidences correspond to the standard deviation after five separate runs of the sampler.

Download table as:  ASCIITypeset image

We attempted fitting a GP model on RVs without conditioning from photometry. For each of these models we sampled λ_e, λ_p, and P_GP from prior distributions ${ \mathcal L }{ \mathcal U }(1,100)$, ${ \mathcal L }{ \mathcal U }(0.1,100)$, and ${ \mathcal L }{ \mathcal U }(10,50)$, respectively (${ \mathcal L }{ \mathcal U }$ denotes the log-uniform distribution, with lower and upper bounds). Indeed, with appropriate amounts of data, observing cadence, and handling of model evidences, planets and GP noise parameters can be inferred from RVs alone (Faria et al. 2016). The GP-only model yields MAP values of λ_e = 30.44 days and P_GP = 15.48 days. The GP + one-planet model returns λ_e = 46.43 days and P_GP = 15.59 days. The returned planet has MAP orbital parameters K_s = 1.1 m s⁻¹ and P = 2.1 days. However, the GP-only model is favored at ${\rm{\Delta }}\mathrm{ln}{ \mathcal Z }\simeq 1$. While the inferred stellar rotation period is consistent with the photometry-derived rotation period, the spot evolution timescale is closer to twice the photometry value, probably due to the sparse sampling of the RV data.

Given that the model evidences are all within $| {\rm{\Delta }}\mathrm{ln}{ \mathcal Z }| \lt 2$ of each other, it is difficult to make decisive inferences through their comparison. The points that we would like to emphasize are that: (1) for this data set, conditioning on high-cadence photometry provides important constrains on the stellar activity model that are otherwise difficult to infer from RVs alone, including a more accurate rotation period and spot evolution timescale; (2) consistent with previous analyses of other stars, the spot decay timescale tends to be longer than the stellar rotation period but within the same order of magnitude; and (3) the Bayesian evidence favors an activity model without planets. The GP's flexibility here is crucial since the activity signal quickly loses coherence, and spots follow various rotation periods depending on latitude. Both of these aspects are addressed by choosing a proper λ_e and a single characteristic P_GP, respectively. In the following section we discuss the temporal sampling of our data, and how observing scheduling can improve the sensitivity of our analysis to short-period signals. In particular, additional high-cadence RVs are necessary to rule out some planets at the K_s ≈ 3 m s⁻¹ level.

To quantify the impact of cadence in our analysis, we perform a simple injection/recovery test as follows. First, we generate a synthetic RV time series by drawing a george GP specified by Equation (4), characteristic of our HD 101501 observations. Hyperparameters are set to the MAP values in Section 4.1. The amplitude is set to 5.06 m s⁻¹, which is the the MAP amplitude of the GP-only fit in Section 4.2. Next, we add a Keplerian component with K_s and P sampled from a grid. Other Keplerian parameters are set to zero. We sample the continuous RV curve at 33 epochs (described below) and add random noise drawn from a zero-mean normal distribution with a standard deviation of 40 cm s⁻¹, which is also the associated uncertainty on each data point. Finally, we fit a GP + sinusoidal model as discussed in Section 3 with three GP hyperparameters fixed, as if they had been predetermined by photometry measurements. We then compare the recovered K_s and P to the actual injected signal. Each full run involves a pair (K_s, P) and new GP draw. The sampler is run only once since we are not interested in precise uncertainties on the model evidence. The procedure is illustrated in Figure 7.

Zoom In Zoom Out Reset image size

The test described above is repeated for several toy observing cadences. We define an N-day cadence as observing for N consecutive nights (separated by one sidereal day), followed by a long gap in time. The gap is drawn from a uniform random distribution between 50 and 80 days. This pattern repeats until the number of exposures totals 33, which was the number of nights HD 101501 was observed by EXPRES at the start of simulations (a couple of additional, recent nights were included in the RV analysis). Each timestamp is then perturbed by a uniformly random variable between −4 and +4 hr to simulate variations in observing scheduling. The 50–80 day gap is about 3–5 × λ_e, chosen to destroy coherence between consecutive bursts of exposures. A variable gap helps avoid unwanted sampling artifacts. We repeat the above test for 2, 5, 10, 20, and 33 day cadences. For example, the 5 day cadence involves timestamps for five consecutive nights, followed by a long gap in time. The pattern repeats six times, followed by three exposures such that the number of data points equals 33. The cadences roughly sample 1/8, 1/4, 1/2, 1, and 2 stellar rotations. As a benchmark we also sample at the identical timestamps of EXPRES data. These tests do not encompass sophisticated modeling of stellar noise or the wide variety of observing strategies and possible Keplerian signals (e.g., eccentric orbits, multiple planets). However, the assumptions made are reasonable for a goal of understanding how consecutive nights of observation relates to stellar activity inferences for stars like HD 101501.

The injected Keplerian is successfully recovered if the following criteria are met: first, the recovered semi-amplitude must be greater than 3σ_K, where σ_K is the standard deviation of K_s draws made by the nested sampler; second, the fractional error on the period must be within 10%. We also check against 1.0027 day⁻¹ aliases of the recovered period (Dawson & Fabrycky 2010; 1 day = 1.0027 sidereal day). Our injection + recovery analysis indicates that high observing cadence, consisting of several nights of consecutive observation, is necessary for the GP framework to identify certain classes of exoplanets orbiting active stars. Figure 8 depicts the recovery of planets for a given cadence (green denotes success, whereas orange denotes convergence on an alias).

Zoom In Zoom Out Reset image size

The 10, 20, and 33 day cadences have better completion across the parameter grid, especially for short periods. The observed cadence also has comparatively good completion, mostly due to there being 76 total timestamps (multiple per night) instead of 33. As expected, long-period/low-amplitude signals are difficult to recover for all cadences. We emphasize that this is not a full Monte Carlo analysis, and that unsuccessful retrievals at K_s > 5 m s⁻¹ and P < 5 days are likely a result of randomness in the GP draw and unfavorable sampling of the orbit. These high-amplitude Keplerians might be retrieved with a more careful analysis (e.g., selection of priors, more thorough posterior sampling, etc.). Planets with periods ≲1.5 days are difficult to distinguish from their 1 day aliases, but are nevertheless recovered by the GP framework. Most modern RV analyses are prone to occasionally identifying aliases over true signals (Dumusque et al. 2017).

Importantly, the results show generally increasing completion with higher cadence. The 2 day cadence fails for nearly all cases except amplitudes ≳5 m s⁻¹, when the planet signal starts to rival the activity signal. The 5 day cadence exhibits improvements, and the 33 day cadence shows the greatest completion of K_s and P combinations. There are diminishing returns going from the 10 day to 20 day cadence, and the 20 day to 33 day cadence. For short-period/low-amplitude planets (K_s ≲ 3 m s⁻¹, P ≲ 5 days), if sampling covers multiple orbital periods within a couple spot lifetimes, then the high-frequency Keplerian can be distinguished from the low-frequency GP/activity component. In this case the actual stellar rotation period may be irrelevant and alternative kernels such as the squared-exponential may be suitable for modeling the activity (Rajpaul et al. 2015). The most challenging planets (K_s ≲ 2 m s⁻¹) are retrieved best by the highest cadences (20 days and 33 days). In these cases, the GP is capable of learning a repeated structure in the activity signal.

Given the above results, we would like to gauge whether the GP framework and EXPRES RVs of HD 101501 are actually sensitive to the full range of possible orbits specified by the priors. We perform a Monte Carlo analysis in which we repeat the injection/recovery procedure, sampled at the actual timestamps of exposures, 10 times (i.e., we generate the bottom-right panel of Figure 8 an additional 10 times, each involving new GP draws). We group correct orbital periods and aliases together as successful recoveries. The results are shown in Figure 9. As expected the longest periods and smallest semi-amplitudes are most challenging to detect, as are periods at nearly 1 day. However, even for modest semi-amplitudes up to ∼3 m s⁻¹, planets with periods less than a couple days are difficult to detect.

Zoom In Zoom Out Reset image size

We explore two additional extensions to the above simulations. First, to what degree is detection reliant on the GP component? Second, how much easier is it to recover a planet with a known ephemeris (i.e., the planet transits)? These investigations involved modifying the underlying model to, in the first case, exclude a GP component and just fit a Keplerian, jitter term, and offset term. For the second case we use the full model, but all Keplerian parameters are fixed to their true values except for semi-amplitude. The results of these tests are shown in the Appendix. Without a GP, the sampler has much greater difficulty identifying planets in the simulated RV data sets. Given the simulated cadences, a GP is necessary to detect most planets with K_s < 5 m s⁻¹. While the sampler does identify Keplerians at P ∼ 15 days, the retrievals might be falsely converging on the stellar activity signal given the similar stellar rotation period. In the future, it would be useful to explore the efficacy of GPs on RV data sets when the planet period is near the stellar rotation period. On the other hand, when the ephemeris is known and fixed a priori, we see a dramatic improvement in recovery when a GP is used. Robust detections are difficult for only the lowest semi-amplitudes and longest periods, limited by factors such as phase coverage, time baseline, and measurement uncertainty.

Several approaches directly incorporate indicators into the GP framework. Indicators may be modeled jointly with RVs as linear combinations of a single, latent GP and its derivative (Rajpaul et al. 2015; Jones et al. 2017; Gilbertson et al. 2020). Another method involves regression on RVs and an indicator time series where both GPs share certain hyperparameters (Suárez Mascareño et al. 2020) or, one can train a GP on an indicator time series and use the results as an initial guess for a subsequent RV analysis (Dumusque et al. 2017).

In Figure 10 we show the RV time series along with CCF FWHM, CCF BIS, and Hα equivalent width, as well as their generalized Lomb–Scargle periodograms (Scargle 1982; Zechmeister & Kürster 2009). The RV and BIS time series have power at the stellar rotation period, with BIS more pronounced. BIS shows a double-peaked structure, split between roughly the GP-derived rotation period and another peak near ∼17.2 days. This pattern is likely due to differential rotation, where long-lived spots at different latitudes exhibit different rotation periods. We attempted joint modeling of RVs and BIS with separate GPs that share λ_e, λ_p, and P_GP. Their likelihood functions were summed together. We sampled the BIS GP amplitude, mean, and jitter from similar distributions as used for the RVs. The zero-planet model is favored over the one-planet model, with $\mathrm{ln}{{ \mathcal Z }}_{0}=-363.98\pm 0.07$ and $\mathrm{ln}{{ \mathcal Z }}_{1}\,=-364.92\pm 0.13$, respectively. The corresponding MAP values of P_GP are 15.87 days and 15.65 days, respectively, similar to the rotation period dervied from RVs alone. The sampler converges to λ_e = 34.01 days and λ_e = 27.47 days for the two models, respectively. The APT light curve changes significantly between rotations, so the actual spot evolution is likely much faster than these estimates.

Zoom In Zoom Out Reset image size

We make a few comments on the periodograms, which may simply be spurious features of the data set. The RV periodogram shows highest power at ∼28 days; however, none of the GP models with planets converged to this orbital period in additional trials where we extended the orbital period prior to 40 days. Rather, they still tended toward the other peak at 15.3 days. Since the indicators do not have significant power at 15.3 days, a Keplerian might be responsible. However, we reiterate that the Bayesian evidence disfavors a planet in all of our fits (light-curve GP conditioning, BIS joint fitting, and freely fitting RVs), and 15.3 days is close to the stellar rotation period.

Our RV analysis warranted use of the most recent seasons of photometry. However, the remarkable 28 season baseline offers an opportunity to obtain more precise constraints on stellar rotation and activity. We again turn to GP regression, which has seen frequent applications toward inferring stellar properties and searching for transits in photometry (Vanderburg et al. 2015; Angus et al. 2018; Barros et al. 2020). We use a celerite GP with a covariance function given by Equation (5) and sample hyperparameters {B, C, L, P_GP}. We find a periodic timescale of ${P}_{\mathrm{GP}}={16.45}_{-0.51}^{+0.60}$ days, and an evolutionary timescale of $L={15.82}_{-3.07}^{+9.04}$ days (uncertainties denote the 16th and 84th percentiles around the median). The MAP P_GP is 16.42 days, which we take as our best estimate of the rotation period P_rot. Technically, this is probably not the equatorial rotation period. However, it is the period that best describes the data, and it might be influenced by the typical latitudes of spots. The GP with MAP hyperparameters is shown in the bottom panel of Figure 1. The evolutionary timescale L is consistent within 1σ of λ_e from our earlier analysis where we fit a george GP to the recent photometry. The periodic timescale P_GP also matches between the two fits. Formally, P_GP has a different definition in Equation (4) than it does in Equation (5), but it has similar meaning and influence on GP behavior. Their similarity gives some assurance that the stellar activity signal in the most recent three seasons shares similar characteristics with those of the whole baseline. The closeness of L and P_rot indicates that the spot distribution changes significantly between consecutive rotations, which also makes the light curve less coherent.

Some of our photometry seasons (e.g., 2001) have very coherent oscillations, while others (e.g., 1998) appear more random. We investigate inter-seasonal variations by fitting a celerite GP to each season individually. The GP provides, in theory, a more reliable estimate of the rotation period than the maximum power of the periodogram. The periodogram is based on a single sinusoid model and most clearly identifies a signal when the phase, amplitude, and period are constant. The quasiperiodic GP can accommodate signals that exhibit small departures from an overall phase, amplitude, and period, for example due to the appearance and disappearance of spots. However, the robustness of the GP is tied to the quantity and cadence of the data (lacking in 2017, for example). Also, if multiple modes are present within a given season, the GP learns a value that maximizes the likelihood, which might not actually be representative of any single mode. In 10 of 28 seasons the best-fit GP has P_GP < 16 days, reaching as low as 12–14 days. It is below 17 days in all cases. The variability in periodic timescales is most likely a sign of differential rotation. In a previous study, Mittag et al. (2017) analyze periodograms of Ca ii H&K and Ca ii IRT line strengths and find power at multiple periods, which they also attribute to differential rotation. We perform a bootstrap by sampling the 28 values of P_GP with replacement 10,000 times, recording the rotational shear α = ΔP/P at each iteration. The difference between maximum and minimum periods, ΔP, assumes that the maximum corresponds to rotation at the poles and the minimum at the equator. We find $\alpha ={0.45}_{-0.19}^{+0.03}$, which is similar to the solar rotational shear of ∼0.4 (Snodgrass & Ulrich 1990).

In order to better understand the evolution of the starspots on the surface of HD 101501, we applied a light-curve inversion algorithm to reconstruct the stellar surfaces. The algorithm light-curve inversion (LI; Harmon & Crews 2000) uses a modified Tikhonov regularization and stellar parameters to converge on a solution and makes no a priori assumptions about starspot shape, number, or size.

For the APT light curve, we divided the light curve into portions lasting approximately one rotation. When necessary to provide more phase coverage, the portion of the light curve used for inversion lasted more than one rotation (this is noted in the lower right corners in Figures 11 and 12). In very few cases, data points were used for overlapping portions of the light curve. For the TESS light curve, there is a gap over 4 days long in the middle of the sector of observation (see Figure 2). We divided the light curve during this time into two rotations and binned the data into bins 0.01 in phase. Both rotations observed with TESS have incomplete phase coverage.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We assigned T_eff = 5500 K and a spot temperature of T_spot = 4500 K, based on Figure 7 of Berdyugina (2005). Combining those temperatures and the appropriate response function for the filters yields a spot-to-photosphere ratio of 0.3452 for the V band and 0.4530 for the TESS observations. In addition to this value, LI also uses limb-darkening coefficients as input. For the V-band APT inversions, quadratic limb-darkening coefficients are used (0.5955 and 0.1488; Claret et al. 2013). For the TESS inversions, quadratic limb-darkening coefficients are also used (0.3876 and 0.2036; Claret 2018). Each inversion used a unique rms error in order to balance the inversion algorithm between overfitting and oversmoothing the data (the average rms for the APT light curves was 0.0013 and both TESS light curves used an rms of 0.0003). A stellar inclination of i = 52° was found from the v sin i given by Brewer et al. (2016), the stellar diameter (Bonneau et al. 2006), and the parallax (Gaia Collaboration et al. 2018).

The reconstructed surfaces can be found in Figures 11 and 12. The reconstructed dark regions on the stellar surface are not necessarily individual starspots, but may be unresolved groups of spots. We cannot differentiate between these possibilities, and we refer to the dark regions as starspots. The starspots of HD 101501 behave similarly to sunspot groups in that they typically evolve on the same timescale as a rotation. On some occasions, the same spot appears to be present for more than one rotation. The latitudes at which the starspots appear are only weakly constrained by the light curve and the limb-darkening coefficients used. Because of this and the short starspot lifetimes, we do not further investigate differential rotation for this star. The variations in the photometric data and the resultant surface reconstructions indicate that at all points of observation, there is evidence of starspots on the surface. Furthermore, the surfaces that result from the inversion of TESS light curves show similar structures, but are not a direct comparison to the V-band reconstructions because the bandpasses are very different. The TESS filter covers a much larger wavelength range, and the contrast between cool starspots and the photosphere is more dramatic at shorter wavelengths.

The Sun-like star HD 101501 presents an interesting case study for understanding stellar activity signatures in RVs and photometry. We present new high-precision RVs of HD 101501 along with simultaneous ground-based photometry with a baseline of 28 years. Several weeks of photometry from TESS are also analyzed. We summarize our findings as follows.

• 1.  
  A GP framework is used to model both the photometry as well as the correlated noise in the RV time series. HD 101501 represents a case study for the GP framework that may be applied to RVs of other EXPRES targets. The RVs, which exhibit variations at a level of ∼10 m s⁻¹, are best explained by an activity-only model. The Bayesian evidence disfavors the presence of a Keplerian signal in our data.
• 2.  
  Through a simple injection and recovery analysis, we explore the space of orbital parameters to which the GP framework is sensitive. We test different cadences in order to understand how observing strategy impacts our ability to detect planets. The lowest cadences, in which exposures are largely spaced in time, contribute very little to the GP retrieval. Higher cadences, in which the star is observed for many consecutive nights, assist the GP framework in separating short-period orbits from the stellar activity signal. These results refine our observing plans and offer important guidance for RV observations of active stars.
• 3.  
  GPs place tight constraints on the stellar rotation period associated with spots, at P_rot ∼ 16.4 days. We use GPs to analyze periodicity in individual photometric observing seasons. The variability from season to season suggests a rotational shear of ∼0.45 and an equatorial rotation period of ∼13 days.
• 4.  
  Reconstructed stellar surfaces show the persistent presence of starspots on the surface of HD 101501 at all times. While starspots are always present, they are observed to change significantly between rotations making it impossible to trace their evolution over many rotations in these data.

Detecting exoplanets around active stars remains a significant challenge. Correlated noise models show great promise for mitigating activity in RVs, especially when combined with simultaneous photometry. GPs have had great success in modeling stellar activity in HARPS, HARPS-N, CARMENES, and ESPRESSO RVs, and this analysis of HD 101501 represents the first application of GPs to EXPRES RVs. The high precision of current spectrographs, optimized observing strategy, and new RV extraction techniques will push exoplanet detection limits in the near future.