A broadband thermal emission spectrum of the ultra-hot Jupiter WASP-18b

Close-in giant exoplanets with temperatures greater than 2,000 K (‘ultra-hot Jupiters’) have been the subject of extensive efforts to determine their atmospheric properties using thermal emission measurements from the Hubble Space Telescope (HST) and Spitzer Space Telescope1–3. However, previous studies have yielded inconsistent results because the small sizes of the spectral features and the limited information content of the data resulted in high sensitivity to the varying assumptions made in the treatment of instrument systematics and the atmospheric retrieval analysis3–12. Here we present a dayside thermal emission spectrum of the ultra-hot Jupiter WASP-18b obtained with the NIRISS13 instrument on the JWST. The data span 0.85 to 2.85 μm in wavelength at an average resolving power of 400 and exhibit minimal systematics. The spectrum shows three water emission features (at >6σ confidence) and evidence for optical opacity, possibly attributable to H−, TiO and VO (combined significance of 3.8σ). Models that fit the data require a thermal inversion, molecular dissociation as predicted by chemical equilibrium, a solar heavy-element abundance (‘metallicity’, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{M/H}}=1.0{3}_{-0.51}^{+1.11}$$\end{document}M/H=1.03−0.51+1.11 times solar) and a carbon-to-oxygen (C/O) ratio less than unity. The data also yield a dayside brightness temperature map, which shows a peak in temperature near the substellar point that decreases steeply and symmetrically with longitude towards the terminators.

Close-in giant exoplanets with temperatures greater than 2,000 K ("ultra-hot Jupiters") have been the subject of extensive efforts to determine their atmospheric properties using thermal emission measurements from the Hubble and Spitzer Space Telescopes 1;2;3 .However, previous studies have yielded inconsistent results because the small sizes of the spectral features and the limited information content of the data resulted in high sensitivity to the varying assumptions made in the treatment of instrument systematics and the atmospheric retrieval analysis 3;4;5;6;7;8;9;10;11;12 .Here we present a dayside thermal emission spectrum of the ultrahot Jupiter WASP-18b obtained with the NIRISS 13 instrument on JWST.The data span 0.85 to 2.85 µm in wavelength at an average resolving power of 400 and exhibit minimal systematics.The spectrum shows three water emission features (at >6σ confidence) and evidence for optical opacity, possibly due to H − , TiO, and VO (combined significance of 3.8σ).Models that fit the data require a thermal inversion, molecular dissociation as predicted by chemical equilibrium, a solar heavy element abundance ("metallicity", M/H = 1.03 +1.11 −0.51 × solar), and a carbon-to-oxygen (C/O) ratio less than unity.The data also yield a dayside brightness temperature map, which shows a peak in temperature near the sub-stellar point that decreases steeply and symmetrically with longitude toward the terminators.
We observed a secondary eclipse of WASP-18b with NIRISS/SOSS 13 as part of the JWST Transiting Exoplanet Community Early Release Science Program 14 .WASP-18b is a 10.4±0.4M J ultra-hot Jupiter on a 0.94 day orbit around a bright (J mag = 8.4) F6V-type star 15 .Our goals were to characterize WASP-18b's atmosphere and demonstrate the capabilities of JWST observations for exoplanets orbiting bright stars.We used the SUBSTRIP96 subarray mode (96 × 2048 pixels) to avoid saturation by minimizing the individual integration times.The SUBSTRIP96 mode covers the first spectral order between 0.85 and 2.85 µm, whereas the SUBSTRIP256 subarray also provides the shorter-wavelength measurements in the second spectral order.The time series spans 6.71 hours and consists of 2720 continuous integrations with 3 groups and 8.88 seconds per integration, delivering an integration efficiency of 67%.We used the F277W filter in the final ten integrations to check for contamination from background stars and found none.We observed for 2.83 hours before the eclipse, and continued for 1.70 hours after the eclipse.The observations captured 107 • of WASP-18b's orbit.Assuming it is tidally locked, the planet rotated by the same angle during the observation.
We analyzed the data using four independent pipelines: NAMELESS, nirHiss, supreme-SPOON, and transitspectroscopy (see Methods).Beginning from the raw uncalibrated data or stage 1 products, we performed custom reductions and extracted 1D spectra from each integration using either a fixed-width aperture (NAMELESS, supreme-SPOON, and transitspectroscopy) or an optimal extraction (nirHiss) technique.We put particular emphasis on the removal of 1/f noise (f is frequency), a signal with power spectrum inversely proportional to the frequency that is introduced through variations of the reference voltage as the detector is being read.Its removal requires careful treatment as the spectral trace covers most of the SUBSTRIP96 subarray (see Methods, Extended Data Fig. 1).Finally, we obtained spectrophotometric light curves by summing the observed flux within 408 spectral bins, each containing five pixel columns on the detector (Extended Data Fig. 2).We also produced a white-light curve by summing the spectrophotometric light curves over all wavelengths.All light curves show a sudden decrease in flux around the 1336th integration, simultaneous to a tilt event from one of the primary mirror's segments 16;17 .In the NIRISS data, this can be independently identified via a small but detectable morphological change in the spectral trace on the detector (Extended Data Fig. 3).We also observed small variations in the spectral trace morphology throughout the time series, mainly of its position and full width at half maximum, that are correlated with the measured flux.We detrended against these morphological changes at the fitting stage.
We analyzed the extracted white and spectrophotometric light curves by fitting the parameters for the secondary eclipse, the partial phase curve, and the systematics using ExoTEP 18 (see Methods).When fitting the white light curve, we allowed the semi-major axis and impact parameter to vary, constrained by Gaussian priors that were derived from an analysis of the full-orbit phase curve observed by the Transiting Exoplanet Survey Satellite (TESS; see Methods).We imposed a uniform prior on the mid-eclipse time and assumed a circular orbit.Those parameters were subsequently fixed when fitting the spectrophotometric bins (see Methods).The maximum signal-to-noise (SNR) ratio for a single pixel spectrophotometric light curve is 617 at 1.14 µm, with the SNR curve closely following the shape of the throughput-weighted stellar spectrum and reaching a minimum of 62 at 2.83 µm.All four reductions yield consistent results (Extended Data Fig. 4), with all resulting thermal emission spectra being consistent at less than one standard deviation on average.The residuals for each spectrophotometric light curve closely follow the expected 1/ √ n (n is the number of events) scaling of Poisson noise when binned in time, and the white light curve bins down to 5 ppm over one hour timescales (see Extended Data Fig. 5).
The secondary eclipse spectrum was created by collating the planet-to-star flux ratio values at mid-eclipse for all the wavelength channels (see Fig. 1).We then multiplied by a PHOENIX stellar spectrum model 19 produced using previously published parameters for WASP-18 (i.e., T eff = 6435 K, log g = 4. 35, and [Fe/H] = 0.1 20;21 ) to convert the dayside secondary-eclipse spectrum into the planet's thermal emission spectrum (see Fig. 1).For clarity, we also computed the brightness temperature spectrum, commonly used in planetary science, by calculating the blackbody temperature corresponding to the observed thermal emission in each wavelength bin (see Fig. 2).This transformation into brightness temperature facilitates identification of the various opacity sources by removing the large average slope caused by the behavior of the Planck emission across the NIRISS/SOSS wavelength range.
The observed brightness temperature spectrum shows strong deviations from a blackbody.It is dominated by the 1.4, 1.9, and 2.5 µm water emission features and a rise in brightness temperature shortwards of 1.3 µm.The rise in brightness is caused by the combined opacities of H − , TiO, and VO and we infer a combined detection significance of 3.8σ for these three species (Fig. 2).All molecular features appear in emission, indicating a thermal inversion (i.e., temperature increases with altitude, see also Fig. 3).The water features are consistent with a solar-composition atmosphere, as predicted by 1D radiative-convective models and 3D general circulation models (GCMs).They are strongly inconsistent with any high C/O or high metallicity scenarios 4 (Fig. 2b).This finding is further strengthened by the lack of detectable CO features at 1.6 and 2.4 µm, which should be the dominant species in a high-metallicity carbon-rich atmosphere (Fig. 2b).Using the free retrieval, we constrain the 3σ upper limit of the CO log mixing ratio to -2.42 (see below, Table 1, Extended Data Fig. 6).
Quantitatively, we infer an atmospheric metallicity of 0.82 +0.59 −0.37 times solar when fitting the NAMELESS reduction to a grid of self-consistent 1D radiative-convective models, and, consistently, 1.19 +1.22 −0.67 times solar when allowing for a free vertical temperature structure in the chemically-consistent retrievals.Both modeling approaches accounted for the thermal dissociation of water in the upper atmosphere and assumed chemical equilibrium.In both cases the best fits are obtained for sub-solar C/O values around 0.03-0.3,with the self-consistent models giving 3σ upper limits of 0.2, while the free-temperature-structure retrievals allows C/O values up to 0.6 at 3σ (Fig. 3), consistent with the upper limit from high-resolution dayside thermal emission observations 22 .We also assessed the effect of disequilibrium chemistry (see Methods) on the observed thermal emission and find the impact to be below 10 ppm due to the short chemical timescales in this hot atmosphere, justifying the assumption of thermochemical equilibrium models.In addition, we performed a free-chemistry atmospheric retrieval 23;24;25 , including the effects of thermal dissociation 9 , and inferred a H 2 O deep atmospheric log mixing ratio of −3.23 +0. 450.29 , consistent with the models assuming chemical equilibrium (Fig. 2c).We identify a strong thermal inversion with a temperature increase of 500 K in the middle atmosphere from 1 bar to 0.01 bar, which corresponds to the pressure range covered by the contribution functions (Fig. 3a).Our bestfit radiative-convective model provides strong evidence that the temperature inversion is caused by the absorption of stellar light by TiO (see Extended Data Fig. 7).At first sight this can seem at odds with high-spectral resolution observations that have detected other species able to create thermal inversions, such as atomic iron 10 , but have had trouble detecting TiO 26 .This tension is easily solved when considering that both TiO and water thermally dissociate in the upper atmospheric layers of ultra-hot Jupiters.Our observations, on the other hand, are sensitive to deeper layers of the atmosphere close to the infrared photosphere, which extends from 0.01 down to 1 bar (see contribution function on Fig. 3a), in the region where molecules such as water and TiO recombine (see Extended Data Fig. 6).Even though our model also predicts that iron can produce a thermal inversion, its near-constant abundances means that inversions due to iron happen at pressures lower than 1 millibar and not where we detect the main thermal inversion.
The precise constraints on the atmospheric metallicity and C/O enable us to investigate possible formation scenarios of WASP-18b.;28 , indicates that accretion of protoplanetary gas, rather than rocky or icy planetesimals, dominated the planet's late-stage formation.The mass-metallicity trend derived from solar system planets 29;30;31;32 predicts that the metallicity decreases as the mass of the planet increases, approaching the composition of the star for the most massive planets.Our finding of solar metallicity, three times lower than that of Jupiter, is consistent with this trend, given WASP-18b's mass of 10.4 M J .The low C/O ratio furthermore disfavors forming WASP-18b beyond the CO 2 ice line followed by an inward migration after the disk has dispersed, as gas accretion in that region would have led to high C/O values 33 .Detailed spectroscopic observations of the 4.5 µm CO feature, which is found within the spectral range of JWST NIRSpec G395H, could lead to a more stringent constraint on the C/O ratio, and, thus, on WASP-18b's formation and migration history.A detailed interpretation of the atmospheric C/O ratio of WASP-18b would require knowledge of the C/O ratio of the host star.This was recently found to be significantly sub-solar (C/O=0.23±0.05) 28based on high-resolution spectroscopy.However, stellar C/O measurements are especially challenging due to stellar model inaccuracies and weak/blended absorption lines 34 , so further confirmation is warranted.
Another possible formation scenario for WASP-18b is through collapse of the disk from gravitational instability 35 with a disk-free migration.This process leads to an atmospheric metallicity and C/O dictated by the local disk composition, and is expected to result in planets with stellar-to-super-stellar metallicities and sub-stellar-to-stellar C/O 36 , in agreement with our results.
;38;39 ).We begin this analysis with the systematics-corrected white-light curve.We performed two independent applications of the method, both enforcing positive flux contribution from visible locations on the planet.We find two brightness map solutions which fit the data similarly well (Fig. 4).We convert brightness maps to brightness temperature maps assuming the star emits as a blackbody at 6432 ± 48 K 20 and R p /R s = 0.09783±0.00028(Table 2).The first solution (blue model) shows a brightness temperature plateau stretching from approximately -40 • to +40 • of longitude relative to the substellar point, with a virtually constant latitudinally averaged brightness temperature of 2781 +25 −13 K.The second solution (red model) shows a more concentrated hot spot at the substellar point with a maximum brightness temperature of 2925 +16 −18 K and a consistent decrease in temperature both eastward and westward of the substellar point.Both solutions consistently reveal a steep temperature drop with longitude toward the terminators, with the inferred brightness temperature falling to 1686 +70 −27 K at the western terminator and 1869 +50 −14 K at the eastern terminator (blue model), and neither shows any significant shift of the brightest region away from the substellar point.This is consistent with what was measured from the HST phase curve of WASP-18b 40 .The high temperatures covering most of the dayside in both solutions, along with the steep decrease in temperature near the limbs, are consistent with the atmospheric retrieval results (Fig. 3d).
Beyond the terminators and leading to the nightside, we infer a continued drop in the thermal emission.Our JWST observations have the sensitivity to probe part of the night side because the planet rotates by 107 • during the time series, providing a view of up to 62.5 • of the night side east of the eastern terminator at the beginning and ending 44.5 • west away from the western terminator.The lack of a hot-spot offset and the large center-to-limb brightness temperature contrast suggest heat transport by winds moving radially away from the sub-stellar point and toward the nightside, rather than redistributing heat to the nightside through the formation of an equatorial jet 41;42 .Lorentz forces are expected to play an important role in atmospheric dynamics of ultra-hot Jupiters, due to their high dayside temperatures 42;43;44 .Thermally ionized alkali metals coupled to an internal dynamo-driven planetary magnetic field interact with the neutral species and are expected to prevent the formation of an eastward equatorial jet 45 .By approximating the effects of the Lorentz force as a locally calculated magnetic drag force in GCMs 46 (see Methods), we find that the observed white-light curve is best explained by an internal planetary field strength larger than 5 G, as this field strength is sufficient to prevent a discernible longitudinal shift of the hot spot from the substellar point (Fig. 4).This is further confirmed by a separate GCM considering spatially uniform drag timescales, for which we find that the case with the highest drag strength (τ drag = 103 s) produces white-light curves that best fit the observations (Fig. 4, Extended Data Fig. 8).Furthermore, self-consistent magnetohydrodynamic (MHD) models of ultra-hot Jupiters considering the response of the magnetic field to the circulation have shown the possibility of time variability in the longitudinal hot spot offset, oscillating between the western and eastern hemispheres over timescales of 10-100 days 47;48;49 , but additional observations are needed to test this possibility.
The large wavelength coverage and high spectral and photometric precision of JWST's NIRISS/SOSS mode present many opportunities for the study and detailed characterization of atmospheric processes through thermal emission spectroscopy.Furthermore, planets with high signal to noise eclipses such as WASP-18b allow for the three-dimensional mapping of their atmospheres to retrieve the temperature structure across the dayside as well as variations in properties such as molecular abundances 38;39 .JWST will enable these measurements for most bright transiting exoplanets, giving rise to the possibility of studying the dynamics and chemistry of a wide range of exoplanets directly from secondary eclipse observations.Methods NIRISS/SOSS Reduction and Spectrophotometric Extraction.We perform four separate reductions of the NIRISS/SOSS eclipse observations of WASP-18b using the NAMELESS, nirHiss1 , supreme-SPOON2 , and transitspectroscopy 3 pipelines for inter-comparison of individual reduction steps 50 .All pipelines are built around the official STScI jwst reduction pipeline 4 with the addition of custom correction steps for systematics such as 1/f noise, zodiacal background, and cosmic rays.Reductions are performed from either the raw uncalibrated data (NAMELESS, nirHiss, and supreme-SPOON) or stage 1 (transitspectroscopy) products up to the extraction of the spectrophotometric light curves.

NAMELESS Reduction
We use the NAMELESS pipeline 50 to reduce the WASP-18 b observations from the uncalibrated data products through spectral extraction.We used the jwst pipeline version 1.6.0,CRDS (Calibration Reference Data System) version 11.16.5, and CRDS context jwst 0977.pmap for the reduction.First, we go through all steps of the jwst pipeline Stage 1, with the exception of the dark current step.We skip the dark subtraction step to avoid introducing additional noise due to the lack of a high-fidelity reference file.After the ramp-fitting step, we go through the assign wcs, srctype, and flat field steps of the jwst pipeline Stage 2; we skip the background step and apply our own custom routine for handling the background.We skip the pathloss and photom steps, as an absolute flux calibration is not needed.We perform outlier detection by computing the product of the second derivatives in the column and row directions for all frames 50 .We divide the frames into windows of 4×4 pixels, where we then compute the local median and standard deviation of the second derivative.We flag all pixels that are ≥ 4σ away from the window median.Furthermore, we flag pixels that show null or negative counts.All flagged pixels are set equal to the local median of their window.We correct for background systematics using the following routine.First, we identified section (x∈ [5,400], y∈[0,20]) as the region of the SUBSTRIP96 sub-array where the contribution from the spectral orders to the counts is minimal.Next, we compute the scaling factor between the median frame and the model background provided on the STScI JDox User Documentation5 within the aforementioned region.We consider the 16 th percentile of the distribution as the scaling value and subtract the scaled background from all integrations.We pay close attention to the 1/f correction for these observations, as the magnitude of the spectral trace variation is highly wavelength-dependent in the secondary eclipse.Therefore, we consider all columns independently when scaling the trace to compute the 1/f noise.Furthermore, we treat this correction in two parts as we observe a tilt event 16;17 around the 1336th integration (Extended Data Fig. 3), possibly due to a sudden movement in one of the primary mirror's segments, resulting in a change in the morphology of the trace that manifests as a sudden decrease in flux.First, we compute the median columns c before and after the tilt event; we use integrations 300-900 and 1350-1900.Then, we define a given column j and row k at an integration i as the sum of the scaled median column m i,j cj,k and the 1/f noise n i,j = c i,j,k − m i,j cj,k .Using the errors i,j,k returned by the jwst pipeline, we minimize the χ 2 to find the value of the 1/f noise that best fits the observed column, such that and We then subtract this value from all columns and integrations.We set the error i,j,k to ∞ for all pixels that have non-zero data quality flags returned by the jwst pipeline such that they are not considered in the fit of the 1/f noise.We also set the errors to ∞ in the region x ∈ [76, 96], y ∈ [530, 1350] of the detector, where a portion of the second spectral order is visible.This is appropriate treatment as the 1/f noise scales independently across each order due to the difference in wavelength coverage.After correction of the 1/f noise, we trace the location of Order 1 on the detector by computing the maximum of the trace convolved with a Gaussian filter for all columns.We further smooth the positions of the trace centroids using a spline function.Finally, we perform a box spectral extraction of the first order using the transitspectroscopy.spectroscopy.getSimpleSpectrumroutine with an aperture diameter of 30 pixels.We follow two steps for trace identification.First, we use the (x,y) position of the trace from the "jwst niriss spectrace 0022.fits"reference file, with y-values offset by ∼25 pixels, to identify Order 2. We mask this region, such that it does not contaminate the trace identification or background routine later on.Next, we use the nirHiss.tracing.maskmethod edges function.This technique identifies the edges of Order 1 using a canny edge detection method from scikit-image, an open-source image processing package 52 .This method uses the derivative of a Gaussian function in order to identify regions with the maximum gradient.From this step, the potential edges are narrowed down to 1pixel curves along the maxima.This results in an image where the outline of Order 1 is presented.We identify the median location along the column from the top and bottom edges of Order 1, and smooth the trace by fitting a 4 th -order polynomial.We use the trace to mask the location of Order 1 when stepping through the nirHiss background routines.For background treatment, we follow a similar method presented in Feinstein et al. ( 2023) 50 , namely we identify a region without significant contamination from the spectral trace, and scale this region to the same region on the model background on the STScI JDox User Documentation.We use the region x ∈ [4, 250] and y ∈ [0, 30] and find an average scaling factor of 0.6007.We apply this scaling factor to the model background and subtract it from all integrations.Next, we remove 1/f noise in a similar manner to transitspectroscopy and scale this 1/f noise treatment to the out-of-eclipse integrations (0 -1250 and 1900 -2500).We remove cosmic rays using the L.A. Cosmic technique 53;54 .Finally, we extract the spectra using the optimal extraction routine, which is a robust means to simultaneously remove additional bad pixels/cosmic ray events while placing non-uniform weighting on each pixel in order to negate distortion produced by the spatial profile 55 .We use a normalized median image to best capture the unique NIRISS/SOSS spatial profile.

supreme-SPOON Reduction
We follow a similar approach with supreme-SPOON as presented in Feinstein et al. (2023) 50 .We start from the raw uncalibrated data files, which we downloaded from the MAST archive, and process them through the supreme-SPOON Stage 1, which performs the detector level calibrations including superbias subtraction, saturation flagging, jump detection, and ramp fitting.As with the previous pipelines, we do not perform any dark current subtraction.supreme-SPOON additionally treats the 1/f noise at the group level.This is done by subtracting a median stack of all in-eclipse integrations, scaled to the flux level of each individ-ual integration via the white light curve, to create a difference image revealing the characteristic 1/f striping.A column-wise median of the n th difference image is then subtracted from the corresponding integration.supreme-SPOON also removes the zodiacal background signal directly before the 1/f correction step via scaling the SUBSTRIP96 SOSS background model provided by the STScI to the flux level of each integration, as is described in Feinstein et al. ( 2023) 50 .We then pass the Stage 1 processed files through supreme-SPOON Stage 2, which performs additional calibrations such as flat fielding and hot pixel interpolation.We extract the stellar spectra at the pixel level using a simple box aperture extraction with a width of 30 pixels centered on the order 1 trace, as the dilution resulting from the order overlap with the second order has been shown to be negligible 56;57 .The y-pixel positions of the trace are determined via the edgetrigger algorithm 57 .We find that the extracted trace positions match with those measured during commissioning and included in the spectrace reference file, and we therefore use the default JWST wavelength solution.
transitspectroscopy Reduction We follow a similar approach adopted by the transitspectroscopy pipeline discussed in Feinstein et al. ( 2023) 50 .This reduction starts from the * rateints.fitsfiles that were processed by the jwst pipeline from STScI.We scaled the zodiacal background model provided on STScI JDox User Documentation to the observed two-dimensional spectra in the box delimited by pixels x ∈ [10, 250] and y ∈ [10, 30].
The scaled background is then subtracted from each integration.In summary, the procedure to remove the 1/f noise is as follows: We take the median of all integrations and subtract it from each integration, which leaves predominantly the 1/f noise.We then take the column-by-column median of this residual noise, considering only the pixels that are 20 pixels away from the center of the trace, and assume it is representative of the structure of the 1/f noise of the images.These values are then subtracted from each column.For the spectral extraction, we used the transitspectroscopy.spectroscopy.getSimpleSpectrumroutine with an aperture width of 30 pixels.We removed the outliers of the extracted spectra caused by cosmic rays or deviating pixel by taking the combined median of all spectra and flagging outlier points that deviate by more than 5σ from this median spectrum.The flagged wavelength bins are then corrected by taking the mean of the neighboring bins.
Spectrophotometric Light Curve Creation.From the aforementioned stellar spectral extraction pipelines, we create and fit models to the spectrophotometric light curves F (t, λ).The light curves are composed of three distinct signals: the planetary flux throughout partial phase curve and eclipse F p , the stellar flux F * , and the systematics S, which we model as a function of time t and wavelength λ via equation 3.
As the main scientific quantity of interest is the planetary signal F p (t, λ, θ), where θ represents the planetary orbital parameters, we aim to properly characterize and correct for the stellar flux and systematics.Light curve fitting is performed in two separate steps: (I) we fit the white light curve (II) we run individual fitting for each spectrophotometric bin.The values of the orbital parameters retrieved from the white light curve are fixed in the spectrophotometric light curves.The following sections describe our treatment of the planetary flux F p (t, λ, θ), the stellar variability F * (t, λ), and the systematics S(t, λ).
Light Curve Component I: Planetary Flux Despite the fact that the main target of these observations was the secondary eclipse of WASP-18 b, we also capture a portion of its phase curve during the before-and after-eclipse baseline.Over the course of the observations, the planet rotates 107 • , revealing significant information on the spatial distribution of its atmosphere.To model the planetary flux in time, we consider a second-order harmonic function where 0 ≤ f (t, θ) ≤ 1 is the time-dependent, visible fraction of the planetary projected disk as a function of the orbital parameters θ and P is the orbital period.The harmonics consist of a term describing the semi-amplitude F n of the planetary flux variation, as well as the time t n where the harmonic reaches its maximum.The visible fraction is computed using the normalized secondary eclipse light curve modeled with the batman python package 58 .The second-order harmonic function provides sufficient precision for the noise floor of the JWST 59 .We fit for the orbital parameters that dictate the shape and duration of the eclipse: the time of superior conjunction T sec (U[2459802.28,2459802.48]), the impact parameter b (N [0.360, 0.026 2 ]), as well as the semimajor axis to stellar radius ratio a/R * (N [3.496, 0.029 2 ]).We assume the orbit to be circular, which is justified by TESS analysis as it finds a strong Bayesian Information Criterion 60 (BIC) and Aikake Information Criterion 61 (AIC) preference for the non-eccentric orbit.Given the close proximity of WASP-18b to its host star, strong tidal interactions lead to the ellipsoidal deformation of the planet and its host.Past studies have shown that this deformation for WASP-18b is of ∼2.5×10 −3 R p , leading to a variation of the flux of order unity ppm, and is thus negligible 40;62 .Finally, near the lower end of the first order (0.85µm), there is also a contribution to the observed flux from the stellar light reflected by WASP-18b.However, the upper limits of the geometric albedo (A g <0.048 63 and A g = 0.025±0.027 64) obtained from TESS correspond to a reflected light contribution of <35 ppm near 0.8µm.We therefore do not consider a specific term for the reflected light and instead assume that this will be fit by the second-order harmonic function.
Light Curve Component II: Stellar Variability We consider three phenomena that can lead to temporal changes in the observed stellar flux: stellar activity, A, ellipsoidal variations, E, and Doppler boosting, D. Stellar activity, generally caused by the presence and movement of star spots on the stellar hemisphere visible to the observer, leads to variations in the observed stellar spectrum on a timescale that is of the order of the stellar rotation period.Past observations of WASP-18 have constrained this period to be P rot ∼ 5.5 days 65 .Despite this relatively short period with respect to stars of similar physical properties (e.g., the effective temperature, T eff , and luminosity, L ) 66 , the star shows abnormally low activity in the UV and X-ray domains, possibly due to tidal interactions with WASP-18b disrupting its dynamo 67;68;69 .As we expect the rotational modulation to be on a relatively long timescale compared to our observations and its amplitude to be quite low, we do not directly fit for this term and instead assume it to be handled by the systematics model.As for the ellipsoidal variation and Doppler boosting, they are both caused by the influence of WASP-18b on its host star.Although the ellipsoidal deformation of WASP-18b leads to a negligible impact on the phase curve, the same is not true for its host.The stellar ellipsoidal effect, with maxima fixed at quadrature when the projected area is at its highest, is found to be of semi-amplitude 172.2 ppm from the TESS analysis.Following previous analyses of HST observations 40 , we consider ellipsoidal variation to be achromatic and fix its amplitude to the TESS value for the full NIRISS/SOSS wavelength range.The Doppler boosting effect is a result of the Doppler shift of the stellar spectral energy distribution as its radial velocity varies throughout its orbit.We fix this amplitude to 21.8 ppm, with the maximum at phase 0.25, as done in Arcangeli et al. (2019) 40 .The observed stellar flux is therefore described as the sum of the ellipsoidal variation and Doppler boosting to the mean stellar flux in time F * (λ).
Light Curve Component III: Systematics Model The white and spectrophotometric light curves show two distinct important systematics: a sudden drop in flux caused by a tilt event shortly after the beginning of full-eclipse as well as high-frequency variations in the flux throughout the observations caused by small changes in the morphology of the trace.We track the trend of these systematics throughout the observations by performing incremental principal components analysis (IPCA) with the open-source scikit-learn 70 package on the processed detector images (Extended Data Fig. 3).The first principal component is the tilt event, which we use to determine the exact integration where it occurs.We handle the tilt event in the white and spectrophotometric light curves by fitting for an offset in flux for the data after the 1336th integration.We also observe two principal components analogous to the y-position and full width at half maximum (FWHM) of the trace in time.We find that these two components are correlated to short-frequency variations in the light curves and therefore detrend linearly against them at the fitting stage.We find that, despite having a lower variance than the y-position, the variation of the FWHM has a larger effect on the light curves when using box extraction.Finally, we fit for a linear trend in time to account for any further potential stellar activity and instrumental effect.
We note that considering a second order polynomial trend or higher in time for the systematics results in significant correlation with the partial phase curve information.However, a linear trend systematics model has been found sufficient to fit past NIRISS/SOSS observations (Feinstein  et al. (2023)  50 , Radica et al. submitted).Furthermore, we find that the curvature around secondary eclipse increases monotonically with wavelength, which is expected from the planetary signal and inconsistent with stellar activity as well as instrumental effects.
Light Curve Fitting Light curve fitting is performed on the extracted spectrophotometric observations using the ExoTEP framework 18 .With the orbital parameter values constrained from the white light curve (T sec = 2459802.381867+0.000092 −0.000089 , a/R * = 3.483 +0.021 −0.020 , b = 0.340 +0.016 −0.018 ), we solely fit for the planetary flux and systematics for each spectrophotometric light curve.We chose a resolution of 5 pixels/bin for our spectrum, corresponding to 408 spectrophotometric bins, to mitigate potential correlation between wavelengths in the atmospheric retrievals as pixels in the spectral direction are not independent.All fits for the 408 bins are then done independently.Retrievals are performed using the Affine Invariant Markov chain Monte Carlo Ensemble sampler emcee 71 , using 20,000 steps and 4 walkers per free parameter.The first 12,000 steps, 60% of the total amount, are discarded as burn-in to ensure that the samples are taken after the walkers have converged.The samples from the white and spectrophotometric light curves are used to produce two science products: the detrended white light curve and dayside thermal emission spectrum.The detrended white light curve is obtained by dividing out the best-fit systematics model and subtracting the stellar variability from the light curve to isolate the planetary signal.For the dayside thermal emission, the median values and uncertainties of F p (T sec , λ) are computed from the samples of the parameters of equation 4.

TESS Phase Curve Analysis
During the TESS Primary Mission, the WASP-18 system was observed in Sectors 2 and 3 (2018 August 22 to October 18).The full-orbit phase curve was analyzed in several earlier publications 63;72 , which reported a robust detection of the planet's secondary eclipse and high signal-to-noise measurements of the planet's phase curve variation and signals corresponding to the host star's ellipsoidal distortion and Doppler boosting.During the ongoing TESS Extended Mission, the spacecraft reobserved WASP-18 in Sectors 29 and 30 (2020 August 26 to October 21).We carried out a follow-up phase curve analysis of all four Sectors' worth of TESS data, following the same methods used previously.
The data from the TESS observations were processed by the Science Processing Operations Center (SPOC) pipeline, which yielded near-continuous light curves at a 2-minute cadence.In addition to the raw extracted flux measurements contained in the simple aperture photometry (SAP) light curves, the SPOC pipeline also produced the pre-search conditioning (PDC) light curves, which have been corrected for common-mode instrumental systematics trends that are shared among all sources on the corresponding detector.We used these PDC light curves for our phase curve analysis.After dividing the light curves into individual segments that are separated by the spacecraft's scheduled momentum dumps, we fit each segment to a combined phase curve and systematics model.The astrophysical phase curve model consists of two components describing the planetary and stellar fluxes: The Doppler boosting D and ellipsoidal distortion E semi-amplitudes are defined as before.Here, the planet's orbital phase is defined relative to the mid-transit time T 0 : φ = 2π(t−T 0 ).The planet's phase curve contribution has a single mode with a semi-amplitude of F atm and oscillates around the average relative planetary flux fp .The transit and secondary eclipse light curves φ t and φ e were modeled using batman with quadratic limb-darkening.In this parameterization, the secondary eclipse depth (i.e., dayside flux) and nightside flux are fp + F atm and fp − F atm , respectively.For the systematics model, we used polynomials in time and chose the optimal polynomial order for each segment individually that minimized the BIC.
We used ExoTEP to calculate the posterior distributions of the free parameters through a joint MCMC fit of all light curve segments.In order to reduce the number of free parameters in the joint fit, we first carried out fits to smaller groups of light curve segments corresponding to each TESS Sector and then divided the light curve segments by the best-fit systematics model.In the final joint fit of the systematics-corrected TESS light curve, no additional systematics parameters were included.We accounted for time-correlated noise (i.e., red noise) by fixing the uncertainty of all datapoints within each Sector to the standard deviation of the residuals, multiplied by the fractional enhancement of the average binned residual scatter from the expected Poisson noise scaling across bin sizes ranging from 30 minutes to 8 hours 72 .In addition to the phase curve parameters described above ( fp , F atm , D, E), we allowed the mid-transit time T 0 , orbital period P , relative planetary radius R p /R * , scaled orbital semi-major axis a/R * , impact parameter b, and quadratic limb-darkening coefficients u 1 and u 2 to vary freely.
The results of our TESS phase curve fit are listed in Table 2.The updated values for the orbital ephemeris and transit shape parameters are statistically consistent with the published results from the previous TESS phase curve analyses 63;72 , while being significantly more precise.We used the median and 1σ uncertainties of a/R * and b as Gaussian priors in the NIRISS/SOSS white light curve fit.We also used the median values of P , R p /R * , D, and E as fixed parameters for the NIRISS/SOSS white and spectrophotometric light curves fit.We obtained a secondary eclipse depth of 357 ± 14 ppm and a nightside flux that is consistent with zero.All three phase curve amplitudes were measured at high signal-to-noise: F atm = 177.5 ± 5.7 ppm, D = 21.8 ± 5.2 ppm, and E = 172.2± 5.6 ppm.
Secondary Eclipse Mapping.To perform eclipse mapping, we use both ThERESA 39 and the methods of Mansfield et al. (2020) 38 , which are separate implementations of the same process introduced by Rauscher et al. (2018) 37 when applied to white light curves, to cross-check our results.First, we generate a basis set of light curves from spherical harmonic maps 73;74 with degree l ≤ l max , then transform these light curves to a new, orthogonal basis set of "eigencurves" using principal component analysis (PCA).Each eigencurve corresponds to an "eigenmap," the planetary flux map which, when integrated over the visible hemisphere at each exposure time, generates the corresponding eigencurve.
We then fit the white light curve with a linear combination of a uniform-map light curve, the N most informative (largest eigenvalue) eigencurves, and a constant offset term to adjust for the fact that the observed planetary flux during eclipse (when the planet is entirely blocked by the star) must be equal to 0 and to allow for adjustments to light-curve normalization.Since the eigenmaps represent differences from a uniform map, it is possible to recover a fit which contains regions of negative planet emission.This is physically impossible, so we impose a positivity constraint on the total flux map.While the eigenmaps are mathematically defined across the entire planetary sphere, our observations only constrain the portion of the planet that is visible during the observation, so we only enforce the flux positivity condition in the visible region of the planet.We test all combinations of l max ≤ 6 and N ≤ 8 using a least-squares minimization and select the optimal values by minimizing the BIC.
We find that the fit with the lowest BIC to the broadband light curve has l max = 5 and N = 5.However, the fit with l max = 2 and N = 5 was only slightly less preferred, so here we explore the inferred brightness distribution from both solutions.Figure 9 shows the resulting light curve for the l max = 5, N = 5 solution after sequential subtraction of each eigencurve.The preference for a fit with five eigencurves is driven by the residuals in ingress and egress, which can be seen by eye and are not sufficiently corrected for with a uniform map.Including the uniform-map light curve and the constant term, the fit thus contained a total of seven free parameters.We used a Markov Chain Monte Carlo (MCMC) procedure to estimate parameter uncertainties.For the analysis following Mansfield et al. (2020) 38 , we test for convergence of the MCMC by ensuring that the chain length is 50 times the autocorrelation timescale, while for the analysis using ThERESA 39 we use the Gelman-Rubin convergence test 75 .
The resulting weights of each eigencurve are then applied to the corresponding eigenmaps to generate a flux map of the planet.We convert the star-normalized flux map to brightness temperature by assuming the planet and star are blackbodies emitting at the NIRISS/SOSS throughputweighted mean wavelength of 1.6 microns 37 .We estimate temperature map uncertainties by computing a subsample of maps from the MCMC posterior distribution and calculating 68.3%, 95.5%, and 99.7% quantiles at each location.
Figure 4 shows the resulting broadband brightness temperature map for the l max = 5, N = 5 case and longitudinal brightness temperature profiles for both the l max = 5, N = 5 and l max = 2, N = 5 cases, calculated by averaging meridian flux at each longitude weighted by cos (latitude) 2 .Additionally, we compare the equatorial slices to predictions from several GCMs (see General Circulation Models section).We note that not all structures on a planetary map will leave an observable signature in a secondary eclipse light curve.When comparing GCMs to secondary eclipse maps, it is important to only compare the components of GCM maps which can be physically accessed with eclipse mapping.Therefore, we use the methods of Luger et al. (2021) 76 to separate each GCM map into the "null space", or components that are inaccessible to eclipse mapping observations, and the "preimage", or components that are accessible through mapping.Figure 4 compares the longitudinal temperature trends from only the preimage of each GCM to the observed map.We find that both map solutions agree on a steep gradient in temperature near the limbs, which is well matched by GCM predictions.Additionally, both maps show a temperature distribution roughly symmetrical in longitude about the substellar point.However, the two maps disagree on the exact shape of the brightness distribution.The l max = 5, N = 5 map shows an extended hot plateau region of roughly constant temperature from -40 • to +40 • in longitude, while the l max = 2, N = 5 map shows a more concentrated hot spot with a steady decrease in temperature away from the substellar point.As these maps both fit the data with similar BIC, the current data do not give us the necessary precision to determine which solution represents the true temperature distribution of WASP-18b.Future observations at higher precision may distinguish between these two modes of solutions.
To test our ability to constrain latitudinal structures, we performed two eclipse-mapping fits: the eigenmapping fit presented above and a fit where the initial basis set of maps is a longitudinal Fourier series that is constant with latitude (see Extended Data Fig. 10).Both methods retrieve similar longitudinal temperature structures, with steep gradients near the limb and a extended hot plateau.However, the constant-with-latitude map is also able to fit the data well, indicating a lack of constraints on latitudinal features within the uncertainties on the data.This is not unexpected as the relatively low impact parameter (b = 0.36) of WASP-18b results in a lower amount of latitudinal information contained in the secondary eclipse signal.Further observations of WASP-18b, or of planets with higher impact parameter, may enable us to pull latitudinal signals out of the noise.
Our eclipse mapping assumes the orbital parameters of the system are precisely known relative to data uncertainties, a safe assumption with Spitzer data 37 .With JWST, data quality may be high enough that uncertainties on orbital parameters impart significant uncertainty on the mapping results.To test this, we ran analyses with impact parameter, orbital semi-major axis, and eclipse time fixed to values ±1σ.In some cases, this led to a "hotspot" model like the red one in Figure 4 being preferred over a "plateau" model like the blue one, which is unsurprising given their similar statistical preference.However, all resulting maps were well within the uncertainties of one of the two models presented.We note that, while the eccentricity is kept fixed to zero throughout this analysis, considering an eccentric orbit would allow to first order for variations in mid-eclipse time and eclipse duration 77 .Past radial velocity measurements have found a small but non-zero eccentricity for WASP-18b of e = 0.0076 ± 0.0010 78 , corresponding to an offset of the time of mid-eclipse of 9 seconds, as well as a difference of 120 seconds between the transit and eclipse durations.These differences in eclipse timing and duration are of the same magnitude as those induced while varying T sec , a/R, and b (8 seconds for the mid-eclipse time and 90 seconds for the eclipse duration).Therefore, performing the mapping considering a circular orbit while varying T sec , a/R, and b is analogous to effects that could be expected from an eccentric orbit.
Atmospheric Retrieval.We perform 1D atmospheric retrievals on the NAMELESS reduction at a resolution of 5 pixels per bin (408 bins) using four techniques with varying levels of physical assumptions: a self-consistent radiative-convective-thermochemical equilibrium grid retrieval (Sc-CHIMERA), a chemical equilibrium with free temperature-pressure profile retrieval (SCARLET), a free chemistry retrieval with thermal dissociation (HyDRA), and a free chemistry retrieval with abundances assumed constant with altitude (POSEIDON).All retrieval methods considered the same PHOENIX stellar spectrum 19 , produced using previously published parameters for WASP-18 (i.e., T ef f = 6435 K, log g = 4.35, and [Fe/H] = 0.1 20;21 ), to convert from model planet flux spectra to F p /F s values.We chose to use a model stellar spectrum instead of the extracted spectrum to avoid the possible introduction of systematic errors in the through the process of absolute flux calibration.
1D Radiative-Convective-Thermochemical Equilibrium (1D-RCTE) Grid Retrieval We use a 1D-RCTE grid retrieval-based method, ScCHIMERA 7;79 , with the opacity sources described in Mansfield et al. (2021) 2 and Glidic et al. (2022) 80 .These 1D-RCTE models assume cloud-free one-dimensional radiative-convective-thermochemical equilibrium using the methods described in Toon et al. (1989)  81 to solve for the net flux divergence across each layer of the atmosphere, and the Newton-Raphson iteration scheme of McKay et al. (1989)  82 to march towards an equilibrium vertical temperature structure.The NASA Chemical Equilibrium with Applications-2 (CEA2) routine 83 is used to compute the thermochemical equilibrium gas and condensate mole fractions for hundreds of relevant species.Opacities are computed with the correlated-K random-overlapresort-rebin method 84 .Input elemental abundances from Lodders & Palme (2009) 85 are scaled to a given metallicity ([M/H]) and carbon-to-oxygen ratio (C/O).
Using WASP-18b's planetary and stellar parameters, we produced a grid of 2730 1D-RCTE models and resulting top-of-atmosphere thermal emission spectra spanning the atmospheric C/O (0.01 -2.0) , [M/H] (-2.0 -2.0, where 0 is solar, 1 is 10×, etc), and heat redistribution (f , 1.0 -2.8, where 1=full, 2=dayside, and 2.67 is the max value, as defined in Arcangeli et al. (2019) 40 ).We additionally include a scale factor, A HS (allowed to vary from 0.5 -2.0), that multiplies the planetary flux by a constant to account for a hot spot area fraction emitting most of the observed flux 86 .The PyMultiNest 87 routine is used to sample the 1D-RCTE spectra via interpolation (and subsequent binning to the data wavelength bins) to obtain posterior probability constraints on the above parameters.We have made public our grid models, including temperature-pressure profiles, molecular abundances, and emission spectra, as well as additional figures showing the posteriors of retrieved parameters and the impact of each parameter on the spectrum.This can be found on Zenodo under the DOI: 10.5281/zenodo.7332105.

Chemical Equilibrium & Free Temperature Retrieval
We use the SCARLET atmospheric retrieval framework 18;88;89;90 to perform a chemical equilibrium retrieval with a free temperaturepressure profile on our retrieved dayside thermal emission spectrum.The SCARLET forward model computes the emergent disk-integrated thermal emission for a given set of molecular abundances, temperature structure, and cloud properties.The forward model is then coupled to the Affine Invariant Markov chain Monte Carlo Ensemble sampler emcee 71 to constrain the atmospheric properties.Due to the high dayside temperature and large pressures probed through thermal emission spectroscopy, we assume the atmosphere is in thermochemical equilibrium.For the equilibrium chemistry, we consider the following species: H 2 , H, H − 91;92 , He, Na, K 93;94 , Fe, H 2 O 95 , OH 96 , CO 96 , CO 2 96 , CH 4 97 , NH 3 98 , HCN 99 , TiO 100 , VO 101 , and FeH 102 .The abundances of these species are interpolated in temperature and pressure using a grid of chemical equilibrium abundances from FastChem2 103 , which includes the effects of thermal dissociation for all the species included in the model.These abundances also vary with the atmospheric metallicity, [M/H] (U[−3, 3]), and carbon-to-oxygen ratio, C/O (U[0, 1]), which are considered as free parameters in the retrieval.As for the temperature structure, we use a free parametrization 104 which here fits for N = 10 temperature points (U[100, 4500] K) with fixed spacing in log-pressure (P = 10 2 -10 −6 bar).Although this parametrization is free, it is regularized by a prior punishing for the second derivative of the profile using a physical hyper-parameter, σ s , with units of kelvin per pres-sure decade squared (K dex −2 ).This prior is implemented to prevent over-fitting and nonphysical temperature oscillations at short pressure scale lengths.For this work, we use a hyper-parameter value of σ s = 1000 K dex −2 , corresponding to a low punishment against second derivatives as we want the retrieval to explore freely the temperature-pressure profile parameter space.We note that further lowering this punishment does not affect the retrieved TP profile.Finally, we fit for an area fraction A HS (U[0, 1]) that is multiplied directly with the thermal emission spectrum, for a total of 14 free parameters.This factor is used to compensate for the presence of a hot spot which, while taking up only a portion of the planetary disk, contributes almost completely to the observed emission 86 .For the retrieval, we use four walkers per free parameter and consider the standard chi-square likelihood for the spectrum fit.We run the retrieval for 25,000 steps and discard the first 15,000 steps, 60% of the total amount, to ensure that the samples are taken after the walkers have converged.Spectra are initially computed using opacity sampling at a resolving power of R = 15,625, which is sufficient to simulate JWST observations 105 , convolved to the instrument resolution and subsequently binned to the retrieved wavelength bins.
Free Chemistry & Free Temperature-Pressure Profile Retrieval We use two independent atmospheric retrieval codes to perform free-chemistry retrievals on WASP-18b's dayside thermal emission spectrum: HyDRA 23;24;25;106 , including the effects of thermal dissociation, and POSEI-DON 107;108 , which here assumes constant-with-depth abundances for all chemical species.
HyDRA consists of a parametric forward atmospheric model coupled to a PYTHON implementation of the MULTINEST 109 Nested Sampling Bayesian parameter estimation algorithm 110 , PYMULTINEST 87 .The inputs to the parametric model are the deep-atmosphere mixing ratios for each of the chemical species included, the temperature-pressure profile parameters (6 free parameters), and a dilution parameter (area fraction) to account for 3D effects on the dayside 86 .Given the high dayside temperatures of WASP-18b, we consider high-temperature opacity sources and the effects of thermal dissociation, as described in Gandhi et al. (2020) 24 .We include opacity due to the molecular, atomic and ionic species with spectral features in the 0.8-2.8µm range which are expected in a high-temperature, H 2 -rich atmosphere: collision-induced absorption (CIA) due to H 2 -H 2 and H 2 -He 111 , H 2 O 96 , CO 96 , CO 2 96 , HCN 112 , OH 96 , TiO 100 , VO 101 , FeH 113 , Na 114 , K 114 and H − 91;92 .The line-by-line absorption cross sections for these species are calculated following the methods in Gandhi & Madhusudhan (2017) 115 , using data from the references listed.The H − free-free and bound-free opacity is calculated using the methods of Bell & Berrington (1987)  91 and John (1988) 92 , respectively.HyDRA includes the effects of thermal dissociation of H 2 O, TiO, VO and H − as a function of pressure and temperature, following the method of Parmentier et al. (2018) 9 .In particular, the abundance profiles of these species are calculated following Equations 1 and 2 of ref. 9 , where the deep abundance of each species (A 0 ) is a parameter in the retrieval, and the α, β and γ parameters are those given in Table 1 of that same work.The abundance profiles of the remaining chemical species are assumed to be constant with depth.
HyDRA uses the parametric temperature-pressure profile of Madhusudhan & Seager  (2009)  116 , which has been used extensively in exoplanet atmosphere retrievals, including ultrahot Jupiters such as WASP-18b 24 .The temperature parameterization is able to capture thermally inverted, non-inverted, and isothermal profiles, spanning the range of possible thermal structures for ultra-hot Jupiters.The HyDRA retrievals also include an area fraction parameter, A HS , which multiplies the emergent emission spectrum by a constant factor to account for the dominant contribution of the hot spot 86 .Given the input chemical abundances, temperature-pressure profile parameters and area fraction, the model thermal emission spectrum is calculated at a resolving power of R ∼15,000, convolved to the instrument resolution, binned to the data resolution, and compared to the data to calculate the likelihood of the model instance.Detection significances are calculated for specific chemical species by comparing the Bayesian Evidences of retrievals which include/exclude the species in question 89;117 .These detection significances factor in the ability of the retrieval to fit the observations with a different temperature profile and/or other chemical species, when the species in question is not included.Since thermal emission spectra are very sensitive to the atmospheric temperature profile, changes in the temperature structure can, in some cases, somewhat compensate for the absence of a particular chemical species, contributing to a lower detection significance.
We additionally use HyDRA to test the sensitivity of the free chemistry retrievals to the limits of the log-normal prior distributions assumed for the chemical abundances.For the HyDRA retrieval including thermal dissociation effects, we test two scenarios: wide, uninformative priors for all 11 species included (log mixing ratio ranging from -12 to -1), and slightly more restricted priors for the refractory species included (log mixing ratio ranging from -12 to -4 for TiO, VO and FeH, -12 to -2 for Na and K, and -12 to -1 for the remaining species).The more restricted prior limits for the refractory species are motivated by the relatively lower abundances expected for these species compared to the volatile species, across a range of metallicities and C/O ratios 118;119;120 .
We find that the atmospheric properties retrieved with HyDRA are consistent within 1σ for the two choices of prior limits.With the wide priors, the HyDRA retrieval infers a H 2 O log mixing ratio of -3.09 +1. 280.32 , with double-peaked posterior probability distributions for the H 2 O and TiO abundances.While the dominant posterior peaks correspond to approximately solar H 2 O and TiO abundances, the second, lower-likelihood peaks correspond to ∼ 30× and ∼ 10 4 × super-solar H 2 O and TiO abundances, respectively.Such an extreme TiO abundance warrants skepticism, and may, for example, be a result of the well-known degeneracy between chemical abundances and the atmospheric temperature gradient (see also Evans et al. (2017)  121 ).The retrieved H 2 O abundance in this case is consistent with the chemical equilibrium retrievals and self-consistent 1D radiativeconvective models described above, though with a larger uncertainty due to the double-peaked posterior distribution.When the restricted priors are used, the low-likelihood, high-abundance peaks in the H 2 O and TiO posterior distributions are no longer present, and the retrieved H 2 O abundance is -3.23 +0. 450.29 , in excellent agreement with the chemical equilibrium retrievals and self-consistent 1D radiative-convective models.We note that the retrieved temperature-pressure profiles and abundances of CO, CO 2 , HCN, OH, H − and FeH are unaffected by the choice of prior limits discussed above.The abundances of Na and K are unconstrained in both cases.While the two choices of prior limits give consistent results, the expectation of chemical equilibrium in the atmosphere of WASP-18b, in addition to the unlikelihood of a ∼ 10 4 × super-solar TiO abundance, motivate the use of the somewhat restricted priors on the refractory chemical abundances.
We note that for either choice of prior, the retrieved deep-atmosphere abundance of VO is significantly higher than that inferred by the chemical equilibrium retrieval (Extended Data Fig. 6).This is due to a difference in the thermal dissociation prescriptions; in the HyDRA retrieval, thermal dissociation results in a significantly depleted VO abundance at the photosphere, while in the chemical equilibrium retrieval thermal dissociation begins at lower pressures.Furthermore, the posterior distribution for the VO abundance peaks at highly super-solar values (∼ 10 4 × solar).Such a high abundance is physically unlikely, and may indicate the presence of further sources of optical opacity not included in the retrieval.
We also conduct free-chemistry retrievals, without the inclusion of thermal dissociation, using POSEIDON.POSEIDON is an atmospheric retrieval code originally designed for the interpretation of exoplanet transmission spectra 107 .We have recently extended POSEIDON to include secondary eclipse emission spectra modelling and retrieval capabilities.For an ultra-hot Jupiter such as WASP-18b, where the dayside can be assumed clear, POSEIDON's emission forward model calculates the emergent flux via a standard single stream prescription without scattering where I layer bot and I layer top are, respectively, the upwards specific intensity incident on the lower layer boundary and the intensity leaving the upper layer boundary, dτ layer is the differential vertical optical depth across the layer, µ = cos θ specifies the ray direction, and B(T layer , λ) is the black body spectral radiance at the layer temperature.Using the boundary condition I deep (µ, λ) = B(T layer , λ), POSEIDON propagates Equation 8 upwards to determine the emergent intensity at the top of the atmosphere.The emergent planetary flux is determined via where W are the Gaussian quadrature weights corresponding to each µ (here taken as second order quadrature over the interval µ ∈ [0, 1]).Finally, the planet-star flux ratio seen at Earth is given by where R p, phot is the effective photosphere radius 122;123 at wavelength λ (evaluated at τ (λ) = 2/3).
Since the calculation of R p, phot requires a reference radius boundary condition to solve the equation of hydrostatic equilibrium, we prescribe r(P = 10 mbar) as a free parameter.
For the WASP-18b POSEIDON retrieval analysis, we calculated emission spectra via opacity sampling and explored the parameter space using MULTINEST via its Python wrapper PY-MULTINEST 87;109 .POSEIDON solves the radiative transfer on an intermediate resolution wavelength grid (here, R = 20,000 from 0.8 to 3.0 µm), onto which high-spectral-resolution (R ∼ 10 6 ), pre-computed cross sections 124 are down-sampled.For WASP-18b, we consider the following opacity sources: H 2 -H 2 125 and H 2 -He 125 CIA, H 2 O 95 , CO 126 , CO 2 127 , HCN 128 , H − 92 , OH 129 , FeH 102 , TiO 100 , VO 101 , Na 130 , and K 130 .We prescribed uniform-in-altitude mixing ratios, defined by a single free parameter for each of the chemical species included.The PYMULTINEST retrievals with POSEIDON use 2,000 live points, and the six-parameter temperature-pressure profile 116 outlined above in the description of HyDRA.POSEIDON accounts for the dominant contribution of the hot spot by prescribing the 10 millibar radius as a free parameter, which is subsequently converted into an equivalent A HS posterior by comparison to the white light planet radius.
We note that both HyDRA and POSEIDON yield consistent retrieval results when thermal dissociation is not considered in the HyDRA retrievals.However, the inclusion of thermal dissociation results in significantly different retrieved H 2 O abundances and temperature-pressure profiles, which are in agreement with the chemical equilibrium retrievals and self-consistent 1D radiative-convective models described above.
Disequilibrium Chemistry Model.To further justify the use of chemical equilibrium models in our analysis of WASP-18b, we produce a grid of disequilibrium chemistry forward models to assess whether disequilibrium effects might strongly shape our observations.We begin by calculating the atmospheric temperature-pressure structure of WASP-18b under radiative-convective equilibrium with the HELIOS 131;132 radiative transfer code.Next, we calculate altitude-dependent mixing ratios of chemical species under this temperature-pressure structure with the VULCAN 133 1-dimensional chemical kinetics code, using an N-C-H-O reaction network that includes ionization and recombination of Fe, Mg, Ca, Na, K, H, and He.We use the current version of VULCAN (VULCAN2 134 ), which includes optional photochemistry and parameterizes the transport flux of chemical species with eddy diffusion, molecular diffusion, thermal diffusion, and vertical advection.We update this code to include the effect of photoionization.Finally, we generate emission spectra with the PLATON radiative transfer code 135 at the resolution and wavelength range of NIRISS/SOSS. 6We modify PLATON to calculate bound-free and free-free H − opacity as a function of altitude; this addition is necessary to assess whether disequilibrium abundance H − opacity could mute spectral features more strongly than predicted by equilibrium chemistry.We furthermore modify PLATON to accept higher-temperature (T > 3000 K) opacity files that we calculate with the HELIOS-K code 136;137 .
For our set of models, we vary the eddy diffusion coefficient, K zz , from 10 7 cm −2 s −1 to 10 13 cm −2 s −1 , holding it constant at all altitudes for a given simulation.We perform this sweep over many orders of magnitude of K zz to understand the maximum effect that disequilibrium chemistry could have on the observed emission spectrum.While K zz is a limited descriptor of vertical mixing and is expected to vary as a function of altitude (e.g., Parmentier et al. (2013)  138 ), we assume that our forward models bracket the expected vertical mixing behavior of this planet.This statement is further motivated by our GCM models, if we approximate K zz = vH for vertical wind velocity v and atmospheric scale height H 139;140 .The minimum dayside-average K zz for our kinematic MHD GCM (see General Circulation Models section) is ∼ 10 8 cm −2 s −1 , and the maximum dayside-average K zz is ∼ 10 9 cm −2 s −1 , well within our VULCAN grid range.Our model grid also toggles the inclusion of molecular diffusion and photochemistry.As input to VULCAN when photochemistry is included, we use a stellar spectrum that is appropriate for WASP-18 from Rugheimer et al. (2013) 141 .The spectrum is constructed by joining synthetic spectra from Kurucz (1979)  142 and UV flux measurements of Piscium HD 222368 by the International Ultraviolet Explorer7 at 300 nm.
We find that our inclusion of disequilibrium chemistry effects-photochemistry, molecular diffusion, thermal diffusion, and eddy diffusion-produces spectra that are not strongly discrepant from spectra computed assuming chemical equilibrium.Indeed, all discrepancies between spectra produced under chemical equilibrium and spectra produced under chemical disequilibrium spectra are less than 10 ppm.This agreement is expected, as chemistry at the pressure levels probed by low-resolution emission spectra are not predicted to be strongly modified by photochemistry or mixing (e.g., Tsai et al. (2021) 134 ).Furthermore, the high temperature of WASP-18b implies that the chemical timescales in the atmosphere are quite short, allowing chemical reactions to occur more quickly than the relevant disequilibrium timescales.
In sum, our grid of chemical disequilibrium forward models indicates that disequilibrium chemistry effects considered here do not strongly affect the emission spectrum of WASP-18b in the NIRISS/SOSS waveband.
General Circulation Models.A suite of General Circulation Models (GCMs) are compared to the retrieved dayside spectrum and dayside temperature map.
We use the SPARC/MITgcm 143 to model the 3D atmospheric structure of WASP-18b.The model solves the primitive equations in spherical geometry using the MITgcm 144 and the radiative transfer equations using a current 1D radiative transfer model 145 .We use the correlated-k framework to generate opacities based on the line-by-line opacities 146 .Our model assumes a solar composition for the elemental abundances 147 and chemical equilibrium gas-phase composition 148 .Our model naturally takes into account the effect of thermal dissociation 9 .We used a time step of 25 s and ran the simulations for about 300 Earth days, averaging all quantities over the last 100 days.
We include additional sources of drag through a Rayleigh drag parametrization with a single constant timescale per model that determines the efficiency with which the flow is damped.The drag is constant over the whole planetary atmosphere.We vary this timescale between models from τ drag = 10 3 to 10 6 s (efficient drag) and a no drag model with τ drag = ∞.Our range of drag strengths cover the transition from a drag-free, wind-circulation case to a drag-dominated circulation.The specific WASP-18b simulations that we use are described in more detail in Arcangeli et  al. (2019)  40 .
The second model we use is the kinematic magnetohydrodynamics (MHD) GCM (described in detail in ref. 46 ) with an updated picket fence radiative transfer scheme 149;150 .Due to the high gravity of this planet, we chose to model the planet from 100 bar to 10 −4 bar over 65 layers at a horizontal resolution of T31 (corresponding to roughly 3 degree resolution at the equator).We use the kinematic MHD prescription described in Perna et al. (2010)  43 , which has been used in models of hot Jupiters HD 209458b and HD 189733b 151 as well as the ultra-hot Jupiter WASP-76b 46;152 .This drag prescription assumes a global dipole magnetic field, generated by an interior dynamo.Because of this geometry, our drag timescale is applied as a Rayleigh drag term solely to the eastwest momentum equation (influencing flow perpendicular to magnetic field lines) and is calculated as where B is the chosen global magnetic field strength, φ is the latitude, ρ is the density, and η is magnetic resistivity.This timescale is calculated locally and often, allowing the timescale to vary by over 10 orders of magnitude throughout the model and respond to changes in atmospheric temperatures.A minimum timescale cutoff (∼ 10 3 s) is applied in locations where τ mag would be less than 1/20 of the the planet's orbit, for numerical stability.We chose to model a range of magnetic field strengths (0 G, 5 G, 10 G, and 20 G) as its true value is not known.Past HST 7 (red points), TESS (see Methods), and Spitzer 153 (gray points) are shown for comparison.We show the best-fit model (blue line) from the SCARLET chemical equilibrium retrieval , extrapolated to the TESS and Spitzer wavelengths considering the same atmospheric parameters.We find that the retrieved spectrum is in good agreement with the past HST observations.The throughput-integrated model is shown for the TESS and Spitzer points (blue points).The white (broadband) light curve (white points) and three example spectrophotometric light curves (blue, green, and orange points at 1.05, 1.72, and 2.77 µm respectively), along with their best fitting models (black line), are shown to scale.The phase variation of the measured planetary flux around the secondary eclipse is clearly visible.b, Planetary thermal emission spectrum of WASP-18b, as computed from the F p /F s spectrum and the PHOENIX stellar spectrum.The shortest wavelengths of the NIRISS/SOSS first order reach the maximum of the planetary spectral energy distribution, thereby enclosing 65% of the total thermal energy emitted by the planet.Blackbody spectra for temperatures T = 2850 (dotted), 2950 (dash-dotted), and 3050 K (dashed) are shown in purple, with the best-fitting blackbody spectrum to the NIRISS data being T = 2950±3 K. | Atmospheric constraints from the chemical equilibrium and free chemistry retrievals.a, Retrieved temperature-pressure profiles with 1 and 2 σ contours for the chemical equilibrium with free temperature-pressure profile (blue), radiative-convective thermochemical equilibrium (1D-RCTE, red), and free chemistry with thermal dissociation (green) retrievals.The retrieved temperature-pressure profiles are consistent between the retrievals and show an inversion in the pressure range that is constrained from the observations, as shown by the contribution functions at 0.85 (dot-dashed gray line), 1.82 (dashed brown line), and 2.83 µm (orange line).The temperature-pressure profile of WASP-18b is above the CaTiO 3 condensation curve 154 (black dashed line) at almost all pressures, which motivates the presence of a temperature inversion caused by TiO as Ti is available in gas form.The dayside average temperature-pressure profile of the τ drag = 10 3 s SPARC/MITgcm (white dashed line) is computed from the viewing angle average of T (P ) 4 and shown for comparison.We also show the posterior probability distributions of the atmospheric metallicity [M/H] (b), C/O ratio (c), and area fraction A HS (d).The area fraction A HS is a scaling factor applied to the thermal emission spectrum to compensate for the possible presence of a concentrated hot spot contributing to most of the observed emission 86 .All methods retrieve metallicities consistent with solar at 1σ.The retrieved C/O 3σ upper limits are of 0.6 and 0.2 for the chemical equilibrium with free temperature-pressure profile and the 1D-RCTE retrievals respectively.Finally, we find the area fraction A HS is consistent with 1 when allowing the temperature-pressure profile to vary freely, indicating the lack of a concentrated hot spot on the dayside contributing to the majority of the observed emission.Extended Data Fig. 6 | Abundance constraints from the free chemistry retrievals.Probability posteriors of the deep abundance of various species considered for the free chemistry and temperature with (blue, HyDRA) and without (red, POSEIDON) thermal dissociation.We also show the median retrieved VMR profiles from the chemical equilibrium with free temperature-pressure profile retrieval (black line, SCARLET).The pressure range probed by the observations (∼0.01-1 bars, see 3) is indicated by the dashed gray lines.The only species independently detected is H 2 O, which is found to be consistent with the retrieved chemical equilibrium abundance when considering the effect of thermal dissociation.All other species considered are found to be unconstrained although consistent with chemical equilibrium predictions.The photosphere as predicted by our radiative-convective model is around 50 millibar, but the retrievals infer the deep molecular abundances.presented in the text.The red model uses a constant-with-latitude Fourier series as the basis set, rather than spherical harmonics, to investigate constraints on latitudinal aspects of the map.Shaded regions denote 1, 2, and 3 σ quantiles.b, Same as a but for the eclipse egress.c, Planetary flux along the equator for the same two models.Note that, regardless of basis functions, we retrieve the same longitudinal structure, giving us confidence in the longitudinal brightness distribution.d, Same as c but along the substellar meridian.Both models fit the data well but find different latitudinal structure, indicating that we are unable to constrain latitudinal variation.

Figure 1 |
Figure1| Dayside thermal emission spectrum of WASP-18b.a, Observed dayside planet-to-star flux ratio spectrum (black points), binned at a fixed resolving power of R = 50 for visual clarity.Past HST 7 (red points), TESS (see Methods), and Spitzer 153 (gray points) are shown for comparison.We show the best-fit model (blue line) from the SCARLET chemical equilibrium retrieval , extrapolated to the TESS and Spitzer wavelengths considering the same atmospheric parameters.We find that the retrieved spectrum is in good agreement with the past HST observations.The throughput-integrated model is shown for the TESS and Spitzer points (blue points).The white (broadband) light curve (white points) and three example spectrophotometric light curves (blue, green, and orange points at 1.05, 1.72, and 2.77 µm respectively), along with their best fitting models (black line), are shown to scale.The phase variation of the measured planetary flux around the secondary eclipse is clearly visible.b, Planetary thermal emission spectrum of WASP-18b, as computed from the F p /F s spectrum and the PHOENIX stellar spectrum.The shortest wavelengths of the NIRISS/SOSS first order reach the maximum of the planetary spectral energy distribution, thereby enclosing 65% of the total thermal energy emitted by the planet.Blackbody spectra for temperatures T = 2850 (dotted), 2950 (dash-dotted), and 3050 K (dashed) are shown in purple, with the best-fitting blackbody spectrum to the NIRISS data being T = 2950±3 K.

Figure 3
Figure3| Atmospheric constraints from the chemical equilibrium and free chemistry retrievals.a, Retrieved temperature-pressure profiles with 1 and 2 σ contours for the chemical equilibrium with free temperature-pressure profile (blue), radiative-convective thermochemical equilibrium (1D-RCTE, red), and free chemistry with thermal dissociation (green) retrievals.The retrieved temperature-pressure profiles are consistent between the retrievals and show an inversion in the pressure range that is constrained from the observations, as shown by the contribution functions at 0.85 (dot-dashed gray line), 1.82 (dashed brown line), and 2.83 µm (orange line).The temperature-pressure profile of WASP-18b is above the CaTiO 3 condensation curve154 (black dashed line) at almost all pressures, which motivates the presence of a temperature inversion caused by TiO as Ti is available in gas form.The dayside average temperature-pressure profile of the τ drag = 10 3 s SPARC/MITgcm (white dashed line) is computed from the viewing angle average of T (P )4 and shown for comparison.We also show the posterior probability distributions of the atmospheric metallicity [M/H] (b), C/O ratio (c), and area fraction A HS (d).The area fraction A HS is a scaling factor applied to the thermal emission spectrum to compensate for the possible presence of a concentrated hot spot contributing to most of the observed emission86 .All methods retrieve metallicities consistent with solar at 1σ.The retrieved C/O 3σ upper limits are of 0.6 and 0.2 for the chemical equilibrium with free temperature-pressure profile and the 1D-RCTE retrievals respectively.Finally, we find the area fraction A HS is consistent with 1 when allowing the temperature-pressure profile to vary freely, indicating the lack of a concentrated hot spot on the dayside contributing to the majority of the observed emission.

lmax = 5 ,Figure 4 |. 1 |. 2 |. 3 |
Figure 4 | Retrieved temperature map of WASP-18b.a, Latitudinally-averaged brightness temperatures (see Methods) of the planet along the equator.The blue and red shaded areas show solutions for l max = 5, N = 5 and l max = 2, N = 5 (see Methods), respectively.Statistically, the blue model is marginally preferred.Dark, medium, and light shading denote the 1, 2, and 3 σ confidence regions respectively, showing the range of model possibilities.Overplotted are several predictions from general circulation models with magnetic-field 152 (green) or uniform 143 (purple) drag timescales (see Methods).The plot only shows longitudes emitting at least 10% of the substellar flux.b, The temperature map of WASP-18b for the l max = 5, N = 5 solution.Along the equator at -90 • , 0 • , and 90 • longitude the temperatures are 1744 K, 3121 K, and 2009 K, respectively.c, The colorbar for the map shown in b.Color represents the brightness temperature and saturation represents the relative contribution to the light curve based on its visibility.

. 5 |
Light curve residuals binned in time.a, Absolute root mean square (RMS) of the residuals as a function of bin size (black line) for the white light curve.The RMS values are plotted against the Poisson noise limit (red dashed line), which goes down as the square root of the number of integrations contained in a single bin.The residuals bin down to ∼5 ppm for bins of 1 hour.The broadband residuals do not follow perfectly the Poisson noise, which is indicative of remaining time correlations.b, Normalized RMS of the 408 spectrophotometric light curves considered in the analysis.We observe that the residuals follow the Poisson noise limit from bin sizes of a single integration up to bins of ∼75 minutes, indicating that there are no time correlations in the residuals.

. 7 | 2 . 10 |
WASP-18b's brightness temperature spectrum fit and source of the thermal inversion.a, The dark blue line indicates the chemical equilibrium median fit to the NIRISS data (black points), with shaded blue regions showing the 1σ and 2σ credible intervals in the retrieved spectrum (medium and light blue, respectively).The spectra are extrapolated to the TESS (visible wavelengths) and the Spitzer/IRAC measurements (3.6 and 4.5 µm) observations (gray points) considering the same atmospheric parameters.b, Best-fit radiative-convective model temperature-pressure profile together with radiative-convective solutions where specific species known to create a thermal inversion are removed from the atmosphere.Absorption by atomic iron contributes to the thermal inversion at pressures lower than 1.0 millibar, whereas TiO is responsible for the thermal inversion seen between 0.1 and 0.01 bar.SiO contributes at pressures lower than 0.1 bar.c, Best-fit radiative-convective brightness temperature spectrum (excluding area fraction) and resulting spectra when removing specific species.As shown by the change from emission to absorption features in the spectra when TiO is removed, the TiO-induced thermal inversion is the one probed by our observations./N = 1.96 χ 2 /N = 1.56 χ 2 /N = 1.11 χ 2 /N = 1.48 b Extended Data Fig. 8 | Light-curve predictions from GCMs.a, GCM light curves compared against the data and the best-fitting eclipse-mapping model.b, Data and the GCM light curves with the eclipse-mapping model subtracted.The GCMs with strong atmospheric drag (RM-GCM 20 G and SPARC/MITgcm τ = 10 3 s) match the data better than their counterparts that have little to no drag.Latitudinal structure in the eclipse map.a, Ingress of the eclipse, with two models overplotted and a 1096 ppm (white-light planet flux at mid-eclipse) uniform planet model subtracted to highlight deviations.The data (small dots) have been binned by a factor of five (dots with error bars) for clarity.The blue model is the eclipse map for l max = 5, N = 5 nirHiss Reduction We use the nirHiss Python open-source data reduction pipeline as described in Feinstein et al. (2023) 50 .To summarize, this pipeline uses Eureka! 51to go from the Stage 0 JWST outputs to Stage 2 calibration, which applies detector-level corrections, produce count rate images, and calibrates individual exposures.From the Stage 2 "calints" FITS files, we use nirHiss to correct for background sources, 1/f noise, cosmic ray removal, trace identification, and spectral extraction.