1.1.

Roman Mission and Exoplanet Coronagraph

The 2.4-m Nancy Grace Roman Space Telescope (Roman; formerly the Wide Field Infrared Survey Telescope, WFIRST)¹ is a NASA mission set to launch in the mid-2020s and operate at the second Sun-Earth Lagrange point (SEL2). Roman will carry two scientific instruments: (1) the Wide Field Instrument, providing imaging resolution comparable to that achieved by the Hubble Space Telescope but with 100 times the field of view and (2) a high-contrast, small field of view Coronagraph Instrument (CGI). The Roman CGI is intended to serve as a pathfinder instrument for future exoplanet imaging missions such as HabEx² by demonstrating key starlight suppression technologies, including the first active wavefront control in space.³ CGI aims to achieve detection limits of ${10}^{-7}$ (minimum requirement) to ${10}^{-8}$ (current best estimate) at separations of $\sim 0.15$ to 1.5 arcsec from bright ($V<6$) host stars, which would enable imaging and low-resolution spectroscopy of giant planets in reflected starlight.⁴ The simulations in the Roman Exoplanet Imaging Data Challenge (Roman EIDC) employed the most recent end-to-end observing scenario models (OS6),⁵ which enabled detections down to $\sim {10}^{-9}$. As of this writing, the Roman mission has passed its preliminary design review and entered phase C (final design and fabrication).

1.2.

Goals and Scope of the Data Challenge

The goals of the Roman EIDC⁶ are as follows:

The scope of the Roman EIDC was to unveil an exoplanetary system hidden in realistic CGI data, which includes high fidelity wavefront control residuals, a realistic EMCCD detector model, and astrophysical contamination. Participating teams were given six simulated imaging epochs for a nearby star (modeled after 47 UMa), calibration files corresponding to OS6, 15 years of simulated precursor radial velocity (RV) data (described in Sec. 2.3), and basic stellar system information (V-band magnitude, distance, mass, spectral type, and proper motion). Participants were not told the number of planets present nor any prior information about background stars, galaxies, or exozodiacal light. From these data, participants were asked to determine the number of planets detected and to find their orbital parameters, masses, radii, and albedos. Our three-planet 47 UMa-analog system was designed so that the planets (b, c, and d) would be in stable orbits and fall in or out of the field of view at various detection levels in the various epochs. In only one instance was the entire planetary system (marginally) detectable in a single epoch. The imaging and RV simulations will be discussed in detail in a future publication.⁷

1.3.

Challenge in Four Steps

The Roman EIDC was rolled out in four increasingly difficult steps to help participants work in a logical manner and avoid overwhelming the teams. Participants submitted the results of each step to receive the next package of simulated data, as follows:

Importantly, the SS images were not provided until step 3 because these data would have undermined the “blind” analyses of the prior epochs by clearly revealing planets that are only marginally detected with the HLC. Figures 1 and 3 show the science grade data products (post-PSF subtraction) for the HLC and SS, respectively.

Fig. 1

Simulated CGI HLC data after applying RDI and roll combination to the co-added images at two epochs, displayed in units of integrated photoelectrons. The faint source corresponding to planet “d” is circled in green. The data analysis presented in this paper was carried out “in house” by the Roman CGI Turnbull Science Investigation Team members.

1.4.

Hack-a-Thons and a Legacy Tutorial Suite

To engage the community prior to launching the Roman EIDC, we held four tutorial “hack-a-thon” events in Baltimore, Los Angeles, New York, and Tokyo. Hack-a-thon participants were given two rehearsal data sets as well as a suite of Jupyter notebooks (mainly written in Python) to begin working with the data. This training material was developed and improved over the duration of the challenge, and it is now available for public use.⁸ Altogether, a diverse group of over 70 people attended the hack-a-thons, a subset of whom ultimately participated in the full data challenge either individually or in teams. In an upcoming paper,⁹ we will report on the Roman EIDC results in detail.

2.1.

HLC Images

The baseline imaging mode of CGI is a Hybrid Lyot Coronagraph (HLC) operating in a bandpass centered at 575 nm (the 546 to 604 nm “band 1” filter of CGI). This configuration provides the smallest inner working angle (IWA) of the available coronagraph modes, $\sim 150\mathrm{mas}$, which results in the best overall sensitivity for detecting exoplanets in reflected starlight. The data challenge files include four epochs of simulated HLC images of the scientific target at time intervals $\mathrm{\Delta }T=0.0$, 0.15, 1.0, and 2.0 years from an initial observation on November 1, 2026. The time-varying residual starlight pattern, including speckles and pointing jitter effects, is based on the “Observing Scenario 6” (OS6) PSF time series prepared by the CGI project’s integrated modeling team in 2018.⁵ The OS6 PSF time series includes alternating observatory roll angles to assist in discriminating astrophysical signals from quasistatic speckles. Therefore, the simulated HLC images are provided as co-added images with a total integration time of 66,000 s ($\sim 18\mathrm{h}$) in each of two roll angles 26 deg apart, along with the science target observation. None of the provided images have been postprocessed, and it was left to the participants to apply either reference differential imaging (RDI) or angular differential imaging to partly subtract the residual starlight pattern.

The result of our “in-house” image reduction applied to two simulated HLC data epochs is shown in Fig. 1. There are two relatively bright point sources: one is a planet (NE quadrant, near the occulting mask) and the other is a background star (SW quadrant, near the edge of the field). In addition, the reduced images reveal diffuse scattered light from a debris disk with a peak surface brightness of $19.5\mathrm{mag}{\mathrm{arcsec}}^{-2}$ in $V$-band (equivalent to 10 “zodi” surface brightness units). Lastly, the source that we highlight for this article is a faint outer planet, circled in green in Fig. 1. This planet “d” source would be indistinguishable from speckle noise if these HLC images were the only evidence for its existence. However, planet d reappears as a much higher signal-to-noise ratio (SNR) point source in the later SS rendezvous images (Sec. 2.2), with motion consistent with a Keplerian orbit and physically consistent with a weak, long-period residual feature in the data challenge RVs (Sec. 2.3).

2.2.

Starshade Images

The SS images were generated with SISTER,¹⁰ an open-source software capable of performing SS simulations with high optical fidelity. We chose the 425- to 552-nm band of the Starshade Rendezvous Probe (SRP),¹¹ see Fig. 2, because it is the closest channel to the HLC simulated images. The simulation uses nominal instrument parameters consistent with Roman end-of-life (EOL) operation and includes epochs consistent with the required solar exclusion angles of the mission.

Fig. 2

The geometry of the SS and Roman telescope in the SRP mission concept.¹¹ The SS’s IWA is the apparent size of the SS’s radius as seen by the Roman telescope. In the 425 to 552 nm passband, the IWA is 72 mas.

2.2.5.

Astrophysical components

The astrophysical components of the SS simulations include the three-planet system itself, plus residual starlight from an imperfect SS, solar glint due to sun’s light scattered through the SS petals, local zodiacal light, exozodiacal dust light, and a background galaxy. The planets’ apparent separations for epoch $\mathrm{\Delta }T=3$ years are 116.6, 271.0, and 488.2 mas, respectively, and for epoch $\mathrm{\Delta }T=4$ years they are 159.6, 182.7, and 502.1 mas. In all cases, the projected distance is significantly larger than the geometric IWA of the SS in this configuration (72 mas).

We assumed a Lambertian phase function to derive the intensity of the reflected light from the planets, which is good enough to fulfill the aims of the data challenge. The flux ratios for epoch $\mathrm{\Delta }T=3$ years are $5\times {10}^{-10}$, $5.6\times {10}^{-10}$, and $6.7\times {10}^{-10}$, for planets b, c, and d, respectively, and for epoch $\mathrm{\Delta }T=4$ years they are $2.4\times {10}^{-9}$, $1.2\times {10}^{-10}$, and $4.7\times {10}^{-10}$, respectively. For comparison, Earth’s flux ratio at quadrature is $1.1\times {10}^{-10}$ at 500 nm. The planet orbital configuration was chosen to provide low flux ratios during the SS epochs to better demonstrate the capability of SS imaging.

The exozodiacal cloud is assumed to have a density five times that of the solar system’s zodiacal cloud. Its dust distribution was derived from an N-body dynamical simulation of dust particles under the gravitational influence of a co-orbiting planet.¹³ The dust grains in this model tend to become trapped in mean motion resonances with a planet. In the data challenge images, the scattered light from these resonant dust structures produces diffuse intensity enhancements that co-orbit with the innermost planet “b.” For the extragalactic background object, we used the Haystacks data cubes¹⁴ and selected the brightest background galaxy from the HST deep field survey. In Sec. 3, we provide the relative contribution of each component for the case of planet d.

Simulated images for the two epochs with the SS without any postprocessing are shown in Fig. 3. Both images show planet d clearly inside the dotted circle. In epoch $\mathrm{\Delta }T=3$ years, planet c is located near $(-0.10,-0.25)\mathrm{arcsec}$ and can be seen as a faint perturbation to the exozodiacal dust, while planet b is masked by a bright clump of the exozodiacal light. In epoch $\mathrm{\Delta }T=4$ years, planet b is clearly visible near (0.10, 0.15) arcsec, while planet c is masked by the exozodiacal light. The images also show how the resonant structure of the exozodiacal revolves around 47 UMa from one epoch to the following one and its impact on planet detection.

Fig. 3

Co-added Roman CGI/SS images of the data challenge target, without any postprocessing, displayed in units of integrated photoelectrons. The point source corresponding to planet “d” is circled in green.

2.3.

Precursor Radial Velocity Data

Synthetic RV data points were created spanning 15 years: first 10 years with RV instrumental precision of $1\mathrm{m}/\mathrm{s}$ similar to instruments such as Keck HIRES and an additional 5 years with $0.3\mathrm{m}/\mathrm{s}$ taking into account the emerging extreme precision RV instruments such as the WIYN NEID spectrograph. A total of 200 data points were created by solving Kepler’s equation and superposing individual signals from each planet. The number of data points and the length of coverage is similar to the archival RV coverage for 47 UMa, except that the data points are spread out nearly evenly (i.e., in the most optimal way). No planet–planet interactions are included in the RV data, since such effects are not observable with RV.

To mimic different sources of noise, each RV data point was passed through a Gaussian filter to randomize the velocity measurements, with the sigma calculated from the quadrature sum of instrumental error and stellar jitter. The jitter value of around $2.1\mathrm{m}/\mathrm{s}$ was estimated using the method from Isaacson and Fischer¹⁵ assuming the stellar mass and $S$-index values of a 47 UMa host star. The total RV as well as the individual phase-folded RV for each planet is shown below. The figure is generated using RV modeling toolkit RadVel ¹⁶ with all of the parameters fixed so that the outermost planet can be properly displayed. Parameters for each planet such as the orbital period (P), semiamplitude (K), and eccentricity (e) are shown in the upper right corner of each phase-folded panel.

3.1.

Orbit Fitting

By combining the two epochs of SS astrometry (Fig. 3) with the residual sinusoidal feature in the RVs [Fig. 4(c)], it is possible to place some preliminary constraints on the Keplerian orbit. In our orbit fitting trials with the open-source RadVel ¹⁶ and orbitize! ¹⁸ packages, we found that, when applying priors on period, eccentricity, and ascending node position angle with the two epochs of SS astrometry, orbitize! constrains the semimajor axis to a 68% confidence interval of 7.4 to 10.5 AU and the inclination angle to a 68% confidence interval of 41 to 56 deg.

With the addition of the weaker point source detections in the HLC images, the fit is improved by the increased time baseline. In Table 1, we list the 68% confidence intervals of the inclination angle and planet mass retrieved from the combined orbit fit, alongside the true values used to generate the simulated data. The planet mass constraints are modest, being limited by the quality of the RV signal and the fact that the imaging astrometry only covers one-fifth of planet d’s 21-year orbital period. Even with all four epochs of astrometry, the confidence interval spans more than a factor of 2 in mass, from 0.24 Jupiter masses ($M_{\text{Jup}}$) to 0.55 $M_{\text{Jup}}$. The other two data challenge planets are situated at smaller orbital radii, with semimajor axes 2.3 AU and 4.2 AU. In those cases, the narrower confidence intervals of the RV semiamplitudes, along with larger orbital period coverage of the imaging data, reduce the final mass uncertainties to within 20% of their true values. We will describe the results for these other data challenge planets, with emphasis on the participants’ analysis, in forthcoming publications. Because those inner planets are brighter and have shorter orbital periods, the HLS data are sufficient for detecting and estimating their complete orbital parameters to similar precision, without SS observations.

Table 1

Estimated orbital inclination and mass of planet d, based on combining RV priors with various subsets of astrometry. These measurements are from an “in-house” analysis of the simulated data; results from the data challenge participants are deferred to forthcoming publications.

3.2.

Photometry and Phase Curve Analysis

One of the aspects of future reflected starlight direct imaging will be the ability to observe the photometric variation over the course of the orbit due to the changing star-planet-observer phase angle. Combined knowledge of the orbit, along with an inferred planet radius, will enable observers to estimate the bulk geometric albedo from a flux ratio time series. As a demonstration, we used our imaging photometry of planet d, in combination with the previous orbit fit, to estimate the planet’s geometric albedo, assuming a Lambertian sphere scatterer.

Before fitting for the planet’s albedo, we first used the posterior mass estimate from the orbit fit and RV semiamplitude and then inferred the planet’s radius from an empirical mass–radius relationship along with the approximate equilibrium temperature of planet d.¹⁹ A simple least-squares optimization finds the albedo that minimizes the distance between the four epochs of photometry and the Lambertian phase curve associated with the orbit posteriors. The result, in Fig. 6, illustrates the four photometric data points, the best-fit flux curve, and the flux curve of the true planet parameters.

Fig. 6

Best fit (red) to the geometric albedo of planet d using photometry and orbit posteriors, assuming a Lambertian sphere scatterer. The true flux curve is shown in blue.

One of the dramatic features shown in Fig. 6 is the improved photometric constraints enabled by the higher SNR SS detections. A slight bias in the SS photometry causes the albedo to be underestimated (0.3 versus the true value of 0.5), but this bias is still much smaller than the overall uncertainty of the albedo estimate. The fit shown in Fig. 6 corresponds only to the best estimated orbit parameters and mass, not the envelope of all photometric phase curves within the confidence intervals of all orbital elements. If we repeat the fit at the minimum and maximum extremes of the 68% confidence intervals of all orbit parameters derived from the data, then the albedo confidence interval spans 0.04 to 0.73. A longer time series with photometry at the SS precision would enable better constraints.