# Characterization of K2-167 b and CALM, a new stellar activity mitigation
method

Zoë L. de Beurs ,1,2 Andrew Vanderburg,3 Erica Thygesen,4,† Joseph E.
Rodriguez,4 Xavier Dumusque, 5 Annelies Mortier,6 Luca Malavolta, 7,8 Lars A.
Buchhave,9 Christopher J. Shallue, 10 Sebastian Zieba, 11,12 Laura Kreidberg
,11 John H. Livingston, 13,14,15 R. D. Haywood,16,17 David W. Latham ,10
Mercedes López-Morales, 11 André M. Silva ,18,19
1Department of Earth, Atmospheric and Planetary Sciences, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA
2NSF Graduate Research Fellow & MIT Presidential Fellow
3 Department of Physics and Kavli Institute for Astrophysics and Space
Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
4 Center for Data Intensive and Time Domain Astronomy, Department of Physics
and Astronomy, Michigan State University, East Lansing, MI 48824, USA
5 Observatoire de Genève, Université de Genève, 51 chemin des Maillettes, 1290
Versoix, Switzerland
6 School of Physics & Astronomy, University of Birmingham, Edgbaston,
Birmingham, B15 2TT, UK
7 Dipartimento di Fisica e Astronomia “Galileo Galilei”, Universitá di Padova,
Vicolo del l’Osservatorio 3, I-35122 Padova, Italy
8 INAF - Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5,
Padova, 35122, Italy
9 DTU Space, National Space Institute, Technical University of Denmark,
Elektrovej 328, DK-2800 Kgs. Lyngby, Denmark
10 Center for Astrophysics $\mid$ Harvard & Smithsonian, 60 Garden St,
Cambridge, MA 02138, USA
11 Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg,
Germany.
12 Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333CA Leiden, The
Netherlands
13 Astrobiology Center, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
14 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan
15 Astronomical Science Program, Graduate University for Advanced Studies,
SOKENDAI, 2-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan
16 Astrophysics Group, University of Exeter, Exeter EX4 2QL, UK
17 STFC Ernest Rutherford Fellow
18 Departamento de Fisica e Astronomia, Faculdade de Ciencias, Universidade do
Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal
19 Instituto de Astrofisica e Ciencias do Espaco, Universidade do Porto, CAUP,
Rua das Estrelas, 4150-762 Porto, Portugal
† Quad Fellow Corresponding author email<EMAIL_ADDRESS>0000-0002-7564-6047
0000-0001-7246-5438 0000-0002-9165-6245 0000-0001-8812-0565
0000-0002-9332-2011 0000-0001-7254-4363 0000-0002-6492-2085
0000-0003-1605-5666 0000-0002-7585-9974 0000-0003-0562-6750
0000-0003-0514-1147 0000-0002-4881-3620 0000-0001-9140-3574
0000-0001-9911-7388 0000-0003-3204-8183 0000-0003-4920-738X

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We report precise radial velocity (RV) observations of HD 212657 (= K2-167), a
star shown by _K2_ to host a transiting sub-Neptune-sized planet in a 10 day
orbit. Using Transiting Exoplanet Survey Satellite (_TESS_) photometry, we
refined the planet parameters, especially the orbital period. We collected 74
precise RVs with the HARPS-N spectrograph between August 2015 and October
2016. Although this planet was first found to transit in 2015 and validated in
2018, excess RV scatter originally limited mass measurements. Here, we measure
a mass by taking advantage of reductions in scatter from updates to the
HARPS-N Data Reduction System (2.3.5) and our new activity mitigation method
called CCF Activity Linear Model (_CALM_), which uses activity-induced line
shape changes in the spectra without requiring timing information. Using the
_CALM_ framework, we performed a joint fit with RVs and transits using
EXOFASTv2 and find $M_{p}=$ $6.3_{-1.4}^{+1.4}$$\,M_{\oplus}$ and $R_{p}=$
$2.33^{+0.17}_{-0.15}$$\,R_{\oplus}$, which places K2-167 b at the upper edge
of the radius valley. We also find hints of a secondary companion at a $\sim$
22 day period, but confirmation requires additional RVs. Although
characterizing lower-mass planets like K2-167 b is often impeded by stellar
variability, these systems especially help probe the formation physics (i.e.
photoevaporation, core-powered mass loss) of the radius valley. In the future,
_CALM_ or similar techniques could be widely applied to FGK-type stars, help
characterize a population of exoplanets surrounding the radius valley, and
further our understanding of their formation.

###### keywords:

exoplanets – techniques: radial velocities – planets and satellites :
terrestrial planets – stars : activity – methods : statistical

††pubyear: 2015††pagerange: Characterization of K2-167 b and CALM, a new
stellar activity mitigation method–A.2

##  1 Introduction

The study of exoplanets has the ultimate goal of understanding the formation
and evolution of planets. In our own solar system, we have a class of dense
smaller rocky planets with iron cores and a class of large gas giants that
have extended H/He envelopes. However, the HARPS spectrograph (Mayor et al.,
2003) and Kepler mission (Borucki et al., 2010) found hundreds of planets
(Mayor et al., 2011; Batalha et al., 2013) that fall between these two
categories and with it brought many new questions: What is the composition of
these in-between planets? Does liquid water exist on their surface? Do they
have gaseous envelopes? If so, do their envelopes mainly consist of H/He? How
did this class of planets form and evolve? And why do we not see any examples
of them in our own solar system? To answer these questions, we need to
characterize these types of systems and constrain their possible compositions
through precise radius and mass measurements.

Precise radial velocity (RV) observations revealed that transiting planets
larger than 1.6 $\,R_{\oplus}$ are typically inconsistent with rocky
compositions (Rogers, 2015), so these intermediate planets are often sub-
divided into two classes: super-Earths ($1\,R_{\oplus}$ $\lesssim$
$R_{p}\lesssim 1.6\,R_{\oplus}$) and sub-Neptunes ($1.6\,R_{\oplus}\lesssim
R_{p}\lesssim 4\,R_{\oplus}$). A gap in the occurrence of planets with radii
between these two populations was predicted by several researchers based on
the effect of photoevaporation (Lopez & Fortney, 2013; Owen & Wu, 2013; Jin et
al., 2014). This gap is often described as the “radius valley” or “radius gap”
and describes a paucity of planets between 1.5 to 2 $\,R_{\oplus}$. This
bimodality in the occurence rates of small planets has been detected (Fulton
et al., 2017) and further refined and characterized using parallax and
asteroseismic measururements (Berger et al., 2018; Van Eylen et al., 2018;
Petigura, 2020).

Recently, Petigura et al. (2022) investigated the location of the radius gap
as a function of stellar mass, metallicity, age, and orbital period and found
that the radius gap location follows $R_{p}\propto P^{m}$ where $m=-0.10\pm
0.03$. Dattilo et al. (2023) also found that the location of the radius gap
varies depending on stellar type for FGK stars. However, the exact underlying
physics that explains the origins of the radius gap is still somewhat debated
(e.g. Ginzburg et al., 2018; Gupta & Schlichting, 2020; Lee et al., 2022) and
determining which mechanism is primarily responsible requires precise
characterization of a larger population of exoplanets (Rogers et al., 2021).
Precise radius and mass measurements of planets at the boundaries of the
radius valley can especially help refine our understanding of the physics by
constraining the bulk planetary composition and inferring either the presence
or absence of a thick volatile envelope.

Unfortunately, the difficulty of modeling stellar variability in RV
observations has limited the study of these important systems. In particular,
spurious RV shifts due to stellar activity signals can mask or mimic planetary
signals. In the past decade, stellar activity has been identified as the
limiting factor in RV precision across the HR diagram (National Academies of
Sciences, Engineering, and Medicine and others, 2018; Haywood et al., 2020),
but it is especially problematic for young stellar hosts, including those
where planets might still be undergoing photoevaporation. Finding new ways to
mitigate stellar activity is critical for our understanding of planetary
composition and evolution.

Stellar activity mitigation methods have traditionally focused on using
indicators of the level of stellar activity, which can either be focused on
measuring the level of magnetic activity on a star (Boisse et al., 2009;
Dumusque et al., 2011; Bonfils et al., 2007; Robertson et al., 2014) or
focused on detecting the presence (or lack thereof) of spots and plagues on
the surface (Figueira et al., 2013; Milbourne et al., 2021; Wise et al.,
2022). The activity indicators that correlate with magnetic activity include
spectral features such as the emission in the core of Ca II H&K lines (Saar et
al., 1998; Meunier & Lagrange, 2013; Saar & Fischer, 2000), the Hα line
(Bonfils et al., 2007; Robertson et al., 2014), and the Ca infrared triplet.
Recently, photospheric unsigned magnetic flux has been shown to be strongly
correlated with RV perturbations caused by solar surface activity (Haywood et
al., 2022) and new methods for deriving proxies for magnetic flux using least
squares deconvolution similarly show some of the strongest correlations among
activity indicators (Lienhard et al., 2023). Beyond these types of activity
indicators, properties in cross-correlation function (CCFs) have been used to
decorrelate stellar activity from RVs. CCFs are computed by cross-correlating
the spectra with a binary mask which contains delta functions at the location
of spectral lines. Conceptually, a CCF represents an average of all line
shapes in a star’s spectrum and is therefore sensitive to line shape changes
that persist in most stellar lines. The line shape changes introduced by
stellar activity that cause asymmetry in the CCF are commonly measured using a
CCF Bisector Inverse Span (e.g. Queloz et al., 2001; Lovis et al., 2011) or
the FWHM of the CCF (e.g. Queloz et al., 2009; Hébrard et al., 2010). Although
the bisector and FWHM have the advantage that they can easily be computed and
provide only two free parameters in RV models, much of the information present
in the CCF is lost by only using these metrics to measure line shape changes.
In addition, these traditional stellar activity indicators have been
insufficient to reach the sub $\rm m\,s^{-1}$ RV measurement precision
necessary to detect earth-twins.

There have been recent developments in using advanced flexible models like
Gaussian Process (GP) regression to model quasi-periodic stellar activity
signals and reveal planet signals (Haywood et al., 2014; Rajpaul et al., 2015;
Jones et al., 2017; Delisle et al., 2018; Barragán et al., 2019; Gilbertson et
al., 2020; Nicholson & Aigrain, 2022; Barragán et al., 2023; Tran et al.,
2023). However, even with the current state-of-the-art stellar variability
mitigation techniques (Zhao et al., 2022), our RV measurement precision must
still improve by an order of magnitude to detect the 10 $\rm cm\,s^{-1}$
signals induced by Earth-mass exoplanets in the habitable zones of Sun-like
stars. To work towards achieving this goal, we can first focus on precise mass
measurements of slightly more massive planets like super-Earths and sub-
Neptunes. This requires us to continue developing new methods to mitigate the
impact of stellar activity on exoplanet mass measurements.

Inspired by the successful implementation of a neural network stellar activity
mitigation technique on solar observations (de Beurs et al., 2022) and other
uses of machine and deep learning methods for modeling stellar activity
(Perger et al., 2023; Colwell et al., 2023), we designed a simplified method
called CCF Activity Linear Model (_CALM_), which can be applied to other stars
with smaller RV datasets observed at nighttime. Our method focuses on
exploiting the information content present in the CCFs and traces the line
shape changes introduced by stellar activity to separate Doppler reflex motion
from signals from the star. We have already successfully applied this method
to RVs for four stars observed with the EXPRES spectrograph and reduced the
scatter by about a factor of two for the more active stars (Zhao et al.,
2022). In this paper, we apply this same method to one star observed by
HARPS-N: K2-167.

K2-167 b is a planet located just above the radius valley and this transiting
sub-Neptune was first found by the _K2_ mission (Vanderburg et al., 2016; Mayo
et al., 2018). The planetary signal was re-detected (Ikwut-Ukwa et al., 2020)
by the Transiting Exoplanet Survey Satellite (_TESS_ ; Ricker et al., 2015),
and later again in its first extened mission (Thygesen et al., 2023). HARPS-N
observations of K2-167, a relatively active, yet slowly rotating F7-type star,
began in late 2015, but in the intervening years, stellar activity had
prevented a confident detection of the planet’s mass. Recently, Bonomo et al.
(2023) analyzed the K2-167 system as part of a larger sample and reported a
mass of $6.5_{-1.5}^{+1.6}$ $\,M_{\oplus}$. Here, we perform a more in-depth
analysis of the RV time series than this previous work (which was included as
part of a catalog paper). Using our _CALM_ framework, we performed a joint fit
of the transits and RVs, constraining the radius and mass, and enabling future
studies of its composition which may provide insight into its formation.

Our paper is organized as follows. In Section 2, we describe the the
spectroscopic observations from HARPS-N and the photometric observations from
_K2_ and _TESS_. In Section 3, we describe our development of the _CALM_
stellar activity mitigation method, how we prevent overfitting on our data,
and introduce a new metric to measure RV scatter. In Section 4, we describe
the transit, spectroscopic, and radial velocity analysis and resulting system
parameters of K2-167. In Section 5, we discuss the implications of these
results and we conclude in Section 6.

##  2 Data

###  2.1 HARPS-N Spectroscopy

K2-167 was observed with the HARPS-N Spectrograph on the 3.58m Telescopio
Nazionale Galileo (TNG). Designed specifically to yield precise radial
velocity measurements, HARPS-N is a vacuum-enclosed cross-dispersed echelle
spectrograph with temperature and pressure stabilization (Cosentino et al.,
2012). We collected 76 precise radial velocity observations of K2-167 (Table
A.1) with typical integration times of 15 minutes111All observations had 15
minute integration times except for one 20 minute exposure on 26 September
2016 and two 30-minute observations on 28 November 2015 and 9 July 2016
respectively. with the HARPS-N spectrograph between August 2015 and October
2016. For the first season of observations (20 August 2015 to 27 December
2015), a ThAr lamp served as a calibration source. For the second season (20
June 2016 to 12 October 2016), a stabilized Fabry Perot was used as the
calibration source. Most of the time, we obtained only one observation per
night, but between July and November 2016, we often obtained two spectra per
night, separated by a few hours, in order to average radial velocity
variations due to granulation.

The 15 minute exposures of K2-167 yielded typical formal measurement
uncertainties of 1.4 $\rm m\,s^{-1}$. For a bright star like K2-167
($m_{v}=8.24$), we expect this uncertainty to be a combination of instrumental
systematics and the photon-limited uncertainty. We performed our analysis on
an older version of the HARPS-N pipeline (DRS 3.7) and the latest version of
the pipeline (DRS 2.3.5) (Dumusque et al., 2021) and found significantly
larger scatter for the older DRS pipeline, which we discuss extensively in
Section 3. There were numerous improvements between this older and new DRS to
correct for systematics and instrumental effects. In addition, in the new DRS,
different stellar masks were used that focus on spectral lines that are less
sensitive to stellar activity and to more closely match the stellar type of
the targets. In the case of K2-167(F7-type star), DRS 3.7 uses a K5-type
binary mask and the DRS 2.3.5 uses a G9-type binary mask. For both of these
pipeline versions, the data reduction was carried out by cross-correlating the
observations with one of these binary masks. For the DRS 3.7, this produces a
161 element cross-correlation function (CCF), whereas for the newer 2.3.5 DRS,
the sampling was reduced and the pipelines yields a 49 element CCF. These
resulting CCFs are used as an input representation for our _CALM_ method to
track the activity-induced line shape changes in the spectra. The DRS measures
the radial velocities of each observation by fitting the CCF with a Gaussian
function. However, a Gaussian is simple and symmetric and therefore unable to
model the small CCF shape changes introduced by stellar activity. Thus, these
radial velocities will include contributions from both Keplerian shifts and
stellar activity signals. The full HARPS-N RV time series from the new DRS is
shown in Figure 1.

Figure 1: RVs from the HARPS-N old DRS (3.7) in blue and the new DRS (2.3.5)
in orange. a) The RVs are plotted over time. There appears to be a slight
upward trend toward the end of the time series. b) Histogram of the old and
new DRS. The new DRS has less scatter than the old DRS.

###  2.2 Photometry

####  2.2.1 K2

K2-167 b was observed by _K2_ , Kepler’s extended mission which achieved
similar precision to the original Kepler mission after applying systematic
corrections (Vanderburg et al., 2016). K2-167 was observed during K2 Campaign
3 from 2014 November 17 to 2015 January 23 in _K2_ ’s long cadence mode (29.4
minute exposure times). We used the light curve produced by Mayo et al.
(2018). In brief, they extracted the light curve for K2-167 from the target
pixel files, which were accessed through the Mikulski Archive for Space
Telescopes (MAST)222mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html. The
light curves were reprocessed to simultaneously fit the known planet transit,
remove mechanical K2 systematics, and fit stellar variability by using the
methods described in Vanderburg & Johnson (2014a) and Vanderburg et al.
(2016). Due to the unusual brightness of the host star, Mayo et al. (2018)
used larger photometric apertures than their standard analysis for this star.
These apertures did not require dilution corrections because K2-167 is
sufficiently isolated from other stars or background objects (Ligi et al.,
2018). To remove high-frequency variability, we divided the simultaneous-fit
light curve by the best-fit spline component. The _K2_ phase-folded transit
light curves are shown in the top panel of Figure 2. For computational
efficiency, we trimmed the light curve away from transit for our global fit,
and only included a baseline of one transit duration on either side of the
full transit (See Section 4.2).

Figure 2: _K2_ (top) and _TESS_ (bottom) light curves phase-folded on the
period and centered at mid-transit from _K2_ Campaign 3 and _TESS_ Sectors 2,
28, 42 . The purple lines shows the best fitting transit models. Figure 3:
Schematic of our _CALM_ pipeline. We use the changing shape of spectral
features (measured by the change in the cross-correlation function or CCF) to
predict and remove stellar activity signals. We use a linear regression model
with 5-25 indexes from the CCF residuals used as inputs. We also
simultaneously fit Keplerians along with the activity corrections to predict
the overall RV signal. The data shown in (b) and (d) are $\Delta$CCFs and the
RVs of K2-167 respectively for the old DRS (3.7).

####  2.2.2 TESS

K2-167 was pre-selected for two-minute cadence observations by the _TESS_
mission (Ricker et al., 2015) and observed by Camera 1 during _TESS_ Sector 2
from 2018 August 22 to September 20, Sector 28 from 2020 July 31 to 2020
August 25, and Sector 42, from 2021 August 20 to 2021 September 16. The _TESS_
light curve was generated by NASA’s Science Processing Operations Center
(SPOC) pipeline using the Lightkurve package (Lightkurve Collaboration et al.,
2018). SPOC received raw data and processed the images to remove systematic
erros and extract photometry (Jenkins et al., 2016). Using the SPOC Transiting
Planet Search (Jenkins, 2002), the light curves are searched for transit
events. The Pre-search Data Conditioned Simple Aperture Photometry (PDCSAP)
flux uses the optimal _TESS_ aperture to extract the flux, corrects for
systematics using PDC, and produces the light curves (Stumpe et al., 2012;
Stumpe et al., 2014; Smith et al., 2012). From these light curves we remove
low frequency trends that correspond to astrophysical variability and remove
remaining systematics based on the out-of-transit photometry using
keplerspline333https://github.com/avanderburg/keplersplinev2 (Vanderburg &
Johnson, 2014b; Shallue & Vanderburg, 2018). The _TESS_ phase-folded transit
light curves are shown in the bottom panel of Figure 2.

####  2.2.3 Spitzer

To perform follow up observations on the K2 transits, one transit of K2-167 b
was observed with the Spitzer InfraRed Array Camera (IRAC; Fazio et al., 2004)
on March 1st 2016 as part of a larger K2 follow-up program (Program ID: 11026;
PI Michael Werner). The observations lasted for 10.16 hours and were taken
with 0.4 second exposures using channel 2 on IRAC (2-4 $\mu$m). The
observations are further described in detail in Duck et al. (2021) who found
that the complete egress was missing due to the uncertainty in scheduling of
the observations. Nonetheless, we reduced the data using BiLinearly
Interpolated Subpixel Sensitivity (BLISS) mapping (Stevenson et al., 2012) and
Pixel Level Decorrelation (Deming et al., 2015; Garhart et al., 2020) and saw
a feature consistent with a partial transit, but with low enough SNR that it
was not constraining to our model. We therefore exclude the Spizter
observations from our final analysis.

##  3 RV Analysis

We performed an extensive series of analyses to determine the mass of K2-167
b. The RV time series is shown in Figure 1 where there appears to be a slight
long-term upward trend in the RVs that may be the result of long-term activity
cycles. For our activity analyses, we first performed a simple RV analysis
where we fitted a Keplerian signal and did not model the stellar variability.
We ran Markov Chain Monte Carlo (MCMC) to sample the parameter space and model
the Keplerian using radvel’s kepler.rv_drive function444We note that we
accounted for the way that radvel used $\omega_{p}$ rather than $\omega_{*}$
in its implementation (Householder & Weiss, 2022) and only report $\omega_{*}$
in this paper. (Fulton et al., 2018). We use the differential evolution MCMC
package edmcmc555https://github.com/avanderburg/edmcmc (Vanderburg, 2021),
which is an implementation of MCMC that incorporates differential evolution
(DE; Ter Braak, 2006) to solve the problem of choosing an appropriate scale
and orientation for the jumping distribution, enabling a significantly faster
fitting process. We computed 5,000 MCMC chains, discard the first 25% as burn-
in, and required that our parameters converged using the Gelman-Rubin
statistic (Gelman & Rubin, 1992), where our Gelman-Rubin threshold < 1.01. In
this analysis, we observed that there is a large difference in RV
scatter666Here we report the remaining RV scatter after removing our modeled
Keplerian signal. In this way, we report only the unmodeled variation in the
RVs as further detailed in Section 3.5. between the two DRS versions (3.95
$\rm m\,s^{-1}$ and 3.01 $\rm m\,s^{-1}$ for the old and new DRS,
respectively) and that the uncertainties on the estimated semi-amplitude (K)
for the Keplerian signal are larger for the old DRS ($2.0_{-0.7}^{+0.7}$ $\rm
m\,s^{-1}$) than the new DRS ($2.0_{-0.5}^{+0.5}$ $\rm m\,s^{-1}$). The
expected scatter due to instrument systematics (photon noise, drift noise, and
wavelenegth calibration errors) for the both DRS versions is $\sim$ 1 $\rm
m\,s^{-1}$, which is significantly less than the scatter that we measured. For
this chromospherically active star K2-167 ($\textrm{log
R}_{\textrm{HK}}^{{}^{\prime}}=-4.66$), stellar variability could be causing a
part of this excess scatter (and the relatively large uncertainty in K). Thus,
we sought to develop new methods that could model this stellar variability by
tracing shape changes in the CCF. For each of our RV analyses, we performed an
independent analysis using only the 3.7 DRS and an independent analysis using
only the 2.3.5 DRS observations. We included both DRS versions in our RV
analyses since we found interesting differences in the effectiveness of our
activity mitigation methods for the two versions as discussed extensively in
Section 3.6.

###  3.1 _CALM_ Method

In order to account for the stellar variability we observe in the RV
measurements of K2-167, we developed a simplified machine learning model
_CALM_ that allows us to separate the stellar activity from the Doppler reflex
motion. This simplified method is inspired by the previous implementation of a
neural network-based stellar activity mitigation method for solar observations
using line shape changes to predict stellar activity signals (de Beurs et al.,
2022). To extend this approach to extrasolar observations, a new simplified
method was required for two primary reasons:

  1. 1.

RV datasets for other stars rarely have as many observations as the solar
dataset from the HARPS-N Solar Telescope used in de Beurs et al. (2022), which
consisted of $\sim$600 days of observations.

  2. 2.

We lack a firm “ground-truth” for the stellar activity signals of other stars.
For the Sun, we can easily remove the RV contributions of the known planets in
the Solar System, but this is not as straighforward for other stars. This is
because there may always be unmodeled planetary signals that fall below the
detection threshold. These unmodeled planetary signals may add noise to our
estimates of the stellar variability signal of a given star.

To overcome these challenges, we made two modifications to our method. First,
we greatly simplified the machine learning models for use on stars other than
the Sun. Instead of using complex and flexible models like our convolutional
neural networks, we use a linear regression model with a significantly smaller
input representation. In particular, the input representation consists of 5-25
CCF activity indicators777We tested a range of input representation sizes as
discussed in Section 3.2.2 and ultimately find that 5 CCF activity indicators
to balances the risk of overfitting and underfitting to the data for the 2.3.5
DRS. For the 3.7 DRS, we find that 25 CCF indexes balances these concerns.
instead of the entire CCF array, which contains 49 elements for the 2.3.5 DRS
and 161 elements for the 3.7 DRS. Essentially, we have sacrificed some
fidelity in our stellar activity corrections in exchange for immediate
applicability to a much larger sample of stars. Eventually, it may be possible
to train a more complex model on the entire ensemble of stars observed, but
our simplified approach is a first step towards this goal.

The second modification we made to our technique for use on stars other than
the Sun involves fitting for planet signals simultaneously. We cannot a priori
remove all planetary signals from the dataset (since the amplitude and period
of such signals are unknown and/or they may be hidden by stellar activity).
Thus, we must fit any remaining planetary signals simultaneously with the
stellar activity correction. Although this increases the computational
difficulty of our method (because we must test many different possible
planetary signals to discover new planets), the increase in computational
difficulty is offset by the simplicity of our linear regression model. An
overview of our simplified approach is shown in Figure 3.

We can describe this simplified approach as simultaneously fitting for $n$
Keplerian signal(s) and a linear regression model with m CCF activity
indicators for stellar activity. The number of CCF activity indicators m is
determined by finding a balance between optimizing goodness-of-fit and
preventing overfitting. In equation form, our model is given as:

$\begin{split}RV={\color[rgb]{0,0,1}w_{1}\cdot CCF_{1}+w_{2}\cdot
CCF_{2}+...+w_{m}\cdot CCF_{m}}\\\ +b_{1}\cdot Kepler_{1}+b_{2}\cdot
Kepler_{2}+...+b_{n}\cdot Kepler_{n}\\\ ={\color[rgb]{0,0,1}RV_{\rm
activity}}+RV_{\rm Keplerians}\end{split}$ (1)

where $CCF_{m}$ is the value of the CCF at index $m$, $w_{m}$ is a weight
parameter for $CCF_{m}$, $Kepler_{n}$ is the $n_{th}$ fitted Keplerian, and
$b_{n}$ is this Keplerian’s weight parameter. This Keplerian weight parameter
$b_{n}$ is equivalent to a planet’s semi-amplitude ($K$), but parametrizing
the model in this way allows us to do a linear least squares regression fit
where we can simultaneously model planetary and stellar variability
contributions to the overall RV signal.

This section describes our new _CALM_ method in detail. Section 3.2 describes
the pre-processing we must apply to run our _CALM_ method; in particular, we
first must prepare the CCF input representation for our model by centering the
data and choosing the most information-rich subset of the CCF to feed into the
model. Before running our models, we check several diagnostics to prevent
overfitting, which are described in Section 3.4. Next, we sample the parameter
space using a differential evolution-based MCMC method described in Section
3.3. Finally, we evaluate our models and compare them with several traditional
activity indicators using a new scatter metric described in Section 3.5 and
describe our RV analysis results in Section 3.6.

Figure 4: Computed residual CCFs ($\Delta$CCFs) for both the older (3.7) and
newer DRS (2.3.5). $\Delta$CCFs are computed by subtracting the mean CCF and
highlight differences in features between CCFs. The model predicted stellar
activity contribution to the RVs is indicated by its color (red = redshifted,
blue = blueshifted).

###  3.2 Preparing the CCF Input Representation

####  3.2.1 Computing and Centering our $\Delta$CCFs

For our _CALM_ model, we have to pre-process our data products into a uniform
format which makes capturing and learning features in our data easier. We
design the input representation such that the model becomes sensitive to shape
changes due to stellar activity, not translational shifts caused by true
Doppler shifts. Simply put, we want the _CALM_ model to find the difference
between a Gaussian fit to the CCF and the true center-of-mass velocity shift.
In this way, our model will be able to disentangle stellar activity signals
from true Keplerian shifts.

To frame the problem in this way, we shift every CCF such that the velocity
predicted by the Gaussian fit to the CCF is at the median velocity across all
observations. We note that we took precautions to prevent improperly centering
the CCFs in a way that planetary reflex motion remains in the input data. This
could cause overfitting and we further describe our precautionary steps in
Section 3.4. Next, we subtract a median CCF from each CCF to compute the
residual CCF888We note that subtracting this median is purely for
visualization purposes and not required for the regression algorithm.
($\Delta$CCF), which are visualized in Figure 4, where they are color-coded
based on how blueshifted or redshifted the RV signals are. There exists a
strong signal in shape changes based on these RV signals for the old DRS and a
hint of a dependence for the new DRS. This difference may be explained by the
new DRS using spectral lines less sensitive to stellar variability and the
difference in stellar type mask used. Overall, we want  our _CALM_ model to
learn how these shape changes relate to the stellar activity contributions to
RVs so that we can model and mitigate the stellar variability using the
$\Delta$CCFs.

After computing the $\Delta$CCFs, we normalize the observations to ensure the
scale of variations across each input parameter is approximately equal. A
subset of indexes across these normalized $\Delta$CCFs are then fed into the
_CALM_ model . We choose to only include a subset of indexes (5 to 25999The
exact number is determined by optimizing the fit while preventing overfitting
as described in Section 3.4 indexes) rather than the entire normalized
$\Delta$CCFs (49 indexes or 160 indexes for DRS 2.3.5 and DRS 3.7
respectively) because we only have $\sim$80 RV observations and the entire
normalized $\Delta$CCFs would provide too many free parameters. We describe
how we select which CCF indexes to use in Section 3.2.2.

####  3.2.2 Selecting the most significant indexes of the $\Delta$CCF

To ensure that our number of free parameters is significantly less than our
number of observations (74), we include only a subset of the CCF rather than
the entire 49-element array of the CCF 101010the entire array of the CCF
contains 161 values for the old DRS (3.7). This was adjusted to 49 in the new
DRS (2.3.5) to correct for oversampling of the CCF.. We choose the subset of
the CCF by checking which indexes provide the most statistically significant
contribution in Equation 1. We check for statistical significance by sampling
the parameter space with Differential Evolution MCMC (described in Section
3.3) and estimating uncertainties on the CCF weights
($w_{1},w_{2},...,w_{n}$). Commonly, when including $\geq 30$ CCF indexes, the
information at these indexes has degeneracies which dilutes the significance,
resulting in the highest significance indexes only having $1-2\sigma$
significance. This means that for these indexes the value of the weights for
$CCF_{i}$ is one to two times its corresponding uncertainty. We then choose
the 5-25 most statistically significant CCF indexes for our final model. This
exact number depends on the DRS version as explained in detail in Section
3.6.1. When only including those “significant indexes” in our model, the
significance increases to $3-5\sigma$. Generally, we find that the estimated
semi-amplitude (K) of the Keplerian signal of K2-167 b is consistent for any
of the number of CCF indexes tested (Figure 5).

Figure 5: Estimated semi-amplitude for various model configurations. (a)
Estimated semi-amplitude (K) as function of number CCF indexes for both the
old DRS (grey) and new DRS (blue). As the number of CCF indexes included
increases, there a slight positive trend in K. (b) Estimated K for K2-167 b in
a two-planet fit model as described in Section 3.7. Generally, the estimates
of K agree within the errors for any of these model configurations.

###  3.3 Running our RV models with Differential Evolution MCMC

Once we have chosen our CCF indexes and the number of Keplerians to include in
our model, we sample our parameter space using differential evolution MCMC
implemented in edmcmc111111https://github.com/avanderburg/edmcmc (Vanderburg,
2021). We implemented Gaussian priors to the $P$ and
$t_{{\color[rgb]{0,0,0}c}}$ parameters based on the EXOFASTv2 transit fit
results (Table 1), which are discussed in Section 4.2. We note that without
including priors on $P$ and $t_{c}$, we recover the period and time of
conjunction in the RVs. However, we include the priors in our analysis since
they allow us to get the tightest constraints on the mass of K2-167 b. We
model the Keplerian signal of K2-167 b by using radvel’s kepler.rv_drive
function (Fulton et al., 2018). After performing the MCMC fit, we determine
convergence of the parameters by using the Gelman-Rubin statistic (Gelman &
Rubin, 1992), where our Gelman-Rubin threshold < 1.01. We computed 5,000 MCMC
chains and discard the first 25% as burn-in before producing our posterior
samples. We compute the median and 68.3% confidence intervals to derive our
final parameter values. The results of our DE-MCMC analysis of the RVs are
described in Section 3.6.

Figure 6: Overfitting diagnostics for stellar activity RV analysis using the
new DRS (2.3.5) HARPS-N $\Delta$CCFs. The three columns in this Figure show
the $\Delta$CCFs, the weight parameters corresponding to the $\Delta$CCFs, and
the phase-folded RVs for three scenarios respectively. For each of the Figures
in the third column (c, f, i), we provide a legend which includes the scatter
of the radial velocities from the HARPS-N pipeline (raw rvs, std) in $\rm
m\,s^{-1}$, the scatter of the activity corrected radial velocities (corr rvs,
std) in $\rm m\,s^{-1}$, and the predicted semi-amplitude of the planet
(planet preds, K) in $\rm m\,s^{-1}$. The first two rows corresponds to two
different overfitting concerns and the last row demonstrates a case where both
of these concerns are addressed. In the first row (a, b, c), we plot the
$\Delta$CCFs that have not been shifted to be centered at the median velocity.
Without centering the $\Delta$CCFs, the algorithm will be able to access
translational shift (i.e. Doppler shifts) information and attribute potential
planet signals to stellar activity signals, significantly attenuating their
amplitude. In the second row(d, e, f), we show a case where the $\Delta$CCFs
are shifted but the number of indexes is too large. This results in
overfitting as seen in the weights plot for the second row. In the third
row(g, h, i), we have only fed in the $\Delta$CCF indexes that are considered
significant (as described in Section 3.2.2) and use the properly shifted
$\Delta$CCFs.

###  3.4 Overfitting Diagnostics

Overfitting is a significant problem in radial velocity analyses and the
state-of the art methods that are widely used to model and predict stellar
variability often do not agree within error. The severity of this problem is
particularly evident in recent work by Zhao et al. (2022) where 11 teams were
invited to use their 22 different methods on ultra precise EXPRES radial
velocities of four stars that encompass various levels of activity. There is a
concerning lack of agreement between these methods across these targets, which
perhaps can be attributed to some of the methods overfitting to noise in the
observations. Recently, (Blunt et al., 2023) demonstrated that Gaussian
process regression, one of the most common and state-of-the-art methods of
modeling stellar variability, can be susceptible to overfitting. We wanted to
prevent overfitting and thus set out to develop various methods to mitigate
the risk of overfitting.

We developed several diagnostic plots as illustrated in Figure 6. We developed
these diagnostics to prevent two common types of overfitting concerns:

  1. 1.

Not properly centering CCFs such that planetary reflex motion remains in the
input data.

  2. 2.

Using too many CCFs indexes such that there are too many free parameters
compared to number of observations

Figure 7: Phase-folded RVs before and after applying our activity correction
model. From the first to second column, we compare the phase-folded RVs for
the old (3.7) and new (2.3.5) DRS. (a) Before applying our linear regression
activity correction, the RV scatter is $\sigma_{\rm raw}$ = $4.02$ $\rm
m\,s^{-1}$ and we detect the semi-amplitude of the planet with 3.0$\sigma$
confidence in the old DRS (3.7) RVS. (b) When we simultaneously fit our _CALM_
activity model based on shape changes in the CCF with a Keplerian, the
detection significance increases to 4.3$\sigma$. c) For the new DRS, the raw
scatter is much lower ($\sigma_{\rm raw}$ = $3.01$ $\rm m\,s^{-1}$) and
without applying an activity model, we get a 4.3$\sigma$ detection. d) With an
activity model, the raw scatter decreases to $\sigma_{\rm raw}$ = $2.41$ $\rm
m\,s^{-1}$and the detection significance remains the same at 4.3$\sigma$. For
a description of the scatter metrics, $\sigma_{\rm raw}$ and $\sigma_{\rm
us}$, see Section 3.5.

To address the first concern, we plot the $\Delta$CCFs in the first column of
Figure 6 to check whether all translation shifts (i.e. planetary reflex
motion) were removed and the CCFs were properly centered. To illustrate how
this problem appears, we intentionally did not center the $\Delta$CCFs in
Figure 6a. We thus observe a sinusoidal-like pattern for the red $\Delta$CCFs
which correspond to more redshifted RVs and the opposite pattern for more blue
$\Delta$CCFs. This is a tell-tale sign of a translation shift still remaining
within the $\Delta$CCFs. We can understand this shape by reminding ourselves
of the definition of a derivative where we define the median CCF as $f(x)$ and
a redshifted CCF as $f(x-h)$ and thus can write

$f^{\prime}(x)=\lim_{h\to 0}\frac{f(x-h)-f(x)}{h}$ (2)

Thus, the shape expected from the difference between the redshifted CCF and
the median CCF is simply the derivative of a Gaussian. This is illustrated in
Figure A.1 where we simulated a few Gaussian curves and added white Gaussian
noise. We show the expected residuals for improperly centered and properly
centered CCFs in Figures A.1a,b and A.1c,d respectively. When $\Delta$CCFs are
properly shifted to the center, this pattern disappears as seen in both the
simulated example in Figure A.1c,d and in the real data in Figure 6d,g.

To address the second concern, we created diagnostic plots that include the
weights for the $\Delta$CCFs as seen in the second column of Figure 6. When
too many CCF weights are included, we tend to see structure in these weights
plots that is caused by overfitting as one weight is compensating for its
neighboring weight. We illustrate this in the second row of Figure 6 where we
show a case where the $\Delta$CCFs are shifted properly but the number of
indexes is too large. Although this result in quite a low RMS of 2.3 $\rm
m\,s^{-1}$, the large number of free parameters results in overfitting on the
noise in the data. We note that even when this overfitting occurs, the
planetary signal remains preserved with a semi-amplitude of K $\sim$ 2.2 $\rm
m\,s^{-1}$.

Finally, we show an example in the third row of Figure 6 where both the
$\Delta$CCFs are properly centered and the number of indexes (5) is small
enough to not result in overfitting. Both the derivative-of-a-Gaussian pattern
in the $\Delta$CCFs and the structure in the weights disappear.

###  3.5 Introducing a new performance metric: unexplained scatter

To measure the original scatter in the data before applying any activity
correction, we compute the raw scatter, $\sigma_{\rm raw}$, which is defined
as the sample mean standard deviation of the RVs minus the expected
measurement error such that:

$\sigma_{\rm raw}=\sqrt{SD(RV_{\rm
HARPS}){\color[rgb]{0,0,0}{}^{2}}-\sigma_{\rm{\color[rgb]{0,0,0}measurement}}^{2}}$
(3)

where $\sigma_{\rm{\color[rgb]{0,0,0}measurement}}$ is the median estimated RV
errors from the DRS which includes photon noise, drift noise, and errors in
the wavelength calibration (Dumusque et al., 2021). Defining the scatter in
this way assumes that we can approximate the source of RV variations as
Gaussian noise and that they can be approximated as independent scatter. Both
of these assumptions are not true but this provides a sufficient approximation
of the noise. We then evaluate the reduction in RV scatter using two scatter
metrics: the activity corrected scatter ($\sigma_{\rm ac}$) and the
unexplained scatter ($\sigma_{\rm us}$). We define $\sigma_{\rm ac}$ as

$\sigma_{\rm ac}=\sqrt{(SD(RV_{\rm HARPS}-RV_{\rm
SA}))^{2}-\sigma_{\rm{\color[rgb]{0,0,0}measurement}}^{2}}$ (4)

where SD is the sample standard deviation, $RV_{\rm HARPS}$ are the raw radial
velocities from the HARPS DRS pipeline, $RV_{\rm SA}$ are the stellar activity
contributions to the radial velocities as predicted by the given activity
model. This type of scatter metric 121212Sometimes the convention is to not
subtract the measurement errors from the instrument ($\sigma_{\rm
measurement}$) in calculating this scatter metric. is sometimes used to
evaluate and compare stellar activity models, but this metric does not account
for the scatter that can be caused by planet signal(s). Focusing on just
reducing this scatter metric could unintentionally incentivize removing planet
signals rather than preserving them with these activity correction methods. To
partially address this concern, we introduce a metric that subtracts modeled
planet signals so we primarily evaluate the activity model correction method
on unexplained scatter rather than any scatter. We define the unexplained
scatter ($\sigma_{\rm us}$) as

$\sigma_{\rm us}=\sqrt{(SD(RV_{\rm HARPS}-RV_{\rm SA}-RV_{\rm
planet}))^{2}-\sigma_{\rm{\color[rgb]{0,0,0}measurement}}^{2}}$ (5)

where $RV_{\rm planet}$ is the contribution of modeled planets to the RV
measurements. We note that there can evidently still remain un-modeled planet
signals so this scatter metric does not fully address the concern of
incentivizing removal of planet signals to reduce RMS. However, this can be a
step towards finding better metrics for evaluating our RV models.

Figure 8: Old DRS (3.7) Periodograms as a function of number of CCF indexes
used in the stellar activity correction model. HARPS-N Raw (a) and Corrected
(b) RVs in Fourier space for 5 CCF indexes. In each panel, the 1.0% and 10.0%
false alarm probabilities (FAP) computed using the bootstrap method are
indicated with green and black dotted lines, respectively. The long-term
activity signals in panel (a) decrease in magnitude after applying the
activity correction in panel (b). The suspected rotation period and half the
rotation period are indicated in grey. The planet period is indicated in
yellow. As we increase the number of included CCF indexes to N=15 (c) and N=25
(d), the peak corresponding to the star’s rotation period ($P_{\rm rot}$)
decreases significantly in magnitude after applying the stellar activity
corrections and a planet signal emerges at $9.978543$ days. Figure 9: New DRS
(2.3.5) Periodograms for 5 CCF indexes. HARPS-N Raw (a) and Corrected (b) RVs
in Fourier space for 5 CCF indexes. The 1.0% and 10.0% false alarm
probabilities (FAP) computed using the bootstrap method are indicated with
green and black dotted lines, respectively. The long-term activity signals in
panel (a) decrease in magnitude after applying the activity correction in
panel (b). The suspected rotation period and half the rotation period are
indicated in grey. The planet period is indicated in yellow. From (a) to (b),
the peak corresponding to a planet signal at $9.978543$ days increases
slightly in magnitude.

###  3.6 MCMC Results

For both the new and old DRS, the MCMC successfully converged to a solution
for the stellar activity model and the planet parameters.

For the old DRS (3.7), we found that _CALM_ stellar activity mitigation method
yields a significant improvement in both the RV scatter (from $\sigma_{\rm
raw}$ = $4.02$ $\rm m\,s^{-1}$to $\sigma_{\rm us}$ = $2.04$$\rm m\,s^{-1}$)
and the precision with which we recover the mass of the transiting planet
(improving a 3.0$\sigma$ significance detection into a 4.3$\sigma$
significance detection; see Figure 7a, b). With our MCMC fit, we derive a mass
of $6.3_{-1.4}^{+1.4}$$\,M_{\oplus}$ and find that the orbital period and
phase of the detection agrees with the known $9.978543$ day period from
transit observations.

For the new DRS (2.3.5), we find that the _CALM_ method does improve the
scatter in the RVs from from $\sigma_{\rm raw}$ = $3.01$ $\rm m\,s^{-1}$ to
$\sigma_{\rm us}$ = $2.41$$\rm m\,s^{-1}$. However, the level of stellar
activity is already sufficiently low in the new DRS that even without an
activity model, we can detect the planet signal at 4.3$\sigma$ significance
(Figure 7c,d). When applying our _CALM_ framework, the detection significance
remains at 4.3$\sigma$. Although this does not provide a notable improvement,
this result does demonstrate that our method effectively preserves planetary
Doppler signals rather than removing them to achieve a lower RV scatter. We
note that after applying our activity mitigation method, the final unexplained
scatter in the new DRS ($\sigma_{\rm us}$ = $2.41$$\rm m\,s^{-1}$) is
marginally higher than in the old DRS ($\sigma_{\rm us}$ = $2.04$$\rm
m\,s^{-1}$). This may be due to a difference in line lists between the old and
new DRS where less activity sensitive lines were used to derive the CCFs and
RVs in the new DRS. Although this may result in cleaner RVs, this may also
make it more challenging to remove remaining stellar activity signals with the
CCFs in the new DRS. In the future, using custom line lists that are
especially sensitive to stellar activity to compute custom CCFs may allow us
to mitigate the stellar activity more effectively using _CALM_.

Figure 10: Phasefolded two planet fit. a) The RV model for K2-167 b is plotted
in black and the binned RVs are plotted in yellow. The data is phase-folded on
K2-167 b’s period of $9.978543$ days. b) RV model for a hypothethical outer
companion is plotted in blue and the binned RVs are plotted in yellow. The
data is phase-folded on the fitted period of 22.01 days.

####  3.6.1 Results in Fourier Space: Activity Signal Peaks Decrease and
Planetary Signal Emerges

For both the old DRS and new DRS, we investigate the behavior of raw and
stellar activity corrected HARPS RVs in the Fourier domain to investigate
which (quasi)-periodic signals are being modeled and removed. To do this, we
compute the Lomb-Scargle Periodogram (Lomb, 1976; Scargle, 1982) of the RVs
before and after applying the _CALM_ activity correction. We used the
periodogram functions in astropy.timeseries (VanderPlas et al., 2012;
VanderPlas & Ivezić, 2015) to implement a generalized Lomb-Scargle
Periodogram, where the periodograms are normalized according to the formalism
in Zechmeister & Kürster (2009). In the periodograms corresponding to the old
DRS (Figure 8) and the new DRS (Figure 9), we indicate the expected planet
period ($9.978543$ days) from transits in yellow. Based on the v$\sin{i}=3.8$
km s-1 derived from the SPC analysis of HARPS-N spectra, we expect the
rotation period of the star to be $\sim 15$ days and half the rotation period
to be at $\sim 7.5$ days as indicated by $P_{\rm rot}$ and $P_{\rm rot}/2$ in
grey in Figures 8, 9.

For the old DRS (3.7), we plot Lomb-Scargle Periodograms of the RVs before and
after applying the activity correction models with increasing number of CCF
indexes in Figure 8. As we increase the number of CCF indexes (N) included in
the _CALM_ model from N=5 (Figure 8b), to N = 15 (Figure 8c), and finally to N
= 25 (Figure 8d), we see that the peaks corresponding to $P_{\rm rot}$,
$P_{\rm rot}/2$, and long-term stellar activity decrease in magnitude. A
strong signal emerges instead at the the expected planet period ($9.978543$
days). Thus, a larger number of CCF indexes (25) clearly more effectively
separates stellar activity signals from Doppler reflex motion and reveals
planetary signals. However, this does not imply that we should infinitely
increase the number of CCF indexes. As discussed in Section 3.4, increasing
the number of CCF indexes can result in too many free parameters compared to
the number of observations and overfitting on noise in the data. We explored
using 5-40 CCF indexes and find that $\sim$ 25 CCF indexes strikes the balance
between preventing overfitting and still modeling the stellar activity signals
effectively.

For the new DRS (2.3.5), we plot periodograms of the RVs before and after
applying the activity correction model with 5 CCF indexes in Figure 9. We do
not include several activity models with increasing number of CCF indexes
since we did not find that increasing the number of CCF indexes made any
notable difference in the periodograms for the new DRS. Although we do see the
magnitude of the peaks corresponding to $P_{\rm rot}$, $P_{\rm rot}/2$, and
long-term stellar activity decrease marginally in magnitude and a slightly
stronger peak corresponding to the planet period ($9.978543$ days) emerge in
9b, this effect is much less pronounced compared to the old DRS. We expect
that this is due to a difference in the use of linelists in the new and old
DRS to compute the CCFs and RVs as discussed in Section 3.6.

###  3.7 Assessing the possibility of a second planetary companion

In our activity corrected periodograms (bottom panels of Figures 8, 9), we
noticed a peak arising at $\sim$ 22 days that could be hinting at a secondary
companion. This period does not correspond to sidereal day, sidereal year, or
lunation period aliases of K2-167 b or the star’s rotation period and aliases.
Thus, we performed a 2 planet fit with our CCF linear activity model where
K2-167 b is modeled as a Keplerian and the possible companion is modeled as a
cosine function to keep free parameters to a minimum. Although our DE-MCMC
model converged (Gelman Rubin values of < 1.009) and the fit looks relatively
convincing as seen in Figure 10, we do not find a significant detection for
this potential second companion. From the corner plot showing the distribution
of possible periods (P2) in Figure A.2, we see that multiple periods between
21.6 and 22.8 could fit the observed data, which may help explain the lack of
significant detection. Additional RV measurements are necessary to constrain
the period and confirm or rule out this companion’s existence.

####  3.7.1 Comparing _CALM_ to traditional activity indicators

For the old DRS, which shows a larger scatter, we compared our CCF activity
indicators with traditional activity indicators such as s-index, H$\alpha$,
NaD, and other CCF metrics such as bisector (BIS), full-width-half-maximum
(FWHM), and contrast. In Figure 11, we evaluate these methods and compare them
with the CCF method by both examining how well the planet signal is preserved
and what the final remaining scatter in the RV measurements is. We use the
scatter metrics defined in detail in Section 3.5 for this comparison. We find
that using any combination of activity indexes and/or BIS, FWHM results in a
higher final unexplained scatter and a less significant detection of K2-167
b’s mass. In particular, using only the BIS, FWHM, and contrast results in the
highest final unexplained scatter of 3.18 $\rm m\,s^{-1}$ with the lowest
significance of detection of 2.92$\sigma$. The other activity indicator models
perform marginally better  and yield detections of 3.16$\sigma$ to
3.25$\sigma$. However, when we include CCFs, we see a significant reduction in
unexplained scatter and increased detection significance. Both the 25 CCF
_CALM_ model and the 25 CCF model with activity indexes (s-index, H$\alpha$,
NaD) yield a $\sigma_{\rm us}=2.0$ $\rm m\,s^{-1}$, which is the lowest
$\sigma_{\rm us}$ we find. The 25 CCF _CALM_ model has the higher detection
significance with 4.29$\sigma$ compared to 3.95$\sigma$ for the CCF with
activity indexes model. Thus, overall  the 25 CCF _CALM_ model leads to the
lowest $\sigma_{\rm us}$ and the highest detection significance among these
methods.

Figure 11: Comparison between using CCFs, CCF metrics (FWHM, bisector,
contrast), and activity indexes (s-index, CaII, NaD, Hα) as activity
correction models for the Old DRS (3.7). On the bottom panel, the corrected
scatter metrics ($\sigma_{\rm ac}$, $\sigma_{\rm us}$) for each of the
activity models are plotted. The activity corrected standard deviation
$\sigma_{\rm ac}$ measures the scatter in the activity corrected RV
measurements. The unexplained scatter ($\sigma_{\rm us}$) measures the scatter
remaining after removing all known contributions (stellar activity, planet
signal, estimated instrumental noise, photon noise) to the RV signal. These
scatter metrics are defined in more detail in Section 3.5. In the top panel,
we include the estimated semi-amplitude of K2-167 b for each activity
correction model and the corresponding error bars and detection significance
to check how well the planet signal is preserved for each activity mitigation
method.

Parameter | Units | Values
---|---|---
Observed Stellar Parameters
$\log{g}$ | Spectroscopic surface gravity (cgs) | $4.10\pm 0.05$
$T_{\rm eff}$ | Effective Temperature (K) | $6083\pm 90$
$[{\rm Fe/H}]$ | Metallicity (dex) | $-0.34\pm 0.09$
$v\sin{i}$ |  Projected rotational velocity ($\rm km\,s^{-1}$) | $3.8\pm 0.5$
Derived Stellar Parameters
$M_{\star}$ |  Mass ($M_{\odot}$) | $1.109^{+0.11}_{-0.091}$
$R_{\star}$ |  Radius ($R_{\odot}$) | $1.554^{+0.051}_{-0.050}$
$L_{\star}$ |  Luminosity ($L_{\odot}$) | $3.08^{+0.16}_{-0.17}$
$\rho_{\star}$ | Density (cgs) | $0.419^{+0.065}_{-0.059}$
$Age$ | Age (Gyr)a | $5.6^{+2.2}_{-2.0}$
$\log{g}$ | Model-derived surface gravity (cgs) | $4.102^{+0.052}_{-0.055}$
$u_{1}$ | Kepler-band linear limb-darkening coeff | $0.342\pm 0.050$
$u_{1}$ |  _TESS_ -band linear limb-darkening coeff | $0.242\pm 0.030$
$u_{2}$ | Kepler-band quadratic limb-darkening coeff | $0.311^{+0.049}_{-0.050}$
$u_{2}$ |  _TESS_ -band quadratic limb-darkening coeff | $0.298\pm 0.029$
Planetary Parameters: | b
$P$ | Period (days) | $9.978543^{+0.000023}_{-0.000020}$
$M_{P}$ |  Mass ($\,M_{\oplus}$) | $7.0^{+1.7}_{-1.7}$
$R_{P}$ |  Radius ($\,R_{\oplus}$) b | $2.33^{+0.17}_{-0.15}$
$T_{C}$ |  Time of conjunction ($\rm{BJD_{TDB}}$)c | $2456979.9329\pm 0.0016$
$a$ | Semi-major axis (AU) | $0.0939^{+0.0029}_{-0.0027}$
$i$ | Inclination (Degrees) | $86.61^{+1.5}_{-0.53}$
$e$ | Eccentricity | $0.30^{+0.16}_{-0.19}$
$\omega_{\star}$ | Argument of periastron (degrees) | $99^{+44}_{-37}$
$e\cos{\omega_{*}}$ |  | $-0.03\pm 0.14$
$e\sin{\omega_{*}}$ |  | $0.26^{+0.17}_{-0.23}$
$M_{P}$ $\sin i$ |  Minimum mass ($\,M_{\oplus}$) | $6.962^{+1.653}_{-1.621}$
$M_{P}/M_{*}$ | Mass ratio | $0.0000173\pm 0.0000044$
K |  RV semi-amplitude ($\rm m\,s^{-1}$) | $1.88^{+0.49}_{-0.48}$
$T_{eq}$ |  Equilibrium temperature (K)d | $1203^{+25}_{-27}$
$R_{P}/R_{\star}$ | Radius of planet in stellar radii | $0.01444^{+0.00077}_{-0.00057}$
$a/R_{\star}$ | Semi-major axis in stellar radii | $13.01\pm 0.64$
$d/R_{\star}$ | Planet/star separation at mid transit | $9.5^{+3.0}_{-2.3}$
$\delta$ |  Transit depth in Kepler $\left(R_{p}/R_{\star}\right)^{2}$ | $0.0002299^{+0.0000086}_{-0.0000085}$
$\delta$ |  Transit depth in _TESS_ $\left(R_{p}/R_{\star}\right)^{2}$ e | $0.0002235^{+0.0000097}_{-0.0000096}$
$T_{14}$ | Total transit duration (days) | $0.1507^{+0.0037}_{-0.0039}$
$b$ | Transit Impact parameter | $0.60^{+0.18}_{-0.37}$
$\langle F\rangle$ |  Incident Flux (109 erg s-1 cm-2) | $0.425^{+0.054}_{-0.045}$

Table 1: Photometric Properties and Median Values and 68% Confidence Interval
for Global EXOFASTv2 Models for the K2-167 system. See Table 3 in Eastman et
al. (2019a) for a detailed description of all parameters.

a) We note that age measurements derived from isochrones have been found to
have large systematic uncertainties (Torres et al., 2010).
b) We note that this is in equatorial earth radii.
c) Time of conjunction is commonly reported as the "transit time"
d) Assumes no albedo and perfect redistribution.
e) The difference in transit depths for Kepler and _TESS_ is is due to limb
darkening effects.

###  3.8 Final RV model

For our final RV model, we use the new DRS (2.3.5) RVs and the 5 CCF activity
indicator _CALM_ activity correction model since this yields the lowest final
unexplained scatter ($\sigma_{\rm us}$). The CCF activity indicators are then
used in the global transit, radial velocity, stellar parameter fit to
decorrelate the RVs as further described in Section 4.2. We do not include the
possible additional planet that we investigated as described in Section 3.7 in
our final RV model.

##  4 K2-167 System Parameters

In addition to performing the RV activity mitigation analysis, we also measure
the spectroscopic properties and use the CALM activity indicators to perform a
global transit and stellar parameter analysis to derive the refined K2-167
system parameters as described in this section.

###  4.1 Spectroscopic Parameters

Using all 74 of the HARPS-N spectra, we measure the spectroscopic properties
of K2-167 using three independent methods: ARES+MOOG131313ARESv2:
http://www.astro.up.pt/ sousasag/ares/,141414MOOG- 2017:
http://www.as.utexas.edu/ chris/moog.html (Sousa, 2014),
CCFpams151515https://github.com/LucaMalavolta/CCFpams (Malavolta et al.,
2017), and Stellar Parameter Classification (SPC) (Buchhave et al., 2014).

First, we measured the spectroscopic parameters using ARES+MOOG, which is a
curve-of-growth method that relies on neutral and ionized iron lines and is
further described by Mortier, et al. (2013). We co-added all the spectra, used
ARES2 (Sousa, 2014; Sousa et al., 2015) to measured equivalent widths of iron
lines, and used MOOG (Sneden, 1973) to determine atmospheric parameters. We
accounted for systematic effects in ARES+MOOG by using the surface gravity
corrected method from Mortier et al. (2014). Our systematic errors were added
to our precision errors in quadrature. These precision errors are intrinsic to
the method. The systematic errors are 60 K for the the effective temperature,
0.1 dex for the surface gravity, and 0.04 dex for the metallicity (Sousa et
al., 2011). This analysis yielded an effective temperature $T_{\rm
eff,MOOG}=6115\pm 69$ K, surface gravity $\log g_{\rm cgs,MOOG}=4.12\pm 0.13$,
and an iron abundance [Fe/H]${}_{\rm MOOG}=-0.32\pm 0.04$.

The CCFpams method (Malavolta et al., 2017), which relies on an empirical
calibration that uses the width of CCFs to determine the effective
temperature, surface gravity, and iron abundance. Using the same systematic
errors for effective temperature, surface gravity, and metallicity as used for
ARES+MOOG from Sousa et al. (2011), this method computed an effective
temperature $T_{\rm eff,CCFpams}=6092\pm{\color[rgb]{0,0,0}68}$ K, surface
gravity $\log g_{\rm cgs,CCFpams}=4.09\pm{\color[rgb]{0,0,0}0.12}$, and an
iron abundance [Fe/H]${}_{\rm CCFpams}=-0.34\pm{\color[rgb]{0,0,0}0.05}$.

We used the Stellar Parameter Classification (SPC) method (Buchhave et al.,
2014), which determines stellar atmospheric parameters by cross-correlating a
stellar spectrum with synthetic spectra from Kurucz (1992) model atmospheres
and then interpolating correlation peaks to find the effective temperature,
surface gravity, and metallicity. We ran SPC on all our HARPS-N spectra. Using
a prior on gravity of log $g_{\rm cgs,SPC}=4.12$, we measured a temperature of
$T_{\rm eff,SPC}=6047\pm 49$K, a surface gravity of log $g_{\rm
cgs,SPC}=4.12\pm 0.10$, an iron abundance $[Fe/H]_{\rm SPC}=-0.37\pm 0.08$,
and the projected rotational velocity v$\sin{i}_{\rm SPC}=3.8\pm 0.5$km s-1.

Lastly, we computed the weighted average across all three methods. The
weighted average yields an effective temperature $T_{\rm eff,avg}=6083\pm
90$K, a surface gravity of log $g_{\rm cgs,avg}=4.10\pm 0.05$, and an iron
abundance $[Fe/H]_{\rm avg}=-0.34\pm 0.09$. We adopt these values and include
$T_{\rm eff,avg}$ and $[Fe/H]_{\rm avg}$ as Gaussian priors in our global
transit and stellar parameter analysis.

###  4.2 Global Transit/Stellar Parameter Analysis

To refine the ephemerides and system parameters of K2-167 b, we used the
fitting software EXOFASTv2 (Eastman et al., 2013; Eastman, 2017; Eastman et
al., 2019b), which allows us to combine transit observations from several
missions with our spectroscopic parameters and archival photometry. In
particular, we simultaneously fit the _TESS_ , K2 photometry, RVs, and stellar
parameters (see Section 2). In addition to including the RVs, we include the
CCF-based activity indicators derived from the _CALM_ model (described in
Section 3.2) and use these to detrend the RVs in the EXOFASTv2 fit. Since this
detrending is identical to performing a linear regression, incorporating the
CCF activity indicators in EXOFASTv2 is identical to using our _CALM_ model
for the RVs. In this fit, we constrain the possible mass, radius, and age of
the host star according to the MESA Isochrones and Stellar Tracks (MIST) and
stellar evolution models (Paxton et al., 2011; Paxton et al., 2013, 2015;
Dotter, 2016; Choi et al., 2016).

We required a Gelman-Rubin statistic (Gelman & Rubin, 1992) of less than 1.01
and at least 1000 independent draws in each parameter to ensure strict
convergence criteria. We imposed Gaussian priors on $[Fe/H]$ and $T_{\rm eff}$
using the weighted averages derived from our spectroscopic analysis described
in Section 4.1. The results of our fit are shown in Figure 2 and Table 1. In
the EXOFASTv2 fit, we find $M_{p}=$ $7.0^{+1.7}_{-1.7}$$\,M_{\oplus}$ and thus
we find close agreement between the _CALM_ estimate of mass
($6.3_{-1.4}^{+1.4}$$\,M_{\oplus}$) and the mass measured by Bonomo et al.
(2023) ($6.5_{-1.5}^{+1.6}$$\,M_{\oplus}$).

###  4.3 Stellar variability signatures in the photometry

We examined both the K2 and TESS light curves to look for stellar rotation
signals. We did not find significant rotation signals in the TESS data. The K2
light curves have higher photometric precision and probe a bluer wavelength
region compared to TESS and thus should be able to probe more stellar
variability. For both the TESS and K2 data, We used the unflattened light
curve that is not corrected for systematics since the systematic corrections
often also remove stellar variability signals. We observe less than 0.1%
variability in the K2 light cuvres. In Figure 12, we examine the K2 light
curve in Fourier space and find no significant signals at the expected
rotation period ($P_{\rm rot}\sim 15$ days). We find that the periodogram is
dominated by K2 instrumental systematics. Significant signal loss is expected
for K2 due to long-term systematic noise (e.g. differential velocity
aberration) for periods greater than 15 days (Van Cleve et al., 2016) and make
it challenging to detect oscillations and rotation signals.

Figure 12: Periodogram of K2 photometry. In grey, the expected stellar
rotation period ($P_{\rm rot}$) of $\sim$ 15 days is indicated and the
instrumental systematics that dominate signals with P $>20$ days are
indicated. The instrumental systematics dominate the Fourier spectrum of the
K2 light curve.  Figure 13: Residual CCFs ($\Delta$CCFs) from the new DRS
(2.3.5) for two example stars, K2-2 (Thygeson et al. in prep) and LTT 3780
(Cloutier et al., 2020). Although the new DRS (2.3.5) did not show a clear
pattern in shape changes for K2-167, we can see that other targets do clearly
exhibit shape changes corresponding to stellar variability and could
potentially benefit from our _CALM_ method.

##  5 Discussion

In this work, we have measured the mass of K2-167 b as $6.3_{-1.4}^{+1.4}$
$\,M_{\oplus}$ by taking advantage of the reduction in scatter resulting from
updates to the HARPS-N Data Reduction System (2.3.5) and our new stellar
activity mitigation method (_CALM_). This measurement agrees within error with
the recently published mass measurement of Bonomo et al. (2023)
($6.5_{-1.5}^{+1.6}$$\,M_{\oplus}$), which was part of a larger statistical
sample of exoplanet systems. Together with our radius measurement
($2.33^{+0.17}_{-0.15}$$\,R_{\oplus}$), we can use our precise mass
measurement to place constraints on K2-167 b’s composition. We compared our
measured mass and radius to mass-radius relations for planets of different
compositions from Zeng et al. (2016). In Figure 14, we demonstrate this newly
determined mass in a mass/radius diagram. The mass and composition of K2-167 b
is not consistent with an Earth-like composition (32.5% Iron core, 67.5%
silicate mantle). Rather, K2-167 b is less dense than earth-like planets and
falls closer to the pure ice (H2O) mass-radius relation from Zeng et al.
(2016).

This leaves several possibilities for the exact composition of K2-167 b. This
planet could be a mixture of rock, iron, and hydrogen/helium (H/He). One
possibility would be that K2-167 b has an earth-like iron core and mantle that
is surrounded by a volatile envelope, which is believed to be a common
composition for super-Earths and sub-Neptunes (Weiss & Marcy, 2014; Rogers,
2015; Dressing et al., 2015). Another possibility is that K2-167 b is a planet
with a much smaller iron core than earth and mostly or entirely consists of
water or methane or some other heavy volatiles.

Despite the fact that this planet was originally detected over 8 years ago,
only recently has it been possible to make such inferences about its
composition. Like many other low-mass planets, K2-167 b was first detected
using the transit method and was validated several years ago but stellar
activity prevented astronomers from measuring its mass until recently. To
ensure the robustness of K2-167 b’s mass measurement, we performed an
extensive RV analysis where we compared our _CALM_ framework to traditional
activity mitigation methods (H-$\alpha$, s-index, FWHM, BIS, etc). Instead of
using these indicators, our method exploits the activity-induced line shape
changes in the spectra without requiring timing information like Gaussian
Process regression. We trace these shape changes using CCFs and find that this
method either outperforms or produces comparable results to traditional
methods.

The fact that different data reduction software versions, which use different
linelists, yield radial velocities with significantly distinct stellar
activity characteristics suggests interesting future directions in terms of
mitigating stellar activity. The effectiveness of spectra-based activity
mitigation methods has already been shown to be highly dependent on linelists
and McWilliam et al. in prep recently found that using chromatic CCFs, which
are CCFs that use only part of the spectral orders, add complementary
information for stellar activity mitigation for the Sun. Thus, in the future,
we plan to generate our our line profiles generated using primarily activity
sensitive lines such that we can apply _CALM_ to these profiles rather than to
CCFs from the HARPS-N DRS, which may greatly increase the effectiveness of
spectra-based activity mitigation techniques like _CALM_.

###  5.1 Prospects for using _CALM_ on other stars with the new DRS

In the future, we plan to apply our _CALM_ method to variety of others stars
of different stellar types and different ages. Although the new DRS (2.3.5)
RVs for K2-167 did not require our method to derive a mass a measurement, this
is not the case for other targets. From preliminary analyses, we have found
that _CALM_ can be effective at predicting and mitigating stellar variability
for several stars with data from the new DRS (2.3.5). In Figure 13, we show
two examples of the residual CCF observations for K2-2 (Thygeson at el.
submitted), a K0-type star, and LTT 3790 (Cloutier et al., 2020), an M4-type
star. Both these targets have planets confirmed with transits and we find that
_CALM_ can reduce the RV scatter while preserving and increasing the signal
strength of the planetary reflex motion.

In addition to using _CALM_ to model magnetic activity, we are hoping to
develop methods to model granulation, which is especially relevant for stars
with high surface gravity. The granulation-induced RV variation is anti-
correlated with flicker, which are low-level photometric variations induced by
granulation (Bastien et al., 2013). Recently, Lakeland et al. (2023) found
that granulation and supergranulation are some of the main remaining limiting
factors in reaching tens of $\rm cm\,s^{-1}$ level for solar observations.
Thus, it is critical that we find ways to model and predict granulation and
supergranulation, especially for solar-like stars.

###  5.2 Applying _CALM_ to young planetary systems

Improved activity mitigation methods are especially interesting for helping to
constrain the mechanisms involved in creating the radius valley. One of the
best ways to distinguish the most commonly proposed mechanisms,
photoevaporation and core-powered mass loss, is by investigating the timescale
on which mass loss takes place. Since photoevaporation is expected to proceed
on timescales shorter than 100 million years, measuring the compositions of
planets younger than this age is of strong interest. Unfortunately, stars this
young typically have very high levels of stellar activity, making their
precise characterization difficult and making them prone to overfitting using
current state-of-the art techniques like Gaussian Processes (Blunt et al.,
2023). Applying _CALM_ or similar techniques to these young stars could open
the door to measuring more masses and densities for transiting planets
orbiting bright young stars (i.e. EPIC 247589423 (Mann et al., 2018), HD
283869 (Vanderburg et al., 2018), V1298 Tau (David et al., 2019), DS Tuc A
(Newton et al., 2019; Benatti et al., 2019)) and provide insight into the
timescales of mass loss for these planets.

###  5.3 Future further characterization of K2-167 b

In the future, additional RV observations of K2-167 b could allow us to use a
larger portion of the average line spectrum and potentially improve our
modeling of the stellar variability. If timed during transit, additional RV
measurements could allow us to measure the spin-orbit angle using the
Rossiter-McLaughlin (R-M) effect (Rossiter, 1924; McLaughlin, 1924). We
compute the expected RM radial velocity deviation $\Delta RV_{\rm RM}$ using
the approximation from Seager & Dotson (2010) that

$\Delta RV_{\rm RM}=(\frac{Rp}{R_{\star}})^{2}\sqrt{1-b^{2}}vsini$ (6)

For K2-167 b, the expected RM radial velocity deviation is $\Delta RV_{\rm
RM}\sim 20.1$ $\rm m\,s^{-1}$and should be measure-able with precision and
high-precision spectographs. R-M measurements could provide insight into the
angular momentum history and thereby dynamical history of the K2-167 system.

For future atmospheric target selection studies, we calculated the
Transmission Spectroscopy Metric (TSM; Kempton et al., 2018) and the Emission
Spectroscopy Metric (ESM; Kempton et al., 2018) for K2-167 b. We find that
$TSM=6.06$ and $ESM=4.12$, which means that K2-167 b is not an ideal target
for transmission and a moderately favorable target for emission.

Figure 14: The mass/radius diagram for small exoplanets. Planet masses and
radii come from the NASA Exoplanet Archive(Akeson et al., 2013), accessed May
7, 2023. The darkness of each of the symbols is directly proportional to the
precision of the mass and radii measurements. The colored lines represent
theoretical mass/radius relations for solid planets of various compositions
from (Zeng et al., 2016) and for cold hydrogen planets from (Seager et al.,
2007). The planets in our solar system are plotted in blue and K2-167 b is
plotted in purple with their corresponding uncertainties.

##  6 Conclusion

In this paper, we have developed a new spectral-based activity mitigation
method called _CALM_ and used it to characterize the K2-167 system, which has
a planet at the edge of the radius valley. In-depth characterization of
systems like this can be especially critical in probing the formation physics
of the radius valley. In this paper, we have specifically obtained the
following results:

(1) We performed a detailed RV analysis, where we compared _CALM_ to various
traditional activity indicators, evaluated various HARPS-N DRS versions, and
tested a range of number of CCF activity indicators within the _CALM_
framework. Across all of these analysis choices, we derive a consistent mass
within uncertainty and find that the period agrees with the transit period of
P = $9.978543$ days. Our final _CALM_ model derives a mass of $M_{p}=$
$6.3_{-1.4}^{+1.4}$$\,M_{\oplus}$.

(2) We find that the effectiveness of _CALM_ at reducing the RV scatter is
highly dependent on the version of the HARPS-N DRS used, which suggests that
this is highly dependent on the choice of linelists. In the future, custom
linelists to compute CCFs may greatly increase the the effectiveness of
spectra based techniques like _CALM_.

(3) We investigated our RV results in Fourier space using periodograms to
understand which periodic signals were responsible for the observed decrease
in RV scatter. Using _CALM_ , we find that both stellar rotation signals and
long-term magnetic activity signals decrease in magnitude while preserving the
planet signal. This effect is most prominent in the old DRS (3.7). In the new
DRS (2.3.5), we do not see the signal strength of activity signals signal
decrease significantly, but the planet signal is preserved nonetheless. We
expect that this difference is the result of the different linelists used in
the different DRS versions.

(4) We performed a joint fit with the RVs and transit using EXOFASTv2. We
found consistent mass estimates and refined the radius to $R_{p}=$
$2.33^{+0.17}_{-0.15}$$\,R_{\oplus}$, which places K2-167 b at the upper edge
of the radius valley. We refined the stellar parameters using three
spectroscopic parameter estimation methods, and computed the weighted average
across these methods, which yields an effective temperature $T_{\rm
eff,avg}=6083\pm 90$K, a surface gravity of log $g_{\rm cgs,avg}=4.10\pm
0.05$, and an iron abundance $[Fe/H]_{\rm avg}=-0.34\pm 0.09$.

(5) We investigated the possibility of a secondary companion at a $\sim$22 day
orbital period, but we found degeneracies in the possible periods and no
significant detection. Further RV follow-up may help refine the period and
confirm or rule out a possible companion.

(6) We introduced a new scatter metric called the unexplained scatter that
accounts for both modelled stellar variability and modelled planetary signals
to asses model performance. We introduce this metric since commonly used
scatter metrics do not account for scatter originating from un-modelled
planets and sometimes inadvertently incentivize removing planet signals rather
than preserving them.

(7) We compared our _CALM_ model to traditional activity indicators (s-index,
H$\alpha$, NaD, bisector of the CCF, FWHM of the CCF, contrast of the CCF) and
find that our models outperform any combination of these indicators for the
old DRS. Our CCF-based _CALM_ model both reduces the unexplained RV scatter
and results in more significant mass detections.

(8) Combining our mass and radius measurements, we place constraints on K2-167
b’s composition and find that it is less dense than would be expectected for
an earth-like composition. This means that it could have either (i) an earth-
like iron core and mantle with a volatile envelope, or (ii) a much smaller
iron core and further consist of water, methane, or other heavy volatiles.

(9) Furthermore, we demonstrate that although _CALM_ only yields a more
significant mass detection for the old DRS (3.7) and not the new DRS (2.3.5),
this is not the case for several other stars (e.g. K-type stars K2-2, M-type
star LTT 3790) where we do reduce the RV scatter and increase the mass
detection significance using the new DRS.

In the future, continued development and broad deployment of stellar activity
mitigation algorithms will have wide-ranging applications in exoplanet
science. Methods like these could unlock the ability to measure masses of
planets orbiting active stars spanning the HR diagram and pave the way towards
precise characterization of larger exoplanet populations in the super-Earth
and sub-Neptune regime.

## Acknowledgements

We acknowledge helpful conversations and feedback from members of Dave
Latham’s Coffee Club. We thank Dr. Sai Ravela for fruitful discussions on the
best ways to choose the model input parameters and suggesting future
directions for refining our process. The HARPS-N project has been funded by
the Prodex Program of the Swiss Space Office (SSO), the Harvard University
Origins of Life Initiative (HUOLI), the Scottish Universities Physics Alliance
(SUPA), the University of Geneva, the Smithsonian Astrophysical Observatory
(SAO), and the Italian National Astrophysical Institute (INAF), the University
of St Andrews, Queen’s University Belfast, and the University of Edinburgh.

ZLD would like to thank the generous support of the MIT Presidential
Fellowship, the MIT Collamore-Rogers Fellowship and to acknowledge that this
material is based upon work supported by the National Science Foundation
Graduate Research Fellowship under Grant No. 1745302. ZLD and AV acknowledge
support from the TESS Guest Investigator Program under NASA grant
80NSSC19K0388. ZLD and AV acknowledge support from the D.17 Extreme Precision
Radial Velocity Foundation Science program under NASA grant 80NSSC22K0848.

AMS would like to acknowledge that this work was supported by FCT - Fundação
para a Ciência e a Tecnologia through national funds by these grants:
UIDB/04434/2020; UIDP/04434/2020". Funded/Co-funded by the European Union
(ERC, FIERCE, 101052347). Views and opinions expressed are however those of
the author(s) only and do not necessarily reflect those of the European Union
or the European Research Council. Neither the European Union nor the granting
authority can be held responsible for them. A.M.S acknowledges support from
the Fundação para a Ciência e a Tecnologia (FCT) through the Fellowship
2020.05387.BD.

R.D.H. is funded by the UK Science and Technology Facilities Council (STFC)’s
Ernest Rutherford Fellowship (grant number ST/V004735/1).

This research has made use of the NASA Exoplanet Archive, which is operated by
the California Institute of Technology, under contract with the National
Aeronautics and Space Administration under the Exoplanet Exploration Program.

## Data Availability

All RVs and activity indicators used in our analyses will be available via
VizieR at CDS. The _TESS_ and _K2_ photometry data are accesible via the
Mikulski Archive for Space Telescopes portal at
https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html.

## References

  * Akeson et al. (2013) Akeson R. L., et al., 2013, Publications of the Astronomical Society of the Pacific, 125, 989
  * Barragán et al. (2019) Barragán O., et al., 2019, MNRAS, 490, 698
  * Barragán et al. (2023) Barragán O., et al., 2023, MNRAS, 522, 3458
  * Bastien et al. (2013) Bastien F. A., Stassun K. G., Basri G., Pepper J., 2013, Nature, 500, 427
  * Batalha et al. (2013) Batalha N. M., et al., 2013, ApJS, 204, 24
  * Benatti et al. (2019) Benatti S., et al., 2019, A&A, 630, A81
  * Berger et al. (2018) Berger T. A., Huber D., Gaidos E., van Saders J. L., 2018, ApJ, 866, 99
  * Blunt et al. (2023) Blunt S., et al., 2023, AJ, 166, 62
  * Boisse et al. (2009) Boisse I., et al., 2009, A&A, 495, 959
  * Bonfils et al. (2007) Bonfils X., et al., 2007, A&A, 474, 293
  * Bonomo et al. (2023) Bonomo A. S., et al., 2023, arXiv e-prints, p. arXiv:2304.05773
  * Borucki et al. (2010) Borucki W. J., et al., 2010, Science, 327, 977
  * Buchhave et al. (2014) Buchhave L. A., et al., 2014, Nature, 509, 593
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Cloutier et al. (2020) Cloutier R., et al., 2020, AJ, 160, 3
  * Colwell et al. (2023) Colwell I., Timmaraju V., Wise A., 2023, arXiv e-prints, p. arXiv:2304.04807
  * Cosentino et al. (2012) Cosentino R., et al., 2012, in Proc. SPIE. p. 84461V, doi:10.1117/12.925738
  * Dattilo et al. (2023) Dattilo A., Batalha N. M., Bryson S., 2023, arXiv e-prints, p. arXiv:2308.00103
  * David et al. (2019) David T. J., et al., 2019, AJ, 158, 79
  * Delisle et al. (2018) Delisle J. B., et al., 2018, A&A, 614, A133
  * Deming et al. (2015) Deming D., et al., 2015, ApJ, 805, 132
  * Dotter (2016) Dotter A., 2016, ApJS, 222, 8
  * Dressing et al. (2015) Dressing C. D., et al., 2015, ApJ, 800, 135
  * Duck et al. (2021) Duck A., Harada C. K., Harrell J., Morris R. R. A., Williams E., Crossfield I., Werner M., Deming D., 2021, AJ, 162, 136
  * Dumusque et al. (2011) Dumusque X., et al., 2011, A&A, 535, A55
  * Dumusque et al. (2021) Dumusque X., et al., 2021, A&A, 648, A103
  * Eastman (2017) Eastman J., 2017, EXOFASTv2: Generalized publication-quality exoplanet modeling code (ascl:1710.003)
  * Eastman et al. (2013) Eastman J., Gaudi B. S., Agol E., 2013, PASP, 125, 83
  * Eastman et al. (2019a) Eastman J. D., et al., 2019a, arXiv e-prints, p. arXiv:1907.09480
  * Eastman et al. (2019b) Eastman J. D., et al., 2019b, arXiv e-prints, p. arXiv:1907.09480
  * Fazio et al. (2004) Fazio G. G., et al., 2004, ApJS, 154, 10
  * Figueira et al. (2013) Figueira P., Santos N. C., Pepe F., Lovis C., Nardetto N., 2013, A&A, 557, A93
  * Fulton et al. (2017) Fulton B. J., et al., 2017, AJ, 154, 109
  * Fulton et al. (2018) Fulton B. J., Petigura E. A., Blunt S., Sinukoff E., 2018, PASP, 130, 044504
  * Garhart et al. (2020) Garhart E., et al., 2020, AJ, 159, 137
  * Gelman & Rubin (1992) Gelman A., Rubin D. B., 1992, Statistical Science, 7, 457
  * Gilbertson et al. (2020) Gilbertson C., Ford E. B., Jones D. E., Stenning D. C., 2020, ApJ, 905, 155
  * Ginzburg et al. (2018) Ginzburg S., Schlichting H. E., Sari R., 2018, MNRAS, 476, 759
  * Gupta & Schlichting (2020) Gupta A., Schlichting H. E., 2020, MNRAS, 493, 792
  * Haywood et al. (2014) Haywood R. D., et al., 2014, MNRAS, 443, 2517
  * Haywood et al. (2020) Haywood R. D., et al., 2020, arXiv e-prints, p. arXiv:2005.13386
  * Haywood et al. (2022) Haywood R. D., et al., 2022, ApJ, 935, 6
  * Hébrard et al. (2010) Hébrard G., et al., 2010, A&A, 512, A46
  * Householder & Weiss (2022) Householder A., Weiss L., 2022, arXiv e-prints, p. arXiv:2212.06966
  * Ikwut-Ukwa et al. (2020) Ikwut-Ukwa M., et al., 2020, AJ, 160, 209
  * Jenkins (2002) Jenkins J. M., 2002, ApJ, 575, 493
  * Jenkins et al. (2016) Jenkins J. M., et al., 2016, in Chiozzi G., Guzman J. C., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9913, Software and Cyberinfrastructure for Astronomy IV. p. 99133E, doi:10.1117/12.2233418
  * Jin et al. (2014) Jin S., Mordasini C., Parmentier V., van Boekel R., Henning T., Ji J., 2014, ApJ, 795, 65
  * Jones et al. (2017) Jones D. E., Stenning D. C., Ford E. B., Wolpert R. L., Loredo T. J., Dumusque X., 2017, arXiv e-prints, p. arXiv:1711.01318
  * Kempton et al. (2018) Kempton E. M. R., et al., 2018, PASP, 130, 114401
  * Kurucz (1992) Kurucz R. L., 1992, in Barbuy B., Renzini A., eds, The Stellar Populations of Galaxies. 149. p. 225
  * Lakeland et al. (2023) Lakeland B. S., et al., 2023, MNRAS,
  * Lee et al. (2022) Lee E. J., Karalis A., Thorngren D. P., 2022, ApJ, 941, 186
  * Lienhard et al. (2023) Lienhard F., Mortier A., Cegla H. M., Cameron A. C., Klein B., Watson C. A., 2023, MNRAS, 522, 5862
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration et al., 2018, Lightkurve: Kepler and TESS time series analysis in Python (ascl:1812.013)
  * Ligi et al. (2018) Ligi R., et al., 2018, AJ, 156, 182
  * Lomb (1976) Lomb N. R., 1976, Ap&SS, 39, 447
  * Lopez & Fortney (2013) Lopez E. D., Fortney J. J., 2013, ApJ, 776, 2
  * Lovis et al. (2011) Lovis C., et al., 2011, arXiv e-prints, p. arXiv:1107.5325
  * Malavolta et al. (2017) Malavolta L., Lovis C., Pepe F., Sneden C., Udry S., 2017, MNRAS, 469, 3965
  * Mann et al. (2018) Mann A. W., et al., 2018, AJ, 155, 4
  * Mayo et al. (2018) Mayo A. W., et al., 2018, AJ, 155, 136
  * Mayor et al. (2003) Mayor M., et al., 2003, The Messenger, 114, 20
  * Mayor et al. (2011) Mayor M., et al., 2011, arXiv e-prints, p. arXiv:1109.2497
  * McLaughlin (1924) McLaughlin D. B., 1924, ApJ, 60, 22
  * Meunier & Lagrange (2013) Meunier N., Lagrange A. M., 2013, A&A, 551, A101
  * Milbourne et al. (2021) Milbourne T. W., et al., 2021, ApJ, 920
  * Mortier, et al. (2013) Mortier, Santos, N. C. Sousa, S. G. Fernandes, J. M. Adibekyan, V. Zh. Delgado Mena, E. Montalto, M. Israelian, G. 2013, A&A, 558, A106
  * Mortier et al. (2014) Mortier A., Sousa S. G., Adibekyan V. Z., Brandão I. M., Santos N. C., 2014, A&A, 572, A95
  * National Academies of Sciences, Engineering, and Medicine and others (2018) National Academies of Sciences, Engineering, and Medicine and others 2018, Exoplanet Science Strategy. National Academies Press
  * Newton et al. (2019) Newton E. R., et al., 2019, ApJ, 880, L17
  * Nicholson & Aigrain (2022) Nicholson B. A., Aigrain S., 2022, MNRAS, 515, 5251
  * Owen & Wu (2013) Owen J. E., Wu Y., 2013, ApJ, 775, 105
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Perger et al. (2023) Perger M., et al., 2023, A&A, 672, A118
  * Petigura (2020) Petigura E. A., 2020, AJ, 160, 89
  * Petigura et al. (2022) Petigura E. A., et al., 2022, AJ, 163, 179
  * Queloz et al. (2001) Queloz D., et al., 2001, A&A, 379, 279
  * Queloz et al. (2009) Queloz D., et al., 2009, A&A, 506, 303
  * Rajpaul et al. (2015) Rajpaul V., Aigrain S., Osborne M. A., Reece S., Roberts S., 2015, MNRAS, 452, 2269
  * Ricker et al. (2015) Ricker G. R., et al., 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Robertson et al. (2014) Robertson P., Mahadevan S., Endl M., Roy A., 2014, Science, 345, 440
  * Rogers (2015) Rogers L. A., 2015, ApJ, 801, 41
  * Rogers et al. (2021) Rogers J. G., Gupta A., Owen J. E., Schlichting H. E., 2021, MNRAS, 508, 5886
  * Rossiter (1924) Rossiter R. A., 1924, ApJ, 60, 15
  * Saar & Fischer (2000) Saar S. H., Fischer D., 2000, ApJ, 534, L105
  * Saar et al. (1998) Saar S. H., Butler R. P., Marcy G. W., 1998, ApJ, 498, L153
  * Scargle (1982) Scargle J. D., 1982, ApJ, 263, 835
  * Seager & Dotson (2010) Seager S., Dotson R., 2010, Exoplanets. The University of Arizona space science series, University of Arizona Press, Tucson
  * Seager et al. (2007) Seager S., Kuchner M., Hier-Majumder C. A., Militzer B., 2007, ApJ, 669, 1279
  * Shallue & Vanderburg (2018) Shallue C. J., Vanderburg A., 2018, AJ, 155, 94
  * Smith et al. (2012) Smith J. C., et al., 2012, PASP, 124, 1000
  * Sneden (1973) Sneden C., 1973, ApJ, 184, 839
  * Sousa (2014) Sousa S. G., 2014, ARES + MOOG: A Practical Overview of an Equivalent Width (EW) Method to Derive Stellar Parameters. pp 297–310, doi:10.1007/978-3-319-06956-2_26
  * Sousa et al. (2011) Sousa S. G., Santos N. C., Israelian G., Mayor M., Udry S., 2011, A&A, 533, A141
  * Sousa et al. (2015) Sousa S. G., Santos N. C., Adibekyan V., Delgado-Mena E., Israelian G., 2015, A&A, 577, A67
  * Stevenson et al. (2012) Stevenson K. B., et al., 2012, ApJ, 754, 136
  * Stumpe et al. (2012) Stumpe M. C., et al., 2012, PASP, 124, 985
  * Stumpe et al. (2014) Stumpe M. C., Smith J. C., Catanzarite J. H., Van Cleve J. E., Jenkins J. M., Twicken J. D., Girouard F. R., 2014, PASP, 126, 100
  * Ter Braak (2006) Ter Braak C. J. F., 2006, Statistics and Computing, 16, 239
  * Thygesen et al. (2023) Thygesen E., et al., 2023, AJ, 165, 155
  * Torres et al. (2010) Torres G., Andersen J., Giménez A., 2010, A&ARv, 18, 67
  * Tran et al. (2023) Tran Q. H., Bedell M., Foreman-Mackey D., Luger R., 2023, arXiv e-prints, p. arXiv:2305.00988
  * Van Cleve et al. (2016) Van Cleve J. E., et al., 2016, PASP, 128, 075002
  * Van Eylen et al. (2018) Van Eylen V., Agentoft C., Lundkvist M. S., Kjeldsen H., Owen J. E., Fulton B. J., Petigura E., Snellen I., 2018, MNRAS, 479, 4786
  * VanderPlas & Ivezić (2015) VanderPlas J. T., Ivezić Ž., 2015, ApJ, 812, 18
  * VanderPlas et al. (2012) VanderPlas J., Connolly A. J., Ivezic Z., Gray A., 2012, in Proceedings of Conference on Intelligent Data Understanding (CIDU. pp 47–54 (arXiv:1411.5039), doi:10.1109/CIDU.2012.6382200
  * Vanderburg (2021) Vanderburg A., 2021, avanderburg/edmcmc: v1.0.0, doi:10.5281/zenodo.5599854, https://doi.org/10.5281/zenodo.5599854
  * Vanderburg & Johnson (2014a) Vanderburg A., Johnson J. A., 2014a, PASP, 126, 948
  * Vanderburg & Johnson (2014b) Vanderburg A., Johnson J. A., 2014b, PASP, 126, 948
  * Vanderburg et al. (2016) Vanderburg A., et al., 2016, ApJS, 222, 14
  * Vanderburg et al. (2018) Vanderburg A., et al., 2018, AJ, 156, 46
  * Weiss & Marcy (2014) Weiss L. M., Marcy G. W., 2014, ApJ, 783, L6
  * Wise et al. (2022) Wise A., Plavchan P., Dumusque X., Cegla H., Wright D., 2022, ApJ, 930, 121
  * Zechmeister & Kürster (2009) Zechmeister M., Kürster M., 2009, A&A, 496, 577
  * Zeng et al. (2016) Zeng L., Sasselov D. D., Jacobsen S. B., 2016, The Astrophysical Journal, 819, 127
  * Zhao et al. (2022) Zhao L. L., et al., 2022, AJ, 163, 171
  * de Beurs et al. (2022) de Beurs Z. L., et al., 2022, AJ, 164, 49

## Appendix

Figure A.1: Understanding CCF Shape Diagnostics. a) Simulated CCFs that are
red-shifted (red), blue-shifted (blue) and the median CCF (grey) are plotted.
These CCFs are not shifted to the center. b) The residual CCF (simulated CCFs
minus the median CCF) are plotted and are equivalent to the derivatives of a
Gaussian. c) A red-shifted (red), blue-shifted (blue) and median CCF (grey)
are plotted. However, in this case  they are all centered at the median RV and
we simulated some white Gaussian noise as a simple approximation to shape
changes induced by stellar activity. d) The residual CCF (simulated CCFs minus
the median CCF) are plotted and do not resemble the derivatives of a Gaussian.

In this appendix, we include the radial velocity measurements included in this
analysis (Table A.1). In addition, we provide a supplemental figure that
visually explains the CCF shape diagnostics used to prevent overfitting in the
RV analysis (Figure A.1. Lastly, we provide an MCMC corner plot for the two
planet RV fit where the the period of a hypothethical companion is poorly
constrained (Figure A.2).

BJD |  RV [$\rm m\,s^{-1}$] |  RVerror [$\rm m\,s^{-1}$]
---|---|---
2457618.62 | -17549.45 | 0.89
2457618.55 | -17548.22 | 1.17
2457270.54 | -17555.49 | 0.79
2457269.56 | -17555.69 | 1.2
2457574.74 | -17552.52 | 0.89
. | . | .
. | . | .
. | . | .
2457565.71 | -17553.78 | 1.29
2457562.71 | -17551.56 | 1.2
2457322.45 | -17548.57 | 0.81
2457356.3 | -17548.43 | 1.41
2457354.36 | -17546.95 | 1.26
Table A.1: Radial velocity observations and estimated instrumental errors for
K2-167 using the HARPS-N DRS 2.3.5. The full table is available online. Figure
A.2: DE-MCMC corner plot for two planet fit. For K2-167 b the free parameters
are semi-amplitude (Amplitude), period (P), Time of periastron passage (Tp).
Eccentricity is frozen at e=0 for this fit. The potential companion is modeled
by a cosine curve and thus the free parameters are semi-amplitude (K2), period
(P2), and phase (phi2). For P2, we can clearly see that the data can be
consistent with multiple possible periods.