Off-axis Aberrations Improve the Resolution Limits of Incoherent Imaging

In a shift-invariant imaging system, the optical field at the image plane may be obtained from the optical field at the object plane via a convolution with the system’s PSF, $\psi$ [1]. In the presence of off-axis aberrations, however, the relation between object field, $U_{o}$, and image field, $U_{i}$, is given by a more general integral transformation. For simplicity, we assume one transverse spatial coordinate for which the aforementioned transformation is given by

$\displaystyle U_{i}(x)=\int_{-\infty}^{\infty}U_{o}(\xi)\psi(x,\xi)\,{\rm d}\xi,$

(1)

where $\xi$ and $x$ are the object and image plane spatial coordinates, respectively. The departure from a shift-invariant imaging system is capture by $\psi(x,\xi)$, which indicates that the system PSF may have a dependence on both $\xi$ and $x$ that is not confined to their difference $x-\xi$. The form for $\psi(x,\xi)$ may be obtained by considering the propagation of a quasi-monochromatic point source with wavenumber $k=2\pi/\lambda$, modeled as a Dirac delta impulse, located at $\xi$, through the imaging system. This impulse is mapped onto a linear phase factor at the pupil plane. There, the linear phase factor is modulated by a Gaussian pupil function of characteristic width $\sigma_{p}$ and a phase aberration function, $\Delta W$. The choice of a Gaussian pupil function, commonly made in related works, is used for mathematical convenience [3, 12, 11]. An inverse Fourier transformation is then performed on the pupil-plane field to give

$\displaystyle\psi(x,\xi)=U_{0}\int_{-\infty}^{\infty}\exp\left(-{\rm i}k\frac{\xi u}{f}\right)\exp\left(-\frac{u^{2}}{4\sigma_{p}^{2}}\right)\exp\left[-{\rm i}k\Delta W(u,\xi)\right]\exp\left({\rm i}k\frac{ux}{f}\right)\,{\rm d}u,$

(2)

where $U_{0}$ is a constant with dimensions of optical field that ensures the intensity-normalization of $\psi$. Furthermore, $u$ is the pupil-plane spatial coordinate, and $f$ is the focal length of the lenses in a $4f$ imaging system (although the results presented in this work are valid for arbitrary other imaging configurations); the Fourier optical details regarding Eq. (2) is described further in Supplement 1. A schematic of the imaging system is shown in Fig. 1. It will be useful to define $\sigma=f\lambda/(4\pi\sigma_{p})$ as the characteristic width of the diffraction-limited [$\Delta W(u,\xi)=0$] PSF. Notice, as is true for the case of a shift-invariant imaging system, the system PSF would explicitly be a function of the difference in coordinates $x-\xi$ if $\Delta W$ were independent of the object location $\xi$. Such on-axis/field-independent aberrations include well-known Zernike terms such as defocus and spherical aberration. Our work focuses instead on examples of $\Delta W$ which depend on $\xi$.

Before proceeding to the analysis of off-axis aberrations, we note that the spatial profile of $\psi$ is affected by the strength of the aberration $\Delta W$. For the purposes of quantifying this strength in later sections, it is useful to use $2\sigma_{p}$ and $2\sigma$ as standard lengths for the pupil coordinate $u$ and the object location $\xi$, respectively. The reason for this choice is due to (1) the Gaussian pupil function strongly attenuating any signal at the pupil plane located at $|u|>2\sigma_{p}$ and (2) the context of quantum-inspired resolution is most relevant for objects located at $|\xi|<2\sigma$; that is, for objects whose features are below twice the diffraction-limited PSF. This regime is henceforth referred to as the sub-Rayleigh regime.

From Eq. (1), it is clear that the PSF of a shift variant system is determined by the aberration function $\Delta W(u,\xi)$. Although there is an infinitude of options; we restrict our analysis to the case where

$\displaystyle\Delta W(u,\xi)=W_{111}\left(\frac{\xi}{2\sigma}\right)\left(\frac{u}{2\sigma_{p}}\right)+W_{220}\left(\frac{\xi}{2\sigma}\right)^{2}\left(\frac{u}{2\sigma_{p}}\right)^{2},$

(3)

which is recognized as the superposition of the two low-order off-axis Seidel aberrations commonly called OAT and Petzval curvature [18]. Notice that Eq. (3) is written in terms of ratios for $\xi$ and $u$ so that $W_{111}$ and $W_{220}$ can be interpreted as the maximal [when $\xi=2\sigma$ and $u=2\sigma_{p}$, as discussed in the paragraph following Eq. (2)] optical path difference caused by OAT and Petzval curvature, respectively. By further defining $T$ and $P$ as

$\displaystyle T\triangleq\frac{W_{111}}{\lambda\cdot(2\sigma)}\quad\text{and}\quad P\triangleq\frac{W_{220}}{\lambda\cdot(2\sigma)^{2}},$

(4)

to be OAT and Petzval curvature strength parameters, respectively (they measure rates, for $u=2\sigma_{p}$, at which the phase difference caused by OAT and Petzval curvature, respectively, increases as a function of $\xi$), one finds that the intensity-normalized shift-variant PSF is given through Eq. (2) by

$\displaystyle\psi(x,\xi;P,T)=\frac{1}{[2\pi g^{2}(\xi;P)]^{1/4}}\exp\left\{-\frac{[x-\xi(1+2\pi T\sigma)]^{2}}{4g^{2}(\xi,P)}+{\rm i}\Phi(x,\xi;P,T)\right\},$

(5)

where

$\displaystyle g(\xi;P)\triangleq\sigma\sqrt{1+4\pi^{2}P^{2}\xi^{4}},$

(6)

is the characteristic width of $\psi$ and

$\displaystyle\Phi(x,\xi;P,T)\triangleq-\frac{1}{2}\left\{\tan^{-1}(2\pi P\xi^{2})+\frac{\pi P\xi^{2}[x-\xi(1+2\pi T\sigma)]^{2}}{g^{2}(\xi;P)}\right\},$

(7)

is the phase of $\psi$. Equation (5) shows why the choice of Eq. (3) was used: OAT and Petzval curvature cause intuitive changes in the PSF’s mean location and width, respectively. Therefore, much insight can be obtained from studying Eq. (3) since, while other off-axis aberrations may induce higher order effects on $\psi$, the effects of OAT and Petzval is expected to provide a useful leading-order description.

OAT and Petzval curvature, among other Seidel aberrations, are common in typical imaging systems and have well-known interpretations particularly in the language of geometrical (ray) optics. For a system with OAT, the plane wave emerging from the pupil associated with each object location $\xi$ is proportionally tilted (posses an additional linear phase). This often arises from unintentional tilts of the optics within the system and causes an incorrect magnification in the image. This magnification factor, according to Eq. (5) and Fig. 1, is $(1+2\pi T\sigma)$; in the case of an object scene with two point sources, this magnification translates to an amplification of the separation in the image field. Petzval curvature, on the other hand, describes an imaging system in which the image is perfect when the intensity is measured along a curved surface. If the intensity is instead measured, as it usually is, along a flat surface, off-axis object points $\xi$ are mapped to a shift-variant PSF whose width increases with $|\xi|$. This width, according to Eqs. (5) and (6), this width is $g(\xi;P)$.

For the aberration-free case ($P=0$ and $T=0$), Eq. (5) reduces to the shift-invariant Gaussian PSF with width $\sigma$. Furthermore, it should be noted that the proceeding comparisons of Fisher information (FI) values between systems where the PSF is given by $\psi(x,\xi;P,T)$ and the aberration-free case $\psi(x,\xi;0,0)$ is fair, since $g(\xi;P)\geq\sigma$. In other words, the presence of OAT and Petzval curvature does not reduce the diffraction-limited spot size and therefore the improved FI is not attributed to the reduction of $\sigma$.

Although this work will primarily be concerned with incoherent imaging, a general framework is provided here for completeness and comprehension. To this end: the density matrix, which represents the cross-spectral density of the optical field at the image plane, for the case of two partially coherent point sources with intensity ratio $A$, a known centroid taken to be at the origin of the coordinate system, separated by a distance $s$ is given by

$\displaystyle\rho(s;P,T)=\frac{|\psi_{+}\rangle\langle\psi_{+}|+A|\psi_{-}\rangle\langle\psi_{-}|+\Gamma|\psi_{+}\rangle\langle\psi_{-}|+\Gamma^{*}|\psi_{-}\rangle\langle\psi_{+}|}{1+A+2\text{Re}[\Gamma d(s;P,T)]},$

(8)

which is expressed over the non-orthogonal, intensity-normalized basis $\{|\psi_{+}\rangle,|\psi_{-}\rangle\}$. These basis kets are defined over position-space through Eq. (1) as

$\displaystyle|\psi_{\pm}\rangle$

$\displaystyle\triangleq\int_{-\infty}^{\infty}\psi_{\pm}(x;P,T)|x\rangle\,{\rm d}x\triangleq\int_{-\infty}^{\infty}\psi\left(x,\pm\frac{s}{2};P,T\right)|x\rangle\,{\rm d}x.$

(9)

Furthermore,

$\displaystyle d(s;P,T)\triangleq\langle\psi_{-}|\psi_{+}\rangle=\int_{-\infty}^{\infty}\psi\left(x,\frac{s}{2};P,T\right)\psi^{*}\left(x,-\frac{s}{2};P,T\right)\,{\rm d}x$

(10)

is the field-overlap integral between the two imaged PSFs and $\Gamma$ is a complex-valued coherence parameter whose values are restricted by the condition $|\Gamma|\leq A$. The normalization factor in Eq. (8) indicates that the density matrix $\rho(s)$ uses the image-plane normalization scheme, in which the unit-trace condition on $\rho$ is ensured by the total number of received photons at the image plane. This choice is in contrast to the object-plane normalization scheme; a detailed discussion of the two options is provided elsewhere.

The act of performing a measurement on the received photons is captured mathematically through computing the modulus-square coefficients over a corresponding set of projections of $\rho$. These coefficients are taken to be the probabilities of finding an image-plane photon in a certain mode (which may be continuous or discrete). In this work, we focus on comparing the measurement schemes of DI and SPADE. For DI, the set of projection modes is the continuous position basis $\{|x\rangle\}$; the probability density of finding a photon in the $x$ position of the image plane (conditioned on a detection event) is given by

	$\displaystyle p_{\text{I}}(x,s;P,T)$	$\displaystyle=\langle x|\rho(s;P,T)|x\rangle$
		$\displaystyle=\frac{|\psi_{+}(x;P,T)|^{2}+|\psi_{-}(x;P,T)|^{2}+2\text{Re}\left[\Gamma\psi_{+}(x;P,T)\psi_{-}(x;P,T)\right]}{1+A+2\text{Re}[\Gamma d(s;P,T)]}.$		(11)

In Section 3, we focus on the case where $\Gamma=0$ and $A=1$ (equally bright incoherent point sources). Figure 2 is provided for the visualization of $p_{\text{I}}(x,s;P,T)$, which is a normalized version of the image-plane intensity distribution, for various levels of OAT and Petzval curvature. The aberration-free case shown in Fig. 2(a) shows the usual intensity pattern of two Gaussian PSFs converging as the separation $s$ vanishes. The effects of non-zero values of $T$ and $P$ are illustrated in Figs. 2(b) - (d): OAT causes the intensity distributions from the two point sources to converge at a different rate due to the magnification factor $(1+2\pi T\sigma)$ and Petzval curvature gives the intensity distributions a $s$-dependent width. Regarding OAT, it is clear from the comparison of Fig. 2(a) and (b) that non-zero values of $T$ allows for the two point sources’ PSFs to be more easily distinguished for smaller values of $s/\sigma$. This effect, which will be detailed later, is responsible for larger values of CFI when using DI for imaging systems with OAT. Petzval curvature, on the other hand, does not have a similarly clear benefit in $s$-estimation whe compared to OAT. The resolution advantages offered by non-zero values of $P$ instead comes from the sensitivity of the width $p_{\text{I}}$ as $s$ approaches zero.

For SPADE, the projection is done over any set of discrete orthonormal modes. A natural choice is the $\sigma$-matched Hermite-Gauss (HG) basis $\{\phi_{q}(x)\}_{q=0}^{\infty}$, where the $q$-th mode is defined as

$\displaystyle\phi_{q}(x)=\frac{1}{(2\pi\sigma^{2})^{1/4}}\frac{1}{\sqrt{2^{q}q!}}H_{q}\left(\frac{x}{\sqrt{2}\sigma}\right)\exp\left(-\frac{x^{2}}{4\sigma^{2}}\right),$

(12)

and $H_{q}$ is the $q$-th physicist’s Hermite polynomial. The probability of finding a photon (conditioned on a detection event) in the $q$-th mode is therefore the projection

	$\displaystyle p_{\text{II}}(q,s;P,T)$	$\displaystyle=\langle\phi_{q}|\rho(s;P,T)|\phi_{q}\rangle$
		$\displaystyle=\frac{|\langle\phi_{q}|\psi_{+}\rangle|^{2}+|\langle\phi_{q}|\psi_{-}\rangle|^{2}+2\text{Re}\left(\Gamma\langle\phi_{q}|\psi_{+}\rangle\langle\psi_{-}|\phi_{q}\rangle\right)}{1+A+2\text{Re}[\Gamma d(s;P,T)]}.$		(13)

The probability densities given by Eqs. (11) and (13) for DI and SPADE, respectively, are needed in the calculation of the CFI of either measurement scheme. For the case where only the separation $s$ is an unknown parameter, the corresponding CFIs for DI and SPADE are given by

$\displaystyle F_{\text{I}}(s;P,T)$

$\displaystyle=\int_{-\infty}^{\infty}\frac{1}{p_{\text{I}}(x,s;P,T)}\left[\frac{\partial}{\partial s}p_{\text{I}}(x,s;P,T)\right]^{2}\,{\rm d}x,$

(14)

and

$\displaystyle F_{\text{II}}(s;P,T)$

$\displaystyle=\sum_{q=0}^{\infty}\frac{1}{p_{\text{II}}(q,s;P,T)}\left[\frac{\partial}{\partial s}p_{\text{II}}(q,s;P,T)\right]^{2},$

(15)

respectively. Although a full background on SPADE has been presented thus far, it is sufficient to consider a simplified version of SPADE known as binary SPADE (BSPADE) when one is interested in the sub-Rayleigh regime. In BSPADE, a photon arriving from the object plane is measured either in the $q=0$ mode or in the combined $q>0$ mode. Intuition for this simplification comes from the fact that, as the object shrinks, most of the photons arriving at the image plane will project into the lowest order modes and the higher order modes contains progressively less information. In this case, the summation in Eq. (15) simplifies and can can be explicitly expressed as

$\displaystyle F_{\text{II}}(s;P,T)\approx\frac{1}{p_{\text{II}}(0,s;P,T)-p_{\text{II}}^{2}(0,s;P,T)}\left[\frac{\partial}{\partial s}p_{\text{II}}(0,s;P,T)\right]^{2}.$

(16)

The quantities $F_{\text{I}}$ and $F_{\text{II}}$ are valuable in that, via estimation theory, their reciprocals provide lower bounds for the variance in the estimation of $s$ when the corresponding measurement scheme is used. In other words, if $\check{s}_{\text{I}}$ and $\check{s}_{\text{II}}$ are estimators for $s$ constructed from data acquired via DI or BSPADE measurements, respectively, then

$\displaystyle\text{Var}(\hat{s}_{\text{I}})\geq F_{\text{I}}^{-1}\quad\text{and}\quad\text{Var}(\hat{s}_{\text{II}})\geq F_{\text{II}}^{-1}.$

(17)

The CFI for DI and BSPADE measurements are compared for various values of $T$ and $P$ in Section 3. There, it is shown that non-zero values of these aberration strength parameters can lead to larger values of $F_{\text{I}}$ and $F_{\text{II}}$.

In addition to the CFI, which depends on the choice of measurement scheme, it is informative to consider the QFI, $Q_{s}$, which provides an upperbound to all possible CFI once the transformation between object and image planes is stipulated. In other words,

$\displaystyle Q_{s}\geq F_{\text{I}}\quad\text{and}\quad Q_{s}\geq F_{\text{II}}.$

(18)

As is standard in deriving the QFI, one first finds the symmetric logarithm derivative (SLD), $L_{s}$, associated with the separation parameter $s$. The SLD is defined implicitly through

$\displaystyle\frac{\partial\rho(s;P,T)}{\partial s}$

$\displaystyle=\frac{\rho(s;P,T)L_{s}+L_{s}\rho(s;P,T)}{2},$

(19)

where $\rho$ is the image-plane density matrix given by Eq. (8); one should note that this prescription of $\rho$ corresponds to the image-plane normalization scheme detailed in []. With $L_{s}$ in hand, the QFI for $s$ is immediately obtained via

$\displaystyle Q_{s}(s;P,T)$

$\displaystyle=\text{Tr}\left[\frac{\partial\rho(s;P,T)}{\partial s}L_{s}\right].$

(20)

Details pertaining to the derivation of the QFI for imaging systems with off-axis aberrations are found in Supplement 1. It should be emphasized that such calculations for the separation QFI for imaging systems in which the PSF is shift-variant require care and themselves constitute a novel result from which interesting discussions may arise. However, to keep the focus on well-known measurement schemes like DI and BSPADE, the QFI will primarily serve as a confirmation of CFI results through Eq. (18).

Although the framework developed in Section 2 is general, we now specialize to the case of two equally bright incoherent point sources. This choice corresponds to the selection of $A=1$ and $\Gamma=0$ in Eqs. (8), (11), and (13). We are interested in calculating the CFI for DI and BSPADE for various values of $P$ and $T$ and comparing them to the case where the imaging system is aberration-free (and therefore shift-invariant), i.e., for $P=0$ and $T=0$.

In the following, the CFI for DI, given by Eq. (14), is calculated numerically. For BSPADE, it can be shown using Eqs. (5), (12), and (13) that

$\displaystyle p_{\text{II}}(0,s;P,T)=\frac{1}{\sqrt{1+P^{2}\pi^{2}(s/2)^{4}}}\exp\left\{-\frac{(s/2)^{2}(1+2\pi T\sigma)^{2}}{4\sigma^{2}[1+P^{2}\pi^{2}(s/2)^{4}]}\right\},$

(21)

which may be inserted into Eq. (16) to obtain $F_{\text{II}}$. Comparisons of $F_{\text{I}}$ and $F_{\text{II}}$ for various values of $T$ and $P$ are displayed in Fig. 3. When considering the $P=0$ case (the imaging system only has OAT and no Petzval curvature), Fig. 3(a) shows that increasing values of $T$ leads to greater $F_{\text{I}}$ for all values of $s/\sigma$: this global improvement in $F_{\text{I}}$ for non-zero values of $T$ is one of the main results of this analysis. For BSPADE, Fig. 3(b) shows that a similar improvement in $F_{\text{II}}$ occurs when $T\neq 0$ in the sub-Rayleigh regime. Additionally, notice that $F_{\text{II}}$ saturates the QFI at exactly $s=0$ for all values of $T$, which indicates that BSPADE continues to be an optimal measurement for $s$ even for shift-variant imagine systems. However, we re-emphasize that the main purpose of this analysis to to show that introducing off-axis aberrations in a system can vastly improve CFI even if one only considers a DI measurement scheme. The introduction of OAT with strength $T=0.4\sigma^{-1}$ leads to very high information (in comparison with the QFI value of $0.25$ in the aberration-free case) well into the sub-Rayleigh regime ($s<2\sigma$). That is, even though $F_{\text{I}}=0$ at exactly $s=0$ always, one in principle can obtain any $F_{\text{I}}$ for non-zero $s/\sigma$ by increasing the strength of OAT even if $s/\sigma$ is small.

Figures 3(c) and (d) consider the case where the imaging system contains Petzval curvature, but no OAT. It can be seen, in a behavior similar to that seen in Figs. 3(a) and (b) that larger values of $P$ lead to larger values of $F_{\text{I}}$ and $F_{\text{II}}$ in the sub-Rayleigh regime, with the latter saturating the QFI near $s=0$. Therefore, like OAT, Petzval curvature can improve the performance of an imaging system regarding two-point separation estimation. However, there are some notable differences between the behaviors of the CFI when the system contains OAT and when it contains Petzval curvature. First, increasing $P$ does not affect the value of the $F_{\text{II}}$ at exactly $s=0$ ($F_{\text{I}}=0$ regardless of $P$ and $T$). In other words, all the curves seen in Fig. 3(d) coincide in the limit of vanishing separation. Second, there is an interval of $s$ within the sub-Rayleigh regime in which DI outperforms BSPADE. For example, comparing the $P=0.4\sigma^{-2}$ case for $F_{\text{I}}$ and $F_{\text{II}}$ in Figs. 3(c) and (d), repsectively, it can be seen that the former reaches a peak of roughly $1.5$ compared to the latter’s $0.75$. Finally, we note that Petzval curvature may lead to the formation of local minima in $F_{\text{I}}$ and $F_{\text{II}}$, which indicate that the corresponding measurement schemes perform optimally neither at exactly $s=0$ nor at $s\rightarrow\infty$, but rather at some intermediate value that, according to Figs. 3(c) and (d), occur in the sub-Rayleigh regime.

For completeness, CFI and QFI curves for the case where both OAT and Petzval curvature are present are shown in Fig. 4 for a fixed value of $T=0.2\sigma^{-1}$ and varying values of $P$. As was observed from Fig. 3, larger values of $T$ raise the maximal attainable information while various values of $P$ lead to CFI curves that peak locally in the sub-Rayleigh regime [although such peaks for the BSPADE CFIs in Fig. 4(b) are not apparent for smaller values of $P$ since the nonzero value of $T$ raises the CFIs at $s=0$].

It is prudent to point out the behavior of the QFI curves in Figs. 3 and 4 as well. When only OAT is present, the QFI remains constant over $s$, with its value given by

$\displaystyle Q_{s}(s;0,T)=F_{\text{II}}(0;0,T)=\frac{(1+2\pi T\sigma)^{2}}{4\sigma^{2}},$

(22)

which indicates that $Q_{s}(s;0,T)$ grows quadratically as a function of $T$. On the other hand, when Petzval curvature is present, the QFI increases from the value given by Eq. (22) at $s=0$ to larger values as $s$ increases. This peculiar behavior, which is had not been seen for incoherent quantum-inspired superresolution studies of shift-invariant systems, is a result of the presence of the second term in the phase, $\Phi$, of the PSF given by Eq. (7). This quadratic (in $x$) phase contribution is due to the presence of Petzval curvature in $\Delta W$, as stipulated in Eq. (3); Supplement 1 provides more insight regarding the QFI’s behavior.

Maximum likelihood estimation (MLE) simulations for the separation $s$ were performed for DI to support the findings from Eq. (14) and Fig. 3(a) and (c). For the $i$-th (out of $M$) iteration of the simulation, the number of photons arriving at the image plane, $N_{i}$, was chosen using Poisson statistics around an average photon number of $N$. That is, the probability mass function for $N_{i}$ is given by

$\displaystyle p(N_{i})=\frac{N^{N_{i}}\exp(-N)}{N_{i}!}.$

(23)

Once $N_{i}$ is chosen, an equivalent number of photon positions in the image plane is chosen using Eq. (11) as the probability density function for $N_{i}$ independent events. Given these $N_{i}$ photon locations, the unbiased MLE estimator, denoted as $\hat{s}$, is used to find an estimate for the separation $s$. The variance of $\hat{s}$, whose average value is $s$, is then calculated over the ensemble of $M$ iterations. In order to compare the simulation results to $F_{\text{I}}$, one must divide the reciprocal of the aforementioned variance by $N$ to compute the per-photon CFI. This process is repeated for different values of $s,T,$ and $P$ and the results of these simulations are shown in Fig. 5 over the sub-Rayleigh regime $(s<2\sigma)$ for $N=2000$ photons and $M=500$ iterations. As is done in Figs. 3(a) and (c), the cases of $P=0$ and $T=0$ were analyzed, respectively. It is evident that there is good agreement between the simulation results and $H_{\text{I}}$.

The primary result of this work is that it is possible to improve upon the two-point resolution of an aberration-free imaging system by intentionally introducing known off-axis aberrations. Particularly, this is true even for DI measurement schemes as seen in Figs. 3 and 4. In fact, whether the improvement is global (over all $s$, as it is when OAT is present) or only over the sub-Rayleigh regime (when only Petzval curvature is present), the improvement can exceed the improvements offer by modal imaging techniques like BSPADE for aberration-free systems. In other words, if a sufficient amount of off-axis aberration can be provided, DI imaging schemes can lead to large $F_{\text{I}}$ even in the sub-Rayleigh regime. Since practical applications of two-point imaging deal with non-zero values of $s/\sigma$, increasing $T$ can effectively allow DI measurement schemes to attain the same information allowed by BSPADE in shift-invariant systems. However, this is not to say that BSPADE (and other modal imaging schemes) are fruitless; it is clear from Figs. 3 and 4 that BSPADE benefits in a similar fashion to DI when off-axis aberrations are introduced. Importantly, at least for OAT and Petzval curvature, BSPADE still (as it did for the aberration-free case) yield a $F_{\text{II}}$ that saturates the QFI in the sub-Rayleigh regime.

The findings in this work are somewhat counter-intuitive as many traditional imaging systems start with substantial aberrations and a significant amount of effort is often needed to minimize such errors towards the goal of an aberration-free (diffraction-limited) system. However, at least in the context of resolving two incoherent and equally bright point sources where the separation is the only unknown parameter, it turns out that the presence of off-axis aberrations can actually improve performance. It may be of interest then, given the number of well-corrected imaging systems in existence, to understand how one may intentionally introduce OAT and Petzval curvature back into these imaging systems in order to take advantage of the larger DI or BSPADE CFI. Of course, it is possible to purposely (re)introduce aberrations by misaligning or adding/subtracting optical elements in the system. However, we should mention that despite their simple geometrical optics interpretations, OAT and Petzval curvature cannot be properly introduced by simply tilting or curving the image plane, respectively. Although these geometrical manipulations of the image surface create effective PSFs that are SV, they do not have the necessary dependence on $\xi$, the object plane coordinate, in order to bring about the results from a genuine $\Delta W$ given by Eq. (3).