EXOPLANET ALBEDO SPECTRA AND COLORS AS A FUNCTION OF PLANET PHASE, SEPARATION, AND METALLICITY

The contrast ratio between an exoplanet and its star depends on many factors, including the physical properties of the planet and star, the system geometry, background sources, and the instrumentation used. A simple way of estimating the required contrast is to assume that the planet reflects its incident starlight isotropically, such that its apparent brightness is constant over its illuminated surface, as a Lambertian surface. Then, the contrast can be expressed as (Sobolev 1975)

$$\begin{equation} C(\alpha) = \frac{2}{3}A_g(\lambda)\left(\frac{R_p}{\it d}\right)^2\left[\frac{\sin \alpha +(\pi -\alpha)\cos \alpha }{\pi }\right], \end{equation} \tag{ 1 }$$

where α is the planet phase angle, R_p is the planet radius in km, d is the planet–star separation in km, and A_g is the geometric albedo, which we formally define in Section 2.2. The geometric albedo generally takes values between 0 and 1 for the fraction of monochromatic light the planet reflects toward the observer at full phase, although it can be larger than 1 for anisotropic scattering atmospheres or surfaces. If one assumes that the planet is at quadrature, α = π/2, and that the geometric albedo at the wavelength of interest is 0.5, then for a Jupiter analog (1 R_J), the contrast at 0.8 AU is ∼3.8 × 10⁻⁸ and at 10 AU is ∼2.4 × 10⁻¹⁰. For a Neptune analog (1 R_N), the contrast at 0.8 AU is ∼4.5 × 10⁻⁹ and at 10 AU is ∼3 × 10⁻¹¹. For an Earth analog (1 R_E) at 1 AU with A_g = 0.3, the contrast is ∼1.2 × 10⁻¹⁰. High contrast is more difficult to achieve closer to the star than further away, however, the d⁻² dependence also means that more distant objects will be fainter in reflected light.

As discussed in Marley et al. (2007), determining atmospheric composition is a central objective for characterizing extrasolar giant planets (EGPs). A census of how composition varies with planet mass and orbital distance will help inform theories of giant planet formation. The giant planets in our solar system motivate us to consider a large range of metallicities in our model exoplanets. As summarized in Fortney et al. (2008b), the Galileo and Cassini spacecraft have shown that the atmospheres of Jupiter and Saturn have enhanced heavy element abundances compared to that of the Sun. For example, Jupiter has ∼2× to 4× the solar abundance of oxygen, carbon, nitrogen, sulfur, and various noble gases (Atreya et al. 2003). Understanding how giant planets formed and the role of processes such as planetesimal bombardment and accumulation during formation (Owen et al. 1999; Gautier et al. 2001a, 2001b; Guillot & Gladman 2000; Alibert et al. 2005), erosion of the heavy-element core (Stevenson 1985; Guillot et al. 2004), direct accretion of metal-rich disk gas (Guillot & Hueso 2006), or chemical fractionation within the planet (Stevenson & Salpeter 1977; Lodders 2004) in enriching gas giants in heavy elements is relevant to understanding both our own solar system and exoplanet systems.

Our goal is to better understand the combined effect of planet–star separation, composition, and planet phase on spectra that are likely to be obtainable in early reflected-light images of exoplanets. For this initial study, we assume simple circular orbits where the planet–star separation is constant and examine how the albedo spectra and colors change with planet phase due to the change in angle from the source to the observer. More complex orbits, such as eccentric orbits, require consideration of how the structure of the exoplanet's atmosphere would evolve as a function of the distance from its parent star. An initial study of the effect of eccentric orbits, using interpolated temperature and pressure profiles corresponding to different planet–star separations, was presented for a 1 M_J EGP in Sudarsky et al. (2005).

As with any planet, exoplanet albedo spectra are very sensitive to the presence and structure of clouds (Marley et al. 1999; Sudarsky et al. 2000). Because the global height and optical thickness of cloud layers are difficult to model on an a priori basis, Marley et al. (1999) considered a few characteristic cloud models in order to test the sensitivity of albedos to clouds. Sudarsky et al. (2000) likewise employed a fairly simple method for computing cloud opacity. In this study, we use the more complex cloud model of Ackerman & Marley (2001) to calculate the height, particle sizes, and consequently optical thickness of water and ammonia clouds. A comparison of the Ackerman & Marley (2001) model with the models used in Marley et al. (1999) can also be found in Ackerman & Marley (2001). The Ackerman & Marley (2001) cloud model provides wavelength-dependent measures of the anisotropic scattering in terms of the asymmetry factor g, and single-scattering albedo $\tilde{\omega }$ of the cloud particles at each temperature and pressure level. In Section 3, we discuss how these parameters are used to approximate anisotropic scattering as well as describe planned modifications to the cloud model to provide additional scattering information.

The Sudarsky et al. (2000, 2005) models generally focused on 1 M_J planets and considered the effects that might be manifested in planets with inclined or highly elliptical orbits (such as transitioning from cloudy to cloud-free and thus bright to dark in reflected light). The effect of geometry on the orbit determination of close-in exoplanets from photometry was also recently presented in detail by Brown (2009). In this work, we consider both cloud-free (0.8 AU) and cloudy (>1 AU) planets. We extend the work of Marley et al. (1999) to model exoplanet observations as a function of phase and Sudarsky et al. (2000, 2005) by considering a greater range of planet masses and metallicities.

We briefly review the basics of imaging exoplanets in the optical and albedo theory in Section 2 as framework for the following results and discussion. While this material is familiar to the solar system planetary science community, it bears reinforcing because of the subtleties between various types of planetary albedos and their relationship to phase. In Section 3, we describe how we generate a series of exoplanet model atmospheres over a range of planet–star separations and metallicities using a one-dimensional radiative–convective model that was previously tailored for use with EGPs (Marley 1997; Fortney et al. 2005, 2006; Marley et al. 2007; Fortney & Marley 2007; Fortney et al. 2008b, 2008a) and brown dwarfs (Marley et al. 1996; Burrows et al. 1997; Marley 1997; Marley et al. 2002; Saumon et al. 2006, 2007). The model was originally used for solar system planetary bodies such as Titan and Uranus (Toon et al. 1977, 1989; McKay et al. 1989; Marley & McKay 1999). These global mean one-dimensional exoplanet models are then used as input to a high-resolution albedo spectra model. We describe our updates to the albedo model that allow us to calculate emergent intensities and thus albedo spectra as a function of phase.

In Section 4, we present results: albedo spectra from 0.35 μm to 1 μm as a function of planet type, planet–star separation, metallicity, and phase. Because of the practical constraints inherent in spacecraft flybys of solar system giant planets, there is relatively little data in the literature reporting planet-averaged phase functions; we do compare Voyager 1 data of Uranus and Neptune from Pollack et al. (1986) with our model phase functions. We also compare the model albedo spectra with observations of solar system giant planets from Karkoschka (1994). Since early exoplanet direct imaging observations will have limited resolution, in Section 5 we consider lower-resolution, coarse spectra derived from our high-resolution results. For the lower resolution cases, we consider R = λ/Δλ = 5 and 15, and we also consider how changes in albedo spectra might present themselves in terms of color–color comparisons using standard filters in the optical. This helps us to evaluate features that might be detectable using a combination of different wide (or narrow) filter bands.

Color–color diagrams were suggested as a comparative analytic tool for direct imaging of exoplanets in an example that considered the color diversity of solar system planets in Traub (2003). Recently, Fortney et al. (2008b) presented detailed color–color comparisons of hot young Jupiters in the infrared. Earlier, Sudarsky et al. (2005) modeled reflected-light albedo spectra and performed color–color analyses in the optical of a 1 M_J EGP using a similar approach but with a different model and implementation than that used in this work. While our 1 M_J results that include clouds tend to be less red at full phase and less blue at new compared with those in Sudarsky et al. (2005), we generally agree with their cloud-free result and extend their approach to include different compositions, Neptune analogs, and comparisons with observed spectra and colors of solar system planets.

To give the reader an idea of the challenges inherent in coronagraphic detection and characterization of gas- or ice-giant exoplanets, Figure 1 shows simulated direct images of Jupiter and Neptune analogs observed in three different 100 nm wide bands around a solar analog at a distance of 10 pc from the observer (Cahoy et al. 2009). The simulations are performed for a 4 m diameter space telescope, where the light from the parent star has been suppressed by a phase-induced amplitude apodization (PIAA) coronagraph (Guyon et al. 2005), and stellar leakage is included as one of several background sources. Importantly, for this figure only, we assume a gray albedo spectrum with a geometric albedo that is a constant 0.3 with wavelength and thus the same for each planet in each bandpass. The simulations are described in greater detail in Cahoy et al. (2009).

Zoom In Zoom Out Reset image size

Figure 1 helps to emphasize the importance of planet–star contrast ratios and angular resolution in a direct imaging observation. It also puts into context several factors that are not specific to a particular instrument or method, but that are important to consider when both interpreting observations and when constructing exoplanet models that correspond to early direct imaging capabilities. For example, the simulated images make clear the relevance of planet size, wavelength-dependent brightness (albedo spectra, one can imagine the effect of a decreasing albedo with wavelength on the images), minimum and maximum detectable separations from the parent star, contributions from background sources like exozodiacal debris disks, and the effect of system geometry and planet phase on the observations. The wavelength-dependent plate scales in this example also give a sense of the level of calibration required to construct coarse albedo spectra even from concurrent measurements. For photon-starved exoplanet observations, the tradeoff between resolution, integration times, and noise is the major challenge in obtaining data in multiple bands.

The term inner working angle (IWA) is often used to describe the minimum detectable separation between a star and planet in a high-contrast direct imaging measurement, where the IWA is generally substantially larger than the diffraction limit (e.g., the 1.22λ/D Rayleigh criterion). For a band-limited coronagraph, the IWA is defined as the separation at which the planet's light is at half of the maximal throughput of the starlight suppression technique at a contrast level sufficient for detection (Cagigal et al. 2009). An instrument that is improved to achieve a smaller IWA can detect exoplanets at closer star–planet separations and thus can also detect exoplanets around stars at greater distances from the observer. The IWA for current observing systems is typically on the order of a few tenths of an arcsecond (Beichman et al. 2010); instruments in development may achieve IWAs on the order of a tenth of an arcsecond. Examples of some current direct imaging methods include high-performance coronagraphs (Kuchner & Traub 2002; Kasdin et al. 2003; Guyon et al. 2005; Trauger & Traub 2007; Mawet et al. 2009, 2010), external occulters (Cash 2006), and optical nullers (Wallace et al. 2000; Martin Levine et al. 2006), as well as post-processing algorithms such as those used by Marois et al. (2006) and Lafreniere et al. (2007).

In practice, for current direct imaging methods the IWA is typically larger than 2λ/D at the required contrast levels for imaging exoplanets in reflected light (the IWA for the PIAA coronagraph simulated in Figure 1 is also ∼2λ/D). Early optical direct images of exoplanets will likely be limited to planet–star separations greater than 2λ/D, since the IWA metric generally does not include external contributing factors such as noise, background sources, coatings, optical throughput, quantum efficiency, etc. There are also platform-specific limitations, for example, primary mirror diameter D is limited on space-based direct imaging missions by either deployment complexity, cost, or launch vehicle fairing size. Ground-based instruments do not have this limitation on D, however, they do have other challenges, such as contrast limitations due to turbulence in Earth's atmosphere. Given the capability of current technology, using a minimum separation of ∼1 AU for our exoplanet models around a Sun-like star seems reasonable for this first set of simulations. In addition, a separation of ∼1 AU is far enough away from a Sun-like star that reflected light will dominate the features in the optical spectrum, even of gas giants (Marley et al. 1999).

The model planets shown in Figure 1 are based on their solar system counterparts and are separated from their parent star by luminosity-scaled equivalents of 1 AU, 2 AU, and 5 AU. We hold the albedos constant versus λ at 0.3 in these simulations to better illustrate the other wavelength-dependent changes. If we were to ignore the wavelength-dependent changes in albedo, the larger-radius gas-giant planets are simply brighter and easier to detect than terrestrial exoplanets. This is particularly true at redder wavelengths where the λ dependence of the IWA "pushes" the planet closer to the IWA, effectively making smaller separations harder to detect at longer wavelengths. As we discuss in Section 4, excluding the effect of photochemical hazes, realistic giant planets get darker at smaller planet–star separations when the temperatures become too warm for bright clouds to form. Depending on the star, the darkening at temperatures that are too high for cloud formation may occur inside of the IWA, making the detection of such planets less likely. For reflected-light direct imaging of cooler exoplanets, Figure 1 shows that it is easier to detect the larger Jupiter and Neptune analogs compared with terrestrial planets, which motivates our use of giant planets as the exoplanet models in this work.

For some direct imaging methods, there may also be limitations on the maximum planet–star separation, or outer working angle (OWA). This can be due to the optical design or elements such as deformable mirrors, detector size, sub-apertures, etc. It is desirable to have an OWA of ∼20λ/D or larger, if possible, so that in addition to exoplanets, exozodical dust and circumstellar debris disk structure can be characterized at much wider separations.

The existence of giant planets at separations larger than ∼10 AU is difficult to account for in standard core accretion models (Pollack et al. 1996; Ida & Lin 2005; Dodson-Robinson et al. 2009). Their presence at very large planet–star separations may be due to a different formation mechanism, such as gravitational fragmentation in the disk (Boss 2000; Kratter et al. 2010) or due to a combination of multiple effects that move planets formed in dense inner regions of a disk out to distant orbits, including outward migration and planet–planet scattering (Veras et al. 2009).

In this work, we model Jupiters and Neptunes out to 10 AU based on the theoretical predictions for giant planet separations and the fact that the contrast requirement for cool giant planets beyond 10 AU is challenging (d⁻² dependence on the incident light, as noted in Section 1.1). Using Equation (1), a 1 R_J planet with geometric albedo of 0.5 has a contrast of ∼1 × 10⁻¹⁰ at 15 AU and ∼6 × 10⁻¹¹ at 20 AU when observed in reflected light at quadrature. Such contrasts are comparable to or fainter than that expected for terrestrial planets at 1 AU, thus motivating our outer model boundary at 10 AU. In Section 3, we discuss the role that planet–star separation plays in controlling the atmospheric structure (clouds) and albedo spectra of our model exoplanets.

Simultaneous images taken with different filters can function as coarse spectra. In Figure 2, we use observed spectra of the outer solar system planets and Titan from Karkoschka (1994) to illustrate how coarse spectra from an early direct imaging observation might be used to differentiate between planet types. The original spectral data from Karkoschka (1994) are binned down to resolution R = λ/Δλ = 5 and R = 15. Even in the R = 5 spectra, one can distinguish between the ice giants (Neptune, Uranus) and the gas giants (Jupiter, Saturn). The ice giants are brightest in the blue and get darker in the red, and the gas giants are not as bright in the blue, due partly to the effects of photochemical products (which we do not model here), but are brighter mid-band and in the red. An important goal of our study is to understand the extent to which such distinctions hold in general for exoplanets and how the phase at which they are observed affects the colors.

Zoom In Zoom Out Reset image size

As we discuss in more detail in Sections 2.2 and 4, coarse spectra can be used in combination with atmospheric models to constrain the temperature, composition, the presence or the absence of clouds in an exoplanet's atmosphere, and the exoplanet's radius (Marley et al. 1999, 2007; Marley & McKay 1999; Fortney et al. 2008b, 2008a). Color–color comparisons of exoplanets can be made (Traub 2003; Sudarsky et al. 2005; Fortney et al. 2008b), and understanding how system geometry and exoplanet composition relate to variability in the observations is important for interpreting these color–color comparisons.

The definition and relationships between geometric albedo A_g, phase α, the phase integral q, spherical albedo A_s, and Bond albedo A_B, and additional parameters of interest such as planetary effective temperature are important to understand in the context of interpreting direct imaging observations. Although this material is not new, we review it since the distinctions between the various types of albedos are important to emphasize. Note that while Marley et al. (1999) presented geometric albedo spectra and tables of Bond albedos for exoplanets in systems with varying stellar types, our focus here is more closely on the combined effects of separation, composition, and planet phase on model gas and ice giant spectra. We do tabulate Bond albedos for our exoplanet models, but only for a solar analog parent star.

Measurements of the planet's brightness at a given wavelength and full phase yield the planet's monochromatic geometric albedo:

$$\begin{equation} A_g(\lambda) = \frac{F_{p}(\lambda, \alpha =0^{\circ })}{F_{\odot,L}(\lambda)}. \end{equation} \tag{ 2 }$$

Geometric albedo A_g(λ) is the ratio of the reflected flux F_p(λ, α = 0°) of an object at full phase (α = 0°) to the flux from a perfect Lambert disk, F_☉,L of the same radius R_p under the same incident flux F_☉ at the same distance d from the star. A planet's phase angle α is the angle between the incident ray from the star to the planet and the reflected ray from the planet to the observer. The scattering angle Θ is related to the phase angle via α = π − Θ.

In this work, we use the term geometric albedo to refer to albedo spectra at full phase. We use the term albedo spectra to refer to spectra observed at different orbital phases.

In order to reduce confusion in later discussions of planet phase functions and particle scattering functions, we use the notation A_g(λ) for geometric albedo here instead of using the other common notation for geometric albedo, p. As noted in Sobolev (1975), for a perfectly reflecting Lambert sphere, the geometric albedo is 2/3, and for a semi-infinite purely Rayleigh scattering atmosphere, it is 3/4 (Dlugach & Yanovitskij 1974).

Flux measurements at different planet phase angles can be used to define the phase function Φ(λ, α):

$$\begin{equation} \frac{F_p(\lambda,\alpha)}{F_\odot (\lambda)} = A_{g}(\lambda)\left(\frac{R_p}{d}\right)^2 \Phi (\lambda,\alpha). \end{equation} \tag{ 3 }$$

As noted for Equation (2), A_g(λ) is defined at α = 0°, F_p(λ, α) is the monochromatic planet flux, F_☉(λ) is the monochromatic stellar flux, R_p is the planet's radius, d is the planet–star separation, and Φ(λ, α) is the planet's phase function at angle α. The phase function is normalized to be 1.0 at full phase. The scattering phase function of the particles in the exoplanet's atmosphere, p(Θ), contributes to the disk-integrated phase function. The observed phase function allows calculation of the phase integral, q:

$$\begin{equation} q(\lambda) = 2 \int _0^{\pi } \Phi (\lambda,\alpha)\sin \alpha d\alpha. \end{equation} \tag{ 4 }$$

The spherical albedo,

$$\begin{equation} A_s(\lambda) =q(\lambda)A_g(\lambda) \end{equation} \tag{ 5 }$$

is the fraction of incident light reflected toward all angles.

The planet's Bond albedo, an incident flux-weighted and wavelength-integrated function of A_s(λ) and Φ(λ, α), normalized by the incident stellar flux, can then be computed:

$$\begin{equation} A_B = \frac{\int _0^{\infty }A_{s}(\lambda)F_{\odot }(\lambda)d\lambda }{\int _0^{\infty }F_{\odot }(\lambda) d\lambda }. \end{equation} \tag{ 6 }$$

It is important to appreciate that since the Bond albedo is weighted by the incident flux, two identical planets with identical geometric albedo spectra will have different Bond albedos if the incident flux is different (e.g., comparing an F star and an M star; Marley et al. 1999). The Bond albedo can be used to find the effective temperature, T_eff of a rapidly rotating planet:

$$\begin{equation} 4\pi R_p^2 \sigma T_{{\rm eff}}^4 = E_{{\rm int}} + \pi \it R_p^2[f(1- A_B)]\frac{\pi \it S}{d^2}, \end{equation} \tag{ 7 }$$

where σ is the Stefan–Boltzmann constant, πS is the incident flux of the parent star at a planet–star separation of 1 AU, d is the planet–star separation, here in AU, and E_int represents the internal energy sources of the planet (e.g., Conrath et al. 1989).

An additional redistribution parameter f may also be introduced to parameterize the global redistribution of energy to solve for a day and nightside temperature. For solar system giant planets, typically f = 1, since Jupiter or Neptune-like planets with rotation periods of hours and effective temperatures below 500 K have sufficiently long radiative time constants such that atmospheric dynamics are likely to redistribute incident flux around the planet. For exoplanets, sometimes different values of f are used due to the diversity of exoplanet system geometries. For example, f = 2 represents the case where the redistribution of the absorbed flux is over the day side only, yielding a higher dayside effective temperature. Here, we use f = 1.

In Section 3 and the Appendix, we describe updates to our albedo spectra model that allow us to generate albedo spectra at different planet phases for a range of input model gas- and ice-giant exoplanet atmospheres, which can then be used to calculate exoplanet colors, phase integrals, spherical albedos, and Bond albedos.

We generate the exoplanet model atmospheres with the one-dimensional iterative radiative–convective atmosphere model previously used in Marley (1997), Burrows et al. (1997), Marley et al. (1999, 2002), Fortney et al. (2005, 2006, 2007), Fortney & Marley (2007), and Fortney et al. (2008b, 2008a). The radiative transfer method was developed by Toon et al. (1977, 1989), implemented by McKay et al. (1989) and Marley & McKay (1999), and used in application to Titan (McKay et al. 1989) and Uranus (Marley & McKay 1999) as well as to EGPs and brown dwarfs, as noted in Section 1. The model computes the deposition of incident flux from $0.268\,\rm \mu m$ to $6.064\,\rm \mu m$ and accounts for thermal emission out to $227\,\rm \mu m$.

In this work, we use elemental abundances from Lodders (2003). We employ tables of chemical equilibrium compositions following Fegley & Lodders (1994), Lodders & Fegley (2002, 2006), and Freedman et al. (2008). The chemistry calculations include "rainout," where refractory species are depleted from the atmosphere due to their condensation into cloud decks (Lodders 1999). The spectra of both brown dwarfs and our solar system's giant planets can only be reproduced when chemistry calculations incorporate this process (Fegley & Lodders 1994; Marley et al. 2002; Burrows et al. 2002).

As noted in Section 1.1, the four gas giants in our solar system show an enhancement of heavy elements compared to solar abundances. The amount of enhancement provides additional clues about planet formation and evolution. Our models include a range of compositions, with metallicities of 1× and 3× solar for the model Jupiters, and 10× and 30× solar for the model Neptunes. To be specific, we use metallicities of [M/H]/[M/H_cosmic] = 0.0, 0.5, 1.0, and 1.5. So although we refer to them as 3× and 30×, the 3× and 30× metallicities are actually 3.16× and 31.6×.

The radiative–convective model uses the correlated-k method for the tabulation of gaseous opacities (Goody et al. 1989). The extensive opacity database we use is fully described in Freedman et al. (2008). We note that we have yet to incorporate the recent tabulation of methane opacities presented in Karkoschka & Tomasko (2010). We also note the general difficulties in determining the accuracy of opacities at the temperatures and pressures of interest, and of reconciling laboratory observations with Earth-based, space-based, and in situ observations.

We use the cloud model of Ackerman & Marley (2001) to describe the location, vertical distribution, and particle sizes of the major cloud-forming species, which are H₂O and NH₃ in this work. Clouds are coupled iteratively into the radiative–convective model using the approach of Ackerman & Marley (2001) as described in Section 1.3. The model computes the vertical variation of particle sizes and number densities given an assumed cloud sedimentation efficiency, f_sed which we discuss further below. Since the cloud model is converged along with the atmospheric thermal structure, the clouds are fully self-consistent with the atmospheric thermal structure of each model planet. Note that we neglect $\rm NH_4SH$, which may produce an intermediate cloud between the water and ammonia clouds in Jupiter's atmosphere (West et al. 2004). Even at Jupiter the relative roles that NH₃ and NH₄SH clouds play in influencing the albedo spectra remain uncertain (Irwin et al. 2005) and thus we treat only the former cloud.

For each model planet, our one-dimensional modeling thus produces a vertical temperature and pressure, T(P), profile, a cloud structure model, and the variation in composition of all major species. We use these profiles and cloud information (wavelength-dependent asymmetry factor g and single-scattering albedo $\tilde{\omega }$ at each temperature and pressure level) as input to the albedo spectral model described in Section 3.2.

3.1.1. Model Temperature and Pressure Profiles

Table 1 summarizes the exoplanet model gas giants (Jupiters) and ice giants (Neptunes) used in this work. The internal heat flows assumed for the Jupiters and Neptunes correspond to effective temperatures (T_int) of 100 K and 50 K, respectively, for an isolated object, and an age of ∼4.5 Gyr. Gravity for Jupiters is 25 ms⁻² and for Neptunes is 10 ms⁻². We compute radiative–convective equilibrium T(P) profiles and low-resolution spectra for exoplanet models at planet–star separations of 0.8 AU, 2 AU, 5 AU, and 10 AU at 1× and 3× solar metallicity for the Jupiters and 10× and 30× solar metallicity for the Neptunes. As noted in Section 2.1, the maximum and minimum planet–star separations were chosen to span the range of exoplanets detectable with early coronagraphic direct imaging methods.

Table 1. Description of the Model Gas and Ice Giants Used as Inputs to the Albedo Spectra Simulation

Separation	Jupiters			Neptunes			Clouds	fsed
	g = 25 m s−2			g = 10 m s−2
	Tint (K)	Teff (K)		Tint (K)	Teff (K)
		1×	3×		10×	30×
0.8 AU	100	274	282	50	285	292	None	...
2 AU	100	123	126	50	134	148	H2O	6
5 AU	100	115	119	50	110	112	H2O, NH3	10
10 AU	100	107	109	50	81	82	H2O, NH3	10

Notes. The separations and cloud types are noted. T_int is a parameterization of the flux from the interior of the planet, due to its cooling with time (temperature in the absence of any incident flux). f_sed is the ratio of the microphysical sedimentation flux to the eddy sedimentation flux.

Download table as:  ASCIITypeset image

Figure 3 shows the resulting T(P) profiles for the models. Solid lines are used for the lower metallicity 1× and 10× cases, and dashed lines are used for the higher metallicity 3× and 30× cases. Condensation curves, courtesy K. Lodders, are shown on each plot for both H₂O (warmer) and NH₃ (cooler). The condensation curves correspond to 1× (solid line) and 3× (dashed line) solar metallicity for the Jupiter models and 10× (solid line) and 30× (dashed line) solar metallicity for the Neptune models.

Zoom In Zoom Out Reset image size

To illustrate the impact of clouds on the albedo spectra and thus the detectability of exoplanets, we chose planet models across the full range of planet–star separations likely to be probed by early direct imaging observations. Our planet–star separations include cloud-free, H₂O, and combined H₂O and NH₃ cloud cases for both the Jupiters and Neptunes. However, for the purposes of comparison and simplicity, we also used the same planet–star separations for both the Jupiters and the Neptunes to limit the number of models to a few representative cases.

While the Ackerman & Marley (2001) approach results in clouds with varying particle sizes as a function of height through the atmosphere, for heuristic purposes the mean or "effective" cloud particle radii can also be computed (see Ackerman & Marley 2001). Table 2 lists the effective radii of the condensate particles for the clouds in each of the exoplanet models. We stress that these are provided here for comparison between the models; the actual radiative transfer calculations do not employ these mean particle sizes, but rather the actual computed particle size profiles. We used the same values for f_sed, the ratio of the microphysical sedimentation flux to the eddy sedimentation flux, for both Jupiter and Neptune models at a given planet–star separation. At 2 AU, we used f_sed = 6, and at 5 AU and 10 AU, we used f_sed = 10. These values of f_sed were selected partly from our experience with fitting L dwarf spectra (Stephens et al. 2009) and partly to allow convergence (very thick clouds modeled with low values of f_sed do not converge well). Ackerman & Marley (2001) found good agreement with the NH₃ cloud in Jupiter's atmosphere using f_sed = 3. Larger values of f_sed lead to clouds of slightly smaller vertical extent and optical depth (Ackerman & Marley 2001). In a future study, we will explore the effect of changing cloud thickness on planet colors. For the current study, all of the clouds have large optical depths, so the impact of varying f_sed is small. We stress that for this study we are primarily interested in the macro effect of clouds and not in detailed changes to the cloud model. A higher value of f_sed was required for the water cloud planets in order to produce properly converged models.

Table 2. Characteristic "Effective Radii" (μm) of Cloud Condensate Particles Predicted by the Ackerman & Marley (2001) Model

Separation	Jupiters				Neptunes
	1 ×		3 ×		1 ×		3 ×
	H2O	NH3	H2O	NH3	H2O	NH3	H2O	NH3
0.8 AU		...	...	...	...	...	...	...
2 AU	69	...	72	...	52	...	53	...
5 AU	102	56	103	61	72	61	73	67
10 AU	100	74	102	87	80	71	80	73

Note. The parameter "effective radius" is used for comparative purposes to summarize the level-dependent distribution of radii actually used in the model.

Download table as:  ASCIITypeset image

Figure 3 shows that at 0.8 AU, both the 1× and 3× Jupiters are just warm enough to be cloud-free. At 2 AU the Jupiters are just warm enough to only condense H₂O clouds, and at 5 AU and 10 AU the Jupiters have both H₂O and NH₃ clouds. Figure 3 also shows that the H₂O and NH₃ condensation curves for the 10× and 30× Neptunes are a bit warmer than their 1× and 3× counterparts. The 0.8 AU Neptunes are also just warm enough to be cloud-free. The 2 AU Neptunes may be just cool enough to condense both H₂O and NH₃. However, in this work we neglect the condensation of NH₃ at 2 AU, for better comparison with the 2 AU Jupiter models that only condense H₂O. The 5 AU and 10 AU Neptune models have both H₂O and NH₃ clouds.

Our 1× Jupiter models appear to be a bit cooler than those shown in Sudarsky et al. (2005). This is probably due to our computing a global mean temperature profile which assumes that the absorbed incident flux is redistributed by atmospheric dynamics evenly around the planet. Sudarsky et al. (2005) instead compute a mean dayside temperature profile without energy redistribution. While the efficiency of energy redistribution is an interesting topic for the very hot Jupiters, at the lower effective temperatures relevant to this study the radiative time constants are long and the incident flux is likely to be evenly distributed, as it is on the solar system Jovian planets, which do not exhibit any significant diurnal thermal differences.

Given the atmosphere models generated by the radiative–convective model we next apply the albedo spectra model. This albedo spectra model has previously been used for a variety of planetary objects and is described in detail in Toon et al. (1977, 1989), McKay et al. (1989), Marley et al. (1999), and Marley & McKay (1999). In this section, we describe the configuration of the albedo spectra model used here. The albedo model takes as input the exoplanet model's gravity and atmospheric temperature, pressure, composition, and cloud properties that are generated within the radiative–convective model discussed in Section 3.1.

The albedo spectra model uses the same chemistry and opacity database as were used for obtaining the thermal structure. In the thermal structure model, opacities were binned into intervals using k-coefficients. For the albedo calculation we use monochromatic opacities over the ∼790 points from 0.35 μm to 1 μm. The opacity varies with altitude through the model and we account for the changing gaseous composition as well as the influence of temperature and pressure on the opacities themselves. In this work, we do not model photochemistry or attempt to account for photochemical hazes or chromophores (see West et al. 2004 for a review), which are responsible for the low albedos for the real Jupiter in the blue and UV (<0.4 μm) and also for its ruddy complexion; this is discussed in more detail in Marley et al. (1999). While we expect that photochemistry will play an important role in exoplanet spectra, we reserve such considerations for future work.

3.2.1. Radiative Transfer Description

Figure 4 shows a schematic of the plane–parallel model geometry with incident angle μ₀ = cos β = sin ηcos(ζ − α) and observed angle μ₁ = cos ϑ = sin ηcos ζ, where ζ and η are planetary coordinates of longitude and co-latitude, respectively.

Zoom In Zoom Out Reset image size

We use a finite plane–parallel atmosphere divided into N levels and N−1 layers. The plane–parallel framework is shown in Figure 4. Emergent intensity is I(τ, μ₁) with optical depth τ_i from level i to the top of the atmosphere and incremental change in optical depth of $\delta _{\tau _i} = \tau _{i+1}-\tau _i$ between levels. Both incident and emergent angles are positive from the surface normal, with incident light from direction μ₀ = cos β and the observer direction μ₁ = cos ϑ. A negative sign is used to denote the direction of incident flux, e.g., −μ₀. With a source at μ₀, the exact expression for the emergent intensities through the N levels is

$$\begin{eqnarray} I(\tau _i, \mu _1, \phi _1) &=& I(\tau _{i+1}, \mu _1, \phi _1)e^{-\delta _{\tau _i}/\mu _1}\nonumber\\ && + \int _{0}^{\delta _{\tau _i}} S(\tau ^{\prime }, \mu _1, \phi _1)e^{-\tau ^{\prime }/\mu _1}\frac{d\tau ^{\prime }}{\mu _1}, \end{eqnarray} \tag{ 8 }$$

where S(τ', μ₁, ϕ₁) is the source function and contains contributions from both direct and diffuse scattering. The incident collimated flux crosses the unit area perpendicular to direction defined by −μ₀ and $\tilde{\omega }$ is the "single-scattering" albedo, a parameter used to represent the amount of scattering taking place versus absorption. $\tilde{\omega }$ is calculated for each layer and at each wavelength by taking the ratio of the scattering coefficients to the total scattering and absorption coefficients in a layer. The asymmetry parameter, g, is similarly calculated for each layer and wavelength using asymmetry contributions and opacities of the layer constituents. The scattering and asymmetry contributions include an approximation to Raman scattering, Rayleigh scattering, and clouds in this model.

While we begin with an exact solution to the radiative transfer in Equation (8), we introduce an approximation by substituting a two-stream solution for the scattered radiation field (intensity) in the source function. In the case of negligible diffuse radiation (when single-scattering dominates, as it tends to for clear atmospheres) this approximation has no effect. As the importance of multiple scattering increases (e.g., when there are clouds and hazes), the influence of this approximation grows. We present our two stream solution in some detail, but stress that the albedo calculation as used here is not simply a "two stream solution" (Toon et al. 1989).

The radiative transfer equation is often solved using approaches that do not preserve the angular dependence of the scattering phase function, such as the two-stream methods summarized in Meador & Weaver (1980) or the Feautrier method in Mihilas (1978). Here, we follow the "source function" approach in Toon et al. (1989) where we first use two-stream quadrature to solve for the diffuse scattered radiation field. We use the resulting two-stream intensity I_t as an approximation to the source function (S_v versus S_s in Equations (1) and (52) of Toon et al. 1989). We are then able to approximate the angular dependence (μ, ϕ) and (μ', ϕ') in the scattering phase function. We use the two-stream "source function" approach here largely because of its compatibility with the existing model. In a future paper, we will compare the approach in this work with other radiative transfer methods that preserve the angular dependence of the scattering phase function, such as the δ−M stream method (Wiscombe 1977). Comparisons with models that include the vector nature of polarized scattering (e.g., Stam et al. 2006; Buenzli & Schmid 2009) instead of the scalar approach taken for speed and simplicity in this work would also be interesting.

Following the two-stream source function method (Toon et al. 1989), we use two-stream quadrature (sometimes called discrete ordinates) to solve the radiative transfer equation for the upward and downward diffuse fluxes. The two-stream quadrature approach is only one of many different two-stream methods to approximate the diffuse fluxes. Other approaches, including the hemispheric constant, Eddington, modified Eddington, modified quadrature, and delta-function methods, are treated in detail in sources such as Meador & Weaver (1980), Liou (1980), and Toon et al. (1989). In the form of the "unified description" of two-stream methods in Meador & Weaver (1980), the Gaussian quadrature parameters γ_i used in this work are

$$\begin{equation} \gamma _1 = \frac{\sqrt{3}}{2}[2-\tilde{\omega }(1+g)], \end{equation} \tag{ 9 }$$

$$\begin{equation} \gamma _2 = \frac{\sqrt{3\tilde{\omega }}}{2}(1-g), \end{equation} \tag{ 10 }$$

$$\begin{equation} \gamma _3 = \frac{1}{2}(1-\sqrt{3}g\mu _0), \end{equation} \tag{ 11 }$$

$$\begin{equation} \gamma _4 = 1-\gamma _3 = \frac{1}{2}(1+\sqrt{3}g\mu _0). \end{equation} \tag{ 12 }$$

The forward-scattering nature of many atmospheric particulates, especially clouds, is not well captured by the two-stream approach as described above. The delta-function technique takes advantage of the similarity principle to split out a fractional forward-scattering contribution f from an original scattering phase function into a delta function, and adjusts the asymmetry factor g, single-scattering albedo $\tilde{\omega }$, and optical depth τ accordingly. This concept is treated in detail in Joseph et al. (1976) who use it along with an Eddington two-stream approach. We use a delta-function approach as described in Liou (1980) along with the two-stream source function method in this work as well.

3.2.2. Boundary Conditions

We use the recursive equation (Equation (8)) to calculate the level intensities from the "ground," N to the top of the atmosphere. Note that in discussion of boundary conditions it is common to refer to the "ground" or "surface." For gas-giant exoplanets, these terms are not accurate since there is no "ground" in the sense that there is one for a terrestrial planet. Instead, we treat the "surface" as a pressure level deep enough that no photons are reflected back out. We set the following boundary condition for the ground level:

$$\begin{equation} I(\tau _N,\mu _1) = R(\mu _1)\frac{F^{\downarrow {\rm total}}(\tau _N)}{\pi }, \end{equation} \tag{ 13 }$$

where F^↓total(τ_N) = F^↓(τ_N) + F^↓direct(τ_N) is the total downward flux onto the surface, both diffuse and direct. We treat the ground layer as a Lambert surface with reflectance R(μ₁) = R. In this work, we use a "black" surface, R = 0. However, in the general case where R ≠ 0, we would need to solve the diffuse radiative transfer equation to get F^↓(τ_N).

We also set boundary conditions at the top of the atmosphere: there is no diffuse flux coming from outside the atmosphere, and any upward flux consists of the reflected direct flux off the Lambert surface plus the downward diffuse flux scattered upward by the Lambert surface:

$$\begin{equation} F_{{\rm top}}^{\downarrow } = 0 \mbox{ and } F_{{\rm Surface}}^{\uparrow } = R\big(F_{{\rm Surface}}^{\downarrow {\rm direct}} + F_{{\rm Surface}}^{\downarrow }\big). \end{equation} \tag{ 14 }$$

These are the boundary conditions used with the two-stream source function method described above.

3.2.3. Scattering

The albedo spectra model uses a wavelength-dependent approximation of anisotropic scattering. The cloud model, which uses a full Mie scattering model (Ackerman & Marley 2001), currently returns the single-scattering albedo $\tilde{\omega }$ and asymmetry parameter g as a function of wavelength for each condensate species and at each temperature and pressure level. In future work, we will update the cloud model to return the angular distribution of the scattering in addition to the asymmetry parameter. With this type of approach we would then use an M-term fit to the full angular distribution of the scattering function returned by the cloud model in conjunction with a δ−M-stream radiative transfer approach. For the moment, we start with the available $\tilde{\omega }$ and g (which are functions of wavelength, temperature, and pressure) as input parameters to a two-term Henyey–Greenstein function (TTHG) as an approximation to the anisotropic scattering (as discussed below, the TTHG function approximates the angular dependence of both forward and backward scattering). As with many previous studies, for our current goal of exploring a wide range of exoplanet atmosphere models that may be observed by generally photon-starved direct imaging observations, our treatment of atmospheric scattering is of a scalar, not vector, nature, meaning that the effects of polarization are not considered. We refer readers interested in the effects of polarization to simulations of Jupiter-like and Rayleigh scattering atmospheres in Stam et al. (2004, 2006), Stam & Hovenier (2005), and Buenzli & Schmid (2009).

We also approximate Raman and Rayleigh scattering. The approximate treatment of Raman scattering is the same as described in Pollack et al. (1986), incorporating Raman scattering into the single-scattering albedo as a function of the stellar spectrum (see Equation (2) in Pollack et al. 1986). We approximate Rayleigh scattering by including an additional term in the scattering phase function as briefly described below.

The scattering phase function p(cos Θ) is a representative function that takes contributions from a number of different scatterers in the atmosphere into account. The scattering phase functions used in this model attempt to represent the cumulative scattering effect of all of the constituents in each layer of the model. As noted earlier, the single-scattering albedo $\tilde{\omega }$ and asymmetry factor g used in the scattering phase function are updated by the cloud model and radiative–convective model at each layer and for each wavelength with weighted contributions from the atmospheric constituents.

The scattering angle Θ can also be written in terms of radiation incoming from (μ', ϕ') that is scattered into (μ, ϕ):

$$\begin{equation} \cos \Theta = \mu \mu ^{\prime } + (1-\mu ^2)^{1/2}(1-\mu ^{\prime 2})^{1/2}\cos (\phi ^{\prime }-\phi). \end{equation} \tag{ 15 }$$

The second term on the right-hand side is sometimes neglected in the case of azimuthal independence or averaging. The scattering phase function may be written in terms of cos Θ or μ', ϕ' and μ, ϕ. The integral of the scattering phase function is typically normalized either to unity or to the single-scattering albedo, $\tilde{\omega }$. Note that the scattering phase function p should not be confused with planet phase function Φ(λ, α). Also note the π difference in directionality between the scattering angle and the planet's phase angle, α = π − Θ; backscattering would be more strongly seen by an observer at full phase, α = 0°, whereas backscattered Θ = 180°.

As shown in Figure 7, for these models, the contribution of scattering to the albedo spectra not only depends on the mixture of gases in a given layer, but also the presence, composition, location, and depth of clouds in the atmosphere. There are many different analytic scattering phase functions that attempt to realistically capture the scattering behavior of atmospheric constituents, such as the Henyey–Greenstein function, as well as approaches to capturing strong forward scattering using delta functions. While there are a variety of scattering phase functions that are used to represent different types of atmospheric constituents, it is common to expand the phase function in terms of Legendre polynomials in cos Θ:

$$\begin{equation} p(\cos \Theta) = \sum _{\ell = 0}^{N}\tilde{\omega _{\ell }}P_{\ell }(\cos \Theta), \end{equation} \tag{ 16 }$$

$$\begin{equation} \tilde{\omega _{\ell }} = \frac{2\ell +1}{2} \int _{-1}^{1} p(\cos \Theta)P_{\ell }(\cos \Theta). \end{equation} \tag{ 17 }$$

Note that despite the similar notation, the single-scattering albedo, $\tilde{\omega }$, is different from the Legendre constants $\tilde{\omega }_{\ell }$. In most two-stream models, the Legendre polynomial expansion of the scattering phase function is not used beyond the first order, P₁. In that case, P₀(cos Θ) = 1 and P₁(cos Θ) = cos Θ. For a second-order expansion, $P_2(\cos \Theta) = \frac{1}{2}(3\cos ^2\Theta -1)$.

For a first-order expansion, azimuthal independence is often assumed, leading to

$$\begin{equation} p(\mu,\mu ^{\prime }) = 1 + 3g\mu \mu ^{\prime }, \end{equation} \tag{ 18 }$$

where g is the asymmetry factor, defined as the first moment of the phase function, and also equal to $\tilde{\omega }_{1}/3$. The asymmetry factor g represents the relative strength of forward scattering.

The Rayleigh scattering phase function is

$$\begin{equation} p(\cos \Theta) = \frac{3}{4}(1+\cos ^2\Theta). \end{equation} \tag{ 19 }$$

The cos²Θ dependence is not included in the first-order expansion in Equation (18). However, continuing the Legendre expansion to second order allows representation of Rayleigh scattering with the cos²Θ term:

$$\begin{equation} p(\cos \Theta) = \tilde{\omega }_0 + \tilde{\omega }_1\cos \Theta + \frac{\tilde{\omega }_2}{2}(3\cos ^2\Theta -1). \end{equation} \tag{ 20 }$$

For both isotropic scattering and Rayleigh scattering, g = 0. We use the original first-order expansion (Equation (18)) for p in the initial radiative transfer calculation of the two-stream fluxes that are to be used in the source function method; the truncation treats the Rayleigh scattering as isotropic (since for Rayleigh, g = 0). However, once we have used the source-function method to plug the fluxes back in, we use a second-order phase function expansion for p and assume azimuthal symmetry such that

$$\begin{equation} p(\mu \mu ^{\prime }) = 1 + 3g\mu \mu ^{\prime } + \frac{g_2}{2}[3\mu ^2(\mu ^{\prime })^2-1]. \end{equation} \tag{ 21 }$$

With this expansion, Rayleigh scattering can be modeled when g = 0 and $g_2 = \frac{1}{2}$. In the model, g goes to zero when the contribution from τ_Rayleigh in the denominator is large, and g₂ goes to 1/2. For the purpose of evaluating the second-order expansion with the two-stream source function, we also define a "mean square angle" $\bar{\mu }_2$:

$$\begin{equation} \bar{\mu }_2^2 = \frac{\int \mu ^2 I d\mu }{\int I d\mu }, \end{equation} \tag{ 22 }$$

where $\bar{\mu }_2$ is treated as a free parameter ranging from 0 to 1 that depends on the two-stream angle approximation used, and that, in combination with the higher order term g₂, is intended to provide an accurate result for the pure Rayleigh case. In the model used, $\bar{\mu }_2 \approx (1/\sqrt{3})^{1/2}$. The weighted asymmetry parameters (g, g₂) in the model balance the relative contribution of Rayleigh scatterers such that the diffuse scattering phase function approaches a Rayleigh scattering phase function when Rayleigh scattering dominates. This capability becomes more important at short wavelengths and in clear (cloud-and haze-free) atmospheres, and is particularly relevant for the direct imaging of exoplanets in the optical.

3.2.4. Two-term Henyey–Greenstein

We use a TTHG as the direct scattering phase function to capture the contributions of very strong forward scattering and moderate backscattering. The TTHG makes use of the asymmetry parameter g that is calculated in the radiative–convective model for the clouds that form within each layer following Ackerman & Marley (2001). As noted above, we plan to update this approach to retain the angular properties from each layer and use an M-term Legendre polynomial to best fit the angular distribution of the scattering phase function for each layer. For the moment, we use a TTHG that is a function of g for each layer and that has a form generally consistent with high forward scattering and moderate backscattering by clouds and aerosols, as shown in Figure 5 and informed by studies of the scattering function in solar system giant planets in Pollack et al. (1986) and Dyudina et al. (2005) and discussion of the scattering properties of clouds and aerosols in Liou (1980). The values of g from the cloud model used in this work, at pressure levels where clouds have formed, typically range between 0.5 at short wavelengths and 0.99 at long wavelengths.

Zoom In Zoom Out Reset image size

The general form of a single Henyey–Greenstein function is

$$\begin{equation} p_{\rm HG} = \frac{1-g^2}{(1 + g^2 -2g\cos \Theta)^{3/2}}. \end{equation} \tag{ 23 }$$

The HG phase function is useful because it both approximates strongly peaked forward scattering and has a simple expansion in Legendre polynomials:

$$\begin{equation} p_{\rm HG} = \sum _{\ell = 0}^{N}(2\ell + 1)g^{\ell }P_{\ell }(\cos \Theta). \end{equation} \tag{ 24 }$$

From its expansion, p_HG has a first moment of g and has a second moment f = g². The ability to expand the HG function in terms of N Legendre polynomials as well as the simple form of its second moment also makes it useful for adjusting two-stream solutions to better include the strong forward-scattering nature of cloud particles and aerosols (Joseph et al. 1976). However, a single-term HG scattering phase function does not capture backscattering. To capture backscattering, TTHG scattering phase functions can be used

$$\begin{equation} p_{\rm TTHG} = b p_{\rm HG}(g_a,\Theta) + (1-b)p_{\rm HG}(g_b,\Theta), \end{equation} \tag{ 25 }$$

where g_a>0 causes the forward peaks and g_b < 0 the backward peaks. The form used here has g_a = g, g_b = −g/2, and b = 1 − g²_b.

3.2.5. Planet Phase Geometry

To compute the albedo spectra, we integrate over emergent intensities resulting from the incident radiation. Previously (McKay et al. 1989; Marley et al. 1999), the emergent intensity calculations assumed an α = 0° system geometry where the source and observer were collinear, such that μ₀ = μ₁. With that restriction, to consider variation with phase, the geometric albedo spectra could then be multiplied by some planet phase function Φ(λ, α) to emulate the change in emergent flux with phase. Such an approach does not capture the changes in the albedo spectra with wavelength that occur for different phases as the paths through the atmosphere and contributions from anisotropic scatterers change. We thus have updated the albedo spectra model to use different μ₀ and μ₁.

To simulate a spherical planet, we cover the illuminated surface of a sphere with many plane–parallel facets, where each facet has different incident and observed angles, as shown in Figure 6. Following the approach of Horak (1950) and Horak & Little (1965), we use two-dimensional planetary coordinates and Chebyshev–Gauss integration to integrate over the emergent intensities and calculate the albedo spectra. The Chebyshev–Gauss angles are the input angles to the intensity calculations. The plane–parallel approach currently does not accurately capture curved geometry (limb effects) at high phase angles, e.g., Chapman (1931a, 1931b); we will include this refinement in a future update.

Zoom In Zoom Out Reset image size

In Figures 7 and 8, we show how the optical depth τ varies throughout the model atmospheres with increasing planet–star separation for our 3× Jupiter and 10× Neptune models. The corresponding albedo spectra are shown in Figure 9. Figures 7 and 8 show pressure P versus wavelength with shaded contours of log τ that correspond to the color scale on the right. Optical depth of τ = 1 is at log τ = 0.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the cloud-free 3× Jupiter case at 0.8 AU, Figure 7 shows that Rayleigh scattering dominates the short wavelengths deep into the atmosphere. Contributions from the alkalis, Na (∼0.59 μm) and K (∼0.78 μm, see Table 3), shape the spectrum and are evident deeper in the atmosphere where they are favored by chemical equilibrium. Thus even though the effective temperature of this model (Table 1) is 290 K, well below the temperature at which we would expect the atmosphere to be devoid of alkali elements, alkalis remain an important opacity source. The reason for this is that the clear atmosphere allows photons to penetrate as deep as 10 bars before they are scattered for the 0.8 AU 3× case (Figure 7, top left panel). At 10 bars the temperature exceeds 800 K and the alkalis are present in the gas and thus are detectable.

At 2 AU, bright and optically thick water clouds form below ∼500 mbar; their presence dominates opacity across the optical. Mie scattering predicts that the extinction efficiency at optical wavelengths will be essentially constant for all particles with radii greater than $1\,\rm \mu m$ and thus the spectral contribution of the clouds appears flat in the figures. It is interesting to note that optical thicknesses can reach ∼10, 000 for H₂O clouds. At 5 AU, ammonia clouds form below ∼50 mbar and water clouds form below ∼1 bar, at a deeper level than that for the warmer 2 AU case. The ammonia clouds also contribute significantly to the optical depth. For the still colder 10 AU Jupiter 3× model, both the NH₃ and H₂O clouds form deeper into the atmosphere, allowing Rayleigh scattering to again dominate the opacities at short wavelengths. Even at long wavelengths, unit optical depth is likely to be reached before encountering the clouds.

In addition to the alkalis (clearly seen at 0.8 AU) and CH₄, opacity due to gaseous H₂O plays a role between ∼0.92 and 0.95 μm in our albedo spectra for all of the cases closer in than 5 AU. Detection of water absorption features in the optical has not been confirmed in Jupiter's reflection spectrum, although absorption features near ∼0.94 μm were noted by Karkoschka (1994) as being present. At the time, H₂O was suggested potentially being the cause, however, later observations suggested that NH₃ could be responsible (Karkoschka 1998). We address the presence of H₂O again later in discussion of Figure 11.

Figure 8 shows the relationship between τ, P, λ and planet–star separation for the 10× Neptune models. Similar to the 0.8 AU 3× Jupiter case, the 0.8 AU 10× Neptune case is dominated by Rayleigh scattering at short wavelengths, only with more pronounced absorption features into the red. The 2 AU case is also dominated by the presence of H₂O clouds that form just below ∼380 mbar, a bit higher than the Jupiter case, implying that the albedo spectra should be a bit brighter as well. At 5 AU, an ammonia cloud forms below ∼100 mbar, and extends down nearly to the H₂O cloud, which forms below ∼470 mbar. Also similar to the 10 AU 3× Jupiter case, the 10 AU 10× Neptune shows both cloud decks forming at substantially higher pressures, below ∼700 mbar for NH₃ and below ∼3.6 bars for H₂O.

In Figure 9, we show albedo spectra for all of the model cases summarized in Table 1. The spectra for our Jupiter-like exoplanet models can be compared with those presented in Marley et al. (1999) and Sudarsky et al. (2000, 2005). For the purpose of comparison, the model Jupiters used here map to the Class III (clear), Class II (water cloud), and Class I (ammonia cloud) nomenclature used in Sudarsky et al. (2000). Our clear and ammonia cloud models are similar to those in Sudarsky et al. (2000), and our water cloud models at 2 AU are a bit brighter. As shown in Figure 10, distinct Na, K, CH₄, and H₂O features are apparent, particularly CH₄ near 0.62, 0.74, and 0.89 μm (see Table 3). As discussed in Section 4.1, the cloud-free 0.8 AU spectra are dominated by Rayleigh scattering at short wavelengths for both 1× and 3× Jupiters and 10× and 30× Neptunes.

Zoom In Zoom Out Reset image size

At any given planet–star separation, the higher metallicity compositions generally have smaller albedos due to the increased opacity of their atmospheres. The presence of relatively high and thick H₂O clouds in the atmosphere of the 2 AU Jupiters results in a higher albedo across the visible. The bright effect of relatively high H₂O clouds is also apparent through mid-band in the 2 AU Neptunes, however, the absorption features are considerably more pronounced at longer wavelengths for the Neptunes than for the Jupiters. At full phase, backscattering from clouds also plays a role. At separations of 5 AU and 10 AU, the albedos show progressively larger contributions by Rayleigh scattering at short wavelengths, consistent with the lower clouds in Figures 7 and 8. At 5 AU and 10 AU, the difference between albedos for 1× and 3× solar abundances of heavy elements for Jupiters is larger at all wavelengths than the difference between albedos for 10× and 30× solar for Neptunes. In the cooler atmospheres at larger planet–star separations the albedo spectra seem to change little with increasing heavy element abundances above about a 10-fold enhancement.

We compare the albedo spectra for our standard 5 AU 3× enhancement Jupiter model at α = 0° with observed data from the real Jupiter in our solar system at near full-phase from Karkoschka (1994) in Figure 11. This figure illustrates how our interpretation of Jupiter might proceed if we were to detect it as an exoplanet. The general agreement in morphology of the spectral features is simply a consequence of the spectrum being primarily (but not exclusively) shaped by methane. As noted in Section 3.2, we do not include the effect of photochemical products such as hazes that would explain the difference between our model and the observed data at short wavelengths (Marley et al. 1999; Sudarsky et al. 2000). Although our 5 AU 3× Jupiter model was not adjusted to fit the data, varying the cloud thickness or f_sed parameter would brighten or darken the spectra.

Zoom In Zoom Out Reset image size

The difference between our 5 AU 3× Jupiter model and the data from Karkoschka (1994) around 0.94 μm is of interest. As noted above Karkoschka (1994, 1998) noticed features near this wavelength region in Jupiter's albedo spectrum. Karkoschka (1994) mentioned that they could potentially be water features, but in the later paper concluded they were more likely to be ammonia. While the features we see in this region of our models are likely not the same as those observed by Karkoschka (1994, 1998), we were able to confirm, by running models both with and without H₂O as shown in Figure 11, that the model features around 0.94 μm are indeed attributable to H₂O. These features would be interesting targets for direct imaging observations of exoplanets since the appearance of water is a sensitive probe of temperature and composition.

In Figure 12, we show the albedo spectra of the 0.8 AU 3× Jupiter model as α progresses from 0° (full) to 180° (new). As expected, the albedo decreases with increasing α. More importantly, Figure 12 also shows the ratio of the albedo spectra at each α increment to the geometric albedo spectrum at α = 0°. The wavelength-dependent changes are apparent as a function of α. There are both broad effects on the spectrum, such as the drop in brightness from shorter to longer wavelengths becoming shallower as the planet phase increases, and narrower effects on the spectrum due to the changing line depths as different volumes of atmosphere are probed at different phases.

Zoom In Zoom Out Reset image size

Figure 13 shows albedo spectra of the Jupiter and Neptune models at 5 AU for a few different phase angles to illustrate the point noted in Section 4.2, that for planet–star separations beyond 2 AU, where it is cooler, there is not as large a difference in the albedo spectra for the 10× and 30× Neptune models as there is between the 1× and 3× Jupiter models, even at different phase angles. This implies that it will be more challenging to use albedo spectra to differentiate between Neptune analogs with large abundances of heavy elements. In future work, we plan to further increase the metallicity of Jupiter analogs to determine to what extent this occurs for the Jupiters as well as Neptunes. We also plan to investigate decreasing the metallicity of the Neptune analogs to examine the range of metallicities that we can differentiate between at cooler, larger planet–star separations. For the Jupiter case, it appears that the difference in metallicity is easier to detect when probing deeper into the atmosphere at smaller phase angles than at larger phase angles.

Zoom In Zoom Out Reset image size

From the albedo spectra as a function of phase, we can also generate phase functions, Φ(λ, α) as in Equation (3). Although we can calculate a phase function for each wavelength individually, we instead present the phase functions for our model exoplanets using the Johnson–Morgan/Cousins UBVRI filter passbands (Fukugita et al. 1995) as shown at the top of Figure 14. Also shown for reference in Figure 14 with the filter responses⁵ is the albedo spectrum of the 0.8 AU 1× Jupiter model. Note that our models extend shortward only to 0.35 μm. Since they do not extend to the shortest wavelengths of the U filter, we do not use U in this work. Figure 14 also shows the solar spectrum used for the parent star,⁶ and the reference spectrum⁷ used to compute color magnitudes. Direct imaging observations will obviously be made with different filter responses, possibly tuned specifically to achieve the particular scientific objectives of the observation. However, the approach presented here would be applicable to any arbitrary filter set as well (see the Appendix).

Zoom In Zoom Out Reset image size

In Figure 15, we first compare the phase function for our 10 AU 10× Neptune model with data points of both Uranus and Neptune from Voyager 1 presented in Pollack et al. (1986). The data points and our phase function are shown with a Lambert phase function for reference. While the error bars are relatively large, and there are two Uranus data points that are closer to the Lambert curve than our model, there is general agreement at higher phases. Since these particular Voyager 1 data were taken with a clear filter,⁸ the model phase function shown is an average of phase functions over all wavelengths. We did not include the phase function for Jupiter used in Dyudina et al. (2005) and Sudarsky et al. (2005) for comparison since the Pioneer data used in the generation of that curve were not disk-integrated.

Zoom In Zoom Out Reset image size

Figure 16 shows phase functions for each model exoplanet and for the B, V, R, and I filters. The phase function is shaped by geometry and the composite of all of the scattering phase functions making contributions to scattering in the exoplanet's atmosphere. For reference, the area beneath a Lambert phase function, or its phase integral q, as in Equation (4), is 3/2, and the area beneath a Rayleigh phase function is 4/3 (q_Lambert>q_Rayleigh). Different line styles represent different metallicities: 1× is solid, 3× is dashed, 10× is dot-dashed, and 30× is dotted. A Lambert phase function is shown in each sub-plot for reference. For the 0.8 AU cloud-free case, it is difficult to distinguish between the metallicities, although the higher metallicity phase functions are generally larger than the lower metallicities. The phase functions start off quite similar to the Lambert phase function at short wavelengths, and become shallower and rise above than the Lambert phase function as the colors move redward.

Zoom In Zoom Out Reset image size

For the 2 AU water cloud case, the phase functions are all just below the Lambert phase function; the difference becomes more pronounced as colors move redward; for the R and I filters, the beginnings of a forward-scattering "toe" are seen at large phase angles. Effects of backscattering would appear at small phase angles as well. For the 5 AU ammonia and water cloud case, the ammonia clouds are higher in the atmosphere than the water clouds, and there is a noticeable difference between the Jupiters and Neptunes, likely because the ammonia clouds are a bit higher for the Jupiters than they are for the Neptunes (also true at 10 AU, see Figures 7 and 8). At shorter wavelengths at 5 AU, the Jupiters fall below the Lambert curve and the Neptunes land closer to it, although the Neptune curves also fall below the Lambert curve as colors get redder. The forward-scattering "toe" becomes more pronounced. At 10 AU, both the ammonia and water clouds condense at much higher pressures (lower altitudes), particularly for the Neptunes. The Neptune phase functions are below both the Jupiter phase functions and the Lambert phase function. We note that our phase functions are not nearly as steep near full phase as are those in Sudarsky et al. (2005). We will further examine the dependence of phase function morphology on the assumptions made in the treatment of the direct and diffuse scattering functions in a future paper. We discuss in the following section the use of color–color diagrams to help distinguish between the different types of exoplanet models, as it would be difficult to observe and confirm the small differences between the phase functions shown in Equation (16).

From the phase functions and the geometric albedos, we can use the equations outlined in Section 2.2 to calculate the phase integral q, spherical and Bond albedos A_s and A_B. Our Bond albedo calculation integrates over 0.35–2.5 μm. Contributions from longer wavelengths are negligible (less than 5%) for all cases since both the albedo spectra and incident flux drop off quickly at longer wavelengths. Table 4 contains geometric albedos calculated for each B, V, R, I color filter for each of our models. Table 5 contains the corresponding phase integrals q. The Bond albedos in Table 6 yield the effective temperatures in Table 1.

4.4.1. Bond Albedos and Implications for Dynamics

We noted in Section 2 that even for the low-resolution spectra of solar system outer planets in Figure 2 generated from Karkoschka's (1994) high-resolution spectra, it was still possible to discern between gas-giant (Jupiter, Saturn) and ice-giant (Uranus, Neptune) planets. This is still the case in Figure 17, where the Jupiter models are generally brighter than the Neptunes. For solar system outer planets, the real Jupiter and Saturn are darker at short wavelengths due to absorption by photochemical hazes, but are brighter in the middle of the optical. The real Neptune and Uranus are brighter at short wavelengths and darker at longer wavelengths, due to stronger absorption features in the red. The brightness of a given planet depends sensitively on the height and optical thickness of the uppermost cloud decks, however, so one could imagine planets with somewhat different cloud structures from those presented here. The actual flux measured by a camera depends on the product of a planet's albedo, cross section, and phase. Within the limitations of the models (which do not include photochemistry) ice giants generally will be more difficult to detect than gas giants at longer wavelengths both because of their lower albedos and their smaller radii.

Zoom In Zoom Out Reset image size

Figure 17 again illustrates the difference between 1× and 3× solar abundances for the Jupiter models which is larger than the difference between 10× and 30× for the Neptune models. In the R = 15 spectra in Figure 17, it is also possible to see the effect of the strong methane absorption in all of the spectra between 0.7 and 0.9 μm. Placement of relatively coarse filters both centered on these features as well as on the continuum around them should make it possible to detect and possibly differentiate between these features in early direct imaging observations.

As mentioned in Section 4.4, to generate illustrative color–color diagrams of exoplanets, we use the Johnson–Morgan/Cousins UBVRI filter passbands (Fukugita et al. 1995) as shown at the top of Figure 14. In Figure 18, we first use the high-resolution solar system outer planet data from Karkoschka (1994) to compute the colors for Jupiter, Neptune, Saturn, Uranus, and the moon Titan (see the Appendix for methods). These serve as a useful reference for comparison with the model Jupiters and Neptunes. We then compare these with the (V − R) and (R − I) colors of the exoplanet models. Different symbols represent different planet–star separations: circle for 0.8 AU, square for 2 AU, triangle for 5 AU, and diamond for 10 AU. Red represents the 1× Jupiters and magenta the 3× Jupiters. Blue represents the 10× Neptunes and cyan the 30× Neptunes. The symbols are strongest at α = 0° and fade with α increasing to 180° in 10° increments.

Zoom In Zoom Out Reset image size

The cloud-free 1× model at 0.8 AU exhibits behavior that is generally consistent with that in Sudarsky et al. (2005), ranging from (R − I) of about −1 to −0.5 and (V − R) of about −0.4 to 0.2 with increasing phase. The 0.8 AU 3× Jupiter model and both Neptune models are also consistent with this behavior. The 1× and 3× Jupiter models at 2 AU and 5 AU at full phase are quite close to the real Jupiter. The 10× and 30× Neptune models also approach the observed Neptune, which is encouraging given that the brightness for the models is a bit darker than the observational data. We also note that for the cloudy Jupiter models, it is difficult to differentiate between the models at planet–star separations wider than 2 AU on the (V − R), (R − I) color–color diagrams.

While additional work remains to be done on both generating and interpreting exoplanet albedo spectra models, another relevant question for color–color comparisons is which colors (parts of the spectrum) are most informative. Making color–color comparisons using the filters that best allow differentiation between planet models also requires further optimization. In Figure 19, we show four additional color–color diagrams using different filter pairings. From this figure, it appears that the (B − V), (R − I) and (B − V), (V − I) diagrams provide the greatest spread between the planet types. These diagrams are also useful for determining how well the models capture the integrated effect of different broad regions of the spectrum in comparison to the solar system outer planet data.

Zoom In Zoom Out Reset image size

It has been suggested by Traub (2003) based on the upper left panel of Figure 18 that color alone will be a definitive test of planet characteristics, for example by distinguishing terrestrial, ice-giant, and gas-giant planets. Figures 18 and 19 suggest that such an approach is not as straightforward as it might appear. Because planetary phase functions vary with wavelength the planet colors also depend upon phase. Furthermore, colors vary as cloud heights vary with atmospheric temperature and composition. Thus for example, at 5 AU in the (V − R) versus (R − I) color space, the distinction between the Neptune and Jupiter-like models is slight, particularly if the phase is unknown. We conclude that a single color–color diagram is likely inadequate to sort planets with certainty and the more photometric bandpasses that are available, the better. Of course, the best characterization of an exoplanet will be with spectral information, preferably with R>15 (Figure 17).

In this work, we have investigated the roles of planet–star separation, heavy element abundance, and planet phase on the reflected-light albedo spectra of Jupiter and Neptune analog exoplanets. Albedo spectra of gas-giant exoplanets are sensitive to both the presence of clouds and the changes in optical depth that occur as a function of planet phase, when different volumes of atmosphere are sampled corresponding to the source–observer geometry. The planet–star separation (temperature) plays a primary role in determining the cloud structure in our model exoplanet atmospheres. The abundance of heavy elements (metallicity) plays a supporting role, with clouds typically extending deeper to higher pressures and higher opacities with increasing metallicity.

Marley et al. (1999) and Sudarsky et al. (2000, 2005) have previously explored the effect of changing separation on albedo spectra of 1 M_J and more massive planets. We expand upon this work, developing a different model implementation to include observational geometry, exploring the effects of changing atmospheric composition and adding Neptune-mass planets. In total, we present 16 different model cases that span a range of parameter space consistent with the expected capabilities of early direct imaging observations of exoplanets in reflected light. The 16 cases include a range of planet sizes (Jupiter and Neptune analogs), planet–star separations (0.8–10 AU), compositions (heavy element abundances between 1× and 30× solar), and planet phases (from 0° to 180° in 10° increments). The exoplanet models capture cloud-free, H₂O, and H₂O + NH₃ cloud cases.

As previously established, clouds play a key role in controlling the brightness of a planet seen in scattered light. Water clouds first form high in the atmosphere for planets found beyond about 1 AU. With increasing planet–star separation, the water clouds are found progressively deeper in the planetary atmosphere until ammonia clouds form above them at planet–star separations ∼5 AU. For cloudless planets, gaseous alkali elements, particularly Na and K, are important opacity sources along with water. Once clouds form, the main absorbing species is methane. However, the reflected-light spectrum of even Jupiter itself is not purely a methane spectrum, for example, it also contains ammonia absorption features (Karkoschka 1994, 1998) and in our model, even water absorption features. The computed energy budgets for our mature model gas- and ice-giant planets also imply that there should be substantial equator to pole temperature gradients and thus atmospheric dynamics that will differ from solar system giant planets.

We find that, as in the solar system, the reflected planetary spectra are sensitive to planetary phase as planets are not simply Lambertian scatterers at all wavelengths. As a result, knowledge of the phase at which the image of an exoplanet was obtained is a key ingredient to the interpretation of the observation. Planet phase will be trivially known for planets that have also been detected by RV observations. It will be more difficult to interpret observations of exoplanets seen in only one or two images (perhaps because the planet has passed inside or outside of the IWA or OWA of the coronagraph) and for which no RV data are available.

Ratios of albedo spectra at taken different planet phases clearly have wavelength-dependent features. Making multiple measurements of exoplanets at different phase angles, taking the ratio, and examining the relative changes in the depth of absorption features can help to characterize the vertical distribution or scattering dependences of atmospheric constituents in an exoplanet's atmosphere. The spectra generated for these exoplanet models can be used to estimate the sensitivity required to detect this type of phase-dependent effect in direct imaging observations.

As has been noted elsewhere (e.g., Marley 1997; Marley et al. 2007), photochemical hazes will add a challenging aspect to any effort to interpret cool giant exoplanet spectra. As the comparison of Jupiter and the model Jupiter in Figure 11 attests, the neglect of photochemical hazes is a first-order discrepancy in the model. A Neptune at 1 AU around a solar-type star, with a rich brew of not only methane, but also water and $\rm H_2S$ in its atmosphere will unquestionably exhibit a complex photochemistry and thermochemistry. Zahnle et al. (2010) have only begun to explore the chemical processes that might take place in such atmospheres, but regardless we can expect that actual spectra may deviate strongly from those shown here. If so this will add new layers of insight to that gained from the study of the scattering spectra of giant exoplanets.

Coarse spectra and color–color diagrams are a first step in characterizing an exoplanet's atmosphere. We use standard UBVRI filters to calculate the colors for each of our exoplanet models at different phase angles as well as for solar system planets Jupiter, Saturn, Saturn's satellite Titan, Neptune, and Uranus. Exoplanet colors are highly dependent on planet–star separation (clouds) and planet phase, and less sensitive to metallicity, although it may be possible to differentiate between high and low metallicity for certain filter pairings. It is important to make multiple observations and to obtain as much information as possible about the geometry of the system in order to interpret the observed spectra.

There are several areas that could be further developed based on this initial work. The models and methods presented here can be used to study systems using different stellar spectra and with different instrument-specific filters. Signal-to-noise ratio estimates can then be computed for spectral features of interest, such as Na, K, CH₄, NH₃, and H₂O at appropriate resolutions to help inform instrument specifications and observing plans. The exoplanet albedo spectra can be combined with model spectra of exozodiacal dust and used as input to instrument performance simulations to investigate the detectability of exoplanets as a function of wavelength in systems with substantial circumstellar debris structure.

Additional exoplanet model "points" in the metallicity and planet–star separation space can be considered. Given the relatively small effect of increasing from 10× to 30× solar metallicity on the albedo spectra for Neptunes at separations larger than 2 AU, and yet the noticeable effect of increasing from 1× to 3× solar metallicity for the Jupiters, it would be interesting to adjust abundances and planet–star separations to determine where the effect begins to diminish.

Cloud formation has a strong effect on albedo spectra. It would be instructive to consider the effect of different sedimentation ratios for these exoplanet models. It would also be interesting to adjust the vertical extent of the clouds to investigate the effect of longer path lengths through the clouds at different phase angles and for different vertical distributions. A similar investigation to compare the absorption features of well-mixed versus stratified or pressure-dependent atmospheric species as a function of phase would also be useful in preparing to characterize exoplanets with direct imaging observations. These types of investigations would benefit from updating the current cloud model to output the full anisotropic phase function in addition to the asymmetry parameter and single-scattering albedo.

In this work, we restricted the system geometry to the simple case of circular orbits with no other variation in the Keplerian elements (such as eccentricity, inclination, or obliquity). Investigation of more complex orbits and the effect of the changing planet–star separation through the orbit on the T(P) structure of the atmosphere would be informative. The time dependence and evolution of an exoplanet's atmosphere through a more complex orbit could have a substantial effect on the observed albedo spectra. The obliquity of a planet and the potential for seasonal change would further complicate the interpretation of albedo spectra.

With the ability to create many facets over an illuminated sphere, we could use models for several types of facets to paint a portrait of a non-uniform exoplanet on the sphere. The contributions from the different facets would then be integrated to generate the albedo spectra. For example, we could create facets that correspond to banded structure, or other variations in cloud topology.