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

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.


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. Spectro-scopic 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 re-quire 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 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 dayto-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 p = 3.23 ± 0.94 M J , M = 3.00 ± 0.21 M ). However, they do not obtain an original measurement of K , 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.

OBSERVATIONS
We observed two transits of KELT-9 b with the highresolutioné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 (∼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∼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 JD UTC as well as HJD UTC ; HARPS-N times are provided in MJD UTC ) to BJD TDB using the Time Utilities 1 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 JD UTC to BJD TDB conversion (see Eastman et al. 2010 for a detailed description of the difference between JD UTC and 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 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.

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.

Least-squares deconvolution of stellar line profiles
In subsequent analysis, K (stellar RV semi-amplitude) and v 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.
Since KELT-9 is a rapidly rotating A0 star (v sin i = 111.4 km s −1 ), we recover our own measurements of K and v 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 and v 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, Y 0 is a logarithmicallysampled, n-element vector of the observed rotationally broadened stellar spectrum and F is a template stellar spectrum of the corresponding T eff . F has the same logarithmic wavelength sampling as the observed spectrum but is not rotationally broadened. We generate F by inputting the VALD3 linelist for a T 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 × 0.01R 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 F.
The deconvolution process recovers the rotational broadening kernel that transforms F into Y 0 . It can be described through matrix multiplication of 1) the cross-correlation between a line mask, 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  Gaudi et al. 2017 homoscedasticity): where Z(v i ) is the m-element vector of the deconvolved line profile (our output of interest). M is an n × 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, Λ is a regularization parameter, and R is the m × m matrix of first-order Tikhonov regularization; we adopt Equation 16 in Donatelli & Reichel 2014 for the form of R. The velocities v i corresponding to the elements in the line profile Z(v i ) are determined by the wavelength shift of each row in the 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).

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: 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, , as and For each observation, we apply least-squares fitting from scipy 2 (Virtanen et al. 2020) to determine the best-fit values of v , , 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 , 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 for each observation and took the 16 th percentile and 84 th percentile samples as the ±1σ errors on the best-fit v .
We also tested fixing the limb-darkening parameter to = 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 is fixed are not significant enough to 2 https://www.scipy.org/index.html affect the resulting stellar orbital parameters (discussed in Section 3.1.3) within errors.

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: where v (t) is the stellar radial velocity as a function of time, K 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 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 emcee 3 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 p ). We also fit three individual systemic velocities that singly apply to their corresponding datasets: v sys,PEPSI , v sys,HARPS−N , and v sys,TRES . Assuming a fixed v 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. 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 L ≡ − χ 2 2 . Here χ 2 is the statistical chi-squared, defined as χ 2 = data−model data errors 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 bestfit model stellar RV curve. Table 3 reports the fitted systemic velocities and stellar RV semi-amplitude from this analysis.

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 p , which we will use later in this work to constrain the dynamical mass of the planet.

Transmission spectrum construction
We constrain the RV semi-amplitude of the planet, K p , and global day-to-nightside winds, v wind , from a line-byline 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β line around 4861.4 A. The outputs from the PEPSI pipeline are continuumnormalized 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α to Hζ (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 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 petitRADTRANS 4 4 https://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 p , R , m p , log g p , and T 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: where λ 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.

Secondary effect modelling
We follow the numerical approach to modelling RME and CLV presented in Casasayas-Barris et al. 2020 (orig- inally presented in Yan et al. 2015 andYan 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 × 0.01R cells that make up the stellar disk. x = x p cos λ + y p sin λ y = −x p sin λ + y p cos λ where In Equations 7-10, λ is the sky-projected obliquity, a is the semi-major axis of the planet's orbit, t 0 is the midtransit time, P is the orbital period of the planet, i 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 p < R < 2.5R p as done in Casasayas-Barris et al. 2020.

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 p ), the radial velocity of the terminator-averaged atmosphere in the planet's restframe (v wind ), the contrast of the planet's Gaussian atmospheric absorption profile (h), the standard deviation of the planet's Gaussian atmospheric absorption profile (σ), 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 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 p , v 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: The strength and width of the planet's atmospheric absorption signal is dependent on the free parameters h and σ, resulting in the following form for the planet's absorption signature as a function of the radial velocities v and orbital phases φ that correspond with the observed flux maps: 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 − χ 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.
The main goal of the model-fitting for the purposes of this work is to recover a value for K p ; v wind , or day-tonightside 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 p by taking the S/N-weighted 50th percentile sample across all chains. The lower and upper 1σ values are the S/N-weighted 16th percentile and 84th percentile samples respectively. See Figure 7 for the best-fit K 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. 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 (λ) and the observed flux map F (λ, φ)). We cross-correlate our PEPSI and HARPS-N transmission spectra between 4950Å to 5400Å (avoiding the broad Hβ 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: where f (λ, 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 leastsquares 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 george 5 (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 lineby-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 5 https://george.readthedocs.io/en/latest/ 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 where m p and M are the masses and K p and K 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 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 p ) and star (K ), orbital inclination (i), and the orbital period (P ): where a p and a are the orbital semi-major axes of the planet and star respectively. We can rewrite M in terms of m p using Equation 14 and apply Equations 16 and 17 to recast the semi-major axes in terms of RV semiamplitudes and orbital period. This allows us to write our quantity of interest, m p , in terms of purely empirical (and predominantly spectroscopic) observables. Upon doing so, we derive the following expression for planet mass: Applying this result to conservation of momentum (Equation 14) yields the following expression for stellar mass: Uncertainties on m p and M are derived according to linear propagation of errors (see Appendix, Section A for the analytic expressions).

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 selfconsistent, purely empirical measurement of KELT-9 b's dynamical mass, 2.17 ± 0.56 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 ± 0.78 M .

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σ by adding in quadrature the errors from two studies of comparison for a parameter. Our measurement of m p = 2.17 ± 0.56 M J agrees with Gaudi et al. 2017and 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 given in Gaudi et al. 2017 and independent measurements of K p (see Figure 7). Yan & Henning 2018 uses K p = 268.7 +6.2 −6.4 km s −1 from their analysis of atmospheric dynamics using Hα as a tracer as well as K = 0.276 ± 0.079 km s −1 from the discovery paper; however our measurement of K is nearly a factor of 1.2 (∼ 0.42σ difference) less than that in the discovery paper. Hoeijmakers et al. 2019 uses a similar value of K 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 crosscorrelated planet absorption track across all fully intransit 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 p to revise the stellar mass, but then update the planet mass according to the value of K 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 p and K 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 andKasper et al. 2021 obtain complementary constraints on K p . In principle, one might expect these dayside measurements to disagree with our terminator measurements due to atmospheric circulation. However, the measurements of K 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σ. 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 andHoeijmakers et al. 2019, but disagrees with Yan & Henning 2018 by 1.1σ.
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 = 2.17 ± 0.33 R 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 = 4 3 πR 3 (assuming a spherical star) with the Gaudi et al. 2017 measurement for R , we obtain a stellar density of ρ = 0.29 ± 0.17 g cm −3 , which is consistent with direct measurement from the lightcurve in Gaudi et al. 2017 of ρ = 0.2702 ± 0.0029 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 (∼55% of the stellar density error) and our empirical mass error (∼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 , which differs from our measurement of M * = 2.11 ± 0.78 M .

Caveats to the orbital motion model
The possibility of an eccentric orbit would undermine our measurement of K 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 ω = 4.77×10 −7 ±7.2×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 ω, 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 ∼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.

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 , K p , a , a 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 and a 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 = 2.17 ± 0.33 R . 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 limbdarkening (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 ,eq and ρ . 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 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 and external constraint on v sin i , 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 planetto-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 β 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 with stellar surface gravity as measured by the broadening of the stellar absorption lines seen in high-resolution spectra. This measurement of log(g ) 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.

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 where α beam is the power law exponent of the emitted flux from the source star as a function of frequency. We can estimate α beam of a blackbody source with effective temperature T eff according to Equation 3  where x = hν kT 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 ∆F F0 = 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.
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 R a is a direct observable from transit lightcurves; we use a R = 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 ± 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 β. 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 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.

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α 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α 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α can energetically escape beyond the planet's Roche lobe, Yan & Henning 2018 estimate a mass loss rate ofṀ ∼ 10 12 g s −1 . Cauley et al. 2019 corroborates this measurement within an order of magnitude, estimatingṀ ∼ 1 × 10 12 g s −1 when using Mg I as a tracer andṀ ∼ 3 × 10 12 g s −1 from Balmer line analysis. Wyttenbach et al. 2020 reportṀ ∼ 10 12.8±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 (∼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 ∼10 10 -10 11 by accounting for the heating efficiency of XUV flux from the KELT-9 host star. Their estimated Hα 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σ range provided in Gaudi et al. 2017 (M p = 2.88 ± 0.84 M J ). As previously noted, 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 ∼33%. Note that recent work by Krenn et al. 2021 comparing the energylimited 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α absorption when adopting a planet mass between 0.80-1.20 M J . This is incompatible with our revised mass.

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 p = 2.17 ± 0.56 M J and the star is M = 2.11 ± 0.78 M . 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.