Integration time adjusted completeness

Abstract. Future, large-scale, exoplanet direct-imaging missions will be capable of discovering and characterizing Earth-like exoplanets. These mission designs can be evaluated using completeness, the fraction of planets from some population that are detectable by a telescope at an arbitrary observation time. However, the original formulation of completeness uses instrument visibility limits and ignores additional integration time and planetary motion constraints. Some of the sampled planets used to calculate completeness may transit in and out of an instrument’s geometric and photometric visibility limits while they are being observed, thereby causing the integration time agnostic calculation to overestimate completeness. We present a method for calculating completeness that accounts for the fraction of planets that leave the visibility limits of the telescope during the integration time period. We define completeness using the aggregate fraction of an orbital period during which planets are detectable, calculated using the specific times that planets enter and leave an instrument’s visibility limits and the integration time. To perform this calculation, we derive analytical methods for finding the planet-star projected separation extrema, times past periastron that these extrema occur, and times past periastron that the planet-star projected separation intersects a specific separation circle. We also provide efficient numerical methods for calculating the planet-star difference in magnitude extrema and times past periastron corresponding to specific values Δmag. Our integration time adjusted completeness shows that, for a planned star observation at 25 pc with 1-day and 5-day integration times, integration time adjusted completeness of Earth-like planets is reduced by 1% and 5% from the integration time agnostic completeness, respectively. Integration time adjusted completeness calculated in this manner also provides a computationally inexpensive method for finding dynamic completeness—the completeness change on subsequent observations.


Introduction
Direct imaging blind search mission schedules can be optimized 1 by maximizing completeness 2the fraction of exoplanets from an assumed planet population that are detectable by a particular instrument at an arbitrary observation time. Completeness is typically parameterized by a limiting planet-star brightness difference (Δmag lim ), the inner working angle (IWA) of the instrument, and the outer working angle (OWA) of the instrument. The original Monte Carlo approach developed by Brown 2 involved creating a cloud of synthetic planets by sampling the underlying Keplerian orbital elements (KOE) and physical parameters of a planet population and determining the fraction of those individually simulated planets within the visible limits of the instrument. Multiplying completeness by the exoplanet occurrence rate gives the expected exoplanet yield for observing a given star. Although completeness is a good metric for predicting instrument performance, the calculation described above only captures an instant in time and does not include whether the time that a planet is within the visible region of the telescope is sufficient to actually make a detection.
The original method of calculating completeness was developed by Brown. 2 Our adaptation of this method in EXOSIMS 3 is extensively outlined in Ref. 1 and involves randomly generating planet KOE, planet geometric albedo (p), and planetary radius (R) sampled from NASA's Exoplanet Program Analysis Group Study Analysis Group 13 (SAG13) 1,4,5 probability density functions, generating a finely binned 2D histogram of Δmag versus s (planet-star separation), and fitting a 2D spline to the histogram bins. In all current methods, completeness is calculated as the double integral of s and Δmag over the joint probability density function as in Eq. (17) of Ref. 1 and Eq. (7) in Ref. 6. The limits of integration define a detectable planet as one in which Δmag < Δmag lim and d i IWA < s < d i OWA, where d i is the distance of the host star from the spacecraft. While the Δmag lim in this paper is used to describe a general upper limit of integration for calculating completeness, this Δmag lim can be formulated as a function of integration time as in Eq. (12) of Ref. 1 by making assumptions about a multitude of instrument parameters and external noise sources. This approach, as conventionally implemented, is a good estimator for completeness, 1 but it does not take planet motion in time into account.
Searches like the Gemini Planet Imager and Nancy Grace Roman Space Telescope (Roman) Coronagraphic Instrument (CGI) 7 are sensitive to larger planets with larger planet-star separations and longer orbital periods. 1 They make use of Brown completeness to plan blind searches. However, future telescope concepts like the Habitable Exoplanet Observatory (HabEx) 8 and the Large UV/Optical/IR Surveyor (LUVOIR) 9 seek to find smaller Earth-like exoplanets, with shorter orbital periods around stars farther away. Their search criteria will result in Brown completeness based yield overestimations due to planet motion out of instrument visibility limits. Figure 1 shows a schematic of a direct imaging observation; it demonstrates how a planet can be within the photometric and geometric visibility limits of the instrument and still not be detected. The red regions of the planet's orbit indicate where the planet is within the visible limits but not detectable because the planet will move out of the region in less time than it takes to detect the planet [case (b) in this figure]. Since there will be some Earth-like exoplanets that are counted toward the completeness score but are not actually within the visibility constraints of the (c) (a) (b) Fig. 1 Schematic of a direct imaging observation. The line represents the projection of a planet's orbit about its parent star (yellow) into the plane of the sky as seen by a distant observer. The arrow indicates the planet's direction of motion. Blue portions of the orbit (a) indicate times when the planet is detectable. Red portions of the orbit (b) indicate times when the planet is within the visibility limits of the instrument, but it is not detectable due to integration time constraints. The dashed portions of the orbit (c) indicate times when the planet is too faint to be observed. The shaded gray regions represent the projected inner and OWAs of the instrument. The planet is unobservable on portions of the orbit intersecting these gray regions. Green dots indicate transitions in and out of instrument visibility limits. instrument long enough to be directly imaged, we need a new method that only counts targets within the instrument visibility long enough to be observed or characterized. Figure 2 demonstrates this phenomenon in the separation versus Δmag phase space commonly used to define instrument contrast curves and completeness. Figure 2(b) shows a shortperiod, fast-moving, Mars-like planet that transitions into and out of the assumed instrument's visible limits before it can be detected. Figure 2(c) shows the effect for a Uranus-like planet that crosses through and exits the visible region in less time than it takes for it to be detected. In these cases, Brown's completeness calculation would include simulated planets in these regions, thus overestimating the overall detection yield. Additionally, Fig. 2(d) shows that the last possible moment a Neptune-like planet could be imaged before it leaves the visible limits of the telescope. This demonstrates a portion of the orbit where the planet is within the visible limits of the telescope but is not detectable.
This paper presents a method, implemented in the EXOSIMS modeling software, 3 for calculating integration time adjusted completeness that allows us to investigate exactly how much the original completeness definition overestimates planet yields. EXOSIMS is an exoplanet direct imaging mission modeling software used to simulate the Roman CGI 7 and future mission concepts such as HabEx 8 and LUVOIR. 9 EXOSIMS simulates populations of planets, observatories, instruments, and underlying dynamics to create full, end-to-end mission simulations. We derive distributions of potential mission yield from ensembles of these mission simulations. EXOSIMS utilizes completeness as a heuristic for target selection and observation scheduling, 1 implementing both Brown's Monte Carlo approach to completeness 2 as well as the analytical completeness formulation from Ref. 6. However, EXOSIMS determines whether actual detections occur by evaluating the effective signal-to-noise ratio (SNR) of individual simulated observations on simulated planets. Keithly et al. 1 discussed the full optimization process for a singlevisit blind search and showed that the ideal yield from an EXOSIMS simulation ensemble (when mission constraints are discounted) for a limited blind search mission of Roman is equivalent to the yield calculated via Brown Completeness. Comparisons between Brown completeness-based exoplanet yield estimations and Monte Carlo of full mission simulations have been done in Ref. 10 and the Standard Definitions and Evaluation Team Final Report. 11 The yields calculated by completeness are higher than the average yield from a Monte Carlo of full missions simulated in EXOSIMS as noted in Table 6 of Ref. 11. These works attribute the difference in yields to the additional mission constraints captured by the full mission simulation approach, inefficiencies in the optimization algorithms, and inefficiencies in the scheduling algorithms implemented in EXOSIMS. 12 The unmentioned assumption is that Brown completeness-based yield estimates are accurate, but this paper demonstrates how the omission of planet motion and integration time reduces the resulting completeness yield estimates.
To know the amount of time that an individual planet spends within the observable limits of an instrument, we need to calculate the times in a planet's orbit when a planet enters or exits these limits. This means that we need methods for calculating when a planet has a given planetstar separation and Δmag. Once we know all of the times when a planet enters or exits an instrument's visibility limits, we know the fraction of time that the planet is able to be detected with that instrument. Averaging over these visibility fractions for a large number of samples evaluates completeness in a fundamentally different way from its original formulation. These same methods also allow for the evaluation of integration time adjusted completeness as well as a new method of calculating dynamic completeness. 13 Section 2 presents the derivation of this chain of calculations and discusses practical aspects of their implementation. Section 3 provides detailed validation of the methods and presents results of the various calculations enabled by them. Finally, in Section 4, we discuss various aspects of the algorithms and results, and we lay out future applications for this methodology.

Methods
In this section, we present our detailed process for calculating integration time adjusted completeness. We do this by finding the time windows in which a planet is within the separation and Δmag visibility limits of an instrument, discounting each time window by the integration time. The general overview of this process is as follows.
1. Calculate locations of apparent intersections between the projection of the 3D orbit into the plane of the sky and the s WA circle about the star. 3. Calculate ν from X and Y for each intersection and extrema. 4. Calculate t from ν of each intersection and extrema. 5. Combine times of s WA and Δmag lim intersections to create time windows between intersections. 6. Identify time windows in which the planet is visible or not visible. 7. Calculate integration time adjusted completeness averaging the orbital fraction of time that a planet is visible discounted by the integration time.
The general equations for Δmag and s used in these derivations are 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 ; 5 7 5 Δmag ¼ −2.5 log 10 p R j̱ r k∕i j 2 ΦðβÞ

;
(1) and 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 ; 5 1 7 Here R is the planet radius, p is the geometric albedo of the planet, and ΦðβÞ is the planet phase function. The other variables are defined in Fig. 3, where β is the star-planet-observer angle (also called the phase angle), r k∕i is the vector from star i to planet k, and ̱r i∕SC is the unit vector from the spacecraft (SC) to the star (̱r is the unit vector of r). The plane of the sky for a given observation lies in the ̱x and ̱ŷ plane, where ̱r i∕SC defines ̱ẑ of the target system. For our purposes, the direction of ̱x is arbitrary, but is typically taken to be a well-defined, inertially fixed direction, such as the ICRS mean equinox or pole direction. Here and throughout this paper, i in a subscript refers to the i'th target star and i (not subscripted) refers to the orbit inclination.

Projected Orbit and Separation Intersection
In this section, we derive an analytical expression for the true anomaly (ν) of s-orbit intersection points between a circle in the plane of the sky and the projection of a 3D Keplerian orbit on the plane of the sky. A practical example of a circle in the plane of the sky is the projected inner or Fig. 3 The orbital path of planet k (black) in a general X Y Z Cartesian coordinate system. r k ∕i (straight blue arrow) describes the location of the planet k (blue circle) relative to star i (yellow circle). The plane of the sky is noted by the translucent red parallelogram entirely in thex − and y − plane. The dashed ellipse (red) is the projection of the planet orbit (black) onto the plane of the sky.ĥ − is the orbit angular momentum vector.ê − is the orbit eccentricity vector.n − is the line of nodes. The star Ã is generally referred to in subscripts as Ã and often referred to in subscripts as the i'th star. This is not to be confused with the variable i, which is the inclination of the planet orbit. ν is the true anomaly of the planet. ω is the argument of periapsis of the planet. Ω is the longitude of the ascending node of the planet.
OWA of the instrument, equal to IWAd i (or OWAd i ) for target star distance d i . An orbit defined by KOE takes the shape of an ellipse in the orbital plane. The perpendicular projection of a 3D elliptical orbit onto the plane of the sky is given by the ̱x and ̱ŷ components of the 3D orbit. However, the resulting expression is not easily solvable for the true anomalies at intersection points. Instead, we can express the orbit projection in the plane of the sky as another ellipse. This simplifies the difficult intersection problem into analytically solvable subproblems and yields a relatively simple solution to the intersection between a circle in the plane of the sky and the projected orbit in the plane of the sky. Here we include the brief outline of this process, with the full procedure detailed in Appendix B.
1. Project the 3D orbital ellipse into a 2D projected ellipse and formulate the separation equation. 2. Find the expected number of s intersection points. 3. Calculate s intersection point coordinates.

Projection of the elliptical orbit
Given the KOE of a 3D orbit and a plane to project it on, we calculate the semimajor and semiminor axes of the projected orbit (a p and b p , respectively) as well as the angle between the projected semimajor axis and ̱x, which we call (θ). We first convert the KOE into orbital radius components in the Cartesian XYZ coordinate system as in Fig. 3: 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 ; 4 8 0 X ¼ rðcosðΩÞ cosðω þ νÞ − sinðΩÞ sinðω þ νÞ cosðiÞÞ Y ¼ rðsinðΩÞ cosðω þ νÞ þ cosðΩÞ sinðω þ νÞ cosðiÞÞ Z ¼ r sinðiÞ sinðω þ νÞ; (3) where r is the orbital radius (magnitude of r k∕i ) given by 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 ; 4 0 0 Ω is the longitude of the ascending node, ω is the argument of periapsis, ν is the true anomaly, and e is the eccentricity of the planet's orbit. Figure 4 shows a schematic view of the projection of the orbit onto the plane of the sky. There are two particularly important points associated with the projected orbit ellipse. First, F is the filled focus (star location) of the orbit and retains the same coordinates in the plane of the sky. When observing a star, F is the center of all working angle circles and the origin of the XYZ coordinate system. Note that in Fig. 4 the projected ellipse is located well below the original 3D ellipse for clarity, but the points F indicated by the orange circle and F indicated by the orange × are in fact coincident.
Second, O is the geometric center of the 3D orbit and is found most efficiently by averaging the XYZ locations of the planet at apoapsis and periapsis: 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 ; 2 2 3 where FO is the line segment (equivalently Euclidean vector) from F to O and r k∕i ðνÞ is the evaluation of Eq. (3) for the given value of ν. As ̱x and ̱ŷ define the plane of the sky, the XYZ coordinates of the geometric center of the projected ellipse (O 0 ) are given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 1 4 3 O 0 ¼ hFO · ̱x; FO · ̱ŷ; 0i: Appendix J provides a proof that any generic 3D ellipse projects to a 2D ellipse, as shown graphically in Fig. 4. The projection of the semimajor axis and semiminor axis from the 3D orbit onto the 2D plane of the sky form conjugate diameters of the projected ellipse. Any two diameters of an ellipse are conjugate diameters if and only if the tangent line to the ellipse at an endpoint of one diameter is parallel to the other diameter. 14 Each pair of conjugate diameters of an ellipse has a corresponding tangent parallelogram, sometimes called a bounding parallelogram. In The semimajor and semiminor axes of the projected ellipse may be found from any conjugate diameters. We take the conjugate diameters A 0 B 0 and D 0 C 0 and draw line QQ 0 through B 0 , perpendicular to D 0 C 0 as shown by the gray line in Fig. 5. Points Q and Q 0 are chosen such that The principal axes IR and ST lie on the bisectors of the angles formed by lines O 0 Q and O 0 Q 0 . From this construction, we calculate the projected ellipse semimajor axis, semiminor axis, and angular offset of the semimajor axis from ̱x. 16 IR and ST are given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 2 0 8 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 1 6 7 Defining the angle between O 0 B 0 and O 0 D 0 as ϕ, the cosine rule applied to triangles O 0 B 0 Q and O 0 B 0 Q 0 yielding E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 1 2 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 8 6 jO 0 Q 0 j 2 ¼ jO 0 B 0 j 2 þ jO 0 D 0 j 2 þ 2jO 0 B 0 jjO 0 D 0 j sin ϕ; where jO 0 B 0 j denotes the length of O 0 B 0 . Inserting these into Eq. (7), we obtain The planet orbits about the star, which is located at focus F (orange circle). The orange dot and orange × denoted as F are the same point in space, but the red projection is shown as offset from the original orbit for clarity. The projection of the original 3D elliptical orbit onto the X Y plane is given by the red ellipse containing points A 0 , B 0 , C 0 , and D 0 (blue diamonds; perpendicular projections of A, B, C, and D).
Point O is the geometric center of the 3D elliptical orbit and projects to O 0 in the plane of the sky. Point P (green circle) is any arbitrary point along the original 3D ellipse and maps to the semimajor axis and semiminor axis components H and K , respectively (green ×). P 0 (pink circle) is the perpendicular projection of P, and H 0 and K 0 (pink ×) are projections of H and K , respectively. The components of these perpendicular projections preserve the ratios of their values to the semimajor and semiminor axes, which, given the equation for an ellipse, can be used to prove the projection of an ellipse is itself an ellipse. The blue lines A 0 B 0 and C 0 D 0 form conjugate diameters of the red ellipse and are the projection of the semimajor axis and semiminor axis of the original 3D ellipse onto the plane of the sky. These conjugate diameters can then be used to find the semimajor axis and semiminor axis of the red projected ellipse IR and ST , respectively. x dr , y dr , z dr are the components of the derotated reference frame (dr ) as defined in the text. respectively. Now that we know all of the parameters necessary to describe the projected ellipse, we can standardize this ellipse into a simpler form to simplify subsequent calculations. We define a new frame (dr as in Fig. 4) as the derotation and geometric centering of the projected orbit such that the semimajor axis of the projected ellipse (O 0 I) is aligned with ̱x dr , the semiminor axis of the projected ellipse (O 0 S) is aligned with ̱ŷ dr , and O 0 is the origin of the dr coordinate system. The dr coordinates of the star location F (x Ã , y Ã ) are given by a simple rotation of the projection of O 0 F onto the ̱ x; ̱ y plane by angle θ: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 1 6 6 x Ã y Ã

Global and local extrema of planet-star separation
Before solving for the true anomalies where the orbit's projected separation s is equal to s WA (general working angle separation WAd i ), we first need to know how many of these s WA -orbit The projected ellipse (red) has semimajor axis IR and semiminor axis T S (purple). Both are calculable from the conjugate diameters that are the projections of the semimajor and semiminor axes of the original 3D ellipse (A 0 B 0 , C 0 D 0 , blue). The line QQ 0 is drawn such that it is perpendicular to the smaller conjugate diameter (C 0 D 0 ) and is bisected by B 0 : jB 0 Qj ¼ jB 0 Q 0 j ¼ jO 0 C 0 j. The semimajor axis of the projected ellipse is the angular bisector of O 0 Q and O 0 Q 0 . 15 Finally, the semiminor axis T S of the projected ellipse is perpendicular to the semimajor axis.
intersections we are looking for. We find the expected number of solutions by finding the s extrema throughout the orbit. We do so by solving for the roots of the derivative of the projected planet-star separation. We start with the general equation for an ellipse: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 6 ; 6 8 7 x e a p 2 þ y e b p 2 ¼ 1; (17) where x e and y e are the coordinates of any point on the ellipse and a p , b p are the semimajor and semiminor axes, respectively. We rewrite this in terms of x e , giving E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 1 1 6 ; 6 1 7 The projected separation is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 1 1 6 ; 5 5 4 Taking the derivative of Eq. (19) with respect to x e , substituting Eq. (18) (and its derivative), and setting it equal to zero, we have E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 6 ; 4 9 8 Isolating the square root term to one side and squaring both sides of the equation gives which is expanded with coefficients of x e , to get the polynomial expression E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 1 1 6 ; 3 6 7 0 This expression is a fourth-order polynomial, which we write in standard quartic form as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 1 1 6 ; 2 7 7 where A 0 , B 0 , C 0 , and D 0 are the constants and functions of a p , b p , x Ã , and y Ã , defined in Appendix F. We now apply the analytical solutions of the general quartic expression in standard form given in Appendix K. Solving this quartic gives us a set of four x e solutions (x e ) corresponding to two global extrema and two local extrema (if the latter exist for a particular orbit). Figure 6 shows a schematic representation of an orbit with four extrema. We use the imaginary components of the solutions, geometry of the ellipse, and magnitude of the higher order terms in the quartic solutions to identify which solutions belong to which extrema. We first use the magnitude of the imaginary components of solutions to determine how many extrema there are and filter out solutions that are not extrema. The algebraic solutions of the quartic polynomial all give values in the first quadrant (quadrants 1 to 4 are numbered counter-clockwise such that quadrant 1 has strictly positive coordinates, see the four quadrants of Fig. 6), so we must use geometry of the problem to determine the proper sign of each extremum's true coordinates x e;g and y e;g [g references solutions 0 to 3 Eq. (91)-Eq. (94)].
All numerical solutions to the quartic have some degree of imaginary component due to accumulation of numerical errors. Quartic solution sets with only two extrema will have two real solutions (solutions with small imaginary components only due to numerical error) and two solutions with large imaginary components (algebraic solutions that are artifacts to be thrown away). The majority of KOE have only two s extrema. The only case in which no extrema exist is in a circular face-on orbit. In this specific case, all solutions to the quartic will be nearly identical and have large imaginary components, and the resulting s extrema will be identical (s min ¼ s max ). We assume that all solutions are real if jIðx e;g Þj < 10 −5 ∀ g ∈ f0;1; 2;3g. We assume only two solutions are real if jIðx e Þj < 10 −5 for only two solutions. We define a new ordered set containing either two or four elements depending on magnitude of the imaginary components as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 4 ; 1 1 6 ; 2 7 2 x R ¼ fjRðx e;g Þj∶jIðx e;g Þj < 10 −5 ∀ g ∈ 0..3g; (24) where the absolute value is due to the algebraic solutions of the quartic being only defined in the first quadrant. We now know the number of expected solutions from the dimension of x R but need to leverage the quartic algebraic solution and geometry to determine which components belong to which quadrant. For orbits with x R containing only two solutions, the first two solutions to the quartic (x R;0 and x R;1 ) produce the largest and smallest x R;g − x Ã , respectively. Due to the shape of an ellipse, these must necessarily produce s min and s max . We define the set y R as the application of Eq. (18) to each element of x R . We define four separation quantities from the possible sign combinations of the coordinate magnitudes in x R and y R as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 1 1 6 ; 1 3 0 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 1 1 6 ; 7 9 s þAE0 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðx R;0 þ x Ã Þ 2 þ ðy R;0 AE y Ã Þ 2 q ; (26) Fig. 6 Planet orbit (black) in the dr frame and planet-star separation extrema. The minimum separation (cyan) always occurs in the same quadrant as the star in the dr frame (orange ×). The maximum separation (red) always occurs in the quadrants opposite the star. The local minimum separation (magenta) and local maximum separation (gold) occur in the same half-plane about the y axis as the star, but in opposite half-planes about the x axis as the star (fourth quadrant). This configuration applies to the vast majority of orbits. A small number of edge cases exist, including circular orbits and some edge on orbits. The blue dots represent the foci of the projected ellipse in the dr frame, the orange dot is the ellipse center, and the purple dashed lines are the semimajor and semiminor axes.
For the orbits with four solutions in x R , the first two elements will always be the global extrema, and the last two elements will be the local extrema, which are given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 4 ; 1 1 6 ; 3 0 4 where s −þ2 and s −þ3 are calculated in the same manner as s −þ0 and s −þ0 . We are able to find the coordinates of all four using the same logic as for finding s min , s max , s lmin , and s lmax . This procedure yields the coordinates of all existing extrema in the dr frame. To find their locations on the projection of orbit in the plane of the sky, we apply the inverse of Eq. (16).

Intersections between a circle and an ellipse
We find s WA -orbit intersections by formulating the circle-ellipse intersections as another quartic, solving this, and assigning the algebraic solutions to intersections. We assign solutions using the number of s extrema, the size of the s WA intersecting circle relative to these s extrema, and the ellipse geometry. To formulate the s WA -orbit intersections as a quartic, we start with Eq. (19) and substitute in Eq. (18). This gives us a separation equation solely as a function of x e and star location. We expand and transform this into a general polynomial of x e with a general s As in Sec. 2.1.2, we divide by the leading coefficient to convert this to the general quartic form: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 6 ; 1 1 6 ; 6 2 1 with A 1 , B 1 , C 1 , and D 1 given in Appendix G. We solve this using the general quartic solution as given in Appendix K. This results in a solution for the x e of intersections in the first quadrant of the dr frame. We always have four algebraic x e solutions that may or may not correspond to actual s WAorbit intersections. The number of s extrema (either two or four) and projected separation relative to these extrema determine how many intersections will occur. (If we expect two intersections, then two of the four algebraic solutions must be real solutions, and the other two are some combination of repeated roots or non-physical imaginary solutions.) In each of these cases, we must handle the assignment of x e and y e solutions to the correct quadrants. Figure 7 shows a schematic representation of two cases corresponding to four total intersections.
For KOE with four extrema and s min < s WA < s lmin , we know that there will be two intersections on the same y side of the ellipse as the star (quadrants 1 or 2 in the dr frame). Of the four quartic solutions to Eq. (36) that we have to choose from, we know that x 0 is one of them. The other solution could either be For KOE with four extrema and s lmax < s WA < s max , we know that there will be two intersections on the opposite x side of the ellipse as the star (quadrants 2 and 3 in the dr frame). In all cases, x 0 and x 1 are the intersection solutions. This is because the first two solutions have the largest term in the quartic solution [Eqs. (91) and (92)]. x 1 is slightly smaller than x 0 because it subtracts the second largest term. The relative magnitudes of x 0 and x 1 determine that y 0 occurs in the same side of the ellipse as the star and y 1 must occur on the opposite side of the ellipse as the star. Therefore, (x 0 , y 0 ) occurs in quadrant 2, and (x 1 , y 1 ) occurs in quadrant 3.
For KOE with four extrema and s lmin < s WA < s lmax , we know that there will be four intersections. Unlike in Sec. 2.1.2 and the rest of this paper where x 0 through x 3 are ordered as in Appendix K, we order x h based off Δx h ¼ jx h − x Ã j. We order the quartic solutions from x 0 to x 3 such that x 0 is where minðfΔx h ∀ h ∈ f0;1; 2;3ggÞ and x 3 is where maxðfΔx h ∀ h ∈ f0;1; 2;3ggÞ. These newly ordered Δx h correspond to those shown in Fig. 7. The (x 2 ; y 2 ) intersection occurs in either quadrant 1 or 2, but it is always above and to the left of the star in the dr frame and has the second largest Δx h component. The (x 3 , y 3 ) intersection occurs in either quadrant 1 or 4, but it always has the largest Δx h component. We resolve the sign of y 3 by testing it in both quadrants. In >99.992% of cases, we assign y 3 to quadrant 1 as in Fig. 7(a), but for a minority of cases, its correct assignment is quadrant 4 as in Fig. 7(b). The (x 1 ; y 1 ) intersection always occurs in quadrant 4. The (x 0 ; y 0 ) intersection occurs in either quadrant 3 or quadrant 4, but it is always below and to the left of the star in the dr frame and has the smallest Δx h component.
For KOE with two extrema and s min < s WA < s max , we know that there will be two circleellipse intersections. The KOE determine where the star is located in the dr frame. The location of the star in the dr frame relative to the vertices of the projected ellipse [ð0; b p Þ, ð0; −b p Þ, ð−a p ; 0Þ, and ða p ; 0Þ] determines the star-vertex separation ordering and subsequently to which quadrants the two intersection solutions belong. Instead of calculating the star-vertex separation for each orbit, we divide the first quadrant into four regions that specify the star location type [types 0 to 3 as indicated in Fig. 8(a)]. This means that any KOE with the star in location type 2 has the associated star-vertex separation ordering.
Using the equidistant lines between ellipse vertices, we divide the first quadrant into four regions (types 0 to 3) as shown in Fig. 8(a). Regions 0 and 2 are divided by the line defined by where a projected ellipse and star is of type 0 if y Ã > y a p −x;b p −y ðx Ã Þ and one of types 1, 2, or 3 otherwise depending on the other equidistant lines. Similarly, regions 1 and 2 are divided by the line defined by and the star-ellipse classification type is calculated similar to type 0. An example star is shown in Fig. 8(a). This star location is type 2 because it is in the region bounded by the three dashed lines (gray, pink, and teal). This type 2 star has the separation ordering indicated in Fig. 8(a) and Table 1. Let us consider a s WA separation circle such that s a p −x Ã ;y Ã < s WA < s x Ã ;b p þy Ã . We, therefore, know that the two intersections must occur in quadrant 2 and quadrant 4 of the dr frame.
The distances s x Ã ;b p þy Ã , s x Ã ;b p −y Ã , s a p þx Ã ;y Ã , and s a p −x Ã ;y Ã in Fig. 8(a) are the distances of the star to each of the ellipse vertices and are calculated by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 0 ; 1 1 6 ; 1 5 6 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 1 ; 1 1 6 ; 1 0 4  Table 1 define orbit types sorted by vertex distances from smallest to largest based on star location type in the first quadrant of the dr frame. Finally, after determining the correct number of intersection solutions and the proper quadrants to which these solutions belong, we rerotate and translate the intersection solution locations back into the 2D projection of the 3D orbit.

Δmag Intersections
In this section, we present our method for calculating the values of ν on a planet's orbit where the planet has a specific value of Δmag lim , called Δmag intersections. As in Sec. 2.1, to compute these solutions, we first need to calculate all Δmag extrema (Δmag min , Δmag max , Δmag lmin , and Δmag lmax ). The process for calculating Δmag extrema is detailed in Appendix C, but briefly included as follows: 1. express Δmag as a polynomial in cosðνÞ; 2. find the values of ν and Δmag for all polynomial roots; and 3. remove invalid and duplicate solutions.
The general process for calculating the true anomalies where a planet has Δmag lim is similar to the extrema-finding process but with some minor modifications. After calculating the Δmag extrema, we determine how many Δmag intersections a given orbit should have with a particular Δmag lim value. If Δmag lim < Δmag min or Δmag lim > Δmag max , there are no intersections. If Δmag min < Δmag lim < Δmag lmin or Δmag lmax < Δmag lim < Δmag max , then there are exactly two intersections. If Δmag lmin < Δmag lim < Δmag lmax , then there are exactly four intersections. When the orbit does not contain local Δmag extrema, there are just two Δmag intersections if Δmag min < Δmag lim < Δmag max .
Knowing the number of solutions to expect, we follow the same steps as above: finding a governing polynomial equation, solving for its roots (represented by the ordered set x), and filtering out the relevant solutions. We additionally throw out any solutions with large errors from the input Δmag lim . The full process outline is included in Appendix D and discussed in depth below.
We start with the definition of Δmag given in Eq. (1). Although this expression contains multiple terms that are fully defined by an orbit's KOE, it is also a function of other planet properties, including the planet's radius, geometric albedo, and phase function. To make the mathematical development presented below tractable, we assume that the quasi-Lambert phase function 17 is a sufficient approximation of any planet's phase function. In general, the quasi-Lambert phase function is a better representation of the Earth's phase function than the Lambert phase function.
This quasi-Lambert phase function is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 2 ; 1 1 6 ; 1 1 5 We take advantage of this function's form by substituting the half-angle formula:  . Note that this expression for β depends on making the approximation that the observer-star and observer-planet vectors are parallel, which introduces minor error given the large distances to even the nearest stars. 18 After making these substitutions, as well as taking Eq. (4) for the planet-star distance term j̱ r k∕i j, simplifying, and collecting all of the non-orbital planet parameters on the left side, we find E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 6 ; 1 1 6 ; 5 3 7 ðe cosðνÞ þ 1Þ 2 ð1 þ sinðiÞ cosðωÞ sinðνÞ þ sinðiÞ sinðωÞ cosðνÞÞ 2 : (46) We now have an expression that isolates terms containing ν and can be decomposed into a numerically solvable polynomial. To improve solving efficiency and determine which subset of planets should have ν solutions for a given Δmag, we need to first calculate the Δmag extrema over the full orbit. The outline of the process for calculating these Δmag extrema is included in Appendix C.
To calculate these Δmag extrema, we first find the derivative of Eq. (46) and multiply by two (for simplification purposes) to get to get an expression in cosðνÞ only. We define x ¼ cosðνÞ, isolate the ffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 − x 2 p term, square both sides, and expand to get an eighth degree polynomial in x of the form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 9 ; 1 1 6 ; 1 0 2 Although eighth degree polynomials do not have analytical solutions, a numerical root solver can determine the eight roots of this function denoted as the ordered set x. We discard any x p > 1 ∀ p ∈ f0: : : 8g, x p < −1 ∀ p ∈ f0: : : 8g, or solutions with large imaginary components. We then calculate the remaining true anomaly solutions by ν 0 ¼ cos −1 ðxÞ and ν 1 ¼ 2π − ν 0 . Of the remaining valid solutions, we identify whether each solution is an extremum by evaluating whether Δmagðν 0 AE δνÞ are both larger or both smaller than Δmagðν 0 Þ. If identified as a potential extremum, the smallest and largest extrema are assigned to Δmag min and Δmag max . The remaining extrema are checked for duplicates, which are identified by solutions with (ν, Δmag) values close to existing extrema. If any solutions are remaining, an additional check and assignment is then made for the local extrema. Through this process, a solution identified in ν 0 has the associated solution in ν 1 removed.
To calculate the Δmag intersections, we apply the same process for turning Eq. (46) into a polynomial as in the Δmag extrema calculation. The full outline for calculating the Δmag intersections is given in Appendix D. By following this process, we arrive at E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 0 ; 1 1 6 ; 5 6 6 (50) where ξ is the collection of constants on the left side of Eq. (46). The coefficients of this polynomial are given in Eq. (76) of Appendix I. We again use a numerical root solver to find the roots of this function x and again discard all x p > 1 ∀ p ∈ f0: : : 8g, x p < −1 ∀ p ∈ f0: : : 8g, and solutions with large imaginary components. We then calculate the true anomalies of all remaining solutions by ν 0 ¼ cos −1 ðxÞ and ν 1 ¼ 2π − ν 0 . We then evaluate Δmag 0 ¼ Δmagðν 0 Þ and Δmag 1 ¼ Δmagðν 1 Þ and remove solutions where jΔmag 0 − Δmag lim j > 0.01 and jΔmag 1 − Δmag lim j > 0.01. (Note that 0.01 is < AE 0.08% error on Δmag.) We do an iterative process of selecting (ν, Δmag), which are unique (not duplicate solutions) and are closest to the expected Δmag lim , until we have the expected number of solutions. In the majority of cases in which solutions exist, there are generally only two viable solutions, the assignment of which is simple. In some cases, the solution selection is more ambiguous as double roots are possible. 1 þ e cosðνÞ að1 − e 2 Þ ¼ 1 X ½cosðΩÞ cosðω þ νÞ − sinðΩÞ sinðω þ νÞ cosðiÞ:

ν from X and Y
In both of these expressions, we substitute angle addition formulas of sinðω þ νÞ and cosðω þ νÞ and subsequently set the two equations equal to each other. We set Eqs. (51) and (52) equal to each other, expand, isolate a cosðνÞ term and sinðνÞ term, solve for ν, and rearrange to get E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 3 ; 1 1 6 ; 1 5 3 We now have an analytical expression for ν solely as a function of the KOE, X, and Y of a particular point on the orbital ellipse.
By combining the original X and Y equations to solve for ν analytically, we have created two potential solutions at ν and ν AE π. One of these is the correct ν value and the other is not. To calculate the correct intersection point, we calculate the separations at both ν and ν AE π to find absolute error in s and use the smaller error of the two.

Calculate t from ν
We have ν for the locations where s extrema, Δmag extrema, s WA -orbit intersections, and Δmag intersections occur, but we need them in terms of time. We calculate eccentric anomaly E of these events directly from ν as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 4 ; 1 1 6 ; 6 1 8 sin ν e þ cos ν ; (54) which gives the corresponding time E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 5 ; 1 1 6 ; 5 6 0 where T is the orbital period of a general planet.

Converting from Star to Star
The planet-star intersection points and visibility ranges can be saved as either true anomalies or as specific times for a reference star. It makes the most sense to store the values as times and use 1M ⊙ as the reference time. The times are scaled to any star mass M by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 6 ; 1 1 6 ; 4 3 7

Calculating Completeness
We calculate integration time adjusted completeness using the aggregated fraction of time that planets are detectable. Using the methods described in Sec. 2 on the instrument's photometric visibility limit (Δmag lim ) and astrometric visibility limits (IWA and OWA), we calculate the specific times that the planet enters or exists the instruments visibility limits. By collecting these times and inspecting intermediate test points, we identify the time windows in which any planet is detectable by the instrument. Given an integration time (t max ) required to reach some SNR consistent with a clear detection (typically a value > 5), we throw out all visibility time windows less than t max and discount all other visibility windows by t max . Dividing the sum of all integration time discounted visibility windows for a planet by its orbital period gives us the fraction of time that a planet is detectable by the instrument. Aggregating the fractions of time that planets are detectable by the instrument gives us integration time adjusted completeness. In mathematical form, this is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 7 ; 1 1 6 ; 2 2 2 Here δt j;k is the j'th time window larger than t max in which the k'th planet is visible, U k is a boolean that indicates that the planet is always visible 0 or at least sometimes visible 1, N is the total number of planets, and T k is the orbital period for the k'th planet. We also calculate the completeness for an Earth-like exoplanet of the Earth-like exoplanets subpopulation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 8 ; 1 1 6 ; 1 1 7 where N È is the total number of Earth-like planets simulated. Similarly, this capability can be extended to the entirety of the Kopparapu classification scheme. 4 Our method uses orders of magnitude fewer exoplanets and an equivalent memory but fills in the Δmag space using an order of magnitude better accounting with strategically calculated true anomalies.

Calculating Dynamic Completeness
The methods presented in this paper can also be used to calculate dynamic completeness. 13 Dynamic completeness extends the original concept of completeness by considering the fraction of planets observable on subsequent observations of the same target. For the original formulation of completeness, which relies on simulating a cloud of planets, this requires propagating every simulated planet along its orbits. Dynamic completeness is frequently used to compute the fraction of planets in a population that are initially undetectable but become detectable upon a second observation some time later.
We define P as the unordered set of all planets, P detected;1 as the set of planets detected at observation time one, and P undetected;1 as the set of planets not detected at observation time one. To find the sets of detected planets at the first observation, we start by generating a random observation time for the k'th planet (t start;k ) between 0 and T k . We then determine which planets are within the instrument's visibility limits at t start;k to create P detected;1 and its complement P undetected;1 .
For a subsequent observation some time later (t wait ), we find the time past periastron at which the observation occurs by [modðt start þ t wait ; T k Þ] and determine which planets are within the instrument's visibility limits at that time. This defines the unordered set of planets detected at observation time two as P detected;2 . Its complement is P undetected;2 , the set of planets not detected at observation time two.
Dynamic completeness for the second observation of a target star is given by the fraction of planets undetected in the initial observation but detected on the second one: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 9 ; 1 1 6 ; 4 0 8 where nðPÞ is a function giving the number of elements in set P and N is the total number of planets. Dynamic completeness of the m'th visit is similarly given by the fraction of previously undetected planets that are detected on that visit: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 6 0 ; 1 1 6 ; 3 1 8 3 Results

ν from s
We apply the methods developed Sec. 2 to calculate the planet-star separation extrema, the true anomalies at which they occur, the times past periastron at which the extrema occur, the true anomalies where the separation circle intersects with the projected orbital ellipse, and the times at which these intersections occur for a planet orbit. Figure 9 shows the projected separation of a sample orbit with global and local extrema calculated via Sec. 2 methods. This figure also demonstrates the ability of our methods to find all true anomalies (and times) when the projected separation takes a specific value (in this case s WA ¼ 1 AU).
We now determine the accuracy of the methods implemented. To do this, we calculate the true anomalies of s WA -orbit intersections for 10 5 orbits using the method in Sec. 2.1.3 and s WA ¼ 1 AU: We then calculate the planet-star separation by substituting these true anomalies back into Eq. (2) and plot the error between these separations and the input s WA in Fig. 10. Of the 10 5 planets orbits simulated, ∼25; 000 orbits produced two or more intersections. The largest intersection error observed is <10 −6 . Machine precision error is ∼10 −16 , and the square root of this is 10 −8 , which indicates that we have approached machine precision in our method.

ν from Δmag
We apply the methods developed Sec. 2 to calculate the Δmag extrema, the true anomalies at which they occur, the times past periastron at which the extrema occur, the true anomalies where the Δmag intersections occur, and the times at which these Δmag intersections occur for a planet orbit. Figure 11 shows the Δmag of a sample orbit with global and local extrema calculated via Sec. 2 methods. This figure also demonstrates the ability to find all true anomalies (and times) when the Δmag takes a specific value (in this case Δmag lim ¼ 25.0).
To check the error in true anomaly produced by the Δmag intersection method in this paper, we compare the ν of Δmag intersections calculated with alternative numerical solving methods. The method presented in this paper finds Δmag intersections of 10 6 planets within ∼6.3 s. The first error checking method uses a cubic spline fit to 300 (ν, Δmag) points along each planet's orbit. We then subtract Δmag lim ¼ 29 (our test point) from the spline and find the roots. This univariate spline root solving method is capable of executing in 419 s on the 510,120 planets in the population that produced two intersections, a rate of 8.1 × 10 −4 s per planet. The univariate spline root solving method is 100× slower than the Δmag intersection method in Sec. 2.2. The separation error between the univariate spline method and the Δmag intersection method is <6 × 10 −5 for >99.98% of planets and <10 −4 for the other 0.02%. These high-error targets are planets with low variation in ΔmagðνÞ orbits. The second-error checking method solves for the (ν, Δmag) intersections using a numerical minimization method on a random subset of 10 4 planets. This method is far more inefficient, taking 1646 s on 10 4 planets with two intersections at a rate of 0.1646 s per planet. Using a minimization function allows us to determine the error in Δmag of an intersection point to within 10 −8 . This numerical minimization method independently confirms the accuracy of solutions to the univariate spline roots method on the limited number of targets tested. We plotted the normalized frequency of true anomaly error of both the numerical minimization method and univariate spline root solving method compared against the Δmag intersection method in Fig. 12. Note that the total number of incidences of a given error is normalized by the bin width and total number of targets, so the non-normalized frequency of the 3 × 10 −11 bin is incredibly large compared with the 5 × 10 −4 bin. Both methods indicate that the resulting true anomalies of the Δmag intersections are within 10 −4 rad of each other.
Although the Δmag intersection method featured in this paper is orders of magnitude faster than the univariate-spline-roots method and just as accurate, the quasi-Lambert phase function required for its use is not the best phase function for all planets; however, it does fit Earth-like planet phase function quite well. At a substantial time cost, the univariate-spline-roots method can be used for planets with any phase function.

Convergence and Validation
We test convergence of our method for calculating completeness by repeatedly calculating completeness on 10 5 planets. We test convergence of Brown's method by repeatedly calculating completeness over a set of a logarithmically increasing number of planets. We use the same SAG13 planet population parameters from Ref. 1 but make the substitution of a quasi-Lambert phase function to define the underlying planet population. Since this technique will be most relevant to future highly capable telescopes, we are using the HabEx IWA of 0.045 arc sec, OWA of 6 arc sec, and upper Δmag limit of 25 on a star that is 10 pc away. 8 This resulting instrument visibility limits and planet population cover a well populated region of the Δmag versus s joint probability density function [see Fig. 1 in Ref. 1 for a comparable example of the joint probability density function]. Convergence of the two completeness methods is demonstrated in Fig. 13.
Both methods presented in Fig. 13 demonstrate similar percent convergence to their own converged means. Completeness repeatedly calculated using the method in Sec. 2.6 with 10 5 planets has a mean of 0.25785 with a standard deviation of 0.0010 for 4.28 × 10 8 planets. In comparison, the Brown completeness method converges to 0.25783. This error in the converged mean of the two methods is consistent with the completeness error identified in Fig. 14. We expect that completeness calculated using the method in this paper will have comparatively better convergence in sparsely populated regions of the Δmag versus s joint probability density function (see Fig. 1

in Ref. 1).
Since completeness for an individual planet is the fraction of time that a planet is visible to the instrument, we validate individual planet visibility windows using test points. To calculate a ground truth fraction of time that a planet is visible, we create 10 5 test points, evenly spaced in time, and calculate both the planet-star separation and planet-star Δmag at each of theose test points. We then determine whether each point is within the visibility limits of the instrument. The fraction of visible points is the fraction of time that a planet is visible. We repeated this for 25,000 planets due to computational cost of the test point method. The histogram of error in completeness for individual planets calculated using the test point method and method in this paper is shown in Fig. 14. The maximum error in completeness is observed to be near 5 × 10 −5 , which is consistent with the converged error between Brown completeness and integration time adjusted completeness.
The test point method that we used to validate integration time adjusted completeness calculations took over 52 h on only 25,000 planets (less than a quarter of the planets needed to compute completeness using integration time adjusted completeness). The computation time for the integration time adjusted completeness calculation on 10 5 planets is 21.34 s with a standard deviation of 0.4 s over 1000 calculations. Brown completeness, including the joint probability density function generation, has a mean execution time of 3.19 s and standard deviation of 0.13 s over 1000 calculations. This execution time of all of these methods scales linearly with the number of planets. Some room for optimization exists in integration time adjusted completeness.

Completeness Versus Integration Time
To evaluate the effects of integration time on completeness, we calculate integration time adjusted completeness for various integration times, star distances, and planet populations (SAG13 and Earth-like planet population). Figure 15 shows the integration time adjusted completeness for the SAG13 planet population 1 and the Earth-like planet subpopulation defined in Appendix L. The decrease in completeness for longer integration times is most prominent for nearby stars for both populations. As expected, longer integration times decrease completeness.
For the assumed observatory parameters, the Brown completeness of an Earth-like population for a target 5 pc away is 0.583. As we increase integration time to 1, 2, and 5 days, the associated integration time adjusted completeness decreases by 0.63%, 1.27%, and 3.15%, respectively. Recalculating Brown completeness for a target at 25 pc in Fig. 15(b) gives 0.583. As we increase integration time to 1, 2, and 5 days on this 25 pc target, the associated integration time adjusted completeness decreases by 0.989%, 1.97%, and 4.92%, respectively. Our integration time adjustment of completeness shows Brown completeness overestimates exoplanet yields for any observation. Brown completeness applied to stars farther away results in a substantial overestimation of exoplanet yield.
Technically, the Δmag lim used as the upper limit for calculating completeness is a function of integration time and will approach a theoretical upper limit as the noise floor is reached. By calculating Δmag lim using Eq. (12) of Ref. 1, we can calculate completeness at each integration time to get the completeness versus integration time curve. We created a 4-m telescope as in Appendix M and evaluated completeness at each integration time to create an example demonstrating how Brown completeness and integration time adjusted completeness vary with integration time in Fig. 16. We can see that the completeness of both methods tracks very closely until they diverge. Brown completeness approaches its asymptotic limit as expected, whereas integration time adjusted completeness approaches a maximum and begins decreasing past ∼0.2 day.

Dynamic Completeness
We compute dynamic completeness for the example test case in Fig. 1 of Ref. 13 using the method described in this paper and replicating the approach of the original work. This computation is only for a second epoch. We replicated the original work by sampling a large number of planets, finding the planets initially visible, propagating these planets to some time past initial observation, and determining what fraction of these planets were not detected the first time but were detected the second time. We replicated Brown's work using both the Lambert phase function (orange line in Fig. 17) to match the original work as well as the quasi-Lambert phase function (red line in Fig. 17).
We compared the dynamic completeness computation time for 1000 dynamic completeness calculations testing 1000 individual points in time past the initial observation. Brown's dynamic completeness had an average execution time of 34.62 s with a standard deviation of 0.99 s. Calculating dynamic completeness using the method presented in this work resulted in an average execution time of 4.757 s with a standard deviation of 0.036 s. The dynamic completeness calculation presented in this paper is ≈7× faster than the traditional Brown completeness method, enabling its use in optimization.

Convergence
Both methods of calculating completeness suffer from sparse sampling of the exoplanet population parameter space. Brown completeness is most affected as it needs planets sampled over Ω, ω, a, e, i, ν, p, and R. The current implementation of our method needs the same parameters except for ν. We are therefore surprised that Brown's method and our method share such similar completeness convergence. We hypothesize that the similar convergence results are because the integration bounds used for calculating completeness are over a dense region of the joint probability density function. If we examined a more sparse region of the joint probability density function such as high Δmag or larger separation (where the larger, lower occurrence rate planets occur), we expect our completeness methods to have marginally better convergence for low numbers of planets. Regardless, integration time adjusted completeness easily reaches completeness errors below 10 −3 (translates to ≈0.8% where C ¼ 0.12, the largest completeness of an observation in Ref. 1), which is a sufficient error for use in optimization.

Reducing Parameter Spaces
The integration time adjusted completeness method implemented in EXOSIMS is not cached for computational efficiency like the Brown completeness is. Although we are able to store the joint probability density function of Δmag versus s with Brown's method, integration time completeness requires the KOE of each simulated planet to be stored. The curse of dimensionality prohibits us from finely sampling the entire subspace of exoplanets, making the caching of integration time adjusted completeness prohibitive, but not all our methods require all 6 KOE, p, and R. So, in the future, it may be possible to create a representative subsampling of planets weighted by occurrence rate based on the parameters needed for different calculations. For example, s extrema, s WA -orbit intersections, and Δmag extrema calculations only require ω, a, e, and i (Note that the ν locations of the Δmag extrema are independent of p and R, but the magnitude of Δmag extrema are dependent upon them). The photometric property p × R could then be sampled and independently combined with individual (ω, a, e, i) combinations. The primary benefit is that, given some instrument parameters and star properties, we could determine the subset of parameters that can physically be observed prior to calculating completeness, thus making the per-star completeness calculations more efficient.

Limitations
The approach we implemented in this paper did not vary the Δmag lim as we varied the integration time. If properly executed, Δmag lim should increase with increased integration time. However, calculating the planet true anomaly intersections with a given Δmag lim is the most expensive calculation, so repeatedly calculating this to convergence is not desirable.
The Δmag lim of coronagraph designs are working angle dependent and therefore separation dependent, but we assume that this limit is a fixed quantity. If there is substantial variation in the Keithly, Savransky, and Spohn: Integration time adjusted completeness limiting Δmag, it is better to use the conservative value with this method. If a higher level of predictive accuracy is desired, these methods could be used to ascertain if any intersection with the Δmag lim exists and subsequently another method could be used on the subset of planets with known intersections.

Impact of the Integration Time Adjusted Completeness
As we showed in this paper, integration time adjusted completeness and Brown completeness converge to the same value to within 0.00002 when evaluated at t max ¼ 0 day. The integration time adjusted completeness for an observation of an Earth-like planet population on a 1-M ⊙ star 5 pc away and integration time of ≈1 day is 0.64% lower than the comparable Brown completeness. For reference, the Roman target list in Table 9 of Ref. 1 has a maximal integration time of 1.71 days but most are <0.6 day. This means that we could expect an average reduction in yield below 0.64% when observing a population of Earth-like planets around Sun-like stars. This integration time completeness adjustment is within the 3.19%, 3σ, margin of error from 1000 Monte Carlo simulations of the cycle 6 Roman in Ref. 1 and cannot be considered statistically significant. However, a ∼4-day integration time on a Sun-like star 25 pc away observing an Earth-like population has a decrease in completeness above the Ref. 1 3σ margin of error. However, the threshold of statistical significance should not deter the widespread use of integration time adjusted completeness as the adoption of this method can shore up the differences between completeness-based yields and simulation-based yields. We can also say that integration time has a muted effect when observing stars that are farther away. This can most likely be attributed to the increase in s min and resulting decrease in the total time-fraction that smaller period planets spend within the visible limits of the telescope. Big planets with the orbital radius of Jupiter will move more slowly and be less affected by integration times.
In Appendix M, we optimized an exoplanet direct imaging mission maximizing single-visit Brown completeness yield and included the resulting Design Reference Mission (DRM) in Table 2. The resulting DRM has a Brown completeness yield of 387.16 exoplanet detections in a single-visit detection survey on average, but the integration time adjusted completeness yield expects only 354.67 exoplanets on average. This means that Brown completeness for the particular instrument parameters optimized over the SAG13 planet population with the particular mission parameters overpredicts the actual exoplanet yield by 32.49 exoplanets on average, 9.16% more than the integration time adjusted completeness yield. When calculating completeness, we are careful to scale the completeness limits of integration by the star's luminosity. For the integration time adjusted completeness method, we additionally scale the planet's orbital periods based off the mass of the host star. The percentage difference above and beyond that expected from Fig. 15 can be found by looking at these two adjustments applied to each star and the target list in Table 2. The average star in the target list has a larger mass and a brighter luminosity than that of the Sun. The brighter luminosity results in a smaller s min and s max , which serves to incorporate more smaller semimajor axis planets into the completeness calculations, and the larger star masses result in shorter periods, meaning that the planet visibility windows all decrease in duration.
Integration time adjusted completeness is crucial for determining the ability of an Earth-like planet to be spectrally characterized. A spectral characterization with a coronagraph could take between a few days and 60 days. Because of the long integration time, the planners of a HabEx use the "characterization completeness" of 10% and maximum integration time of 60 days as a filter on stars to consider observing. 8 The calculation of this critical filter could be substantially improved by considering integration time adjusted completeness if the spectrum of the planet must be taken all at once (i.e., observation spanning multiple weeks) and not spread across multiple epochs. Just because a planet can be detected does not mean that it can be spectrally characterized.

Dynamic Completeness and Computation Cost
The greatest benefit of using the methods in this paper to calculate completeness is the marginal additional cost of calculating dynamic completeness. Generally, dynamic completeness requires the computation of true anomalies from time, which is prohibitive for >10 5 orbits at >100 different times in the future. With our methods, calculating visibility windows allows us to use boolean operations to compute dynamic completeness about 7× faster than Brown's dynamic completeness method. The computational cost of calculating completeness using the method described in this paper is an order of magnitude larger than Ref. 2, but it can use orders of magnitude fewer planets to do so. Unlike Brown's method, within the computation time of our method, we also get additional desirable information about the detectable planets such as their s extrema and Δmag extrema.

Revisiting the Same Exoplanet
Another limitation present in the planning of exoplanet direct-imaging missions is the telescope keep-out angles. Figure 3 of Ref. 1 contains a keep-out map for a subset of the DRM created in that paper. The smallest percent of time that a target star is visible for the Roman is nominally 28%. Due to symmetry of the keep-out region, this translates into two separate time windows of visibility of ∼51 days. Instead of considering the integration time an integration time input, it can also be considered a revisit time input for determining the probability of being able to observe a planet twice in the same target-star visibility window.

Exoplanet Classification
The underlying methods in this paper are used to find locations along a planet's orbit where s and Δmag intersections occur. The methods in this paper, with some modification, can also be used to probabilistically classify an exoplanet subtype. 4 If an exoplanet is detected with a particular (s, Δmag) and an uncertainty region of s AE σ s and Δmag AE σ Δmag , then the methods in this paper can be applied to each of these four bounding lines. By finding the average orbital time-fraction that exoplanets of each type spend in the bounding uncertainty box, we can find the probability that the exoplanet belongs to a specific exoplanet subtype. This requires additional work beyond the scope of this paper.

Conclusion
We have demonstrated an accurate method for calculating integration time adjusted completeness and its adaptation to calculating dynamic completeness. In the process, we also created fast and accurate methods for calculating the true anomalies where a planet's orbit has specific values of projected separation Δmag and their extrema. We demonstrated how to use these methods to calculate a more accurate integration time adjusted completeness using the fractions of time that a planet is detectable by an instrument. We demonstrate that traditional methods of calculating completeness overestimate the number of observable planets because they do not subtract the integration time used in observing the target. For a Sun-like star at 25 pc with 1-day and 5-day integration times, integration time adjusted completeness of Earth-like planets is reduced by 1% and 5%, respectively. We applied integration time adjusted completeness to a target list optimized using the Brown completeness method and found that Brown completeness overestimated yields by 9.61%. We also demonstrated that our methods can be used to quickly calculate dynamic completeness for determining the optimal time to revisit a target star.

Appendix A: Common Notation
There are many common variable forms used in this paper that are simply summarizable. Any variable x refers to a 3D vector in X, Y, and Z coordinates of something. We use x to indicate an array of variables, specifically used when referring to multiple roots of a polynomial. jABj refers to the length of line segment AB. Line segment AB is treated as a vector. C 0 is a projection of point C.
There are multiple different subscripts with different meanings used in this paper. A subscript with x i refers to the index of a host star (out of some whole, non-descript, star catalog).
A subscript x k refers to an individual planet out of the whole large set of planets. Subscripts of x min , x max , x lmin , and x lmax are descriptors on the individual variable x that indicate that the variable is associated with the minimum, maximum, local minimum, or local maximum, respectively. When solving for the distance between points inside an ellipse and the vertices or covertices of that ellipse with semimajor axis a, we use the shorthand notation of the form s aþx;y . [This is specifically the distance between point (x, y) and (−a, 0). ] We denote coefficients of the polynomials of the four methods that solve the quartic in this paper as A # .
In Fig. 4, we reference many points on the 3D elliptical orbit and 2D projection of this orbit into the plane of the sky. Points on the original 3D ellipse are labeled B, H, P, C, K, O, D, A, B 0 , H 0 , P 0 , C 0 , K 0 , O 0 , D 0 , and A 0 .
In Appendix K, we use p 0 to p 11 as intermediate constants for simplifying the full quartic expression. P, D, and Δ (by itself and only in this section) are intermediate constants derived from quartic coefficients for determining the sign of the quartic solutions.   This results in a total of 84,210 planet-star intersections. After calculating and identifying true anomalies of intersections using the methods described in this paper, we evaluated the planet-star separation of each orbit at the true anomalies and calculated error from the input s WA .
(a) (b) Fig. 11 (a) The Δmag of a planet (black line) plotted against ν and (b) time past periastron. Δmag extrema are indicated by diamonds, where the maximum is red, local maximum is yellow, local minimum is magenta, and minimum is teal. Separation minimum and maximum are indicated by horizontal lines. For an input Δmag lim ¼ 25 (green line), we calculated the specific true anomalies (green dots) where the planet's Δmag intersects this line. The specific planet KOE are Fig. 13 Convergence of % error in completeness for varying numbers of planets used in its computation for Brown's method (green) and the method presented in this work with t max ¼ 0 day (purple). Completeness at the maximum number of planets is assumed to be the converged value of completeness of the respective methods. The converged mean of Brown's method is 0.25783 compared with the converged mean of the method in this paper of 0.25785. The standard deviation of the method in this paper is 0.0010. Fig. 12 A histogram of νðΔmagÞ error between the root solving method described in this paper and a root solver of a cubic spline fit method (purple) for 510, 120 planet orbits that produce two intersections (totaling 1,020,240 datapoints). We also compare the error between the root solving method described in this paper and an optimization method (blue) for 10 4 planet subset of the planet orbits that produce intersections. The optimization method investigates fewer planets because it is orders of magnitude more computationally expensive, but it is a fundamentally different approach to validating our method than a root solver. Keithly, Savransky,

Appendix B: True Anomalies of s Intersection Process
The process outlined here is used in Sec.   i. For all solutions (all real and only 2 real): A. Calculate s min from x 1 and assign ðx min ; y min Þ of quadrant 1. B. Calculate s max from x 0 and assign ðx max ; y max Þ of quadrant 1.
iii. Assign solution signs to proper quadrants. 9. Find circle and projected derotated centered ellipse intersection points. iv. Inside outer separation bounds, s min < s < s max . v. Inside local min/max separation bounds, s lmin < s < s lmax . vi. Outside outer separation bounds (no intersections).
(t) Calculate intersection solutions for planets with s min , s max , s lmin , and s lmax .
• Sort Δx k from min to max and rearrange the associated x k to match that order.
ix. Two intersections on opposite x side of ellipse as the star, s lmax < s < s max .
• y 1 is opposite sign of jyj.
(u) Calculate intersection solutions for planets with only s min and s max .
x. Calculate projected ellipse quadrant separation bounds s x;bþy , s x;b−y , s aþx;y , and s a−x;y . xi. Identify star location type and quadrant order.
• Type 0 occurs where s xþa;y < s x;yþb . Smallest to largest order: s x;b−y , s a−x;y , s xþa;y , s x;yþb .
• Type 1 occurs where s x;yþb < s a−x;y . Smallest to largest order: s x;b−y , s x;yþb , s a−x;y , s xþa;y .
• Type 2 occurs where s a−x;y < s x;yþb and s x;yþb < s xþa;y and s x;b−y < s a−x;y . Smallest to largest order: s x;b−y , s a−x;y , s x;yþb , s xþa;y .
• Type 3 occurs where s a−x;y < s x;b−y . Smallest to largest order: s a−x;y , s x;b−y , s x;yþb , s xþa;y . 8 Appendix C: ν from Δmag Extrema process Here we present the full outline of our process for calculating true anomalies where the planet's orbit has a Δmag extrema. These steps are discussed in Sec. 2.2.

1.
Substitute into Eq. (1) components and expand until all terms are functions of KOE.

2.
Manipulate the Δmag equation until all ν terms are isolated on the right side.

6.
Isolate all square root terms to general form Square both sides. 8.
Find polynomial coefficients of x. 9.
Remove solutions that are not extrema.
Remove duplicate solutions.

Appendix D: ν from Δmag lim Intersection Process
Here we present the full outline of our process for calculating true anomalies where the planet's orbit has a specified Δmag, which we call an intersection. These steps are discussed in Sec. 2.2.

3.
Substitute into Δmag equation components and expand until all terms are functions of KOE.

4.
Manipulate the equation until all ν terms are isolated on the right side.
Replace all cosðνÞ with x.

8.
Isolate all square root terms to general form A 2 ffiffiffiffiffiffiffiffiffiffiffiffi ffi Square both sides. 10. Expand.

11.
Find polynomial coefficients of x.
Verify that there are only two viable solutions. 20.
Assign solutions to ν.
10 Appendix E: Projected Ellipse Semimajor Axis and Semiminor Axis, and θ We write the analytical expressions for a p , b p , and θ using full expansions of Eqs. (13)- (15). These expressions are far too long to be practically conveyed, so we use the following intermediary parameters here and only here to simplify these expressions: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 6 1 ; 1 1 6 ; 5 9 4 K 0 ¼ eð1 − e 2 Þ K 1 ¼ sinðΩÞ cosðωÞ þ sinðωÞ cosðΩÞ cosðiÞ K 4 ¼ − sinðΩÞ sinðωÞ cosðiÞ þ cosðΩÞ cosðωÞ The simplified expression for the semimajor axis of the projected ellipse is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 6 2 ; 1 1 6 ; 3 4 2 The simplified expression for the projected ellipse semiminor axis is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 6 3 ; 1 1 6 ; 2 4 3 The simplified expression for the angle between the projected ellipse semimajor axis and X axis is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 6 4 ; 1 1 6 ; 1 3 3 Let us suppose that we have some ellipse on an arbitrary plane indicated by the black ellipse in Fig. 4. We will say AB and CD are the principal axes of some this 3D ellipse. AB and CD intersect at point O, the geometric center of this 3D ellipse. We will say point P is any point on the ellipse with H being the projection of point P onto axis AB and K being the projection of point P onto axis CD. By the definition of an ellipse, we have E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 7 7 ; 1 1 6 ; 6 5 3

Appendix K: Quartic Solution
The quartic function has a known analytical solution. 19 We start with a set of useful simplifying terms E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 7 9 ; 1 1 6 ; 4 1 6 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 0 ; 1 1 6 ; 3 5 4 (80) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 1 ; 1 1 6 ; 3 1 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 2 ; 1 1 6 ; 2 7 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 3 ; 1 1 6 ; 2 3 6 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 4 ; 1 1 6 ; 1 9 5 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 5 ; 1 1 6 ; 1 6 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 6 ; 1 1 6 ; 1 0 3 Theoretically, we always have four solutions. Theoretically, we can use equations of coefficients to determine how many roots are real, complex, and double roots. The general form of the quartic defines the expressions Δ, P, D 2 , R, and Δ 0 , which are E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 9 5 ; 1 1 6 ; 1 9 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 9 6 ; 1 1 6 ; 1 1 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 9 7 ; 1 1 6 ; 8 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 9 9 ; 1 1 6 ; 7 0 1 After evaluation of these constants, we can determine E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 1 0 0 ; 1 1 6 ; 6 7 7 solution types ¼ However, the use of this theoretical classification does not work when numerical errors are introduced. The accumulation of numerical errors causes the solutions to be improperly classified for a sufficiently large number of cases to prohibit the use of these classifications. The Δ expression rarely evaluates to Δ ¼ 0 when it should do so quite frequently. Numerical rounding errors can also frequently result in the evaluation Δ < 0 when in reality Δ > 0 and vice versa.

Appendix L: Earth-Like Subpopulation
We use a definition of Earth-like exoplanet similar to that used in Ref. 8. We use a more conservative planetary radius range of 0.9R È ≤ R ≤ 1.4R È . We use the flux-at-planet range of 0.3586 ≤ L plan ≤ 1.1080. To classify a planet, we use the time-averaged incident flux on the planet calculated by The resulting Brown completeness optimized integration times for a DRM are included in Table 2. The summed Brown completeness for this DRM is 68.89, whereas the same summed integration time adjusted completeness is 63.11, a difference of 5.78. We can multiply either of these numbers by the SAG13 exoplanet occurrence rate of 5.62 to get the expected number of exoplanets detected in a single-visit blind search. These are 387. 16 and 354.67 on average, respectively. This means that 32.49 planets are lost in a mission by not taking into account integration times when optimizing the mission schedule.              Laboratory. He received his PhD from Princeton University in 2011 and was a postdoctoral researcher at Lawrence Livermore National Laboratory, where he worked on integration, testing and commissioning of the Gemini Planet Imager. His research focuses on the optimization of ground and space-based exoplanet imaging surveys, control of autonomous optical systems, and advanced image processing.
Corey Spohn received a Bachelor of Science in both physics and engineering science and mechanics from Virginia Tech in May 2018. During his time at Virginia Tech he focused heavily on astrophysics while conducting research on anti-frosting surfaces, inter-droplet ice bridging, and flying snakes. For his senior design project, he investigated human-robot dynamics in a crowded panic situation by having humans play the children's game sharks and minnows with robots in a motion capture studio. He began working towards a PhD in the fall of 2018 and joined the Space Imaging and Optical Systems Lab in December 2018 where he is currently researching exoplanet mission scheduling.