DECIPHERING THERMAL PHASE CURVES OF DRY, TIDALLY LOCKED TERRESTRIAL PLANETS

Data from the Kepler telescope indicate that ∼50%–100% of nearby cool stars host a rocky planet (Dressing & Charbonneau 2013; Morton & Swift 2014). If we were able to characterize even a fraction of these planets we could vastly expand our understanding of processes fundamental for terrestrial planets, including planet formation, atmospheric escape, photochemistry, and atmospheric dynamics. The observational best-case scenario is a transiting planet, whose orbit we happen to view edge-on, so that the planet periodically passes in front of and behind its star. Broadly speaking, to characterize such a planet we would want to determine its atmospheric composition, temperature structure, and atmospheric mass.¹ The composition reflects how the planet formed, how its atmosphere subsequently evolved (e.g., via degassing from the interior or atmospheric escape), and the chemical state of its atmosphere. The temperature structure indicates the dynamical regime of the atmosphere and, if retrievable down to the surface, whether the planet could be habitable. The atmospheric mass reflects the planet's atmospheric evolution and also determines its habitability (by controlling whether water can exist as a liquid).

The most mature techniques for characterizing transiting planets are transit spectroscopy, in which starlight is measured as it filters through a planet's atmosphere, and emission spectroscopy, in which a planet's thermal emission is measured just before the planet is occulted by its star. In theory, high-resolution transit and emission spectra both contain enough information to uniquely constrain atmospheric composition, temperature structure, atmospheric mass, and planetary mass (Madhusudhan & Seager 2009; Benneke & Seager 2012; Line et al. 2012; Lee et al. 2012; Wit & Seager 2013). In practice, it is difficult to comprehensively characterize even hot Jupiters with any single spectroscopic technique due to measurement error and observational degeneracies (e.g., Burrows 2013; Hansen et al. 2014; Griffith 2014).

It is therefore desirable to seek additional methods for characterizing terrestrial planets that complement high-resolution spectroscopy. One simple approach is to observe a planet's broadband thermal phase curve, which is the net infrared flux the planet emits as it orbits its star. Before the planet passes in front of its star we observe flux emitted from the planet's nightside, and just before the planet passes behind its star we observe flux from the planet's dayside. The resulting phase curve can then be used to infer five pieces of information: the planet's average thermal emission, the location of hot and cold spots and the flux emitted at the hot and cold spots (Cowan & Agol 2008). This technique has already been applied to hot Jupiters. For example, Knutson et al. (2007) were able to infer equatorial superrotation on HD 189733b from the fact that its hot spot is shifted eastward of the substellar point, consistent with the theoretical prediction of Showman & Guillot (2002). It will be more challenging to measure thermal phase curves of smaller and cooler planets, but it should be possible to perform such measurements using next-generation instruments like the James Webb Space Telescope (JWST; Deming et al. 2009).

Although thermal phase curve measurements of terrestrial planets will soon be technically feasible, more work is needed to determine how they can be fully exploited. A natural starting point is to assume that planets accessible to near-future observations will be tidally locked (even though planets could also be trapped in higher-order spin resonances; see Section 5). For a tidally locked planet the phase curve depends largely on the atmospheric redistribution of energy between dayside and nightside. Many researchers have therefore proposed using phase curves to characterize tidally locked planets (Seager & Deming 2009; Cowan & Agol 2011; Selsis et al. 2011; Menou 2012a; Yang et al. 2013; Mills & Abbot 2013; Perez-Becker & Showman 2013; Yang & Abbot 2014; Kataria et al. 2014). At the same time, these results have also shown that phase curves are sensitive to multiple atmospheric parameters, which makes them difficult to interpret. For example, a small phase curve amplitude is compatible with: (1) a massive atmosphere because thicker atmospheres transport heat more effectively (Selsis et al. 2011); (2) an atmosphere containing large amounts of H₂, which has a higher heat capacity than high mean-molecular-weight gases and therefore loses heat more slowly as air is advected to the nightside (Menou 2012a); (3) relatively weaker absorption of shortwave radiation, so that stellar energy is deposited at higher pressures before being reemitted to space (Burrows et al. 2010; Heng et al. 2011); and (4) a low magnetic drag in ionized atmospheres, which allows higher wind speeds and thus more efficient heat transport (Rauscher & Menou 2012).

In this paper we disentangle how different atmospheric parameters affect phase curves and show how the phase curve amplitude can be used to constrain atmospheric mass. We focus on the phase curve amplitude because we find that, for many terrestrial planets, hot/cold spot offsets will be small (see Section 3). We focus on atmospheric mass because it will be difficult to infer from either transit or emission spectroscopy, as can be seen from the following argument. Following Lecavelier des Etangs et al. (2008), the maximum pressure that can be probed in transit is $p_{{\rm max}} = 0.56 \times g/\kappa _{{\rm min}} \times \sqrt{H/(2 \pi a)}$, where g is the acceleration of gravity, κ_min is the opacity per unit mass in the most transparent part of the spectrum, H is the scale height, and a is the planetary radius.² If we assume that Rayleigh scattering dominates the transit spectrum up to 0.75 μm and that this is the most transparent part of the spectrum, then, for an N₂ atmosphere, κ_min ∼ 2.59 × 10⁻⁶ m² kg⁻¹ (Table 5.2, Pierrehumbert 2011). For an Earth analog with a = a_⊕, g = 10 m s⁻² and H = 8 km, we find that p_max = 0.3 bar. In reality it would be even harder to probe an atmosphere this deeply in transit due to clouds and hazes (Fortney 2005) or atmospheric refraction (Betremieux & Kaltenegger 2013; Misra et al. 2014b). Emission spectroscopy generally probes deeper into an atmosphere than transit spectroscopy. Atmospheric mass can then be constrained using the fact that pressure-broadening widens molecular absorption features at higher pressures. Nevertheless, pressure and molecular abundances are largely degenerate in their effect on emission spectra, which complicates the interpretation of emission spectra. For example, von Paris et al. (2013) estimate that emission spectroscopy of a cool Earth analog with a 1 bar atmosphere could place an upper bound on the surface pressure of about 5.6 bar. Obtaining the upper bound would require a low-resolution spectrum (λ/Δλ = 20) with a signal-to-noise ratio (S/N) of 10 (their Table 3). To estimate how much observation time this would require on JWST, we use the S/N maps in Belu et al. (2011) as a guideline. We estimate that, for a cool M dwarf at 5 pc, S/N = 10 at this spectral resolution would require ∼14 days of continuous monitoring.³ This amounts to observing roughly every eclipse the planet makes during JWST's five-year mission lifetime. The pressure-dependent formation of O₂ dimers offers another method for measuring atmospheric mass (Misra et al. 2014a). However, being able to detect the dimer spectral signature requires an atmosphere with O₂ concentrations similar to Earth's. The observation time necessary would again amount to ∼12 days of continuous observation, or almost all available transits over JWST's mission lifetime. Such long-term and detailed observations could be feasible for high-priority observation targets, but even in those cases it would be desirable to have an independent and less time-consuming way of estimating the atmospheric mass.

In the following sections, we first analyze the dynamical and radiative equations relevant for terrestrial planet atmospheres. We specifically consider atmospheres that are "dry" (i.e., condensation is negligible) and tidally locked. For such atmospheres we identify six nondimensional parameters that could influence phase curves (Section 2). Next, we use an idealized general circulation model (GCM) to numerically test which of the six parameters' phase curves actually are sensitive to (Section 3). Except for hot and rapidly rotating planets with optically thick atmospheres, we do not find significant hot/cold spot offsets. We therefore focus on how the phase curve amplitude can be used to constrain an atmosphere's properties. We find that, except for hot and rapidly rotating planets, the phase curve amplitude is mainly sensitive to two nondimensional parameters. We then show how the phase amplitude can be combined with information from transit or emission spectroscopy to constrain the atmospheric mass of a terrestrial planet (Section 4). We estimate that one measurement of a nearby super-Earth's phase amplitude with JWST, taken over half the planet's orbit, could constrain the atmospheric mass to within a factor of two.

We adopt a basic, yet comprehensive, model for the phase curve of a terrestrial planet. We assume the planet is tidally locked and in a circular orbit. This regime is particularly relevant for planets orbiting smaller main-sequence stars, that is, K and M dwarfs. At minimum, the phase curve of such a planet is set by atmospheric fluid dynamics, radiative transfer, and surface–atmosphere exchange of energy and momentum. As is standard for planetary atmospheres, we model the atmospheric fluid dynamics using the primitive equations (Vallis 2006). The primitive equations assume hydrostatic equilibrium and that horizontal length scales are much larger than vertical ones, both of which tend to be excellent approximations for large-scale motions. We focus on atmospheres cooler than 1000 K, for which magnetic effects should be negligible (Menou 2012b). We use bulk aerodynamic formulae for the surface exchanges of energy and momentum. We model the radiative transfer as two-band (shortwave and longwave) gray radiation. We neglect scattering and assume that longwave and shortwave opacities increase linearly with pressure. The linear dependency approximates the effects of pressure broadening and continuum absorption in a well-mixed atmosphere (Robinson & Catling 2014).

We assume that the thermodynamics are dry. This is a natural starting point for a theoretical investigation, but our results should apply to a wide range of actual atmospheres. First, we expect that many terrestrial planets will be dry because post-formation delivery of volatiles via planetesimals and comets is a stochastic process (Morbidelli et al. 2000). In addition, for planets hotter than Earth, volatiles can be lost via atmospheric escape (the so-called moist greenhouse; Kasting 1988). On tidally locked planets, volatiles can also become cold-trapped on the nightside (Leconte et al. 2013; Menou 2013). Moreover, the dry regime is a useful approximation even for atmospheres like Earth's with moderate amounts of a condensing substance (Schneider 2006). We therefore expect that insight gained in the dry regime will carry over to the moist case. For example, if the atmospheric dynamics were insensitive to one parameter in the dry regime (e.g., surface friction), this suggests that moist atmospheres could be similarly insensitive. Finally, observations will be able to control for cases in which our analysis no longer applies. For example, condensation and cloud formation would lead to anomalously high bond albedos and could also reverse the expected day–night phase curve pattern (Yang et al. 2013).

The equations of our assumed model, shown in Appendix A, contain twelve dimensional parameters. The parameters are: stellar constant L_*, planetary albedo α, rotation rate Ω, planetary radius a, surface gravity g, specific heat capacity c_p, specific gas constant R, shortwave and longwave opacities at some reference pressure κ_SW and κ_LW, the Stefan–Boltzmann constant σ_SB, surface pressure p_s, and surface drag coefficient C_D. For a single gas species the specific gas constant is R ≡ k_B/(m_pM), where k_B is the Boltzmann constant, m_p is the mass of a proton, and M is the molecular weight of a gas molecule. For multiple species in a well-mixed dry atmosphere one can similarly assign bulk values of R and c_p (Caballero 2014). The opacities κ_SW and κ_LW are defined at a reference pressure, p₀. The choice of p₀ is arbitrary and one could set it equal to p_s, so it does not provide an additional dimensional parameter. To better compare our choices of κ_SW and κ_LW with previous work, we keep p₀ and p_s distinct. We also note that L_*, α, and σ_SB are not independent. The product L_*(1 − α) only appears in the stellar forcing term of the radiative equations, and σ_SB only appears in the radiative equations (Appendix A). We account for the degeneracy among L_*, α, and σ_SB by defining a characteristic temperature, T_eq = [L_*(1 − α)/(4σ_SB)]^1/4, which is the equilibrium emission temperature of a spatially homogeneous planet. This reduces the number of dimensional parameters to ten.

Following Frierson (2005), we use the Buckingham–Pi theorem to express ten dimensional parameters measured in four different units (mass, length, time, and temperature) as only six nondimensional parameters (Buckingham 1914). There is no unique choice for these nondimensional parameters; we form them using characteristic scales that we consider most appropriate for relatively slowly rotating tidally locked atmospheres. Our choice of scales nevertheless leads to nondimensional parameters that are well known in the literature. As a characteristic velocity scale we choose the speed of gravity waves, c ∼ NH, where N is the Brunt–Väisälä frequency and H ≡ RT_eq/g is the scale height. Adjustment via gravity waves is key in setting the day–night temperature gradients, and hence phase curves, of relatively slowly rotating planets (Perez-Becker & Showman 2013; Showman et al. 2013). The Brunt–Väisälä frequency is given by N² = g/T(g/c_p + dT/dz). The lapse rate, dT/dz, is a priori unknown for any atmosphere. To place an upper bound on the velocity scale we assume an isothermal atmosphere, dT/dz ∼ 0, so $c_{{\rm wave}} = \sqrt{R/c_p} \times \sqrt{gH}$. This amounts to assuming that gravity waves are very fast (on the order of the speed of sound, $c_{{\rm sound}} = \sqrt{\gamma _{a} R T_{eq}} = \sqrt{\gamma _{a} g H}$, where γ_a is the adiabatic index; cf. Heng & Kopparla 2012). As a characteristic length scale we choose the planetary radius a. We note that another potential length scale is given by the equatorial Rossby deformation radius, $L_{Ro}=\sqrt{a c_{{\rm wave}}/(2\Omega)}$, which is the maximum distance that equatorial waves can travel poleward under the influence of rotation. For slowly rotating planets the Rossby radius exceeds the planetary radius, L_Ro > a, and equatorial waves can propagate planetwide. We estimate that our choice of a as the length scale is valid for planets with orbital period $\gtrsim \mathcal {O}(6)$ days.⁴ From conservation of mass, we choose a vertical velocity scale c_wave × p_s/a. The remaining scales and nondimensionalized equations are shown in Appendix A.

We arrive at the following six nondimensional parameters:

$$\begin{eqnarray*} &&\left(\frac{R}{c_p}, \frac{2 \Omega a }{c_{{\rm wave}}}, \frac{a}{c_{{\rm wave}}} \frac{g \sigma _{SB} T_{eq}^3}{c_p p_s},\frac{\tau _{SW}}{\tau _{LW}},\frac{\kappa _{LW} p_0}{g} \left[\frac{p_s}{p_0}\right]^2, \frac{C_D a g}{R T_{eq}} \right)\\ &&\quad \equiv \left(\frac{R}{c_p}, \frac{a^2}{L_{Ro}^2}, \frac{t_{{\rm wave}}}{t_{{\rm rad}}},\gamma,\tau _{LW}, C_D \frac{a}{H} \right). \end{eqnarray*}$$

The six parameters are related to physical processes as follows: the adiabatic coefficient R/c_p controls the lapse rate and is also identical to the ratio 2/(2 + n), where n is the degrees of freedom of a gas (Pierrehumbert 2011). The nondimensional Rossby radius $a^2/L_{Ro}^2$ governs the latitudinal extent over which equatorial waves can transport energy and momentum (Matsuno 1966; Showman & Polvani 2011; Leconte et al. 2013). We emphasize that instead of $a^2/L_{Ro}^2$ one could choose different scales and arrive at, for example, a Rossby number or a nondimensionalized Rhines scale (Showman et al. 2010, 2013). Which scale to choose depends on the processes under consideration, and our anticipation of wave adjustment processes naturally leads to $a^2/L_{Ro}^2$. Our results support our analysis, and we find that phase curves are largely insensitive to planetary rotation when $a^2/L_{Ro}^2 \lesssim 1$, that is, as long as waves propagate planetwide (Section 3). The ratio t_wave/t_rad compares the time it takes for waves to redistribute energy across the planet, t_wave ≡ a/c_wave, to the atmosphere's radiative cooling time, $t_{{\rm rad}} \equiv c_p p_s/(g \sigma _{SB} T_{eq}^3)$ (Showman et al. 2013; Perez-Becker & Showman 2013). The atmospheric shortwave and longwave optical thicknesses at the surface are τ_SW ≡ κ_SWp₀/g × (p_s/p₀)² and τ_LW ≡ κ_LWp₀/g × (p_s/p₀)², and their ratio is γ. We note that the precise forms of τ_SW and τ_LW depend on details such as pressure broadening and scattering (e.g., Pierrehumbert 2011; Heng et al. 2014). Our definition of γ is equal to the more commonly used ratio of shortwave to longwave opacities (e.g., Guillot 2010) when shortwave and longwave opacities have the same pressure dependency. The influence of surface friction and surface heating on the atmosphere is governed by C_Da/H.

Two atmospheres governed by the equations that we assume are guaranteed to be dynamically similar (identical dynamics in statistical equilibrium) if their six nondimensional parameters are identical (also see Section 5). We note that only the nondimensionalized dynamics will be similar; the physical values of, for example, temperature gradients or wind speeds could be quite different. We also note that dynamical similarity does not depend on how we nondimensionalize the equations, that is, our particular choice of characteristic scales and nondimensional parameters. Nondimensionalization therefore allows us to identify atmospheres that one might consider distinct based on their dimensional parameters, but that turn out to be dynamically similar.

We test this idea in an idealized GCM. The model is based on the GFDL Flexible Model System (Held & Suarez 1994) and was subsequently modified by Frierson et al. (2006) and Merlis & Schneider (2010). This model has already been used to simulate the atmospheres of Earth (Frierson et al. 2006), Jupiter (Liu & Schneider 2011), hot Jupiters (Heng et al. 2011), tidally locked terrestrial planets (Merlis & Schneider 2010; Mills & Abbot 2013), and non-synchronously rotating terrestrial planets (Kaspi & Showman 2014). For our simulations, we remove moisture and replace the model's convective parameterization with an instantaneous dry convection scheme (Manabe et al. 1965). We run all simulations for at least 1000 days with a spatial resolution of either 64 × 128 × 30 or 48 × 96 × 20 grid points (latitude×longitude×vertical, corresponding to T42 or T31 spectral resolution). We use model time steps between 30 and 1200 s. We vary the time step because hot atmospheres require smaller time steps for numerical stability, whereas colder atmospheres can be integrated using longer time steps but also take longer to reach equilibrium. We expect this behavior, given that $c_{{\rm wave}} \propto \sqrt{T}$, so that hotter atmospheres are more likely to violate the Courant–Friedrichs–Lewy criterion. We consider a simulation equilibrated once the global-averaged radiative imbalance between incoming stellar and outgoing longwave radiation has fallen below at least 1% of the incoming stellar radiation. We note that the GCM simulates additional higher-order physics, and therefore contains additional parameters, which we did not include in the derivation of the six nondimensional parameters. In particular, the model contains a full Monin–Obukhov surface boundary layer scheme which self-consistently computes the depth of the boundary layer, diffusion of surface fluxes, and surface drag. This means the drag coefficient, C_D, is computed by the model instead of being a fixed parameter.⁵ For example, for a neutrally stratified boundary layer C_D = [k_vk/log (z/z₀)]², where k_vk is the von Karman constant, z is the height of the lowest model layer and z₀ is the roughness length. Because z and z₀ only enter into this equation logarithmically we modify C_D by adjusting k_vk. Similarly, the GCM requires additional parameters for its numerical algorithms. For example, the momentum equations are implemented using numerical dissipation via horizontal ∇⁸ hyperdiffusion. The hyperdiffusivity is chosen to damp the smallest resolved scale on a timescale of 12 hr, which sets a dissipative timescale. The dynamical core also uses a Robert–Asselin time filter, which is controlled by another nondimensional parameter. Our assumption, which we test in Section 3, is that the equations described in Appendix A capture the most important physics simulated by the GCM.

Although the GCM uses standard latitude–longitude coordinates, we adopt a tidally locked coordinate system to present our numerical results. Tidally locked planets in relatively long-period orbits tend to exhibit a strong symmetry about the axis connecting the substellar and antistellar points. Figure 1 shows the surface temperature in our first reference simulation, which is a cool, slowly rotating, Earth-sized planet with (a, Ω, T_eq) = (a_⊕, 2π/[50 days], 283 K) (see Table 1). The surface temperature is symmetric because the reference simulation is in a slowly rotating dynamical regime, $a^2/L_{Ro}^2 \ll 1$. The symmetry is not perfect, but it captures the dominant spatial variability. We therefore define a tidally locked coordinate system with a tidally locked latitude θ_TL and longitude λ_TL, where θ_TL is the angle away from the terminator and λ_TL is the angle about the substellar point (see Figure 1(b), Appendix B).

Zoom In Zoom Out Reset image size

To compute phase curves we follow Cowan & Agol (2008) and assume an edge-on viewing geometry (see Appendix C). We normalize the disk-integrated fluxes by F_rock ≡ 2/3 × L_*(1 − α), which is the dayside-averaged observer-projected flux emitted by a planet without an atmosphere. A bare rock will therefore have a nondimensional phase curve, F/F_rock, that varies between zero and one. On the other hand, a planet with efficient heat transport will have a constant phase curve equal to⁶ $F/F_{{\rm rock}} = \sigma _{SB} T_{eq}^4/F_{{\rm rock}} = 3/8$. We also define the phase curve peak-to-trough amplitude as the normalized difference between the phase curve maximum and minimum, (F_max − F_min)/F_rock.

First, we test whether our model (with six nondimensional parameters) captures the main physics of dry, tidally locked atmospheres. We consider our model adequate if different GCM simulations produce identical climates when their nondimensional parameters are identical. Our reference case is a cool, slowly rotating, Earth-sized planet (Table 1). Figure 2 compares the reference case to simulations in which we change dimensional parameters, but keep the six nondimensional parameters fixed (parameter choices are shown in Table 1 a). We find that the nondimensional surface temperature, T_s/T_eq, differs less than 1% between the reference simulation and the simulations with fixed nondimensional parameters. The nondimensional meridional wind velocity, v/c_wave, in the upper troposphere is more variable, and differs by ≲ 3% over most of the model domain. The largest deviation in meridional wind is ∼30% near the antistellar point. The deviation partly arises because we project simulated wind speeds into a tidally locked coordinate system, which mixes the wind vector components. If we instead consider the total wind speed, $\sqrt{u^2+v^2}/c_{{\rm wave}}$, the deviation is <10%. Moreover, v → 0 at the antistellar point, requiring longer averaging periods, and wave breaking occurs on the nightside, creating small-scale structure and numerical dissipation. Both effects can lead to deviations from dynamical similarity. For comparison we show some simulations in which we vary the nondimensional parameters (colored curves in Figure 2). In these simulations surface temperature and wind velocities change up to 300% compared with the reference, which demonstrates that the dynamical similarity predicted by the nondimensionalization is not trivial. We conclude that the above six nondimensional parameters are sufficient to capture the most important dynamics of the idealized GCM simulations.

Zoom In Zoom Out Reset image size

Table 1. (a) Parameters for the Simulations in Figure 2, (b) Parameters for the Simulations in Figures 3 and 4(a)

Dimensional Parameters and Nondimensional Parameters
Reference:
a = a⊕, Ω = 2π/(50 days), Teq = 283K, $R=R_{N_2}$, $c_p=c_{p,N_2}$,	$(\frac{R}{c_p}, a^2/L_{Ro}^2, \frac{t_{{\rm wave}}}{t_{{\rm rad}}},\gamma,\tau _{LW}, C_D \frac{a}{H})$
ps = 1 bar, g = 10 m s−2, κLW = 10−4 m2 kg−1, κSW = 0	=(0.29, 0.12, 5.1 × 10−3, 0, 1, 1.4)
(a) In 16 simulations, we vary the dimensional parameters while keeping all nondimensional parameters fixed.
For comparison, in the bottom three simulations, we allow the nondimensional parameters to change.
(below: relative to reference)
(R, cp, CD) × 2, Ω × 21/2, ps × 2−3/2, κLW × 23	Same as reference
(R, cp, CD) × 2−1, Ω × 2−1/2, ps × 23/2, κLW × 2−3	⋅⋅⋅
(R, cp, g) × 2, Ω × 21/2, ps × 2−1/2, κLW × 22	⋅⋅⋅
(R, cp, g) × 2−1, Ω × 2−1/2, ps × 21/2, κLW × 2−2	⋅⋅⋅
(R, cp, g) × 5, Ω × 51/2, ps × 5−1/2, κLW × 52	⋅⋅⋅
(R, cp, g) × 5−1, Ω × 5−1/2, ps × 51/2, κLW × 5−2	⋅⋅⋅
(Ω, g, κLW) × 2, a × 2−1	⋅⋅⋅
(Ω, g, κLW) × 2−1, a × 21	⋅⋅⋅
$\Omega \times (\frac{5}{4})^{1/2}, (T_{eq},g) \times \frac{5}{4}, p_s \times (\frac{5}{4})^{7/2}, \kappa _{LW} \times (\frac{5}{4})^{6}$	⋅⋅⋅
$\Omega \times (\frac{4}{5})^{1/2}, (T_{eq},g) \times \frac{4}{5}, p_s \times (\frac{4}{5})^{7/2}, \kappa _{LW} \times (\frac{4}{5})^{6}$	⋅⋅⋅
$(a,R,c_p) \times \frac{3}{2}, (\Omega,p_s) \times (\frac{3}{2})^{-1/2}, \kappa _{LW} \times \frac{3}{2}$	⋅⋅⋅
(a, R, cp) × 2, (Ω, ps) × 2−1/2, κLW × 2	⋅⋅⋅
$(a,R,c_p) \times \frac{2}{3}, (\Omega,p_s) \times (\frac{2}{3})^{-1/2}, \kappa _{LW} \times \frac{2}{3}$	⋅⋅⋅
$(a,R,c_p) \times \frac{1}{2}, (\Omega,p_s) \times (\frac{1}{2})^{-1/2}, \kappa _{LW} \times \frac{1}{2}$	⋅⋅⋅
(R, cp, CD) × 5, Ω × 51/2, ps × 5−3/2, κLW × 53	⋅⋅⋅
$(R,c_p,C_D) \times \frac{1}{5}, \Omega \times (\frac{1}{5})^{1/2}, p_s\times (\frac{1}{5})^{-3/2}, \kappa _{LW} \times (\frac{1}{5})^3$	⋅⋅⋅
κLW × 2	τLW × 2
Ω = 2π/(2 days)	$a^2/L_{Ro}^2\times 25$
ps/10, κLW × 102	twave/trad × 10
(b) This illustrates how we vary one nondimensional parameter at a time while keeping other nondimensional
parameters fixed (also see Appendix D). The dimensional parameters remain within the constraints shown in Table 2.
(below: relative to reference)
cp × 1.5, (R, g) × 1.51/2, ps × 1.5−1/2, κLW × 1.53/2	R/cp × 0.82
(R, cp, CD) × 0.75, ps × 0.75−3/2, κLW × 0.753	$a^2/L_{Ro}^2 \times 1.15$
(R, cp, CD) × 4.64, ps × 4.64−3/2, κLW × 100	$a^2/L_{Ro}^2 \times 0.46$
ps × 6.32−1, g × 2.5, CD × 2.5−1, κLW × 100	twave/trad × 15.8
ps × 0.1, κLW × 100	twave/trad × 10
ps × 0.141, κLW × 50	twave/trad × 7.1
ps × 0.5, κLW × 4	twave/trad × 2
ps × 21/2, κLW × 0.5	twave/trad × 0.71
ps × 3.162, κLW × 0.1	twave/trad × 0.32
ps × 2, g × 2.5−1, CD × 2.5, κLW × 0.1	twave/trad × 0.2
κSW = 5 × 10−5	γ = 0.5
κSW = 10−4	γ = 1
κLW × 100	τLW × 100
κLW × 50	τLW × 50
κLW × 10	τLW × 10
κLW × 0.5	τLW × 0.5
κLW × 0.2	τLW × 0.2
κLW × 0.1	τLW × 0.1
CD × 10, (g, ps) × 2.5, κLW × 0.4	CDa/H × 25
CD × 10	CDa/H × 10
CD × 0.1	CDa/H × 0.1
CD × 0.1, (g, ps) × 2.5−1, κLW × 2.5	CDa/H × 0.04

Notes. The drag coefficient C_D is not a fixed parameter in the model, so the values shown for C_D(a/H) is only approximate. Symbols for dimensional parameters are defined in Table 2, nondimensional parameters are defined in Section 2.

Download table as:  ASCIITypeset image

Next, we explore how sensitive phase curves are to each of the nondimensional parameters. We consider different reference simulations and vary their nondimensional parameters one at a time to see how this affects the resulting phase curves. For the reference simulations we consider different scenarios where a, Ω, and T_eq are fixed. We do so because a, Ω, and T_eq are relatively easily constrained for a transiting planet. To vary one nondimensional parameter at a time we first find all possible transformations of the dimensional parameters that only modify a given nondimensional parameter (see Appendix D). For example, to decrease t_wave/t_rad we could increase p_s but adjust κ_SW and κ_LW such that γ and τ_LW remain constant. We then use these transformations to vary each nondimensional parameter over its largest range compatible with fixed (a, Ω, T_eq) and the constraints in Table 2.

Table 2. Maximal Range of Dimensional Values we Consider

Dimensional Parameter	Symbol	Unit	Minimum Value	Maximum Value
Planetary radius	a	a⊕	1	2
Rotation rate	Ω	days−1	2π/50	2π/2
Equilibrium temperature	Teq	K	100	600
Surface gravity	g	10 m s−2	$\frac{2}{5} \times (a/a_{\oplus})$	$\frac{5}{2} \times (a/a_{\oplus})$
Specific heat capacity	cp	J kg−1 K−1	820	14230
Specific gas constant	R	J kg−1 K−1	190	4157
Surface pressure	ps	bar	10−2	10
Longwave opacitya	κLW	m2 kg−1	10−5	10−2
Shortwave opacitya	κSW	m2 kg−1	0	10−2
Surface drag coefficient	CD, via kvk	–	× 0.1	× 10

Notes. For a given reference simulation, we fix (a, Ω, T_eq) and change the remaining dimensional parameters such that only one nondimensional parameter varies at a time. C_D is not a fixed parameter, so we vary the von-Karman constant k_vk to increase and decrease C_D by an order of magnitude. We vary R and c_p, but require that R/c_p stays within the range of diatomic and triatomic gases (0.22 ≲ R/c_p ⩽ 0.29; Section 2.3.3., Pierrehumbert 2011). The maximum value of R and c_p corresponds to H₂, the minimum value to CO₂. In addition, we require that shortwave optical depth does not exceed longwave optical depth (γ ⩽ 1). ^aOpacities are defined at a reference pressure of p₀ = 1 bar.

Download table as:  ASCIITypeset image

We find that, for most planets, only the phase curve peak-to-trough amplitude is robustly sensitive to changes in the nondimensional parameters. We start with the above reference scenario of a cool, slowly rotating planet. The dashed curve in Figure 3 shows the phase curve of the reference simulation. The planet's thermal flux is phase-locked with the incoming stellar radiation, that is, there is no hot spot phase offset. The phase-locking arises because all stellar radiation is absorbed at the ground (γ = 0), while a significant part of this energy can also escape directly from the ground to space without being advected (τ_LW = 1). We note that the cold spot offset is larger than the hot spot offset; however, it would be difficult to detect the cold spot offset because the phase curve is approximately constant near the antistellar point (dashed curve in Figure 3). Next, we vary each nondimensional parameter while keeping the other nondimensional parameters fixed (Table 1b). The gray lines in Figure 3 show how the phase curve varies in response to changes in the nondimensional parameters. We find that the phase curve generally stays phase-locked with the stellar radiation. This only changes once the atmosphere becomes optically thick (τ_LW ≫ 1), but in those cases the offset would again be hard to detect because the planet's thermal emission essentially does not vary (approximately constant curves in Figure 3). In contrast to the negligible hot/cold spot offsets, the phase curve amplitude is much more sensitive to changes in the nondimensional parameters (Figure 3).

Zoom In Zoom Out Reset image size

We explore other reference simulations to see when hot/cold spot offsets become significant. We find that significant hot/cold spot offsets only occur when the atmosphere is optically thick, τ_LW ≫ 1, the planet has a high rotation rate, $a^2/L_{Ro}^2 \gtrsim 1$, and is relatively hot, t_wave/t_rad ≳ 0.01. Table 3 summarizes our results. While a cool and slowly rotating planet with τ_LW = 10 shows a hot spot offset of up to 82°, the phase curve in that case is almost constant and the offset therefore not detectable (second-to-bottom row in Table 3). Only in the hottest and most rapidly rotating scenario with τ_LW ≫ 1 that we consider do we find a large hot spot offset of 26° which would also be detectable (bottom row in Table 3). If future observations found a large hot spot offset, this would therefore not only imply that the planet has an atmosphere, but also that the atmosphere would have to be optically thick. Many terrestrial planets, however, should have small hot spot offsets (Table 3). For the rest of this paper we therefore focus on the phase amplitude, and how it could be used to characterize the atmosphere of a planet.

Table 3. We Explore a Broad Range of Atmospheric Scenarios

Reference Parameters						Reference	Ref. Hot	Amplitude Sensitivity to
$\frac{R}{c_p}$	$\frac{a^2}{L_{Ro}^2}$	$\frac{t_{{\rm wave}}}{t_{{\rm rad}}}$	γ	τLW	$C_D \frac{a}{H}$	Amplitude	Spot Offset	$\frac{R}{c_p}$	$\frac{a^2}{L_{Ro}^2}$	$\frac{t_{{\rm wave}}}{t_{{\rm rad}}}$	γ	τLW	$C_D \frac{a}{H}$
0.29	0.12	5.1 × 10−3	0	1	1.4	0.24*	0°	+0	+0	$\color{red}{+0.31}$	+0	$\color{red}{+0.49}$	$\color{red}{+0.1}$
								−0.02	−0	−0.08	$\color{red}{-0.11}$	$\color{red}{-0.24}$	−0
0.29	0.3	5.1 × 10−3	0	1	1.4	0.26	3°	+0	+0	$\color{red}{+0.32}$	+0	$\color{red}{+0.48}$	+0.09
								−0.02	−0.01	−0.07	$\color{red}{-0.11}$	$\color{red}{-0.25}$	−0
0.29	0.6	1 × 10−2	0	1	2.8	0.36	3°	+0	+0.02	$\color{red}{+0.34}$	+0	$\color{red}{+0.42}$	+0.06
								−0.03	−0.05	$\color{red}{-0.11}$	$\color{red}{-0.16}$	$\color{red}{-0.32}$	−0.02
0.29	4.1	6.6 × 10−2	0	1	1.3	0.8**	7°	+0	+0.01	$\color{red}{+0.16}$	+0	$\color{red}{+0.1}$	+0.09
								−0.03	$\color{red}{-0.16}$	$\color{red}{-0.34}$	$\color{red}{-0.13}$	$\color{red}{-0.51}$	$\color{red}{-0.1}$
0.29	10	7.6 × 10−4	0	1	8.1	0.19	0°	+0	+0	$\color{red}{+0.32}$	+0	$\color{red}{+0.48}$	+0.04
								−0.03	−0.03	−0	−0.07	−0.14	−0.03
0.29	0.12	5.1 × 10−3	0.5	10	1.4	0.01	$\color{red}{82^\circ}$	+0	+0	+0	+0	$\color{red}{+0.71}$	+0
								−0	−0	−0.01	−0.01	−0.01	−0
0.29	4.1	6.6 × 10−2	0.5	10	1.4	0.5	$\color{red}{26^\circ}$	+0	+0.04	$\color{red}{+0.13}$	+0.03	$\color{red}{+0.37}$	+0.06
								−0.05	$\color{red}{-0.17}$	$\color{red}{-0.31}$	−0	−0.01	−0.08

Notes. The six columns on the left show the nondimensional parameters for the reference simulations, the center two columns show the phase curve amplitude and the hot spot offset in the reference simulation, and the six columns on the right show the maximum/minimum change in phase curve amplitude in response to each nondimensional parameter. The red font in the column "Ref. Hot Spot Offset" emphasizes entries with values larger than 10°. The red font in the six rightmost columns emphasizes entries whose absolute value is ⩾0.1. The top row, *, corresponds to the cool slowly rotating scenario in Figures 3 and 4(a), and the fourth row from the top, **, corresponds to the hot and rapidly rotating scenario in Figure 4(b).

Download table as:  ASCIITypeset image

We find that, except for hot and rapidly rotating planets, the phase amplitude is primarily sensitive to the ratio of dynamical to radiative timescales, t_wave/t_rad, and the optical depth, τ_LW. We again start with the reference scenario of a cool, slowly rotating planet. The dashed line in Figure 4(a) shows the reference phase amplitude, (F_max − F_min)/F_rock, and the vertical lines show how sensitive the phase amplitude is to changes in each nondimensional parameter. The phase amplitude is far more sensitive to t_wave/t_rad and τ_LW than to any of the other nondimensional parameters. For example, when we vary t_wave/t_rad (primarily by changing g and p_s, while adjusting other parameters; see Table 1b), the phase amplitude varies between 0.2 and 0.6. In contrast, when we vary the nondimensional Rossby radius, $a^2/L_{Ro}^2$ (by changing R, and thus c_wave, while adjusting other parameters), the phase amplitude varies by less than 0.01. We emphasize this does not mean that the atmospheric dynamics or phase curve are insensitive to $a^2/L_{Ro}^2$ in general. As we show in the next paragraph, $a^2/L_{Ro}^2$ can affect phase curves when $a^2/L_{Ro}^2 \gtrsim 1$. However, for this particular scenario, once (a, Ω, T_eq) are known then $a^2/L_{Ro}^2$ is already constrained to be much smaller than one. The remaining observational uncertainty in $a^2/L_{Ro}^2$ barely affects our interpretation of the planet's phase curve amplitude. Figure 4(a) therefore shows that, for cool, slowly rotating planets with known (a, Ω, T_eq), a phase curve measurement contains essentially no information about the parameters R/c_p, $a^2/L_{Ro}^2$, γ, and C_Da/H. On the other hand, a measurement of the phase amplitude would constrain the combination of t_wave/t_rad and τ_LW.

Zoom In Zoom Out Reset image size

We explore other reference simulations to determine whether there are regimes in which the phase amplitude is sensitive to other parameters. Similar to our result for hot spot offsets, we find that phase amplitude only becomes sensitive to $a^2/L_{Ro}^2$, γ, and C_Da/H for large, hot, and rapidly rotating planets. Specifically, a planet has to have both $a^2/L_{Ro}^2 \,{\gtrsim}\, 1$ and t_wave/t_rad ≳ 0.01 for additional nondimensional parameters to affect the phase amplitude. Our results are summarized on the right-hand side of Table 3. We find that, in all scenarios, the phase amplitude is most sensitive to t_wave/t_rad and τ_LW. We also find that $a^2/L_{Ro}^2$ and t_wave/t_rad both have to be large for the phase amplitude to become sensitive to additional parameters; a large value of $a^2/L_{Ro}^2$ by itself is not sufficient (third-to-bottom row in Table 3). Together with our above result that hot spot offsets also require $a^2/L_{Ro}^2 \gtrsim 1$ and t_wave/t_rad ≳ 0.01, this suggests that a regime shift occurs in the atmospheric dynamics near this threshold. Figure 4(b) shows the scenario in which phase amplitude is most sensitive to additional parameters (** in Table 3). This scenario corresponds to a super-Earth with (a, Ω, T_eq) = (2a_⊕, 2π/[2 days], 600 K). In this scenario the phase amplitude is additionally sensitive to variations in $a^2/L_{Ro}^2$, γ, and C_Da/H (Figure 4(b)). For such a planet, a measurement of the phase amplitude would be degenerate with multiple atmospheric parameters, although the hot/cold spot offsets could provide additional information (Table 3). Many terrestrial planets, however, will have phase amplitudes that are, to good approximation, only sensitive to t_wave/t_rad and τ_LW (Table 3).

For planets whose phase curves primarily depend on the ratio of dynamical to radiative timescales, t_wave/t_rad, and the optical depth, τ_LW, we consider how a phase curve could constrain an atmosphere's properties. Expanded in terms of dimensional quantities, the two nondimensional parameters are $a g \sigma _{SB} T_{eq}^{5/2}/(p_s R c_p^{1/2})$ and $\kappa _{LW} p_s^2/(g p_0)$. The most important unknowns are the longwave opacity at a reference pressure, κ_LW, and the surface pressure, p_s, because the other dimensional parameters are relatively easy to constrain. For a transiting planet, the planetary radius, a, would be known. One can constrain the equilibrium temperature, T_eq, because a planet's broadband thermal emission averaged over one orbit is equal to $\sigma _{SB} T_{eq}^4$. The specific gas constant and heat capacity, R and c_p, vary most significantly between H₂-dominated atmospheres and high mean-molecular-weight atmospheres. Because R also sets the atmospheric scale height, a transit spectrum could be sufficient to distinguish between an H₂ and a high mean-molecular-weight atmosphere. If one can determine whether an atmosphere is H₂-dominated or not, the detailed value of R and c_p is secondary; for example, t_wave/t_rad only varies by a factor of two between a pure N₂ atmosphere and a pure CO₂ atmosphere. The surface gravity, g, can be constrained via radial-velocity or transit-timing measurements. Moreover, interior models indicate that bulk compositions ranging from water ice to iron would only change the bulk density, and thus g, by a factor of ∼2 (Seager et al. 2007). In contrast, κ_LW and p_s can change the values of t_wave/t_rad and τ_LW by several orders of magnitude.

This means a planet's phase amplitude can be used to characterize longwave opacity, κ_LW, and surface pressure, p_s. To evaluate the feasibility of doing so, we estimate the observable phase amplitude signal and the precision possible with JWST. For the signal we assume an optimistic scenario similar to that assumed by Yang et al. (2013). Specifically, we assume a cool super-Earth with $(a,T_{eq},\Omega,g,R,c_p) = (2a_\oplus, 300\ {\rm K}, 2\pi /[10 \ \mathrm{days}],20\, \mathrm{m\, s}^{-2}, R_{N_2}, c_{p,N_2})$, orbiting a GJ1214-like star with (a_*, T_*) = (0.2a_☉, 3000 K). We assume the star is 5 pc away, and the phase curve is observed between transit and secondary eclipse, for 5 days total. This planet would have $a^2/L_{Ro}^2 \sim 1$. That is roughly the regime for which our findings start to apply, i.e., the planet's phase amplitude largely depends on only two nondimensional parameters. Cooler and/or smaller planets would be even more solidly in the slowly rotating and cool regime, but more difficult to observe. We performed simulations that explore the phase amplitude as a function of κ_LW and p_s. For a given simulation we compute the normalized phase amplitude (F_max − F_min)/F_rock. We multiply this amplitude by the planet–star contrast of a bare rock, F_rock/F_*, to get the phase amplitude relative to the stellar flux, (F_max − F_min)/F_*. To compute F_rock/F_* we approximate the planetary and stellar emission as blackbody radiation. Following Yang et al. (2013), the planet–star contrast in some band [λ₁, λ₂] is then

$$\begin{eqnarray*} &&\frac{F_{{\rm rock}}}{F_*} = \left(\frac{a}{a_*} \right)^2 \frac{\int _{\lambda _1}^{\lambda _2} B(T_{{\rm rock}},\lambda) d\lambda }{\int _{\lambda _1}^{\lambda _2} B(T_*,\lambda) d\lambda }. \end{eqnarray*}$$

Here B is the Planck function, and T_rock is the dayside-averaged observer-projected temperature of a bare rock (Appendix C). We assume 16.5 ⩽ λ ⩽ 19.5 μm, which corresponds to the F1800W filter on JWST's Mid-Infrared Instrument. We choose this band because it avoids the 15 μm CO₂ absorption feature, but other spectral window regions would be similarly suitable. For a bare rock the phase amplitude would be (F_max − F_min)/F_* = (F_rock − 0)/F_* = 358 ppm, while an atmosphere with perfect day–night heat transport would have a phase amplitude of 0 ppm.

To estimate the precision possible with JWST, we assume that the observational error is dominated by stellar photon noise. In the photon noise limit the precision is $\sigma /F_* = 1/\sqrt{N}$, where the number of stellar photons N is

$$\begin{eqnarray*} &&N = \pi \left(\frac{D}{2} \right)^2 \Delta t \left(\frac{a_*}{d} \right)^2 \int _{\lambda _1}^{\lambda _2} \frac{B(T_*,\lambda)}{E(\lambda)} d\lambda. \end{eqnarray*}$$

Here D is the diameter of JWST's mirror (=6.5m), Δt is the length of observation, d is the distance between observer and star, and E(λ) = hc/λ is the energy per photon. We first make an optimistic estimate for the precision. Recent measurements with the Hubble Space Telescope almost reached the photon noise limit (Kreidberg et al. 2014; Knutson et al. 2014), and we assume JWST will do similarly well. To account for imperfect instrument throughput and detector efficiency we degrade the photon-limited precision by a factor of 1/3 (Figure 3, Glasse et al. 2010). We find that, over a 12 hr integration, JWST should be able to measure the planet–star flux ratio with a precision of 12 ppm.⁷ Because the phase amplitude is the difference between two fluxes, the 3σ uncertainty interval for the phase amplitude is⁸ $\sqrt{2}\times 3 \times 12$ ppm = 51 ppm. For a pessimistic estimate we repeat the previous calculation, but additionally impose a noise floor of 40 ppm. This floor represents unexpected instrumental systematics, zodiacal light or other noise sources. With this noise floor, a 12 hr integration would only reach 10% of the photon noise limit. Given that Spitzer measurements were able to reach ∼30% of the photon noise limit (Figure 3, Cowan et al. 2012), we consider this estimate very pessimistic. In this case a 2σ (3σ) measurement of the phase amplitude would have a precision of $\sqrt{2}\times 2\times 40= 113$ (170) ppm.

Figure 5 shows our simulation results and optimistic 1σ precision estimate. We find that thin atmospheres (p_s ⩽ 0.2 bar) have phase amplitudes close to 358 ppm for small and moderate values of κ_LW (⩽10⁻³ m² kg⁻¹). Taking into account measurement uncertainties of ∼50 ppm, these atmospheres would be difficult to distinguish from bare rocks. Similarly, thick atmospheres (p_s > 5 bar) tend to have phase amplitudes close to zero. A phase curve would constrain atmospheric mass most effectively between those two limits. Figure 5 also shows that, between those limits, the phase amplitude is sensitive to both κ_LW and p_s. We find that any observed phase amplitude would be compatible with both a thin atmosphere that is opaque to longwave radiation (large κ_LW) and a thick atmosphere that is transparent to longwave radiation (small κ_LW). Nevertheless, if transit or emission spectroscopy could determine the concentration of greenhouse gases in an atmosphere, and therefore κ_LW, the phase amplitude would yield the value of p_s. As an example we highlight a simulation with a 1 bar atmosphere and κ_LW = 4 × 10⁻⁴ m² kg⁻¹ (black square in Figure 5). The exact value of this atmosphere's phase amplitude is 152 ppm. The observed phase amplitude would therefore be 152 ± 51 ppm with 3σ confidence, which constrains the surface pressure to 0.7 ⩽ p_s ⩽ 1.3 bar. Even using our pessimistic precision estimate, we find that JWST would be able to constrain surface pressure to 0.5 ⩽ p_s ⩽ 2.3 bar, albeit only with 2σ confidence (152 ± 113 ppm). We further note that our pessimistic precision estimate would only place a lower bound on the surface pressure, p_s ⩾ 0.2 bar, with 3σ confidence (152 ± 170 ppm). These values are the most precise constraints that the phase curve amplitude can place on surface pressure, because we assumed the other dimensional parameters (κ_LW, R, c_p, etc.) are already well characterized via transit or emission spectroscopy. Observational uncertainties in the other dimensional parameters would increase the uncertainty in the inferred surface pressure. Nevertheless, our results show that a phase amplitude measurement can place meaningful bounds on a planet's atmospheric mass, while the necessary observation time is competitive with the time required to constrain atmospheric mass via transit or emission spectroscopy (see Section 1).

Zoom In Zoom Out Reset image size

Dimensional analysis is a crucial tool in comparative planetology and the study of exoplanets (Golitsyn 1970; Mitchell & Vallis 2010; Showman et al. 2010; Read 2011; Potter et al. 2013; Del Genio 2013; Mitchell et al. 2014). Our approach highlights the utility of the Buckingham–Pi theorem for the study of planetary atmospheres (cf. Frierson 2005). We show that the primitive equations coupled to the two-stream equations are governed by a fairly small set of nondimensional parameters. These nondimensional parameters also encapsulate the atmospheric dynamics of an idealized GCM. Our analysis reveals basic dimensional degeneracies, which could allow modelers to sample large parameter spaces more efficiently (Figure 2, Appendix D). It is straightforward to expand our analysis to include additional physics, for example, moist thermodynamics, multi-band radiation, non-hydrostatic atmospheres, chemical disequilibrium or magnetohydrodynamics. It follows from the Buckingham–Pi theorem that each additional independent physical parameter will introduce another nondimensional parameter.

Our analysis suggests that the dynamics of gaseous planets could be even easier to understand than the dynamics of terrestrial planets. Gaseous planets do not have a distinct surface, so we suppose that their atmospheric dynamics are to first order independent of the bottom boundary. A range of modeling studies tend to support this assumption (Heng et al. 2011; Menou 2012a; Kataria et al. 2013). This means the dynamics are insensitive to the surface friction/heating parameter, C_Da/H. We furthermore need to replace the surface pressure, p_s, with a new characteristic pressure. We note that for terrestrial planets p_s denotes the depth of the dynamically active part of the atmosphere, where winds are driven by gradients in the stellar forcing. For a gaseous planet we analogously use the photon deposition depth, that is, the pressure where stellar radiation is absorbed, $p_{D} \sim \sqrt{g p_0/\kappa _{SW}}$ (see Equation (117), Heng et al. 2014). The longwave optical depth τ_LW, which was previously defined at the surface, now becomes the optical depth at the level of photon deposition, τ_LW = κ_LWp₀/g × (p_D/p₀)² = κ_LW/κ_SW = γ⁻¹. The last step shows that τ_LW and γ cease to be independent degrees of freedom. This is also consistent with the Buckingham–Pi theorem because we should lose one nondimensional parameter in the limit p_s → ∞. The dynamics of gaseous planets then depend on only four nondimensional parameters (where p_s is now replaced by p_D):

$$\begin{eqnarray*} &&\left(\frac{R}{c_p}, \frac{a^2}{L_{Ro}^2}, \frac{t_{{\rm wave}}}{t_{{\rm rad}}},\gamma \right). \end{eqnarray*}$$

This set is sufficiently small for easy numerical exploration. Moreover, observations that cannot be explained by a model with the above four parameters would strongly point to the importance of additional physics, for example, breakdown of well-mixed gaseous opacities (via chemical disequilibrium), non-gray radiative effects, condensation or clouds.

Next, we discuss the observational effort necessary for phase curve observations. While there were initial attempts to monitor phase curves discontinuously with Spitzer (Harrington et al. 2006), more recent Spitzer observations tend to cover the whole course of an orbit to account for long-term instrumental drift (Knutson et al. 2012). Similar continuous observations with JWST would become very time-consuming for planets with period ≳10 days. Although we assumed in Section 4 that the planet is observed for half an orbit, our results indicate that full phase curve coverage is not necessary as long as hot spot offsets are small or negligible (Figure 3). In theory this means that one only needs to observe a planet's thermal emission near primary and secondary eclipse, which would greatly reduce the required observation time. In practice, it will be challenging to relate the observed thermal fluxes from distinct observation periods, but precise characterization of JWST's instrumental drift might still permit discontinuous observation strategies.

Finally, we discuss further physics that might influence our conclusions. Any model necessarily only approximates the dynamics of real atmospheres. The dynamical similarities predicted by our analysis will break down if processes neglected in our model become significant or if the neglected processes are singular perturbations. We do not expect condensation to be a singular perturbation, given that dry models are able to reproduce many aspects of Earth's atmospheric dynamics (Schneider 2006). For a planet with a condensing substance, our method should provide an upper bound on the atmospheric mass. That is because latent heat transport, and for true Earth analogs ocean heat transport, would increase the day–night energy transport. Clouds would similarly reduce the phase curve amplitude, by reducing the dayside brightness temperature while not strongly affecting the nightside brightness temperature (Yang et al. 2013). If this effect is strong enough it can even reverse the expected day–night phase curve pattern. Detection of a inverted phase curve pattern would therefore be a telltale sign of a condensing atmosphere. In all other cases, a planet with a condensing substance should have a reduced phase curve amplitude and thus resemble a dry planet with larger atmospheric mass.

For the radiative transfer, we assume that shortwave and longwave opacities both increase linearly with pressure due to pressure broadening and collision-induced absorption. This assumption breaks down if atmospheric opacities are set by different mechanisms, for example, if shortwave radiation was absorbed by dust (which is insensitive to pressure broadening). We therefore explored simulations in which the shortwave opacity is independent of pressure. We find that reducing the pressure dependency has an effect qualitatively similar to increasing the ratio of shortwave to longwave optical depths, γ, in our standard simulations. Specifically, the resulting increase in shortwave absorption at higher altitudes creates stratospheric inversions, but only has a limited impact on thermal phase curves (∼10% − 15% decrease in phase amplitude compared to our standard simulations with γ = 1; see Figure 4).

We also assume that planets are tidally locked, but the effect of atmospheric thermal tides (Cunha et al. 2014) and trapping in higher-order resonances (Makarov et al. 2012) could result in non-synchronous orbital states. Non-synchronous rotation would introduce an additional nondimensional parameter into our model that compares the length of the day–night cycle with the planet's cooling timescale. We do not expect that a small deviation from synchronous rotation would be a singular perturbation of this parameter. Yang et al. (2014) show that non-synchronously rotating planets with sufficiently slow rotation and/or short cooling timescale smoothly approach the tidally locked regime (their Figure 1(a)). Such planets rotate non-synchronously but in an instantaneous sense still appear tidally locked. Cunha et al. (2014) show that thermal tides are generally not compatible with synchronous rotation. Nevertheless, the deviation from synchronous rotation is small for planets in close orbits with zero eccentricity around small stars, such as M dwarfs (Table 2 in Cunha et al. 2014). Unless synchronous rotation represented a singular limit (see above), the potential effect of thermal tides should then not greatly affect the types of planets we consider here. Many terrestrial planets might never reach synchronous rotation and instead get trapped in a higher-order orbital resonance, like Mercury did in our solar system (Makarov et al. 2012). The probability of being trapped decreases for planets with lower eccentricity (Figure 6 in Makarov et al. 2012), which means that our assumption of synchronous rotation is at least consistent with the fact that also we do not consider non-zero eccentricities. It is beyond the scope of this article to investigate higher-order spin states, but phase curves at visible wavelengths could provide a consistency check for applying our method to planets with optically thin atmospheres (Fujii et al. 2014).

We use dimensional analysis to find a set of six nondimensional parameters that captures the main atmospheric dynamics of dry, tidally locked terrestrial planets in an idealized GCM. We use the GCM to investigate the sensitivity of thermal phase curves to each of the nondimensional parameters. Except for hot and rapidly rotating atmospheres that are optically thick in the longwave, we do not find significant hot spot offsets. On the other hand, the phase curve amplitude remains sensitive to changes in the atmospheric parameters across a large range of atmospheric scenarios. Focusing on the phase amplitude, we find that the phase amplitude of many terrestrial planets is sensitive to only two nondimensional parameters. The main unknowns in the two nondimensional parameters are the surface pressure and the longwave opacity. The longwave opacity can be constrained by transit or emission spectroscopy, in which case the phase amplitude would constrain the surface pressure and atmospheric mass. As an example, we estimate that a broadband phase curve near 18 μm with JWST, taken over a single half-orbit, could be sufficient to constrain the surface pressure of a cool super-Earth to within a factor of two. Constraints like the one we propose will be crucial for understanding atmospheric evolution, in particular atmospheric escape. Moreover, constraining atmospheric mass is important for characterizing the surface conditions of potentially habitable planets.

We are grateful to Nick Cowan, Laura Kreidberg, Feng Ding, Wendy Zhang, and Jacob Bean for their insightful comments on an early draft, and Kevin Heng for thoughtfully reviewing our work. We additionally thank Nick Cowan, Laura Kreidberg, and Sean Mills for discussing JWST precision estimates with us. D.D.B.K. was supported by NSF DMS-0940261, which is part of the Mathematics and Climate Research Network. D.S.A. acknowledges support from an Alfred P. Sloan Research Fellowship. This work was completed with resources provided by the University of Chicago Research Computing Center.

The primitive equations expressed in standard latitude–longitude pressure coordinates (θ, λ, p) are

$$\begin{eqnarray*} \frac{D\mathbf {u}}{Dt} &=& {-}2 \Omega \sin \theta \mathbf {k} \times \mathbf {u} - \nabla \phi -g \frac{\partial \mathbf {F}_{m} }{\partial p}, \\ \frac{\partial \phi }{\partial p}&=& {-}\frac{R T}{p}, \\ \nabla \cdot \mathbf {u} & + & \frac{\partial \omega }{\partial p}= 0, \\ \frac{DT}{Dt}&=& \frac{R T \omega }{c_p p} + \frac{g}{c_p} \frac{\partial F_{{\rm rad}} }{\partial p} + \frac{g}{c_p} \frac{\partial F_{{\rm sens}}}{\partial p} + \frac{g}{c_p} \frac{F_{{\rm conv}}}{\partial p}. \end{eqnarray*}$$

Here u = (u, v) is the horizontal wind velocity, (D/Dt) = (∂/∂t) + u · ∇ + ω(∂/∂p) is the material derivative, k is a unit vector pointing along the axis of planetary rotation, ϕ is the geopotential, T is temperature, ω ≡ (Dp/Dt) is the vertical velocity (expressed as a change in pressure), and dimensional parameters are defined in Section 2. From the top, these equations express conservation of momentum, the hydrostatic approximation, conservation of mass, and conservation of energy. Although the mass conservation equation looks like it assumes incompressibility, it does not. The primitive equations are compressible, but mass conservation can be written in the above simple form using p coordinates and the hydrostatic equation (Vallis 2006, p.79). The forcing terms are

$$\begin{eqnarray*} F_{{\rm rad}}&=&F_{LW}^{\uparrow }-F_{LW}^{\downarrow }-F_{SW}^{\downarrow }, \\ \frac{\partial F_{LW}^{\uparrow }}{\partial p}&=& \frac{2 \kappa _{lw}}{g} \left(F_{LW}^{\uparrow } -\sigma _{SB} T^4\right), \\ \frac{\partial F_{LW}^{\downarrow }}{\partial p}&=& {-}\frac{2 \kappa _{lw}}{g} \left(F_{LW}^{\downarrow }-\sigma _{SB} T^4\right), \\ \frac{\partial F_{SW}^{\downarrow }}{\partial p} & = & {-}\frac{2 \kappa _{sw}}{g} F_{SW}^{\downarrow },\\ F_{SW}^{\downarrow }|_{p=0} & = &\left\lbrace \begin{array}{@{}lc}(1-\alpha) L_* \cos \theta \cos \lambda & \mathrm{if }\ 270^\circ \le \lambda \le 90^\circ \\ 0 & \mathrm{elsewhere}, \end{array}\right.\\ \mathbf {F}_{m} |_{p_s}&=& \rho _s C_D |\mathbf {u}_s| \mathbf {u}_s, \\ F_{{\rm sens}} |_{p_s}&=& \rho _s c_p C_D |\mathbf {u}_s| (T_s-T|_{p_s}). \end{eqnarray*}$$

Here u_s is the surface wind velocity, T_s is the surface temperature, $T|_{p_s}$ is the near-surface air temperature, and ρ_s is the atmospheric density at the surface ($\rho _s = p_s R^{-1} T|_{p_s}$, using the ideal gas law). The convective heat flux F_conv instantaneously adjusts an unstable lapse rate toward the dry adiabat while conserving dry enthalpy

$$\begin{eqnarray*} &&\frac{{d} }{{d} t} \int _0^{p_s} c_p T \frac{dp}{g} = 0. \end{eqnarray*}$$

We neglect scattering in the radiative equations. We assume the hemi-isotropic closure, which is why the radiative equations contain a factor of two (Section 6.4 in Heng et al. 2014). We also assume that opacities increase with p due to pressure broadening and/or collision-induced absorption, so κ_sw = κ_SW(p/p₀) and κ_lw = κ_LW(p/p₀). We define nondimensional shortwave and longwave optical depths as $d\hat{\tau }_{sw}/dp = 2\kappa _{sw}(p)/g$ and $d\hat{\tau }_{lw}/dp = 2\kappa _{lw}(p)/g$. We integrate to find the total optical depths, $\tau _{SW}\equiv \hat{\tau }_{sw}(p_s)=\kappa _{SW} p_0/g \times (p_s/p_0)^2$ and $\tau _{LW}\equiv \hat{\tau }_{lw}(p_s) = \kappa _{LW} p_0/g \times (p_s/p_0)^2$. The surface temperature, T_s, is determined by energy balance

$$\begin{eqnarray*} \mathcal {C} \frac{{d} T_s}{{d} t} & = & F_{SW}^{\downarrow } |_{p_s} - \left(\sigma _{SB} T_s^4 - F_{LW}^{\downarrow }|_{p_s}\right) - F_{{\rm sens}} |_{p_s}. \end{eqnarray*}$$

Here $\mathcal {C}$ is the surface thermal inertia. For a tidally locked planet the stellar forcing does not depend on time. As a first approximation we ignore internal atmospheric variability and assume that T_s is time-independent. This means that $\mathcal {C}$ does not enter the list of dimensional quantities (see below for a discussion of when this assumption is valid).

We form the following nondimensional quantities, marked with the hat symbol: $({D}/{Dt}) = ({1}/{t_{{\rm wave}}}) ({D}/{D\hat{t}})$, $\nabla = \hat{\nabla }/a$, $p=p_s\hat{p}$, $\mathbf {u} = c_{{\rm wave}} \hat{\mathbf {u}}$, $\omega = c_{{\rm wave}} p_s/a \times \hat{\omega }$, $\phi = g H \hat{\phi }$, $T = T_{eq} \hat{T}$, $F_{{\rm rad}} = \sigma _{SB} T_{eq}^4 \hat{F}_{{\rm rad}}$, $F_{{\rm conv}} = \sigma _{SB} T_{eq}^4 \hat{F}_{{\rm conv}}$, $\mathbf {F}_{m} = p_s C_D c_{{\rm wave}}^2 /(R T_{eq}) \times \hat{\mathbf {F}}_{m}$, $F_{{\rm sens}} = p_s c_p C_D c_{{\rm wave}} /R \times \hat{F}_{{\rm sens}}$.

The primitive and radiative equations in nondimensional form are

$$\begin{eqnarray*} \left(\frac{D\hat{\mathbf {u}}}{D\hat{t}} \right) &=& {-}\frac{a^2}{L_{Ro}^2} (\sin \theta \mathbf {k} \times \hat{\mathbf {u}}) - \frac{c_p}{R} (\hat{\nabla }\hat{\phi }) -\frac{C_D a}{H} \left(\frac{\partial \hat{\mathbf {F}}_{m} }{\partial \hat{p}}\right), \\ \left(\frac{\partial \hat{\phi }}{\partial \hat{p}}\right)&=& {-} \left(\frac{\hat{T}}{\hat{p}}\right), \\ (\hat{\nabla } \cdot \hat{\mathbf {u}}) & + & \left(\frac{\partial \hat{\omega }}{\partial \hat{p}}\right)= 0, \\ \left(\frac{D\hat{T}}{D\hat{t}}\right)&=& \frac{R}{c_p} \left(\frac{\hat{T} \hat{\omega }}{\hat{p}}\right) + \frac{t_{{\rm wave}}}{t_{{\rm rad}}} \left(\frac{\partial \hat{F}_{{\rm rad}} }{\partial \hat{p}}\right) + \frac{C_D a}{H} \left(\frac{\partial \hat{F}_{{\rm sens}} }{\partial \hat{p}}\right)\\ &&+ \frac{t_{{\rm wave}}}{t_{{\rm rad}}} \left(\frac{\partial \hat{F}_{{\rm conv}}}{\partial \hat{p}}\right), \\ \left(\frac{\partial \hat{F}_{LW}^{\uparrow }}{\partial \hat{\tau }_{lw}}\right)&=& \hat{F}_{LW}^{\uparrow } - \hat{T}^4 (0\le \hat{\tau }_{lw}\le \tau _{LW}), \\ \left(\frac{\partial \hat{F}_{LW}^{\downarrow }}{\partial \hat{\tau }_{lw}}\right)&=&{-} \left(\hat{F}_{LW}^{\downarrow }- \hat{T}^4\right) (0\le \hat{\tau }_{lw}\le \tau _{LW}), \\ \left(\frac{\partial \hat{F}_{SW}^{\downarrow }}{\partial \tau _{sw}}\right) & = &{-}\hat{F}_{SW}^{\downarrow } (0\le \hat{\tau }_{sw}\le \gamma \tau _{LW}). \end{eqnarray*}$$

The surface energy budget in nondimensional form is

$$\begin{eqnarray*} \frac{t_{{\rm rad},s}}{t_{{\rm rad}}} \left(\frac{{d} \hat{T_s}}{{d} \hat{t}} \right) & = & \hat{F}_{SW}^{\downarrow } |_{p_s} - \left(\hat{T}_s^4 - \hat{F}_{LW}^{\downarrow }|_{p_s} \right) - \frac{t_{{\rm rad}}}{t_{{\rm wave}}} \frac{C_D a}{H} \hat{F}_{{\rm sens}} |_{p_s}, \end{eqnarray*}$$

where $t_{{\rm rad},s} = \mathcal {C}/\sigma _{SB} T_{eq}^3$ is the surface radiative timescale. Our assumption that T_s is time-independent, and t_{rad, s} can be ignored, therefore breaks down when t_{rad, s} ≳ t_rad. In our reference simulations we assume a heat capacity equivalent to that of a well-mixed water layer with depth 3 m. This means the surface thermal inertia, $\mathcal {C}$, is actually large enough that t_{rad, s} ∼ t_rad. To check if this affects our results we recomputed the reference simulations in Figure 4 with $\mathcal {C}$ reduced by a factor of 300, such that t_{rad, s} ≪ t_rad. We note that this heat capacity is far less than realistic values for $\mathcal {C}$. At individual grid points, the time-averaged surface temperature changes up to 2.5% in the cool Earth-size scenario and up to 6% for the hot super-Earth. However, reducing $\mathcal {C}$ affects the phase curve amplitude of both runs by less than 0.5%. This confirms that our choice of six nondimensional parameters captures the dominant dynamics in the GCM simulations (Figure 2).

A standard geographic coordinate system is defined via the radial distance from a planet's center, r, the latitude, θ, which is the angle away from the equator, and the longitude, λ, which is the angle about the planet's north pole. Atmospheres of fast rotating planets are approximately symmetric around the axis of rotation due to conservation of angular momentum, so their time-averaged properties are often displayed as averages over λ. Here we make use of the approximate symmetry of slowly rotating tidally locked planets about the axis connecting the substellar and antistellar points (Figure 1). We define the tidally locked latitude, θ_TL, as the angle away from the terminator, and the tidally locked longitude, λ_TL, as the angle about the substellar point. We choose (θ, λ) = (0, 0) to coincide with the substellar point, and (θ_TL, λ_TL) = (0, 0) to coincide with the north pole (see top row in Figure 1). For example, in tidally locked coordinates λ_TL = 0 and 90° ⩾ θ_TL ⩾ −90° defines the arc that connects substellar and antistellar points via the north pole.

To translate between standard and tidally locked coordinates we first transform both spherical coordinate systems into Cartesian coordinates, so that the north pole lies at (x, y, z) = (0, 0, r) and the substellar point lies at (x, y, z) = (r, 0, 0)

$$\begin{eqnarray} x & = & r \cos \theta \cos \lambda, \nonumber\\ y & = & r \cos \theta \sin \lambda, \nonumber \\ z & = & r \sin \theta, \end{eqnarray} \tag{ B1 }$$

and

$$\begin{eqnarray} x & = & r \sin \theta _{TL}, \nonumber\\ y & = & r \cos \theta _{TL} \sin \lambda _{TL}, \nonumber \\ z & = & r \cos \theta _{TL} \cos \lambda _{TL}. \end{eqnarray} \tag{ B2 }$$

By combining Equations (B1) and (B2) we can express θ_TL and λ_TL in terms of θ and λ:

$$\begin{eqnarray} \theta _{TL} & = & \sin ^{-1}(\cos \theta \cos \lambda), \nonumber\\ \lambda _{TL} & = & \tan ^{-1}\left(\frac{\sin \lambda }{\tan \theta } \right). \end{eqnarray} \tag{ B3 }$$

To plot GCM output in tidally locked coordinates we first express the GCM output in terms of θ_TL and λ_TL, and then linearly interpolate the output onto an evenly spaced (θ_TL, λ_TL) grid.

Transforming GCM wind velocities into tidally locked coordinates is slightly more complicated. Horizontal winds are defined as (u, v) ≡ (rcos θ(Dλ/Dt), r(Dθ/Dt)). We analogously define wind velocities in a tidally locked coordinate system as

$$\begin{eqnarray*} u_{TL} & \equiv & r \cos \theta _{TL} \frac{\mathrm{D} \lambda _{TL}}{\mathrm{D} t} \\ & = & r \cos \theta _{TL} \left(\frac{\partial \lambda _{TL}}{\partial \lambda } \frac{\mathrm{D} \lambda }{\mathrm{D} t} + \frac{\partial \lambda _{TL}}{\partial \theta } \frac{\mathrm{D} \theta }{\mathrm{D} t} \right) \\ & = & \cos \theta _{TL} \left(\frac{\partial \lambda _{TL}}{\partial \lambda } \frac{u}{\cos \theta } + \frac{\partial \lambda _{TL}}{\partial \theta } v \right), \end{eqnarray*}$$

and

$$\begin{eqnarray*} v_{TL} & \equiv & r \frac{\mathrm{D} \theta _{TL}}{\mathrm{D} t} \\ & = & r \left(\frac{\partial \theta _{TL}}{\partial \lambda } \frac{\mathrm{D} \lambda }{\mathrm{D} t} + \frac{\partial \theta _{TL}}{\partial \theta } \frac{\mathrm{D} \theta }{\mathrm{D} t} \right) \\ & = & \frac{\partial \theta _{TL}}{\partial \lambda } \frac{u}{\cos \theta } + \frac{\partial \theta _{TL}}{\partial \theta } v. \end{eqnarray*}$$

We evaluate the partial derivatives (∂λ_TL/∂λ, ∂λ_TL/∂θ, ∂θ_TL/∂λ, ∂θ_TL/∂θ) using Equations (B3). The resulting expressions are long and lead to little insight, so we omit them here.

The area-averaged and observer-projected flux from a planet as seen by a distant observer is

$$\begin{eqnarray*} F(\xi) & = & \frac{\int _{-\pi /2}^{\pi /2} \int _{-\xi -\pi /2}^{-\xi +\pi /2} F_{LW}^{\uparrow }|_{p=0} \cos (\lambda +\xi) \cos ^2(\theta) d\lambda d\theta }{\int _{-\pi /2}^{\pi /2} \int _{-\xi -\pi /2}^{-\xi +\pi /2} \cos (\lambda +\xi) \cos ^2(\theta) d\lambda d\theta }, \end{eqnarray*}$$

where ξ is the phase angle, i.e., the angle between the observer's line of sight and the substellar point (ξ = 0 at secondary eclipse, ξ = π at transit), and $F_{LW}^{\uparrow }|_{p=0}$ is the outgoing thermal flux at the top of atmosphere. This equation is expressed in standard latitude–longitude coordinates, and assumes that the orbit is viewed edge-on (Cowan & Agol 2008). The planet's total flux as seen by a distant observer is πa² × F(ξ).

The maximum thermal flux emitted by a planet corresponds to the flux emitted by the dayside of a bare rock, which we call F_rock. We compute F_rock by setting the outgoing thermal flux equal to the incoming stellar flux at every point, $F_{LW}^{\uparrow }|_{p=0} = L_*(1-\alpha) \cos (\theta) \cos (\lambda)$, so F_rock = 2/3 × L_*(1 − α). We define the dayside-averaged observer-projected temperature of a bare rock as $\sigma _{SB} T_{{\rm rock}}^4 = F_{{\rm rock}}$. This temperature is related to the equilibrium temperature of a planet with effective heat transport via T_rock = (8/3)^1/4T_eq.

We illustrate how to find the set of all transformations on the dimensional parameters that only vary one nondimensional parameter. Our approach is similar to the commonly used technique in dimensional analysis of finding a complete set of nondimensional parameters via matrix methods (Price 2003). The method is general, but we illustrate it assuming (a, Ω, T_eq) are fixed to directly explain our parameter choices in Table 1b. To transform the remaining dimensional parameters, we consider multiplying each of them by a different constant,

$$\begin{eqnarray*} &&(R^{\prime },c_p^{\prime },g^{\prime }, p_s^{\prime },\kappa _{SW}^{\prime },\kappa _{LW}^{\prime },C_D^{\prime }) = \\ &&\quad (C^{v_1} R,C^{v_2} c_p, C^{v_3} g, C^{v_4} p_s, C^{v_5}\kappa _{SW}, C^{v_6}\kappa _{LW}, C^{v_7} C_D), \end{eqnarray*}$$

where a prime denotes a transformed parameter, C is a constant, and (v₁, v₂, ...) are different exponents. At the same time we want to keep all nondimensional parameters except one fixed. As an example we consider all transformations that multiply R/c_p by a factor of C. This means

$$\begin{eqnarray*} \left(\frac{R}{c_{p}} \right)^{\prime } = C \times \left(\frac{R}{c_{p}} \right) & \Rightarrow & C^{v_1} C^{-v_2} \left(\frac{R}{c_{p}} \right) = C \left(\frac{R}{c_{p}} \right) \\ & \Rightarrow & C^{v_1-v_2} = C^1 \\ \left(\frac{a^2}{L_{Ro}^2} \right)^{\prime } = 1 \times \left(\frac{a^2}{L_{Ro}^2} \right) & \Rightarrow & C^{-v_1+v_2/2} \left(\frac{a^2}{L_{Ro}^2} \right) = C^0 \left(\frac{a^2}{L_{Ro}^2} \right) \\ & \Rightarrow & C^{-v_1+v_2/2} = C^0 \\ \left(\frac{t_{{\rm wave}}}{t_{{\rm rad}}} \right)^{\prime } = 1 \times \left(\frac{t_{{\rm wave}}}{t_{{\rm rad}}} \right) & \Rightarrow & C^{-v_1-v_2/2+v_3-v_4} \left(\frac{t_{{\rm wave}}}{t_{{\rm rad}}} \right) \\ && = C^0 \left(\frac{t_{{\rm wave}}}{t_{{\rm rad}}} \right) \\ & \Rightarrow & C^{-v_1-v_2/2+v_3-v_4} = C^0 \\ & \vdots \end{eqnarray*}$$

We consider the exponents of C in the resulting equations. They form a linear system of equations and can be written in matrix form Av = (1, 0, 0, ...) as

$$\begin{eqnarray*} \left(\begin{array}{rcrrcrc}1 & -1 & 0 & 0 & 0 & 0 & 0 \\ -1 & 1/2 & 0 & 0 & 0 & 0 & 0 \\ -1 & -1/2 & 1 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 & 0 \\ 0 & 0 & -1 & 2 & 0 & 1 & 0 \\ -1 & 0 & 1 & 0 & 0 & 0 & 1 \\ \end{array} \right) \left(\begin{array}{c}v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_5 \\ v_6 \\ v_7 \end{array}\right) & = & \left(\begin{array}{c}1\\ 0\\ 0\\ 0\\ 0\\ 0 \end{array}\right) \end{eqnarray*}$$

The matrix columns correspond to the exponents of (R, c_p, g, p_s, κ_SW, κ_LW, C_D), and the matrix rows correspond to each nondimensional parameter. For example, the first row of the matrix corresponds to R/c_p, and has only two non-zero entries (corresponding to the exponents with which R and c_p appear in R/c_p).

This system of equations is underdetermined, that is, it has infinitely many solutions (there are six rows/six equations, but seven columns/seven unknowns). All solutions can be expressed as a particular solution v to the equation Av = (1, 0, 0, ...), plus any vector that lies in the kernel (nullspace) of A: {v + x, where x satisfies Ax = 0}. We express the set of all solutions as

$$\begin{eqnarray*} \lbrace \mathbf {v} + k \cdot \mathbf {x} \rbrace & = & \left(-1, -2, -1-k, 1-k, -3+k, -3+k, k \right), \end{eqnarray*}$$

where k is an arbitrary number, and each entry of this vector corresponds to the power to which C is being raised for each dimensional parameter. For example, to increase R/c_p by C, the first dimensional parameter, R, is multiplied by C⁻¹, the second dimensional parameter, c_p, is multiplied by C⁻², etc.

To vary R/c_p over the largest possible range compatible with the dimensional constraints in Table 2 amounts to finding the largest and smallest values of C that are still consistent with Table 2, while k is allowed to take on any value. We first solve this problem by inspection, and then compare our parameter choices to the values we get from numerical optimization. For our reference case of a cool, Earth-sized planet, we find k = 0 and $1/\sqrt{1.5} \lesssim C \le 1$ (Table 1b).

We note that one can use this method similarly to find the set of all transformations that leave all nondimensional parameters invariant (see Table 1 a, Figure 2). This set of transformations is simply given by the kernel of A.