2.1.

Planet Populations

To calculate completeness for each target for a specific instrument, we first generate a joint probability density function of $\mathrm{\Delta }\mathrm{mag}$ and planet–star separation projected onto the image plane, $s$, using Monte Carlo as in Ref. 3 for an assumed planet population. We use two planet populations in this paper; one derived from Kepler’s detections from Q1-Q6 data¹⁹ (Kepler-like)²⁰ and another derived from NASA’s Exoplanet Program Analysis Group (ExoPAG) Study Analysis Group 13 (SAG13).^21,22 The Kepler-like and SAG13 populations both share the same eccentricity and albedo distributions but differ in their occurrence rates, semimajor axis, and planetary radii distributions. The final output necessary for calculating completeness is the joint probability distribution of $\mathrm{\Delta }\mathrm{mag}$ and $s$ which is found by sampling the respective planet population and fitting a rectilinear bivariate spline. We use the same sampling methods and distributions for both calculating the joint probability distributions and generating planets in our simulated universes.

2.1.1.

Kepler-like planetary semimajor axis and radii

For the Kepler-like planet population, we adopt a power-law distribution for semimajor axis ($a$) modified from the power law fit in Ref. 23 to be of the form

Eq. (1)

$$f_{\overline{a}}(a)=\frac{a^{-0.62}}{a_{\text{norm}}}\exp (-\frac{a^2}{a_{\text{knee}}^2}),$$

where $-0.62$ is adopted from Ref. 24. In this model, we include an exponential decline in semimajor axis past a “semimajor axis knee” ($a_{\text{knee}}$). This “knee” is motivated by a lack of giant planet discoveries by the Gemini Planet Imager, which indicates the RV-based power-law planet occurrence rate from Ref. 23 approaches 0 between 10 and 100 AU.²⁵ Using giant planet radial velocity data, Ref. 26 fit many period break points to the giant planet occurrence rate power-law model. We assume an intermediate knee equivalent to $a_{\text{knee}}=10$ AU. The normalization factor is given by integrating the non-normalized distribution over a specific $a$ range

Eq. (2)

$$a_{\text{norm}}={\int }_{a_{\min }}^{a_{\max }}{a}^{-0.62}\exp (-\frac{a^2}{a_{\text{knee}}^2})\mathrm{d}a,$$

where we consider values of $a$ range in $a_{\min }=0.1$ AU to $a_{\max }=30$ AU, again based on the paucity of wide-separation planets discovered to date. The per-simulation average distribution of $\mathbf{a}$ is shown in the histogram above Fig. 1(a). We note that $a_{\min }<\min (s_{\min ,i})$ for WFIRST and this planet population. Since the closest target list star has distance $\min (d_i)=2.63$ pc and has the smallest observable planet star separation $\min (s_{\min ,i})\approx \mathrm{IWA}\times \min (d_i)$, then WFIRST’s inner working angle (IWA) of 0.15 arcsec means the smallest planet–star separation observable by WFIRST is 0.394 AU. By assuming a maximum planet eccentricity $(e_{\max })=0.35$, the smallest observable semimajor axis $a_{\min }$ from $(s_{\min })\approx a_{\min }(1+e_{\max })=0.292\mathrm{AU}$.

Fig. 1

The universe of planets generated over all simulations for (a) Kepler-like and (b) SAG13 planet populations. (c) and (d) The populations of detected planets for simulations run with completeness calculated using the Kepler-like planet population observing a universe of (a) the Kepler-like planets and (b) SAG13 universe of planets. (e) and (f) The populations of detected planets for simulations run with completeness calculated using the SAG13 planet population observing a universe of (a) the Kepler-like planets and (b) SAG13 universe of planets. Overlay text on (a) and (b) shows average occurrences per grid-space averaged over the Monte Carlo of simulations. Overlay text on (e) and (f) shows average detections per grid-space averaged over the Monte Carlo of simulations. These plots reference RpvsSMAdetectionsDATA_WFIRSTcycle6core_CKL2_PPKL2_2019_04_05_19_34_.txt.

We base our Kepler-like planetary radii ($R$) based on Fig. 7 in Ref. 19. We define the bin counts in Fig. 7 of Ref. 19 as ${\mathbf{R}}_{85}$. These bins range from $R_{\min }=1R_{\oplus }$ to $R_{\max }=22.6R_{\oplus }$ and are based solely on data from planets with periods $>0.8$ days, but $<85$ days. To properly tune the overall Kepler-like planet occurrence rates (${\eta }_{\mathrm{KL}}$) starting with ${\mathbf{R}}_{85}$, we first use Eq. (2) evaluated at $a_{\min }=a_{0.8}$ to $a_{\max }=a_{85}$ to get $a_{85,\text{norm}}$. Here, $a_{0.8}$ and $a_{85}$ are the semimajor axes corresponding to periods of 0.8 and 85 days around a Sun mass star. We then scale ${\mathbf{R}}_{85}$ by Eq. (2) and $a_{85,\text{norm}}$ to get the adjusted planetary radii occurrence rates:

Eq. (3)

$${\mathbf{R}}_{\mathrm{vals}}={\mathbf{R}}_{85}\frac{a_{\text{norm}}}{a_{85,\text{norm}}}.$$

We multiply the last five bins in ${\mathbf{R}}_{\mathrm{vals}}$ by 2.5 to account for longer orbital baseline data and more closely match the larger period orbit distributions available from radial velocity surveys at the time when this distribution was first derived.^23,27 Even with this multiplication, the right-hand side histogram of Fig. 1(a) shows the characteristic drop-off in planetary radius clearly observable in Fig. 7 of Ref. 19. The Kepler-like planet occurrence rate per star is ${\eta }_{\mathrm{KL}}=\sum {\mathbf{R}}_{\mathrm{vals}}=2.375$ over the specific $a$ and $R$ ranges discussed in this section.

When we simulate a universe of Kepler-like planets, we first calculate the number of planets to generate for each planet radius bin ($N_q$) for all bins in ${\mathbf{R}}_{\mathrm{vals}}$. This is given as

Eq. (4)

$$N_q=N_p\lceil \frac{R_{\mathrm{vals},q}}{\sum {\mathbf{R}}_{\mathrm{vals}}}\rceil .$$

In this paper, $N_p$ is the number of planets to generate. In our Monte Carlo completeness calculations, we generate ${10}^8$ planets. We take $N_q$ samples from a log-uniform distribution over the $q$’th bin range in ${\mathbf{R}}_{\mathrm{vals}}$ to get a set of planets radii (${\mathbf{R}}_q$).¹⁹ Finally, we randomly select $N_p$ radii from the collective nine sets of planetary radii to get $\mathbf{R}={\bigcup }_{q=1}^9{\mathbf{R}}_q$. At the same time, we use an inverse transform sampler to generate $\mathbf{a}=\{a_k\forall k\in 1..N_p\}$ based on Eq. (1).

2.1.2.

SAG13 planetary radii and semimajor axis

The SAG13 planet population represents a significant update on the Kepler-like population, incorporating all available Kepler data circa 2017. The occurrence rate model presented in Ref. 28 is a power-law model fit as a function of $\ln P$ and $\ln R$ but does not include an occurrence rate turnover indicated in giant planet occurrence rates from radial velocity studies²⁶ and the lack of detections of large $a$ planets by subsequent Gemini Planet Imager survey results.²⁵ We need an occurrence rate model that is in terms of linear Keplerian orbital elements and decreases occurrence rates for large $a$ planets so our extrapolation of the model beyond the bounds in Ref. 28 does not produce substantially more large $a$ planets than expected. We include a full derivation of the SAG13 occurrence rate model implemented^22,29 in EXOSIMS in Sec. 7.8 but summarize the steps, important components, and results here.

We start with the original occurrence rate model fit from Ref. 28 included as Eq. (37) and convert from $(\ln P,\ln R)$ to $(P,R)$ resulting in Eq. (38). We reconstructed the occurrence rate grid from the linear model in Fig. 2(a). When directly compared to the occurrence rate grid for G-type stars in Ref. 28, we see a maximum deviation of 0.74 ($100\times$ planets per bin occurrence) and maximum percent deviation of 13%. We then convert Eq. (38) from period ($P$) space to $a$ space using Eqs. (39) and (40) assuming a solar mass star to arrive at Eq. (41). We compare the analytical integral over $a$ space of this new occurrence rate model to arrive at Fig. 2(b) which shows identical occurrence rates in each of their counterparts’ grid spaces in Fig. 2(a). We now append a semimajor axis “knee” to decay occurrence rates of larger $a$ planets to 0, which results in Eq. (42). We calculate the overall SAG13 exoplanet occurrence rate per planet (${\eta }_{\mathrm{SAG}13}$) by performing the double integral over the entire $a$ and $R$ space as in Eq. (43). The SAG13 model parameters in Table 1 are split at $3.4R_{\oplus }$, so we use the notation $i=0$ for $R_{\min }=0.666R_{\oplus }\le R<3.4R_{\oplus }$ and $i=1$ for $3.4R_{\oplus }\le R<17.086R_{\oplus }=R_{\max }$. An intermediate simplification of the double integral in the calculation of ${\eta }_{\mathrm{SAG}13}$ is

Eq. (5)

$${\eta }_{\mathrm{SAG}13}=\frac{{\gamma }_0{K}_0}{{\alpha }_0}[R_{\max }^{{\alpha }_0}-{(3.4R_{\oplus })}^{{\alpha }_0}]+\frac{{\gamma }_1{K}_1}{{\alpha }_1}[{(3.4R_{\oplus })}^{{\alpha }_1}-R_{\min }^{{\alpha }_1}],$$

where ${\gamma }_i$, ${\beta }_i$, and ${\alpha }_i$ are constants from Table 1. The ${\alpha }_i$ in the denominator and $({\alpha }_i-1)$ to ${\alpha }_i$ exponent are a result of the indefinite integral over $R$. The four $R$ terms are the result of evaluating the indefinite integral over $R_{\min }$ to $3.4R_{\oplus }$ to $R_{\max }$. Finally, $K_i$ is an intermediate calculation representing the marginalization of Eq. (42) over $a$ written as

Eq. (6)

$$K_i={\int }_{a_{\min }}^{a_{\max }}{\left(2\pi \sqrt{\frac{a^3}{\mu }}\right)}^{({\beta }_i-1)}\left(3\pi \sqrt{\frac{a}{\mu }}\right)\exp (-\frac{a^3}{a_{{\text{knee}}^3}})\mathrm{d}a.$$

Here, $\mu$ is the gravitational parameter assuming a solar mass star. We use the same $a_{\text{knee}}=10$ AU in the SAG13 planet population that we use for the Kepler-like planet population. The average number of planets generated per star in the EXOSIMS implemented SAG13 universe over the range specified is ${\eta }_{\mathrm{SAG}13}=5.62$.

Table 1

Parametric fit parameters for the SAG13 planet population implemented in EXOSIMS.28

Parameter	[i=0,i=1]
${\beta }_i$	[0.26, 0.59]
${\alpha }_i$	$[-0.19,-1.18]$
${\gamma }_i$	[0.38, 0.73]

Fig. 2

We calculate the SAG 13 exoplanet occurrence rate cumulative distribution by integrating over the linear occurrence rate model derived from Ref. 28 and included in Eq. (38) over $R$ and $P$ of each bin in (a). The purple numbers in (a) represent $100\times$ the per-bin per star exoplanet occurrence rate and can be directly compared to the planet occurrence rate grid per star in Ref. 28, which shows good agreement with the largest deviation between implementations of 0.74%. The total SAG 13 occurrence rate per star found using Eqs. (38) and (43) over the plot range in (a) is 1.95 planets per star. The black percentages in (a) are calculated by dividing the $100\times$ per bin occurrence (purple number) by the occurrence rate over the entire grid range and represents the probability a single planet will be in that bin. The exoplanet occurrence rate cumulative distribution in (b) is found using Eqs. (41) and (43) and integrating over the $R$ and $a$ ranges in (b). The $P$ range in (a) and the $a$ range in (b) are mapped using Eq. (39) assuming a Sun mass star. The $100\times$ and percentage occurrence rates per bin between (a) and (b) are identical. We show the joint probability distribution of $R$ and $a$ from Eq. (44) implemented in EXOSIMS in (c) which is extrapolated beyond the range of (a) and includes the semimajor axis “knee.” The apparent difference in color gradient between (a) or (b) and (c) is due to the logarithmic scale of (c); the top right bins are orders of magnitude larger in area than the bottom left bins.

In order to generate planets with $R$ and $a$ from the occurrence rate model, we still need a joint probability density function $f_{{\overline{R}}_p,\overline{a}}(R,a)$ as well as a single variable probability density function and conditional probability density function. By normalizing Eq. (42) by ${\eta }_{\mathrm{SAG}13}$, we get the joint probability distribution in Eq. (44). From here, we calculate the probability distribution of $R$ for the SAG13 planet population by marginalizing Eq. (44) over $a$ only. This gives us the planetary radius distribution

Eq. (7)

$$f_{\overline{R}}(R)=\frac{{\gamma }_i{K}_i}{{\eta }_{\mathrm{SAG}13}}{R}^{({\alpha }_i-1)}.$$

This probability density function inherits the parametric conditions of the joint probability density function. Since $f_{\overline{R}}(R)$ is a marginalization over the joint probability density function, the integral of $f_{\overline{R}}(R)$ over the range $R_{\min }\le R\le R_{\max }$ is 1. The distribution of generated planets in $\mathbf{R}$ is shown in the right-side histogram of Fig. 1(b). We now calculate the conditional probability distribution $f_{\overline{a}|\overline{R}}(a)$ using $f_{\overline{R},\overline{a}}(R,a)$, $f_{\overline{R}}(R)$, and Bayes rule. The conditional probability distribution is a simple division of the joint probability density function by the marginalized probability density function to get

Eq. (8)

$$f_{\overline{a}|{\overline{R}}_p}(a)=\frac{1}{K_i}{\left(2\pi \sqrt{\frac{a^3}{\mu }}\right)}^{({\beta }_i-1)}\left(3\pi \sqrt{\frac{a}{\mu }}\right)\exp (-\frac{a^3}{a_{\text{knee}}^3}).$$

We now have the tools to sample semimajor axis and planetary radius.

To generate a planet from the probability distributions in this section, we first sample Eq. (7) using a single variable inverse transform sampler to get a set of planetary radius, ${\mathbf{R}}_p$. We also use a single variable inverse transform sampler to sample Eq. (8) and get $\mathbf{a}$. The $a_{\min }$ and $a_{\max }$ used in the SAG13 population are the same as in the Kepler-like planet population. SAG13’s $a$ frequency distribution is included above Fig. 1(b). The two-dimensional (2-D) contour plot and associated grid values in Fig. 1(b) show $\mathbf{R}$ and $\mathbf{a}$ interdependence.

2.1.5.

Joint probability density function

The final remaining parameters needed to describe planets at the time of an observation are the physical location of the planets relative to their host stars (${\underset{\_}{r}}_{k/i}$), which is defined by Eq. (1) in Ref. 31. We assume that the direction of the orbit eccentricity vector is uniformly distributed in space, so that the orbit inclination is sinusoidally distributed between 0 and $\pi$, while the remaining Keplerian elements (longitude of the ascending node, argument of periapsis, and mean anomaly) are all uniformly distributed from 0 to $2\pi$.³¹

We sample all of the parameters aforementioned for each planet which is sufficient information to calculate a projected planet–star separation $(s_k)$ and difference in magnitude between the host star and the planet $(\mathrm{\Delta }\mathrm{mag})$. We calculate the projected planet–star separation as

Eq. (10)

$$s_k=|{\underset{\_}{r}}_{k/i}-{\underset{\_}{r}}_{k/i}·\frac{{\underset{\_}{r}}_{i/SC}}{|{\underset{\_}{r}}_{i/SC}|}|.$$

We calculate $\mathrm{\Delta }{\mathrm{mag}}_k$ using Eq. (4) from Ref. 7 and adopt their use of a Lambert phase function.

Sampling all of the parameters described above for $N_p={10}^8$ planets allows us to calculate the joint probability density of projected separation and star–planet magnitude difference, $f_{\overline{s},\overline{\mathrm{\Delta }\mathrm{mag}}}(s,\mathrm{\Delta }\mathrm{mag})$. As in Eq. (7) in Ref. 7, we assume independence between all parameters, except for semimajor axis and planet radius in the SAG13 population as well as the planet albedo dependence on planet radius. We bin these planets into a $1000\times 1000$ grid and fit a rectilinear bivariate spline to the resulting 2-D histogram. This rectilinear bivariate spline consists of integrable high order polynomials with total volume under the surface of 1. The resulting $f_{\overline{s},\overline{\mathrm{\Delta }\mathrm{mag}}}(s,\mathrm{\Delta }\mathrm{mag})$ densities for the two populations are shown in Fig. 3.

Fig. 3

Joint probability density functions of projected separation and $\mathrm{\Delta }\mathrm{mag}$ based on (a) Kepler-like and (b) SAG13 planet populations. The Kepler-like distribution produces larger orbital radii (and therefore projected separations) than SAG13 for the same $a_{\text{knee}}$ values due to the use of the quadratic and cubic semimajor axis rollover functions [see Eqs. (1) and (6)].

2.2.

Star Catalog

We need a list of target stars, along with their positions on the sky, distance ($d_i$), and apparent brightness in B and V bands for our calculations of completeness and integration time. We derive the list of targets stars from the EXOCAT-1 star catalog discussed in Ref. 32, which contains a variety of targets out to a distance of $\sim 30\mathrm{pc}$. Some of these stars are missing photometric values, including the luminosity, absolute magnitude, V band bolometric correction, and the apparent VBHJK magnitudes, all of which we attempt to fill in by interpolating a table of standard stars by spectral type in Ref. 33.

From the $N=2396$ targets, we have an initial set of target stars $\mathbf{I}$ which we pare down through a series of filters. There are many stars in the EXOCAT-1 star catalog with “not a number” data entries which can propagate through our equations so we remove targets with these kinds of data entries. Not all of these entries are absolutely necessary so some targets may have been unnecessarily filtered. Some of these star systems are binary systems which we do not account for in our equations so we omit them from $\mathbf{I}$. In addition, some stars have low apparent visual magnitudes so we filter any target stars we know would take excessively long to achieve a reasonable $\mathrm{\Delta }\mathrm{mag}$ on.

Missing data filter removes 429 targets with any fields containing a “not a number” value within the set of IPAC fields {hip_name, st_spttype, parx, st_vmag, st_j2m, st_h2m, st_vmk, st_dist, st_bmv, st_mbol, st_lbol, coords, st_pmra, st_pmdec, wds_sep} as well as whenever we could not fill in a missing photometric value of a star.

Binary star filter removes 164 targets using the Washington Double Star catalog (filters targets with companion stars within $<10\mathrm{arcsec}$).³⁴

Integration time cut-off filter removes 1436 targets with integration times $t_i>30\text{days}$, where $t_i$ is calculated assuming local zodiacal light to be $Z=23.0$ mag $\mathrm{arc}{\sec }^{-2}$, exozodiacal light with magnitude $EZ=22.0$ mag ${\mathrm{arcsec}}^{-2}$, $\mathrm{\Delta }{\mathrm{mag}}_0=22.5$ [used on the $t_i(\mathrm{\Delta }\mathrm{mag})$ equations in Refs. 4 and 35], and a working angle ${\mathrm{WA}}_0=0.28\mathrm{arcsec}$ (see Sec. 7). This WA results in reasonable core mean intensity and throughout combination based on inspection of the instrument parameter files in Figs. 4(a) and 4(c).

Fig. 4

System throughput in the FWHM region of the planet PSF core (a) from /WFIRST_cycle6/G22_FIT_565/G22_FIT_565_thruput.fits. Intensity transmission of extended background sources such as zodiacal light including the pupil mask, Lyot stop, and polarizer (b) from /WFIRST_cycle6/G22_FIT_565/G22_FIT_565_occ_trans.fits. Core mean intensity as a function of wavelength and working angle. Black lines and dot represent detection mode wavelength (565 nm) and IWA (0.3 arcsec) for integration time calculation (c) from /WFIRST_cycle6/G22_FIT_565/G22_FIT_565_mean_intensity.fits. Area of the FWHM region of the planet PSF in ${\mathrm{arcsec}}^2$ (d) from /WFIRST_cycle6/G22_FIT_565/G22_FIT_565_area.fits.

After paring down $\mathbf{I}$ using a missing data filter, binary star filter, and integration time cut-off filter, $\mathbf{I}$ is reduced down to 651 targets (the filters are not mutually exclusive). Our initial filtering of target stars reduces the number of degrees of freedom in the optimization process and increases computation time. Larger telescopes will result in larger target lists with less candidate stars removed via the initial integration time filter. The statically optimized target list is included in Table 7. This optimized target list contains mainly FGK-type stars and 8 A-type stars. All of these targets are located within 20 pc.

2.3.

Calculating Completeness

Completeness is calculated for the $i$’th target by integrating over the joint probability density function of $s$ and $\mathrm{\Delta }\mathrm{mag}$¹⁵

Fig. 5

Quantum efficiency of the detector as a function of wavelength. White dot represents detection mode wavelength (565 nm). From /WFIRST_cycle6/QE/Basic_Multi2_std_Si.fits

The theoretical maximum completeness for the $i$’th target ($c_{\infty ,i}$) is found by integrating Eq. (11) to the upper limit of $\mathrm{\Delta }{\mathrm{mag}}_i$ at $t_{\infty }$ (this assumes an infinite observing time is available). Using the WFIRST parameters, we see a universal upper limiting $\mathrm{\Delta }{\mathrm{mag}}_i$ of 23.137. It is important to note that the inclusion of speckle residuals means that integrating past a certain point will not produce any improvements in the achieved SNR, meaning there is a specific time for every target past which it makes no sense to integrate further.

While the equations derived thus far are sufficient to perform continuous integration time optimization, gradient-based solvers, such as the ones we employ as follows, all perform significantly better with analytical gradient functions. To expedite the optimization process in Sec. 2.4, we calculate the derivative of completeness with respect to integration time as a function of integration time. As in Ref. 15, the derivative of Eq. (12) with respect to time is

2.4.

Optimization Process

Our goal is to maximize the number of unique detections throughout a blind-search survey. We evaluate the reward potential of a target using the completeness metric and the cost as the integration time required to achieve that completeness. In our optimization formulation, we seek to maximize the summed completeness that fits within the rigid mission time constraints coupled with an additional time cost per observation. Maximizing summed completeness with fixed overhead per observation in a time constrained mission is the full nonlinear optimization problem we aim to solve in this section.

Our optimization process is broken down into three major steps:

Both step 1 and step 2 enable the solving of the full nonlinear optimization problem in step 3. Step 3 can be seeded with any “good” initial guess (an initial solution sufficiently close to an optimal solution). A “good” initial guess for SLSQP can even be an infeasible solution (i.e., in the case of an unexpected reduction in total observing time, we can seed step 3 with the previous solution to step 3). Optimization step 3 is fundamentally important because it allows individual adjustment of each $t_i$, individual removal of targets, and individual addition of targets (something Ref. 4 cannot do) while conforming to all time-related constraints of the full nonlinear optimization problem.

2.4.4.

Step 3: SLSQP minimization

In step 3, we formulate the SLSQP optimization process with an initial solution seeded with ${\mathbf{x}}_1^{*}{\mathbf{t}}_0$ or ${\mathbf{x}}_2^{*}{\mathbf{t}}_2^{*}$, whichever produces a larger summed completeness. In practice, we could take any $\mathbf{xt}$ as a good initial guess and resolve with new mission time constraints based on new information. This seeded solution should prove sufficiently close to an optimal solution such that the $\mathbf{c}(\mathbf{t})$ sigmoid-like inflection point is exceeded. We now replace our previous assumed $f_{Z0}$ with a more optimistic surface brightness. We calculate zodiacal light surface brightness every $\approx 1/3$ of a day for a year, interpolating the lookup tables from Ref. 37 as shown in Fig. 6. For our optimization, we specifically use the per annum minimum, ${\mathbf{f}}_{Z,\min }$ for all targets excluding $f_{Z,i}$ in keep-out regions shown in Fig. 7(a). This is crucial as the $f_{Z,i}$ intensity of targets with 0 deg heliocentric ecliptic latitude, ($b$), has local zodiacal light minimum within 56 deg of the observatory’s antisolar point (see Fig. 6), a region not visible to WFIRST due to solar panel pointing requirements in Table 3 (Algorithm 3).

Fig. 6

Local zodiacal light intensity interpolant [$f_Z(l,b,\lambda )$ of Eq. (19)] at target star-spacecraft heliocentric ecliptic longitude [$l({\widehat{r}}_{i/SC})$] relative to the Sun-spacecraft heliocentric ecliptic longitude [$l({\widehat{r}}_{\odot /SC})$] and target star heliocentric ecliptic latitude [$b({\widehat{r}}_{i/SC})$] relative to spacecraft heliocentric ecliptic latitude [$b({\widehat{r}}_{\odot /SC})$] and $\lambda =565\mathrm{nm}$. The minimum zodiacal light intensity [$f_{Z,\min }(b)$] for each $b$ is indicated by densely packed red squares, which appear to form lines caused using a linear interpolant and the coarseness of the underlying datapoints. Our implementation makes observations only when target stars are coincident with black squares, the minimum $f_Z$ with keep-out restrictions ($b<37\mathrm{deg}$). Note Ref. 37 is significant to three decimal places but the interpolant introduces machine precision numbers, the red and black dots have functionally equivalent values for $b>37\mathrm{deg}$. Only one quadrant is shown as the data are assumed to be reflection symmetric.

Fig. 7

Keep-out map of WFIRST/targets over the first mission year showing visible times of targets (white) and different sources of keep-out occlusion including the Sun (yellow), Earth (blue), Moon (gray), and Mercury/Venus (red). The horizontal histogram shows percentage of time each target is visible (a). The keep-out map filter is implemented in EXOSIMS at the step 2 shown in Fig. 8. The histogram and cumulative distribution of visibility of all targets is shown (b). Minimum target visibility is 28%.^38,39

Algorithm 3

SLSQP optimization.

Input:$\mathbf{I}$, ${\mathbf{f}}_Z$, $\mathbf{t}$, $T_{\mathrm{OH}}$, $T_{\text{settling}}$, and $T_{\max }$

Output:${\mathbf{t}}_3^{*}$, the integration times to spend on each star $$\begin{gathered} {\mathbf{t}}^{*}=\underset{\mathbf{t}}{\min }-\sum _{i=0}^{N-1}{c}_i(t_i) \\ \mathrm{s.t.} \\ {t}_i<T_{\max },\forall i\in \mathbf{I} \\ -t_i<0,\forall i\in \mathbf{I} \\ \sum _{i\in \mathbf{I}}{x}_i(T_{\mathrm{OH}}+T_{\text{settling}})+t_{0,i}<T_{\max }. \end{gathered}$$

The output ${\mathbf{t}}_3^{*}$ is an optimal allocation of integration times to each target accounting for the nonlinear overhead time assignment. This optimal allocation encodes our observation priority over all possible targets and means each target in the target list is equivalently important to observe to achieve the expected summed completeness. In the validation section, we choose the minimum zodiacal light intensity target selection metric to complement our optimization assumptions. The drawback is, if variations occur in the total mission time for any reason, high reward targets might not be observed. There are technically infinite “near-optimal” solutions. Equivalent maximum summed completeness target lists can be achieved for a range of different numbers of target stars. This is directly caused by the increasing number of approximately equivalent target stars at further stellar distances. These integration times are used as inputs to our validation process discussed in Sec. 2.5 and implemented in EXOSIMS, which runs a Monte Carlo of full mission survey simulations, observing each target for the integration times prescribed in ${\mathbf{t}}_3^{*}$ and bookkeeping all dynamic aspects of the mission.

2.5.

Validation

We validate the ability to schedule our optimized integration times from Sec. 2.4 and listed in Table 7 of Sec. 8, via a Monte Carlo of full survey simulations using the EXOSIMS code base, detailed in Ref. 16. This software allows us to bookkeep dynamic aspects of the mission while scheduling observations on randomly generated planetary systems about a real set of target stars. Using EXOSIMS, we account for exoplanet, solar system planet, and observatory orbit propagation as well as recalculating look vectors, keep-out regions, and local zodiacal light noise.

At the start of each survey simulation, we randomly generate planets around stars; drawn from the Kepler-like or SAG13 planet populations discussed in Sec. 2.1 using sampling methods described in Ref. 16. To define the observable times of individual stars, we combine solar system planet locations with instrument-specific keep-out angles and observatory look vectors for each target in the star catalog throughout the duration of the mission. These solar system planet locations are taken from the Horizons Ephemeris System furnished by NASA’s Navigation and Ancillary Information Facility (NAIF)^40,41 based on an assumed mission start modified Julian date (60,634 for WFIRST). The EXOSIMS framework default observatory orbit is a quasiperiodic, stable, halo orbit at the second Earth–Sun Lagrange point with period of $\approx 180$ days based on Ref. 42. We stitch the observatory position along the orbit using a one-dimensional (1-D) spline interpolant and propagate the observatory along the interpolant throughout MET assuming a general observatory start location on the Halo orbit when the Earth is at the Earth–Sun equinox (60,575.25 MJD for WFIRST). The starting location of the observatory for a single-visit blind search survey coronagraph mission has negligible impact on yield. We discuss the calculation of keep-out regions in Sec. 2.5.1. In each simulation, we incrementally filter available targets, simulate observations and their outcomes, and propagate orbits as shown in Fig. 8.

Fig. 8

EXOSIMS survey simulation simplified flowchart depicting major filtering steps discussed in Sec. 2.2, calculating integration time, selecting the next target, making the detection observation, making the characterization observation under conditions outlined in Sec. 7, and advancing time. EXOSIMS is additionally capable of strictly adhering to predefined observing blocks, but this functionality was not used for the results presented here.

At the start of each main loop shown in Fig. 8 (steps 1 to 9), filters remove targets with too long of integration times ($t_i>30$ day); targets currently in keep-out regions of planetary bodies; previously observed stars not currently scheduled for revisits; and filters targets not observable within the nearest time constraint. In this paper, we do not revisit targets, so we filter out any targets that have been observed during the mission. We then use an intelligent method of choosing the next target at the current mission time, $t_c$. Several selection metrics have been studied in Refs. 15, 17, and 43, but we choose to observe targets at their minimum zodiacal light intensity. Here, our method of selecting targets for observation advances time by the smallest amount $\min (t_c-{\mathbf{t}}_{f_{Z,\min }})$. ${\mathbf{t}}_{f_{Z,\min }}$ are the mission times when targets in $\mathbf{I}$ have their next local zodiacal light minimum, ${\mathbf{f}}_{Z,\min }$. This identifies target star $i$ which we proceed to observe for the precalculated time $t_i$ from the method described in Sec. 2.4. We divide the observations into discrete time intervals and calculate the signal and noise at each interval, calculate the total achieved SNR, and propagate planets around the star as well as the observatory and solar system planets. Splitting up observations into time intervals enables us to approximate the achieved SNR via Riemann sum using the new location of the observatory, planets, and simulated exoplanets along with the associated zodiacal light, planet phase angle, and planet position in the instrument working angle. We divide our observations into two equal intervals over the duration of the integration time, $t_i$, which are all strictly $<2$ days in this case. We then advance the mission time by $t_i+T_{\mathrm{OH}}+T_{\text{settling}}$ and check if the mission is over.

The EXOSIMS framework relies upon probabilistic planet generation and random draws. However, our Python implementation is capable of not only replicating results but also fully reproducing each survey simulation run by resetting the simulation from the simulation’s random seed. EXOSIMS also keeps track of all inputs required to replicate the simulation.

2.5.2.

Local zodiacal light

The local zodiacal light intensity ($f_{Z,i}$) is the largest, known, time varying noise source we account for and manifests itself in the background noise count rate introduced in Eq. (12) and explicitly written in Sec. 7. Variations in local zodiacal light can vary the summed completeness by up to 10%.^17,43,44

In Sec. 2.4, for step 1 for our integration time optimization, we use a static $Z_0$ of 23.0 mag ${\mathrm{arcsec}}^{-2}$ as done in Ref. 4. However, in step 3, we optimize our final integration times using ${\mathbf{f}}_{Z,\min }$ (or equivalently ${\mathbf{Z}}_{\max }$) which implies specific times of the year when these observations can be made, as shown by the black lines in Fig. 6. The black points in Fig. 6 are curved at low $l$, and this is directly caused by the large anti-solar keep-out region which marginally inhibits observation at $f_{Z,\min }$.

There is a distinction between the $f_{Z,i}$ used in the integration time optimization and the $f_{Z,i}$ used in the Monte Carlo validation. When we evaluate whether a detection has been made at the top of step 7 in Fig. 8, we calculate the planet SNR using $f_{Z,i}$ based on ${\underset{\_}{r}}_{i/SC}(t_c)$ where $t_c$ is the current time in the simulation. We calculate the interpolant in Fig. 6 using Eq. (19), a 2-D linear interpolation of intensity from Table 17 of Ref. 37 [$f_{\beta }(l,b)$], a quadratic interpolation of wavelength dependence from Table 19 in Ref. 37 {$f_{\lambda }[{\log }_{10}(\lambda )]$} and applying a Sun color correction of ${\mathcal{F}}_0(\lambda )$ to get

Eq. (19)

$$f_Z(l,b,\lambda )=\frac{f_{\lambda }[{\log }_{10}(\lambda )]f_{\beta }(l,b)}{h{c}^{\prime }{\mathcal{F}}_0(\lambda )}.$$

Here, $h$ and $c^{\prime }$ are the Planck’s constant and speed of light in a vacuum, and $f_{\lambda }[{\log }_{10}(\lambda )]$ is a quadratic 1-D interpolant of Ref. 37 data in Fig. 9 of Sec. 7. Since the Table 17 from Ref. 37 is in the geocentric ecliptic coordinates, but the spacecraft will physically be located on an Earth–Sun L2 orbit, we assume the table’s $b=0\mathrm{deg}$ and $l=0\mathrm{deg}$ is coincident with ${\underset{\_}{r}}_{\odot /SC}$, and the additional distance of the spacecraft from the Sun has no effect. Further discussion of these two assumptions is included in Sec. 7.

Fig. 9

Local zodiacal light correction factor from Table 19 of Ref. 37 with region of high accuracy between 200 nm (red dashed line) and $2.0\mu \mathrm{m}$ (blue dashed line) and decreasing accuracy at larger wavelengths due to infrared and scattering parity in contribution.

Since the statically optimized integration time assumes a fixed zodiacal light, which coincides with a specific time of year and our actual observation takes into account the local zodiacal light conditions, our target selection method has an impact on the zodiacal light of the target we observe. In previous work, we tested different target selection metrics incorporating completeness, integration time, and deviations from the maximum or minimum zodiacal light intensity.^15,17,42 For a statically optimized target list, we found the strategy observing at minimum zodiacal light to be most effective. However, this assumes we have all the observing time allotted. If we included instrument degradation due to radiation or possibly early mission terminations, other selection metrics which preferentially select higher performing targets may become more attractive.

Using the 2-D interpolation of Fig. 6 combined with instrument specific filters in Sec. 3 and optimization process in Sec. 2.4, we can plot the histogram of minimum (black dots in Fig. 6) and maximum (at edge of solar keep-out) local zodiacal light intensity in Fig. 10. Depending on the time of year observations are made, the variation in zodiacal light intensity changes (excluding when targets are in keep-out) by upward of two magnitudes. We also see using an estimated $Z_0$ of 23.0 mag ${\mathrm{arcsec}}^{-2}$ is an overly optimistic approximation of the local zodiacal light noise. Realistically, for our filtered set of target stars, the appropriate ${\mu }_{Z_{\min }}\approx 22.79$ mag ${\mathrm{arcsec}}^{-2}$ and ${\mu }_{Z_{\max }}\approx 21.59$ mag ${\mathrm{arcsec}}^{-2}$. Since $f_{Z,\min }$ and $f_{Z,\max }$ in Fig. 10 include the instrument-specific keep-out regions for WFIRST, increasing the maximum solar keep-out could decrease $f_{Z,\min }$ and decreasing the minimum solar keep-out angle would increase $f_{Z,\max }$.

Fig. 10

A histogram of zodiacal light intensity in mag ${\mathrm{arcsec}}^{-2}$ ($Z$), for stars at the minimum observed intensity ($Z_{\max }$), and stars at maximum observed intensity ($Z_{\min }$), taking into account WFIRST keep-out regions. $Z_0$ corresponds to the optimistic zodiacal light intensity in mag ${\mathrm{arcsec}}^{-2}$ used by Stark in Ref. 4 and Brown in Ref. 3. Dashed lines represent target list ${\mu }_{Z_{\min }}=22.79$ mag ${\mathrm{arcsec}}^{-2}$ and ${\mu }_{Z_{\max }}=21.59$ mag ${\mathrm{arcsec}}^{-2}$.

2.5.4.

Convergence

We determine the required number of simulations to ensure the accuracy of our results by executing 10,000 simulations of a generic input specification similar to that used in Sec. 3. By assuming the true mean yield is equivalent to the mean yield over 10,000 simulations (${\mu }_{\det ,\mathrm{10,000}}$), we can demonstrate convergence of ${\mu }_{\det ,i}$ to ${\mu }_{\det ,\mathrm{10,000}}$, where $i$ in this section is the number of simulations in the ensemble. We calculate the mean yield from a subset of ensembles incrementally for each $i$ using Eq. (38) from Ref. 47. We then calculate the percent deviation between ${\mu }_{\det ,i}$ and ${\mu }_{\det ,\mathrm{10,000}}$ shown in Fig. 11. Figure 11 shows the mean number of unique detections converges for a sufficiently large number of simulations.

Fig. 11

Absolute percent error from ${\mu }_{{\det }_{\mathrm{10,000}}}$ convergence combined with $1\sigma$, $2\sigma$, and $3\sigma$ confidence intervals.

By assuming the mean number of detections per run are normally distributed, we show the $1\sigma$, $2\sigma$, and $3\sigma$ confidence intervals of the mean in Fig. 11. In reality, the ensembles are some form of gamma distribution because the numbers are all positive and priors are exponentially distributed. This normal distribution assumption fits better for high yield telescopes. The absolute percent errors for varying confidence intervals for 100 and 1000 simulations are presented in Table 4. Absolute percent errors represent the uncertainty in the mean number of unique detections. Excluding Table 4 and Fig. 11, all results in this paper are derived from single simulations or ensembles of 1000 simulations. We use the $3\sigma$ confidence interval to make comparisons and determine whether a result is significant. We can say a general mean number of unique detections is accurate to within $\pm 3.19\%$ at $3\sigma$ for 1000 simulations.

Table 4

Absolute percent error confidence intervals for 100 and 1000 simulations. This table references data created using runs from Dean22May18RS09CXXfZ01OB01PP01SU01 and file convergenceDATA_Dean22May18RS09CXXfZ01OB01PP01SU01_2019_04_09_01_23_.txt.

# Sims	Confidence interval	Absolute percent error (%)
1000	$1\sigma$	1.16
1000	$2\sigma$	2.33
1000	$3\sigma$	3.19
100	$1\sigma$	3.45
100	$2\sigma$	6.95
100	$3\sigma$	9.58

3.1.

Completeness and Planned Observations

By applying the optimization process from Sec. 2.4 to the target list filtered in Sec. 2.2 and assuming a Kepler-like planet population with per observation $T_{\mathrm{OH}}=T_{\text{settling}}=0.5$ day and maximum observing time of $T_{\max }=91.3125$ day (3 months), we get the integration times in Table 7 in Sec. 8. The planned summed completeness of this target list is $\sum c_i(t_{3,i})=2.35$. Since completeness is the probability of detecting planets from a population around a star, multiplying $\sum c_i$ by the population occurrence rate (${\eta }_{\mathrm{KL}}$ from Sec. 2.1) gives the expectation value of planets detected, in this case equal to 5.58 detections. We calculate the ultimate completeness¹⁰ for all 651 potential targets by evaluating $\sum c_i$ at an infinite integration time ($t_{\infty }$) to get $\sum c_i(t_{\infty })=9.01$ and an expectation value of 21.40 exoplanets detected. We note here that the theoretical maximum summed completeness, if there were no $\mathrm{\Delta }\mathrm{mag}$ or WA constraints, for all 651 target stars is 651. The summed ultimate completeness of only the targets in Table 7 is $\sum c_i$ $\forall i\in \mathbf{I}=2.99$ with an expectation value of 7.10 exoplanets detected. The ratio $2.35/9.01\times 100=26.1\%$ is a measure of planned target list yield to the maximum theoretical summed completeness observing all 651 targets with $t_{\infty }$. The ratio $2.35/2.99\times 100=78.6\%$ is a measure of the summed completeness of the planned target list to the maximum theoretical summed completeness of only the target stars observed with $t_{\infty }$. These ratios represent the fraction of planet phase space about all scheduled targets that could be probed in the assumed, finite, total mission time. We can conclude from these results that, in the short amount of time allocated for WFIRST, our optimized target list is capable of observing 78.6% of the summed ultimate completeness it could possibly gather.

There is a non-negligible difference between the planned completeness, $c_{3,i}$, and the completeness actually observed $c_{t_{\mathrm{obs}},i}$. The planned observations and Monte Carlo simulations entirely based on the Kepler-like population are shown in Fig. 12. Note the loss of one observation between the planned and actual observations, which we attribute to accumulation of machine precision errors and our strict adherence to the $T_{\max }$ upper bound on observing time as well as the allowance of characterizations in the implemented survey. The loss of an observation is evident by the single red square without an associated blue circle in Fig. 12. The $\sum c_i$ actually observed in this particular simulation of the Monte Carlo was 2.33. For each observation made in a survey simulation, observed completeness (blue circles) coincident with planned completeness (red squares) indicates each simulated observation occurs under optimal conditions for a single visit. Because we are observing targets solely at the local zodiacal light minimum and do not modify integration times, all observations made are optimal. If targets were observed at suboptimal zodiacal light levels, the $c_{t_{\mathrm{obs}},i}$ would be below the $c_{0,i}$ as shown in Ref. 17.

Fig. 12

Completeness as a function of integration times calculated using the Kepler-like distribution, EXOCAT-1 star catalog, Nemati SNR model,³⁵ and Leinert Table Zodiacal Light.³⁷ Black lines show completeness versus integration time for the 10 highest planned completeness targets, the median completeness target, and the lowest completeness target. Red squares indicate planned integration time and planned completeness based on $\epsilon$ maximizing summed completeness. Blue dots indicate observation integration time and observation completeness of the simulated observation. White dots represent completeness at $\mathrm{\Delta }{\mathrm{mag}}_{\lim }$ and are plotted for the 10 highest, median, and lowest completeness planned targets. Blue diamonds show the completeness and integration time of the maximum $c_i/t_i$ point for the 10 highest, median, and lowest completeness targets. This plot references data generated using C0vsT0andCvsTDATA_WFIRSTcycle6core_CKL2_PPKL2_2019_10_07_11_58_.txt with specific target stars, integration times, and completeness included in Table 7.

From the combined time varying limit in Eq. (12) applied to Eq. (11), we get the black sigmoid-like lines in Fig. 12, specifically plotting $c_i(t_i)$ for the top 10, median, and lowest completeness optimized targets in the target list. The median and lowest $c_i(t_i)$ lines are characteristic for the majority of similar low completeness targets and the addition of targets will typically be below the lowest completeness target in this list. The $c_i(t_i)$ and completeness side histogram shows a clustering of targets at lower completeness values, which can be attributed to the increased number of targets at larger $d_i$. The upper limit of completeness lines is consistent with the theoretical maximum completeness values.

Demonstrating the importance of including overhead and settling times in observations are the max $c_i/t_i$ diamonds in Fig. 12. These are universally located at some small integration time ($t_i<{10}^{-3}\text{day}$, $\forall i\in \mathbf{I}$) and small completeness meaning any optimized target list maximizing $\sum \frac{c_i(t_i)}{t_i}$ without overhead constraints will have strictly nonzero $t_i>0$, $\forall i\in \mathbf{I}$. This means optimizing summed completeness without $T_{\mathrm{OH}}$ and $T_{\text{settling}}$ results in an observation target list of length $N$ (651 in the case of WFIRST), which cannot be executed under realistic conditions. Similarly, observing at $\mathrm{\Delta }{\mathrm{mag}}_{\lim }$ would result in suboptimal individual target completeness and also leaves a substantial amount of unobserved phase space around each target.

Calculating completeness using the SAG13 planet population shows the summed ultimate completeness of all 651 targets is 13.538. The summed ultimate completeness of the observed targets is 3.641. The minimum completeness of these targets has increased due to the increased likelihood of larger $\mathrm{\Delta }\mathrm{mag}$ planets at larger $s$ in the SAG13 distribution in Fig. 3(b). As in the Kepler-like optimized target list, we see a single target is not being observed due to strict enforcement of observing time constraints. The major difference from Figs. 12 and 13 is the change in shape of the $c_i(t)$ lines. The SAG13 $c_i(t)$ lines have a lower ultimate completeness, but more targets have this limit. In addition, each of these targets is observed for shorter integration times, which mostly have a higher theoretical maximum completeness and general shift toward higher completeness at lower integration times. There is additionally a larger separation of lower completeness targets in Fig. 13 compared to Fig. 12. In general, completeness calculated using the SAG13 population is larger than completeness calculated using the Kepler-like population.

Fig. 13

Completeness as a function of integration times calculated using the SAG13 distribution, EXOCAT-1 star catalog, Nemati SNR model,³⁵ and Leinert Table Zodiacal Light.³⁷ Black lines show completeness versus integration time for the 10 highest planned completeness targets, the median completeness target, and the lowest completeness target. Red squares indicate planned integration time and planned completeness based on $\epsilon$ maximizing summed completeness. Blue dots indicate observation integration time and observation completeness of the simulated observations. White dots represent completeness at $\mathrm{\Delta }{\mathrm{mag}}_{\lim }$ and is plotted for the 10 highest, median, and lowest completeness planned targets. Blue diamonds show the completeness and integration time of the maximum $c_i/t_i$ point for the 10 highest, median, and lowest completeness targets. This plot references data generated using C0vsT0andCvsTDATA_WFIRSTcycle6core_CSAG13_PPSAG13_2019_10_07_14_29_.txt.

3.2.

Sky Distribution of Completeness, Integration Times, and Targets

We are able to take the optimized target list included in Table 7 and bin the heliocentric ecliptic coordinates $(l,b)$ of each target, into triangular regions of approximately equivalent size and approximately isotropic distribution on the sky. When we sum integration time for all targets in each bin and normalize by bin area, we get the skymap distribution shown in Fig. 14. Since Fig. 12 shows all integration times are between 0.1 and 2 days, we can conclude the $(l,b)=(-140\mathrm{deg},0\mathrm{deg})$ and $(l,b)=(20\mathrm{deg},50\mathrm{deg})$ bins have the highest concentrations of observing time. By summing the bins over heliocentric ecliptic latitudes ($b$), we see a large disparity in $\sum t_i$ versus $l$, target count versus $l$, and $\sum c_i$ versus $l$. Since WFIRST is planned to be on an L2 halo orbit and has a Sun-orbital period of $\sim 365.25$ days, the Leinert local zodiacal light³⁷ is fixed in this rotating frame, and the time-distribution of stars is uneven, the optimally scheduled mission will have preferential observing times consistent with the distribution in Fig. 14. This result is important for optimally distributing limited CGI time under the constraints of a 5-year mission shared with multiple other instruments. Using this distribution, we can create preferentially distributed observing blocks for the CGI.

Fig. 14

The distributions of a Kepler-like optimized target list including (a) a skymap divided into approximately evenly sized triangular bins with isotropic sky distribution showing the time/area density of observations, (b) a histogram of total sky time versus heliocentric ecliptic longitude ($l$), (c) a histogram of target counts versus $l$, and (d) a histogram of summed completeness versus $l$.

3.3.

Detected Planet Properties

From our ensembles of survey simulations, we can look at how $a$ and $R$ of the detected planets are distributed. The top row of Fig. 1 is the average distribution of $R$ and $a$ for all generated planets in a universe. In each subplot of Fig. 1, the top left and bottom right show the number of simulations used to generate the resulting distribution. Each 2-D contour plot is normalized such that the integral over the area is 1 for the individual ensemble so the color scale can be shared. The white number in each gridded region shows the average number of planets generated or detected per simulation in that bin over the ensemble of simulations. The number in the top right of Figs. 1(a) and 1(b) is the sum total of all planets generated in the ensemble of universes. Subplots of Figs. 1(c)–1(f) show the average distribution of detected planets over an ensemble. The number in the top right of Figs. 1(c)–1(f) is the sum total of all planets detected in the ensemble. These summations have been tabulated as average yields in Table 5. Subplots of Figs. 1(c) and 1(d) use a target list optimized for the Kepler-like planet population to observe universes of Kepler-like and SAG13 simulated planets. Similarly, subplots of Fig. 1(e) and 1(f) use a target list optimized for the SAG13 planet population to observe universes of Kepler-like and SAG13 simulated planets.

Table 5

Summary of overfitting average unique detection yield and average characterizations from four Monte Carlo ensembles with optimized target list integration times calculated for different planet populations observing universes of different planets WFIRSTCompSpecPriors_WFIRSTcycle6core_3mo_405_19.

We can do an analysis on the kinds of planets WFIRST is expected to detect in a blind search survey. Since our universe is randomly generated, we show the distribution of generated $R$ versus $a$ planets for all stars for the Kepler-like and SAG13 planet populations in Fig. 1. The implemented planet generation rate for the universe of the Kepler-like planets is ${\eta }_{\mathrm{KL}}\approx 2.377$, consistent to two decimal places with the planet population model in Sec. 2.1. Here, $\eta$ is simply calculated by dividing the sum total of planets in the ensemble of universes by the number of simulations in the ensemble (1000 from Sec. 2.5.4) and number of target stars (651 from Sec. 2.2). The ${\eta }_{\mathrm{SAG}13}\approx 5.618$ for SAG13 is also consistent to within two decimal places with the model in Sec. 2.1. The generated planet populations are consistent with the limits presented in Sec. 2.1. In Sec. 2.1, we showed the smallest observable planet–star separation observable to be 0.292 AU and each of the four detected planets plots shows all planets detected have $a>0.292$ AU. A distinctive feature of the planet $a$ generation is the “knee” applied at 10 AU which can be seen by the sharp drop-off in both populations.

There is a marginal difference in the $y$ axis of the contour plots between Kepler-like and SAG13 populations, as SAG13 generates smaller planet radii than Kepler-like, so direct comparisons cannot be made between individual gridspace averages across universes. We also note the SAG13 universe generates an order of magnitude more large $R$ and large $a$ planets as indicated by the 24.29 and 221.19 grids from Kepler-like and SAG13, respectively [Figs. 1(a) and 1(b) gridspace (4,4)]. Comparisons between planet populations used to calculate completeness observing the same planet population universe indicates that optimizing with the Kepler-like population results in marginally smaller planet $R$ detections. Specifically for the Kepler-like universe’s largest average detection bin, we see that optimizing with the Kepler-like planet population results in 1.35 detected exoplanets on average and optimizing with the SAG13 planet population results in 1.11 detected exoplanets on average, which is significant, given our convergence results above. For the SAG13 universe’s largest averaged detection bin, we see optimizing with the Kepler-like planet population results in the 3.51 detections and optimizing with SAG13 yields 3.55 detections but scaling by the number of detected planets gives $3.54=3.51\times 16.266/16.101$ which is nearly identical. From these observations, we can conclude observing a universe of SAG13 planets with integration times optimized using completeness from a Kepler-like planet population result in more detections of small $R$ planets. Each gridspace in the bottom two rows of Fig. 1(d) is greater than or equal to each gridspace in the bottom two rows of Fig. 1(f). We can specifically point to the 0.81 and 0.69 grid spaces that most evidently confirms this observation.

We can draw several conclusions by inspecting the average population of planets detected. Our simulations show WFIRST will not detect planets with $a<a_{\text{mercury}}$. We also see WFIRST is not capable of detecting planets with $R<R_{\oplus }$. We also note that WFIRST is not sensitive to planets beyond $a_{\text{saturn}}$. The majority of planets detectable by WFIRST are cold-Jovians and Super Earths. Based on Fig. 1, we can conclude that optimizing integration times with the more pessimistic Kepler-like planet population universally biases detections toward smaller planets.

3.4.

Overfitting

We have chosen to use the Kepler-like and SAG13 planet populations to optimize target lists in this paper, both of which are created based on the known population of exoplanets. Since a motivation for the WFIRST CGI is to observe new exoplanets in an unexplored region of space, we must investigate how yield for a target list of integration times optimized for one planet population changes when observing a universe full of planets based on another planet population.

Optimizing a target list using completeness based on a planet population and observing that same planet population results in the highest yield in given in Table 5. However, observing a universe of Kepler-like planets with a target list optimized for SAG13 planets results in a 5.06% decrease in exoplanets detected, a greater decrease than observing SAG13 planets with a target list optimized using Kepler-like planets (a 1.01% decrease). A possible explanation for this difference is our inclusion of the rare characterization. The characterization part of Table 5 does not indicate this is the case and shows more planets are characterized on average when optimized with the wrong planet population range. A characterization observation is triggered whenever a planet is detected with sufficiently small $\mathrm{\Delta }\mathrm{mag}$ and separation such that the immediate reobservation could achieve an $\mathrm{SNR}>10$ with a newly calculated $t_i<30$ days and all other detection observation filters in Sec. 2.5. Since the observation would be carried out with a different instrument, the changed parameters are included in Sec. 7.

The conclusion we can draw from this exercise in overfitting is that optimizing with the more pessimistic Kepler-like planet population yields more robust detections and more characterizations if the planet population is actually SAG13. This warrants reinvestigation with varying margins placed on $t_{\mathrm{obs},i}$.

6.1.

Symbols

Throughout this paper, we use some common subscript conventions and notation discussed here. We use a generic variable $X$ to describe these conventions. The subscript $k$, such as in $X_k$ or $X_{x,k}$, is referring to a $X$ or $X_x$ of an individual planet. The subscript $i$, such as in $X_i$ or $X_{x,i}$, is referring to a $X$ or $X_x$ of a individual target star. The subscript min, such as in $X_{\min }$ or $X_{x,\min }$, describes that this parameter is a minimum of $X$ or $X_x$; this extends to the max subscript. We use bold variables such as $\mathbf{X}$ when we are referring to a set or collection of something. We use underlines such as $\underset{\_}{X}$ to identify position vectors and $\widehat{X}$ to identify unit vectors.

We make use of common variables differentiated by their subscript or boldness to differentiate between what they specifically refer to. All $ϵ$ parameters are related to efficiency of something. Lower case $a$ refers to a semimajor axis-related parameter. Lower case $p$ refers to an albedo-related parameter. Lower case $c$ refers to a completeness. Lower case $r$ refers to a direction or position vector. $t$ variables are all time related. $s$ refers to a planet–star separation related value. $a$ refers to a semimajor axis-related value. $R$ refers to a planet radius-related value.

7.1.

Zero-Magnitude Flux

The zero-magnitude flux, ${\mathcal{F}}_0$, used in the calculation of spectral flux density, $C_{{\mathcal{F}}_0}$, in Eqs. (12), (13), and (19), is given as