<figure><img src="image_1.png" /></figure>

# The Radiative Effects of Photochemical Hazes on the Atmospheric Circulation and Phase Curves of Sub-Neptunes

Maria E. Steinrueck<sup>1,2,14</sup> (1, 2), Vivien Parmentier<sup>3,4</sup> (1, 2), Laura Kreidberg<sup>2</sup> (1, 2), Peter Gao<sup>5</sup> (1, 2), Eliza M.-R. Kempton<sup>6</sup> (1, 2),

Michael Zhang<sup>1,14</sup> (1, 2), Kevin B. Stevenson<sup>7</sup> (1, 2), Isaac Malsky<sup>8</sup> (1, 2), Michael T. Roman<sup>9,10</sup> (1, 2), Emily Rauscher<sup>11</sup> (1, 2), Matej Malik<sup>6</sup> (1, 2), Roxana Lupu<sup>12</sup> (1, 2), Tiffany Kataria<sup>8</sup> (1, 2), Anjali A. A. Piette<sup>13</sup> (1, 2), Jacob L. Bean<sup>1</sup> (1, 2), and Matthew C. Nixon<sup>6</sup> (1, 2)

<sup>1</sup>Department of Astronomy & Astrophysics, University of Chicago, Chicago, IL 60637, USA

<sup>2</sup>Max-Planck-Institut für Astronomie, 69117 Heidelberg, Germany

<sup>3</sup>Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France

<sup>4</sup>Atmospheric, Ocean, and Planetary Physics, Department of Physics, Oxford University, OX1 3PU, UK

<sup>5</sup>Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road, NW, Washington, DC 20015, USA

<sup>6</sup>Department of Astronomy, University of Maryland, College Park, MD 20742, USA

<sup>7</sup>JHU Applied Physics Laboratory, 11100 Johns Hopkins Rd, Laurel, MD 20723, USA

<sup>8</sup>Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA

<sup>9</sup>Facultad de Ingeniería y Ciencias, Universidad Adolfo Ibáñez, Santiago, Chile

<sup>10</sup>School of Physics and Astronomy, University of Leicester, Leicester, UK

<sup>11</sup>Department of Astronomy and Astrophysics, University of Michigan, Ann Arbor, MI 48109, USA

<sup>12</sup>Eureka Scientific, Inc., Oakland, CA 94602, USA

<sup>13</sup>School of Physics & Astronomy, University of Birmingham, Edgbaston, Birmingham, UK

Received 2024 November 4; revised 2025 March 28; accepted 2025 April 3; published 2025 May 16

## Abstract

Measuring the atmospheric composition of hazy sub-Neptunes like GJ 1214b through transmission spectroscopy is difficult because of the degeneracy between mean molecular weight (MMW) and haziness. It has been proposed that phase-curve observations can break this degeneracy because of the relationship between MMW and phase-curve amplitude. However, photochemical hazes can strongly affect phase-curve amplitudes as well. We present a large set of general circulation model simulations of the sub-Neptune GJ 1214b that include photochemical hazes with varying atmospheric composition, haze opacity, and haze optical properties. In our simulations, photochemical hazes cause temperature changes of up to 200 K, producing thermal inversions and cooling deeper regions. This results in increased phase-curve amplitudes and adds a considerable scatter to the phase-curve amplitude-metallicity relationship. However, we find that if the haze production rate is high enough to significantly alter the phase curve, the secondary eclipse spectrum will exhibit either emission features or strongly muted absorption features. Thus, the combination of a white-light phase curve and a secondary eclipse spectrum can successfully distinguish between a hazy, lower-MMW and a clear, high-MMW scenario.

*Unified Astronomy Thesaurus concepts: Exoplanet atmospheric dynamics (2307); Exoplanet atmospheres (487); Exoplanet atmospheric structure (2310); Extrasolar gaseous giant planets (509); Extrasolar gaseous planets (2172)*

## 1. Introduction

Sub-Neptunes are among the most common planets found around other stars (A. W. Howard et al. 2012). Given that no such planet exists in our own solar system, planet formation models historically did not account for them. The formation (J. J. Fortney et al. 2013; E. J. Lee & E. Chiang 2016; M. Lambrechts et al. 2019; L. Zeng et al. 2019; R. Burn et al. 2024), evolution (J. E. Owen & Y. Wu 2013; S. Ginzburg et al. 2018; C. Mordasini 2020), and interior structure (e.g., L. A. Rogers & S. Seager 2010; N. Nettelmann et al. 2011; L. A. Rogers et al. 2011; D. Valencia et al. 2013; R. Luque & E. Pallé 2022; M. C. Nixon et al. 2024) of these planets are currently a subject of debate. Measuring the atmospheric composition of sub-Neptunes is expected to provide constraints on the history of sub-Neptunes and to break degeneracies in

interior structure models (e.g., J. L. Bean et al. 2021b; B. Benneke et al. 2024; M. C. Nixon et al. 2024).

Transmission spectroscopy of sub-Neptunes has proven to be challenging owing to the degeneracy between aerosols (photochemical hazes and condensate clouds) and mean molecular weight (MMW), both of which can reduce the amplitude of spectral features (B. Benneke & S. Seager 2012; M. R. Line & V. Parmentier 2016). This degeneracy is particularly limiting for observations taken in a narrow bandpass such as the WFC3 G141 grism of the Hubble Space Telescope (HST; e.g., H. A. Knutson et al. 2014; L. Kreidberg et al. 2014; I. J. M. Crossfield & L. Kreidberg 2017; J. Brande et al. 2024). The wider wavelength coverage and higher precision of JWST are able to mitigate this degeneracy in some cases, due to covering multiple absorption bands of the same species and resolving subtle differences in spectral shape between those two scenarios (e.g., B. Benneke & S. Seager 2012, 2013; C. Mai & M. R. Line 2019; C. Piaulet-Ghorayeb et al. 2024). However, aerosols remain a challenge even with JWST (e.g., N. L. Wallack et al. 2024). Furthermore, the benchmark sub-Neptune GJ 1214b, which has the highest transmission spectra metric (E. M. R. Kempton et al. 2018) of all known sub-Neptunes and thus should be the easiest target

<figure><img src="image_2.png" /></figure>

<footnote n="14"><sup>14</sup> 51 Pegasi b fellow.</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
for atmospheric characterization, has an unexpectedly flat transmission spectrum (L. Kreidberg et al. 2014; P. Gao et al. 2023; K. Ohno et al. 2024; E. Schlawin et al. 2024). The extremely low feature amplitude of its spectrum cannot be explained with a high MMW alone and therefore requires high-altitude aerosols.

Most sub-Neptunes that are promising targets for atmospheric characterization are in the temperature range of $300 \text{ K} < T_{\text{eq}} < 900 \text{ K}$. For these temperatures, models predict that photochemical hazes are likely to form (C. V. Morley et al. 2013, 2015; Y. Kawashima & M. Ikoma 2018, 2019; D. Adams et al. 2019; P. Lavvas et al. 2019) and to be the most important opacity source in the atmosphere (P. Gao et al. 2020). Laboratory experiments also suggest that photochemical hazes readily form across a range of atmospheric compositions and temperatures relevant to sub-Neptunes (C. He et al. 2018, 2020; S. M. Hörst et al. 2018). For GJ 1214b in particular, general circulation models (GCMs; B. Charnay et al. 2015a) and microphysics models (K. Ohno & S. Okuzumi 2018; K. Ohno et al. 2020) suggest that it may be hard to lift condensate clouds to high-enough altitudes to explain the flatness of the spectrum. Thus, photochemical hazes have become a favored explanation for the featureless atmosphere of GJ 1214b in particular (e.g., P. Gao et al. 2023).

Phase curves and emission spectra are alternative ways of characterizing the atmospheres of sub-Neptunes that could be used to break the degeneracies from transmission spectroscopy. In particular, phase curves have been proposed as an alternative measurement of the MMW of sub-Neptune atmospheres (J. L. Bean et al. 2021a). This is because atmospheric dynamics predicts a relationship between phase-curve amplitude and MMW, with an increase in the phase-curve amplitude with increasing MMW (T. Kataria et al. 2014; X. Zhang & A. P. Showman 2017). E. M. R. Kempton et al. (2023) recently applied this approach to the benchmark sub-Neptune GJ 1214b, utilizing JWST MIRI LRS to observe the first thermal phase curve of a sub-Neptune. Emission spectra are also of interest for planets with strong aerosol coverage because they probe deeper into the atmosphere than transmission spectra (J. J. Fortney 2005).

However, 1D models show that photochemical hazes can strongly impact the temperature structure (C. V. Morley et al. 2015; A. A. A. Piette & N. Madhusudhan 2020; P. Lavvas & A. Arfaux 2021). The absorption of stellar photons by hazes results in a thermal inversion at low pressures and a cooling of the deep atmosphere. 3D simulations of hot Jupiters show that for sufficiently high haze production rates hazes can strongly affect atmospheric circulation, increasing or decreasing the strength of the equatorial jet depending on the assumed haze refractive index, and increase the phase-curve amplitude (M. E. Steinrueck et al. 2023). It thus seems likely that the MMW-phase-curve amplitude relation for sub-Neptunes, which so far has been mainly studied under the assumption of clear atmospheres (with one study including radiative feedback of condensate clouds; B. Charnay et al. 2015b), will be affected by photochemical hazes. However, the radiative feedback of hazes in GCMs of sub-Neptunes has not been studied systematically.

While condensate clouds in sub-Neptune atmospheres have been studied using GCMs with multiple approaches (B. Charnay et al. 2015a, 2015b; D. A. Christie et al. 2022), photochemical hazes have not been included in sub-Neptune

GCMs until recently. As part of the effort to interpret the phase-curve observations of GJ 1214b, we included a large grid of GCM simulations of that planet with horizontally homogeneous hazes in the observational paper (E. M. R. Kempton et al. 2023). These were the first published sub-Neptune GCMs that considered the radiative effects of a haze layer. The models showed that hazes can have a first-order effect on the phase curve and that both a high metallicity and highly scattering aerosols were required to match the observed phase curve. In this publication, we are conducting an in-depth analysis of this large and rich set of simulations of GJ 1214b to better understand how hazes affect the atmospheres of sub-Neptunes and what conclusions can be drawn for the broader sub-Neptune population.

The remainder of this paper is structured as follows: In Section 2, we describe our modeling approach. Section 3 describes temperature structure and atmospheric circulation of the simulations, while Section 4 focuses on observable trends in the 5–12 µm region. Further, we briefly discuss how our models compare to the recent JWST MIRI LRS observations of GJ 1214b, including the original (E. M. R. Kempton et al. 2023) data reduction and an alternative data reduction presented in a second paper from our collaboration (I. Malsky et al. 2025), in Section 5. Finally, we discuss possible observing strategies for future observations of sub-Neptunes in Section 6 and summarize our conclusions in Section 7.

## 2. Methods

We use SPARC/MITgcm (A. P. Showman et al. 2009; T. Kataria et al. 2013), which couples the dynamical core of A. Adcroft et al. (2004) to the wavelength-dependent radiative transfer code of M. S. Marley & C. P. McKay (1999), to simulate the atmosphere of sub-Neptune GJ 1214b. SPARC/MITgcm has in the past mainly been used to simulate the atmospheres of hot Jupiters (e.g., V. Parmentier et al. 2013; T. Kataria et al. 2016; M. E. Steinrueck et al. 2019, 2023; S.-M. Tsai et al. 2023), but it has also been successfully applied to hot Neptunes (N. K. Lewis et al. 2010) and sub-Neptunes (T. Kataria et al. 2014).

### 2.1. Atmospheric Dynamics

We use the dynamical core of A. Adcroft et al. (2004) to solve the primitive equations on a cubed-sphere grid. Our numerical setup uses a fourth-order Shapiro filter (R. Shapiro 1970) to damp numerical fluctuations at the grid scale that could otherwise lead to instabilities. Further, we include a pressure-dependent linear drag, $(dv/dt)_{\text{drag}} = -k_v v$, in the deep atmosphere ($p > p_{\text{drag,top}} = 10$ bars), which ensures convergence independent of initial conditions given long-enough simulation run times (B. Liu & A. P. Showman 2013) and stabilizes the simulation. The drag strength is given by $k_v = k_F(p - p_{\text{drag,top}})/(p_{\text{bottom}} - p_{\text{drag,top}})$, with the bottom boundary of the simulation $p_{\text{bottom}} = 200$ bars and the maximum drag strength $k_F = 10^{-5} \text{ s}^{-1}$. The most important model parameters are summarized in Table 1 for parameters shared by all simulations and in Table 2 for parameters depending on the assumed composition of the atmosphere.

### 2.2. Radiative Transfer

The radiative transfer code used in our simulations (M. S. Marley & C. P. McKay 1999) assumes a plane-parallel
<table><caption>**Table 1**
Model Parameters Shared by All Simulations</caption><thead><tr><th>Parameter</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Radius<sup>a</sup></td><td>1.7469282 × 10<sup>7</sup></td><td>m</td></tr><tr><td>Gravity<sup>a</sup></td><td>10.65</td><td>m s<sup>-2</sup></td></tr><tr><td>Rotation period<sup>a,b</sup></td><td>1.58</td><td>days</td></tr><tr><td>Semimajor axis<sup>a</sup></td><td>0.0149</td><td>au</td></tr><tr><td>Interior flux<sup>c</sup></td><td>4.59 × 10<sup>-2</sup></td><td>W m<sup>-2</sup></td></tr><tr><td>Horizontal resolution</td><td>C32<sup>e</sup></td><td></td></tr><tr><td>Vertical resolution<sup>d</sup></td><td>60</td><td>layers</td></tr><tr><td>Lower pressure boundary<sup>d</sup></td><td>1.75 × 10<sup>-7</sup></td><td>bars</td></tr><tr><td>Upper pressure boundary</td><td>200</td><td>bars</td></tr></tbody></table>

**Notes.**

<sup>a</sup> R. Cloutier et al. (2021).

<sup>b</sup> Assumed to be equal to the orbital period due to tidal locking.

<sup>c</sup> Corresponding to an intrinsic temperature of 30 K, based on E. D. Lopez & J. J. Fortney (2014).

<sup>d</sup> With the exception of the 3000× solar clear-atmosphere simulation, which has 53 layers and a lower pressure boundary of 2 × 10<sup>-6</sup> bars.

<sup>e</sup> Equivalent to a resolution of 128 × 64 on a longitude–latitude grid.

atmosphere and uses the two-stream approximation. Wave-length dependence is taken into account with the correlated-k method with 11 wavelength bins (T. Kataria et al. 2013) within the simulations and 196 wavelength bins for post-processing simulations to obtain spectra. We used molecular opacities based on R. S. Freedman et al. (2008, 2014) and R. E. Lupu et al. (2014), with the more recent updates described in M. S. Marley et al. (2021). For the opacity calculation, atmospheric abundances at thermochemical equilibrium were calculated over a grid of temperatures (75–6000 K) and pressures (10<sup>-6</sup>–3000 bars) based on the models in E. Gharib-Nezhad et al. (2021) and the Sonora atmosphere model (M. S. Marley et al. 2021), extended here to metallicities up to 3000× solar. While the atmosphere of GJ 1214b is not expected to be in chemical equilibrium, studies show that for planets in the temperature range of GJ 1214b, at least for hydrogen-dominated atmospheres, the impact of transport-induced disequilibrium chemistry on the temperature structure is smaller than the possible effects of hazes (M. E. Steinrueck et al. 2019; M. Zamyatina et al. 2023). Furthermore, the current data have so far been insensitive to disequilibrium species.

For the incoming starlight, we used a stellar spectrum based on a PHOENIX model (T. O. Husser et al. 2013) using the stellar parameters of GJ 1214b from R. Cloutier et al. (2021). At the bottom boundary of the model, we include an upward radiative flux corresponding to an intrinsic temperature of 30 K (E. D. Lopez & J. J. Fortney 2014).

The white-light phase curves were calculated by integrating the planetary flux weighted by the MIRI filter response at each phase and dividing it by the similarly integrated stellar flux, following Equation (1) of D. Charbonneau et al. (2005). Details about our post-processing procedure can be found in V. Parmentier et al. (2016).

## 2.3. Haze Model

We assume a haze profile that is static and horizontally uniform on isobars. For a first exploration of the radiative effects of photochemical hazes in sub-Neptunes, this is justified because a tracer-based model with radiatively active hazes as in M. E. Steinrueck et al. (2023) is much more computationally

expensive and would not allow an as-thorough exploration of the parameter space. Furthermore, for GJ 1214b-like conditions, haze particle radii are expected to vertically vary by more than two orders of magnitude (Y. Kawashima & M. Ikoma 2019; P. Lavvas et al. 2019). It is thus likely that the vertical variations in particle radius, and thus opacity, are larger than the horizontal ones.

The nominal haze profile was derived from a 1D simulation with CARMA (O. B. Toon et al. 1979; R. P. Turco et al. 1979; A. S. Ackerman et al. 1995), a bin-based microphysical model for clouds and hazes. In particular, the simulation setup of D. Adams et al. (2019) with a haze production rate of 10<sup>-12</sup> g cm<sup>-2</sup> s<sup>-1</sup>, a metallicity of 100× solar, and spherical hazes was used. This haze production rate is in the typical range of predicted haze formation rates for GJ 1214b based on photolysis mass fluxes of haze precursor molecules in photochemical models (Y. Kawashima & M. Ikoma 2019; P. Lavvas et al. 2019). We choose the 100× solar case because it is in the intermediate range of metallicities that we consider. Throughout this work, we use the same haze profile independent of metallicity. We do not scale haze opacity as a function of metallicity because neither experimental nor theoretical studies show evidence for a monotonous relationship between haze production rate and metallicity (C. He et al. 2018; S. M. Hörst et al. 2018; P. Lavvas et al. 2019). Number density and mean radius from the microphysics model are shown in the left panel of Figure 1.

We then calculated haze optical depth, single-scattering albedo, and asymmetry parameter as a function of pressure for the particle size distribution from the microphysics model with a Mie scattering code (C. F. Bohren & D. R. Huffman 1983), assuming the same pressure grid as used in the microphysics model. Because the haze opacity varies weakly within each bin, we consider it as a continuum opacity that can just be added to the k-coefficients within each bin. We took the wavelength average of each of these three quantities within each of the 11 wavelength bins used in the GCM and interpolated the profiles onto the vertical grid of the GCM. For the few GCM layers above the low-pressure boundary of the microphysics model (1 µbar), quantities were extrapolated linearly in pressure.

To simulate different haze production rates, we scaled the optical depth profile by a factor of 10 or 100, while leaving the single-scattering albedo and asymmetry parameter unchanged. This approach neglects that particle sizes can vary between different metallicities and haze production rates. However, it provides a simple and fast way of exploring parameter space. It also makes it easier to isolate changes in atmospheric dynamics and radiative transfer due to the gas-phase composition. Multiple studies (D. Adams et al. 2019; Y. Kawashima & M. Ikoma 2019; P. Lavvas et al. 2019; P. Gao et al. 2023) have explored the effects of varying metallicity and haze production rate on particle sizes in 1D. Our study aims to complement this work by focusing on changes in atmospheric dynamics while leaving the particle size distribution constant.

For the optical properties of the photochemical hazes, we explored three different scenarios. In the first scenario, we assumed refractive indices of soot based on the compilation of P. Lavvas & T. Koskinen (2017) in the Mie calculation. This is a frequently used assumption in simulations of hot exoplanet atmospheres (e.g., C. V. Morley et al. 2015; D. Adams et al. 2019; P. Lavvas et al. 2019; M. E. Steinrueck et al. 2021), as soots can survive much higher temperatures than most other
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Summary of the haze properties assumed in the simulations. The left panel shows number density and mean radius of the haze distribution that was used as input in the GCM simulations. The middle panel shows the differential optical depth for soot (solid) and tholin (dashed) hazes at 1 µm (representative of the incoming starlight) and 8 µm (representative of the outgoing thermal radiation) for the nominal haze opacity case. The pressure levels where the optical depth reaches 0.01 (circle), 0.1 (star), and 1 (square) are highlighted. For soot at 1 µm, the optical depth reaches 0.01 at a pressure of 6 × 10⁻⁷ bars, outside the pressure range shown. The right panel shows the wavelength-averaged single-scattering albedo weighted by the stellar flux. For the “maximally reflective” case, the optical depth is identical to the soot case, but the single-scattering albedo is fixed to 0.9999.</figcaption></figure>

<table><caption>**Table 2**
Composition-dependent Model Parameters</caption><thead><tr><th>Composition</th><th>Mean Molecular Weight (g mol<sup>-1</sup>)</th><th>Specific Heat Capacity (J kg<sup>-1</sup> K<sup>-1</sup>)</th><th>Specific Gas Constant (J kg<sup>-1</sup> K<sup>-1</sup>)</th><th>Hydrodynamic Time Step (s)</th><th>Radiative Time Step (s)</th></tr></thead><tbody><tr><td>1 × solar</td><td>2.24</td><td>13,000</td><td>3714</td><td>25</td><td>50</td></tr><tr><td>100 × solar</td><td>4.38</td><td>6474</td><td>1898</td><td>25</td><td>50</td></tr><tr><td>300 × solar</td><td>8.50</td><td>3500</td><td>978</td><td>10</td><td>20</td></tr><tr><td>3000 × solar</td><td>28.5</td><td>1300</td><td>292</td><td>10</td><td>20</td></tr></tbody></table>

candidate haze compositions. The second scenario assumes refractive indices of Titan-like tholins (B. N. Khare et al. 1984) as an example of more scattering hazes with a stronger wavelength dependence. As a third scenario, motivated by the observed high albedo of GJ 1214b, we also consider hazes with the extinction cross section and asymmetry parameter equal to the soot scenario but with the single-scattering albedo set to 0.9999. We dub this scenario “maximally reflective hazes.” These assumptions are identical to what was assumed for “maximally reflective hazes” in E. M. R. Kempton et al. (2023). We note that compared to scattering and absorption calculated self-consistently from a realistic refractive index, this approach may overestimate extinction at long wavelengths, due to the absorption cross section dropping off ∝λ⁻¹ while the scattering cross section is ∝λ⁻⁴ in the Rayleigh scattering regime. However, our goal here is to provide an extreme limiting case for illustration, and this caveat is unlikely to qualitatively change the conclusions.

The optical depth profiles for two representative wavelengths and the single-scattering albedo profile averaged across the full wavelength range, weighed by the stellar spectrum, are shown in the middle and right panels of Figure 1. Furthermore, Figure 18 illustrates the wavelength dependence of these properties at a range of pressures.

Recently, additional optical properties from experiments specifically aiming to simulate haze formation on exoplanets were published (L. Corrales et al. 2023; C. He et al. 2024). These sets of optical constants still are qualitatively similar to the B. N. Khare et al. (1984) tholins in that they feature a strong UV opacity slope, an opacity window region with a very low extinction coefficient in the optical to near-infrared, moderate extinction coefficients at near-infrared wavelengths longer than the opacity window, and somewhat higher extinction coefficients in the mid-infrared. It thus seems likely that these optical property sets will produce qualitatively similar results to the tholin case, and we thus do not run separate simulations for them.

## 2.4. Simulations

Simulations were initiated from a state of rest, with a globally uniform initial pressure-temperature profile derived from a 1D radiative transfer calculation with HELIOS (M. Malik et al. 2017, 2019). All simulations were run for 1000 Earth days. The results presented in this paper are based on the time average of the last 100 simulation days. We ran a total of 27 simulations. This includes all 24 simulations listed in Extend Data Table 3 in E. M. R. Kempton et al. (2023), as
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Dayside-averaged (solid) and nightside-averaged (dashed) temperature profiles for all clear-atmosphere simulations. Temperatures were weighted by the cosine of the angle of incidence (dayside)/the angle from the antistellar direction (nightside) during averaging.</figcaption></figure>

well as three additional simulations with tholins (300× solar, haze scale factor 100×; 3000× solar, haze scale factor 10×; 3000× solar, haze scale factor 100×).

# 3. Results

## 3.1. Temperature Structure

For clear-atmosphere models there is a well-established trend of less efficient heat redistribution and thus higher day-to-night temperature contrasts with higher MMW (T. Kataria et al. 2014; B. Charnay et al. 2015a; X. Zhang & A. P. Showman 2017). This behavior is once again borne out by our clear-atmosphere models (Figure 2). Introducing hazes can dramatically change the temperature structure of the atmosphere. The magnitude of these changes increases with enhanced haze opacity (Figure 3). If the haze opacity increases sufficiently, all models form a dayside thermal inversion and exhibit a colder temperature profile in the deep atmosphere. However, the shape of the thermal inversion differs substantially depending on the assumed haze optical properties.

Soot hazes, being highly absorbing across the whole wavelength spectrum, lead to strong heating at low pressures, increasing dayside temperatures by up to 200 K at pressures ≈1 mbar. At the same time, the atmosphere cools compared to the clear-atmosphere case at pressures ≳1–5 mbar owing to the antigreenhouse effect (e.g., G. E. Thomas & K. Stamnes 2002, p. 449). The strength of the thermal inversion strongly depends on the haze opacity: for the nominal haze opacity, the thermal inversion remains moderate, with a temperature increase of ≈40 K between 2 and 0.3 mbar for the 300× solar case, with a large nearly isothermal region of the atmosphere above the inversion. For the 10× increased haze opacity, however, the thermal inversion spans over three orders of magnitude in pressure (between 20 mbar and 10 µbar for the 300× solar case), with a temperature increase of over 300 K. This agrees well with the findings of C. V. Morley et al. (2015).

In the deep atmosphere ($p \gtrsim 0.1$ bars), where atmospheric heat transport is efficient, the nightside temperature profiles closely follow the dayside ones. Thus, nightside temperatures at high and intermediate pressures are cooler for higher haze opacities. This cooling effect is strongest for the soot hazes compared to other haze types. At very low pressures ($p \lesssim 0.1$ mbar), the nightside temperature profiles in the soot simulations become isothermal and are warmer than in the clear-atmosphere profile, which continues to decrease with height.

For the maximally reflective hazes, there is practically no absorption. Thus, the main way in which these hazes affect the temperature structure is by scattering incoming starlight at low pressures, preventing it from reaching deeper layers. The resulting dayside temperature profile is close to the haze-free temperature profile for the nominal haze opacity and exhibits a moderate to large thermal inversion for larger haze opacities. Due to the large amount of light reflected at low pressures, temperature profiles with maximally reflective hazes have cooler daysides at almost all pressures. However, compared to the soot case, the deeper region below the isothermal and inverted-temperature regions cools less.

Nightside temperatures in the deep atmosphere are again mainly set by the amount of dayside cooling. At low pressures ($p \lesssim 0.1$ mbar), the temperature profile becomes warmer and more isothermal with increasing haze opacity. This is likely due to a combination of radiative effects—due to the added haze opacity, radiation can escape less easily—and advection of heat from the dayside due to the increased strength of the equatorial jet at high haze opacities (see Section 3.2).

With a single-scattering albedo much higher than that of soot but less than the maximally reflective case, the temperature changes in the tholin simulations exhibit some characteristics from each of these two extreme cases. The amount of cooling of deeper layers is comparable to or slightly higher than the maximally reflective case. Given the high single-scattering albedo (reaching up to 0.95 between 0.1 and 1 mbar; see Figure 1), this appears to be mainly driven by reflection of incoming starlight. However, there are also small to moderate temperature increases at low pressures in all tholin simulations, due to absorption by tholin haze. In most cases, the additional heating is not enough to change the structure of the atmosphere at low pressures (i.e., in the region above the cooling-driven thermal inversion, temperature still decreases with height). Only in simulations with 100× increased haze opacity does a more extended thermal inversion due to heating and a temperature peak at low pressures form that is comparable to the temperature structure in the soot simulations. We further note that the temperature minimum due to cooling is located at higher pressures than in the maximally reflective case. This is likely caused by the reduced extinction cross section of tholins compared to soot and our maximally reflective hazes (for which the extinction cross section of soot was assumed) at low pressures (see Figures 1 and 18).

## 3.2. Atmospheric Circulation

The atmospheric circulation in the clear-atmosphere models is consistent with previous work and depends strongly on metallicity (Figure 4). In the solar-metallicity case, an equatorial jet and two midlatitude jets (one per hemisphere) form, similar to previous simulations of GJ 1214b (T. Kataria et al. 2014; B. Charnay et al. 2015a) and GJ 436b (N. K. Lewis
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Same as Figure 2, but for a selection of simulations with photochemical hazes compared to the clear simulation with the same metallicity.</figcaption></figure>

et al. 2010). The midlatitude jets are broader and peak at lower pressures than the equatorial jet. The jets are not well separated from each other—while there are distinct peaks, there are no regions of low zonal-mean zonal velocity or westward flow in between those peaks. In the 100× solar simulation, a single broad equatorial jet forms, resembling the results from B. Charnay et al. (2015a). The zonal velocity pattern in the 300× solar case looks similar to the 100× solar case, but the wind speeds are lower. This is consistent with the finding of X. Zhang & A. P. Showman (2017) that a higher MMW leads to lower wind speeds. The highest-metallicity simulation, 3000× solar, features three well-separated jets (one equatorial jet and one high-latitude jet in each hemisphere). X. Zhang & A. P. Showman (2017) suggest that a high MMW leads to narrower jets (and thus a total higher number of jets on the planet for sufficiently high MMW), due to the smaller Rossby

deformation radius. In their CO₂ and N₂ simulations using Newtonian cooling, high-latitude jets can be seen at similar locations to those in our 3000× solar simulation. However, the relative strength of their high-latitude jets compared to their equatorial jet is low (about 25%), while in our simulations they reach more than 50% of the equatorial jet speed, suggesting that full wavelength-dependent radiative transfer and detailed atmospheric composition matter. Qualitatively, the zonal-mean circulation is also in between the H₂O-rich and the 50% H₂O–50% CO₂ cases from T. Kataria et al. (2014) and displays similarity to the H₂O simulation of B. Charnay et al. (2015a).

For simulations including photochemical hazes, the atmospheric circulation transitions to a single superrotating equatorial jet if the haze opacity is sufficiently high. Figure 5 illustrates this transition for the 300× solar case for all three different optical properties. A similar behavior can be observed
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Zonal-mean zonal velocity for the clear-atmosphere simulations. The contours highlight where the zonal-mean zonal velocity reaches 50% and 75% of the peak velocity within the simulation.</figcaption></figure>

for most other metallicities.<fnref n="15" /> In most cases, the jet is relatively narrow, confined to ±20° of latitude. The jet speed typically increases with increasing haze opacity. For soot hazes, the transition to a narrow equatorial jet happens at a lower haze opacity than for tholins and maximally reflective hazes, likely because of their stronger impact on temperature structure and heating rates.

The mechanism for the transition to a more focused and faster jet likely differs between the soot and the maximally reflective hazes. For the soot hazes, it seems likely that strongly enhanced thermal forcing in the regions of jet formation causes the transition. For the maximally reflective hazes, however, the net heating decreases owing to the large fraction of radiation being reflected away from the planet. A different mechanism thus must cause the jet to accelerate and narrow. One possibility is that, due to the cooling of the atmosphere, the Rossby deformation radius shrinks and causes the jet to become confined to a smaller latitude range. The acceleration then might happen because there is a narrower latitude range at which the jet can transport heat from the dayside to the nightside. For the tholin case, both mechanisms could be at work. A detailed investigation of the mechanisms driving the atmospheric circulation changes, however, is outside the scope of our work.

The increased jet speed also leads to more efficient heat advection from the dayside to the nightside at low pressures. This is one of the factors leading to an increase in nightside temperatures compared to the haze-free case at very low pressures ($p \lesssim 10^{-4}$ bars; see Figure 3). In the case of maximally reflective cases, the advection of heat leads to a somewhat decreased day-night temperature contrast at low pressures. In the soot case, however, the day-night temperature contrast at low pressures increases substantially because the heating on the dayside greatly outweighs the increased temperature advection. For tholins, the day-night temperature difference decreases slightly in some simulations, while in other simulations it increases at some pressures and decreases at other pressures.

For higher metallicities, higher haze opacities are required for the atmospheric circulation to transition to a single jet. For

example, for the 100× and 300× solar tholin and maximally reflective simulations, the single narrow jet already occurs at 10× the nominal haze opacity. For the 3000× solar case, the corresponding simulations feature an equatorial jet and two high-latitude jets (Figure 6). Only the simulations with 100× the nominal haze opacity have a single equatorial jet for this metallicity. One possible explanation is that atmospheres with higher metallicity have a larger abundance of strongly absorbing gas species such as CO₂, H₂O, and CO. Thus, the haze opacity makes up a lower fraction of the total opacity at the same haze opacity.

# 4. Observable Trends in the 5–12 µm Region

Now we turn our attention to how the changes in temperature structure and atmospheric circulation described above translate to observable changes to phase curves and emission spectra. We focus on examining the 5–12 µm region covered by JWST MIRI LRS for our examination of observable trends because it is close to the peak of thermal emission of GJ 1214b and other planets in the temperature region where hazes are expected to be most dominant for sub-Neptunes. Furthermore, E. M. R. Kempton et al. (2023) estimated that 50%–60% of the thermal flux of this planet is emitted within the MIRI LRS bandpass.

## 4.1. Phase Curves

We start by focusing on properties that can be extracted directly from white-light phase-curve observations. We first turn our attention to how photochemical hazes affect the relationship between phase-curve amplitude and MMW (or, for our parameterization of atmospheric composition, equivalently, metallicity; X. Zhang & A. P. Showman 2017). The white-light phase-curve amplitudes in the MIRI LRS bandpass are shown in Figure 7. We first point out that the clear-atmosphere simulations reproduce the expected trend of increasing phase-curve amplitude with increasing metallicity. When hazes are added, the phase-curve amplitude generally increases.<fnref n="16" /> For any given combination of haze opacity and haze optical properties, the phase-curve amplitude also increases monotonously with

<footnote n="15">¹⁵ The solar-metallicity case shows a somewhat more complex behavior; however, the focus of our work is on higher metallicities.</footnote>

<footnote n="16">¹⁶ The solar-metallicity simulations with maximally reflective hazes are an exception to this: In this particular case, the phase-curve amplitude drops ever so slightly.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Similar to Figure 4, but for all hazy simulations with a metallicity of 300× solar. The color scale is identical to Figure 4.</figcaption></figure>

increasing metallicity. However, looking at the full set of simulations with varying haze opacities and haze properties, hazes add a considerable spread to the trend, with the difference between the clear-atmosphere amplitude and the amplitude of highest haze opacities with the same metallicity reaching up to 0.4. This spread is enough that it is no longer possible to uniquely determine the metallicity or MMW from the phase-curve amplitude alone. However, the most extreme cases (solar vs. 3000× solar) are still clearly separated.

When looking at the relationship between phase-curve amplitude and haze opacity, it becomes obvious that for low and nominal haze opacities (0.1× and 1× the nominal haze opacity) the phase-curve amplitude typically remains unchanged. When increasing the haze opacity enough to substantially affect the dayside thermal structure, the phase-curve amplitude increases. For tholins and maximally reflective

hazes, this transition happens between 1× and 10× our nominal haze opacity. For soots, however, even the nominal haze opacity is enough to substantially increase the phase-curve amplitude. This is consistent with Section 3.1, where we found that for soot a thermal inversion already forms for the nominal haze opacity.

One might hope that the degeneracy between phase-curve amplitude, haze opacity, and haze optical properties could be resolved by considering additional characteristics of the white-light phase curve such as secondary eclipse depth and the offset of the phase-curve maximum with respect to secondary eclipse. However, these additional observables cannot fully break the degeneracy. This is demonstrated in Figure 8, which shows all simulations, but information on haziness has been removed. There is a large scatter in secondary eclipse depth in addition to the scatter in amplitude. Some simulations with different
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Similar to Figure 4, but for hazy simulations with a metallicity of 3000× solar. The color scale is identical to Figure 4.</figcaption></figure>

metallicities are seen very close to each other in the amplitude–eclipse depth and amplitude–offset diagrams. Depending on the haze optical properties, the secondary eclipse depth can either increase for primarily absorbing hazes like soot or decrease for primarily reflecting hazes like tholins or our “maximally reflective” hazes (Figure 9). In the case of reflecting hazes, it is thus possible to distinguish between the case of a clear atmosphere with higher metallicity and a hazy atmosphere with a lower metallicity with the same phase-curve amplitude. Different metallicities with different levels of haze opacity, however, can still look very similar to each other. However, for soot hazes, the secondary eclipse depth increases with increasing haze opacity, as it would for a clear atmosphere with increasing metallicity. It is thus impossible to tell the

difference between these two scenarios based on the white-light phase-curve amplitude and secondary eclipse depth alone.

Further, we note that the phase-curve offset (Figure 10) is unlikely to provide useful information on metallicity or haze opacities in most cases. The phase-curve offset is strongly correlated with the phase-curve amplitude (Figure 8, right panel) and thus does not provide any complementary information to the phase-curve amplitude. In addition, the offset is already small (≈25°) in the 100× solar clear-atmosphere simulation. While it further decreases to ≈15° at 3000× solar (clear atmosphere), the difference of 10° between these two models is comparable to typical observational uncertainties of the phase-curve offset. For example, the 1σ uncertainty on the GJ 1214b phase curve in the Malsky et al.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** White-light phase-curve amplitudes in the JWST MIRI LRS bandpass of all simulations included in this study, plotted against metallicity (left panel) and haze opacity (right panel). The phase-curve amplitudes generally increase with increasing metallicity and with increasing haze opacity. We define the phase-curve amplitude as $(F_{max} - F_{min})/F_{max}$, where $F_{max}$ and $F_{min}$ are the maximum and minimum values of the flux of the planet. The gray horizontal lines and shaded regions indicate the observed values from the Kempton et al. and Malsky et al. data reductions, including their $1\sigma$ uncertainty intervals.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** White-light phase-curve amplitudes in the JWST MIRI LRS bandpass of all simulations included in this study, plotted against secondary eclipse depth (left panel) and eastward phase-curve offsets (right panel). Information on hazes has been removed to demonstrate that differences in haze opacity and haze optical properties lead to a spread in phase-curve amplitudes, secondary eclipse depth, and phase-curve offset that make it impossible to determine metallicity from the white-light phase curve alone. The color scale represents metallicity.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** White-light secondary eclipse depths in the JWST MIRI LRS bandpass of all simulations included in this study, plotted against metallicity (left panel) and haze opacity (right panel). Generally, absorbing hazes like soots increase the eclipse depth, while more scattering hazes like tholins and maximally reflective hazes decrease the eclipse depth.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Eastward white-light phase-curve offsets in the JWST MIRI LRS bandpass of all simulations included in this study, plotted against metallicity (left panel) and haze opacity (right panel).</figcaption></figure>

reduction is ±4°. Furthermore, observational publications often underestimate the uncertainty on phase-curve offsets because they do not take pipeline-induced uncertainties into

account (T. J. Bell et al. 2021). Thus, even in the clear-atmosphere case, it is already hard to distinguish between different metallicities based on the phase-curve offset. The
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Dayside and nightside emission spectra of the haze-free simulations.</figcaption></figure>

additional spread introduced by the different levels of haze opacities further ruins the already slim prospects of determining metallicity based on the phase offset.

All solar-metallicity simulations, clear and hazy, exhibit a much larger eastward phase-curve offset than all other simulations. All of the solar-metallicity hazy models further can be clearly distinguished from the clear-atmosphere case by having an even larger phase-curve offset. The larger phase-curve offset is driven by a stronger eastward jet. In the case of the solar metallicity, the jet speed is very sensitive to even the nominal haze opacity, and the resulting increase in the hot spot offset is larger than the effect of probing higher layers (with a smaller hot spot offset), due to the added haze opacity. However, we caution that all of these models also exhibit a very low phase-curve amplitude. For low-amplitude phase curves, the offset generally is harder to measure observationally and the observational uncertainty will be larger. Past studies have also shown that the phase-curve offset can be strongly

affected by inhomogeneous cloud coverage (V. Parmentier et al. 2021), an effect not included in our models.

We thus conclude that a band-integrated phase curve alone may be insufficient to determine the metallicity or MMW of the atmosphere in the presence of hazes. Next, we turn to emission spectra. We will further discuss potential observational strategies for determining atmospheric composition for hazy atmospheres in Section 6.

## 4.2. Emission Spectra

Before we examine the effect of photochemical hazes on dayside and nightside spectra, we provide a brief overview of the clear-atmosphere emission spectra (Figure 11). Contribution functions for the clear-atmosphere cases are shown in Figure 19 in Appendix B. Water and methane are responsible for the most prominent absorption feature visible at all metallicities, spanning from approximately 5 to 9 µm. Water here dominates the spectrum below ≈7 µm, and methane dominates above ≈7 µm, with both absorption bands smoothly blending into each other. With increasing metallicity, CO₂ gains increasing importance (see, e.g., Figure 2 in P. Gao et al. 2023), first by increasing absorption in the 15 µm band (which almost, but not quite, reaches into the MIRI LRS wavelength range) and then by also adding an absorption feature between 9 and 10 µm in the 3000× solar spectrum.

With increasing metallicity, the dayside emission increases at almost all wavelengths, while the change in the nightside spectra is more complicated. On the nightside, the emission drops dramatically in the H₂O and CH₄ band as metallicity increases. However, in the spectral window beyond 9 µm, the emission first increases substantially for the 100× and 300× solar simulations, with an emitted flux close to the dayside flux in this spectral region. These spectral windows are probing deeper layers, where atmospheric circulation is efficient at transporting heat and day-night gradients are small. Only with the highly increased CO₂ abundance in the 3000× solar simulation does the nightside flux drop off in this wavelength region, falling below the flux of the solar-metallicity simulation.

### 4.2.1. Dayside Spectra

Even though the change in the band-integrated emitted flux due to hazes varies dramatically between the three sets of haze optical properties we examined, we identify a common pattern in how the spectra change with increasing haze opacity for all three haze properties. First, the amplitude of absorption features decreases until the spectrum is almost featureless (i.e., probing an isothermal portion of the atmosphere). Then, emission features appear (Figure 12). This behavior is consistent with previous studies on the effect of photochemical hazes on emission spectra of sub-Neptunes in 1D (C. V. Morley et al. 2015) and hot Jupiters in 1D (P. Lavvas & A. Arfaux 2021; A. Arfaux & P. Lavvas 2022) and 3D (M. E. Steinrueck et al. 2023). However, the overall spectral shape and where this transition happens depend on haze properties and metallicity.

For soot hazes, a thermal inversion shows up at lower haze opacities—even at the nominal haze opacity—compared to the less absorbing hazes, resulting in emission features. The only exception to this is the 3000× solar model, which is almost completely featureless from 5 to 10 µm and exhibits small absorption features at shorter and longer wavelengths. Because
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Dayside emission spectra at secondary eclipse for multiple haze opacities and haze properties, for a metallicity of 300× solar (top row) and 3000× solar (bottom row). With increasing haze opacity, spectral features seen as absorption features in the haze-free case gradually turn into emission features. Gray lines indicate spectra for a planet emitting blackbody radiation at multiple temperatures (indicated by the labels) for comparison. The JWST MIRI LRS secondary eclipse spectrum of E. M. R. Kempton et al. (2023) and the alternative data reduction of I. Malsky et al. (2025) are shown as red circles and purple triangles, respectively, with error bars.</figcaption></figure>

the thermal inversion is driven by additional heating through absorption by soot hazes (see Section 3.1), the emitted planetary flux increases substantially in molecular absorption bands such as the water and methane bands from around 6 to 8 µm.

For the maximally reflective hazes, as described in Section 3.1, the thermal inversion is driven entirely by cooling in deeper layers, due to incoming starlight being reflected at high altitudes. Thus, the overall dayside emission drops off dramatically with increasing haze opacity. At the nominal haze opacity, there is a decreased amplitude of absorption features at all metallicities, but no emission features are visible yet. For 10× the nominal haze opacity, the water/methane feature from 5 to 9 µm shows up in emission for all metallicities except for the 3000× solar case. Increasing the haze opacity even further to 100× leads to a drastic decrease in the overall flux and even stronger emission features. For this haze opacity, the 3000× solar case also exhibits emission features.

We note that because the emissivity of the maximally reflective hazes effectively is zero owing to their high single-scattering albedo, no photons are emitted from the haze particles (e.g., J. Taylor et al. 2021). Consequently, the emission effectively emerges from regions with significant

gas-phase opacity, and thus from much deeper regions of the atmosphere than in the soot case. Therefore, even though the optical depth (which includes scattering and absorption) reaches unity at $\approx 10^{-5}$ bars in the 100× haze opacity case, the emission features indeed probe the thermal inversion at $\approx 10^{-3}$ to $10^{-4}$ bars in the 300× solar and 3000× solar-metallicity cases.

For tholin hazes, the dayside spectrum transitions from absorption to emission features at similar haze opacities to those for the maximally reflective hazes. However, due to the additional heating due to absorption by hazes at low pressures, the emitted flux does not drop nearly as dramatically with increasing haze opacity as in the maximally reflective case. In fact, in the center of the water and methane bands, the emitted flux stays roughly the same or slightly increases, while it drops by less than 30% outside the band even in the 100× haze production cases. As with the other haze properties, for the 3000× solar simulation, a higher haze opacity (100×) is needed to produce emission features, again consistent with the fact that at haze opacities $< 100$× the nominal one, no thermal inversion forms. This again confirms the picture emerging in previous sections that for higher metallicities a larger haze
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Nightside emission spectra for multiple haze opacities and haze properties, for a metallicity of 300× solar (top) and 3000× solar (bottom). With increasing haze opacity, the nightside flux drops off owing to a combination of cooling of the deep atmosphere and emission emerging from higher altitudes. Gray lines indicate spectra for a planet emitting blackbody radiation at multiple temperatures (indicated by the labels) for comparison. The JWST MIRI LRS nightside spectrum of E. M. R. Kempton et al. (2023) and the alternative data reduction of I. Malsky et al. (2025) are shown as red circles and purple triangles, respectively, with error bars.</figcaption></figure>

opacity is needed to alter the atmospheric structure than for lower metallicities.

### 4.2.2. Nightside Spectra

The spectral shape of the nightside emission spectra (Figure 13) remains similar to the clear-atmosphere spectrum except for very high haze opacities. However, the emitted flux is reduced for all

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Models that provide a relatively good match to the observed JWST MIRI LRS white-light phase curve of GJ 1214b, compared to the data reductions of E. M. R. Kempton et al. (2023) and I. Malsky et al. (2025). Generally, models with a metallicity $≥100×$ solar and tholins or maximally reflective hazes provide the closest qualitative match.</figcaption></figure>

hazy simulations, with a larger reduction in emitted flux for higher haze opacities. Two factors contribute to this reduction in flux: the cooling of deep layers of the atmosphere, and the additional haze opacity resulting in the emitted radiation emerging from lower pressure levels with lower temperatures.

When comparing identical haze opacities, the reduction in flux is smallest for tholins and largest for soots. Tholins have a lower extinction cross section than soots at most wavelengths; thus, the outgoing emission probes deeper, hotter layers than for soot or maximally reflective hazes. Because we assumed an extinction cross section identical to soot for the maximally reflective hazes, the outgoing emission should emerge from roughly similar regions in the atmosphere for soot and maximally reflective hazes. The lower nightside flux in the soot case thus can be attributed to the stronger cooling of the deep atmosphere. As the haze opacity increases, the region probed by the outgoing emission shifts to higher regions in the atmosphere, at some point reaching the region in which the temperature profile is close to isothermal (Figure 3). This explains why the emission spectra are almost featureless for the 100× haze production cases.

## 5. Comparison to GJ 1214b Phase-curve Observations

In this section, we compare our models to the recent phase-curve observations of GJ 1214b with JWST MIRI/LRS (Figure 14). A comparison between most of our models and the observed phase curve of GJ 1214b was already presented in E. M. R. Kempton et al. (2023). They conclude that both a high metallicity ($≥100×$ solar) and a high (10× nominal) haze opacity with hazes that are more reflective than tholins are needed to qualitatively match the white-light phase curve. These models roughly reproduce the observed dayside flux and phase-curve offset, with observations falling in between the tholin and “maximally reflective” models, and also produce a large phase-curve amplitude. None of the models, however, predict nightside fluxes as low as in the data reduction presented in that paper. It is worth noting that E. M. R. Kempton et al. (2023) stressed that the nightside flux was the most uncertain aspect of the observational analysis.
One aspect not discussed in E. M. R. Kempton et al. (2023) is that the dayside spectra of almost all models that qualitatively match the phase curve (i.e., that have a haze opacity high enough to cool down the dayside flux substantially) exhibit quite prominent emission features. The only exception to that is the 3000× solar model with 10× haze opacity, which shows an almost featureless spectrum with very subtle emission features. In contrast, the secondary eclipse spectrum reported in E. M. R. Kempton et al. (2023) shows subtle absorption features. The observed shape of the nightside spectrum matches models with a metallicity ≳100× solar and ≈10× the nominal haze opacity qualitatively, though the flux between 6 and 8 µm dips noticeably below the values predicted by these models.

I. Malsky et al. (2025) present an alternative data reduction that includes a previously unidentified potential systematic (a sudden drop in flux near the phase-curve minimum). Including this additional systematic leads to a much higher nightside flux. The phase-curve amplitude and flux at the phase-curve maximum are now close to a range of models, especially the 100× and 300× solar models with tholins and 10× the nominal haze opacity. We note that, in terms of amplitude and maximum flux, the phase curve also lies between the 300× solar models with maximally reflective hazes with 1× and 10× the nominal haze opacity and is relatively close to the 3000× solar model with tholins and 10× the nominal haze production. However, all of the models predict a phase-curve offset that is much smaller than the observed offset using the I. Malsky et al. (2025) reanalysis.

The shape of the dayside spectrum also changes with the reanalysis, becoming closer to a blackbody spectrum (Figure 12). This could be indicative of hazes on the dayside. However, the tension between our models predicting substantial emission features and the lack of emission features in the observed spectrum remains. The shape of the nightside spectrum matches the models that provide the closest match to the white-light phase curve relatively well when ignoring the likely unreliable data points with $\lambda > 10.5$ µm. Past 10.5 µm, the data exhibited considerable red noise and the phase-curve fits in these wavelength bins were poor (see I. Malsky et al. 2025).

To utilize the spectral information from other phases as well, we proceed to look at the spectral phase curves of the models we identified as qualitatively matching the white-light phase curve best (Figures 15 and 16). For the simulations with tholins among this set (Figure 15), generally the models overpredict the emitted flux at wavelengths < 8.5 µm for a substantial portion of the phase curve (centered at the dayside). At wavelengths >9.5 µm, in turn, the 3000× solar phase curve matches the observed maximum and minimum fluxes (especially the Kempton et al. reduction) well but is shifted toward a smaller phase offset than observed, while the lower-metallicity simulations exhibit a lower phase-curve amplitude and higher nightside fluxes than observed. The pattern of overpredicting flux within the water and methane absorption bands while matching the phase curve better in channels outside these absorption bands can be attributed to the emission features exhibited by the models but not seen in the data. Seeing this discrepancy with observations arise in the spectral phase curve as well demonstrates that the thermal inversion dominates the model spectra throughout the entire dayside. Thus, the mismatch to the observed spectral shape is not just limited to a narrow phase region near the secondary eclipse.

The simulations with maximally reflective hazes that provide a reasonable match to the white-light phase curve, in turn, exhibit a closer match to the observed spectral phase curve at wavelengths < 9 µm but underpredict the flux at almost all phases at longer wavelengths. For the 100× and 300× solar simulations, the phase-curve amplitude also is smaller than observed at these wavelengths, while the 3000× solar phase curve provides a somewhat better match to observations in the longest-wavelength channels. The pattern of having lower-than-observed flux at wavelengths < 9 µm can also be attributed to the presence of emission features in most model-predicted spectra, which is not seen in the observations. An exception to this behavior is the 300× solar simulation with nominal haze opacity, which overpredicts the flux in all channels at all phases except for the 7.75 and 8.25 µm channels, in which the emitted flux falls below the Malsky et al. reduction but above the Kempton et al. reduction for the phases in which primarily the nightside is visible.

We further note that the trend of higher phase-curve amplitudes with higher metallicities, when the haze properties are held constant, holds true even in the narrow wavelength channels.

To summarize, some of our models are able to match multiple aspects of the observed phase curve of GJ 1214b well, especially the low dayside flux, the shape of the nightside spectrum, and, in the case of the I. Malsky et al. (2025) reanalysis, the phase-curve amplitude. However, remaining tensions are the observed lack of emission features on the dayside, the low nightside flux in the case of the Kempton data reduction, and the high observed phase-curve offset in the case of the Malsky reanalysis. One simulation (3000× solar, maximally reflective hazes, 10× haze opacity) only shows very subtle emission features and thus qualitatively provides a somewhat better match to the secondary eclipse spectrum and spectral phase curve. Yet this model overall overestimates the albedo of the planet and thus emits too little flux at almost all phases when integrated across the MIRI bandpass. Our model may be too simplistic to fully match the observations. For example, simultaneously including condensate clouds could lower the haze production rate that is required to explain the low dayside flux (I. Malsky et al. 2025) and result in a more blackbody-like spectrum. Using a more consistent haze profile derived from microphysics models run for each atmospheric composition and haze production rate (rather than applying the same haze profile to all atmospheric compositions and scaling the haze optical depth) could also change the vertical temperature structure and thus the strength of the thermal inversion. Spatially inhomogeneous hazes may also be possible (M. E. Steinrueck et al. 2021), and including a dynamic feedback between atmospheric transport of hazes, atmospheric circulation, and temperature structure could further impact the emission spectra and phase curve (M. E. Steinrueck et al. 2023). Finally, we note that the haze optical properties remain highly uncertain. While we explored multiple limiting cases, exploring a range of possible optical properties, such as those presented in recent studies (L. Corrales et al. 2023; C. He et al. 2024), may also improve the match to observations.

## 6. Recommended Observing Strategy

As discussed in Section 4.1, photochemical hazes add a considerable spread to the phase-curve amplitude–MMW/ phase-curve amplitude–metallicity trend, making it impossible to distinguish between a high-metallicity clear atmosphere and
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Spectral phase curves for a selection of models with tholin hazes that provide a relatively good match to the observed white-light phase curve. The models substantially overpredict the emitted flux for a large fraction of the phase curve in spectral regions that show emission features in the secondary eclipse spectrum (6–8.5 µm). The observed spectral phase curves are plotted in the background in light gray with filled circles and solid lines representing the E. M. R. Kempton et al. (2023) reduction and open circles and dashed lines representing the I. Malsky et al. (2025) reduction. We note that, due to the relatively low resolution of our post-processing, the spectral bins of the models do not line up exactly with the spectral bins of the observations (see Appendix C for more details).</figcaption></figure>

a lower-metallicity atmosphere with soot hazes. At the same time, intermediate<fnref n="17" /> haze opacities can lead to an isothermal atmosphere on the dayside of the planet in pressure regions

probed by emission spectroscopy, and thus a blackbody-like emission spectrum, again prohibiting the determination of the composition of the atmosphere. Here we argue that a combination of a band-integrated phase curve and a secondary eclipse spectrum can break the degeneracy between metallicity and haze opacity and thus is a promising observing strategy for sub-Neptunes.

In Section 4.1, we observed that generally the phase-curve amplitude will only be affected by hazes if the haze opacity is large enough to affect the temperature structure of the atmosphere. In Section 4.2, we found that for these cases the dayside emission spectrum shows dramatically reduced

<footnote n="17">¹⁷ We use the term “intermediate” to mean intermediate within our chosen grid of parameters. However, when comparing to the predicted photochemical mass flux of haze precursors in photochemistry simulations, assuming haze yield rates similar to the solar system, and assuming that the haze opacity roughly scales with haze production rate, the nominal (1× haze opacity) scenario would be considered a typical assumption and the 10× haze opacity scenario would be considered high, while the 100× haze opacity scenario would be considered extremely high—comparable to the surprisingly high haze production rates inferred for GJ 1214b by P. Gao et al. (2023) and K. Ohno et al. (2024).</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Same as Figure 15, but for a selection of models with maximally reflective hazes that provide a relatively good match to the observed white-light phase curve. The models provide a good qualitative match at wavelengths < 9.5 µm (Kempton reduction)/<9.5 µm (Malsky reduction) but underestimate the flux and, in some channels, amplitude in the longest-wavelength channels.</figcaption></figure>

absorption features, an almost featureless spectrum, or emission features. Thus, the dayside spectrum can inform the interpretation of the phase curve: If there are prominent absorption features present, with an amplitude comparable to what is expected for a clear atmosphere, it is safe to use the clear-atmosphere phase-curve amplitude–metallicity relationship to deduce the metallicity. However, if emission features are present or the dayside spectrum looks blackbody-like, the phase curve should be compared to GCM simulations that include photochemical hazes for interpretation. An example, comparing the emission spectra of two atmospheres with nearly

identical phase curves, is shown in the top row of Figure 17. The emission spectrum of the higher-metallicity, clear-atmosphere model exhibits easily visible absorption features, while for the hazy, lower-metallicity atmosphere the same spectral features show up in emission.

Vice versa, a white-light phase curve adds complementary information to the secondary eclipse spectrum and can be crucial for narrowing down the atmospheric composition. As we demonstrated in Section 4.2, intermediate (10×) haze opacities can lead to almost featureless emission spectra that contain little information on the atmospheric composition. In
<figure><img src="image_19.png" /><figcaption>Example for similar secondary eclipse spectra</figcaption></figure>

<figure><img src="image_20.png" /></figure>

**Figure 17.** This figure demonstrates the advantages of our recommended observing strategy of combining a white-light phase curve with a secondary eclipse spectrum. The top row shows an example of two models with different atmospheric composition with nearly identical phase curves but differing emission spectra, while the bottom row features an example of two models with nearly identical emission spectra but vastly differing white-light phase curves. The left column shows the white-light phase curve in the MIRI bandpass, and the right column shows the secondary eclipse spectrum.

these cases, a phase curve can distinguish between a higher- metallicity and a lower-metallicity scenario. An example is shown in the bottom row of Figure 17.

For high haze opacities, our models predict strong emission features. In this scenario, it should be possible to obtain constraints on atmospheric composition and haze production rate from the emission features in the secondary eclipse spectrum. However, a caveat to this is that our models seem to overpredict the strength of emission features when compared to the emission spectrum of GJ 1214b. A more thorough examination of the scenario of high haze production rates with more complex haze models (e.g., considering the change in particle size distribution, and thus opacity structure, with changing metallicity and haze production rate) may thus be necessary to correctly interpret scenarios with hints of high haze opacity.

Overall, we conclude that the combination of a white-light phase curve and a secondary eclipse spectrum is promising for characterizing the atmospheres of sub-Neptunes. For some sub- Neptunes, like GJ 1214b, the signal-to-noise ratio is good enough that a meaningful dayside spectrum can be obtained from a single secondary eclipse or two secondary eclipses. For

these planets, a secondary eclipse spectrum will automatically be included with a phase-curve observation.<fnref n="18" /> However, for fainter targets, it may be necessary to stack multiple secondary eclipses to get a reasonable emission spectrum. For some of these cases, a band-integrated phase curve could still be obtained in a single visit. We suggest that in such cases obtaining additional secondary eclipses (but not full phase curves, which would be much more time intensive) is the best observing strategy. This will be particularly important for planets for which transmission spectroscopy has not been able to determine the composition of the atmosphere, for example, due to high-altitude aerosols (like GJ 1214b) or due to contamination from starspots, which generally has a lower impact on emission spectra.

<footnote n="18">¹⁸ For the very few targets that are bright enough that a spectral phase curve can be obtained, the spectral phase curve and in particular the nightside spectrum can be highly valuable in determining the composition of the atmosphere as well. Nightside spectra often will still exhibit spectral features even in the case of a featureless secondary eclipse spectrum, as evidenced by both our models (Figure 13) and the observations of GJ 1214b. However, even for the best-case scenario of GJ 1214b, the nightside spectrum and spectral phase curves have significant uncertainty. The number of sub-Neptunes for which a spectral phase curve is feasible thus is severely limited.</footnote>
## 7. Conclusion

We examined how photochemical hazes impact the global temperature structure, atmospheric circulation, phase curves, and emission spectra of sub-Neptunes with a variety of atmospheric compositions and haze properties using the example of GJ 1214b. To do this, we utilized a set of simulations with a GCM with wavelength-dependent radiative transfer. We assumed horizontally uniform hazes with a vertical profile derived from a microphysics model. Scattering and absorption from hazes are taken into account. We explored three different assumptions for the haze optical properties: soot, tholin, and “maximally reflective” hazes, which use the extinction cross section of soot but have a single-scattering albedo artificially set to 0.9999. We assumed that the haze vertical profile does not vary with atmospheric composition. To simulate adjusting the haze production rate, we scaled the optical depth by a constant factor. Our main conclusions can be summarized as follows:

1. The effect of photochemical hazes on the 3D temperature structure of sub-Neptunes strongly depends on the assumed haze optical properties. On the dayside, for the highly absorbing soot hazes, the atmosphere heats up at low pressures compared to the clear-atmosphere case, while deeper layers experience strong cooling. The result is a strong thermal inversion. For the other extreme case we considered, “maximally reflective hazes,” the deeper layers of atmosphere cool, while the temperature at low pressures remains unchanged. A smaller thermal inversion forms owing to this cooling. For tholin hazes, which reflect a large fraction of the starlight but still absorb some of it, the changes are qualitatively between both of the other cases, with some warming at low pressures.

2. The nightside temperatures are strongly determined by the amount of dayside cooling in deeper layers of the atmosphere. Thus, soot simulations tend to have the coldest nightsides at a given haze production rate.

3. Photochemical hazes change the atmospheric circulation of sub-Neptunes. For high-enough haze opacities, photochemical hazes tend to produce a single narrow equatorial jet, even for atmospheric compositions for which clear-atmosphere simulations predict a broad equatorial jet or multiple jets. The higher the haze production rate is, the faster the jet speed.

4. Hazes tend to reduce the amplitude of absorption features on the dayside for moderate levels of haze production opacity and produce emission features for high haze production rates. With increasing metallicity, a higher haze opacity is needed for a thermal inversion to form and thus for emission features to appear in the dayside spectrum.

5. Generally, phase-curve amplitudes stay the same for low haze opacities and increase for higher haze opacities. Soot hazes result in larger increases to the phase-curve amplitude than tholins or highly reflective hazes.

6. Hazes significantly increase the scatter in the phase-curve amplitude–MMW relationship. While the relationship generally holds true for hazy simulations if haze opacity and optical properties are held constant, neither haze production rates nor optical properties are currently known, and they may vary strongly between planets. It

is thus no longer possible to estimate the MMW from the phase-curve amplitude alone.

7. However, when combining white-light phase curves with secondary eclipse spectra, the degeneracy between haziness and metallicity can be disentangled: if the haze abundance is high enough to affect the phase-curve amplitude, the secondary eclipse spectrum will show signatures of hazes such as a close-to-isothermal temperature profile or emission features. Thus, phase curves continue to be a promising technique to characterize the atmospheres of sub-Neptunes for which transmission spectroscopy has proven to be challenging.

However, a few caveats remain. While our models qualitatively explain many features of the GJ 1214b phase curve (dayside flux, shape of nightside spectra, high phase-curve amplitude), some tensions remain: Almost all models that match the white-light dayside flux show strong emission features (with one model showing subtle emission features), when the data show subtle absorption features or blackbody spectra. Furthermore, none of our models can produce a nightside flux as low as in the Kempton reduction or a phase-curve offset as large as in the Malsky reduction. It is likely that a more realistic haze model is needed to fully explain the observations (see discussion at the end of Section 5).

To fully utilize the potential of phase curves in the study of sub-Neptunes, we will need robust measurements of the planet-to-star flux at better than 100 ppm precision. This is particularly true for MIRI LRS, which covers the wavelength range where the emission from many hazy sub-Neptunes such as GJ 1214b peaks and thus is the instrument of choice. However, MIRI data have the most complex instrument systematics of any JWST instrument. While data reduction for a relatively high signal-to-noise ratio hot Jupiter phase curve led to consistent results across multiple teams (T. J. Bell et al. 2024), for the more challenging case of GJ 1214b there are currently two independent analyses of the phase curve that do not agree, particularly on the nightside flux. We are still in the early days of JWST data analysis, and history has shown that for HST and Spitzer observations it took several years to build consensus on the best way to reduce and analyze the data. The high scientific potential of JWST phase-curve observations motivates a continued push to improve our understanding of the instruments and their behavior in coming years.

## Acknowledgments

We thank Channon Visscher for providing equilibrium chemistry abundance tables. M.S. and M.Z. acknowledge support from the 51 Pegasi b Fellowship, funded by the Heising-Simons Foundation. This work is based in part on observations made with the NASA/ESA/CSA James Webb Space Telescope. These observations are associated with program #1803. Support for this program was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. R.L. was supported by NASA XRP grant 80NSSC22K0953 and STScI grant JWST-AR-01977.007-A. A portion of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). Computations were performed on the
HPC system Vera at the Max Planck Computing and Data
Facility (MPCDF).

Facilities: JWST.

*Software:* MITgcm (A. Adcroft et al. 2004), Numpy (C. R. Harris et al. 2020), SciPy (P. Virtanen et al. 2020), Astropy (Astropy Collaboration et al. 2022), matplotlib (J. D. Hunter 2007).

## Appendix A
### Haze Optical Properties

The differential optical depth and single-scattering albedos of the haze profiles we use are shown in Figure 18. An alternative visualization of the haze properties and also the asymmetry parameter can be found in Figure 1 of I. Malsky et al. (2025).

<figure><img src="image_21.png" /><figcaption>**Figure 18.** Differential optical depth (as a proxy for extinction cross section; top row) and single-scattering albedo (bottom row) for the assumed haze profile for soot (left) and tholin (right) hazes at different pressures.</figcaption></figure>
Appendix B
Contribution Functions of Clear-atmosphere Simulations

The contribution functions of the clear-atmosphere simulations are shown in Figure 19.

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Wavelength-dependent contribution functions for the dayside (left column) and nightside (right column) emission spectra. Shown are the clear-atmosphere simulations.</figcaption></figure>
# Appendix C
## Bins for Spectral Phase Curve

Due to the fairly low resolution of our post-processing (196 custom wavelength bins, with the correlated-k method being used to calculate the outgoing flux within each bin), the wavelength bins of the models do not line up exactly with the edges of the wavelength bins used in the observational spectral phase curves. For transparency, we show the wavelength bins used to generate Figures 15 and 16 in Table 3. We found that the figures only changed marginally when choosing to display a slightly different set of the post-processing wavelength bins for display in the plots.

<table><caption>**Table 3**
Wavelength Bins for Spectral Phase Curves</caption><thead><tr><th>Plot Title<br>(μm)</th><th>λ<sub>obs,0</sub><br>(μm)</th><th>λ<sub>obs,1</sub><br>(μm)</th><th>λ<sub>model,0</sub><br>(μm)</th><th>λ<sub>model,1</sub><br>(μm)</th></tr></thead><tbody><tr><td>5.25</td><td>5.0</td><td>5.5</td><td>5.08</td><td>5.45</td></tr><tr><td>5.75</td><td>5.5</td><td>6.0</td><td>5.45</td><td>5.98</td></tr><tr><td>6.25</td><td>6.0</td><td>6.5</td><td>5.98</td><td>6.45</td></tr><tr><td>6.75</td><td>6.5</td><td>7.0</td><td>6.45</td><td>7.00</td></tr><tr><td>7.25</td><td>7.0</td><td>7.5</td><td>7.00</td><td>7.62</td></tr><tr><td>7.75</td><td>7.5</td><td>8.0</td><td>7.62</td><td>8.00</td></tr><tr><td>8.25</td><td>8.0</td><td>8.5</td><td>8.00</td><td>8.54</td></tr><tr><td>8.75</td><td>8.5</td><td>9.0</td><td>8.54</td><td>9.20</td></tr><tr><td>9.25</td><td>9.0</td><td>9.5</td><td>9.20</td><td>9.50</td></tr><tr><td>9.75</td><td>9.5</td><td>10.0</td><td>9.50</td><td>9.80</td></tr><tr><td>10.25</td><td>10.0</td><td>10.5</td><td>9.80</td><td>10.40</td></tr></tbody></table>

## ORCID iDs

Maria E. Steinrueck https://orcid.org/0000-0001-8342-1895
Vivien Parmentier https://orcid.org/0000-0001-9521-6258
Laura Kreidberg https://orcid.org/0000-0003-0514-1147
Peter Gao https://orcid.org/0000-0002-8518-9601
Eliza M.-R. Kempton https://orcid.org/0000-0002-1337-9051

Michael Zhang https://orcid.org/0000-0002-0659-1783
Kevin B. Stevenson https://orcid.org/0000-0002-7352-7941
Isaac Malsky https://orcid.org/0000-0003-0217-3880
Michael T. Roman https://orcid.org/0000-0001-8206-2165
Emily Rauscher https://orcid.org/0000-0003-3963-9672
Matej Malik https://orcid.org/0000-0002-2110-6694
Roxana Lupu https://orcid.org/0000-0003-3444-5908
Tiffany Kataria https://orcid.org/0000-0003-3759-9080
Anjali A. A. Piette https://orcid.org/0000-0002-4487-5533
Jacob L. Bean https://orcid.org/0000-0003-4733-6532
Matthew C. Nixon https://orcid.org/0000-0001-8236-5553

## References

Ackerman, A. S., Toon, O. B., & Hobbs, P. V. 1995, JGR, 100, 7121

Adams, D., Gao, P., de Pater, I., & Morley, C. V. 2019, ApJ, 874, 61

Adcroft, A., Campin, J.-M., Hill, C., & Marshall, J. 2004, MWRv, 132, 2845

Arfaux, A., & Lavvas, P. 2022, MNRAS, 515, 4753

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Bean, J. L., Kempton, E. M. R., Fu, G., et al. 2021a, JWST Proposal, 1803

Bean, J. L., Raymond, S. N., & Owen, J. E. 2021b, JGRE, 126, e06639

Bell, T. J., Crouzet, N., Cubillos, P. E., et al. 2024, NatAs, 8, 879

Bell, T. J., Dang, L., Cowan, N. B., et al. 2021, MNRAS, 504, 3316

Benneke, B., Roy, P.-A., Coulombe, L.-P., et al. 2024, arXiv:2403.03325

Benneke, B., & Seager, S. 2012, ApJ, 753, 100

Benneke, B., & Seager, S. 2013, ApJ, 778, 153

Bohren, C. F., & Huffman, D. R. 1983, Absorption and Scattering of Light by Small Particles (New York: Wiley)

Brande, J., Crossfield, I. J. M., Kreidberg, L., et al. 2024, ApJL, 961, L23

Burn, R., Mordasini, C., Mishra, L., et al. 2024, NatAs, 8, 463

Charbonneau, D., Allen, L. E., Megeath, S. T., et al. 2005, ApJ, 626, 523

Charnay, B., Meadows, V., & Leconte, J. 2015a, ApJ, 813, 15

Charnay, B., Meadows, V., Misra, A., Leconte, J., & Arney, G. 2015b, ApJL, 813, L1

Christie, D. A., Mayne, N. J., Gillard, R. M., et al. 2022, MNRAS, 517, 1407

Cloutier, R., Charbonneau, D., Deming, D., Bonfils, X., & Astudillo-Defru, N. 2021, AJ, 162, 174

Corrales, L., Gavilan, L., Teal, D. J., & Kempton, E. M. R. 2023, ApJL, 943, L26

Crossfield, I. J. M., & Kreidberg, L. 2017, AJ, 154, 261

Fortney, J. J. 2005, MNRAS, 364, 649

Fortney, J. J., Mordasini, C., Nettelmann, N., et al. 2013, ApJ, 775, 80

Freedman, R. S., Lustig-Yaeger, J., Fortney, J. J., et al. 2014, ApJS, 214, 25

Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504

Gao, P., Piette, A. A. A., Steinrueck, M. E., et al. 2023, ApJ, 951, 96

Gao, P., Thorngren, D. P., Lee, E. K. H., et al. 2020, NatAs, 4, 951

Gharib-Nezhad, E., Marley, M. S., Batalha, N. E., et al. 2021, ApJ, 919, 21

Ginzburg, S., Schlichting, H. E., & Sari, R. 2018, MNRAS, 476, 759

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

He, C., Hörst, S. M., Lewis, N. K., et al. 2018, AJ, 156, 38

He, C., Hörst, S. M., Lewis, N. K., et al. 2020, PSJ, 1, 51

He, C., Radke, M., Moran, S. E., et al. 2024, NatAs, 8, 182

Hörst, S. M., He, C., Lewis, N. K., et al. 2018, NatAs, 2, 303

Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15

Hunter, J. D. 2007, CSE, 9, 90

Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Kataria, T., Showman, A. P., Fortney, J. J., Marley, M. S., & Freedman, R. S. 2014, ApJ, 785, 92

Kataria, T., Showman, A. P., Lewis, N. K., et al. 2013, ApJ, 767, 76

Kataria, T., Sing, D. K., Lewis, N. K., et al. 2016, ApJ, 821, 9

Kawashima, Y., & Ikoma, M. 2018, ApJ, 853, 7

Kawashima, Y., & Ikoma, M. 2019, ApJ, 877, 109

Kempton, E. M. R., Bean, J. L., Louie, D. R., et al. 2018, PASP, 130, 114401

Kempton, E. M. R., Zhang, M., Bean, J. L., et al. 2023, Natur, 620, 67

Khare, B. N., Sagan, C., Arakawa, E. T., et al. 1984, Icar, 60, 127

Knutson, H. A., Dragomir, D., Kreidberg, L., et al. 2014, ApJ, 794, 155

Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Natur, 505, 69

Lambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, A&A, 627, A83

Lavvas, P., & Arfaux, A. 2021, MNRAS, 502, 5643

Lavvas, P., & Koskinen, T. 2017, ApJ, 847, 32

Lavvas, P., Koskinen, T., Steinrueck, M. E., García Muñoz, A., & Showman, A. P. 2019, ApJ, 878, 118

Lee, E. J., & Chiang, E. 2016, ApJ, 817, 90

Lewis, N. K., Showman, A. P., Fortney, J. J., et al. 2010, ApJ, 720, 344

Line, M. R., & Parmentier, V. 2016, ApJ, 820, 78

Liu, B., & Showman, A. P. 2013, ApJ, 770, 42

Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1

Lupu, R. E., Zahnle, K., Marley, M. S., et al. 2014, ApJ, 784, 27

Luque, R., & Pallé, E. 2022, Sci, 377, 1211

Mai, C., & Line, M. R. 2019, ApJ, 883, 144

Malik, M., Grosheintz, L., Mendonça, J. M., et al. 2017, AJ, 153, 56

Malik, M., Kitzmann, D., Mendonça, J. M., et al. 2019, AJ, 157, 170

Malsky, I., Rauscher, E., Stevenson, K., et al. 2025, AJ, 169, 221

Marley, M. S., & McKay, C. P. 1999, Icar, 138, 268

Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85

Mordasini, C. 2020, A&A, 638, A52

Morley, C. V., Fortney, J. J., Kempton, E. M. R., et al. 2013, ApJ, 775, 33

Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2015, ApJ, 815, 110

Nettelmann, N., Fortney, J. J., Kramm, U., & Redmer, R. 2011, ApJ, 733, 2

Nixon, M. C., Piette, A. A. A., Kempton, E. M. R., et al. 2024, ApJL, 970, L28

Ohno, K., & Okuzumi, S. 2018, ApJ, 859, 34

Ohno, K., Okuzumi, S., & Tazaki, R. 2020, ApJ, 891, 131

Ohno, K., Schlawin, E., Bell, T. J., et al. 2025, ApJL, 979, L7

Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105

Parmentier, V., Fortney, J. J., Showman, A. P., Morley, C., & Marley, M. S. 2016, ApJ, 828, 22

Parmentier, V., Showman, A. P., & Fortney, J. J. 2021, MNRAS, 501, 78
Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91

Piaulet-Ghorayeb, C., Benneke, B., Radica, M., et al. 2024, ApJL, 974, L10

Piette, A. A. A., & Madhusudhan, N. 2020, ApJ, 904, 154

Rogers, L. A., Bodenheimer, P., Lissauer, J. J., & Seager, S. 2011, ApJ, 738, 59

Rogers, L. A., & Seager, S. 2010, ApJ, 716, 1208

Schlawin, E., Ohno, K., Bell, T. J., et al. 2024, ApJL, 974, L33

Shapiro, R. 1970, RvGSP, 8, 359

Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564

Steinrueck, M. E., Koskinen, T., Lavvas, P., et al. 2023, ApJ, 951, 117

Steinrueck, M. E., Parmentier, V., Showman, A. P., Lothringer, J. D., & Lupu, R. E. 2019, ApJ, 880, 14

Steinrueck, M. E., Showman, A. P., Lavvas, P., et al. 2021, MNRAS, 504, 2783

Taylor, J., Parmentier, V., Line, M. R., et al. 2021, MNRAS, 506, 1309

Thomas, G. E., & Stamnes, K. 1999, Radiative Transfer in the Atmosphere and Ocean (Cambridge: Cambridge Univ. Press),

Toon, O. B., Turco, R. P., Hamill, P., Kiang, C. S., & Whitten, R. C. 1979, JAtS, 36, 718

Tsai, S.-M., Steinrueck, M., Parmentier, V., Lewis, N., & Pierrehumbert, R. 2023, MNRAS, 520, 3867

Turco, R. P., Hamill, P., Toon, O. B., Whitten, R. C., & Kiang, C. S. 1979, JAtS, 36, 699

Valencia, D., Guillot, T., Parmentier, V., & Freedman, R. S. 2013, ApJ, 775, 10

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Wallack, N. L., Batalha, N. E., Alderson, L., et al. 2024, AJ, 168, 77

Zamyatina, M., Hébrard, E., Drummond, B., et al. 2023, MNRAS, 519, 3129

Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, PNAS, 116, 9723

Zhang, X., & Showman, A. P. 2017, ApJ, 836, 73