Abstract
This paper concerns optimization and analysis of telescope-coronagraph systems for direct imaging of exoplanets. In this paper, the coronagraph system, with arbitrary pupil geometry, is theoretically considered, and the governing equation for the pupil design is derived. The method of moments is applied to solve the generalized energy-concentration eigenvalue problem to obtain the optimal pupil apodization and complete sets of orthonormal basis functions for arbitrary pupil geometries. The method yields eigenvalues indicating the fraction of starlight energy encircled in the area of the focal-plane mask (FPM), where starlight can be occulted and/or nulled. In other words, a higher eigenvalue implies less leakage/spillover of light outside of the FPM region and into the planet-discovery zone. Thus, a higher eigenvalue supports better starlight suppression for a given type of coronagraph. This methodology is useful for semi-quantitatively ranking different modes of perturbation with respect to energy spillage in the focal plane independent of coronagraph design details. A model-order–reduction-based sensitivity analysis is conducted to investigate the coupling between different pupil modes induced by aberrations. A pupil mode recovery scheme is presented to offer a theoretically rigorous and computationally efficient approach to reconstruct the optimal pupil mode under an arbitrary phase perturbation. The reconstruction coefficients and recovery-effectiveness factors are derived theoretically and demonstrated numerically. Several numerical examples, including the LUVOIR A and B pupils, are provided to validate and demonstrate the applicability of the proposed methods. The reported methodology enables model-order reduction based on degree of focal-plane energy concentration and reconstruction of optimal pupil apodization vis-á-vis phase aberrations using a precomputed basis set. These features should enhance computational efficiency for coronagraph design and sensitivity analysis.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. INTRODUCTION
Humans have long been observing the night sky and wondering if other planets like Earth, abounding with life, exist elsewhere in the universe. The National Academies of Sciences, Engineering, and Medicine’s recently published Pathways to Discovery in Astronomy and Astrophysics for the 2020’s (aka “Astro2020”) reflects this intellectual yearning by recommending a top-priority NASA flagship “IR/O/UV” mission to search for biosignatures from exoplanets, as well as a major technology-maturation program to enable such a mission [1]. Thus far, Earth is the only place known to harbor life; hence, observations of extrasolar systems resembling our own Earth–Sun system are of utmost interest. If one were to observe the solar system from 10 parsecs away, Earth would be approximately 10 billion times dimmer than the Sun in the visible and separated by only 0.1 arcsec from the Sun (at quadrature). Hence, starlight suppression at extraordinarily small angular separations from the parent star is a key capability in the search for Earth-size exoplanets and life. Coronagraphy is a major approach to achieving such a capability and enabling “exoEarth” observations.
A coronagraph instrument attenuates stellar light in the center of the telescope’s field of view while transmitting off-axis, planetary light to reveal faint features around the star. Figure 1 depicts a typical coronagraph. A plane wave representing starlight arrives upon the telescope pupil ($a$) at normal incidence. A real telescope produces wavefront aberrations and amplitude modifications (e.g., due to structural obscurations and movements, thermal changes, coating nonuniformities, and optical aberrations), depicted as squiggly lines in the figure. These variations and obscurations diffract light to off-axis positions in the image, imposing noise that can hamper planet detection. Therefore, a coronagraph instrument usually includes a set of deformable mirrors (DMs) to correct wavefront aberrations and manipulate wavefronts to enhance starlight suppression. Modern coronagraphs usually employ two DMs, one in a conjugate plane of the telescope pupil to manipulate wavefront phase, and one in a nonconjugate plane to add amplitude manipulation via the Talbot effect. Subsequent optics focus the light onto a focal plane ($b$), where a focal-plane mask (FPM) attenuates the starlight (focused on-axis), while transmitting the off-axis planetary light. An opaque disk with a radius equal to a few diffraction-limited widths of the stellar image is an example of a simple FPM. For a sense of typical scale, an FPM radius of 40 µm corresponds to 2.5 diffraction widths ($F\lambda \sim \lambda /D$, where $F$ is the focal ratio and $D$ is the pupil diameter) for a $f/30$ focal ratio at 0.55 µm wavelength. The FPM inevitably diffracts some amount of starlight; a downstream spatial filter residing in a pupil-conjugate plane ($c$) commonly called the “Lyot stop” attenuates the residual starlight. The Lyot stop can be a simple circular aperture that blocks light diffracted by the FPM at the outer perimeters of the pupil image. Finally, subsequent optics focus the remaining light (containing planetary light and residual starlight) onto a camera to form the observational image ($d$). Note that, in principle, the entrance-aperture, focal, and Lyot-stop planes can all have masks with complex (real and imaginary) transmission functions to manipulate the amplitude and phase of light for maximal attenuation of starlight and minimal attenuation of planetary light.
In practice, attenuating starlight to achieve contrast of $1 \times {10^{- 10}}$ while preserving adequate planetary throughput only $0.1^{\prime \prime}$ off axis is a daunting technological challenge. In the laboratory, Trauger and Traub [2] set a breakthrough milestone in 2007, achieving $6 \times {10^{- 10}}$ contrast (planet-to-star flux ratio) averaged over a circular segment in the image planet between 4 and 10 $\lambda /D$ off axis. Nonetheless, the experiment used monochromatic laser light to simulate starlight; thus, it did not fully demonstrate the capability of observing continuum starlight over a 10%–20% bandwidth to enable spectroscopy, which is necessary in searching for spectral biosignatures. Moreover, the circular segment covered only 30% of the full annular field of view between 4 and 10 $\lambda /D$ [79 of 263 ${(\lambda /D)^2}$ in the image area], thereby increasing the required mission observation time to search a given sky region. To date, Seo et al. [3] published the best overall coronagraphic performance demonstrated in the laboratory: $4 \times {10^{- 10}}$ contrast over a full annular region between 3 and 8 $\lambda /D$ [area $\approx 170\;{(\lambda /D)^2}$], with 10% broadband light. Thus, contrast performance within an order of magnitude from requirement has been achieved in controlled laboratory environments, and designing coronagraphs that can achieve requisite performance in realistic space-telescope environments is a major challenge toward realizing an exoEarth mission.
Within this grand challenge, the need for a large telescope primary mirror and its required stability is a key driver of future mission designs. In addition to determining the photon-collecting area and, thus, achievable signal-to-noise ratio (SNR), a telescope’s entrance-pupil (primary mirror) diameter determines the coronagraph’s achievable “inner working angle” (IWA), which, in turn, affects the number of exoEarths detectable by the mission [4]. The IWA indicates the innermost angular separation from the target star at which a planet can be detected, and it is commonly defined as the innermost angle at which the coronagraph’s optical throughput is 50%. As the diffraction-limited spatial resolution of an optical system is proportional to $\lambda /D$, the IWA scales inversely with aperture diameter. For instance, the above-envisioned exoEarth 10 parsecs away would require an IWA $\le 0.1$ arcsec, and stellar systems farther away require yet smaller IWA. Thus, the larger the telescope, the smaller the IWA (and the larger the photon-collecting area), and the deeper a field can be observed to detect more planets. Also, for a given target planet, a smaller IWA allows detection over a larger range of planetary orbital phase, thereby increasing detection probability as well.
These facts elucidate the fundamental importance of telescope size. However, mitigating deformations of the reflecting surface, which directly produces wavefront errors (WFEs), becomes increasingly challenging as the telescope mirror size and, thereby, structural flexibility increases. Indeed, exoEarth observation requires extraordinary (picometer-level) stability on pupil wavefronts. System analyses indicate that holding WFE (incident on the coronagraph) below $\approx 10\;{\rm pm}$ RMS over a timescale of approximately 10 min is required for spatial frequencies corresponding to the target detection zone in the image plane [5]. The above-mentioned requirement applies to pupil wavefronts immediately before the coronagraph’s FPM. Possible implementations of telescope stabilization system(s) and coronagraph wavefront-correction mechanisms (usually involving DMs) can relax the requirement on the surface stability of upstream optics (including the telescope). Despite such mitigations, WFE involving the primary mirror’s mechanical stability is a major noise term affecting exoEarth detection [6].
The James Webb Space Telescope’s (JWST) successful deployment and operation in space opened a new era of unprecedentedly large space observatories enabled by segmented, deployable mirrors. Such a telescope architecture circumvents limits on the launchable telescope size set by the rocket fairing size. While a segmented primary mirror potentially offers essential benefits for coronagraphy, rigid-body movements and deformations of segments create additional degrees of freedom that need to be mechanically stabilized and included in the analysis of wavefront variations impacting coronagraph performance. Regardless of telescope architecture, stabilizing wavefronts arriving at the coronagraph FPM to a tenth of the hydrogen atom’s diameter is a challenging problem in the co-design and co-engineering of thermo-mechanical, optical, and wavefront-correction systems.
A key element in the co-design cycle is the analysis of the impact that telescope instabilities (vibrations and deformations) have on coronagraph performance, or inversely, deriving telescope requirements based on performance objectives. Such analyses entail conversion of telescope instabilities into optical wavefront aberrations at the coronagraph’s entrance pupil and subsequent propagation through the coronagraph (with main optical planes shown in Fig. 1). The resulting two-dimensional (2D) contrast map informs the coronagraph’s planet-detection capability. Potier et al. recently [7] reported such an analysis estimating sensitivities of mission exoEarth yield to telescope vibrations, whereas Laginja et al. [8] developed an inverse approach to derive telescope tolerance given the contrast requirement. In addition, Tajdaran et al. [9] used a model-order reduction (MOR) method to convert telescope motions into WFE at the exit pupil. MOR is a common feature of all three studies mentioned above. In the case of Potier et al. and Tajdaran et al., the large number of vibrational modes and temporal sampling required MOR to achieve computational feasibility. In Laginja et al., the inverse approach required matrix inversion based on a finite basis set. In particular, Potier et al. employed principal-component analysis (PCA) to derive spatial phase-aberration modes in order of contribution to the phase variance of wavefront changes due to telescope vibrations. They found that a small number accounted for the vast majority of phase variance, enabling reduction of 2D phase-aberration maps to $\le 10$ modal coefficients and, thus, greatly enhancing computational speed. Tajdaran et al. used a conceptually similar approach to derive modes of telescope vibration ranked with respect to phase variance; this enabled performing analyses of the telescope opto-mechanical system with adequate computation speed. Laginja et al. used Zernike polynomials.
In this paper, we report a new approach to derive complete sets of orthonormal basis functions for arbitrary pupil geometries offering some inherent advantages compared to existing bases, such as Zernike polynomials, Bessel–Fourier series, and data-derived approaches (e.g., PCA). Zernike polynomials are widely utilized as a good basis set for representing low-order optical aberrations (e.g., those due to polishing errors and misalignments). However, the unequal distribution of nodal lines over the unit disk produces ringing near the perimeter of the unit disk, especially for high-order Zernike terms. Moreover, other effects, such as those generated by telescope vibrations, are generally not efficiently represented as a superposition of Zernike modes. Bessel–Fourier modes are another well-established basis for the disk domain. However, they are derived with either zero-amplitude or zero-derivative (Dirichlet or Neumann) boundary conditions, which do not necessarily represent the physical system at hand. Furthermore, both Zernike and Bessel–Fourier bases are defined over a circular support area, which does not conform with noncircular apertures, e.g., segmented pupils using an assembly of hexagonal panels. Data-derived approaches circumvent these issues. The data-derived approaches utilized thus far (e.g., by Potier et al. and Tajdaran et al.) have shown enabling enhancements in computational efficiency. However, these modal bases are ranked with respect to phase variances in the telescope’s exit pupil (or, equivalently, the coronagraph’s entrance pupil), not directly ranked with respect to energy distribution/concentration in the focal plane, where the FPM and images of the star and candidate planet(s) reside. Equal variances of phase aberrations in different modes in the pupil produce unequal effects in the focal plane, because different spatial frequencies in the pupil plane diffract light into different locations in the focal plane, and some interfere with planet detection more than others. Therefore, we seek a basis set ranked with respect to energy distribution/concentration in the focal plane.
This paper describes a methodology to solve an energy-concentration eigenvalue problem in deriving a complete set of orthonormal basis functions for arbitrary pupil geometries. The eigenvalues directly correspond to encircled energy in a prescribed focal-plane center zone (e.g., inside the IWA), thereby indicating how much light can be occulted or nulled by the FPM inside the center zone, and how much spills out of the zone to interfere with planet detection. As such, we derive basis sets that are quantitatively ranked with respect to energy concentration in the focal plane’s occulting/nulling zone. The mathematical methodology can be viewed as an extension of Slepian’s seminal work on prolate-spheroidal wave-functions (PSWFs) to arbitrary domains [10–15]. As Soummer et al. [16] have shown, the eigenvalue problem we solve is one that represents a major type of coronagraph design, the apodized-pupil Lyot coronagraph (APLC), and the lowest-order eigenfunction corresponds to the optimal pupil apodization for an APLC, which concentrates the most amount of light into the FPM’s occulting area.
We posit that this method provides efficient modal decomposition for analyzing coronagraph performance or sensitivity because energy concentration in the focal plane is important for the effectiveness of every coronagraph that utilizes a FPM. In addition, we derive a perturbative method to estimate the sensitivity of coronagraph performance to wavefront aberrations. This method also provides intuitive insights into which aberration modes are most degrading of coronagraph performance. A pupil mode recovery scheme is proposed to recompute the optimal eigenmode when the wavefront is arbitrarily perturbed. The recovery effectiveness is also derived as a function of the perturbation intensity and the number of eigenmodes used in the reconstruction.
This paper is organized as follows. The coronagraph pupil design method, the method of moments (MoM) [17–20], is detailed in Section 2. The MOR-based sensitivity analysis [21] is described in Section 3, in which the pupil mode recovery scheme is also described. In Section 4, numerical examples are presented to validate and demonstrate the application of the proposed methods. Section 5 concludes this paper.
2. CORONAGRAPH PUPIL DESIGN
Consider the starlight incident from the normal direction upon the pupil with a pupil function ${\Psi _A}({\boldsymbol \rho})$ defined in a 2D aperture region ${\boldsymbol \rho} \in {\cal P} \subset {\mathbb{R}^2}$. Physically, ${\Psi _A}$ stands for the reflection/transmission coefficient of the reflective/transmissive pupil aperture. The shapes of the apertures can be different in different telescopes, but the areas of the apertures are denoted uniformly as $S$. The image of the starlight on the focal plane ${\boldsymbol k} \in {\cal I} \subset {\mathbb{R}^2}$ is then the Fourier transform of the pupil function. With the application of an FPM with amplitude transmission $\Omega$, the transmitted light ${\Psi _B}({\boldsymbol k})$ can be expressed as
By assuming $M = 1$, Eq. (5) cannot be satisfied precisely, but only approximately. Physically, using a flat-top mask results in a small leakage of light due to the incomplete cancellation between the incident light and the light scattered by the mask. Mathematically, a factor $\lambda$ can be introduced into Eq. (5) to yield the eigenvalue problem
In Eq. (6), the integration over the image plane can be calculated analytically as
To solve Eq. (8) with the MoM [20], a triangulation of the pupil plane ${\cal P}$ is first constructed by discretizing the plane into $T$ small triangular elements ${\cal P} = \bigcup\nolimits_{e = 1}^T {\Delta ^e}$ with a total of ${N_{\cal P}}$ nodes. At node $i$ of the triangulation, a linear basis function ${l_i}({\boldsymbol \rho})$ is defined such that
By solving the GEVP (11), a total of ${N_{\cal P}}$ eigenvalues ${\lambda _n}$, $n = 1,2, \ldots ,{N_{\cal P}}$, with their corresponding eigenvectors ${\{\psi \} ^n}$, can be obtained, and ${N_{\cal P}}$ pupil functions can be expressed as
Due to the mathematical properties of the Slepians [15], they have distinct eigenvalues between 0 and 1. From Eq. (16), the angular dependence of the PSWFs leads to degenerate eigenvalues with a multiplicity of 2 when $m \ne 0$. The eigenvalues can be arranged in a descending order
It is also proved that only the first ${N^{{2D}}}$ eigenvalues are close to unity, and the eigenvectors are spatially concentrated in ${\cal P}$ and bandlimited in ${\cal I}$. For the purpose of the coronagraph design, it is desirable to choose the eigenvector corresponding to the largest eigenvalue as the pupil function since it yields the least residual energy [16]3. PERTURBATION ANALYSIS
The GEVP (11) obtained from the MoM discretization of Eq. (8) can be expressed in matrix form as
where $[\Psi]$ is a square matrix with its $n$th column being the $n$th pupil mode ${\{\psi \} ^n}$, and $[\Lambda]$ is a diagonal matrix with its $n$th diagonal entry being the $n$th eigenvalue ${\lambda _n}$ arranged in a descending order as in Eq. (17).The fact that $[\Lambda]$ is a diagonal matrix indicates that there is no energy coupling between two eigenmodes ${\{\psi \} ^n}$ and ${\{\psi \} ^m}$ where $n \ne m$. This is the result under conditions that the pupil is perfectly aligned in a plane such that light exiting the pupil remains as a plane wave. Unfortunately, in real scenarios, static aberrations and dynamic disturbances introduce phase variations across the pupil. In this section, the effect of such phase variations is studied using perturbation theory to investigate the resulting energy coupling between different eigenmodes. The resultant sensitivity matrices can then be used to adjust the shape of the pupil plane to compensate for the phase difference and restore contrast.
A. Model-Order Reduction
With a general FPM $M({\boldsymbol k})$, the integration kernel in Eq. (7) can be expressed as
Denoting the number of significant eigenvalues as $L$, the eigenvalue matrix can be rewritten as
Following the conventions introduced in Eqs. (24)–(26), the eigenvector matrix can be expressed in block form as
Conditions (37) and (38) constitute a model with a reduced order since the number of significant eigenvalues, $L$, is much smaller than the total number of eigenvalues, ${N_{\cal P}}$, in the original system (22) and (23).
B. Sensitivity Analysis
In this section, the sensitivity analysis is performed when a small perturbation is introduced to DoF $w$ in the $z$ direction that is perpendicular to the pupil plane. This perturbation can represent, for example, a piston motion of a mirror panel or a DM actuator, which results in a small change in the phase of the propagated light. Since in the GEVP (19), both matrices $[A]$ and $[B]$ are calculated on the transverse plane as in Eqs. (12) and (13), their values will not change to the perturbation. The only change caused by the perturbation is the eigenvector matrix $[\Psi]$, which now becomes
where $[P]$ denotes a diagonal perturbation matrix and Here, $k = 2\pi /\lambda$ and $\lambda$ denote the wavenumber and wavelength of the incoming starlight, respectively, and ${z^\prime} = z/\lambda$ stands for the electrical length of the perturbation on panel $w$. Note that even though only phase perturbations are considered in this paper, amplitude perturbations can be easily incorporated in Eq. (41), for instance, using $\Delta p = a{{e}^{- {i}2\pi {z^\prime}}}$, where a denotes the amplitude perturbation. All theories, formulations, discussions, and conclusions remain valid in such a case. In Eq. (39), each column of the perturbed eigenvector matrix $[\tilde \Psi]$ represents the complex starlight intensity reflected by the pupil when the $w$th panel/actuator has a small perturbation ${z^\prime}$ from the remaining DoFs. The objective of the sensitivity analysis is to find the perturbed eigenvalue matrix, $[\tilde \Lambda]$, or its MOR version, $[{\tilde \Lambda _L}]$, and the sensitivity matrix, $[\frac{{\partial \tilde \Lambda}}{{\partial {z^\prime}}}]$, due to the motion of DoF $w$. This is governed by the following modified GEVP: With Eqs. (22), (23) and (39), this GEVP can be rewritten as which leads to the solution of the perturbed eigenvalue matrixApparently, Eq. (45) is a similarity transformation, and $[\tilde \Lambda]$ and $[\Lambda]$ are similar matrices. In Eq. (45),
is called the un-perturbation matrix for a reason that will be explained in the next section. Similarly, Note that $[Q]$ and ${[Q]^*}$ are symmetric matrices. To obtain Eq. (45), ${[\Psi]^{- 1}} = [\Psi {]^T}[B]$ from Eq. (22) and ${[P]^{- 1}} = [P{]^*}$ from Eq. (40) are used, where $*$ denotes the complex conjugate. It is easy to verify that when $z^\prime $ is an integer, $[P] = [Q] = [I]$ and $[\tilde \Lambda] = [\Lambda]$. Also, $[\tilde \Lambda]$ is complex valued, and its off-diagonal elements are in complex-conjugate pairs such that ${\tilde \lambda _{\textit{ij}}} = \tilde \lambda _{\textit{ji}}^*$.The following perspective helps us understand the physical interpretation of $[\tilde \Lambda]$. When an operator operates on a vector, one can compute the result by decomposing the vector onto the eigenvector basis set and stretching the projection along each basis by the corresponding eigenvalue. In the unperturbed case, the eigenvalue matrix $[\Lambda]$ is diagonal and the eigenvectors ${\{\psi \} _i}$ are orthonormal. In the perturbed case, the altered eigenvectors ${\{\tilde \psi \} _i}$ are no longer orthonormal, and we now have nonzero off-diagonal elements in the eigenvalue matrix $[\tilde \Lambda]$. The off-diagonal elements involve projections along directions that are linear combinations of $\{\tilde \psi \}$ basis vectors. The imaginary part of the eigenvalue represent a phase shift in the superposition of basis vectors. Hence, the off-diagonal elements are coupling coefficients spreading energy between modes. Since higher-order modes have higher energy leakage, off-diagonal elements that couple low-order modes with high-orders ones are particularly deleterious to coronagraph contrast.
From Eq. (45), the sensitivity matrix can be calculated as
Note that the trace of a commutator is always zero, so the trace of ${[{\partial _{{z^\prime}}}\tilde \Lambda]_{{z^\prime} = 0}}$ is zero. This simply reflects the fact that coupling energy of a mode onto itself is physically meaningless. Also note that if $[{D_w}] \approx [I]$, i.e., local orthogonality is largely preserved, then the sensitivity is very low. This mathematical fact can provide some physical insights regarding the effect of different types of perturbations. For example, perturbations that affect the entire pupil in an approximately uniform manner (or, on the other extreme, confined to an almost infinitesimal area of the pupil) largely preserve orthogonality and cause only small energy leaks, and thereby small contrast degradation, whereas local perturbations of a significant area tend to be most harmful. The piston of the entire pupil plane is an extreme example of a perfect preservation of orthogonality, causing no contrast degradation. On the other hand, movements of a single mirror panel (a local disturbance with significant area) or uncorrelated movements of multiple panels cause relatively large contrast degradation. In addition, if the pupil geometry has a particular symmetry, one can potentially apply group theory to determine null elements of $[\tilde \Lambda]$ without actually computing $[\tilde \Lambda]$. In the section below on numerical examples, we provide further physical interpretations elucidated by computational results.
C. MOR-Based Sensitivity Analysis
Expressions (45) and (51) can be used directly to calculate the perturbed eigenvalue and sensitivity matrices, respectively. However, the calculation of the full matrices requires knowledge of all eigenvalue–eigenvector pairs and becomes computationally expensive if ${N_{\cal P}}$ is large. In the coronagraph analysis, only the first few eigenmodes are of interest. In this section, the dimensions of these matrices are reduced to account for the first $L$ significant modes only.
To this end, the perturbed eigenvalue and sensitivity matrices are partitioned as
From Eq. (24), the eigenvalue matrix can be approximated as
It should be pointed out that the accuracy of the approximations in Eqs. (55) and (58) depends on the moduli of the insignificant eigenvalues that have been discarded. The larger $L$ is chosen, the more accurate these approximations become. Also, since $[B]$ is sparse, $[P]$ is diagonal, and as will be seen in the next section, $L \ll {N_{\cal P}}$ in practice, the computational and storage costs in the calculation of these approximations are very cheap.
D. Eigenmode Recovery
It has been shown in the preceding sections that under an incident waveform perturbation, the pupil functions $[\tilde \Psi] = [P][\Psi]$ are no longer orthonormal. Consequently, the residual energy (18) will increase, as can be calculated using Eq. (45) or Eq. (55). To compensate for the waveform perturbation, a combination of eigenmodes $[\Psi]$ can be used to construct a desired pupil mode that minimizes the residual energy. This is possible because the eigenmodes obtained from the GEVP (19) are orthonormal basis functions that the pupil aperture supports, which can be used to express an arbitrary function on the pupil. To this end, consider the following perturbed system:
The objective here is to find $[{\Psi ^\prime}]$ such that the eigenvalues $[\Lambda]$ of the above GEVP are the same as those in the unperturbed system (19). Comparing Eq. (61) with Eq. (19), it is apparent that if $[P][{\Psi ^\prime}] = [\Psi]$, the eigenvalues remain the same.Expressing $[{\Psi ^\prime}]$ in terms of the linear combination of $[\Psi]$ as
yields Obviously, in the above expression, $[P]$ is the perturbation matrix and ${[Q]^*}$ is the un-perturbation matrix. From Eq. (62), it can be seen that the $i$th column of $[{\Psi ^\prime}]$, ${\{{\psi ^\prime}\} _i}$, is expressed as the linear combination of all columns of $[\Psi]$, ${\{\psi \} _j}$ ($j = 1,2, \ldots ,{N_{\cal P}}$), with their respective weights $q_{\textit{ji}}^*$, which is the element in the $j$th row and $i$th column of matrix ${[Q]^*}$Using Eq. (22) and noting ${[P]^{- 1}} = [P{]^*}$ yields the following: which is the same expression defined in Eq. (46). The perturbation matrix $[P]$ in Eq. (65) can be made more general to permit different perturbations to multiple panels at the same time. On panels $w \in {\cal W}$, we generalize $\Delta p$ in Eq. (41) to permit different phase perturbations $\Delta {p_w} = {{e}^{- {i}2\pi z_w^\prime}}$, where $z_w^\prime $ can take different values on different panels. Defining $[{E_w}] = {\rm diag}\{{O_1},{O_2}, \cdots ,{I_w}, \cdots ,{O_W}\}$, ${[P]^*}$ can be rewritten as and the un-perturbation matrix ${[Q]^*}$ can be expressed asSince in GEVP (19), ${\{\psi \} _1}$ has the largest eigenvalue ${\lambda _1}$, it is always desirable to obtain a pupil mode ${\{{\psi ^\prime}\} _1}$ such that under a given perturbation $[P]$, the perturbed mode $[P]\{{\psi ^\prime}{\} _1} = \{\psi {\} _1}$. Such a pupil mode can be obtained by combining all eigenmodes $[\Psi]$ using the reconstruction coefficients given by the first column of ${[Q]^*}$, $\{q\} _1^*$
In Eq. (68), all eigenmodes contribute to the un-perturbation process, as can be seen from the inner summation from 1 to ${N_{\cal P}}$ in the second term. When all eigenmodes contribute, the un-perturbed system will have an un-perturbed eigenvalue ${\lambda _1}$ that is the theoretical maximum supported by the pupil aperture and the flat-top mask. When only the first few eigenmodes are used in the construction of the un-perturbed pupil mode, the corresponding eigenvalue will be smaller that ${\lambda _1}$. When only the first $L$ eigenmodes contribute, the un-perturbed pupil mode can be expressed as
Since ${[Q]^*} = [Q{]^{- 1}}$, from Eq. (63) it is clear that $[P][\Psi] = [\Psi][Q]$. For the first eigenmode, it can be shown that
4. NUMERICAL EXAMPLES
In this section, several examples are given to validate and demonstrate the proposed methods in the calculation of the optimal pupil function and the perturbation and un-perturbation behaviors.
A. Circular Pupil
As the first example, a circularly shaped pupil aperture with a radius of 6.75 m is considered to illustrate the eigenvalue and eigenvector supported by an aperture with classical shape and, therefore, with an analytical solution. The accuracy of the pupil function calculation using the MoM will be demonstrated by comparing the numerical solution with the analytical solution. The radius of the flat-top mask plane is chosen as $K = 0.4\pi\; {\rm rad/m}$ such that $N$, which stands for the spatial frequency in terms of number of cycles across the aperture diameter, is 2.7. With these settings, the planar Shannon number ${N^{{2D}}}$ is 17.99, indicating that the first 18 eigenvalues are close to unity. By solving the GEVP (11) with the MoM, the first 24 eigenvalue–eigenvector pairs (“eigen pairs” hereafter) are presented in Fig. 3. Clearly, for every $m \ne 0$ in Eq. (16), there are two degenerate modes for each ${R_{\textit{mn}}}$, as can be seen from this figure. It is also very clear from this figure that the eigenvalues decrease rapidly as their indices increase. As the eigenvalues decrease, their corresponding eigenvectors shown in the color plots present an increasing “leakage” of energy to the outside of the circular aperture indicated by the gray circle. In the design of a coronagraph, it is desirable that all the incident starlight gets shielded by the Lyot stop that has the same shape as the pupil aperture. As a result, only the first few eigen pairs are of interest in the design.
To demonstrate the accuracy of the numerical solution, the eigenvectors obtained by the MoM are compared to the analytical solutions, the generalized PSWFs in the radial direction, which is the Slepians that satisfy Eq. (15). The comparisons between analytically and numerically obtained PSWFs are presented in Fig. 4. In particular, solutions of ${R_{00}}$, ${R_{01}}$, and ${R_{02}}$ are shown and compared in Fig. 4(a); ${R_{11}}$, ${R_{12}}$, and ${R_{13}}$ are shown and compared in Fig. 4(b); ${R_{22}}$, ${R_{23}}$, ${R_{33}}$, and ${R_{34}}$ are shown and compared in Fig. 4(c); and ${R_{44}}$, ${R_{55}}$, ${R_{66}}$, and ${R_{77}}$ are shown and compared in Fig. 4(d). It is very clear from these figures that the numerical and analytical solutions for all of these PSWFs are on top of each other, indicating excellent accuracy of the proposed numerical solver.
B. LUVOIR-A Telescope Coronagraph
When a pupil aperture with an irregular shape is encountered, there is generally no analytical solution available. In such a case, the proposed numerical solver becomes extremely useful. In this section, the Large UltraViolet Optical InfraRed (LUVOIR) Surveyor is considered and modeled. LUVOIR is a mission concept for a highly capable, multi-wavelength space observatory with ambitious scientific goals, which was studied by NASA [5]. With its large telescope aperture, one of LUVOIR’s major goals is to characterize a wide range of exoplanets that might be habitable or even inhabited. In this section, the LUVOIR surveyor structure A (LUVOIR-A) telescope aperture is considered, and in the next section, structure B (LUVOIR-B) is considered.
The LUVOIR-A pupil aperture is made of 120 identical hexagonal panels, which have an edge length of 0.7061 m and are separated from the neighboring panels by a 6-mm gap for them to move freely in the direction perpendicular to the aperture. The outer diameter of the entire structure is about 15.0 m, and the inner diameter is about 13.5 m. The shadows from the telescope supporting structures are modeled as 0.235-m-wide gaps. Shown in Fig. 5(a) is the pupil structure and panels with ID labeled from 1 to 120. The triangulation of the pupil plane results in 10,998 triangular elements and 7118 nodes on the plane, as shown in Fig. 5(b). The zoomed-in view in Fig. 5(d) shows the 6-mm gaps between adjacent panels. In this example, the radius of the flat-top mask is once again chosen as $K = 0.4\pi$, resulting in ${N^{{\rm 2D}}} = 18.97$. With the MoM, the first 24 eigenmodes along with their corresponding eigenvalues are calculated and presented in Fig. 6. In this figure, only the eigenmode distributions on the pupil are shown. It can be seen that the small gaps between the panels do not interfere with the eigenmode distributions but the wider shadows can split the eigenmodes. It is also seen that the eigenvalues decrease rapidly along with a decreased energy concentration of the eigenmodes on the pupil. The distribution of all 7118 eigenvalues calculated using the MoM is presented in Fig. 7(a) in log scale. Clearly, the eigenvalues decrease very fast. Out of the total of 7118 eigenvalues, there are only 60 eigenvalues larger than ${10^{- 4}}$, as indicated by the blue horizontal and vertical dotted lines in this figure. Also note that after about 200 eigenvalues, the decrease of eigenvalues flattens out, as indicated by the black vertical dashed line. This serves as the foundation as the MOR-based sensitivity analysis.
Next, the sensitivity analysis of the pupil eigenmode under small perturbation is performed. While the numerical values are not presented, the observations are summarized here. When small perturbations are introduced to different panels, the resulting perturbed eigenvalue matrices $[{\tilde \Lambda _L}]$ are fully populated and complex valued, indicating the fact that the energy of an eigenmode ”spills out” to other eigenmodes under perturbations. The intensity of each matrix element indicates the coupling strength between two eigenmodes. The purely imaginary matrix ${[{\partial _{{z^\prime}}}{\tilde \Lambda _L}]_{{z^\prime} = 0}}$ is antisymmetric. Also, when a panel located at the “hotspot” of the an eigenmode is perturbed, a larger change in the perturbed eigenvalue matrix and a stronger energy coupling are observed. This example provides insight and intuitive understanding of the effects of panel perturbation based on the perturbed panel locations and the eigenmode distributions.
C. LUVOIR-B Telescope Coronagraph
In this section, LUVOIR telescope mirror structure B is investigated, and its pupil coronagraph is modeled and simulated. Shown in Fig. 5(c) is the geometry of the LUVOIR-B pupil aperture, which is made of 55 hexagonal panels. The edge length of each panel is 0.5514 m. The edges of its outermost panels are trimmed down to a circle with a radius of 3.35 m. Between every two adjacent panels there is a 6-mm gap as shown in Fig. 5(d). This aperture is discretized into a total of 5376 triangular elements and 3463 nodes, resulting in ${N_{\cal P}} = 3463$ unknowns in the GEVP (11). The radius of the mask plane is chosen as $K = 0.8\pi$ rad/m with a flat-top mask function $M = 1$ in this example. With these settings, the Shannon number ${N^{{2D}}}$ is 17.50, indicating that the first 18 eigenvalues are close to unity. By solving the GEVP (11), the first 24 eigen pairs are presented in Fig. 8. Due to the circular symmetry, degenerate eigenmodes present in pairs with a 90° rotation. The slight differences between the degenerate eigenvalues (for example, see degenerate modes 7 and 8, 18 and 19, and 20 and 21 in Fig. 8) are caused by the relative locations of the segment gaps with respect to the specific eigenmodes. The distribution of all 3463 eigenvalues calculated using the MoM is presented in Fig. 7(b) in log scale. Out of the total of 3463 eigenvalues, there are only 53 eigenvalues larger than ${10^{- 4}}$, as indicated by the blue horizontal and vertical dotted lines in this figure.
To better understand the effects of the segment gaps, the GEVP of a circular pupil aperture with a radius of 3.35 m is solved, and the resulting eigenmode distributions are identical to those shown in Fig. 8. It is interesting and worthwhile to notice the eigenvalue differences of the corresponding eigenmodes in these two cases. As summarized in Table 1, the eigenvalues of the LUVOIR-B pupil are always a little smaller that those of the circular pupil, which is caused by the energy leakage due to the tiny 6-mm gaps in the LUVOIR-B pupil. Specifically, the surface area of the circular pupil is $35.26\;{\rm m}^2$ and that of the LUVOIR-B pupil is $34.82\;{{\rm m}^2}$, representing a 1.25% area reduction due to the segment gaps. The energy leakage due to the gaps, taking the first eigenmode as an example, can be calculated as follows: (1) the residual energy of the circular pupil, according to Eq. (18), is ${(1 - {\lambda _1})^2} = (1 - {0.999991)^2} = 8.1 \times {10^{- 11}}$; (2) the residual energy of the trimmed LUVOIR-B pupil is ${(1 - {\lambda _1})^2} = (1 - {0.987543)^2} = 1.6 \,\times\def\LDeqbreak{} {10^{- 4}}$. As a result, the energy leakage due to the segment gaps is $1.6 \times {10^{- 4}}/8.1 \times {10^{- 11}} = 1.98 \times {10^6}$, or about 2 million times. This can be observed more clearly in Fig. 9, where the images of the first eigenmodes of the circular and LUVOIR-B pupil apertures are presented. These images are calculated by Fourier transforming the eigenmodes to the corresponding Fourier space. In these figures, the white circles in the center indicate the location and size of the opaque flat-top masks, which are placed to shield as much direct starlight as possible. Clearly, in Fig. 9(a) the intensity of the image outside the white circle is very small compared to that in Fig. 9(b), in which the hexagonal pattern of high intensity speckles is present. The presence of such a pattern results from the segment gaps, which requires DM action to mitigate.
As demonstrated in the preceding section, when the incident starlight is no longer in phase, the original pupil eigenmode loses its effectiveness in shielding the starlight and results in a larger energy leakage. To compensate for the perturbed starlight wavefront, several eigenmodes are needed to recover a good starlight blockage. As indicated by Eqs. (67) and (68), in principle, the weighted sum of all eigenmodes is needed to recover a certain eigenmode under perturbations. The combination weights can be calculated using the perturbation intensities $\Delta {p_w}$ on panel $w$ and the local orthogonal matrix ${D_w}$. While $\Delta {p_w}$ changes in different cases, ${D_w}$ does not. In this section, the elements of the local orthogonal matrix ${D_w}$ are examined to demonstrate how important each eigenmode is in the reconstruction process. Since the first eigenmode has the largest eigenvalue, it is always the desired mode to be recovered. Thus, only the elements in the first column of ${D_w}$ are presented. Shown in Fig. 10 are the local orthogonal matrix elements ${d_{w,j1}}$ calculated on panels $w = 5$, 20, and 28 of the LUVOIR-B pupil aperture. The subscript $j$ denotes the indices of the eigenmodes with their corresponding eigenvalues arranged in a descending order. Very important observations can be made from this figure. First, most eigenmodes are needed if the recovery is to be made perfect since they all have nonzero ${d_w}$ values. Second, for the LUVOIR-B pupil, only the first 200 eigenmodes contribute significantly to the mode recovery. As can be seen, in all three cases under consideration, the ${d_w}$ values of eigenmodes 201 to the last are in general smaller that those of the first 200 modes. This observation correlates well with Fig. 7(b), where the eigenvalues decrease logarithmically until around the 200th eigenvalue, as indicated by the black vertical dashed line. After the 200th eigenvalue, such a decrease flattens out. Third, the ${d_w}$ values on panel 5 are significantly smaller than those on panel 20, which are smaller than those on panel 28. This can be explained by referring to the panel locations shown in Fig. 5(c). Clearly, panel 5 is located on the edge of the pupil where the first eigenmode, which is to be recovered, has the smallest values. Panel 20 is closer to the center, and panel 28 is the center panel of the pupil aperture, which has the largest eigenmode values. It can then be concluded that the importance of each panel is different, which depends on their location and the values of the desired eigenmode on that panel. Lastly, the number of nonzero ${d_w}$ values on panel 28 is much smaller than those on the other two panels, due to the fact that all eigenmode distributions on the center panel 28 are symmetrical, as can be seen in Fig. 8. Such symmetry results in a better local orthogonality between different eigenmodes.
Finally, the eigenmode reconstruction
is performed and the ERF on for the first eigenmode, ${{\rm ERF}_1}$, is calculated in different cases and presented in Figs. 11 and 12. Figure 11 presents the changes of ${{\rm ERF}_1}$ when a single panel, $w = 5$, 20, or 28, is perturbed by a displacement ${z^\prime}$. In these figures, the horizontal axis stands for the total number of eigenmodes, starting from the first one, that are used to reconstruct the desired eigenmode ${\{\psi \} _1}$. Clearly, for each ${z^\prime}$ case, the more the eigenmodes are used in the reconstruction, the larger the ERF becomes. Only when all eigenmodes are used can a 100% recovery can be achieved. As ${z^\prime}$ increases, the ERF decreases in all cases, which indicates that more eigenmodes are needed to have a better eigenvalue recovery and a smaller energy leakage. It is also noted by comparing Figs. 11(a)–11(c) that with the same displacement ${z^\prime}$, the ERF reduces most significantly on panel 28, followed by panel 20, and least significantly on panel 5. This is due to the same reason explained earlier, that panel 28 sits in the center of the pupil aperture where the first eigenmode has the largest value. Figure 12 is presented to demonstrate the recovery capability of formulation (73) when multiple panels are disturbed with different displacements. In this demonstration, panels 5, 20, and 28 are perturbed simultaneously, with different sets of perturbation values ${z^\prime}$ labeled in the figure. Similar observations to those in Fig. 11 can be made. It is also observed that when multiple panels are perturbed with the same ${z^\prime}$, the resulting ERF decreases compared to the case in which only one panel is perturbed. This agrees with the physical intuition where uncorrelated movements of multiple panels at the same time will increase energy leakage.5. CONCLUSION
In this paper, the governing equation for coronagraph pupil design and optimization has been derived rigorously (in a fashion similar to that of Soummer et al. [16]), which leads to a GEVP for the calculation of the pupil modes. The MoM has been proposed to solve the GEVP. The development and implementation of this numerical method offer a powerful approach to calculating the optimal pupil-apodization function concerning a pupil aperture with an arbitrary shape. Based on the MOR technique, the sensitivity analysis has been performed when a small perturbation is introduced to a pupil panel or a deformable mirror pixel. The closed forms of the perturbed eigenvalue and the sensitivity matrices have been derived analytically, which permit an efficient calculation of these matrices and provide better physical and engineering insights. An eigenmode recovery scheme has been investigated based on the pupil mode analysis. A recovery formulation has been derived to provide the eigenmode reconstruction coefficients when multiple panels are perturbed simultaneously and with different displacements. The reconstruction effectiveness can then be calculated based on the intensities of the perturbations and the number of eigenmodes used in the reconstruction. This offers a theoretically rigorous and computationally efficient approach to explore which pupil panels or deformable mirror pixels have more significant effects under perturbations and which pupil eigenmodes are more important, both qualitatively and quantitatively, in a recovery process. The ability to recompute optimal solutions via-á-vis aberrations using precomputed basis functions can potentially facilitate real-time coronagraph optimizations as well. In future efforts, we will investigate the optimization and implementation of the FPM function, and will develop a numerical scheme that optimizes the pupil-apodization and FPM functions concurrently.
Funding
Jet Propulsion Laboratory.
Acknowledgment
Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data and pseudo-code underlying the results presented in this paper are not publicly available at this time but may be obtained from the first author (Prof. Su Yan, su.yan@howard.edu) upon request.
REFERENCES
1. NASEM, Pathways to Discovery in Astronomy and Astrophysics for the 2020s (National Academies, 2021).
2. J. T. Trauger and W. A. Traub, “A laboratory demonstration of the capability to image an earth-like extrasolar planet,” Nature 446, 771–773 (2007). [CrossRef]
3. B.-J. Seo, K. Patterson, K. Balasubramanian, B. Crill, T. Chui, D. Echeverri, B. Kern, D. Marx, D. Moody, C. M. Prada, G. Ruane, F. Shi, J. Shaw, N. Siegler, H. Tang, J. Trauger, D. Wilson, and R. Zimmer, “Testbed demonstration of highcontrast coronagraph imaging in search for earth-like exoplanets,” Proc. SPIE 11117, 111171V (2019). [CrossRef]
4. C. C. Stark, A. Roberge, A. Mandell, and T. D. Robinson, “Maximizing the exoearth candidate yield from a future direct imaging mission,” Astrophys. J. 795, 122 (2014). [CrossRef]
5. LUVOIR Mission Concept Study Team, LUVOIR Final Report (NASA, 2019).
6. B. Nemati, H. P. Stahl, M. T. Stahl, G. J. Ruane, and L. J. Sheldon, “Method for deriving optical telescope performance specifications for earth-detecting coronagraphs,” J. Astron. Telesc. Instrum. Syst. 6, 039002 (2020). [CrossRef]
7. A. Potier, G. Ruane, C. Stark, P. Chen, A. Chopra, L. Dewell, R. J. Parramon, A. Nordt, L. Pueyo, D. Redding, A. E. Riggs, and D. Sirbu, “Adaptive optics performance of a simulated coronagraph instrument on a large, segmented space telescope in steady state,”J. Astron. Telesc. Instrum. Syst. 8, 035002 (2022). [CrossRef]
8. I. Laginja, R. Soummer, L. M. Mugnier, L. Pueyo, J.-F. Sauvage, L. Leboulleux, L. Coyle, and J. S. Knight, “Analytical tolerancing of segmented telescope co-phasing for exo-earth high-contrast imaging,” J. Astron. Telesc. Instrum. Syst. 7, 015004 (2021). [CrossRef]
9. K. Tajdaran, L. D. Dewell, M. S. Jacoby, A. A. Nordt, and J. Z. Lou, “Line-of-sight and wavefront error dynamic stability during coronagraphic imaging for a 6.7-meter inscribed diameter UVOIR segmented telescope with non-contact pointing and vibration isolation,” Proc. SPIE 12180, 121802L (2022). [CrossRef]
10. D. Slepian and H. O. Pollak, “Prolate spheroidal wave functions, Fourier analysis and uncertainty–I,” Bell Syst. Tech. J. 40, 43–63 (1961). [CrossRef]
11. H. J. Landau and H. O. Pollak, “Prolate spheroidal wave functions, Fourier analysis and uncertainty–II,” Bell Syst. Tech. J. 40, 65–84 (1961). [CrossRef]
12. H. J. Landau and H. O. Pollak, “Prolate spheroidal wave functions, Fourier analysis and uncertainty–III: The dimension of the space of essentially time- and band-limited signals,” Bell Syst. Tech. J. 41, 1295–1336 (1962). [CrossRef]
13. D. Slepian, “Prolate spheroidal wave functions, Fourier analysis and uncertainty–IV: Extensions to many dimensions; generalized prolate spheroidal functions,” Bell Syst. Tech. J. 43, 3009–3057 (1964). [CrossRef]
14. D. Slepian, “Prolate spheroidal wave functions, Fourier analysis, and uncertainty–V: The discrete case,” Bell Syst. Tech. J. 57, 1371–1430 (1978). [CrossRef]
15. F. J. Simons and D. V. Wang, “Spatiospectral concentration in the Cartesian plane,” Int. J. Geomath. 2, 1–36 (2011). [CrossRef]
16. R. Soummer, L. Pueyo, A. Ferrari, C. Aime, A. Sivaramakrishnan, and N. Yaitskova, “Apodized pupil Lyot coronagraphs for arbitrary apertures. II. Theoretical properties and application to extremely large telescopes,” Astrophys. J. 695, 695–706 (2009). [CrossRef]
17. R. F. Harrington, Field Computation by Moment Methods (Macmillan, 1968).
18. W. C. Chew, J.-M. Jin, E. Michielssen, and J. M. Song, eds., Fast and Efficient Algorithms in Computational Electromagnetics (Artech House, 2001).
19. S. Yan and J.-M. Jin, “Self-dual surface integral equations for electromagnetic scattering from IBC objects,” IEEE Trans. Antennas Propag. 61, 5533–5546 (2013). [CrossRef]
20. S. Yan, P. Chen, M. I. Wade, T. L. Gill, and J. T. Trauger, “Method of moments based coronagraph pupil design for exoplanet exploration,” in International Applied Computational Electromagnetics Society (ACES) Symposium [online] (2021).
21. S. Yan, P. Chen, M. I. Wade, T. L. Gill, and J. T. Trauger, “Model order reduction based sensitivity analysis of coronagraph pupil aperture,” in Proceedings of IEEE Antennas and Propagation Symposium (Denver, 2022).
22. J.-M. Jin, The Finite Element Method in Electromagnetics, 3rd ed. (Wiley, 2014).
23. S. Zhang and J.-M. Jin, Computation of Special Functions (Wiley, 1996).