Weak Seasonality on Temperate Exoplanets Around Low-mass Stars

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 $\ll1$ for planets around very-low-mass M dwarfs and $\gtrsim1$ 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.


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. 2018Fletcher et al. , 2020 and Uranus (see a review by Fletcher 2021). Even Jupiter, which has a negligible obliquity of ∼ 3 • , 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 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 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 locked 1 . 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.

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 τ rad ∼ p g cp 4σT 3 eq , 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 σ 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 τ rad ∼ cp 4κσT 3 eq , where κ is the mean IR opacity. Here, T 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: where r p is the planetary radius, a is the semimajor axis, e is the eccentricity, A is the planetary bond albedo, T is the stellar effective temperature, and R is the stellar radius. The annual-mean stellar flux F with non-zero eccentricity has been well known (e.g., Williams & Pollard 2002). Using Kepler's third law P 2 orbit ≈ 4π 2 GM a 3 where P orbit is the orbital period, M is the stellar mass 2 and G is the gravitational constant, the ratio between 1 Some close-in planets may still have chances to maintain nonzero obliquity through perturbations from additional companions in the system (Su & Lai 2021).
2 When the planetary mass Mp is a nontrivial fraction of the total mass of the system, one should use the orbital period and the radiative timescale is where κ and c p are in SI unit, T , R and M are the effective temperature, radius and mass of the sun which are adopted as 5772 K, 6.957 × 10 8 m and 1.9891 × 10 30 kg, respectively. For a particular set of atmospheric parameters of κ, c p and A, P orbit /τ 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 solarsystem 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 orbit /τ rad .
Given the complications mentioned above, we compute P orbit /τ rad as a function of varying κ/c p and stellar mass. Results are shown in Figure 1 as a function of stellar mass and κ/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 κ ref = 5 × 10 −4 m 2 kg −1 and c p,ref = 13000 Jkg −1 K −1 . This κ ref results in a mean thermal photosphere at about 0.5 bar for Jovian gravity which is reasonable for Jovianlike 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. 3 At a given κ/c p , P orbit /τ rad spans about two orders of magnitude when the stellar mass varies from 0.18 to 1.7 M . For κ/c p near the reference value, P orbit /τ rad is much less than 1 for low-mass M dwarfs and > 1 with stellar masses 0.6 M with moderate bond albedo and eccentricity. At sufficiently high κ/c p which represents high metallicity in the atmosphere, P orbit /τ 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.
A small P orbit /τ 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 where T RE = T RE + ∆T RE cos ωt is the time-varying radiative equilibrium temperature, T RE is the timeindependent component and ω = 2π/P orbit . The timedependent temperature is where sin α = − ωτ rad √ 1+ω 2 τ 2 rad and α is a phase lag of the temperature perturbation relative to the irradiation. In the presence of a finite τ rad , temperature varies with an amplitude a factor of 1/ 1 + ω 2 τ 2 rad of that of the radiative equilibrium, and the phase lag α is in between −90 • and 0. A larger P orbit /τ rad leads to a stronger temperature variation and a smaller magnitude of phase lags α. 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, κ ref , c p,ref and a few values for A. It appears that the seasonal temperature variations are negligible and there is a nearly −90 • phase lag near the photosphere for Jupiter-like planets around stars with masses 0.6 M .

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 of database of stellar properties is not expected to change our qualitative results and conclusions. 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.

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-ofatmosphere (TOA) flux F (τ = 0, t) and photospheric temperature T (τ = 1, t) in the form of where ∆F v is the magnitude of the irradiation variation. Physical quantities are the real parts of Eq. (5). The relative amplitude |F |/∆F v and phase lag α f of the TOA thermal flux perturbation, as well as the phase lag of the temperature perturbation α T are determined only by P orbit /τ rad . The amplitude of temperature perturbations |T | depends on T eq and ∆F v too, but the ratio to a reference amplitude |T |/|T ref | depends on P orbit /τ rad only, where T ref is a perturbation with an extremely large P orbit /τ rad which is arbitrarily set to 10 4 . As an example, we assume κ = κ ref , c p = c p,ref , visible opacity κ v = 2 × 10 −4 m 2 kg −1 , zenith angle of the irradiation µ v = 1, A = 0.5 and e = 0. Time evolution The logarithm of the ratio between the orbital period and the radiative timescale at τ = 1 level, log(P orbit /τ rad ), as a function of stellar mass and κ/cp relative to the reference value (the reference infrared opacity is κ ref = 5 × 10 −4 m 2 kg −1 and specific heat is 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 orbit /τ rad = 1. Discontinuity in the plots arises from the segmental fits in Eker et al. (2018). -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 κ is assumed to be 5 × 10 −4 m 2 kg −1 and the heat capacity cp is assumed to be 13000 Jkg −1 K −1 .
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 (τ = 0) and T (τ = 1) as a function of stellar mass exhibit the same general trend as that in Section 2. These relative amplitudes are 0.1 at low stellar mass ( 0.6 M ), quantitatively similar to those using relaxation method ( Figure 2). However, they increase only up to 0.6 at M = 1.7 M , 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 • at M < 0.2 M . The phase lags in the radiative transfer calculation increase (i.e., the amplitudes decrease) more rapidly with increasing M when M 1 M than that in the relax-ation calculation, but their upper bound is lower at the highest M . 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 κ v > κ. 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.
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 ab-  The picket fence model considers two different values for the thermal opacity κ 1 and κ 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 κ 1 and κ 2 . Therefore, the total (wavelength-integrated) thermal emission through bands with κ 1 and κ 2 are approximated as (1 − χ)B and χB, respectively, where χ is the fractional band width associated with κ 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 κ 1 , κ 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 κ 1 = 5 × 10 −4 m 2 kg −1 , κ 2 = 5 × 10 −3 m 2 kg −1 , χ = 0.3 and three visible bands with κ v = [2 × 10 −4 , 2 × 10 −3 , 2 × 10 −2 ] m 2 kg −1 and partition coefficients of visible bands = [0.7, 0.2, 0.1]. As expected, the temperature response is weak at high optical depth and smoothly increases with decreasing τ 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 • from low to high τ 1 and this crossing occurs at higher τ 1 for cases with higher stellar mass. After the −90 • crossing, the phase lag in the case with 0.3 M quickly converges to nearly −90 • at high optical depth; the one with 1 M first shows oscillations then converges to a value < −100 • at depth; the one with 1.7 M continuously decreases with increasing τ 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 semigrey case in Figure 4, but overall the |F |/∆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 ( 30 • ), whereas they show a maximum of about 60 • in cases with the stellar mass of about 0.6 M . 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 ∼ 0.3, P orbit /τ 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 • and −15 • at 1 mbar and 0.1 mbar, respectively. See Fletcher et al. (2018Fletcher et al. ( , 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.

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 , Rauscher 2017.
The stellar irradiation pattern in the presence of nonzero obliquity and eccentricity is described as follow. The latitude of the substellar point φ ss in our model is given by sin φ ss = sin ψ cos f, where ψ 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): The time evolution of E is given by the Kepler's equation  Different bond albedo A are assumed for different curves, and eccentricity is assumed to be zero for all models. The infrared opacity κ is assumed to be 5 × 10 −4 m 2 kg −1 and the heat capacity cp is assumed to be 13000 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 orbit /τ rad = 10 4 in the upper panel. The bottom panel shows the phase lag of temperature anomalies at τ = 1 relative to the irradiation anomaly. The star-planet distance is d = a(1 − e cos f ). Equation (6) states that the maximum Northern excursion of the substellar latitude (φ ss = ψ) occurs when the planet is at the periapsis (f = 0). Of course, Equation (6) where H = cos −1 (− tan φ tan φ ss ) is the length from starrise to noon at latitude φ 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 π, respectively. Note that d, H and φ 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 κ = 5 × 10 −4 m 2 kg −1 , a visible opacity of κ v = 2 × 10 −4 m 2 kg −1 , and specific heat c p = 13000 Jkg −1 K −1 , the same as those in Section 3.1. We utilize the numerical package twostr (Kylling et al. 1995) to integrate the radiative transfer equations, and the mean zenith angle of irradiation is assumed to be 1/ √ 3 following the default value in 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 eq = 250 K, but with different obliquity and around stars with different masses. The top and middle rows in Figure 7 show the top-ofatmosphere (TOA) thermal flux as a function of time and latitude for cases with stellar masses of M = 0.4, 1 and 1.6 M and two obliquities of 30 • and 60 • . The stellar effective temperature T and radius R are found  using relations in Eker et al. (2018). With the given T eq , the orbital periods for cases with 0.4, 1 and 1.6 M 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 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 • . Cases with an obliquity of 60 • 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 and 60 • 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 • 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 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.

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 τ rad (see their Figure  1). In this study, we show that their regime IV and V with τ rad P orbit can mostly apply to gaseous planets around low-mass stars, while their regime II and III with τ rad P 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 ψ c 54 • 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 ψ c , the smaller the meridional gradient of the longitudinal-annual-mean irradiation. Planets with an obliquity sufficiently close to the critical value ψ 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 X. Tan All models assume an equilibrium temperature of 250 K, a bond albedo A = 0.3, a thermal opacity of κ = 5×10 −4 m 2 kg −1 , a visible opacity of κv = 2 × 10 −4 m 2 kg −1 and the specific heat of cp = 13000 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 • , a zero eccentricity but with three stellar masses of 0.4 (left), 1 (middle) and 1.6 M (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 for cases with M = 0.4 and 1.6 M , respectively. Middle row: same as those at the top row except that the obliquity is 60 • . 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 , respectively. close to ψ 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 ψ c around moderate-mass stars to low-mass stars would be interesting to investigate. On the other hand, planets with obliquities much different from ψ 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 grav-ity 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 hydrogendominated 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(Fletcher et al. , 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 lowand 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/EXOCON-DENSE, PI: R.T. Pierrehumbert).
The plane-parallel, two-stream approximation of radiative transfer equations for the diffuse, azimuthally averaged IR intensity I are (e.g., Liou 2002Chapter 4.6, or Komacek et al. 2017 for the exact setup as here): where I + is the upward intensity, I − is the downward intensity in the IR, τ is IR optical depth, µ is the mean IR zenith angle associated with a semi-isotropic hemisphere of radiation, and B(τ ) is the Planck function at the local temperature T (τ ). The diffusive thermal fluxes are obtained as F ± = πI ± . We use the hemi-isotropic closure for the thermal band with µ = 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 where F 0 is the top-of-atmosphere (TOA) irradiative flux density at a given location with µ v being the mean local zenith angle of the irradiation, and τ v is the optical depth in the visible band. The total atmospheric net flux is 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 Perturbations of temperature and thermal fluxes relative to the radiative equilibrium as a function of τ are denoted as T , F + and F − . The perturbation of the Planck function is approximated as B = 4σT 3 RE T . We write the perturbation of the irradiation as the following periodic oscillation: in which the change of µ 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 = F + − F − : where β = κ v /(κµ v ) and κ v is the opacity at the visible band. For analytic simplicity, we have neglected the variation of T RE with τ in Eq. (A5) and adopt T RE ∼ T eq below. Applying boundary conditions of vanishing perturbations at τ → ∞ and F − = 0 at τ = 0, the solution to Eq. (A5) is with constants The TOA outgoing thermal flux perturbation is F (τ = 0, t) = (C +Z)∆F v exp(iωt), and the temperature perturbation at an arbitrary τ is For this particular problem we may naturally define the radiative timescale as |T |/|dT th /dt|, where dT th /dt = κ cp dF dτ is the heating rate caused solely by the perturbations of thermal radiation F given by Eq. (A6). Panel (a) in Figure  8 shows the radiative timescales as a function of τ for a few different η = P orbit /τ rad . All cases assume parameters the same as those used in Section 3.1, including an equilibrium temperature of T eq = 250 K, a bond albedo A = 0.3, a thermal opacity of κ = 5 × 10 −4 m 2 kg −1 , a visible opacity of κ v = 2 × 10 −4 m 2 kg −1 , µ v = 1, and the specific heat of c p = 13000 Jkg −1 K −1 . The radiative timescales are generally small at low τ and increase rapidly with τ 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 τ rad ∼ p g cp 4σT 3 eq , or τ rad ∼ τ cp 4σκT 3 eq (represented as the grey dashed line in panel (a)) if assuming a constant κ, that has been commonly used in literature (e.g., Showman & Guillot 2002) is a rather good approximation for those calculated by |T |/|dT th /dt| near the thermal photosphere. The value at τ = 1, cp 4σκT 3 eq , which is plotted as the dotted grey line in panel (a), crosses the profiles calculated by |T |/|dT th /dt| in between τ = 1 to 1.4 for different cases. These comparisons strongly support the use of τ rad ∼ cp 4σκT 3 eq 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 , at phase 0 (ωt = 2nπ 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 η exhibit a mode with a single long vertical wavelength. A higher η = 5 shows a mixture of the long mode and a short mode. The case with η = 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 as a function of τ for these cases. The phase lags of all cases are relatively large at low τ and decrease with increasing τ . These phase lags all go across −90 • from low to high τ and this crossing occurs at higher τ for higher η. After the −90 • crossing, the phase lag with η = 0.1 quickly converges to nearly −90 • at high optical depth; the one with η = 1 converges to a slightly smaller value at depth; the one with η = 5 first shows oscillations then converges to a value < −100 • at depth; the one with η = 12.6 continuously decreases with increasing τ .
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 η = 0.1, 5 and 12.6. In the case with a low η = 0.1, the long vertical mode appears to oscillate left and right. In the case with a high η = 12.6, the short vertical mode appears to propagate downward with time. The case with a medium η = 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 (−π, π). For the case with η = 12.6 (or other cases with η 10), the downward propagation of the short mode implies that the phase lag will continuously decrease with increasing τ rather than oscillating in between (−π, π). In the phase lag for η = 12.6 shown in panel (c), with increasing τ , once the phase lag jumps to a value that is 2π greater than the adjacent value, it is subtracted by 2π to ensure continuity.
Understanding of these behaviours can be obtained from limiting cases of the solution Eq. (A7). With a general η/µ, at phase 0, temperature perturbation is positive at low optical depth and encounters a first transition to a negative value at layers with τ 1 which also comes with a phase lag transition to < −90 • . The critical τ where this transition occurs is larger for case with higher η. 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 • 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 τ 1. When η/µ 1, cos γ ∼ 1, sin γ ∼ η/µ, k ∼ (1 + iη/(2µ))/µ, |Z| 1; because |k| ∼ 1/µ > β for parameters used in Figure 8, T is dominated by terms of exp (−βτ ) at τ 1: Therefore, when η/µ 1, the real part of T does not change sign and T approaches a phase lag of nearly −90 • at τ 1, corresponding to the long mode shown in Figure 8. In the other limit of η/µ 1, γ ∼ π/2, k ∼ (ηµ) −1/2 (1/ √ 2 + i/ √ 2). As |k| < β, at high optical depth, T is dominated by in which the nonzero imaginary part of k generates oscillations as a function of τ , corresponding to the downward propagation of short vertical modes. The vertical wavelength is larger with larger η. For a moderate η/µ, as long as |k| > β, the long mode given by Eq. (A8) still holds at τ 1 with a phase lag converging to values < −90 • . But at medium τ , the term with exp(−kτ ) 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 η = 5 in Figure 8. With the set of parameters used here, the transition where |k| < β occurs at η 10, explaining the transition from η = 5 to η = 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 cp 4σκT 3 eq 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 β = κ v /(κµ v ) which can be achieved by strong visible opacity, |k| < β may be always satisfied and there will be short vertical modes at τ 1 with arbitrary η. In such cases the calculated radiative timescales |T |/|dT th /dt| show overall smaller values for smaller η.

B. ANALYTIC SOLUTION OF NON-GREY PERTURBATIONS
The radiative transfer equations with absorption only for the (picket fence) two characteristic thermal bands are: where χ is the fractional band width associated with κ 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 κ v,j where j = 1, 2, .... The perturbation of the irradiation in this case is treated as where β j = κ v,j /(κ 1 µ v ) and j j = 1. Writing the perturbations of the net thermal flux in band 1 and 2 asF 1 = F + 1 − F − 1 andF 2 = F + 2 − F − 2 , together with equations (A3) and (B2), the perburbation form of equations (B1) is organized as Where R = κ 1 /κ 2 . The time-dependent flux perturbations can be written asF 1 (τ 1 , t) = F 1 (τ 1 ) exp(iωt) and F 2 (τ 1 , t) = F 2 (τ 1 ) exp(iωt), and the above equation set becomes a set of coupled nonhomogeneous second-order differential equations with constant coefficients (again, we assume that T RE is a constant over τ and T RE = T 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 with constants where η = P 4σT 3 eq κ 1 /c p = P orbit /τ rad . The general solution to the coupled equations (B4) is a linear combination of exponential functions on the form F 1 = a exp(kτ 1 ) and F 2 = b exp(kτ 1 ). By substituting these characteristic forms to equations (B4), we obtain four values of k. Since perturbations should vanish when τ 1 → ∞, only the set of k with negative real parts is retained, which is The constants are b 1,2 = a 1,2 B 2 k 2 1,2 with a 1,2 indeterminate. Then we seek particular solutions to the full equation set (B3) in the form of F * 1 exp(iωt) and F * 2 exp(iωt), with Substituting equations (B7) to (B3), we obtain coefficients as follows Z 1,j = A 2 β 2 j (1 − B 1 β 2 j − B 2 β 2 j ) (A 1 B 1 − A 2 B 2 )β 4 j − (A 1 + B 1 )β 2 j + 1 (B8) and Z 2,j = B 2 β 2 j B 1 β 2 j − 1 (Z 1,j − 1).
The complete solution of the time-independent part of the net thermal flux is F 1 = a 1 exp(k 1 τ 1 ) + a 2 exp(k 2 τ 1 ) + F * 1 , F 2 = b 1 exp(k 1 τ 1 ) + b 2 exp(k 2 τ 1 ) + F * 2 , in which a 1,2 are yet to be determined. We make use the upper boundary conditions F − 1 = 0 and F − 2 = 0 to determine a 1,2 , yielding a 1 = ∆F v j j Z 2,j Finally, the desired solution for thermal fluxes is that in Eqs. (B10) times exp(iωt). Temperature perturbation as a function of time and optical depth τ 1 is derived using which is given by T (τ 1 , t) = κ 1 iωc p ∂F tot ∂τ 1 = κ 1 iωc p exp(iωt)× {a 1 k 1 exp(k 1 τ 1 ) + a 2 k 2 exp(k 2 τ 1 ) − ∆F v j [ j Z 1,j β j exp(−β j τ 1 )] 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.