EXPRES. III. Revealing the Stellar Activity Radial Velocity Signature of Eridani with Photometry and Interferometry

Rachael M. Roettenbacher, Samuel H. C. Cabot, Debra A. Fischer, John D.

Monnier, Gregory W. Henry, Robert O. Harmon, Heidi Korhonen, John M.

Brewer, Joe Llama, Ryan R. Petersburg, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Rachael M. Roettenbacher, Samuel H. C. Cabot, Debra A. Fischer, John D. Monnier, Gregory W. Henry, et al. EXPRES. III. Revealing the Stellar Activity Radial Velocity Signature of Eridani with Photometry and Interferometry. The Astronomical Journal, 2022, 163, 10.3847/1538-3881/ac3235. insu-03705395

HAL Id: insu-03705395

https://insu.hal.science/insu-03705395

Submitted on 27 Jun 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

EXPRES. III. Revealing the Stellar Activity Radial Velocity Signature of ϵ Eridani with Photometry and Interferometry

Rachael M. Roettenbacher<fnref n="1" />,<fnref n="2" /> ID, Samuel H. C. Cabot<fnref n="2" /> ID, Debra A. Fischer<fnref n="2" /> ID, John D. Monnier<fnref n="3" /> ID, Gregory W. Henry<fnref n="4" /> ID,
Robert O. Harmon<fnref n="5" />, Heidi Korhonen<fnref n="6" /> ID, John M. Brewer<fnref n="7" /> ID, Joe Llama<fnref n="8" /> ID, Ryan R. Petersburg<fnref n="9" /> ID, Lily L. Zhao<fnref n="2" /> ID,
Stefan Kraus<fnref n="10" /> ID, Jean-Baptiste Le Bouquin<fnref n="11" /> ID, Narsireddy Anugu<fnref n="12" /> ID, Claire L. Davies<fnref n="10" /> ID, Tyler Gardner<fnref n="3" /> ID,
Cyprien Lanthermann<fnref n="13" />,<fnref n="14" /> ID, Gail Schaefer<fnref n="14" /> ID, Benjamin Setterholm<fnref n="3" /> ID, Catherine A. Clark<fnref n="15" /> ID, Svetlana G. Jorstad<fnref n="16" /> ID,
Kyler Kuehn<fnref n="8" />,<fnref n="17" /> ID, and Stephen Levine<fnref n="8" /> ID

Received 2021 July 31; revised 2021 September 30; accepted 2021 October 20; published 2021 December 16

Abstract

The distortions of absorption line profiles caused by photospheric brightness variations on the surfaces of cool,
main-sequence stars can mimic or overwhelm radial velocity (RV) shifts due to the presence of exoplanets. The
latest generation of precision RV spectrographs aims to detect velocity amplitudes ~10 cm s⁻<fnref n="1" />, but requires
mitigation of stellar signals. Statistical techniques are being developed to differentiate between Keplerian and
activity-related velocity perturbations. Two important challenges, however, are the interpretability of the stellar
activity component as RV models become more sophisticated, and ensuring the lowest-amplitude Keplerian
signatures are not inadvertently accounted for in flexible models of stellar activity. For the K2V exoplanet host ϵ
Eridani, we separately used ground-based photometry to constrain Gaussian processes for modeling RVs and
TESS photometry with a light-curve inversion algorithm to reconstruct the stellar surface. From the reconstructions
of TESS photometry, we produced an activity model that reduced the rms scatter in RVs obtained with EXPRESS
from 4.72 to 1.98 m s⁻<fnref n="1" />. We present a pilot study using the CHARA Array and MIRC-X beam combiner to
directly image the starspots seen in the TESS photometry. With the limited phase coverage, our spot detections are
marginal with current data but a future dedicated observing campaign should allow for imaging, as well as allow
the stellar inclination and orientation with respect to the debris disk to be definitively determined. This work shows
that stellar surface maps obtained with high-cadence, time-series photometric and interferometric data can provide
the constraints needed to accurately reduce RV scatter.

*Unified Astronomy Thesaurus concepts: Planet hosting stars (1242); Radial velocity (1332); Starspots (1572)*

*Supporting material: machine-readable tables*

**1. Introduction**

Radial velocity (RV) surveys for exoplanets have not yet
been able to detect planets with similar masses and radii to
those of the Earth in Earth-like orbits around Sun-like stars.
However, the latest generation of spectrographs are designed to
reduce instrumental error sources with the goal of isolating the
stellar signals that obstruct the detection of low-amplitude
velocity signals. With extreme-precision RV (EPRV) surveys

of solar analogs, features such as starspots on the stellar surface
produce temporal variations in the shapes of line profiles
that add time-correlated variations to the center-of-mass RV
measurements and must be properly accounted for.

A number of phenomena contribute to the absorption line
profile signatures caused by stellar activity that result in RV shifts.
The convective envelope of cool, Sun-like stars is composed of
cells, or granules, in which hot stellar material rises and then falls
as it cools. Localized, strong magnetic fields suppress convection
and manifest as bright regions, including faculae, plages, and
networks, and dark starspots on the stellar surface. As these
features rotate in and out of view, they create rotationally-
modulated absorption line signatures that lead to periodic RV
signatures. Faculae, which contribute RV amplitudes <1 m s⁻<fnref n="1" />,

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms
of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further
distribution of this work must maintain attribution to the author(s) and the title
of the work, journal citation and DOI.

<footnote n="1">¹ Yale Center for Astronomy and Astrophysics, Yale University, 46 Hillhouse Avenue, New Haven, CT 06511, USA; rachael.roettenbacher@yale.edu</footnote>

<footnote n="2">² Department of Astronomy, Yale University, 52 Hillhouse Avenue, New Haven, CT 06511, USA</footnote>

<footnote n="3">³ Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="4">⁴ Tennessee State University, Center of Excellence in Information Systems, Nashville, TN 37203, USA</footnote>

<footnote n="5">⁵ Department of Physics and Astronomy, Ohio Wesleyan University, Delaware, OH 43015, USA</footnote>

<footnote n="6">⁶ European Southern Observatory, Alonso de Córdova 3107, Vitacura, Casilla 19001, Santiago, Chile</footnote>

<footnote n="7">⁷ Department of Physics and Astronomy, San Francisco State University, 1600 Holloway Avenue, San Francisco, CA 94132, USA</footnote>

<footnote n="8">⁸ Lowell Observatory, 1400 West Mars Hill Road, Flagstaff, AZ 86001, USA</footnote>

<footnote n="9">⁹ Department of Physics, Yale University, 217 Prospect Street, New Haven, CT 06511, USA</footnote>

<footnote n="10">¹⁰ Astrophysics Group, Department of Physics & Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, UK</footnote>

<footnote n="11">¹¹ Institut de Planetologie et d'Astrophysique de Grenoble, Grenoble F-38058, France</footnote>

<footnote n="12">¹² Steward Observatory, Department of Astronomy, University of Arizona, 933 North Cherry Avenue, Tucson, AZ, 85721, USA</footnote>

<footnote n="13">¹³ Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, B-3001, Leuven, Belgium</footnote>

<footnote n="14">¹⁴ The CHARA Array of Georgia State University, Mount Wilson Observatory, Mount Wilson, CA 91203, USA</footnote>

<footnote n="15">¹⁵ Department of Astronomy and Planetary Science, Northern Arizona University, Flagstaff, AZ 86011, USA</footnote>

<footnote n="16">¹⁶ Institute for Astrophysical Research, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA</footnote>

<footnote n="17">¹⁷ Australian Astronomical Optics, Faculty of Science and Engineering, Macquarie University, Macquarie Park, NSW 2113, Australia</footnote>
are bright with respect to the photosphere and are most apparent when near the limb of the stellar disk (Haywood et al. 2016). Plages and networks are also bright with respect to the photosphere, but have a more significant RV amplitude contribution on the order of a few meters per second, as seen on the Sun (Milbourne et al. 2019). Starspots are dark features against the photosphere and can cause line profile distortions that contribute a wide range of RV amplitudes—from less than 1 m s⁻¹ to several kilometers per second for large starspots (Roettenbacher et al. 2015; Haywood et al. 2016). Even for chromospherically quiet stars, these signatures can overwhelm the ~10 cm s⁻¹ RV signature of an Earth analog in the habitable zone.

Ongoing efforts to characterize and isolate stellar activity signals include modeling stellar activity with flexible correlated noise models (e.g., Gaussian processes (GPs); Haywood et al. 2014; Gilbertson et al. 2020); statistically identifying stellar activity (e.g., PCA or F-statistic; Davis et al. 2017; Holzer et al. 2021); estimating RV variations from a spot model applied to photometry (Aigrain et al. 2012; Dumusque et al. 2014); extracting signatures of stellar activity from the cross-correlation function used to measure RVs (Collier Cameron et al. 2021); and using Doppler imaging of young stars to filter stellar activity from line profiles (Heitzmann et al. 2021). Several additional or related methods are compared by Dumusque et al. (2017) and L. L. Zhao et al. (2021, submitted). These methods are based on photometric and/or spectroscopic data and aim to account for RV scatter due to photospheric activity.

There are a number of ways to reconstruct the stellar surfaces to resolve some of the surface structure contributing to RV signatures. Photometric light-curve inversion uses one or more light curves of a star showing rotational variability to reconstruct the stellar surface (Harmon & Crews 2000; Savanov & Strassmeier 2008; Luo et al. 2019). This method makes no assumptions on the spot size, shape, or number, and has been shown to be reliable as compared to simulations and interferometric images (Harmon & Crews 2000; Roettenbacher et al. 2017). While the method gives starspot longitudes, degeneracies only allow for the determination of relative starspot latitudes when multiple features are present, and degeneracies will remain between the size of a starspot and its absolute latitude (e.g., Harmon & Crews 2000). Doppler imaging is a method to reconstruct relatively large surface spots using high-resolution spectra that can better determine starspot latitude (e.g., Vogt et al. 1987; Rice et al. 1989), but a degeneracy between the hemispheres remains. Doppler imaging requires a rotational velocity of ~10 km s⁻¹, restricting the stars to which method can be applied (example applications of the method include those in Korhonen et al. 2021 and Şenavcı et al. 2021). Interferometric aperture synthesis imaging is a third method that can be used to reconstruct the stellar surface. Because interferometric imaging allows for stars and their spots to be imaged as they appear on the sky, degeneracies in the starspot latitudes of other techniques are resolved. With observations spanning a stellar rotation period, the stellar inclination and the star's orientation on the sky —the position angle of the axis of rotation—can be measured. Resolving the surfaces of stars is currently only possible for bright stars that have angular diameters ~2 mas or more and relatively large starspots (e.g., Roettenbacher et al. 2016, 2017; Martinez et al. 2021; Parks et al. 2021).

Knowing the brightness inhomogeneities, such as starspots, present at the time of RV observations can provide a way to

separate their impact on line profiles and thus on RV signatures from those of the planets. This has provided motivation for solar telescopes that measure RVs for disk-integrated spectra of the Sun (Collier Cameron et al. 2021; Dumusque et al. 2021). Here, we use a similar approach to analyze the star ϵ Eridani (HD 22049, TIC 118572803). The closest K2 dwarf to the Sun, ϵ Eri is at a distance of 3.220 ± 0.004 pc (Gaia Collaboration et al. 2021). It is a bright (V = 3.73, H = 1.75; Ducati 2002), main-sequence star with a radius of 0.74 ± 0.01 R☉ derived from interferometry (di Folco et al. 2007; Baines & Armstrong 2012). It is known to be active with variable starspots and detected global activity cycles (Metcalfe et al. 2013). The star has a rotation period of approximately 11 days, detected with a variety of techniques, including MOST photometry and RV variations (Giguere et al. 2016), modulation of Ca II H&K measurements (Hempelmann et al. 2016), and photometry (Lanza et al. 2014). ϵ Eri is also a known exoplanet host star (Hatzes et al. 2000; Mawet et al. 2019) with an RV-detected 0.8 MJup planet in a 7.4 yr orbit.

Here, we utilize stellar surface images to isolate the associated RV signatures in EPRV data. In Section 2, we present the spectroscopic, interferometric, and photometric observations used in this work. Spectroscopic observations include an extensive, multi-decade baseline archival data set, as well as new, high-precision RV measurements. In Section 3, we detail our GP analysis, infer attributes of the stellar activity, and confirm orbital parameters of the known planetary companion. In Section 4, we reconstruct the stellar surface with a light-curve inversion algorithm and model interferometric data obtained on two nights. In Section 5, the surface reconstructed from a light curve is combined with a disk model that simulates the stellar spectrum. We discuss the simulated stellar spectra and show how they successfully account for a significant portion of scatter in contemporaneous RV measurements. The model’s efficacy, additional considerations, and future work are discussed in Section 6.

## 2. Observations

ϵ Eri has been extensively observed owing to its brightness and proximity to the Sun. In this section, we describe the observations that were used for this work, including brief descriptions of the archival data.

### 2.1. *EXPRES* Spectroscopy

High-resolution spectra of ϵ Eri were acquired with the EXtreme PREcision Spectrograph (EXPRES; Jurgenson et al. 2016) commissioned at the 4.3 m Lowell Discovery Telescope (LDT; Levine et al. 2012). EXPRES is an optical spectrograph optimized for wavelengths 380–780 nm, reaching a typical resolving power of R ~ 137,500. Additional specifications can be found in works detailing the RV pipeline, instrument performance verification, and first science results (Blackman et al. 2020; Brewer et al. 2020; Petersburg et al. 2020). EXPRES attains ~30 cm s⁻¹ RV precision for spectra of slowly-rotating, main-sequence FGK-stars when the signal-to-noise ratio reaches 250 at 550 nm. While the lowest rms RV scatter is observed around chromospherically inactive stars, we observed ϵ Eri as an interesting case study for characterizing and mitigating RV jitter in a moderately active star.

We obtained 164 RVs of ϵ Eri on 39 distinct nights between 2019 August 15 and 2020 November 13 (see Table 1). Besides
<table><caption>**Table 1**
EXPRES RVs</caption><thead><tr><th>Reduced Julian Date<br>(RJD = JD - 2,400,000.0)</th><th>v<br>(m s<sup>-1</sup>)</th><th>σ<sub>v</sub><br>(m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>58,710.983924</td><td>-9.762</td><td>0.482</td></tr><tr><td>58,710.985348</td><td>-10.127</td><td>0.466</td></tr><tr><td>58,710.986526</td><td>-9.149</td><td>0.692</td></tr><tr><td>58,716.996130</td><td>-14.496</td><td>0.387</td></tr><tr><td>58,716.997549</td><td>-10.840</td><td>0.380</td></tr><tr><td>...</td><td>...</td><td>...</td></tr></tbody></table>

(This table is available in its entirety in machine-readable form.)

one seasonal gap, $\epsilon$ Eri was observed every 3–10 nights. Stellar activity dominates the RVs, yielding an rms scatter of 6.6 m s⁻¹ with an average single-measurement uncertainty of 35 cm s⁻¹.

## 2.2. Archival Spectroscopy

Mawet et al. (2019) conducted an extensive analysis of archival $\epsilon$ Eri RVs spanning 30 yr, coupled with direct imaging of the system. They placed tight constraints on the orbit of the 7 yr planet first identified by Hatzes et al. (2000). The study made use of over 450 RVs obtained with the following instruments (their respective observatories are listed immediately after them): High Resolution Echelle Spectrometer/Keck, Levy Automated Planet Finder/Lick, Hamilton/Lick, Coudé Echelle/La Silla, and High Accuracy Radial Velocity Planet Searcher/La Silla. We refer the reader to Mawet et al. (2019) and Zechmeister et al. (2013) for details on these RVs. Our analysis of the archival data additionally includes CHIRON RVs acquired in 2014 (Giguere et al. 2016).

## 2.3. TESS Photometry

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) observed $\epsilon$ Eri during Sector 31 (2020 October 21 through 2020 November 19). The 2 minute cadence simple aperture photometry (SAP) light curves were obtained through the Barbara A. Mikulski Archive for Space Telescopes (MAST).

To account for scattered light and other systematic issues, we typically remove co-trending basis vectors (CBVs) from SAP light curves (following the procedure used in Roettenbacher & Vida 2018; Cabot et al. 2021). However, we did not remove any CBVs from the Sector 31 light curve, as the CBVs provided through MAST for Sector 31 (with a file creation date of 2020 December 12) add trends to the data or increase the noise. The shape of the SAP light curve is consistent with the nearly contemporaneous, ground-based light curve described in Section 2.4, suggesting stellar variability is likely the dominant signature. We removed data flagged as bad by the TESS pipeline (a nonzero quality flag). There were no usable data provided for the SAP light curve between JD 2,459,155.8948312 and JD 2,459,158.8670810 due to scattered light from the Moon and the transition between orbits of the satellite. $\epsilon$ Eri was observed with Camera 1, which was strongly impacted by scattered light at the time of these removed data points according to the Sector 31 Data Release Notes (DR47).<fnref n="18" /> We removed a total of 2376 observations from the original 18,314 observations.

<table><caption>**Table 2**
Strömgren $(b + y)/2$ Differential Photometric Data of $\epsilon$ Eri with the APT T8</caption><thead><tr><th>Modified Julian Date<br>(MJD = JD - 2,400,000.5)</th><th>$(b + y)/2$ Differential<br>Magnitude</th><th>Trend Removed</th></tr></thead><tbody><tr><td>56,576.813</td><td>-1.98170</td><td>-1.97666</td></tr><tr><td>56,576.901</td><td>-1.98095</td><td>-1.97666</td></tr><tr><td>56,577.815</td><td>-1.98250</td><td>-1.97661</td></tr><tr><td>56,584.791</td><td>-1.97605</td><td>-1.97630</td></tr><tr><td>56,584.875</td><td>-1.97405</td><td>-1.97629</td></tr><tr><td>...</td><td>...</td><td>...</td></tr></tbody></table>

(This table is available in its entirety in machine-readable form.)

<table><caption>**Table 3**
Strömgren $(b + y)/2$ Differential Photometric Data of $\epsilon$ Eri with the APT T4</caption><thead><tr><th>Modified Julian Date<br>(MJD = JD - 2,400,000.5)</th><th>$(b + y)/2$ Differential<br>Magnitude</th><th>Trend Removed</th></tr></thead><tbody><tr><td>59,152.758</td><td>-1.96490</td><td>-1.96735</td></tr><tr><td>59,152.820</td><td>-1.96360</td><td>-1.96735</td></tr><tr><td>59,152.844</td><td>-1.96520</td><td>-1.96735</td></tr><tr><td>59,152.878</td><td>-1.96450</td><td>-1.96735</td></tr><tr><td>59,152.941</td><td>-1.96405</td><td>-1.96735</td></tr><tr><td>...</td><td>...</td><td>...</td></tr></tbody></table>

(This table is available in its entirety in machine-readable form.)

## 2.4. Automated Photoelectric Telescope Photometry

$\epsilon$ Eri was observed with two Automatic Photoelectric Telescopes (APTs) at Fairborn Observatory (Henry 1999). Six hundred observations were obtained with the T8 (0.8 m) telescope from 2013 October 11 through 2020 February 20 and 64 observations on the T4 (0.75 m) telescope from 2020 October 30 through 2020 November 30. Differential photometric observations were obtained in the Strömgren *b* and *y* bandpasses, but are presented as a combined $(b + y)/2$ bandpass. The stars HD 22243 and HD 23281 were used as comparison stars to ensure the variable signature was that of $\epsilon$ Eri.

We removed the long-term trends in both of these data sets. To do so, we smoothed the data with a Gaussian kernel with an FWHM of 50 days and removed that signature to leave only the signature of the starspots that were rotating in and out of view. These data are included in Tables 2 and 3.

## 2.5. MIRC-X Interferometry

Two sets of interferometric observations were obtained on 2020 November 2 (JD 2,459,155.77) and 2020 November 5 (JD 2,459,158.87) at the Center for High Angular Resolution Astronomy (CHARA) Array (ten Brummelaar et al. 2005) using the Michigan InfraRed Combiner-eXeter (MIRC-X; Anugu et al. 2020). The CHARA Array consists of six 1 m telescopes in a Y-shaped array with nonredundant baselines (B) extending from 34 to 330 m (ten Brummelaar et al. 2005). The snapshot observations were obtained in the H band (spanning 1.5–1.8 µm; $\lambda/2B \approx 0.5$ mas) in the grism mode (R ~ 190). The observations of $\epsilon$ Eri were followed by those of the calibration star, HD 26912 (uniform disk angular diameter measured in H band with no limb darkening, $\theta_{UD,H} = 0.285 \pm 0.026$ mas; Chelli et al. 2016).

<footnote n="18">¹⁸ https://archive.stsci.edu/missions/tess/doc/tess_drn/tessSector_31_drn47_v02.pdf</footnote>
We reduced the interferometric observations with the standard MIRC-X reduction pipeline (version 1.3.5)<fnref n="19" /> with the default reduction parameters, and we set the number of coherent co-adds (ncoh) to 10, the flux threshold to 5, and the signal-to-noise threshold to 3. We calibrated the data with a modified version of the calibration software for the previous Michigan InfraRed Combiner (Monnier et al. 2012), which allowed for the removal of bad data that the automated pipeline would not properly flag.

### 3. GP Analysis

GPs have become a frequently-used tool in analyses of RV time series since the case studies by Haywood et al. (2014) and Rajpaul et al. (2015). Specifically, they are often used as a flexible model of correlated structure attributed to stellar activity. GPs are advantageous because of their analytically tractable likelihood function, and simple parameterization through their covariance function (further details surrounding GPs may be found in Rasmussen & Williams 2006).

We recently applied GPs to high-precision EXPRESS RVs of the bright, Sun-like star HD 101501 and demonstrated that the combination of high-amplitude stellar activity and sparse observing cadence inhibits the detection of low-mass planets (Cabot et al. 2021). We showed that high-cadence RVs drastically improve the detectable parameter space for planets around active stars, and that simultaneous photometry provides important constraints on stellar activity that are otherwise difficult to infer from RVs alone. We used the same framework in the following analysis. We opted to use the celerite quasiperiodic covariance kernel (Foreman-Mackey et al. 2017) with hyperparameters {$B, C, L$, and $P_{GP}$}, corresponding to the amplitude of the covariance, the weighting of the sinusoidal term, the decay parameter, and the recurrence timescale, respectively. The covariance matrix was constructed for pairs of timestamps ($t_i, t_j$) and has the form

$$K_{ij} = \frac{B}{2+C} e^{-|t_i - t_j|/L} \left[ \cos \frac{2\pi|t_i - t_j|}{P_{GP}} + (1+C) \right], \quad (1)$$

where all the hyperparameters are positive (Foreman-Mackey et al. 2017). celerite covariance matrices may be inverted with reduced complexity, which makes their GPs appropriate for the extensive RV data set considered here. The kernel is markedly faster than the commonly-used quasiperiodic kernel (e.g., Haywood et al. 2014), which may be implemented with the `george` package (Ambikasaran et al. 2015). In contrast to the `george` quasiperiodic kernel, the `celerite` quasiperiodic kernel is not mean-square differentiable (Rasmussen & Williams 2006), and its covariance decreases faster for a fixed decay timescale parameter (denoted $L$ for the `celerite` kernel, and $\lambda_e$ for the `george` kernel). Our fitting process made use of the nested sampler `PyMultinest` (Feroz & Hobson 2008; Feroz et al. 2009; Buchner et al. 2014; Feroz et al. 2019). The sampling parameters included 4000 live points (except for our GP and 2-Planets model, which used 6000 live points), a sampling efficiency of 0.6, and an evidence tolerance of 0.5. For each model, the sampler was run three times in order to obtain a median and standard deviation on the log-evidence

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Ground-based APT light curve of $\epsilon$ Eri. Top: Differential $(b + y)/2$ photometry (in magnitudes) from the T8 (black circles) and T4 (gray diamonds, latest data set) telescopes. The long-term trends are shown in red for the T8 data and in blue for the T4 data. Bottom: The same photometry as above, but with the long-term trends removed. The signature that remains is assumed to be rotational variation and not overall brightening or dimming trends.</figcaption></figure>

In $\mathcal{Z}$. We also confirmed that the inferred parameters are consistent across the three runs (within $1\sigma$ uncertainties).

We note that Mawet et al. (2019) excluded a GP from their model after showing that it is statistically disfavored by the Bayesian information criterion, and that it does not significantly impact the derived parameters. Our study focused on characterizing stellar activity, as well as on searching for Keplerian signals that have amplitude less than that of the stellar activity signal. In particular, we were interested in the effect of rotationally-modulated signals from starspots, which can be spatially identified with a light-curve inversion algorithm (see Section 4.1). Therefore, we retained a quasiperiodic GP in our model.

### 3.1. Photometry Preconditioning

The eight-year data set of ground-based APT photometry, described in Section 2.4 and shown in Figure 1, was used to obtain posterior distributions on the GP hyperparameters. Broad, log-uniform ($\mathcal{L}U$) or uniform ($U$) priors were assigned to the GP hyperparameters $B, L, P_{GP}$, and $C$. The fit involved a jitter term $s$ added in quadrature to all uncertainties, as well as a global offset $\gamma$. The jitter and offset terms were given a broad $\mathcal{L}U$ and broad $U$ prior, respectively. The TESS photometry did not span enough rotations to constrain the GP parameters, and was not used in this analysis. Posterior samples are shown in Appendix A, along with distribution medians and 16% to 84% confidence intervals. The periodic timescale $P_{GP}$ is very well constrained at $11.4 \pm 0.2$ days. The decay timescale $L$ is found to be $40_{-10}^{+20}$ days. The constraints on these hyperparameters are based on the physical process of spots and faculae evolving on the rotating stellar surface. The parameters $P_{GP}$ and $L$ approximately correspond to the stellar rotation period at the typical latitude of spots and the typical spot lifetime, respectively. While not imposed explicitly as a prior, the results indicate the decay timescale is longer than the periodic timescale.

There is no accurate, analytical model for the influence of stellar activity on RV measurements; however, RV studies searching for exoplanets have successfully used GPs regressed

<footnote n="19">¹⁹ https://gitlab.chara.gsu.edu/lebouquij/mircx_pipe</footnote>
<table><caption>**Table 4**
GP Analysis Results</caption><thead><tr><th>Parameter</th><th>Units</th><th>GP and 1-Planet<br>(Free)</th><th>GP and 1-Planet<br>(Preconditioned)</th><th>GP and 2-Planets<br>(Preconditioned)</th></tr></thead><tbody><tr><td><i>B</i></td><td>(m s<sup>-1</sup>)<sup>2</sup></td><td>59<sup>+17</sup><sub>-8</sub></td><td>69<sup>+11</sup><sub>-8</sub></td><td>70<sup>+17</sup><sub>-8</sub></td></tr><tr><td>ln <i>C</i></td><td>...</td><td>-11 ± 7</td><td>-6.7</td><td>-6.7</td></tr><tr><td><i>L</i></td><td>m s<sup>-1</sup></td><td>21<sup>+9</sup><sub>-5</sub></td><td>36.5</td><td>36.5</td></tr><tr><td><i>P</i><sub>GP</sub></td><td>days</td><td>11.8 ± 0.4</td><td>11.4</td><td>11.4</td></tr><tr><td><i>K</i><sub>s,b</sub></td><td>m s<sup>-1</sup></td><td>10<sup>+1</sup><sub>-2</sub></td><td>10<sup>+1</sup><sub>-2</sub></td><td>10<sup>+1</sup><sub>-3</sub></td></tr><tr><td><i>φ</i><sub>0,b</sub></td><td>rad</td><td>0.8<sup>+0.4</sup><sub>-0.3</sub></td><td>0.8 ± 0.3</td><td>0.8<sup>+0.5</sup><sub>-0.3</sub></td></tr><tr><td><i>P</i><sub>b</sub></td><td>days</td><td>2650 ± 50</td><td>2650 ± 50</td><td>2650 ± 60</td></tr><tr><td><i>ω</i><sub>b</sub></td><td>rad</td><td>3 ± 2</td><td>3 ± 2</td><td>3 ± 2</td></tr><tr><td><i>e</i><sub>b</sub></td><td>...</td><td>0.01<sup>+0.06</sup><sub>-0.01</sub></td><td>0.01<sup>+0.05</sup><sub>-0.01</sub></td><td>0.01<sup>+0.06</sup><sub>-0.01</sub></td></tr><tr><td><i>K</i><sub>s,c</sub></td><td>m s<sup>-1</sup></td><td>...</td><td>...</td><td>0.1<sup>+0.7</sup><sub>-0.1</sub></td></tr><tr><td><i>φ</i><sub>0,c</sub></td><td>rad</td><td>...</td><td>...</td><td>3 ± 2</td></tr><tr><td><i>P</i><sub>c</sub></td><td>days</td><td>...</td><td>...</td><td>80<sup>+700</sup><sub>-80</sub></td></tr><tr><td><i>ω</i><sub>c</sub></td><td>rad</td><td>...</td><td>...</td><td>3 ± 2</td></tr><tr><td><i>e</i><sub>c</sub></td><td>...</td><td>...</td><td>...</td><td>0.03<sup>+0.27</sup><sub>-0.03</sub></td></tr><tr><td>ln <i>Z</i></td><td>...</td><td>-3289.24 ± 0.02</td><td>-3287.76 ± 0.10</td><td>-3287.12 ± 0.13</td></tr><tr><td>ln <i>L</i><sub>MAP</sub></td><td>...</td><td>-3227.4</td><td>-3231.5</td><td>-3220.8</td></tr></tbody></table>

**Note.** Results of our *celerite* GP analysis on the full RV data set. The three models correspond to GP and 1-Planet in which all GP hyperparameters were left free; GP and 1-Planet where three GP hyperparameters were fixed to their MAP values from the light-curve analysis; and GP and 2-Planets, where again three GP hyperparameters were fixed. The columns contain the median of the marginalized distribution of each parameter of interest, and the uncertainties correspond to the 16th and 84th percentiles. Values in parentheses “(O)” denote the MAP values. The bottom rows contain the log-evidences returned by the nested sampler and the log-likelihood of the MAP vector. The value and uncertainty of each ln *Z* represent the median and standard deviation of three separate sampler runs, respectively. Jitter and offset parameters for each time series are omitted from the table. Note that for the 2-Planets model, a significant number of samples involve the sampler effectively swapping between planet b and planet c. Rather than enforcing a prior to maintain ordering between the planets, we simply filtered samples for the statistics reported above for this model. Orbital parameters correspond to samples with *P*<sub>c</sub> < 2000 days.

to the RV measurements (Haywood et al. 2014) to model the activity component. Since RVs are usually sparser than photometry, it is useful to first regress a GP to photometry and determine the best-fit hyperparameters (as we have done above), and then regress another GP to the RV measurements where most of the hyperparameters are fixed to the photometry-based best-fit values (Haywood et al. 2014). We adopted the same strategy of Haywood et al. (2014) in our RV analysis by fixing three GP hyperparameters to their maximum a posteriori (MAP) values from the photometry fit: log *C* = -6.7, *L* = 36.5 days, and *P*<sub>GP</sub> = 11.4 days. These fixed hyperparameters are listed in Table 4 for the “preconditioned” models. The amplitude hyperparameter *B* was left free, since the photometry and RVs naturally concern different units. Kosiarek & Crossfield (2020) demonstrated that the best-fit hyperparameters for quasiperiodic GPs are similar when regressed separately to contemporaneous RVs and photometry of the Sun. While this result supports the notion of preconditioning with photometry, a variety of factors such as spot distributions, instrument systematics, or non-contemporaneous time series could plausibly result in different best-fit hyperparameters. To check the efficacy of preconditioning GP hyperparameters on the photometry fit, we explored a GP and 1-Planet model that left all GP hyperparameters free (Table 4, first column). The orbital parameters are largely unaffected. Interestingly, the inferred periodic timescale *P*<sub>GP</sub> = 11.8 ± 0.4 days is very close to the photometry preconditioning best-fit value (*P*<sub>GP</sub> = 11.4 days), demonstrating a strong periodic component in the RV data set arising from stellar rotation, and that the constraint proposed by Kosiarek & Crossfield (2020) is unnecessary in our case. The decay timescale *L* = 21<sup>+9</sup><sub>-5</sub> days is roughly half of the best-fit GP photometry value (*L* = 36.5 days), allowing the

GP to more rapidly evolve between subsequent stellar rotations; however, the two estimates overlap at their 1σ upper and lower limits.

## 3.2. RV Analysis

We proceeded to model over 30 yr of RVs, with a focus on characterizing stellar activity, obtaining tight constraints on the orbit of the known planet *ε* Eri b, and searching for additional planets. Mawet et al. (2019) showed that secular acceleration of *ε* Eri has negligible impact on derived parameters, and eventually excluded it from their model. Additionally, they partitioned the Hamilton/Lick RVs (Fischer et al. 2014) into four distinct data sets on account of instrument upgrades. We repeated these decisions in our analysis. In total, we jointly analyzed 11 different RV data sets, which include distinct time series for each Lick epoch and the two cameras used in conjunction with Coudé Echelle (see Figure 2). Therefore, our entire RV model consists of the following: a GP activity component, which contributes one free hyperparameter *B* corresponding to the covariance amplitude; a systemic velocity offset *γ*<sub>k</sub> and jitter term *s*<sub>k</sub> added in quadrature with each data point’s uncertainty, for each RV data set *k* ∈ {1, 2, 3...11}; and five orbital elements {*K*<sub>s</sub>, *φ*<sub>0</sub>, *P*, *ω*, and *e*} for each Keplerian component, corresponding to the semi-amplitude, the phase of the first epoch, the orbital period, the longitude of the periastron, and the eccentricity, respectively. The prior on *B* is a log-uniform probability density function distribution ranging from 0.1 to 2500, or *L*U(0.1, 2500) (m s<sup>-1</sup>)<sup>2</sup>. The priors on all *γ*<sub>k</sub> and *s*<sub>k</sub> are *L*U(-30, 30) m s<sup>-1</sup> and *L*U(0.01, 20) m s<sup>-1</sup>, respectively. The semi-amplitude *K*<sub>s,b</sub> prior is *L*U(0.1, 20) m s<sup>-1</sup>, and the orbital period *P*<sub>b</sub> prior is *L*U(2, 10,000) days.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Best-fit GP and 1-Planet model that was fit to the combination of archival $\epsilon$ Eri RVs with new EXPRES RVs (dark blue scatter points). The archival data set is identical to that analyzed by Mawet et al. (2019), with the addition of CHIRON RVs (Giguere et al. 2016). For clarity, the MAP offset values ($\gamma_k$) have been subtracted from each RV data set. A zoomed-in panel of the EXPRES RVs is shown on the right. The GP mean (red line) predominantly tracks the 7 yr planetary signal (black line), with deviations owing to stellar activity. The quasiperiodic activity signal is more clearly seen in the zoomed-in panel. The period of oscillations is usually close to the periodic timescale hyperparameter, which was fixed to the 11.4 day stellar rotation period. By choosing a quasiperiodic kernel, the GP can accommodate small variations in period and amplitude, as well as gradual changes in the activity signal's phase. The characteristic timescale for these variations is set by the parameter $L$. The GP $1\sigma$ and $2\sigma$ confidence intervals are depicted as lightly shaded regions around the GP mean. CES+LC and CES+VLC correspond to data obtained with the Coudé Echelle Spectrograph with the Long Camera and Very Long Camera, respectively. Lick data are color-coded by their corresponding upgrade epoch.</figcaption></figure>

The fit results for the parameters of interest are summarized in the second column of Table 4. Our fitted orbital parameters for the 7 yr planet agree with those of Mawet et al. (2019) within $1\sigma$ uncertainties. The fitted RV curve is shown in Figure 2, where $1\sigma$ and $2\sigma$ confidence regions are derived analytically from the GP posterior predictive distribution, with the hyperparameters fixed at the MAP values. We confirm that the orbit is consistent with a circular one (Mawet et al. 2019). We fit an additional model of two planets with preconditioning of GP hyperparameters (Table 4, third column). We fixed the GP hyperparameters $L$, $P_{GP}$, and $\log C$ to photometry-derived values in order to reduce the fit dimensionality, and adopted the same priors on the second Keplerian parameters as used for the first Keplerian component. The Bayesian evidence prefers the two-planet model marginally, but insufficiently to justify the additional Keplerian parameters. The second planet's semi-amplitude is consistent with zero.

In Section 5, we more intensively study the activity-induced signal in the recent EXPRES RVs. The MAP orbital solution in the one-planet, preconditioned model (Table 4, second column) serves as the RV baseline.

## 4. Photospheric Brightness Variations

Here, we model the spotted stellar surface of $\epsilon$ Eri with two different techniques: light-curve inversion of TESS photometry and spot modeling with interferometry from the CHARA Array using the MIRC-X beam combiner. Both techniques assume the star is spherical, as distortions are not expected for this star with $v \sin i = 2.93 \pm 0.5$ km s$^{-1}$ (Giguere et al. 2016).

### 4.1. Light-curve Inversion Reconstructions

While ground-based photometry (Section 2.4) is available for $\epsilon$ Eri, we focus here only on the 2 minute cadence TESS light curve, which has significant overlap with both the EXPRES observations and our two MIRC-X observations.

The Sector 31 TESS light curve encompasses just over two rotations of $\epsilon$ Eri based on our $P_{rot} = 11.4$ days. At the gap in data between the two orbits that make up Sector 31, we split the data into two rotation periods (2,459,144.5196172–2,459, 155.8934423 and 2,459,158.8670810–2,459,169.9489975). The observations were averaged in 100 equally sized bins in phase across each of the two rotation periods. We reconstructed each rotation separately with an algorithm Light-curve Inversion (LI; Harmon & Crews 2000). LI makes no *a priori* assumptions of the starspot shape, number, or size and uses a modified Tikonhov regularizer to reconstruct the stellar surface (for more details on LI see Harmon & Crews 2000; Roettenbacher et al. 2011, 2013). We provided the algorithm input parameters of effective temperature $T_{eff} = 5100$ K (Giguere et al. 2016), starspot temperature $T_{spot} = 4100$ K (based on Berdyugina 2005), quadratic limb-darkening coefficients $a = 0.4258$ and $b = 0.1936$ (Claret 2018), and stellar inclination $i = 70^\circ$ (Giguere et al. 2016), where equator-on viewing corresponds to $i = 90^\circ$. We estimated the spot-to-photosphere brightness ratio by integrating the radiation of blackbodies with the temperatures of the starspot and the photosphere over the spectral response function of the TESS bandpass (as in Roettenbacher et al. 2013). For $T_{spot} = 4100$ K and $T_{eff} = 5100$ K, the brightness ratio is 0.40, i.e., the spot has 40% of the brightness of the photosphere. LI finds the optimum reconstruction that simultaneously fits the prescribed spot-to-photosphere brightness ratio and the prescribed rms difference between the observed and model light curves. We required the rms difference for these observations to be 0.0005 (a value comparable to that used for applying LI to a Kepler light curve; Roettenbacher et al. 2016). Because of the small amplitude of the TESS light curve, this rms allows for the systematic deviations between the observed and reconstructed light curves shown in Figure 3. Reducing the rms value to improve the fit resulted in surfaces with features that are characteristic of overfitting (e.g., elongated dark features that are not consistent with structures analogous to starspots; Harmon & Crews 2000). The light curve of $\epsilon$ Eri evolves on a timescale shorter than the stellar rotation
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Top: Binned TESS light curve (blue solid line) and LI-reconstructed light curve (red dashed line) of ϵ Eri. The binned light curve is normalized to the maximum value for each rotation. Two stellar rotations were reconstructed using LI with the gap in TESS data separating the rotations. The vertical black dashed lines indicate the times of the interferometric MIRC-X observations discussed in Section 4.2. Bottom: LI-reconstructed surfaces of the binned TESS light curve. The visible pole is marked with a white line and the surface rotates counterclockwise around this pole. The surface temperature ranges from a spot temperature of T<sub>spot</sub> = 4100 K to a photospheric temperature of T<sub>phot</sub> = 5100 K. The surfaces are shown at the times of the MIRC-X observations.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Left: LI-reconstructed pseudo-Mercator surface of ϵ Eri for the first rotation observed by TESS (ending before TJD 2,459,157). At phase 0.0, as viewed by TESS, the center of the star is at longitude 90°. As time increases, the longitude of the center of the star decreases. Right: LI-reconstructed pseudo-Mercator surface of ϵ Eri for the second rotation observed by TESS (starting after TJD 2,459,157).</figcaption></figure>

period, so we emphasize that LI assumes that the light curve does not evolve in the light curve being inverted.

In Figure 3, we present the two LI-reconstructed surfaces from Sector 31 at the time of the MIRC-X observations. Each of the two surface maps shown in Figure 3 features one prominent starspot. The star rotates 98° from left to right between the two MIRC-X observations, so the spot seen on the left surface map is the same as the spot near the right limb in the right surface map. We also present Figure 4, which shows pseudo-Mercator projections of the entire stellar rotation from each portion of the TESS light curve. In this presentation, JD 2,459,144.5196172 is considered phase 0.0, which corresponds to the longitude of 90° as the center of the star, as seen by the

observer. As time increases, the longitude at the center of the star decreases because the star is assumed to rotate counterclockwise as seen from above the visible pole, while the stellar longitude also increases in the counterclockwise direction. A second starspot is visible in both reconstructions.

## 4.2. Interferometric Models of Starspots

For the two sets of interferometric data, we individually fit a stellar surface to each one. The limb-darkening coefficient, the starspot location, and the starspot size were fit separately to the data, as described below.
For limb darkening, we assumed power-law limb darkening, which is defined as

$$I(\mu) = I_0 \mu^{\alpha}, \qquad (2)$$

where $\mu$ is the cosine of the angle between the observer and the normal to the stellar surface, $I$ is the intensity, $I_0$ is the intensity at the center of the stellar surface, and $\alpha$ is the limb-darkening coefficient.

For both sets of data, the starspot is defined to have a brightness of 61% of that of the photosphere. This value was estimated based on the ratio between the estimated spot and photospheric temperatures in the H band, the bandpass in which our MIRC-X observations were obtained. The difference between the spot and photospheric temperatures is the same as was used with LI and was estimated based on Berdyugina (2005).

To construct a model of the star, we assumed the H-band, limb-darkened angular diameter measured by Baines & Armstrong (2012), $\theta_{LD} = 2.153$ mas. While Baines & Armstrong (2012) used a linear limb-darkening law (Milne 1921) and a model limb-darkening coefficient from Claret et al. (1995), we used a power-law limb-darkening law, as described above, which was shown to be an appropriate model for interferometric observations by Lacour et al. (2008) and described by Hestroffer (1997). To determine the best value of $\alpha$ for limb darkening, we fit a surface without spots to the data. We selected $\alpha$ with the fit that had the lowest reduced $\chi^2$ value for the combination of the visibilities, closure phases, and triple amplitudes with the visibilities being weighted 10 times more strongly than the other parameters to allow the closure phases and triple amplitudes to contribute to but not dominate the fit. The visibility measurements were favored because the shape of the visibility curve more strongly constrains the limb-darkening parameter than the closure phases, which more strongly constrain surface asymmetries. For the 2020 November 2 data, we found a best-fit limb-darkening coefficient of $\alpha = 0.27 \pm 0.02$. For the 2020 November 5 data, we found $\alpha = 0.27 \pm 0.01$. Errors were determined with 100 bootstraps. For example, for each bootstrap, from the 985 visibility data points for 2020 November 2, a point was randomly chosen 985 times with replacement, and this set, as well as similarly selected closure phases and triple amplitudes, was used to fit for $\alpha$.

We then assumed the stellar diameter and the limb darkening and performed a grid search for the location of a simple, circular starspot with varying size to obtain the best fit for its location, which we show in Figure 5. Ideally, to choose the best-fit spot location, we would select the location with the lowest reduced $\chi^2$ value for the combination of the visibilities, closure phases, and triple amplitudes with the closure phases being weighted 10 times more strongly than the other parameters in this case because the closure phases are more sensitive to surface asymmetries than the other observables. Here, the weighting was chosen to allow the visibilities and triple amplitudes to contribute to but not dominate the fit. While this best-fit method is valid for the night of 2020 November 2, when the starspot was near the center of the stellar disk and at nearly its largest contribution to the light curve, the starspot was on the limb on 2020 November 5 and at a much lower contribution to the light curve. The best-fit spot location for 2020 November 5 is likely an artifact, but there is a local minimum that is consistent with the expected location of

the spot based on its previous location and the star having rotated through 98° between the observations. The local minimum spot location is shown in Figure 5. The details of our recovery tests for simulated starspots using the uv plane coverage of the 2020 November 5 data are in Appendix B. We note that the limited data available for these models makes these spot detections marginal, though consistent with the photometric observations.

For the 2020 November 2 data, we performed 100 bootstraps, as described above, for the location of the starspot of the same size and darkness used for the best-fit data set to illustrate the quality of our fit. The starspot locations were not easily quantified to error bars on the starspot location; therefore, we show the location of the starspots in Figure 9 in Appendix B. The locations of the starspots in these bootstraps are roughly consistent with the region surrounding the minimum of the $\chi^2$ surface of the best-fit data set (see Figure 10 also in Appendix B).

### 4.3. Reconstruction and Model Differences

While the TESS light curve and, consequently, the LI surface reconstructions of Figure 3 show more starspot evolution than is detectable in the interferometric models of Figure 5, we note that the TESS light curve is missing nearly all observations between the two interferometric observations. As a result, the reconstruction of the evolving spot that is visible at the time of both MIRC-X observations is more informed by the prior rotation (2020 November 2) and the following rotation (2020 November 5). The LI surfaces give the impression that the morphology of the spot has changed in the 3 days between the MIRC-X observations. While we assume some evolution has occurred, the differences between the spots as seen in Figure 4 are potentially misleading and likely due to noise artifacts of the reconstructions. Unfortunately, we are unable to resolve structure in the starspots with the interferometric observations for comparison.

A notable difference between the LI reconstructions and the interferometric models is the locations of the spots. LI has no ability to either constrain or determine the position angle of the rotational pole in the plane of the sky. Interferometric observations, however, do provide the opportunity to constrain the star's orientation on the sky, both the position angle and the inclination, given sufficient data. Consistent with our LI reconstructions, we assumed an inclination of $i = 70^\circ$ (Giguere et al. 2016) for the interferometric observations. We investigated the position angle of the rotational pole of $\epsilon$ Eri by comparing the interferometric models to the surfaces reconstructed with LI. The starspot longitudes taken from LI are reliable, but the starspot latitude may not be reliable due to a lack of sufficient information from the input light curve(s). Because we only used a single bandpass of data for these reconstructions, we assumed that the starspots are not necessarily located in the latitude where they appear on the LI surfaces.

To approximate the position angle of $\epsilon$ Eri in the plane of the sky, we superimposed the interferometric models with the LI-reconstructed surfaces rotated in the plane of the sky. For each snapshot interferometric data set, there are two orientations for which the modeled starspot and the corresponding LI-reconstructed starspot are aligned in longitude; however only one orientation is consistent with both interferometric models. Using the reconstructions and models for both 2020 November
<figure><img src="image_9.png" /></figure>

**Figure 5.** Top: MIRC-X observations of ϵ Eri from 2020 November 2 (left) and 2020 November 5 (right). The observed data are shown in black with one baseline of each triangle shown in the closure phase and triple amplitude panels. The model data for a limb-darkened, spotted surface are plotted in red. Bottom: The spotted, limb-darkened (power-law) stellar surface models used to generate the interferometric data in red. The visible pole is marked with a white line and the surface rotates counterclockwise around this pole. As in Figure 3, the surface temperature ranges from a spot temperature of T<sub>spot</sub> = 4100 K to a photospheric temperature of T<sub>phot</sub> = 5100 K. The surface models are presented such that they are on the plane of the sky where north (toward the celestial pole) is up and east is to the left. Because of the sparse uv plane coverage, these starspot detections are marginal. More data are needed for future studies aiming to image the star and determine the inclination and position angle of the rotational pole.

<figure><img src="image_10.png" /></figure>

**Figure 6.** Left: The best-fit interferometric fit for the 2020 November 2 data set with the LI reconstruction and its latitude and longitude grid overlaid, assuming *i* = +70°. The LI reconstruction is rotated ~335° east of north in order to have the LI and MIRC-X starspots lie on the same line of longitude. The visible pole is marked with a white line and the surface rotates counterclockwise around this pole. Right: Same as the left figure, but for the 2020 November 5 data set. This position angle determination is tentative due to the limited interferometric data.

2 and 5, we found that the data suggest that the position angle of the rotation axis of the star is ~335°, as presented in Figure 6. This assumes that the stellar inclination is *i* = +70°, but the method used to determine the inclination in Giguere et al. (2016) and LI both cannot make a distinction between *i* = +70° and *i* = −70°. Therefore, a position angle of 155° is also possible if *i* = −70°. The position angle suggested in this analysis is tentative and must be confirmed with future

data sets. Further details on this determination are found in
Appendix B.

**5. Spectroscopic Simulations**

Sector 31 TESS photometry spans nearly 2.4 stellar rotations of ϵ Eri when it was also observed with EXPRESS. From these photometric observations, we reconstructed the stellar surface
(Figures 4) and, here, we discuss our method for identifying and isolating the RV signature of the starspots at any phase of rotation.

## 5.1. Integrated Spectrum over a Simulated Disk

We simulated the stellar spectrum over time as the surface inhomogeneities reconstructed in the previous section rotated into and out of view for the time baseline of the TESS Sector 31 observations. This timeframe overlaps with 32 (of 164 total) spectra of $\epsilon$ Eri obtained with EXPRES between 2020 October and November. Our simulation is similar to the SOAP software (Boisse et al. 2012), which models perturbations to the spectral line profile induced by circular spots and plages, while accounting for limb darkening, stellar rotation, and stellar geometry. Additional physics, including an improved limb darkening and a model of convective blueshift suppression, were incorporated into SOAP2 by Dumusque et al. (2014). In a similar vein, 2D simulations of stellar disks have been used in studies of exoplanet transits (Casasayas-Barris et al. 2019), with the purpose of investigating spectroscopic artifacts from center-to-limb variations and the Rossiter-McLaughlin effect. The premise is the same, in that the exoplanet occults a region of the star with its own local properties (e.g., white-light flux and Doppler shift). Our simulation code uses the physics and methodology used in SOAP2, but with EXPRES spectra for $\epsilon$ Eri instead of the solar spectrum.

A high-fidelity template spectrum of $\epsilon$ Eri was obtained by simultaneously fitting all 164 EXPRES spectra of the star with a B-spline regression similar to what is done by SERVAL (Zechmeister et al. 2018). The B-spline is crucial for providing a smooth, continuous function to which we may apply arbitrarily small Doppler shifts. The disk model itself comprises an 80 × 80 pixel grid representing the visible surface of the star. Corner pixels beyond 1 $R_*$ remain empty throughout the following steps. Each pixel was assigned (1) a finely-sampled ($R \sim 800,000$) spectrum derived from the spectral template; (2) a flux weight (i.e., the relative contribution of flux to the integrated spectrum due to limb darkening; the effect of spots is discussed below); and (3) a local velocity determined by stellar rotation. We approximated the spectra for various values of $\mu$ by assigning flux weights according to an appropriate limb-darkening law and used the appropriately Doppler-shifted EXPRES template spectrum. To generate a disk-integrated spectrum at each moment in time, we co-added the Doppler-shifted flux-weighted spectra from every pixel in the grid.

Critical to this study, the model is based on the LI surface reconstructions of the TESS light curve described in the previous section. First, the LI surface was interpolated from the existing spatial structure (approximately equal-area rectangular zones) onto a uniformly-spaced latitude/longitude grid (Figure 4). We found 90 latitude divisions and 180 longitude divisions were sufficient for resolving the finest details in the LI surface. The rotation of the star was simulated by transforming the latitude and longitude coordinates, taking into account stellar inclination. Finally, the surface was mapped from 3D spherical coordinates to a 2D projection (via interpolation) on the pixel grid described above. At a given time step, we multiplied the flux weights of each pixel (determined by a quadratic limb-darkening law) by the relative brightness of the projected stellar surface, and then integrated the spectrum. In this way, we accounted for spots downweighting local contributions to the integrated spectrum. As a consistency check, we summed the white-light flux of the projected disk at

30 equally-spaced intervals during a single rotation and recovered the relative brightness variations in the TESS light curve. The effects of resolution were explored, and we found the results did not change appreciably when more pixels were used to model the disk. RVs were obtained from the integrated spectrum via cross-correlation with the original stellar template.

Our LI surface model for the TESS Sector 31 data is shown in Figure 7 and compared to both the actual EXPRES data, the GP model described in Section 3, and an $FF'$ model (Aigrain et al. 2012) discussed in the next section. The Keplerian component of the known planet has been subtracted, but this represents only a marginal RV trend of $\sim 18 \text{ cm s}^{-1}$ over the $\sim 20$ day timespan of the RVs considered. The RVs from the LI model contain relative velocity variations that arise from perturbations to the spectral line profile caused by the simulated spot. The actual EXPRES RVs are also relative velocities; however, the relative velocity offset in the EXPRES data changes as new data are acquired and will only asymptotically approach a constant value after all of the signals (i.e., the known planet in a 7 yr orbit and the photospheric contributions) have been well sampled over all phases. The temporal baseline of the EXPRES RVs, which is slightly over 1 yr, is not sufficient to reach that constant offset, so we derived and removed the best-fit offset between the simulated and observed RVs.

The subset of EXPRES RVs taken during the TESS Sector 31 observations has an rms scatter of $4.72 \text{ m s}^{-1}$. If we adopt corrections from the GP model, this scatter is reduced to $0.76 \text{ m s}^{-1}$. However, the GP model, by definition, is conditioned on the data to which it is fit. Therefore, the residuals reflect the measurement uncertainties added in quadrature with the jitter term $s_k$, added to account for intra-night scatter. Intra-night $s_k$ can arise from p-mode oscillations, granulation, or underestimation of formal uncertainties. This illustrates one of our primary concerns: even when the hyperparameters are conditioned on photometric observations, the GP's flexibility (which is based on its parameterization and likelihood) allows it to conform to the RV measurements very closely, typically within $<1 \text{ m s}^{-1}$ from individual measurements or the mean RV of a given night.

The LI surface model, however, is completely independent of the EXPRES RVs. Aside from the constant RV offset discussed above, the LI surface model does not represent a fit to the RV data, but was derived strictly from the TESS photometry and the stellar spectrum template. After subtracting the RVs derived with the LI model and interpolated to the time of the EXPRES RV measurements, the residual rms of the EXPRES RVs decreases from $4.72$ to $1.98 \text{ m s}^{-1}$. This is a 58% reduction in the RV scatter for $\epsilon$ Eri.

# 6. Discussion

## 6.1. Efficacy of the LI Surface RV Model

We have shown that the LI surface stellar activity model presented in Section 5 accounts for a significant portion of scatter in the EXPRES RV measurements. The success of this model marks an important step in robustly separating Keplerian and activity-related components in an RV time series; however, it is important to address the remaining $\sim 2 \text{ m s}^{-1}$ rms scatter and potential avenues for improving the model. For example, Dumusque et al. (2014) accounted for additional physics in their spot models, including inhibition of convective blueshift in regions affected by spots. Qualitatively, the net effect of
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Comparison of three models of activity-based RV variations. The red curve and shaded regions depict the GP model and confidence intervals (also shown in Figure 2). The orange curve depicts the RV variations obtained by rotating the LI-reconstructed stellar surface and integrating the stellar spectrum over a pixelated disk. For reference, snapshots of the disk model are shown for several selected timestamps of RV exposures. The features on the surface are responsible for the RV variations in the model, and a faint dashed line connects each snapshot to the corresponding point on the RV model curve. The gray curve represents an *FF'* model (Aigrain et al. 2012), which serves as a benchmark for the LI surface model. The bottom panel shows the residuals for the LI surface model, GP model, and *FF'* model in orange, red, and gray, respectively. The rms scatter of residuals after subtracting the GP is 0.76 m s⁻¹, which is largely a reflection of the EXPRESS measurement uncertainties and intra-night scatter. The rms scatter of residuals after subtracting the LI surface model and best-fit offset is 1.98 m s⁻¹, compared to the 2.20 m s⁻¹ rms from subtracting the *FF'* model. All three are a reduction from the original rms scatter of 4.72 m s⁻¹ of this subset of RVs, which were selected based on their overlap with TESS photometry. In both the plot and the rms calculations, the marginal Keplerian contribution has been subtracted.</figcaption></figure>

accounting for convective blueshift in an equatorial spot is to break the symmetry in its corresponding RV signal, and push the RV signal toward more positive velocities. We experimented with a simple implementation of the convective blueshift effect by assigning a constant velocity offset inside the active region (e.g., Dumusque et al. (2014) adopted 350 m s⁻¹), but this addition did not improve the model fit. It is worth noting that $\epsilon$ Eri is a K2 dwarf, and the interplay between convective cells and magnetic fields in active regions may differ from that in the Sun. It might also help to more carefully model spectra at different values of $\mu$. For example, Casasayas-Barris et al. (2019) generated synthetic spectra for regions extending from the disk center to the limb. The change in effective temperature may have a minor impact on the cross-correlation and RV inference and it would be useful to retrieve faculae and accurate sizes for high-latitude features in the surface reconstructions; however, this will require longer-baseline interferometric measurements than are presently available. The most prudent next step is to improve the cadence of the RVs to allow for more dynamic modeling of the surface and to extend the time baseline of observations for both the RVs and the space-based photometry.

It is also worthwhile to explore additional physics. The full impact of stellar activity itself involves many components, including features on the rotating surface, as well as granulation, *p*-mode oscillations, and magnetic activity cycles (Fischer et al. 2016, and references therein). While the amplitudes and timescales of different activity sources have been studied, the details needed to model their impact on RVs are not well understood. Nevertheless, the growing body of literature

adopting quasiperiodic GPs as an activity model supports the premise that, on timescales of several days, rotationally-modulated signals often have the strongest amplitude effect on RVs, and are the most likely to generate false-positive Keplerian candidates. There is also strong precedent for inferring rotationally-modulated RV variations from photometry (Aigrain et al. 2012; Haywood et al. 2014).

## 6.2. Comparison to GPs

We modeled the full RV time series in Section 3 for the purpose of independently measuring orbital parameters of $\epsilon$ Eri b, evaluating evidence for additional planets, and characterizing stellar rotation and typical spot lifetimes. Our model employed a GP via the *celerite* implementation. GPs are being applied as flexible stellar activity models by other groups as well in the current era of extreme-precision spectroscopy (e.g., Faria et al. 2020). Our application of GPs closely resembles that of Haywood et al. (2014); however, one may simultaneously model activity indicators as done by Rajpaul et al. (2015). A number of even more advanced models are described in contributions to Zhao et al. (2021, in preparation), and a new GP model for inferring parameters governing the distribution of starspots, with potential applications to RV data sets, has been presented by Luger et al. (2021).

A particular advantage of the LI surface model is its derivation from an independent, photometric data set, which eliminates the possibility of Keplerian signals being absorbed into the model. We previously showed that photometry-preconditioned GPs are prone to conform to low-amplitude variations from a planet if that
Keplerian signal is rejected under the Bayesian evidence comparison (i.e., adding five additional parameters to the model does not sufficiently improve its log-likelihood; Cabot et al. 2021). This effect is especially true for sparse RV data sets. Another advantage is that the LI surface model has a direct, clearly interpretable correspondence to a physical characteristic of the star (i.e., resolved surface features), while GPs typically do not. One may use specific, physically-motivated GP kernels (e.g., quasiperiodic ones) and find a repeating structure that corresponds to long-lived surface features, but it is a degenerate problem to invert the GP model and resolve surface features or other specific qualities of the star. In our case, the GP clearly does not correspond to the effects of surface features at all times and deviates at a >2σ level from the LI surface model at most times (Figure 7). GP models that fit multiple, contemporaneous time series at once (e.g., Rajpaul et al. 2015; Gilbertson et al. 2020) in implementations such as pyaneti (Barragán et al. 2022) were not tested in this study, but may make for interesting future comparison against the LI model. Compared to the photometry-preconditioned GP used here, they may have reduced flexibility.

At present, GPs remain a useful model for stellar activity, which can greatly assist with identifying and constraining Keplerian signals. However, accurately modeling RV variations with alternative and complementary data sets is an important goal, both for robustness against inadvertently removing Keplerian signals and for the interpretability of the activity model.

### 6.3. Comparison to FF'

The FF' technique (Aigrain et al. 2012) models RV variations based on contemporaneous photometry and provides a useful benchmark for our LI surface model. FF' exploits the geometry of a spot moving across the stellar surface to remove explicit dependence on the rotation period, stellar inclination, and spot latitude. However, it is accurate only to first order in the presence of multiple spots, and neglects limb darkening and spot projection effects. Extensions of the FF' technique include those of Rajpaul et al. (2015) and Giguere et al. (2016), which are not explored here. The full FF' model is

$$ \Delta RV = \frac{\Psi(t)}{\Psi_0} \left(1 - \frac{\Psi(t)}{\Psi_0}\right) \frac{R_*}{f} + \left(1 - \frac{\Psi(t)}{\Psi_0}\right)^2 \frac{\delta V_c \kappa}{f}, \quad (3) $$

or equivalently

$$ \Delta RV = \Delta RV_{\text{rot}} + \Delta RV_c, \quad (4) $$

by denoting the left term as the rotation-related component and the right term as the convection-related component. In the above equation, $Ψ(t)$ is the light curve, and the constants $Ψ_0$, $δV_c$, $f$, $κ$, and $R_*$ represent the disk's flux if no spots are present, the convective blueshift inhibition within a magnetized region, the reduction in flux for a spot at the center of the disk, the ratio between the areas of the magnetized region and spot surface, and the stellar radius, respectively. Two parameters may be estimated directly from the light curve: $Ψ_0 ≈ Φ_{max} + σ_Ψ$, where $σ_Ψ$ is the standard deviation of the light curve and $Φ_{max}$ is the light-curve maximum; and $f ≈ (Ψ_0 - Φ_{min})/Ψ_0$, where $Φ_{min}$ is the light-curve minimum (Aigrain et al. 2012).

The binned TESS light curves (Figure 3) were interpolated onto an oversampled grid and subsequently smoothed with a Savitzky-Golay filter, which allowed us to compute a smooth

time derivative of the light curve $Ψ(t)$. The window length was approximately 6% of an 11 day TESS observing window. We manually varied the window size, and found the final residual rms scatter changed up to ~20 cm s⁻¹. The adopted window size yielded the lowest residual rms scatter. We optimized $δV_cκ$ and a global model offset to fit the model to the subset of EXPRES data in Figure 7. The convection term in Equation (4) did not improve the fit, so $δV_cκ$ was fixed to 0 in the best-fit model. The convection term was significantly smaller than the rotation term in a case study of HD 189733 (Aigrain et al. 2012), a star with a similar rotation period to $ϵ$ Eri; however, suppression of convective blueshift may become the dominant process for slower rotators (e.g., Haywood et al. 2014).

The residual rms for the best-fit FF' model is 2.20 m s⁻¹. The FF' model has considerably more fine structures than the LI surface model (Figure 7), which is due to either noise in the TESS light curve or smaller surface features unresolved by the inversion. However the largest features in both models are similar, which is expected since coincidentally only one spot is visible at most times, and the spots are small as compared to the stellar disk, although the two models deviate by a few meters per second midway through the first rotation, which may be due to the noncircular morphology of the visible spot. In addition to moderately outperforming FF' by a ~20 cm s⁻¹ reduction in rms, the LI surface model is more easily interpreted since it reveals the correspondence between specific spots and their RV perturbations. The remaining scatter in both sets of residuals indicates that more accurate modeling is needed for certain spot distributions (e.g., high-latitude features, or an odd-numbered multipole component), or that other features (e.g., faculae) or physics (e.g., granulation) are responsible for a significant portion of the activity signal.

### 6.4. Imaging Sun-like Stars

While LI reconstructions are valuable in revealing the stellar surface's RV contribution, the surfaces are affected by the degeneracies of the light-curve inversion method. The longitude of starspots is well constrained by the light curve; however, the starspot latitude is not. With this method, information on constraining stellar latitudes comes from limb darkening. Starspots at different latitudes will impact the light curve differently at different wavelengths because limb darkening is different in different bandpasses. For our $ϵ$ Eri inversions, we used a single-bandpass light curve from TESS. Improving the latitude constraints of starspots being reconstructed with LI will require simultaneous light curves in a range of photometric filters.

As stated in Section 1, reconstructing a stellar surface with Doppler imaging will provide more latitude information. Although the method is being utilized for searching for hot Jupiters orbiting young stars (Heitzmann et al. 2021), Doppler imaging is not an appropriate imaging method for main-sequence stars, like $ϵ$ Eri, because the stars' slow rotation does not provide sufficient spatial resolution.

To date, interferometric aperture synthesis imaging, which unambiguously provides latitude information, has not been performed on main-sequence stars. However, our models of two data sets from the CHARA Array with the MIRC-X beam combiner indicate that detailed images of bright, spatially large main-sequence stars are possible with sufficient uv plane coverage and prominent spots.

Between the two interferometric observations, $ϵ$ Eri rotated in approximately a quarter rotation with the starspot present on
both nights of the interferometric data being in the same spot. Our models potentially constrain stellar orientation. The interferometric data suggests that the star's rotation axis is oriented along the position angle ~335° east of north. This was determined by rotating the LI reconstruction such that the MIRC-X observations would have the same longitude (see Section 4.3 and Appendix B). As we have only a small number of interferometric observations, confirmation of this position angle requires further observations.

The starspot was observed approximately 20° of latitude away from the starspot in the LI reconstructions. This is in line with the expectation that surfaces reconstructed with LI will not necessarily place the starspots at the appropriate stellar latitudes. LI favors smaller spots, which the regularizer identifies as “smoother”; in combination with only one bandpass with limited limb-darkening constraints, LI will reconstruct the surface with starspots at the sub-Earth latitude. For the 70° inclination of ϵ Eri, the sub-Earth latitude is +20°.

While the starspot latitude appears to be lower than what LI predicts, we did not use this information to inform the simulated spectra, as our interferometric model for 2020 November 5 was not independently determined, and was informed by the previous observation and the LI surface due to the low contrast between the starspot and the limb (for more details, see Appendix B). However, these observations serve as a proof of concept that starspots can be interferometrically detected on a main-sequence star.

Because the position angle was not independently determined for both nights of observation, we do not definitively state the star's orientation on the sky. We also assumed an inclination of 70° after Giguere et al. (2016). The debris disk around ϵ Eri has been detected in a number of different studies that report inclinations ranging from about 20° to 30° (Greaves et al. 2014; MacGregor et al. 2015; Chavez-Dagostino et al. 2016; Booth et al. 2017; Holland et al. 2017). All of those studies give position angles of ~0°, except for Holland et al. (2017), who give 61° ± 3° east of north. The inclination of the plane of the planet is consistently found to be ~30° (Hatzes et al. 2000; Benedict et al. 2006; Reffert & Quirrenbach 2011). The longitude of the ascending node of the orbit was found to be 254° ± 7° by Benedict et al. (2006) and 282° ± 20° by Reffert & Quirrenbach (2011), both of whom are in disagreement regarding the debris disk position angle measurements, which are defined by the major axis. To confirm whether the stellar equatorial plane is aligned with the debris disk or with the planetary orbit, more MIRC-X observations are required to interferometrically measure both the star's position angle and its inclination, from which an improved understanding of the system's orientation and evolution could be derived in addition to a better stellar image for our spectroscopic analysis.

## 6.5. Conclusions and Future Work

We have shown that if available, a light-curve inversion image of the stellar surface can provide crucial information for disentangling the signature of stellar activity from that of planets. Obtaining a more accurate image will result in more accurate simulated RVs.

We recommend both RV and imaging observations should be of high cadence. Cabot et al. (2021) emphasized the importance of high-cadence RV observations for modeling stellar activity with GPs and improving the RV precision for detecting planets. Obtaining high-cadence observations for complementary imaging is also vital to having an accurate

image of the stellar surface from which RVs will be modeled. Here, we obtained two short observations of ϵ Eri with the MIRC-X beam combiner. More observations throughout a night would provide denser uv plane coverage, and observations across the stellar rotation would allow for the entire stellar surface to be imaged in a method analogous to Doppler imaging and light-curve inversion. The method developed here will work best when the imaging and RV observations are obtained contemporaneously to ensure that the same stellar surface evolution is being observed. As seen for ϵ Eri in Figure 3, the surface structures of main-sequence stars can change on timescales shorter than the stellar rotation.

While each high-cadence data set investigating the stellar activity and potential for planetary companions is valuable, combining the data sets as we describe here allows for the strengths of EPRV and stellar surface imaging to be harnessed in a way otherwise only accessible for the Sun. The unique value of solar studies is prior detailed knowledge of surface features coupled with EPRV measurements. We are working toward extending this to other stars by using high-cadence photometry to reconstruct the stellar surface and verify that current interferometric capabilities can detect starspots on a main-sequence star. Interferometric stellar images will provide unambiguous prior information about surface features to better understand the impact of photospheric activity on RVs.

We are grateful to A. E. Szymkowiak for his role in the development of EXPRES; J. Ennis and A. Labdon for their role in the development of MIRC-X; and P. Gatkine, C. Trujillo, L. H. Wasserman, and M. West for contributing observations. We thank the anonymous referee for thoughtful comments and insights that improved this paper. These results made use of the LDT at Lowell Observatory. Lowell is a private, nonprofit institution dedicated to astrophysical research and public appreciation of astronomy and operates the LDT in partnership with Boston University, the University of Maryland, the University of Toledo, Northern Arizona University, and Yale University. Lowell Observatory sits at the base of mountains sacred to tribes throughout the region. We honor their past, present, and future generations, who have lived here for millennia and will forever call this place home. Support for the design and construction of EXPRES is provided by the National Science Foundation (NSF) MRI-1429365 and ATI-1509436 and Yale University. We gratefully acknowledge support to carry out this research from NSF 2009528, NSF 1616086, NASA 17-XRP17 2-0064, the Heising-Simons Foundation, and an anonymous donor in the Yale alumni community. The CHARA Array is supported by the NSF under grant Nos. AST-1636624 and AST-1715788. Institutional support has been provided by the GSU College of Arts and Sciences and the GSU Office of the Vice President for Research and Economic Development. MIRC-X received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant No. 639889). This research has made use of the Jean-Marie Mariotti Center Aspro service.<fnref n="20" /> The APT photometric data were supported by NASA, NSF, Tennessee State University, and the State of Tennessee through its Centers of Excellence program. This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by NASA's Science Mission Directorate. R.M.R. acknowledges support from the Yale Center for

<footnote n="20">²⁰ Available at http://www.jmmc.fr/aspro.</footnote>
Astronomy and Astrophysics Prize Postdoctoral Fellowship, the Heising-Simons 51 Pegasi b Postdoctoral Fellowship, and NASA EPRV 80NSSC21K1034. J.M.B. was supported in part by NASA XRP 80NSSC21K0571. S.K. acknowledges support from an ERC Starting Grant (grant agreement 639889). S.G.J. acknowledges partial support from the NASA TESS GI grant 80NSSC21K0243.

Facilities: TESS, DCT, CHARA, Fairborn APTs.

## Appendix A
### Posterior Distribution of the GP Fit to APT Photometry

Posterior draws of the GP hyperparameters are shown in Figure 8, based on the fit to APT photometry. Median values are listed with uncertainties corresponding to the 16th and 84th percentiles. The MAP values of $P_{GP}$, $L$, and $C$ were used in the preconditioned RV fits.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Corner plot showing posterior distribution draws of GP hyperparameters (celerite quasiperiodic covariance kernel), after burn-in by the nested sampler. The 16%, 50%, and 84% quantiles are marked with vertical dashed lines in the marginalized histograms, which are also used to define the median values and uncertainties printed above each column.</figcaption></figure>
## Appendix B
Determining Starspot Location and Position Angle and Starspot Recovery Tests

In Section 4.2, we discussed the interferometric models of $\epsilon$ Eri. In order to show the robustness of the location of the starspot in the 2020 November 2 data, we performed 100 bootstraps where the visibility, closure phase, and triple amplitude data were chosen with replacement, as described in Section 4.2 for the limb-darkening parameter, $\alpha$. We plot the results in Figure 9 with each bootstrap represented by a gray circle, of the approximate size of the starspot used in the fit, and the best-fit spot location shown in red. These bootstrap locations are approximately consistent with the region around the minimum of the $\chi^2$ surface shown in Figure 10, described below.

As mentioned in Section 4.3, we determined that the data suggest the position angle of $\epsilon$ Eri is ~335° east of north. This was found by superimposing the LI surfaces on top of the interferometric models and rotating the LI surfaces such that the LI starspot and the interferometric model starspot were aligned at the same longitude. Because only one starspot was visible, there were two orientations that put the starspots at the same longitude. For 2020 November 2, we found that the best-fit interferometric-modeled and LI-reconstructed starspots were aligned in longitude if the LI reconstruction was rotated ~75° and ~335° east of north in the plane of the sky.

However, the best-fit model for 2020 November 5 and the ~98° rotation of the star based on the $P_{rot} = 11.4$ days between the observations are not consistent. In Figure 10, we show the reduced $\chi^2$ surfaces for both 2020 November 2 and 5. For 2020 November 2 there is a global minimum (represented by the black region within the reduced $\chi^2$ space), but for 2020 November 5, there are multiple local minima, and the lowest reduced $\chi^2$ value is again represented in black.

Because we expected the starspot to be on the limb from the LI reconstructions, we investigated the ability of our interferometric spot-fitting model algorithm to detect starspots of varying sizes for a star with the angular diameter and limb darkening of $\epsilon$ Eri with the uv plane coverage and telescope configuration of the 2020 November 5 observations.

We aimed to recover a starspot in different locations and with varying contributions to the overall brightness of the star (reductions in light from 0.1% to 6.4%). We assigned an angular diameter $\theta_{LD} = 2.153$ mas and a limb-darkening coefficient $\alpha = 0.27$, as above. We used a circular starspot that had a spot-to-photosphere brightness ratio of 0.61 and allowed its position and size to vary, but we required that the starspot had to stay within the star. In Table 5, we include the test and recovered spot parameters.

In Figures 11 and 12, we show that we were able to reasonably accurately recover the location of the starspot when

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Outline of the model stellar surface with the starspot location of each bootstrap indicated by a gray circle. Each circle is the size of the starspot used in the fit. The thicker red circle is the location of the best-fit starspot location. The stellar surface is oriented as in Figure 5.</figcaption></figure>

it reduced the brightness of the star by at least 0.4%. The reduced $\chi^2$ surfaces included in these figures are a combination between the visibilities, closure phases, and triple amplitudes with the closure phases weighted 10 times more than the other observables, as described above. For spots that reduce the stellar brightness by less than 0.4%, the best-fit solution is an artifact. However, when examining the reduced $\chi^2$ space of the fit, the actual spot location was found to be a local minimum. The starspot removing the smallest amount of light from the stellar surface explored here, 0.1%, is analogous to our estimations for the starspot on the limb.

Because the recovery tests indicate that the actual location of the starspot on 2020 November 5 may be recoverable as a local minimum, we considered these regions of the reduced $\chi^2$ space. Extrapolating from the two possible position angles found for just 2020 November 2, ~75° was rejected because the 2020 November 5 spot location was not in a local minimum, but the position angle of ~335° placed the starspot in a local minimum, on the lower right of the right side of Figure 10. As a result, we included the starspot in this orientation and the associated model data in Figure 5. More data are required, however, to confirm this orientation.

Because we cannot verify the position angle and because the LI reconstructions of $\epsilon$ Eri indicate the presence of another starspot that was out of view on both 2020 November 2 and 5, we did not attempt to inform the photometric surfaces with the interferometric spot location.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Left: Reduced $\chi^2$ surface for detecting a starspot in the 2020 November 2 interferometric data. The color gradient was chosen to particularly emphasize the location of the global and local minima. The minimum within the black region is the location of the model starspot shown in Figure 5. Right: Reduced $\chi^2$ surface for detecting a starspot in the 2020 November 5 interferometric data. The black region is the global minimum, which our recovery tests show is likely to be an artifact, but the minimum on the lower right of the plot is aligned with our estimate of the starspot location given the position angles possible for the 2020 November 2 data and the LI reconstructions. The spot at this location and its associated interferometric observations are shown in Figure 5.</figcaption></figure>

<table><caption>**Table 5**
Interferometric Starspot Recovery Test Results</caption><thead><tr><th>Test Identifier</th><th>Assigned Distance from Center</th><th>Assigned Position Angle (deg E of N)</th><th>Assigned Brightness Decrease (%)</th><th>Model Distance from Center</th><th>Model Position Angle (deg E of N)</th><th>Model Brightness Decrease (%)</th></tr></thead><tbody><tr><td>A0</td><td>0.75</td><td>45</td><td>0.10</td><td>0.74</td><td>61</td><td>0.13</td></tr><tr><td>A1</td><td>0.75</td><td>225</td><td>0.10</td><td>1.01</td><td>356</td><td>0.10</td></tr><tr><td>B0</td><td>0.75</td><td>90</td><td>0.20</td><td>0.75</td><td>65</td><td>0.35</td></tr><tr><td>B1</td><td>0.75</td><td>270</td><td>0.20</td><td>0.93</td><td>152</td><td>0.12</td></tr><tr><td>C0</td><td>0.75</td><td>135</td><td>0.40</td><td>0.79</td><td>133</td><td>0.62</td></tr><tr><td>C1</td><td>0.75</td><td>315</td><td>0.40</td><td>0.95</td><td>331</td><td>0.30</td></tr><tr><td>D0</td><td>0.75</td><td>180</td><td>0.80</td><td>0.83</td><td>173</td><td>0.93</td></tr><tr><td>D1</td><td>0.75</td><td>0</td><td>0.80</td><td>0.83</td><td>358</td><td>0.91</td></tr><tr><td>E0</td><td>0.75</td><td>225</td><td>1.60</td><td>0.74</td><td>223</td><td>1.50</td></tr><tr><td>E1</td><td>0.75</td><td>45</td><td>1.60</td><td>0.76</td><td>45</td><td>1.98</td></tr><tr><td>F0</td><td>0.75</td><td>270</td><td>3.20</td><td>0.75</td><td>270</td><td>3.33</td></tr><tr><td>F1</td><td>0.75</td><td>90</td><td>3.20</td><td>0.75</td><td>90</td><td>3.39</td></tr><tr><td>G0</td><td>0.75</td><td>315</td><td>6.40</td><td>0.58</td><td>313</td><td>7.19</td></tr><tr><td>G1</td><td>0.75</td><td>135</td><td>6.40</td><td>0.63</td><td>133</td><td>7.19</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>**Figure 11.** First and third columns: Best-fit model for the interferometric spot recovery test. The surface temperature ranges from a spot temperature of $T_{\text{spot}} = 4100$ K to a photospheric temperature of $T_{\text{phot}} = 5100$ K. Second and fourth columns: Reduced $\chi^2$ surface for the starspot size indicated by the best-fit model (lower values are white and higher values are dark blue). The red crosses × indicate where the center of the best-fit spot is located. The green crosses × indicate where the center of the spot is positioned. Each row is labeled on the left and corresponds to the appropriately labeled row of test and recovered spot parameters in Table 5.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 12.** First and third columns: Best-fit model for the interferometric spot recovery test. Second and fourth columns: Reduced $\chi^2$ surface for the starspot size indicated by the best-fit model (lower values are white and higher values are dark blue). The red crosses × indicate where the center of the best-fit spot is located. The green crosses × indicate where the center of the spot is positioned. Each row is labeled on the left and corresponds to the appropriately labeled row of test and recovered spot parameters in Table 5.</figcaption></figure>

ORCID iDs

Rachael M. Roettenbacher https://orcid.org/0000-0002-9288-3482

Samuel H. C. Cabot https://orcid.org/0000-0001-9749-6150

Debra A. Fischer https://orcid.org/0000-0003-2221-0861

John D. Monnier https://orcid.org/0000-0002-3380-3307

Gregory W. Henry https://orcid.org/0000-0003-4155-8513

Heidi Korhonen https://orcid.org/0000-0003-0529-1161

John M. Brewer https://orcid.org/0000-0002-9873-1471

Joe Llama https://orcid.org/0000-0003-4450-0368

Ryan R. Petersburg https://orcid.org/0000-0003-2168-0191

Lily L. Zhao https://orcid.org/0000-0002-3852-3590

Stefan Kraus https://orcid.org/0000-0001-6017-8773

Jean-Baptiste Le Bouquin https://orcid.org/0000-0002-0493-4674

Narsireddy Anugu https://orcid.org/0000-0002-2208-6541

Claire L. Davies https://orcid.org/0000-0001-9764-2357

Tyler Gardner https://orcid.org/0000-0002-3003-3183

Cyprien Lanthermann https://orcid.org/0000-0001-9745-5834

Gail Schaefer https://orcid.org/0000-0001-5415-9189

Benjamin Setterholm https://orcid.org/0000-0001-5980-0246

Catherine A. Clark https://orcid.org/0000-0002-2361-5812

Svetlana G. Jorstad https://orcid.org/0000-0001-6158-1708

Kyler Kuehn https://orcid.org/0000-0003-0120-0808

Stephen Levine https://orcid.org/0000-0002-1050-3539

References

Aigrain, S., Pont, F., & Zucker, S. 2012, MNRAS, 419, 3147

Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O'Neil, M. 2015, *ITPAM*, 38, 252

Anugu, N., Le Bouquin, J.-B., Monnier, J. D., et al. 2020, AJ, 160, 158

Baines, E. K., & Armstrong, J. T. 2012, ApJ, 744, 138

Barragán, O., Aigrain, S., Rajpaul, V. M., & Zicher, N. 2022, MNRAS, 509, 866

Benedict, G. F., McArthur, B. E., Gatewood, G., et al. 2006, AJ, 132, 2206

Berdyugina, S. V. 2005, LRSP, 2, 8

Blackman, R. T., Fischer, D. A., Jurgenson, C. A., et al. 2020, AJ, 159, 238

Boisse, I., Bonfils, X., & Santos, N. C. 2012, A&A, 545, A109

Booth, M., Dent, W. R. F., Jordán, A., et al. 2017, MNRAS, 469, 3200

Brewer, J. M., Fischer, D. A., Blackman, R. T., et al. 2020, AJ, 160, 67

Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125

Cabot, S. H. C., Roettenbacher, R. M., Henry, G. W., et al. 2021, AJ, 161, 26

Casasayas-Barris, N., Pallé, E., Yan, F., et al. 2019, A&A, 628, A9
Chavez-Dagostino, M., Bertone, E., Cruz-Saenz de Miera, F., et al. 2016, MNRAS, 462, 2285

Chelli, A., Duvert, G., Bourgès, L., et al. 2016, A&A, 589, A112

Claret, A. 2018, A&A, 618, A20

Claret, A., Diaz-Cordoves, J., & Gimenez, A. 1995, A&AS, 114, 247

Collier Cameron, A., Ford, E. B., Shahaf, S., et al. 2021, MNRAS, 505, 1699

Şenavcı, H. V., Kılıçoğlu, T., Işık, E., et al. 2021, MNRAS, 502, 3343

Davis, A. B., Cisewski, J., Dumusque, X., Fischer, D. A., & Ford, E. B. 2017, ApJ, 846, 59

di Folco, E., Absil, O., Augereau, J. C., et al. 2007, A&A, 475, 243

Ducati, J. R. 2002, yCat, 22, 37

Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132

Dumusque, X., Borsa, F., Damasso, M., et al. 2017, A&A, 598, A133

Dumusque, X., Cretignier, M., Sosnowska, D., et al. 2021, A&A, 648, A103

Faria, J. P., Adibekyan, V., Amazo-Gómez, E. M., et al. 2020, A&A, 635, A13

Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, OJAp, 2, 10

Fischer, D. A., Anglada-Escude, G., Arriagada, P., et al. 2016, PASP, 128, 066001

Fischer, D. A., Marcy, G. W., & Spronck, J. F. P. 2014, ApJS, 210, 5

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1

Giguere, M. J., Fischer, D. A., Zhang, C. X. Y., et al. 2016, ApJ, 824, 150

Gilbertson, C., Ford, E. B., Jones, D. E., & Stenning, D. C. 2020, ApJ, 905, 155

Greaves, J. S., Sibthorpe, B., Acke, B., et al. 2014, ApJL, 791, L11

Harmon, R. O., & Crews, L. J. 2000, AJ, 120, 3274

Hatzes, A. P., Cochran, W. D., McArthur, B., et al. 2000, ApJL, 544, L145

Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517

Haywood, R. D., Collier Cameron, A., Unruh, Y. C., et al. 2016, MNRAS, 457, 3637

Heitzmann, A., Marsden, S. C., Petit, P., et al. 2021, MNRAS, 505, 4989

Hempelmann, A., Mittag, M., Gonzalez-Perez, J. N., et al. 2016, A&A, 586, A14

Henry, G. W. 1999, PASP, 111, 845

Hestroffer, D. 1997, A&A, 327, 199

Holland, W. S., Matthews, B. C., Kennedy, G. M., et al. 2017, MNRAS, 470, 3606

Holzer, P. H., Cisewski-Kehe, J., Zhao, L., et al. 2021, AJ, 161, 272

Jurgenson, C., Fischer, D., McCracken, T., et al. 2016, Proc. SPIE, 9908, 99086T

Korhonen, H., Roettenbacher, R. M., Gu, S., et al. 2021, A&A, 646, A6

Kosiarek, M. R., & Crossfield, I. J. M. 2020, AJ, 159, 271

Lacour, S., Meimon, S., Thiébaut, E., et al. 2008, A&A, 485, 561

Lanza, A. F., Das Chagas, M. L., & De Medeiros, J. R. 2014, A&A, 564, A50

Levine, S. E., Bida, T. A., Chylek, T., et al. 2012, Proc. SPIE, 8444, 844419

Luger, R., Foreman-Mackey, D., & Hedges, C. 2021, AJ, 162, 124

Luo, T., Liang, Y., & IP, W.-H. 2019, AJ, 157, 238

MacGregor, M. A., Wilner, D. J., Andrews, S. M., Lestrade, J.-F., & Maddison, S. 2015, ApJ, 809, 47

Martinez, A. O., Baron, F. R., Monnier, J. D., Roettenbacher, R. M., & Parks, J. R. 2021, ApJ, 916, 60

Mawet, D., Hirsch, L., Lee, E. J., et al. 2019, AJ, 157, 33

Metcalfe, T. S., Buccino, A. P., Brown, B. P., et al. 2013, ApJL, 763, L26

Milbourne, T. W., Haywood, R. D., Phillips, D. F., et al. 2019, ApJ, 874, 107

Milne, E. A. 1921, MNRAS, 81, 361

Monnier, J. D., Che, X., Zhao, M., et al. 2012, ApJL, 761, L3

Parks, J. R., White, R. J., Baron, F., et al. 2021, ApJ, 913, 54

Petersburg, R. R., Ong, J. M. J., Zhao, L. L., et al. 2020, AJ, 159, 187

Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., & Roberts, S. 2015, MNRAS, 452, 2269

Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for Machine Learning (Cambridge, MA: MIT Press)

Reffert, S., & Quirrenbach, A. 2011, A&A, 527, A140

Rice, J. B., Wehlau, W. H., & Khokhlova, V. L. 1989, A&A, 208, 179

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, Proc. SPIE, 9143, 914320

Roettenbacher, R. M., Harmon, R. O., Vutisalchavakul, N., & Henry, G. W. 2011, AJ, 141, 138

Roettenbacher, R. M., Kane, S. R., Monnier, J. D., & Harmon, R. O. 2016, ApJ, 832, 207

Roettenbacher, R. M., Monnier, J. D., Fekel, F. C., et al. 2015, ApJ, 807, 23

Roettenbacher, R. M., Monnier, J. D., Harmon, R. O., Barclay, T., & Still, M. 2013, ApJ, 767, 60

Roettenbacher, R. M., Monnier, J. D., Korhonen, H., et al. 2017, ApJ, 849, 120

Roettenbacher, R. M., & Vida, K. 2018, ApJ, 868, 3

Savanov, I. S., & Strassmeier, K. G. 2008, AN, 329, 364

ten Brummelaar, T. A., McAlister, H. A., Ridgway, S. T., et al. 2005, ApJ, 628, 453

Vogt, S. S., Penrod, G. D., & Hatzes, A. P. 1987, ApJ, 321, 496

Zechmeister, M., Kürster, M., Endl, M., et al. 2013, A&A, 552, A78

Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12