Simulating Reflected Light Exoplanet Spectra of the Promising Direct Imaging Target, υ Andromedae d, with a New, Fast Sampling Method Using the Planetary Spectrum Generator

We use the PSG to examine these sampling techniques and to calculate the spectra in this paper. The PSG is a flexible radiative transfer suite that allows users to implement targeted observing scenarios through the integration of a range of spectroscopic, atmospheric, and instrument databases. PSG enables users to synthesize a broad range of spectra through a user-friendly web interface to these models and databases. The full description of the tool is beyond the scope of this paper, but an overview of PSG is given in Villanueva et al. (2018). Different modules allow the user to specify a particular scene, either directly through the web interface or through the application program interface, which enables queued runs of multiple simulations. This is possible by creating scripts that can modify configuration (labeled "config" from now on) files that PSG uses as input for a simulation. Each simulation run on the online interface of PSG also allows a user to download the config file associated with the simulation. We provide several config files in the appendix as a reference for the reader. Users have the ability to customize a number of different parameters in a simulation, with modules that enable customization of orbital properties and observational geometry, surface and atmosphere properties, and instrument properties. As a visual example, Figure 1 displays the module setup that was used (in addition to simulations using scripts that called the API) in the Jupiter spectra validation simulations in Section 3.

Zoom In Zoom Out Reset image size

Users may also upload settings for specific observational targets using preloaded templates or a lookup function where PSG extracts orbital parameters from the NASA Exoplanet Archive (Akeson et al. 2013) for known exoplanets. The atmosphere and surface module also contains a number of preset atmospheric templates the user can choose, including many specifically designed for or relevant to exoplanets. Two examples are a set of templates that compute temperature/pressure profiles (Parmentier & Guillot 2014) and line-by-line abundances (Kempton et al. 2017) for gas-giant exoplanets with varying chemistry, and another set that produces templates for terrestrial exoplanets (Turbet et al. 2017).

A common method of integrating the flux for these planets involves sampling the emergent flux from a sphere using many plane-parallel facets, where facets correspond to different pairs of incidence and emergence angles. This has been used in numerous studies (Cahoy et al. 2010; Webber et al. 2015; MacDonald et al. 2018) examining reflected light spectra of directly imaged exoplanets. In these and many other studies, the method has relied upon the ability to formulate the planetary albedo in a manner that can employ common numerical integration techniques to solve the integral for the flux. In this case, Chebyshev–Gauss quadrature using Chebyshev polynomials of the second kind is used to yield an exact result for polynomials of degree 2n – 1 or less by choosing suitable Chebyshev and Gaussian nodes/angles and weights at which to evaluate the flux (Horak 1950; Horak & Little 1965). These are available online for any given number of nodes or can be calculated using formulas given in Webber et al. (2015). Precision to the real value can be improved by increasing the number of nodes at which flux is evaluated, and numerical techniques for improving precision are also available (Dehghan et al. 2005). A visual example of the distribution of nodes using Chebyshev–Gauss integration using a 10 × 10 (100 total) set of nodes on a sphere at full phase is given in Figure 2(a).

Zoom In Zoom Out Reset image size

Chebyshev–Gauss integration is easily implementable in PSG by calling the API and was used in order to calculate reflected light spectra and validate accuracy and efficiency of the subsampling method described in Section 2.3. Because simulations use a finite beam size (for example, see the annular beam size in white around Jupiter in the left portion of Figure 1), the integrated flux calculated using Chebyshev–Gauss integration needs to be normalized by the beam size, which with a circular beam is just the area of the spherical cap. In addition, Chebyshev–Gauss points across the longitude need to be normalized to their limb-to-limb chord length, which is dependent on the latitude of interest and is largest at the equator. In order to validate that Chebyshev–Gauss integration was being implemented correctly, we calculated fluxes for different scenes in order to compare to theoretical limits and to other data and simulations. We validated our implementation of the Chebyshev–Gauss integration with PSG/API by running both it and the subsampling method we developed, for a perfect reflecting surface (albedo = 1) and iterating over the number of sampling points. We compared the albedo using our sampling method to the analytic Lambert scattering phase function as given in Madhusudhan & Burrows (2012) and compared all the values at zero phase, where all three methods converged to the theoretical 2/3 limit as the number of sampling points increased (see Figure 3).

Zoom In Zoom Out Reset image size

In order to enable a more rapid but similarly accurate integration of disk-integrated flux, we develop a numerical algorithm that creates a 2D matrix of regions of similar outgoing flux based on the stellar incidence angle and observer–atmosphere angle. The disk is then divided across this 2D array by selecting the number of subsamples that are used to encapsulate regions and define contributions/weighting functions for each eigenvalue. Internally, PSG divides the sampled disk into a map of 140 × 140 pixels and computes incidence and emission angles for each pixel (19,600 sets). Without subsampling, PSG computes a single set of effective incidence and emission angles from these 19,600 sets, which is then used to compute a single radiative transfer calculation. When subsampling is enabled, PSG creates histograms of incidence/emission angles, with N defining the number of bins between 0° and 90°. It then identifies the different possible combinations, and the relative occurrence of each incidence/emission combination (among the 19,600 sets), ultimately establishing the weight for those supersets. PSG then computes radiative transfer calculations using those supersets and adds the spectra for each case yet weighted by the numerically computed weights. For instance, for a symmetric simple case (subobserver angle = substellar angle = 0), we are in full symmetry, and only the diagonal elements of the matrix have a weight. So instead of running 100 (10 × 10) simulations, simulations only need to be run for the diagonal elements (10). This is shown in middle image on the right side of Figure 2. The full disk can be divided into 10 concentric rings because each of these regions shares a common incidence/observer angle. In the more general cases, a particular subsample choice leads to a set of distinct incidence/observer angle regions; based on the number of angles subdivisions, PSG then chooses how many simulations are needed and the weight for each subregion. For example, because these angles range from 0 to 90, a subsample choice of N = 5 would lead to a subdivision with bins of 18° and 22 distinct radiative transfer regions, as shown in the rightmost image in Figure 2.

Using this technique, PSG is able to accurately subsample planets' heterogeneous geometries in single simulations—but with an important implicit assumption of homogeneous atmospheric/surface properties in each subsample region. The ability to leverage the symmetry of observational properties in this averaging scheme is based upon the assumption that the radiative properties of a region vary smoothly across the region—versus potential spatial variations due to inhomogeneous clouds or spikes in surface albedo. While the number of subsamples chosen using this method can be increased in order to capture the potential effects of inhomogeneities, this comes at a higher computational expense. Instead, inhomogeneous atmospheric/surface properties can also be simulated using GlobES (Global Exoplanet Spectra, https://psg.gsfc.nasa.gov/apps/globes.php) or via the API to call simulations for different regions and then weighting them appropriately.

In order to validate and compare the performance of the subsample method to the Chebyshev–Gauss method, we performed radiative transfer simulations for a Jupiter-like atmosphere employing both methods at different sampling resolutions. We used profiles for atmospheric temperature, pressure, vertical mixing ratio, and cloud density from MacDonald et al. (2018, subsequently referred to as M18); we also used the same opacity catalog as M18 to ensure that no deviations resulted from different opacity assumptions. We produced simulated spectra and examined comparisons for both a cloudy and cloud-free Jupiter simulation, in order to examine how the differences between the two methods propagated into a realistic simulation. Additional details of these validation simulations are given in Section 3.

Results of sampling resolution tests for the Jupiter simulation with clouds are shown in Figure 4, and results for the cloud-free simulations were also very similar. As is evident from the figure, albedo calculations for both the Chebyshev–Gauss and subsample method converge to very similar values with a relatively small number of nodes/samples. The accuracy of the simulations using the PSG subsample method is within 5%, 3%, and 0.2% of the integrated flux respectively for the one-sample, four-sample, and six-sample simulations relative to the 10-sample simulation. The accuracy of the simulations using the Chebyshev–Gauss method is within 5%, 1.5%, and 1.5% of the integrated flux respectively for the 1 × 1 (1 node), 4 × 4 (16 node), and 6 × 6 (36 node) simulations relative to the 10 × 10 (100 node) simulation. The accuracy of particular Chebyshev–Gauss and subsample node selections relative to each other is within 3% for the integrated flux for the 10 × 10 (100 node) and 10-subsample case. These differences in the average and integrated flux between the two methods are small—and are likely due to differences in which they sample the atmosphere for radiative transfer calculations. For example, in both the subsample and Chebyshev–Gauss methods implemented in PSG, there is a finite region that is used to calculate the flux for each region (the subsample region and the annular beam of the Chebyshev–Gauss point) before summing or weighting. This discrete region will produce small variations in the flux that become less significant with a larger number of samples (particularly for the Chebyshev–Gauss case as sampling increases toward the limb). In addition to that, the subsample method's averaging algorithm will inherently converge to a limit with a large number of samples, which is also evident from the comparison in Figure 4.

Zoom In Zoom Out Reset image size

Finally, the subsample method was also used in the validation case shown in Figure 5, where PSG was used to simulate a disk-integrated Jupiter albedo spectrum from Karkoschka (1994). As discussed in the validation section (Section 3), the simulation was able to reproduce the Jupiter spectra to a great degree (with albedo variations of ∼0.026 on average; see the difference as a function of wavelength in the right panel of Figure 5), with small differences likely driven by inhomogeneous variations in atmospheric profile and cloud and haze properties as a function of location on Jupiter in both latitude and longitude. Based on these tests, the subsample method appears to be able to accurately simulate reflectance spectra, similar to the Chebyshev–Gauss method.

Zoom In Zoom Out Reset image size

The computational time for the subsample method is less expensive than the time required for Chebyshev–Gauss sampling in PSG (see Table 1). While this is to be expected given that the technique was developed to take advantage of symmetries in the emergent flux from a planet, we also ran comparisons of the computational expense of the two methods for several different cases. We ran these simulations on a personal laptop workstation with the following computer specifications: Macbook Pro with 3.1 GHz Dual-Core Intel Core i7 processor. Computer code calling the API to PSG for both methods was written in Python. We calculated the time comparisons for the previous simulations (a cloud-free Jupiter case, a cloudy Jupiter case, and a case where a haze was added to the cloudy Jupiter profile) using different subsamples. In all three cases the relative computational expense was generally the same. The cloudy Jupiter simulations are given in Figure 4 and the corresponding computational expense of the runs (as well as the other cases) are given in Table 1. The lower computation time required for the subsample method approaches approximately an order of magnitude as a greater number of samples are used, where the computational expense of the subsample method is $O(n\mathrm{log}(n))$ while that of Chebyshev–Gauss is O(n²). Given the similar accuracy of the subsample method and significant gain in computational expense, we use the subsample method in simulations described in the rest of the paper.

υ And d is an attractive target for the Roman mission and other future efforts due to its favorable observational properties. The 2015 Roman Science Definition Team report (Spergel et al. 2015) lists some of the most promising targets for the observatory from the sample of known radial-velocity planets at that time. In Tables 2–7 of the report, observational parameters for some of the most favorable targets are listed, including separation from the host star, relative contrast to the primary, and integration time required to obtain an S/N of 5 at 565 nm with a 10% bandpass using the Coronagraph Instrument (CGI) Hybrid Lyot Coronagraph (HLC). While the integration times made simple assumptions for planetary albedo and noise contributions, the relative values give a useful measure of the potential efficacy of observations of different planets. Due to the υ And A system's distance and host star's brightness, the integration time required for observations of υ And d to reach the required S/N is one of the shortest among all of the listed planets. The angular separation for υ And d is fairly small at 01805 for mean orbital separation but is larger than all of the planets with shorter integration times in Tables 2–7 from Spergel et al. (2015), which makes it an advantageous target with respect to inner working angle considerations relative to brighter planets. The contrast ratio for υ And d is also favorable and is better (using the assumptions listed in the report) than other common targets of simulation, such as 47 UMa c and υ And e. Indeed, all planets in Tables 2–7 from Spergel et al. (2015) with better contrast ratios are at shorter angular separation <016 from their host star than υ And d. Some early simulations (described in Section 5.2) of the S/N also reinforce the observational favorability of υ And d. Simulated S/N predictions using Jupiter-like spectra for all the then-known RV planets at a number of different bandpasses also demonstrated the relatively high S/N that could be obtained for υ And d (Lacy et al. 2019). This is especially true at the bluer wavelengths relevant to the current imaging and spectroscopic capabilities of Roman (though Lacy et al. 2019 found that at redder wavelengths, the assumed spectra for υ And d yielded less favorable S/N). This combination of observational qualities motivates us to simulate a range of different υ And d spectra that reflect different atmospheric states (given in Table 2), in order to help plan for potential future observational programs.

Table 2.  υ And d Parameters Used for Spectra Simulations

υ And d Parameters
1Observation Derived
Planetary Mass—${10.25}_{-3.3}^{+0.7}$ MJup
Planetary Radius—1.02 RJup
Semimajor Axis—2.53 au
Orbital Eccentricity—0.316
Orbital Inclination—23758
Model Parameters
Teff (max, min, mean), Tint = 0–260 K, 188 K, 215 K
Teff,mean inc. 2 Tint (10.25Mj , 6.95Mj )—319 K, 270 K
Planetary Gravity (g)—244.23 m s−2
3Metallicity ([Fe/H]star = 0.131)—1/5/10/15×
Water-cloud Particle Size—0.1/1 μm
Haze Properties—no haze/Tholin Haze
1 McArthur et al. (2010); Deitrick et al. (2015)
2 Marley et al. (2018)
3 Gonzalez & Laws (2007)

Download table as:  ASCIITypeset image

There have been some earlier simulations of the spectra of υ And d. Lacy et al. (2019) was the first study to specifically examine the S/N for υ And d assuming both a Jupiter- and Neptune-like spectra, but earlier work by Sudarsky et al. (2000, 2003) did try to produce generic spectra for the class of planets that υ And d may fit into. These works used atmospheric structure and radiative transfer modeling to produce broad extrasolar giant-planet classifications based on qualitative similarities in the composition and spectra of objects within broad effective temperature ranges. In Sudarsky et al. (2003), υ And d was modeled as a "Class II" "water class" planet characterized by tropospheric water clouds. Simulations propagated the "Class II" planet spectra using the system and orbital constraints at that time. While the spectra produced were for a Jupiter-fixed orbital radius, phase-averaged planet, Sudarsky et al. (2003) did acknowledge that the substantial eccentricity of the planet means it may be too warm in its outer atmosphere to possess water clouds at periastron and cross over into the water-cloud-free "Class III" planet regime. This potential variation in planetary and atmospheric parameters due to υ And d's eccentric orbit is just one of the interesting facets of the system that suggest more detailed modeling may be warranted.

Planets such as 47 Uma c and υ And e may be more similar to and are often modeled as Jupiter analogs, but υ And d represents a giant exoplanet that exists in an atmospheric phase space not seen in the solar system (for example, see Figure 1 of Lupu et al. 2016). Its orbit places it in the conventional liquid water habitable zone of its system (Buccino et al. 2006). The variation in separation due to eccentricity may mean that the planet experiences significant variation in the existence and location/morphology of clouds in its atmosphere, which may modify its albedo. The planet may transition from a relatively cloud-free "Class III" giant planet at periastron to a water-cloud-possessing "Class II" planet at apastron, which would result in a shift from a relatively low-albedo average reflectance to a high albedo average reflectance (indeed, Sudarsky et al. 2000 suggest planets with high water clouds may be brighter in the visible than Jupiter analogs); brightening would occur as the planet's orbital separation approaches a maximum, which would also make it less likely to be interior to the inner working angle of the Roman CGI coronagraphs. The variation in temperature and consequently cloud location as a function of phase may therefore produce a measurable effect on the planetary spectra. Finally, as is the case for gas giants in the solar system, υ And d may possess hazes that alter its spectra. Given the location of υ And d and the common presence of moons and subsequently the infall of their material on their giant-planet hosts in the solar system, such hazes may contain albedo-influencing chromophores that may be either exogenously or endogenously sourced. This suggests that there are a number of different types of potential atmospheres that υ And d may possess, which guides our simulations of the planets' spectra in the rest of this section.

In order to simulate spectra of potential atmospheres of υ And d, we first obtain an atmospheric temperature and pressure structure for the planet. This is based on effective temperature and planetary gravity parameters that are dependent on planetary radius. Because there are no observational constraints for the radius, we use models that explore the mass–radius relation for giant planets in order to produce an estimate for the radius. While the radius is a function of a number of different variables, including age of the system and core size (Fortney et al. 2007, 2008), most models of the mass–radius relation find that the radius of evolved gas giants in the range of Jupiter's mass are not likely to exhibit large radius variations with mass (Bashi et al. 2017). Using relations (Fortney et al. 2007, 2008) for evolved planets with masses and orbital separations relevant to υ And d, we use a value of 1.02 R_Jupiter for the planet. This yields a value of g ∼ 245 m s⁻², which is an order of magnitude larger than Jupiter. With these values we can use the analytical fits provided by M18 in order to produce a temperature and pressure profile for υ And d. M18 created these by fitting a number of self-consistent T/P models using the methods described in Fortney et al. (2008). While the models use a maximum of 100 m s⁻² for the planetary gravity, we test lower values of gravity between that value and the calculated value for υ And d and find it does not result in a significant difference in the profile and potential locations of cloud formation. Due to the substantial eccentricity and subsequent variation in orbital separation from its host star that υ And d possesses, it is possible that the atmospheric profile for the planet may vary from periastron to apastron. A calculation of effective temperature (T_eff) indicates a variation of ∼70 K from periastron to apastron with a mean T_eff of 215 K. In this case, we do not include the potential contribution of internal heating (T_int) to T_eff given the uncertainty of υ And d's internal properties. The value of T_eff for the solar system giant planets is approximately 10–20 K greater than their T_eq—such an increase would not change the T_eff of υ And d into a significantly different condensation regime for apastron and mean separation. For periastron, such an increase would potentially result in a cloud-free atmosphere at larger metallicity values. In order to explore the potential effect of this variation on the atmosphere, we plotted a T/P profile for the planet (with several metallicity values for the planet as a function of the stellar metallicity for υ And A) for the corresponding effective temperatures at periastron, apastron, and its mean orbital separation (listed in Table 2). The profiles are given in Figure 7, with condensation curves for water and ammonia overlaid. While the profiles at the orbital extrema are unlikely to exactly match the actual profile given the latency inherent in the structure of the atmosphere due to past conditions, these profiles give general guideposts for the atmospheric parameters that may control observables. We do also include simulations (given in Table 2) with modeled T_int values for the constrained mass and age of υ And d. Those values are based on the closest values that correspond to υ And d in the Sonora-Bobcat tables (Marley et al. 2018) and are simulated for the reference mass value but also the upper and lower bound on mass for the planet. We then use these T_int and T_eff values in order to simulate spectra for these cases and examine potential observations by Roman.

Zoom In Zoom Out Reset image size

We also examined the general effect of gravity and metallicity on the profiles. As discussed above, in general there was no significant variation of where profiles intersected with condensation curves as planetary gravity was varied. At significantly higher values (2× the value for υ And d's observation derived mass) of gravity than the value calculated, the periastron and apastron profiles converged toward the mean separation profile with a temperature difference of less than 10 K above approximately 1 bar. The intersection of the water condensation curve with the temperature profile did not vary significantly, and the main difference was that at around the 50 mbar level and above, temperatures were actually higher for the apastron case versus the periastron case, though the difference was less than 5 K and both were close to the mean separation value (the periastron temperature decreased with higher g). At lower values of g, the difference in temperatures at the orbital extrema was more pronounced, as apastron temperatures decreased and periastron increased (for example, at Jupiter-like values of g—unrealistic given observed parameters—stratospheric temperatures increased by about 20 K at periastron, and apastron temperatures decreased by a smaller amount). At these lower values of g, the periastron profile did not intersect with the water condensation curve, suggesting a potential lack of water clouds.

Metallicity had a more significant effect on the profiles for observed metallicities in the solar system gas giants (Wong et al. 2004; Fletcher et al. 2009). Figure 7 shows profiles for the three orbital distances at different metallicities, with increased metallicity denoted by increasingly transparent lines (temperature/pressure profiles for additional cases with modeled T_int values are given in Appendix A). There are some pronounced effects of increased metallicity—notably that at metallicities near or greater than Jupiter-like enhancement values, the periastron T/P curve may suggest a water-cloud-free atmosphere. Indeed, even the mean separation T/P profile at 15× metallicity is close to a water-cloud-free case, suggesting that more enhanced metallicity values would lead to a planet where water clouds may only be favored as the planet approaches apastron (tests of much higher metallicity values, >25, show that such a planet would still likely possess a thin layer of clouds at mean separation, but that such clouds would become thinner). The actual process of formation and destruction of these clouds is beyond the scope of this paper but would be dependent on the stellar irradiation environment, and sedimentation and production timescales associated with cloud formation. These differences also result in variations of the intersection of T/P profiles with the water condensation curves of over an order of magnitude in pressure for potential deeper clouds. The variation of these profiles and subsequently the atmospheric features with realistic metallicity values necessitates models that consider a range of different atmospheric properties.

In our simulations of υ And d's atmosphere, we first used a thermochemical equilibrium chemistry code in order to compute the molecular abundances of the likely dominant molecules in the atmosphere. We used the package GGchem (Woitke et al. 2018), an open-source computer code that can determine the chemical composition of gases in thermochemical equilibrium down to 100 K. GGchem allows users to choose elements and their associated abundances, sources for equilibrium constants, and temperature and pressure settings when setting up a model. This enables chemistry calculations for atmospheres of varying metallicity with a given temperature–pressure profile for different vertical layer representations. The software also includes an equilibrium condensation code (which we do not use in this study) in addition to the gas-phase equilibrium chemistry component. The software is especially useful for υ And d as it can simulate equilibrium chemistry at temperature ranges (∼160–700 K) relevant to the planet (shown in Figure 7).

We ran GGchem for the T/P profiles we obtained for different metallicity values using a 120-layer atmosphere that is evenly distributed in pressure log space from 10⁻⁷ to 10³ bars. However, our radiative transfer calculations, and consequently our simulated spectra, are based on a maximum pressure limit of 10 bars (which is often taken as the troposphere boundary or atmosphere surface for the solar system gas-giant planets; Seiff et al. 1998; Encrenaz 2004). Once the model is run using a selected list of species and their abundances, an output of the most common molecules at a particular T/P value is produced. We then screened the molecules we would use for the rest of our calculations by only choosing the compounds and elements that exceeded greater than 1 part per million at the 1 bar level for the 1× or 5× metallicity cases. This leads to an atmosphere that is composed of six dominant molecules that comprise the vast majority (>99.999%) of the total abundance: H₂, He, CH₄, H₂O, H₂S, and NH₃. While this threshold does capture the vast majority of the modeled atmosphere composition that is likely to influence spectra, we do note that this excludes Alkali elements and compounds such as sodium and potassium, which may exhibit spectral features for higher metallicity cases.

The metallicity cases that were modeled were the 1×, 5×, and 15× scenarios and were chosen to encompass a range of values that reflect the current understanding of values for Jupiter and Saturn (Atreya et al. 2016). Increased granularity of metallicity values, a wider range of values, and the inclusion of more realistic metal abundances that are less simplistically correlated would all be valuable improvements for future study. These metallicity values and the appropriate T/P profiles for each orbital position were then used to determine the existence and morphology of clouds. As discussed above, the latency in atmospheric properties between different orbital positions may mean some of the input parameters do not instantaneously affect the atmospheric structure. However, as physically motivated fiducial models, such external parameter-driven atmospheric structure can help inform how spectra may potentially change over a sequence of observations.

Cloud thickness, vertical location, and other properties were all chosen based upon physically motivated assumptions taken from modeled profiles. Cloud extent was chosen based upon intersections of T/P profiles with condensation curves, as is visible in Figure 7. Total cloud mass was then chosen in order to obey mass conservation based upon the integrated condensible column above the deepest saturation layer using the equilibrium chemistry calculations of the abundance of the volatile of interest. Given the relevant T/P range, this entails integrating water vapor abundance in condensation regions and then converting to water ice, with a very small, decreasing amount of water vapor retained in the profile. Cloud thickness structure was set with a layer of maximum cloud condensate abundance ratio located where supersaturation was the highest. The condensate mixing ratio in adjacent layers was then set by prescribing gradual attenuation in layers below the maximum layer and a steeper attenuation in layers above (with end points in the condensation range). Variations in the attenuation were tested and do have some effects on spectra, but are generally at levels that do not influence the broader interpretation relevant to near-term observations (for example, the narrowest versus broadest maximum layer schemes had differences less than ∼several percent but one intermediate scheme had larger ∼5%–10% differences in a small portion of the spectra). Finally, given the exploratory nature of these simulations, two particle sizes were tested for the cloud particles, 0.1 and 1 μm. They were chosen as initial exploratory values based upon the upper atmosphere particle sizes measured by the Galileo Probe Nephelometer and to test smaller particle sizes motivated by the evidence of refractory clouds on hot Jupiters (West et al. 2004; Lee et al. 2016). Additional particle sizes and more realistic distributions of sizes will be important variables to test in future work. The water-ice cloud properties and scattering model details are given in the config files included in Appendix C, and PSG allows the user to choose from a number of options for the aerosol properties. In the simulations described in this paper, the refractive indices for the water-ice clouds were taken from Massie & Hervig (2013) and Gordon et al. (2017) while the Mie scattering implementation used was from Bohren & Huffman (1983) and used 20 angles and 200 size bins from 0.005 to 20 μm. The hypothetical effect of a haze layer was the last additional component included in the simulations, with a thin, high layer of tholin haze composed of 0.1 and 1 μm particles (particle size based on the approximate peak of a number of modeled haze particle distributions—see Figure 8 of Gao et al. (2021; specific details on haze particles used are available in the attached config files in Appendix C).

Modeled geometric albedo spectra from a number of selected simulations for the T_int = 0 case are given in Figures 8 and 9. Current Roman CGI filters are overlaid in olive green, yellow, orange, and red. Figure 8 shows the spectra for models that all simulated a 5× stellar metallicity atmosphere for υ And d. These spectra were produced for the planet at three separate orbital positions—periastron, mean separation distance from the host star, and apastron. These different orbital distances corresponded to different effective temperatures and consequently atmospheric profiles. For the periastron case, the atmosphere was too warm for water clouds to condense in the portion of the atmosphere that was simulated. The result is a cloud-free atmosphere that produces a relatively low albedo across the entire visible spectral range, with slightly higher reflectivity in the bluer wavelengths due to Rayleigh scattering. Both the mean separation and apastron cases produced atmospheres with water-cloud condensation—in both cases with high water clouds that significantly increased albedo through the entire visible range. For these cloudy atmospheres, we tested both of the cloud particle sizes described above and, for the apastron case, also included a model with the tholin haze. Apastron cases are given with the dashed lines, while mean separation cases are denoted with the nonblack solid lines. Within those cases, particle size is distinguished by the color of the lines (green is used for 0.1 μm, and blue is used for 1 μm). The justification for the choice of particle sizes is given in Section 5.5. The haze has the expected effect of reducing albedo in the bluer wavelengths (by more than a factor of 2 at 0.3 μm) and having little effect toward redder wavelengths. Unsurprisingly, particle size also has a significant effect on the spectra. While the 1 μm particle size cases tend to have less variation in albedo in the spectral range examined, the 0.1 μm cases exhibit more significant variation from 0.3 to 1 μm, with peak and trough albedo location varying by wavelength due to a combination of both the particle size and the vertical position and extent of the clouds. The brightest albedo spectra across the entire wavelength range are in the case of the planet at apastron with water clouds composed of 1 μm particles. Spectra with the same cloud particle size, except with the planet at mean separation, possess albedos approximately 10%–20% less bright due to changes in cloud height, extent, and thickness. A summary table that details the cloud properties for these simulations is given in Appendix B. This effect of planet orbital distance is more complicated in the 0.1 μm case, as the albedo is higher at bluer wavelengths at apastron and at redder wavelengths for the mean separation case due to the different cloud structure. Absorption signatures, mostly due to methane, are equally prominent in just about all cases, with the small variations due to the vertical location of the cloud decks. Variations between the separate cases at a particular orbital distance (periastron/mean/apastron separations are 1.73/2.53/3.33 au, respectively) may be distinguishable, but at least in these models, some of the most diagnostic wavelengths may be bluer than current Roman filters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows a comparison of the model geometric albedo of the atmosphere assuming different metallicity values for when the planet is at apastron. The left panel shows how the spectra change with varying metallicity while assuming cloud particle sizes of 1 μm, while the right panel shows the same set of cases but with a cloud particle size of 0.1 μm. While particular metallicity cases may exhibit broad patterns of brighter or darker continuum spectra across the entire wavelength range relative to each other, the relationship between metallicity and overall continuum brightness is not a linear one. This is because the relationship of a particular metallicity for the atmosphere and its corresponding T/P profile does not translate to cloud formation with respect to cloud height/extent or density in a simple pattern. Because these cloud properties (see Appendix B and Massie & Hervig 2013; Gordon et al. 2017 for details) all have significant and complicated effects on spectra, continuum brightness can vary in nonintuitive ways. For example, in the 1 μm cases, the continuum brightness of the thicker cloud deck in the 1× metallicity case is greater than that of the 15× metallicity case. However, the 5× metallicity case is the brightest with respect to continuum albedo. The only relatively intuitive relationship is that of the strength of the absorption signatures, which increase with increasing metallicity due to higher molecular abundances. In the 0.1 μm cloud cases, the nonintuitive cloud effects for varying metallicity cases are even more pronounced, as the peak albedo for the 5× metallicity case is offset (as is the entire spectra) from the 1× and 15× metallicity cases. This difference (and perhaps the differences in the 1 μm case) may be due to the location of where the cloud layers in each case become opaque. For example, the 1× and 15× metallicity cases experience a maximum cloud fraction at 0.058 and 0.024 bars respectively, while the 5× metallicity case has a maximum cloud fraction at 0.069 bars. The deeper atmosphere cloud appears to account for the difference in the 0.1 μm curve peaks—artificially setting the highest fraction at 0.344 bars for the 5× metallicity case instead leads to a similar morphology to 1× and 15× cases, as is demonstrated by the "shifted" blue line in the 0.1 μm case. The effect of this demonstration of the impact of maximum cloud fraction and its corresponding effect on opacity and the spectra are also included in the 1 μm case, where the 5× metallicity case goes from the brightest continuum albedo (approximately 0.75) to the darkest one (approximately 0.65) in this contrived scenario. With the current Roman CGI filters (Kasdin et al. 2020) overlaid in these plots as well, it is apparent that depending on the S/N described in Section 5.7, it may be possible to resolve some of these properties based upon the spectra. However, given the dramatic effect that such variations in atmospheric properties and cloud morphology can have on albedo, a more detailed study of the different atmosphere scenarios that are considered plausible is highly warranted.

While the spectra in Figures 8 and 9 were examples of the planet's spectra if it were viewed at full phase, a realistic simulation of υ And d's expected spectra requires consideration of the orbital phase and illumination state of the planet. υ And d is somewhat unique in this regard, in that its three-dimensional orbital structure has been explored due to the existence of mass constraints from RV and astrometry measurements of the system (McArthur et al. 2010; Deitrick et al. 2015). We refer to spectra that incorporate existing constraints on inclination and eccentricity along with phase variation as "phase and illumination-appropriate" spectra. Realistic simulations of the planet's spectrum also need to take into account the field-of-view-dependent throughput (Kasdin et al. 2020) of the instrument used to observe the system. For the purposes of this study, we examine simulated Roman CGI observations of υ And d given current instrument specifications, including two of the baseline observing modes: spectroscopy in Band 3 (675–785 nm) with a Shaped Pupil Coronagraph (SPC) and broadband imaging in Band 1 (546–604 nm) with a Hybrid Lyot Coronagraph (HLC) (Mennesson et al. 2020). ⁵

The relative point-spread function transmission of the Roman CGI Band 3 SPC mode, mapped over the coronagraph field of view, is illustrated in Figure 10(a). In Figure 10(b), we show the sky-projected offset of υ And d as a function of time, color-coded by the phase function of a Lambertian sphere. Finally, in Figure 10(c), we show the combined effect of the Lambertian phase function modulated by the field-point-dependent coronagraph throughput. Together, these indicate the most favorable dates for observing υ And d after the expected 2026 commissioning phase of the Roman Space Telescope span roughly from the beginning of 2028 to the beginning of 2029 (approximately 90°–150° in phase from υ And d's periastron). Using these orbital positions and times as a guide, we then simulate spectra of υ And d for a number of our atmosphere models with the appropriate phase-dependent illumination and instrument-specific throughput with PSG.

Zoom In Zoom Out Reset image size

The most observationally favorable portion of the orbit occurs during the planet's approach to apastron, which is why we typically use the apastron-specific atmosphere models in the simulations examining the impact of phase and illumination. The mean separation cases plotted using geometric albedo in Figure 8 are also potentially a useful model that can be examined given potential latency as the planet moves from periastron toward apastron during the observationally optimal portion of its orbit. However, while the optimal observations may be a little past mean separation, the two cases differ by such a small enough amount that it does not substantially affect our conclusions. We also restrict our simulations in this case to the haze-free scenarios, given the relatively uncertain nature of potential haze existence and properties. Using PSG and the config file included in Appendix C, we simulate the different coronagraph instrument scenes (both the SPC and HLC modes; Mennesson et al. 2020) for our atmosphere model, the appropriate orbital and viewing configuration, and the relevant instrument parameters. This includes the current coronagraph throughput as a function of separation as well as the other optical loss factors and detector noise characteristics. The phase- and illumination-specific albedo spectra are given in Figures 11 (for the T_int = 0 case) and 12 (for the cases where T_int was modeled based upon Marley et al. 2018). There are a number of key takeaways from these simulations that are relevant to potential means of extracting atmospheric properties.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the phase- and illumination-specific albedo spectra for the T_int = 0 case, the 1 and 0.1 μm cloud particle size cases exhibit different continuum flux values in their spectra at the same phase (see Table 3 for Roman-band-averaged values). This is most apparent at the two most extreme phase cases plotted in Figure 11, where differences between the cases are even apparent in the regions overlaid with the Roman filter bandpasses. This also then relates to another distinguishing pattern between the two cases—that the magnitude of total flux variation as a function of phase for a particle size is also significantly different in the two cases. As the planet moves from a phase of 90°–150°, the difference in albedo in some of the Roman filter regions is greater for the 1 μm case versus the 0.1 μm case. The Roman-averaged band albedos given in Table 3 illustrate this—for example, the 825 nm filter captures the larger variation for models with a 1 μm cloud particle size for a phase of 90°–150° from periastron versus the 0.1 μm particle size case (a difference of 0.178 versus 0.103).

The last prominent, potentially useful characteristic of the simulated spectra in Figure 11 is the potential for a color–color change as the planet moves in its orbit. The color–color ratios in the Roman bands for these simulations, as normalized to the 575 nm band, are given in Table 4. As an example, for the 1 μm case, the relative ratio of flux in the bluest Roman filter versus the reddest filter flips as the planet moves from a phase of 90°–150° (shown by the change in the cell text color from vermilion to sky blue in the 1 μm cases). This flip does not occur in the 0.1 μm case and may be an additional potentially useful diagnostic of the atmosphere depending on the achievable S/N for the observation. While these differences and diagnostics in the spectra are relevant to a change in only the particle size assumed in the model calculations, additional variation of model parameters is likely to produce both degeneracies between potential diagnostics and additional pathways to extracting information about the atmosphere.

We also examined the phase- and illumination-specific albedo spectra for cases where T_int corresponds to a planet mass of 10.25 M_j (corresponding to a T_eff,mean = 319 K) and 6.95 M_j (corresponding to a T_eff,mean = 270 K) (Marley et al. 2018). The resulting simulated spectra for a selected number of these cases are given in Figure 12 (a baseline config file used for those cases is given in Appendix C). That figure shows the simulated spectra for a phase of 110° from periastron and in a wavelength region corresponding to the SPC spectroscopy modes. The error bars correspond to 1σ shot noise at R = 25 for a 400 hr observation. The 400 hr observation time is based on the exposure time required to achieve previously stated detection limits for Roman (Bailey et al. 2019) and the approximately 2000 hr/3 months of time CGI will nominally receive (including calibration overheads) during its technology demonstration phase (Akeson et al. 2019). Based on the spectra shown in Figure 12, it is unlikely that spectroscopic observations would be able to extricate variations in metallicity between 1× and 15× the stellar value for the T_eff,mean = 319 K case. However, for cases with a lower T_int and consequently T_eff,mean, the presence of a water-cloud deck that would not exist at the higher temperature cases would be distinguishable. Thus, observations of υ And d using Roman should be able to put a bound on the value of T_int and T_eff,mean for the planet. This is especially true when considering that filter observations would enable shorter observation times and higher S/N for the same albedo values.

Finally, with these Roman scene spectra in hand and additional simulated observations for the Roman mission, we examined the optimal time to observe υ And d and its flux during the potential period of the prime Roman mission from 2026–2031 (Akeson et al. 2019). The left image of Figure 13 displays the PSF-normalized intensity for υ And d for both the SPC and HLC modes as a function of time and Roman observing windows. υ And d peaks in flux very close to the first winter observing window, exhibiting a normalized intensity of ∼5 × 10⁻⁹ that makes it a viable target for both the Band 3 SPC and Band 1 HLC observing modes. Note that the Band 1 HLC mode has a smaller inner working angle (150 mas, versus 190 mas for the Band 3 SPC), so the higher PSF throughput at the planet's angular separation results in a larger normalized intensity, despite a slightly lower albedo in Band 1 (546–604 nm). The maximum normalized intensity values, 5 × 10⁻⁹ for Band 1 HLC (with error <3%) and 3 × 10⁻⁹ for Band 3 SPC (with error 5%–10% over the spectral range), are within the detection limits of Roman CGI under current best estimates of instrument performance (Bailey et al. 2019).

Zoom In Zoom Out Reset image size

We simulated Roman CGI observations of υ And d using two software tools: (i) the PSG and (ii) a purpose-built Python script for simulating CGI SPC zero-order deviation (ZOD) prism spectroscopy data called zodprism (Groff et al. 2021; Poberezhskiy et al. 2021). Both tools indicate that υ And d is a viable target for Roman CGI, under the assumptions of orbit parameters and albedo given in Table 2. The right-hand plot of Figure 13 shows one simulated PSG scene assuming a peak intensity observing window, using the config file included in Appendix C. When a noise model corresponding to the Roman CGI SPC spectroscopy mode is included, PSG predicts reaching an S/N of 10 in the central wavelength bin of the R = 50 spectrum in approximately 400 hr. This time-to-S/N estimate was cross-checked with the instrument team's zodprism script, which incorporates a model of the coronagraph point-spread function and wavelength-resolved residual starlight pattern (Krist et al. 2016), the slit mask and ZOD prism dispersion profile for the spectroscopy mode, the system throughput, and the detector noise parameters (Groff et al. 2020). The integration time needed to reach the same S/N as the Band 1 HLC imaging mode would be at least 10× less, due to the combination of higher source intensity (Figure 13) and wider signal bandpass.