Boundary diffraction wave integrals for diffraction modeling of external occulters

An occulter is a large diffracting screen which may be flown in conjunction with a telescope to image extrasolar planets. The edge is shaped to minimize the diffracted light in a region beyond the occulter, and a telescope may be placed in this dark shadow to view an extrasolar system with the starlight removed. Errors in position, orientation, and shape of the occulter will diffract additional light into this region, and a challenge of modeling an occulter system is to accurately and quickly model these effects. We present a fast method for the calculation of electric fields following an occulter, based on the concept of the boundary diffraction wave: the 2D structure of the occulter is reduced to a 1D edge integral which directly incorporates the occulter shape, and which can be easily adjusted to include changes in occulter position and shape, as well as the effects of sources---such as exoplanets---which arrive off-axis to the occulter. The structure of a typical implementation of the algorithm is included.


Introduction
One of the major goals of the coming decades is to directly image a terrestrial exoplanet. Direct imaging allows spectral information about the planet's atmosphere to be extracted from the spectrum of reflected light, including the potential presence of biomarkers, gases whose presence and abundance may indicate the presence of life. Detecting planets directly is a difficult task, however, for two reasons: intensity ratio (contrast) and angular resolution. An Earth-twin which orbits a Sun-twin at 1AU emits 10 10 times less flux than the parent star in the visible band; if the system was located at 10 parsecs from Earth, the angular separation of the two objects would be 100mas, which for most proposed space telescopes puts the separation under 4λ /D. Specialized methods are required to remove this flux at these small working angles.
One proposed method of suppressing the starlight is an occulter, a spacecraft with a shaped edge flown in front of the telescope to block the starlight before it arrives at the telescope. The occulter size (tens of meters) and distance (tens of thousands of kilometers) are chosen so the angular extent of the occulter is smaller than some desired angle, on the order of 100 milliarcseconds, so exoplanets in orbit about the star will still be visible. The edge is shaped to control diffraction, with the form chosen to suppress the light to a factor of 10 10 across the telescope aperture and over a wide spectral passband.
A major advantage of this approach, compared to an internal coronagraph, is that it significantly loosens the tolerances on the optical quality and thermal stability of the telescope optics, as well as eliminating the need for a wavefront control system for active correction. Conversely, the occulter itself must be built, deployed, and flown in formation to its own set of tolerances, and demonstrating this capability is an area of active research [1]. These tolerances include limits on permanent errors in the occulter shape due to manufacturing and deployment errors; transient shape deformations from thermal fluctuations and vibration; and position and orientation changes of the occulter relative to the telescope-target axis from formation flying [2,3].
The scale of an occulter-and the distance it must maintain with respect to its targetprecludes optical testing on the ground. All of these effects must be modeled in order to build an error budget and verify that an occulter can satisfy these requirements. The capture of the large dynamic range in the resulting field-10 5 or greater in amplitude-drives the accuracy of the propagator used for the models, while the modeling of time-varying thermal shape deformations and closed-loop formation flying in broadband light drive the speed of the calculation. Some of the most in-depth simulated observations can take from hours to days to run, and the ability to quickly and accurately model propagation past an occulter is thus a major enabling factor in the characterization of occulter system performance. Section 2 gives an overview of the current techniques; the new technique presented in this work is described in Sections 3 and 3.1, and a suggested implementation scheme, along with a pseudocode overview of the algorithm and some computational results, are given in Section 4.

Current occulter modeling techniques
Occulter designs take the general form of a solid central disk surrounded by N p identical tapering structures (petals) around the edge, giving the whole structure the general appearance of a flower. (An example is shown in Fig. 1.) This shape is a result of their provenance as (0,1)valued approximations to apodizers; the general method to design an occulter is to determine a smooth apodization profile A(r) which can provide the necessary starlight suppression, and then convert it into a binary shape with a sufficient number of petals so that performance is not undermined. The actual form of A(r) can be determined by optimization [4,5] or selecting an existing functional form [6,7]. For an occulter designed in this manner, the region Ω in which the occulter is opaque can be written directly: with (r, θ ) being polar coordinates in the plane of the occulter. Each of the [. . . , . . .] is a set of points which defines the outline of a single petal; denotes the union of these sets, which sets the boundary for the entire occulter. (We denote this boundary as ∂ Ω, as it will be used later.) For a typical occulter, the occulter-telescope distance z = 5 × 10 7 m; the wavelength λ = 5 × 10 −7 m; the occulter radius R = 25m; and the maximum excursions in the plane of the telescope |ξ | and |η| ≤ 3m. These values place the occulter well within the Fresnel regime; that is, the region in which it is valid to approximate the exponent in the propagation integral with the first two terms of a power series (the Fresnel approximation). This approximation is suitable when [8].
Further, a standard propagation integral computes the electric field past a finite aperture; the occulter, on the other hand, is a finite obstruction with an effectively infinite aperture. Here, we can turn to Babinet's theorem [9], which notes that field propagated through an aperture and the field propagated past its complement sum to the field that would have resulted were no obstruction present at all.
If a plane wave with amplitude A and wavelength λ is normally-incident on the occulter, as would be the case for an occulter properly aligned with a target star, the electric field at a distance z past the occulter would be written with the Fresnel approximation and Babinet's principle as where ξ and η are Cartesian coordinates at the downstream location. The "1−" term before the integral results from the use of Babinet's principle. Errors in (for example) occulter shape will modify Ω, but Eq. (4) will hold regardless. Computing the integral for points in the plane of a telescope aperture is the key step in modeling the performance of an occulter. One could attempt this integral by brute force: discretizing the shape of the occulter over a sufficiently fine grid to capture all of the necessary structure, and propagating this to the telescope aperture with direct integration or a matrix Fourier transform [10]. This approach is not an easy task, as the occulter will be tens of meters across, with tolerances at the level of tens of microns [2], requiring a grid up to 10 6 × 10 6 or more to capture the full shape of the profile, not including array padding.
A better approach is to take advantage of the structure of the occulter: in the absence of errors, it has an N p -fold symmetry in the azimuthal direction, and most errors introduce small perturbations to that. Taking advantage of the symmetry leads to three main approaches to modeling propagation past an external occulter: 1. Integrals in r: The Bessel function expansion [4] does the integral over θ explicitly, reducing the 2D integral to a series of single integrals in r whose contributions fall off exponentially fast near the optical axis, which makes computation extremely quick in this region. The series can be modified to incorporate some shape errors [11] at the expense of slower convergence. Unfortunately, it is slower at locations far from the optical axis, and requires a good deal of effort to incorporate many types of errors, such as occulter tilt, into the model. Aside from the propagation modeling, the first, radially-independent term in this series is used in optimization approaches [4,5] to design occulter shapes.
2. Integrals in θ : The Dubra-Ferrari integral [12] starts one step back from the Fresnel form of Eq. (4), at the Raleigh-Sommerfeld diffraction integral, and evaluates the integral over r to produce a pair of single integrals in angular coordinates, one of which holds for points inside the geometric extent of the starshade, and one which holds outside. This approach can incorporate most errors straightforwardly and spends the same amount of time to determine the field any point in the telescope aperture plane. Because of the manner in which the integrals are segmented, however, the integrands are more difficult to determine correctly for points physically outside the extent of the starshade.
3. Perturbations to a nominal shape: The slit approach [13] assumes the electric field for an occulter with no errors has been calculated by one of the above methods, and approximates the difference between the shape with and without errors as a series of long, thin boxes whose Fresnel transforms can be calculated quickly and exactly. This method is very fast but necessarily approximate.
In this paper, we propose a fourth method-the boundary-diffraction-wave approach-which seems to out-perform the three extant methods for general purpose use. It maintains the constant-time property of the Dubra-Ferrari integral with a 2-3× improvement in computational speed; see Sec. 4.1 for further discussion. It incorporates all errors on the shape of the occulter natively, and can be written straightforwardly to include position and orientation errors, as well as model waves from off-axis sources, such as an exoplanetary system or the finite diameter of the star. It can also take advantage of the structure of the integral to evaluate the field at multiple wavelengths with little additional computation. The approach of defining the occulter solely in terms of its boundary is a good fit for tolerancing the sorts of errors that are expected for an occulter. Manufacturing errors can be applied by specifically deforming desired edge regions; errors in petal placement can be included by translating and rotating only the boundary points associated with that petal.

The boundary diffraction wave formulation
Most diffraction approaches to propagation through a finite aperture find their basis in the Huygens-Fresnel principle, which expresses the electric field at a point following the aperture as the sum over spherical waves emitting from every point within the aperture. However, a different representation was introduced by Maggi and Rubinowicz which splits the field into two parts, a part which propagates according to geometrical optics and a part which incorporates the effect of diffraction from the boundary. This boundary-diffraction-wave (BDW) representation was codified by the work of Miyamoto and Wolf [14,15], who showed that the boundary integral could be represented for any incident field entirely by a line integral around a vector potential. We will review their analysis briefly here.
Miyamoto and Wolf [14] begin by writing the Kirchhoff formulation (see e.g. [9]) for the diffraction integral from the aperture. Consider a volume of space, bounded by a surface S; the field at any point P within the volume may be expressed as an integral over all surface points Q on S: Here we let s be the distance between P and Q, n be a vector normal to S at point Q, and ∇ Q is the gradient evaluated at point Q. In the case of the occulter, the volume is the space z ≥ 0, and S stretches over the z = 0 plane, as well as having a component at infinity which is assumed to vanish. (Showing that this term vanishes requires some minor additional assumptions; see Sec. 8.3.2 of [9]). It can then be shown [14] that for V of this form, there exists a vector potential W, such that V = ∇ × W, and for plane waves incident on an aperture, that associated vector potential is Here r is the vector from the origin to P; Q is a point within the aperture on S, and r is the vector from the origin to Q. k = 2π/λ is the wavenumber and λ is the wavelength under consideration. The various s-variables become and we assume the form of the plane wave in the vector direction p to be with amplitude A. A coordinate diagram is shown in Fig. 2. Next, we consider for the moment the case of an opaque screen with an opening over the region Ω, a subregion of S with boundary ∂ Ω. We recall here that S is the entire bounding surface of the half-plane z ≥ 0, including the z = 0 plane; Ω is the section of it on which the occulter lies. (We are not considering the occulter case at the moment, but its complement; we will return to the occulter shortly.) In this case, the electric field from Eq. (7) across the aperture is: Our approach is to reduce this integral to a line integral, using Stokes' Theorem. To do this, we consider singularities of the vector potential. In general, W will have singularities somewhere on S for a given P, as otherwise the field in the half-space past the aperture will be zero. In the specific case of the plane wave, the lone singularity occurs at 1 +ŝ · p = 0, which occurs only whenŝ is parallel to the propagation direction of the plane wave; thus, the singularity will fall in the aperture only for points P which fall in the beam as dictated by geometric optics. (Fig. 3 shows an example of this.) Miyamoto and Wolf [15] show that when Stokes' Theorem is applied, we get two distinct cases for the field at point P, depending on the singularity location: U ap (P) = A exp (ikp · r) +U (B) (P) for points inside the extent of the aperture (13) where U (B) (P) is the counterclockwise line integral about the edge of Ω: Here, is a unit vector in the direction tangent to the edge at any point, and d is a differential element of the boundary. This formulation gives the field following an aperture, but in our case we have the opposite: the occulter is opaque, and the surroundings are open. Here we return to Babinet's theorem. As we started with the plane wave in Eq. (11), this implies our occulter field U(P) satisfies U 0 (P) = U(P) +U ap (P), This equation is the general form of the occulter field, in the boundary-diffraction-wave formulation.

Vector potentials for occulter propagation
From this suitably general representation, we can make appropriate substitutions for the specific case of an occulter; these will simplify computation of electric fields later. We define the vector from the origin to P as r = (ξ , η, z); r , the vector from the origin to Q, is defined as r = (x, y, 0). (See again Fig. 2.) We can write the s-variables explicitly as We consider first the case where the plane wave is normally-incident on the occulter, as starlight would be when the system is correctly aligned. Given this, our vector terms can be written as and the potential in Eq. (6) becomes More generally, an off-axis source at an angle of ψ 1 off-axis and ψ 2 from the x-axis will have p = (sin ψ 1 cos ψ 2 , sin ψ 1 sin ψ 2 , cos ψ 1 ), (27) p · r = x sin ψ 1 cos ψ 2 + y sin ψ 1 sin ψ 2 . (28) Before the expression for the vector potential becomes too complicated, we will define three intermediate variables ( f , g, and h) which we can substitute into Eq. (19); these choices will greatly simplify subsequent computations.
Substituting these variables gives: − f cos ψ 2 + g sin ψ 2 cos ψ 1 , a similar relation to Eq. (26). This equation is still exact; for ψ 1 1, the usual case for objects in an exoplanetary system, we can note that and so the potential becomes v( f , g) = ( f sin ψ 2 + g cos ψ 2 cos ψ 1 , (40) − f cos ψ 2 + g sin ψ 2 cos ψ 1 , This form proves to be particularly useful, as ( f 2 + g 2 ) and v · can be calculated without loss of precision, as neither will be of order z for the small ψ 1 case. A typical occulter might have z = 5 × 10 7 m, R = 25m, and |ξ | and |η| ≤ 3m, and k ≈ 10 7 m −1 ; a representative exoplanet target for this occulter might have ψ 1 = 5×10 −7 rad. We note that exponential terms of order kz have been separated from smaller terms; these should be calculated independently, as kz ∼ 10 14 , and combining terms will lose precision in evaluating the exponent.

Efficient implementation
Evaluating the field at each downstream point requires the determination of two parts: the geometric field and the boundary field. We assume, to begin, that we are given a set of N points representing the edge of the occulter which form a simple closed curve, as well as a list of M downstream points and a list of L wavelengths at which the field is to be determined. We also assume that the edge is sufficiently well-sampled that the section of the edge between explicitlyspecified points is well-modeled by a linear segment, as this allows us to use the midpoint rule for numerically approximating the integral. Given the slowly-changing shape over the majority of the petal, this is a good assumption; small regions such as petal tips may be specified more finely.
The geometric field is still a plane wave of form Eq. (11) outside the aperture, and is 0 everywhere inside; finding the geometric field thus becomes a problem of checking whether a downstream point falls behind it or not. We do this by treating the edge as a polygon and finding the winding number of the downstream point (ξ 0 , η 0 ); this can be done efficiently with an O(N) routine such as polywind in Numerical Recipes [16], and holds regardless of the complexity of the occulter shape. This approach also speeds up multiband calculations, as the geometric extent of the shadow determined this way is wavelength-independent. Lateral errors in the occulter-telescope position may be included by adding a ∆ξ and ∆η to all points at the telescope plane, and that should be done first. Note that in the case of an off-axis source, the geometric shadow is shifted laterally as well; in this case, the winding number should be calculated around (ξ 0 − z sin ψ 1 cos ψ 2 , η 0 − z sin ψ 1 sin ψ 2 ). We then have: for points outside the shadow = 0 otherwise.
In some cases-for example, when the telescope aperture remains close to the center of the occulter and the perturbations to the ideal shape of the occulter are small-the winding number may be able to be determined independently. This determination should be done if possible, as even an efficient algorithm can increase runtime significantly for boundaries containing a large number of points.
The boundary field can be evaluated efficiently with midpoint-rule quadrature running around the edge. First, for each of the N line segments running counterclockwise around the occulter edge, we derive the midpoint of line segment j going from (x j , y j , 0) to (x j+1 , y j+1 , 0): and the vector pointing along the segment with (x j+1 , y j+1 ) equal to (x 1 , y 1 ) when j = N. This derivation need only be done once at the beginning, as it holds for all M points and L wavelengths. (We note that, while the midpoint rule has been used for the quadrature, this is certainly not a requirement, and higher order methods could be used; Eq. (39) will still hold. Using a higher-order method may come at the cost of runtime.) For each point in the downstream field, we calculate f , g, and h for all N segments following Eq. (29): and we then need only two intermediate terms: to give us: Note that all of the steps prior to Eq. (49) are wavelength-independent; different values of k may be iterated over in the final step. It is advisable to keep terms of order kz (e.g. kz cos ψ 1 or kh) separate when evaluating to maintain numerical precision. A typical implementation might take the form shown in Algorithm 1. The for-loops over j in particular are well-suited to vectorization where this is supported.

Computational results
In the standard formulation of the Dubra-Ferrari integral (DF), points outside the geometric shadow are calculated with a different set of integrals which require searches for nearest-and further-neighbor points along the boundary for each point at the telescope aperture. To simplify direct comparison between the Dubra-Ferrari approach and the boundary-diffraction-wave approach for this paper, we choose to select a region which lies entirely in the geometric shadow of the occulter. This selection also fixes the winding number to 1 at all points under consideration, letting us bypass an explicit calculation. This choice is not an unreasonable assumption; the telescope will generally be entirely in the geometric shadow for a properly aligned occultertelescope system. Both approaches have similar formulations, with the primary difference coming in the fact that the Dubra-Ferrari integral uses angular coordinates θ centered about (ξ 0 , η 0 ), and so the angular distance between occulter edge points must be rederived for each point at the telescope aperture. This calculation turns out to be one of the major sources of overhead during the DF calculation. To compare the runtimes of the Dubra-Ferrari and BDW implementations, we consider a closed occulter boundary specified by 164000 individual points, corresponding to the occulter shown in Fig. 1. We select subsets of these points, and run Matlab implementations of both algorithms on this subset. Fig. 4 shows the time to implement the propagation as a function of number of boundary points; the BDW algorithm remains 2-3× faster over a wide range of polygon sizes. Fig. 5 shows the intensity profile in a region around the telescope aperture, plotted on a base-10 logarithmic scale. The plots of this intensity profile created with Dubra-Ferrari and BDW are not visually distinguishable; within the telescope aperture, the maximum difference in intensity is 4.8 × 10 −15 .

Conclusion
We present a new approach to modeling fields following planet-finding occulters, based on the boundary-diffraction-wave formulation laid out by Miyamoto and Wolf. It has proved to be simple to modify the shape and recalculate with this approach to include various types of occulter errors, as well as fast to evaluate, improving on the current best method by 2-3× with no loss in precision. An efficient implementation has been laid out, and we hope it proves useful to others in the field.