# Weak seasonality on temperate exoplanets around low-mass stars

Xianyu Tan  Atmospheric, Oceanic and Planetary Physics, Department of Physics,
University of Oxford, OX1 3PU, UK
<EMAIL_ADDRESS>
###### Abstract

Planets with non-zero obliquity and/or orbital eccentricity experience
seasonal variations of stellar irradiation at local latitudes. The extent of
the atmospheric response can be crudely estimated by the ratio between the
orbital timescale and the atmospheric radiative timescale. Given a set of
atmospheric parameters, we show that this ratio depends mostly on the stellar
properties and is independent of orbital distance and planetary equilibrium
temperature. For Jupiter-like atmospheres, this ratio is $\ll 1$ for planets
around very-low-mass M dwarfs and $\gtrsim 1$ when the stellar mass is greater
than about 0.6 solar mass. Complications can arise from various factors,
including varying atmospheric metallicity, clouds and atmospheric dynamics.
Given the eccentricity and obliquity, the seasonal response is expected to be
systematically weaker for gaseous exoplanets around low-mass stars and
stronger for those around more massive stars. The amplitude and phase lag of
atmospheric seasonal variations as a function of host stellar mass are
quantified by idealized analytic models. At the infrared emission level in the
photosphere, the relative amplitudes of thermal flux and temperature
perturbations are negligible, and their phase lags are closed to $-90^{\circ}$
for Jupiter-like planets around very-low-mass stars. The relative amplitudes
and phase lags increase gradually with increasing stellar mass. With a
particular stellar mass, the relative amplitude and phase lag decrease from
low to high infrared optical depth. We also present numerical calculations for
a better illustration of the seasonal behaviors. Lastly, we discuss
implications for the atmospheric circulation and future atmospheric
characterization of exoplanets in systems with different stellar masses.

###### Subject headings:

Exoplanet atmospheres — Planetary atmospheres

##  1\. introduction

Planetary seasonality generally arises from the change of stellar irradiation
at certain latitudes over the orbital timescale in the precence of non-zero
planetary obliquity or orbital eccentricity. Solar-system planets with non-
zero obliquity exhibit significant seasonal variations, which have been
observed on Saturn (see reviews by Fletcher et al. 2018, 2020) and Uranus (see
a review by Fletcher 2021). Even Jupiter, which has a negligible obliquity of
$\sim 3^{\circ}$, shows evidence of an eccentricity season (Nixon et al.
2010). The orbital configurations of exoplanets are more diverse, and we
expect richer seasonal behaviors in the atmospheres of exoplanets.

Previous work has studied the climate and seasonal dynamics of exoplanets and
the implications for observation and habitability. A branch of these studies
has focused on modeling the climate and atmospheric dynamics of terrestrial
planets with Earth-like atmospheres (e.g., Williams & Kasting 1997; Williams &
Pollard 2002; Gaidos & Williams 2004; Ferreira et al. 2014; Shields et al.
2016; Guendelman & Kaspi 2019; Kang et al. 2019; Kang 2019; Palubski et al.
2020). The second category contains investigations of variability, atmospheric
general circulation and potential observable signatures of giant planets with
non-zero obliquity (Langton & Laughlin 2007; Rauscher 2017; Ohno & Zhang
2019a, b) and high eccentricities (Langton & Laughlin 2008; Lewis et al. 2010;
Kataria et al. 2013; Lewis et al. 2014; Mayorga et al. 2021). Some others do
not directly address the climate and seasonal dynamics but try to identify
means by which we could measure or constrain planetary obliquity of exoplanets
(e.g., Kawahara & Fujii 2010; Cowan et al. 2013; Schwartz et al. 2016; Luger
et al. 2021).

Studies of the seasonality of exoplanets have mostly focused on specific
targets or systems around certain types of stars. In this work, we show that,
although the seasonal variation of incoming stellar flux is mainly controlled
by the local latitude, obliquity and eccentricity, given a type of atmosphere,
the characteristic atmospheric response (which may be represented as
nondimensional quantities related to the variation of stellar flux) to such
seasonal variations primarily depends on the stellar properties and is
independent of the orbital distance and planetary equilibrium temperature.
Variations of the mean thermal opacity and specific heat introduce anomalies
of the seasonal response, but the seasonal response is expected to be
systematically weaker on exoplanets around low-mass stars and stronger on
those around higher-mass stars.

Throughout this study, we demonstrate this point for planets with moderate
equilibrium temperatures and thick hydrogen-dominated atmospheres, partly
because the gaseous planets are more likely to be characterized in the near
future via upcoming missions such as the Nancy Grace Roman Space Telescope
(Spergel et al. 2013) and mission concepts such as LUVOIR (Bolcar et al. 2016)
or HabEx (Mennesson et al. 2016) that aim to measure reflected lights of cold
and temperate exoplanets by direct imaging. We expect that a large fraction of
temperate gaseous planets should have sufficient orbital distances to avoid
significant obliquity tidal damping by the host stars if they acquired non-
zero obliquity during their formation and evolution. Our analysis does not
apply to close-in planets that may have been tidally locked111Some close-in
planets may still have chances to maintain non-zero obliquity through
perturbations from additional companions in the system (Su & Lai 2021). .
Atmospheres of terrestrial planets likely have more diversities, including
varying atmospheric mass, composition and surface conditions, which
sufficiently complicate the simple trend that holds for thick hydrogen-
dominated atmospheres.

This paper is organized as follows. In Section 2, we carry out the essential
argument using comparisons between the orbital and atmospheric radiative
timescales. In Section 3, we further present idealized analytic and numerical
radiative transfer calculations to support results in Section 2 and also show
certain quantitative differences. Finally, we discuss limitations and
implications of our results in Section 4.

##  2\. Timescale comparison

The characteristic atmospheric response to the seasonally varying irradiation
at certain latitudes can be intuitively understood by the comparison between
the orbital timescale and the radiative timescale (Ohno & Zhang 2019a), the
latter of which is the characteristic timescale it would take for an
atmosphere to relax its temperature back to the radiative equilibrium via
thermal radiation. The amplitude of the response is expected to be smaller
when the radiative timescale is much longer than the orbital timescale, and
vice versa.

With the cool-to-space approximation (e.g., Andrews et al. 1987; Showman &
Guillot 2002), the radiative timescale is
$\tau_{\rm{rad}}\sim\frac{p}{g}\frac{c_{p}}{4\sigma T_{\mathrm{eq}}^{3}}$,
where $p$ is pressure where the infrared optical depth is around unity, $g$ is
surface gravity, $c_{p}$ is specific heat at constant pressure, and $\sigma$
is the Stefan-Boltzmann constant. This formula ceases to be valid at pressures
that significantly deviate from the thermal photospheric pressure (Showman &
Guillot 2002). Near the mean infrared (IR) photosphere where the mean IR
optical depth is 1, the radiative timescale is then
$\tau_{\rm{rad}}\sim\frac{c_{p}}{4\kappa\sigma T_{\mathrm{eq}}^{3}}$, where
$\kappa$ is the mean IR opacity. Here, $T_{\mathrm{eq}}$ is the planetary
equilibrium temperature defined as a homogeneous blackbody temperature the
planet would have in order to globally radiate away the annual-mean stellar
irradiation:

$\begin{split}4\pi r^{2}_{p}\sigma T_{\mathrm{eq}}^{4}=\langle
F_{\star}\rangle,\\\ \langle
F_{\star}\rangle=\frac{1-A}{\sqrt{1-e^{2}}}\frac{\sigma
T^{4}_{\star}R^{2}_{\star}}{a^{2}}\pi r^{2}_{p},\end{split}$ (1)

where $r_{p}$ is the planetary radius, $a$ is the semimajor axis, $e$ is the
eccentricity, $A$ is the planetary bond albedo, $T_{\star}$ is the stellar
effective temperature, and $R_{\star}$ is the stellar radius. The annual-mean
stellar flux $\langle F_{\star}\rangle$ with non-zero eccentricity has been
well known (e.g., Williams & Pollard 2002). Using Kepler’s third law $P_{\rm
orbit}^{2}\approx\frac{4\pi^{2}}{GM_{\star}}a^{3}$ where $P_{\rm orbit}$ is
the orbital period, $M_{\star}$ is the stellar mass222When the planetary mass
$M_{\rm p}$ is a nontrivial fraction of the total mass of the system, one
should use $P_{\rm orbit}^{2}=\frac{4\pi^{2}}{G(M_{\star}+M_{\rm p})}a^{3}$.
and $G$ is the gravitational constant, the ratio between the orbital period
and the radiative timescale is

$\begin{split}\frac{P_{\rm orbit}}{\tau_{\rm{rad}}}=&\leavevmode\nobreak\
\frac{\sqrt{8}\pi\sigma\kappa
T_{\star}^{3}R_{\star}^{3/2}}{c_{p}\sqrt{GM_{\star}}}\left(\frac{1-A}{\sqrt{1-e^{2}}}\right)^{\frac{3}{4}}\\\
\approx&\leavevmode\nobreak\ 5.94\left(\frac{\kappa}{5\times
10^{-4}}\right)\left(\frac{c_{p}}{13000}\right)^{-1}\left(\frac{T_{\star}}{T_{\odot}}\right)^{3}\\\
&\times\left(\frac{R_{\star}}{R_{\odot}}\right)^{\frac{3}{2}}\left(\frac{M_{\star}}{M_{\odot}}\right)^{-\frac{1}{2}}\left(\frac{1-A}{\sqrt{1-e^{2}}}\right)^{\frac{3}{4}},\end{split}$
(2)

where $\kappa$ and $c_{p}$ are in SI unit, $T_{\odot}$, $R_{\odot}$ and
$M_{\odot}$ are the effective temperature, radius and mass of the sun which
are adopted as 5772 K, $6.957\times 10^{8}$ m and $1.9891\times 10^{30}$ kg,
respectively.

For a particular set of atmospheric parameters of $\kappa$, $c_{p}$ and $A$,
$P_{\rm orbit}/\tau_{\rm{rad}}$ increases rapidly with increasing stellar
luminosity and decreases with stellar mass. This ratio is primarily determined
by properties of the host star (i.e., stellar mass, effective temperature and
radius) and also by the eccentricity to some extent. This suggests that the
stellar host is essential for the planetary seasonality for all planets that
have a certain type of atmosphere in the system regardless their orbital
distance and effective temperature. Below we discuss certain complications.
The mean thermal opacity increases with increasing metallicity in the
atmosphere (e.g., Freedman et al. 2014), which could span more than one order
of magnitude in atmospheres of gaseous planets (e.g., Kreidberg et al. 2014;
Line et al. 2021). Extra complications come from condensation. For instance,
water condenses out in the mean thermal photospheres of all solar-system
gaseous planets. This reduces significant thermal opacity from water vapor,
but at the same time loads thermal and visible opacity from cloud particles
which can be highly uncertain and inhomogeneous. The change of specific heat
is not expected to be significant as long as the atmosphere is dominated by
hydrogen and helium. All major planetary atmospheres in the solar system
exhibit non-zero bond albedo, and their major sources include Rayleigh
scattering and scattering due to clouds and aerosoles (Conrath et al. 1989).
For exoplanets, the bond albedo is highly uncertain and may vary from
negligible to a fraction of one (e.g., Schwartz & Cowan 2015). In contrast to
solar-system planets, the eccentricities of exoplanets can span a wide range
(see http://exoplanets.org/) which could increase the variation of $P_{\rm
orbit}/\tau_{\rm{rad}}$.

Given the complications mentioned above, we compute $P_{\rm
orbit}/\tau_{\rm{rad}}$ as a function of varying $\kappa/c_{p}$ and stellar
mass. Results are shown in Figure 1 as a function of stellar mass and
$\kappa/c_{p}$ relative to a reference value, and each panel contains results
with a set of eccentricity $e$ and bond albedo $A$. We adopt Jupiter-like
reference values of $\kappa_{\rm ref}=5\times 10^{-4}\;\rm m^{2}kg^{-1}$ and
$c_{p,{\rm ref}}=13000\;{\rm Jkg^{-1}K^{-1}}$. This $\kappa_{\rm ref}$ results
in a mean thermal photosphere at about 0.5 bar for Jovian gravity which is
reasonable for Jovian-like conditions (Sromovsky et al. 1998). The relations
between stellar radius, effective temperature and stellar mass are adopted
from the fits in Eker et al. (2018) that are based on a large sample of stars
in the solar neighbourhood.333Since this is an order-of-magnitude estimation,
a slight change of database of stellar properties is not expected to change
our qualitative results and conclusions. At a given $\kappa/c_{p}$, $P_{\rm
orbit}/\tau_{\rm{rad}}$ spans about two orders of magnitude when the stellar
mass varies from 0.18 to 1.7 $M_{\odot}$. For $\kappa/c_{p}$ near the
reference value, $P_{\rm orbit}/\tau_{\rm{rad}}$ is much less than 1 for low-
mass M dwarfs and $>1$ with stellar masses $\gtrapprox 0.6\;M_{\odot}$ with
moderate bond albedo and eccentricity. At sufficiently high $\kappa/c_{p}$
which represents high metallicity in the atmosphere, $P_{\rm
orbit}/\tau_{\rm{rad}}$ is generally close to or larger than 1, suggesting
that metal-rich atmospheres tend to hold strong seasonality regardless the
stellar mass. Nevertheless, seasonality on gaseous planets around low-mass
stars should be statistically weaker than those orbiting more massive stars.

Figure 1.—  The logarithm of the ratio between the orbital period and the
radiative timescale at $\tau=1$ level, $\log(P_{\rm orbit}/\tau_{\rm{rad}})$,
as a function of stellar mass and $\kappa/c_{p}$ relative to the reference
value (the reference infrared opacity is $\kappa_{\rm ref}=5\times
10^{-4}\;{\rm m^{2}kg^{-1}}$ and specific heat is $c_{p,{\rm ref}}=13000\;{\rm
Jkg^{-1}K^{-1}}$). The reference values are appropriate for Jovian atmosphere.
Different bond albedo $A$ and eccentricity $e$ are assumed for results in
different panels as labeled above. Stellar radius and effective temperature as
a function of stellar mass are adopted from Eker et al. (2018). Black lines
represent contours of $P_{\rm orbit}/\tau_{\rm{rad}}=1$. Discontinuity in the
plots arises from the segmental fits in Eker et al. (2018).  Figure 2.—
Relative amplitude (upper panel) and phase lag (bottom panel) of the
temperature anomalies relative to the variation of the radiative equilibrium
temperature as a function of stellar mass. Different bond albedo $A$ and
eccentricity $e$ are assumed for different curves as labelled. The infrared
opacity $\kappa$ is assumed to be $5\times 10^{-4}\;{\rm m^{2}kg^{-1}}$ and
the heat capacity $c_{p}$ is assumed to be $13000\;{\rm Jkg^{-1}K^{-1}}$.

A small $P_{\rm orbit}/\tau_{\rm{rad}}$ leads to a small amplitude and a
significant phase lag of the temperature variations relative to the
irradiation. At a minimum level of complexity, we consider the response of the
local temperature to time varying irradiation crudely as a relaxation to the
radiative equilibrium over the radiative timescale as

$\frac{dT}{dt}=\frac{T_{\rm RE}(t)-T}{\tau_{\rm{rad}}},$ (3)

where $T_{\rm RE}=\overline{T}_{\rm RE}+\Delta T_{\rm RE}\cos{\omega t}$ is
the time-varying radiative equilibrium temperature, $\overline{T}_{\rm RE}$ is
the time-independent component and $\omega=2\pi/P_{\rm orbit}$. The time-
dependent temperature is

$T=\overline{T}_{\rm RE}+\frac{\Delta T_{\rm
RE}}{\sqrt{1+\omega^{2}\tau_{\rm{rad}}^{2}}}\cos{(\omega t+\alpha)},$ (4)

where
$\sin{\alpha}=-\frac{\omega\tau_{\rm{rad}}}{\sqrt{1+\omega^{2}\tau_{\rm{rad}}^{2}}}$
and $\alpha$ is a phase lag of the temperature perturbation relative to the
irradiation. In the presence of a finite $\tau_{\rm{rad}}$, temperature varies
with an amplitude a factor of $1/\sqrt{1+\omega^{2}\tau_{\rm{rad}}^{2}}$ of
that of the radiative equilibrium, and the phase lag $\alpha$ is in between
$-90^{\circ}$ and 0. A larger $P_{\rm orbit}/\tau_{\rm{rad}}$ leads to a
stronger temperature variation and a smaller magnitude of phase lags $\alpha$.
Figure 2 shows the amplitude as a function of stellar mass in the upper panel
and the phase lag in the lower panel assuming $e=0$, $\kappa_{\rm ref}$,
$c_{p,{\rm ref}}$ and a few values for $A$. It appears that the seasonal
temperature variations are negligible and there is a nearly $-90^{\circ}$
phase lag near the photosphere for Jupiter-like planets around stars with
masses $\lessapprox 0.6\;M_{\odot}$.

##  3\. Radiative and temperature responses to time varying irradiation

The timescale comparison has laid out the essence of our argument. In the
following we take a step further and present one-dimensional (1D) analytic and
two-dimensional (2D) numerical calculations of radiative transfer and
atmospheric thermal structure responding to the time-varying irradiation.
These calculations demonstrate some quantitative differences to those
calculated by simple relaxation, but the qualitative conclusions remain the
same. In addition, the radiative transfer models can capture the pressure-
dependent seasonal variations in a more natural way.

In Section 3.1 and 3.2, we perform 1D analytic calculations in which we do not
attempt to realistically model the geometry and magnitude of the seasonal
irradiation but merely capture the characteristic response to a periodically
varying irradiation. Therefore, the models are linearized in which the
response is normalized to the magnitude of the irradiation variation. One may
consider the calculations as for a latitudinal band of a rapidly rotating
planets (such that the diurnal cycle can be negligible) with non-zero
obliquity.

Seasonality at certain latitudes depends critically on the magnitude of the
seasonal variation of irradiation which depends on the planetary obliquity,
eccentricity and the latitude. In many cases, the seasonal variations are not
sinusoidal. To capture the geometry and magnitude of the seasonal variations,
we carry out 2D numerical calculations as a function of latitude and pressure
in Section 3.3.

###  3.1. An analytic semi-grey model

We start with a semi-grey framework in which the stellar irradiation is
represented by a single visible band with a characteristic constant opacity
and the atmospheric thermal IR fluxes are represented with another single band
with a characteristic constant IR opacity. Here we consider a plane-parallel,
two-stream approximation which has been widely used in modeling planetary
atmospheres. With a steady irradiation and boundary conditions, the structure
at radiative equilibrium can be straightforwardly found (e.g., Guillot 2010;
Heng et al. 2012 in the context of hot Jupiters). We do not concern the
equilibrium structure but are interested in the perturbations around the
equilibrium. The variation of irradiation is parameterized as a sinusoidally
varying flux perturbation with an arbitrary amplitude assuming zero
eccentricity. The equation set and analytic solution of perturbations of
temperature and thermal flux are given in the Appendix A.

The seasonal variations of irradiation are not sinusoidal in many cases. In
the presence of a non-zero eccentricity and obliquity, the seasonal variation
at any latitude is not sinusoidal; even with a circular orbit, the variation
is not sinusoidal at some latitudes. However, the seasonal variation of
irradiation may be approximated as a linear combination of sinusoids with
different frequencies and amplitudes. The thermal response can then be the
linear combination of solutions to these sinusoids because of the linearity of
our problem.

One may express perturbations of the top-of-atmosphere (TOA) flux
$F^{\prime}(\tau=0,t)$ and photospheric temperature $T^{\prime}(\tau=1,t)$ in
the form of

$\begin{split}\frac{F^{\prime}(\tau=0,t)}{\Delta
F_{v}}=\frac{|F^{\prime}|}{\Delta F_{v}}\exp(i\alpha_{f})\exp(i\omega t),\\\
T^{\prime}(\tau=1,t)=|T^{\prime}(\tau=1)|\exp(i\alpha_{T})\exp(i\omega
t),\end{split}$ (5)

where $\Delta F_{v}$ is the magnitude of the irradiation variation. Physical
quantities are the real parts of Eq. (5). The relative amplitude
$|F^{\prime}|/\Delta F_{v}$ and phase lag $\alpha_{f}$ of the TOA thermal flux
perturbation, as well as the phase lag of the temperature perturbation
$\alpha_{T}$ are determined only by $P_{\rm orbit}/\tau_{\rm{rad}}$. The
amplitude of temperature perturbations $|T^{\prime}|$ depends on $T_{\rm eq}$
and $\Delta F_{v}$ too, but the ratio to a reference amplitude
$|T^{\prime}|/|T^{\prime}_{\rm ref}|$ depends on $P_{\rm
orbit}/\tau_{\rm{rad}}$ only, where $T^{\prime}_{\rm ref}$ is a perturbation
with an extremely large $P_{\rm orbit}/\tau_{\rm{rad}}$ which is arbitrarily
set to $10^{4}$.

As an example, we assume $\kappa=\kappa_{\rm ref}$, $c_{p}=c_{p,{\rm ref}}$,
visible opacity $\kappa_{v}=2\times 10^{-4}\;{\rm m^{2}kg^{-1}}$, zenith angle
of the irradiation $\mu_{v}=1$, $A=0.5$ and $e=0$. Time evolution of the
seasonal irradiation (with a reversed sign) and the TOA thermal flux
variations for cases with different stellar mass are shown in Figure 3. Both
the flux magnitude and the phase lag decrease with decreasing stellar mass.
More results of relative amplitudes and phase lags for thermal flux and
temperature perturbations are shown in Figure 4 as a function of stellar mass
for cases with a few different bond albedos. The relative amplitudes and phase
lags of both $F^{\prime}(\tau=0)$ and $T^{\prime}(\tau=1)$ as a function of
stellar mass exhibit the same general trend as that in Section 2. These
relative amplitudes are $\lessapprox 0.1$ at low stellar mass ($\lessapprox
0.6\;M_{\odot}$), quantitatively similar to those using relaxation method
(Figure 2). However, they increase only up to $\lesssim 0.6$ at
$M_{\star}=1.7\;M_{\odot}$, and this upper bound is obviously smaller than
that (above 0.95) in the relaxation calculation. All phase lags from two types
of models are nearly $-90^{\circ}$ at $M_{\star}<0.2\;M_{\odot}$. The phase
lags in the radiative transfer calculation increase (i.e., the amplitudes
decrease) more rapidly with increasing $M_{\star}$ when $M_{\star}\lessapprox
1\;M_{\odot}$ than that in the relaxation calculation, but their upper bound
is lower at the highest $M_{\star}$. The difference between results from the
radiative transfer and relaxation calculations at high stellar mass presumably
arises from the fact that communications between the deeper layers where there
is a growing thermal inertia and the photosphere is naturally taken into
account in the radiative transfer calculation while not in the relaxation
calculation.

The mean visible opacity is usually smaller than the mean IR opacity in
temperate atmospheres (Freedman et al. 2014). However, enhanced absorption of
irradiation may be possible. We have calculated cases with
$\kappa_{v}>\kappa$. The relative amplitudes and phase lags of the TOA thermal
flux variations are generally larger than those shown in Figure 4. This is
reasonable because in this case energy is absorbed and re-emitted at a lower
thermal optical depth where the thermal inertia is smaller.

We have focused on the TOA flux and temperature variations near the thermal
photosphere. The amplitude of the perturbations decreases exponentially with
increasing optical depth. They are difficult to observe and may be easily
affected by other sources such as atmospheric dynamics or interior convection.
Thus, the deep perturbations in our calculations are not discussed in the main
text but are referred to Appendix A for details.

Figure 3.— Time evolution of seasonal irradiation (with a reversed sign) and
TOA thermal flux variation for cases with different stellar mass. The dots
mark the peak fluxes at their first period and illustrate phase lags between
different cases. Parameters used in these cases are $\kappa=5\times
10^{-4}\;{\rm m^{2}kg^{-1}}$, $c_{p}=13000\;{\rm Jkg^{-1}K^{-1}}$,
$\kappa_{v}=2\times 10^{-4}\;{\rm m^{2}kg^{-1}}$, $\mu_{v}=1$, $A=0.5$ and
$e=0$.

Figure 4.— _Left:_ Relative amplitude (upper panel) and phase lag (bottom
panel) of the top-of-atmosphere (TOA) thermal flux anomalies relative to the
variation of the irradiation as a function of stellar mass. Different bond
albedo $A$ are assumed for different curves, and eccentricity is assumed to be
zero for all models. The infrared opacity $\kappa$ is assumed to be $5\times
10^{-4}\;{\rm m^{2}kg^{-1}}$ and the heat capacity $c_{p}$ is assumed to be
$13000\;{\rm Jkg^{-1}K^{-1}}$. _Right_ : Relative amplitude of the temperature
anomalies at optical depth of one relative to that calculated using a large
ratio of $P_{\rm orbit}/\tau_{\rm{rad}}=10^{4}$ in the upper panel. The bottom
panel shows the phase lag of temperature anomalies at $\tau=1$ relative to the
irradiation anomaly.

###  3.2. An analytic non-grey model: picket fence in the thermal emission and
arbitrary bands in the visible absorption

We are also interested in the seasonal responses at different layers because
these might provide observational signatures to better constrain planetary
seasonality. Non-grey (wavelength dependent) effects in the absorption and
emission affect the radiative heating and cooling rates and therefore the
amplitudes and phase lags at different layers. To capture the non-grey effects
in the upper atmospheric layers, we consider a model with the picket fence
opacity structure in the thermal band and arbitrary bands in the visible. This
setup enables improvements in the radiative transfer calculation while
remaining analytically trackable. We still consider absorption only.

The picket fence model considers two different values for the thermal opacity
$\kappa_{1}$ and $\kappa_{2}$, and the characteristic band width associated
with them are assumed to be much smaller than the width over which the Planck
function varies. The Planck function is considered constant over a
characteristic band width that contains both $\kappa_{1}$ and $\kappa_{2}$.
Therefore, the total (wavelength-integrated) thermal emission through bands
with $\kappa_{1}$ and $\kappa_{2}$ are approximated as $(1-\chi)B$ and $\chi
B$, respectively, where $\chi$ is the fractional band width associated with
$\kappa_{2}$. For more detailed descriptions of the picket fence models, see
Chandrasekhar (1935); King (1955) and Parmentier & Guillot (2014). In the
visible absorption, we consider partitioning of irradiative energy into
multiple channels with different opacities. The radiative transfer equations
and analytic solution are given in the Appendix B.

We perform an example to show the optical-depth dependent response. In
principle, the values of $\kappa_{1}$, $\kappa_{2}$ and visible opacities may
be chosen according to the full opacities for specific cases with certain
atmospheric temperature, composition, and stellar spectrum (Parmentier et al.
2015; Lee et al. 2021). Here we are interested in a general conceptual
illustration, and so in the following results we do not tune the opacities to
match specific cases. Figure 5 shows relative amplitudes on the left and phase
lags on the right for temperature variations as a function of optical depth
for a set of cases with different stellar mass and with $\kappa_{1}=5\times
10^{-4}\;\rm m^{2}kg^{-1}$, $\kappa_{2}=5\times 10^{-3}\;\rm m^{2}kg^{-1}$,
$\chi=0.3$ and three visible bands with $\kappa_{v}=[2\times 10^{-4},2\times
10^{-3},2\times 10^{-2}]\;\rm m^{2}kg^{-1}$ and partition coefficients of
visible bands $\epsilon=[0.7,0.2,0.1]$. As expected, the temperature response
is weak at high optical depth and smoothly increases with decreasing
$\tau_{1}$. Planets around higher mass stars have overall larger temperature
perturbations due to the longer orbital period compared to the radiative
timescale. The phase lags of all cases are relatively large at low optical
depth and decrease with increasing optical depth. These phase lags all go
across $-90^{\circ}$ from low to high $\tau_{1}$ and this crossing occurs at
higher $\tau_{1}$ for cases with higher stellar mass. After the $-90^{\circ}$
crossing, the phase lag in the case with 0.3 $M_{\odot}$ quickly converges to
nearly $-90^{\circ}$ at high optical depth; the one with 1 $M_{\odot}$ first
shows oscillations then converges to a value $<-100^{\circ}$ at depth; the one
with 1.7 $M_{\odot}$ continuously decreases with increasing $\tau_{1}$. More
discussion of the deep perturbations is referred to Appendix A.

The left panels in Figure 6 show amplitudes and phase lags of TOA thermal
fluxes as a function of stellar mass using atmospheric parameters as those
used in Figure 5. The TOA fluxes show the same trend as that in the semi-grey
case in Figure 4, but overall the $|F^{\prime}|/\Delta F_{v}$ is slightly
higher and the phase lags magnitude is smaller. This is because some portion
of the irradiation is absorbed and re-emitted from lower pressures where the
responsive timescales are shorter. On the right column, we display temperature
perturbation amplitudes and phase lags at different optical depth. At higher
stellar masses, the phase lag differences between different optical depth are
relatively small ($\lessapprox 30^{\circ}$), whereas they show a maximum of
about $60^{\circ}$ in cases with the stellar mass of about 0.6 $M_{\odot}$.
Then the phase lag differences decline again at lower stellar mass.

Saturn has the most detailed measurements of seasonal variability among
gaseous planets thanks to the longevity of the Cassini mission. Assuming that
the atmospheric parameters are similar to those of Jupiter with $A\sim 0.3$,
$P_{\rm orbit}/\tau_{\rm{rad}}$ is about 4.6 for Saturn according to Eq. (2),
suggesting a somewhat strong seasonality. Indeed, observations and detailed
radiative transfer calculations show obvious seasonal hemispheric asymmetry in
Saturn’s stratosphere. The amplitudes of temperature variations increase with
decreasing pressure. For instance, the amplitude of the south pole is about 40
K, 10 K and 4 K at pressure between 1 to 5 mbar, 110 mbar and 440 mbar,
respectively. The phase lag is about $-30^{\circ}$ and $-15^{\circ}$ at 1 mbar
and 0.1 mbar, respectively. See Fletcher et al. (2018, 2020) for comprehensive
reviews. Here we do not attempt to model the case for Saturn, but these
observations are in agreement to our qualitative argument.

Figure 5.— _Left:_ Relative amplitudes of the temperature variations relative
to that calculated using a large ratio of $P_{\rm
orbit}/\tau_{\rm{rad}}=10^{4}$ as a function of optical depth from a non-grey
model with three different stellar mass. These calculations assume infrared
opacities of $\kappa_{1}=5\times 10^{-4}\;\rm m^{2}kg^{-1}$,
$\kappa_{2}=5\times 10^{-3}\;\rm m^{2}kg^{-1}$, $\chi=0.3$ and three visible
bands with $\kappa_{v}=[2\times 10^{-4},2\times 10^{-3},2\times 10^{-2}]\;\rm
m^{2}kg^{-1}$ and $\epsilon=[0.7,0.2,0.1]$. The bond albedo $A$ and the
eccentricity $e$ are assumed zero. _Right:_ The phase lags of the temperature
variations relative to the irradiation variation as a function of optical
depth for cases shown on the left panel. The grey dashed line represents a
$-90^{\circ}$ phase lag.

Figure 6.— _Left:_ Relative amplitudes (upper panel) and phase lags (bottom
panel) of the TOA thermal flux variations relative to the variation of the
irradiation as a function of stellar mass from a set of non-grey cases.
Different bond albedo $A$ are assumed for different curves, and the
eccentricity is assumed to be zero for all models. Parameters are the same as
those used in Figure 5. _Right:_ Relative amplitudes (upper panel) and phase
lags (bottom panel) of temperature variations at three optical depth
$\tau_{1}=0.01,0.1$ and 1 as a function of stellar mass for the a subset of
non-grey cases.

###  3.3. Numerical two-dimensional semi-grey calculations

To capture the geometry and magnitude of the seasonal irradiation at different
latitudes, we turn to numerical integration of radiative transfer equations
performed over a latitudinal grid. At each latitude, the irradiation is
treated as the longitudinal mean which is reasonable for temperate, rapidaly
rotating exoplanets (Showman et al. 2015; Rauscher 2017).

The stellar irradiation pattern in the presence of non-zero obliquity and
eccentricity is described as follow. The latitude of the substellar point
$\phi_{\rm ss}$ in our model is given by

$\sin\phi_{\rm ss}=\sin\psi\cos{f},$ (6)

where $\psi$ is the planet’s obliquity and $f$ is the true anomaly. The
relation between $f$ and time $t$ is linked by the eccentric anomaly $E$
(Murray & Dermott 1999):

$\tan{\frac{f}{2}}=\sqrt{\frac{1+e}{1-e}}\tan{\frac{E}{2}}.$ (7)

The time evolution of $E$ is given by the Kepler’s equation

$\frac{dE}{dt}=\frac{2\pi}{P_{\rm orbit}}\frac{1}{1-e\cos{E}}.$ (8)

The star-planet distance is $d=a(1-e\cos{f})$. Equation (6) states that the
maximum Northern excursion of the substellar latitude ($\phi_{\rm ss}=\psi$)
occurs when the planet is at the periapsis ($f=0$). Of course, Equation (6) is
only one of the possible configurations, but here we show this particular case
for illustrations of seasonal behaviors in the presence of both non-zero
obliquity and eccentricity. The longitudinal averaged stellar flux at the top
of the atmosphere at certain latitude $\phi$ is (e.g., Liou 2002)

$\begin{split}F_{v}(\phi,t)=&\frac{(1-A)\sigma
T^{4}_{\star}(R_{\star}/d)^{2}}{\pi}\times\\\ &(\sin{\phi}\sin{\phi_{\rm
ss}}H+\cos{\phi}\cos{\phi_{\rm ss}}\sin{H}),\end{split}$ (9)

where $H=\cos^{-1}(-\tan{\phi}\tan{\phi_{\rm ss}})$ is the length from star-
rise to noon at latitude $\phi$ measured in radiance. When the argument in the
inverse cosine function is larger than 1 or smaller than -1, $H$ is set to 0
or $\pi$, respectively. Note that $d$, $H$ and $\phi_{\rm ss}$ are functions
of time $t$.

The 2D numerical model integrates a set of radiative transfer equations
(Equations 9, A1 and A3) as a function of latitude and time. We adopt a
thermal opacity of $\kappa=5\times 10^{-4}\;{\rm m^{2}kg^{-1}}$, a visible
opacity of $\kappa_{v}=2\times 10^{-4}\;{\rm m^{2}kg^{-1}}$, and specific heat
$c_{p}=13000\;{\rm Jkg^{-1}K^{-1}}$, the same as those in Section 3.1. We
utilize the numerical package ${\it twostr}$ (Kylling et al. 1995) to
integrate the radiative transfer equations, and the mean zenith angle of
irradiation is assumed to be $1/\sqrt{3}$ following the default value in ${\it
twostr}$. The model is initialized with a uniform temperature-pressure profile
that is in a radiative equilibrium with respect to the global-annual-mean
irradiation. At the lower boundary layer, we apply a net upward heat flux with
an internal temperature of 100 K, appropriate for Jupiter-like planets. The
simulations are usually integrated for over $10^{5}$ Earth days after which
they are in statistical equilibrium.

We first show results for planets with a circular orbit ($e=0)$, a prescribed
bond albedo $A=0.3$ and a fixed equilibrium temperature $T_{\rm eq}=250$ K,
but with different obliquity and around stars with different masses. The top
and middle rows in Figure 7 show the top-of-atmosphere (TOA) thermal flux as a
function of time and latitude for cases with stellar masses of
$M_{\star}=0.4$, 1 and 1.6 $M_{\odot}$ and two obliquities of $30^{\circ}$ and
$60^{\circ}$. The stellar effective temperature $T_{\star}$ and radius
$R_{\star}$ are found using relations in Eker et al. (2018). With the given
$T_{\rm eq}$, the orbital periods for cases with 0.4, 1 and 1.6 $M_{\odot}$
and a zero eccentricity are about 28, 381, 1427 Earth days, respectively,
based on Equation (1). In each panel, the dashed line is the substellar
latitude as a function of time. As is easily visualized, for cases with both
obliquities, the TOA thermal flux shows stronger hemispheric asymmetry and
higher flux variation over the seasonal cycles around stars with higher
masses. Magnitude of the phase lag between the thermal flux and irradiation
(represented by the dash lines) increases with decrease stellar mass. The case
with 0.4 $M_{\odot}$ exhibits negligible seasonal variation and hemispheric
asymmetry, maintaining a nearly constant hot-equator and cold-poles
configuration over the seasons for the case with an obliquity of $30^{\circ}$.
Cases with an obliquity of $60^{\circ}$ are at the regime in which the annual-
mean irradiation at the poles are greater than that at the equator (Ward
1974), and therefore the one with 0.4 $M_{\odot}$ and $60^{\circ}$ obliquity
shows a nearly constant cold-equator and hot-poles configuration.

We also present similar results with a moderate eccentricity of $e=0.3$ and an
obliquity of $30^{\circ}$ in the bottom row of Figure 7. The time variation of
the substellar latitude is not sinusoidal as a result of the non-zero
eccentricity. The non-zero eccentricity also amplifies the seasonal variation
in the specific configuration given by Equation (6). In cases with higher
stellar masses, the seasonal flux variations are more pulsive than those shown
in cases with $e=0$ (first row). Despite the larger amplitude of the seasonal
variation of irradiation, the thermal flux in the case with 0.4 $M_{\odot}$
still exhibits a rather weak seasonal variability. Overall, the picture that
planets around lower-mass stars tend to exhibit weaker seasonal response than
those around higher-mass stars holds well in the presence of non-zero
eccentricity.

Figure 7.— Top-of-atmosphere thermal flux as a function of latitude and time.
The dash line in each panel is the substellar latitude as a function of time.
All models assume an equilibrium temperature of 250 K, a bond albedo $A=0.3$,
a thermal opacity of $\kappa=5\times 10^{-4}\;{\rm m^{2}kg^{-1}}$, a visible
opacity of $\kappa_{v}=2\times 10^{-4}\;{\rm m^{2}kg^{-1}}$ and the specific
heat of $c_{p}=13000\;{\rm Jkg^{-1}K^{-1}}$. These models have been integrated
over $10^{5}$ Earth days and have been in statistical equilibrium. _Top row:_
results of cases with an obliquity of $30^{\circ}$, a zero eccentricity but
with three stellar masses of 0.4 (left), 1 (middle) and 1.6 $M_{\odot}$
(right). The orbital periods are about 28, 381 and 1427 Earth days,
respectively. The stellar effective temperatures are 3564 K and 7001 K, and
the stellar radius are 0.34 and 1.9 $R_{\odot}$ for cases with $M_{\star}=0.4$
and 1.6 $M_{\odot}$, respectively. _Middle row_ : same as those at the top row
except that the obliquity is $60^{\circ}$. _Bottom row:_ same as thoes at the
top row except that the eccentricity is $e=0.3$. Accordingly, the orbital
periods are about 33, 450 and 1687 Earth days for stellar mass of 0.4, 1 and
1.6 $M_{\odot}$, respectively.

##  4\. Discussion

In this work, we have shown that given a set of atmospheric parameters of
gaseous planets, the atmospheric response to the seasonal variation of stellar
irradiation due to non-zero planetary obliquity and/or orbital eccentricity
primarily depends on the stellar properties. Despite complications could arise
from diverse atmospheric conditions, gaseous planets around low-mass stars
should exhibit statistically weaker seasonality than those around more massive
stars. We carried out timescale comparisons in Section 2 for the essential
argument. In Section 3, we showed further demonstrations using both analytic
and numerical calculations of radiative transfer. This finding may help
interpret future atmospheric characterization of exoplanets around different
types of stars via directly imaging missions.

There are interesting implications for the atmospheric general circulation of
planets in systems of different stellar mass. Using 2D shallow-water dynamical
models, Ohno & Zhang (2019a) summarized five circulation regimes depending on
the planetary obliquity and the atmospheric radiative timescale
$\tau_{\rm{rad}}$ (see their Figure 1). In this study, we show that their
regime IV and V with $\tau_{\rm{rad}}\gg P_{\rm orbit}$ can mostly apply to
gaseous planets around low-mass stars, while their regime II and III with
$\tau_{\rm{rad}}\lesssim P_{\rm orbit}$ likely apply to those around more
massive stars. As such, our work is complementary to that in Ohno & Zhang
(2019a).

A critical obliquity of $\psi_{c}\approxeq 54^{\circ}$ is well-known, below
which the annual-mean irradiation is maximized at the equator and above which
the annual-mean irradiation is maximized at the poles (Ward 1974). It also
suggests that, the closer the planetary obliquity to $\psi_{c}$, the smaller
the meridional gradient of the longitudinal-annual-mean irradiation. Planets
with an obliquity sufficiently close to the critical value $\psi_{c}$ around
low-mass stars may result in not only weak seasonal variation but also small
horizontal temperature variation. In this regime, the driving force of the
global circulation that is solely raised by the irradiation difference can be
weak. Other forms of the driving force, such as that powered by convective
perturbations (Williams 1978; Showman et al. 2019) could potentially shape the
circulation and drive multiple strong zonal jets. Planets with an obliquity
sufficiently close to $\psi_{c}$ but around more massive stars could host
significant seasonal hemispheric asymmetry of the temperature structure
(regimes II and III in Ohno & Zhang 2019a). The circulation patterns should be
constrained by the temperature difference and exhibit seasonal hemispheric
asymmetry. This circulation transition of planets with obliquities close to
$\psi_{c}$ around moderate-mass stars to low-mass stars would be interesting
to investigate. On the other hand, planets with obliquities much different
from $\psi_{c}$ around stars of any mass could sustain large meridional
temperature gradients, regardless whether the temperature patterns vary
strongly over the season. Their atmospheric circulation should be largely
shaped by the differential irradiation and their seasonal variations as shown
in Ohno & Zhang (2019a).

Planetary rotation plays a key role in the circulation as well. Near low
latitudes, horizontally propagating gravity waves and thermally driven
overturning circulations (analog to the Hadley circulation on Earth) could be
efficient to remove temperature anomalies. The latitudinal width of this
tropical zone is sensitive to rotation, and is expected to the wider for
slower rotators (Held & Hou 1980). At mid-to-high latitudes where rotation
dominates the large-scale dynamics, meridional heat transport efficiency is
expected to be weaker than that at low latitudes and decreases with increasing
rotation rate (Kaspi & Showman 2015). We expect that the meridional
temperature structure of rapid rotators may be able to hold closer to that
determined purely by seasonal irradiation, while atmospheric dynamics of slow
rotators might shape the seasonality further in addition to that determined
just by the radiative response. These may have consequence on cloud formation,
haze and chemistry transport, which could have effects on reflective lights
probed by direct imaging of exoplanets.

Our analytic models are vastly idealized, and multiple factors could introduce
complications. First, as mentioned above, atmospheric dynamics can transport
heat horizontally and vertically. For rapidly rotating planets such as Saturn,
dynamical effects on the seasonal pattern is not dominant (Fletcher et al.
2020). But this can be important for slow rotators. Second, as alluded in
Section 2, variations of opacity due to varying atmospheric metallicity and
temperature contribute to the scattering on the simple trend predicted for
hydrogen-dominated atmospheres. Third, the presence and variations of clouds
and hazes represent another major uncertainty in the atmospheres. They
contribute to the energy budget of both visible and thermal bands, and
sometimes can be dominant. Their coupling to the atmospheric dynamics makes
quantitative predictions more challenging. Finally, the presence of the
internal heat in the gaseous planets when the internal heat is comparable or
greater than the irradiation could disrupt the trend predicted in this study.
The onset of convection could eliminate the seasonal patterns, which is the
case for Saturn (Fletcher et al. 2018, 2020). Another example is Jupiter.
Despite that there is an equator-to-pole insolation difference, its thermal
emission does not show an equator-to-pole pattern which may be attributed to
the interior convection of Jupiter (Ingersoll & Porco 1978).

Future work should include quantitative predictions of seasonal dynamics for
planets around different stellar types using more realistic models. Such
models should include global-scale three-dimensional dynamics which can
capture vertically propagating waves, baroclinic instability and their
interactions with the jet streams, as well as proper representation of
radiative transfer that represents pressure dependent thermal damping.
Parameterization of chemical reactions and haze formation, their dynamical
transport and effects on the energy budget are interesting to include. Post-
processing using realistic radiative transfer would be necessary to identify
observational signatures of the seasonal dynamics and would be useful for
future observational interpretations and strategies. One important aspect is
that the observational signature of seasonal variation should be sensitive to
pressure because the radiative timescale varies with pressure. A combination
of low- and high-resolution spectroscopy will be a valuable tool for
diagnosing seasonal variability as they can measure layers at both relatively
high and low pressures. The study of seasonal responses over the planet
population around different types of stars will be somewhat analog to the day-
to-night temperature variations of the hot Jupiter population, which similarly
invoke comparison of the radiative timescale to other key timescales in the
system (Perez-Becker & Showman 2013; Komacek & Showman 2016; Komacek et al.
2017, see also the latest data compilation in Wong et al. 2021). We expect
trends in the planet population but also significant deviations due to
intrinsic complications in their atmospheres.

We thank discussion with A.P. Showman, R.T. Pierrehumbert, V. Parmentier and
X. Zhang. We thank the reviewer for constructive and encouraging comments.
This study is supported by European Research Council Advanced Grant (grant
agreement No. 740963/EXOCONDENSE, PI: R.T. Pierrehumbert).

## References

  * Andrews et al. (1987) Andrews, D. G., Holton, J. R., & Leovy, C. B. 1987, Middle atmosphere dynamics, Vol. 40 (Academic press)
  * Bolcar et al. (2016) Bolcar, M. R., Feinberg, L., France, K., Rauscher, B. J., Redding, D., & Schiminovich, D. 2016, in Space telescopes and instrumentation 2016: Optical, infrared, and millimeter wave, Vol. 9904, International Society for Optics and Photonics, 99040J
  * Chandrasekhar (1935) Chandrasekhar, S. 1935, MNRAS, 96, 21
  * Conrath et al. (1989) Conrath, B. J., Hanel, R. A., & Samuelson, R. E. Thermal structure and heat balance of the outer planets., ed. S. K. Atreya, J. B. Pollack, & M. S. Matthews, 513–538
  * Cowan et al. (2013) Cowan, N. B., Fuentes, P. A., & Haggard, H. M. 2013, Monthly Notices of the Royal Astronomical Society, 434, 2465
  * Eker et al. (2018) Eker, Z., Bakış, V., Bilir, S., Soydugan, F., Steer, I., Soydugan, E., Bakış, H., Aliçavuş, F., Aslan, G., & Alpsoy, M. 2018, Monthly Notices of the Royal Astronomical Society, 479, 5491
  * Ferreira et al. (2014) Ferreira, D., Marshall, J., O’Gorman, P. A., & Seager, S. 2014, Icarus, 243, 236
  * Fletcher (2021) Fletcher, L. N. 2021, arXiv e-prints, arXiv:2105.06377
  * Fletcher et al. (2018) Fletcher, L. N., Greathouse, T. K., Guerlet, S., Moses, J. I., & West, R. A. Saturn’s Seasonally Changing Atmosphere. Thermal Structure Composition And Aerosols, ed. K. H. Baines, F. M. Flasar, N. Krupp, & T. Stallard, 251–294
  * Fletcher et al. (2020) Fletcher, L. N., Sromovsky, L., Hue, V., Moses, J. I., Guerlet, S., West, R. A., & Koskinen, T. 2020, arXiv e-prints, arXiv:2012.09288
  * Freedman et al. (2014) Freedman, R. S., Lustig-Yaeger, J., Fortney, J. J., Lupu, R. E., Marley, M. S., & Lodders, K. 2014, The Astrophysical Journal Supplement Series, 214, 25
  * Gaidos & Williams (2004) Gaidos, E. & Williams, D. 2004, New Astronomy, 10, 67
  * Guendelman & Kaspi (2019) Guendelman, I. & Kaspi, Y. 2019, The Astrophysical Journal, 881, 67
  * Guillot (2010) Guillot, T. 2010, A&A, 520, A27
  * Held & Hou (1980) Held, I. M. & Hou, A. Y. 1980, Journal of the Atmospheric Sciences, 37, 515
  * Heng et al. (2012) Heng, K., Hayek, W., Pont, F., & Sing, D. K. 2012, Monthly Notices of the Royal Astronomical Society, 420, 20
  * Ingersoll & Porco (1978) Ingersoll, A. P. & Porco, C. C. 1978, Icarus, 35, 27
  * Kang (2019) Kang, W. 2019, ApJ, 877, L6
  * Kang et al. (2019) Kang, W., Cai, M., & Tziperman, E. 2019, Icarus, 330, 142
  * Kaspi & Showman (2015) Kaspi, Y. & Showman, A. P. 2015, The Astrophysical Journal, 804, 60
  * Kataria et al. (2013) Kataria, T., Showman, A. P., Lewis, N. K., Fortney, J. J., Marley, M. S., & Freedman, R. S. 2013, The Astrophysical Journal, 767, 76
  * Kawahara & Fujii (2010) Kawahara, H. & Fujii, Y. 2010, The Astrophysical Journal, 720, 1333
  * King (1955) King, I. J. I. F. 1955, ApJ, 121, 711
  * Komacek & Showman (2016) Komacek, T. D. & Showman, A. P. 2016, The Astrophysical Journal, 821, 16
  * Komacek et al. (2017) Komacek, T. D., Showman, A. P., & Tan, X. 2017, ApJ, 835, 198
  * Kreidberg et al. (2014) Kreidberg, L., Bean, J. L., Désert, J.-M., Line, M. R., Fortney, J. J., Madhusudhan, N., Stevenson, K. B., Showman, A. P., Charbonneau, D., McCullough, P. R., Seager, S., Burrows, A., Henry, G. W., Williamson, M., Kataria, T., & Homeier, D. 2014, ApJ, 793, L27
  * Kylling et al. (1995) Kylling, A., Stamnes, K., & Tsay, S.-C. 1995, Journal of Atmospheric Chemistry, 21, 115
  * Langton & Laughlin (2007) Langton, J. & Laughlin, G. 2007, The Astrophysical Journal Letters, 657, L113
  * Langton & Laughlin (2008) —. 2008, The Astrophysical Journal, 674, 1106
  * Lee et al. (2021) Lee, E. K. H., Parmentier, V., Hammond, M., Grimm, S. L., Kitzmann, D., Tan, X., Tsai, S.-M., & Pierrehumbert, R. T. 2021, MNRAS, 506, 2695
  * Lewis et al. (2014) Lewis, N. K., Showman, A. P., Fortney, J. J., Knutson, H. A., & Marley, M. S. 2014, The Astrophysical Journal, 795, 150
  * Lewis et al. (2010) Lewis, N. K., Showman, A. P., Fortney, J. J., Marley, M. S., Freedman, R. S., & Lodders, K. 2010, The Astrophysical Journal, 720, 344
  * Line et al. (2021) Line, M. R., Brogi, M., Bean, J. L., Gandhi, S., Zalesky, J., Parmentier, V., Smith, P., Mace, G. N., Mansfield, M., Kempton, E. M.-R., et al. 2021, Nature, 598, 580
  * Liou (2002) Liou, K.-N. 2002, An introduction to atmospheric radiation, Vol. 84 (Academic press)
  * Luger et al. (2021) Luger, R., Foreman-Mackey, D., Hedges, C., & Hogg, D. W. 2021, AJ, 162, 123
  * Mayorga et al. (2021) Mayorga, L., Robinson, T. D., Marley, M. S., May, E., & Stevenson, K. B. 2021, arXiv preprint arXiv:2105.08009
  * Mennesson et al. (2016) Mennesson, B., Gaudi, S., Seager, S., Cahoy, K., Domagal-Goldman, S., Feinberg, L., Guyon, O., Kasdin, J., Marois, C., Mawet, D., et al. 2016, in Space telescopes and instrumentation 2016: Optical, infrared, and millimeter wave, Vol. 9904, International Society for Optics and Photonics, 99040L
  * Murray & Dermott (1999) Murray, C. D. & Dermott, S. F. 1999, Solar system dynamics (Cambridge university press)
  * Nixon et al. (2010) Nixon, C. A., Achterberg, R. K., Romani, P. N., Allen, M., Zhang, X., Teanby, N. A., Irwin, P. G. J., & Flasar, F. M. 2010, Planet. Space Sci., 58, 1667
  * Ohno & Zhang (2019a) Ohno, K. & Zhang, X. 2019a, The Astrophysical Journal, 874, 1
  * Ohno & Zhang (2019b) —. 2019b, The Astrophysical Journal, 874, 2
  * Palubski et al. (2020) Palubski, I. Z., Shields, A. L., & Deitrick, R. 2020, The Astrophysical Journal, 890, 30
  * Parmentier & Guillot (2014) Parmentier, V. & Guillot, T. 2014, Astronomy & Astrophysics, 562, A133
  * Parmentier et al. (2015) Parmentier, V., Guillot, T., Fortney, J. J., & Marley, M. S. 2015, Astronomy & Astrophysics, 574, A35
  * Perez-Becker & Showman (2013) Perez-Becker, D. & Showman, A. P. 2013, ApJ, 776, 134
  * Rauscher (2017) Rauscher, E. 2017, The Astrophysical Journal, 846, 69
  * Schwartz & Cowan (2015) Schwartz, J. C. & Cowan, N. B. 2015, Monthly Notices of the Royal Astronomical Society, 449, 4192
  * Schwartz et al. (2016) Schwartz, J. C., Sekowski, C., Haggard, H. M., Pallé, E., & Cowan, N. B. 2016, Monthly Notices of the Royal Astronomical Society, 457, 926
  * Shields et al. (2016) Shields, A. L., Barnes, R., Agol, E., Charnay, B., Bitz, C., & Meadows, V. S. 2016, Astrobiology, 16, 443
  * Showman & Guillot (2002) Showman, A. P. & Guillot, T. 2002, A&A, 385, 166
  * Showman et al. (2015) Showman, A. P., Lewis, N. K., & Fortney, J. J. 2015, The Astrophysical Journal, 801, 95
  * Showman et al. (2019) Showman, A. P., Tan, X., & Zhang, X. 2019, ApJ, 883, 4
  * Spergel et al. (2013) Spergel, D., Gehrels, N., Breckinridge, J., Donahue, M., Dressler, A., Gaudi, B., Greene, T., Guyon, O., Hirata, C., Kalirai, J., et al. 2013, arXiv preprint arXiv:1305.5422
  * Sromovsky et al. (1998) Sromovsky, L. A., Collard, A. D., Fry, P. M., Orton, G. S., Lemmon, M. T., Tomasko, M. G., & Freedman, R. S. 1998, J. Geophys. Res., 103, 22929
  * Su & Lai (2021) Su, Y. & Lai, D. 2021, Monthly Notices of the Royal Astronomical Society, 509, 3301
  * Ward (1974) Ward, W. R. 1974, J. Geophys. Res., 79, 3375
  * Williams & Kasting (1997) Williams, D. M. & Kasting, J. F. 1997, Icarus, 129, 254
  * Williams & Pollard (2002) Williams, D. M. & Pollard, D. 2002, International Journal of Astrobiology, 1, 61
  * Williams (1978) Williams, G. P. 1978, Journal of the Atmospheric Sciences, 35, 1399
  * Wong et al. (2021) Wong, I., Shporer, A., Zhou, G., Kitzmann, D., Komacek, T. D., Tan, X., Tronsgaard, R., Buchhave, L. A., Vissapragada, S., Greklek-McKeon, M., Rodriguez, J. E., Ahlers, J. P., Quinn, S. N., Furlan, E., Howell, S. B., Bieryla, A., Heng, K., Knutson, H. A., Collins, K. A., McLeod, K. K., Berlind, P., Brown, P., Calkins, M. L., de Leon, J. P., Esparza-Borges, E., Esquerdo, G. A., Fukui, A., Gan, T., Girardin, E., Gnilka, C. L., Ikoma, M., Jensen, E. L. N., Kielkopf, J., Kodama, T., Kurita, S., Lester, K. V., Lewin, P., Marino, G., Murgas, F., Narita, N., Pallé, E., Schwarz, R. P., Stassun, K. G., Tamura, M., Watanabe, N., Benneke, B., Ricker, G. R., Latham, D. W., Vanderspek, R., Seager, S., Winn, J. N., Jenkins, J. M., Caldwell, D. A., Fong, W., Huang, C. X., Mireles, I., Schlieder, J. E., Shiao, B., & Noel Villaseñor, J. 2021, AJ, 162, 256

##  Appendix A Analytic solution of semi-grey perturbations

The plane-parallel, two-stream approximation of radiative transfer equations
for the diffuse, azimuthally averaged IR intensity $I$ are (e.g., Liou 2002
Chapter 4.6, or Komacek et al. 2017 for the exact setup as here):

$\begin{split}\mu\frac{dI^{+}}{d\tau}=I^{+}-B(\tau),\\\
-\mu\frac{dI^{-}}{d\tau}=I^{-}-B(\tau),\end{split}$ (A1)

where $I^{+}$ is the upward intensity, $I^{-}$ is the downward intensity in
the IR, $\tau$ is IR optical depth, $\mu$ is the mean IR zenith angle
associated with a semi-isotropic hemisphere of radiation, and $B(\tau)$ is the
Planck function at the local temperature $T(\tau)$. The diffusive thermal
fluxes are obtained as $F^{\pm}=\pi I^{\pm}$. We use the hemi-isotropic
closure for the thermal band with $\mu=0.5$, which assumes that the thermal
radiation is isotropic in each hemisphere which is reasonable for atmospheres
without much scattering. The stellar irradiation in the visible band is
included simply as

$F^{-}_{v}=F_{0}(1-A)\exp{\left(-\frac{\tau_{v}}{\mu_{v}}\right)},\quad
F^{+}_{v}=0,$ (A2)

where $F_{0}$ is the top-of-atmosphere (TOA) irradiative flux density at a
given location with $\mu_{v}$ being the mean local zenith angle of the
irradiation, and $\tau_{v}$ is the optical depth in the visible band. The
total atmospheric net flux is $F_{\rm tot}=F^{+}_{v}-F^{-}_{v}+F^{+}-F^{-}$.
For simplicity, we consider absorption only and explicitly omit scattering.
Effects of scattering in the energy balance is implicitly included as non-zero
bond albedo $A$. The response of temperature to the time-dependent radiative
fluxes is

$\frac{dT}{dt}=\frac{g}{c_{p}}\frac{dF_{\rm
tot}}{dp}=\frac{\kappa}{c_{p}}\frac{dF_{\rm tot}}{d\tau}.$ (A3)

Perturbations of temperature and thermal fluxes relative to the radiative
equilibrium as a function of $\tau$ are denoted as $T^{\prime}$, $F^{\prime+}$
and $F^{\prime-}$. The perturbation of the Planck function is approximated as
$B^{\prime}=4\sigma T^{3}_{\rm RE}T^{\prime}$. We write the perturbation of
the irradiation as the following periodic oscillation:

$F^{\prime-}_{v}(t)=\Delta F_{v}\exp{(i\omega
t)}\exp{\left(-\frac{\tau_{v}}{\mu_{v}}\right)},$ (A4)

in which the change of $\mu_{v}$ is neglected. Substituting these
perturbations to Eqs. (A1), (A2) and (A3) and subtracting the radiative
equilibrium state, we obtain the equation for the perturbation of the net
thermal flux $F^{\prime}=F^{\prime+}-F^{\prime-}$:

$\begin{split}&\mu^{2}\frac{\partial^{3}F^{\prime}}{\partial\tau^{2}\partial
t}=\frac{\partial F^{\prime}}{\partial t}-\\\ &8\pi\sigma T^{3}_{\rm
RE}\mu\frac{\kappa}{c_{p}}\left(\frac{\partial^{2}F^{\prime}}{\partial\tau^{2}}-\Delta
F_{v}\exp{(i\omega t)}\exp{(-\beta\tau)}\beta^{2}\right),\end{split}$ (A5)

where $\beta=\kappa_{v}/(\kappa\mu_{v})$ and $\kappa_{v}$ is the opacity at
the visible band. For analytic simplicity, we have neglected the variation of
$T_{\rm RE}$ with $\tau$ in Eq. (A5) and adopt $T_{\rm RE}\sim
T_{\mathrm{eq}}$ below. Applying boundary conditions of vanishing
perturbations at $\tau\rightarrow\infty$ and $F^{\prime-}=0$ at $\tau=0$, the
solution to Eq. (A5) is

$F^{\prime}=\left[C\exp{(-k\tau)}+Z\exp{(-\beta\tau)}\right]\Delta
F_{v}\exp{(i\omega t)},$ (A6)

with constants

$\begin{split}&C=-Z\frac{k+\beta\mu k}{\beta+\beta\mu k},\quad
Z=\frac{\beta^{2}}{\beta^{2}+i(\mu\beta^{2}-1/\mu)/\eta},\\\
&k=(\mu^{2}+\eta^{2})^{-1/4}\mu^{-1/2}\exp{(i\gamma/2)},\\\
&\cos\gamma=\frac{\mu}{\sqrt{\mu^{2}+\eta^{2}}},\quad\sin\gamma=\frac{\eta}{\sqrt{\mu^{2}+\eta^{2}}},\\\
&\eta=P_{\rm orbit}4\sigma T_{\mathrm{eq}}^{3}\kappa/c_{p}=\frac{P_{\rm
orbit}}{\tau_{\rm{rad}}}.\end{split}$

The TOA outgoing thermal flux perturbation is
$F^{\prime}(\tau=0,t)=(C+Z)\Delta F_{v}\exp(i\omega t)$, and the temperature
perturbation at an arbitrary $\tau$ is

$\begin{split}&T^{\prime}(\tau,t)=\frac{\kappa\Delta F_{v}}{i\omega
c_{p}}\exp(i\omega t)\times\\\ &[-kC\exp(-k\tau)-\beta
Z\exp(-\beta\tau)+\beta\exp(-\beta\tau)].\end{split}$ (A7)

Figure 8.— Analytic calculations of semi-grey perturbations with different
$\eta=P_{\rm orbit}/\tau_{\rm{rad}}$. All cases assume an equilibrium
temperature of $T_{\rm eq}=250$ K, a bond albedo $A=0.3$, a thermal opacity of
$\kappa=5\times 10^{-4}\;{\rm m^{2}kg^{-1}}$, a visible opacity of
$\kappa_{v}=2\times 10^{-4}\;{\rm m^{2}kg^{-1}}$, $\mu_{v}=1$, and the
specific heat of $c_{p}=13000\;{\rm Jkg^{-1}K^{-1}}$. _Panel (a):_ radiative
timescale calculated as $|T^{\prime}|/|dT_{\rm th}/dt|$ where $dT_{\rm
th}/dt=\frac{\kappa}{c_{p}}\frac{dF^{\prime}}{d\tau}$ and $F^{\prime}$ is
given by Eq. (A6) as a function of optical depth $\tau$ for a few cases with
different $\eta$. The grey dotted line remarks the value of
$\frac{c_{p}}{4\sigma\kappa T_{\mathrm{eq}}^{3}}$ and the grey dashed line
represents a linear relation $\frac{\tau c_{p}}{4\sigma\kappa
T_{\mathrm{eq}}^{3}}$. _Panel (b):_ temperature perturbations, that is the
real part of $T^{\prime}$, at phase 0 ($\omega t=2n\pi$ where $n$ is an
integer) as a function of $\tau$ with different $\eta$. These perturbations
are greatly amplified and normalized so it is possible to visualize. The
profiles with $\eta=0.1,1$ and 5 are multiplied by $\exp{(0.39\tau)}$ and the
profile with $\eta=12.6$ is multiplied by $\exp{(0.27\tau)}$. _Panel (c):_
phase lags of $T^{\prime}$ as a function of $\tau$. _Panel (d), (e) and (f):_
temperature perturbations as a function of $\tau$ at different phases for
cases with different $\eta$. These perturbations are amplified and normalized
as mentioned for panel (b).

For this particular problem we may naturally define the radiative timescale as
$|T^{\prime}|/|dT_{\rm th}/dt|$, where $dT_{\rm
th}/dt=\frac{\kappa}{c_{p}}\frac{dF^{\prime}}{d\tau}$ is the heating rate
caused solely by the perturbations of thermal radiation $F^{\prime}$ given by
Eq. (A6). Panel (a) in Figure 8 shows the radiative timescales as a function
of $\tau$ for a few different $\eta=P_{\rm orbit}/\tau_{\rm{rad}}$. All cases
assume parameters the same as those used in Section 3.1, including an
equilibrium temperature of $T_{\rm eq}=250$ K, a bond albedo $A=0.3$, a
thermal opacity of $\kappa=5\times 10^{-4}\;{\rm m^{2}kg^{-1}}$, a visible
opacity of $\kappa_{v}=2\times 10^{-4}\;{\rm m^{2}kg^{-1}}$, $\mu_{v}=1$, and
the specific heat of $c_{p}=13000\;{\rm Jkg^{-1}K^{-1}}$. The radiative
timescales are generally small at low $\tau$ and increase rapidly with $\tau$
toward the thermal photosphere. But they stop increasing (some even decrease)
and then converge to some values at greater depth. Interestingly, despite
having the same atmospheric parameters, the profiles of the radiative
timescale can quantitatively differ depending on the frequency of the seasonal
irradiation. With this set of atmospheric parameters, the expression
$\tau_{\rm{rad}}\sim\frac{p}{g}\frac{c_{p}}{4\sigma T_{\mathrm{eq}}^{3}}$, or
$\tau_{\rm{rad}}\sim\frac{\tau c_{p}}{4\sigma\kappa T_{\mathrm{eq}}^{3}}$
(represented as the grey dashed line in panel (a)) if assuming a constant
$\kappa$, that has been commonly used in literature (e.g., Showman & Guillot
2002) is a rather good approximation for those calculated by
$|T^{\prime}|/|dT_{\rm th}/dt|$ near the thermal photosphere. The value at
$\tau=1$, $\frac{c_{p}}{4\sigma\kappa T_{\mathrm{eq}}^{3}}$, which is plotted
as the dotted grey line in panel (a), crosses the profiles calculated by
$|T^{\prime}|/|dT_{\rm th}/dt|$ in between $\tau=1$ to $1.4$ for different
cases. These comparisons strongly support the use of
$\tau_{\rm{rad}}\sim\frac{c_{p}}{4\sigma\kappa T_{\mathrm{eq}}^{3}}$ in
Section 2.

We discuss some interesting features of the deep perturbations by the seasonal
irradiation. These perturbations are small and difficult to observe. In
addition, temperature anomalies could be dominated by other sources including
dynamics and diabatic effects from condensation/evaporation in the deep
layers. Therefore, the following discussion is mostly for theoretical interest
within our framework. Panel (b) in Figure 8 shows the temperature
perturbations, the real parts of $T^{\prime}$, at phase 0 ($\omega t=2n\pi$
where $n$ is an integer) for the same set of cases. These perturbations are
greatly amplified as well as normalized so that it is possible to visualize.
The oscillations at the deep layers of cases with low and moderate $\eta$
exhibit a mode with a single long vertical wavelength. A higher $\eta=5$ shows
a mixture of the long mode and a short mode. The case with $\eta=12.6$ is
dominated by the mode with short vertical wavelengths, showing many
alternating perturbations at depth. Panel (c) shows the phase lags of
$T^{\prime}$ as a function of $\tau$ for these cases. The phase lags of all
cases are relatively large at low $\tau$ and decrease with increasing $\tau$.
These phase lags all go across $-90^{\circ}$ from low to high $\tau$ and this
crossing occurs at higher $\tau$ for higher $\eta$. After the $-90^{\circ}$
crossing, the phase lag with $\eta=0.1$ quickly converges to nearly
$-90^{\circ}$ at high optical depth; the one with $\eta=1$ converges to a
slightly smaller value at depth; the one with $\eta=5$ first shows
oscillations then converges to a value $<-100^{\circ}$ at depth; the one with
$\eta=12.6$ continuously decreases with increasing $\tau$.

Time evolution of the deep layers is depicted in panels (d), (e) and (f) in
Figure 8, which show the temperature perturbations at different phases for
cases with $\eta=0.1,5$ and 12.6. In the case with a low $\eta=0.1$, the long
vertical mode appears to oscillate left and right. In the case with a high
$\eta=12.6$, the short vertical mode appears to propagate downward with time.
The case with a medium $\eta=5$ shows a mixture of the long mode that
oscillates left and right and a short mode that propagates downward. Note that
the numerical inversion of the phase function generates phase lags within the
range of ($-\pi,\pi$). For the case with $\eta=12.6$ (or other cases with
$\eta\gtrsim 10$), the downward propagation of the short mode implies that the
phase lag will continuously decrease with increasing $\tau$ rather than
oscillating in between ($-\pi,\pi$). In the phase lag for $\eta=12.6$ shown in
panel (c), with increasing $\tau$, once the phase lag jumps to a value that is
$2\pi$ greater than the adjacent value, it is subtracted by $2\pi$ to ensure
continuity.

Understanding of these behaviours can be obtained from limiting cases of the
solution Eq. (A7). With a general $\eta/\mu$, at phase 0, temperature
perturbation is positive at low optical depth and encounters a first
transition to a negative value at layers with $\tau\gtrsim 1$ which also comes
with a phase lag transition to $<-90^{\circ}$. The critical $\tau$ where this
transition occurs is larger for case with higher $\eta$. This remarks the
transition from a somewhat direct response to irradiation to a thermally
diffusive response at higher optical depth where the irradiation cannot reach
and the oscillation is maintained by the diffusive thermal flux. The phase lag
that is smaller than $-90^{\circ}$ cannot be obtained in the relaxation method
discussed in Section 2, which represents a qualitative difference between the
relaxation and radiative transfer calculations. Whether the deep oscillation
is dominated by the long or short mode then depends on properties of the
solution at $\tau\gg 1$.

When $\eta/\mu\ll 1$, $\cos{\gamma}\sim 1$, $\sin{\gamma}\sim\eta/\mu$,
$k\sim(1+i\eta/(2\mu))/\mu$, $|Z|\ll 1$; because $|k|\sim 1/\mu>\beta$ for
parameters used in Figure 8, $T^{\prime}$ is dominated by terms of
$\exp{(-\beta\tau)}$ at $\tau\gg 1$:

$T^{\prime}\sim\frac{\kappa\Delta F_{v}}{i\omega
c_{p}}\beta(1-Z)\exp(-\beta\tau)\exp(i\omega t).$ (A8)

Therefore, when $\eta/\mu\ll 1$, the real part of $T^{\prime}$ does not change
sign and $T^{\prime}$ approaches a phase lag of nearly $-90^{\circ}$ at
$\tau\gg 1$, corresponding to the long mode shown in Figure 8.

In the other limit of $\eta/\mu\gg 1$, $\gamma\sim\pi/2$,
$k\sim(\eta\mu)^{-1/2}(1/\sqrt{2}+i/\sqrt{2})$. As $|k|<\beta$, at high
optical depth, $T^{\prime}$ is dominated by

$T^{\prime}\sim-\frac{\kappa\Delta F_{v}}{i\omega
c_{p}}kC\exp(-k\tau)\exp(i\omega t),$ (A9)

in which the nonzero imaginary part of $k$ generates oscillations as a
function of $\tau$, corresponding to the downward propagation of short
vertical modes. The vertical wavelength is larger with larger $\eta$.

For a moderate $\eta/\mu$, as long as $|k|>\beta$, the long mode given by Eq.
(A8) still holds at $\tau\gg 1$ with a phase lag converging to values
$<-90^{\circ}$. But at medium $\tau$, the term with $\exp(-k\tau)$ can still
contribute short modes in addition to the long mode, and the typical behavior
is a mixture of a long and short modes as shown for the case with $\eta=5$ in
Figure 8. With the set of parameters used here, the transition where
$|k|<\beta$ occurs at $\eta\gtrsim 10$, explaining the transition from
$\eta=5$ to $\eta=12.6$ in Figure 8.

The physical intuition is that in order to excite the short vertical modes,
the oscillation period of the irradiation has to exceed the radiative
timescale slightly below the thermal photosphere which is typically a factor
of several of $\frac{c_{p}}{4\sigma\kappa T_{\mathrm{eq}}^{3}}$ for parameters
used for Figure 8. The significant vertical variation of radiative timescale
near the thermal photosphere generates a large vertical phase lag variation
over a short vertical distance. Only when the irradiation varies slower than
the local radiative timescale, the short vertical phase variation has time to
propagate downward diffusively without being disrupted by the irradiation.

In cases with a higher $\beta=\kappa_{v}/(\kappa\mu_{v})$ which can be
achieved by strong visible opacity, $|k|<\beta$ may be always satisfied and
there will be short vertical modes at $\tau\gg 1$ with arbitrary $\eta$. In
such cases the calculated radiative timescales $|T^{\prime}|/|dT_{\rm th}/dt|$
show overall smaller values for smaller $\eta$.

##  Appendix B Analytic solution of non-grey perturbations

The radiative transfer equations with absorption only for the (picket fence)
two characteristic thermal bands are:

$\begin{split}\mu\frac{dI_{1}^{+}}{d\tau_{1}}&=I_{1}^{+}-(1-\chi)B(\tau_{1}),\\\
-\mu\frac{dI_{1}^{-}}{d\tau_{1}}&=I_{1}^{-}-(1-\chi)B(\tau_{1}),\\\
\mu\frac{dI_{2}^{+}}{d\tau_{2}}&=I_{2}^{+}-\chi B(\tau_{2}),\\\
-\mu\frac{dI_{2}^{-}}{d\tau_{2}}&=I_{2}^{-}-\chi B(\tau_{2}),\end{split}$ (B1)

where $\chi$ is the fractional band width associated with $\kappa_{2}$ with
respect to the total band width. Similar to Appendix A, we solve for the
structure of perturbations around the radiative equilibrium. Then we consider
partitioning of irradiative energy into multiple channels in the visible band
with different opacities $\kappa_{v,j}$ where $j=1,2,...$. The perturbation of
the irradiation in this case is treated as

$\begin{split}F^{\prime-}_{v}(t)&=\sum_{j}\Delta
F_{v}\epsilon_{j}\exp{(i\omega
t)}\exp{\left(-\frac{\tau_{v,j}}{\mu_{v}}\right)}\\\ &=\Delta
F_{v}\exp{(i\omega
t)}\sum_{j}\epsilon_{j}\exp(-\beta_{j}\tau_{1}),\end{split}$ (B2)

where $\beta_{j}=\kappa_{v,j}/(\kappa_{1}\mu_{v})$ and
$\sum_{j}\epsilon_{j}=1$.

Writing the perturbations of the net thermal flux in band 1 and 2 as
$\tilde{F}^{\prime}_{1}=F^{\prime+}_{1}-F^{\prime-}_{1}$ and
$\tilde{F}^{\prime}_{2}=F^{\prime+}_{2}-F^{\prime-}_{2}$, together with
equations (A3) and (B2), the perburbation form of equations (B1) is organized
as

$\begin{split}\mu\frac{d^{3}\tilde{F}^{\prime}_{1}}{d\tau_{1}^{2}dt}=&\frac{1}{\mu}\frac{d\tilde{F}^{\prime}_{1}}{dt}-8\pi\sigma
T_{\rm RE}^{3}(1-\chi)\frac{\kappa_{1}}{c_{p}}\\\
&\left(\frac{d^{2}\tilde{F}^{\prime}_{1}}{d\tau^{2}_{1}}+\frac{d^{2}\tilde{F}^{\prime}_{2}}{d\tau^{2}_{1}}-\Delta
F_{v}\exp(i\omega
t)\sum_{j}[\epsilon_{j}\beta_{j}^{2}\exp(-\beta_{j}\tau_{1})]\right)\\\
\mu\mathcal{R}^{2}\frac{d^{3}\tilde{F}^{\prime}_{2}}{d\tau_{1}^{2}dt}=&\frac{1}{\mu}\frac{d\tilde{F}^{\prime}_{2}}{dt}-8\pi\sigma
T_{\rm RE}^{3}\chi\mathcal{R}\frac{\kappa_{1}}{c_{p}}\\\
&\left(\frac{d^{2}\tilde{F}^{\prime}_{1}}{d\tau^{2}_{1}}+\frac{d^{2}\tilde{F}^{\prime}_{2}}{d\tau^{2}_{1}}-\Delta
F_{v}\exp(i\omega
t)\sum_{j}[\epsilon_{j}\beta_{j}^{2}\exp(-\beta_{j}\tau_{1})]\right),\end{split}$
(B3)

Where $\mathcal{R}=\kappa_{1}/\kappa_{2}$. The time-dependent flux
perturbations can be written as
$\tilde{F}^{\prime}_{1}(\tau_{1},t)=F^{\prime}_{1}(\tau_{1})\exp(i\omega t)$
and $\tilde{F}^{\prime}_{2}(\tau_{1},t)=F^{\prime}_{2}(\tau_{1})\exp(i\omega
t)$, and the above equation set becomes a set of coupled nonhomogeneous
second-order differential equations with constant coefficients (again, we
assume that $T_{\rm RE}$ is a constant over $\tau$ and $T_{\rm RE}=T_{\rm
eq}$). We first seek general solutions to the homogeneous equation set and
then obtain particular solutions for the full set. The complete solution is
the combination of the general solutions and the particular solutions together
with boundary conditions. The homogeneous part of Equations (B3) can be
written as

$\begin{split}\mathcal{A}_{1}\frac{d^{2}F^{\prime}_{1}}{d\tau_{1}^{2}}=F^{\prime}_{1}+\mathcal{A}_{2}\frac{d^{2}F^{\prime}_{2}}{d\tau_{1}^{2}}\\\
\mathcal{B}_{1}\frac{d^{2}F^{\prime}_{2}}{d\tau_{1}^{2}}=F^{\prime}_{1}+\mathcal{B}_{2}\frac{d^{2}F^{\prime}_{1}}{d\tau_{1}^{2}},\end{split}$
(B4)

with constants

$\begin{split}&\mathcal{A}_{1}\equiv\mu^{2}-i\eta\mu(1-\chi),\quad\mathcal{A}_{2}\equiv
i\eta\mu(1-\chi),\\\
&\mathcal{B}_{1}\equiv\mu^{2}\mathcal{R}^{2}-i\eta\mu\chi\mathcal{R},\quad\mathcal{B}_{2}\equiv
i\eta\mu\chi\mathcal{R},\end{split}$

where $\eta=P4\sigma T_{\mathrm{eq}}^{3}\kappa_{1}/c_{p}=P_{\rm
orbit}/\tau_{\rm{rad}}$. The general solution to the coupled equations (B4) is
a linear combination of exponential functions on the form
$F^{\prime}_{1}=a\exp(k\tau_{1})$ and $F^{\prime}_{2}=b\exp(k\tau_{1})$. By
substituting these characteristic forms to equations (B4), we obtain four
values of $k$. Since perturbations should vanish when
$\tau_{1}\rightarrow\infty$, only the set of $k$ with negative real parts is
retained, which is

$k_{1,2}=-\left(\frac{-(\mathcal{A}_{1}+\mathcal{B}_{1})\pm\sqrt{(\mathcal{A}_{1}-\mathcal{B}_{1})^{2}+4\mathcal{A}_{2}\mathcal{B}_{2}}}{2(\mathcal{A}_{2}\mathcal{B}_{2}-\mathcal{A}_{1}\mathcal{B}_{1})}\right)^{1/2}.$
(B5)

The constants are

$b_{1,2}=a_{1,2}\frac{\mathcal{B}_{2}k_{1,2}^{2}}{\mathcal{B}_{1}k_{1,2}^{2}-1},$
(B6)

with $a_{1,2}$ indeterminate.

Then we seek particular solutions to the full equation set (B3) in the form of
$F^{\ast}_{1}\exp(i\omega t)$ and $F^{\ast}_{2}\exp(i\omega t)$, with

$\begin{split}F^{\ast}_{1}=\Delta
F_{v}\sum_{j}[\epsilon_{j}\mathcal{Z}_{1,j}\exp(-\beta_{j}\tau_{1})]\\\
F^{\ast}_{2}=\Delta
F_{v}\sum_{j}(\epsilon_{j}\mathcal{Z}_{2,j}\exp(-\beta_{j}\tau_{1})].\end{split}$
(B7)

Substituting equations (B7) to (B3), we obtain coefficients as follows

$\mathcal{Z}_{1,j}=\frac{\mathcal{A}_{2}\beta_{j}^{2}(1-\mathcal{B}_{1}\beta_{j}^{2}-\mathcal{B}_{2}\beta_{j}^{2})}{(\mathcal{A}_{1}\mathcal{B}_{1}-\mathcal{A}_{2}\mathcal{B}_{2})\beta_{j}^{4}-(\mathcal{A}_{1}+\mathcal{B}_{1})\beta_{j}^{2}+1}$
(B8)

and

$\mathcal{Z}_{2,j}=\frac{\mathcal{B}_{2}\beta_{j}^{2}}{\mathcal{B}_{1}\beta_{j}^{2}-1}(\mathcal{Z}_{1,j}-1).$
(B9)

The complete solution of the time-independent part of the net thermal flux is

$\begin{split}F^{\prime}_{1}=a_{1}\exp(k_{1}\tau_{1})+a_{2}\exp(k_{2}\tau_{1})+F^{\ast}_{1},\\\
F^{\prime}_{2}=b_{1}\exp(k_{1}\tau_{1})+b_{2}\exp(k_{2}\tau_{1})+F^{\ast}_{2},\end{split}$
(B10)

in which $a_{1,2}$ are yet to be determined. We make use the upper boundary
conditions $F^{\prime-}_{1}=0$ and $F^{\prime-}_{2}=0$ to determine $a_{1,2}$,
yielding

$\begin{split}&a_{1}=\frac{\Delta
F_{v}\sum_{j}\left(\epsilon_{j}\mathcal{Z}_{2,j}\frac{\beta_{j}\mathcal{R}\mu+1}{\beta_{j}\mathcal{R}\mu}-\epsilon_{j}\mathcal{Z}_{1,j}\frac{\beta_{j}\mu+1}{\beta_{j}\mu}\frac{\mathcal{B}_{2}k^{2}_{2}}{\mathcal{B}_{1}k^{2}_{2}-1}\frac{1-k_{2}\mu\mathcal{R}}{\mathcal{R}(1-k_{2}\mu)}\right)}{\frac{\mathcal{B}_{2}k_{1}}{\mathcal{B}_{1}k^{2}_{1}-1}\frac{1-k_{1}\mu\mathcal{R}}{\mu\mathcal{R}}-\frac{\mathcal{B}_{2}k^{2}_{2}}{\mathcal{B}_{1}k^{2}_{2}-1}\frac{1-k_{2}\mu\mathcal{R}}{\mathcal{R}(1-k_{2}\mu)}\frac{1-k_{1}\mu}{k_{1}\mu}}\\\
&a_{2}=\frac{\Delta
F_{v}\sum_{j}\left(\epsilon_{j}\mathcal{Z}_{2,j}\frac{\beta_{j}\mathcal{R}\mu+1}{\beta_{j}\mathcal{R}\mu}-\epsilon_{j}\mathcal{Z}_{1,j}\frac{\beta_{j}\mu+1}{\beta_{j}\mu}\frac{\mathcal{B}_{2}k^{2}_{1}}{\mathcal{B}_{1}k^{2}_{1}-1}\frac{1-k_{1}\mu\mathcal{R}}{\mathcal{R}(1-k_{1}\mu)}\right)}{\frac{\mathcal{B}_{2}k_{2}}{\mathcal{B}_{1}k^{2}_{2}-1}\frac{1-k_{2}\mu\mathcal{R}}{\mu\mathcal{R}}-\frac{\mathcal{B}_{2}k^{2}_{1}}{\mathcal{B}_{1}k^{2}_{1}-1}\frac{1-k_{1}\mu\mathcal{R}}{\mathcal{R}(1-k_{1}\mu)}\frac{1-k_{2}\mu}{k_{2}\mu}}.\end{split}$
(B11)

Finally, the desired solution for thermal fluxes is that in Eqs. (B10) times
$\exp(i\omega t)$. Temperature perturbation as a function of time and optical
depth $\tau_{1}$ is derived using

$\frac{dT}{dt}=\frac{g}{c_{p}}\frac{dF_{\rm
tot}}{dp}=\frac{\kappa_{1}}{c_{p}}\frac{dF_{\rm tot}}{d\tau_{1}},$ (B12)

which is given by

$\begin{split}T(\tau_{1},t)&=\frac{\kappa_{1}}{i\omega c_{p}}\frac{\partial
F^{\prime}_{\rm tot}}{\partial\tau_{1}}=\frac{\kappa_{1}}{i\omega
c_{p}}\exp(i\omega t)\times\\\
&\\{a_{1}k_{1}\exp(k_{1}\tau_{1})+a_{2}k_{2}\exp(k_{2}\tau_{1})-\Delta
F_{v}\sum_{j}[\epsilon_{j}\mathcal{Z}_{1,j}\beta_{j}\exp(-\beta_{j}\tau_{1})]\\\
&+b_{1}k_{1}\exp(k_{1}\tau_{1})+b_{2}k_{2}\exp(k_{2}\tau_{1})-\Delta
F_{v}\sum_{j}[\epsilon_{j}\mathcal{Z}_{2,j}\beta_{j}\exp(-\beta_{j}\tau_{1})]\\\
&+\Delta
F_{v}\sum_{j}[\epsilon_{j}\beta_{j}\exp(-\beta_{j}\tau_{1})]\\}.\end{split}$
(B13)

Likewise, these solutions can be expressed as a combination of amplitudes and
phase lags.

The perturbations in the deep layers in the non-grey situation should have
qualitatively similar properties as those in the semi-grey cases shown in
Appendix A. This is because at high optical depth both cases collapse into the
diffusive limit which may be adaptively described by the semi-grey framework.