Granulation in Color: The Chromatic Stellar Flicker of LTT-9779 and WASP-121

Introduction

Granulation is the surface phenomenon associated with near-surface convective energy transport in stellar photospheres (e.g., [53]). On the Sun, this phenomenon was first reported through visual observations by Herschel [29] and others, and later captured in detail through drawings [60] and short-exposure photography [32, 33]. Solar granulation subsequently became a key observable signature of photospheric convection. Granulation produces stochastic brightness variations as hot plasma rises, cools radiatively near the stellar surface, and descends through darker intergranular lanes (Dravins 1982; [53]). Stars also support acoustic oscillations excited by granulation, contributing a stochastic variability on the timescale of minutes ([44]). In disk-integrated observations, these phenomena introduce stellar variability with characteristic amplitudes and timescales. For low-activity stars like the Sun, granulation and oscillations are among the dominant contributors to disk-integrated photometric variability on timescales of minutes to hours, while magnetic activity becomes increasingly important on rotational and longer timescales ([8, 61]). These convective processes can hinder the detection and characterization of exoplanets ([2] Cegla [12]; Meunier & Lagrange [49] by raising the stellar noise floor and creating a stochastic stellar “flicker”.

Flicker affects exoplanet characterization through both spectroscopic and photometric channels. In high-resolution spectroscopy, it alters spectral-line shapes, asymmetries, and wavelength shifts, introducing time- and wavelength-correlated radial-velocity (RV) variations that can obscure the mass measurement of Earth-like planets orbiting Sun-like stars (e.g., [27, 23, 10] (Schrijver & Zwaan 2000)). Simulations of RV surveys indicate that granulation can increase the uncertainty of Earth-analogue mass measurements by approximately a factor of three compared with a case without stellar variability ([46]). In transit observations, neglecting time-correlated stellar flicker can introduce errors of up to 10% in the inferred planet-to-star radius ratio for an Earth-sized planet orbiting a Sun-like star ([66]). The wavelength dependence of this variability can also affect transit depths measured in individual spectral channels and, consequently, the inferred planetary transmission spectrum. Recent transit studies with JWST’s NIRISS/SOSS instrument ([16]) have found signals attributed to granulation with characteristic amplitude of ∼50-100 ppm in the light curves of main sequence stars ([11, 14, 63, 54]). With the extreme precision of space telescopes and ground-based RV spectrographs, accounting for stellar variability has become an essential aspect of exoplanet discovery and characterization.

Additionally, granulation and acoustic oscillations are important observable phenomena for constraining fundamental stellar parameters and understanding the physics of stellar interiors ([51, 1]). Granulation amplitudes and timescales depend on the fundamental stellar parameters; effective temperature ($T_{\mathrm{eff}}$), surface gravity ($\log g$), and metallicity ([M/H]), as well as magnetic activity ([69, 6, 45, 18, 48, 39, 35, 15, 9, 13, 67, 59]). The intensity of the granulation pattern scales with $T_{\mathrm{eff}}$ and [M/H] while the number of granules on the surface of the star, $N_{\mathrm{granules}}$, scales with $\log g$. The disk-integrated amplitude of stellar flicker scales directly with the intensity contrast of the granulation pattern and inversely with the square root of the number of independently varying granules $1/\sqrt{N_{\mathrm{granules}}}$. In other words, the observed effect decreases for cooler stars because they have a stronger surface gravity and more granules.

Cranmer et al. [15] suggested that increased surface magnetism would reduce the magnitude of stellar flicker because magnetic fields suppress convection. Using three-dimensional radiative-MHD simulations, Kostogryz et al. [40] studied the wavelength-dependence of granulation for K2, G2, and F3 stellar models, including G2 models with different metallicities. They showed that magnetic fields decrease the intensity contrast of granulation by suppressing convective temperature fluctuations while also enhancing the intensity contrast from localized magnetic structures, especially near the limb. In disk-integrated spectra, these effects compete with each other to modulate the observed variability depending on stellar type, metallicity, and wavelength. Thus, chromatic granulation provides a promising diagnostic of stellar photospheric properties, surface magnetism, and near-surface convection.

The vast majority of stellar granulation research has focused on single-bandpass light curves of evolved stars, with only recent work by Sreenivas et al. [64] taking advantage of TESS and Kepler light curves to perform a two-channel chromatic analysis of granulation and oscillation and demonstrating evidence of the predicted decrease in granulation amplitude with wavelength ([38, 47]). Studying the wavelength-dependence of convective processes can help us clarify their driving mechanisms and provide a systematic improvement to the study of exoplanet transmission spectra. Here, we aim to detect and characterize spectral granulation signatures in continuous high-cadence time-series spectroscopy from JWST NIRISS/SOSS. This work starts a series on “Granulation in Color” (AR #11527 (PIs: J. Valenti & A. Glidden) which will use the unprecedented spectral sensitivity and photometric precision of JWST to compare wavelength-dependent stellar flicker with 3D radiative MHD models.

This paper is laid out as follows: in § 2, we describe the target selection process and how we reduced and conditioned the JWST observations in order to study their power spectra. In § 3, we present the results of Bayesian inference on the power spectra, which lead to measurements of the granulation timescale and amplitude. We further compare the chromatic variation in amplitude to theoretical predictions from stellar granulation simulations. Finally, in § 4, we place these results in the context of previous studies of solar and stellar granulation.

Methods

Target Selection and Data Reduction

LTT-9779 and WASP-121

We chose these targets from the list of publicly available JWST transit observations compiled in the Transiting Exoplanets List of Space Telescope Spectroscopy (TrExoLiSTS) database ([52]). Starting with the full list of all publicly available transmission spectrum observations from JWST, we narrowed the search to F- and G-type stars as their flicker amplitudes are greater than those of cooler stars. From those observations we selected our final targets as those with the longest baselines – NIRISS/SOSS phase curve observations of LTT-9779 (21-hour baseline with 16.5s cadence) and WASP-121 (37-hour baseline with 38.5s cadence). Instrumental and stellar parameters relevant to this study are listed in Table 1.

Parameter	LTT-9779	WASP-121
$T_{\mathrm{eff}}$ [K]	5445 ± 84	6628 ± 66
$M_{\star}$ [$M_{\odot}$]	1.03 ± 0.04	1.385 ± 0.016
$R_{\star}$ [$R_{\odot}$]	0.95 ± 0.01	1.461 ± 0.005
M/H [dex]	0.25 ± 0.08	0.17 ± 0.05
$\log g$ [cgs]	4.43 ± 0.31	4.251 ± 0.003
Cadence [s]	16.5	38.5
Baseline [hr]	20.6	36.9

The NIRISS/SOSS observations of WASP-121 and LTT-9779 analyzed in this work are both phase curves obtained as part of the JWST GTO Program 1201 ([42]). The observations of WASP-121 were taken in October 2023 using the SUBSTRIP256 subarray and published in Splinter et al. [63]. The observations of LTT-9779 were taken in July 2022 using the SUBSTRIP256 subarray and published in Coulombe et al. (2025). We downloaded the observations from the Mikulski Archive for Space Telescopes (MAST) and reduced the data independently using the exoTEDRF software (Feinstein et al. [19]; Radica et al. [57]; Radica [56]).

Data reduction with exoTEDRF consists of three stages. Stage 1 is detector-level processing which performs bias and dark subtraction, 1/f noise correction, and cosmic ray rejection. Stage 2 is the spectroscopic processing where flat-fielding, background subtraction, bad pixel masking, and a principle component analysis (PCA) are performed. Eigenvalues from the PCA analysis are used to de-trend the data. The first component of the PCA corresponds to the astrophysical components of the observation and are not removed, so we used only components 2, 3, and 4. These components correspond to drifts in the position of the trace (components 2 and 3) and the characteristic beating pattern induced by the thermal cycling of the telescope that is common to all NIRISS/SOSS observations (component 4). Stage 3 of the exoTEDRF reduction is where the 1D spectrum is extracted from along the detector with a simple box extraction using exoTEDRF’s built-in optimization to determine the ideal box width (30 pixels in our case). This software also has a fourth stage which will perform light curve fitting, but we perform this part ourselves. Further details of this pipeline and its procedure can be found in the exoTEDRF tutorial² and are discussed in (Radica et al. [57]; Radica [56]).

We trim the noisier wavelength regions near the edges of the detector (blue-ward of 0.62 $\mu$m and red-ward of 2.82$\mu$m). We then model the resulting light curves by accounting for planetary transit and eclipses, phase variations in the planet’s thermal emission and reflected light, stellar rotational modulation, and instrumental systematics. The exoplanet transit and eclipse models are calculated with batman (Kreidberg [41]) and an optimized set of parameters is found using scipy.optimize. Details about the transit model can be found in Appendix A. We subtract the best-fitting model from each light curve to obtain the residual time series (Figure 1) which we use to investigate stellar granulation and oscillations in the presence of observational noise and any remaining instrumental systematics.

Solar Observations

To provide a solar reference for interpreting the granulation and oscillation signatures measured in LTT-9779 and WASP-121, we also analyze Sun-as-a-star observations from the VIRGO experiment on board the Solar and Heliospheric Observatory (SOHO). In particular, we use VIRGO/PMO6 Total Solar Irradiance (TSI) measurements (Fröhlich et al. [24]) obtained at a cadence of 1 min between February 1996 and October 2018. We divide the more than 8000 days of observations into 48-hour segments, calculate the power spectral density of each segment, and average the resulting spectra. This procedure reduces the stochastic scatter in the solar power spectrum while using individual time series with durations comparable to the JWST observations.

Power Spectrum Analysis

We calculate the discrete fast Fourier transform of the de-trended light curves with numpy.fft.fft and convert that Power Spectrum (PS) into a Power Spectral Density (PSD):

where $\nu$ is frequency, $\Delta t$ is the time cadence between exposures, and N is the number of exposures. For the chromatic analysis, we chose 6 wavelength bins and summed the light curve along wavelength for each bin. Bin edges were defined by eye such that there was similar scatter in the residual (de-trended) time series for each bin.

We model the power spectrum, $P(\nu)$, as a sum of Harvey-like profiles [28, 65], a white-noise term, and a Gaussian envelope representing the power excess due to solar-like oscillations:

where $\nu$ is the linear frequency, $P'$ is a frequency-independent white-noise component, and $N_{\mathrm{components}}$ is the number of Harvey-like granulation components, taken to be one or two in this work. For each component, $A_i$ sets the lower-frequency plateau of the power spectral density, $B_i$ determines the turnover frequency and is related to the characteristic timescale, and $C_i$ controls the high-frequency slope of the profile. The oscillation component, $P_{\mathrm{osc}}$, is modeled as a Gaussian envelope:

where $H_{\mathrm{osc}}$ is the peak power spectral density of the envelope, $\nu_{\max}$ is its central frequency, and $\sigma_{\mathrm{env}}$ is the envelope width. We relate model parameters $A_i$ and $B_i$ to their physical counterparts in the time domain, $\sigma_{\mathrm{gran}}$ (RMS) and $\tau_{\mathrm{gran}}$ (s), with

and

where $\xi_{C_i}$ is a normalization constant which satisfies

For $C_i=2$ or $C_i=4$, analytic solutions for $\xi_{C_i}$ work out to

and

respectively [35]. For other values of $C_i$, this integral must be evaluated numerically.

The exponent $C$ controls the high-frequency slope of the Harvey-like profile: $C=2$ is routinely applied in the literature and corresponds to an instantaneously appearing and gradually decaying pattern of granulation variability, while $C=4$ is related to a symmetric rise and decay which exhibits a steeper slope in the power spectrum. Based on a statistical comparison of one-versus two-component fits and a variety of treatments of $C$ (see Appendix § 8), we focus our results and discussion on three models. Summarized in Table 2: a one-component model with $C=2$ (hereafter $\mathcal{M}1$), a one-component model with $C=4$ ($\mathcal{M}2$), and a two-component model with $C \in\{2,4\}$ (hereafter $\mathcal{M}3$). Although $\mathcal{M}3$ includes two Harvey-like components, we do not interpret the second, longer-timescale component as a secure detection of an additional granulation signal. Instead, we treat it as a contribution that accounts for excess low-frequency power that may arise from residual instrumental systematics, imperfect subtraction of planetary signals, or other stellar variability.

Model	$N_{\mathrm{components}}$	$C$
$\mathcal{M}1$	1	2
$\mathcal{M}2$	1	4
$\mathcal{M}3$	2	High-frequency = 2, Low-frequency = 4

MCMC Setup

We infer the power spectrum parameters using Bayesian Inference. We use the Goodman & Weare [26] affine-invariant Markov Chain Monte Carlo (MCMC) with the emcee [22] ensemble sampler to compute the marginalized probability functions. In the minimization we use the likelihood function in the form of the exponential distribution, which is appropriate for power spectra fitting (see equations (9) and (10) of Larsen et al. [43]). Models are run with 50 walkers, 60000 steps, and 30000 burn-in steps. As shown in Table 3 all parameters have uniform priors with wide bounds with the exception of $\nu_{\max}$ and $\sigma_{\mathrm{env}}$ which we restrict to reasonable ranges to avoid spurious solutions far from their expected values.

Parameter	LTT-9779	WASP-121
$\log\sigma_{\mathrm{noise}}$	$\mathcal{U}(4, 7)$	$\mathcal{U}(4, 7)$
$\log A_1$	$\mathcal{U}(5, 8)$	$\mathcal{U}(5, 8)$
$\log B_1$	$\mathcal{U}(2.5, 4)$	$\mathcal{U}(2.5, 4)$
$\log A_2$	$\mathcal{U}(4, 8.0)$	$\mathcal{U}(4, 8)$
$\log B_2$	$\mathcal{U}(B_1, 4.5)$	$\mathcal{U}(B_1, 4.5)$
$H_{\mathrm{osc}}$	$\mathcal{U}(4, 7)$	$\mathcal{U}(4, 7)$
$\nu_{\max}$	$\mathcal{U}(2800, 4200)$	$\mathcal{U}(1300, 2100)$
$\sigma_{\mathrm{env}}$	$\mathcal{U}(100, 500)$	$\mathcal{U}(10, 200)$

Results

Solar Power Spectrum

As a benchmark for our analysis, we fit the solar power spectral density with $\mathcal{M}1$, which includes a single Harvey-like granulation component with $C=2$, a Gaussian oscillation envelope, and a white-noise term. Before fitting, we mask a prominent instrumental peak and two associated harmonics which are caused by the time series sampling frequency. The fitted spectrum is shown in Figure 2, and the inferred parameters are listed in Table 4. We adopt this model as the solar reference without performing the additional model comparisons applied to the JWST targets.

Parameter	Sun	LTT-9779	WASP-121
$\mathcal{M}1$
$\sigma_{\mathrm{gran}}$ [ppm]	$43.6^{+1.5}_{-1.4}$	$55.3^{+3.8}_{-3.5}$	$62.8^{+3.9}_{-3.6}$
$\tau_{\mathrm{gran}}$ [min]	$3.4^{+0.5}_{-0.4}$	$5.6^{+1.5}_{-1.1}$	$8.5^{+2.0}_{-1.5}$
$\nu_{\max}$ [$\mu$Hz]	$3129^{+97}_{-76}$	$3377^{+164}_{-191}$	$1647^{+80}_{-67}$
$\sigma_{\mathrm{env}}$ [$\mu$Hz]	$372^{+99}_{-82}$	$395^{+174}_{-136}$	$83^{+80}_{-42}$
$\mathcal{M}2$
$\sigma_{\mathrm{gran}}$ [ppm]		$50.0^{+3.4}_{-3.2}$	$57.2^{+3.6}_{-3.3}$
$\tau_{\mathrm{gran}}$ [min]		$3.2 \pm0.4$	$4.6 \pm0.5$
$\nu_{\max}$ [$\mu$Hz]		$3325^{+161}_{-197}$	$1633^{+48}_{-61}$
$\sigma_{\mathrm{env}}$ [$\mu$Hz]		$466 \pm156$	$109^{+81}_{-48}$
$\mathcal{M}3$
$\sigma_1$ [ppm]		$53.8^{+4.1}_{-4.0}$	$59.3 \pm4.3$
$\tau_1$ [min]		$4.8^{+1.5}_{-1.3}$	$6.3^{+2.1}_{-1.8}$
$\sigma_2$ [ppm]		$14.0^{+15.9}_{-11.3}$	$26.6^{+11.5}_{-8.9}$
$\tau_2$ [min]		$33.0^{+54.2}_{-20.0}$	$38.1^{+48.5}_{-16.0}$
$\nu_{\max}$ [$\mu$Hz]		$3382^{+166}_{-191}$	$1649^{+53}_{-69}$
$\sigma_{\mathrm{env}}$ [$\mu$Hz]		$387^{+179}_{-137}$	$82^{+82}_{-44}$
Predicted
$\sigma_{\mathrm{gran}}$ [ppm]		$42 \pm2$	$57 \pm5$
$\tau_{\mathrm{gran}}$ [min]		$3.2 \pm0.4$	$6 \pm1$
GP
$\sigma_{\mathrm{gran}}$ [ppm]		$66 \pm4$	$59 \pm4/67 \pm5$
$\tau_{\mathrm{gran}}$ [min]		$3.3^{+0.5}_{-0.3}$	$4.8 \pm0.6/5.6^{+1.1}_{-0.8}$

Bandpass-Integrated Stellar Power Spectra

While most instrumental systematics are removed in the data processing step, residual instrumental signatures appear as strong peaks in the power spectra (see Figure 3). Identified by Albert et al. [3], these peaks at 3.40 and 3.63 minutes (204 and 218 seconds or 4902 and 4587 $\mu$Hz) are due to slight forces exerted on the instrument when thermostats are cycled on and off. For both targets, we mask the primary signature as well as harmonic signatures at $2\nu_{\mathrm{Inst}}$ and $\frac{1}{2}\nu_{\mathrm{Inst}}$. Additional harmonic signatures may be present but are difficult to identify above the noise.

Table 4 summarizes the parameters inferred from the three models. The Bayesian information criterion favors $\mathcal{M}1$ for both JWST targets, although the Akaike information criterion provides little discrimination between $\mathcal{M}1$ and $\mathcal{M}3$ for WASP-121 (see Appendix § 8). Model $\mathcal{M}2$ has the same granulation-profile shape as the GP kernel used by Coulombe et al. [14]; Splinter et al. [63], while $\mathcal{M}3$ includes an additional component to account for excess low-frequency power that is not necessarily attributable to granulation. Figure 3 shows the fits obtained with each models for LTT-9779 and WASP-121.

The timescale of the second component in $\mathcal{M}3$ is poorly constrained. Following the interpretation outlined in Sect. 2.2, we compare only the first, shorter-timescale component with the granulation components of the single-component models.

Chromatic Measurements

The chromatic granulation amplitude modeling are shown in Figure 4 and Table 5 and Table 6. For each model, the characteristic timescale or timescales were fixed to the corresponding posterior medians inferred from the bandpass-integrated PSD. The reported amplitude uncertainties are therefore conditional on these fixed values. Both stars show an overall decrease in granulation amplitude towards longer wavelengths across all three model prescriptions.

$\lambda$ (µm)	$\Delta\lambda$ (µm)	$\sigma_1$ (ppm)	$\sigma_2$ (ppm)	$\sigma_{\mathrm{total}}$ (ppm)
$\mathcal{M}1$
0.735	0.115	$59.0^{+4.9}_{-4.7}$
0.975	0.125	$54.3^{+4.8}_{-4.6}$
1.200	0.100	$52.6^{+4.6}_{-4.3}$
1.400	0.100	$52.2^{+4.7}_{-4.5}$
1.610	0.110	$49.9^{+5.5}_{-5.3}$
2.270	0.550	$45.2^{+4.9}_{-4.8}$
$\mathcal{M}2$
0.735	0.115	$54.0^{+4.5}_{-4.3}$
0.975	0.125	$49.6^{+4.4}_{-4.2}$
1.200	0.100	$47.5^{+4.1}_{-3.9}$
1.400	0.100	$47.2^{+4.2}_{-4.0}$
1.610	0.110	$45.1^{+4.9}_{-4.7}$
2.270	0.550	$41.4^{+4.5}_{-4.4}$
$\mathcal{M}3$
0.735	0.115	$57.9^{+5.1}_{-4.9}$	$16.0^{+16.8}_{-13.0}$	$74.0^{+15.7}_{-11.4}$
0.975	0.125	$50.1^{+5.3}_{-5.1}$	$29.3^{+11.7}_{-11.5}$	$78.9^{+12.1}_{-11.1}$
1.200	0.100	$49.8^{+5.0}_{-4.9}$	$26.8^{+13.8}_{-14.9}$	$76.2^{+13.4}_{-12.9}$
1.400	0.100	$51.7^{+4.8}_{-4.7}$	$11.0^{+15.0}_{-8.7}$	$63.5^{+13.7}_{-8.8}$
1.610	0.110	$50.3^{+5.6}_{-5.6}$	$8.1^{+16.6}_{-6.3}$	$60.3^{+13.8}_{-8.5}$
2.270	0.550	$40.9^{+6.2}_{-6.4}$	$27.8^{+14.6}_{-17.4}$	$67.8^{+13.8}_{-13.6}$

$\lambda$ (µm)	$\Delta\lambda$ (µm)	$\sigma_1$ (ppm)	$\sigma_2$ (ppm)	$\sigma_{\mathrm{total}}$ (ppm)
$\mathcal{M}1$
0.735	0.115	$66.6^{+4.6}_{-4.5}$
0.975	0.125	$58.5^{+4.4}_{-4.6}$
1.200	0.100	$55.2^{+4.7}_{-4.0}$
1.400	0.100	$56.0^{+4.1}_{-4.8}$
1.610	0.110	$52.3^{+5.2}_{-4.6}$
2.270	0.550	$50.7^{+5.0}_{-4.9}$
$\mathcal{M}2$
0.735	0.115	$59.7^{+3.3}_{-4.9}$
0.975	0.125	$53.2^{+4.0}_{-4.1}$
1.200	0.100	$49.8^{+3.2}_{-4.5}$
1.400	0.100	$50.9^{+4.0}_{-4.0}$
1.610	0.110	$47.6^{+4.0}_{-4.6}$
2.270	0.550	$46.9^{+4.4}_{-4.4}$
$\mathcal{M}3$
0.735	0.115	$65.9^{+7.7}_{-1.2}$	$11.5^{+2.5}_{-20.8}$	$78.1^{+4.8}_{-16.7}$
0.975	0.125	$58.1^{+5.9}_{-3.5}$	$14.1^{+4.4}_{-18.7}$	$72.1^{+4.6}_{-15.5}$
1.200	0.100	$56.0^{+5.4}_{-3.3}$	$5.2^{+3.4}_{-14.1}$	$64.6^{+3.9}_{-13.0}$
1.400	0.100	$56.5^{+4.7}_{-3.9}$	$6.0^{+2.7}_{-12.0}$	$63.7^{+3.1}_{-12.2}$
1.610	0.110	$53.7^{+6.6}_{-3.2}$	$5.7^{+3.0}_{-11.8}$	$61.0^{+4.8}_{-11.1}$
2.270	0.550	$51.9^{+3.5}_{-6.3}$	$3.4^{+8.2}_{-3.4}$	$56.5^{+9.0}_{-3.6}$

For WASP-121, $\mathcal{M}1$ and the first component of $\mathcal{M}3$ yield similar amplitudes, while $\mathcal{M}2$ gives amplitudes approximately 6 ppm lower and follows a similar wavelength dependence. For LTT-9779, $\mathcal{M}2$ also yields systematically lower amplitudes than $\mathcal{M}1$. The first component of $\mathcal{M}3$ shows a less smooth wavelength dependence, with a modest local enhancement near 1.4–1.5 $\mu$m. However, this structure is not clearly resolved within the reported uncertainties.

Discussion

Oscillation and Granulation Scaling Relations

We use solar scaling relations to predict the oscillation and granulation parameters of our targets. The frequency of maximum oscillation power, $\nu_{\max}$, is predicted using the scaling relation of Kjeldsen & Bedding [38]:

where $M$, $R$, and $T_{\mathrm{eff}}$ are the stellar mass, radius, and effective temperature, respectively, and the subscript $\odot$ denotes solar reference values. For the characteristic granulation timescale and RMS amplitude, we use the relations given by Gilliland et al. [25]:

and:

When calculating predictions from the scaling relations, we use the median parameters measured from $\mathcal{M}1$ fits of the solar data

(also found in Table 4):

and

From these we predict a granulation amplitude of $42 \pm2\mathrm{ppm}$ with timescale of $3.2 \pm0.4$ minutes for LTT-9779. For WASP-

121, we predict an amplitude of $\sigma_{\mathrm{gran}} = 57 \pm5\mathrm{ppm}$ and timescale of $\tau_{\mathrm{gran}} = 6 \pm1$ minutes.

Interpretation of Timescale

The timescale we measure for the Sun is $3.4^{+0.5}_{-0.4}$ minutes, which agrees with the $\sim3.7$ minute timescale measured in the VIRGO PSD by Gilliland et al. [25], but is notably shorter than the 5-6 minute lifetime of large individual solar granules (Title et al. 1989; Hirzberger et al. [30]. On the other hand, Kallinger et al. [35] fit the VIRGO PSD and measured $\tau_{\odot}$ to be 6.3 minutes.

We emphasize that the timescale inferred from the fit does not necessarily correspond to the full lifetime of a granule. Instead, it may characterize a particular phase of granular evolution during which the brightness changes most strongly. Fits obtained with different values of $C$ may therefore probe different phases of a granule’s lifetime.

Oscillations

The scaling-relation predictions provide a consistency check on whether the fitted excess power can be attributed to solar-like oscillations. Using the stellar parameters in Table 1 and Equation (4), we predict $\nu_{\max} = 3664 \pm143\ \mu\mathrm{Hz}$ for LTT-9779, compared with our measured value of $3382^{+166}_{-191}\ \mu\mathrm{Hz}$. For WASP-121, the predicted frequency is $1888 \pm77\ \mu\mathrm{Hz}$, whereas the measured value is $1649^{+55}_{-69}\ \mu\mathrm{Hz}$.

The interpretation of these measurements is complicated by narrow peaks of excess power around 2200–2900 $\mu\mathrm{Hz}$, which may be associated with a subharmonic at half the frequency of the NIRISS instrumental beat signature. Our Gaussian model describes the collective power of the modes of oscillation modes rather than individual modes. Masking potentially instrumental peaks can therefore also remove stellar oscillation power, reducing the information available to constrain the central frequency and width of the oscillation envelope.

For LTT-9779, the measured $\nu_{\max}$ is broadly consistent with the scaling-relation prediction, supporting a stellar origin for the fitted excess power. For WASP-121, the measured frequency is lower than predicted, and possible contamination by instrumental features makes the interpretation less secure. We therefore

@suppress regard the identification of solar-like oscillations in WASP-121 as tentative.

Granulation Parameters in GP and Power Spectrum Analysis

Recent investigations have analyzed granulation signatures directly in the light curve using Gaussian Processes (GP) in order to enhance transit modeling [55, 7, 11, 14, 63, 54]. Here we compare the granulation amplitudes and characteristic timescales inferred from our analysis with both scaling-relation predictions and published GP-based measurements. To do this comparison we must first translate the parameters reported by the GP analyses, which used the angular-frequency SHO of the celerite package Foreman-Mackey et al. [21], into the linear frequency space of the Harvey model. We also attempt to reconcile the fact that the timescale and amplitude measured by a GP are larger than predicted from scaling relations and measured for the Sun (see discussion in Pass et al. [54]).

Timescale

We start by examining the timescales. Part of the apparent discrepancy comes from a difference between the timescale definitions adopted by and in our Harvey-like power-spectrum model. In the GP approach, the exact form of the simple harmonic oscillator (SHO) profile used by celerite is:

which, for $Q=1/\sqrt{2}$, can be simplified to:

with the timescale defined as

If we compare the denominator of eqn. 8 to the Harvey profile denominator in eqn. 2, we can see that $\omega/\omega_0 = B\nu$. With the definition that linear frequency $\nu$ corresponds to angular frequency $\omega$ as $\omega=2\pi\nu$, we arrive at:

In the following discussion, we have taken the reported (GP-derived) timescales from Coulombe et al. [14] and Splinter et al. [63] and divided them by $2\pi$, calculating the following timescales:

Amplitude

There is also a significant difference between the GP-derived granulation signal and both the scaling relation predictions and our results. We find a resolution by again examining the analytic forms of the simple harmonic oscillator used in our method and the GP method, arriving at the following relation (derivation provided in Appendix C):

We use eqn. (10) to convert the Coulombe et al. [14] measurement of 93 $\pm$ 6ppm to 66 $\pm$ 4ppm and the Splinter et al. [63] measurements of 83 $\pm$ 5ppm and 95 $\pm$ 7ppm to 59 $\pm$ 4ppm and 67 $\pm$ 5ppm for order 1 and order 2, respectively. Scaling relation predictions and converted GP measurements are listed alongside our measurements in Table 4.

Wavelength Dependence of Granulation

The wavelength dependence of the granulation amplitude provides an additional test of the stellar origin of the measured variability. Whereas the bandpass-integrated amplitude describes the overall photometric fluctuations, its spectral dependence probes how the contrast between bright and dark surface structures changes with wavelength. Comparing this dependence with stellar-atmosphere simulations therefore allows us to assess whether the observed variability is consistent with granulation and to explore its sensitivity to stellar properties.

We measure a wavelength-dependent decrease in amplitude for both stars (see Table 5 and Table 6). To interpret these trends, we compare our measurements with the granulation-induced amplitude spectra of Kostogryz et al. [40], calculated from three-dimensional magnetohydrodynamic simulations using the MURaM code. We use an F3 model with [M/H] = 0 as a proxy for WASP-121 and G2 models with [M/H] = 0 and +0.5 as proxies for LTT-9779. These are the closest available MURaM-simulated models to the stellar parameters of our targets. However, their surface gravities and effective temperatures differ from those of the observed stars. Because granulation-induced variability amplitude depends strongly on surface gravity, this mismatch limits the interpretation of differences in absolute amplitude. Instead, we focus on a qualitative comparison of the spectral shapes. A quantitative comparison will require simulations matched to the surface gravity, effective temperature, and metallicity of each target.

For LTT-9779, the predicted simulated amplitudes are significantly lower than what we measure but reproduce the general spectral trend. Notably, the enhanced metallicity model (purple) appears more similar to that of the spectral profile we measure, demonstrating the promise of characterizing stellar atmospheres from spectral granulation profiles.

For WASP-121 we do not have an enhanced metallicity simulation available. The simulated profile shows agreement with the overall level of signal strength albeit with a much steeper slope than what we measure.

Summary & Conclusions

We analyzed the bandpass-integrated and wavelength-resolved power spectra of JWST/NIRISS/SOSS observations of LTT-9779 and WASP-121 to characterize their granulation signatures. Our main findings are as follows.

• We measured granulation amplitudes and characteristic timescales from bandpass-integrated timeseries of the G7V-type star LTT-9779 and the F6V-type star WASP-121 using Harvey-like power spectrum models. Comparing one- and two-component models as well as different treatments of C, we can see that there is general agreement between model results, with slightly smaller timescales and amplitudes measured for $\mathcal{M}2$ compared to $\mathcal{M}1$ and $\mathcal{M}3$. For $\mathcal{M}3$, we interpret the second longer-timescale component as a contribution from excess low-frequency power (Figure B.1) rather than as a secure detection of an additional granulation component. Under this interpretation, we recommend taking the measurements from the first component of $\mathcal{M}3$ as the best approximation of granulation for these observations. This gives us a granulation amplitude of $\sigma_{\mathrm{gran}} = 53.8^{+4.8}_{-4.0}$ppm with characteristic timescale $\tau_{\mathrm{gran}} = 4.8^{+1.5}_{-1.3}$min for LTT-9779 and granulation amplitude of $\sigma_{\mathrm{gran}} = 59.3^{+4.3}_{-4.3}$ppm with characteristic timescale $\tau_{\mathrm{gran}} = 6.3^{+2.1}_{-1.8}$min for WASP-121.

• We identified an excess-power envelope consistent with solar-like oscillations in LTT-9779, with $\nu_{\max} = 3382^{+166}_{-191}\ \mu\mathrm{Hz}$, broadly consistent with the scaling-relation prediction. For WASP-121, the fitted value, $\nu_{\max} = 1649^{+55}_{-69}\ \mu\mathrm{Hz}$, is lower than predicted, and we regard the oscillation identification as tentative. Potential instrumental features and the masking of affected frequency intervals limit the interpretation of these envelopes.

• We compared Coulombe et al. [14] and Splinter et al. [63] GP measurements with our results and found that in both cases, the GP measures greater amplitudes and longer timescales than we obtain. We believe this is not a physical contradiction or methodological inaccuracy, but a confusion driven by inconsistent definitions of timescale and signal amplitude. Following the discussion of this point in the text, we scale the GP measured timescales by $\frac{1}{2\pi}$ and their associated amplitudes ($a_{\mathrm{gran}}$) by $\frac{1}{\sqrt{2}}$. These conversion factors bring the GP-based measurements closer to our results and the scaling relation predictions.

• Finally, we investigated the wavelength-dependence of granulation amplitude, finding that both stars’ granulation signal decreases with wavelength in ways that are consistent with theoretical models for similar stars. From Kostogryz et al. [40] we expect granulation amplitude to increase for more active and higher metallicity stars, which we will investigate in future work as part of the “Granulation in Color” archival grant.

These results demonstrate the potential of JWST time-series spectrophotometry to characterize the spectral dependence of stellar granulation and inform its treatment in exoplanet light-curve analyses.

Acknowledgements. This work is supported by the Volkswagen Foundation. N.K., V.V., and L.G. acknowledge financial support through German Aerospace Center (DLR) grant “PLATO Data Center” #50OP1902 and European Research Council (ERC) Synergy Grant “WHOLE SUN” #810218. L.G. is supported in part by the joint ANR–DFG project BUTTERFLY under the Agence nationale de la recherche (ANR grant 23-CE46-0009-01) and the Deutsche Forschungsgemeinschaft (DFG grant 530101854). Data presented in this publication makes use of observations taken with the NASA/ESA/CSA JWST from GO # 1201 (PI D. Lafrenière), accessed through the Mikulski Archive for Space Telescopes (MAST) hosted at the Space Telescope Science Institute. We thank Wolfgang Finsler for providing the calibrated VIRGO/PMO6 observations. In addition to the software cited in the text, we gratefully acknowledge use of the open-source codes astropy [5], chromatic, corner [20], matplotlib [31], numpy [70], and scipy [71]. This analysis is associated with JWST AR #11527 (PIs: J. Valenti & A. Glidden) with support provided by NASA through a grant from the Space Telescope Science Institute, operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-03127. WCW is supported by Volkswagen Foundation grant no. 9E126. R.L. is funded by the European Union (ERC, THIRDSTEE, 101164189) and acknowledges support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. I.R. and J.B.-P. acknowledge support from the European Research Council (ERC) under the European Union’s Horizon Europe programme (ERC Advanced Grant SPOTLESS; no. 101140786), from the Spanish grant PID2024-158486OB-C31 funded by MCIU/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, from the programme Unidad de Excelencia María de Maeztu CEX2020-001058-M financed by MCIN/AEI/10.13039/501100011033 and by the MaX-CSIC Excellence Award MaX4-SOMMA-ICE, and from the Generalitat de Catalunya/CERCA programme. J.B.-P. also acknowledges financial support from Spanish grant PID2020-120375GB-I00, funded by MCIU/AEI, Consolidación 2022 CNS2022-136050 and PRE2022-101942 funded by MCIU/AEI/10.13039/501100011033 and ESF+.

Light Curve Fitting

We model the light curves as the sum of a planetary component and a systematic component, following the general form of the models constructed by Coulombe et al. [14] and Splinter et al. [63]. The planetary component contains the transit and eclipse models along with phase curve components for the planetary reflection and emission. WASP-121 also requires model components for stellar variability and ellipsoidal variations of the planetary radius, an effect which is observable for gas giants but not for smaller (more rigid) planets. The systematic portion of the model accounts for visit-long changes to the baseline and “tilt” offsets caused by slight movements of JWST’s primary mirror [58, 4]. The LTT-9779 light curve has two of these events, while WASP-121 has a single event.

Taking all of this into account, we assemble the following model for both targets’ light curves:

where $F(t)$ is the total normalized flux, $A(t)$ is the astrophysical component, and $S(t)$ is the systematics model. The specifics of $A(t)$ and $S(t)$ vary between targets, so we start with LTT-9779:

where $F_\ast(t)$ and $F_p(t)$ are the stellar and planetary flux models. For LTT-9779 the $>45$-day rotation period of this star is significantly longer than the $\sim1$-day observational window, so we set $F_\ast(t) = 1.0$ and model $F_p(t)$ as:

where $\tau(t)$ is the transit model, $\epsilon(t)$ is the eclipse model, $f_p$ is the average depth between the two eclipses, and $\Phi(t)$ is the thermal phase curve. We model $\Phi(t)$ as a 2nd order (i.e., up to $n = 2$) Fourier series:

where $C_n$ and $D_n$ are coefficients for $n = 1$ and $n = 2$ (4 coefficients total), and $\theta(t)$ is the orbital phase, defined as:

with eclipse midpoint $T_{0,\mathrm{ecl}}$ and orbital period calculated from the difference in eclipse midpoints, $P_{\mathrm{ecl}}$.

For the systematics model:

where $B(t)$ is the baseline and $J_{1,2}$ represents the first and second tilt events. The out-of-transit baseline is modeled with a 3rd order polynomial:

where $a_n$ are free coefficients. The tilt events are modeled as:

where $j_{1,2}$ are the tilt amplitudes, $T_{\mathrm{tilt},1,2}$ represents the time of the tilt events and $H(t - T_{\mathrm{tilt},1,2})$ is the Heaviside function that applies a flux offset to the baseline for all $t > T_{\mathrm{tilt},1,2}$.

For WASP-121 the model looks very similar, with a few key differences. Our astrophysical model must now account for ellipsoidal variations of the planet’s radius (relevant for gas giants like WASP-121 but not smaller more rigid planets like LTT-9779) as well as stellar rotation:

with

with rotation amplitude $A_\ast$, stellar rotation period $P_\ast$ (fixed at 1.13 days; Sing et al. [62]), mid-transit time $t_0$, and stellar phase offset $\phi_\ast$. The planetary contribution is the combined effect of a transit, eclipses, ellipsoidal variation, and thermal phase:

$\Phi(t)$ is a second-order Fourier series and the ellipsoidal variations function $\Omega(t)$ is defined as:

where $i$ is the orbital inclination and $R_{\mathrm{sub}}/R_p$ is the ratio of the planetary radius at the substellar point (the line connecting the star and planet) to the planetary radius along the polar axis.

Finally, the systematics and baseline model for WASP-121 is very similar to that used for LTT-9779:

where $B(t)$ is a linear function with a slope and offset, and $J(t)$ is the model for a single tilt event.

In the transit model, we fix eccentricity to zero and fit for the transit midpoint $T_{0,\mathrm{tra}}$, planet radius $R_p/R_\ast$, scaled semi-major axis $a/R_\ast$, inclination $i$, and the quadratic limb darkening coefficients $u_1$ and $u_2$. The eclipse model adds two additional free parameters per eclipse: the eclipse midpoint $T_{0,\mathrm{ecl}}$ and the depth $f_p$. The minimization is first applied to the BI light curve for both targets, and then achromatic transit parameters are fixed for subsequent fits to the spectroscopic light curves as shown in Figures A.1 and A.2. While this approach to fitting planetary phase curves is relatively low-complexity and inappropriate for inferring physical parameters, it is appropriate for the purpose of removing the dominant planetary signatures to study the stellar noise.

Model Selection

We compare models which have different C (the Harvey Profile slope) and $N_{\mathrm{components}}$. The standard Harvey profile parametrization using $C = 2$ corresponds in the time domain to a simple harmonic oscillator, while $C = 4$ corresponds to a symmetric rising and decaying pulse [50, 37]. Several studies have instead set $C = 4$ or have it as a free parameter (e.g., [50, 36, 35, 27]). In order to determine the best parameterization of $N$ and $C$, we fit a suite of different models to either star’s BI PSD. We vary $N$ between one-component and two-components, and vary $C$ to be fixed at 2, fixed at 4, and in the case of two-components, all permutations of 2 and 4 ({2, 2}, {4, 4}, {2, 4}, and {4, 2}); a total of 8 models. We then compare both the Bayesian Information Criteria (BIC) and Akaike Information Criteria (AIC) for the entire suite of models. We performed a series of tests in order to determine how many components and what values of $C$ are justified by these data. Table B.1 shows the results of this test, from which we conclude that a single component with $C=2$ ($\mathcal{M}_1$) is the best choice for measuring and reporting results for these systems. Comparing all results across targets, the next best model appears to be $N=2$ with $C=\{2, 4\}$ ($\mathcal{M}_3$).

Solar Comparison

These model comparisons are not performed for the TSI data, as the Solar PSD is unambiguously best fit by a single component with $C=2$. Figure B.1 shows that, while a rise in power at lower frequencies is a genuine physical signature in the TSI data, the amplitude and frequency at which the “plateau” occurs is highly dependent on the length of the window from which the power spectrum is calculated. This illustrates the sensitivity of the inferred low-frequency structure to the observing window and limits direct comparisons with the shorter JWST observations. Based on this comparison, we do not expect to see any genuine evidence of a second granulation component on the timescale of these JWST observations. We have computed and shown results of a “second component” for LTT-9779 and WASP-121 but this is not a granulation component – rather, it is likely to be from some combination of stellar activity, residual planetary signatures, instrumental noise, or other sources of excess noise. The solar comparison alone does not establish the physical origin of an additional component in the JWST spectra. We therefore treat the longer-timescale component of $\mathcal{M}3$ as a phenomenological description of excess low-frequency power, without identifying it as a separate granulation signal.

Derivation of $\sigma_{\mathrm{GP}}$

The time series variance of the power spectrum used in celerite is defined as

Excluding negative frequencies and substituting eqn. 8 for $P(\omega)$, we can rewrite this as:

Model	$k$	AIC	$\Delta$AIC	$p(\mathrm{model})_{\mathrm{AIC}}$	BIC	$\Delta$BIC	$p(\mathrm{model})_{\mathrm{BIC}}$
LTT-9779 ($n_{\mathrm{data}}=2016$)
N=1, C=2 (M1)	6	57330.08	0.00	0.5716	57363.74	0.00	0.9747
N=2, C={2, 4} (M3)	8	57332.53	2.45	0.1679	57377.41	13.67	0.0010
N=2, C={2, 2}	8	57332.96	2.88	0.1353	57377.84	14.10	0.0008
N=2, C={4, 2}	8	57334.71	4.63	0.0565	57379.58	15.85	0.0004
N=2, C={4, 4}	8	57334.75	4.67	0.0555	57379.62	15.88	0.0003
N=1, C=4 (M2)	6	57337.60	7.52	0.0133	57371.26	7.52	0.0227
WASP-121 ($n_{\mathrm{data}}=1096$)
N=2, C={2, 4} (M3)	8	32462.48	0.00	0.3410	32502.48	9.07	0.0105
N=2, C={2, 2}	8	32463.34	0.86	0.2219	32503.34	9.93	0.0068
N=1, C=2 (M1)	6	32463.41	0.93	0.2140	32493.41	0.00	0.9745
N=2, C={4, 4}	8	32464.51	2.03	0.1235	32504.51	11.10	0.0038
N=2, C={4, 2}	8	32464.95	2.47	0.0992	32504.95	11.54	0.0030
N=1, C=4 (M2)	6	32476.45	13.97	0.0003	32506.46	13.04	0.0014

Solving this integral gives us:

Now we return to the definition of $a_{\mathrm{gran}}$, the measured parameter reported as granulation amplitude in Coulombe et al. [14] and Splinter et al. [63], which is

Comparing these two equations, we find the following:

This equation tells us how to relate the time series variance to the GP-reported amplitude, $a_{\mathrm{gran,GP}}$: