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

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\pm0.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\pm28$ 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\pm38$ K and a measured dayside brightness temperature of $3237\pm59$ K (assuming zero geometric albedo), TOI-1518b is a promising candidate for future emission spectroscopy to probe for a thermal inversion.


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. * 51 Pegasi b Fellow 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 , HD 2685b , 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 nearultraviolet 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;. 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.

OBSERVATIONS AND SYSTEM CHARACTERIZATION
This section describes our analysis of the available TESS photometry of TOI-1518, as well as highresolution 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.

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 fullframe 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) 1 .
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. , 2014Smith 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σ. However, there were residual longterm 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σ 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 starplanet 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.

Full-orbit Phase-curve Model
We fit the full-orbit phase curve with a composite flux model for the planet ψ p and the star ψ (e.g., Wong et al. 2020b,c;Wong et al. 2021): Here, A atm , A ellip , and A 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 δ 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φ)). 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.
The transits and secondary eclipse light curves (λ t and λ 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 atm cos(π + δ). Likewise, the hemisphere-averaged nightside flux is given by D n = f p − A atm cos(δ). 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: 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 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 . 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 ), planet-star radius ratio R p /R , 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 atm , A ellip , and δ to vary. We also include a uniform per-point uncertainty parameter σ k for each light-curve segment as a free parameter in order to ensure a reduced χ 2 value of one and retrieve realistic uncertainties on the astrophysical parameters. The median values of σ 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 = 0.0987 ± 0.0017 and well-constrained transit-shape parameters: b = 0.9103 ± 0.0065 and a/R = 4.231 ± 0.064. We detect the secondary eclipse with a depth of ∼380 ppm and a significant atmospheric phase-curve modulation with a semiamplitude of roughly 160 ppm. There is a nearly 5σ 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 ω 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σ upper limit of e < 0.01 (formally, e = 0.0031 +0.0047 −0.0022 ); the inclusion of e and ω as free parameters is strongly disfavored by the Bayesian Information Criterion (∆BIC = 16). The corresponding e cos ω and e sin ω values, which relate to offsets in the secondary eclipse timing and duration, respectively, are 0.0007 +0.0016 −0.0012 and −0.0005 +0.0030 −0.0061 . 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σ 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.

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) 2 . 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 .
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 (λ 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 gband (λ 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 (λ 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 long ; λ 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 (λ 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; λ cut−on = 500 nm). The light-curve data are available at ExoFOP-TESS. 3 The raw ground-based transit light curves are shown in the Appendix.
The follow-up light curves confirm that the TESSdetected 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 inter-polating 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 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 groundbased transit datasets translates to large relative uncertainties on the measured transit depth, exceeding 10% across all five visits. Nevertheless, we obtain R p /R values that are consistent with the measurement from fitting the TESS light curve alone at better than the 2σ level. Similarly, the five ground-based transit depths are mutually consistent to within 2σ, indicating an achromatic transit.

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 phasecurve 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.
The results of our joint fit are listed in Table 1. Figure 3 shows the binned, phase-folded, and systematicscorrected TESS light curve alongside the best-fit phasecurve 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 ± 0.000011 is measured to ∼1 s precision. We obtain a planet-star radius ratio of R p /R = 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 = 4.291 +0.057 −0.061 . The secondary eclipse depth is measured to more than 12σ significance: 364±28 ppm. The atmospheric phase-curve modulation has a semiamplitude of 160.4 ± 6.7 ppm. No significant phase shift in the planet's phase curve is measured, indicating that the location of maximum brightness on  (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).  the dayside hemisphere is well-aligned with the substellar point. The derived nightside flux is 43 ± 27 ppm. The ellipsoidal distortion signal from the host star is detected at 5.2σ significance, with a semiamplitude of 31.3 ± 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 , and f p , in addition to the typical degeneracy between b and a/R .

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 .

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 ∼300 s exposures. The extracted spectra have a signal-to-noise (S/N) of ∼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 O 2 and H 2 O 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 highresolution 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).

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 linelist (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 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 rota-tional 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 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 ± 6 km s −1 , which agrees with expectations for this spectral type (Nielsen et al. 2013). A fitted [Fe/H] of −0.1 ± 0.12 is low for a star hosting a hot Jupiter (Fischer & Valenti 2005); only ∼ 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 eff = 6910 ± 445 K and log g = 3.97 ± 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 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 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. (2017Stassun et al. ( , 2018a. We took the B T V T magnitudes from Tycho-2, the BV i magnitudes from APASS, the JHK S magnitudes from 2MASS, the W1-W4 magnitudes from WISE, the GG BP G 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 µm (see Figure 6).
We performed a fit using Kurucz stellar atmosphere models, with the free parameters being the effective temperature (T 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 χ 2 of 20.3 (with 12 degrees of freedom) and best-fit parameters T eff = 7300 ± 100 K, [Fe/H] = 0.0 ± 0.2, log g = 4.1 ± 0.2, and A V = 0.05 ± 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 bol = (6.52 ± 0.31) × 10 −9 erg s −1 cm −2 . Taking F bol and T eff together with the Gaia DR2 parallax (4.398 ± 0.033 mas) gives a stellar radius of R = 1.950 ± 0.048 R . In addition, we use R together with log g to obtain an empirical mass estimate of M = 1.79 ± 0.26 M , which is consistent with that calculated via the empirical relations of Torres et al. (2010) -M = 1.70 ± 0.12 M .

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-oftransit 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 ∼ 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σ 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σ error of 75 m s −1 , i.e. less than 2σ 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 M J , well within expectations for hot Jupiters.
To determine the systemic velocity, we compute the weighted mean of the measured RVs, −14.79 ± 0.06 km s −1 , which must be corrected for an instrumental offset of −0.87 ± 0.16 km s −1 , found from standard stars. We arrive at a systemic velocity V sys of −13.94 ± 0.17 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.

THE SPECTROSCOPIC TRANSIT
In this section, we describe the methods used to analyze the spectroscopic transit observation from EX-PRES. Cross-correlation was performed with the X-COR pipeline, previously used for atmospheric detections in  WASP-121b Ben-Yami et al. 2020) and MASCARA-2b (Hoeijmakers et al. 2020). Crosscorrelation 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 andPino 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 crosscorrelation 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.

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 out and then divided each individual spectrum by F out . Interstellar medium features were removed through division by F 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 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: 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 .

Spin-Orbit Misalignment
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 out . While F out is a good template for the out-oftransit 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 λ. 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 The orthogonal distance x ⊥ is determined by the position of the planet: x ⊥ (t) = x p (t) cos(λ) − y p (t) sin(λ) (7) x p (t) = a R sin(2πφ) (8) We used the physical model of Cegla et al. (2016) and the emcee sampler . Free parameters included the above four as well as Vsys, which returned a posterior distribution that was very similar to its prior Gaussian distribution. The parameters a/R and ip were constrained by Gaussian priors derived from the results of our TESS light-curve fit (Table 1).
Therefore, we can obtain independent constraints on a/R , λ, 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 < 12, 0 < λ < 2π and 0 < i p < π. We define Gaussian priors for the rotation speed and global offset: v sin i ∼ N (µ = 80, σ = 50) km s −1 , V sys ∼ N (µ = −14.5, σ = 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 = 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 λ and i p . Next, we rerun the MCMC using photometricallyderived priors on a/R 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 λ = 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 3 M J 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.

K p − V sys Analysis
Closer inspection of Figure 8 shows a faint, white trail spanning approximately ±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 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 sys . It is useful to determine K p and V 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 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 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.

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 H 2 . 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.

Detections
We detect Fe in the atmosphere of TOI-1518b at the 5.2σ level. We also report evidence of Fe + at the 3.4σ level. The PHOENIX model, which contains both species in addition to other atoms and ions, yields an enhanced atmospheric detection at 5.9σ confidence, while a combined Fe/Fe + model yields a 5.4σ detection. The Doppler shadow correction removes an artifact that otherwise biases detection significances. The K p and V 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 sys = −17 +3 −2 km s −1 . For Fe the values are K p = 157 +68 −44 km s −1 and V sys = −16 +2 −4 km s −1 , and for Fe + they are K p = 178 +41 −62 km s −1 and V sys = −18 +3 −3 km s −1 . Uncertainties correspond to the range of K p and V sys within a 1σ 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 sys found here is offset by about 3 km s −1 at the ∼ 1−2σ level. This blueshift may  Table 1, we predict a planetary RV semiamplitude of K p = 2πa sin i p /P = 217.4 ± 6.2 km s −1 . This value is higher than the K p measured from cross-correlation, but still consistent to within the 1σ 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 δ , which is equal to the peak value of the stacked CCF over all K p and V 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 ± 0.8) × 10 −4 and (1.5 ± 0.4) × 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 sc , of length hundreds of kilometers for hot Jupiters) (Madhusudhan et al. 2014). The excess absorption beyond the transmission spectrum continuum (R p /R ) 2 is approximately δ ≈ 2R p H/R 2 ; in other words, H ≈ δR p /2(R p /R ) 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 , which is typical for the H − continuum of an UHJ. We also take H sc ∼ 880 km, estimated from the measured T eq and log g, as well as taking the mean molecular weight as µ = 2.3 for an H 2 -dominated atmosphere; however µ may be affected by H 2 dissociation on the planet's dayside. While the mass is highly uncertain, we take the posterior median value of 1.4 M J in order to estimate log g. The resultant pressures corresponding to the absorption are P ∼ 6 × 10 −4 bar for Fe and P ∼ 2 × 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σ, 4.7σ, and 5.2σ 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σ when using exposures from only the first half of the transit. It drops to 1-2σ 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.

Temperature and Circulation
From the stellar radius, we can use the values of R p /R and a/R from our photometric analysis to straightforwardly compute the planet's radius and orbital semimajor axis: R p = 1.875 ± 0.053 R J and a = 0.0389 ± 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 hemisphereaveraged brightness temperature T p via the following relation (e.g., Shporer 2017): Here, the stellar and planetary flux spectra are given by F λ (T eff ) and F λ (T p ), respectively, and τ (λ) is the transmission function of the TESS bandpass. For simplicity, we assume that the planet's emission spectrum is wellmodeled 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 {T eff , [Fe/H], log g} that smoothly interpolates these values. The planet's brightness temperature can then be fit for using an MCMC routine, with Gaussian priors for T eff , [Fe/H], log g, and R p /R derived from the SED and TESS light-curve fits.
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 consis-tent with zero (e.g., Shporer et al. 2019;Wong et al. 2020b;Wong et al. 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 B and the recirculation efficiency . In this parametrization, 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σ) and = 0.5 ± 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.

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.

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 eq 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(Hoeijmakers et al. , 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σ), 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 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 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.

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): Here, the pre-factor α ellip is a function of the linear limbdarkening and gravity-darkening coefficients u and g for the host star: Similar to our treatment of the quadratic limbdarkening 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 ± 0.05 and g = 0.12 ± 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 ellip , a/R , i p , M , u, and g. We obtain a photometric mass estimate of M p = 4.8 +1.3 −1.1 M J . This value is significantly (2.3σ) larger than the RV-derived mass upper limit of 2.3 M J .
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 planetfacing 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.

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;. 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). 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 atm , and stellar ellipsoidal distortion amplitude A ellip are given in parts-per-million. The phase offset in the planetary phase curve δ is provided in degrees. Note the significant correlations between the impact parameter b, scaled semimajor axis a/R , radius ratio R p /R , and f p -a consequence of the grazing nature of the planetary transit. 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 -, and R-band observations is 7 minutes; a shorter 4-minute bin size is used for the higher-precision R 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.