# KELT-9 as an eclipsing double-lined spectroscopic binary: a unique and self-
consistent solution to the system

Anusha Pai Asnodkar The Ohio State University, McPherson Laboratory, 140 W 18th Ave., Columbus, OH 43210, USA  Ji Wang (王吉) The Ohio State University, McPherson Laboratory, 140 W 18th Ave., Columbus, OH 43210, USA  B. Scott Gaudi  The Ohio State University, McPherson Laboratory, 140 W 18th Ave., Columbus, OH 43210, USA  P. Wilson Cauley Laboratory for Atmospheric and Space Physics, University of Colorado Boulder, Boulder, CO 80303  Jason D. Eastman Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA  Ilya Ilyin  Leibniz-Institut for Astrophysics Potsdam (AIP), An der Sternwarte 16, D–14482 Potsdam, Germany  Klaus Strassmeier Leibniz-Institute for Astrophysics Potsdam (AIP), An der Sternwarte 16, D–14482 Potsdam, Germany  Thomas Beatty  Department of Astronomy and Steward Observatory, University of Arizona, Tucson, AZ 85721 Anusha Pai Asnodkar<EMAIL_ADDRESS>
(Received August 20, 2021; Revised October 19, 2021; Accepted to AJ October
21, 2021)

###### Abstract

Transiting hot Jupiters present a unique opportunity to measure absolute
planetary masses due to the magnitude of their radial velocity signals and
known orbital inclination. Measuring planet mass is critical to understanding
atmospheric dynamics and escape under extreme stellar irradiation. Here, we
present the ultra-hot Jupiter system, KELT-9, as a double-lined spectroscopic
binary. This allows us to directly and empirically constrain the mass of the
star and its planetary companion, without reference to any theoretical stellar
evolutionary models or empirical stellar scaling relations. Using data from
the PEPSI, HARPS-N, and TRES spectrographs across multiple epochs, we apply
least-squares deconvolution to measure out-of-transit stellar radial
velocities. With the PEPSI and HARPS-N datasets, we measure in-transit planet
radial velocities using transmission spectroscopy. By fitting the circular
orbital solution that captures these Keplerian motions, we recover a planetary
dynamical mass of $2.17\pm 0.56\ \mathrm{M_{J}}$ and stellar dynamical mass of
$2.11\pm 0.78\ \mathrm{M_{\odot}}$, both of which agree with the discovery
paper. Furthermore, we argue that this system, as well as systems like it, are
highly overconstrained, providing multiple independent avenues for empirically
cross-validating model-independent solutions to the system parameters. We also
discuss the implications of this revised mass for studies of atmospheric
escape.

exoplanets, exoplanetary atmospheres

††journal: AAS Journals††facilities: LBT (PEPSI), TNG (HARPS-N), Fred L.
Whipple Observatory (TRES)††software: scipy (Virtanen et al., 2020),
petitRADTRANS (Mollière et al., 2019), SME (Valenti & Piskunov, 1996, 2012),
Time Utilities (Eastman, 2012), emcee (Foreman-Mackey et al., 2013), george
(https://george.readthedocs.io/en/latest/)

##  1 Introduction

Absolute mass is a critical but elusive quantity throughout the field of
observational astronomy. Most empirical constraints on mass rely on the
analysis of dynamical information of bodies interacting gravitationally, which
are not always attainable from two-dimensional projections on the plane of the
sky. Stellar eclipsing double-lined spectroscopic binaries (SB2s) are a
classic case in which dynamical masses can be directly measured. Spectroscopic
observations of two eclipsing stars, thus stars of known inclination, yield
orbital velocities by harnessing the Doppler effect for light. Newton’s
fundamental law of gravitation can be combined with knowledge of the system’s
orbital motion to recover the dynamical masses of both stars purely
empirically. These measurements calibrate the stellar evolution models that
are typically used to determine the masses of stars and indirectly the planets
they host (Popper, 1980; Harmanec, 1988; Andersen, 1991; Torres et al., 2010;
Stevens et al., 2018).

Mass is no less essential for the characterization of transiting planet
systems. For example, determining the interior composition of terrestrial
planets or measuring atmospheric escape on highly irradiated planets both
require knowledge of the planet’s mass. Transiting planets and their host
stars are in an eclipsing orbital configuration that, in principle, enables a
direct measurement of their masses. Transiting hot Jupiters (HJs), in
particular, are particularly suitable for this task. Due to their relatively
large mass and proximity to their host stars, HJs impart a stronger reflex
motion through their gravitational influence on their host stars than other
classes of planets. Thus they are the most accessible class of planets for
transit and radial velocity (RV) techniques. The RV signature of the host star
in principle can be constrained spectroscopically for the most massive and
tightly bound HJs, even without high-precision radial velocity instruments.

However, because of their extremely small planet/star flux ratios, most HJs
systems are essentially eclipsing single-lined spectroscopic binaries (SB1s).
Such systems do not allow for a unique solution to the masses and radii of the
individual objects (Seager & Mallén-Ornelas, 2003); rather there is a one-
dimensional degeneracy between the mass and radius of the primary. Thus any
inferences about the mass of a transiting planet requires an external
constraint on stellar mass, typically from stellar evolution models or
empirical scaling relations between the properties of main-sequence stars
(Torres et al., 2010; Dotter, 2016; Choi et al., 2016; Duck et al., ).
However, the reliability of the estimates and uncertainties of these semi-
empirical measurements has recently been called into question (Tayar et al.,
2020). Systematic uncertainties in stellar properties that may deviate from
the representative population used to calibrate evolutionary tracks can
propagate into an incorrect estimation of planet mass and other planetary
parameters. Furthermore, Tayar et al. 2020 show that current planetary
parameter uncertainties are often disproportionately underestimated relative
to the uncertainties on the host star properties from which the planetary
parameters are derived. Thus, obtaining masses by purely empirical methods is
important because it provides a concrete check on the semi-empirical practices
that are commonly employed.

If the radial velocity of transiting planet can be also measured, then the
system essentially becomes an eclipsing SB2, and thus an empirical measurement
of the system parameters, including the masses of the planet and star, can be
made. In this work, we apply the classical techniques used to analyze double-
lined spectroscopic binaries to the KELT-9 system (Gaudi et al., 2017), using
observations from the The Potsdam Echelle Polarimetric and Spectroscopic
Instrument (PEPSI) spectrograph (Strassmeier et al., 2015) on the Large
Binocular Telescope and the HARPS-N spectrograph on Telescopio Nazionale
Galileo (TNG). KELT-9 b is an ultra-hot Jupiter (UHJ) and the hottest known
planet to-date ($T_{\mathrm{eq}}=4050$ K). UHJs generally have orbital periods
below 3 days, yielding orbital velocities in the range of hundreds of km s-1;
this is approximately 3 orders of magnitude greater than the orbital velocity
of the host star.

Furthermore, the planet’s radial velocity can be measured from its atmospheric
absorption signature observed during transit. Snellen et al. 2010 is a
pioneering work in this field, not only for providing the first measurement of
winds on an exoplanet, but also for using this technique to empirically
constrain the absolute masses of the planet and its host. CO was adopted as a
tracer of HD 209458 b’s absorption signature to simultaneously chart the
planet’s orbital motion and measure its day-to-nightside winds. As a result,
the Doppler velocity shift of the planet’s atmospheric absorption features
over the course of a transit span a broader range of radial velocities (RVs)
than that of the relatively stationary stellar spectral features. Transmission
spectroscopy harnesses this distinction in velocity space to disentangle the
planet’s spectrum from its host star, allowing for a self-consistent
measurement of the planet’s absolute mass.

Yan & Henning 2018 adopt this method to measure the masses of KELT-9 b and its
host ($M_{\mathrm{p}}=3.23\pm 0.94\ \mathrm{M_{\mathrm{J}}}$,
$M_{\star}=3.00\pm 0.21\ \mathrm{M_{\odot}}$). However, they do not obtain an
original measurement of $K_{\star}$, the stellar orbital velocity; instead
they adopt the value given in Gaudi et al. 2017 measured from TRES stellar
spectra. In our work, we utilize the original TRES data in Gaudi et al. 2017
as well as additional RVs from higher precision spectrographs, PEPSI (LBT) and
HARPS-N (TNG), to constrain the dynamical masses of the KELT-9 system through
a fully self-contained analysis. In combination with the Hipparcos parallax
and SED fitting in Gaudi et al. 2017 that allow us to constrain KELT-9 b’s
radius, these observations will provide a purely empirical measurement of the
system parameters.

The KELT-9 system is notably overconstrained. A key focus of this work is to
emphasize that complementary observations (discussed in Section 5.2) will
contribute complete model-independent solutions to the system parameters. All
of these constraints can be combined to determine the system parameters to
higher accuracy. Improved precision of planetary parameters is critical for
science cases pertinent to KELT-9 b, such as atmospheric escape. A tight
constraint on mass will enable a more precise understanding of KELT-9 b’s
atmospheric escape since the mass loss rate measurement is dependent on
planetary parameters, i.e. mass or surface gravity.

Like Snellen et al. 2010, we use a fully self-consistent technique analogous
to the study of double-lined eclipsing binaries to obtain the absolute masses
of the KELT-9 system. In section 2, we describe our new datasets from the
PEPSI spectrograph as well as the archival HARPS-N observations we use in our
analysis. In section 3, we outline the techniques we use to recover stellar
and planetary orbital velocities and subsequently determine the absolute
masses of the planet and star. In section 4, we report our resulting mass
constraints and compare our findings with previous literature. In section 5,
we discuss the ways in which the properties of the KELT-9 system are
potentially empirically over-constrained. We consider contributions from
current and future complementary observations. When combined, we argue that
these will lead to a complete solution of system parameters. We also discuss
the implications for atmospheric escape. Finally, we present our conclusions
in 6.

##  2 Observations

Table 1: Datasets used to measure stellar RV: columns provide the name of the
instrument, the date of the starting observation, the start and end times of
observing, and the number of observations.

Instrument | Night |  $t_{\mathrm{start}}$ (UTC) |  $t_{\mathrm{end}}$ (UTC) | $N_{\mathrm{obs}}$
---|---|---|---|---
PEPSI | 2018-07-03 | 04:07:34.4 | 11:16:41.3 | 82
PEPSI | 2019-06-22 | 05:19:33.7 | 11:29:06.5 | 65
PEPSI | 2021-06-28 | 05:55:41.3 | 11:48:23.5 | 62†
HARPS-N | 2017-07-31 | 20:59:04.0 | 05:19:10.0 | 49
HARPS-N | 2018-07-20 | 21:20:24.0 | 05:09:58.0 | 46
TRES | 2014-2016 |  |  | 60†

†Out-of-transit only.

We observed two transits of KELT-9 b with the high-resolution échelle
spectrograph PEPSI (Strassmeier et al., 2015) on the Large Binocular Telescope
(LBT; two 8.4-m mirrors, effective aperture of 11.8 m) in Arizona (see Table 1
for the specific nights of observations). PEPSI has a blue arm (nominally
3830–5440 Å) and a red arm (nominally 5440–9070 Å) with six cross-dispersers
for full optical coverage. In this work, we use high-resolution data from the
blue arm taken with cross-disperser 3 ($\sim$4750–5430 Å, R=50,000)
exclusively because of negligible telluric contamination (Cauley et al.,
2019). The PEPSI pipeline produces wavelength-calibrated 1D spectra of each
order which are then continuum-normalized, corrected for solar barycentric
motion, and stitched into a single 1D spectral vector.

Our PEPSI dataset taken on 2018-07-03 (hereafter PEPSI 2018) was originally
presented with an analysis of KELT-9 b’s Balmer and metal lines in Cauley et
al. 2019. In the blue arm, the spectra were taken with exposure times between
220 s and 387 s (depending on fluctuations in observing conditions) to
approximately maintain a constant signal-to-noise ratio (SNR) of 210 in the
continuum; in practice, the SNR ranged between 182 to 219 across observations.
Additionally, we present a new dataset taken on 2019-06-22 (hereafter PEPSI
2019) for which observations began during the transit. The exposure times of
the blue arm spectra were between 214 s and 270 s with continuum SNRs ranging
between 286 and 321.

We utilized two archival HARPS-North (HARPS-N) datasets from the Italian
Center for Astronomical Archives (IA2) Facility to increase the number of
stellar RV samples for our measurement of the host star’s orbital motion.
HARPS-N is a high resolution (R$\sim$115,000) optical spectrograph (wavelength
range between 3874 Å and 6909 Å ) on the Telescopio Nazionale Galileo (TNG) in
La Palma, Spain. The first night of HARPS-North observations are from
2017-07-31 (hereafter HARPS-N 2017) and the other are from 2018-07-20
(hereafter HARPS-N 2018); both were originally presented in Hoeijmakers et al.
2019. The SNR of a given observation in these datasets ranges between 35-140,
depending on the order, at exposure times of 600 s. We retrieved the 1D,
order-stitched spectra from the IA2 Archive Facility.

Three of the datasets (PEPSI 2018, HARPS-N 2017, HARPS-N 2018) include
observations taken immediately before, during, and immediately after transit.
PEPSI 2019 only includes observations during and immediately after transit. We
converted all observation timings from their respective timing systems (PEPSI
times are provided in both $\mathrm{JD_{UTC}}$ as well as
$\mathrm{HJD_{UTC}}$; HARPS-N times are provided in $\mathrm{MJD_{UTC}}$) to
$\mathrm{BJD_{TDB}}$ using the Time
Utilities111https://astroutils.astronomy.osu.edu/time/utc2bjd.html online
software tool (Eastman, 2012) to make them comparable with our revised
ephemerides of the KELT-9 system (see Table 2). This is a crucial step for
precision radial velocity measurements (Eastman et al., 2010), especially of
atmospheric dynamics. The $\mathrm{JD_{UTC}}$ to $\mathrm{BJD_{TDB}}$
conversion (see Eastman et al. 2010 for a detailed description of the
difference between $\mathrm{JD_{UTC}}$ and $\mathrm{BJD_{TDB}}$) is not
accounted for in previous literature, e.g. Cauley et al. 2019, and yields a
difference up to 4.4 minutes in our PEPSI datasets. Note that the discovery
paper (Gaudi et al., 2017) ephemeris, which is commonly adopted in KELT-9 b
literature, is also in the $\mathrm{BJD_{TDB}}$ timing system.

Table 2: KELT-9 system parameters: stellar and planetary parameters generally from the discovery paper as well as updated measurements of orbital configuration and ephemerides from this work. Parameter | Units | Symbol | Value | Source
---|---|---|---|---
Stellar parameters: |  |  |  |
Stellar mass | $M_{\odot}$ | $M_{\star}$ | $2.52^{+0.25}_{-0.20}$ | Gaudi et al. 2017
Stellar radius | $R_{\odot}$ | $R_{\star}$ | $2.362^{+0.075}_{-0.063}$ | Gaudi et al. 2017
Stellar density | g $\mathrm{cm^{-3}}$ | $\rho_{\star}$ | $0.2702\pm 0.0029$ | Gaudi et al. 2017
Effective temperature | K | $T_{\mathrm{eff}}$ | $10170\pm 450$ | Gaudi et al. 2017
Projected rotational velocity | km s-1 | $v\sin{i}$ | $111.4\pm 1.3$ | Gaudi et al. 2017
Planetary parameters: |  |  |  |
Planet mass | $M_{J}$ | $m_{\mathrm{p}}$ | $2.88\pm 0.84$ | Gaudi et al. 2017
Planet radius | $R_{J}$ | $R_{\mathrm{p}}$ | $1.891^{+0.061}_{-0.053}$ | Gaudi et al. 2017
Semi-major axis | AU | $a$ | $0.03462^{+0.00110}_{-0.00093}$ | Gaudi et al. 2017
Eccentricity |  | $\varepsilon$ | $0$ | Gaudi et al. 2017
Spin-orbit alignment | degree | $\lambda$ | $-84.8\pm 1.4$ | Gaudi et al. 2017
Orbital inclination | degree | $i_{\mathrm{orbit}}$ | $86.79\pm 0.25$ | Gaudi et al. 2017
Ephemeris: |  |  |  |
Mid-transit time | $\mathrm{BJD_{TDB}}$ | $T_{0}$ | 2458566.436560 $\pm$ 0.000048 | This work
Time of secondary eclipse | $\mathrm{BJD_{TDB}}$ | $T_{\mathrm{S}}$ | 2458584.950546 $\pm$ 0.000048 | This work
Orbital period | days | $P$ | 1.48111890 $\pm$ 0.00000016 | This work
Ingress/egress transit duration | days | $\tau$ | $0.012808^{+0.000027}_{-0.000026}$ | This work
Total transit duration | days | $T_{14}$ | $0.15949\pm 0.00011$ | This work

##  3 Methods

###  3.1 Stellar orbital properties

To measure the dynamical mass of KELT-9 b, we need to know the orbital
properties of the planet and its host star. We first recover the out-of-
transit stellar velocities by least-squares deconvolution (LSD) of the stellar
spectra across our datasets to fit for the stellar orbital velocity and
systemic radial velocity as measured by each instrument.

####  3.1.1 Least-squares deconvolution of stellar line profiles

In subsequent analysis, $K_{\star}$ (stellar RV semi-amplitude) and
$v_{\mathrm{sys},\ i}$ (systemic velocity measured by instrument $i$) are
crucial values. The systemic radial velocity of the KELT-9 system has been a
point of controversy in previous literature (Gaudi et al., 2017; Hoeijmakers
et al., 2019; Borsa et al., 2019); see Table 3 for a compilation of the
different literature values. The standard procedure to recover stellar
velocities is to centroid the cross-correlation function (CCF) profiles of the
observations, where the CCF is the observed stellar spectra cross-correlated
with a template spectrum corresponding to the star’s effective temperature,
with a Gaussian fit. This method becomes increasingly imprecise for fast
rotators due to significant rotational broadening.

Table 3: KELT-9 revised system parameters Parameter |  Units |  Value |  Source
---|---|---|---
$v_{\mathrm{sys,PEPSI}}$ |  km s-1 |  $-17.86\pm 0.044$ |  This work
$v_{\mathrm{sys,HARPS-N}}$ |  km s-1 |  $-17.15\pm 0.11$ |  This work
|  |  $-17.74\pm 0.11$ |  Hoeijmakers et al. 2019
|  |  $-19.819\pm 0.024$ |  Borsa et al. 2019
$v_{sys,\mathrm{TRES}}$ |  km s-1 |  $-18.97\pm 0.12$ |  This work
|  |  $-20.567\pm 0.1$ |  Gaudi et al. 2017
$K_{\star}$ |  km s-1 |  $0.23\pm 0.060$ |  This work
|  |  $0.276\pm 0.079$ |  Gaudi et al. 2017
|  |  $0.293\pm 0.032$ |  Borsa et al. 2019
$K_{\mathrm{p}}$ |  km s-1 |  $239.07^{+5.83}_{-5.79}$ |  This work
|  |  $268.7^{+6.2}_{-6.4}$ |  Yan & Henning 2018
|  |  $234.24\pm 0.90$ |  Hoeijmakers et al. 2019
|  |  $241.5^{+3}_{-2}$ |  Pino et al. 2020
$M_{\star}$ |  $\mathrm{M_{\odot}}$ |  $2.11\pm 0.78$ |  This work
|  |  $2.52^{+0.25}_{-0.20}$ |  Gaudi et al. 2017
|  |  $3.00\pm 0.21$ |  Yan & Henning 2018
|  |  $1.978\pm 0.023$ |  Hoeijmakers et al. 2019
$m_{\mathrm{p}}$ |  $\mathrm{M_{J}}$ |  $2.17\pm 0.56$ |  This work
|  |  $2.88\pm 0.84$ |  Gaudi et al. 2017
|  |  $3.23\pm 0.94$ |  Yan & Henning 2018
|  |  $2.44\pm 0.70$ |  Hoeijmakers et al. 2019
$\rho_{\star}$ |  g $\mathrm{cm^{-3}}$ |  $0.29\pm 0.17$ |  This work
|  |  $0.2702\pm 0.0029$ |  Gaudi et al. 2017

Since KELT-9 is a rapidly rotating A0 star ($v\sin{i}=111.4$ km s-1), we
recover our own measurements of $K_{\star}$ and $v_{\mathrm{sys},i}$ by
performing LSD on our time-resolved stellar spectra to recover the rotational
broadening kernel (with in-transit observations affected by the Rossiter-
McLaughlin effect) of the star at each time of observation. The LSD procedure
allows for tunable regularization of the recovered line profile, a feature
which cross-correlation does not accomplish. The rotational kernel is centered
on the star’s radial velocity at the time of observation, which we can fit for
using an analytical kernel defined according to the star’s $v\sin{i}$. These
velocities over time can be fit with the orbital RV equation to determine
$K_{\star}$ and $v_{\mathrm{sys},i}$. This procedure was previously adopted in
Borsa et al. 2019.

For the first step, we refer to the LSD procedure provided in Kochukhov et al.
2010 with the modification for regularization given in Wang et al. 2017. In
our application of this technique, $\mathbf{Y}^{0}$ is a logarithmically-
sampled, $n$-element vector of the observed rotationally broadened stellar
spectrum and $\mathbf{F}$ is a template stellar spectrum of the corresponding
$T_{\mathrm{eff}}$. $\mathbf{F}$ has the same logarithmic wavelength sampling
as the observed spectrum but is not rotationally broadened. We generate
$\mathbf{F}$ by inputting the VALD3 linelist for a $T_{\mathrm{eff}}=10170$ K
star into the IDL software Spectroscopy Made Easy (SME) at 21 different limb-
darkening angles (Valenti & Piskunov, 1996, 2012). We treat the stellar disk
as a pixelated grid of $0.01R_{\star}\times 0.01R_{\star}$ cells to: 1)
interpolate in limb-darkening angle to generate the corresponding spectrum for
each cell, 2) add up the spectra of each cell, and 3) continuum-normalize to
generate the disk-integrated stellar spectrum $\mathbf{F}$.

The deconvolution process recovers the rotational broadening kernel that
transforms $\mathbf{F}$ into $\mathbf{Y}^{0}$. It can be described through
matrix multiplication of 1) the cross-correlation between a line mask,
$\mathbf{M}$, and the observed spectrum and 2) the inverse of the
autocorrelation of the line mask modified by regularization. This is
mathematically represented by the following equation (analogous to equation 1
in Wang et al. 2017 assuming homoscedasticity):

$\mathbf{Z}(v_{i})=(\mathbf{M}^{\mathrm{T}}\cdot\mathbf{M}+\Lambda\mathbf{R})^{-1}\cdot\mathbf{M}^{\mathrm{T}}\cdot\mathbf{Y}^{0}$
(1)

where $\mathbf{Z}(v_{i})$ is the m-element vector of the deconvolved line
profile (our output of interest). $\mathbf{M}$ is an $n\times m$ line mask
constructed by expanding the template spectrum into its corresponding Toeplitz
matrix; we adopt the definition provided in Donati et al. 1997. For
regularization, $\Lambda$ is a regularization parameter, and $\mathbf{R}$ is
the $m\times m$ matrix of first-order Tikhonov regularization; we adopt
Equation 16 in Donatelli & Reichel 2014 for the form of $\mathbf{R}$. The
velocities $v_{i}$ corresponding to the elements in the line profile
$\mathbf{Z}(v_{i})$ are determined by the wavelength shift of each row in the
$\mathbf{M}$ matrix relative to the template spectrum and should be linearly-
spaced since we constructed a template spectrum that is logarithmically-spaced
in wavelength.

We apply this deconvolution algorithm to all of our out-of-transit
observations since the in-transit kernels contain a deficiency in the
rotationally-broadened line profile due to the Rossiter-McLaughlin effect. As
the planet transits, it blocks a region of the stellar disk and the
contribution from this portion of the stellar disk is not added to the
integrated stellar disk spectrum. This manifests as a deficiency in the
rotational broadening kernel at the radial velocity of that portion of the
stellar disk based on stellar rotation. We pooled together our PEPSI data with
publicly available HARPS-N data and TRES data from the discovery paper (Gaudi
et al., 2017) to increase our sample size and out-of-transit coverage (see
Table 1).

(a)

(b)

Figure 1: (a) Example of least-squares deconvolution of a stellar line profile
with an analytic rotational kernel fit to centroid stellar radial velocity.
(b) Stellar radial velocity curve of the KELT-9 system derived from fitting a
circular orbital solution to the out-of-transit stellar radial velocities. The
data points are shifted by the best-fit systemic velocities of their
corresponding datasets. The inset in the top panel zooms in on the range of
velocities along the y-axis spanned by the best-fit orbital solution.

####  3.1.2 Analytic rotational kernel fit to stellar line profiles

Upon generating the observational line profiles, we fit each one with an
analytic rotational kernel model to determine their centers and recover an
orbital RV curve of the star. We use the analytical expression for a
rotational broadening kernel given in Gray 2005 (Equation 18.14), which is
adapted for our application as follows:

$\mathbf{G}(v_{i},t)=c_{1}\Big{[}1-\Big{(}\frac{v_{i}-v_{\star}(t)}{v\sin{i}}\Big{)}^{2}\Big{]}^{1/2}+c_{2}\Big{[}1-\Big{(}\frac{v_{i}-v_{\star}(t)}{v\sin{i}}\Big{)}^{2}\Big{]}$
(2)

where $\mathbf{G}(v_{i})$ is the analytical rotational kernel defined over a
range of velocities $v_{i}$; $v_{\star}(t)$ is the centroid of the rotational
kernel and represents the radial velocity of the star at the time $t$ of a
given observation; $v\sin{i}$ is the width of the kernel and represents the
sky-projected rotational velocity of the star (we used the value from Gaudi et
al. 2017 of 111.4 km s-1 to ensure the same value was used across
observations; the residuals in the left panel of Figure 1 suggest this value
sufficiently matches the data); and $c_{1}$ and $c_{2}$ are constants defined
in terms of $v\sin{i}$ and the linear limb-darkening coefficient of the star,
$\epsilon$, as

$c_{1}=\frac{2(1-\epsilon)}{\pi v\sin{i}(1-\epsilon/3)}$ (3)

and

$c_{2}=\frac{\epsilon}{2v\sin{i}(1-\epsilon/3)}$ (4)

For each observation, we apply least-squares fitting from scipy
222https://www.scipy.org/index.html (Virtanen et al., 2020) to determine the
best-fit values of $v_{\star}$, $\epsilon$, a multiplicative scaling factor
(rescales the analytical kernel to match the scaling of the empirical
deconvolved line profiles), and an additive offset (the empirical deconvolved
line profiles may not have a baseline centered at 0 due to a lack of flux
conservation between the template and observed spectra). Of these parameters,
we are most interested in $v_{\star}$, which provides measurements of the
star’s radial velocity over time. We estimated radial velocity errors by
bootstrapping the residuals of the flat region of the deconvolved kernel
($|$velocity$|\ >v\sin{i}$ in the left panel of Figure 1), adding the samples
to the best-fit model kernel, and refitting the line profile. We repeated this
bootstrapping procedure to obtain 1000 resampled values of $v_{\star}$ for
each observation and took the 16 percentile and 84 percentile samples as the
$\pm 1\sigma$ errors on the best-fit $v_{\star}$.

We also tested fixing the limb-darkening parameter to $\epsilon=0.3356$ based
on the stellar parameters that best describe KELT-9 in the limb-darkening
tables provided by Claret 2017. Upon refitting the rotational broadening
profiles, we find that the differences in the recovered stellar RVs when
$\epsilon$ is fixed are not significant enough to affect the resulting stellar
orbital parameters (discussed in Section 3.1.3) within errors.

####  3.1.3 Stellar orbital solution

To obtain the stellar orbital parameters, we fit a circular orbital RV
solution to our time-resolved measurements of stellar velocity, which has the
following form:

$v_{\star}(t)=K_{\star}\sin\Big{(}\frac{t-t_{0}}{p}\Big{)}+v_{\mathrm{sys},i}$
(5)

where $v_{\star}(t)$ is the stellar radial velocity as a function of time,
$K_{\star}$ is the stellar RV semi-amplitude (stellar orbital velocity for a
circular orbit), $t_{0}$ is the mid-transit time, $P$ is the orbital period,
and $v_{\mathrm{sys},i}$ is the systemic velocity of a given instrument. We
apply the Markov Chain Monte Carlo (MCMC) method to sample the parameter space
using the emcee333https://emcee.readthedocs.io/en/stable/ code (Foreman-Mackey
et al., 2013). We fit three common parameters that apply across all datasets:
mid-transit time ephemeris ($t_{0}$), period ($p$), and stellar RV semi-
amplitude ($K_{\mathrm{p}}$). We also fit three individual systemic velocities
that singly apply to their corresponding datasets: $v_{\mathrm{sys,PEPSI}}$,
$v_{\mathrm{sys,HARPS-N}}$, and $v_{\mathrm{sys,TRES}}$. Assuming a fixed
$v_{\mathrm{sys},i}$ across all datasets for a given instrument is valid since
all three instruments are wavelength-calibrated against a ThAr reference which
guarantees instrumental RV stability.

Table 4: Stellar orbital solution MCMC priors. Parameter | Units | Prior
---|---|---
$t_{0}$ | $\mathrm{BJD_{TDB}}$ | Gaussian prior
|  | Center: 2458566.436560
|  | Amplitude: 0.000048
$p$ | days | Gaussian prior
|  | Center: 1.48111890
|  | Amplitude: 0.00000016
$K_{\star}$ | km $\mathrm{s}^{-1}$ | Linearly uniform prior
|  | Lower bound: -1
|  | Upper bound: 0
$v_{\mathrm{sys,PEPSI}}$ | km $\mathrm{s}^{-1}$ | Linearly uniform prior
|  | Lower bound: -100
|  | Upper bound: 0
$v_{\mathrm{sys,HARPS-N}}$ | km $\mathrm{s}^{-1}$ | Linearly uniform prior
|  | Lower bound: -100
|  | Upper bound: 0
$v_{\mathrm{sys,TRES}}$ | km $\mathrm{s}^{-1}$ | Linearly uniform prior
|  | Lower bound: -50
|  | Upper bound: 50

We apply linearly uniform priors on all model parameters except $t_{0}$ and
$p$ for which we use Gaussian priors centered on the mid-transit time
ephemeris and period respectively as provided in Table 2 ($T_{0}$ and $P$)
with their corresponding errors as the standard deviations of the Gaussians.
See Table 4 for priors. The priors $T_{0}$ and $P$ are updated ephemerides for
the KELT-9 b system refitted with EXOFASTv2 (Eastman et al., 2019) to include
recent TESS observations (Ricker et al., 2015) in addition to the follow up
lightcurves and TRES RVs from the discovery paper. The updated ephemerides are
critical for this analysis since the observations span multiple years; the
error propagation using the original ephemerides can be as large as 4.8
minutes for the PEPSI 2019 dataset. Introducing the TESS 2019 data lightcurves
to the global fit improves the precision of the ephemerides and eliminates the
issue of ”stale” ephemerides by spanning a broad temporal baseline.

We evaluate the goodness-of-fit for a given model using the log-likelihood
$\ln{\cal L}\equiv-\frac{\chi^{2}}{2}$. Here $\chi^{2}$ is the statistical
chi-squared, defined as
$\chi^{2}=\sum\Big{(}\frac{\mathrm{data}-\mathrm{model}}{\mathrm{data\
errors}}\Big{)}^{2}$, of the model. To ensure the log-prior is on a comparable
scale as the log-likelihood, we scale the priors by the number of elements in
the observed stellar RV curve before summing the log-prior and log-likelihood
in the log-posterior. We run the emcee sampler with 10 walkers until the chain
length is at least 100 times the estimated autocorrelation time and the
estimated autocorrelation time has changed by less than 1%, checking every 500
steps. Under these criteria, the posterior distributions of all parameters
appear sufficiently Gaussian or converged. See the right panel of Figure 1 for
the observed and best-fit model stellar RV curve. Table 3 reports the fitted
systemic velocities and stellar RV semi-amplitude from this analysis.

###  3.2 Planet orbital properties

To obtain the orbital properties of the planet, we generate transmission
spectra which feature the planet’s atmospheric absorption track during transit
as well as secondary effects due to the geometry of the planet’s transit
across a non-uniform stellar disk. We simultaneously model both of these
signatures and apply a Bayesian framework for fitting the data with MCMC; the
most relevant outputs of this procedure is the orbital velocity of the planet,
$K_{\mathrm{p}}$, which we will use later in this work to constrain the
dynamical mass of the planet.

####  3.2.1 Transmission spectrum construction

We constrain the RV semi-amplitude of the planet, $K_{\mathrm{p}}$, and global
day-to-nightside winds, $v_{\mathrm{wind}}$, from a line-by-line analysis of
the atmospheric absorption signature seen in our transmission spectra of
KELT-9 b from PEPSI. To extract the planet’s atmospheric absorption signature,
we chose to focus on 6 Fe II lines in the wavelength range between 4915–5400 Å
to avoid telluric contamination and broad features such as the H$\beta$ line
around 4861.4 Å. The outputs from the PEPSI pipeline are continuum-normalized
stellar spectra, which we interpolate onto the same logarithmically-spaced
wavelength grid constructed such that absorption features are as well-sampled
as they are by the original wavelength grid. This enables uniform-spacing in
velocity across all observations. We constructed an empirical combined stellar
spectrum by taking the average of the out-of-transit spectra, which were
identified using updated ephemerides of the KELT-9 b system (see Table 2). To
remove the stellar component, we divide all spectra by the combined stellar
spectrum. The residual transmission spectra contain the planet’s signature as
well as secondary geometric effects caused by the Rossiter-McLaughlin Effect
(RME) (Rossiter, 1924; McLaughlin, 1924; Queloz et al., 2000; Ohta et al.,
2005; Gaudi & Winn, 2007) and center-to-limb variation (CLV) (Stenflo, 2015;
Yan et al., 2015; Yan & Henning, 2018), which are both consequences of non-
uniformity across the stellar disk (see the Doppler Shadow in Figure 2, a
byproduct of RME). These features are distinguishable because they span
different regions of velocity space and contribute opposing signs to the flux
map.

Upon inspection of our transmission spectra in the wavelength range between
4915–5400 Å for planet atmospheric signatures, we noticed that Fe II lines
presented the strongest atomic absorption signatures (apart from Balmer lines,
which are significantly broader as well). Previous studies of KELT-9 b’s
atmosphere with transmission spectroscopy revealed a diverse array of Balmer
lines ranging from H$\alpha$ to H$\zeta$ (Yan & Henning, 2018; Cauley et al.,
2019; Wyttenbach et al., 2020). Due to the extreme heating of KELT-9 b’s
atmosphere, traces of heavy metals, both neutral and ionized, have also been
found either by directly investigating the atomic absorption lines (Cauley et
al., 2019) or by harnessing the cross-correlation technique to boost the
absorption signal (Hoeijmakers et al., 2019). Like Cauley et al. 2019, we find
that the strongest metal features in our selected wavelength range are
produced by Fe II; although signatures of Fe I, Ti II, and Mg I can be
observed without cross-correlation in our PEPSI datasets, we have chosen to
focus on the Fe II lines to ensure that our measurement of $K_{\mathrm{p}}$ is
reliably derived from narrow absorption features that have the highest signal-
to-noise.

For each of our six selected Fe II lines, we focus on a tightly restricted
span of wavelengths that completely encapsulates the planet absorption and
Doppler shadow signatures for the observations that are fully in-transit
(between 2nd and 3rd contact). Thus we have a flux map of the fully in-transit
observations for each Fe II line. Wavelength is on the x-axis and orbital
phase is on the y-axis. The values in the map correspond to the continuum-
normalized transmission spectrum fluxes during transit.

We then generate a template spectrum of the Fe II species the planet’s
atmosphere using
petitRADTRANS444https://petitradtrans.readthedocs.io/en/latest/ (Mollière et
al., 2019), a radiative transfer code for modelling transmission and emission
spectra of planetary atmospheres. We note that at the time of writing, the Fe
II opacities used by petitRADTRANS are in air, unlike most other species which
are given in vacuum; hence we did not need to perform any wavelength
correction to match the template spectrum wavelengths with our observations.
The planetary parameters we adopt for $R_{\mathrm{p}}$, $R_{\star}$,
$m_{\mathrm{p}}$, $\log{g_{\mathrm{p}}}$, and $T_{\mathrm{eq}}$ are provided
in Table 2. Additional parameters that are necessary for generating a
transmission spectrum in petitRADTRANS are available in Table 5. Abundances
were chosen to be close to solar (Palme et al., 2014), while the pressure
structure was constructed to encompass the region of the atmosphere that most
strongly affects line formation in transmission spectra. It is not essential
to accurately model the specific shape of lines in the resulting template
spectrum of the planet’s atmosphere, as the line centers are the most relevant
quantity for our analysis of the planet’s orbital velocity through the Doppler
effect.

With the wavelengths of the Fe II line centers from the template spectrum, we
can convert our flux maps from wavelength space to velocity space with the
simplified Doppler effect for non-relativistic speeds:

$v\approx c\frac{\lambda-\lambda_{0}}{\lambda_{0}}$ (6)

where $\lambda_{0}$ is the line center corresponding to the Fe II line that
generates the absorption feature in a given flux map. We shift our flux maps
to the rest-frame of the star-planet system by subtracting the systemic
velocity calculated in Section 3.1.3 from the velocity grid of each flux map.
The resulting flux maps of each Fe II line is presented in Figure 3.

Figure 2: 2D map of transmission spectra over the course of KELT-9 b’s transit for the PEPSI 2018 dataset; the blue track is formed by the planet’s atmospheric absorption while the red track is the Doppler shadow from the Rossiter-McLaughlin effect.Top panel displays fully in-transit observations. Middle panel shows the best-fit model from MCMC sampling, with the Doppler shadow and CLV determined from numerical modelling of the planet’s transit using SME stellar models, while the planet absorption track is a uniform Gaussian signal shifted in velocity according to the best-fit orbital motion of the planet, systemic velocity, and best-fit day-to-nightside winds. The bottom panel shows the residuals (data - model). Table 5: Additional petitRADTRANS parameter inputs Parameter | Units | Value
---|---|---
Equilibrium temperature | K | 4050
Internal temperature | K | 100
Pressure range | bar |  $10^{-10}$ – $10^{2}$
Reference pressure | bar | $10^{-9}$
Infrared atmospheric opacity |  | 0.01
Ratio between optical and IR opacity |  | 0.4
Abundances: |  |
Fe II |  | 0.00125
$\mathrm{H_{2}}$ |  | 0.748
He |  | 0.250

####  3.2.2 Secondary effect modelling

Figure 3: Expanded version of Figure 2 displaying the six Fe II absorption
lines chosen for fitting $K_{\mathrm{p}}$ (PEPSI 2018 dataset).

We follow the numerical approach to modelling RME and CLV presented in
Casasayas-Barris et al. 2020 (originally presented in Yan et al. 2015 and Yan
et al. 2017, modified for planet radius interpolation in Casasayas-Barris et
al. 2020). We obtain linelists corresponding to the properties of the KELT-9
host star from VALD3 (Pakhomov et al., 2017). Known stellar properties from
Gaudi et al. 2017 are input into the IDL software Spectroscopy Made Easy (SME)
to model the host star spectrum using the VALD3 (Pakhomov et al., 2017)
linelists at 21 different limb-darkening angles (Valenti & Piskunov, 1996,
2012). We model the planet’s transit across a pixelated grid of
$0.01R_{\star}\times 0.01R_{\star}$ cells that make up the stellar disk. The
sky-projected position of the planet is computed according to Equations 7 and
8, which are analogous to Equations 7, 8, and 10 in Collier Cameron et al.
2010 (with the sign error in their Equation 10 corrected; this error was
originally noted in Eastman et al. 2019):

$\displaystyle x$
$\displaystyle=x_{\mathrm{p}}\cos{\lambda}+y_{\mathrm{p}}\sin{\lambda}$ (7)
$\displaystyle y$
$\displaystyle=-x_{\mathrm{p}}\sin{\lambda}+y_{\mathrm{p}}\cos{\lambda}$ (8)

where

$\displaystyle x_{\mathrm{p}}$ $\displaystyle=a\sin{\frac{2\pi(t-t_{0})}{P}}$
(9) $\displaystyle y_{\mathrm{p}}$
$\displaystyle=a\cos{\frac{2\pi(t-t_{0})}{P}}\cos{i_{\mathrm{orbit}}}$ (10)

In Equations 7–10, $\lambda$ is the sky-projected obliquity, $a$ is the semi-
major axis of the planet’s orbit, $t_{0}$ is the mid-transit time, $P$ is the
orbital period of the planet, $i_{\mathrm{orbit}}$ is the orbital inclination
of the planet, and $t$ is the time of a given observation; the system
parameters related to these equations can be found in Table 2. We generate a
grid of orbital phases that Nyquist samples the phases of the observations and
for each phase in the grid, we identify which cells are unobscured by the
planet’s transit. We interpolate in limb-darkening angle to generate a
spectrum for each unobscured cell and shift the spectrum in velocity according
to the cell’s radial velocity as determined by its location on the stellar
disk and the $v\sin{i}$ of the star. Then we add up the velocity-shifted
spectra of each unobscured cell, continuum-normalize, and divide by the out-
of-transit model spectrum (the sum of the spectra of all cells interpolated in
limb-darkening angle and shifted in radial velocity, then continuum-
normalized) to generate the disk-integrated transmission spectrum (excluding
the planet’s absorption signature) for each phase. In this manner, we produce
a 2D map of spectra over different orbital phases that encapsulate both the
Rossiter-McLaughlin effect and the center-to-limb variation that perturb the
transmission spectra during the planet’s transit. We generate a grid of such
models for a range of planet radii spanning
$0.7R_{\mathrm{p}}<R<2.5R_{\mathrm{p}}$ as done in Casasayas-Barris et al.
2020.

####  3.2.3 Line-by-line analysis of Fe II lines

We apply a procedure analogous to the Markov chain Monte Carlo (MCMC) fitting
of the observed transmission spectra presented in Casasayas-Barris et al. 2020
using emcee. To avoid additional asymmetries and velocity offsets from
equatorial jets and rotation, we exclude observations during ingress and
egress and instead focus on fully in-transit observations such that any
velocity shifts from atmospheric dynamics can be treated as a constant
blueshifted offset from day-to-nightside winds. We fit for the following free
parameters: the effective planet radius factor in RME/CLV modelling ($f$), the
RV semi-amplitude of the planet ($K_{\mathrm{p}}$), the radial velocity of the
terminator-averaged atmosphere in the planet’s restframe
($v_{\mathrm{wind}}$), the contrast of the planet’s Gaussian atmospheric
absorption profile ($h$), the standard deviation of the planet’s Gaussian
atmospheric absorption profile ($\sigma$), and the mid-transit time ($t_{0}$).
See Figure 4 for a sample corner plot.

We generate model flux maps according to the input parameters as follows. The
parameter $f$ corresponds to a model map of the RM and CLV effects from the
grid of models generated in Section 3.2.2 by interpolating between models in
planet radius according to the value of $f$, then interpolating the resulting
map in wavelength (then converted to velocity according to the reference Fe II
line of the data map being fitted and applying Equation 6) and orbital phase
to match the observations, i.e. the single line data maps generated in Section
3.2.1.

Recall that $t_{0}$ is a free parameter in this fitting procedure (this is not
done in Casasayas-Barris et al. 2020 or Yan et al. 2017; we have added this
extra free parameter so that we can include the uncertainty of mid-transit
time, which is strongly correlated with the parameter $v_{\mathrm{wind}}$, in
our analysis. We inject a Gaussian signal in our model maps that varies in
velocity with time but remains constant in amplitude and width to represent
the planet’s absorption signature. The parameters $K_{\mathrm{p}}$,
$v_{\mathrm{wind}}$, and $t_{0}$ determine the center of the planet’s
atmospheric absorption signal as a function of orbital phase (which is, once
again, dependent on the free parameter $t_{0}$) in our models as follows:

$v_{\mathrm{p}}(\phi)=K_{\mathrm{p}}\sin(2\pi\phi)+v_{\mathrm{wind}}$ (11)

The strength and width of the planet’s atmospheric absorption signal is
dependent on the free parameters $h$ and $\sigma$, resulting in the following
form for the planet’s absorption signature as a function of the radial
velocities $v$ and orbital phases $\phi$ that correspond with the observed
flux maps:

$T_{v}(\phi)=1+h\mathrm{e}^{\frac{v-v_{\mathrm{p}}(\phi)}{2\sigma}}$ (12)

This is analogous to Equation 1 in Yan & Henning 2018.

We adopt a Bayesian framework for sampling the parameter space with MCMC. We
apply linearly uniform priors on all model parameters except $t_{0}$, for
which we use a Gaussian prior centered on the mid-transit time value provided
in Table 2 ($T_{0}$, propagated to the epoch of the corresponding dataset
being fitted) with a standard deviation that matches the propagated mid-
transit time error. See Table 6 for priors. As before, we use
$-\frac{\chi^{2}}{2}$ as the log-likelihood of a given model and scale the
priors by the number of elements in the observed flux map (which should be the
same as the number of elements in the model flux map) before summing the log-
prior and log-likelihood in the log-posterior. We run the emcee sampler with
10 walkers until the the estimated autocorrelation time is 1.5% of the chain
length and the estimated autocorrelation time has changed by less than 1% ,
checking every 500 steps. Under these criteria, the posterior distributions of
all parameters appear sufficiently Gaussian or converged.

Table 6: Planet orbital solution MCMC priors (same for both line-by-line and cross-correlation analysis). Parameter | Units |  Prior
---|---|---
$f$ |  |  Uniform prior
|  |  Lower bound: 0.7
|  |  Upper bound: 2.5
$K_{\mathrm{p}}$ | km $\mathrm{s}^{-1}$ |  Uniform prior
|  |  Lower bound: 100
|  |  Upper bound: 350
$v_{\mathrm{wind}}$ | km $\mathrm{s}^{-1}$ |  Uniform prior
|  |  Lower bound: -50
|  |  Upper bound: 50
$h$ |  |  Uniform prior
|  |  Lower bound: -1
|  |  Upper bound: 0
$\sigma$ | velocity pixels |  Uniform prior
|  |  Lower bound: 0
|  |  Upper bound: 1000
$t_{0}$ | $\mathrm{BJD_{TDB}}$ |  Gaussian prior
|  |  Center: 2458566.436560 propagated to epoch of dataset
|  |  Width: 0.000048 propagated (according to error on orbital period) to epoch of dataset

The main goal of the model-fitting for the purposes of this work is to recover
a value for $K_{\mathrm{p}}$; $v_{\mathrm{wind}}$, or day-to-nightside winds,
will be the topic of a future work. The last 25% of samples are extracted from
the 10 walkers for a given Fe II line and aggregated into one chain for each
Fe II line. We weight each sample from all six chains according to the signal-
to-noise ratio (S/N) of the Fe II line flux map it corresponds to. Signal of a
given flux map is the 50th percentile sample of $h$ in the flux map’s chain.
Noise is the standard deviation of the residuals, i.e. best-fit model
subtracted from the flux map, where the best-fit model is a model generated by
the 50th percentile parameters from the chain of the flux map. We combine the
weighted chains to determine a representative value for the quantity of
interest. We recover a representative value of $K_{\mathrm{p}}$ by taking the
S/N-weighted 50th percentile sample across all chains. The lower and upper
1$\sigma$ values are the S/N-weighted 16th percentile and 84th percentile
samples respectively. See Figure 7 for the best-fit $K_{\mathrm{p}}$ for each
line analyzed in both the PEPSI datasets as well as the representative value
(dashed black line) from combining the MCMC analysis of all the lines.
Literature values are provided for comparison.

Figure 4: Example of a corner plot from MCMC fitting of the Fe II $\lambda
4923.9$ Å flux map.

####  3.2.4 Planet orbital properties: Cross-correlation analysis of Fe II
lines

(a)

(b)

Figure 5: Cross-correlated flux maps before artifact correction in the top
panel and after artifact correction in the bottom panel for the (a) PEPSI 2018
dataset and (b) HARPS-N 2017 dataset.

Since the SNR for the HARPS-N data is significantly lower than the PEPSI data,
individual lines are noisier in the HARPS-N observations. Cross-correlation is
a technique that can boost the SNR by adding up contributions from matching
signals between an observed and template spectrum across multiple features.
For comparison with the publicly available archival HARPS-N observations of
KELT-9 b, we cross-correlate our transmission spectra with a template Fe II
template spectrum (we will call it $f(\lambda)$ and the observed flux map
$F(\lambda,\phi)$). We cross-correlate our PEPSI and HARPS-N transmission
spectra between 4950 Å to 5400 Å (avoiding the broad H$\mathrm{\beta}$
feature) with the Fe II template generated in petitRADTRANS from 3.2.1 over a
sufficient range of velocities to span the atmospheric absorption signature.
The cross-correlation function we adopt is defined as follows:

$CCF(v,\phi)=\sum_{\lambda_{i}=\lambda_{\mathrm{min}}}^{\lambda_{\mathrm{max}}}f(\lambda_{i},v)F(\lambda_{i},\phi)$
(13)

where $f(\lambda,v)$ is the template spectrum shifted by velocity $v$.
Likewise, we cross-correlate our grid of RME/CLV models from 3.2.2 with the
same Fe II template. Consequently, we obtain CCF maps analogous to Figure 2
that can be fit with the same MCMC procedure as given in 3.2.3. See Figure 6
for observed and model CCF maps and the corresponding residuals for all 4
transmission spectroscopy datasets.

Our cross-correlated flux maps display distinct features that are not
attributable to the planet’s atmospheric absorption track or the RME/CLV
perturbations. In particular, our cross-correlation procedure reproduces the
artifacts described in Hoeijmakers et al. 2019 as ”stellar pulsations” in the
HARPS-N datasets (see Figure 5). We model these artifacts as individual
Gaussian features which have time-variable amplitudes, standard deviations,
and centroids that gradually change across observations. We fit these features
using non-linear least-squares optimization and subtract the model from the
observed cross-correlation flux map. Since these artifacts overlap in
velocity-phase space with the atmospheric absorption track, we first subtract
a preliminary fit of the planet’s absorption and the Doppler shadow before
removing the artifacts.

We also adopt a correlated noise model in our MCMC fitting procedure of the
cross-correlated maps using the Gaussian process regression library
george555https://george.readthedocs.io/en/latest/ (see panels in the second
row of Figure 6). Without this added component to the model, the parameter
errors are notably underestimated. We assume a 2-dimensional Matérn-3/2
covariance kernel, resulting in one amplitude parameter and two length
parameters in addition to our original model parameters to fit with MCMC. This
additional step does not yield parameter posterior distributions that are
statistically significantly different for the line-by-line analysis. For our
cross-correlated data, however, including the Gaussian process to model the
covariance structure of the data accurately captures the uncertainty from
correlated noise.

###  3.3 Measuring planet and stellar masses

As in Snellen et al. 2010, we treat the system as a double-lined spectroscopic
binary to determine the planet and stellar masses. In previous sections, we
demonstrated that: 1) the absorption lines of the planet’s atmosphere during
transit capture the planet’s orbital velocity while 2) the out-of-transit
observations can be deconvolved with a template stellar spectrum to retrieve
the stellar orbital velocity. From conservation of momentum, we know that

$m_{\mathrm{p}}K_{\mathrm{p}}=M_{\star}K_{\star}$ (14)

where $m_{\mathrm{p}}$ and $M_{\star}$ are the masses and $K_{\mathrm{p}}$ and
$K_{\star}$ are the RV semi-amplitudes (or orbital velocity in the case of a
circular orbit) of the planet and star respectively. Kepler’s 3rd law of
orbital motion states

$P^{2}=\frac{4\pi^{2}}{G(M_{\star}+m_{\mathrm{p}})}(a_{\star}+a_{\mathrm{p}})^{3}$
(15)

Assuming KELT-9 b is on a circular orbit based on an estimate of its
circularization timescale Borsa et al. 2019, we can relate the semi-major axes
with observables we have previously measured, namely the RV semi-amplitudes of
the planet ($K_{\mathrm{p}}$) and star ($K_{\star}$), orbital inclination
($i$), and the orbital period ($P$):

$\displaystyle\frac{K_{\mathrm{p}}}{\sin{i}}$ $\displaystyle=\frac{2\pi
a_{\mathrm{p}}}{P}$ (16) $\displaystyle\frac{K_{\star}}{\sin{i}}$
$\displaystyle=\frac{2\pi a_{\star}}{P}$ (17)

where $a_{\mathrm{p}}$ and $a_{\star}$ are the orbital semi-major axes of the
planet and star respectively. We can rewrite $M_{\star}$ in terms of
$m_{\mathrm{p}}$ using Equation 14 and apply Equations 16 and 17 to recast the
semi-major axes in terms of RV semi-amplitudes and orbital period. This allows
us to write our quantity of interest, $m_{\mathrm{p}}$, in terms of purely
empirical (and predominantly spectroscopic) observables. Upon doing so, we
derive the following expression for planet mass:

$m_{\mathrm{p}}=\frac{PK_{\star}}{2\pi
GK_{\mathrm{p}}}\Big{(}\frac{K_{\mathrm{p}}}{\sin{i}}\Big{)}^{3}\Big{(}1+\frac{K_{\star}}{K_{\mathrm{p}}}\Big{)}^{2}$
(18)

Applying this result to conservation of momentum (Equation 14) yields the
following expression for stellar mass:

$M_{\star}=\frac{P}{2\pi
G}\Big{(}\frac{K_{\mathrm{p}}}{\sin{i}}\Big{)}^{3}\Big{(}1+\frac{K_{\star}}{K_{\mathrm{p}}}\Big{)}^{2}$
(19)

Uncertainties on $m_{\mathrm{p}}$ and $M_{\star}$ are derived according to
linear propagation of errors (see Appendix, Section A for the analytic
expressions).

Figure 6: Same as Figure 3 except displaying CCF maps for each dataset.

##  4 Results

###  4.1 Summary of findings

We present the revised planetary and stellar parameters resulting from our
analysis in Table 3. Most notable is our measurement of the planetary mass,
the first self-consistent, purely empirical measurement of KELT-9 b’s
dynamical mass, $2.17\pm 0.56\ \mathrm{M_{J}}$, by treating the system as a
double-lined spectroscopic binary. Another byproduct of this analysis is the
mass of the star, which we measure to be $2.11\pm 0.78\ \mathrm{M_{\odot}}$.

###  4.2 Comparison with previous literature

We now compare our measurement of the dynamical mass of KELT-9 b against
previous literature (see Table 3). Throughout the following discussion, we
define $1\sigma$ by adding in quadrature the errors from two studies of
comparison for a parameter. Our measurement of $m_{\mathrm{p}}=2.17\pm 0.56\
\mathrm{M_{J}}$ agrees with Gaudi et al. 2017 and Yan & Henning 2018, and
Hoeijmakers et al. 2019.

The discrepancy in mass measurements across the literature can be largely
attributed to the use of $K_{\star}$ given in Gaudi et al. 2017 and
independent measurements of $K_{\mathrm{p}}$ (see Figure 7). Yan & Henning
2018 uses $K_{\mathrm{p}}=268.7^{+6.2}_{-6.4}$ km s-1 from their analysis of
atmospheric dynamics using $\mathrm{H\alpha}$ as a tracer as well as
$K_{\star}=0.276\pm 0.079$ km s-1 from the discovery paper; however our
measurement of $K_{\star}$ is nearly a factor of 1.2 ($\sim 0.42\sigma$
difference) less than that in the discovery paper. Hoeijmakers et al. 2019
uses a similar value of $K_{\mathrm{p}}$ as ours that agrees within errors,
which they obtained from fitting Fe II lines as well. Their procedure slightly
differs from ours in regards to order of operations: instead of fitting their
cross-correlated planet absorption track across all fully in-transit
observations simultaneously for a given dataset, they fit each observation
with an independent Gaussian absorption feature and fit a single circular
orbital solution to the individual planet RVs they obtain across both of their
datasets. However, they use the planet mass from Gaudi et al. 2017 along with
their value of $K_{\mathrm{p}}$ to revise the stellar mass, but then update
the planet mass according to the value of $K_{\star}$ from Gaudi et al. 2017
and their revised stellar mass. Thus their analysis is not self-consistent.
Our analysis, on the other hand, is self-consistent, because we obtain
empirical values of $K_{\mathrm{p}}$ and $K_{\star}$ exclusively by applying
the same routines and wavelength-consistent template linelists on all the data
we present in this work.

Studies of KELT-9 b’s neutral iron emission such as Pino et al. 2020 and
Kasper et al. 2021 obtain complementary constraints on $K_{\mathrm{p}}$. In
principle, one might expect these dayside measurements to disagree with our
terminator measurements due to atmospheric circulation. However, the
measurements of $K_{\mathrm{p}}$ from both Pino et al. 2020 (displayed in
Figure 7) and Kasper et al. 2021 (not included in Figure 7 since this work did
not consolidate their multiple constraints from different nights of
observations into one measurement) agree within $1\sigma$. This strengthens
our claim that the planet’s orbital parameters can be procured from the
planet’s transmission absorption signature without secondary effects from
atmospheric processes like circulation or condensation.

We also find that our measurement of stellar mass agrees with Gaudi et al.
2017 and Hoeijmakers et al. 2019, but disagrees with Yan & Henning 2018 by
1.1$\sigma$.

Since stellar density is a direct observable from transit lightcurves (Seager
& Mallén-Ornelas, 2003), we can use our revised stellar mass measurement and
empirical measurements of stellar radius from literature to compute a stellar
density and compare with the stellar density derived from the lightcurve as
provided in Gaudi et al. 2017. The stellar radius provided in Gaudi et al.
2017 is directly constrained to be $R_{\star}=2.17\pm 0.33\
\mathrm{R_{\odot}}$ from the Hipparcos parallax, effective temperature,
bolometric flux of the star from integrating its spectral energy distribution
(SED), and interstellar extinction. Using our measurement of the host star’s
dynamical mass and $V=\frac{4}{3}\pi R_{\star}^{3}$ (assuming a spherical
star) with the Gaudi et al. 2017 measurement for $R_{\star}$, we obtain a
stellar density of $\rho_{\star}=0.29\pm 0.17\ \mathrm{g\ cm^{-3}}$, which is
consistent with direct measurement from the lightcurve in Gaudi et al. 2017 of
$\rho_{\star}=0.2702\pm 0.0029\ \mathrm{gcm^{-3}}$. This suggests that our
dynamical mass analysis agrees with independent metrics from previous studies.
Note that the error on our empirical measurement of stellar density has nearly
equivalent contributions from the error on the empirical radius ($\sim$55% of
the stellar density error) and our empirical mass error ($\sim$45% of the
stellar density error).

We note that the spherical star assumption is not necessarily valid for fast
rotators. Ahlers et al. 2020 goes an additional step to account for effects
from gravity-darkening in their measurement of stellar radius using TESS
transit lightcurves of KELT-9 b. Gravity-darkening is a phenomenon by which
effective temperature varies across a stellar surface due to a rapidly
rotating star’s oblateness perturbing the star’s hydrostatic equilibrium near
its equator (Ahlers et al., 2020). While Ahlers et al. 2020 do provide an
equatorial radius and oblateness parameter that can, in principle, be combined
with our mass measurement to recover a stellar density, we do not perform this
step for comparison because the stellar equatorial radius they derive assumes
a prior on the stellar mass based on (Gaudi et al., 2017) of
$M_{*}=2.52^{+0.25}_{-0.20}~{}M_{\odot}$, which differs from our measurement
of $M_{*}=2.11\pm 0.78\ \mathrm{M_{\odot}}$.

##  5 Discussion

###  5.1 Caveats to the orbital motion model

The possibility of an eccentric orbit would undermine our measurement of
$K_{\mathrm{p}}$ and by extension the KELT-9 b’s mass. Our model presumes that
KELT-9 b and its host star obey circular orbital motion based on the planet’s
estimated circularization timescale Borsa et al. 2019. Furthermore, we obtain
an empirical constraint of $e\cos{\omega}=4.77\times 10^{-7}\pm 7.2\times
10^{-5}$ from the ephemerides of the primary transit and secondary eclipse;
the precision of this constraint is possible due to the inclusion of TESS
observations in our global fit for the system ephemerides. With such a low
value for $e\cos{\omega}$, it is unlikely that the planet’s orbit is
significantly eccentric enough to affect our empirical measurement of the
planet’s mass within errors.

We also assume that the dominant effect of atmospheric dynamics on the
planet’s absorption signature between 2nd and 3rd contact is a constant RV
offset due to day-to-nightside winds. Ehrenreich et al. 2020 has shown that
the wind component of WASP-76 b’s atmospheric absorption signature changes
from 0 km s-1 to $\sim$11 km s-1 over the course of the planet’s transit. This
would imply that our model of the planet’s atmospheric absorption signature
may not sufficiently capture the effect of KELT-9 b’s atmospheric dynamics.
However, we note that unlike Ehrenreich et al. 2020, our analysis does not
include observations taken during ingress or egress in the procedure for
fitting the orbital motion of the planet. Ingress and egress are expected to
contribute the strongest asymmetries in the absorption signal due to eastward
equatorial jets and rotation, yielding an additional redshift on the leading
limb during ingress and a blueshift on the trailing limb during egress
(assuming the planet’s rotational axis and orbital axis are aligned);
Ehrenreich et al. 2020 demonstrated this effect empirically. Since we exclude
observations that capture a partial limb of the planet, these asymmetries
should not manifest as strongly in our analysis.

Even if we were to include the entire transit in our analysis, the absorption
signal during ingress and egress is not very strong in our observations. When
we align all in-transit observations between 1st and 4th contact in the
planet’s restframe and additionally remove the measured offset component for
each dataset, the absorption feature is visibly centered around 0 km s-1 over
the course of the entire transit. This suggests that the wind speed does not
appear to vary significantly over the course of the transit. Furthermore, the
asymmetry observed by Ehrenreich et al. 2020 shows a constant RV offset during
the second half of the transit, which is where KELT-9 b’s absorption signal is
strongest and will contribute the most to our measured day-to-nightside wind
measurement.

###  5.2 Complementary constraints on system parameters

Among its numerous significant attributes, one noteworthy characteristic of
the KELT-9 system is that its physical properties are empirically over-
constrained. In particular, additional empirical constraints are possible due
to the gravity-darkening signature during the primary transit that originates
from rotational flattening of the host star. The geometry of such systems
introduces mathematical complexities that make quantitative analyses
analytically intractable. Rather, numerical methods are required to obtain the
system parameters; these are beyond the scope of this work. Therefore, we
proceed to qualitatively describe the ways in which a complete, model-
independent solutions to the system parameters can be achieved.

The first method treats the system as a double-lined eclipsing binary (SB2) as
done in this work. We have determined the masses of the system from
spectroscopic observables in this work by treating the KELT-9 system as an
eclipsing double-lined spectroscopic binary; thus all the orbital elements
($K_{\star}$, $K_{\mathrm{p}}$, $a_{\star}$, $a_{\mathrm{p}}$, etc.) of the
system are known. Furthermore, one can apply the assumption that the orbit is
circular to get the stellar and planetary semi-major axes $a_{\star}$ and
$a_{\mathrm{p}}$ respectively. Determining the radii of the host star and
planet is complicated somewhat by the fact the star is oblate and gravity-
darkened. However, the photometric gravity-darkening signature during the
primary transit (as presented in Ahlers et al. 2020) uniquely constrains the
inclination of the stellar rotational axis as well as the stellar oblateness
parameter. These two parameters then completely specify the geometry of the
planet’s transit across the oblate stellar disk, taking into account the
planet’s orbital inclination (angular offset along the line-of-sight between
the plane of the planet’s orbit and the plane of the sky) and stellar
inclination (angular offset along the line-of-sight between the stellar
rotation axis and the plane of the sky), and spin-orbit obliquity (angular
offset in the plane-of-sky between the stellar rotation axis and the planet’s
orbital axis). The chord transited by the planet is related to the equatorial
stellar radius by these angular quantities and the oblateness of the star.
From this, the equatorial stellar radius can be related to the orbital semi-
major axis by the FWHM transit duration and orbital period. The planet’s
radius is consequently obtained from the ratio of planet radius to stellar
equatorial radius provided in Ahlers et al. 2020. Thus a complete mass-radius
solution to the KELT-9 system can be obtained.

A second approach is to use just the host star RVs in combination with an
empirical stellar radius from an SED and parallax. The latter part of this
approach is provided in Gaudi & Winn 2007 using KELT-9’s Hipparcos parallax,
yielding $R_{\star}=2.17\pm 0.33\ \mathrm{R_{\odot}}$. This empirical radius
is more appropriately thought of as an effective radius since it assumes a
spherical star with a uniformly illuminated sky-projected disk. This empirical
radius can be translated into a measurement of the true stellar equatorial
radius by determining the equatorial radius (given the oblateness parameter
from Ahlers 2016) that yields an equivalent integrated disk flux as the
spherical case. Numerically integrating the oblate star disk flux as done in
Ahlers 2016 requires accounting for: 1) the star’s longitudinal variation in
flux from limb-darkening (Equation 3 of Ahlers 2016), 2) latitudinal variation
in flux from gravity-darkening (Equation 2 of Ahlers 2016), and 3) inclination
of the stellar rotational axis, which affects both the sky-projected disk
shape of the oblate star and the range of latitudes visible along the line-of-
sight (relevant for the latitude-dependence of the gravity-darkening profile).
As before, the planet’s radius is deduced from the ratio of planet radius to
stellar equatorial radius given in Ahlers et al. 2020. The stellar mass can be
determined from the stellar radial parameters in combination with stellar
density, a transit observable as given in Seager & Mallén-Ornelas 2003. The
relation for stellar density will require modifications for the geometry of an
oblate and inclined spheroid star transited by a planet on an inclined and
oblique circular orbit. These nuances directly effect the relations between
transit observables, i.e. depth and duration, and the system parameters of
interest, i.e. $a/R_{\star,\mathrm{eq}}$ and $\rho_{\star}$. Lastly, the
planet mass can be related to the stellar mass by the stellar RV semi-
amplitude as given in Equation 9 of Hoeijmakers et al. 2019 for the complete
solution to the system. A revised measurement of the empirical stellar radius
with a Gaia parallax would significantly improve the precision of this
constraint.

The third approach combines the single-lined spectroscopy measurement of
$K_{\star}$ with the photometric gravity-darkening signature during primary
transit, which provides constraints on stellar surface gravity. Since the host
star is an oblate spheroid, the surface gravity varies spatially as given by
Equation 10 of Barnes 2009. The transit signature maps stellar surface
brightness as the planet crosses the stellar disk. For a given $R_{\star}$ and
external constraint on $v\sin{i_{\star}}$, one can uniquely determine limb-
darkening coefficients, oblateness, stellar inclination, orbital inclination,
and spin-orbit obliquity as done in Ahlers et al. 2020. For a given stellar
equatorial radius, the system can be uniquely solved. Stellar mass is related
to stellar radius using the stellar density, a transit observable (taking into
account the nuances in the case of an oblate star as described above). Then as
before, planet radius can be obtained from the planet-to-star radius ratio
from the gravity-darkening transit signature and the planet mass comes from
Equation 9 of Hoeijmakers et al. 2019. The only missing component to solve the
system is the stellar radius. Ahlers et al. 2020 uses a gravity-darkening
correction on the SED from Gaudi et al. 2017, similar to our description in
the second approach. To distinguish this third approach from the second, we
note that stellar equatorial radius is related to the gravity-darkening
exponent $\beta$ and stellar inclination by plugging in Equation 10 of Barnes
2009 into Equation 9 of the same work. These quantities are not degenerate
because their effects are not scaled versions of each other. Thus, it is
possible to determine a unique solution for the equatorial radius based on the
shape of the transit lightcurve, and thus all the other geometric parameters
parameters (stellar oblateness, inclinations, spin-orbit obliquity, limb-
darkening coefficients, and the gravity-darkening exponent).

The fourth and final approach we present combines the single-lined
spectroscopy measurement of $K_{\star}$ with stellar surface gravity as
measured by the broadening of the stellar absorption lines seen in high-
resolution spectra. This measurement of $\log(g_{\star})$ does not account for
the fact that surface gravity varies spatially across the surface of an oblate
star. Instead it serves more as an effective surface gravity assuming a
spherical star. Under this assumption, Section 2.3.1 of Stevens et al. 2018
presents a derivation of system parameters from spectroscopic stellar RV semi-
amplitude and spectroscopic stellar surface gravity. This constraint is the
weakest because it does not account for the oblate geometry of the star.
Additionally, spectroscopically determining the surface gravity of hot, rapid
rotators is challenging since their absorption features are dominated by
rotational broadening.

We hope these four approaches guide future follow-up work to cross-validate
and tighten constraints on the KELT-9 system and other transiting two-body
systems in which the primary is a rapid rotator.

Figure 7: $K_{\mathrm{p}}$ measurements from Fe II line-by-line and cross-
correlation analysis with MCMC sampling errors.

###  5.3 Independent metrics of mass in the literature

With sufficient photometric precision, the BEaming, Ellipsoidal, and
Reflection (BEER) algorithm (Faigler & Mazeh, 2011) enables the detection of
short-period massive planets, both transiting and non-transiting. This
technique encompasses three distinct photometric signatures, of which two
(Doppler beaming and ellipsoidal variations) depend on the mass of the planet.
We provide order-of-magnitude estimates of these signals based on our revised
mass constraints.

Doppler beaming is the periodic variation in flux due to the orbital motion of
a star hosting a companion (Loeb & Gaudi, 2003). The fractional amplitude in
flux modulations is given in Loeb & Gaudi 2003 as

$\frac{\Delta F}{F_{0}}=\frac{3-\alpha_{\mathrm{beam}}}{c}K_{\star}$ (20)

where $\alpha_{\mathrm{beam}}$ is the power law exponent of the emitted flux
from the source star as a function of frequency. We can estimate
$\alpha_{\mathrm{beam}}$ of a blackbody source with effective temperature
$T_{\mathrm{eff}}$ according to Equation 3 in Loeb & Gaudi 2003 as
$\alpha_{\mathrm{beam}}(\nu)\approx\frac{\mathrm{e}^{x}(3-x)-3}{\mathrm{e}^{x}-1}$
where $x=\frac{h\nu}{kT_{\mathrm{eff}}}$. Adopting the frequency corresponding
to the $V$-band wavelength of 551 nm as a representative frequency, we obtain
a Doppler beaming signal on the order $\frac{\Delta F}{F_{0}}=2.17$ ppm. This
estimated signal falls far below the photometric precision of present-day
instruments like TESS to be discernable.

Ellipsoidal variations are flux perturbations that arise as a consequence of
tidal effects induced by a companion. Equation 2 of Faigler & Mazeh 2011
estimate ellipsoidal variations of amplitudes

$A_{\mathrm{ellip}}=\alpha_{\mathrm{ellip}}\frac{m_{\mathrm{p}}\sin{i_{\mathrm{orbit}}}}{M_{\star}}\Big{(}\frac{R_{\star}}{a}\Big{)}^{3}\sin{i_{\mathrm{orbit}}}$
(21)

where $\alpha_{\mathrm{ellip}}$ is defined in Equation 4 of Faigler & Mazeh
2011 as

$\alpha_{\mathrm{ellip}}\approx 0.15\frac{(15+u)(1+g)}{3-u}$ (22)

The parameter $u$ is the linear limb-darkening coefficient of the star, which
we take to be $u=0.3356$ based on Claret 2017. We adopt edge cases of $g=0.3$
and $g=1$ for the gravity-darkening coefficient (not to be mistaken with the
gravity-darkening exponent from Ahlers et al. 2020) as suggested in Faigler &
Mazeh 2011. The ratio $\frac{R_{\star}}{a}$ is a direct observable from
transit lightcurves; we use $\frac{a}{R_{\star}}=3.153$ from Gaudi et al.
2017. The resulting estimate of the ellipsoidal variations signal based on our
mass ratio ranges between 35.0 and 53.8 ppm.

Wong et al. 2020 presents photometric harmonics in the phase curve of the
KELT-9 system by phase-folding TESS lightcurves and binning at 30-minute
intervals. They measure periodic variations attributable to ellipsoidal
variations with an amplitude of $39.6\pm 4.5$ ppm. This measurement is within
the range of the ellipsoidal variation signal amplitudes we estimate based on
our revised masses of the system. In fact, it is marginally easier to explain
Wong et al. 2020’s measurement with our planet-star mass ratio than with the
original value in Gaudi et al. 2017, which yield ellipsoidal variation signal
amplitudes between 38.9 to 60.0 ppm for the same range of values for $\beta$.
Note that numerous uncertainties play into our estimate of ellipsoidal
variations, such as the uncertainty of the limb- and gravity-darkening
coefficients of this particular star. In particular, metallicity and
evolutionary effects on the limb-darkening coefficient are stronger after the
onset of convection ($\log{T_{\mathrm{eff}}}>3.9$) (Claret, 2017).
Furthermore, the star is significantly distorted by rotation and the planet’s
orbit is nearly orthogonal to this distortion. This geometry limits
gravitational distortion of the star by the companion planet and weakens the
signal from ellipsoidal variations.

###  5.4 Consequences for atmospheric escape

As the hottest planet known to-date, KELT-9 b is a unique target for studying
the hydrodynamic escape of planetary atmospheres driven by extreme stellar
irradiation (Fossati et al., 2018; García Muñoz & Schneider, 2019; Krenn et
al., 2021). Observations of atmospheric loss on KELT-9 b use light element
tracers such as Balmer lines (Yan & Henning, 2018; Cauley et al., 2019;
Wyttenbach et al., 2020) and Mg I (Cauley et al., 2019) to probe mass loss
rates according to empirically-validated species densities and a slew of
assumptions surrounding elemental abundances and equilibrium conditions.
Notably, these studies generally assume thermodynamic and ionization
equilibrium as well as solar abundances of species; the former may be
discounted by NLTE studies of KELT-9 b’s upper atmosphere (García Muñoz &
Schneider, 2019; Wyttenbach et al., 2020; Fossati et al., 2021), while the
latter may be invalidated by recent work on atmospheric retrieval of hot
Jupiters (Giacobbe et al., 2021).

We summarize the state of the field to-date regarding the mass loss rate of
KELT-9b. Gaudi et al. 2017 provided an initial estimate of KELT-9 b’s mass
loss rate between $10^{10}$ – $10^{13}$ g s-1 based on the equation for
energy-limited mass loss rate provided in Equation 22 of Murray-Clay et al.
2009, which scales inversely with the mass of the planet. Yan & Henning 2018
empirically constrain KELT-9 b’s mass loss rate by identifying an excess in
H$\alpha$ absorption depth relative to the photometric transit depth in the
continuum. By estimating 1) the number density of hydrogen from model fits to
the H$\alpha$ line profile (these models depend on planetary parameters, most
notably planet mass) and 2) the contribution from the altitude regime of the
planet where H$\alpha$ can energetically escape beyond the planet’s Roche
lobe, Yan & Henning 2018 estimate a mass loss rate of $\dot{M}\sim 10^{12}$ g
s-1. Cauley et al. 2019 corroborates this measurement within an order of
magnitude, estimating $\dot{M}\sim 1\times 10^{12}$ g s-1 when using Mg I as a
tracer and $\dot{M}\sim 3\times 10^{12}$ g s-1 from Balmer line analysis.
Wyttenbach et al. 2020 report $\dot{M}\sim 10^{12.8\pm 0.3}$ g s-1 from Balmer
line analysis as well.

Reconciling observations of KELT-9 b’s atmospheric escape with theory
currently faces unresolved challenges. HJ atmospheric escape is typically
modelled as hydrodynamic escape due to heating from Lyman continuum absorption
in the X-ray and extreme-ultraviolet (XUV). Fossati et al. 2018 notes that
ionizing XUV fluxes in the wavelength regime relevant to heating (and thus
escape) are weaker in hotter intermediate mass stars, such as the stellar host
of the KELT-9 system, than their cooler ($\sim$8000-8500 K) counterparts, such
as the host of another UHJ WASP-33 b. XUV flux is driven by coronal heating,
which is related to stellar magnetic activity. A convective envelope is
necessary for the interactions that generate magnetic activity. The surface
convective zone vanishes in hotter stars and consequently they emit less XUV
flux. Fossati et al. 2018 estimate KELT-9 b’s mass loss rate is $\sim$1010 \-
1011 by accounting for the heating efficiency of XUV flux from the KELT-9 host
star. Their estimated H$\mathrm{\alpha}$ transit depth of 0.7% does not agree
with the observed 1.8% in Yan & Henning 2018. They propose that one way of
bridging this gap is by adopting a planetary mass on the lower end of the
1$\sigma$ range provided in Gaudi et al. 2017 ($M_{\mathrm{p}}=2.88\pm 0.84\
\mathrm{M_{J}}$). As previously noted, $\dot{M}$ scales inversely with planet
mass in the most simplified energy-limited case, which balances gravitational
potential and thermal heating from the host star’s X-ray and extreme UV (XUV)
radiation. This would bring the planetary mass required to explain
observations of KELT-9 b’s atmospheric loss in closer agreement with our
revised empirical mass. We estimate that our revised mass increases the
estimated energy-limited mass loss rate (see Fossati et al. 2018 for a
discussion of this estimate for KELT-9 b and its nuances) by $\sim$33%. Note
that recent work by Krenn et al. 2021 comparing the energy-limited
approximation against hydrodynamic simulations of atmospheric escape shows
that the energy-limited approximation is not suited for estimating UHJ mass-
loss rates (or planets in extreme temperature or mass regimes in general; UHJs
satisfy both of these conditions).

García Muñoz & Schneider 2019 expands upon the implications of Fossati et al.
2018 by proposing that Balmer continuum absorption in the near-ultraviolet
(NUV) is the dominant source of heating. However, their models best match
observations of KELT-9 b’s H$\alpha$ absorption when adopting a planet mass
between 0.80–1.20 $M_{J}$. This is incompatible with our revised mass.

##  6 Conclusion

We have presented the analysis of spectroscopic data of KELT-9b in the context
of an eclipsing double-lined spectroscopic binary. This has enabled us to
direcly and empirically obtain the dynamical mass the KELT-9 b and its host
star. Using multi-epoch spectroscopic observations of the system, we find that
the dynamical mass of the planet is $m_{\mathrm{p}}=2.17\pm 0.56\
\mathrm{M_{J}}$ and the star is $M_{\star}=2.11\pm 0.78\ \mathrm{M_{\odot}}$.
Our planet mass measurement generally agrees with previous literature (Gaudi
et al., 2017; Hoeijmakers et al., 2019). We also obtain a purely empirical
measurement of stellar density, a direct observable from transit lightcurves,
that agrees with the value in the discovery paper; this suggests that our
analysis is trustworthy. Our methodology can be applied to many Hot Jupiter
systems, thereby enabling the direct and empirical measurement of their
planets and host stars.

We note that the KELT-9 system is empirically overconstrained due to the
unique geometric information provided by the in-transit gravity-darkening
signature of rapid rotator systems. We present a framework for obtaining a
complete solution to the system parameters in three purely observational ways.

An order-of-magnitude estimate shows that this revised planetary mass is large
enough to induce ellipsoidal variations observable with TESS phase curves of
the KELT-9 system. Furthermore, this result is especially crucial for studies
of atmospheric escape, which depend on mass for empirical measurement and
theoretical modelling; as of now, the high levels of atmospheric escape seen
on KELT-9 b have not been reconciled with its high mass. Our purely empirical
confirmation of the planet’s mass motivates further exploration of this
conundrum surrounding the substantial mass loss on KELT-9 b.

We thank the anonymous referee for taking the time to provide insightful
comments which improved the quality of this paper. A.P.A would like to thank
the David G. Price Fellowship in Astronomical Instrumentation for funding her
work this year. Work by B.S.G. and J.W. was partially supported by the Thomas
Jefferson Chair for Space Exploration endowment from the Ohio State
University. This work is based on observations made with the Large Binocular
Telescope. The LBT is an international collaboration among institutions in the
United States, Italy and Germany. LBT Corporation partners are: The University
of Arizona on behalf of the Arizona Board of Regents; Istituto Nazionale di
Astrofisica, Italy; LBT Beteiligungsgesellschaft, Germany, representing the
Max-Planck Society, The Leibniz Institute for Astrophysics Potsdam, and
Heidelberg University; The Ohio State University, representing OSU, University
of Notre Dame, University of Minnesota and University of Virginia. This paper
includes data collected by the TESS mission. Funding for the TESS mission is
provided by the NASA’s Science Mission Directorate. We thank Dr. George Zhou
for providing TRES observations of the KELT-9 system. We also thank Dr.
Marshall Johnson, Dr. Francesco Borsa, Dr. Jens Hoeijmakers, Dr. Lorenzo Pino,
and Dr. Luca Fossati for contributing their valuable expertise regarding high-
resolution spectroscopy and the KELT-9 system.

## References

  * Ahlers (2016) Ahlers, J. P. 2016, ApJ, 832, 93, doi: 10.3847/0004-637X/832/1/93
  * Ahlers et al. (2020) Ahlers, J. P., Johnson, M. C., Stassun, K. G., et al. 2020, AJ, 160, 4, doi: 10.3847/1538-3881/ab8fa3
  * Andersen (1991) Andersen, J. 1991, A&A Rev., 3, 91, doi: 10.1007/BF00873538
  * Barnes (2009) Barnes, J. W. 2009, ApJ, 705, 683, doi: 10.1088/0004-637X/705/1/683
  * Borsa et al. (2019) Borsa, F., Rainer, M., Bonomo, A. S., et al. 2019, A&A, 631, A34, doi: 10.1051/0004-6361/201935718
  * Casasayas-Barris et al. (2020) Casasayas-Barris, N., Pallé, E., Yan, F., et al. 2020, A&A, 640, C6, doi: 10.1051/0004-6361/201935623e
  * Cauley et al. (2019) Cauley, P. W., Shkolnik, E. L., Ilyin, I., et al. 2019, AJ, 157, 69, doi: 10.3847/1538-3881/aaf725
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Claret (2017) Claret, A. 2017, A&A, 600, A30, doi: 10.1051/0004-6361/201629705
  * Collier Cameron et al. (2010) Collier Cameron, A., Bruce, V. A., Miller, G. R. M., Triaud, A. H. M. J., & Queloz, D. 2010, MNRAS, 403, 151, doi: 10.1111/j.1365-2966.2009.16131.x
  * Donatelli & Reichel (2014) Donatelli, M., & Reichel, L. 2014, Journal of Computational and Applied Mathematics, 272, 334, doi: https://doi.org/10.1016/j.cam.2013.08.015
  * Donati et al. (1997) Donati, J. F., Semel, M., Carter, B. D., Rees, D. E., & Collier Cameron, A. 1997, MNRAS, 291, 658, doi: 10.1093/mnras/291.4.658
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8, doi: 10.3847/0067-0049/222/1/8
  * (14) Duck, A., Gaudi, B. S., Rodriguez, J., & Eastman, J.
  * Eastman (2012) Eastman, J. 2012, Time Utilities.  http://ascl.net/1206.012
  * Eastman et al. (2010) Eastman, J., Siverd, R., & Gaudi, B. S. 2010, PASP, 122, 935, doi: 10.1086/655938
  * Eastman et al. (2019) Eastman, J. D., Rodriguez, J. E., Agol, E., et al. 2019, arXiv e-prints, arXiv:1907.09480.  https://arxiv.org/abs/1907.09480
  * Ehrenreich et al. (2020) Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, Nature, 580, 597, doi: 10.1038/s41586-020-2107-1
  * Faigler & Mazeh (2011) Faigler, S., & Mazeh, T. 2011, MNRAS, 415, 3921, doi: 10.1111/j.1365-2966.2011.19011.x
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fossati et al. (2018) Fossati, L., Koskinen, T., Lothringer, J. D., et al. 2018, ApJ, 868, L30, doi: 10.3847/2041-8213/aaf0a5
  * Fossati et al. (2021) Fossati, L., Young, M. E., Shulyak, D., et al. 2021, arXiv e-prints, arXiv:2106.11263.  https://arxiv.org/abs/2106.11263
  * García Muñoz & Schneider (2019) García Muñoz, A., & Schneider, P. C. 2019, ApJ, 884, L43, doi: 10.3847/2041-8213/ab498d
  * Gaudi & Winn (2007) Gaudi, B. S., & Winn, J. N. 2007, ApJ, 655, 550, doi: 10.1086/509910
  * Gaudi et al. (2017) Gaudi, B. S., Stassun, K. G., Collins, K. A., et al. 2017, Nature, 546, 514, doi: 10.1038/nature22392
  * Giacobbe et al. (2021) Giacobbe, P., Brogi, M., Gandhi, S., et al. 2021, arXiv e-prints, arXiv:2104.03352.  https://arxiv.org/abs/2104.03352
  * Gray (2005) Gray, D. F. 2005, The Observation and Analysis of Stellar Photospheres, 3rd edn. (Cambridge University Press), doi: 10.1017/CBO9781316036570
  * Harmanec (1988) Harmanec, P. 1988, Bulletin of the Astronomical Institutes of Czechoslovakia, 39, 329
  * Hoeijmakers et al. (2019) Hoeijmakers, H. J., Ehrenreich, D., Kitzmann, D., et al. 2019, A&A, 627, A165, doi: 10.1051/0004-6361/201935089
  * Kasper et al. (2021) Kasper, D., Bean, J. L., Line, M. R., et al. 2021, ApJ, 921, L18, doi: 10.3847/2041-8213/ac30e1
  * Kochukhov et al. (2010) Kochukhov, O., Makaganiuk, V., & Piskunov, N. 2010, A&A, 524, A5, doi: 10.1051/0004-6361/201015429
  * Krenn et al. (2021) Krenn, A. F., Fossati, L., Kubyshkina, D., & Lammer, H. 2021, arXiv e-prints, arXiv:2105.05858.  https://arxiv.org/abs/2105.05858
  * Loeb & Gaudi (2003) Loeb, A., & Gaudi, B. S. 2003, ApJ, 588, L117, doi: 10.1086/375551
  * McLaughlin (1924) McLaughlin, D. B. 1924, ApJ, 60, 22, doi: 10.1086/142826
  * Mollière et al. (2019) Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67, doi: 10.1051/0004-6361/201935470
  * Murray-Clay et al. (2009) Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23, doi: 10.1088/0004-637X/693/1/23
  * Ohta et al. (2005) Ohta, Y., Taruya, A., & Suto, Y. 2005, ApJ, 622, 1118, doi: 10.1086/428344
  * Pakhomov et al. (2017) Pakhomov, Y., Piskunov, N., & Ryabchikova, T. 2017, in Astronomical Society of the Pacific Conference Series, Vol. 510, Stars: From Collapse to Collapse, ed. Y. Y. Balega, D. O. Kudryavtsev, I. I. Romanyuk, & I. A. Yakunin, 518.  https://arxiv.org/abs/1710.10854
  * Palme et al. (2014) Palme, H., Lodders, K., & Jones, A. 2014, Solar System Abundances of the Elements, Vol. 2 (Elsevier Ltd.), 15–36
  * Pino et al. (2020) Pino, L., Désert, J.-M., Brogi, M., et al. 2020, ApJ, 894, L27, doi: 10.3847/2041-8213/ab8c44
  * Popper (1980) Popper, D. M. 1980, ARA&A, 18, 115, doi: 10.1146/annurev.aa.18.090180.000555
  * Queloz et al. (2000) Queloz, D., Eggenberger, A., Mayor, M., et al. 2000, A&A, 359, L13.  https://arxiv.org/abs/astro-ph/0006213
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003, doi: 10.1117/1.JATIS.1.1.014003
  * Rossiter (1924) Rossiter, R. A. 1924, ApJ, 60, 15, doi: 10.1086/142825
  * Seager & Mallén-Ornelas (2003) Seager, S., & Mallén-Ornelas, G. 2003, ApJ, 585, 1038, doi: 10.1086/346105
  * Snellen et al. (2010) Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049, doi: 10.1038/nature09111
  * Stenflo (2015) Stenflo, J. O. 2015, A&A, 573, A74, doi: 10.1051/0004-6361/201424685
  * Stevens et al. (2018) Stevens, D. J., Gaudi, B. S., & Stassun, K. G. 2018, ApJ, 862, 53, doi: 10.3847/1538-4357/aaccf5
  * Strassmeier et al. (2015) Strassmeier, K. G., Ilyin, I., Järvinen, A., et al. 2015, Astronomische Nachrichten, 336, 324, doi: 10.1002/asna.201512172
  * Tayar et al. (2020) Tayar, J., Claytor, Z. R., Huber, D., & van Saders, J. 2020, arXiv e-prints, arXiv:2012.07957.  https://arxiv.org/abs/2012.07957
  * Torres et al. (2010) Torres, G., Andersen, J., & Giménez, A. 2010, A&A Rev., 18, 67, doi: 10.1007/s00159-009-0025-1
  * Valenti & Piskunov (1996) Valenti, J. A., & Piskunov, N. 1996, A&AS, 118, 595
  * Valenti & Piskunov (2012) —. 2012, SME: Spectroscopy Made Easy.  http://ascl.net/1202.013
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wang et al. (2017) Wang, J., Prato, L., & Mawet, D. 2017, ApJ, 838, 35, doi: 10.3847/1538-4357/aa6345
  * Wong et al. (2020) Wong, I., Shporer, A., Kitzmann, D., et al. 2020, AJ, 160, 88, doi: 10.3847/1538-3881/aba2cb
  * Wyttenbach et al. (2020) Wyttenbach, A., Mollière, P., Ehrenreich, D., et al. 2020, A&A, 638, A87, doi: 10.1051/0004-6361/201937316
  * Yan et al. (2015) Yan, F., Fosbury, R. A. E., Petr-Gotzens, M. G., Zhao, G., & Pallé, E. 2015, A&A, 574, A94, doi: 10.1051/0004-6361/201425220
  * Yan & Henning (2018) Yan, F., & Henning, T. 2018, Nature Astronomy, 2, 714, doi: 10.1038/s41550-018-0503-3
  * Yan et al. (2017) Yan, F., Pallé, E., Fosbury, R. A. E., Petr-Gotzens, M. G., & Henning, T. 2017, A&A, 603, A73, doi: 10.1051/0004-6361/201630144

##  Appendix A Error propagation of stellar and planetary mass from RV
observables

We apply linear propagation of errors to recover the following expression for
the propagated uncertainty of the planet mass:

$\begin{split}\delta m_{\mathrm{p}}=&\Bigg{[}\Big{(}\frac{\partial
m_{\mathrm{p}}}{\partial P}\Big{)}^{2}\delta P^{2}+\Big{(}\frac{\partial
m_{\mathrm{p}}}{\partial i}\Big{)}^{2}\delta i^{2}+\Big{(}\frac{\partial
m_{\mathrm{p}}}{\partial K_{\mathrm{p}}}\Big{)}^{2}\delta
K_{\mathrm{p}}^{2}+\Big{(}\frac{\partial m_{\mathrm{p}}}{\partial
K_{\star}}\Big{)}^{2}\delta K_{\star}^{2}\Bigg{]}^{1/2}\end{split}$

We present the analytic form of each partial derivative term in Equation A
below:

$\displaystyle\frac{\partial m_{\mathrm{p}}}{\partial P}$
$\displaystyle=\frac{K_{\star}(K_{\mathrm{p}}+K_{\star})^{2}}{2\pi
G\sin^{3}{i}}$ $\displaystyle\frac{\partial m_{\mathrm{p}}}{\partial i}$
$\displaystyle=\frac{-3K_{\star}(K_{\mathrm{p}}+K_{\star})^{2}P\cos{i}}{2\pi
G\sin^{2}{i}}$ $\displaystyle\frac{\partial m_{\mathrm{p}}}{\partial
K_{\mathrm{p}}}$
$\displaystyle=\frac{K_{\star}(K_{\mathrm{p}}+K_{\star})P}{\pi G\sin^{3}{i}}$
$\displaystyle\frac{\partial m_{\mathrm{p}}}{\partial K_{\star}}$
$\displaystyle=\frac{(K_{\mathrm{p}}+K_{\star})(K_{\mathrm{p}}+3K_{\star})P}{2\pi
G\sin^{3}{i}}$

Analogously, the propagated stellar mass uncertainty is:

$\begin{split}\delta M_{\star}=&\Bigg{[}\Big{(}\frac{\partial
M_{\star}}{\partial P}\Big{)}^{2}\delta P^{2}+\Big{(}\frac{\partial
M_{\star}}{\partial i}\Big{)}^{2}\delta i^{2}+\Big{(}\frac{\partial
M_{\star}}{\partial K_{\mathrm{p}}}\Big{)}^{2}\delta
K_{\mathrm{p}}^{2}+\Big{(}\frac{\partial M_{\star}}{\partial
K_{\star}}\Big{)}^{2}\delta K_{\star}^{2}\Bigg{]}^{1/2}\end{split}$

with partial derivative terms of the form

$\displaystyle\frac{\partial m_{\mathrm{p}}}{\partial P}$
$\displaystyle=\frac{K_{\mathrm{p}}(K_{\mathrm{p}}+K_{\star})^{2}}{2\pi
G\sin^{3}{i}}$ $\displaystyle\frac{\partial m_{\mathrm{p}}}{\partial i}$
$\displaystyle=\frac{-3K_{\mathrm{p}}(K_{\mathrm{p}}+K_{\star})^{2}P\cos{i}}{2\pi
G\sin^{2}{i}}$ $\displaystyle\frac{\partial m_{\mathrm{p}}}{\partial
K_{\mathrm{p}}}$
$\displaystyle=\frac{(K_{\mathrm{p}}+K_{\star})(3K_{\mathrm{p}}+K_{\star})P}{2\pi
G\sin^{3}{i}}$ $\displaystyle\frac{\partial m_{\mathrm{p}}}{\partial
K_{\star}}$
$\displaystyle=\frac{K_{\mathrm{p}}(K_{\mathrm{p}}+K_{\star})P}{\pi
G\sin^{3}{i}}$