Detection of Solar-like Oscillations and Granulation at Infrared Wavelengths using JWST

Introduction

The surfaces of cool stars exhibit low-level variability due to processes arising from convection [see 44 for a review]. Granulation—the rising and falling plasma in convective cells—produces brightness fluctuations on timescales from minutes to hours. The granulation timescale, $\tau_{\mathrm{gran}}$, is related to the size of the granules, and has been shown to depend on stellar parameters [86], namely surface gravity [60, 5], $T_{\mathrm{eff}}$ [52, 50, 67], and metallicity [60, 69].

As well as producing granulation, near-surface convection also drives the stochastic excitation and damping of pressure-mode oscillations. These oscillations can be summarized by two global properties: the frequency maximum power, $\nu_{\max}$ [12], and the average frequency separation between consecutive oscillation modes of the same angular degree, $\Delta\nu$ [83]. These scale with stellar properties such as $\log g$, $T_{\mathrm{eff}}$ and [Fe/H] [51, 52], and empirical scaling relations thus provide precise constraints on stellar density, radius, and age [51, 8, 6]. This has been demonstrated for hundreds of thousands of stars using space-based data from Kepler and TESS (The Transiting Exoplanet Survey Satellite, [68]) (see [13, 30, 31, 48]).

Email: [email protected]

Amplitudes of stellar granulation and oscillations directly probe the physics of convective energy transport, which is poorly understood [41, 42]. Traditional 1D models using Mixing Length Theory (MLT, [10]) and non-local MLT (e.g. W. Unno 1967; D. O. Gough 1977a,b; M.-A. Dupret et al. 2006) have been used to derive amplitude scaling relations based on energy equipartition (e.g. P. Goldreich & D. A. Keeley 1977a,b; J. Christensen-Dalsgaard & S. Frandsen 1983a,b), which show broad agreement with observations [51, 45, 50]. Detailed 3D convection models of individual stars can predict amplitudes from first principles (e.g. [58, 82, 72, 70]), but are computationally expensive. Scaling relations derived from the results of 3D models offer the best test for the largest samples of stars, but still often require ad-hoc correction factors (e.g. [21]) to match observations.

A powerful test of convection models would be the measurement of granulation and oscillation amplitudes at different wavelengths. Both scaling relations and 3D models predict that amplitudes of granulation and oscillations decrease at longer wavelengths [51, 4, 55, 72], consistent with the assumption that convection can be approximated as small temperature fluctuations [61, 4]. However, this assumption has only been tested observationally over the small wavelength difference between Kepler and TESS bandpasses for stars other than the sun (e.g. [80]), or between the three visible filters from the SPM (Sun Photometers) on the VIRGO (Variability of solar IRradiance and Gravity Oscillations, C. Fröhlich et al. [29] mission, for the Sun (E. Michel et al. [61]; J. Ballot et al. [4]; M. N. Lund et al. [56]).

The James Webb Space Telescope (JWST) offers the opportunity to study stellar variability due to convection at infrared wavelengths for the first time. JWST regularly obtains continuous light curves of stars to study exoplanet atmospheres that have excellent long-term stability and photometric precision (see e.g. J. L. Bean et al. [7]; L. Alderson et al. [3]; Z. Rustamkulov et al. [71]; A. D. Feinstein et al. [26]; L.-P. Coulombe et al. [19]; T. J. Bell et al. [9]). Granulation is expected to be a significant source of variability for some of these JWST observations (B. Morris & D. Huber, in prep), and has been reported as a likely source of correlated non-instrumental signal for some solar-type stars (L.-P. Coulombe et al. [20]; J. Splinter et al. [79]; J. K. Teske et al. [84]; R. Hu et al. 2024; S. Boucher et al. [11]). However, no confirmed detection of oscillations and granulation has been made to date.

Here, we report the detection of solar-like oscillations and granulation in WASP-121 using 36.9 hours of NIRISS/SOSS (R. Doyon et al. [23]; L. Albert et al. [2]) observations obtained as part of the NIRISS Exploration of the Atmospheric diversity of Transiting exoplanets (NEAT) Guaranteed Time Observations Programme (GTO 1201; PI Lafrenière). WASP-121 is a late F-type dwarf star, hosting an ultra-hot Jupiter on a 1.27 day orbit (L. Delrez et al. [22]). Scaling relations (T. M. Brown et al. [12]) predict a frequency of maximum power of $\nu_{\max} \simeq1730\ \mu\mathrm{Hz}$, and smoothing the power spectrum, which is generated using the NIRISS/SOSS Order 1 integrated white light curve (centered at $\sim1.5\ \mu\mathrm{m}$), reveals an oscillation excess that is consistent with this value. Comparison with a Kepler star with high SNR and similar stellar properties demonstrates a close match, further strengthening the detection. We also detected at least five evenly spaced peaks, which is a clear signature of solar-like oscillations. This is the first detection of solar-like oscillations at infrared wavelengths in a star, including our own Sun.

Asteroseismology

We used the residuals of the NIRISS/SOSS Order 1 integrated lightcurve to search for signatures of stellar variability. We first converted the time series to a Lomb-Scargle periodogram, and then converted power into power spectral density (PSD). We used the asteroseismic scaling relations H. Kjeldsen & T. R. Bedding [51]; T. M. Brown et al. [12] as follows, to predict the global seismic parameters, $\nu_{\max}$ and $\Delta\nu$:

Here $M_\star$, $R_\star$ and $T_{\mathrm{eff}}$ are the stellar mass, radius and effective temperature adopted from D. K. Sing et al. [77], who used high precision spectroscopy, with values of $M_\star= 1.35 \pm0.09\,M_\odot$, $R_\star= 1.46 \pm0.03\,R_\odot$, and $T_{\mathrm{eff}} = 6628 \pm140\mathrm{K}$. This yielded a predicted $\Delta\nu$ and $\nu_{\max}$ of $\sim89\ \mu\mathrm{Hz}$ and $\sim1830\ \mu\mathrm{Hz}$, respectively.

Comparison to a known oscillator with similar physical parameters is a useful way to verify low-SNR detections of oscillations. We used a similar Kepler star, KIC 9139163 (HD 176071), which was part of the Kepler Legacy Sample (M. N. Lund et al. [57]; V. Silva Aguirre et al. [76]) – one of the most precisely characterised samples of asteroseismic targets observed to-date. KIC 9139163 is similar to WASP-121 in $T_{\mathrm{eff}}$, mass, radius and [Fe/H] ($M_\star= 1.42 \pm0.06\ M_\odot$, $R_\star= 1.58 \pm0.02\ R_\odot$, [Fe/H] $= 0.15 \pm0.09$, and $T_{\mathrm{eff}} = 6400 \pm84\ \mathrm{K}$), and has a measured $\nu_{\max}$ of $1730\ \mu\mathrm{Hz}$ – close to our predicted value for WASP-121.

The PSD of WASP-121 is shown in Figure 1. The solid black curve shows these data with a Gaussian filter of width $4\Delta\nu$, which confirms an oscillation envelope in the region of the predicted $\nu_{\max}$, shown via the orange band. Applying the same smoothing to the KIC 9139163 data (the solid blue curve) confirms that the $\nu_{\max}$ is similar for both stars. The WASP-121 observations spanned only 36.9 hours (approximately 1.5 days), whereas the KIC 9139163 data span a period of over 4 consecutive years, which explains the relatively low signal-to-noise ratio (SNR) of the WASP-121 detection.

The inset of Figure 1 shows the WASP-121 power spectrum in linear frequency space, cropped around the predicted $\nu_{\max}$. Five peaks are clearly visible around the between 1500 and 1900$\mu$Hz, which are evenly spaced by a frequency value which is similar to the predicted $\Delta\nu$. The PSD average over 15 1.5-day segments of the KIC 9139163 data is shown on the same plot in blue, for comparison.

To measure the $\nu_{\max}$ and $\Delta\nu$, we used pySYD (A. Chontos et al. [15]) to fit a Harvey-like background component and white noise component to remove the lower-frequency granulation signals. With this method, we measured a $\nu_{\max}$ of 1615 $\pm$ 67 $\mu$Hz with an amplitude of 5.21 $\pm$ 2.16 ppm. The right hand panel of Figure 2 compares these measurements to stars observed by Kepler (M. Sayeed et al. [75]), which show the well-known relation that stars with lower $\nu_{\max}$ oscillate with larger amplitudes. The amplitude of WASP-121 is slightly lower than the average amplitude of Kepler stars with similar $\nu_{\max}$ and $T_{\mathrm{eff}}$, such as with KIC 9139163 (A = 6.19 ppm $\pm$ 1.46 ppm), in agreement with expectations for infrared measurements.

We also compared the échelle diagrams of the two stars. The WASP-121 échelle shows only a few peaks, however the locations of the visible vertical ridges compare well to the échelle of the Kepler star, as seen in the right-hand panel of Figure 2. Mode identification is not possible with the JWST échelle alone, because hot F-type dwarfs typically exhibit blending between the $l=0$ and $l=2$ modes in frequency space, making mode identification challenging (see e.g. T. R. White et al. [91]). However, by comparing to the well-resolved Kepler star, whose dipole and quadrupole modes are detectable individually, due to the 4-year Kepler dataset, we identified the $l=0/2$ and $l=1$ mode ridges on the left and right hand side of the WASP-121 échelle respectively.

Granulation

Granulation, the surface signature of convection, introduces a broad low-frequency background in observations. We detect granulation in WASP-121, which appears as a slope with negative gradient at frequencies lower than the solar-like oscillations. Many studies use methods similar to pySYD to account for this ‘background’, fitting one or more Harvey-like profiles (J. Harvey [40]) plus a white-noise term to the power spectrum (D. Huber et al. [45]). et al. 2011; A. Chontos et al. [16]; J. R. Larsen et al. [53]). However, frequency-domain fits can be unreliable for short or noisy datasets. We therefore used a Gaussian Process (GP) to model oscillations and granulation simultaneously in the time domain. We found granulation amplitudes and frequencies consistent with predictions from literature (e.g. T. Kallinger et al. [50]).

Similarly to L.-P. Coulombe et al. [20] and J. Splinter et al. [79], and following D. Foreman-Mackey et al. [27], we described the variability as the sum of two damped simple harmonic oscillator (SHOTerm) components implemented in celerite2: one for the granulation background and the other the p-mode oscillation envelope. The power spectral density of each term is

where $S_0$ is the amplitude, $\omega_0$ the characteristic frequency, and $Q$ the quality factor. For the background component, we fix $Q = \frac{1}{\sqrt{2}}$ (D. Foreman-Mackey et al. [27]; F. Pereira et al. [64]). For the oscillations, $Q$ was allowed to vary. The full time-domain covariance of the kernels is

We note that the lower-frequency supergranulation is challenging to measure due to the limited length of the dataset. We therefore recommend that with JWST data of similar length, only the median-timescale granulation is likely to be measurable with this method, especially with simultaneous detection of solar-like oscillations. A future paper will investigate the feasibility of measuring the granulation in datasets where oscillations are not detectable. A decline in the smoothed power spectrum of WASP-121 is indeed visible at several hundred $\mu\mathrm{Hz}$ (Figure 1), consistent with the expected granulation frequencies (T. Kallinger et al. [50]). To robustly predict the frequency of granulation, we used the empirical scaling relation of T. Kallinger et al. [50],

together with our measured $\nu_{\mathrm{max}}$ of 1615 $\mu\mathrm{Hz}$, yielding $\nu_{\mathrm{gran}} \approx439$ $\mu\mathrm{Hz}$.

Using this estimate as a prior, we calculated the maximum a posteriori estimate of the GP model, and used the outputs as starting values to sample the likelihoods with PyMC3. For the white noise component, the mean $\mu_w$ was assigned a normal prior $\mathcal{N}(0, 10)$, while its variance $S_w^2$ followed a log-normal prior $\ln(\mathcal{N}(y_{\mathrm{err}}, 10))$, centered on the mean reported measurement uncertainty. For the granulation component, the amplitude $A_{\mathrm{gran}}$ was given a log-normal prior $\ln(\mathcal{N}(\sigma_y, 50))$ based on the standard deviation of the light curve. We note that the generous prior width was chosen here such that the photon noise should not bias the measurement of the amplitude, which is further ensured by the simultaneous fitting of the white-noise term. The characteristic frequency $\nu_{\mathrm{gran}}$ has a normal prior $\mathcal{N}(439, 50)$ $\mu\mathrm{Hz}$, and the quality factor $Q_{\mathrm{gran}}$ was fixed at $1/\sqrt{2}$. For the oscillation component, the amplitude $A_{\mathrm{osc}}$ followed a log-normal prior $\ln(\mathcal{N}(\sigma_y^2, 1))$, the frequency of maximum power $\nu_{\mathrm{max}}$ a normal prior $\mathcal{N}(1615, 100)$ $\mu\mathrm{Hz}$, and the quality factor $Q_{\mathrm{osc}}$ a log-normal prior $\ln(\mathcal{N}(3, 1))$.

We applied the same GP approach to 30, 30-day sections of Kepler data for KIC 9139163, taking the mean and 1-sigma values of each GP parameter. Although granulation in Kepler stars is commonly modelled with a two-component mixture (T. Kallinger et al. [50]), doing so here would not provide a like-for-like comparison with the JWST analysis. The JWST time series has a model of a $\sim1.13$ day rotation period removed, using a semi-sinusoidal stellar activity model. We therefore applied a 2-term Savitzky–Golay filter (A. Savitzky & M. J. E. Golay [74]) to the Kepler data to remove variability on comparable timescales. After this pre-processing, we expected that the remaining Kepler variability is dominated by granulation and the oscillation envelope.

Similar granulation properties are determined for both stars. The GP in both the time and frequency domains is shown in Figure 3, for both WASP-121 (right), and for one of the 30-day data sections of KIC 9139163 (left). In general, the amplitudes of the variability for both stars are similar, and the two-component GP model describes the convective processes well in both stars. For WASP-121, we measured $\nu_{\mathrm{max}} = 1619^{+72.71}_{-78.22}$ $\mu\mathrm{Hz}$, $\Delta\nu= 82.147 \pm0.09$ $\mu\mathrm{Hz}$, granulation timescale $\tau_{\mathrm{gran}} = 9.79^{+1.26}_{-0.22}$ min, and the granulation amplitude $A_{\mathrm{gran}} = 65.96 \pm4.88$ ppm. For KIC 9139163, the corresponding values were $\nu_{\mathrm{max}} = 1686^{+68.82}_{-74.36}$ $\mu\mathrm{Hz}$, $\Delta\nu= 81.2 \pm0.09$ $\mu\mathrm{Hz}$, $\tau_{\mathrm{gran}} = 8.33^{+0.31}_{-0.28}$ min, and $A_{\mathrm{gran}} = 60.5 \pm2.28$ ppm.

To validate the astrophysical nature of the oscillation and granulation signals, we considered the relationship between their amplitudes. This comparison was made without the correction from F. Pereira et al. [64], such that the ratio

was described by

For the Order 1 NIRISS/SOSS data, this yielded a value of R=2.78. Using the equation and coefficients from tbl:2 of T. Kallinger et al. [50] gives a predicted R=3.10, which agrees well with the measured value.

Testing the Chromaticity of Granulation

The different observed amplitudes of granulation and oscillations in KIC 9139163 and WASP-121 may be attributed to the slight difference in stellar parameters, instrumental effects, or brightness. Therefore, these measurements alone are not sufficient to test for the wavelength dependence of the amplitudes. Such a test requires simultaneous data of the same star, made at different wavelengths.

To test whether the amplitudes of the signals observed vary with wavelength, we performed an identical analysis on the Order 2 data. Although Order 2 NIRISS/SOSS data is known to be of lower SNR ([79]), oscillations were tentatively observed.

We repeated the GP analysis for Order 2 data, and found that the measured amplitude of granulation is lower in the Order 1 data than the Order 2 data, with $A_{\mathrm{O1}} = 71.9^{+4.9}_{-4.56}$ against $A_{\mathrm{O2}} = 84.1^{+6.8}_{-5.8}$, in alignment with expectations from the literature ([58, 55, 80]). These predict that the amplitudes of convective processes should decrease with wavelength, because the contrast between intergranular lanes and granule surfaces induces intensity differences set by the Planck function, which grows less sensitive to temperature at longer wavelengths.

The high signal-to-noise of the granulation allows us to measure its amplitude as a function of wavelength in multiple wavelength bins. To do this, we binned the 2048 individual pixel-level lightcurves, obtained as detailed in J. Splinter et al. [79] (again without the GP model removed), into 5 wavelength bins, centered at 1.05 $\mu$m, 1.44 $\mu$m, 1.84 $\mu$m, 2.24 $\mu$m, and 2.63 $\mu$m. The light curves and smoothed power spectra for each bin are shown in Figure 4. We found that splitting the data into these 5 bins reduced the SNR, and therefore we could not statistically verify the presence of the oscillation power excess at all wavelengths. The excess can be seen visually in the first two bluest wavelength bins, however becomes indistinguishable from the white noise at redder wavelengths.

Therefore, we fit the GP with one kernel representing the granulation, fixing $\nu_{\mathrm{gran}} = 410\ \mu\mathrm{Hz}$ from the integrated Order 1 data, to avoid any shorter-timescale variation affecting the estimation of granulation timescale. We applied the GP methodology to the 5 binned spectral time series from Order 1, in addition to the white light Order 2 data, this time fixing frequency at the value measured from the integrated Order 1 data, and only allowing amplitude to vary.

Figure 5 shows the GP granulation amplitude as a function of wavelength. The fifth and sixth bins yielded large error bars due to the higher RMS, so should be interpreted with caution. Nevertheless, a clear decrease in amplitude is observed with wavelength.

These measurements confirm that the amplitudes of granulation decrease with wavelength. In order to compare this result to predictions, we calculated expected relative amplitudes assuming a blackbody following the methodology from M. N. Lund [55] for each wavelength bin of the Order 1 data, and for the integrated data.

In order to compare the observed amplitude decrease to theoretical predictions, we calculated the bolometric correction factor following M. N. Lund [55]. The bolometric amplitude of a signal must be divided by this factor to obtain the predicted amplitude at a given bandpass, such that

where $A_{\mathrm{bol}}$ is the bolometric amplitude, and $A_{\mathrm{bp}}$ is the amplitude in a given bandpass. This bolometric factor allows us to build a simple wavelength-dependent model of the predicted amplitude dependence of a convective signal, like granulation.

If we model the atmosphere of a star as a blackbody (Planck) function, then it follows

In flux density units this is given by

The derivative of this with respect to $T_{\mathrm{eff}}$ is:

We use the above equations, along with the definition of a bolometric correction factor from M. N. Lund (2019), to define:

Using the spectral response function for the NIRISS Order 2 spectral trace, and the $T_{\mathrm{eff}}$ from D. K. Sing et al. [77], we calculate $C_{\mathrm{bol}}$ for each wavelength bin used in our analysis. These values are used to plot the blue curve in the right panel of Figure 5. As a test of how the amplitude scaling affects the results, we exclude the Order 2 data and repeat the calculation and comparison for only the 5 Order 1 wavelength bins. This resulted in a slightly reduced deviation between measurement and prediction, however the decrease in amplitude with wavelength remained.

We also calculate values for the bolometric correction factor across the Order 1 wavelengths for a range of $T_{\mathrm{eff}}$ between 4000 and 8000K. We present these factors in Table 1 for use by the community, and plot them in the left-hand panel of Figure 5. The predictions made using this model, alongside the observations presented in this paper, reinforce the need to model the stellar variability independently in each spectral channel.

$T_{\mathrm{eff}}$ [K]	$c_{1.05}$	$c_{1.44}$	$c_{1.84}$	$c_{1.975}$	$c_{2.17}$	$c_{0.86}$
4000	1.1408	1.4451	1.7376	1.9758	2.1702	0.8600
4444	1.2331	1.5502	1.8475	2.0853	2.2768	0.9323
4889	1.3215	1.6484	1.9483	2.1843	2.3723	1.0032
5333	1.4057	1.7402	2.0408	2.2741	2.4581	1.0724
5778	1.4859	1.8258	2.1259	2.3558	2.5356	1.1399
6222	1.5623	1.9059	2.2044	2.4305	2.6058	1.2055
6667	1.6349	1.9809	2.2769	2.4989	2.6698	1.2691
7111	1.7039	2.0510	2.3440	2.5617	2.7283	1.3308
7556	1.7696	2.1168	2.4063	2.6196	2.7820	1.3905
8000	1.8321	2.1786	2.4642	2.6732	2.8314	1.4483

We find that the decrease in granulation amplitude broadly matches the predictions from blackbody models, which serves as observational evidence that the changes to stellar surfaces due to convective processes can be approximated as small adiabatic temperature fluctuations ([72]). The data appear to follow a slightly different functional form to the predictions, with larger amplitudes at $\approx1.5\ \mathrm{\mu m}$, possibly due to variability of spectral features due to granulation which have been predicted in 3D models ([93, 78] (D. Dravins & H.-G. Ludwig 2023)).

We further validated the astrophysical nature of the stellar signals by measuring the correlation of total flux across wavelengths. A Pearson correlation test between consecutive spectral light curves produced correlation coefficients of 0.57, 0.38, 0.25 and 0.18, all with P-values $< 10^{-25}$. This confirms that the temporal parameter of the observed variability is, on average, consistent at all wavelengths. To isolate the granulation signals from oscillation signals, each time series was smoothed and high-pass filtered around a frequency of 1000$\mathrm{\mu}$Hz. We then fitted a model of the form $y=kx$ between the flux of the bluest band and each subsequent band. This resulted in values of K = 0.663 $\pm$ 0.085, 0.391 $\pm$ 0.107, 0.405 $\pm$ 0.163 and 0.337 $\pm$ 0.279. These provide simplified estimates for granulation amplitudes, which decrease with wavelength, in agreement with predictions and with our GP results. The difference in absolute values between this method and the GP fit can likely be attributed to the coarse signal isolation method, and the different treatment of white noise.

Implications

The wavelength dependence of convection-driven variability has important implications across a range of astrophysics. Additional comparisons between observations in the infrared and theoretical scaling relations ([72]) will shed light on how well convection can be modeled as simple perturbations that only depend on temperature, and to what level 3D models accurately predict the time and wavelength dependence of granulation. Improved constraints on 3D models directly translate into 1D stellar structure model grids ([54, 87, 59, 62, 49]), which are widely used to study stellar populations within and outside of our galaxy.

The chromaticity of stellar variability is also important for transit spectroscopy studies. Unmodelled granulation can systematically bias planet radius measurements ([1, 66]), and the resulting error floor on planet radius can reach 10%, for an Earth-sized planet ([81]). If planets transit on chords that have different brightnesses due to granulation, this can bias the interpretation of temporal atmospheric variations made from e.g. phase-curve analysis ([14]). Incorrectly modelling these effects at multiple wavelengths could lead to wavelength-dependent biases in the results [66, 28] because transit depths measured at shorter wavelengths could suffer from more significant granulation noise floors than longer wavelengths [73]. Observations made at longer wavelengths are predicted to suffer from 10–20% transit depth error floors as a result of stellar variability (B. Morris & D. Huber, in prep).

Finally, the wavelength dependence of oscillations is critical for estimating asteroseismic yields of future space-based missions. For example, the Nancy Grace Roman Telescope [63] will obtain continuous near-infrared time-domain observations in the Galactic Bulge, which is expected to detect oscillations in hundreds of thousands of red giants [89, 46]. The results presented here support the assumptions made to predict those yields, and thus open up the opportunity to analyze variations in asteroseismic yields due to astrophysical effects such as stellar activity and close binaries [33, 34, 32].

Conclusion

In this paper, we used JWST/NIRISS SOSS observations of an exoplanet host star to presented the first detections of stellar granulation and solar-like oscillations at infrared wavelengths. By modelling the time-series data with a Gaussian-process framework, we have estimated granulation timescale and oscillation envelope parameters for WASP-121, which match predictions well. These results demonstrate the suitability of JWST for resolving near-surface variability processes in the infrared. We also demonstrate, for the first time using observations, the chromaticity of convective processes in sun-like stars. This finding is crucial for our understanding of 3D stellar structure, mitigating stellar variability in multi-wavelength exoplanet observations, and for predicting the yields of future infrared asteroseismology surveys.

ACKNOWLEDGMENTS

This work was supported by the Australian Research Council (DP 250102562). D.H. acknowledges support from the National Aeronautics and Space Administration (80NSSC25K7155, 80NSSC25K0177) and NASA’s Interdisciplinary Consortia for Astrobiology Research (NNH19ZDA001N-ICAR) under award number 19-1CAR19 20-0041. DL would like to acknowledge financial support from the Canadian Space Agency. Some/all of the data presented in this article were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via doi:10.17909. We acknowledge and appreciate the comments provided by the anonymous reviewer, which strengthened the quality of the manuscript.

Facilities: JWST (NIRISS)

Software: The work presented in this letter has made use of Astropy, [65], Matplotlib [47], NumPy [39], pandas [90] (The pandas development team 2020), and Celerite [27].