ExoVista: A Suite of Planetary System Models for Exoplanet Studies

The first module of exoVista, load_stars, loads the target list. Currently, load_stars reads the ∼8000 stars within 50 pc from the HabEx and LUVOIR master target list (Stark et al. 2019). The HabEx/LUVOIR master target list was generated by combining the Hipparcos (Perryman et al. 1997; van Leeuwen 2007) and Gaia TGAS/DR2 (Gaia Collaboration 2016; Marrese et al. 2019) target lists within 50 pc. All Hipparcos distances were updated with Gaia distances. The Washington Double Star (WDS) catalog was referenced to obtain simple binary parameters (separation and differential magnitude of the nearest companion); no other binary catalogs were referenced. Spectral types and photometry were obtained from the original Hipparcos release and updated with SIMBAD.

We remove all luminosity class III and IV stars, leaving only main-sequence (MS) and subgiant stars, along with unclassified luminosity classes that we assume are predominantly MS. We then estimate stellar angular diameters using the B − V analytic fits from Boyajian et al. (2014); we calculate physical stellar diameters from the angular diameter and distance. We use the MS stellar properties table from Pecaut & Mamajek (2013) to estimate stellar mass and effective temperature by interpolating B − V. We then use the estimated mass and stellar radius to estimate $\mathrm{log}g$. Finally, we estimate luminosity from the absolute V band magnitude and a color correction using Equation (9) of Torres (2010) and the color correction factors in footnote 10 of Pecaut & Mamajek (2013).

We emphasize that this target list is not intended to be accurate on a star-by-star basis. In addition to containing occasionally inaccurate stellar parameters, many M stars are missing due to the ∼8th magnitude limit of the Hipparcos catalog, and we expect a non-negligible fraction of these stars to have unknown binary companions. The current input catalog is intended to roughly represent the population of nearby MS and subgiant stars brighter than 8th magnitude. This target list could be substantially improved in accuracy and completeness, but a thorough analysis of every star in this list is beyond the scope of this project. While extending the target list to include more M stars should have a negligible impact on future direct imaging mission simulations that target potentially Earth-like planets, it would significantly expand the applicability of this tool to transit survey simulations.

We note that exoVista is not limited to the target list generated by load_stars. Users may supply their own target list as long as the format of the data structure is consistent and a minimum set of parameters are provided.

load_stars returns an n-element array of structures, s, where n is the number of stars. Each entry of s contains a stellar ID number, Hipparcos ID, Tycho ID, R.A., decl., UBVRIJHK photometric magnitudes, absolute V magnitude, distance, spectral type, bolometric luminosity, effective temperature, angular diameter, physical diameter, mass, $\mathrm{log}g$, and WDS separation and differential magnitude of the nearest companion star. Stellar spectra are determined later in the generate_scene module.

After the target list is loaded, we use the generate_planets module to distribute planets among stars. generate_planets modifies the structure array s to include a system midplane inclination, i, and position angle (PA). We draw inclinations randomly from a distribution uniform in $\cos i$ and position angles uniformly from [0, π).

We then randomly assign planets to stars. To do so, we use the occurrence rate map in mass-semimajor axis space for FGK stars from Dulz et al. (2020). Dulz et al. (2020) extrapolate the Kepler and RV occurrence rates to smaller masses and longer periods while using dynamical stability arguments to limit the number of planets at large distances. Importantly, we apply the occurrence rate map in $(M,{a}_{{L}_{\star }})$ space, where M is planet mass and ${a}_{{L}_{\star }}$ is the semimajor axis multiplied by the square root of the stellar bolometric luminosity. By doing so, the occurrence rate maps scale with the habitable zone such that η_Earth is independent of spectral type. We treat the occurrence rate map as a probability distribution. For all n stars we perform a Monte Carlo draw on every bin of the occurrence rate map to find the total number of planets to distribute to all stars in each bin. As a result, when running exoVista the number of planets generated may vary from run to run and be greater or less than the expectation value. We note that the occurrence rates of Dulz et al. (2020) have significant uncertainties at large periods. These uncertainties must be kept in mind when applying the astrometric or direct imaging techniques to the exoVista data set. To help estimate the impact of this uncertainty, we provide additional occurrence rate maps that correspond to the 1σ upper and lower limits.

After determining the total number of planets to generate for each $(M,{a}_{{L}_{\star }})$ bin, we randomly assign those planets to stars. We assume that planet mass and semimajor axis are randomly spaced logarithmically within each $(M,{a}_{{L}_{\star }})$ bin. We use planet mass to determine radius via the relationship of Chen & Kipping (2017), given by

$$\begin{eqnarray}R=\left\{\begin{array}{ll}1.008{\left(M/1.00\right)}^{0.279} & \mathrm{for}\,M\lt 2.04{M}_{\mathrm{Earth}}\\ 1.230{\left(M/2.04\right)}^{0.589} & \mathrm{for}\,2.04\leqslant M\leqslant 131.6{M}_{\mathrm{Earth}}\\ 14.31{\left(M/131.6\right)}^{-0.044} & \mathrm{for}\,M\gt 131.6{M}_{\mathrm{Earth}}\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where R is the radius of the planet. Currently, we assume a singular mass-to-radius conversion with no spread in the distribution. A realistic distribution of radii for each mass, which we leave for future work, would lead to a broader range of possible planet fluxes and transit depths. As long as that distribution is symmetric, we do not expect it to significantly impact future mission studies in the limit of large target sample sizes.

By default, we set all orbital eccentricities to zero, as was also assumed by Dulz et al. (2020) when estimating the impact of dynamical stability constraints on occurrence rates. However, the user can specify a nonzero range and the orbital eccentricities will propagate through the rest of the code self-consistently (e.g., when determining stability, the debris disk distribution, and integrating orbits). Eccentric orbits would lead to greater spacing between simulated planets as well as larger peak stellar radial velocity signals near periapse passage. We also assume planetary orbital inclinations are uniformly distributed from 0° to 5° by default, relative to the system midplane. All other orbital angles are randomly distributed between [0, 2π). Our orbital assumptions effectively exclude the presence of eccentric, highly inclined hot Jupiters. While such exoplanets are rare, studies using the radial velocity and transit methods, which are biased toward finding such planets, should not expect the dynamical aspects of this population to be represented within the exoVista data set.

After adding planets, we check all planet pairs to determine if they are dynamically stable. To do so, we take the same approach as Dulz et al. (2020) and calculate the mutual Hill radius, given by

$$\begin{eqnarray}&&{\rm{\Delta }}=2\left(\displaystyle \frac{{a}_{\mathrm{out}}-{a}_{\mathrm{in}}}{{a}_{\mathrm{out}}+{a}_{\mathrm{in}}}\right){\left(\displaystyle \frac{3{M}_{\star }}{{M}_{\mathrm{out}}+{M}_{\mathrm{in}}}\right)}^{1/3},\end{eqnarray} \tag{ 2 }$$

where a_in and M_in are the semimajor axis and mass of the inner planet, a_out and M_out are the semimajor axis and mass of the outer planet, and M_⋆ is the mass of the host star. We require Δ > 6 and delete all planets that do not meet this criterion. We then iteratively add more planets until the desired number per $(M,{a}_{{L}_{\star }})$ bin is generated. If any bin has not achieved its desired number of planets in 50 iterations, the creation of planets ends. We note that the occurrence rates of Dulz et al. (2020) adopted Δ > 9 to ensure Gyr of stability. Adopting the same criteria here greatly slows the planet distribution process. Because we are generating planetary systems that are plausible, and not necessarily strictly stable for many Gyr, we relaxed this criterion to substantially speed up the distribution of planets. We note that we do not currently attempt to reproduce the empirically measured planet multiplicity, nor do we include mean motion resonant configurations of exoplanets.

We then randomly assign each planet a wavelength-dependent geometric albedo based on simple rules applied to each $(M,{a}_{{L}_{\star }})$ bin. For the geometric albedos, we currently adopt the geometric albedos included in the Haystacks project (Roberge et al. 2017), which includes models of Venus, Earth, Mars, Jupiter, Saturn, Uranus, and Neptune. Briefly, the terrestrial planet albedos were derived from reflectance calculations by the Virtual Planet Laboratory using the Spectral Mapping Atmospheric Radiative Transfer model (Meadows & Crisp 1996). The disk-integrated reflectance of Earth was derived assuming an equatorial view averaging over the daily rotation, but can be applied to a broad range of viewing angles with ∼30% accuracy. The albedos of the giant planets were derived from the 300–900 nm albedo spectra of Karkoschka (1998). The 900–2500 nm SpeX observations of Rayner et al. (2009) were used to cover the near-infrared after dividing by the solar spectrum and scaling to match the albedo spectra of Karkoschka (1998). Further details of how these albedos were determined, as well as plots of their wavelength dependence, can be found in Roberge et al. (2017) and references therein. We emphasize that the exoVista model albedos are not intended to be used to fit to observations. Rather, they are intended to loosely represent the variety of exoplanet albedos that may occur in nature and serve as a known input for exoplanet spectral retrieval studies. We assign each albedo file an integer probability from 1 to 100 that determines the relative likelihood of drawing that planet in the given $(M,{a}_{{L}_{\star }})$ bin.

We treat the phase space of exoEarth candidates (EECs) separately to more easily control the distribution of planets within this region. We adopt the same definition of EECs as in the HabEx and LUVOIR final reports (The LUVOIR Team 2019; Gaudi et al. 2020). All EECs are assigned some form of an Earth-like atmosphere (Archean, Hazy Archean, Proterozoic High-O₂, Proterozoic Low-O₂, and Modern Earth) based on the approximate relative timescales of these phases during Earth's history (Arney et al. 2016).

The result of generate_planets is a modification to the structure array s. In addition to containing the midplane inclination and PA for each system, generate_planets adds an additional structure array p to each entry of s. p is a 30 element array containing each existing planet's parameters: mass, radius, semimajor axis, eccentricity, inclination, longitude of ascending node, argument of periastron, initial mean anomaly, and assigned albedo file. For a target list of ∼8000 stars, generate_planets typically takes a total of ∼10 s to run.

After the planets are distributed around each star, we use the module generate_disks to quasi-self-consistently generate parameters defining the debris disk of each star based on the underlying planetary system. Here we describe the debris disk generation process in detail, as this is a novel model developed specifically for exoVista.

Observations of debris disks have revealed several common aspects about their geometry that our models attempt to replicate. For a summary of observed debris disk geometries, see Hughes et al. (2018). First, observed disks, which are hundreds to thousands of times denser than the solar system zodiacal cloud, commonly appear as circumstellar rings. These rings have been shown to be relatively well fit with radial Gaussians at mm wavelengths, suggesting a well-constrained population of larger grains that track their planetesimal parent bodies. At short wavelengths we observe smaller grains, which are significantly impacted by radiative forces. The small grains create an extended halo of dust beyond the circumstellar ring to large distances, as they are blown onto eccentric orbits with large semimajor axes by radiation pressure. Radiative drag forces also cause these small dust grains to spiral inward from the circumstellar ring. For the massive, collisional disks we currently observe, only a small fraction of the dust spirals inward, such that the disk density interior to the ring is a much smaller fraction of the ring itself. However, theoretical studies suggest a much larger fraction of dust will spiral inward in less massive disks where collisions are rare, so much so that the disk can appear to be continuous all the way from the ring into the sublimation distance when viewed at visible wavelengths (e.g., Wyatt 1999; Kuchner & Stark 2010; Löhne et al. 2012; Kennedy & Piette 2015).

Observations have also revealed that debris disks commonly have a cold component at tens of astronomical units, and a warm component at a few astronomical units (Hughes et al. 2018). By default, we distribute two disk components (warm and cool) to each system, though a maximum of three is possible if desired by the user. To distribute the dust, we first determine plausible locations for the circumstellar ring of parent bodies based on the underlying planetary system. To do so, we determine all semimajor axes that are not within 3 Hill radii of a planet. For the inner warm component, we find the most stable region (defined as the largest width of stable semimajor axes in log space) between 0.5 and 5 au. For the outer cool component, we repeat the process from 5 to 50 au. If the user desires a third cold disk component, we repeat the process again from 50 to 500 au.

We model each disk component using a piece-wise combination of analytic models as shown in Figure 1. This includes a Gaussian parent body ring of half-width Δr centered at circumstellar distance r₀ (yellow line in Figure 1). The width of the Gaussian ring ranges from [0.05, 0.3]r₀, and is assigned based on the smaller of the range maximum (0.3) and the maximum width of the stable region at r₀. Only stable regions wider than the range minimum (0.05) are considered. Within the Gaussian ring we adopt a Dohnanyi size distribution (Dohnanyi 1969), given by

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{ds}}\propto {s}^{-3.5}.\end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

Outside of r₀ + Δr (blue line in Figure 1), we model the disk as an r^−1.5 power-law dust halo. Although the outer halo portion should have a non-Dohnanyi size distribution (Thebault et al. 2014), the expected size distribution is not well described analytically, nor well constrained observationally; we currently assume a Dohnanyi size distribution for the halo as well. The size distribution of our optically thin disk models mostly impacts the color of the disk, as the surface brightness was normalized to that expected for one "zodi." A steeper-than-Dohnanyi size distribution, with an enhanced population of small grains, would appear bluer, whereas a flatter size distribution would appear redder, which could impact exoplanet spectral retrievals.

Interior to r₀ − Δr (green line in Figure 1) we model the disk as inward-migrating, colliding dust. We use Equation (3.36) from Wyatt (1999), which accounts for the collisional destruction of dust grains via

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)=\displaystyle \frac{{\rm{\Sigma }}({r}_{0})}{1+4{\eta }_{0}(1-\sqrt{r/\left({r}_{0}-{\rm{\Delta }}r\right)})},\end{eqnarray} \tag{ 4 }$$

where Σ is the surface density, η₀ = t_PR(r = r₀ − Δr)/t_coll(r =r₀ − Δr), t_PR(r) is the Poynting–Robertson (PR) drag time evaluated at circumstellar distance r (Wyatt & Whipple 1950), and t_coll(r) is the collisional timescale evaluated at r. We approximate the collision time for each disk as

$$\begin{eqnarray}&&{t}_{\mathrm{coll}}=\displaystyle \frac{{t}_{\mathrm{orbit}}}{4\pi \left(z\times {10}^{-7}\right)},\end{eqnarray} \tag{ 5 }$$

where z is the number of zodis assigned to the disk and t_orbit is the orbital period. We calculate the PR drag time (Wyatt & Whipple 1950) individually for each grain size in the disk, such that the inward-migrating portion can have a significantly non-Dohnanyi size distribution for massive disks. As a result, the inner regions of the disk can exhibit a radial color gradient.

Systems that host massive planets can dynamically eject dust (e.g., Moro-Martín & Malhotra 2002). To include this effect, we multiply the dust distribution by a factor of ${(r/{r}_{\mathrm{truncate}})}^{3}$, where r_truncate is set to 1.1 times the apastron distance of the outermost planet interior to the belt that has a mass >100 Earth masses (red line in Figure 1). This truncation factor is motivated by the modeled optical depth reduction of inward-migrating Kuiper Belt dust, due to Jupiter and Saturn (Kuchner & Stark 2010).

The density of each analytic piece is scaled to ensure continuity at the transitions. All analytic pieces of a given disk component are randomly assigned the same scale height H/r₀ ranging from [0.03, 0.2]. The vertical variation of each disk component's density is modeled by

$$\begin{eqnarray}&&n\left(r,z\right)=\displaystyle \frac{n\left(r\right)}{r}{e}^{-0.5{\left(z/\left(r\left(H/{r}_{0}\right)\right)\right)}^{2}},\end{eqnarray} \tag{ 6 }$$

where n(r) is the radial density distribution, z is the distance normal to the disk component's midplane, and the 1/r factor results from the normalization of the vertical distribution.

Finally, each disk component is assigned a scattering phase function (SPF). Observations have revealed that the shape of debris disks' SPFs can be successfully described using a linear combination of Henyey–Greenstein (HG) functions (Stark et al. 2014; Sai et al. 2015; Matthews et al. 2017; Milli et al. 2017; Chen et al. 2020), expressed as

$$\begin{eqnarray}&&p\left(g,\theta \right)=\sum _{i}\displaystyle \frac{{w}_{i}}{4\pi }\displaystyle \frac{1-{g}_{i}^{2}}{{[\,1+{g}_{i}^{2}-2{g}_{i}\cos \theta \,]}^{3/2}},\end{eqnarray} \tag{ 7 }$$

where θ is the scattering phase angle, w_i and g_i are the weight and asymmetry parameters of the ith HG function, respectively, ∑_i w_i = 1, and g_i ranges from −1 for perfect backscattering to 1 for perfect forward scattering (Henyey & Greenstein 1941). Unfortunately, because we cannot access the full range of scattering angles of any observed debris disk, we cannot know for sure that fits to observations remain valid at all scattering angles (Hedman & Stark 2015). Critically, this lack of knowledge creates an uncertainty in the normalization of measured SPFs, making their application to models challenging.

In lieu of a measured, properly normalized SPF for debris disks, we chose to use the next best thing: the SPFs of Saturn's optically thin G and D ringlets, which serve as a good proxy for debris disks, have been measured over nearly the full range of scattering angles and can be properly normalized with reasonable accuracy (Hedman & Stark 2015). Hedman & Stark (2015) showed that a linear combination of three HG functions with w_i = [0.754, 0.151, 0.095] and g_i = [0.995, 0.585, 0.005] fit the measured SPFs well. However, simply adopting this single best fit would not properly sample the variety of SPFs that has been observed in debris disks (Hughes et al. 2018). Therefore, we opt to randomly vary the weights and asymmetry parameters of each of the three HG functions. Specifically, we randomly and uniformly select values of 0.85 ≤ g₁ ≤ 0.98, 0.5 ≤ g₂ ≤ 0.7, and − 0.2 ≤ g₃ ≤ 0.2 for the asymmetry parameters, and 0.6 ≤ w₁ ≤ 0.9, 0.1 ≤ w₂ ≤ (1.0 − w₁), and w₃ =1.0 − w₁ − w₂ for the weights. This variation produces the SPFs shown in Figure 2, which sample a variety of slopes and values near θ = 90°, where most debris disk dust is observed. This variety will alter the degree of forward scattering in individual disks, as well as the apparent albedo, while maintaining an average apparent albedo consistent with Hedman & Stark (2015). While we allow the SPF to vary between disk components, we do not currently vary it with dust grain size.

Zoom In Zoom Out Reset image size

The density of the inner warm disk is drawn randomly from the free-form distribution that best fits the LBTI HOSTS survey results, which has a median of 3 zodis of dust (Ertel et al. 2020). The outer cool disk is randomly assigned a density within a factor of 5 of the inner disk, in a basic effort to roughly correlate coeval disks. The uncertainties on the LBTI HOSTS best-fit distribution are significant and can impact studies of future direct imaging missions, as brighter disks make it more difficult to detect the faint planets. To allow the user to constrain this impact, we include alternative distributions that represent the 1σ upper and lower limits from the LBTI HOSTS results (Ertel et al. 2020). Importantly, we introduce an additional normalization factor to account for the SPF. We normalize the density of the disks such that a solar system twin model with our average SPF, inclined by 60° from face-on, produces a V band surface brightness of 22 mag arcsec⁻² at a separation of 1 au and scattering angle of 90°. This defines 1 "zodi" of dust using the same convention as the HabEx and LUVOIR Final Reports (The LUVOIR Team 2019; Gaudi et al. 2020).

We note that our disk models currently ignore debris disk asymmetries, including clumpy mean motion resonant structures (e.g., Jackson & Zook 1989) as well as warps and offsets due to secular perturbations from massive planets (e.g., Dermott et al. 1995). These asymmetries could have significant effects on future directing imaging methods by mimicking or obfuscating planets (e.g., Defrère et al. 2012). While future versions of exoVista may include some of these effects, we emphasize that the impact of these effects on observations must currently be examined separately.

generate_disks modifies the structure array s to include a new three-element array structure disk for each element of s. Each entry of disk includes all of the disk's parameters described above. For a target list of ∼8000 stars, generate_disks typically takes a total of a few seconds to run.

Finally, after all planet and disk parameters are generated for each star, we use the module generate_scene to model the astrophysical scenes. We first define a wavelength array for the star and planets (0.3–1.0 μm with spectral resolution R = 300 by default). We also define the wavelength array for the disk, which defaults to the same range as the star and planets but with reduced spectral resolution R_disk = 10. We can use a lower spectral resolution for the disk because the disk's spectral features are broad at visible wavelengths. The lower spectral resolution allows us to save substantial disk space for the output products. While the user can modify the wavelength array for any object, we emphasize that exoVista does not currently handle thermal emission.

Next, we define the range of dust grain sizes to include in the debris disk model. We set the maximum grain size to 100× the longest wavelength desired, which ensures that contributions from large grains do not contribute appreciably to the image. We set the minimum grain size to the blowout size, below which dust grains are removed from the system on short timescales by radiation pressure. While the blowout size is believed to vary with spectral type, we currently set the blowout size to 0.5 μm for all stars, as the blowout size of individual disks is poorly understood and depends strongly on the dust composition and porosity (Lebreton et al. 2012). We also set the dust grain size resolution, s/Δs. By default s/Δs = 5, which we have tested to produce adequate results using a wide range of grain compositions and size distributions.

We then loop over each star and generate its model outputs. For each star, we first retrieve an approximate stellar spectrum (Castelli & Kurucz 2004) based on the stellar effective temperature, luminosity, and log(g). We normalize the stellar spectral model to reproduce the star's observed V band flux. We do not normalize the stellar spectrum to the bolometric luminosity or radius, which are approximate quantities derived from analytic fits to B − V for main-sequence stars. We also do not attempt to fit the stellar spectrum to any additional observed photometry, and thus reddening due to extinction is not included.

Next, we distribute dust grains based on the geometry and parameters of each disk component and generate a contrast image of each disk component. We do this on a pixel-by-pixel basis over the full image of the disk, solving the flux in each pixel to a predetermined tolerance. By default, the pixel scale is 2 mas and the field of view is 05. The tolerance per pixel η is set to 0.05/z and limited to a range of [0.0005, 0.05]. This tolerance is roughly equivalent to 1/3 of a Δmag = 26.5 point source observed with a 4 m telescope, the assumed 3σ noise floor for HabEx. We limit η to a minimum of 0.0005 under the assumption that disks brighter than one hundred zodis will not be probed via direct imaging for planets with Δmag = 26.5, consistent with the yield simulations for the HabEx mission (Gaudi et al. 2020). Because pixels near the star dominate the run time, we set η = 0.1 for all pixels within 15 mas of the star, equivalent to 3.5λ/D for a 15 m telescope at a wavelength of 300 nm; we do not expect any future direct imaging mission to routinely access circumstellar distances this small and are justified in relaxing the precision of disk flux calculations at these locations.

To calculate the disk flux in each pixel in the image, we generate an array of coordinates subsampling $(x^{\prime} ,y^{\prime} ,z^{\prime} )$ space within that pixel, where $(x^{\prime} ,y^{\prime} )$ sample the east–west and north–south directions, respectively, $z^{\prime}$ samples along the line of sight in the image, and $z^{\prime} =0$ at the distance to the target star. Samples in the $(x^{\prime} ,y^{\prime} )$ dimensions are uniformly spaced, whereas samples in the $z^{\prime}$ dimension are spaced $\propto z{{\prime} }^{2}$ to achieve better resolution at the disk midplane and small circumstellar distances. We set the maximum value of $z^{\prime}$ along our line of sight to ${r}_{0,\max }{\eta }^{-2/7}$, where ${r}_{0,\max }$ is the maximum value of r₀ of all disk components. This ensures that $z^{\prime}$ extends to a distance at which the surface brightness falls to a fraction η of the brightness at ${r}_{0,\max }$, based on the r^−1.5 dust halo density distribution and an r⁻² illumination factor.

For each pixel, we start by subsampling the pixel using n_xy = 1 points in the $(x^{\prime} ,y^{\prime} )$ dimensions and n_z = 128 points in the $z^{\prime}$ dimension. We first resolve the disk in the $z^{\prime}$ dimension by doubling n_z until the total flux in the pixel varies by less than the set tolerance. We then resolve the disk in the $(x^{\prime} ,y^{\prime} )$ dimensions by doubling n_xy until the total flux varies by less than the set tolerance. The final number of points used in a given pixel is equal to ${n}_{z,f}\times {n}_{{xy},f}^{2}$, where n_z,f and n_xy,f are the final resolutions in the $z^{\prime}$ and $(x^{\prime} ,y^{\prime} )$ dimensions, respectively.

At each subsampled point, we calculate the disk density of all three disk components based on the parameters generated in the module generate_disks. We first transform the subsampled $(x^{\prime} ,y^{\prime} ,z^{\prime} )$ points to the disk coordinates (x, y, z), then evaluate the piece-wise analytic radial functions illustrated in Figure 1. If a massive planet is present, we include the truncation model shown in red. We then multiply by the vertical variation given by Equation (6).

Next, we multiply this density by the dust grain size distribution, scattering cross section, and scattering efficiency, and then integrate over all grain sizes. Dust grain scattering efficiencies are determined via Mie theory calculations. Performing Mie theory calculations on the fly would greatly slow the calculations, so we precalculated scattering efficiencies for specific grain sizes at the default size resolution s/Δs = 5. For each discrete value of s representing a bin of sizes of width Δs, we calculated the scattering efficiencies by averaging over 100 subsampled grain sizes assuming the Dohnanyi size distribution given by Equation (3) within the subsample. This subsampling is necessary to reduce "ringing" artifacts from Mie theory. We note that the assumption of a Dohnanyi distribution over our subsampled Δs does not limit us to a Dohnanyi distribution overall; tests show reasonable agreement between our subsampled method and a high-resolution calculation even for size distributions as steep as ${dN}/{ds}\propto {s}^{-5.5}$. Currently, scattering efficiencies have only been calculated for astronomical silicates (Li & Draine 2001); future versions may include additional dust compositions.

Finally, we multiply by the SPF of the dust, as given by the linear combination of three HG functions generated by generate_disks, the differential volume of each subsampled $(x^{\prime} ,y^{\prime} ,z^{\prime} )$ point, and sum over all points to determine the disk contrast per pixel (flux of disk per pixel divided by stellar flux). While the flux cube would have the star's spectral features imprinted on it, and thus need to be calculated at a higher spectral resolution, the contrast cube does not, allowing for lower spectral resolution, quicker calculations, and substantially smaller file sizes.

We note that in a handful of anomalous cases, a small number of pixels in high-zodi disks can converge on the final flux very slowly as we increase n_z and n_xy , oscillating about a final value. This behavior eventually dominates available memory and overall run time. To mitigate this issue, we limit n_z ≤ 32768 and n_xy ≤ 32, regardless of the desired tolerance. To track this issue, exoVista outputs a map of the approximate precision achieved in each pixel set equal to the absolute value of the fractional difference from the previous calculation. Users can reference this map to approximate the numerical noise in the debris disk image at all wavelengths.

Next, we evolve the orbits of the star and planets, recording the barycentric and heliocentric positions and velocities of all bodies as a function of time over the interval ${t}_{\max }$, resolved into timesteps of Δt. By default, ${t}_{\max }=5$ yr and Δt = 10 days. We use a Bulirsch-Stoer n-body integrator to calculate positions and velocities of all bodies to a tolerance of 10⁻¹² (Stark & Kuchner 2008). We include only gravitational forces and neglect the gravitational influence of the debris disk on the rest of the system. Because all disks are assumed to be azimuthally symmetric, we do not include the dust particles in the integration, as the disk's appearance is assumed to remain fixed.

We then calculate each planet's reflected light contrast as a function of wavelength for each time step/position. To do so, we load each planet's geometric albedo spectrum at its native resolution, typically much higher than the default R = 300, calculate the stellar spectrum at the same wavelengths, and determine the planet's reflected light flux at the higher native resolution. We then bin down to the desired output resolution, taking care to calculate the planet's flux at the exact transition wavelength between the lower resolution spectral bins. To convert to contrast, we divide by the stellar flux at the lower output resolution, such that users can retrieve the planet spectrum calculated at high fidelity simply by multiplying by the output stellar flux.

Finally, we produce an output scene file in .fits format. This file includes all of the stellar parameters generated by load_stars, the stellar position, velocity, and flux spectrum as a function of time, all of the planet parameters generated by generate_planets, the position, velocity, and contrast spectrum of all planets as a function of time, and the disk contrast data cube.

The generate_scene module is the most numerically taxing step of the exoVista code. Whereas previous modules execute in seconds, this step requires ∼1 day of execution on a 32 core machine to generate ∼8000 astrophysical scenes. The vast majority of the time is spent generating the debris disk image. Typical run times are ∼4 minutes per planetary system using a single core.

The outputs of exoVista are .fits files that adhere to a specific format. The format is described in detail in the public documentation. To simplify the use of the data sets we have included a sample module called load_scene in both IDL and Python. load_scene provides an example of how to convert the data products into an astrophysical data cube at a given time. load_scene opens the desired .fits file, interpolates all positions, velocities, and spectra to the desired epoch, and returns these quantities along with a disk flux data cube, allowing the user to easily generate an image by convolving with a user-supplied point-spread function (PSF).

Figure 3 illustrates 28 random V band scenes generated by exoVista via the load_scene module. Each image is individually normalized, measures 03 across, and neglects the stellar contribution. We modeled point sources using a 10 m telescope with an Airy pattern PSF. For illustrative purposes, to highlight the fidelity of the disk model, we did not convolve the disk with the PSF. A movie showing the systems in motion over a 5 yr time span is available for download. ² Some disks show central clearings associated with massive planets due to the truncation feature of our disk model. Others show ring-like features associated with the parent body belt. Some systems exhibit bright inner planets, as evidenced by bright Airy rings extending to large distances. In highly inclined systems, these bright inner planets can create a stroboscopic effect as they quickly orbit from crescent to gibbous phase. The scenes shown in Figure 3 can be used to study a wide range of aspects of a future direct imaging mission, from how well we can model and subtract exozodi, to planet–planet confusion and orbit determination, to survey optimization.

We can also use the outputs of exoVista to investigate other exoplanet detection methods. Figures 4 and 5 plot the barycentric radial velocity and astrometric position of a single star from our models, which clearly exhibit the signatures of several planets within the system. We note that exoVista does not currently include stellar noise, the gravitational influence of companion stars, or photometric asymmetries within the planetary system. As a result, these simulated signals are relatively crude and any detailed study of an RV or astrometric survey applied to these models would need to include such effects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 shows the position information plotted for multiple planets in a single near-edge-on system, each shown with a separate color. We plot the stellar radius at the center (black) to show that the innermost planet in this system (red) is predicted to transit. By comparing the planet radius with the stellar radius, users can generate transit light curves from the exoVista models. While the majority of transiting planets in the data set should be reliably detected by interpolating the default time step of 10 days, grazing transits may require a finer time step.

Zoom In Zoom Out Reset image size

ExoVista can also be used to perform rudimentary TTV analyses. Figure 7 shows variations in the transit time of a modeled planet that transits, plotted as the difference of the ingress time from the average. Here we have simply interpolated the planet's position to estimate the moment of ingress. Detailed TTV analyses may require a finer time step in the data set.

Zoom In Zoom Out Reset image size

The exoVista source code is publicly available for download. ³ The exoVista code is written in a combination of IDL and C. The provided documentation describes the installation process for several use case scenarios.

We expect that most users of exoVista may be more interested in the output fits-formatted models than the source code. Therefore, we have pregenerated planetary systems for 2258 stars within 30 pc from the LUVOIR/HabEx master target list, as a preliminary exoVista "universe." Because each realization of a planetary system is random, and a single instance does not capture the stochastic nature of these calculations, we repeated this process a total of 12 times. All 12 sets of 2258 planetary systems, which we refer to as the DEC21 data set, are available for download via the Goddard EMAC web server. Users can access the files via a browser-based file viewer. ⁴ At this same location we provide a README.txt file that gives instructions on how to search for planetary systems that meet some parametric criteria and then retrieve those matching data sets. The data files can be read using the IDL or Python versions of the load_scene module in the exoVista source code distribution.