1. INTRODUCTION

The challenge of discovering habitable planets beyond our solar system has motivated astronomers to develop a diverse repertoire of exoplanet detection techniques. Broadly speaking, transit photometry, radial velocity, gravitational microlensing, and astrometry methods all monitor perturbations to the brightness, position, and spectrum of a prospective host star over time to infer the presence and dynamics of a faint orbiting companion [1]. While these methods have enjoyed great success in detecting exoplanets, contributing over 5500 confirmed discoveries to date [2,3], they fundamentally rely on indirect observations, which provide limited information about more detailed planetary features. Remotely characterizing atmospheric composition, weather patterns, surface temperature, and surface gravity is crucial for understanding extrasolar chemical environments and identifying potential biosignatures [4,5].

By comparison, direct imaging techniques aspire to spatially observe/resolve orbiting exoplanets, providing more comprehensive planetary data [6,7]. However, direct imaging faces two compounding obstacles: photon shot noise and diffraction. Exoplanets are extremely faint compared to their host stars, with relative brightness factors ranging from ${10^{- 5}}$ for Hot Jupiters to ${10^{- 11}}$ for exo-Earths in the habitable zone. Moreover, the distance between an exoplanet and its host star often falls below the optical resolution capabilities of current space-based telescopes, residing in the so-called “sub-diffraction regime” [8]. Under these circumstances, light from the exoplanet overlaps with prominent diffraction features of the host star. This overlap, combined with the overwhelming shot noise generated by the bright star, effectively renders the exoplanet undetectable with conventional telescopes. Developments in coronagraphy techniques have succeeded in nulling an on-axis point-like star so that only light from the exoplanet reaches the detector [9–14]. In this way, state-of-the-art coronagraphs suppress photon shot noise intrinsic to measuring classical states of light, thereby enhancing the signal-to-noise ratio (SNR) of exoplanet signatures [15,16].

Inspired by new insights in passive super-resolution imaging [17–20], we recently reported the quantum information limits for exoplanet detection and localization [21]. Our findings revealed that these limits are achieved by a direct-imaging coronagraph that exclusively rejects the fundamental mode of the telescope. In contrast, current state-of-the-art coronagraphs discard information-bearing photons in higher-order modes [22], resulting in sub-optimal performance over the sub-diffraction regime, as illustrated in Fig. 1. Quantum-optimal coronagraphs, however, preserve information at sub-diffraction star–planet separations, where an abundance of exoplanets is expected to reside given current statistical models [23–26].

 

Fig. 1. Comparisons between (a) the Chernoff exponent and (b) the Fisher information per photon as a function of the normalized star–planet separation ${r_\Delta}/\sigma$ for various coronagraph designs relative to the quantum limits (black) of exoplanet detection and localization at high contrasts. The perfect coronagraph (dashed yellow) is the focus of this work and is known to theoretically saturate both quantum limits in [21]. State-of-the-art systems, such as the phase-induced amplitude apodization complex mask coronagraph (PIAACMC) [11] and the vortex coronagraph [9], are notably sub-optimal in the sub-diffraction regime. Spatial mode demultiplexing (SPADE) [17] is an alternative quantum-optimal measurement shown for reference, though it does not qualify as a direct imaging coronagraph.

Download Full Size | PDF

Emerging spatial mode demultiplexing (SPADE) coronagraph platforms, such as photonic lanterns [27,28], offer an alternative sensing paradigm that involves measuring optical power in different spatial modes. Although this approach provides a promising direction for exoplanet spectroscopy, it departs from the imaging paradigm, which we argue remains desirable in coronagraphy. Images immediately provide the context and composition of a scene, thereby relaxing the amount of prior information needed to interpret and process measurements. Consider, for example, a prospective host star surrounded by multiple luminous sources (e.g., systems of exoplanets or exozodiacal dust). The exitance distribution around the star is impossible to uniquely determine from modal power measurements alone, without additionally measuring the relative phase between modes or sorting multiple non-degenerate mode bases [29–31]. In comparison, a direct imaging coronagraph readily provides information on the spatial distribution of luminous sources around a star without added measurement complexity.

In this work, we propose a direct imaging coronagraph using a spatial mode sorter implemented with a multi-plane light converter (MPLC) [32–34]. To the best of our knowledge, this is the first experimental verification of a coronagraph design that theoretically saturates the quantum limits for exoplanet discovery tasks. By applying a maximum likelihood estimator (MLE) to the images collected with our benchtop setup, we localize an artificial exoplanet at sub-diffraction separations from an artificial host star with a contrast ratio of 1000:1. Denoting the Rayleigh diffraction limit of our imaging system as $\sigma$, our experiment achieves an absolute localization error below $0.03\sigma$ over the separation range $[0, 0.6]\sigma$ with the MLE.

2. METHODS

A. Experimental Design

Figure 2(a) illustrates the working principles of a quantum-optimal coronagraph design based on two cascaded mode sorters. The first mode sorter decomposes the incident optical field into a point spread function (PSF)-adapted transverse spatial mode basis [35] in order to isolate and eliminate photons residing in the fundamental mode. The second sorter inverts the mode decomposition, coherently recombining light in the residual modes to form an image of the exoplanet on a detector array. This scheme can be viewed as spatial mode filtering.

 

Fig. 2. (a) Conceptual design of a quantum-optimal direct-imaging coronagraph based on spatial mode sorting. The forward mode sorter demultiplexes the field incident on the image plane into a PSF-adapted basis. Light in the fundamental mode is rejected by a mask at the sorting plane, while light in the remaining modes propagates freely to an inverse mode sorter. The inverse mode sorter coherently recombines (multiplexes) the remaining modes to form an image at the detector. (b) A simplified schematic of the experimental implementation of a quantum-optimal coronagraph using an MPLC mode sorter and non-reciprocal polarization elements. The MPLC is comprised of a sequence of diffractive phase elements designed to (de)multiplex Fourier–Zernike modes $\{{{\boldsymbol \psi}_0},{{\boldsymbol \psi}_1},{{\boldsymbol \psi}_2},{{\boldsymbol \psi}_3}\}$, where ${{\boldsymbol \psi}_0}$ is the PSF of the imaging system. A maximum likelihood estimator is applied on the captured image to localize the position of a sub-diffraction exoplanet.

Download Full Size | PDF

Our experimental setup, shown in Fig. 2(b), emulates the working principles of Fig. 2(a) by double-passing the optical field through a single mode sorter implemented on a three-plane MPLC, as detailed in [36] and [37]. During the forward pass, the MPLC spatially demultiplexes the optical field in the Fourier–Zernike modes $\{{{\boldsymbol \psi}_{{nm}}}(\vec r)\}$ (defined in Appendix A), which constitute a PSF-adapted basis for circular apertures, and focuses light on each mode to spatially resolve Gaussian spots on the sorting plane. The spot corresponding to the fundamental mode is directed to the opening of a pinhole mirror and absorbed at a beam dump. The remaining modes reflect off the pinhole mirror and are sent backward through the mode sorter. The unitary nature of spatial mode sorting reverses the mode transformation during the backward pass. Non-reciprocal polarization elements split the optical path for the forward (pre-nulling) and backward (post-nulling) pass, sending the filtered field to a detector. Through this process, the field at the detector plane is identical to the field at the focal plane minus the optical contribution from the fundamental mode. The $4f$ imaging system used in our setup is characterized by a circular aperture diameter $D = 400\;{\unicode{x00B5}\text{m}}$ and a focal length $f = 200\;\text{mm}$ operating at wavelength $\lambda = 532\;\text{nm}$, yielding a Rayleigh resolution of $\sigma = 1.22\lambda f/D = 324\;{\unicode{x00B5}\text{m}}$ on the object plane.

We focus on demonstrating exoplanet localization at sub-Rayleigh star–planet separations. The regime where quantum-optimal coronagraphs offer the greatest theoretical advantage over existing high-performance coronagraph designs. To sample this separation regime, we align the coronagraph to a bright on-axis point source (artificial star) and vertically step the position of a second dim point source (artificial exoplanet) ${\vec r_e} = (0,{y_e})$ over the discrete domain ${y_e} \in {\mathscr Y} = [- h:\Delta : + h]$, with endpoints $h = 0.85\sigma$ and sampling step size $\Delta = 0.0215\sigma$.

Light from a sub-diffraction exoplanet couples predominantly to a small set of lower-order Fourier–Zernike modes. We therefore configured the MPLC to sort a truncated basis $\{{{\boldsymbol \psi}_0},{{\boldsymbol \psi}_1},{{\boldsymbol \psi}_2},{{\boldsymbol \psi}_3}\} = \{{{\boldsymbol \psi}_{00}},{{\boldsymbol \psi}_{1 - 1}},{{\boldsymbol \psi}_{20}},{{\boldsymbol \psi}_{22}}\}$ depicted in Fig. 2(b), where ${{\boldsymbol \psi}_0}(\vec r)$ is the fundamental mode of the imaging system. Collectively, these modes contain a majority of the energy in the field generated by a sub-Rayleigh companion traversing the $y$-axis, as shown in Fig. 3(a). We define the nominal region of support for this truncated basis to be $|{y_e}| \le 0.6\sigma$. For any exoplanet located in this region, ${\gt}{95}\%$ of the optical energy couples to the truncated mode basis.

 

Fig. 3. (a) Theoretical photon detection probabilities (coupling strengths) in each mode as a function of the off-axis source location for the truncated set of Fourier–Zernike modes in the absence of cross-talk. The nominal “region of support” is given by off-axis positions in the domain $[0,0.6]\sigma$ within which ${\gt}95\%$ of the light energy couples to the truncated mode basis. (b) Distribution of the classical Fisher information in each mode as a function of the off-axis position. The classical Fisher information of the truncated mode basis remains within $10\%$ of the quantum Fisher information limit over the domain ${\sim}[0,0.4]\sigma$. Outside this domain, the information concentrates in higher-order modes that are ignored by our mode sorter.

Download Full Size | PDF

 

Fig. 4. (a) Cross-talk calibration data (dotted) depicting the relative optical intensity measured in each target region of the sorting plane as a function of off-axis source position. Solid curves correspond to theoretically predicted mode intensities under a constrained least-squares fit of the unitary cross-talk matrix to the calibration data. (b) A visual representation of the least-squares fitting algorithm used to estimate the experimental cross-talk matrix under the constraint that $\Omega$ be a unitary matrix (further details in Appendix D). (c) Example calibration images of the optical intensity distribution measured at the sorting plane for various off-axis source locations. The total intensity of each mode was taken to be the integrated intensity within its designated circular region on the detector.

Download Full Size | PDF

The mode set used in our experimental setup was chosen specifically for the 1D localization of a sub-diffraction exoplanet displaced along the $y$-axis. However, configuring the setup for 2D localization by augmenting/modifying the truncated basis to include modes that are sensitive to displacement along the $x$-axis is straightforward. For example, we could have chosen to sort the modes $\{{{\boldsymbol \psi}_{00}},{{\boldsymbol \psi}_{1 - 1}},{{\boldsymbol \psi}_{11}},{{\boldsymbol \psi}_{20}}\}$, which affords sub-diffraction sensitivity to horizontal displacement through the ${{\boldsymbol \psi}_{11}}$ mode. Expanding the basis to more than four modes was found to significantly degrade the cross-talk of the mode sorter due to the limited number of phase masks available on our programmable MPLC. In principle, introducing more masks would allow one to sort more modes and temper modal cross-talk, enabling access to larger star–planet separations and higher contrasts.

B. Measurement Model

We construct a statistical forward model for the images captured with our experimental coronagraph. To circumvent erroneous interference effects that would arise if both the star and exoplanet were simultaneously illuminated by the coherent laser source in our setup, we synthesize a measurement of the star–planet scene by adding measurements of each body illuminated independently, so ${\boldsymbol Y} = {{\boldsymbol X}_s} + {{\boldsymbol X}_e}$, where ${{\boldsymbol X}_s}$ and ${{\boldsymbol X}_e}$ are random vectors denoting images captured with our coronagraph when illuminating the artificial star and exoplanet, respectively. Importantly, the statistics of the synthetic measurement ${\boldsymbol Y}$ are identical to those expected for an incoherent star–planet pair.

The ${{\boldsymbol X}_s}$ and ${{\boldsymbol X}_e}$ are themselves constructed from independent realizations of a single-shot measurement. Specifically, ${{\boldsymbol X}_s} = \sum\nolimits_{i = 1}^{{N_s}} {{\boldsymbol X}^{(i)}}(\vec 0)$ and ${{\boldsymbol X}_e} = {\boldsymbol X}({\vec r_e})$, where ${\vec r_e}$ is the position of the exoplanet, and the choice of ${N_s}$ sets the star–planet contrast. Meanwhile, the single-shot measurement ${\boldsymbol X}({\vec r_0}) \in {\mathbb{N}^M}$ is a random vector denoting the number of photons measured at each pixel of an $M$-pixel detector array over a fixed exposure time $T$ when imaging a single point source located at position ${\vec r_0}$. The number of photons measured in each pixel is statistically independent Poisson random variables such that

By using the additive property of Poisson random variables, the complete statistical model for the synthetic measurement is

C. Contrast Sensitivity Limit and Image Quality

The cross-talk matrix $\Omega$ of the mode sorter determines the performance of our coronagraph in two critical ways. First, energy leakage into and out of the fundamental mode sets a heuristic limit on the maximum accessible star–planet contrast. Second, interference with residual light in the fundamental mode causes distortions to the image of the exoplanet. Figure 4 shows calibration measurements made at the sorting plane of our system for determining the phase and amplitude of each entry in the cross-talk matrix $\Omega$ depicted in Fig. 5

 

Fig. 5. (a) Amplitude and (b) phase of the experimental cross-talk matrix recovered from the best-fit calibration method (depicted in Fig. 4) under the approximation that $\Omega$ be unitary. This approximation ignores scattering into higher-order modes beyond the span of the truncated mode set (see Appendix E for details). The first row of the amplitude matrix corresponds to leakage from the fundamental mode into higher-order modes (imperfect star nulling), while the first column of the amplitude matrix corresponds to leakage from the higher-order modes into the fundamental mode (lost exoplanet light). The phase components of the experimental cross-talk matrix are artifacts of the Fourier–Zernike mode definitions and the unitary approximation of $\Omega$.

Download Full Size | PDF

The maximum star–planet contrast accessible with our benchtop coronagraph is set by the leakage of the fundamental mode into higher-order modes. This leakage causes light from the on-axis star to corrupt an otherwise pristine signal from the exoplanet. The throughput of the coronagraph is given by $P({\vec r_0}) = {{\textbf 1}^{\intercal}}{\boldsymbol q}({\vec r_0})$, which represents the probability that a photon emitted by a point source at position ${\vec r_0}$ arrives at the detector after propagating through the system. The leakage out of the fundamental mode is given by the throughput of the on-axis star ${P_s} = P(\vec 0)$. Equivalently, the throughput for an exoplanet located at ${\vec r_e}$ is given by ${P_e} = P({\vec r_e})$.

Under shot-noise-limited conditions, the SNR of the exoplanet signature is

 

Fig. 6. End-to-end throughput of the coronagraph as measured at the image plane (black scattered dots). The throughput predicted by the experimental model (solid black curve) makes use of the MPLC cross-talk matrix of Fig. 5. For reference, we also plot the ideal throughput (solid red curve) achieved by a perfect coronagraph that exclusively rejects the fundamental mode. The inset depicts the throughput on a log scale near the optical axis. When the source is on-axis, the measured throughput approximates the leakage from the fundamental mode predicted by the cross-talk matrix ${\sim}{10^{- 3}}$.

Download Full Size | PDF

The entries of the cross-talk matrix corresponding to imperfect isolation of the fundamental mode (i.e., the first row/column of $\Omega$) also introduce distortions to the final image of the exoplanet. To see this, first consider an ideal cross-talk matrix of the form

 

Fig. 7. Comparison between theoretical, experimental, and simulated direct imaging measurements of the artificial exoplanet for various sub-diffraction star–planet separations along the $y$-axis under a 1000:1 star–planet contrast. In each image, the ground truth location of the exoplanet coincides with the white cross-hairs, while the star is located at the origin. All images are max-normalized to illustrate dominant spatial features. THEORETICAL: These images correspond to the perfect coronagraph operating in the absence of background or modal cross-talk. EXPERIMENTAL: These are images of the artificial exoplanet collected with our benchtop coronagraph after subtraction of the background. SIMULATED: These images were generated by numerically propagating the optical field through a digital twin of the MPLC device used in our experimental setup. Noticeable distortions and asymmetries in the experimental and simulated images for equidistant off-axis ${\pm}{y_e}$ exoplanet locations arise from modal cross-talk matrices depicted in Figs. 5 and 14, respectively.

Download Full Size | PDF

Conversely, if the mode sorter does not perfectly isolate the fundamental mode [i.e., $\Omega$ does not assume the form of Eq. (6)], then interference between the residual light in the fundamental mode and orthogonal modes cause distortions to the image of the exoplanet. The precise nature of these distortions depends on the specific amplitudes and phases of each off-diagonal entry in the cross-talk matrix, given the target basis that the mode sorter is configured to (de)multiplex. In practice, Eq. (7) must be approached with progressively higher fidelity to gain access to more extreme star–planet contrasts and lower image distortion. In Appendix E, we perform a simulated analysis of the cross-talk matrix for the MPLC deployed in our experimental setup.

3. RESULTS

A. Exoplanet Localization

In Fig. 7, we compare theoretical, experimental, and simulated images of an artificial star–planet system captured with our coronagraph. For separations $|{y_e}| \gt 0.1\sigma$, there are strong qualitative similarities between the theoretical and experimental image intensity profiles, indicating that the exoplanet signal dominates over the background noise. However, at separations below this threshold, where the vast majority of exoplanet photons are discarded with the fundamental mode, the inherent noise in our experimental system becomes more prominent than the exoplanet signal. Additionally, the asymmetry observed between images of the exoplanet positioned at equal distances above and below the optical axis (${\pm}{y_e}$) emerges due to the cross-talk between modes ${{\boldsymbol \psi}_1}$ and ${{\boldsymbol \psi}_2}$.

For each exoplanet location ${y_e} \in {\mathscr Y}$, we compiled repeated synthetic measurements $\{{{\boldsymbol Y}^{(i)}}({y_e})\}$ for $i = 1, \ldots ,(\ell = 100)$ with an exposure time of $T = 100\;\text{ms} $. To localize the exoplanet, we employ a MLE,

We note that the MLE of the exoplanet position implicitly assumes a priori knowledge of the star–planet contrast $b$ via the distribution on ${\boldsymbol Y}$. In the idealized case where the cross-talk matrix takes the form of Eq. (6), such that all light from the star is eliminated, knowledge of the contrast is not required because the measurement only contains contributions from the exoplanet ${\boldsymbol Y} \propto {{\boldsymbol X}_e}$. In practice, any true implementation of a mode sorter will likely suffer from some amount of cross-talk that prohibits complete extinction of the fundamental mode. Consequently, the star–planet contrast becomes a nuisance parameter for exoplanet localization and detection tasks. This reality applies to any coronagraphy system with leakage out of the fundamental mode, not specifically our design. A possible solution is to jointly estimate the exoplanet position and the contrast, though this extension lies beyond the scope of the current work.

Figure 8(a) shows a map of the log-likelihood ${\mathscr L}({y_e};{y_e^\prime}) = \log P({\boldsymbol Y}({y_e})|{y_e^\prime})$ averaged over all experimental trials. The likelihood map exhibits a peak ridge line (maximum likelihood) that corresponds closely with the ground-truth position of the exoplanet. However, for exoplanets residing near the optical axis and outside the nominal region of support for the truncated mode basis, the likelihood exhibits weak curvature, indicating greater estimator uncertainty. Figure 8(b) shows the MLEs obtained for all repeated measurements over the exoplanet translation scan, reinforcing salient features of the log-likelihood map. The MLE accurately localizes the exoplanet within the nominal region of support $|{y_e}| \le .6\sigma$. Outside of this domain, the estimator experiences an increasing bias and uncertainty as photons are lost to modes beyond the support of the truncated mode set.

 

Fig. 8. (a) Map of the log-likelihood ${\mathcal L}(\check{y} _e;{y_e})$ over sub-diffraction exoplanet locations. The peak likelihood ridge line (dark red regions) nearly coincides with perfect localization (black diagonal). However, certain regions near the optical axis and outside the support of the truncated mode basis exhibit weakly peaked likelihoods, indicating a greater variance in the MLE. (b) Scatter plot of the MLEs obtained from the complete measurement ensemble collected during the exoplanet position scan. The color bar indicates the relative frequency of an estimate at each exoplanet position. Obvious outliers were removed for visual clarity.

Download Full Size | PDF

B. Statistical Performance Analysis

To further quantify the performance of our coronagraph, we analyze the statistical error and imprecision across the MLE realizations $\{\check{y} _e^{(i)}({y_e}):{y_e} \in {\mathscr Y}\}$. Given a ground truth exoplanet position ${y_e}$, we denote the mean and variance of the MLE to be ${\bar y_e}({y_e}) \equiv {\mathbb{E}_{{\boldsymbol Y}|{y_e}}}[{{\check{y}} _e}]$ and $\sigma _e^2({y_e}) = {\mathbb{V}_{{\boldsymbol Y}|{y_e}}}[{{\check{y}} _e}]$, respectively. The unbiased empirical estimators of the MLE’s mean and variance are given by the statistical averages,

 

Fig. 9. (a) Average error of the exoplanet MLE as a function of the ground truth position (black curve). Green and blue shaded regions demarcate where the absolute localization error is below $\sigma /50$ and $\sigma /25$, respectively. (b) Statistical standard deviation of the MLE as a function of the exoplanet position as compared to the Cramér–Rao lower bound (CRLB) for the measurement model of Eq. (3). Uncertainty bars on the imprecision estimator were calculated using jackknife resampling [38]. The red shaded region indicates viable bounds of the CRLB given uncertainty in the experimental model parameters. In particular, we show the upper and lower bounds of the CRLB associated with the left and right sides of the FWHM, respectively, for ${\lambda _0}$ in Fig. 12 (Appendix D).

Download Full Size | PDF

The variance of the MLE is subject to the Cramér–Rao lower bound (CRLB) for a biased estimator,

4. CONCLUSION

In summary, we present an experimental quantum-optimal direct imaging coronagraph design using spatial mode (de)multiplexing, which has inherent advantages for exoplanet exploration relative to indirect observation techniques. Furthermore, we demonstrate high-accuracy localization of an artificial exoplanet at sub-diffraction distances from its host star, which provides new capability to search and study exoplanets that may be inaccessible (i.e., in terms of inner working angle) to current state-of-the-art coronagraphs.

The performance of our experimental system at deeply sub-diffraction separations is presently limited by a combination of modal cross-talk, detector dark noise, and ambient background light. These experimental non-idealities can be mitigated comprehensively with higher-fidelity detectors, mode-sorters, and further optical design optimizations, enabling operation at deeper sub-diffraction separations and higher contrast levels. For large star–planet separations, our system performance is limited by the truncation of the mode basis, which may be improved by extending the number of phase planes comprising the MPLC. Furthermore, recent techniques for in situ optimization of MPLC phase masks [39] could overcome model mismatch between our experimental system and its digital counterpart, enabling faithful numerical performance analyses. Additionally, while we employ a synthetic measurement construction in this work, follow-up studies involving two simultaneously radiating sources would provide more robust validation of our coronagraph design. Overall, we believe our work provides a compelling basis for further exploration and development of mode-sorting solutions for high-contrast astronomical imaging tasks.

Looking forward, even-order coronagraphs may be implemented using spatial mode filtering techniques similar to those presented here to accommodate the finite size of the host star. Previous research have shown that shot noise induced by an extended host star may be progressively tempered by judiciously filtering/attenuating higher-order optical modes [22]. Furthermore, extending functionality to broadband sources will be critical for analyzing the spectrum of the exoplanet, which embeds key planetary science information. Broadband sources introduce mode mismatch, such that the cross-talk matrix becomes wavelength-dependent. For sub-diffraction exoplanets, the constraints on acceptable levels of mode mismatch become more stringent [40,41], ostensibly limiting the efficacy of spatial mode sorting methods. However, recent numerical [42] and experimental [43] work suggests that MPLCs may be used to simultaneously sort spectral and spatial modes, providing a common platform for exoplanet spectroscopy and/or nulling of a broadband star.

APPENDIX A: FOURIER–ZERNIKE MODES

For an imaging system with a circular pupil of radius $R$, focal length $f$, and operating wavelength $\lambda$, let $({X_a},{Y_a})$ and $({X_b},{Y_b})$ denote the coordinate space of the pupil plane and focal plane, respectively. We define the dimensionless pupil plane and focal plane coordinate vectors as

(A1a)$$\vec u \equiv \frac{1}{R}({X_a},{Y_a}),$$

(A1b)$$\vec r \equiv \frac{R}{{\lambda f}}({X_b},{Y_b}).$$

The Zernike modes ${{\boldsymbol {\tilde \psi}}_{\textit{nm}}}(\vec u)$ constitute a PSF-matched basis over a circular pupil. In polar coordinates, they are given by

(A2a)$${{\boldsymbol {\tilde \psi}}_{\textit{nm}}}(u,\theta) \equiv {R_{\textit{nm}}}(u){\Theta _m}(\theta)\text{circ}(u),$$

(A2b)$${R_{\textit{nm}}}(u) \equiv \sum\limits_{j = 0}^{(n - |m|)/2} \frac{{{{(- 1)}^j}\sqrt {n + 1} (n - j)!}}{{j![(n + m)/2 - j]![(n - m)/2 - j]!}}{u^{n - 2j}},$$

(A2c)$${\Theta _m}(\theta) \equiv \left\{{\begin{array}{*{20}{l}}{\sqrt 2 \cos (|m|\theta)}&{(m \gt 0)}\\1&{(m = 0)}\\{\sqrt 2 \sin (|m|\theta)}&{(m \lt 0)}\end{array}}, \right.$$

(A2d)$$\text{circ}(u) \equiv \left\{{\begin{array}{*{20}{c}}{1,}&{u \le 1}\\{0,}&{u \gt 1}\end{array}} ,\right.$$

where the radial index range is $n = 0,1,2, \ldots ,\infty$. The angular index range for a given radial index is $m \in {S_n} = \{- n, - n + 2, \ldots ,n - 2,n\}$. These modes are defined to satisfy the orthonormality condition,

$$\int {\boldsymbol {\tilde \psi}}_{\textit{nm}}^*(u,\theta){{\boldsymbol {\tilde \psi}}_{n^\prime m^\prime}}(u,\theta)u\text{d}u\text{d}\theta = {\delta _{{nn^\prime}}}{\delta _{{mm^\prime}}}.$$

The Fourier transform of the Zernike modes over the pupil is found in [44] to be

(A3)$${{\boldsymbol \psi}_{\textit{nm}}}(r,\phi) = {i^{n + 2|m|}}\sqrt {n + 1} \frac{{{J_{n + 1}}(2\pi r)}}{{\sqrt \pi r}}{\Theta _m}(\phi),$$

which we refer to as the “Fourier–Zernike” modes throughout the main text. In this work, we sort the truncated mode basis $\{{{\boldsymbol \psi}_0},{{\boldsymbol \psi}_1},{{\boldsymbol \psi}_2},{{\boldsymbol \psi}_3}\}$ corresponding to modes $\{{{\boldsymbol \psi}_{00}},{{\boldsymbol \psi}_{1, - 1}},{{\boldsymbol \psi}_{2,0}},{{\boldsymbol \psi}_{2,2}}\}$. The amplitude of these modes is shown in Fig. 2(b).

APPENDIX B: MEASUREMENT MODEL

Here, we expand on the single-shot measurement model of Eq. (3). The post-nulling optical intensity distribution on the detector is given in Eq. (2) to be

$${\boldsymbol q}({\vec r_0}) = |\Psi \Omega C{\Omega ^\dagger}{\Psi ^\dagger}{\boldsymbol u}({\vec r_0}{)|^2},$$

which provides a matrix description of forward propagating the optical field through the MPLC, nulling the fundamental mode, back-propagating the field to an image plane, and measuring the intensity on a detector array. Here, ${\boldsymbol u}({\vec r_0}) = \text{vec}[{\boldsymbol \psi}(\vec r - {\vec r_0})]$ is the vectorized form of the shifted PSF satisfying ${\boldsymbol u}{({\vec r_0})^\dagger}{\boldsymbol u}({\vec r_0}) = 1$. $\Psi$ is the truncated change-of-basis matrix for which the columns are the vectorized Fourier–Zernike modes ${\psi _k} = \text{vec}[{{\boldsymbol \psi}_k}(\vec r)]$,

(B1)$$\Psi = \left[{\begin{array}{*{20}{c}}|&{}&|\\{{\psi _0}}& \cdots &{{\psi _{K - 1}}}\\|&{}&|\end{array}} \right],$$

satisfying the orthonormality condition ${\Psi ^\dagger}\Psi = {I_K}$, where ${I_K}$ is the identity matrix on ${\mathbb{C}^K}$. The coronagraph matrix is

(B2)$$C = \pm \left[{\begin{array}{*{20}{c}}0&\quad {{{\textbf 0}^ \intercal}}\\{\textbf 0}&\quad {{I_{K - 1}}}\end{array}} \right],$$

which serves to null the port assigned to the fundamental mode. The sign of the coronagraph matrix depends on whether the system is transmissive $(+)$ at the sorting plane, as shown in Fig. 2(a), or reflective $(-)$ at the sorting plane, as in our experimental configuration depicted in Fig. 2(b). While we treat the optical field as a scalar quantity (as the field is always in a definite linear polarization state throughout our experiment), details on the experimental polarization manipulation for path splitting are illustrated in Fig. 10.

 

Fig. 10. Unfolded beam path illustrating polarization control for path splitting of forward- (black arrows) and backward- (red arrows) propagating beams. Through this manipulation, we implement the design of Fig. 2 using a single-mode sorter.

Download Full Size | PDF

APPENDIX C: FUNCTION VECTORIZATION

Let $\Delta x,\Delta y$ be the horizontal and vertical pixel pitch dimensions for a specific camera, and let $\{({x_i},{y_i})\} _{i = 1}^M$ be the set of points corresponding to the center of each pixel. The $\text{vec}(\cdot)$ operation takes a well-behaved function $f(x,y)$ over the plane and converts it to a discrete column vector ${\textbf f}$, whose entries are given by

(C1)$${{ f}_i} = \int_{{x_i} - \Delta x/2}^{{x_i} + \Delta x/2} \int_{{y_i} - \Delta y/2}^{{y_i} + \Delta y/2} f(x,y)\text{d}x\text{d}y.$$

APPENDIX D: CALIBRATION

This section addresses our procedures for characterizing all system-level parameters $\{\Omega ,{\lambda _D},{\lambda _0},{\boldsymbol s}\}$ taken to be known quantities in the measurement model.

A. Cross-Talk Matrix $\Omega$

Figure 4 shows calibration measurements made at the sorting plane of the MPLC for characterizing modal cross-talk. For each position ${y_i} \in \mathscr Y$ of the point source, where $i = 1, \ldots ,(P = |\mathscr Y|)$, we measure the intensity in each mode channel ${{\textbf v}_i} \in {\mathbb{R}^K}$. We estimate the cross-talk matrix from these calibration measurements by optimizing a least-squares objective function $\mathscr{F}(\Omega)$ under the constraint that $\Omega$ be unitary,

(D1)$$\check{\Omega} = \mathop {\text{argmin}}\limits_{\Omega : \Omega \in U(K)} \underbrace {||V - \tilde V(\Omega)||_F^2}_{\mathscr{F}(\Omega)},$$

where $|| \cdot {||_F}$ is the Frobenius-norm, and $\odot$ is the Hadamard element-wise product. In formulating this optimization problem, we define the following terms:

• • $V \in {\mathbb{R}^{K \times P}}$: measurements of intensity in each mode channel. The ${i^{\text{th}}}$ column of $V$ is the measurement vector ${{\boldsymbol v}_i}$ collected for the ${i^{\text{th}}}$ position of the point source.
• • ${\tilde V}(\Omega) = ({\Omega ^\dagger}Z) \odot {({\Omega ^\dagger}Z)^*} \in {\mathbb{R}^{K \times P}}$: forward model of the expected intensities measured in each mode channel under the unitary cross-talk matrix $\Omega$.
• • $Z \in {\mathbb{R}^{K \times P}}$: mode expansion coefficients for each off-axis source location in the calibration scan. The ${i^{\text{th}}}$ column of $Z$ is the vector ${{\boldsymbol z}_i} = {\Psi ^\dagger}{\boldsymbol u}({y_i}) \in {\mathbb{C}^K}$ for ${y_i} \in \mathscr Y$.
• • ${\Gamma _\Omega} = (\frac{{\partial \mathscr{F}}}{{\partial {\Omega ^*}}})(\Omega) \equiv \frac{1}{2}\big[\frac{{\partial \mathscr{F}}}{{\partial {\Omega _R}}} + i\frac{{\partial\mathscr{F}}}{{\partial {\Omega _I}}}\big]$: gradient of the objective function with respect to the cross-talk matrix in the Euclidean space ${\mathbb{R}^{2K \times 2K}}$, where ${\Omega _R} = \Re \Omega$ and ${\Omega _I} = \Im \Omega$. The explicit Euclidean gradient under the squared error loss function is

  (D2)$$\left(\frac{{\partial \mathscr F}}{{\partial {\Omega ^*}}}\right)(\Omega) = - 2Z {[(V - \tilde V(\Omega)) \odot ({\Omega ^\dagger}Z)]^\dagger}.$$
• • $G(\Omega) = {\Gamma _\Omega}{\Omega ^\dagger} - \Omega \Gamma _\Omega ^\dagger$: Riemannian gradient of the objective function over the Lie group $U(K)$.

We solve this optimization using the constrained descent algorithm detailed in [45] and [46]. The algorithm is briefly summarized as follows. Instantiate ${\Omega _0} = I$ with a learning rate $\mu \gt 0$. Then, at iteration $t$, we update the cross-talk matrix as

(D3)$${\Omega _{t + 1}} = \textsf{expm}\,(-\mu{G_t}){\Omega _t},$$

where ${G_t} = G({\Omega _t})$, and $\textsf{expm}(\cdot)$ is a numerical matrix exponentiation function. The iterations are performed until convergence. In Fig. 4(a), we show the measured mode intensity data alongside the theoretical mode intensities after applying the best-fit cross-talk matrix shown in Fig. 5.

B. Dark Noise Rate ${\lambda _D}$ and Structured Background ${\boldsymbol s}$

The background contribution term employed in our single-shot measurement model arises from a combination of undesired stray light reflections and the dark noise rate of the pixels in the camera. We thus decompose the background into two terms as

(D4)$${\lambda _B}{{\boldsymbol p}_B} = {\boldsymbol s} + {\lambda _D}{\textbf 1},$$

where ${\boldsymbol s} \in {\mathbb{R}^M}$ is a structured background rate observed across our experimental measurements, and ${\lambda _D} \in \mathbb{R}$ is the spatially uniform dark click rate of each pixel operating at room temperature. Furthermore, ${\textbf 1} \in {\mathbb{R}^M}$ is a vector of ones, and we have imposed the constraint ${{\textbf 1}^\top}{{\boldsymbol p}_B} = 1$, such that ${{\boldsymbol p}_B}$ is a discrete probability distribution. Figure 11 shows calibration data for estimating the dark noise rate ${\lambda _D}$ and the structured background rate ${\boldsymbol s}$ of Eq. (D4). An estimate of the dark noise rate $\check{\lambda}_D = 246$ [photons/$T$] was found by fitting a Poisson distribution to a histogram of photon counts recorded by our camera operating in a dark environment at room temperature. The structured background rate was estimated by averaging $n={ 10^3}$ images of the on-axis source as viewed through the coronagraph,

(D5)$$\check{s} = \frac{1}{n}\sum\limits_{i = 1}^n [{{\boldsymbol X}^{(i)}}(\vec 0) - {\lambda _D}{\textbf 1}],$$

where we have assumed that the photons arriving at the detector due to imperfect nulling of the fundamental mode (i.e., modal cross-talk) are negligible compared to the dark noise rate ${\lambda _0}{\boldsymbol q}(\vec 0) \ll {\lambda _D}$.

 

Fig. 11. (a) Poisson fit to the detector dark-noise histogram. The dark rate of the detector was found to be ${\lambda _D} = 246$ [photons $\cdot {T^{- 1}}$] per pixel. (b) Structured background profile ${\boldsymbol s}$ induced by stray light and undesired reflections in the experimental setup.

Download Full Size | PDF

C. Signal Photon Rate ${\lambda _0}$

The mean photon rate entering the pupil from the star–planet scene is approximated using the best linear unbiased estimator (BLUE) [47]. Let the random variable $Y_j^{(i)}({\vec r_e})$ represent the photons collected in the ${j^{\text{th}}}$ pixel of the ${i^{\text{th}}}$ measurement for a given exoplanet location ${\vec r_e}$. In accordance with Eq. (3), each pixel follows a distribution:

$$Y_j^{(i)}({\vec r_e}) \sim \text{Poiss}({\Lambda _0}{p_j}({\vec r_e}) + {\Lambda _B}{p_{{B_j}}}),$$

which we approximate as a non-random term linear in ${\Lambda _0}$, with zero-mean additive Gaussian noise,

(D6)$$Y_j^{(i)}({\vec r_e}) = {\Lambda _0}{p_j}({\vec r_e}) + {\Lambda _B}{p_{{B_j}}} + {\varepsilon _j},$$

(D7)$${\varepsilon _j}({\Lambda _0},{\vec r_e}) \sim {\cal N}(0,{\Lambda _0}{p_j}({\vec r_e}) + {\Lambda _B}{p_{{B_j}}}).$$

Defining the residuals,

(D8)$$r_j^{(i)}({\Lambda _0},{\vec r_e}) = Y_j^{(i)}({\vec r_e}) - {\Lambda _0}{p_j}({\vec r_e}) - {\Lambda _B}{p_{{B_j}}},$$

the BLUE is found by optimizing the sum of squared residuals weighted by the variance of the photon arrivals in each pixel,

(D9)$${\check{\Lambda}_0} = \mathop {\text{argmax}}\limits_{{\Lambda _0}} \sum\limits_{{y_e} \in {\mathscr Y}} \sum\limits_{i = 1}^\ell \sum\limits_{j = 1}^M \frac{{{{(r_j^{(i)}({\Lambda _0},{y_e}))}^2}}}{{\mathbb{V}[{\varepsilon _j}({\Lambda _0},{y_e})]}}.$$

We numerically solve this optimization problem over the entire collected data set and estimate the mean photon rate entering the pupil from a single point source to be $\check{\lambda}_0=\check{\Lambda}_0/N=3.17\times10^{6}$ [photons/T].

D. Magnification and Alignment Correction

Figure 13 shows the experimental coronagraph throughput as a function of the scanned source position. Comparing the measured throughput to the theoretical throughput curve of our system revealed a scaling $m$ and shift $a$ between the true source displacements ${y_e}$ relative to the diffraction limit and the displacements computed from our system specifications $y_e^\prime = m({y_e} - a)$. The scaling and shift terms arise from experimental errors in magnification and alignment of the optical axis, respectively. After accounting for these discrepancies, we see that the throughput of the curve of our transformed (scaled and shifted) coordinates align with theoretical predictions.

 

Fig. 12. Histogram over estimates of the scene photon emission rate ${\lambda _0}$ from which we compute the best linear unbiased estimator (dashed red) $\check{\lambda}_{0} = 3.17 \times {10^6}$ [Photons $\cdot {T^{- 1}}$]. The edges of the FWHM are used in Fig. 9 to determine the range of the Cramér–Rao Bounds associated with our experiment.

Download Full Size | PDF

 

Fig. 13. Alignment and magnification discrepancies in the $4f$ system are accounted for by fitting the experimental throughput to the theoretically predicted values.

Download Full Size | PDF

 

Fig. 14. (a) Amplitude and (b) phase of the complex cross-talk matrix calculated from numerical simulations. (c) Simulated target and output fields of the mode sorter at the sorting plane after propagating the Zernike modes through the MPLC. The entries of the cross-talk matrix are the pair-wise projections (overlap integrals) between target modes and output modes.

Download Full Size | PDF

APPENDIX E: SIMULATED MPLC CROSS-TALK MATRIX

In this section, we address the complex cross-talk matrix of the MPLC mode sorter in greater detail. Historically, literature in spatial mode sorting has placed emphasis on defining the cross-talk matrix in terms of the power leakage between modes (amplitude-only characterization). However, a complete characterization of the cross-talk requires the estimation of both the amplitude and phase. Direct experimental characterization of the complex cross-talk matrix is critical to forthcoming high-performance mode-sorting applications, such as coronagraphy and spatial mode filtering. In this section, we study the cross-talk matrix of the MPLC used in our experimental design with numerical field propagation simulations.

The goal of an MPLC is to passively map a given set of input modes $\{{\psi _i} \in {\mathbb{C}^M}\} _{i = 1}^K$ into a desired set of output modes $\{{\phi _i} \in {\mathbb{C}^M}\} _{i = 1}^K$. For concreteness, let $\Psi ,\Phi \in {\mathbb{C}^{M \times K}}$ be matrices whose columns are (vectorized) orthonormal input and output modes represented in the pixel basis. In general, the number of pixel degrees of freedom $M$ is much greater than the number of modes $K$ that one wishes to sort. Additionally, we take the modes to be orthonormal, such that ${\Psi ^\dagger}\Psi = {\Phi ^\dagger}\Phi = {I_K}$. The cross-talk matrix $\Omega \in {\mathbb{C}^{K \times K}}$ describes the interferometric connection between $\Psi$ and $\Phi$ within the MPLC transfer matrix,

$$T = \Phi {\Omega ^\dagger}{\Psi ^\dagger}.$$

Under this decomposition of the transfer matrix $T$, the action of the MPLC on an arbitrary incident field can be understood as follows: first, project the arriving field onto the basis of input modes via ${\Psi ^\dagger}$. Then, interfere light in the input mode channels according to the cross-talk matrix $\Omega$. Finally, expand the mixed channels back into the pixel basis via $\Phi$. Note that if $\Omega = {I_K}$, then the MPLC perfectly maps the input modes to the output modes in a bijective fashion. In the most general case, however, $\Omega$ admits a singular value decomposition, which comprehensively handles interferometric mixing between mode channels as well as loss through the device.

A simple way to calculate the elements of the complex cross-talk matrix in simulation is by forward propagating each input mode independently and taking the projection between the resulting output field at the sorting plane and the target output mode. That is, we perform the computation

$${\Omega _{\textit{ij}}} = (\phi _i^\dagger T{\psi _j}{)^*},\quad T = {H_\tau}{D_\tau} \cdots {H_1}{D_1},$$

where we compute the transfer matrix via alternating the application of phase masks (diagonal matrices ${D_i}$) inter-spaced by plane-to-plane light propagation (unitary matrices ${H_i}$). Figures 14(a) and 12(b) show the simulated cross-talk of our MPLC device. We note that this simulated cross-talk matrix is not unitary as the forward-propagated modes span a slightly different subspace compared to the target mode set. This is visually evident in Fig. 14(c), where residual speckles in the output fields are observed to lie outside the designated regions of support for the target modes. Contrast this matrix to the experimental cross-talk matrix displayed in Fig. 5, which we evaluated using the methods described in Appendix D(A). Importantly, for calibrating the experimental cross-talk matrix, we assumed loss was negligible and forced the requirement that $\Omega$ be unitary. Despite this reduction in the degrees of freedom, the unitary cross-talk matrix calibrated for our benchtop implementation of the coronagraph provides an experimental model that aligns reasonably well with our measured results.

APPENDIX F: NUMERICAL COMPARISONS TO OTHER CORONAGRAPHS

As points of comparison, we numerically evaluate the performance of the experimental coronagraph against the PIAACMC, the vortex (charge-2) coronagraph, and the perfect coronagraph. The perfect coronagraph is the idealized version of our experimental implementation in the absence of mode truncation or cross-talk. We assign the same background illumination profile espoused in our experiment ${\Lambda _B}{{\boldsymbol p}_B}$ to the measurement model of each coronagraph. Then, the distinguishing factor between coronagraphs manifests in the term ${\Lambda _0}{\boldsymbol p}({\vec r_e})$, which now depends on the singular-value decomposition of the field transformation performed by a given coronagraph,

(F1)$${\boldsymbol q}({\vec r_0}) = |U\Sigma {V^\dagger}{\boldsymbol u}({\vec r_0}{)|^2},$$

where $U,V \in {\mathbb{C}^{M \times M}}$ are unitary matrices, and $\Sigma \in {\mathbb{R}^{M \times M}}$ is a diagonal matrix of singular values. Figure 15 shows the Cramér–Rao lower bounds associated with each coronagraph. We see that the idealized version of our experimental setup (Experiment 4-Mode) outperforms both the PIAACMC and vortex coronagraph for over sub-diffraction exoplanet locations in the range $|{y_e}| \lt = 0.4\sigma$ that corresponds to where the CFI of the truncated mode basis begins to depart from the QFI limit, as shown in Fig. 3(b).

 

Fig. 15. Comparison of the CRLBs for unbiased estimators of the exoplanet position under measurements with different coronagraphs. Each coronagraph measurement model includes the same empirical background noise characterized in our experimental system.

Download Full Size | PDF

APPENDIX G: CRAMÉR–RAO LOWER BOUND FOR THE BIASED MLE

The maximum likelihood estimator is known to be an unbiased estimator in the asymptotic limit of many repeated measurements (here, the term “measurement” formally refers to the detection of a photon). While this assumption is increasingly valid for larger star–planet separations, wherein the number of detected exoplanet photons is large, it becomes invalid for small star–planet separations, as the vast majority of exoplanet photons are discarded by the coronagraph. This motivates the analysis of Eq. (11) in the main text. The expectations ${\mathbb{E}_{{\boldsymbol Y}|{y_e}}}[\cdot]$ are taken over the joint distribution on the measurement outcome ${\boldsymbol Y}$ given by

(G1a)$$P({\boldsymbol Y}|{y_e}) = \prod\limits_{j = 1}^M P({Y_j}|{y_e}),$$

(G1b)$$P({Y_j} = n|{y_e}) = \frac{{{e^{- {\tau _j}({y_e})}}{\tau _j}{{({y_e})}^n}}}{{n!}},$$

(G1c)$${\tau _j}({y_e}) = {\Lambda _0}{p_j}({y_e}) + {\Lambda _B}{p_{{B_j}}}.$$

Finally, $I({y_e})$ is the total classical Fisher information defined as

(G2)$$I({y_e}) = \sum\limits_{\mathscr{D}} P({\boldsymbol Y}|{y_e})\left[\frac{d}{{d{y_e}}}\log P({\boldsymbol Y}|{y_e})\right]^2,$$

where $\mathscr D$ is the event space of possible measurement outcomes for ${\boldsymbol Y}$. In Fig. 16(a), we plot the bias of the MLE, which we numerically computed by Monte–Carlo sampling the probabilistic measurement model of Eq. (3) to determine the expectation of the MLE ${\mathbb{E}_{{\boldsymbol Y}|{y_e}}}[\check{y}_{e}({y_e})]$. We see that near the optical axis, the MLE is clearly biased, as it no longer operates in the asymptotic regime. In Fig. 16(b), we show the corresponding biased CRLB for the MLE of the exoplanet position.

 

Fig. 16. Numerical evaluation of the MLE bias and the associated CRLB curve under the experimental model. This CRLB curve is equivalent to that used in Fig. 9(b).

Download Full Size | PDF

Funding

National Science Foundation Graduate Research Fellowship Program (DGE-2137419).

Acknowledgment

ND acknowledges support from the National Science Foundation Graduate Research Fellowship. IO and SG acknowledge that this research was supported by Raytheon and recognize the contributions of Jaime Bucay and Mark Meisner for their insights.

Disclosures

The authors declare no conflicts of interest.

Data availability

Our codebase for figure generation, coronagraph simulation, and measurement processing conducted in this work can be found in the project GitHub repository [48]. Raw measurement data files are hosted on Zenodo and are free to download at [49].

REFERENCES

1. J. T. Wright and B. S. Gaudi, Exoplanet Detection Methods (Springer, 2012).

2. R. L. Akeson, X. Chen, D. Ciardi, et al., “The NASA exoplanet archive: data and tools for exoplanet research,” Publ. Astron. Soc. Pac. 125, 989 (2013). [CrossRef]  

3. J. Schneider, C. Dedieu, P. Le Sidaner, et al., “Defining and cataloging exoplanets: the exoplanet.eu database,” Astron. Astrophys. 532, A79 (2011). [CrossRef]  

4. T. Currie, B. Biller, A.-M. Lagrange, et al., “Direct imaging and spectroscopy of extrasolar planets,” arXiv, (2023). [CrossRef]  

5. W. A. Traub and B. R. Oppenheimer, Direct Imaging of Exoplanets (University of Arizona, 2010).

6. S. Seager and D. Deming, “Exoplanet atmospheres,” Annu. Rev. Astron. Astrophys. 48, 631–672 (2010). [CrossRef]  

7. B. A. Biller and M. Bonnefoy, Exoplanet Atmosphere Measurements from Direct Imaging (Springer International Publishing, 2018), pp. 2107–2135.

8. B. Biller, “Detecting and characterizing exoplanets with direct imaging: past, present, and future,” Proc. Int. Astron. Union 8, 1–11 (2013). [CrossRef]  

9. G. Foo, D. M. Palacios, and G. A. Swartzlander, “Optical vortex coronagraph,” Opt. Lett. 30, 3308–3310 (2005). [CrossRef]  

10. E. Mari, F. Tamburini, G. A. Swartzlander, et al., “Sub-Rayleigh optical vortex coronagraphy,” Opt. Express 20, 2445–2451 (2012). [CrossRef]  

11. O. Guyon, F. Martinache, R. Belikov, et al., “High performance PIAA coronagraphy with complex amplitude focal plane masks,” Astrophys. J. Suppl. Ser. 190, 220 (2010). [CrossRef]  

12. C. Aime, R. Soummer, and A. Ferrari, “Total coronagraphic extinction of rectangular apertures using linear prolate apodizations,” Astron. Astrophys. 389, 334–344 (2002). [CrossRef]  

13. R. Soummer, C. Aimé, and P. Falloon, “Stellar coronagraphy with prolate apodized circular apertures,” Astron. Astrophys. 397, 1161–1172 (2003). [CrossRef]  

14. R. Soummer, “Apodized pupil LYOT coronagraphs for arbitrary telescope apertures,” Astrophys. J. 618, L161 (2004). [CrossRef]  

15. A. Zurlo, Direct Imaging of Exoplanets (Springer, 2024).

16. O. Guyon, E. A. Pluzhnik, M. J. Kuchner, et al., “Theoretical limits on extrasolar terrestrial planet detection with coronagraphs,” Astrophys. J. Suppl. Ser. 167, 81 (2006). [CrossRef]  

17. M. Tsang, R. Nair, and X.-M. Lu, “Quantum theory of superresolution for two incoherent optical point sources,” Phys. Rev. X 6, 031033 (2016). [CrossRef]  

18. M. Tsang, “Resolving starlight: a quantum perspective,” Contemp. Phys. 60, 279–298 (2019). [CrossRef]  

19. K. K. Lee, C. N. Gagatsos, S. Guha, et al., “Quantum-inspired multi-parameter adaptive Bayesian estimation for sensing and imaging,” IEEE J. Sel. Top. Signal Process. 17, 491–501 (2023). [CrossRef]  

20. S. Prasad, “Quantum limited super-resolution of an unequal-brightness source pair in three dimensions,” Phys. Scripta 95, 054004 (2020). [CrossRef]  

21. N. Deshler, S. Haffert, and A. Ashok, “Achieving quantum limits of exoplanet detection and localization,” arXiv, (2024). [CrossRef]  

22. R. Belikov, D. Sirbu, J. B. Jewell, et al., “Theoretical performance limits for coronagraphs on obstructed and unobstructed apertures: how much can current designs be improved?” Proc. SPIE 11823, 118230W (2021). [CrossRef]  

23. A.-M. Lagrange, F. Philipot, P. Rubini, et al., “Radial distribution of giant exoplanets at solar system scales,” Astron. Astrophys. 677, A71 (2023). [CrossRef]  

24. R. B. Fernandes, G. D. Mulders, I. Pascucci, et al., “Hints for a turnover at the snow line in the giant planet occurrence rate,” Astrophys. J. 874, 81 (2019). [CrossRef]  

25. J. Chen and D. Kipping, “Probabilistic forecasting of the masses and radii of other worlds,” Astrophys. J. 834, 17 (2016). [CrossRef]  

26. B. Ning, A. Wolfgang, and S. Ghosh, “Predicting exoplanet masses and radii: a nonparametric approach,” Astrophys. J. 869, 5 (2018). [CrossRef]  

27. Y. Xin, N. Jovanovic, G. Ruane, et al., “Efficient detection and characterization of exoplanets within the diffraction limit: nulling with a mode-selective photonic lantern,” Astrophys. J. 938, 140 (2022). [CrossRef]  

28. Y. Xin, D. Echeverri, N. Jovanovic, et al., “Laboratory characterization of a mode-selective photonic lantern for exoplanet characterization,” Proc. SPIE 12680, 126800I (2023). [CrossRef]  

29. Y. J. Kim, M. P. Fitzgerald, J. Lin, et al., “Coherent imaging with photonic lanterns,” Astrophys. J. 964, 113 (2024). [CrossRef]  

30. M. Tsang, “Quantum limit to subdiffraction incoherent optical imaging,” Phys. Rev. A 99, 012305 (2019). [CrossRef]  

31. M. Tsang, “Quantum limit to subdiffraction incoherent optical imaging. II. A parametric-submodel approach,” Phys. Rev. A 104, 052411 (2021). [CrossRef]  

32. J. Carpenter, N. K. Fontaine, B. R. M. Norris, et al., “Spatial mode sorter coronagraphs,” in Conference on Lasers and Electro-Optics Pacific Rim (CLEO-PR) (Optica Publishing Group, 2020).

33. N. K. Fontaine, R. Ryf, H. Chen, et al., “Laguerre-Gaussian mode sorter,” Nat. Commun. 10, 1865 (2019). [CrossRef]  

34. G. Labroille, B. Denolle, P. Jian, et al., “Efficient and mode selective spatial mode multiplexer based on multi-plane light conversion,” Opt. Express 22, 15599–15607 (2014). [CrossRef]  

35. J. R. Řehaček, Z. Hradil, B. Stoklasa, et al., “Multiparameter quantum metrology of incoherent point sources: towards realistic superresolution,” Phys. Rev. A 96, 062107 (2017). [CrossRef]  

36. I. Ozer, M. R. Grace, and S. Guha, “Reconfigurable spatial-mode sorter for super-resolution imaging,” in Conference on Lasers and Electro-Optics (CLEO) (IEEE, 2022), pp. 1–2.

37. I. Ozer, M. R. Grace, P.-A. Blanche, et al., “Adaptive super-resolution imaging without prior knowledge using a programmable spatial-mode sorter,” arXiv, (2024). [CrossRef]  

38. C.-F. J. Wu, “Jackknife, bootstrap and other resampling methods in regression analysis,” Ann. Statist. 14, 1261–1295 (1986). [CrossRef]  

39. J. C. A. Rocha, U. G. Būtaitė, J. Carpenter, et al., “Self-configuring high-speed multi-plane light conversion,” arXiv, (2025). [CrossRef]  

40. S. Prasad, “Quantum limited source localization and pair superresolution in two dimensions under finite-emission bandwidth,” Phys. Rev. A 102, 033726 (2020). [CrossRef]  

41. M. R. Grace, “Quantum limits and optimal receivers for passive sub-diffraction imaging,” Ph.D. thesis (The University of Arizona, 2022).

42. Y. Zhang, H. Wen, A. Fardoost, et al., “Simultaneous sorting of wavelengths and spatial modes using multi-plane light conversion,” arXiv, (2020). [CrossRef]  

43. M. Mounaix, N. K. Fontaine, D. T. Neilson, et al., “Time reversed optical waves by arbitrary vector spatiotemporal field generation,” Nat. Commun. 11, 5813 (2020). [CrossRef]  

44. G. M. Dai, “Zernike aberration coefficients transformed to and from Fourier series coefficients for wavefront representation,” Opt. Lett. 31, 501–503 (2006). [CrossRef]  

45. T. E. Abrudan, J. Eriksson, and V. Koivunen, “Steepest descent algorithms for optimization under unitary matrix constraint,” IEEE Trans. Signal Process. 56, 1134–1147 (2008). [CrossRef]  

46. T. Abrudan, J. Eriksson, and V. Koivunen, “Optimization under unitary matrix constraint using approximate matrix exponential,” in 39th Asilomar Conference on Signals, Systems and Computers (2005), pp. 242–246.

47. A. C. Aitken, “On least squares and linear combination of observations,” Proc. R. Soc. Edinburgh 55, 42–48 (1936). [CrossRef]  

48. N. Deshler, “Perfect-coronagraph-experiment,” GitHub, 2025, https://github.com/I2SL/Perfect-Coronagraph-Experiment.

49. N. Deshler, “Experimental demonstration of a quantum-optimal coronagraph using spatial mode sorters: Zenodo repository,” Zenodo, 2025, https://zenodo.org/records/15028674.