Jupiter as an Exoplanet: Insights from Cassini Phase Curves

Due to its proximity to Earth, Jupiter of the Solar System serves as a unique case study for gas-giant exoplanets. In the current study, we perform fits of ab initio, reflective, semi-infinite, homogeneous model atmospheres to 61 phase curves from 0.40 to 1.00 $\mu$m, obtained from the Cassini spacecraft, within a Bayesian framework. We reproduce the previous finding that atmospheric models using classic reflection laws (Lambertian, Rayleigh, single Henyey-Greenstein) provide poor fits to the data. Using the double Henyey-Greenstein reflection law, we extract posterior distributions of the single-scattering albedo and scattering asymmetry factors and tabulate their median values and uncertainties. We infer that the aerosols in the Jovian atmosphere are large, irregular, polydisperse particles that produce strong forward scattering together with a narrow backscattering lobe. The near-unity values of the single-scattering albedos imply that multiple scattering of radiation is an important effect. We speculate that the observed narrow backscattering lobe is caused by coherent backscattering of radiation, which is usually associated with Solar System bodies with solid surfaces and regolith. Our findings demonstrate that precise, multi-wavelength phase curves encode valuable information on the fundamental properties of cloud/haze particles. The method described in this Letter enables single-scattering albedos and scattering asymmetry factors to be retrieved from James Webb Space Telescope phase curves of exoplanets.


INTRODUCTION
Named after the Roman god of the sky and thunder, Jupiter is the most massive planet of our Solar System, has been studied for several centuries and was the subject of close scrutiny by recent space missions (Porco et al. 2004;Bolton et al. 2017). Since gas giants were the first exoplanets to be discovered (Mayor & Queloz 1995) and their atmospheres remain the most well-characterised (Deming & Seager 2017), Jupiter holds a special place among Solar System planets as a unique case study for exoplanets and their atmospheres.
The visible and near-infrared phase curves of Jupiter are of particular interest, because they quantify the fraction of sunlight reflected by the Jovian atmosphere as a function of the orbital phase angle. Unlike for ∼ 1000 K hot Jupiters, the Jovian phase curves are not contaminated by the thermal (infrared) emission of Jupiter and can be safely assumed to comprise predominantly reflected sunlight. Using data from the Cassini space mission, Dyudina et al. (2016) and Mayorga et al. (2016) previously showed that the Jovian phase curves are "cuspy" and peak more sharply towards zero phase angle than classic laws of reflection (Lambertian, Rayleigh). Figure 1 shows an example, where we reproduce the findings of Dyudina et al. (2016) and Mayorga et al. (2016). Neither study performed fits for fundamental physical param-  Examples of failed fits of model phase curves of semi-infinite, homogeneous atmospheres to the 0.50 µm Cassini phase curve of Jupiter assuming the single Henyey-Greenstein and Rayleigh reflection laws. The Lambertian reflection law is not a fit as it does not possess any free parameters (i.e., ω = 1 and Ag = 2/3).
eters (single-scattering albedo, scattering asymmetry factor) within a Bayesian framework, as such numerical calculations are computationally expensive.
Recently, Heng et al. (2021) reported first-principles, fully analytical solutions of reflected light phase curves for any reflection law. In the current Letter, we exploit this develop-ment to fit ab initio 1 models of semi-infinite, homogeneous atmospheres to a set of 61 Cassini phase curves of Jupiter (Li et al. 2018) from wavelengths of λ = 0.40-1.00 µm (with bin sizes of 0.01 µm) within a Bayesian framework (Foreman-Mackey et al. 2013). We demonstrate that the double Henyey-Greenstein (DHG) reflection law (Henyey & Greenstein 1941;Kattawar 1975; provides reasonable fits to the Jovian phase curves, which indicates the scattering of light by large, irregular particles (McGuire & Hapke 1995) that are poorly described by Mie theory (Mie 1908;Kitzmann & Heng 2018).
In Section 2, we describe our methods and data. In Section 3, we demonstrate that the set of fits yields inferred values (with uncertainties) of the single-scattering albedo and scattering asymmetry factors at each wavelength (Table 1), which provide important input for workers studying the chemistry of Jovian clouds and hazes. In Section 4, we compare our work to previous studies and discuss its implications.

Background theory
Building on the classic work of Chandrasekhar (1960), Sobolev (1975) and Hapke (1981), Heng et al. (2021) derived novel analytical solutions for the geometric albedo A g and integral phase function Ψ of semi-infinite atmospheres for any law of reflection (scattering phase function P ). For context, some of the key findings of Heng et al. (2021) are concisely reviewed here.
Let incident sunlight with intensity I impinge upon an atmosphere at a zenith angle µ . Let the intensity of light reflected by the atmosphere be I 0 , which is generally a formal solution of the radiative transfer equation (Chandrasekhar 1960;Hapke 1981). The reflection coefficient is (Dlugach & Yanovitskij 1974;Sobolev 1975;Seager 2010) The quantity ρ/π is sometimes termed the "bidirectional reflection distribution function" (BRDF; e.g., Dyudina et al. 2016).
The geometric albedo may be expressed in terms of the reflection coefficient (Dlugach & Yanovitskij 1974;Sobolev 1975), where µ = cos θ and θ is the polar angle in the frame of reference of the planet. In the observer-centric coordinate 1 There are no "tuning parameters" beyond one's choice of the reflection law (scattering phase function P ), since the formulae for Ag and Ψ are formal solutions of the radiative transfer equation (Heng et al. 2021). The absence of parameters to finetune implies the ability to extract or retrieve fundamental physical parameters associated with the scatterers, which may take the form of atoms, molecules, ions or aerosols (clouds/hazes). system, the flux (divided by I ) observed at any orbital phase angle α is (Sobolev 1975), where Θ and Φ are the observer's latitude and longitude, respectively. If we define F 0 ≡ F (α = 0 • ), then the integral phase function is (Sobolev 1975) The approach of Hapke (1981) is followed in assuming that single scattering is described by any scattering phase function P , but multiple scattering occurs isotropically. Building on this rich body of work, Heng et al. (2021) showed that closed-form solutions for A g and Ψ may be derived for any P . Within this formalism, the phase curve is A g Ψ. In other words, both the shape and normalisation of the reflected light phase curve are computed from first principles.

Bayesian framework for data fitting
To compute A g Ψ requires the specification of 4 free parameters: g 1 , g 2 , f and the single-scattering albedo ω. We perform fits of  the possibility that the data uncertainties have been underestimated by adding (δA g Ψ) 2 to the variance, where ln δ is formally part of the fit (Hogg et al. 2010) 4 . Uniform prior distributions are assumed: 0 ≤ ω ≤ 1, 0 ≤ g 1 ≤ 1, −1 ≤ g 2 ≤ 0, 0 ≤ f ≤ 1 and −10 < ln δ < 1. The joint posterior distributions of the parameters are plotted using the corner routine written in the Python programming language (Foreman-Mackey 2016).

Cassini data of Jupiter
The global images of Jupiter recorded by the Cassini spacecraft are used to generate the light curves at visible and near-infrared wavelengths (Li et al. 2018). The Cassini spacecraft observed Jupiter from October 2000 to March 2001 en route to Saturn. Among the dozen scientific instruments of Cassini, the Imaging Science Subsystem (ISS; Porco et al. 2004) and the Visual and Infrared Mapping Spectrometer (VIMS; Brown et al. 2004) conducted observations of Jupiter at visible and near-infrared wavelengths. The ISS is an imager with two cameras and multiple filters onboard the Cassini spacecraft (Porco et al. 2004;Li et al. 2018). We mainly use the ISS observations because of their better spatial resolution and more complete coverage of phase angle compared to the VIMS observations. To our knowledge, the Cassini ISS provides the best coverage of phase angles among all available global images of Jupiter. The selected high-spatial-resolution ISS global images recorded have a range of phase angles from about 0 • to 140 • . This limit of 140 • comes from the need to prevent stray light from entering the camera while pointing close to the Sun. The ISS global images have observational gaps in phase angle, which are filled by a polynomial function (Li et al. 2018) using a least-squares method (Bevington & Robinson 2003). In addition, the ISS global images were recorded at discrete wavelengths and some observational gaps in wavelength exist. We interpolate the phase curves from the ISS wavelengths to all wavelengths from 0.40-1.00 µm (0.01 µm bin size) by using ground-based spectra measured at a phase angle of 6.8 • (Karkoschka 1998). Each phase curve has 141 data points (α = 0 • -140 • with 1 • resolution).
There are two dominant sources of uncertainties in our processing of Jupiter's light curves: 1. Calibrating the Cassini ISS global images; 2. Filling in the observational gaps in phase angle and wavelength. Calibration was performed using the Cassini ISS CALibration (CISSCAL) software (http://ciclops.org/sci/cisscal.php). The calibration uncertainties (i.e., uneven bit-weighting, bias subtraction, 2-Hz noise, dark current in the ISS cameras, bright/dark pixel pair artifacts from anti-blooming mode, flat-field artifacts) are discussed in detail in West et al. (2010) and Knowles et al. (2020). All of these uncertainties typically 4 https://emcee.readthedocs.io/en/stable/tutorials/line/ amount to a few percent of the calibrated radiance (Knowles et al. 2020). We apply the correction factors provided by the ISS calibration team to account for the aforementioned calibration uncertainties. The remaining calibration uncertainties should be less than one percent of the calibrated radiance. To account for the uncertainties associated with filling in the observational gaps in phase angle and wavelength, we use the fit residuals, which are the differences between the observed and fitted values at the ISS observed phase angles and wavelengths, to estimate the uncertainties associated with filling in the observational gaps. By fitting for ln δ, we formally allow for the possibility that the uncertainties have been underestimated (Hogg et al. 2010).  2018), albeit within a Bayesian framework: the Jovian phase curves are generally too cuspy near α = 0 • and are poorly fitted by the single Henyey-Greenstein and Rayleigh reflection laws. By definition, the Lambertian reflection law always has ω = 1 and A g = 2/3 (Sobolev 1975); it possesses no free parameters and a fit is not performed.

Fitting for fundamental physical parameters
In Table 1, the median values of ω, g 1 , g 2 , f and ln δ, along with their 1-σ uncertainties, are reported. Figure 2 shows examples of fits to data at 0.50 µm (the peak of the solar spectrum) and 0.89 µm (the MT3 methane absorption band; Table VIII of Porco et al. 2004). Generally, the reduced chisquare value of the fits is always much less than unity. Typically, we have δ e −4 ≈ 2%, but in some cases we have δ ∼ e −2.3 ≈ 10%. Figure 3 shows the retrieved values of ω, g 1 , g 2 and 1 − f as functions of wavelength. It is striking that sharp decreases in ω and g 1 , as well as an increase in g 2 (which is a negative quantity), coincide with the MT2 (0.723-0.730 µm) and MT3 (0.882-0.898 µm) methane absorption filters. This behavior is consistent with a reduction in the strength of scattering. It is less noticeable for the MT1 (0.615-0.623 µm), but it has been previously noted that methane absorption in this filter is weak (Li et al. 2006).
(6) Figure 3 shows a visualisation of the average scattering phase function, which is inconsistent with Rayleigh scattering. The particles exhibit strong forward scattering (f ∼ 1) together with a narrow backscattering lobe. In other words, most of the power is in the Henyey-Greenstein scattering phase function that quantifies forward scattering.  −0.16 . Uniform prior distributions were adopted: 0 ≤ ω ≤ 1, −1 ≤ g ≤ 1, 0 ≤ B0 ≤ 10, 0 ≤ h ≤ 10 and −10 ≤ ln δ ≤ 1. As there are hints of bimodality in some of the posterior distributions, we performed this particular fit for 10,000 Monte Carlo steps.

Interpretation
A detailed interpretation of the inferred values of the DHG parameters in the context of scattering theory and atmospheric chemistry is beyond the scope of the current Letter, since the particles are likely to be non-spherical and chemically heterogeneous; see, e.g., Zhang et al. (2013) and Guerlet et al. (2020) for recent discussions. However, several general statements may be made. Firstly, the general inference of f ∼ 1 implies that the narrow backscattering lobe is a higherorder correction to a single Henyey-Greenstein reflection law describing large particles that produce forward scattering. Secondly, the general inference that g 1 > 0 indicates the presence of large particles with sizes that are somewhat larger than the range of wavelengths probed. Thirdly, the nonmonotonic variation of g 1 with wavelength indicates the presence of a size distribution-the particles are polydisperse. If the particles were monodisperse, then g 1 would monotonically decrease with increasing wavelength. Fourthly, the high inferred values of ω ∼ 1 indicate that multiple scattering of sunlight is a significant effect (Hapke 1981) and are consistent with previous cloud models of the Jovian atmosphere (Sromovsky & Fry 2002 Dyudina et al. (2016) previously performed fits to Pioneer and Cassini data of Jupiter using the DHG reflection law. Instead of fitting A g Ψ to data, these authors directly fitted P to the reflection coefficient; see Figures 1 and 2 of Dyudina et al. (2016). Inferred values of the DHG fitting parameters are reported in Table 2 of Dyudina et al. (2016), but the associated uncertainties and posterior distributions are not reported as no formal uncertainty estimates were performed. Figures  3 and 4 of Dyudina et al. (2016) show comparisons of the phase curves of Jupiter to numerical calculations of A g Ψ, but no fits of A g Ψ to data are shown. Mayorga et al. (2016) and Li et al. (2018) fitted the phase curves with polynomials, but these fits do not allow ω, g 1 and g 2 to be extracted.

Coherent backscattering in the Jovian atmosphere?
The physics behind the cuspy profiles of the Jovian phase curves remain unelucidated. In Solar System bodies with solid surfaces and regolith, cuspy phase curves are caused by a combination of shadow hiding and coherent backscattering (Hapke et al. 1998). The same effects have been observed for the rings of Saturn (Déau et al. 2013).
Coherent backscattering is the phenomenon of multiply scattered light within a non-uniform medium interfering constructively to produce a brightness peak at zero phase angle (Hapke et al. 1993;Hapke 2002). It has been studied in a variety of contexts within physics (see Akkermans et al. 1988 and references therein).
To explore the possibility that coherent backscattering is the cause of the cuspy nature of the Jovian phase curves near α = 0 • , we multiply ρ by a simplified function that is given by equations (28) and (29) of Hapke (2002), where B 0 (magnitude of coherent backscattering) and h (ratio of λ/4π to transport mean free path) are two additional fitting parameters of the model. In order to keep the number of parameters at five (including ln δ), we adopt the single Henyey-Greenstein reflection law (f = 1, g 1 = g). Instead of fitting all 61 Jovian phase curves with the modification stated above, we focus on one example: 0.50 µm. The fit to data, residuals and posterior distributions of ω, g, B 0 and h are shown in Figure 4. Future studies should elucidate, from first principles, if coherent backscattering is a significant effect in gaseous atmospheres with large, irregular aerosols.

Implications for exoplanets
It is unknown if the large, irregular, polydisperse scatterers present in the Jovian atmosphere are also present in hot Jovian atmospheres. If they are, a possible obstacle with detecting the cuspy profiles of reflected light phase curves at zero phase angle is that no information is available during secondary eclipse. At the time of writing, the most robust detection of a reflected light phase curve is of the hot Jupiter Kepler-7b (Demory et al. 2013). Figures 2 and 3 of Demory et al. (2013) show that the duration of the secondary eclipse corresponds to about 5% of the orbit, which implies that information on any potentially cuspy profile is unavailable across −9 • α 9 • . Fits performed with data from 0 • ≤ α ≤ 9 • excluded indicate that this does not alleviate the difficulty with fitting the single Henyey-Greenstein reflection law to the Jovian phase curves (not shown). Cuspy profiles in the reflected light phase curves of gas-giant exoplanets remain to be detected.
There is a key difference between the interpretation of reflected light from Jupiter versus hot Jupiters. Since its rotational and orbital periods are about 10 hours and 12 years, respectively, Jupiter is a fast rotator. To lowest order, a fast rotator has longitudinal symmetry in terms of insolation and may be interpreted as a homogeneous sphere. The fact that all 61 Cassini phase curves considered in the current study peak at α = 0 • is strongly consistent with a homogeneous sphere. By contrast, tidally locked hot Jupiters do not possess longitudinal symmetry and are interpreted as inhomoge-neous spheres, which produce phase curves that do not peak at α = 0 • -an expectation that is borne out by both visible/optical photometric observations (Demory et al. 2013;Hu et al. 2015;Shporer & Hu 2015) and general circulation models (Oreshenko et al. 2016;Parmentier et al. 2016;Roman & Rauscher 2017).
The current investigation suggests that precise, multiwavelength phase curves that will be procured by the James Webb Space Telescope encode valuable information on the fundamental properties of clouds and hazes. Retrieving for the single-scattering albedos and scattering asymmetry factors using an ab initio model will motivate further studies into the microphysics and chemistry of clouds and hazes in exoplanetary atmospheres. In this context, the method described in the current Letter provides a crucial bridge between observations and simulations.
KH acknowledges financial support from the European Research Council (via a Consolidator Grant to KH; grant number 771620), the PlanetS National Center of Competence in Research (NCCR) and the Center for Space and Habitability (CSH), as well as a honorary professorship from Warwick University. KH is grateful to his senior Ph.D student Chloe Fisher for teaching him how to use emcee and corner and to his postdoctoral research associate Brett Morris for advice on data fitting. Note: Uncertainties are stated for 1σ (16th, 50th and 84th percentile).