Disk-integrated brightness of the sunlit portion of Earth (and occasionally the Moon) observed by EPIC serves as an ideal proxy for point-source measurements, analogous to observations of distant exoplanets (Jiang et al., 2018). The 388 nm radiance (Figure 1) is predominantly influenced by highly reflective clouds, with minimal contribution from the darker surface background, thus serving as an effective proxy for cloud distribution. In contrast, the 780 nm channel captures significant signals from both clouds and surface features. Fourier analysis using Lomb–Scargle periodograms (Lomb, 1976; Scargle, 1982; Hans et al., 1999) reveals distinct periodic signals in these channels, capturing the temporal variability driven by Earth’s rotation, orbital configuration, and the spatial distribution of clouds and surface features.

The qualitative results (Jiang et al., 2018) show that the variability of Earth’s light curve is largely influenced by the distribution of landmasses, oceans, and clouds. In the UV and visible bands, the reflected radiance is dominated by clouds and enhanced whenever clouds are present in the EPIC FOV. Surface contributions (i.e., positive landmass and negative ocean signals) are mixed with cloud signals in the visible wavelengths and become more prominent in the near-infrared channels. Earth appears brighter during the Southern Hemisphere (SH) summer due to its orbital eccentricity (∼0.017), which brings it closer to the Sun during that season. Additionally, the prominent feature in the Fourier spectra is the 24 h (1 day) peak (Figure 2), reflecting Earth’s 24 h rotation period. A few notable features with periods shorter than 24 h appear across all wavelengths at 12, 8, and 6 h, with the 12 h peak especially strong in the UV (Figures 2b,d). These signals likely result from periodic patterns of the Earth’s surface and cloud systems, which repeat daily when the Earth rotates. The period of 12 h could be a combined signature from clouds, oceans, and landmass patterns. For example, Asia and North America enter the FOV roughly 12 h apart in local time, and the Western Pacific faces EPIC around 02:00–03:00 UTC, followed by the Atlantic Ocean around 14:00–15:00 UTC. The 8-h and 6-h components may arise from the spatial variation of the cloud system and the continents, such as the ∼8 h difference between western Australia and northeastern Africa, and the ∼6 h interval between South America and southern Africa as they sequentially enter the EPIC FOV.

A quantitative analysis of the individual contributions from these components is essential for a comprehensive understanding of the light curve. While the spectral properties of known surface types can be estimated under different stellar spectral energy distributions (Shields et al., 2013), the specific surface composition of exoplanets remains unknown. Therefore, to maintain generality in the absence of such information, all types of land surfaces are treated as equivalent and represented with an “average” land surface (Fan et al., 2019). Acknowledging that future exoplanet observations may span different wavelength ranges, the reflectance time series is scaled to have zero mean and unit variance prior to further analysis. This normalization ensures equal weighting across spectral channels, which is particularly important for multivariate techniques such as singular value decomposition (SVD) and Gradient Boosted Regression Trees (GBRT; Friedman, 2001), which are sensitive to differences in data scale. The input dataset consisted of 10 reflectance channels along with 2 fractional coverage labels (land and cloud fraction, Figure 3) at each time step. In contrast, unnormalized reflectance values can be retained to preserve physical brightness information from different wavelengths. Additionally, EPIC Level 2 products, which are derived based on well-characterized properties of Earth (e.g., the distribution of surface types, the composition of the surface and atmosphere, and the spectral signatures of various spatial features), provide an opportunity to construct physically grounded models of planetary reflectance. The spatially resolved EPIC-view land cover (NASA/LARC/SD/ASDC 2017) and cloud data (NASA/LARC/SD/ASDC 2018), classified as the more specific fractional contributions of six distinct spatial features (i.e., ocean, desert, snow/ice, vegetation, and low and high clouds) were incorporated as additional indicators, resulting in a time series of 10 reflectance channels and 6 fractional coverage labels at each time step (Gu et al., 2021).

The GBRT results (Figure 3) as well as the eigenvalues of PCs (Figure 4a) suggest that PC1 and PC2 are the most important components with distinct contributions to land and cloud features, while the contributions from other PCs are relatively negligible (Figure 3; Fan et al., 2019). For land fraction, PC2 is dominant with a weight of 0.88, whereas PC1 contributes only 0.03. The variation in land/ocean fraction is thus primarily captured by PC2, showing a strong correlation of 0.91. Both PCs contribute substantially to cloud fraction, with weights of 0.28 for PC1 and 0.49 for PC2. The GBRT results also show that PC1 is strongly associated with low cloud, showing a flat spectral shape from UV to NIR wavelengths, consistent with cloud spectra, and a prominent peak corresponding to low clouds in observed data (∼70% relative importance) (Gu et al., 2021). The time series of PC1 is highly correlated with low cloud fraction (R = −0.8). PC2 captures information about land and ocean surfaces, with a spectrum that increases in reflectance with wavelength, matching their spectral behavior (Figure 4). GBRT shows peak importance for land/ocean labels (>84%), except for snow/ice, which may be due to its relatively low coverage on Earth’s surface (∼10%), confinement to high latitudes (resulting in smaller weight for disk integration), and/or slow temporal change. PC2 is strongly correlated with land/ocean contrast (R = −0.96). PC4 contains high cloud information, showing a flat spectrum with oxygen absorption features (688 and 764 nm). While GBRT indicates high cloud dominance (>71%), with a strong correlation (R = −0.82) between high cloud fraction and PC4, PC3 remains unclear but may relate to atmospheric Rayleigh scattering or ozone absorption due to its steep UV slope. PC5 shows peaks for vegetation/desert labels, suggesting it reflects the vegetation-desert contrast. Its spectrum (Figure 4b) is consistent with the VRE, and in experiments with only vegetation or desert, PC5 decreases significantly, providing weak evidence of its link to vegetation.

Given the strong correlation between PC2 and surface features (Figure 3), and the comparable importance of PC1 and PC2 for clouds, it is inferred that clouds consist of two components: one independent of surface type, and the other associated with surface features. This interpretation aligns with previous studies that incorporate surface type into cloud albedo calculations (Thompson and Barron, 1981; Biasiotti et al., 2022). Specifically, PC2 exhibits smoother and more consistent daily changes (Figures 2g,h), while PC1 shows greater variability across days (Figures 2e,f). PC2 also tends to peak during summer, when the sunlit side of the Northern Hemisphere, with extensive landmass, is visible to both the Sun and the spacecraft. Its consistent diurnal variation likely reflects the combined effect of Earth’s rotation and the asymmetry in surface distribution, supporting the conclusion that the cloud variability is temporally correlated with surface features.

PC2 primarily captures the surface information about Earth with a strong linear correlation, thus making it possible to reconstruct a two-dimensional (2D) surface map. To facilitate generalization for future observations of Earth-like exoplanets, a set of simplified assumptions was adopted (Fan et al., 2019; Fan and Yung, 2020): the spectral properties of reflective surfaces were treated as unknown, the incoming solar flux was assumed to be uniform and known, and the entire planetary surface was modeled as a Lambertian reflector. Although Earth’s oceans are known to be strongly non-Lambertian, this assumption allows the method to be applied more broadly when surface characteristics are uncertain. The viewing geometry was assumed to be known and derived from DSCOVR navigation data over the 2-year observation period. Under these assumptions, the surface reconstruction becomes a linear regression problem (Fan et al., 2019). The reconstructed map (Figure 5a) shows PC2 values are positively correlated with land fraction. Coastlines are defined using the median PC2 value, aligning with the assumption that the total land coverage is unknown. Compared with the global land/ocean map of Earth (Figure 5b), the reconstructed map successfully captures all major continental features.

Building upon the framework (Fan et al., 2019; Fan and Yung, 2020), further studies have proposed enhancements to improve mapping fidelity and flexibility. Since global mapping is inherently ill-posed, Aizawa et al. (2020) applied sparse modeling techniques to improve the spatial resolution of inferred maps. Besides, Kawahara (2020) formulated a single inverse problem that unifies spin-orbit tomography and NMF (nonnegative matrix factorization)-based spectral unmixing, enabling the simultaneous retrieval of both spectral and geographical information from disk-integrated exoplanet light curves. Moreover, to account for temporal variability, Kawahara and Masuda (2020) developed a dynamic mapping approach by extending the Bayesian framework of Farr et al. (2018), utilizing analytical solutions to the Bayesian inverse problem with (multivariate) Gaussian priors and their isomorphic representations, thereby enabling the modeling of time-varying geography.

Since exoplanets will remain spatially unresolved point sources for the foreseeable future, quantifying the contribution of each spatial component to the integrated light curve is crucial for characterizing their habitability. Although PC2 is primarily interpreted (Fan et al., 2019) in terms of land and cloud fractions, Earth, as a well-observed and well-understood planet, provides a unique opportunity for more detailed modeling beyond such general correlations. Leveraging the availability of high-resolution data on Earth’s surface composition and atmospheric constituents, it becomes feasible to simulate its reflective signatures through physical models that account for diverse surface types, cloud structures, and atmospheric scattering effects.

As discussed above, the EPIC Level 2 products provide the opportunity to clarify 6 spatial features. By using a two-dimensional moving average applied to the pixel-level reflectance data within the phase space defined by the solar and viewing zenith angles, the reflected intensity for each spatial feature at a given wavelength is calculated. This enables the construction of a temporally varying Earth image model, which can then be used to reconstruct DSCOVR observations, analyze their relationship with the spatial features, and interpret the physical meaning of the principal components (Gu et al., 2021).

Assumed that the entire sunlit disk is uniformly covered by each of the six spatial features individually, the corresponding disk-averaged reflection spectra for a synthetic Earth (Figures 6a,b) provides characteristic spectral profiles associated with reflecting surfaces (Gu et al., 2021). For each time step, the reflectance contribution of a given spatial feature was calculated by dividing its total reflectance by the disk-integrated reflectance and a daily average of these contributions was then applied (Figure 6c). The slopes in the first three UV channels result from Rayleigh scattering and ozone absorption. Absorption features from the oxygen A and B bands at 764 nm and 688 nm influence the spectral shape, leading to minimal spectral differences among spatial features. Oxygen absorption is weaker in high clouds compared to low clouds (or the surface) due to shorter optical path lengths, suggesting the possibility of distinguishing high and low clouds in the disk-integrated spectrum. The spectral shape of ice/snow reflectance closely matches that of low clouds, making it difficult to identify their contribution to individual light curves due to smaller coverage. Vegetation is the only feature showing a monotonic increase in reflectance at the longest wavelengths, consistent with the VRE (i.e., 680, 688, and 764 nm). In the first four UV channels (Figure 6c), ocean and cloud contributions dominate (around 40%), with atmospheric scattering and ozone absorption influencing the reflectance (Jiang et al., 2018), making it hard to distinguish surface features. At longer wavelengths, low clouds contribute up to ∼60%, while high clouds are more significant at 688 and 764 nm due to stronger oxygen absorption in the lower atmosphere. Vegetation shows higher contributions in the NIR channels, consistent with the VRE. The seasonal cycle, with more vegetation and less ocean in the northern summer, is also apparent in the time series.

The eigenvalues (Figure 4a) of PC1 and PC2 indicate that they together explain over ∼97% of the total variance, with values nearly an order of magnitude larger than that of PC3 (Gu et al., 2021). PC3–6 form a secondary group, potentially capturing weaker yet meaningful spatial signals. Components beyond PC6, each contributing less than ∼0.02% of the variance, are treated as noise and excluded from further analysis. Despite slightly lower eigenvalues, the synthetic single-point light curves successfully reproduce all characteristic spectral features, except for minor switching between the observed and reconstructed eigenvectors for PC3 and PC4 (Figure 4b). The result is promising, especially considering that the reconstruction was achieved using only six spatially distinct yet temporally invariant surface and atmospheric components.

Although a non–radiative transfer (non-RT) Earth model was employed (Gu et al., 2021) and successfully demonstrated that information related to low clouds, surface type distribution, and high clouds is encoded in different PCs of Earth’s light curves, the absence of radiative transfer processes in such models still limited their ability to assess the influence of atmospheric composition and scattering. Previous studies (e.g., Tinetti et al., 2006a; Tinetti et al., 2006b; Robinson et al., 2011; Feng et al., 2018; Batalha et al., 2019) have attempted to model Earth’s light curves using RT models, yet these efforts lack validation against observations extending beyond a few days. Consequently, the development of a comprehensive RT model of Earth, rigorously tested against long-term and diverse observations, remains both critical and pressing.

To investigate the influence of atmospheric radiative transfer on Earth’s disk-integrated reflectance, the Earth Spectrum Simulator (ESS) is developed (Gu et al., 2022). The ESS model incorporates the two-stream-exact-single-scattering (2S-ESS) line-by-line RT model that combines precise single-scattering calculations with a two-stream approximation for multiple scattering and offers a balance between computational efficiency and accuracy (Spurr and Natraj, 2011) and has been widely applied in Earth remote sensing of greenhouse gases and aerosols (Xi et al., 2015; Zhang et al., 2015; 2016; Zeng et al., 2017; 2018; Huang et al., 2020; Natraj et al., 2021), thus accounting for surface reflection, cloud and aerosol scattering, and gaseous absorption.

In the 2S-ESS model, surface reflectance was simulated through BRDFs (Bidirectional Reflectance Distribution Functions, Figure 7a), and then land surfaces were categorized into seven types (NASA/LARC/SD/ASDC, 2017). For snow- and ice-covered pixels, a Lambertian surface albedo was assumed. Additionally, the Rayleigh scattering and molecular absorption were simulated considering only the influence of O₂, H₂O, and O₃, due to their dominant influence in the narrow spectral bands of EPIC. The vertical water vapor and temperature were assumed spatially invariant while surface pressure varied spatially and temporally. The surface pressure (Figure 7c) data was derived from MERRA-2 (Modern-Era Retrospective analysis for Research and Application, Version 2). The 3D ozone distributions were derived (Figure 7d), and all atmospheric fields were sampled on the 5th, 15th, and 25th days of each month to reduce computational demands. Aerosol optical properties (Figure 7e) were calculated using Mie scattering theory, except for dust, where nonsphericity is considered (Meng et al., 2010; Crisp et al., 2021). Clouds (Figure 7f) were further classified as either liquid or ice, with one type per pixel and fixed vertical distributions and their optical properties were based on Mie theory (Hess et al., 1998) for liquid clouds and a cirrus habit mixture (Baum et al., 2011) for ice clouds. For pixels with mixed cloud types at degraded resolution, multiple RT simulations were performed and averaged based on cloud fraction to obtain the final reflectance.

Generally, the ESS can reproduce the major features observed accurately. In the shortest wavelength channels, strong Rayleigh scattering results in generally bright images, with the 340 nm channel being the brightest (Figure 8), which is consistent with the brightness variations caused by the enhanced ozone absorption between 340 nm and 317 nm. In the longer wavelength channels, surface information becomes more clearly visible, except for the two oxygen absorption channels (688 nm and 764 nm). In particular, due to the VRE, the reflectance of vegetated pixels increases significantly in the 780 nm channel, forming a sharp contrast with desert pixels, whose reflectance gradually increases across the 551 nm, 680 nm, and 780 nm channels. These spectral differences among different surface types benefit from the detailed modeling of land surface BRDF. In terms of atmospheric scattering features, the use of MAIAC data enables effective simulation of the backscattering peak observed under EPIC viewing geometry; meanwhile, the adoption of the Cox–Munk surface model allows the simulation to capture the feature of ocean glint (Figure 7b).

Specifically, the ESS model successfully reproduces the time-averaged spectra with errors generally below 5%, and within 7% overall, except at the two oxygen absorption bands (688 nm and 764 nm), where simulated reflectance exceeds observations by up to 19.53%. The model exhibits slightly larger temporal variability than observations. In the 317–443 nm range, the simulated reflectance tends to be slightly darker than observed, while simulations at longer wavelengths generally appear brighter. A pronounced seasonal cycle is captured in the Southern Hemisphere, with higher reflectance during summer, likely due to the closer Earth–Sun distance and the presence of the Antarctic ice sheet. A smaller enhancement is also observed in the Northern Hemisphere summer, which may be attributed to increased land surface reflectance and cloud cover driven by the South Asian monsoon (Jiang et al., 2018).

The comparison between the simulated and observed PCs (Figure 4) includes an additional subsampled dataset of approximately 300 time points, selected to align with the simulation output. This subset captures the main temporal variability of the light curves, as indicated by the relatively small amplitudes of higher-order singular values (Figure 4a; Gu et al., 2021). At the same time points, the singular values from the simulations are generally larger than those from the observations, suggesting stronger temporal variability in the simulated light curves. Notably, PC3 and PC4 appear to be interchanged between the simulation and observation (Figure 4b). The simulation successfully reconstructs PC1-5; in particular, the first two PCs show near-perfect agreement, demonstrating that the ESS model accurately reproduces the two most dominant signals, namely, low cloud and surface reflectance patterns, which is in line with the previous findings (Gu et al., 2021). The observed PC4 and PC3 also correspond well to the simulated PC3 and PC4, respectively. Since observed PC4 is primarily associated with high clouds (Gu et al., 2021), and ESS tends to simulate higher reflectance for high clouds, the corresponding eigenvalue increases due to the enhanced variance contribution from high cloud coverage.

Although the details of cloud features show some differences, mainly because the simulation uses the DSCOVR composite cloud dataset, which does not exactly match the instantaneous cloud distribution observed by EPIC, the ESS model still serves as a valuable tool for generating light curves of Earth-like exoplanets and provides an important platform for combining with global climate models (GCMs) to support exoplanet characterization studies with its high fidelity, robustness, and efficiency.

In the context of measuring exoplanet complexity, two key aspects must be considered: 1) how to formally quantify complexity, and 2) how to use remotely detectable signals to assess the complexity of a given exoplanet. To address the first, the formal epsilon machine (EM) framework has been in development for several decades (Crutchfield, 2012), and was applied to the challenge of exoplanet characterization (Bartlett et al., 2022). EMs are discrete mathematical objects (finite state machines) that can compactly represent both simple random and simple deterministic processes (as contrasted with minimal Turing machine programs, which ascribe high complexities to simple random processes). They also have the capacity to reveal the underlying phase-space structure of nonlinear systems (Brodu, 2011; Packard et al., 1980; Crutchfield et al., 1986). EMs are particularly well suited to the analysis of planetary light curves, as they are designed to fulfill several key criteria: 1) reproduce deterministic features of the input time series; 2) reproduce random stochastic features of the input time series; 3) serve as the optimal predictor of the generating process; and 4) represent the internal structure of the generating process in the most compact form (Brodu, 2011; Crutchfield, 2012). These kinds of algorithms take a time series as input and produce a minimal optimized model (the EM) capable of reproducing time series that are, subject to parameter tuning and data quality, etc., statistically equivalent to the original. The size of the resulting model, particularly the information content of its state-space distribution, defines a metric known as statistical complexity.

To address the second challenge, it is hypothesized that planets exhibiting the highest complexity are the most likely to host biospheres, with Earth serving as the initial reference. Even if lifeless planets with high complexity are identified, the approach remains valuable as a characterization tool, given the largely unknown space of exoplanet types, which could vary widely in complexity. For example, the Moon exhibits very low complexity, while the Sun and Jupiter display somewhat richer structures and variability. Earth, with its active atmosphere, hydrosphere, and biosphere, appears to be the most complex body in the Solar System (Kasting and Siefert, 2002; Dyke et al., 2011; Olejarz et al., 2021), although this has yet to be quantitatively confirmed. By using EMs derived from reflectance time series, exoplanet complexity levels can be estimated, which could serve as a potential biosignature.

As an initial statistical assessment, the Shannon entropy was calculated for each proxy exoplanet reflectance time series (Bartlett et al., 2022). The reflectance values were discretized and normalized, and the entropies were computed using the standard formula based on the frequency of occurrence for each reflectance value. These entropies were then averaged over the ten EPIC wavelength channels. On the complexity-entropy plane (Figure 9a), groups are observed along a monotonic slope of increasing complexity and entropy despite inherent noise and data sampling issues. The simplest synthetic Earths (cloudless, mono-surface worlds) exhibit the smallest Shannon entropy, while the highest entropy is observed in the most qualitatively and quantitatively complex case (the original Earth). Within the cloudy worlds, there is a relatively small distinction between the mono-surface and multi-surface groups, but it is promising to observe any distinction despite the cloud cover’s masking effect.

To further explore the potential for distinguishing living from non-living planets, the gas giant, Jupiter, which is generally assumed to be sterile and known for its complex atmosphere that exhibits both periodic and stochastic features, was compared with Earth (Bartlett et al., 2022; Figure 9b). For each wavelength considered, Earth exhibited higher complexity and entropy than Jupiter, with an average difference of approximately 1.25 bits for complexity and 1 bit for entropy, although these values vary depending on the wavelength. Even though the lower image resolution of the Cassini data may have affected the information content of the Jupiter data, these results still suggest that Earth’s planetary complexity exceeds that of Jupiter, which is consistent with the hypothesis that planetary complexity correlates with the presence of life. However, further investigation is required to fully understand the causal influences leading to this observed correlation.

When analyzing an object using discrete data, it’s possible that important details in the original measurements are lost during the discretization (binning) step. In contrast, methods applicable to continuous time series data may reveal important, perhaps subtle features that could aid in planetary characterization. This was the motivation for Segal et al. (2024) to employ an approximation for the entropy of a continuous variable to DSCOVR EPIC reflectance time series. Specifically, they used the Vasicek estimate of the differential entropy (Vasicek, 1976). PCA applied to the spatially unresolved multiwavelength reflectance data for each synthetic Earth used by Bartlett et al. (2022) and the average differential entropy (Figure 9c) was derived. This was found to be approximately proportional to the joint differential entropy (JDE). Figure 9d shows that the average JDE of eigencolors, obtained via PCA, reveals a new level of distinctions between the synthetic Earth time series data. Unlike the original wavelength channels, eigenvectors are linearly independent which means that cross-wavelength mutual information, such as that introduced by cloud cover, is only accounted for once (Segal et al., 2024). This is supported by the average Pearson correlation coefficient among the 10 eigencolors, which is on the order of 10^–17 for both Earth and the reconstructed worlds. In contrast, when the correlation is computed across the original EPIC narrow band wavelengths, the average value is approximately R = 0.57 for the original Earth data and decreases to R = 0.08 for the cloud-removed version. These results are consistent with previous findings that reflectance variability due to clouds is across all 10 EPIC bands (Jiang et al., 2018).