© ESO 2022

1 Introduction

We can directly observe the disk-integrated light from exoplanet atmospheres by measuring their phase curves, or the apparent brightness of the star-planet system as a function of orbital phase (see review by Parmentier & Crossfield 2018). Phase curve observations of close-in, hot exoplanets can reveal hints of atmospheric heat redistribution (Knutson et al. 2007), typically observed as a hot spot offset in longitude from the synchronously rotating sub-stellar point (Cowan & Fujii 2018; Madhusudhan 2019). By searching for variations in the phase curve with time, one can search for signatures of non-steady-state dynamics in exoplanet atmospheres (e.g. Heng & Showman 2015).

One technique for fitting the phase curves of hot exoplan-ets was presented by Cowan & Agol (2008), which essentially applies Fourier’s theorem to phase curves. It demonstrates that a mixture of a few sinusoidal components can fit most phase curves with a relatively small number of free parameters. Cowan & Agol (2008) study the information loss as one takes disk-integrated observations of exoplanet atmospheres, showing that the light curve can only sufficiently constrain roughly five free parameters. Knutson et al. (2007); Hu et al. (2015) fit the more demanding light curves of HD 189733 b and Kepler-7 b with a “beach-ball” model, with longitudinal segments at different temperatures.

Several works have used spherical harmonics to describe the maps of hot exoplanets. For example, Cowan et al. (2013) derived phase curves for spherical harmonic basis maps, Rauscher et al. (2018) fit for the temperature map of HD 189733 b with a spherical harmonic expansion and a dimensionality reduction technique, and Luger et al. (2019) derive closed-form solutions for phase curves and occultations of planetary flux maps represented with spherical harmonics. One strength of these approaches is that any planetary temperature map can be represented at sufficiently large degrees, and one possible weakness is that the flexible fits may correspond to implausible climate scenarios.

Matsuno (1966) and Longuet-Higgins (1968) previously showed that a reduced set of two-dimensional fluid equations (the so-called “shallow water system” or Laplace’s tidal equations) on a rotating sphere may be reduced to the quantum harmonic oscillator equation. The shallow water framework has been applied to the atmospheres of neutron stars (Heng & Spitkovsky 2009), the solar tachocline (Gilman 2000), the atmosphere of the Earth (Gill 1980), the top layer of the core of the Earth (Braginsky 1998) and the behaviour of submerged islands in the ocean of the Earth (Longuet-Higgins 1967). In the absence of rotation, the solutions to the shallow water framework are the familiar spherical harmonics. With rotation, the solutions are the so-called “parabolic cylinder functions”, which are combinations of exponential functions and Hermite polynomials. Heng & Workman (2014) further showed that when stellar heating (forcing) and drag/friction (both hydrodynamic and mag-netohydrodynamic) are added to the system, one still obtains the parabolic cylinder functions in certain limits. For brevity, we will refer to these functions as the h_mℓ generalised spherical harmonic basis and use them to describe exoplanet photospheric temperature maps.

There are several advantages to using the h_mℓ basis functions to describe an exoplanet’s 2D temperature map: (1) the input arguments to the special functions, α and ω_drag, are physically-meaningful parameters which represent properties of the atmosphere; (2) it is equally straightforward to fit the 2D basis function description of the temperature map to general circulation model (GCM) simulations and to observations of real exoplanets, allowing theory and observation to be compared in the same mathematical language; (3) even when limited to only a few free parameters, the h_mℓ basis can reproduce a diverse range of phase curve phenomena including offset hots pots and chevron-shaped hot regions. We will show in this work that some of the physically-meaningful parameters can be constrained using fits of the h_mℓ basis to GCMs, which reduces the total number of free parameters in fits to phase curve observations.

In Sect. 2, we outline the h_mℓ basis functions used to compute the temperature map and the phase curves from the temperature map. In Sect. 3, we develop formalism for measuring Bond albedos and redistribution efficiencies from 2D temperature maps. In Sect. 4, we outline a technique for constraining the latitudinal temperature distributions of hot exoplanets by fitting the h_mℓ basis functions to the 2D temperature maps produced by GCMs. In Sect. 5, we fit phase curve observations from Spitzer using priors extracted from the fits to the GCM temperature maps, and compare results obtained from fitting the Spitzer phase curves and the GCMs. We discuss the results in Sect. 6 and conclude in Sect. 7.

2 Phase curve model

Mathematically, h_mℓ is the perturbation to the height of the fluid in the so-called “shallow water system”, which is a proxy for the temperature perturbation (Cho et al. 2008). Therefore, we define a 2D temperature map as a function of latitude and longitude T (θ,ϕ): $$T(\theta ,\varphi )=\overline{T}\left(1+{\displaystyle \sum _{m,\mathcal{l}}^{{\mathcal{l}}_{\max }}{h}_{m\mathcal{l}}(\theta ,\varphi )}\right)$$(1)

The prefactor $\overline{T}$ on the right hand side acts as a constant scaling term for the mean temperature field, on which the h_mℓ terms are a perturbation. A plausible initial guess for the value of this prefactor is $\overline{T}=T_{\star }\sqrt{R_{\star }/a}f\text{'}{(1-A_B^{\text{'}})}^{1/4}$ for stellar effective temperature T_⋆, and normalized semimajor axis $a/R_{\star }.A_B^{\text{'}}$ is the Bond albedo and ƒ′ may be interpreted either geometrically (see Sect. 3.1) or as a greenhouse factor. We mark the ƒ′ and $A_B^{\text{'}}$ variables with primes here to denote that these two terms will be redefined in Sect. 3, and thus we do not report the values of the primed variables in this work. We note that the two terms in the factor $f\text{'}{(1-A_B^{\text{'}})}^{1/4}$ are perfectly degenerate.

For efficient sampling of all plausible temperature maps given the degenerate parameterization of the two-parameter term $f\text{'}{(1-A_B^{\text{'}})}^{1/4}$ in Eq. (1), we take the following conventions in generating temperature maps: (1) we set $A_B^{\text{'}}$ in Eq. (1) to zero always, and (2) we allow ƒ′ to vary uniformly from zero to unity, which extends beyond the range normally given to this parameter in the literature. The net effects of these choices are that (1) the two free parameters ƒ′ and $A_B^{\text{'}}$ are sampled as a single free parameter which describes $f\text{'}{(1-A_B^{\text{'}})}^{1/4}$, and (2) the parameters provided to the sampler $(f\text{'}\in [0-1],A_B^{\text{'}}=0)$ need not be tracked during sampling. Instead, we infer the Bond albedo and the heat redistribution efficiency for a given temperature map by computing the temperature map with the parameter $f\text{'}{(1-A_B^{\text{'}})}^{1/4}$, and then integrate outgoing flux from the planetary atmosphere to compute the Bond albedo and redistribution efficiency, as outlined in the Sect. 3.4.

The h_ml (α, ω_drag) terms are defined by Eq. (258) of Heng & Workman (2014), $$\begin{array}{l}{h}_{m\mathcal{l}}=\frac{C_{m\mathcal{l}}}{{\omega }_{\text{drag}}^2{\alpha }^4+m^2}{e}^{-{\tilde{\mu }}^2/2}[\mu m{H}_{\mathcal{l}}\cos (m\varphi )\hfill \\ +\alpha {\omega }_{\text{drag}}(2\mathcal{l}{H}_{\mathcal{l}-1}-\tilde{\mu }{H}_{\mathcal{l}})\sin (m\varphi )],\hfill \end{array}$$(2)

where α is a dimensionless fluid number that is constructed from the Rossby, Reynolds and Prandtl numbers. ω_drag is the dimensionless drag frequency (normalised by twice the angular velocity of rotation), $\mu =\cos \theta ,\tilde{\mu }=\alpha \mu ,\text{and}{H}_{\mathcal{l}}(\tilde{\mu })$ are the physicist’s Hermite polynomials: $$\begin{array}{l}{H}_0=1\hfill \\ H_1=2\tilde{\mu },\hfill \\ H_3=8{\tilde{\mu }}^3+12\tilde{\mu }\hfill \\ H_4=16{\tilde{\mu }}^4-48{\tilde{\mu }}^2+12\hfill \end{array}$$

The h_mℓ basis is derived without any special reference point in longitude. As a result, the hottest point on the atmosphere need not be the sub-stellar point. We parameterize shifts in the prime meridian of the h_mℓ basis from the sub-stellar point using a simple phase offset parameter, Δϕ, which describes a rotation of the coordinate system about the pole. It is important to note that the phase offset parameter Δφ is not equivalent to the traditional “hotspot offset” parameter often used in the literature, which measures the difference in longitude between the hottest longitude and the sub-stellar longitude. The phase offset Δϕ will vary with ω_drag, though in practice, drag values approaching and above of ω_drag ≳ 3 correspond to maps where the hottest point on the map is practically indistinguishable from the sub stellar point. In the remainder of this paper, we will report the phase offset Δϕ parameter with the implicit caveat that the distinction between phase offset and hotspot offset is negligible for the large drag frequencies discussed here.

For a phase curve fit with ℓ_max = 1, for example, the free parameters in the fit could include {C₁₁, α, ω_drag, ƒ′, Δϕ}. In practice however, α and ω_drag control the latitudinal concentration of heat in the planet’s atmosphere, so these two parameters are unconstrained by the phase curve alone. In Sect. 4, we will determine reasonable values for α and ω_drag based on fits with the h_mℓ basis to the 2D temperature maps produced by GCMs. These constraints can act as physically-motivated priors in the infrared phase curve fits which follow in Sect. 5.

The generalised spherical harmonic h_mℓ basis defines an “alphabet” with only a few free parameters which we can use to succinctly describe exoplanet temperature fields, shown in Fig. 1. Each panel represents an individual generalized spherical harmonic expansion term, and we reconstruct temperature maps of exoplanets by evaluating linear combinations of these basis maps.

We have fully defined the temperature map of the planet with the parameters above, and now we must construct the Planck function ℬ_λ(T(θ, ϕ)), integrate ℬ_λ within the bandpass of observations, and then disk-integrate to construct a thermal phase curve. We integrate the ratio of flux from the planet and the host star (Cowan & Agol 2011a): $$\frac{F_{\text{p}}}{F_{\star }}=\frac{1}{\pi I_{\star }}{\left(\frac{R_{\text{p}}}{R_{\star }}\right)}^2{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_{-\xi -\pi /2}^{-\xi +\pi /2}I}}(\theta ,\varphi )\cos (\varphi +\xi ){\sin }^2(\theta )\text{d}\varphi \text{d}\theta$$(3)

where the intensity I is given by $$I={\displaystyle \int \lambda }{\mathcal{F}}_{\lambda }{\mathcal{B}}_{\lambda }(T(\theta ,\varphi ))\text{d}\lambda$$(4)

for a filter response function ℱ_λ. The above relation neglects any reflected light. We assume the spectrum of each host star is accurately described by PHOENIX model spectra with the nearest effective temperature and surface gravity (Husser et al. 2013).

We have developed an efficient Python package for computing this triple integral for time series phase curve observations, written with Cython, theano and JAX for maximum speed and readability, called kelp, which is publicly available¹.

Fig. 1 Four montages of the first several terms in the generalized spherical harmonic expansion of the temperature map in the hmℓ basis. Each montage corresponds to a different combination of α and ωdrag, labeled in the upper right of each quadrant. Each map shows the temperature perturbation (purple to yellow is cold to hot) as a function of latitude and longitude (shown in Moll weide projections such that the sub stellar longitude is in the center of the plot). The m = 0 terms are always zero. The even-ℓ terms are asymmetric about the equator and therefore do not represent typical GCM results in equilibrium, so we keep all even-ℓ power coefficients fixed to zero in the subsequent fits (this approach assumes edge-on orbits, for exceptions see Cowan et al. 2017). In the upper two quadrants, when ωdrag is set to a smaller value, the chevron shape becomes more pronounced. In the two quadrants on the right, when the α is set to a larger value, the concentration of heat near the equator is more pronounced.

3 From 0D to 2D: Bond Albedos and Redistribution Efficiencies

3.1 Background

Zero-dimensional (0D) “box models” are used ubiquitously in exoplanet science. As a prime example, the “equilibrium temperature” of an exoplanet is standard fare in publications (e.g. Seager 2010). It originates from considering the stellar irradiation received by, and subsequently re-radiated χπ steradians of, a spherical exoplanet. The stellar luminosity is $4\pi R_{\star }^2{\sigma }_{\text{SB}}{T}_{\star }^4$ with R_⋆ being the stellar radius and σ_SB the Stefan-Boltzmann constant. The star is approximated as a spherical black body with a uniform photospheric temperature T_⋆. The stellar flux fills a celestial sphere with a surface area of 4πa², where a is the distance between the star and exoplanet. The exoplanet presents a cross section of πR² for intercepting this stellar flux, but absorbs only a fraction 1 – A_B of it with A_B being its Bond albedo. Therefore, the box model of luminosity received is $$L=4\pi R_{\star }^2{\sigma }_{\text{SB}}{T}_{\star }^4(1-A_{\text{B}})\frac{\pi R^2}{4\pi a^2}$$(5)

Assuming the exoplanet to radiate like a black body, the box model of luminosity re-radiated by the exoplanet is $$L^{\prime }=\chi \pi R^2{\sigma }_{\text{SB}}{T}^4$$(6)

where T is a characteristic temperature associated with the exoplanet. If we assume that its flux is re-radiated over χπ = 4π steradians, then one recovers the standard expression for the equilibrium temperature (e.g., Seager 2010), $$T_{\text{eq}}=\frac{T_{\text{irr}}}{\sqrt{2}}{(1-A_{\text{B}})}^{1/4},$$(7)

where $T_{\text{irr}}=T_{\star }\sqrt{R_{\star }/a}$ is the irradiation temperature. The equilibrium temperature is often quoted assuming A_B = 0 (e.g., López-Morales & Seager 2007).

More generally, one may write ƒ = 1/χ as a redistribution factor and derive an expression for the he mispherically-averaged day side temperature of the exoplanet, $$T_{\text{d}}=T_0{f}^{1/4}$$(8)

where we define² T₀ ≡ T_irr(1 − A_B)^1/4. By analogy, one may write the he mispherically- averaged night side temperature as T_n = T₀g^1/4. By assuming a linear relationship between ƒ and g and applying the boundary conditions that ƒ = g = 1/4 for full heat redistribution and ƒ = 2/3 and g = 0 for no redistribution, one obtains $$T_{\text{n}}=T_0{\left(\frac{2}{5}-\frac{3f}{5}\right)}^{1/4}$$(9)

One may also rewrite the preceding expressions in terms of a redistribution efficiency ε = 8/5 − 12 ƒ/5 that is bounded between 0 and 1, which yields Eqs. (4) and (5) of Cowan & Agol (2011b).

The day side and night side brightness temperatures of an exo-planet may be directly inferred from data. If one observes the secondary eclipse of an exoplanet purely in thermal emission, then the day side brightness temperature may be extracted in a model-independent way (Seager 2010). Extracting the night side brightness temperature requires a thermal phase curve to be measured as it is the difference in flux between the top of the transit and the bottom of the secondary eclipse (Schwartz et al. 2017). For example, the night side brightness temperature of 55 Cancri e was measured using a Spitzer phase curve (Demory et al. 2016a). It is worth noting that, in extracting the temperature map, Eq. (4) already accounts for thermal emission within a Spitzer bandpass. For a black body, the temperature and brightness temperature are equal.

Equations (8) and (9) may be inverted to solve for the Bond albedo and redistribution factor/efficiency, if the day side and night side (brightness) temperatures are known empirically. However, it is possible to obtain negative Bond albedos by using Eqs. (8) and (9). This artefact arises from approximating a two-dimensional (2D) situation with a 0D box model and is tied to the choice of boundary condition made for g, which we will next demonstrate. It affects the inference of A_B and ƒ or ε from all thermal phase curves of exoplanets.

3.2 General Box Models

Starting again from T_n = T₀g^1/4, we apply more general boundary conditions on g: $$\begin{array}{l}\text{Full redistribution}:f=\text{g}=\frac{1}{4}\hfill \\ \text{No}\text{redistribution}:f=f_0,\text{g}=0.\hfill \end{array}$$(10)

This yields a general expression for g in terms of ƒ and ƒ₀, $$g=\frac{f-f_0}{1-4f_0}$$(11)

Similarly, more general boundary conditions on ε are: $$\begin{array}{l}\text{Full}\text{redistribution}:f=\frac{1}{4},\epsilon =1\hfill \\ \text{No}\text{redistribution}:f=f_0,\epsilon =0,\hfill \end{array}$$(12)

which yields $$\epsilon =\frac{4(f-f_0)}{1-4f_0}=4g$$(13)

The general expression for the day side temperature is $$T_{\text{d}}=T_0{f}^{1/4}=T_0{\left(\frac{\epsilon }{4}-\epsilon f_0+f_0\right)}^{1/4}$$(14)

The general expression for the night side temperature is $$T_{\text{n}}=T_0{\left(\frac{f-f_0}{1-4f_0}\right)}^{1/4}=\frac{T_0{\epsilon }^{1/4}}{\sqrt{2}}$$(15)

In the limit of ƒ₀ = 2/3 (Hansen 2008), we recover Eq. (4) of Cowan & Agol (2011b) for the dayside temperature.

If we define the contrast between the night side and day side hemispheric temperatures as C = T_n/T_d, then the redistribution factor may be expressed as $$f=\frac{f_0}{1-C^4(1-4f_0)}$$(16)

3.3 Why Box Models May Produce Negative Bond Albedos

One may rewrite the expression for the day side temperature in Eq. (14) to obtain an expression for the Bond albedo, $$A_{\text{B}}=1-\frac{1}{f}{\left(\frac{T_{\text{d}}}{T_{\text{irr}}}\right)}^4$$(17)

If we demand that A_B < 0, then we obtain the following inequality, $$\frac{T_{\text{d}}}{T_{\text{irr}}}>f^{1/4}$$(18)

In the high-temperature limit, we expect C ≈ 0 and ƒ ≈ ƒ₀. Therefore, negative Bond albedos are obtained when $T_{\text{d}}/T_{\text{irr}}\gtrsim f_0^{1/4}$

It only remains to derive the value of ƒ₀. The standard value of ƒ₀ = 2/3 is obtained from adding an irradiation term to the expression for the intensity associated with Milne’s solution (Hansen 2008). Section 2.2 of Cowan & Agol (2011b) discusses how other choices may be made. In the context of our box models, we may write $$L^{\prime }=R^2{\displaystyle {\int }_{-\pi /2}^{\pi /2}{\displaystyle {\int }_0^{\pi }F}}\sin \theta \text{d}\theta \text{d}\varphi ,$$(19)

where the integration from ϕ = −π/2 to π/2 occurs only over the dayside hemisphere. If the day side flux of the exoplanet is F = σ_SBT⁴ and T is assumed to be constant, then we recover L'′= 2πR²σ_SB T⁴, χ = 2 and ƒ₀ = 1/2. If we account approximately for anisotropic emission over the day side hemisphere (caused by anisotropic stellar irradiation; Hansen 2008) by writing F = σ_SBT⁴ sin θ, then we obtain ƒ₀ = 2/π ≈ 0.64. This assumes that the radiative times cale is much shorter than the dynamical timescale (López-Morales & Seager 2007).

Table 8 of Wong et al. (2021) lists the dayside and irradiation temperatures of several hot Jupiters obtained from analysing Spitzer data. Specifically, it allows us to estimate T_d/T_irr values of 0.96 ± 0.04 for HAT-P-7b, 0.94 ± 0.05 for KELT-1b, 0.92 ± 0.02 for WASP-18b and 0.83 ± 0.02 for WASP-33b.

Figure 2 demonstrates that these estimates of T_d/T_irr satisfy the condition for negative Bond albedos stated in Eq. (18). Unsurprisingly, Table 8 of Wong et al. (2021) lists a negative Bond albedo for HAT-P-7b.

Generally, if the empirically-obtained value of T_d/T_irr exceeds about 0.9 (with this exact threshold depending on the assumed value of f₀, then the use of Eqs. (8) and (9) results in a negative value for the Bond albedo. This artifact originates from the global non-conservation of energy.

Fig. 2 Lower limit on Td/Tirr (curve) for producing negative Bond albedos. For comparison, empirically obtained values of Td/Tirr from Wong et al. (2021) are shown; the uncertainties on these values (not shown for plot clarity) are between 2% and 5%.

3.4 General Expression for Bond Albedo

We propose more general expressions for extracting A_B and ε directly from a two-dimensional temperature map, which have been previously proposed in Keating et al. (2019). Generally, the expression for L in Eq. (5) holds regardless of one’s assumption on f₀. However, the expression for L′in Eq. (6) needs to be generalised for application to infrared phase curves, $$L^{\prime }=R^2{\displaystyle {\int }_{-\pi }^{\pi }{\displaystyle {\int }_0^{\pi }F}}\sin \theta \text{d}\theta \text{d}\varphi .$$(20)

A general expression for inference of the Bond albedo is $$A_{\text{B}}=1-{\left(\frac{a}{R_{\star }}\right)}^2\frac{{\displaystyle {\int }_{-\pi }^{\pi }{\displaystyle {\int }_0^{\pi }F}}\text{s}\stackrel{.}{\text{in}}\theta \text{d}\theta \text{d}\varphi }{\pi \sigma T_{\star }^4}.$$(21)

The two-dimensional flux distribution F(θ, ϕ) may be inferred from an infrared phase curve, where the temperature T(θ, ϕ) is not assumed to be constant.

The redistribution efficiency is simply the ratio of the night-side to day side fluxes, $$\epsilon =\frac{{\displaystyle {\int }_{\pi /2}^{\pi }{\displaystyle {\int }_0^{\pi }F}}\sin \theta \text{d}\theta \text{d}\varphi +{\displaystyle {\int }_{-\pi }^{-\pi /2}{\displaystyle {\int }_0^{\pi }F}}\sin \theta \text{d}\theta \text{d}\varphi }{{\displaystyle {\int }_{-\pi /2}^{\pi /2}{\displaystyle {\int }_0^{\pi }F}}\sin \theta \text{d}\theta \text{d}\varphi }$$(22)

3.5 Limitations of Two-Dimensional Models

While more general than 0D box models, the 2D temperature-map approach has the following limitations:

• The atmospheric temperature generally varies in all three dimensions, implying that thermal phase curves measured at different wavelengths sample radial/vertical variations in temperature (Stevenson et al. 2014). In other words, a thermal phase curve at a specific wavelength probes a 2D surface; at other wavelengths, other surfaces are probed. These surfaces are generally not isobaric;

• The spectral energy distribution of the exoplanet atmosphere may deviate strongly from a blackbody, especially far away from the black body peak where the intensity is proportional to λ⁻⁴ (Rayleigh-Jeans law). This implies that the brightness temperature is not the true local temperature.

In the limit of sparse data, it is plausible to combine the analysis of infrared phase curves, measured in two different Spitzer channels, to infer a 2D temperature map. In the era of the James Webb Space Telescope (JWST), the need for a three-dimensional approach will be urgently needed as it will affect the inference of Bond albedo values from JWST phase curves.

4 2D Temperature maps from GCMs

The incentive is strong to infer 2D temperature maps for exoplanets, despite the lack of constraints from phase curve observations. Planets such as KELT-9 b have dayside atmospheres as hot as a K star and night side temperatures as hot as an M star, and the distribution of temperature with latitude may be similarly extreme, spanning ~ 15OO K from pole to equator (see below). These tremendous temperature variations suggest that different chemistry may be important within the same planetary hemisphere. For example, this has important implications in translating what we learn from the phase curve into expectations for observations at the planet’s terminator, which is probed during transmission spectroscopy.

The disk-integration that occurs when constructing a phase curve in Eq. (3) effectively discards information about the distribution of temperatures with latitude. In some cases this information can be recovered for the day sides of exoplanets on slightly inclined orbits during secondary eclipse (see, e.g.: Majeau et al. 2012; de Wit et al. 2012). In general, however, the distribution of temperatures with latitude is not constrained by the phase curve.

In this section, we explore whether we can use three-dimensional GCMs to inform inferred temperature maps of real hot Jupiters from phase curve observations. We introduce a set of GCMs in the next subsection, then outline a procedure for fitting them in the following subsection. Finally, we summarize the results in the limit of a planet with very high irradiation.

4.1 GCM Simulations

Given that phase curves of exoplanets do not constrain the latitudinal distribution of heat in an exoplanet atmosphere, we will “observe” some ubiquitous trends in model atmospheres which set the expectations for the latitude distribution of heat in real exoplanet atmospheres. Perna et al. (2012) studied heat redistribution and energy dissipation in hot Jupiter atmospheres using three-dimensional circulation models with dual-band or “double grey” radiative transfer. The authors used the GCMs to infer the effects of temperature inversions and stellar irradiation on hotspot offsets from the substellar point, among other quantities. We summarize their grid of 16 GCMs in Table 1. The greenhouse parameter γ is the ratio of the mean infrared to optical/visible opacity and controls whether a temperature inversion is present (Guillot 2010); γ = 0.5 and 2 correspond to temperature-pressure profiles without and with inversions, respectively.

4.2 Fits to The GCM Temperature Maps

4.2.1 Phase Offset

The temperature map formulation by Heng & Workman (2014) is constructed without any special definitions for an origin in longitude, ϕ. As a result, the temperature maps need a constant offset of π/2 to behave correctly in the limit of ω_drag ≫ 1, which ensures that the sub stellar point is also the hottest point on the map.

Similarly, we introduce an arbitrary phase offset as a fitting parameter, Δϕ, which allows the temperature map to be rotated with respect to the substellar point. While the h_mℓ basis could in principle fit any temperature map at sufficiently large spherical harmonic degree, maps with significant hotspot offsets would require ℓ_max ≫ 1 in order to approximate the offset temperature distribution, incurring many free C_mℓ parameters and making the inference intractable. Therefore, we introduce a single arbitrary rotation offset angle Δϕ to the fits of the temperature maps to concisely describe hotspot offsets at small ℓ_max.

4.2.2 Empirical limits on ℓ_max

As noted in the previous section, introducing the coordinate system offset angle Δϕ of the temperature map fits allows us to reproduce the GCM temperature maps with small ℓ_max. Choosing ℓ_max = 1 requires five free parameters: {ƒ, C₁₁,α, ω_drag, Δϕ}. The ℓ = 2 maps are not important for fitting mean phase curves of exoplanets as they are not symmetric across the planetary equator, and therefore should not have significant power for an atmosphere near equilibrium (see Fig. 1). ℓ_max = 3 has seven parameters, and so on. We first fit the 2D GCM temperature maps with the h_mℓ basis using ℓ_max = 3, and find that the three C_ml terms for ℓ = 3 are constrained to be very small compared to the power in C₁₁. We then inspect the quality of the fits for ℓ_max = 1 and ℓ_max = 3, and find that the ℓ_max = 1 maps are indeed strikingly similar to the ℓ_max = 3 maps, and thus we reason that ℓ_max = 1 is likely sufficient to approximate temperature maps for phase curve analyses, requiring up to five free parameters.

We note that previous studies have suggested that about five free parameters can be constrained by fits to a 1D phase curve (Cowan & Agol 2011a; Hu et al. 2015). It is encouraging that we arrive at a similar conclusion about the information content of phase curves as studies that take very different approaches to constructing phase curves, but as we will discuss in the next subsection, if we narrow our focus to the hottest exoplanets we can collapse the parameter space further.

4.2.3 Trends Towards High Equilibrium Temperature

Next, we look for trends with irradiation in the physically meaningful fitting parameters with which we can constrain the 2D temperature maps of real exoplanets. We fit the 2D temperature maps of the Perna et al. (2012) GCMs for the parameters which control the 2D temperature map in the h_mℓ basis with ℓ_max = 1, {ƒ, C_mℓ, α, ω_drag, Δϕ}, using Markov chain Monte Carlo (Foreman-Mackey et al. 2013). The GCM temperature fields and the maximum-likelihood h_mℓ maps are shown in Fig. 3. Comparison between the GCM results and the maximum-likelihood h_mℓ reconstructions confirms that the major features of the map are reproduced by the h_mℓ model at small ℓ_max, and more subtle features are reproduced at large ℓ_max. When α and ω_drag are allowed to vary, the concentration of heat near the equator is easily reproduced by the h_mℓ basis.

Remarkably, the four hottest temperature maps in Fig. 3 can all be fit with similar α and ω_drag parameters, α = 0.6 and ω_drag = 4.5. This is demonstrated by the posterior distributions for α and ω_drag for each of the four hottest GCM maps in Fig. 4. In other words, the latitudinal temperature distribution across maps with varying equilibrium temperature is quite similar when phrased in terms of the h_mℓ “alphabet”. The chevron shapes that are produced by this combination of α and ω_drag will be familiar to many GCM practitioners – it is a shape reproduced by many GCMs and has an elucidated theoretical basis (e.g.: Showman & Polvani 2011; Tsai et al. 2014).

At ℓ_max = 3 we are unable to reproduce some of the subtle features of the map, such as the equatorial jet that is noticeable in the night side atmosphere, or the “point” that forms at the leading edge of the day side atmosphere. These features could be reproduced by going to higher ℓ_max, but since we are focused here on jointly inferring properties from GCMs and phase curves, we assume that such fine features will remain undetectable with phase curve photometry, and neglect these small features. In the bottom row of Fig. 3, we show a linear h_mℓ solution to the GCM temperature maps with ℓ_max = 50, which better reproduces the precise shape of the dayside hotspot but not the jet, for example.

To summarize, our experiments fitting the h_mℓ maps to the GCMs have shown that: (1) we must introduce a phase offset that was not accounted for in the original definition of the h_mℓ basis; (2) ℓ_max = 1 or 3 is likely sufficient to capture the global-scale temperature inhomogeneities that are routinely in agreement between different GCMs and likely to be constrained by the phase curves; (3) the hottest maps all demand a common combination of α = 0.6 and ω_drag = 4.5, collapsing our parameter space from five to three free parameters.

5 Spitzer Observations

Next, we apply the h_mℓ basis functions in the context of phase curve observations from the Spitzer Space Telescope. Spitzer photometry at 3.6 and 4.5 μm provides coverage of spectral regions where the majority of the planetary flux is thermal emission rather than reflected light (more on this in Sect. 6.6). The Spitzer phase curves are disk-integrated measurements of the planet-to-star flux ratio, and thus constrain the longitudinal distribution of flux from the planet. In this section, we present original photometric extractions for seven of the eight targets presented, and interpret those phase curves with the h_mℓ basis.

Fig. 3 Comparison between the 2D temperature maps of hot Jupiter GCMs in the upper row (Perna et al. 2012), and the 2D temperature maps represented with the hmℓ basis in the middle and lower rows with ℓmax = 1 and 50, respectively. The horizontal axis is longitude, the vertical axis is latitude, and the center of each map is the substellar point. For model input parameters see Table 1; for an outline of the linear algorithm used to fit maps with ℓmax = 50, see Appendix A.

Fig. 4 Posterior correlations between α and ωdrag from fits of the hmℓ basis maps with ℓmax= 1 to the four hottest GCM temperature maps. The posterior distributions overlap at a few standard deviations for the values we adopt in the remainder of this work, α = 0.6 and ωdrag = 4.5 (marked with black lines).

5.1 Reduction

We retrieved archival spitzer/IRAC observations from the Spitzer Heritage Archive³, as published in the references listed in Table 2. The reduction of these AORs follows the procedure of Demory et al. (2016a), and we refer the reader to details in that reference. We model the IRAC intra-pixel sensitivity (Ingalls et al. 2016) using a modified implementation of the BLISS (BiLinearly-Interpolated Sub-pixel Sensitivity) mapping technique (Stevenson et al. 2012).

The baseline photometric model includes the pixel response function’s FWHM along the x and y axes, which significantly reduces the level of correlated noise as shown in previous studies (e.g. Lanotte et al. 2014; Demory et al. 2016b; Demory et al. 2016a; Gillon et al. 2017; Mendonça et al. 2018). The model does not include time-dependent parameters. We vary the sampling parameters with Markov chain Monte Carlo (MCMC) in a framework already presented in Gillon et al. (2012). We run two chains of 200 000 steps each to determine the phase-curve properties at 3.6 and 4.5 μm based on the entire dataset. We trimmed the ramp in the first 30 min of the first AOR of the HAT-P-7 Channel 1 visit.

The photometry of KELT-9 b in Channel 1 is being studied in a forthcoming work by Beatty et al. (in prep.), and therefore is excluded from this work. The photometry of HD 189733 is affected by stellar variability at the mmag level requiring a special analysis unique to this object, including time-dependent parameters which we exclude from the other analyses. As such, we adopt the fully-reduced Spitzer phase curve from Knutson et al. (2012) for subsequent analysis in both Spitzer channels.

Fig. 5 Spitzer phase curves of the eight exoplanets (gray, binned in black) fit with the hmℓ basis (red), plotted in order of equilibrium temperature (low to high from top left to bottom right). Each planet is plotted with the IRAC Channel 1 phase curve in the first and third row and the Channel 2 phase curve in the second and fourth row (see row labels on the right side). Maximum likelihood parameters for each fit are listed in Table 3. Residuals are plotted beneath each fit with a histogram of residuals on the left.

5.2 Analysis

In the previous section we determined that the h_mℓ basis can be reduced to three fitting parameters: the redistribution factor (or greenhouse parameter) ƒ or redistribution efficiency ε, the power in the first spherical harmonic term C₁₁, and the phase offset parameter Δϕ. In practice, these parameters respectively set: the baseline temperature field, the amplitude of the temperature variations from the day to the night side, and the hotspot offset.

Each phase curve can be fit using the h_mℓ basis with the free parameters: {C₁₁,Δϕ,ε}. We compute posterior distributions on each parameter with the No U-Turn Sampler using a STAN-like tuning schedule, implemented by (Foreman-Mackey et al. 2021). α and ω_drag are varied with tight Gaussian priors: α ~ 𝒩(0.6,0.1) and ω_drag ~ 𝒩(4.5, 1.5). By allowing these parameters to have non-zero variance, we are slightly exploring the degeneracies between parameters, especially between α and C₁₁ which are highly correlated in fits to the phase curve with ℓ_max = 1. We compute a fixed eclipse model for the planets using the best-fit planetary parameters from each paper in Table 2 using the Mandei & Agol (2002) model and implemented by Kreidberg (2015). The eclipse model incurs no new free parameters – we assume that all of the flux stemming from the phase curve signal vanishes during totality of secondary eclipse, so that the flux ratio F_p/F_⋆ = 0 in eclipse.

We allow the thermal phase curves to be contaminated by sinusoidal contributions from ellipsoidal variations and Doppler beaming with unknown amplitudes but fixed phases and periods. Typical amplitudes for these sinusoidal terms are small compared to the phase variations. The most significant ellipsoidal variations affected the phase curves of WASP-103 and WASP-19.

In Fig. 5, a typical draw from the posterior distribution is plotted in red for each of the Spitzer phase curves shown in gray, and binned in black. The upper row contains observations from the 3.6 μm Channel 1, and the lower row contains 4.5 μm Channel 2 observations. The phase curves are plotted in order of increasing equilibrium temperature from left to right. The fits in red are quite good approximations to the behavior of the light curves in black. A notable exception not shown here is WASP-12, which has been claimed to have a special phase curve shape at 4.5 μm due to co-orbiting material between the planet and the star (Bell et al. 2019).

One notable phenomenon in Fig. 5 is the variation in phase curve amplitude between Spitzer channels, particularly for WASP-19 b. The difference in amplitude and phase offset between the two channels is likely because different wavelengths probe different depths and therefore different temperatures in the planet’s atmosphere. We discuss this more in Sect. 6.

Posterior distributions and parameter correlations are shown for WASP-43 b in Fig. 6. The spherical harmonic power term C₁₁ is correlated with the dimensionless fluid parameter α, which has a tight prior.

Figure 7 shows the day side and night side temperatures derived from the Spitzer photometry of the hot Jupiters. We confirm that day side temperatures are proportional to the equilibrium temperatures as found by Schwartz & Cowan (2015). We also independently confirm the observation of Keating et al. (2019); Beatty et al. (2019) which suggested that the night side temperatures of most hot Jupiters are similar to 1100 K, irrespective of equilibrium temperature. It has been suggested that silicates and hydrocarbon hazes can produce this trend (Gao et al. 2020). We note that this trend is broken by KELT-9 b, which is hot enough to dissociate H₂ molecules on its dayside, affecting the global heat transport (see, e.g.: Bell & Cowan 2018; Wong et al. 2020; Roth et al. 2021; Mansfield et al. 2020).

The best-fit h_mℓ basis parameters C₁₁ and Δϕ for both the fits to the GCM temperature maps and the Spitzer phase curves are plotted in Fig. 8. The C₁₁ values from the fits to the hottest GCMs agree well with the results from the Spitzer phase curves. The Spitzer phase curves show a roughly linear increase with equilibrium temperature given by $$C_{11,\text{Ch}1}=(0.00019\pm 0.00008)T_{\text{eq}}+(0.02\pm 0.17)$$(23) $$C_{11,\text{Ch}2}=(0.00010\pm 0.00007)T_{\text{eq}}+(0.12\pm 0.14)$$(24)

for planetary equilibrium temperatures between 1200 K and 3100 K. The phase offset parameter Δϕ is close to zero for all Spitzer phase curves except for the two most extreme temperature cases, while all GCMs below 1500 K indicate significantly nonzero phase offsets. The rather extreme phase offset of HD 189733 b is consistent with the range observed in GCM temperature maps of similar equilibrium temperatures. The simulations and the observations are in good agreement in terms of both C₁₁ and Δϕ where the equilibrium temperatures overlap near 2000 K.

Fig. 6 Posterior distributions for the three derived parameters {Δϕ, C11,ε} which determine the amplitude and shape of the phase curve for WASP-43 b as observed in Spitzer IRAC 1 and 2 (3.6, 4.5 μm), with priors on {α, ωdrag}. The global α parameter is degenerate with the C11 specific to each bandpass, which introduces a cross-bandpass correlation between the two C11 parameters.

Fig. 7 Hemispherically-averaged day side and night side temperatures Td, Tn in Spitzer Channel 1 (red squares) and Channel 2 (black circles) as a function of planetary equilibrium temperature. These temperatures are derived by averaging the day side and night side temperature maps derived from the full phase curve fit. The hottest planet is KELT-9 b, which clearly breaks the apparent “uniform night side temperature” trend among cooler hot Jupiters noted by Keating et al. (2019). Uncertainties are smaller than the points.

Fig. 8 Maximum-likelihood parameters for the spherical harmonic power coefficient C11 (upper) and hotspot offset (lower) for the full sample of GCM fits and Spitzer phase curve fits in two filter bands. Blue points are the results of fits to GCMs, circles and squares correspond to inverted and non-inverted atmospheres. Red squares and black circles represent fits to the Spitzer phase curves in filters IRAC 1 and 2, respectively.

6 Discussion

6.1 Comparison With Bell et al. (2021)

Recently Bell et al. (2021) produced re-reduced Spitzer phase curves of many targets, including seven planets in common with this analysis. We compare their results with the ones presented in this work in Fig. 11, showing excellent agreement on the day side temperatures, indicating that the day side temperature inferences are not strongly model dependent. The two hottest planets show significantly hotter temperatures in the Bell et al. (2021) analysis than in the analysis presented here, which may be attributed to differences in the reduction process.

6.2 Temperature Forecasts For hot Jupiters

By tightly constraining the α and ω_drag parameters we have essentially fit the phase curves of each planet with a nearly identical temperature map with different relative scalings and rotations. Figure 8 shows roughly monotonic trends for the retrieved values of Δϕ and C₁₁ as a function of the equilibrium temperature for our analysis of both GCMs and Spitzer data. These trends suggest a well-behaved sequence of Spitzer -derived temperature maps, which we show in Fig. 9. This sequence in turn suggests, to first order, a roughly universal climate state for hot Jupiters, consistent with theoretical expectations (Showman & Polvani 2011; Tsai et al. 2014).

In other words, one unified map can explain the phase variations of all of the hot Jupiters, if we rescaled all of the derived maps by the three parameters (ƒ, Δϕ and C₁₁). We propose that most hot Jupiters should have temperature maps like Fig. 10. We can test this hypothesis with upcoming observations from the CHEOPS space telescope, which will measure a precise and time-resolved phase curve of KELT-9 b in the optical, among other planets. Deviations from this “mean hot Jupiter” temperature map could be indicative of variations in the atmospheric structure with depth, and time-varying atmospheric conditions, for example.

6.3 When To Use h_mℓ to fit a phase curve

The h_mℓ basis is a generalized spherical harmonic expansion representation (parabolic cylinder functions) of the temperature map of a celestial body based on solutions to the shallow water system. In this work, we have shown that spherical harmonic degrees as small as ℓ_max = 1 are sufficient to approximate the temperature maps of highly irradiated, giant exoplanets, for certain values of α and ω_drag.

The shallow water framework is applicable for systems where the horizontal dimension of dynamics in the fluid far exceeds the radial one (hence the “shallow” part of the name). This is true for giant exoplanets, for example, where a few pressure scale heights in the outer parts of the atmosphere interact with light, and the bulk of the planetary radius is insulated from stellar irradiation.

The shallow water system at small spherical harmonic degrees would struggle to approximate the temperature maps of planets without atmospheres. Such worlds can have T_d ≫ 0 K and T_n ≈ 0 K when synchronously rotating, and the discontinuity near the terminator could only be successfully reproduced by the h_mℓ basis at large spherical harmonic degrees. We note that this challenge is not unique to the h_mℓ basis, but rather a shortcoming of spherical harmonic expansion maps in general. As such, we recommend that the h_mℓ basis should be applied only to planets that are likely to have atmospheres.

Fig. 9 Maximum-likelihood temperature maps for several of the hottest exoplanets, in order of increasing equilibrium temperature from left to right, where the horizontal axes are planetary longitude and vertical axes are latitude. The horizontal axis is the longitude centered on the sub stellar point, the vertical axis is latitude from pole to pole. The upper row shows the Spitzer IRAC Channel 1 (3.6 μm) inferences, the lower row shows IRAC Channel 2 (4.5 μm). The longitudinal temperature variations are constrained by the Spitzer phase curve observations while the latitudinal temperature variations are fixed at the parameterization that best describes the GCMs.

6.4 Discrepancies Between Bandpasses

Differences between the inferred parameters across multiple filter band passes can be attributed to variations in the typical altitude probed by the photosphere of the planet at each wavelength. Table 3 shows the inferred phase curve parameters using the h_mℓ basis for Spitzer IRAC Channels 1 and 2 for eight planets. The day side and night side temperatures, for example, can be significantly different between the two Spitzer band-passes. This is the case for WASP-18 b, which has a day-side integrated temperature nearly 250 K cooler in Channel 1 than Channel 2, corresponding to a difference of ≈ 10× the uncertainty on either measurement. Such differences between band passes can occur when variable opacity with wavelength determines the depth of the photosphere. The 3D nature of the planet’s atmosphere is often described with a piecewise, non-linear temperature-pressure profile with altitude (see, e.g.: Kitzmann et al. 2020), and the phase curve brightness temperatures presented here represent samples from the T–P profile weighted by the contribution function for a given bandpass.

Fig. 10 The “universal” temperature map for hot Jupiters: scaled versions of this temperature map are consistent with both the temperature distributions in GCMs (which constrain the latitudinal distribution) and the phase curves of the hottest exoplanets (which constrain the longitudinal distribution). This map has α = 0.6 and ωdrag = 4.5 and an arbitrary scaling for ƒ, C11 and Δϕ = 0.

Fig. 11 Hemispherically-averaged day side and night side temperatures from analysis with the hmℓ framework compared with a more traditional analysis by Bell et al. (2021), showing good agreement between techniques.

6.5 Reporting Bond Albedos and Heat Redistribution From 2D Temperature Maps

As discussed in Sect. 3, Cowan & Agol (2011b) devised a simple parameterization for the Bond albedo A_B and heat redistribution efficiency parameter ε. This “0D” approach constructs a planet with a day side and night side temperature, and the relationship between those temperatures is dictated by the energy balance in the star-planet system, which is parameterized by A_B and ε. We show in Sect. 3 that the Cowan & Agol (2011b) relations for these parameters can be manipulated to compute A_B and ε given the temperature map of the planet.

Since the temperature maps inferred in this work are 2D, we outline a more general relationship between the full temperature map, A_B, and ε in Sect. 3.4. This 2D, energy-balanced approach sets the thermal emission of a planet through its entire surface equal to the radiation incident on the planet’s day side. We adopt the resulting expression for A_B throughout this work.

We also redefine the heat redistribution efficiency as the ratio of the flux emitted by the night side of the planet and the flux emitted through the day side. This definition preserves the limits ε ∈ [0,1], as in Cowan & Agol (2011b), and is straightforward to compute from both GCM outputs and h_mℓ maps. The Bond albedos and heat redistribution parameters derived in this work are plotted in Fig. 12.

Fig. 12 Posterior distributions for the Bond albedos and heat redistribution efficiencies from the Spitzer phase curves of several hot exoplanets, colored by their equilibrium temperatures, in the style of Fig. 3 of Schwartz et al. (2017). Section 3 defines AB and ε for the 2D temperature maps inferred from the hmℓ basis.

6.5.1 On The Definition Of The Bond Albedo

The Bond albedo is defined as the spherical albedo A_s(λ) integrated over all wavelengths, weighted by the stellar spectrum $$A_B=\frac{{\displaystyle \int I_{\star }}{A}_s\text{d}\lambda }{{\displaystyle \int I_{\star }}\text{d}\lambda },$$(25)

where I_⋆(λ, T_⋆) is the incoming stellar radiation. In Sect. 3.4, we express the Bond albedo assuming the brightness temperature of the planet in any finite bandpass will balance the incoming stellar radiation over all wavelengths; in other words, the planet is assumed to be a blackbody⁴ which reflects light with some efficiency 1 – A_B. A real planet does not emit as a black body and has spectral features, and thus its brightness temperatures in finite band passes are not equivalent to the equilibrium temperature. As a result of these subtleties, (1) the A_B inferred from a finite bandpass may be negative, and (2) A_B estimates inferred from brightness temperatures in distinct band passes may have different values⁵.

Fig. 13 Fraction of light reflected from the planetary surface at secondary eclipse normalized by the total flux from both reflection and thermal emission, assuming constant Ag(λ) = 0.1, after Fig. 4 of Schwartz & Cowan (2015). Black body emission is assumed for both the planet and the star. See Appendix B for details.

6.5.2 On The Definition of The Heat Redistribution Efficiency

The heat redistribution efficiency ε as defined for 2D temperature maps in Eq. (22) is distinct from the definition in Cowan & Agol (2011b). For example, consider a planet that has a uniform dayside temperature T_d and uniform nightside temperature T_n. The 0D derivation yields $${\epsilon }_{0\text{D}}=\frac{8T_{\text{n}}^4}{3T_{\text{d}}^4+5T_{\text{n}}^4},$$(26)

whereas the 2D approach reduces to $${\epsilon }_{2\text{D}}=\frac{T_{\text{n}}^4}{T_{\text{d}}^4}.$$(27)

These expressions are equivalent when ε_0D = ε_2D = 1 and T_d = T_n. For imperfect recirculation efficiencies, the appearance of the day side and night side temperatures in the denominator of Eq. 26 ensures that the 0D and 2D approaches yield different efficiencies. It can be shown that the 0D efficiencies are larger than the 2D efficiencies by up to ≲ 0.24, where the maximal difference occurs for T_d + T_n = 8(T_d – T_n).

6.6 Reflected Light

One caveat of the analysis presented here is that we make the usual assumption that the reflected light component of the phase curve in the Spitzer bandpass is negligible. Figure 13 shows the expected contribution of reflected light to the observed, reflected and thermal emission spectrum of the planet, assuming A_g =0.1, similar to typical hot Jupiters in the Kepler bandpass (Heng & Demory 2013). Details of the calculation are provided in Appendix B.

The contribution of reflected light to the total planetary flux is on the order of a few percent for the warm Spitzer band passes for all targets in this work. Reflected light contributes significant flux at optical wavelengths for cooler exoplanets such as HD 189733 b and WASP-43 b. Joint inferences are required for the thermal emission and reflectance properties of cooler exoplanets at optical wavelengths (see, e.g.: Heng et al. 2021). Ultra-hot Jupiters such as KELT-9 b are dominated by thermal emission even at blue-optical wavelengths. Phase curves of cooler giant exoplanets from the JWST in the NIRSpec/Prism mode, for example, will be a mixture of thermal emission and reflected light that varies with wavelength. A similar conclusion on reflected light was reached by Keating & Cowan (2017).

6.7 Comparison With Traditional Spherical Harmonics

include in Appendix C calculations which indicate that the h_mℓ basis reproduces GCM temperature maps more accurately than the traditional spherical harmonic $Y_{\mathcal{l}}^m$ basis. We conclude that the h_mℓ basis functions provide better fits to GCM temperature maps than similar $Y_{\mathcal{l}}^m$ basis functions, which may be relevant when self-consistently computing hotspot offsets and phase curve amplitudes directly from GCM temperature maps and indirectly from phase curves, as is the subject of this work.

7 Conclusion

• GCM temperature maps and real photometric phase curves can be compared using the same language, which is a set of modified spherical harmonic basis functions which we call the h_mℓ basis (Sect. 2)

• The “0D” derivation for the Bond albedo A_B and the heat redistribution efficiency ε popularized by Cowan & Agol (2011b) can be generalized when considering 2D temperature maps (Sect. 3)

• There are only three free thermal parameters required to fit each phase curve in a given bandpass, (C₁₁, Δϕ, ƒ}

• We can fit GCM outputs (Sect. 4) and Spitzer phase curves (Sect. 5) with the same {α, ω_drag}, and therefore there exists one “universal” hot Jupiter map, which after scaling and rotation can matches existing hot Jupiter observations

• Bond albedos, heat redistribution efficiencies, phase offsets, and day- and night-side temperatures are computed for several hot Jupiters with Spitzer phase curves (Table 3)

• We confirm that reflected light from planets will contaminate the bluest wavelengths accessible to the James Webb Space Telescope, and the contamination is most severe for cooler planets

• Comparisons of the h_mℓ basis with the traditional $Y_{\mathcal{l}}^m$ basis suggest that GCMs are more accurately reproduced with the h_mℓ basis, though the resulting phase curves are usually indistinguishable.

• We offer code that computes the h_mℓ maps and their phase curves efficiently called kelp, which is publicly available with implementations in Cython, theano and JAX⁶.

Acknowledgements.

This manuscript was carefully and thoughtfully reviewed and improved by Nicolas Cowan. This work has been carried out in the framework of the PlanetS National Centre of Competence in Research (NCCR) supported by the Swiss National Science Foundation (SNSF). We gratefully acknowledge software testing and discussions with Matthew Hooton and Adrien Deline. KH acknowledges partial financial support by a European Research Council Consolidator Grant (number 771620, project EXOKLEIN). KJ acknowledges financial support by a Hans Sigrist Foundation Ph.D Fellowship. We thank Russell Deitrick, Elsie Lee and Nathan Mayne for providing GCM output from previous publications. This research has made use of NASA’s Astrophysics Data System. This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA Explorer Program. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). The original description of the VizieR service was published in A&AS 143, 23. CP acknowledges financial support from a scholarship awarded by the Fonds de Recherche Québécois – Nature et Technologie (FRQNT B1X and International Internships Program – Programme de stages internationaux; Québec), a Technologies for Exo-Planetary Science (TEPS) PhD scholarship and the Vanier scholarship from the Natural Sciences and Engineering Research Council (NSERC) of Canada. We gratefully acknowledge the open source software that made this work possible: astropy (Astropy Collaboration et al. 2013, 2018), ipython (Pérez & Granger 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), batman (Kreidberg 2015), emcee (Foreman-Mackey et al. 2013), PyMC3 (Salvatier et al. 2016), theano (Theano Development Team 2016), pymc3-ext (Foreman-Mackey et al. 2021), corner (Foreman-Mackey 2016), cython (Behnel et al. 2011), JAX (Bradbury et al. 2018).

Appendix A h_mℓ linear solve

In this Appendix we outline how to reproduce a given GCM temperature map, for example, with the h_mℓ basis, as done in the bottom row of Figure 3. We outline here the linear solution for the C_mℓ terms that produce h_mℓ maps at arbitrary ℓ_max, following the linear solve technique of Luger et al. (2019).

For all m, ℓ, one can construct the h_mℓ perturbation map from Equation (2), setting each C_mℓ to a positive constant. Each 2D perturbation map is then flattened into a column vector, and each column vector is stacked with the others to produce a design matrix; this is known in the starry framework as a “pixelization matrix,” P. We stack onto this design matrix a column vector of ones to represent the mean temperature. One can then linearly solve for Q in $$AQ=P^{\top }w,$$(A.1)

where w is a weight matrix which can be assigned to sin θ in order to reproduce equatorial features and down-weight the polar regions, for example, and the matrix A is $$A=P^{\top }{w}^{2}P+\in I,$$(A.2)

where I is the identity matrix and ϵ is a small positive constant. The vector of unnormalized C_mℓ coefficients c′ which reproduce a given temperature map T_input can be found by multiplying Q by the temperature map, $$c^{\prime }=Q{T}_{\text{input}}.$$(A.3)

Finally we normalize each c′ by the mean of the temperature map $c^{\prime }=c^{\prime }/{\overline{\text{T}}}_{\text{input}}$ The resulting c vector represents each C_m,ℓ power coefficient to be passed through the h_mℓ framework in order to produce the h_mℓ representation of T_input.

Appendix B Reflected light

The observed flux from an exoplanet is the sum of the thermal emission spectrum and the reflected starlight, $$F_p(\lambda ,\xi )=F_{\text{refl}}(\lambda ,\xi )+F_{\text{therm}}(\lambda ,\xi ).$$(B.1)

The thermal emission component is written in Equation (3). The reflected light component is defined as a function of the geometric albedo at a given wavelength $$F_{\text{refl}}(\lambda ,\xi )=A_{\text{g}}(\lambda ){\left(\frac{R_p}{a}\right)}^2\Psi (\xi )F_{\ast }(\lambda ),$$(B.2)

where Ψ(ξ) is the integral phase function and F_⋆(λ) is the stellar emission spectrum.

The ratio of the reflected light spectrum observed at secondary eclipse (ξ = 0) to the total planetary flux is given by $$\frac{F_{\text{refl}}}{F_{\text{p}}}(\lambda ,0)={\left(1+\frac{F_{\text{therm}}(\lambda ,R_{p}0)}{A_g(\lambda ){\left(\frac{R_p}{a}\right)}^2{F}_{\ast }(\lambda ),}\right)}^{-1},$$(B.3) $$={\left(1+{\left(\frac{a}{R_{\ast }}\right)}^2{\displaystyle \iint \frac{{\mathcal{B}}_{\lambda }(T_{\text{p}})\cos (\varphi )\sin \text{d}\varphi \text{d}\theta }{\pi A_g(\lambda ){\mathcal{B}}_{\lambda }(T_{\ast })}}\right)}^{-1}.$$(B.4)

In Figure 13, we compute each curve by solving the equation above with the corresponding planetary and stellar parameters, assuming A_g(λ) = 0.1 is constant with wavelength. Note that the thermal emission component is computed using the temperature map inferred from the Spitzer Channel 2 phase curves, though as discussed in Section 6.4, the observed, brightness temperature of the planet is a function of wavelength.

Appendix C Comparison with traditional spherical harmonics

In Figure C.1, we select a model with three free parameters in both the h_mℓ and $Y_{\mathcal{l}}^m$ representations to solve linearly (see our Appendix A for the mathematical approach). The differences between h_mℓ and $Y_{\mathcal{l}}^m$ can be as large as 400 K peak-to-peak. The χ² between both the h_mℓ and $Y_{\mathcal{l}}^m$ representations compared with the GCM temperature map is labelled on each h_mℓ and $Y_{\mathcal{l}}^m$ map. We note that the χ² is significantly smaller for each h_mℓ representation than the corresponding $Y_{\mathcal{l}}^m$ representation with the same number of free parameters. This indicates that the GCM temperature maps are more accurately reproduced by the h_mℓ basis than the $Y_{\mathcal{l}}^m$ basis, even though the phase curves produced by either the h_mℓ or $Y_{\mathcal{l}}^m$ solutions are nearly indistinguishable.

Fig. C.1 The four hottest GCM temperature maps (first column), the best-fit linear solutions with a low-order hmℓ representation (second column), the fits with low-order ordinary spherical harmonic representation (third column, labeled $Y_{\mathcal{l}}^m$), the difference between the two representations (fourth column), and the phase curves produced by the two representations of the GCM map as observed in the CHEOPS bandpass (fifth column).

References

All Tables

All Figures

Fig. 1 Four montages of the first several terms in the generalized spherical harmonic expansion of the temperature map in the hmℓ basis. Each montage corresponds to a different combination of α and ωdrag, labeled in the upper right of each quadrant. Each map shows the temperature perturbation (purple to yellow is cold to hot) as a function of latitude and longitude (shown in Moll weide projections such that the sub stellar longitude is in the center of the plot). The m = 0 terms are always zero. The even-ℓ terms are asymmetric about the equator and therefore do not represent typical GCM results in equilibrium, so we keep all even-ℓ power coefficients fixed to zero in the subsequent fits (this approach assumes edge-on orbits, for exceptions see Cowan et al. 2017). In the upper two quadrants, when ωdrag is set to a smaller value, the chevron shape becomes more pronounced. In the two quadrants on the right, when the α is set to a larger value, the concentration of heat near the equator is more pronounced.

In the text

	Fig. 2 Lower limit on Td/Tirr (curve) for producing negative Bond albedos. For comparison, empirically obtained values of Td/Tirr from Wong et al. (2021) are shown; the uncertainties on these values (not shown for plot clarity) are between 2% and 5%.
In the text

Fig. 3 Comparison between the 2D temperature maps of hot Jupiter GCMs in the upper row (Perna et al. 2012), and the 2D temperature maps represented with the hmℓ basis in the middle and lower rows with ℓmax = 1 and 50, respectively. The horizontal axis is longitude, the vertical axis is latitude, and the center of each map is the substellar point. For model input parameters see Table 1; for an outline of the linear algorithm used to fit maps with ℓmax = 50, see Appendix A.

In the text

	Fig. 4 Posterior correlations between α and ωdrag from fits of the hmℓ basis maps with ℓmax= 1 to the four hottest GCM temperature maps. The posterior distributions overlap at a few standard deviations for the values we adopt in the remainder of this work, α = 0.6 and ωdrag = 4.5 (marked with black lines).
In the text

Fig. 5 Spitzer phase curves of the eight exoplanets (gray, binned in black) fit with the hmℓ basis (red), plotted in order of equilibrium temperature (low to high from top left to bottom right). Each planet is plotted with the IRAC Channel 1 phase curve in the first and third row and the Channel 2 phase curve in the second and fourth row (see row labels on the right side). Maximum likelihood parameters for each fit are listed in Table 3. Residuals are plotted beneath each fit with a histogram of residuals on the left.

In the text

	Fig. 6 Posterior distributions for the three derived parameters {Δϕ, C11,ε} which determine the amplitude and shape of the phase curve for WASP-43 b as observed in Spitzer IRAC 1 and 2 (3.6, 4.5 μm), with priors on {α, ωdrag}. The global α parameter is degenerate with the C11 specific to each bandpass, which introduces a cross-bandpass correlation between the two C11 parameters.
In the text

Fig. 7 Hemispherically-averaged day side and night side temperatures Td, Tn in Spitzer Channel 1 (red squares) and Channel 2 (black circles) as a function of planetary equilibrium temperature. These temperatures are derived by averaging the day side and night side temperature maps derived from the full phase curve fit. The hottest planet is KELT-9 b, which clearly breaks the apparent “uniform night side temperature” trend among cooler hot Jupiters noted by Keating et al. (2019). Uncertainties are smaller than the points.

In the text

Fig. 8 Maximum-likelihood parameters for the spherical harmonic power coefficient C11 (upper) and hotspot offset (lower) for the full sample of GCM fits and Spitzer phase curve fits in two filter bands. Blue points are the results of fits to GCMs, circles and squares correspond to inverted and non-inverted atmospheres. Red squares and black circles represent fits to the Spitzer phase curves in filters IRAC 1 and 2, respectively.

In the text

Fig. 9 Maximum-likelihood temperature maps for several of the hottest exoplanets, in order of increasing equilibrium temperature from left to right, where the horizontal axes are planetary longitude and vertical axes are latitude. The horizontal axis is the longitude centered on the sub stellar point, the vertical axis is latitude from pole to pole. The upper row shows the Spitzer IRAC Channel 1 (3.6 μm) inferences, the lower row shows IRAC Channel 2 (4.5 μm). The longitudinal temperature variations are constrained by the Spitzer phase curve observations while the latitudinal temperature variations are fixed at the parameterization that best describes the GCMs.

In the text

	Fig. 10 The “universal” temperature map for hot Jupiters: scaled versions of this temperature map are consistent with both the temperature distributions in GCMs (which constrain the latitudinal distribution) and the phase curves of the hottest exoplanets (which constrain the longitudinal distribution). This map has α = 0.6 and ωdrag = 4.5 and an arbitrary scaling for ƒ, C11 and Δϕ = 0.
In the text

	Fig. 11 Hemispherically-averaged day side and night side temperatures from analysis with the hmℓ framework compared with a more traditional analysis by Bell et al. (2021), showing good agreement between techniques.
In the text

	Fig. 12 Posterior distributions for the Bond albedos and heat redistribution efficiencies from the Spitzer phase curves of several hot exoplanets, colored by their equilibrium temperatures, in the style of Fig. 3 of Schwartz et al. (2017). Section 3 defines AB and ε for the 2D temperature maps inferred from the hmℓ basis.
In the text

	Fig. 13 Fraction of light reflected from the planetary surface at secondary eclipse normalized by the total flux from both reflection and thermal emission, assuming constant Ag(λ) = 0.1, after Fig. 4 of Schwartz & Cowan (2015). Black body emission is assumed for both the planet and the star. See Appendix B for details.
In the text

Fig. C.1 The four hottest GCM temperature maps (first column), the best-fit linear solutions with a low-order hmℓ representation (second column), the fits with low-order ordinary spherical harmonic representation (third column, labeled $Y_{\mathcal{l}}^m$), the difference between the two representations (fourth column), and the phase curves produced by the two representations of the GCM map as observed in the CHEOPS bandpass (fifth column).

In the text