Starshade Imaging Simulation Toolkit for Exoplanet Reconnaissance

Starshade Imaging Simulation Toolkit for Exoplanet Reconnaissance (SISTER) is a versatile tool designed to provide accurate models of the images of exoplanet systems when observed with a starshade positioned to block the light from the host star. SISTER allows one to control a set of observational parameters including: (1) the starshade design, position, orientation, and glint properties; (2) the telescope and optical system pupil, aberrations, bandpass, and throughput including a detector model; (3) the exoplanetary system, including stellar distance and spectral type, parallax and proper motion, planet size, reflection properties, orbital parameters, and exozodiacal dust; and (4) background objects. Additionally, there is a substantial library of built-in plotting software added, but the simulations may be stored on disk and plotted with any other software. We describe SISTER’s algorithms, its operational modules, and how it can be used to generate starshade optical simulations with a high degree of fidelity. We include some imaging examples. © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI. [DOI: 10.1117/1.JATIS.7.2.021217]


Introduction
Starshades have emerged as a promising concept for direct imaging and spectral characterization of exoplanets. With their high throughput, small inner working angle (IWA), unlimited outer working angle (OWA), and broad spectral bandwidth, a starshade transforms a moderate-to large-aperture telescope into a powerful exoplanet imaging instrument. The starshade concept is illustrated in Fig. 1.
In parallel with technology development and mission studies, design reference missions [17][18][19] have established the starshade's potential rich scientific returns. To date, these studies have relied on simplistic representations of starshade performance, e.g., spatially stationary imaging point spread functions (PSFs), representative stellar leakage profiles, and average solar glint characteristics, to suit the specific study needs. Oakley and Cash 20 studied the ability of the New Worlds Observer starshade to determine the nature of the surface and clouds and the diurnal rotation rate of Earth-like planets. They implemented highly detailed planetary models and used the *Address all correspondence to Sergi R. Hildebrandt, srh@caltech.edu ZODIPIC 21-23 model to generate the exozodiacal light distribution while accounting for starshade throughput. However, since the precise nature of the PSF was not required for their study, they did not perform a detailed convolution with a spatially variant PSF. Likewise, the Exo-S mission concept study 1,11 made use of throughput determined by detailed propagations but did not implement a PSF convolution with localized PSFs.
Hu et al. 24 developed a detailed convolution code to produce realistic images that could be used to test image processing algorithms and the impact of starshade shape errors, e.g., petal displacements. Starting with astronomical images from the Haystacks model 25 with a 0.1-AU pixel scale, they used Fresnel diffraction to calculate the telescope pupil field for all pixels that were significantly affected by the presence of the ideal or perturbed starshade. They used direct Fourier convolution for pixels beyond the starshade "region of influence" such that the imaging process was stationary.
In this paper, we describe the Starshade Imaging Simulation Toolkit for Exoplanet Reconnaissance (SISTER), a versatile suite of MATLAB-based software for accurately simulating images of exoplanet systems observed with a starshade (Fig. 2). SISTER allows one to  An input astrophysical scene gets convolved with the optical response of the telescope and the starshade at different wavelengths to produce a broadband image or spectral-image data cube. An example of an astrophysical scene, together with some telescope perturbation, is provided in Sec. 4. The resulting intensity at the image plane for a broadband simulation can also be found in Sec. 4. The diffraction code is reviewed in Sec. 3.1. A few examples with starshade perturbations are shown in Sec. 3.8. An example of a static starshade is provided in Secs. 3.4 and 3.5, whereas Secs. 3.6.2-3.7 deal with a spinning starshade. The detector simulator is described in Sec. 3.10. An example of a spectral simulation can be found elsewhere. 26 For an extensive set of examples and simulations, see Ref. 27. control a set of observational parameters including: (1) the starshade design, position, orientation, and glint properties; (2) the telescope and optical system pupil, aberrations, bandpass, and throughput including a detector model; (3) the exoplanetary system, including stellar distance and spectral type, parallax and proper motion (PPM), planet size, reflection properties, orbital parameters, and exozodiacal dust; and (4) background objects. SISTER has recently been used to generate the starshade simulations of the Starshade Exoplanet Data Challenge. (The simulations and corresponding documentation are publicly available in Ref. 28.) 29 This paper describes SISTER's algorithms, its operational modules, and how it can be used to generate starshade optical simulations with a high degree of fidelity. The MATLAB code for all examples is available along with the software suite and a detailed handbook in Ref. 27.
This paper is organized as follows. First, in Sec. 2, we describe the astronomical components that are modeled, from the host stars and their spectra to planetary and orbital models, zodiacal light, and background sources. Next, in Sec. 3, we describe the instrumental components, including the optical model which, in a nutshell, is the convolution of the system PSF with the astronomical scene, with added detector and instrument noise including solar glint from the edge of the starshade. We also discuss in this section the implementation approach that enables SISTER to work efficiently on a modern laptop. Section 4 includes an example demonstrating imaging of the components that make up the astronomical scene.

Overview
The input astrophysical image, also called the scene, is generated with a combination of SISTER's internal building blocks and externally generated images, e.g., extracted from Haystacks 25 data cubes, or any other image with the correct irradiance units and spatial scale, as discussed below. SISTER's built-in scene components include the host star and its planets, the exo-zodiacal emission, an exo-Kuiper belt, our local zodiacal light, and a sky model for ground telescopes.
The scene is defined on a regular equatorial system grid with the pixel size expressed in mas. The grid spacing assumes true arc-length since the scene's field of view (FOV) is small enough to be considered a flat surface (of the order of 1 arc sec 2 ). The pixel pitch is generally a few mas. The energy units of the scene are spectral irradiance at the telescope aperture, also called flux density F λ or F ν , in units of W∕m 2 ∕μm and Jy, respectively (see also Appendix B). When using F ν in SISTER, the data are still ordered by wavelength.

Host Star
The stellar spectrum and other star characteristics may either be user-defined in SISTER or read from a star in ExoCat. 30 ExoCat includes 2347 stars taken from the Hipparcos Catalogue with measured parallaxes corresponding to distances of ≤30 pc. This sample is nearly complete down to V ¼ 8, corresponding to stars brighter than 0.5L ⊙ . ExoCat also includes several stellar properties that are used by SISTER: parallax, proper motion, spectral type, luminosity, V-band photometry, and estimated stellar masses and radii. Figure 3 shows the distribution of the stars as a function of distance and V magnitude. Among them, we highlight β CVn, 31 a solar type star at a distance of 8.44 pc, that we use to run some examples of SISTER later in this paper. β CVn has also been considered in some starshade mission studies such as the starshade rendezvous mission with the Roman Space Telescope 32 and HabEx. 12 A recent analysis based on ancillary radial velocity data and dynamical stability shows that β CVn could harbor several terrestrial planets in its Habitable Zone. 33 More precise catalogs exist, e.g., Gaia, 34 that are not provided with SISTER although one may pick up any star by setting its properties, such as stellar type, distance to Earth, apparent V magnitude, stellar mass, radius, coordinates, and proper motion. SISTER's handbook 35 provides full details on how to set up the host star properties.
Stars and planets are treated as point sources. SISTER presently has no built-in accommodation for the finite-angular diameter of the stars. For most stars, the effect leads to a small reduction in the diameter of the starshade shadow, which could impact formation flying requirements. For example, a star with a 1-mas diameter would reduce the shadow diameter for HabEx (telescope-starshade separation = 76.6 Mm) from its design value of 6 to 5.6 m. A star with a large angular diameter, e.g., alpha Cen with d star ¼ 8 mas, would reduce the shadow size to below the 4-m telescope diameter, severely compromising the starshade's ability to suppress the starlight. The finite stellar diameter can be approximated by modeling the star as a set of closely spaced point sources with identical spectra while preserving the star's total flux. This function will be included in a future release of SISTER.

Stellar Spectra
If an external stellar spectrum is not provided, SISTER finds the closest stellar spectrum to one of these main types: A0V, A5V, F0V, F5V, G0V, G5V, K0V, K5V, M0V, and M5V. SISTER uses the spectral irradiance normalized to V ¼ 0 with respect to Vega's flux. The stellar spectra are stored from 200 to 2000 nm in the original steps of 0.5 nm. Figure 4 shows stellar spectra for a range of stellar types, including the solar spectrum. 36 In the future, we will incorporate the SSFL library of 131 stellar spectra. 37 Combining the information from ExoCat and the stellar spectra, Fig. 5 shows the number of photons per unit area and time over two passbands: 425 to 552 nm "blue" and 615 to 800 nm  "green" that we use later when producing some simulation examples. Figure 5 can be used to estimate the photon flux arriving at the telescope by multiplying any value on it by the collecting area of the telescope A C and the planet's average flux ratio over the observing passband Φ. In addition to these factors, one has to include other instrumental factors, such as optical transmission, quantum efficiency of the detector, and any other photon losses in the system, some of which are wavelength dependent. SISTER provides these factors for Roman and HabEx.

Planets
SISTER simulates the reflection of starlight by exoplanets; it does not include infrared selfradiating exoplanets. The relative position of the exoplanet, the host star, and the observer, as well as the particular physical properties of the planet, determines how much radiation from the exoplanet arrives at the starshade-telescope location.
The electromagnetic radiation arriving from an exoplanet due to reflected light is summarized as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 3 2 3 where F p is the flux of the planet, F s is the flux of the host star, and Φ is the flux ratio between the star and the planet. All of them depend on the wavelength of the radiation. Moreover, Φ also depends on the actual distance of the planet to the star r p−s that varies with time, the planet's radius R p , phase angle α, phase function f p , and geometric albedo g p , i.e., the planet's albedo observed in superior conjunction (f p ¼ 1, α ¼ 0) by a distant observer. The phase function and geometric albedo also depend on the wavelength. SISTER has compiled a few geometric albedos for the solar system planets from publicly available data 25 and some examples of archean atmospheres for Venus, Earth, and Mars. 38 Figure 6 shows a summary of the geometric albedos available in SISTER, 25,38 and Fig. 7 shows the phase function choices available in SISTER, where ω is the degree of forward scattering asymmetry, 39 although the user may use any external phase function. 40 SISTER can generate the Keplerian orbit (a two-body orbital solution defined by six Keplerian elements) and the projected orbit on the image plane, including the associated phase angle α following well-known expressions (see Ref. 41 and SISTER's handbook 35 for more details). The gravitational interactions among the planets are neglected. If the user sets some, or all, of the Keplerian elements and the phase angle, SISTER checks that the equations for orbital motion and phase angle can be made consistent and lets the user know the adopted values. Otherwise, it lets the user know where the conflict occurred. For simplicity, SISTER includes the physical information for all planets in the solar system. When setting the planet's properties in SISTER, one may modify any of the parameters while keeping the remaining unchanged. For instance, one may move Jupiter closer to the host star or define a super-Earth with a larger radius than Earth. One can use Eq. (1) with the stellar photon flux from Fig. 5 and the spectra of Fig. 6 to estimate exoplanet photon rates per unit area at the telescope's aperture. One can also consider planets with a relative incident stellar flux that is the same for stars of different luminosity. This may be accomplished by rescaling their semimajor axes by the square root of the luminosity of the star. Figure 8 derives the photon flux for an exo-Earth with a semimajor axis of 1 AU rescaled by the square root of the luminosity, integrated over the 425-to 552-nm and 615-to 800-nm passbands. Finally, for reference, Table 1 lists an averaged planet flux ratio Φ assuming the average orbital properties of the solar system planets observed in superior conjunction (f p ¼ 1, α ¼ 0) by a distant observer and integrated over the respective spectra. We would like to remark that these photon rates per unit area are usually converted into photon rates per detector pixel, although that depends on the specific characteristics of the instrument under consideration. One needs to take into account the collecting area of the telescope, as well as other instrumental factors, such as optical transmission, quantum efficiency of the detector, and any other photon losses in the system, some of which are wavelength dependent. Finally, one also needs to incorporate the angular size of the detector's pixels and some region of interest of the full PSF response [e.g., 1 full width at half maximum (FWHM)]. The values for each of these factors for starshade rendezvous and HabEx mission concepts can be found in their respective studies. 12,32 Let us remind the reader that, for the cases of Roman and HabEx, SISTER incorporates all of these factors, including any wavelength dependence.

Local Zodiacal Light and Scattered Starlight
Local zodiacal light is introduced following the STScI model. 42 The spectrum of the zodiacal light follows that of the Sun between 0.2 and 10 μm because the scattering strength depends weakly on the wavelength. Given the spatial distribution of the surface brightness of the local zodiacal light and the very small angular areas being imaged (∼1 arc sec 2 ), the surface brightness is effectively a spatially uniform background that affects the Poisson noise level of the simulations. Figure 9 shows the local zodiacal light brightness in units of V magnitude per arc sec 2 for different values of the heliocentric Ecliptic coordinates derived by SISTER based on the HST data. 42 SISTER's implementation exactly matches the data at the same coordinate grid points and linearly interpolates the HST data for any other coordinate choice (Appendix B.3 provides more technical details). In Fig. 10, we estimate the photon rate arriving at the telescope aperture from local zodiacal light for Roman taking into account the camera's projected pixel scale of 21.85 mas. 43 The median value of the local zodiacal brightness in the regions shown in Fig. 9 is 23.02 V mag∕arc sec 2 , which is typical of the value used by several exposure time calculators, 12,32 and it corresponds to 1.44 × 10 −2 photons∕s in the 425-to 552-nm Roman passband and 2.23 × 10 −2 photons∕s in the 615-to 800-nm band. In practice, <10% of incident photons reach the detector due to optical losses and detector losses. The contribution of reflected starlight from the Milky Way (maximum value of ∼20.5 V mag∕arc sec 2 ) on the starshade is weaker than that of local zodiacal due to the low reflectance of the starshade's surface (<5%). However, it is still possible to have galactic (individual) background stars.

Exozodiacal Dust and Debris Disk
The exozodiacal emission is one of the least well-known factors in an exoplanetary system. Its spatial distribution and composition may be quite complex. Moreover, the population of its brightness distribution is not well constrained. In the case of significant asymmetries, such as a blob structure or a hole, these should orbit around the host star as well. SISTER does not generate a data cube of simulated data for an arbitrary exodust emission, although the capability to interface or incorporate some other public tools may be added in the future.
However, one may generate externally a data cube of exodust scattered light and merge it with internally generated SISTER components (e.g., the simulations generated for the Roman CGI data challenge with SISTER; see Ref. 44). Since the generation of exozodiacal dust templates is not trivial, in general, we suggest the use of ZODIPIC, 21,22 which can produce  an exozodiacal dust template that resembles that of the solar system. The result is a smooth dust distribution, but it may be sufficient in some circumstances. SISTER scales the output of ZODIPIC to the stellar spectrum of the host star assuming that the scattering strength of the exozodiacal emission follows the host star spectrum from 0.2 to 10 μm. 42 SISTER uses a version of ZODIPIC that was used in the Haystacks model 25 that applies a Henyey-Greenstein phase scattering function. The software can be found in Ref. 45. In Sec. 4, we show an example.

Parallax and Proper Motion
PPM is computed in SISTER simulations. Both effects may help to disentangle exoplanets from background objects when observations at different epochs are considered. The model assumes that the Earth's orbit is circular and coplanar with the ecliptic and uses the J2000 equinox.
Starshade observations of a given target will likely be spaced several months apart due to scheduling constraints on the line of sight relative to the Sun and the number of targets in the observational program. After a few months, PPM moves many targets by a distance larger than the FWHM of the PSF, as shown in Fig. 11. Figures 11(a) and (b) identify the FWHM for VRI photometric passbands 46 (V ¼ 551 nm, R ¼ 658 nm, and I ¼ 806 nm). The motion increases the likelihood that background objects can be disambiguated from planets with a pair of observations.

Extra Galactic Background
The Haystacks Project 25 provides, in addition to the present and archean solar system models, data for a deep extragalactic field. The original Haystacks field has an angular extension of 36 arc sec × 36 arc sec, much larger than the actual FOV used for most simulations with SISTER, and it has a native pixel scale of 10 mas. SISTER takes care automatically to crop and resize the Haystacks extragalactic background data to match the pixel scale and FOV selected on the astrophysical scene. SISTER allows one to choose to center the extragalactic background image at any location on the Haystacks data. It also provides a way to add background stars to its simulated data. In the current version of SISTER, there is no specific option to add background stars.

Sky Model for Ground Telescope
For ground telescopes, SISTER includes the total sky radiance and transmittance based upon the Cerro Paranal Sky Model Calculator 47,48 (see the SISTER's handbook 35 for an example and more details). Briefly, SISTER considers three scenarios for the sky radiance and transmittance depending on the moon phase. (The files were provided by Stefan Kimensberger at Universität Innsbruck, Austria.) Each of them includes scattered moonlight, local zodiacal light, scattered starlight, emission lines of the upper atmosphere, airglow (residual continuum), and molecular emission of lower atmosphere (the latter essentially zero below 1.7 μm). Figure 12 shows the atmospheric transmittance and total radiance of the ELT sky model with half moon phase.
However, SISTER does not provide a convolution of the astronomical scene that takes into account the effects of the turbulence of the atmosphere into the optical response of the telescope. The latter depends on several factors that vary on a case by case basis. If the telescope's PSF is known, it is possible to postprocess SISTER's simulation to the actual telescope's turbulent PSF.

Instrumental Components
Imaging simulations in SISTER are performed in two steps. First, a set of basis PSFs is generated. These determine the imaging properties over a grid of source positions and wavelengths. The basis computation can be quite computer intensive although it only needs to be done once for a given starshade-telescope system. We utilized image plane symmetries where possible and minimized the number of computations by judicious choice of starshade and source parameters to minimize the computation time while preserving accuracy. Once the basis is set up, the second step is to define the scene and call SISTER to generate the image by convolving the scene with the set of basis PSFs. In particular, to obtain the results of this work, we derived 206 independent PSF bases. In a standard dual-core laptop, the computation of a single PSF basis would take on average 16 h and 800 Mb of memory. In practice, we used a cluster with 32 cores. An imaging simulation, like the one shown in Section 4 that aims at 1% errors in the optical response, takes about 90 s to complete. In general, SISTER prints out timely information on the progress of both the imaging simulation and PSF basis construction, which may be paused and resumed at any time (see SISTER's handbook 35 for full details).
The main goal of this section is to derive the PSF response for a set of instrumental and optical parameters and determine the parameter values that lead to a subpercent PSF error for a source position where the starshade transmittance τ ≥ 0.25, which includes an area on the image plane well inside the nominal IWA (see Fig. 13). This range of locations should cover any practical application of the starshade simulations. Higher precision comes at the price of longer run times. We recognize that the SNR for exoplanet observations will almost always be low (i.e., SNR ≲ 100); thus once the noise-free PSF is shown to be accurate at levels below 1%, we can expect that the images will in general be limited by photometric and detector noise.

Diffraction Algorithm
SISTER uses the boundary diffraction wave (BDW) algorithm 49 to compute the electric field of the scene at the telescope pupil. BDW represents the electric field by a combination of geometrical optics and a wave that diffracts from the boundary of the aperture. The algorithm is based on the Kirchhoff formulation 50 and efficiently computes a 1D integral around the edge of the starshade. It is similar to the 1D algorithm of Dubra and Ferrari 51 which was first implemented for starshades by Cash. 52 The efficiency is an important characteristic as the starshade is a broadband imaging system with a non-stationary point-spread function.
A starshade's performance is measured in terms of the normalized intensity, defined as the ratio of intensity averaged over an image plane resolution element (typically the PSF FWHM) to the peak intensity of the target star when the starshade is not present. Because we are mainly concerned with planets located at the angles ≥IWA, where the throughput is ∼unity, the normalized intensity is identical to instrument contrast, which is more commonly used to express starshade performance. Harness et al. 53 provided a formal definition of the contrast. Since exo-Earths have a flux ratio of ∼10 −10 relative to their host star, it is important for the starshade to achieve ∼10 −10 instrument contrast (often simply called "contrast") for a high signal-to-noise ratio observation. This also sets the requirement on the accuracy of the algorithm that predicts starshade performance.
The algorithm has been validated to contrast levels of at least 10 −10 both in the laboratory and with respect to independently written codes. Harness et al. 2 performed experiments at flight-like Fresnel numbers using 1/1000 scale starshade masks. Harness et al. 54 designed the masks using his independently written version of a 1D algorithm and verified the performance against a 2D gridded Fresnel propagation to contrast levels better than 10 −11 . The SISTER BDW algorithm evaluated the mask and predicted results consistent with this level as well. The laboratory results have a measured contrast floor of 10 −10 at the IWA (corresponding to the tips of the inner starshade) and a floor of 2 × 10 −11 at the OWA. 53 The laboratory performance was limited by polarization-dependent propagation through the few-micron thick gaps in the masks, which was significant at the laboratory scale but will be well below 10 −12 contrast at the full flight scale. These results validate the accuracy of the diffraction algorithm to at least a contrast of 10 −10 .

Starshade Mission Concepts
This work emphasizes the starshade rendezvous 32 and HabEx 12 missions. These two missions give us enough variety to demonstrate how SISTER can provide optical simulations with errors below 1%. An important difference between these missions is the nature of their pupils. Although Habex has a circular, unobstructed pupil, SRM has a large secondary mirror and six non-radial secondary struts. Another important difference is the bandpass; SRM operates in two distinct bandpasses, 435 to 552 nm and 615 to 800 nm, with the starshade position shifted to accommodate each bandpass. Habex, on the other hand, is designed to work from 200 to 1800 nm, with a broadband starshade that can accommodate a bandpass of 400 to 1000 nm in a single position. 12 To make the comparison between both missions more useful and to limit the number of cases considered, we compare the same passbands for both cases: 425 to 552 nm and 615 to 800 nm. In the case of starshade rendezvous, the IWA for the 425-to 552-nm band is 72 mas, whereas for the 615 to 800 band, it is 104 mas. In the case of HabEx, the IWA is 70 mas for both bands. Both starshades have 24 petals and are spinning about the axis of their central disk. An analogous analysis to the one presented here would also find the adequate set of parameters for any other spinning or non-spinning starshade mission.

Methodology
In the following sections, we change parameter values, e.g., the number of points used to represent the outline of the starshade, until we observe convergence of the PSF peak response, overall photometry, and astrometry to a precision better than 1%. The most extreme value of the parameter under consideration will be called the "reference" value. The reference value is usually set based on some experience with the analysis or by computational limits. In any case, we show consistently that the adequate values are less extreme than the reference ones.
The convergence depends on the relative position of the source with respect to the center of the starshade. For instance, it is clear that, for sources located far from the starshade, the diffraction effects will be much less than when close to the IWA. In our study, we consider all of these possible locations, including the central one where the star is located. However, for the sake of simplicity in the implementation of SISTER's algorithms, we choose a single adequate value for each parameter that is valid for all relative positions of the source and the starshade.
In the context of direct imaging of exoplanetary systems, we chose three tests that will be applied to each parameter analyzed in this study. These tests evaluate the maximum difference between the PSF and its reference value, the integrated photometry of the PSF compared with its reference value, and the astrometric position of the PSF compared with its reference value. We do not add any instrumental noise to the PSF response, and our results will in most applications be conservative because the 1% criterion applies to the worst case position in the field, which generally occurs within the radius of the IWA.
The first test estimates the maximum difference of the PSF response between a value of one parameter and its reference value. We define the maximum percentage difference with respect to the reference PSF's peak, or simply PSF peak test, as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 3 2 7 where PSF PAR is the PSF response corresponding to some value of the parameter under consideration and PSF REF is the one associated with the most extreme value of the parameter considered. In this expression, MaxDiffðSÞ ¼ s ∈ S∕absðsÞ ≥ absðsÞ ∀ s ∈ S, for an array S, and max is the maximum value of the array. MaxDiff allows one to preserve the sign information of the differences. We note that this test is the maximum difference between the PSFs rather than the difference between the maxima. It is thus a simple means of testing changes in the PSF morphology even when the photometry and astrometry are unaffected. For the photometric test, since our PSF response is known for any of the parameter values studied, we choose a PSF fitting method. Let us assume that we have an input point-like source S at some location of the image plane ðx 0 ; y 0 Þ with some peak amplitude normalized to 1; then the corresponding optical response with the reference PSF is simply given by On the other hand, if one estimates the photometry of the same input source with a PSF corresponding to some different values of the instrumental parameter, one would write S ¼ bPSF PAR ðx 0 ; y 0 Þ, where b in general is different from 1. Since there's no background in our simulations and the images are noiseless, the best estimate for b can be obtained with minimum least squares. Thus we define the percentage difference with respect to the reference PSF's photometry, (b − 1), or simply PSF photometry test, as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 7 3 5 where · means the dot product of two arrays with the same number of elements, A · B ¼ Σ i A i B i , i ¼ 1; : : : ; N, and N is the number of elements of A. Finally, for the astrometric test, we compare the PSF centroids corresponding to different values of each parameter and the corresponding reference values. To give a better sense of the relative importance of the centroid differences, the absolute difference between the centroids is then compared with the PSF FWHM of the reference case far from the starshade center, i.e., the PSF response associated with the telescope itself that we will call hereafter "stationary" because it does not vary with distance far from the starshade. In summary, we define the percentage difference with respect to the reference PSF's centroid as a fraction of the FWHM of the stationary reference PSF: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 5 7 8 The reference PSF FWHM values at mid-band are 40.6 and 58.9 mas for Roman's 425-to 552nm and 615-to 800-nm bands, respectively, including the effect of its pupil, and 25.7 and 37.2 mas for HabEx's 425-to 552-nm and 615-to 800-nm bands, respectively. In the following sections, we also show the same results for Δ LOG PSF ¼ log 10 ½absðΔPSFÞ to be able to see the differences when they are small or change significantly between different parameter values. We note that ΔPSF contains information about the sign of the differences, whereas Δ LOG PSF does not.
The analyses that we perform contain a large number of combinations that are all accessible in SISTER's website (see Ref. 45) together with a MATLAB script that can reproduce any of the results of this study. In particular, we show explicit results for two starshade mission concepts: the starshade rendezvous 32 and the starshade for HabEx, 12 although a very similar procedure would be applied to any other case. We analyze the parameters that constitute the fundamental building blocks of the scene simulation. The diffraction calculation for a given starshade orientation and source position is governed by the number of points at the diffracting aperture N PE and the number of points in the pupil, which is tied to the PSF spatial extent N λ∕D . The starshade spins about its axis, and this is represented by a finite number of closely spaced rotations n pos . The starshade diffraction pattern evolves rapidly with distance from the line of sight and with wavelength, especially within the IWA. The required spacings between source points and wavelengths to allow an accurate interpolation of the PSF at any source point and wavelength are given by Δr and Δλ. Specifically, these parameters are as follows: • the number of sampling points N PE along each edge of the starshade petal; • the number of starshade rotations n pos spanning the 15 deg angle between petals; • the radial spatial extent N λ∕D of the PSF expressed relative to the ratio of the wavelength and telescope diameter; • the spacing Δr between adjacent source points in mas; and • the spacing Δλ between adjacent monochromatic wavelength channels in nanometers.
We do not treat the number of pixels sampling the pupil N PIX separately because we have seen that given a PSF spatial extent N λ∕D it suffices to choose N PIX ≥ 2N λ∕D , the Nyquist sampling criterion. A lower value introduces aliasing effects on the image plane. Higher values of N PIX without increasing N λ∕D do not significantly improve the precision of the simulations, so we set N PIX ¼ 2N λ∕D þ 1, so there is a well-defined center in the pupil array.
When plotting the results, we show a summary figure to provide a quick overview of the maximum difference for each of the three cases as a function of the parameter value, at any distance from the starshade and plots for each of the three cases showing the radial dependence of the difference for each parameter value. Additionally, in Appendix C, we show actual PSFs and PSF differences for several parameter values.
Before dealing with each of the five parameters listed above, we begin by looking at two general effects. First, we show how the signal is affected in the case of a non-spinning starshade when the source is located at different positions between or behind the starshade petals. The result also applies to a spinning starshade and justifies the inclusion in the analysis of a large set of radial distances that extend well beyond the IWA. Second, since the starshade rendezvous has an obstructed pupil with a secondary and struts, we show how photometry and astrometry are affected by the struts in the pupil and demonstrate that, within our precision goal, a radial PSF basis is sufficient. We then derive the set of adequate parameters for a circularly symmetric pupil.

Relative Position of a Source with Respect to the Starshade Petals
This section shows the effect of having a point source located at some angle between two starshade petals (Fig. 14). This is an important instrumental effect for exoplanet imaging that also has a direct effect on how the PSF basis of a spinning starshade is generated. To be able to speak about the angular distance between a source and a starshade petal, the starshade must be nonspinning. The reference PSF at any radius is chosen as the one associated with a source located along a radial axis that goes from the starshade's center through one of its tips. The other cases have a relative angle between the source and the tips of the starshade α ss . We considered α ss ¼ 0 deg The results shown in this section were derived for HabEx only to avoid any confusion between the effects due to α ss and those from Roman's non-circular pupil (see next section). Since HabEx's starshade has 24 petals, the relative position between a source and the starshade petals becomes identical after a rotation of 15 deg. Finally, given the relevance of this effect, we analyzed cases spaced by 1 mas from on axis to a radial distance of 230 mas, well beyond the 70-mas IWA. Figure 15 shows a panel showing the results for one starshade mission and one passband. Table 2 provides a summary of the results for both starshade missions and passbands, whereas all of the individual results can be found in Ref. 45. The effect of the relative position of a source with respect to the tips of the starshade is noticeable for almost any value of 0 deg < α ss < 15 deg at distances close to the IWA. The differences with respect to the reference case become negligible close to the starshade axis and far from it. In particular, the differences become subpercent for r ≳ 125 mas, which is significantly farther than the IWA ¼ 70 mas. As mentioned before, these results have a direct implication in precise exoplanet imaging and when generating a PSF basis. As expected, the panel shows that, for α ss ¼ 15 deg, the differences are compatible with numerical precision. Also notice that, for some values of α ss , the curves seem to be missing. This is because the differences should depend on modðα ss ; 7.5 degÞ. In these cases, the two curves get overplotted, and only the last (greater) value of α ss can be seen.

Non-Circularly Symmetric Pupil
In the case of the Roman Space Telescope, the 2.36-m diameter pupil contains a secondary mirror with a diameter of 0.715 m 55 that covers 10.2% of the collecting area of the primary mirror. The secondary is supported by 6 struts with 120 deg symmetry 56 that cover a non-negligible 7.2% of the collecting area of the primary mirror; see Fig. 16. SISTER constructs the PSF basis along a single radial direction and then rotates the PSFs by an azimuthal angle α to account for the starshade's circular symmetry. When the pupil contains struts or some noncircular symmetry, this process rotates the pupil by α, erroneously simulating a telescope with non-symmetric sidelobes that rotate with for each azimuthal field angle.
In this section, we show that making the assumption that the pupil rotates with the telescope gives rise to subpercent errors, except when the source is very close to the starshade center, where it becomes ≤10%. Figure 17 shows a panel showing the results for the Roman starshade and a wavelength of 425 nm. Results for other bands are similar. The pupil rotation error only affects the fidelity of flight formation simulations and with 3% to 10% errors, whereas astrophysical imaging is more accurate (recall Fig. 13). Figure 17 shows that, at and beyond the IWA, pupil rotation errors lead to peak errors ≤1%, whereas photometric and astrometric errors are below 0.3%. In the future, we may optimize SISTER to deal with the actual two-dimensional response across the image plane, if more accurate results are necessary. As for HabEx, we consider a circular pupil since it has an unobstructed telescope. It also allows one to compare our results for both circular and non-circular pupils. Table 3 provides a summary of the results for both starshade missions and passbands, whereas all of the individual results can be found in Ref. 45.   16 Roman's pupil with its secondary and six supporting struts with a high resolution. The diameter of the secondary is 0.303 times the diameter of the primary aperture and covers 10.2% of the collecting area. The struts have a 120-deg symmetry and cover 7.2% of the collecting area. Table 2 Results of the precision tests for a non-spinning starshade with a PSF that is derived for different rotation angles α ss of its petals with respect to the nominal configuration. We show the results for one starshade loci, two passbands, and six sets of radial distances from the center of the starshade. Each cell triplet shows the log 10 of the average differences in percentage for the peak/photometry/centroid tests, respectively, where a reference basis PSF REF   Here we define the five modeling parameters mentioned in Sec. 3, beginning with N PE , the number of points sampling each petal edge from its base to its tip. A denser sampling of locus points should lead to a more accurate representation of starshade diffraction, but at the expense of computation time that increases linearly with the number of points. This parameter is at the core of any optical response, and we thus start our study with it, following the methodology described in Sec. 3.3. The total number of points used to define the locus of the starshade is

Number of Starshade Rotations that Represent a Spinning Starshade
Simulating images for a spinning starshade requires the consideration of different starshade positions as we illustrate in Fig. 19. Even if individual exposures with the detector are shorter than the time it takes one petal to move to the position of the next petal, the scientific data require total   static starshade, and the PSF response is calculated along a radial axis that goes from the center of the starshade through a petal's tip (see Fig. 19). The results derived on this section have N PE ¼ 4000 points along each petal edge.

PSF spatial extent
In this section, we pay attention to the spatial extent of the PSF that allows one to fulfill our precision goal. The PSF response to the pupils considered in our work is not band limited.
As it is well known, the PSF of an unobstructed, circular pupil is the Airy pattern. In this case, the encircled energy extends to distances far from the PSF's center. For small angles α ¼ Nðλ∕DÞ ≪ 1, where λ is the monochromatic wavelength under consideration and D is the diameter of the aperture, the fraction of the encircled energy with respect to the total energy EE is given by the expression EE Airy ðαÞ ¼ 1 − J 2 0 ðπNÞ − J 1 ðπNÞ 2 , where J 0 is the Bessel function of the first kind of order 0, and J 1 is the one of order 1. 50 For instance, 99% of the encircled energy of an Airy disk is contained within a radius of 21λ∕D, whereas 99.5% is contained within 40λ∕D. In practice, the PSF does not maintain its diffraction-limited behavior at large N because of aberrations and scatter in the system. One could thus think of setting the PSF extent to some value, for instance N ¼ 5 or similar, and rescale the PSF library accordingly. However, in this work, we are not interested in deriving a PSF response that is adapted to some practical postprocessing scenario. Our aim is to emulate as accurately as possible what will happen in a real experiment and provide optical simulations with subpercent differences with respect to an ideal reference case. In this work, we have two different types of pupils: (i) Roman: obstructed with a secondary mirror and the struts that support it and (ii) HabEx: unobstructed. Even in the case of HabEx, with a circular unobstructed pupil, we could not rely on the analytical result of EE Airy because the PSF arrays that we use in the simulations are square not circular. Hence, EE is slightly greater than in the circular case. And, even if the change in EE is small, the effect on N is noticeable. In summary, we derived the value of EE for different values of α ¼ Nðλ∕DÞ using the intensity of the diffraction pattern of an obstructed circular pupil, 57,58 numerically integrating it over a square grid, making sure the result converged for large values of N. We also verified that we recovered the analytical results if EE is integrated over a circle. We found that the encircled energy contained within a distance of 40λ∕D for a square aperture is 99.9% for the unobstructed case (HabEx) and 99.5% for the obstructed case (Roman). We also found that 99% of the encircled energy for a square aperture is attained at 8.4λ∕D for the unobstructed case (HabEx) and 22.0λ∕D for the obstructed case (Roman).
Our goal is to be able to simulate the PSF response with a subpercent precision. Thus to reach 99% of the encircled energy, the relative error with respect to the 40λ∕D case has to be 0.5% (1 to 99/99.5) in the case of Roman and 0.9% for HabEx. Fig. 19 Spinning starshade as viewed by SISTER. Since the detector integration times are significantly longer than the spinning period of the starshade (nominally 1/3 rpm), SISTER derives the optical response for each of a series of multiple starshades with petals that have been rotated by angle with respect to a nominal configuration (top left), which are incoherent and co-adds them. In the figures above, the number of starshade configurations n pos evenly spans the angle between two adjacent petals. Our analysis shows that, for an ideal starshade, n pos ¼ 7 provides subpercent precision compared with many more rotations.
We then simulated different PSFs with a spatial extension given by 2N λ∕D × 2N Figure 29. They show that 99% of the encircled energy for a square aperture of an unobstructed case (HabEx) is contained within a distance of 8.8 λ∕D, compared with 8.4 λ∕D from the previous numerical analysis. And, in the obstructed case of Roman, our simulations derive 21.5 λ∕D, whereas the numerical analysis derives 22.0 λ∕D. We consider this agreement satisfactory inasmuch as we choose slightly greater values for the adequate choice for the simulations to leave some additional margin: 10 λ∕D for HabEx and 23 λ∕D for Roman; see also Sec. 3.7. One also notice in our results the differences between obstructed and unobstructed telescopes.
Given a value of the PSF spatial extent, the pupil needs to be sampled with a number of pixels N PIX that avoids aliasing effects in the image plane: N PIX ¼ 2 Ã N λ∕D . We choose, however, N PIX ¼ 2 Ã N λ∕D þ 1 to have a well-defined center of the pupil's array. We have seen how N PIX < 2 Ã N λ∕D produces aliasing and how increasing N PIX without increasing N does improve the precision of the simulations.
In summary, we need N λ∕D ¼ 23 for Roman and N λ∕D ¼ 10 for HabEx to safely achieve a subpercent precision in the encircled energy of the telescope's PSF. Figure 20 shows the adequate sampling of Roman's pupil.

PSF spacing
This section shows the result of comparing a PSF basis for a spinning starshade with a PSF that is derived for point sources separated by some distance (Δr or "PSF spacing"). Figure 21 shows how the PSF changes at different distances from the starshade's axis. The reference basis PSF REF is built with Δr ¼ 0.5 mas, i.e., every half mas starting at the starshade's center. The other PSF responses PSF PAR have been derived with Δr ¼ 1;2; 3;4; 5;6; 7;8; 9, and 10 mas. For any radial distance that is not an exact multiple of Δr, the PSF response is approximated by five different interpolation methods. (i) Assigning the PSF response corresponding to the closest lower multiple of Δr: "previous," (ii) assigning the PSF response corresponding to the closest upper multiple of Δr: "next," (iii) assigning the PSF response corresponding to the closest multiple of Δr: "nearest," (iv) by means of a spatial linear interpolation between the PSF response corresponding to the closest lower and upper multiple of Δr: "linear," and (v) by means of a spatial spline interpolation between the PSF response corresponding to the closest lower and upper multiple of Δr: "spline." We also considered two other non-linear interpolations: Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) and the Modified Akima Spline method (Makima), which can handle very high changes or nearly constant behavior better than spline. In our case, PCHIP did not provide better results than spline and took about two orders of magnitude longer to interpolate a PSF basis with Δr ¼ 5 mas compared with the spline interpolation. Makima took even longer than PCHIP. We report here the best results that were obtained with the spline interpolation method.

Wavelength spacing
In this section, we study the PSF response of a spinning starshade across a passband with monochromatic wavelengths that are separated by some value Δλ. Figure 22 shows how the PSF response changes at four different monochromatic wavelengths in the case of starshade rendezvous and for two different source locations. Let us call this set of wavelengths the "pivot" wavelengths. For any other wavelength, we combine the PSF response at these pivot wavelengths into a new PSF response. Since the PSF changes with wavelength both in intensity and spatial extent, our analysis includes a peak intensity interpolation and a spatial interpolation. In the case of the peak intensity interpolation, we found that it is enough to consider two methods. (i) The PSF response at some wavelength is identified with the PSF response of the closest pivot PSF with a wavelength that is equal or higher than the wavelength under consideration. We label this case "upper" peak approximation. (ii) The two closest pivot PSFs are linearly combined. Specifically, the linear factor is given by the wavelength difference of the actual wavelength with the lower and upper pivot wavelengths. When merging both pivot PSFs, given the complexity of the spatial structure of the PSF close or inside the IWA, we considered a total of four spatial interpolation schemes, although they can be classified as two different approaches. (i) There is no modification of the spatial extent of the pivot PSFs, that is, the two pivot PSF responses are directly co-added and the only interpolation is the peak intensity one. (ii) The pivot PSFs are radially stretched/contracted with a factor that is directly proportional to the quotient between the wavelength under consideration and the lower and upper pivot wavelengths. We label this type of transformation "conformal." For the sake of completeness, we consider three specific interpolation methods when implementing the conformal transformation: bilinear, Lanczos 2, and Lanczos 3. (The conformal transformation and its methods were implemented in our analysis using the MATLAB function "imresize.") In total, we have 2 × 4 different interpolation schemes for each Δλ. To simplify the amount of calculation involved in this section, we did not consider the full extension of both passbands: 425 to 552 and 615 to 800 nm, but three representative subsets for each passband instead: a lower, a middle, and an upper sub-band. Specifically, we performed the analysis just described on the sub-bands: 425 to 435, 483 to 493, 542 to 552, 615 to 625, 702 to 712, and 790 to 800 nm. The "reference" basis PSF REF was built with Δλ ¼ 0.5 nm. The other PSF responses PSF PAR were derived for Δλ ¼ 1;2; 3;4; 5;6; 7;8; 9, and 10 nm. Each image has been normalized to its peak intensity (i.e., the color scale is linear from 0 to 1) (see also Fig. 13). The effects of diffraction are noticeable for sources inside the IWA, affecting both photometry and astrometry.
After performing all of these combinations, we conclude that the case that gave the least differences between PSF PAR and PSF REF was the linear interpolation for the peak intensity and no spatial interpolation. This is the case for which we report the results in Table 11.

Summary of Adequate Parameter Values
We performed extensive modeling in several iterations for the set of parameters described above. In a first iteration, we selected nominal values for four fixed parameters and studied the sensitivity to the fifth. In subsequent iterations, we reset the parameters to a value from the earlier iterations that led to subpercent photometric and astrometric errors.
The results are summarized in Table 4. See Appendix A for detailed tables of sensitivities to each parameter, plots showing the worst case sensitivities for each parameter as a function of working angle, and plots of the sensitivity for each parameter value. We performed the analysis described above for different values of the PSF pixel pitch: 1, 3, 9, and 12 mas. We did not  Table 4 Results of the precision tests for a spinning starshade. We show the results for two starshade loci, two passbands, and five model parameters. Each cell triplet shows the log 10 of the average differences in percentage for the peak/photometry/centroid tests, respectively, where the reference basis for each respective parameter is built from N PE ¼ 16;000 points, n pos ¼ 37 angles, N λ∕D ¼ 40, Δr ¼ 0.5 mas, and Δλ ¼ 0.5 nm. For instance, for Roman, the 425-to 552-nm passband, a radial distance of 65 to 75 mas from the starshade center, and a PSF built with N PE ¼ 4000 points along the edge of the petal compared with the reference built with N PE ¼ 16;000 points, the results in this table are −1.2∕ − 1.6∕ − 1.5, which correspond to averaged differences of ∼0.063%, 0.025%, and 0.0032%, respectively. According to Eq. (4), the difference in the centroid test for the chosen example corresponds to 1 μas. The on-axis response is the same for any Δr . We show the results for Δr ¼ 1 mas and 1 to 6 mas that provide a precise response near the axis. observe any significant deviation from conclusions in any of them. Thus, in practice, the PSF pixel pitch can be chosen to match the astrophysical scene's pixel pitch without loss of accuracy. These results apply to starshade designs used over a range of Fresnel numbers spanning ∼8 to 21, where the Fresnel number is defined as F ¼ r 2 ∕ðλZÞ for a starshade of diameter r, a wavelength of λ, and telescope-starshade separation Z. The range is set by the limits of the SRM and Habex bands shown in this table. The results will likely hold over a much larger range of Fresnel numbers, though we do not show this explicitly in this work. The parameters from Table 4 were used for the imaging simulations described in Sec. 4.  Table 4 that needs to be updated is the number of rotations used to derive the starlight response within the precision goal. SISTER automatically spins the perturbed starshade into different petal positions until the precision goal is reached. This value depends on each perturbation case and is usually 3 to 5 times n pos . The PSF response is then stored on disk, and the next simulation is as fast as the ideal starshade case. Of course, if desired, one can compute the full basis functions for any non-ideal starshade locus.

Starshade with Perturbation
Two other common perturbations are handled: formation flying and starshade tilts. Formation flying misalignments, 5,59 which are required to be <1 m radially, lead to leakage of starlight at the 10 −11 contrast level. SISTER allows the user to adjust the lateral and longitudinal positions of the starshade, treating all such cases as instantaneous. There is no alignment control mechanism in the code. Additionally, the starshade tilts about its horizontal axis. This case is handled by adjusting the axial component of each locus point and computing the on-axis contrast, just as with a spatially perturbed starshade. Figure 23 shows the nominal on-axis response, a tilted starshade, one with some lateral misalignment (1 m to the left), and a non-ideal starshade with a 10 −10 contrast, which is the result of some random fabrication errors along the petal edges with 400-μm RMS. This figure is derived for starshade rendezvous and 425 nm.

Solar Glint
Sunlight will glint, via both reflection and diffraction, from the shape-defining edges of the starshade petals. Both Lambertian 3,60-63 and specular 62-64 edge characteristics have been studied, with nearly razor-sharp specular edges generally appearing to be fainter overall. From measurements of coated edged and modeling of optimized designs, the expected visual magnitude of the average glint lobe brightness around the IWA is ∼30 for the SRM 425 to 552 band, ∼29 for the SRM 615 to 800 band, and ∼30 over the HabEx range of 400 to 1000 nm. 65 The location and brightness of the glint are a function of the position of the Sun and the orientation of the starshade. The Sun's position is defined by its angle from the starshade normal ϕ and its angle from observational East α, as shown in Fig. 24. The glint is localized to the regions where the sharp, specular edges are broadside to the Sun-starshade-telescope plane (Fig. 25). This results in a two-lobe pattern with some polarization dependence.
SISTER makes use of laboratory measurements of the glint from flight-like edge samples. 3,65 The measurements are analyzed in the form of scatter distribution functions (SDFs) as shown in Figure 25(a), where the abscissa represents the angle of the Sun α relative to the orientation of the edge θ with θ − α ¼ 0 defining the specular angle. The ordinate axis is the Sun's angle from the starshade normal ϕ. The SDF shows that the etched amorphous metal edges are highly specular Fig. 24 Orientation of the Sun with respect to the scene. The Sun is ϕ deg from the starshade normal and oriented α deg from the East axis. Fig. 25 The glint brightness along each small section of the starshade perimeter for a solar angle of ϕ and angle relative to the edge of θ − α is determined from the measured SDF. (a) The logscale SDF for an amorphous metal edge that has been coated with an antireflection coating. 65 A perfectly specular edge would appear as a single vertical line at θ − alpha ¼ 0. In this example, the Sun is located at ϕ ¼ 60 deg from the starshade normal and below the starshade. (b) Edges oriented perpendicular to the Sun-starshade-telescope plane appear brightest as shown in high resolution. (c) The local flux at the starshade edges is convolved with the telescope PSF and added to the imaged scene. and that the flux is brightest at low ϕ where diffraction dominates. Separate SDFs are measured for polarization parallel and perpendicular to the edge. SISTER reads the SDFs and the locus of points constituting the outline of the starshade. Given the position of the Sun relative to the starshade and telescope and the local normals of the starshade edges, the code determines where on the (θ − α; ϕ) SDF to sample the scatter. Edges oriented for specular reflection appear brighter, as shown in the middle panel of Fig. 25. The edges are categorized as "leading" and "trailing." Leading edges are defined as having the Sun-starshade-telescope angle <180 deg. The light can reflect and diffract into the telescope. Trailing edges require the light to diffract toward the telescope, with no possibility of reflection. The code assumes that these edges, as well as all edges out to a radius of 7.5 m, are shadowed using structure on the Sun-facing side of the starshade and do not contribute to the glint. The contributions of all remaining edge segments are convolved with the telescope PSF at each wavelength being simulated. The PSF accounts for the defocus of the image since the starshade is not distant enough to be in the far-field of the telescope. In the case of a spinning starshade, the petals are rotated to different positions, and the resulting scatter maps are averaged.

Solar glint model validation
To validate the accuracy of the SISTER calculations, we generated a test starshade with a triangular shape having a 1-m long edge oriented for specular reflection of Sunlight. We assumed that the edge was perfectly conducting and had a negligible terminal radius of curvature so that only diffraction was present. We applied the Sommerfeld diffraction equations 50 for the "S" (parallel) polarization and generated an SDF for a specular edge. This SDF was substituted for the experimental measurements. With the Sun positioned at ϕ ¼ 80 deg (10 deg behind the plane of the starshade), SISTER then calculated a glint lobe. We calibrated the magnitude of the glint lobe by modeling a planet of known magnitude. We used SISTER to generate a planet well to the side of the starshade (so as not to be attenuated by the starshade). We chose the host star to be the Sun at a distance of 10 pc (V ¼ 4.83) and chose the planet to have a flux ratio of 8.55 × 10 −9 (delta mag = 20.17). This created planet with a visual magnitude of V ¼ 25.
Comparing the planet with the glint lobes, we found that the integrated light of the glint lobe was V ¼ 24.63.
We then compared this with the analytical calculation for the 1-m edge at the distance of the starshade. Assuming only that the Sun had a visual apparent magnitude of −26.76, 66 we obtained the predicted glint result of V ¼ 24.64. We thus conclude that errors in the SISTER glint code are ΔV ¼ 0.01 mag or 1%.

Detector
SISTER implements the same general detector model used in Roman's performance studies. 67 It allows one to deal with the main sources of noise: read noise, clock induced charge, dark current, and shot noise and to correctly handle the gain of an EMCCD detector. From a noiseless, optical SISTER simulation, SISTER derives a realization of the detector model at each detector pixel producing a noise image that can be added to the noiseless simulation. In this way, SISTER may quickly produce several independent noise realizations from the same input simulation that can be used in the statistical analysis of starshade images. However, effects such as cosmic rays, brighter-fatter response, or correlated noise are not included in SISTER. Since SISTER provides the simulated scenes without noise, one can run any more elaborated detector noise simulator through them. (For instance, the starshade data in the Roman Exoplanet Imaging Data Challenge Ref. 44 were generated from SISTER noiseless simulations to which we added the output noise from the EMCCD noise simulator from Ref. 68.) Section 4 shows an example.

Limitations and Future Development
Although solar glint is expected to be a significant source of instrument noise, there are other known sources of scattered light that are under development but have not yet been fully implemented in SISTER. These include reflections of bright bodies, e.g., Earthshine, and the Milky Way; leakage of starlight and Sunlight through micrometeoroid holes; glint of Sunlight through gaps in the starshade edges; and solar illumination of exhaust plumes from the stationkeeping thrusters. The expected brightness of these sources is addressed by Hu et al. 69 The plan for bright body reflections is to model the reflection and scatter properties of telescope-facing surfaces, determine the starshade orientation in space, and reflect the light from any significant bright bodies in the starshade-telescope hemisphere into the telescope. The image includes a defocus term representative of the distance between the telescope and starshade. The telescope-facing surface is represented with a three-part model: the nominal surfaces, which include tilted, flat petals and a conical interior section that covers the spacecraft; a model of the surface ripples in the optical shield; and a model of the bidirectional reflectance distribution function of the surface material. The bright bodies include the crescent Earth, crescent Moon, Venus, Jupiter, and a spatially and spectrally coarse model representing the Milky Way. The defocus term is already present in SISTER and is used in the solar glint model.
The optical shield is designed to have a high enough optical density to restrict Sunlight leakage to a negligible level. However, micrometeoroids impact and perforate both sides of the starshade. Although we estimate that over a 5-year mission this will be a minor effect, 69 we nonetheless plan to include a model that accounts for the nature of the hole pattern, which is expected to be a small penetration hole followed by two larger holes formed from the cone of material vaporized with the impact in the first layer. The model will account for diffraction and reflection at the holes, scattering between the layers, and a final scattering out of the telescopefacing holes.
Approaches to modeling solar glint through edge gaps and solar illumination of exhaust plumes are still being evaluated.

Example
Using the nominal SISTER modeling parameters listed in Table 4, we simulate the imaging of a scene containing a nearby star with three planets and an exozodical disk. The starshade is assumed to have a lateral displacement of 1 m, shape imperfections, and solar glint. The simulation also includes an extragalactic background object and local zodiacal light. We choose β CVn 31 as the reference star and a 1-day long integration using the 615-to 800-nm band of the SRM mission. β CVn is a G0V star of absolute V ¼ 4.64 at 8.44 pc, with a mass of 1.02 solar masses and a luminosity 1.15 times that of the Sun. Its proper motion is significant: ð−704.75; 292.74Þ mas/year and would help disambiguate exoplanet candidates from other background objects.
We placed three planets in circular orbit around the star: an Earth-twin at 1 AU, a sub-Neptune having a Neptune-like albedo at 2.28 AU, and a Jupiter-twin at 5.2 AU. All are assumed to have a Lambertian phase function with semimajor axis and phase angle shown in Table 5. The system meets the dynamical stability criterion of Eqs. (23) and (35) in Gladman. 70 Even though a recent analysis based on radial velocity ancillary data has ruled out the presence of a Jupiter-like planet within 10 AU of β CVn, 33 we include it here at 5.2 AU for the sake of showing the relative intensity of these three major types of planets at three representative distances from the host star. The system is inclined at 35 deg. The planets' apparent separations in Table 5 are determined by the system inclination, planet phase angles, and SMA. The simulation uses nominal SRM instrument parameters consistent with end-of-life (EOL) operation. 71 The optical throughput is the product of reflection losses in the telescope (0.81), in the coronagraph instrument (0.60, exclusive of coronagraph masks, which are not used), and the starshade dichroic filter (0.9) for a net throughput of 0.44. The telescope pupil includes the secondary mirror central obscuration and struts and is 2.36-m in diameter.
The starshade is assumed to have a 1-m lateral formation flying displacement (apparent motion is 0.8 mas) along the positive declination axis. This places the starshade at the limit of the required alignment tolerance.
The starshade is assumed to be imperfect; to simulate possible manufacturing and deployment errors, we displace the petals in-plane as 24 independent rigid bodies with radial and azimuthal offsets such that the root-mean-square displacements of all points along the petals is 0.56 mm. This creates a starshade that under ideal alignment has an instrument contrast of 1.4 × 10 −10 at the IWA. The contrast is substantially worse when combined with the 1-m offset.
Solar glint is computed under the assumption that the starshade edges are coated with an antireflection coating similar to one that has been tested in the laboratory. 65 The Sun is 60 deg from the starshade normal, which is a median value among planned observations, and is oriented at 45 deg from horizontal.
For the extragalactic background object, we used the Haystacks data cubes 25 and selected the brightest background galaxy from the HST deep field survey.
The exozodiacal dust cloud is generated using ZODIPIC as described in Sec. 2.6. β CVn is a solar type star, and the debris disk from ZODIPIC, which is based on solar system data, is consistent with what might be the actual debris disk in β CVn. The forward scattering function uses a value of g ¼ 0.2, which is an average value found in the solar system and nearby systems. 72 The cloud is assumed to have a density 5 times that of the solar system zodiacal cloud. However, we do not include any particular structure (dust rings or clumps) produced by the three-planet system. This is a complex task in itself that would require intensive N-body simulations and is out of the scope of SISTER. As mentioned in Sec. 2.6, if the user has such a dust template at different wavelengths, SISTER will simulate the corresponding image.
The detector noise follows the model described in Sec. 3.10 for its EOL. The particular values of the different parameters come from on-line Roman Space Telescope parameter tables. 73 The readout noise is 100e − per read out frame, EM gain is 1000, dark current at EOL is 0.77 e − ∕pix∕h, and clock induced charge noise is 0.02 e − ∕pix∕frame. The detector quantum efficiency is based on laboratory measurements of Roman's detector and is included in SISTER's distribution. We also include a current best estimate, a reduction of 25%, for its degradation at EOL. The pixel scale in the Roman coronagraph instrument is 21.85 mas∕pixel. We assumed that all pixels are identical; there are no hot pixels. The simulation does not include contamination from cosmic rays or sources of noise not intrinsic to the detector electronics.
Simulated images are shown in Fig. 26. Each subpanel shows a separate component, except the lower right panel, which is the combined scene consisting of the laterally displaced, shapedistorted starshade, solar glint, background galaxy, planets, exozodiacal light, local zodiacal background, shot noise, and detector noise. The average local zodiacal background has been subtracted to keep the combined image on the same scale as the other subpanels. The brightness scale is the logarithm of the average number of detected photons (the actual number in the combined image) for a 1-day long observation. A circle of 24 blue dots at the center of each subpanel indicates the radial extent of the starshade petal tips (104 mas), corresponding to the IWA. Table 6 shows the peak counts at any position and peak counts at the IWA for each component.
The planet contrasts in Table 6 do not exactly match the planet flux ratios in Table 5 because we are comparing peak counts rather than integrated counts and because the starshade partially attenuates planets that are close to the IWA, recall Fig. 13, including the exo-Earth at 97 mas. At this location, the transmittance is ∼0.64, accounting for the contrast and flux ratio difference.
This simulation demonstrates that bright exozodiacal disks are the dominant sources of background. As noted above, the perturbed starshade combined with the lateral offset diffracts significant light at the IWA, well above the diffraction expected from an aligned starshade. This is a reminder of the importance of these simulations in capturing effects that could be missed when considering diffraction, reflection, and other scattering effects in isolation.
We have described the starshade imaging simulation toolbox known as SISTER. The main conclusion from this paper is that we have identified the set of model parameters, shown in Table 4, resulting in subpercent errors of the PSF modeling over the image plane. Although higher accuracy is possible using different values of the model parameters, exoplanet observations will be challenged with noise that, in the majority of the situations, will be larger than the SISTER modeling precision. The complete set of sensitivity tables in Appendix A serves as a guide to the accuracy obtainable for reduced or greater parameter values. The SISTER code and documentation are open source and available in Ref. 27.

Appendix A: SISTER Parameter Sensitivities
In this appendix we show the numerical results for the starshade modeling parameters N PE , the number of points sampling each starshade petal's edge, n pos , the different number of positions of  the starshade that emulate a spinning starshade, N λ∕D , the spatial extent of the PSF basis in units of λ∕D, Δr, the angular distance between two consecutive PSF elements of the PSF basis, and Δλ, the wavelength spacing between two consecutive monochromatic PSF responses. First, Figure 27 shows a panel with the results for N PE , for one starshade mission and one passband. Table 7 provides a summary of the results for both starshade missions and passbands, whereas all of the individual results can be found in Ref. 45. When analyzing the results, it is important to recall the value of the transmittance at different distances from the starshade's center (Fig. 13). Taking into account the results of the analysis, a value of N PE ¼ 4000 is enough to provide optical simulations with differences that are subpercent with respect to the reference PSF response, or with even higher values of N PE . This conclusion applies to the case of a nonspinning starshade. Later, we confirm that the same choice of N PE ¼ 4000 is adequate for the case of a spinning starshade. Figure 28 shows a panel with the results for n pos , for one starshade mission and one passband. Table 8 provides a summary of the results for both starshade missions and passbands, whereas all of the individual results can be found in Ref. 45. Even though n pos ¼ 4 already shows differences that are subpercent with respect to the reference case, we choose n pos ¼ 7 as the adequate value of this parameter because it clearly shows differences well below the percent level at all relevant distances, helping to ensure that the combination with other instrumental parameters is still subpercent while still providing reasonable times for the generation of the PSF basis. Figure 29 shows a panel with the results for N λ∕D , for one starshade mission and one passband. Table 9 provides a summary of the results for both starshade missions and passbands,  Table 7 provides a numerical summary for Roman, HabEx, and the 425-to 552 nm and 615-to 800-nm passbands. If the curves in the summary plot become smaller as N PE increases, one can conclude that the PSF response converges to values with a difference that is smaller than the results obtained in our analysis. Table 7 Results of the precision tests for a non-spinning starshade with a PSF that is derived with different samplings along the edge of the petal N PE with respect to the nominal configuration. We show the results for two starshade loci, two passbands, and six sets of radial distances from the center of the starshade. Each cell triplet shows the log 10 of the average differences in percentage for the peak/photometry/centroid tests, respectively, where a reference basis PSF REF built with N PE ¼ 16;000 is compared with four other PSF responses with N PE ¼ 1000, 2000, 4000, and 8000 points along the edge of the petal. For instance, for Roman, the 425-to 552-nm passband, a radial distance of 65 to 75 mas from the starshade center, and a PSF built with N PE ¼ 1000 points along the edge of the petal compared with the PSF REF , the results on this table are −1.4∕ − 1.5∕ − 2.3, which correspond to averaged differences of ∼0.040%, 0.032%, and 0.0050%, respectively. According to Eq. (4), the difference in the centroid test for the chosen example corresponds to 2 μas. The rows in bold face correspond to the parameter choice to build the standard PSF: 4000 points along the edge of the petal.  Table 8 provides a numerical summary for Roman, HabEx, and the 425-to 552-nm and 615-to 800-nm passbands. If the curves become smaller as n pos increase, one can conclude that the PSF response converges to values with a difference that is smaller than the results obtained in our analysis.  whereas all of the individual results can be found in Ref. 45. When analyzing the results, it is important to recall the value of the transmittance at different distances from the starshade's center (Fig. 13). Let us remark that for each value of N λ∕D , the total energy of the PSF changes because the spatial extension is different, but the spatial features do not change. This is why the centroid               remains the same for different values of N λ∕D , but the photometry and peak tests are "flat." The results show that for Roman and N λ ≥ 13, the differences are subpercent. For HabEx, N λ ≥ 7 provides differences that are just subpercent. Since instrumental parameters are set altogether, we choose N λ ¼ 16 such that the differences are ≲0.1% for both Roman and HabEx. Figure 30 shows a panel with the results for Δr, for one starshade mission and one passband. Table 10 provides a summary of the results for both starshade missions and passbands, whereas all of the individual results can be found in Ref. 45. Figure 31 shows a panel with the results for Appendix, for one starshade mission and one passband. Table 11 provides a summary of the results for both starshade missions and passbands, whereas all of the individual results can be found in Ref. 45. 7 Appendix B: Photometric Notes 7.1 Spectral Irradiance Units SISTER scenes must be given in terms of spectral irradiance, also named flux density. There are two ways of expressing it with units. One option is to provide the spectral irradiance per unit wavelength, and the other one is per unit frequency. In the first case, a standard unit is W∕m 2 ∕μm, or with other choices in the wavelength unit, e.g., nanometers. In the case of choosing the frequency, the common choice is Janskys (1 Jy ¼ 10 −26 W∕m 2 ∕Hz), a non-SI unit. This is typical practice in visible astronomy, 25 e.g., the spectral irradiance of the Sun 36 at 10 pc from Earth and at 425 nm is 23.29 Jy, corresponding to 3.78 × 10 −10 W∕m 2 ∕μm. The conversion between both systems of units is provided by the relationship F λ dλ ¼ −F ν dν, where the minus sign reminds us that increasing frequency corresponds to decreasing wavelength and viceversa λν ¼ c.

Solar Spectrum
In the last decades, there have been different compilations of the solar irradiance. SISTER uses by default the one prepared by Willmer 36 in 2018 (Willmer18). This data analysis sets the solar apparent magnitude to −26.76 in the V-band and its absolute magnitude to V ¼ 4.81. However, since other solar data have been used quite often and are still used, we provide the option of choosing another data set. 35 The other two common data sets are the one prepared by the American Society for Testing and Materials 74 in 2000 (AM0) and the one from Wehrli 75 compiled in 1985 (WMO85). Figure 32 shows the three with a comparison of their ratios for the range of wavelengths covered by SISTER.

Local Zodical Light
The units of the local zodiacal light are the standard V magnitudes (λ ¼ 551 nm and Δλ ¼ 88 nm) per square arc sec. The pixel size of the camera's detector is taken into account when deriving the contribution on the simulations (e.g., 21.85 mas for Roman 43 or 11.7 mas for HabEx 12 ). For some lines of sight, the local zodiacal light is considered too bright to carry any observations. In those cases, SISTER will generate a black image. If SISTER is dealing with Keplerian orbits at different epochs, it is expected that for some of them the simulation will be a black image, whereas for other epochs it will be a regular simulation. When translating the surface brightness into actual spectral flux, SISTER uses the solar spectrum as explained in Sec. 2.3.
The equation used to derive the photon rate per camera's pixel at the telescope's aperture in Fig. 10 is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 5 9 7 r γ ¼ 10 −0.4ΔV LZL × pix mas 1000 where r γ is measured in photons/s/pix, ΔV LZL ¼ 26.76 − V LZL is the magnitude difference between the apparent solar magnitude and the surface brightness of the local zodiacal light V LZL in units of V mag∕arc sec 2 , pix mas is the side length of the camera's pixel, SIS is the solar irradiance spectrum, 36 and ðλ 0 ; λ 1 Þ are the lower and upper limits of the passband under consideration.
We finish this section reminding the reader that the contribution of the local zodiacal light is twofold. On the one hand, it adds a constant background to the total number of counts of the simulation, which can be estimated by computing the median value of the output simulation in a region far enough from the host star and subtracted off. On the other hand, the local zodiacal light contributes to the noise level with a Poisson distribution term.

Appendix C: Individual PSF Differences
In this section, we show a few examples of the individual PSF differences between the reference PSFs used for each instrumental parameter studied in Sec. 3. These PSF differences together with Eqs. (2)-(4) provide all of the results on the parameter sensitivities. There are over 3 million such images that have been analyzed individually. Each panel in Figs. 33 and 34 shows the PSF differences for the starshade rendezvous mission and the monochromatic wavelength 425 nm. The source is (i) at an angular distance of 5 mas from the starshade axis and (ii) at 65 mas from the starshade axis, which is slightly inside the IWA, 72 mas for one specific choice of the instrumental parameter.