# Obliquity of an Earth-like Planet from Frequency Modulation of Its Direct-imaged Lightcurve: Mock Analysis from General Circulation Model Simulation

Yuta Nakagawa, Takanori Kodama, Masaki Ishiwatari, Hajime Kawahara, Yasushi Suto, Yoshiyuki O. Takahashi, George L. Hashimoto, Kiyoshi Kuramoto, Kensuke Nakajima, Shin-Ichi Takehiro, et al.

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

## ► To cite this version:

Yuta Nakagawa, Takanori Kodama, Masaki Ishiwatari, Hajime Kawahara, Yasushi Suto, et al.:: Obliquity of an Earth-like Planet from Frequency Modulation of Its Direct-imaged Lightcurve: Mock Analysis from General Circulation Model Simulation. The Astrophysical Journal, 2020, 898, 10.3847/1538-4357/ab9eb8. insu-03678865

HAL Id: insu-03678865

https://insu.hal.science/insu-03678865v1

Submitted on 13 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Obliquity of an Earth-like Planet from Frequency Modulation of Its Direct-imaged Lightcurve: Mock Analysis from General Circulation Model Simulation

Yuta Nakagawa<fnref n="1" />, Takanori Kodama<fnref n="2" />, Masaki Ishiwatari<fnref n="3" />,<fnref n="4" />, Hajime Kawahara<fnref n="5" />,<fnref n="6" />, Yasushi Suto<fnref n="1" />,<fnref n="6" />, Yoshiyuki O. Takahashi<fnref n="4" />,<fnref n="7" />, George L. Hashimoto<fnref n="8" />, Kiyoshi Kuramoto<fnref n="3" />,<fnref n="4" />, Kensuke Nakajima<fnref n="9" />, Shin-ichi Takehiro<fnref n="10" />, and Yoshi-Yuki Hayashi<fnref n="4" />,<fnref n="7" />

Received 2019 November 29; revised 2020 June 19; accepted 2020 June 19; published 2020 July 28

## Abstract

Direct-imaging techniques of exoplanets have made significant progress recently and will eventually enable monitoring of photometric and spectroscopic signals of Earth-like habitable planets. The presence of clouds, however, would remain as one of the most uncertain components in deciphering such direct-imaged signals of planets. We attempt to examine how the planetary obliquity produces different cloud patterns by performing a series of general circulation model simulation runs using a set of parameters relevant for our Earth. Then we use the simulated photometric lightcurves to compute their frequency modulation that is due to the planetary spin-orbit coupling over an entire orbital period, and we attempt to see to what extent one can estimate the obliquity of an Earth twin. We find that it is possible to estimate the obliquity of an Earth twin within the uncertainty of several degrees with a dedicated 4 m space telescope at 10 pc away from the system if the stellar flux is completely blocked. While our conclusion is based on several idealized assumptions, a frequency modulation of a directly imaged Earth-like planet offers a unique methodology to determine its obliquity.

*Unified Astronomy Thesaurus concepts:* Exoplanet atmospheres (487); Exoplanet surface composition (2022); Hydrodynamics (1963); Radiative transfer (1335); Oblique rotators (1144); Biosignatures (2018); Direct imaging (387); Stellar photometry (1620); Earth atmosphere (437); Atmospheric circulation (112)

## 1. Introduction

Direct imaging of Earth-like planets is a quite challenging but indispensable technique to revolutionize our understanding of planets in the near future. The amplitude modulation of a photometric lightcurve from a *color-changing dot* is sensitive to its surface pattern and thus would reveal the presence of lands, oceans, clouds, and even vegetation on the surface of the planets (e.g., Sagan et al. 1993; Ford et al. 2001; Cowan et al. 2009; Oakley & Cash 2009; Fujii et al. 2010, 2011; Rushby et al. 2019; Suto 2019). Indeed, continuous monitoring of oblique planets over their orbital periods may even enable one to reconstruct their two-dimensional surface map (Kawahara & Fujii 2010, 2011; Fujii & Kawahara 2012; Farr et al. 2018). The feasibility of the mapping has recently been tested using continuous Earth observations by the Deep Space Climate Observatory orbiting at an altitude of 150 km (Jiang et al. 2018; Fan et al. 2019; Aizawa et al. 2020).

The lightcurve carries complementary information for the planet as well. The autocorrelation analysis of the photometric variation roughly provides the rotation period of the planet (Pallé et al. 2008). The obliquity can also be inferred from a simultaneous fitting of the spin vector and planet surface (e.g., Kawahara & Fujii 2010; Schwartz et al. 2016; Farr et al. 2018). Such dynamical parameters of the planet are of interest for a general circulation modeling of Earth-like

planets (e.g., Kaspi & Showman 2015; Deitrick et al. 2018; Komacek et al. 2019).

Strictly speaking, an apparent photometric period observed by a distant observer is not necessarily identical to the true spin rotation period that is due to the planetary orbital motion. This is related to why a sidereal day of our Earth $P_{\text{spin}}$ is approximately 365.24/366.24 × 24 ≈ 23.934 hr, which corresponds to the true spin frequency $f_{\text{spin}}$ ≈ 1.00274 (day⁻<fnref n="1" />), instead of the $f_{\text{spin,heliocentric}} = 1$ (day⁻<fnref n="1" />). The difference between the observed and true spin rotation frequencies, $f_{\text{obs}}$ and $f_{\text{spin}}$, is time dependent and sensitive to the geometrical configuration of the system, including the planetary obliquity, $\zeta$, the inclination of the planetary orbital plane for the observer, $i$, and the observer's direction (the orbital phase angle $\Theta_{\text{eq}}$ measured from the ascending node, for instance).

Thus the corresponding frequency modulation of the periodicity in the lightcurve may reveal those parameters, through the presence of the large-scale inhomogeneity of the surface. We emphasize that the frequency modulation signal is much less sensitive to the specific distribution pattern of the surface than the amplitude modulation is. Kawahara (2016, hereafter K16) proposed a novel idea to measure the planetary obliquity from the frequency modulation and successfully demonstrated its feasibility using a static cloud-subtracted Earth model.

<footnote n="1">¹Department of Physics, The University of Tokyo, Tokyo 113-0033, Japan; suto@phys.s.u-tokyo.ac.jp</footnote>

<footnote n="2">²Laboratoire d'astrophysique de Bordeaux, Université de Bordeaux, CNRS, Allée Geoffroy Saint-Hilaire, F-33615 Pessac, France</footnote>

<footnote n="3">³Division of Earth and Planetary Sciences, Hokkaido University, Sapporo 060-0810, Japan</footnote>

<footnote n="4">⁴Center for Planetary Science, Kobe University, Kobe 650-0047, Japan</footnote>

<footnote n="5">⁵Department of Earth and Planetary Sciences, The University of Tokyo, Tokyo 113-0033, Japan</footnote>

<footnote n="6">⁶Research Center for the Early Universe, School of Science, The University of Tokyo, Tokyo 113-0033, Japan</footnote>

<footnote n="7">⁷Department of Planetology, Kobe University, Kobe 657-8501, Japan</footnote>

<footnote n="8">⁸Department of Earth Sciences, Okayama University, Okayama 700-8530, Japan</footnote>

<footnote n="9">⁹Department of Earth and Planetary Sciences, Kyushu University, Fukuoka 819-0395, Japan</footnote>

<footnote n="10">¹⁰Research Institute for Mathematical Sciences, Kyoto University, Kyoto 606-8502, Japan</footnote>
<figure><img src="image_4.png" /></figure>

$$P_{\text{obs}} = P_{(\text{lspin} + (R \to R'))}$$

$$2\pi f_{\text{spin}} P_{\text{obs}} = 2\pi + 2\pi f_{\text{orb}} P_{\text{obs}}$$

$$\therefore f_{\text{obs}} \equiv \frac{1}{P_{\text{obs}}} = f_{\text{spin}} - f_{\text{orb}}$$

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

$$P_{\text{obs}} = P_{(\text{lspin} - (R' \to R))}$$

$$2\pi f_{\text{spin}} P_{\text{obs}} = 2\pi - 2\pi f_{\text{orb}} P_{\text{obs}}$$

$$\therefore f_{\text{obs}} = f_{\text{spin}} + f_{\text{orb}}$$

**Figure 1.** Schematic illustration for the periodicity in the photometric lightcurve of a planet observed from the direction $i = 0^{\circ}$ relative to the normal vector of the orbital plane. Panels (a) and (b) indicate prograde ($\zeta = 0^{\circ}$) and retrograde ($\zeta = 180^{\circ}$) planets, respectively, in a *geocentric* frame. The host star moves around the planet from $S$ to $S'$ after one heliocentric day of the planet. The illuminated and visible part of the planet for a face-on observer (shaded region) changes accordingly, and the reflective point at each epoch also moves from $R$ to $R'$.

The basic principle of frequency modulation can be understood from Figure 1. For a perfectly prograde planet ($\zeta = 0^{\circ}$), the illuminated and visible part of the planet viewed from a face-on observer ($i = 0^{\circ}$) moves in the same direction as the planetary spin (panel (a)). The reflective point, at which the reflected flux of the star is maximal on the planetary surface, moves accordingly, and thus it takes slightly more than one spin rotation period $P_{\text{spin}}$ for the observer to see the exact same part of the planet. Therefore, the observed photometric variation frequency becomes $f_{\text{obs}} = f_{\text{spin}} - f_{\text{orb}}$. Applying the same argument, one can easily understand that $f_{\text{obs}} = f_{\text{spin}} + f_{\text{orb}}$ for a perfectly retrograde planet ($\zeta = 180^{\circ}$), as illustrated in panel (b) of Figure 1.

In general, the photometric variation frequency $f_{\text{obs}}$ is not constant and varies according to the mutual geometry between the star and the planet, leading to a frequency modulation of the photometric lightcurve of the planet. Figure 2 illustrates an example of the time-dependent frequency modulation for a $\zeta = 90^{\circ}$ planet viewed by a distant observer at $i = 0^{\circ}$. In this case, the motion of the reflective point on the planetary surface changes the direction relative to the planetary spin axis in a time-dependent fashion, resulting in the frequency modulation of the observed period.

When the star is located in S1 (and also in S3), the reflective point on the planetary surface moves along the constant longitude, and the planet exhibits a nearly identical illuminated and visible part of its surface after one spin rotation period. This implies that $f_{\text{obs}} \approx f_{\text{spin}}$. In contrast, when the star is located at S2 (S4), the reflective point after one spin rotation period moves slightly westward (eastward), leading to $f_{\text{obs}} \approx f_{\text{spin}} + f_{\text{orb}}$ ($f_{\text{obs}} \approx f_{\text{spin}} - f_{\text{orb}}$).

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Same as Figure 1 but for $\zeta = 90^{\circ}$ and $i = 0^{\circ}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Schematic configuration of the system in a *geocentric* frame. The directions of the observer and the planetary spin vector do not vary in time, while the direction of the star is time dependent.</figcaption></figure>

While the above frequency modulation is basically determined by the geometrical configuration of the system characterized by $\zeta$, $i$, and $\Theta_{\text{eq}}$ as mentioned above (see also Figure 3), the most important uncertain factor in modeling the lightcurve is the time-dependent cloud pattern. A planet completely covered by thick, homogeneous clouds, for instance, does not exhibit any photometric variation, and thus one cannot probe the surface information at all. In the case of our Earth, approximately 50%–60% of the surface is covered by clouds on average. Thus, it is not clear to what extent the interpretation from the frequency modulation of the lightcurve is affected or even biased by the properties and time-dependent distribution pattern of clouds.

Since the planetary obliquity is supposed to sensitively change the cloud pattern among others, a feasibility study of
the obliquity measurement from the frequency modulation requires a self-consistent modeling of clouds over the entire surface of a planet. This is why we perform the general circulation model (GCM)<fnref n="11" /> simulation and analyze the simulated lightcurves for different planetary obliquities.

The rest of the paper is organized as follows. Section 2 describes the basic model of the frequency modulation in the lightcurve, the GCM simulation of the Earth with different obliquities, and radiation transfer to simulate lightcurves. Section 3 shows the analysis method of the frequency modulation and the result of the frequency modulation signal extracted from simulated lightcurves. Finally, Section 4 is devoted to the summary and conclusion of the present paper.

## 2. Computational Methods

### 2.1. Basic Strategy to Estimate the Planetary Obliquity from Photometric Variation

For simplicity, we consider a star-planet system in a circular orbit, which is schematically illustrated in Figure 3. In order to compute the photometric variation of the planet, it is convenient to define a geocentric frame in which the planet is located at the origin. The stellar orbit defines the xy-plane, and the star orbits around the z-axis in a counterclockwise manner. The unit vector of the planetary spin is on the yz-plane and expressed as (0, sin ζ, cos ζ) in terms of the planetary obliquity ζ. Thus the direction of the x-axis corresponds to that of the vernal equinox.

The unit vector toward a distant observer is given by (cos Θeq sin i, −sin Θeq sin i, cos i), where i is the inclination, and Θeq is the phase angle measured clockwise from the x-axis (i.e., the vernal equinox).

In this frame, the location of the star on the orbit is specified by its phase angle Θ(t) measured from the observer's projected direction. Because we consider a circular orbit below, Θ(t) = 2πf<sub>orb</sub>t (mod 2π).

K16 computed the frequency modulation based on a maximum-weighted longitude approximation and derived the following formula for f<sub>obs</sub> in the case of a circular orbit:

$$f_{\text{model}} = f_{\text{spin}} + \epsilon_{\zeta}(\Theta)f_{\text{orb}}, \quad (1)$$

where $\epsilon_{\zeta}(\Theta)$ is the modulation factor,<fnref n="12" />

$$\epsilon_{\zeta}(\Theta) = \frac{-\cos \zeta [1 + \cos \Theta \sin i] + \sin \zeta \cos i \sin(\Theta - \Theta_{\text{eq}})}{[\cos(\Theta - \Theta_{\text{eq}}) + \sin i \cos \Theta_{\text{eq}}]^2 + [\cos \zeta \sin(\Theta - \Theta_{\text{eq}}) - \cos \zeta \sin i \sin \Theta_{\text{eq}} - \sin \zeta \cos i]^2} \quad (2)$$

We apply the maximum-weighted longitude approximation and derive a general formula for noncircular orbits in the Appendix. In the present analysis, however, we focus on a circular orbit and adopt Equation (2) for the frequency modulation template.

Following K16, we use the pseudo-Wigner distribution to estimate the frequency modulation of the photometric variation of a given lightcurve. The pseudo-Wigner distribution is the

Fourier transform of the autocorrelation of the data, emphasizing the periodicity near the time of interest and reducing the cross terms and noises. Further detail will be described in Section 3.

### 2.2. GCM Simulation of the Earth with Different Obliquities

We would like to emphasize that the main purpose of the present paper is to examine the feasibility of the planetary obliquity measurement through the frequency modulation of the lightcurve. The cloud covering pattern and fraction are important factors that would degrade the measurement. On the other hand, the precise modeling of the climate is not supposed to be essential for the feasibility. Therefore, various assumptions and limitations of our current GCM simulation described below need to be clarified and understood, but they do not change the main conclusion of the present paper.

We use the GCM code DCPAM5 (the Dennou Club Planetary Atmospheric Model), which has been developed by GFD-Dennou Club<fnref n="13" /> for planetary climate modeling. DCPAM5 has been developed with the aim of being able to calculate an atmospheric condition of various terrestrial planets, using general formulae as much as possible, by excluding properties and modules specific to the Earth (e.g., Noda et al. 2017). DCPAM5 employs the primitive equation system assuming that the vertical component of the equation of motion is hydrostatic.

#### 2.2.1. Setup and Subgrid Physical Processes

We set the computational grids of 32 × 64 × 26 corresponding to latitudinal, longitudinal, and vertical directions, respectively. We carry out calculations in the region up to about 6 mbar, which includes the whole troposphere and a part of the stratosphere. The vertical extent of the model domain is enough for our study to express the generation and motion of clouds because clouds are generated and advected in the troposphere. Our simulation resolves the typical Hadley cell with ~5 × 10 grids, and thus reproduces the global meridional circulation observed on Earth reasonably well.

We use some parameterized physical processes. In the shortwave (visible and near-infrared, corresponding to the range of incident stellar flux) radiation process, we take account of absorption by H₂O and CO₂, absorption and scattering by clouds, and the Rayleigh scattering. In the longwave (mid- and

far-infrared, corresponding to the range of planetary thermal emission) radiation process, we take account of absorption by H₂O and CO₂ molecules and clouds. The level-2.5 closure scheme of Mellor & Yamada (1982) is used for turbulent diffusion. The methods of Beljaars & Holtslag (1991) and Beljaars (1995) are used for surface flux calculation. Moist convection is parameterized by the relaxed Arakawa-Schubert scheme described in Moorthi & Suarez (1992). Large-scale condensation (nonconvective condensation) is parameterized by the scheme of Le Treut & Li (1991). The amount of cloud water is calculated by integrating a time-dependent equation

<footnote n="11">¹¹ “General climate model” is also referred to as GCM. The two terms are often used interchangeably, but sometimes “general circulation model” more specifically implies a part of modules in the “general climate model.” In this sense, our model may be referred to as “an atmospheric general circulation model,” but we do not distinguish between them in the present paper.</footnote>

<footnote n="12">¹² Equation (2) is the correct version of Equation (13) in K16, which contains a couple of typos in signs.</footnote>

<footnote n="13">¹³ http://www.gfd-dennou.org/, and http://www.gfd-dennou.org/library/dcpam/.</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 4. Annual mean cloud column density (g m⁻²) of GCM experiments with six individual obliquities ζ.</figcaption></figure>

including condensation, evaporation, advection, turbulent diffusion, and sedimentation of cloud water. The extinction rate of cloud water is assumed to be proportional to the amount of cloud water, and extinction time is given as an external parameter. The bucket model of Manabe (1969) is used for soil moisture calculation. We use a slab ocean model and its depth to 60 m, the value of Rose (2015).

Our simulation is intended to produce a simulated lightcurve for an Earth twin but with different obliquity ζ. Thus we basically adopt the known parameters of the Earth, except for its obliquity. For simplicity, we set the orbital eccentricity and the orbital period to be $e = 0$ and $P_{orb} = 365.0$ day.

We solve surface temperature and sea ice concentration directly from our simulation, instead of adopting the observed value for Earth with $ζ = 23°.44$, since those values change with the different values of ζ. We use observational data of surface geological properties, neglecting that the change of climate also affects those parameters. Surface albedo is calculated at each grid point according to the surface geological properties, land moisture, and temperature.<fnref n="14" /> Because our GCM does not include the microphysics of cloud formation, cloud parameters are fixed to those for the Earth: the effective radii of water and ice cloud particles are set to be 10 µm and 50 µm, respectively. The lifetimes of water and ice clouds are chosen to be 3240 s and 8400 s, respectively.

### 2.2.2. Initial Conditions

In the present paper, we consider six different values for the obliquity: $ζ = 0°$, $30°$, $60°$, $90°$, $150°$, and $180°$. The simulation runs for $ζ = 0°$, $150°$, and $180°$ start from the isothermal atmosphere of temperature $T_{init} = 280$ K and surface pressure $p_s = 10^5$ Pa with initially vanishing specific humidity and wind speed. Then we evolve those three models for 20 yr so that they reach equilibrium. We call this process the relaxation run.

For $ζ = 30°$, $60°$, and $90°$, we first run the case of $ζ = 15°$ with exactly the same initial conditions mentioned above for 20 yr as the relaxation run. We adopt the final result of the relaxation run with $ζ$ as the initial condition for the next model with $ζ + Δζ$. We set $Δζ = 5°$. The system with $ζ + Δζ$ becomes almost in equilibrium 10 yr after the final epoch of the relaxation run for $ζ$; the annual mean of the total atmospheric energy is constant within the level of 0.1% for the last 5 yr. Thus we stop the relaxation run after 10 yr. We repeat the process up to $ζ = 90°$. The incremental procedure is mainly to save the computation time. In the retrograde runs for $ζ = 150°$ and $180°$, we skip the incremental procedure and made sure that the results reach the equilibrium state after the 20 yr relaxation run directly from the isothermal atmosphere.

The simulated data that we analyze below are computed for an additional one year after each relaxation run. We output the physical parameters every 3 hr for the entire period, which is the required time resolution for detecting the rotation frequency (corresponding to 24 hr) and its modulation from the simulated lightcurve.

### 2.2.3. Climate of Earths with Different Obliquities

Figure 4 shows the annual mean cloud column density distribution of planets with different obliquities. The results for $ζ ≤ 30°$ show cloud belts on the equator and midlatitudes. The clouds around the equator are generated by the Hadley circulation. This circulation also produces subtropical highs, which are shown as the partially cloudless continents around the latitude $λ = 20°–30°$.

The cloud patterns for $ζ = 150°$ and $180°$ are very similar to those for $ζ = 30°$ and $0°$, respectively. This is due to the symmetry with respect to the stellar location for the cases of $ζ$ and $180° - ζ$. The results for $ζ = 60°$ and $90°$ have different cloud patterns because of their atmospheric circulation from the day-side pole to the equator. The present result is roughly consistent with that shown in Williams & Pollard (2003), but a quantitative comparison is beyond the scope of this paper. As we mentioned earlier, however, the precise modeling of the

<footnote n="14">¹⁴ The model codes and related data for the GCM experiments are available at http://www.gfd-dennou.org/library/dcpam/sample/.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 5. Schematic configuration of the scattering geometry.</figcaption></figure>

climate is not the focus of this work. We plan to make a further comparison elsewhere.

## 2.3. Simulated Lightcurves

### 2.3.1. Scattering Model and Radiative Transfer through the Planetary Atmosphere

The total flux of the scattered light from the planet $F(\lambda)$ at wavelength $\lambda$ is computed by integrating the intensities $I$ over the illuminated (I) and visible (V) region of the planetary surface:

$$F(\lambda) = \int_{I \cap V} I(\vartheta_0, \vartheta_1, \varphi; \lambda) \cos \vartheta_1 dS \frac{1}{D_{\text{obs}}^2}, \quad (3)$$

where $\cos \vartheta_1 dS$ is the projected area element of the planetary surface viewed by the observer located at a distance of $D_{\text{obs}}$.

The location of each planetary surface area element is specified by the three angles ($\vartheta_0$, $\vartheta_1$, and $\varphi$), as illustrated in Figure 5. Then the intensity $I$ from the planetary surface area element is given by

$$I(\vartheta_0, \vartheta_1, \varphi; \lambda) = F_{*,p}(\lambda) \cos \vartheta_0 f(\vartheta_0, \vartheta_1, \varphi; \lambda), \quad (4)$$

where $F_{*,p}$ is the incident flux and $f$ is the BRDF (bidirectional reflectance distribution function) that characterizes the scattering properties of the planetary surface.

Because $f$ includes the entire radiative effects of atmosphere, clouds, and solid/liquid planetary surface, we need to perform a numerical radiative transfer calculation through the planetary atmosphere. For that purpose, we compute $f$ using the public code `libRadtran` (Mayer & Kylling 2005; Emde et al. 2016), which solved the radiative transfer based on various detailed models of optical properties of Earth's atmosphere, clouds, aerosols, lands, and ocean.<fnref n="15" />

The `libRadtran` code provides several different options for specific models. We choose the following options:

1. We choose REPTRAN (Gasteiger et al. 2014) for optical properties of the planetary atmosphere.

<table><caption>Table 1
Photometric Bands of Our Mock Observation</caption><thead><tr><td>Band Number</td><td>MODIS Central Wavelength</td></tr></thead><tbody><tr><td>1</td><td>0.469 µm</td></tr><tr><td>2</td><td>0.555 µm</td></tr><tr><td>3</td><td>0.645 µm</td></tr><tr><td>4</td><td>0.858 µm</td></tr><tr><td>5</td><td>1.240 µm</td></tr><tr><td>6</td><td>1.640 µm</td></tr></tbody></table>

2. We compute optical properties of clouds according to Hu & Stamnes (1993). We adopt 10 µm for the effective radius of water cloud particles, as assumed in our GCM simulation.

3. We select the Ross–Li BRDF model (Wanner et al. 1995) for land scattering. We adopt three Ross–Li parameters that are required in `libRadtran` from a remote sensing project of Earth called MODerate resolution Imaging Spectroradiometer (MODIS; Salomonson et al. 1989). More specifically, we choose their data set “snow-free gap-filled MODIS BRDF model parameters.” In doing so, we employ the data in March, neglecting the annual variation. Also, we sample the three parameters at the center of each grid on the planetary surface (32 × 64), instead of averaging over the entire grid. We adopt the above approximation just for simplicity.

4. Since the above particular data set does not have sufficient information for Antarctica, we assume the Lambert scattering and employ the ice albedos of (0.948, 0.921, 0.891, 0.837, 0.562, 0.233), corresponding to the six MODIS bands from 1 to 6 described below. These values of ice albedo are picked from the data “snow-free gap-filled MODIS BRDF model parameters” at (N69°20, W39°35). This approximation is not serious because the ice albedos do not change so much depending on the area.

5. We select the ocean reflection BRDF model of Nakajima (1983) that is implemented in `libRadtran`. We choose 4 m s⁻¹ for the wind speed at 10 m above the ocean. Further detail can be found in Fujii et al. (2010).

6. Finally, we solve the radiative transfer equation through the atmosphere under a plane-parallel approximation. We choose the DIScrete-Ordinate-method Radiative Transfer model (Stamnes et al. 1988).

We use the GCM outputs of water cloud density, ice cloud density, temperature, air density, and vapor mixing ratio as the input vertical profiles of atmosphere and clouds for `libRadtran`. While our GCM simulations distinguish between ice cloud and water cloud, we regard the ice cloud as a water cloud in `libRadtran` so as to reduce the computational cost. For simplicity, we ignore the radiative transfer outside the region of GCM simulation ($z \sim 0-30$ km), including effects due to the upper atmosphere of the planet, exozodiacal dust, and the interstellar medium.

We compute the intensity in six photometric bands centered at the wavelengths of the MODIS bands (Table 1) but with an expanded bandwidth of $\Delta\lambda = 0.1$ µm.

The MODIS project selected their photometric bands so as to characterize the reflection properties of Earth's surface by remote sensing. Figure 6(a) shows examples of an effective albedo (reflectance) spectrum for different components of the Earth's surface: soil, vegetation, and ocean. Three bands (1–3) roughly correspond to the visible colors blue, green, and red,

<footnote n="15">¹⁵ We use the `libRadtran` version 2.0.1. URL: http://www.libradtran.org/doku.php.</footnote>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** (a) Reflectance of the Earth's surface components through the Earth's atmosphere. (b) Transmittance of the Earth's atmosphere and positions of the MODIS photometric bands. Effective albedo and transmittance of the atmosphere of the Earth. Numbers from 1 to 6 above the panels indicate numbers of bands in the MODIS project (see Table 1). For an atmospheric profile, the US Standard Atmosphere (Anderson 1986) without clouds is used.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Effective albedo of clouds calculated with libRadtran from our GCM outputs with the obliquity $\zeta = 30^{\circ}$.</figcaption></figure>

respectively. Figure 6(a) exhibits a clear difference among the three components, ocean, soil, and vegetation. Incidentally, the MODIS project chooses three near-IR bands that correspond to observational windows of Earth's atmosphere (Figure 6(b)).

As we have already emphasized, the cloud distribution is the most important ingredient in our mock simulation. In order to examine the dependence on their properties, we generate a simple cloud distribution as follows. Our GCM result for $\zeta = 30^{\circ}$ indicates that the simulated cloud distribution has typical column densities of $0.040_{-0.025}^{+0.050} \text{ kg m}^{-2}$. Thus we redistribute all of the clouds homogeneously within 0.0–0.3, 0.5–1.0, and 3.0–8.0 km, which roughly correspond to the typical heights for mist, lower clouds, and middle clouds for Earth.

Figure 7 shows the resulting effective albedos for those mock clouds, indicating that the albedos are mainly determined by the column density and are fairly insensitive to the height of the clouds.

### 2.3.2. Simulated Images and Lightcurves of an Earth Twin

Before performing the frequency modulation analysis, let us present examples of the apparent images and lightcurves from our mock observation.

Figures 8(a) and (b) show the images of an Earth twin in January and July, respectively, with different obliquities viewed by a distant observer at $i = 0^{\circ}$. Plotted from left to right are input surface distribution, illuminated and visible part of the cloudless Earth with atmosphere, illuminated and visible part of the Earth with both cloud and atmosphere from our GCM simulation, and the corresponding cloud distribution. The arrows indicate the incident direction of the starlight.

The input surface distribution (the left images) is computed from the intensity of land alone, neglecting the contribution of the ocean reflection. The land is assumed to be covered by the US Standard Atmosphere (Anderson 1986), and land scattering
<figure><img src="image_13.png" /><figcaption>**Figure 8.** (a) Images of an Earth twin from our GCM simulations with different obliquities viewed by a distant observer at *i* = 0° in January. From left to right, we plot the input surface images, illuminated and visible part of the cloudless Earth with atmosphere, illuminated and visible part of the Earth with both cloud and atmosphere from our GCM simulation, and the corresponding cloud distribution. The orange arrows show the direction of stellar illumination. We adopt the RGB flux ratio to be the intensity ratio of bands 3:2:1 (0.645 µm:0.555 µm:0.469 µm) and apply the gamma correction with *γ* = 1/2.2 so as to roughly represent the apparent colors. (b) Same as (a) but in July.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 8.** (Continued.)</figcaption></figure>

is approximated by a Lambertian distribution. Since those images are just for reference, we assume the geometric configuration with ($\vartheta_0, \vartheta_1$) = (0°, 0°).

The different surface components are illustrated in orange, green, and blue for continents, vegetation, and oceans,

respectively. In the left images, one may identify North and South America, Eurasia, Africa, and Antarctica. The images of the cloudless Earth exhibit relatively well the colors of the surface below the atmosphere, and they also show an oceanic glint (oceanic mirror reflection) in the illuminating direction
band 1 (0.469 µm), Δλ = 0.1 µm, January

<figure><img src="image_15.png" /><figcaption>**Figure 9.** Examples of simulated lightcurves of an Earth twin with different obliquities for an observer at i = 0°. The photon counts N1,0 and N4,0 correspond to the set of parameters Dobs = 10 pc, Dtel = 4 m, texp = 3 hr, and Δλ = 0.1 µm and are scaled as Equation (5). The quoted error bars consider the photon shot noise alone.</figcaption></figure>

band 4 (0.858 µm), Δλ = 0.1 µm, January

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

(Sagan et al. 1993; Robinson et al. 2010). Those signatures of the surface components are significantly degraded by the clouds, but one may still identify the presence of the Sahara desert for ζ < 60° in Figure 8(a), for instance. Although one may not identify the Sahara desert for ζ = 90° in January (Figure 8(a)), the Sahara desert appears in the visible and illuminating part in July (Figure 8(b)). Thus it can be still used as a frequency modulation indicator for part of a year.

Figures 8(a) and (b) reconfirm that the cloud distribution weakens the surface information in photometric monitoring, but it still indicates that the diurnal variation and possibly its frequency modulation detection are feasible if there exists a good tracer of the global planetary surface like the Sahara desert.

Mock photometric monitoring of the images presented in Figures 8(a) and (b) generates the corresponding simulated lightcurves. Throughout the analysis in what follows, we consider an observer located at i = 0° for simplicity. Since we output results of our GCM simulations every three hours, we construct simulated lightcurves from those discrete snapshots. Then we ignore the change of the lightcurve during the three hours, and we construct mock lightcurves sampled every three hours. While this approximate method significantly affects the lightcurve variation on a timescale less than three

hours, the variation around a planetary spin period (24 hr) of interest to us is hardly affected.

Figure 9 shows an example of one-week lightcurves in January for Earth twins with different obliquities; the left and right panels correspond to those in bands 1 and 4, respectively. We assume that the star–planet system is located at a distance Dobs away from the telescope of diameter Dtel and exposure time texp. In an idealized case where both the light from the host star and other instrumental noises are completely neglected, the photon counts at band i with bandwidth Δλ are scaled as

$$N_i(t) = N_{i,0}(t) \left( \frac{D_{\text{obs}}}{10 \text{ pc}} \right)^{-2} \left( \frac{D_{\text{tel}}}{4 \text{ m}} \right)^2 \times \left( \frac{t_{\text{exp}}}{3 \text{ hr}} \right) \left( \frac{\Delta \lambda}{0.1 \ \mu\text{m}} \right). \quad (5)$$

The photon counts in Figure 9 correspond to N1,0 (left panel) and N4,0 (right panel) for bands 1 and 4 in Equation (5). In practice, we compute Nᵢ(t) from snapshots every three hours, assuming texp = 3 hr.

The simulated lightcurves for ζ ≤ 60° exhibit a kind of diurnal periodicity, which does not reflect the surface information directly, but comes mainly from the cloud pattern
correlated with the surface distribution. As $\zeta$ increases ($\zeta \ge 90^{\circ}$), the diurnal periodicity is not easy to identify. As we mentioned above, the Sahara desert played an important role as a tracer of the planetary rotation, and the annual-average cloud pattern is also correlated to the distribution of the surface components. This is why the diurnal periodicity is more visible for photometric monitoring of the Northern Hemisphere in the case of the Earth. Although in the case of $\zeta = 60^{\circ}$, the northern part of South America takes the role as well, it eventually moves out of the visible and illuminated region as $\zeta$ increases.

### 3. Time-Frequency Analysis of Simulated Lightcurves and Parameter Estimation

Given the simulated lightcurves, we perform the frequency modulation analysis following K16. In practice, we use a numerical code juwvid to compute the pseudo-Wigner distribution, which is publicly available from the website.<fnref n="16" />

Our time-frequency analysis proceeds as follows. First we compute $N_i(t)$ from our simulated lightcurves every three hours over an orbital period of one year. Then we sample $N_{i,obs}(t)$ from the Poisson distribution with the expectation value of $N_i(t)$. In other words, we consider the shot noise alone in the analysis below. In total, we have $N_{data} = 2920$ (= 1 yr/3 hr) data points, and we duplicate the data points with the period of 1 yr.

We divide each lightcurve into 73 segments consisting of 40 consecutive data points (i.e., 3 hr × 40 ≈ 5 days). Then we compute the mean $\mu$ and standard deviation $\sigma$ of $N_i(t)$ in each segment and convert to the normalized lightcurve $s(t) \equiv (N_i(t) - \mu)/\sigma$. Finally we compute the pseudo-Wigner distribution:

$$g(f, t) = \int_{-\infty}^{\infty} H(\tau)z(t + \tau/2)z^*(t - \tau/2)e^{-2\pi i \tau} d\tau, \quad (6)$$

where

$$z(t) = \frac{1}{\pi} \int_{0}^{\infty} \tilde{s}(w)e^{iwt}dw \quad (7)$$

is the analytic signal of $s(t)$, with $\tilde{s}(w)$ being the Fourier transform of the normalized lightcurve $s(t)$ in the present case. We choose the window function $H(\tau)$ as the following Hamming window function:

$$h(\tau; T_w) = \begin{cases} 0.54 + 0.46 \cos(2\pi\tau/T_w) & \text{for } |\tau| \le T_w/2 \\ 0 & \text{otherwise.} \end{cases} \quad (8)$$

In practice, we adopt $T_w = 0.25$ yr for the window width of the Hamming window function.

The pseudo-Wigner distribution is an appropriate time-frequency distribution for extracting the instantaneous frequency (e.g., Cohen 1995), as explained below. Let us consider a single mode signal $z = A(t) e^{i\psi(t)}$ with an instantaneous phase $\psi(t)$, where $A(t) \in \mathbb{R}$ is the amplitude of the mode. The ideal time-frequency representation is a delta function $\rho(f, t) = A(t)^2 \delta_D(f - f_{ins}(t))$, where $f_{ins}(t)$ is the instantaneous frequency defined by

$$f_{\text{ins}}(t) \equiv \frac{1}{2\pi} \frac{d\psi(t)}{dt}. \quad (9)$$

Then, the inverse Fourier transform of $\rho$ can be written as

$$\begin{align} \hat{\rho}(\tau, t) &= A(t)^2 e^{2\pi i f_{\text{ins}}(t)\tau} = A(t)^2 e^{i\tau\psi'(t)} \nonumber \\ &\approx A(t)^2 e^{i\psi(t+\tau/2)-i\psi(t-\tau/2)} \nonumber \\ &= z(t + \tau/2)z^*(t - \tau/2), \tag{10} \end{align}$$

where we use the linear approximation $\psi'(t) \approx [\psi(t + \tau/2) - \psi(t - \tau/2)]/\tau$ in the last two terms.

Performing the Fourier transform of Equation (10) with the time window, we obtain the pseudo-Wigner distribution. Because the linear approximation is valid only for the linear frequency modulation such as $f_{ins}(t) \approx at + b$ ($a, b$ are constant values), the width of the window should be chosen to be comparable to the scale of the nonlinear feature of the frequency modulation. The derivative of Equation (10) by $\tau$ at $\tau = 0$ provides

$$\begin{align} \frac{d}{d\tau}[z(t + \tau/2)z^*(t - \tau/2)]|_{\tau=0} &= iA(t)^2\psi'(t) \nonumber \\ &= 2\pi i \int_{-\infty}^{\infty} f\rho(f, t) df. \tag{11} \end{align}$$

Also, the mode amplitude is rewritten as

$$|z(t)|^2 = A(t)^2 = \int_{-\infty}^{\infty} \rho(f, t) df. \quad (12)$$

Then, the instantaneous frequency is formally estimated by the weighted form as

$$f_{\text{ins}}(t) = \frac{1}{2\pi} \frac{d\psi(t)}{dt} = \frac{\int_{-\infty}^{\infty} f\rho(f, t) df}{\int_{-\infty}^{\infty} \rho(f, t) df}. \quad (13)$$

In practice, one can estimate the peak value of the pseudo-Wigner distribution as an instantaneous frequency to avoid the effect of noise. In this expression, we need a complex-valued signal with a nonnegative frequency component of the signal. That is why we convert a real-valued signal $s(t)$ to the analytic signal $z(t)$ in Equation (7).

We calculate the pseudo-Wigner distribution $g(f, t)$ over the range $f_{min} < f < f_{max}$ using Equation (6). Specifically we choose $f_{min} = 0.98$ (day$^{-1}$) and $f_{max} = 1.02$ (day$^{-1}$) throughout the analysis. Since our lightcurves are sampled every 3 hr, the corresponding frequency resolution is not good enough to determine the value of $f_{spin}$ precisely. Therefore we adopt a nonuniform fast Fourier transform scheme (Greengard & Lee 2004) following K16, and we achieved the frequency resolution of $\delta f = (f_{max} - f_{min})/N_f$ after applying an appropriate smoothing of the lightcurves. We choose $N_f = 1024$ in what follows, and the resulting resolution $\delta f \approx 4 \times 10^{-5}$ (day$^{-1}$) is better than the modulation amplitude detected in Figure 16 by a factor of 100.

#### 3.1. Single-band Analysis

Consider first the frequency modulation for single-band lightcurves. Figure 10 is similar to Figure 9 but plots simulated noiseless (without shot noise) lightcurves in the photometric bands 1–6 for ($\zeta, i$) = (0°, 0°).

<footnote n="16">¹⁶ https://github.com/HajimeKawahara/juwvid</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 10.** Same as Figure 9 for $\zeta = 0^{\circ}$, but in photometric bands 1–6 without photometric noise. Solid and dotted symbols indicate the lightcurves with and without clouds, respectively.</figcaption></figure>

As clearly indicated by Figure 10, the apparent diurnal variation in each band originates from the cloud pattern that is correlated with the land-ocean distribution. These surface-correlated clouds were also found in Earth observations by the Deep Space Climate Observatory as the second component of the principal component analysis (Fan et al. 2019). While our analysis did not directly identify the component, it appears to be imprinted in the diurnal variation in a single band. The amplitude of the single-band lightcurves is basically determined by the cloud albedo (Figure 7) multiplied by the incident solar flux. This is why the amplitude of the diurnal variation with clouds in Figure 10 is relatively large around the visible wavelengths (bands 2 and 3) and declines sharply in the near-infrared (bands 4–6).

We note that Figure 10 also indicates the anticorrelation of the lightcurve modulation between cloudless and cloudy cases. For a cloudless case, the photometric variation is mainly due to the land component that has larger albedos (Figures 6(a) and 7). Since clouds are much brighter, however, the photometric variation of a cloudy case is sensitive to the location of clouds, which tend to avoid the continent, in particular desert regions, and rather form preferentially above the ocean. Thus the locations of lands and clouds are anticorrelated, leading to the anticorrelation illustrated in Figure 10. This also explains that the periodic signature of the lightcurve for a cloudy case is weaker for redder bands, because lands become brighter in redder bands and compensate for the variation due to clouds.

The corresponding color map for the pseudo-Wigner distribution on the time-frequency plane (Figure 11) clearly illustrates the above trend that redder bands have weaker

signals. The color indicates the absolute value of the time-frequency distribution density $g(f, t)$, whose maximum value is normalized as unity. Since Figure 10 is for $\zeta = 0^{\circ}$, the period for the apparent diurnal variation should be constant and does not show any frequency modulation. The tiny frequency modulation of $\sim 0.001 \text{ day}^{-1}$ visible in Figure 11 is simply due to the time-dependent inhomogeneous distribution of clouds.

Consider next the time-frequency representation of the band 1 lightcurves for different obliquities (Figure 12). We adopt band 1 because it produces the clearest ridge on the time-frequency representation in Figure 11. The dashed lines show the model frequency modulation $f_{\text{model}}(t)$, Equation (1). The signature of the frequency modulation from the single-band lightcurves is not strong and is barely identifiable only for $\zeta \le 30^{\circ}$. Though the amplitude of frequency modulation is zero for $\zeta = 0^{\circ}$, the signature of the *constant* apparent frequency is clearly visible. This obliquity dependence reflects the specific distribution pattern of land and ocean on the Earth. As shown in the $\zeta = 150^{\circ}$ image in Figure 8(a), the illuminated and visible part in winter is dominated by Antarctica, and there is no significant diurnal variation in the lightcurve. On the other hand, in summer, Antarctica is almost invisible, and parts of Africa and South America generate the diurnal variation signal instead.

## 3.2. Multiband Analysis

As shown in Section 3.1, the single-band analysis does not properly extract the information of the correct frequency modulation, due to the anticorrelation between lands and
<figure><img src="image_18.png" /><figcaption>Figure 11. Time-frequency representation corresponding to the noiseless lightcurves of Figure 10.</figcaption></figure>

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

<figure><img src="image_20.png" /><figcaption>Figure 12. Time-frequency representation of band 1 noiseless signal for different obliquities ($\zeta$ = 0°, 30°, 60°, 90°, 150°, and 180°), corresponding to the results with clouds (upper panel) and without clouds (lower panel) plotted in Figure 11. Thick dashed lines are the prediction based on the maximum-weighted longitude approximation model, $f_{model}(t)$; see Equations (1) and (2) in the main text. Thin blue points indicate $f_{data,max}(t)$, the frequency corresponding to the maximum value of g ($f, t$) over a range $f_{min} < f < f_{max}$ at each epoch; see Equation (16).</figcaption></figure>

clouds. In order to detect the diurnal period that is due to the planetary surface distribution, we need to remove the time-dependent cloud pattern as much as possible.

As inferred from the wavelength dependence of albedos for land and clouds, bands 1 and 4 are mainly sensitive to clouds and clouds+land, respectively (see Figures 6(a) and 7). Thus

the difference of the photon counts $N_1(t)$ and $N_4(t)$ roughly removes the contribution from clouds.

For definiteness, we choose the following linear combination of bands 1 and 4:

$$C_{1-4}(t) = N_1(t) - \alpha_{1-4}N_4(t), \quad (14)$$
<figure><img src="image_21.png" /><figcaption>Figure 13. Simulated noiseless C1-4 lightcurves for different obliquities (ζ = 0°, 30°, 60°, 90°, 150°, and 180°). We adopt the same set of parameters as in Figure 9.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 14. Pseudo-Wigner distribution g(f, t) of the noiseless C1-4. Thick dashed lines indicate fmodel(t), while thin blue points indicate fdata,max(t), the frequency corresponding to the maximum value of g(f, t) over a range fmin < f < fmax at each epoch; see Equation (16). Due to the quality of the data, the values of fdata,max(t) are not robust for ζ = 150° and 180° and are discontinuous.</figcaption></figure>

$$ \alpha_{1-4} = \frac{F_{*1}\lambda_1\Delta\lambda}{F_{*4}\lambda_4\Delta\lambda} \sim 1.12. \qquad (15) $$

The above combination is derived assuming that the albedo of clouds is roughly independent of the wavelength, and thus the contribution of the clouds is canceled, at least partially as shown in Figure 13. While the cloud effect may be removed more efficiently by combining other bands appropriately, it is beyond the scope of the present paper. Thus we perform the frequency modulation analysis using Equation (14) in what follows.

Figure 13 shows an example of simulated noiseless lightcurves using C1-4(t), and Figure 14 is the corresponding

time-frequency representation. A comparison between Figures 12 and 14 clearly indicates that the multiband analysis suppresses the time-dependent cloud effect and significantly improves the frequency modulation signal.

We note that the amplitude of the frequency modulation signature depicted in Figure 14 sensitively depends on the value of ζ, reflecting the specific surface distribution on Earth. As indicated in Figures 8(a) and (b), the Southern Hemisphere, especially around the South Pole, of the Earth is occupied by Antarctica and ocean, in an approximately axisymmetric manner. Thus the diurnal variation of the Southern Hemisphere (for example, viewed from the direction of i = 0° if ζ = 180°) is difficult to detect. This also applies to the ζ = 150° case in
<figure><img src="image_23.png" /><figcaption>**Figure 15.** Pseudo-Wigner distribution for oblique Earth twins from the shot-noise-limited photometric monitor. The top, middle, and bottom panels are for the space telescope apertures $D_{tel} = 2, 4$, and $6$ m, and the left, center, and right panels are for the planetary obliquities $\zeta = 30^{\circ}, 60^{\circ}$, and $90^{\circ}$. We assume that $D_{obs} = 10$ pc, $t_{exp} = 3$ hr, and $\Delta\lambda = 0.1$ µm.</figcaption></figure>

which the frequency modulation signal is clear only in summer, as described at the end of Section 3.1.

In contrast, the Northern Hemisphere is roughly divided into two major distinct components: the Eurasian continent and the Pacific Ocean. This large-scale inhomogeneity, in particular the Sahara desert, acts as a good tracer of an asymmetric surface pattern, yielding a relatively large amplitude signal of the frequency modulation (see Figures 8(a) and (b)). This is why a clear frequency modulation signal in the case of $\zeta \le 90^{\circ}$ can be detected for an observer located at $i = 0^{\circ}$.

### 3.3. Feasibility of the Obliquity Estimate

All of the pseudo-Wigner distributions above (Figures 11, 12, and 14) are based on noiseless data. Now we are in a position to examine to what extent one can estimate the planetary obliquity from the long-term photometric monitoring via the frequency modulation method. For that purpose, we assume a dedicated space mission with a telescope aperture of $D_{tel} = 2, 4$, and $6$ m. Again we consider idealized cases in which the photometric noise comes from the photon shot noise alone, and we generate a set of $C_{1-4}(t)$ lightcurves from the photon counts $N_1(t)$ and $N_4(t)$ obeying the Poisson statistics. Examples of the resulting frequency modulation are presented in Figure 15.

The model frequency modulation is determined by the five parameters ($\zeta, f_{spin}, \Theta_{eq}, i, f_{orb}$) that are listed in Table 2: the planetary obliquity $\zeta$, the planetary spin frequency $f_{spin}$, the angle of the vernal equinox measured from the location of the observer projected on the orbital plane $\Theta_{eq}$, the observer's inclination $i$, and the orbital frequency of the planet $f_{orb}$. Among them, $i$ and $\Theta_{eq}$ simply specify the location of the observer relative to the system and are not so interesting. The remaining three parameters, $\zeta, f_{spin}$, and $f_{orb}$, are important because they characterize the star-planet system.

<table><caption>**Table 2**
Initial and Fixed Parameter Sets for Best-fit Search</caption><thead><tr><th>Initial Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Obliquity ζ</td><td>15<i>i</i>° (<i>i</i> = 1, 2, ..., 12)</td></tr><tr><td>Θ<sub>eq</sub></td><td>60<i>j</i>° (<i>j</i> = 1, 2, ..., 5)</td></tr><tr><td>Spin frequency f<sub>spin</sub></td><td>366 (yr<sup>-1</sup>)</td></tr><tr><td>Fixed parameter</td><td>value</td></tr><tr><td>Orbital inclination <i>i</i></td><td>0°</td></tr><tr><td>Orbital frequency f<sub>orb</sub></td><td>1 (yr<sup>-1</sup>)</td></tr></tbody></table>

In order to estimate $\zeta$, which cannot be estimated otherwise and thus is of primary interest to us, we need to perform eventually a joint analysis of the five parameters in a Bayesian fashion. In the present study, however, we would like to examine the feasibility of the determination of $\zeta$ and $f_{spin}$, assuming that $i$ and $f_{orb}$ are known, for simplicity. The precise spectroscopic and astrometric data would determine $i$ and $f_{orb}$. Also, $f_{spin}$ may be estimated from the photometric data on relatively short timescales apart from the uncertainty of $\epsilon(\Theta)f_{orb}$ in Equation (1).

Under a similar assumption, K16 attempted to find the best-fit values for $\zeta, f_{spin}$, and $\Theta_{eq}$ by minimizing

$$R_1(\Theta_{eq}, \zeta, f_{spin}) = \sum_{j=1}^{N_{data}} |f_{data,max}(t_j) - f_{model}(t_j; \Theta_{eq}, \zeta, f_{spin})|^2, \quad (16)$$

where $f_{model}(t)$ is the frequency derived from the maximum-weighted longitude approximation, Equation (1), and $f_{data,max}(t)$ corresponds to the maximum value of $g(f, t)$ over a range $f_{min} < f < f_{max}$ at each epoch $t$. We tried the same fitting, but
<figure><img src="image_24.png" /><figcaption>Figure 16. Spin parameters ($\zeta, f_{spin}$) estimated from the normalized Gaussian model $\tilde{g}_{model}$. Black crosses show the input spin values ($\zeta, f_{spin}$) = (30°, 366 yr$^{-1}$) and (60°, 366 yr$^{-1}$) for the left and right panels, respectively. We plot the best-fit values from 1000 realizations of shot-noise-limited observation. The shot noise is assumed from the telescope diameter of $D_{tel}$ = 4 m and 6 m for the top and bottom panels, respectively.</figcaption></figure>

the result is not robust against the shot noise, especially when the frequency modulation signal is weak.

Therefore, we empirically improve the fit by taking account of the distribution around the $f_{data,max}(t)$ as well. More specifically, we construct a Gaussian weighted model $\tilde{g}_{model}(f, t)$ for the time-frequency distribution:

$$ \tilde{g}_{model}(f, t) = \exp\left[-\frac{(f - f_{model}(t; \Theta_{eq}, \zeta, f_{spin}))^2}{2\sigma_f^2}\right], \quad (17) $$

where $\sigma_f$ is a new fitting parameter that is introduced to account for the finite width of the frequency distribution around $f_{model}$. Then we minimize the following quantity,

$$ R_2(\Theta_{eq}, \zeta, f_{spin}, \sigma_f) = \sum_{i=1}^{N_f} \sum_{j=1}^{N_{data}} \times \left| \frac{g_{data}(f_i, t_j)}{g_{data}(f_{data,max}(t_j), t_j)} - \tilde{g}_{model}(f_i, t_j) \right|^2, \quad (18) $$

to find the best-fit $\zeta, \Theta_{eq}, f_{spin}$, and $\sigma_f$. The value of $\sigma_f$ should be roughly equal to $1/T_w$ because the time-frequency representation of a signal $z(t) = e^{if_{ins}t}$ based on the pseudo-Wigner distribution has a dispersion corresponding to the Fourier transform of the window function $\tilde{h}(f - f_{ins}; T_w)$, and this

dispersion is flattened due to the noise and nonlinear frequency modulation.

In practice, we use the Levenberg-Marquardt algorithm mpfit (Markwardt 2009) to find the best-fit parameters. This algorithm is a practical and fast algorithm of the least squares method for nonlinear functions. We fit the time-frequency distribution for $\zeta = 30°$ and 60°. Table 2 summarizes our initial parameters in addition to the fixed orbital parameters that we assume to be a priori known.

Figure 16 shows the distribution of the best-fit estimates on the $\zeta-f_{spin}$ plane from 1000 different realizations. The black cross symbols indicate the input values, ($\zeta, f_{spin}$) = (30°, 366 yr$^{-1}$) and (60°, 366 yr$^{-1}$), for the left and right panels, respectively. The top and bottom panels show the results based on the shot-noise-limited observations with $D_{tel}$ = 4 m and 6 m, respectively. The numbers in each panel denote the mean and $1\sigma$ estimated from 1000 realizations.

The systematic offsets of $(\Delta\zeta)_{sys} \approx 3°$ and $(\Delta f_{spin})_{sys} \approx 0.03$ yr$^{-1}$ result most likely from the specific pattern of the continents on the Earth. Indeed, the previous simplified analysis by K16 also found a similar level of systematic offset of the planetary obliquity (approximately several degrees; see Figure 8 of K16). K16 added noises empirically into his mock data, neglecting the time-dependent cloud distribution that we compute here.
The fact that the systematic offsets between the two analyses are similar indicates, therefore, that they should be ascribed to the specific surface pattern of the Earth itself. Indeed Eurasia, North Africa, and South America are distributed roughly from northeast to southwest. This latitudinal pattern is consistent with the positive systematic offset of the obliquity exhibited in Figure 16. Since the amplitude of the systematic offset would depend on the specific pattern of the planetary surface to some extent, it is difficult to predict it a priori, but it is important to bear in mind that it could amount to several degrees, much larger than the statistical uncertainty shown in Figure 16.

## 4. Summary and Conclusion

The direct imaging of Earth-like planets is very challenging but will provide ground-breaking data sets for astronomy, planetary science, and biology, if successful eventually. One notable example is the reconstruction of surface components (e.g., Sagan et al. 1993; Ford et al. 2001; Fujii et al. 2010, 2011; Kawahara & Fujii 2011; Fujii & Kawahara 2012; Suto 2019), and it may be even possible to measure the planetary obliquity through the frequency modulation of the photometric lightcurve of future direct-imaged Earth-like planets, as proposed by Kawahara (2016).

We have examined the feasibility of the methodology by creating simulated lightcurves of our Earth–Sun systems but with different planetary obliquities. First, we performed the GCM simulation for those systems with particular emphasis on the time-dependent cloud distribution. Second, we computed the scattered light in six photometric bands by solving the radiation transfer of the incident starlight through the clouds and atmosphere, taking into account the scattering due to the different surface components under the parameterized bidirectional reflectance distribution function models (Nakajima 1983; Wanner et al. 1995). Third, the resulting light from the planet was mock-observed every three hours over the orbital period of one year, and simulated lightcurves were constructed by combining the different photometric bands so as to suppress the effect of the time-dependent cloud pattern. Finally, we computed the frequency modulation of the lightcurves using the pseudo-Wigner distribution and attempted to estimate the planetary obliquities for cases dominated by photon shot noise.

We found that the frequency modulation signal is crucially dependent on the presence of a large-scale inhomogeneity on the planetary surface. Indeed, this is the case for the Northern Hemisphere of our Earth; in particular, the Sahara desert turned out to be a useful tracer of the planetary spin rotation. The Southern Hemisphere, on the other hand, is relatively featureless, and the frequency modulation signal is weak.

As a result, we found that a dedicated 4 m space telescope at 10 pc away from the system in the face-on view relative to the observer can estimate the planetary obliquity within the uncertainty of several degrees in principle (in the shot-noise-limited case). Although this conclusion is based on several idealized assumptions at this point, we believe that it is very encouraging for the future exploration of the direct imaging of Earth-like planets.

We thank an anonymous referee for numerous constructive suggestions and comments that significantly improved the early manuscript of the paper. This work is supported by the Japan Society for the Promotion of Science (JSPS) Core-to-Core Program “International Network of Planetary Sciences,” the

European Research Council under the European Union’s Horizon 2020 research and innovation program (Grant Agreement 679030/WHIPLASH), JSPS KAKENHI Grant Nos. JP18H01247 and JP19H01947 (Y.S.), JP17K14246 and JP18H04577 (H.K.), JP19H01947 (M.I.), and JP17H06457 (K. K. and Y.O.T), and the Astrobiology Center Program of National Institutes of Natural Sciences (grant No. AB311025). Numerical computations and analyses were partly carried out on the PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. Y.N. acknowledges the support by a fellowship from the Advanced Leading Graduate Course for Photon Science at the University of Tokyo.

Software: DCPAM5 (http://www.gfd-dennou.org/library/dcpam), libRadtran (Mayer & Kylling 2005; Emde et al. 2016), REPTRAN (Gasteiger et al. 2014), mpfit (Markwardt 2009), juwvid (https://github.com/HajimeKawahara/juwvid).

## Appendix

### Behavior of the Frequency Modulation Factor $\epsilon_{\zeta}(\Theta)$ for an Eccentric Orbit

The modulation factor, Equation (2), first derived by K16 assumes a circular orbit for simplicity. We compute a generalized expression for an eccentric orbit, and we present the effect of eccentricity on the frequency modulation based on the maximum-weighted longitude approximation.

For an eccentric orbit, it is more convenient to consider a geocentric frame where the x-axis is the direction toward the periapsis as shown in Figure A1, instead of the vernal equinox (see Figure 3). In this frame, the spin vector is no longer on the yz-plane, and we introduce a new parameter $\beta$, which denotes the azimuthal angle of the planetary spin measured from the y-axis. Similarly, the location of the observer is specified by the phase angle from the periapsis $\Theta_{per}$, and $\Theta$ is now the azimuthal angle measured clockwise from the periapsis, that is, the true anomaly. The frame reduces to that shown in Figure 3 for $e \to 0$, $\beta \to 0^{\circ}$, $\Theta_{per} \to \Theta_{eq}$, and $\Theta \to \Theta - \Theta_{eq}$.

Following K16, we assume that the longitude of the reflective point on the planetary surface, $\hat{\phi}_{M}$, traces faithfully the observable periodicity of the planetary scattered light. Then the observed frequency $f_{obs}$ is given in terms of $\hat{\phi}_{M}$ as

$$ 
\begin{aligned}
f_{\text{obs}}(t) &= -\frac{1}{2\pi} \frac{d\hat{\phi}_{M}}{dt} = -\frac{1}{2\pi} \frac{d\Theta}{dt} \frac{\partial\hat{\phi}_{M}}{\partial\Theta} \\
&= f_{\text{spin}} + \frac{1}{2\pi} \frac{d\Theta}{dt} \epsilon_{\zeta}(\Theta),
\end{aligned} 
\quad (A1) $$

<figure><img src="image_25.png" /><figcaption>Figure A1. Schematic configuration of a geocentric frame for eccentric orbits.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure A2.** Samples of frequency modulation applied to eccentric orbits: (a) $\zeta = 30^{\circ}$ and (b) $\zeta = 60^{\circ}$. Each line shows the frequency modulation specific to the spin configuration of the planet. Solid lines are eccentric cases ($e = 0.2$), and dashed lines are circular cases ($e = 0.0$).</figcaption></figure>

where

$$ \epsilon_{\zeta}(\Theta) \equiv -\frac{\partial(\hat{\phi}_{\mathrm{M}} + \Phi)}{\partial\Theta} = -\frac{\kappa'(\Theta)}{1 + \kappa(\Theta)^2}, \quad (\mathrm{A}2) $$

$$ \kappa(\Theta) \equiv \tan(\hat{\phi}_{\mathrm{M}} + \Phi), \quad (\mathrm{A}3) $$

$$ \Phi \equiv 2\pi f_{\mathrm{spin}} t. \quad (\mathrm{A}4) $$

For a circular orbit, $d\Theta/dt$ is equal to $f_{\mathrm{orb}}$. For $e \neq 0$, however, it cannot be written explicitly in terms of $t$, but is expressed as

$$ \frac{d\Theta}{dt} = 2\pi f_{\mathrm{orb}} (1 - e^2)^{-3/2} (1 + e \cos \Theta)^2. \quad (\mathrm{A5}) $$

In the geocentric frame, unit vectors toward the star and the observer, $e_s$ and $e_o$, are given as $e_s = (\cos\Theta, \sin\Theta, 0)$ and $e_o = (\cos\Theta_{\mathrm{per}}\sin i, -\sin\Theta_{\mathrm{per}}\sin i, \cos i)$, respectively. Thus the unit vector toward the reflective point, $e_M$, is

$$ e_M = \frac{e_s + e_o}{|e_s + e_o|} = \frac{1}{L} \begin{pmatrix} \cos \Theta + \cos \Theta_{\mathrm{per}} \sin i \\ \sin \Theta - \sin \Theta_{\mathrm{per}} \sin i \\ \cos i \end{pmatrix}, \quad (\mathrm{A6}) $$

where $L \equiv |e_s + e_o| = \sqrt{2 + 2 \cos(\Theta + \Theta_{\mathrm{per}}) \sin i}$.

Consider a point on the planetary surface specified by the latitude $\lambda$ and longitude $\phi$ in the rest frame of the planet. The

in the geocentric frame as

$$ e_R = R_z(\beta) R_x(-\zeta) \hat{S}(\Phi) e_R', \quad (\mathrm{A7}) $$

where $\hat{S}(\Phi)$ is a spin rotation operator ($\phi \to \phi + \Phi$), and $R_i$ is the rotation matrix counterclockwise around the $i$-axis. Note that $R_z(\beta)$ is required for a noncircular case.

We apply the generic transformation, Equation (A7), to compute the component of the reflective point in the planetary frame. Then we obtain

$$ \begin{aligned} e_M'(\hat{\phi}_M + \Phi) &= \hat{S}(\Phi)e_M'(\hat{\phi}_M) = R_x(\zeta)R_z(-\beta)e_M \\ &= \frac{1}{L} \begin{pmatrix} \cos(\Theta - \beta) + \sin i \cos(\Theta_{\mathrm{per}} + \beta) \\ \cos \zeta \{\sin(\Theta - \beta) - \sin i \sin(\Theta_{\mathrm{per}} + \beta)\} - \sin \zeta \cos i \\ \sin \zeta \{\sin(\Theta - \beta) - \sin i \sin(\Theta_{\mathrm{per}} + \beta)\} + \cos \zeta \cos i \end{pmatrix}. \end{aligned} \quad (\mathrm{A8}) $$

The ratio of the $x$ and $y$ components in Equation (A8) yields

$$ \tan(\hat{\phi}_M + \Phi) = \frac{\cos \zeta \{\sin(\Theta - \beta) - \sin i \sin(\Theta_{\mathrm{per}} + \beta)\} - \sin \zeta \cos i}{\cos(\Theta - \beta) + \sin i \cos(\Theta_{\mathrm{per}} + \beta)}. \quad (\mathrm{A9}) $$

Therefore, Equation (A2) reduces to

$$ \epsilon_{\zeta}(\Theta) = \frac{-\cos \zeta \{1 + \sin i \cos(\Theta + \Theta_{\mathrm{per}})\} + \sin \zeta \cos i \sin(\Theta - \beta)}{[\cos(\Theta - \beta) + \sin i \cos(\Theta_{\mathrm{per}} + \beta)]^2 + [\cos \zeta \{\sin(\Theta - \beta) - \sin i \sin(\Theta_{\mathrm{per}} + \beta)\} - \sin \zeta \cos i]^2}. \quad (\mathrm{A10}) $$

surface-normal unit vector at the point is $e'_R(\phi, \lambda) = (\cos \phi \cos \lambda, \sin \phi \cos \lambda, \sin \lambda)$. One can transform $e'_R$ to $e_R$

Finally, we obtain the eccentric frequency modulation in terms of the true anomaly $\Theta$:

$$ f_{\mathrm{orb}} = f_{\mathrm{spin}} + \frac{f_{\mathrm{orb}} (1 - e^2)^{-3/2} (1 + e \cos \Theta)^2 [-\cos \zeta \{1 + \sin i \cos(\Theta + \Theta_{\mathrm{per}})\} + \sin \zeta \cos i \sin(\Theta - \beta)]}{[\cos(\Theta - \beta) + \sin i \cos(\Theta_{\mathrm{per}} + \beta)]^2 + [\cos \zeta \{\sin(\Theta - \beta) - \sin i \sin(\Theta_{\mathrm{per}} + \beta)\} - \sin \zeta \cos i]^2}. \quad (\mathrm{A11}) $$
The above equation reproduces Equation (2) for $e \to 0$, $\beta \to 0^\circ$, $\Theta_{\text{per}} \to \Theta_{\text{eq}}$, and $\Theta \to \Theta - \Theta_{\text{eq}}$.

In the main part of the paper, we consider the circular orbit alone just for simplicity, but the effect of *e* is important as well. To show this, we plot Equation (A11) for the Earth-like planet viewed from $i = 0^\circ$ for $e = 0$ and 0.2 in Figure A2. The horizontal axis indicates the time in units of $P_{\text{orb}}$ that is numerically computed from the true anomaly. The left and right panels correspond to the obliquity of $\zeta = 30^\circ$ and $\zeta = 60^\circ$. Different curves indicate cases for $\beta = 0^\circ, 90^\circ, 180^\circ$, and $270^\circ$.

The dashed ($e = 0$) and solid ($e = 0.2$) lines clearly exhibit different amplitudes and phases of the frequency modulation. Thus the effect of eccentricity biases the estimates of $\zeta$ and $\Theta_{\text{eq}}$ if the formula for $e = 0$ is used in fitting the data. Since it is likely that the orbit of direct-imaging targets is precisely determined prior to its monitoring, one can use Equation (A11) as a template using the estimated value of $e$.

## ORCID iDs

Yuta Nakagawa https://orcid.org/0000-0002-1004-685X

Takanori Kodama https://orcid.org/0000-0001-9032-5826

Masaki Ishiwatari https://orcid.org/0000-0001-7490-4676

Hajime Kawahara https://orcid.org/0000-0003-3309-9134

Yasushi Suto https://orcid.org/0000-0002-4858-7598

Yoshiyuki O. Takahashi https://orcid.org/0000-0003-4060-7379

George L. Hashimoto https://orcid.org/0000-0002-3821-6881

Kiyoshi Kuramoto https://orcid.org/0000-0002-6757-8064

Kensuke Nakajima https://orcid.org/0000-0002-0333-2364

Shin-ichi Takehiro https://orcid.org/0000-0002-0857-3992

Yoshi-Yuki Hayashi https://orcid.org/0000-0003-4745-4520

## References

Aizawa, M., Kawahara, H., & Fan, S. 2020, ApJ, 896, 22

Anderson, G. P. 1986, Optical Physics Division, U.S. Air Force Geophysics Laboratory, Environmental Research Papers, NO. 954, Tech. Rep. AFGL-TR-86-0110

Beljaars, A. C. M. 1995, QJRMS, 121, 255

Beljaars, A. C. M., & Holtslag, A. A. M. 1991, JApMe, 30, 327

Cohen, L. 1995, Time-frequency Analysis, Prentice Hall Signal Processing Series (Englewood Cliffs, NJ: Prentice-Hall), http://ci.nii.ac.jp/ncid/BA25293241

Cowan, N. B., Agol, E., Meadows, V. S., et al. 2009, ApJ, 700, 915

Deitrick, R., Barnes, R., Bitz, C., et al. 2018, AJ, 155, 266

Emde, C., Buras-Schnell, R., Kylling, A., et al. 2016, GMD, 9, 1647

Fan, S., Li, C., Li, J.-Z., et al. 2019, ApJL, 882, L1

Farr, B., Farr, W. M., Cowan, N. B., Haggard, H. M., & Robinson, T. 2018, AJ, 156, 146

Ford, E. B., Seager, S., & Turner, E. L. 2001, Natur, 412, 885

Fujii, Y., & Kawahara, H. 2012, ApJ, 755, 101

Fujii, Y., Kawahara, H., Suto, Y., et al. 2010, ApJ, 715, 866

Fujii, Y., Kawahara, H., Suto, Y., et al. 2011, ApJ, 738, 184

Gasteiger, J., Emde, C., Mayer, B., et al. 2014, JQSRT, 148, 99

Greengard, L., & Lee, J.-Y. 2004, SIAMR, 46, 443

Hu, Y. X., & Stamnes, K. 1993, JCLI, 6, 728

Jiang, J. H., Zhai, A. J., Herman, J., et al. 2018, AJ, 156, 26

Kaspi, Y., & Showman, A. P. 2015, ApJ, 804, 60

Kawahara, H. 2016, ApJ, 822, 112

Kawahara, H., & Fujii, Y. 2010, ApJ, 720, 1333

Kawahara, H., & Fujii, Y. 2011, ApJL, 739, L62

Komacek, T. D., Jansen, M. F., Wolf, E. T., & Abbot, D. S. 2019, ApJ, 883, 46

Le Treut, H., & Li, Z.-X. 1991, CIDy, 5, 175

Manabe, S. 1969, MWRv, 97, 739

Markwardt, C. B. 2009, in ASP Conf. Ser. 411, Astronomical Data Analysis Software and Systems XVIII, ed. D. A. Bohlender, D. Durand, & P. Dowler (San Francisco, CA: ASP), 251, arXiv:0902.2850

Mayer, B., & Kylling, A. 2005, ACP, 5, 1855

Mellor, G. L., & Yamada, T. 1982, RvGSP, 20, 851

Moorthi, S., & Suarez, M. J. 1992, MWRv, 120, 978

Nakajima, T. 1983, JQSRT, 29, 521

Noda, S., Ishiwatari, M., Nakajima, K., et al. 2017, Icar, 282, 1

Oakley, P. H. H., & Cash, W. 2009, ApJ, 700, 1428

Pallé, E., Ford, E. B., Seager, S., Montañés-Rodríguez, P., & Vazquez, M. 2008, ApJ, 676, 1319

Robinson, T. D., Meadows, V. S., & Crisp, D. 2010, ApJL, 721, L67

Rose, B. E. J. 2015, JGRD, 120, 1404

Rushby, A. J., Shields, A. L., & Joshi, M. 2019, ApJ, 887, 29

Sagan, C., Thompson, W. R., Carlson, R., Gurnett, D., & Hord, C. 1993, Natur, 365, 715

Salomonson, V. V., Barnes, W. L., Maymon, P. W., Montgomery, H. E., & Ostrow, H. 1989, ITGRS, 27, 145

Schwartz, J. C., Sekowski, C., Haggard, H. M., Pallé, E., & Cowan, N. B. 2016, MNRAS, 457, 926

Stamnes, K., Tsay, S.-C., Jayaweera, K., & Wiscombe, W. 1988, ApOpt, 27, 2502

Suto, Y. 2019, in How to Search for Possible Bio-signatures on Earth-Like Planets: Beyond a Pale Blue Dot, ed. A. Yamagishi, T. Kakegawa, & T. Usui (Singapore: Springer), 441

Wanner, W., Li, X., & Strahler, A. H. 1995, JGRD, 100, 21077

Williams, D. M., & Pollard, D. 2003, IJAsB, 2, 1