# TOI-1518b: A Misaligned Ultra-hot Jupiter with Iron in its Atmosphere

Samuel H. C. Cabot Yale University, 52 Hillhouse Avenue, New Haven, CT 06511,
USA  Aaron Bello-Arufe National Space Institute, Technical University of
Denmark, Elektrovej, DK-2800 Kgs. Lyngby, Denmark  João M. Mendonça National
Space Institute, Technical University of Denmark, Elektrovej, DK-2800 Kgs.
Lyngby, Denmark  René Tronsgaard National Space Institute, Technical
University of Denmark, Elektrovej, DK-2800 Kgs. Lyngby, Denmark  Ian Wong 51
Pegasi b Fellow  Department of Earth, Atmospheric and Planetary Sciences,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  George Zhou
Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA  Lars A. Buchhave National Space Institute, Technical
University of Denmark, Elektrovej, DK-2800 Kgs. Lyngby, Denmark  Debra A.
Fischer Yale University, 52 Hillhouse Avenue, New Haven, CT 06511, USA  Keivan
G. Stassun Vanderbilt University, Department of Physics & Astronomy, 6301
Stevenson Center Ln., Nashville, TN 37235, USA  Victoria Antoci National Space
Institute, Technical University of Denmark, Elektrovej, DK-2800 Kgs. Lyngby,
Denmark  Stellar Astrophysics Centre, Department of Physics and Astronomy,
Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark  David Baker
Physics Department, Austin College, Sherman, TX 75090, USA  Alexander A.
Belinski Sternberg Astronomical Institute, M.V. Lomonosov Moscow State
University, 13, Universitetskij pr., 119234, Moscow, Russia  Björn Benneke
Department of Physics and Institute for Research on Exoplanets, Université de
Montréal, Montreal, QC, Canada  Luke G. Bouma Department of Astrophysical
Sciences, Princeton University, NJ 08544, USA  Jessie L. Christiansen NASA
Exoplanet Science Institute – Caltech/IPAC Pasadena, CA 91125 USA  Karen A.
Collins Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA  Maria V. Goliguzova  Sternberg Astronomical
Institute, M.V. Lomonosov Moscow State University, 13, Universitetskij pr.,
119234, Moscow, Russia  Simone Hagey University of Saskatchewan, Saskatchewan,
Canada  Jon M. Jenkins NASA Ames Research Center, Moffett Field, CA, 94035
Eric L. N. Jensen Dept. of Physics & Astronomy, Swarthmore College, Swarthmore
PA 19081, USA  Richard C. Kidwell Jr  Space Telescope Science Institute,
Baltimore, MD, USA  Didier Laloum  Société Astronomique de France, 3 Rue
Beethoven, 75016 Paris, France  Bob Massey Villa ’39 Observatory, Landers, CA
92285, USA  Kim K. McLeod Department of Astronomy, Wellesley College,
Wellesley, MA 02481, USA  David W. Latham Center for Astrophysics $|$ Harvard
& Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Edward H. Morgan
Department of Physics and Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  George Ricker
Department of Physics and Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Boris S.
Safonov Sternberg Astronomical Institute, M.V. Lomonosov Moscow State
University, 13, Universitetskij pr., 119234, Moscow, Russia  Joshua E.
Schlieder NASA Goddard Space Flight Center, 8800 Greenbelt Rd, Greenbelt, MD
20771, USA  Sara Seager Department of Physics and Kavli Institute for
Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Department of Earth, Atmospheric and Planetary
Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Department of Aeronautics and Astronautics, MIT, 77 Massachusetts Avenue,
Cambridge, MA 02139, USA  Avi Shporer Department of Physics and Kavli
Institute for Astrophysics and Space Research, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Jeffrey C. Smith NASA Ames Research
Center, Moffett Field, CA, 94035  SETI Institute, Mountain View, CA 94043, USA
Gregor Srdoc  Kotizarovci Observatory, Sarsoni 90, 51216 Viskovo, Croatia
Ivan A. Strakhov  Sternberg Astronomical Institute, M.V. Lomonosov Moscow
State University, 13, Universitetskij pr., 119234, Moscow, Russia  Guillermo
Torres Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA  Joseph D. Twicken NASA Ames Research Center, Moffett
Field, CA, 94035  SETI Institute, Mountain View, CA 94043, USA  Roland
Vanderspek Department of Physics and Kavli Institute for Astrophysics and
Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139,
USA  Michael Vezie  Department of Physics and Kavli Institute for Astrophysics
and Space Research, Massachusetts Institute of Technology, Cambridge, MA
02139, USA  Joshua N. Winn Department of Astrophysical Sciences, Princeton
University, NJ 08544, USA

###### Abstract

We present the discovery of TOI-1518b — an ultra-hot Jupiter orbiting a bright
star ($V=8.95$). The transiting planet is confirmed using high-resolution
optical transmission spectra from EXPRES. It is inflated, with $R_{p}=1.875\pm
0.053\,R_{\rm J}$, and exhibits several interesting properties, including a
misaligned orbit (${240.34^{+0.93}_{-0.98}}$ degrees) and nearly grazing
transit ($b=0.9036^{+0.0061}_{-0.0053}$). The planet orbits a fast-rotating F0
host star ($T_{\mathrm{eff}}\simeq 7300$ K) in 1.9 days and experiences
intense irradiation. Notably, the TESS data show a clear secondary eclipse
with a depth of $364\pm 28$ ppm and a significant phase curve signal, from
which we obtain a relative day–night planetary flux difference of roughly 320
ppm and a 5.2$\sigma$ detection of ellipsoidal distortion on the host star.
Prompted by recent detections of atomic and ionized species in ultra-hot
Jupiter atmospheres, we conduct an atmospheric cross-correlation analysis. We
detect neutral iron (${5.2\sigma}$), at ${K_{p}={157^{+68}_{-44}}}$ km s-1 and
${V_{\rm sys}={-16^{+2}_{-4}}}$ km s-1, adding another object to the small
sample of highly irradiated gas-giant planets with Fe detections in
transmission. Detections so far favor particularly inflated gas giants with
radii $\gtrsim 1.78$$R_{\rm J}$; although this may be due to observational
bias. With an equilibrium temperature of $T_{\rm eq}=2492\pm 38$ K and a
measured dayside brightness temperature of $3237\pm 59$ K (assuming zero
geometric albedo), TOI-1518b is a promising candidate for future emission
spectroscopy to probe for a thermal inversion.

Exoplanets — Hot Jupiters — Exoplanet atmospheres — Spectroscopy

††software: Tapir (Jensen, 2013), AstroImageJ (Collins et al., 2017), molecfit
(Smette et al., 2015), radvel (Fulton et al., 2018), Lightkurve (Lightkurve
Collaboration et al., 2018), Helios-K (Grimm et al., 2021), FastChem (Stock et
al., 2018), batman (Kreidberg, 2015), emcee (Foreman-Mackey et al., 2013), SME
(Valenti & Piskunov, 1996), exoplanet (Foreman-Mackey et al., 2021), astropy
(Astropy Collaboration et al., 2018)

##  1 Introduction

Transiting exoplanets — those that pass directly between their host stars and
an observer — offer a wealth of information about their systems. The transit
itself is detectable through the minuscule fraction of starlight occulted by
the planet, which is well within the sensitivity of many current ground- and
space-based telescopes. The Kepler (Borucki et al., 2010) and K2 (Howell et
al., 2014) missions together yielded thousands of transiting exoplanet
candidates, some of which are among the most notable and well-characterized to
date. Planets found by surveys such as HATNet (Bakos et al., 2004), KELT
(Pepper et al., 2007), and WASP (Pollacco et al., 2006) orbit some of the
brightest stars, and have hence been popular targets for atmospheric
characterization. Today, the frontier lies with the Transiting Exoplanet
Survey Satellite (TESS; Ricker et al. 2014), which is searching for planets
transiting bright stars across the entire sky.

The science drivers behind exoplanet transit observations are several-fold.
Newly discovered systems improve our baseline understanding of exoplanet
populations and distributions (Howard et al., 2012; Fressin et al., 2013;
Fulton et al., 2017), and how their properties may be linked to system
architecture and formation scenarios (Lissauer et al., 2011; Fabrycky et al.,
2014; Millholland et al., 2017; Weiss et al., 2018). The presence of
additional, non-transiting exoplanets can be inferred from transit timing
perturbations (Holman & Murray, 2005; Ballard et al., 2011). The host star’s
obliquity can be probed by the Rossiter-McLaughlin effect (Winn et al., 2010;
Triaud, 2018), which also relates to formation pathways (Dawson & Johnson,
2018, and references therein). Finally, transits enable the study of exoplanet
atmospheres based on the excess absorption of starlight from high-altitude
species (e.g. Seager & Sasselov, 1998; Charbonneau et al., 2002; Snellen et
al., 2010; Sing et al., 2016). These latter investigations require dedicated
spectroscopic followup.

In this study, we present the confirmation of the TESS transiting planet
candidate TOI-1518b, a highly irradiated gas giant planet possessing iron
vapor in its atmosphere. Of exoplanets discovered by TESS, this is the first
high-resolution detection of an atmospheric species. Several TESS candidates
have been confirmed as hot Jupiters so far, including HD 202772Ab (Wang et
al., 2019), HD 2685b (Jones et al., 2019), TOI-150b (Cañas et al., 2019), HD
271181b (Kossakowski et al., 2019), TOI-172b (Rodriguez et al., 2019),
TOI-564b, and TOI-905b (Davis et al., 2019). However, TOI-1518b is unique due
to its close-in orbit (1.9 day period) and high level of irradiation from its
F-type host star.

The new planet falls within the category of ultra-hot Jupiters (UHJs), which
have equilibrium temperatures exceeding 2000 K (Fortney et al., 2008;
Parmentier et al., 2018). Many UHJs contain vaporized metals, both neutral and
ionized, in their upper atmospheres (e.g. Hoeijmakers et al., 2018; Casasayas-
Barris et al., 2018). These metals and molecules containing them are
recognized as strong sources of opacity in the optical and near-ultraviolet
regions (Fortney et al., 2008; Lothringer et al., 2020). UHJs often exhibit
thermal inversions (Haynes et al., 2015; Evans et al., 2017); however, the
exact species responsible for the inversions are debated (Fortney et al.,
2008; Lothringer et al., 2018; Gandhi & Madhusudhan, 2019). High-resolution
spectroscopy has become a common method for detecting important species in UHJ
atmospheres and also serves as a means of probing winds (Louden & Wheatley,
2015; Casasayas-Barris et al., 2019) and extended atmospheres (Yan & Henning,
2018).

Our paper is organized as follows. In Section 2 we analyze the TESS photometry
of TOI-1518. We reproduce the detection of a planet candidate, obtain
constraints on its orbital parameters, and report a robust detection of the
secondary eclipse and phase-curve modulations. We also present high-resolution
spectroscopic observations of the system during transit. This spectroscopic
transit is analyzed in Section 3, from which we measure the Rossiter-
Mclaughlin effect and obtain further constraints on the orbit and host star.
Section 3 also includes a review of the cross-correlation method for
atmospheric characterization, the results of which are presented in Section 4.
Finally, we discuss TOI-1518b in the context of previously studied UHJs in
Section 5.

##  2 Observations and System Characterization

Figure 1: The normalized Presearch Data Conditioning Simple Aperture
Photometry (PDCSAP) light curve of TOI-1518 generated by the SPOC pipeline.
The scheduled momentum dumps are indicated by the vertical blue dashed lines.
The red points denote flux ramps and regions of severe systematics that were
trimmed prior to our light curve fits. The orbital phase curve modulations are
discernible in the raw photometry.

This section describes our analysis of the available TESS photometry of
TOI-1518, as well as high-resolution optical spectra of the system. We measure
the system parameters by simultaneously modeling the transit, secondary
eclipse, and full-orbit phase curve. We also fit spectral lines to determine
properties of the star. As detailed below, radial velocity (RV) measurements
of the system provide some broad constraints. However, the deduced parameters
have large uncertainties owing to the rapid rotation speed of the star.

###  2.1 TESS Photometry

The star TIC 427761355 (also designated as BD+66 1610) was observed by Camera
3 of the TESS instrument during Sectors 17 and 18 (UT 2019 Oct 7 to Nov 27).
The Quick Look Pipeline (QLP; Huang et al., 2020) detected a likely transit
signal in the photometry and flagged the companion as a candidate transiting
exoplanet with parameters characteristic of a close-in hot Jupiter. The system
was released as a TESS Object of Interest (TOI) with the designation TOI-1518.
The full-frame images (FFIs) were processed by the Science Processing
Operations Center (SPOC; Jenkins et al., 2016) and made publicly available on
the Mikulski Archive for Space Telescopes (MAST)111https://mast.stsci.edu/.

We obtained the TESS-SPOC HLSP light curves (Caldwell et al., 2020) for
TOI-1518 from MAST. The SPOC data includes two versions of the photometry at
the standard 30-minute cadence: (1) the Simple Aperture Photometry (SAP) light
curve, i.e., the raw photometry extracted from the SPOC pipeline-derived
photometric aperture (Twicken et al., 2010; Morris et al., 2020), and (2) the
Presearch Data Conditioning SAP (PDCSAP) light curve, which has been corrected
for common-mode systematics trends shared by other sources on the detector
(i.e., co-trending basis vectors, or CBVs), while preserving the key
astrophysical signals of interest (Stumpe et al., 2012, 2014; Smith et al.,
2012).

The PDCSAP light curve is considerably cleaner than the SAP photometry, and in
this paper, we present the analysis of the PDCSAP light curve. For
completeness, we carried out an analogous analysis of the SAP light curve.
Systematics were modeled using linear combinations of the CBVs, similar to the
detrending methodology in the SPOC pipeline. We obtained results that are
statistically consistent with the main PDCSAP-derived values to within
$1\sigma$. However, there were residual long-term systematics trends even
after detrending with the CBVs, which led to a roughly $10\%$ increase in
residual scatter from the best-fit light-curve model when compared to the
PDCSAP analysis.

Our analysis methodology closely mirrors the techniques utilized in the
extensive previous work on TESS phase curves (e.g., Shporer et al., 2019; Wong
et al., 2020a, b, c); consult those references for a detailed description of
the data processing and light-curve fitting. The full PDCSAP light curve of
TOI-1518 is shown in Figure 1. Each TESS Sector consists of two spacecraft
orbits, separated by a pause in science observations for data downlink.
Momentum dumps are scheduled during each spacecraft orbit to reset the onboard
reaction wheels. In Sectors 17 and 18, these occurred twice per spacecraft
orbit and are indicated in the plot by vertical blue dashed lines. The
momentum dumps induce small discontinuities in the photometry, as well as
occasional short-term flux ramps. We therefore divide the light curve into
individual segments separated by the momentum dumps and model the remaining
systematics within each segment separately. Significant ramps are trimmed
prior to the final fit; the trimmed points are shown in Figure 1 in red. The
last data segment of Sector 18 is not included in our analysis due to severe
residual systematics. We also apply a 16-point-wide moving median filter to
the light curve after masking the transits and remove $3\sigma$ outliers. The
final light curve contains 1,845 points, divided among 11 segments.

Visual inspection of Figure 1 reveals coherent flux modulations synchronized
to the planet’s orbit, indicative of a phase curve. To examine the harmonic
content of the TESS photometry in more detail, we trim the transits and
secondary eclipses from the light curve (after correcting for instrumental
systematics; see Section 2.2) and generate the Lomb–Scargle periodogram. The
result is plotted in Figure 2. We find a very strong signal at the orbital
frequency, as well as another significant periodicity at the first harmonic of
the orbital period (i.e., two maxima per orbital period).

The phase curve of a star–planet system formally contains contributions from
both the planet and the host star (see review in Shporer 2017). Close-in
exoplanets are tidally-locked, with fixed dayside and nightside hemispheres;
as the planet rotates, the viewing geometry changes, resulting in a periodic
modulation of the observed atmospheric flux that varies as the cosine of the
orbital phase. Massive orbiting companions can also raise a tidal bulge on the
host star’s surface, resulting in a periodic flux modulation that comes to
maximum at quadrature (i.e., a signal with a leading-order term at the first
harmonic of the cosine); this is typically referred to as ellipsoidal
distortion. Lastly, the mutual star–planet gravitational interaction causes
Doppler shifting of the star’s spectrum, producing a modulation in the total
system flux within the bandpass that can sometimes by detected in visible-
light photometry. This so-called Doppler boosting signal has the same phase
alignment as the RV signal, i.e., the sine of the orbital phase.

###  2.2 Full-orbit Phase-curve Model

We fit the full-orbit phase curve with a composite flux model for the planet
$\psi_{p}$ and the star $\psi_{\star}$ (e.g., Wong et al., 2020b, c, 2021):

$\displaystyle\psi_{p}(t)$
$\displaystyle=f_{p}-A_{\mathrm{atm}}\cos(\phi+\delta),$ (1)
$\displaystyle\psi_{\star}(t)$
$\displaystyle=1-A_{\mathrm{ellip}}\cos(2\phi)+A_{\mathrm{Dopp}}\sin(\phi).$
(2)

Here, $A_{\mathrm{atm}}$, $A_{\mathrm{ellip}}$, and $A_{\mathrm{Dopp}}$
indicate the semiamplitudes of the planet’s atmospheric brightness modulation,
the star’s ellipsoidal distortion signal, and the Doppler boosting,
respectively; the signs are assigned so that the measured amplitudes are
positive under normal circumstances. The variables $f_{p}$ and $\delta$
signify the average relative brightness of the planet across its orbit and the
phase shift in the planet’s phase curve, respectively.

We note that the stellar ellipsoidal distortion signal contains additional
higher-order terms (e.g., Morris, 1985; Shporer, 2017). The second-highest
amplitude is expected at the second harmonic of the cosine (i.e.,
$\cos(6\phi)$). However, there is no significant power precisely at that
harmonic in the Lomb–Scargle periodogram (Figure 2); the weak signal around
1.6 d-1 is centered at a slightly higher frequency than the second harmonic
and is likely attributable to low-level residual systematics in the light
curve. Indeed, when fitting for the second-harmonic amplitudes in the light-
curve analysis, we do not measure any significantly nonzero amplitudes.
Therefore, we do not include any higher-order terms of the ellipsoidal
distortion when generating the final set of phase-curve fit results.

Figure 2: Lomb–Scargle periodogram of the detrended TESS PDCSAP light curve of
TOI-1518, with the transits and secondary eclipses removed. Significance
thresholds are indicated by the horizontal lines. The red vertical lines
denote the first three harmonics of the orbital period. There are clear
signals at the orbital frequency and at the first harmonic, corresponding to
the planetary atmospheric brightness modulation and stellar ellipsoidal
distortion, respectively.

The transits and secondary eclipse light curves ($\lambda_{t}$ and
$\lambda_{e}$) are modeled using batman (Kreidberg, 2015). The secondary
eclipse depth (i.e., total dayside hemisphere flux) is related to the phase-
curve parameters via the expression
$D_{d}=f_{p}-A_{\mathrm{atm}}\cos(\pi+\delta)$. Likewise, the hemisphere-
averaged nightside flux is given by
$D_{n}=f_{p}-A_{\mathrm{atm}}\cos(\delta)$. To accurately model the 30-minute
exposures during transit and secondary eclipse, we use an oversampling factor
of 60, i.e., averaging the flux from 30-second subexposures at each timestamp.

Any remaining systematics trends in each light curve segment $k$ are detrended
using generalized polynomials in time:

$S_{N}^{\\{k\\}}(t)=\sum\limits_{j=0}^{N}c_{j}^{\\{k\\}}(t-t_{0})^{j},$ (3)

where $t_{0}$ is the first timestamp of the segment, and $N$ is the order of
the detrending polynomial, which in the final joint fit is set to the order
that minimizes the Bayesian information criterion (BIC) for each segment. The
optimal polynomial orders for the 11 light-curve segments included in our
analysis are 2, 0, 0, 1, 0, 3, 3, 2, 3, 1, and 3. The total astrophysical-
plus-systematics light-curve model, normalized to unity, is

$F(t)=\frac{\psi_{\star}(t)\lambda_{t}(t)+\psi_{p}(t)\lambda_{e}(t)}{1+f_{p}}\times
S_{N}^{\\{k\\}}(t).$ (4)

To obtain an initial set of results from the TESS photometry, we jointly fit
all 11 light-curve segments using the affine-invariant Markov chain Monte
Carlo (MCMC) sampler emcee (Foreman-Mackey et al., 2013). The free
astrophysical parameters in our fit that are unconstrained by any priors
include the transit ephemeris (mid-transit time $T_{c}$ and orbital period
$P$), transit shape parameters (impact parameter $b$ and scaled semimajor axis
$a/R_{\star}$), planet–star radius ratio $R_{p}/R_{\star}$, and the phase-
curve parameters. The predicted Doppler boosting amplitude assuming the RV-
derived mass (see Section 2.8) is roughly 2 ppm — significantly smaller than
the uncertainties on the phase-curve amplitudes. Therefore, we do not fit the
Doppler signal, while allowing $f_{p}$, $A_{\mathrm{atm}}$,
$A_{\mathrm{ellip}}$, and $\delta$ to vary. We also include a uniform per-
point uncertainty parameter $\sigma_{k}$ for each light-curve segment as a
free parameter in order to ensure a reduced $\chi^{2}$ value of one and
retrieve realistic uncertainties on the astrophysical parameters. The median
values of $\sigma_{k}$ range from 147 to 190 ppm across the 11 segments.

The low cadence of the photometry and the grazing nature of the planetary
transit mean that the stellar limb darkening is not well constrained by the
light curve. We employ the standard quadratic limb-darkening law and apply
Gaussian priors to each coefficient. The median values are set to the values
from Claret (2018), interpolated for the measured stellar parameters (see
Section 2.7) of TOI-1518: $u_{1}=0.28$ and $u_{2}=0.23$; the width of the
Gaussian is generously set to 0.05, which is several times larger than the
corresponding range of coefficient values spanned by the stellar parameter
uncertainty regions.

From our preliminary fit to the full TESS light curve, we find that the
transit is grazing, corresponding to a planet–star radius ratio of
$R_{p}/R_{\star}=0.0987\pm 0.0017$ and well-constrained transit-shape
parameters: $b=0.9103\pm 0.0065$ and $a/R_{\star}=4.231\pm 0.064$. We detect
the secondary eclipse with a depth of $\sim$380 ppm and a significant
atmospheric phase-curve modulation with a semiamplitude of roughly 160 ppm.
There is a nearly $5\sigma$ detection of the ellipsoidal distortion signal
from the host star, with a semiamplitude of around 30 ppm.

To probe for deviations from a circular orbit, we also carry out a separate
light-curve fit with the orbital eccentricity $e$ and argument of periastron
$\omega$ as additional free parameters. From the photometry, the orbital
eccentricity is mostly constrained by the timing of the secondary eclipse
relative to the mid-transit time and, to a much lesser extent, the relative
durations of the transit and secondary eclipse. We obtain a tight $2\sigma$
upper limit of $e<0.01$ (formally, $e=0.0031_{-0.0022}^{+0.0047}$); the
inclusion of $e$ and $\omega$ as free parameters is strongly disfavored by the
Bayesian Information Criterion ($\Delta\mathrm{BIC}=16$). The corresponding
$e\cos\omega$ and $e\sin\omega$ values, which relate to offsets in the
secondary eclipse timing and duration, respectively, are
$0.0007_{-0.0012}^{+0.0016}$ and $-0.0005_{-0.0061}^{+0.0030}$. We therefore
conclude that the orbit of TOI-1518b is consistent with circular.

Due to the relatively short timespan contained within each segment, there is a
possibility of small correlations between the coefficients in the detrending
polynomials and the phase-curve parameters. To examine the effect of our
choice of polynomial orders, we experiment with allowing only polynomials up
to first order (i.e., no curvature in the systematics model). The results from
the corresponding joint fit agree well with the aforementioned values. In
particular, the measured secondary eclipse depth, atmospheric brightness
modulation amplitude, and stellar ellipsoidal distortion amplitude are
statistically consistent at much better than the $1\sigma$ level. Therefore,
we conclude that the optimized polynomial orders listed above, which include
orders as high as 3, do not bias the astrophysical parameters in any
significant way.

###  2.3 Ground-based Light Curves

We acquired ground-based time-series follow-up photometry of TOI-1518 as part
of the TESS Follow-up Observing Program
(TFOP)222https://tess.mit.edu/followup. We used the TESS Transit Finder, which
is a customized version of the Tapir software package (Jensen, 2013), to
schedule our transit observations. The photometric data were extracted using
AstroImageJ (Collins et al., 2017).

A full transit was observed from Adams Observatory at the Austin College
(Sherman, TX, USA) 0.6 m telescope on UT 2020 January 5 in I-band
($\lambda_{\rm eff}=806$ nm). A nearly full in-transit portion of a transit
was observed from the Whitin Observatory (Wellesley, MA, USA) 0.7 m telescope
on UT 2020 January 6 in Sloan $g^{\prime}$-band ($\lambda_{\rm eff}=475$ nm).
A full transit was observed from the private Observatory of the Mount (Saint-
Pierre-du-Mont, France) 0.2 m telescope on UT 2020 January 08 in R-band
($\lambda_{\rm eff}=647$ nm). A full transit was observed from the Kotizarovci
Observatory (Viskovo, Croatia) 0.3 m telescope on UT 2020 January 12 in the
Baader R 610 nm longpass band ($R_{\rm long}$; $\lambda_{\rm cut-on}=610$ nm).
A full transit was observed from the Villa ’39 observatory (Landers, CA, USA)
0.36 m telescope on UT 2020 January 24 in B-band ($\lambda_{\rm eff}=442$ nm).
An egress was observed from the University of Saskatchewan Observatory
(Saskatoon, SK, Canada) 0.3 m telescope on UT 2020 March 23 using an Astrodon
Clear with Blue Blocking Filter (CBB; $\lambda_{\rm cut-on}=500$ nm). The
light-curve data are available at ExoFOP-
TESS.333https://exofop.ipac.caltech.edu/tess The raw ground-based transit
light curves are shown in the Appendix.

The follow-up light curves confirm that the TESS-detected event occurs on
target relative to known Gaia stars. We analyze the five transit observations
with full event coverage (i.e., excluding the UT 2020 March 23 egress-only
light curve) by fitting each time series with batman. The mid-transit time,
orbital period, impact parameter, and scaled semimajor axis are constrained by
Gaussian priors based on the results of the TESS phase-curve fit (Section
2.2). Similar to our treatment of the TESS-band transit modeling, the limb-
darkening coefficients are constrained by priors derived by interpolating the
tabulated values in Claret et al. (2013) for the appropriate bandpass to the
measured stellar parameters and uniformly applying a Gaussian width of 0.05.
In the case of the non-standard $R_{\rm long}$ filter used for the UT 2020
January 12 observation, we approximate the bandpass with the Cousins I-band.
The systematics trends in every transit light curve are modeled as a linear
combination of the airmass and the width of the target’s point-spread-
function, along with a constant offset for normalization.

The comparatively low signal-to-noise of the ground-based transit datasets
translates to large relative uncertainties on the measured transit depth,
exceeding 10% across all five visits. Nevertheless, we obtain
$R_{p}/R_{\star}$ values that are consistent with the measurement from fitting
the TESS light curve alone at better than the $2\sigma$ level. Similarly, the
five ground-based transit depths are mutually consistent to within $2\sigma$,
indicating an achromatic transit.

###  2.4 Joint Photometric Analysis

To leverage the additional time baseline and complementary constraints on
transit geometry provided by the follow-up transit light curves, we carry out
a joint analysis of the TESS photometry and ground-based observations. The
orbital ephemeris, transit-shape, and phase-curve parameters are allowed to
freely vary, while the limb-darkening coefficients remain constrained by the
previously-defined priors. The astrophysical light curve and instrumental
systematics are simultaneously modeled for all six datasets in the MCMC
analysis.

QLP/Atlas Parameters | Symbol | Units | Value
---|---|---|---
Right Ascension | RA | — | 23 29 04.224
Declination | Dec | — | +67 02 05.377
V-band Magnitude | $V$ | mag. | $8.952$
Transit and Orbital Parameters |  |  |
Orbital Period | $P$ | days | $1.902603\pm 0.000011$
Mid-transit Time | $T_{c}$ | BJDTDB | $2458787.049255\pm 0.000094$
Radius Ratio | $R_{p}/R_{\star}$ | — | $0.0988^{+0.0015}_{-0.0012}$
Impact Parameter | $b$ | — | $0.9036_{-0.0053}^{+0.0061}$
Scaled Semimajor Axis | $a/R_{\star}$ | — | $4.291_{-0.061}^{+0.057}$
Orbital Eccentricity | $e$ | — |  $<0.01$ ($2\sigma$)
Orbital Inclination† | $i_{p}$ | deg. | $77.84^{+0.23}_{-0.26}$
Phase-curve Parameters |  |  |
Average Relative Planetary Flux | $f_{p}$ | ppm | $204\pm 27$
Planetary Phase-curve Amplitude | $A_{\mathrm{atm}}$ | ppm | $160.4\pm 6.7$
Planetary Phase-curve Offset | $\delta$ | deg. | $-0.7\pm 2.2$
Stellar Ellipsoidal Distortion Amplitude | $A_{\mathrm{ellip}}$ | ppm | $31.3\pm 6.0$
Secondary Eclipse Depth† | $D_{d}$ | ppm | $364\pm 28$
Nightside Flux† | $D_{n}$ | ppm | $43\pm 27$
Dayside Brightness Temperature† | $T_{d}$ | K | $3237\pm 59$
Nightside Brightness Temperature† | $T_{n}$ | K | $1700^{+700}_{-1200}$
Stellar Parameters |  |  |
Effective Temperature | $T_{\rm eff}$ | K | $7300\pm 100$
Metallicity | [Fe/H] | — | $-0.1\pm 0.12$
Surface Gravity | $\log g$ | — | $4.1\pm 0.2$
Projected Rotational Speed | $v\sin i$ | km s-1 | $85.1\pm 6.3$
Stellar Mass | $M_{\star}$ | $M_{\odot}$ | $1.79\pm 0.26$
Stellar Radius | $R_{\star}$ | $R_{\odot}$ | $1.950\pm 0.048$
RV Parameters |  |  |
RV Semiamplitude | $K_{s}$ | m s-1 | $<281\,{(2\sigma)}$
Systemic Velocity | $V_{\rm sys}$ | km s-1 | $-13.94\pm 0.17$
Planetary Parameters |  |  |
Planet Mass | $M_{p}$ | MJ | $<2.3\,{(2\sigma)}$
Planet Radius | $R_{p}$ | $R_{\rm J}$ | $1.875\pm 0.053$
Orbital Semimajor Axis | $a$ | au | $0.0389\pm 0.0011$
Equilibrium Temperature | $T_{\rm eq}$ | K | $2492\pm 38$
Table 1: Parameters for the TOI-1518 (TIC 427761355) planetary system.
Relevant observing information is obtained from the TESS Quick Look Pipeline
(QLP) and Atlas parameters. The V-band magnitude is obtained from the TESS
input catalog (Stassun et al., 2018b). The transit and phase curve parameters
are simultaneously obtained from a joint fit of the full-orbit TESS light
curve and ground-based full-transit photometric datasets (Section 2.4).
Derived parameters (i.e., quantities not directly fit for in the light-curve
analysis) are indicated by the superscript †. The stellar parameters are
determined by fitting a co-added high-resolution spectrum with a stellar model
using Spectroscopy Made Easy and by a model fit to the broadband SED (Section
2.7). The RV parameters are measured from FIES radial velocities (Section
2.8).  Figure 3: Top panel: systematics-corrected and phase-folded TESS light
curve of TOI-1518, binned in 30-minute intervals, with the best-fit phase-
curve model plotted in red. Middle panel: zoomed-in view of the phase-curve
modulations and secondary eclipse. The atmospheric brightness modulation and
ellipsoidal distortion signals are plotted separately in the solid and dashed
blue lines. Bottom panel: corresponding residuals from the best-fit model.
Figure 4: Zoomed-in views of the primary transit (left) and secondary eclipse
(right) of TOI-1518b. The light curves are binned in 3-minute intervals. Note
the difference in vertical scale between the two plots. The difference in out-
of-occultation baselines primarily reflects the planetary phase-curve
modulation. The bottom panels show the corresponding residuals from the best-
fit model.

The results of our joint fit are listed in Table 1. Figure 3 shows the binned,
phase-folded, and systematics-corrected TESS light curve alongside the best-
fit phase-curve model. Close-up views of the primary transit and secondary
eclipse portions of the light curve are provided in Figure 4. The secondary
eclipse and phase-curve modulations are clearly discernible. The detrended
ground-based transit light curves are plotted in the Appendix.

The orbital period of $1.902603\pm 0.000011$ is measured to $\sim$1 s
precision. We obtain a planet–star radius ratio of
$R_{p}/R_{\star}=0.0988^{+0.0015}_{-0.0012}$, which is marginally more precise
than the value derived from the TESS light curve alone. Likewise, we find
slightly-improved values for the impact parameter and scaled semimajor axis:
$b=0.9036^{+0.0061}_{-0.0053}$, $a/R_{\star}=4.291^{+0.057}_{-0.061}$. The
secondary eclipse depth is measured to more than $12\sigma$ significance:
$364\pm 28$ ppm. The atmospheric phase-curve modulation has a semiamplitude of
$160.4\pm 6.7$ ppm. No significant phase shift in the planet’s phase curve is
measured, indicating that the location of maximum brightness on the dayside
hemisphere is well-aligned with the substellar point. The derived nightside
flux is $43\pm 27$ ppm. The ellipsoidal distortion signal from the host star
is detected at $5.2\sigma$ significance, with a semiamplitude of $31.3\pm 6.0$
ppm. All of the phase-curve parameters are statistically identical to the
values that we obtain from fitting the TESS light curve independently. The
planet’s atmospheric brightness modulation and the star’s ellipsoidal
distortion signal are plotted separately in the middle panel of Figure 3.

The full set of marginalized two-parameter posteriors for the fitted
astrophysical quantities (excluding the limb-darkening coefficients) is
plotted in the Appendix. As expected, due to the grazing nature of the
transits and secondary eclipses, there are significant correlations between
$b$, $R_{p}/R_{\star}$, and $f_{p}$, in addition to the typical degeneracy
between $b$ and $a/R_{\star}$.

###  2.5 SPP Speckle Interferometry

TOI-1518 was observed using speckle interferometry on 2020 October 26 with the
SPeckle Polarimeter (SPP; Safonov et al. 2017) on the 2.5 m telescope at the
Sternberg Astronomical Institute of Lomonosov Moscow State University (SAI
MSU). The spectral band has a central wavelength of $880$ nm and a FWHM of 70
nm. The detector has a pixel scale of 20.6 mas px-1, and the angular
resolution was 89 mas. The detection limit for faint stellar companions is
provided in Figure 5. We did not detect any companion brighter than this
limit, e.g., 6.5 mag at $1^{\prime\prime}$.

Figure 5: SPP $5\sigma$ contrast curve for TOI-1518 with autocorrelation
function (ACF) inset. The observations were obtained at
$\lambda_{\mathrm{c}}=880$ nm (FWHM$=$70 nm).

###  2.6 EXPRES Spectroscopy

EXPRES is an ultra-stable optical spectrograph recently commissioned at the
Lowell Discovery Telescope (Levine et al., 2012). It is designed for extreme-
precision radial velocity surveys (see Jurgenson et al., 2016; Blackman et
al., 2020; Petersburg et al., 2020; Brewer et al., 2020, for details about the
instrument specifications and reduction pipeline) and also has the capacity
for atmospheric characterization (see, for example, the recent study of ultra-
hot Jupiter MASCARA-2b by Hoeijmakers et al. 2020). One transit of TOI-1518b
was observed on the night of 2020 August 2, involving 41 $\sim$300 s
exposures. The extracted spectra have a signal-to-noise (S/N) of $\sim$20–40
for pixels in the continuum. Orders were continuum normalized (Petersburg et
al., 2020), and subsequently stitched together to form one-dimensional
spectra. Telluric absorption from O2 and H2O in Earth’s atmosphere was
corrected with molecfit (Smette et al., 2015) in the geocentric rest-frame
using similar fitting parameters as Allart et al. (2017). Indeed, telluric
modeling with molecfit has become a frequent step in high-resolution optical
atmosphere studies (e.g. Casasayas-Barris et al., 2019), and is advantageous
over empirical models for resolving some atmospheric spectral features
(Langeveld et al., 2021).

###  2.7 Spectroscopic Modeling

Before analyzing the transit, we used Spectroscopy Made Easy (SME 423; Valenti
& Piskunov 1996) to infer stellar parameters from the high-resolution spectra.
The analysis closely follows that of Brewer et al. (2016), including the
choice of fitting parameters and wavelength segments. The model made use of a
VALD3 line-list (Ryabchikova et al., 2015), an ATLAS9 atmospheric model
(Kurucz, 1993; Heiter et al., 2002), and a Gaussian convolution instrument
profile with $R=137,000$. Microturbulence was fixed at 0.85 km s-1, and
macroturbulence was scaled to $T_{\rm eff}$ following the parametrization of
Brewer et al. (2016). However, the fit was largely insensitive to these
parameters since the broadening is completely dominated by stellar rotation.
The rotational broadening also prevents a robust fit to abundances of
individual species. We opt to solve for a global [M/H] with the assumption of
a solar abundance pattern for individual elements.

The true uncertainties on effective temperature ($T_{\rm eff}$), metallicity
([Fe/H]), and rotation speed ($v\sin i$) are difficult to gauge (Piskunov &
Valenti, 2017). The Levenberg-Marquardt optimization algorithm involves
computing a curvature matrix at the minimum of the objective function, the
inverse of which is the covariance matrix. The square root of the diagonal
elements are the formal uncertainties on the parameters, assuming that the
dominant source of uncertainty is from measurement errors (i.e. Poisson
statistics on the spectrum). The actual uncertainty is dominated by systematic
effects and model errors, as opposed to measurement errors. Piskunov & Valenti
(2017) describe a method to incorporate model errors. It involves measuring
the sensitivity of each spectral pixel to changes in the parameters and
estimating the change necessary to reduce the fit residuals to zero. The
cumulative distribution function (CDF) of these parameter perturbations is
then calculated. The central region of each CDF gives an estimate of the model
error. Piskunov & Valenti (2017) discuss this method in greater detail, and we
adopt it for our analysis.

We find that TOI-1518 is a rapidly rotating F0 star with $v\sin i=85\pm 6$ km
s-1, which agrees with expectations for this spectral type (Nielsen et al.,
2013). A fitted [Fe/H] of $-0.1\pm 0.12$ is low for a star hosting a hot
Jupiter (Fischer & Valenti, 2005); only $\sim 4\%$ of planet hosts have [Fe/H]
near $-0.1$. However, the uncertainties on [Fe/H] are large due to the
widening and blurring of spectral lines (a consequence of the rapid rotation),
so the star may be more metal-rich than the best-fit value suggests. The best
fit effective temperature and surface gravity are $T_{\rm eff}=6910\pm 445$ K
and $\log g=3.97\pm 0.62$, respectively. More detailed investigation of the
stellar spectrum might warrant modeling non-LTE effects in the deepest lines
and calibrating line positions and $\log gf$ values. However, these
considerations are most important for cooler stars with total rotational
broadening $\lesssim 10$ km s-1 (Brewer et al., 2016), and their impact on
TOI-1518 is reduced due to the rotation speed. Measurements of $v\sin i$ and
[Fe/H] are listed in Table 1. However, we opt to report the better constrained
measurements of $\log g$ and $T_{\rm eff}$ from our spectral energy
distribution modeling (see below). Our inferred $v\sin i$ is used to analyze
the Rossiter-McLaughlin (RM) effect (Rossiter, 1924; McLaughlin, 1924) in
Section 3.2.

As an independent determination of the stellar parameters, we performed an
analysis of the broadband spectral energy distribution (SED) of the star
together with the Gaia DR2 parallaxes (adjusted by $+0.08$ mas to account for
the systematic offset reported by Stassun & Torres, 2018), following the
procedures described in Stassun & Torres (2016); Stassun et al. (2017, 2018a).
We took the $B_{T}V_{T}$ magnitudes from Tycho-2, the $BVi$ magnitudes from
APASS, the $JHK_{S}$ magnitudes from 2MASS, the W1–W4 magnitudes from WISE,
the $GG_{\rm BP}G_{\rm RP}$ magnitudes from Gaia, and the NUV magnitude from
GALEX. Together, the available photometry spans the full stellar SED over the
wavelength range 0.2–22 $\mu$m (see Figure 6).

Figure 6: Spectral energy distribution of TOI-1518. Red symbols represent the
observed photometric measurements, where the horizontal bars represent the
effective width of the passband. Blue symbols are the model fluxes from the
best-fit Kurucz atmosphere model (black).

We performed a fit using Kurucz stellar atmosphere models, with the free
parameters being the effective temperature ($T_{\rm eff}$), metallicity
([Fe/H]), surface gravity ($\log g$), and extinction ($A_{V}$); the extinction
was restricted to the maximum line-of-sight value from the dust maps of
Schlegel et al. (1998). The resulting fit (Figure 6) has a $\chi^{2}$ of 20.3
(with 12 degrees of freedom) and best-fit parameters $T_{\rm eff}=7300\pm 100$
K, $\mathrm{[Fe/H]}=0.0\pm 0.2$, $\log g=4.1\pm 0.2$, and $A_{V}=0.05\pm
0.05$. The relatively low $A_{V}$ may be surprising considering the low
galactic latitude; however, this $A_{V}$ is consistent with the 3D dust maps
for this system’s position from Green et al. (2019).

Integrating the (unreddened) model SED gives the bolometric flux at Earth,
$F_{\rm bol}=(6.52\pm 0.31)\times 10^{-9}$ erg s-1 cm-2. Taking $F_{\rm bol}$
and $T_{\rm eff}$ together with the Gaia DR2 parallax ($4.398\pm 0.033$ mas)
gives a stellar radius of $R_{\star}=1.950\pm 0.048$ R⊙. In addition, we use
$R_{\star}$ together with $\log g$ to obtain an empirical mass estimate of
$M_{\star}=1.79\pm 0.26$ M⊙, which is consistent with that calculated via the
empirical relations of Torres et al. (2010) — $M_{\star}=1.70\pm 0.12$ M⊙.

###  2.8 FIES Spectroscopy

Starting on June 14th 2020 and ending on February 3rd 2021, we monitored
TOI-1518 with the Nordic Optical Telescope (NOT; Djupvik & Andersen, 2010)
using the FIber-fed Echelle Spectrograph (FIES; Telting et al. 2014). This was
done in order to constrain the out-of-transit Doppler motion of the star,
although the high rotation rate of the star broadens the spectral lines and
makes it difficult to measure. The FIES high-resolution fiber reaches $R\sim
67,000$ and covers wavelengths from 3760 Å to 8840 Å with no gaps below 8200
Å. We obtained 22 spectra, which we extract as described in Buchhave et al.
(2010) and assign wavelengths using ThAr calibrations taken immediately before
and after each exposure. The SNR per resolution element ranges from 49 to 141,
measured in the 5500 Å spectral order. We did not include RVs from the EXPRES
spectra when constraining the Doppler motion, as this would require an extra
instrumental offset parameter for a single night of data.

To extract the radial velocities from the FIES spectra, we perform a least-
squares deconvolution (LSD) analysis to derive the spectroscopic broadening
profiles from each observation (Donati et al., 1997). We deconvolve each
spectrum against a synthetic non rotating spectral template generated via the
ATLAS9 library (Castelli & Kurucz, 2003), and fit the resulting line profiles
with a kernel incorporating the rotational, instrumental, and macroturbulent
components of the line broadening function, similar to the recent analysis of
HAT-P-70 by Zhou et al. (2019). The extracted RVs are listed in Table 2.

One point is excluded from the analysis, since it overlaps with the transit.
Using the radvel package (Fulton et al., 2018), we model the orbit as circular
with no other planets in the system; the stipulation of a circular orbit is in
line with the results of our TESS light-curve fit, which indicated a $2\sigma$
upper limit on orbital eccentricity of 0.01 (Table 1). We define Gaussian
priors for period and time of conjunction (using the values and uncertainties
from Table 1), as well as a broad, uniform prior on the RV semiamplitude
$K_{s}$. We sample the parameter space with an MCMC analysis using the default
radvel setup and let the software run until it determines that the chains are
well-mixed.

The $K_{s}$ posterior distribution peaks near its median at 152 m s-1 with a
$1\sigma$ error of 75 m s-1, i.e. less than $2\sigma$ significance. We derive
a 95% upper limit of 281 m s-1. Adopting the stellar mass determined in
Section 2.7 and the orbital inclination determined in Section 3.2, this
corresponds to a planetary upper mass limit of 2.3 MJ, well within
expectations for hot Jupiters.

To determine the systemic velocity, we compute the weighted mean of the
measured RVs, $-14.79\pm 0.06\,\mbox{km s${}^{-1}$}$, which must be corrected
for an instrumental offset of $-0.87\pm 0.16\,\mbox{km s${}^{-1}$}$, found
from standard stars. We arrive at a systemic velocity $V_{\textrm{sys}}$ of
$-13.94\pm 0.17\,\mbox{km s${}^{-1}$}$. The derived RVs are displayed in
Figure 7, with the posterior distribution of $K_{s}$ visualized along with the
phase-folded velocities. The observations provide generally good sampling of
the orbital phase, and have mean cadence of $11.2$ days between adjacent
observations; we do not expect the RV signature to arise from sampling
artifacts or aliases. More data is needed though to determine if the scatter
in the RVs could be caused by one or more additional planets in the system.

Time (BJD) | Phase |  $v$ (km s-1) |  $\sigma_{v}$ (km s-1)
---|---|---|---
2459014.69572 | 0.65 | -14.64 | 0.31
2459021.70842 | 4.34 | -15.19 | 0.29
2459036.66525 | 12.20 | -15.28 | 0.36
2459037.65405 | 12.72 | -15.00 | 0.30
2459038.66970 | 13.25 | -15.34 | 0.30
2459039.72700 | 13.81 | -14.77 | 0.25
2459093.63742 | 42.14 | -14.29 | 0.49
2459095.66605 | 43.21 | -14.87 | 0.24
2459105.52107 | 48.39 | -14.64 | 0.24
2459119.56182 | 55.77 | -14.70 | 0.26
2459123.52136 | 57.85 | -14.43 | 0.27
2459132.56838 | 62.60 | -14.63 | 0.25
2459133.54446 | 63.12 | -14.67 | 0.29
2459134.53120 | 63.63 | -14.64 | 0.27
2459167.49441 | 80.96 | -15.03 | 0.30
2459169.40210 | 81.96 | -14.92 | 0.29
2459182.36116 | 88.77 | -14.71 | 0.38
2459202.62976 | 99.43 | -14.54 | 0.43
2459236.39450 | 117.17 | -14.95 | 0.22
2459247.35102 | 122.93 | -14.23 | 0.28
2459248.35228 | 123.46 | -14.73 | 0.24
2459249.34239 | 123.98 | -15.01 | 0.28
Table 2: Radial velocities of TOI-1518 extracted from FIES spectra. Columns
correspond to the timestamp of the exposure, orbital phase, velocity, and
uncertainty on velocity. Figure 7: Out-of-transit RVs measured with FIES.
Upper panel shows the full RV time series. Lower panel shows the same RVs
phase-folded from $T_{c}$ with the known orbital period. The $K_{s}$ posterior
distribution is visualized as shaded, purple curves in the background (darker:
higher density). The last observation (gray) overlaps with the transit and has
therefore been excluded from the fit. While the data have large uncertainties,
the amplitude of the velocity variation is consistent with a planetary
companion of $M_{p}<2.3~{}M_{\rm J}$.

##  3 The Spectroscopic Transit

In this section, we describe the methods used to analyze the spectroscopic
transit observation from EXPRES. Cross-correlation was performed with the
X-COR pipeline, previously used for atmospheric detections in WASP-121b (Cabot
et al., 2020; Ben-Yami et al., 2020) and MASCARA-2b (Hoeijmakers et al.,
2020). Cross-correlation has become a standard approach for exoplanet
atmospheric analyses at high-resolution (e.g., Snellen et al., 2010; Brogi et
al., 2012; Birkby et al., 2013). This method relies on resolving the orbital
motion of the planet via its Doppler shift on absorption lines (or more
recently emission lines, as shown by Nugroho et al. 2017 and Pino et al.
2020). While individual lines are generally low-S/N, their contributions may
be stacked by cross-correlating an atmospheric model with the data. Then, one
can analyze the resultant cross-correlation function (CCF). This technique has
led to a slew of molecular detections in the near-infrared (NIR), as well as
atomic and ion detections in the optical, starting with KELT-9b (Hoeijmakers
et al., 2018). Please see Madhusudhan (2019) and Ben-Yami et al. (2020) for
more examples of recent atmospheric detections at high-resolution. We briefly
discuss the relevant methods in the following subsection. We then turn our
attention to the RM effect and atmospheric signals present in the CCFs.

###  3.1 Detrending and Cross-Correlation

The most prominent features in the time-series spectra of TOI-1518b are
absorption lines originating in the stellar photosphere, as well as telluric
lines caused by Earth’s atmosphere. As mentioned above, we corrected tellurics
by fitting and dividing each spectrum by a molecfit model. The spectra were
then linearly interpolated onto a common 0.01 Å wavelength grid in the
barycentric rest-frame. We observed a significant narrow sodium absorption
component in the original spectra, which is likely due to the interstellar
medium. Next, we co-added all out-of-transit spectra into a master $F_{\rm
out}$ and then divided each individual spectrum by $F_{\rm out}$. Interstellar
medium features were removed through division by $F_{\rm out}$ since we opted
to not correct for the RV motion of the star (Casasayas-Barris et al., 2018).
Since stellar lines are significantly broadened from rotation, the RV motion
has negligible effect on the planet’s transmission spectrum. Remaining
broadband variations in the spectra were removed by a high-pass Gaussian
filter with a standard deviation of 75 pixels. We restricted our analysis to
the region $4000-6800$ Å. The S/N falls off at bluer wavelengths, and redder
wavelengths suffer from particularly severe telluric absorption. Throughout
the analysis, about $1\%$ of the data were masked to avoid particularly low
S/N pixels on the blue edge of the spectrum and within Balmer lines.

Cross-correlation was performed between each transmission spectrum and a
continuum-subtracted PHOENIX stellar model (Husser et al., 2013). The model
parameters were selected from a grid and chosen to be close to the inferred
parameters: $T_{\rm eff}=7000$ K, $\log g=4.0$ and [Fe/H] $=0.0$. The CCF is
essentially a sliding dot product between the observed spectra and the model
template. It is defined as a function of time $t$ and velocity $v$:

${\rm CCF}(v,t)=\frac{\sum_{i}{f(i|t)m(i|v)w(i)}}{\sum_{i}m(i|v)w(i)}.$ (5)

Here, the observed spectrum $f(i|t)$ corresponds to the flux in pixel $i$ at
time $t$. The PHOENIX stellar template, denoted by $m(i|v)$, has been Doppler
shifted by some velocity $v$ and is interpolated onto the observed wavelength
grid. The weighting term $w(i)$ is chosen to be the inverse time variance of
each pixel, so as to downweight contributions from pixels previously in the
cores of stellar or telluric lines. The CCF velocities are a grid spanning
$-500$ to $+500$ km s-1 in increments of 2 km s-1.

###  3.2 Spin-Orbit Misalignment

Figure 8: Cross-correlation function between the PHOENIX stellar template and
individual transmission spectra. Top Panel: CCF annotated with the start of
ingress and egress. The Doppler shadow (dark) and atmospheric trail (light)
form a “V” shape with a vertex at about $-30$ km s-1. Middle Panel: Doppler
shadow model as described in the text. Bottom Panel: Corrected CCF where the
Doppler shadow model has been subtracted.

Although we have isolated the planetary atmospheric transmission spectrum,
there are residuals at former locations of stellar lines that arise from the
division by $F_{\rm out}$. While $F_{\rm out}$ is a good template for the out-
of-transit stellar spectrum, the stellar line profiles during transit are
distorted because the planet occults part of stellar disk. The projected
location of the planet against the stellar disk changes throughout the
transit, dependent on its impact parameter $b$ and projected obliquity
$\lambda$. The star has a projected rotation speed $v\sin i$, and the flux
emitted at each point on the star’s surface is Doppler shifted by some local
velocity. The transit removes part of the integrated stellar flux, and breaks
the symmetry between each side of the rotating star. This phenomenon is known
as the Rossiter-McLaughlin effect. It is observed by the apparent “Doppler
shadow” in the CCFs (Collier Cameron et al., 2010a), where a dark trail traces
the local velocity of the occulted stellar region.

We model the shadow in a similar fashion as Hoeijmakers et al. (2020) and show
the steps in Figure 8. First, we fit a double-Gaussian profile (sum of two
Gaussians) to the Doppler shadow in each CCF row and record the inner
profile’s fitted mean, standard deviation, and amplitude. The inner profile
models the core of the Doppler shadow, whereas the outer profile models
positive wings on either side that result from normalizing the spectra. The
second Gaussian’s mean was fixed to that of the first, and the standard
deviation was fixed to 18 km s-1. A third degree polynomial is then fit to the
means as a function of time, and then evaluated at the times of each exposure.
This step was repeated for the remaining fitted parameters. Finally, the
Doppler shadow was modeled as a series of double-Gaussian profiles, with
parameters determined by the above polynomials. The polynomials ensure that
the model smoothly varies in time. While this is not a sophisticated physical
model of the shadow, it is effective at correcting the CCF so that the Doppler
shadow does not adversely affect the atmospheric analysis. Serendipitously,
the shadow and planetary signal do not overlap except for a small window at
the start of transit. This configuration is only possible when the planet’s
path is roughly parallel to the projected stellar rotation axis and the
transit takes place near the limb of the star. Nevertheless, it is still
important to model out the Doppler shadow to correctly interpret the S/N of
the atmospheric signal.

The path traced out by the Doppler shadow provides additional constraints on
the transit geometry (Collier Cameron et al., 2010b; Bourrier et al., 2015;
Cegla et al., 2016). The portion of the stellar disk occulted by the planet
has a local velocity

$v_{\star}(t)=x_{\perp}(t)v\sin i.$ (6)

The orthogonal distance $x_{\perp}$ is determined by the position of the
planet:

$x_{\perp}(t)=x_{p}(t)\cos(\lambda)-y_{p}(t)\sin(\lambda)$ (7)
$x_{p}(t)=\frac{a}{R_{\star}}\sin(2\pi\phi)$ (8)
$y_{p}(t)=-\frac{a}{R_{\star}}\cos(2\pi\phi)\cos(i_{p}).$ (9)

Therefore, we can obtain independent constraints on ${a}/{R_{\star}}$,
$\lambda$, $v\sin i$, and $i_{p}$ from the light curve and spectrum fitting
(note the distinction between $i_{p}$ and stellar inclination $i$, the latter
of which we do not investigate here; however it also may be probed by
considering differential rotation (Cegla et al., 2016)). We run an MCMC
routine that samples these parameters and fits the path of the shadow
described by the polynomial fit described above. As an initial check, we use
uniform priors: $2<{a}/{R_{\star}}<12$, $0<\lambda<2\pi$ and $0<i_{p}<\pi$. We
define Gaussian priors for the rotation speed and global offset: $v\sin
i\sim\mathcal{N}(\mu=80,\>\sigma=50)$ km s-1, $V_{\rm
sys}\sim\mathcal{N}(\mu=-14.5,\>\sigma=2)$ km s-1. The results are not
strongly dependent on the choice of prior for the global offset, owing mainly
to the large rotation speed. The sampler includes 15 walkers with 50,000 steps
each. We set the uncertainty on each point equal to the standard deviation of
the Gaussian profile. We assume that the difference between each data point
and the model is independent and normally distributed. We discard the first
5,000 steps and thin the chains by a factor of 40 (approximately the
autocorrelation time).

From this initial analysis, we obtain a scaled semimajor axis
$a/R_{\star}={2.95^{+0.95}_{-0.72}}$. The inclination is in better agreement
with Table 1, at $i_{p}={76.1^{+3.3}_{-4.9}}$ degrees. We also note a strong
correlation between $\lambda$ and $i_{p}$. Next, we rerun the MCMC using
photometrically-derived priors on $a/R_{\star}$ and $i_{p}$ in order to
establish a tighter constraint on obliquity. The final results of our MCMC
analysis, listed in Table 3, show that TOI-1518b is a highly-misaligned,
retrograde planet, with $\lambda={240.34^{+0.93}_{-0.98}}$ degrees. Indeed,
close-in gas giants around hot stars are commonly misaligned (Winn et al.,
2010). Companions with mass $\gtrsim 3$ MJ around hot stars are less likely to
be found in retrograde orbits (Hébrard et al., 2011; Triaud, 2018), but the
RV-derived mass of TOI-1518b is below this threshold.

RM Parameter | Symbol | Units | Value
---|---|---|---
Scaled Semimajor Axis | $a/R_{\star}$ | - | ${{4.272^{+0.058}_{-0.057}}}$
Proj. Obliquity | $\lambda$ | deg. | ${240.34^{+0.93}_{-0.98}}$
Orbital Inclination | $i_{p}$ | deg. | ${77.92\pm 0.24}$
Proj. Rot. Speed | $v\sin i$ | km s-1 | ${74.4\pm 2.3}$
Table 3: Rossiter McLaughlin (RM) parameters, inferred by fitting the path
traced by the Doppler shadow in Section 3.2. We used the physical model of
Cegla et al. (2016) and the emcee sampler (Foreman-Mackey et al., 2013). Free
parameters included the above four as well as $V_{\rm sys}$, which returned a
posterior distribution that was very similar to its prior Gaussian
distribution. The parameters $a/R_{\star}$ and $i_{p}$ were constrained by
Gaussian priors derived from the results of our TESS light-curve fit (Table
1).

###  3.3 $K_{p}-V_{\rm sys}$ Analysis

Closer inspection of Figure 8 shows a faint, white trail spanning
approximately $\pm 50$ km s-1. This feature is a signature of the planet’s
atmosphere. Throughout the transit, the planet’s apparent radial velocity
changes as it moves towards and then away from the observer, given by

$v_{p}(t)=-K_{p}\sin(2\pi(t-T_{c})/P),$ (10)

where $K_{p}$ is the semiamplitude of the planet’s radial velocity. Because
the planet orbits close in, the change in velocity is of order tens of km s-1.
The CCF at each time $t$ peaks when the PHOENIX model template is Doppler
shifted by the planet’s velocity, and features in the model line up with
features in the actual transmission spectrum. The result is a trail in the
CCFs that traces out a small portion of a sinusoidal curve. The planetary
signal may be further enhanced by aligning and co-adding CCF rows, thus
stacking the peaks and improving the signal’s S/N. The slope of the CCF trail
near transit is completely determined by $K_{p}$ through Equation 10. It is
also offset from 0 by the systemic velocity $V_{\rm sys}$. It is useful to
determine $K_{p}$ and $V_{\rm sys}$ by sampling values from a grid and
attempting to shift and stack the CCFs for each combination of values (Brogi
et al., 2012). The signal is maximized at the correct set of values.

The CCF trail only appears if the cross-correlation template contains features
present in the planet’s transmission spectrum. The trail in Figure 8 indicates
that the atmosphere contains neutral and/or ionized species present in the
PHOENIX spectrum. The absorption line positions and relative strengths are
unique to each species. Therefore, we can cross-correlate with a model
template containing only one species, and then perform the $K_{p}-V_{\rm sys}$
analysis to search for an atmospheric signal. If the stacked CCF contains a
sufficiently high significance peak, then we confirm the presence of that
species in the atmosphere of the planet. Here, we define detection
significance (S/N) as the number of standard deviations that the CCF peak lies
away from mean of all values, for all combinations of $K_{p}$ and $V_{\rm
sys}$. Many species of interest are present in the stellar spectrum and have a
Doppler shadow in their CCFs. Therefore, after cross-correlating with each
model template, we scale the shadow model obtained in Section 3.2 by a best-
fitting constant value and subtract it from the the CCF.

###  3.4 Transmission Spectrum Model

During a planet transit, a fraction of the stellar light is filtered by the
planetary atmosphere. To compute the high-resolution transmission spectra of
the planet’s atmosphere, we first need to calculate the opacities of the
elements in the atmosphere. In this work, the Fe and Fe+ opacities were
computed using the HELIOS-K software (Grimm et al. 2021). Our models for Fe
and Fe+ make use of the line-list tables from Kurucz (2018). The lines for
both Fe and Fe+ were computed assuming Voigt profiles, 0.032 cm-1 spectral
resolution, and a fixed line cutoff of 100 cm-1. To calculate the transmission
spectra, we developed our code based on the simple formalism presented in
Gaidos et al. (2017) and Bower et al. (2019). Our model computes the effective
tangent height in an atmosphere that was discretised in 200 annuli. The model
included some simplifications due to the unknown composition of the atmosphere
of TOI-1518b and a weakly constrained planet bulk density: we assumed a
surface gravity of $\log g=3$ and an atmosphere in chemical equilibrium. The
chemical calculations were done with the open-source code FastChem (Stock et
al. 2018), assuming solar metallicities. We include in our model the H-
bound–free and free–free absorption from John (1988). As shown in Kitzmann et
al. (2018), the H- continuum in UHJs is generally between 1 mbar and 10 mbar.
Each high-resolution transmission spectrum includes Fe or Fe+ along with H-
continuum absorption and scattering by H and H2. We generated a grid of high-
resolution transmission spectra assuming isothermal atmospheres ranging from
2000 to 4000 K in steps of 500 K. Following subtraction of the continuum with
a sliding maximum filter and convolution with a Gaussian filter to match the
EXPRES instrumental resolution, these models serve as cross-correlation
templates.

##  4 Atmospheric Characterization

###  4.1 Detections

Figure 9: Atmospheric detections in TOI-1518b and their $K_{p}-V_{\rm sys}$
maps. The top-left corner of each panel indicates the cross-correlation
template, and the bottom right corner lists properties of the peak value,
including S/N of the detection, the average absorption depth after co-adding
CCF rows ($\langle\delta\rangle$), and the maximal value of $K_{p}$, and
$V_{\rm sys}$. In all panels the Doppler shadow has been corrected per the
methods in Section 3.2. For reference, we show results from cross-correlation
with the PHOENIX spectrum used to model the Doppler shadow, revealing there
are species common to both the planet’s atmosphere and star. The subsequent
panels show results from cross-correlating with templates containing Fe and
Fe+. In each panel, the white dotted lines indicate the $K_{p}$ and $V_{\rm
sys}$ with the highest signal.

We detect Fe in the atmosphere of TOI-1518b at the ${5.2}\sigma$ level. We
also report evidence of Fe+ at the ${3.4}\sigma$ level. The PHOENIX model,
which contains both species in addition to other atoms and ions, yields an
enhanced atmospheric detection at $5.9\sigma$ confidence, while a combined
Fe/Fe+ model yields a ${5.4}\sigma$ detection. The Doppler shadow correction
removes an artifact that otherwise biases detection significances. The $K_{p}$
and $V_{\rm sys}$ corresponding to the peak value are consistent across the
various templates. For the PHOENIX model we find $K_{p}={163^{+49}_{-30}}$ km
s-1 and $V_{\rm sys}={-17^{+3}_{-2}}$ km s-1. For Fe the values are
${K_{p}={157^{+68}_{-44}}}$ km s-1 and ${V_{\rm sys}={-16^{+2}_{-4}}}$ km s-1,
and for Fe+ they are ${K_{p}={178^{+41}_{-62}}}$ km s-1 and ${V_{\rm
sys}={-18^{+3}_{-3}}}$ km s-1. Uncertainties correspond to the range of
$K_{p}$ and $V_{\rm sys}$ within a $1\sigma$ contour around the peak. Because
we only sample a small portion of the planet’s orbit, only loose constraints
on the semiamplitude $K_{p}$ are possible. The $V_{\rm sys}$ found here is
offset by about $3$ km s-1 at the $\sim 1-2\sigma$ level. This blueshift may
indicate winds in the upper atmosphere of the planet (Miller-Ricci Kempton &
Rauscher, 2012; Casasayas-Barris et al., 2019). Using values in Table 1, we
predict a planetary RV semiamplitude of $K_{p}=2\pi a\sin i_{p}/P=217.4\pm
6.2$ km s-1. This value is higher than the $K_{p}$ measured from cross-
correlation, but still consistent to within the $1\sigma$ uncertainties.

Equation 5 involves a normalization term in the denominator that allows the
CCF to return a physically meaningful quantity (Hoeijmakers et al., 2019). The
CCF peak is a weighted average of the depths of individual lines in the
transmission spectrum of the planet. In practice, the average depth depends on
the weighting used for low S/N pixels ($w(i)$) and the wavelength range of the
cross-correlation; it also does not correspond to the depth of any particular
line. However, it provides an order-of-magnitude estimate of typical
absorption depths, and hence the altitude of the species in the exoplanet’s
atmosphere. We refer to the average absorption depth as
$\langle\delta\rangle$, which is equal to the peak value of the stacked CCF
over all $K_{p}$ and $V_{\rm sys}$ combinations. As shown by Hoeijmakers et
al. (2019), Fe lines probe much deeper in the atmosphere than Fe+ lines under
chemical equilibrium. While Fe+ lines are stronger in the optical, they are
fewer in number; Fe+ absorption is generally much stronger in the near
ultraviolet (e.g. Sing et al., 2019). We find average absorption depths of
${{(3.6\pm 0.8)\times 10^{-4}}}$ and ${{(1.5\pm 0.4)\times 10^{-3}}}$ for Fe
and Fe+ respectively (note, the significance of the Fe+ signal only indicates
evidence of the species, but we can still proceed with using the signal to
learn about the planet).

Per Equation 5, the average absorption depth depends on the absolute depths of
lines in the data, as well as the relative (but not absolute) depths of lines
in the model. The results above are of the same order of magnitude as those
for KELT-9b (Hoeijmakers et al., 2019). The height of the atmosphere ($H$)
extends 5–10 scale heights ($H_{\rm sc}$, of length hundreds of kilometers for
hot Jupiters) (Madhusudhan et al., 2014). The excess absorption beyond the
transmission spectrum continuum $(R_{p}/R_{\star})^{2}$ is approximately
$\delta\approx 2R_{p}H/R_{\star}^{2}$; in other words, $H\approx\delta
R_{p}/2(R_{p}/R_{\star})^{2}$. For order of magnitude estimates, we use values
in Table 1 and assume the base of the atmosphere has a pressure of 0.01 bar
(Kitzmann et al., 2018), which is typical for the H- continuum of an UHJ. We
also take $H_{\rm sc}\sim 880$ km, estimated from the measured $T_{\rm eq}$
and $\log g$, as well as taking the mean molecular weight as $\mu=2.3$ for an
H2-dominated atmosphere; however $\mu$ may be affected by H2 dissociation on
the planet’s dayside. While the mass is highly uncertain, we take the
posterior median value of $1.4$ MJ in order to estimate $\log g$. The
resultant pressures corresponding to the absorption are ${P\sim{6\times
10^{-4}}}$ bar for Fe and ${P\sim{2\times 10^{-7}}}$ bar for Fe+.
Interestingly, the blueshift is similar between both Fe and Fe+ signals,
suggesting that high-velocity winds might be fairly consistent across various
depths in the atmosphere.

The 4000 K Fe model returns the highest-significance detection. The Fe
detection significances are 4.2$\sigma$, 4.7$\sigma$, and 5.2$\sigma$ for
temperatures of 2000, 3000, and 4000 K, respectively. The cross-correlation
signal also decreases significantly during the second half of transit. The Fe
detection significance is 4.6$\sigma$ when using exposures from only the first
half of the transit. It drops to 1–2$\sigma$ if only exposures from the second
half are used. This variability could trace differential chemistry between the
morning and evening terminators. For example, Ehrenreich et al. (2020) infer a
lack of neutral Fe vapor on the dayside terminator of WASP-76b based on the
changing Doppler shift of the cross-correlation peak in each of their
exposures. Hoeijmakers et al. (2020) observe slightly stronger Fe absorption
in the second half of a transit of MASCARA-2b, which they suggest could be due
to different temperatures or chemistry between terminators. In the case of
TOI-1518b, additional transits would help improve our confidence that the
observed variability is indeed of physical origin.

###  4.2 Temperature and Circulation

From the stellar radius, we can use the values of $R_{p}/R_{\star}$ and
$a/R_{\star}$ from our photometric analysis to straightforwardly compute the
planet’s radius and orbital semimajor axis: $R_{p}=1.875\pm 0.053$ $R_{\rm J}$
and $a=0.0389\pm 0.0011$ au. We also utilize the stellar parameters from the
SED fit to further characterize the planet’s atmosphere. The relative flux of
the planet $D$ in the TESS bandpass, assuming no reflected starlight (i.e.,
zero geometric albedo), is related to the hemisphere-averaged brightness
temperature $T_{p}$ via the following relation (e.g., Shporer, 2017):

$D=\left(\frac{R_{p}}{R_{\star}}\right)^{2}\frac{\int
F_{\lambda}(T_{p})\tau(\lambda)\lambda d\lambda}{\int
F_{\lambda}(T_{\mathrm{eff}})\tau(\lambda)\lambda d\lambda}.$ (11)

Here, the stellar and planetary flux spectra are given by
$F_{\lambda}(T_{\mathrm{eff}})$ and $F_{\lambda}(T_{p})$, respectively, and
$\tau(\lambda)$ is the transmission function of the TESS bandpass. For
simplicity, we assume that the planet’s emission spectrum is well-modeled by a
blackbody function.

For the stellar spectrum, following the technique described in Wong et al.
(2020c), we use PHOENIX stellar models (Husser et al., 2013) and calculate the
integrated stellar flux in the denominator of Equation (11) for a grid of
stellar parameters in the vicinity of the values derived from the SED fit. We
then construct an empirical polynomial function in
$\left\\{T_{\mathrm{eff}},\mathrm{[Fe/H]},\log g\right\\}$ that smoothly
interpolates these values. The planet’s brightness temperature can then be fit
for using an MCMC routine, with Gaussian priors for $T_{\mathrm{eff}}$,
$\mathrm{[Fe/H]}$, $\log g$, and $R_{p}/R_{\star}$ derived from the SED and
TESS light-curve fits.

We use the secondary eclipse depth and nightside flux (Table 1) to calculate
the corresponding dayside and nightside brightness temperatures of TOI-1518b:
$T_{d}=3237\pm 59$ K and $T_{n}=1700^{+700}_{-1200}$ K. The extremely high
dayside temperature makes TOI-1518b among the hottest exoplanets hitherto
discovered, comparable to other UHJs such as WASP-18b ($3100\pm 49$ K; Wong et
al. 2020b) and WASP-33b ($3105\pm 95$ K; von Essen et al. 2020).

We note that any reflected light off the dayside atmosphere (i.e., nonzero
geometric albedo) would decrease the contribution of the planet’s thermal
emission to the measured secondary eclipse, resulting in a lower inferred
dayside brightness temperature. However, at these high temperatures, all known
condensate species are expected to be in the vapor phase across the dayside
hemisphere, making reflective clouds unlikely (e.g., Helling et al., 2019).
This is supported by emission spectrum modeling of other UHJs spanning optical
and thermal infrared wavelengths, which break the degeneracy between short-
wavelength reflectivity and planetary thermal emission and indicate geometric
albedos consistent with zero (e.g., Shporer et al., 2019; Wong et al., 2020b,
2021).

In the broader context of atmospheric circulation, the measured dayside and
nightside brightness temperatures reflect the amount of absorbed insolation
and the efficiency of day–night heat transport. We can use the simple thermal
balance model outlined in Cowan & Agol (2011) to simultaneously constrain the
Bond albedo $A_{\mathrm{B}}$ and the recirculation efficiency $\epsilon$. In
this parametrization, $\epsilon$ ranges from 0 (no recirculation) to 1
(uniform global temperature). To properly propagate the uncertainties on the
stellar and orbital parameters, we use the methodology described in Wong et
al. (2020c). Due to the highly-uncertain nightside brightness temperature, we
retrieve very poor constraints: $A_{B}<0.2$ ($2\sigma$) and $\epsilon=0.5\pm
0.3$. Higher signal-to-noise is required to construct a more precise picture
of the atmospheric heat budget. This may be achieved either by including
additional visible-wavelength photometry of the system from the TESS Extended
Mission or by obtaining full-orbit phase-curve observations at infrared
wavelengths, where the planet–star contrast ratio is significantly higher.

##  5 Discussion and Conclusions

As there have been only a handful of previous detections of iron in UHJs,
TOI-1518b adds an important additional data point in our efforts to understand
the dynamics and thermal structure in highly irradiated atmospheres. We make a
few concluding remarks about the planet below, and then compare it to other
recently characterized UHJs.

Planet |  $T_{\rm eq}$ (K) |  $R_{p}$ ($R_{\rm J}$) |  $\log g$ (cgs) | Fe (Transmission/Emission) | Reference
---|---|---|---|---|---
TOI-1518b | $2492\pm 38$ | $1.875\pm 0.053$ | $<3.229$ | Y/- | this study
KELT-9b | $4050\pm 180$ | $1.783\pm 0.009$ | $3.30^{+0.11}_{-0.15}$ | Y/Y | H18, H19, P20
MASCARA-2b | $2260\pm 50$ | $1.83\pm 0.07$ | $<3.467$ | Y/- | CB19, S20, H20
WASP-121b | $2358\pm 52$ | $1.865\pm 0.044$ | $2.973\pm 0.017$ | Y/- | D16, C19
WASP-76b | $2228\pm 122$ | $1.854\pm 0.077$ | $2.806\pm 0.034$ | Y/- | E20
WASP-189b | $2641\pm 34$ | $1.619\pm 0.021$ | $3.274^{+0.048}_{-0.042}$ | -/Y | A18, C20, L20, Y20
WASP-33b | $2710\pm 50$ | $1.679^{+0.019}_{-0.030}$ | $3.297^{+0.043}_{-0.041}$ | -/Y | Y19, N20
WASP-19b | $2372\pm 60$ | $1.392\pm 0.040$ | $2.616^{+0.065}_{-0.070}$ | -/- | W16, Se21
TOI-1431b | $2181\pm 95$ | $1.546\pm 0.063$ | $4.148^{+0.043}_{-0.041}$ | -/- | S21, A21
Table 4: Summary of recent high-resolution spectroscopy iron detections,
comparing TOI-1518b to known transiting ultra-hot Jupiters. Values and
uncertainties for equilibrium temperature and planet radius are reported in
the references. Surface gravity was calculated from available parameters, if
not reported explicitly. References: H18 (Hoeijmakers et al., 2018), H19
(Hoeijmakers et al., 2019), P20 (Pino et al., 2020) CB19 (Casasayas-Barris et
al., 2019), H20 (Hoeijmakers et al., 2020), S20 (Stangret et al., 2020), D16
(Delrez et al., 2016), C19 (Cabot et al., 2019), E20 (Ehrenreich et al.,
2020), A18 (Anderson et al., 2018), C20 (Cauley et al., 2020), L20 (Lendl et
al., 2020), Y20 (Yan et al., 2020), Y19 (Yan et al., 2019), N20 (Nugroho et
al., 2020), W16 (Wong et al., 2016), Se21 (Sedaghati et al., 2021), S21
(Stangret et al., 2021), A21 (Addison et al., 2021).

###  5.1 TOI-1518b In the Context of Other Iron Detections

Alkali metals (Na and K) have been detected in transmission for numerous hot
Jupiters (e.g. Sing et al., 2016). Over the past two years, Fe has also become
an increasingly common detected species, albeit mostly in UHJs with $T_{\rm
eq}\gtrsim 2000$ K (Parmentier et al., 2018). Fe traces winds in the upper
atmosphere through the systemic velocity offset of the cross-correlation peak
and is also a potential non-oxide contributor to thermal inversions
(Lothringer et al., 2018). In the literature, Fe has been detected in
transmission in the following exoplanets: KELT-9b (Hoeijmakers et al., 2018,
2019), WASP-121b (Cabot et al., 2019), MASCARA-2b (Stangret et al., 2020;
Hoeijmakers et al., 2020), WASP-76b (Ehrenreich et al., 2020), and TOI-1518b
(this study). Fe has been detected in emission in KELT-9b (Pino et al., 2020),
WASP-189b (Lendl et al., 2020), and WASP-33b (Yan et al., 2020). These targets
are listed in Table 4.

Interestingly, Cauley et al. (2020) do not detect Fe in transmission in
WASP-189b, despite it being one of the brightest and hottest systems and the
fact that Fe is detected in emission (however, the observations were made
under poor weather conditions). We note that, although Ca+ was found in
transmission in WASP-33b (Yan et al., 2019), and Fe in emission (Yan et al.,
2020), there has been no claim of Fe in transmission. Fe may be especially
difficult to detect in WASP-33b due to stellar pulsations. We acknowledge a
few additional recent studies, including the non-detection Fe in WASP-19b
(Sedaghati et al., 2021) which is listed in Table 4 (however this target is
considerably fainter than the others, at $V=12.3$), a recent transmission
spectroscopy study of HD149026b (Ishizuka et al., 2021) (however the Fe signal
was only at $2.8\sigma$), and a non-detection in TOI-1431b (which orbits a
relatively bright $V=8.0$ star; this target is listed in Table 4).

While the statistical sample is small, Fe detections seem to favor
particularly inflated UHJs, potentially with a cutoff around $1.7-1.8$ $R_{\rm
J}$. One explanation is that Fe detections require particularly large
atmospheric scale heights in order for the atoms to imprint sufficiently deep
absorption lines on top of the continuum of the transmission spectrum.
However, the surface gravity, which is inversely proportional to scale height,
does not show a discernible relationship to Fe detections. For example, Fe was
detected in transmission in KELT-9b, whose large mass yields a similar $\log
g$ as WASP-189b. The $\log g$ of TOI-1518b is less than $3.229$ at $95\%$
confidence. There are a few bright targets with $R_{p}<1.7$ $R_{\rm J}$ that
are without detailed, cross-correlation atmospheric analyses, and do not have
reported detections of Fe in transmission: MASCARA-1b (Talens et al., 2017),
KELT-7b (Bieryla et al., 2015), and KELT-17b (Zhou et al., 2016). As more gas
giants are detected and characterized, it will be interesting to see if such a
trend between Fe detection and planetary radius continues to hold.

###  5.2 Photometric Mass Measurement and Caveats

In our analysis of the TESS photometry, we obtain a strong detection of the
ellipsoidal distortion component of the phase-curve variability. This signal
is driven by the tidal response of the stellar surface to the mutual
star–planet gravitational interaction, which in turn depends on the mass ratio
between the two components. It follows that the measured amplitude of the
ellipsoidal distortion signal can be used to obtain an independent estimate of
the planet’s mass.

The ellipsoidal distortion of the star is formally modeled as a series of
cosine terms, with the semiamplitude of the leading term (at the first
harmonic of the orbital phase) related to fundamental parameters of the system
via the following expression (e.g., Morris, 1985; Shporer, 2017):

$A_{\mathrm{ellip}}=\alpha_{\mathrm{ellip}}\frac{M_{p}}{M_{\star}}\left(\frac{R_{\star}}{a}\right)^{3}\sin^{2}i_{p}.$
(12)

Here, the pre-factor $\alpha_{\mathrm{ellip}}$ is a function of the linear
limb-darkening and gravity-darkening coefficients $u$ and $g$ for the host
star:

$\alpha_{\mathrm{ellip}}=\frac{3}{20}\frac{(u+15)(g+1)}{3-u}.$ (13)

Similar to our treatment of the quadratic limb-darkening coefficients in the
TESS phase-curve analysis (Section 2.2), we construct Gaussian priors for $u$
and $g$ using values interpolated from the coefficients listed in Claret
(2017): $u=0.41\pm 0.05$ and $g=0.12\pm 0.05$. We then use Equations (12) and
(13) to construct the posterior for $M_{p}$ through Monte Carlo sampling of
the distribution of values for $A_{\mathrm{ellip}}$, $a/R_{\star}$, $i_{p}$,
$M_{\star}$, $u$, and $g$. We obtain a photometric mass estimate of
$M_{p}=4.8^{+1.3}_{-1.1}$ MJ. This value is significantly ($2.3\sigma$) larger
than the RV-derived mass upper limit of 2.3 MJ.

This discrepancy between the phase-curve-derived and RV-derived masses may be
attributable to oversimplifications in the stellar tidal response formalism.
Gomel et al. (2021) found a discrepancy of up to 30% between the amplitudes of
the ellipsoidal distortion derived from the analytic expressions of Morris
(1985) and those derived numerically. More fundamentally, the classical theory
of stellar ellipsoidal distortion from which Equations (12) and (13) are
derived makes several key assumptions: (1) steady-state approximation, which
assumes that the star is in hydrostatic balance and ignores fluid inertia and
the possibility of dynamical tides, (2) equatorial orbit of the companion, and
(3) no effects from stellar rotation. The last two assumptions in particular
are ostensibly invalid in the case of the TOI-1518 system, which contains a
hot Jupiter on a misaligned orbit around a rapidly-rotating star (see Section
3.2). The fast rotation of the star and the resulting rotational bulge,
combined with the spin-orbit misalignment, mean that the tidal bulge raised by
the planet traverses regions of the stellar surface with significantly
different surface gravities. This is expected to directly affect the tidal
response of the star and the corresponding amplitude of the ellipsoidal
distortion signal.

Another possible contributor to an unexpected first harmonic phase-curve
modulation is the variable stellar irradiation experienced by the planet. This
scenario was explored in detail for the case of KELT-9 — a similarly
misaligned system with an ultra-hot Jupiter around a rapidly-rotating star —
where it was found to be the primary source of the unusual phase alignment of
the measured first harmonic photometric modulation (Wong et al., 2020c). In
short, the rapid stellar rotation induces variations in the effective
temperature of the planet-facing hemisphere, which cause the planetary thermal
emission to change in response to the time-varying insolation. The three-
dimensional orientation of TOI-1518’s rotation axis is not known from the
available data, preventing us from being able to directly model the relative
phasing of this additional irradiation signal (as was done for the KELT-9
system). Nevertheless, we do expect some level of photometric variability at
the first harmonic that is due to the planet’s variable dayside temperature,
which may bias the photometric mass estimate.

The previous discussion serves as a cautionary tale about the reliability of
photometric mass measurements derived from the ellipsoidal distortion signal.
The complexities of the stellar tidal response and the possibility of
additional contributions from the planet’s thermal emission mean that many
systems are susceptible to significant discrepancies between the measured and
expected first harmonic amplitudes. Future RV monitoring of this system will
improve the precision of the planet’s mass.

###  5.3 Conclusion

TESS continues to find numerous transiting exoplanet candidates. As these
planets are confirmed, some are bound to become interesting case studies for
atmospheric characterization. In this paper, we reported the confirmation of
an ultra-hot Jupiter on a close-in, highly misaligned orbit around TOI-1518.
The stellar, planetary, and orbital parameters derived from fitting the TESS
light curve, ground-based transit photometry, and spectral energy distribution
are listed in Table 1. The photometry displays a clear secondary eclipse
signal, as well as phase-synchronized modulations in flux attributed to the
day–night brightness contrast of the planet and the tidal distortion of the
host star. In addition, we searched for neutral and ionized Fe in the
companion’s atmosphere through high-resolution transmission spectroscopy. We
detected Fe at high confidence, and also found evidence for Fe+. TOI-1518b is
highly inflated, which makes it amenable to intensive atmospheric
characterization. The equilibrium temperature of TOI-1518b is in the regime
where the planet might exhibit a thermal inversion (Fortney et al., 2008;
Lothringer et al., 2018; Malik et al., 2019; Gandhi & Madhusudhan, 2019).
This, combined with the brightness of the host star, makes TOI-1518b an
attractive target for follow-up emission spectroscopy (Pino et al., 2020;
Nugroho et al., 2020; Yan et al., 2020).

This paper includes data collected by the TESS mission. Funding for the TESS
mission is provided by the NASA Explorer Program. Resources supporting this
work were provided by the NASA High-End Computing (HEC) Program through the
NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the
production of the SPOC data products. We acknowledge the use of TESS High
Level Science Products (HLSP) produced by the Quick-Look Pipeline (QLP) at the
TESS Science Office at MIT, which are publicly available from the Mikulski
Archive for Space Telescopes (MAST). This work used data from the EXtreme
PREcision Spectrograph (EXPRES) that was designed and commissioned at Yale
with financial support by the U.S. National Science Foundation under
MRI-1429365 and ATI1509436 (PI D. Fischer). We gratefully acknowledge support
for telescope time using EXPRES at the LDT from the Heising-Simons Foundation
and an anonymous Yale donor. We acknowledge support from U.S. National Science
Foundation grant 2009528. This research made use of Lightkurve, a Python
package for Kepler and TESS data analysis (Lightkurve Collaboration, 2018).
This paper is partially based on observations made with the Nordic Optical
Telescope, operated by the Nordic Optical Telescope Scientific Association at
the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto
de Astrofisica de Canarias. K.K.M. gratefully acknowledges support from the
New York Community Trust’s Fund for Astrophysical Research. I.W. is supported
by a Heising-Simons 51 Pegasi b postdoctoral fellowship. A.A.B., B.S.S. and
I.A.S. acknowledge the support of Ministry of Science and Higher Education of
the Russian Federation under the grant 075-15-2020-780 (N13.1902.21.0039).
This paper is partially based on observations made at the CMO SAI MSU with the
support by M.V. Lomonosov Moscow State University Program of Development. VA
was supported by a research grant (00028173) from VILLUM FONDEN. Funding for
the Stellar Astrophysics Centre is provided by The Danish National Research
Foundation (Grant agreement no.: DNRF106). This research made use of exoplanet
(Foreman-Mackey et al., 2021) and its dependencies (Agol et al., 2020; Kumar
et al., 2019; Astropy Collaboration et al., 2013, 2018; Luger et al., 2019;
Salvatier et al., 2016; Theano Development Team, 2016). We also acknowledge
very useful input from an anonymous referee, which improved the clarity and
structure of the manuscript.

## References

  * Addison et al. (2021) Addison, B. C., Knudstrup, E., Wong, I., et al. 2021, arXiv e-prints, arXiv:2104.12078
  * Agol et al. (2020) Agol, E., Luger, R., & Foreman-Mackey, D. 2020, AJ, 159, 123
  * Allart et al. (2017) Allart, R., Lovis, C., Pino, L., et al. 2017, A&A, 606, A144
  * Anderson et al. (2018) Anderson, D. R., Temple, L. Y., Nielsen, L. D., et al. 2018, arXiv e-prints, arXiv:1809.04897
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Bakos et al. (2004) Bakos, G., Noyes, R. W., Kovács, G., et al. 2004, PASP, 116, 266
  * Ballard et al. (2011) Ballard, S., Fabrycky, D., Fressin, F., et al. 2011, ApJ, 743, 200
  * Ben-Yami et al. (2020) Ben-Yami, M., Madhusudhan, N., Cabot, S. H. C., et al. 2020, ApJ, 897, L5
  * Bieryla et al. (2015) Bieryla, A., Collins, K., Beatty, T. G., et al. 2015, AJ, 150, 12
  * Birkby et al. (2013) Birkby, J. L., de Kok, R. J., Brogi, M., et al. 2013, MNRAS, 436, L35
  * Blackman et al. (2020) Blackman, R. T., Fischer, D. A., Jurgenson, C. A., et al. 2020, AJ, 159, 238
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977
  * Bourrier et al. (2015) Bourrier, V., Lecavelier des Etangs, A., Hébrard, G., et al. 2015, A&A, 579, A55
  * Bower et al. (2019) Bower, D. J., Kitzmann, D., Wolf, A. S., et al. 2019, A&A, 631, A103
  * Brewer et al. (2016) Brewer, J. M., Fischer, D. A., Valenti, J. A., & Piskunov, N. 2016, ApJS, 225, 32
  * Brewer et al. (2020) Brewer, J. M., Fischer, D. A., Blackman, R. T., et al. 2020, arXiv e-prints, arXiv:2006.02303
  * Brogi et al. (2012) Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2012, Nature, 486, 502
  * Buchhave et al. (2010) Buchhave, L. A., Bakos, G. Á., Hartman, J. D., et al. 2010, ApJ, 720, 1118
  * Cañas et al. (2019) Cañas, C. I., Stefansson, G., Monson, A. J., et al. 2019, ApJ, 877, L29
  * Cabot et al. (2019) Cabot, S. H. C., Madhusudhan, N., Hawker, G. A., & Gandhi, S. 2019, MNRAS, 482, 4422
  * Cabot et al. (2020) Cabot, S. H. C., Madhusudhan, N., Welbanks, L., Piette, A., & Gandhi, S. 2020, MNRAS, 494, 363
  * Caldwell et al. (2020) Caldwell, D. A., Tenenbaum, P., Twicken, J. D., et al. 2020, Research Notes of the American Astronomical Society, 4, 201
  * Casasayas-Barris et al. (2018) Casasayas-Barris, N., Pallé, E., Yan, F., et al. 2018, A&A, 616, A151
  * Casasayas-Barris et al. (2019) —. 2019, A&A, 628, A9
  * Castelli & Kurucz (2003) Castelli, F., & Kurucz, R. L. 2003, in IAU Symposium, Vol. 210, Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, A20
  * Cauley et al. (2020) Cauley, P. W., Shkolnik, E. L., Ilyin, I., et al. 2020, Research Notes of the American Astronomical Society, 4, 53
  * Cegla et al. (2016) Cegla, H. M., Lovis, C., Bourrier, V., et al. 2016, A&A, 588, A127
  * Charbonneau et al. (2002) Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L. 2002, ApJ, 568, 377
  * Claret (2017) Claret, A. 2017, A&A, 600, A30
  * Claret (2018) Claret, A. 2018, A&A, 618, A20
  * Claret et al. (2013) Claret, A., Hauschildt, P. H., & Witte, S. 2013, A&A, 552, A16
  * Collier Cameron et al. (2010a) Collier Cameron, A., Bruce, V. A., Miller, G. R. M., Triaud, A. H. M. J., & Queloz, D. 2010a, MNRAS, 403, 151
  * Collier Cameron et al. (2010b) Collier Cameron, A., Guenther, E., Smalley, B., et al. 2010b, MNRAS, 407, 507
  * Collins et al. (2017) Collins, K. A., Kielkopf, J. F., Stassun, K. G., & Hessman, F. V. 2017, AJ, 153, 77
  * Cowan & Agol (2011) Cowan, N. B., & Agol, E. 2011, ApJ, 729, 54
  * Davis et al. (2019) Davis, A. B., Wang, S., Jones, M., et al. 2019, arXiv e-prints, arXiv:1912.10186
  * Dawson & Johnson (2018) Dawson, R. I., & Johnson, J. A. 2018, ArXiv e-prints, arXiv:1801.06117
  * Delrez et al. (2016) Delrez, L., Santerne, A., Almenara, J. M., et al. 2016, MNRAS, 458, 4025
  * Djupvik & Andersen (2010) Djupvik, A. A., & Andersen, J. 2010, Astrophysics and Space Science Proceedings, 14, 211
  * Donati et al. (1997) Donati, J. F., Semel, M., Carter, B. D., Rees, D. E., & Collier Cameron, A. 1997, MNRAS, 291, 658
  * Ehrenreich et al. (2020) Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, arXiv e-prints, arXiv:2003.05528
  * Evans et al. (2017) Evans, T. M., Sing, D. K., Kataria, T., et al. 2017, Nature, 548, 58
  * Fabrycky et al. (2014) Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146
  * Fischer & Valenti (2005) Fischer, D. A., & Valenti, J. 2005, ApJ, 622, 1102
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Foreman-Mackey et al. (2021) Foreman-Mackey, D., Savel, A., Luger, R., et al. 2021, exoplanet-dev/exoplanet v0.4.4, , , doi:10.5281/zenodo.1998447.  https://doi.org/10.5281/zenodo.1998447
  * Fortney et al. (2008) Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419
  * Fressin et al. (2013) Fressin, F., Torres, G., Charbonneau, D., et al. 2013, ApJ, 766, 81
  * Fulton et al. (2018) Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109
  * Gaidos et al. (2017) Gaidos, E., Kitzmann, D., & Heng, K. 2017, MNRAS, 468, 3418
  * Gandhi & Madhusudhan (2019) Gandhi, S., & Madhusudhan, N. 2019, arXiv e-prints, arXiv:1903.11603
  * Gomel et al. (2021) Gomel, R., Faigler, S., & Mazeh, T. 2021, MNRAS, 501, 2822
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93
  * Grimm et al. (2021) Grimm, S. L., Malik, M., Kitzmann, D., et al. 2021, ApJS, 253, 30
  * Haynes et al. (2015) Haynes, K., Mandell, A. M., Madhusudhan, N., Deming, D., & Knutson, H. 2015, ApJ, 806, 146
  * Hébrard et al. (2011) Hébrard, G., Ehrenreich, D., Bouchy, F., et al. 2011, A&A, 527, L11
  * Heiter et al. (2002) Heiter, U., Kupka, F., van’t Veer-Menneret, C., et al. 2002, A&A, 392, 619
  * Helling et al. (2019) Helling, C., Gourbin, P., Woitke, P., & Parmentier, V. 2019, A&A, 626, A133
  * Hoeijmakers et al. (2018) Hoeijmakers, H. J., Ehrenreich, D., Heng, K., et al. 2018, Nature, 560, 453
  * Hoeijmakers et al. (2019) Hoeijmakers, H. J., Ehrenreich, D., Kitzmann, D., et al. 2019, A&A, 627, A165
  * Hoeijmakers et al. (2020) Hoeijmakers, H. J., Cabot, S. H. C., Zhao, L., et al. 2020, arXiv e-prints, arXiv:2004.08415
  * Holman & Murray (2005) Holman, M. J., & Murray, N. W. 2005, Science, 307, 1288
  * Howard et al. (2012) Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15
  * Howell et al. (2014) Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398
  * Huang et al. (2020) Huang, C. X., Vanderburg, A., Pál, A., et al. 2020, Research Notes of the American Astronomical Society, 4, 204
  * Husser et al. (2013) Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
  * Ishizuka et al. (2021) Ishizuka, M., Kawahara, H., Nugroho, S. K., et al. 2021, arXiv e-prints, arXiv:2102.00211
  * Jenkins et al. (2016) Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in Proc. SPIE, Vol. 9913, Software and Cyberinfrastructure for Astronomy IV, 99133E
  * Jensen (2013) Jensen, E. 2013, Tapir: A web interface for transit/eclipse observability, Astrophysics Source Code Library, , , ascl:1306.007
  * John (1988) John, T. L. 1988, A&A, 193, 189
  * Jones et al. (2019) Jones, M. I., Brahm, R., Espinoza, N., et al. 2019, A&A, 625, A16
  * Jurgenson et al. (2016) Jurgenson, C., Fischer, D., McCracken, T., et al. 2016, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, EXPRES: a next generation RV spectrograph in the search for earth-like worlds, 99086T
  * Kitzmann et al. (2018) Kitzmann, D., Heng, K., Rimmer, P. B., et al. 2018, ApJ, 863, 183
  * Kossakowski et al. (2019) Kossakowski, D., Espinoza, N., Brahm, R., et al. 2019, MNRAS, 490, 1094
  * Kreidberg (2015) Kreidberg, L. 2015, PASP, 127, 1161
  * Kumar et al. (2019) Kumar, R., Carroll, C., Hartikainen, A., & Martin, O. A. 2019, The Journal of Open Source Software, doi:10.21105/joss.01143.  http://joss.theoj.org/papers/10.21105/joss.01143
  * Kurucz (1993) Kurucz, R. 1993, ATLAS9 Stellar Atmosphere Programs and 2 km/s grid. Kurucz CD-ROM No. 13. Cambridge, 13
  * Kurucz (2018) Kurucz, R. L. 2018, in Astronomical Society of the Pacific Conference Series, Vol. 515, Workshop on Astrophysical Opacities, 47
  * Langeveld et al. (2021) Langeveld, A. B., Madhusudhan, N., Cabot, S. H. C., & Hodgkin, S. T. 2021, MNRAS, 502, 4392
  * Lendl et al. (2020) Lendl, M., Csizmadia, S., Deline, A., et al. 2020, arXiv e-prints, arXiv:2009.13403
  * Levine et al. (2012) Levine, S. E., Bida, T. A., Chylek, T., et al. 2012, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8444, Status and performance of the Discovery Channel Telescope during commissioning, 844419
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python, Astrophysics Source Code Library, , , ascl:1812.013
  * Lissauer et al. (2011) Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, ApJS, 197, 8
  * Lothringer et al. (2018) Lothringer, J. D., Barman, T., & Koskinen, T. 2018, ApJ, 866, 27
  * Lothringer et al. (2020) Lothringer, J. D., Fu, G., Sing, D. K., & Barman, T. S. 2020, ApJ, 898, L14
  * Louden & Wheatley (2015) Louden, T., & Wheatley, P. J. 2015, ApJ, 814, L24
  * Luger et al. (2019) Luger, R., Agol, E., Foreman-Mackey, D., et al. 2019, AJ, 157, 64
  * Madhusudhan (2019) Madhusudhan, N. 2019, ARA&A, 57, 617
  * Madhusudhan et al. (2014) Madhusudhan, N., Knutson, H., Fortney, J. J., & Barman, T. 2014, Protostars and Planets VI, 739
  * Malik et al. (2019) Malik, M., Kitzmann, D., Mendonça, J. M., et al. 2019, AJ, 157, 170
  * McLaughlin (1924) McLaughlin, D. B. 1924, ApJ, 60, 22
  * Miller-Ricci Kempton & Rauscher (2012) Miller-Ricci Kempton, E., & Rauscher, E. 2012, ApJ, 751, 117
  * Millholland et al. (2017) Millholland, S., Wang, S., & Laughlin, G. 2017, ApJ, 849, L33
  * Morris et al. (2020) Morris, R. L., Twicken, J. D., Smith, J. C., et al. 2020, Kepler Data Processing Handbook: Photometric Analysis, Kepler Data Processing Handbook (KSCI-19081-003), ,
  * Morris (1985) Morris, S. 1985, ApJ, 295, 143
  * Nielsen et al. (2013) Nielsen, M. B., Gizon, L., Schunker, H., & Karoff, C. 2013, A&A, 557, L10
  * Nugroho et al. (2020) Nugroho, S. K., Gibson, N. P., de Mooij, E. J. W., et al. 2020, ApJ, 898, L31
  * Nugroho et al. (2017) Nugroho, S. K., Kawahara, H., Masuda, K., et al. 2017, AJ, 154, 221
  * Parmentier et al. (2018) Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, A&A, 617, A110
  * Pepper et al. (2007) Pepper, J., Pogge, R. W., DePoy, D. L., et al. 2007, PASP, 119, 923
  * Petersburg et al. (2020) Petersburg, R. R., Ong, J. M. J., Zhao, L. L., et al. 2020, AJ, 159, 187
  * Pino et al. (2020) Pino, L., Désert, J.-M., Brogi, M., et al. 2020, ApJ, 894, L27
  * Piskunov & Valenti (2017) Piskunov, N., & Valenti, J. A. 2017, A&A, 597, A16
  * Pollacco et al. (2006) Pollacco, D. L., Skillen, I., Collier Cameron, A., et al. 2006, PASP, 118, 1407
  * Ricker et al. (2014) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Proc. SPIE, 914320
  * Rodriguez et al. (2019) Rodriguez, J. E., Quinn, S. N., Huang, C. X., et al. 2019, AJ, 157, 191
  * Rossiter (1924) Rossiter, R. A. 1924, ApJ, 60, 15
  * Ryabchikova et al. (2015) Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr, 90, 054005
  * Safonov et al. (2017) Safonov, B. S., Lysenko, P. A., & Dodin, A. V. 2017, Astronomy Letters, 43, 344
  * Salvatier et al. (2016) Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Computer Science, 2, e55
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525
  * Seager & Sasselov (1998) Seager, S., & Sasselov, D. D. 1998, ApJ, 502, L157
  * Sedaghati et al. (2021) Sedaghati, E., MacDonald, R. J., Casasayas-Barris, N., et al. 2021, MNRAS, arXiv:2103.12858
  * Shporer (2017) Shporer, A. 2017, PASP, 129, 072001
  * Shporer et al. (2019) Shporer, A., Wong, I., Huang, C., et al. 2019, AJ, 157, 178
  * Sing et al. (2016) Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2016, Nature, 529, 59
  * Sing et al. (2019) Sing, D. K., Lavvas, P., Ballester, G. E., et al. 2019, AJ, 158, 91
  * Smette et al. (2015) Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77
  * Smith et al. (2012) Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000
  * Snellen et al. (2010) Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049
  * Stangret et al. (2020) Stangret, M., Casasayas-Barris, N., Pallé, E., et al. 2020, arXiv e-prints, arXiv:2003.04650
  * Stangret et al. (2021) Stangret, M., Pallé, E., Casasayas-Barris, N., et al. 2021, arXiv e-prints, arXiv:2104.12414
  * Stassun et al. (2017) Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136
  * Stassun et al. (2018a) Stassun, K. G., Corsaro, E., Pepper, J. A., & Gaudi, B. S. 2018a, AJ, 155, 22
  * Stassun & Torres (2016) Stassun, K. G., & Torres, G. 2016, AJ, 152, 180
  * Stassun & Torres (2018) —. 2018, ApJ, 862, 61
  * Stassun et al. (2018b) Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018b, AJ, 156, 102
  * Stock et al. (2018) Stock, J. W., Kitzmann, D., Patzer, A. B. C., & Sedlmayr, E. 2018, MNRAS, 479, 865
  * Stumpe et al. (2014) Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100
  * Stumpe et al. (2012) Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985
  * Talens et al. (2017) Talens, G. J. J., Albrecht, S., Spronck, J. F. P., et al. 2017, A&A, 606, A73
  * Telting et al. (2014) Telting, J. H., Avila, G., Buchhave, L., et al. 2014, Astronomische Nachrichten, 335, 41
  * Theano Development Team (2016) Theano Development Team. 2016, arXiv e-prints, abs/1605.02688.  http://arxiv.org/abs/1605.02688
  * Torres et al. (2010) Torres, G., Andersen, J., & Giménez, A. 2010, A&A Rev., 18, 67
  * Triaud (2018) Triaud, A. H. M. J. 2018, The Rossiter-McLaughlin Effect in Exoplanet Research, 2
  * Twicken et al. (2010) Twicken, J. D., Clarke, B. D., Bryson, S. T., et al. 2010, in Proc. SPIE, Vol. 7740, Software and Cyberinfrastructure for Astronomy, 774023
  * Valenti & Piskunov (1996) Valenti, J. A., & Piskunov, N. 1996, A&AS, 118, 595
  * von Essen et al. (2020) von Essen, C., Mallonn, M., Borre, C., et al. 2020, A&A, 639, A34
  * Wang et al. (2019) Wang, S., Jones, M., Shporer, A., et al. 2019, AJ, 157, 51
  * Weiss et al. (2018) Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48
  * Winn et al. (2010) Winn, J. N., Fabrycky, D., Albrecht, S., & Johnson, J. A. 2010, ApJ, 718, L145
  * Wong et al. (2020a) Wong, I., Benneke, B., Shporer, A., et al. 2020a, AJ, 150, 104
  * Wong et al. (2020b) Wong, I., Shporer, A., Daylan, T., et al. 2020b, AJ, 160, 155
  * Wong et al. (2016) Wong, I., Knutson, H. A., Kataria, T., et al. 2016, ApJ, 823, 122
  * Wong et al. (2020c) Wong, I., Shporer, A., Kitzmann, D., et al. 2020c, AJ, 160, 88
  * Wong et al. (2021) Wong, I., Kitzmann, D., Shporer, A., et al. 2021, AJ, in press, arXiv:2106.02610
  * Yan & Henning (2018) Yan, F., & Henning, T. 2018, Nature Astronomy, 2, 714
  * Yan et al. (2019) Yan, F., Casasayas-Barris, N., Molaverdikhani, K., et al. 2019, A&A, 632, A69
  * Yan et al. (2020) Yan, F., Pallé, E., Reiners, A., et al. 2020, A&A, 640, L5
  * Zhou et al. (2016) Zhou, G., Rodriguez, J. E., Collins, K. A., et al. 2016, AJ, 152, 136
  * Zhou et al. (2019) Zhou, G., Bakos, G. Á., Bayliss, D., et al. 2019, AJ, 157, 31

Figure 10 depicts a corner plot containing all the astrophysical parameters
fitted for in our joint analysis of the TESS light curve and ground-based
full-transit photometry (Section 2.4); for clarity, the limb-darkening
coefficients for each dataset are not shown. The values of the average
relative planetary flux $f_{p}$, planetary atmospheric brightness modulation
amplitude $A_{\mathrm{atm}}$, and stellar ellipsoidal distortion amplitude
$A_{\mathrm{ellip}}$ are given in parts-per-million. The phase offset in the
planetary phase curve $\delta$ is provided in degrees. Note the significant
correlations between the impact parameter $b$, scaled semimajor axis
$a/R_{\star}$, radius ratio $R_{p}/R_{\star}$, and $f_{p}$ — a consequence of
the grazing nature of the planetary transit.

Figure 10: Corner plot of parameters involved in the joint TESS and ground-
based light-curve fit.

Figure 11 shows the full-transit light curves collected as part of ground-
based followup observations, as described in Section 2.3. Each light curve is
labeled with the filter used. The best-fit transit model from the joint TESS
and ground-based photometric fit is shown in the bottom panels.

Figure 11: Ground-based light curves of TOI-1518 with full coverage of the
primary transits, collected as part of the TESS Follow-up Program. Top: the
photometry at the native time resolution (gray points) and binned (colored
points). Each panel is labeled by the respective bandpass. The binning
interval for the $B$-, $g^{\prime}$-, and $R$-band observations is 7 minutes;
a shorter 4-minute bin size is used for the higher-precision $R_{\rm long}$\-
and $I$-band transits. Bottom: binned, systematics-corrected light curves,
with the best-fit transit model from the joint TESS and ground-based
photometric fit (Table 1) plotted in black.