Coronagraphic Data Post-processing Using Projections on Instrumental Modes

Directly observing exoplanets with coronagraphs is impeded by the presence of speckles from aberrations in the optical path, which can be mitigated in hardware with wavefront control as well as in post-processing. This work explores using an instrument model in post-processing to separate astrophysical signals from residual aberrations in coronagraphic data. The effect of wavefront error (WFE) on the coronagraphic intensity consists of a linear contribution and a quadratic contribution. When either of the terms is much larger than the other, the instrument response can be approximated by a transfer matrix mapping WFE to detector plane intensity. From this transfer matrix, a useful projection onto instrumental modes that removes the dominant error modes can be derived. We apply this projection to synthetically generated Roman Space Telescope hybrid Lyot coronagraph (HLC) data to extract"robust observables,"which can be used instead of raw data for applications such as detection testing. The projection improves planet flux ratio detection limits by about 28% in the linear regime and by over a factor of 2 in the quadratic regime, illustrating that robust observables can increase sensitivity to astrophysical signals and improve the scientific yield from coronagraphic data. While this approach does not require additional information such as observations of reference stars or modulations of a deformable mirror, it can and should be combined with these other techniques, acting as a model-informed prior in an overall post-processing strategy.


INTRODUCTION
Specialist high-contrast techniques are required to directly observe faint astrophysical objects near brighter objects, such as exoplanets, brown dwarfs, or circumstellar disks orbiting much brighter central stars.High contrast observations are essential for answering scientific questions involving binary and planetary system population statistics, planet and disk formation and evolution, planetary atmospheres, and planet habitability and the search for biosignatures (Traub & Oppenheimer 2010).Measuring these exoplanet signals is difficult because they often lie at small angular separations from their host star and can be many orders of magnitude fainter.One major obstacle for high contrast observations is photon noise from the light of the central star.Practical matters such as detector saturation aside, if the star is orders of magnitude brighter than its companion, the photon noise associated with the outer lobes of star's point-spread-2.CORONAGRAPHIC SIGNALS

Data Formation
The model used in this work assumes the light through the instrument is monochromatic.With a discrete representation of the optical planes of an instrument, a coronagraph can be modeled as a linear operator C, a constant 2D matrix transforming the electric field vector at the pupil plane, E s , into the electric field vector at the detector plane, E det .If E s0 is the electric field vector of the central source (star) at the pupil plane in the absence of aberrations, and ∆E s a vector of small perturbations to that electric field, representing wavefront aberrations (which can be variable in time), then the electric field vector at the detector plane, assuming that the star is the only source of light, is (1) The intensity measured is the element-wise norm squared of the detector plane electric field (here, x indicates the element-wise complex conjugate of x and • indicates the element-wise product): (2) The vector of the pupil plane electric field of a binary companion is given by where c is the flux ratio between the planet and the star, k is the pupil plane wave vector indicating the companion's location, and x is the pupil plane coordinate vector.Namely, the planet's pupil-plane electric field is the star's electric field, but tilted and scaled by the square root of the flux ratio.The detector plane intensity for the planet can be expressed as (4) The total intensity on the detector plane from the star and the planet is the sum of Eqs. 2 and 4.However, we can make two simplifying assumptions.The first assumption is that the flux of the planet is small relative to the flux of the star, such that c ≪ 1.The second assumption is that the magnitude of the wavefront error is small relative to the total magnitude of the electric field, namely ∆E s (t) ≪ E s0 , which implies ∆E p (t) ≪ E p0 .This is true if we are both in the "small phase regime" (when there is much less than one wave of wavefront error) and the fractional amplitude error is much less than one.These assumptions imply that the last two terms of Eq. 4 are small relative to the other terms, so we can approximate the total intensity as (5) The first term |CE s0 | 2 is the residual starlight not blocked by the coronagraph in the case of no aberrations.The second term 2 Re (CE s0 ) • C∆E s is linear in the wavefront aberration and corresponds to the interference of the aberration, propagated to the focal plane, with the underlying residual starlight from the coronagraph, analogous to speckle pinning (Bloemhof 2002;Perrin et al. 2003).The third term |C∆E s | 2 is the quadratic term, corresponding to the norm squared of the wavefront error propagated to the focal plane.The last term |CE p0 | 2 is the nominal off-axis signal of interest.
Whether the effect of wavefront errors at some location in the detector plane are dominated by the linear term or the quadratic term depends on the attenuation of starlight by the coronagraph and the level of the propagated wavefront error at that location.If the propagated wavefront aberrations are smaller in complex amplitude than the residual starlight after the coronagraph with no aberrations, the linear term is dominant.When a coronagraph is not used, this corresponds to the speckle pinning regime, in which the aberrations primarily interfere with the wings of the telescope's PSF (Bloemhof 2002).The same phenomenon occurs with a coronagraph; however, as the amplitude of the PSF wings are reduced by the coronagraph, the range of WFE over which this occurs is much more limited.Otherwise, if the propagated wavefront aberrations have relatively larger magnitudes, the quadratic term is dominant.For a given location in the focal plane, the local point of transition between the linear and quadratic regimes occurs when This point of transition is different for each pixel, and also depends on the coronagraph design as well as the "nominal" wavefront (whether it is flat, as is typical for ground-based coronagraphs, or the wavefront corresponding to a dark hole, as is planned for space-based coronagraphs).For this work, we use as an example the Hybrid Lyot Coronagraph (HLC) of the Coronagraph Instrument of the Roman Space Telescope.With the dark hole presented in Section 5.1, which has an average raw contrast (residual stellar intensity divided by unocculted peak intensity) of 5.6 × 10 −9 , the point at which |C∆E s | > |2CE s0 | for 50% of the pixels in the dark hole region occurs at roughly 0.1 waves root-mean-square (RMS) of phase error, on average.This means that wavefront errors less than 0.1 waves RMS will primarily be in the linear regime, while wavefront errors larger than 0.1 waves RMS will primarily be in the quadratic regime, although this is somewhat dependent on the form of the wavefront's spatial power spectral density (PSD) that we use in Section 5.4.1.
In this work, robust observables are only formulated for WFE that is predominantly linear or predominantly quadratic throughout the entire focal plane.However, it may be possible to obtain robust observables for when both terms have comparable contributions, a topic that is left for future work.

Linear Regime
From Equation 5, if we then assume that the linear error term is dominant, then we can drop the quadratic contribution such that the detector plane intensity is approximately The contribution of the wavefront error to the intensity can be expressed as a linear transformation A l acting on the wavefront error: The transfer matrix A l can be calculated semi-analytically from the coronagraph operator and the unaberrated electric field, as derived from Equation 6: The indices i and j label the input basis vectors used to represent the wavefront error, and the index k labels the detector pixel.
It is desirable to reduce the term dependent on WFE, A l ∆E s (t), relative to the terms containing astrophysical signals of interest.This can be achieved by left-multiplying the measured intensities by a matrix K l , that projects out the dominant modes of A l .The following Section 3 describes the process of calculating A l and finding from it an appropriate K l .The observables obtained using projection matrix K l are given by When the wavefront errors are in the linear regime, this projection is expected to suppress the contribution of wavefront errors to the measured data.As long as the measurements retain most of the astrophysical signal, then the projection will boost its SNR.

Quadratic Regime
In the quadratic-dominated regime, we can drop the linear contribution in Eq. 5, such that the detector plane intensity is approximately In a discrete numerical model, the contribution of the quadratic term to each pixel labeled k in the detector plane can be expressed as The indices m, n, i, and j label the input basis vectors used to represent the wavefront error (expressed here in terms of perturbation to the complex electric field), and the index k labels the detector pixel.The quantity M with elements M kij is a 3-tensor containing the second order partial derivative matrix (Hessian) of each pixel intensity with respect to the wavefront error, and relates each pairwise combination of pupil basis vectors to its effect on each detector plane pixel k.Each entry can be calculated semi-analytically from the coronagraph operator using the following formula derived from Equation 11: Assuming there are N pix pixels of interest on the detector, and N basis vectors are used to represent the wavefront error, then, through a remapping, the 3-tensor M of size (N pix × N × N ) can be expanded into a matrix acting on the space of all pairwise combinations of pupil basis vectors.Since Hessians are symmetric because partial derivatives commute (M kij = M kji ), the ordering of each pair of segments does not matter, and the derivatives corresponding to the same pair of original basis vectors can be consolidated into the same entry.This results in an input vector space of size N +1 The 3-tensor M can thus be represented as a (N pix × N +1 2 ) matrix A q of second derivatives, acting on a vector β of perturbations defined for each pairwise combination ∆E si ∆E sj .This results in the following expression for the quadratic term: The projection is similar to the linear case: the detector intensities can be left-multiplied by a matrix K q that projects out the dominant quadratic error modes of M. The observables with the appropriate projection K q are given by O q = K q I tot,q . (14) 3. RESPONSE MATRICES AND ROBUST OBSERVABLES

Calculating the Response Matrix
This section details the numerical calculation of instrument response matrices and the projection matrices.In this work, the response matrix is calculated with the wavefront aberrations represented in the Zernike basis.In this basis, ∆E Zn is the coefficient of the aberration induced by the n th Noll ordered Zernike polynomial (Noll 1976), and N is the total number of polynomials chosen to construct the response matrix: We define N pix as the total number of detector pixels of the optical model and N basis as the number of Zernike modes to include.The coronagraph operator C is the N pix × N basis matrix that, when applied to a vector of Zernike coefficients, gives the perturbation they induce in the focal plane electric field.This operator is typically either already part of the optical model, or obtainable by propagating Zernike modes through the optical model and using finite differences to populate its columns.Given the operator C and the initial unaberrated focal plane electric field, we can calculate both A l and A q using Equations 8 and 12.Note that the term ( k C kj E s0j ) in Equation 8 is simply the initial unaberrated focal plane electric field at pixel k.For more complicated models without simple analytical solutions (such as those that include distortion), automatic differentiation, in which arbitrary exact derivatives can be computed without finite differences, may be useful (Pope et al. 2021).
The linear transfer matrix poses no computational problems, as its size is N pix × N basis .However, for the quadratic transfer matrix A q , the size of the input dimension quickly becomes computationally burdensome for high N basis .For the example system shown in Section 5, a N basis of 528 results in a A q matrix of width N basis +1 2 = 139, 656 (the number of pairwise combinations of pupil basis vectors), and length 5,476 (the number of detector pixels of the model).This A q matrix, when represented as (non-complex) doubles, is over 6GB in size.As explained in Section 3, the calculation of the projection matrix involves a singular value decomposition (SVD) of the response matrix.Since calculating the SVD of a matrix of this size is too computationally expensive, we restrict our quadratic transfer matrix to only include the first N redu = 100 Zernikes, which results in an A q with a width of only N redu +1 2 = 5050.This model is valid only in a smaller area closer to the central star -namely within ∼ 5 λ/D, where λ is the wavelength and D the telescope diameter.However, in Appendix A, we explore using an approximation of the quadratic transfer matrix that can extend the area of applicability while circumventing impractical computational costs.

Calculating the Projection Matrix
Once a response matrix A is obtained, a singular value decomposition of A = USV T is performed, revealing its singular modes and corresponding singular values.Then a choice of the number of modes to project out (N m ) is made.The remaining N pix − N m modes are kept in the post-processing projection K. Accordingly, K is the subset of U that contains the m + 1 th and higher left singular modes of A. A pseudo-code summary of the process to find K is given in Algorithm 1.The optimal N m depends on the signal of interest.For the point-source companion signals explored in this work, N m is chosen as the cutoff that results in the best detection limit at the separation of interest.

Xin et al.
Algorithm1 Calculate projection matrix K Input: Transfer matrix A Input: Cutoff mode (number of modes to project out), Nm Input: Indices of detector plane pixels in region of interest, idx Output: If the linear and quadratic projections are used in the appropriate regimes to increase SNR, they could, for example, allow for a binary signal detection with a deeper flux ratio than using the raw intensity data.Detection tests can be performed on both projected and unprojected data to quantify this effect.

DETECTION TESTING
Detections are typically claimed from a statistical hypothesis test (see e.g.Kasdin & Braems 2006;Jensen-Clem et al. 2017;Ceau et al. 2019).A test statistic T is calculated from the data and compared to a threshold ξ.A detection is claimed if T ≥ ξ, and a lack of a detection is claimed otherwise.The fraction of real companions detected is the true positive rate (TPR).A false positive occurs if there is no companion in the data, but the detection test incorrectly claims a detection.The rate at which this occurs is the false positive rate (FPR).
As the detection threshold ξ is decreased, detecting real companions becomes more likely, but false detections also become more likely (Jensen-Clem et al. 2017).Varying the threshold and plotting the TPR as a function of the FPR results in a receiver operating characteristic (ROC) curve, an example of which is shown in Figure 6.ROC curves characterize the performance of a detection scheme and are used in the determination of flux ratio detection limits.
This work uses a simple Delta Reduced Chi-squared (∆χ 2 r ) statistic, or the difference in the reduced χ 2 of the data assuming it contains only noise, and the reduced χ 2 of the data assuming it contains noise and the companion signal.The formula for calculating this test statistic from the data is given by Equation 16(the bars indicate vector norm, the divisions are element-wise): In this formula, y is the data, which is the synthetically generated realizations of I tot , with or without a planet.Meanwhile, x is the unaberrated model of the planet signal I p0 = |CE p0 | 2 (assuming that it is known, such as through a maximum-likelihood-estimation).The estimated uncertainty of the data is denoted by σ, and ν is the degrees of freedom (the number of data elements minus the number of free parameters; a binary system's three free parameters are the flux ratio, separation, and position angle).This use of this test statistic is motivated by an assumption that the noise is Gaussian and uncorrelated, under which this quantity is related to the relative log-probabilities of the data containing both the planet signal and noise, versus containing only noise.The noise being uncorrelated and Gaussian is generically not the case.However, the effects of the correlation and non-Gaussianity of the injected noise on the resulting test statistic distributions is properly simulated and captured by the Monte Carlo methods used in this work.

EXAMPLE: NANCY GRACE ROMAN SPACE TELESCOPE HYBRID LYOT CORONAGRAPH
In this section, the use of robust observables with the Hybrid Lyot Coronagraph of the Roman Space Telescope is analyzed through simulation.However, this approach could also be applied to other coronagraphs, as long as the exposure times are short enough such that wavefront error has not been averaged out.The optical model of CGI is shown in Figure 1 (Kasdin et al. 2020).The optical elements corresponding to the HLC mode (the relevant mode for this work) are depicted in the top panel.

Optical Model
The HLC operates around a dark hole state, which is obtained using focal plane wavefront control (with deformable mirrors) to measure and minimize the electric field in the detector plane.Such focal plane wavefront control significantly suppresses the amount of starlight in the dark hole, and allows for much deeper raw contrasts than with just a flattened wavefront.Before an observation, the dark hole is generated using high order wavefront sensing and control loop on a bright reference star.Then, the DM shapes are fixed, and the telescope slews to the target star for the observation.During the observation, wavefront errors accrue as a result of instrumental disturbances and drifts.This work aims to mitigate the effects of those wavefront errors in post-processing.Note that as a result of the dark-hole generation, the nominal electric field E s0 is not a flat wavefront, but the pupil plane electric-field obtained at the end of the dark-hole generation sequence.
A Lightweight Space Coronagraph Simulator (LSCS)1 derived from the HLC model in the Fast Linearized Coronagraph Optimizer (FALCO; Riggs et al. 2018) toolbox is used for the following simulations.The LSCS relies on the HLC numerical model and focal-plane wavefront control algorithm included in FALCO to first generate the initial dark hole electric field.The numerical model in FALCO is also used to calculate C from the finite-difference sensitivities of the focal plane electric field to pupil plane phase error expressed in the Zernike basis (we have made the assumption that the matrix transformation is approximately linear in phase, valid when the phase error is much less than a wave).Although we use finite-differences to calculate C, one could also construct it by multiplying together all the matrix transformations of the optical model.These simulations are conducted at a single wavelength of 546 nm.
The average raw contrast of the initial dark hole is 5.6 × 10 −9 .The LSCS model takes in Zernike coefficients for phase aberrations, calculates their effect on the focal plane electric field, and adds them to the initial dark hole electric field to obtain the focal plane electric field in the presence of wavefront errors.The intensity can be calculated as the norm-square of the focal plane electric field.Detector and photon noise are not simulated.Since the default LSCS models only the first 136 Zernikes, FALCO is first used to extend the LSCS model to 528 Zernikes in order for the entire dark hole to be sampled.This results in using 528 Zernikes to sample the entire dark hole, or a N basis of 528.The LSCS models a detector that is 74 × 74 pixels, with 3 pixels per λ/D, for a total pixel number of N pix = 5476.The number of pixels defined to be in the dark hole is N DH = 2608.This model does not consider the effects of amplitude errors, and only analyzes phase errors, which, from end-to-end modeling of Roman CGI, are expected to be the dominant form of dynamic aberrations (Krist et al. 2023).However, for a system where dynamic amplitude errors are comparable to dynamic phase errors, both should be included.

Response Matrices
The Zernike coefficient drift values from the Observing Scenario simulations (OS 9; Krist 2020), based on physical modeling of the telescope, indicate that the WFE expected on Roman will fall within the linear regime of this dark hole.However, the level of wavefront error may end up being higher than currently expected.Additionally, on ground based telescopes, wavefront error from adaptive optics residuals is typically in the quadratic regime.Therefore, for illustrative purposes, both a linearly-dominated noise model and quadratically-dominated noise model are examined.
The matrices A l and A q are calculated according to Section 3.1.The linear matrix includes all Zernikes present in the optical model, and thus has an input dimension of N basis = 528.The quadratic matrix includes only the first 100 Zernikes, and thus has an input dimension of N redu +1 2 = 5, 050.The relevant dimensions of the objects used in this analysis are listed in Table 5.2.Note that the cutoff number N m is a variable to optimized over.From both the linear and the quadratic transfer matrix, model-based projection matrices with a range of cutoff modes are calculated according to Algorithm 1.To rule out the effect of dimensionality alone on the dataset, random projection matrices of the same size are also generated.This is done by taking the SVD of a matrix the same size as the A matrices, but populated with values drawn uniformly from -1 to 1, and then removing the same number of dominant modes as is done with A. These matrices are applied to synthetically generated data to quantify their effect on the detectability of binary companion signals.

Synthetic Data Generation
FALCO is used to generate a library of off-axis PSFs corresponding to the dark hole state, which can be injected as binary companions.These off-axis PSFs do not incorporate any WFE that is added on top of the dark hole state.However, the effect of WFE on the off-axis signal is expected to be much, much smaller than its effect on the on-axis stellar signal, so not modeling the effects of wavefront error on the off-axis signal should have a negligible impact on the data.
The optical system is first initialized in the dark hole state.Two noise models are considered, one in the linear regime, and one in the quadratic regime.Each dataset thus consists of 20 instantaneous frames of independent noise realizations.For each frame, the spatial PSD given in Equation 17 is used to generate the wavefront error.
In this equation, n z is Noll-ordered index of the Zernike coefficient.The normalization parameter a is chosen to be 10 nm for the linear regime, and 130 nm for the quadratic regime.The power law exponent b is chosen to be -2.These PSDs correspond to an average wavefront error (calculated over 100 realizations) of about 7 nm (0.013 waves) RMS for the linear regime data, and about 110 nm (0.2 waves) RMS for the quadratic regime data.As discussed in Section 2.1, the linear-quadratic transition occurs at approximately 0.1 waves RMS.Although 110 nm RMS of dynamic wavefront error is unrealistically high for the Roman HLC, we include this regime for demonstration purposes, as this level of WFE would be relevant on ground-based telescopes.
The resulting 528 Zernike coefficients are propagated through the LSCS to calculate the resulting dark hole intensities.In order to create data with an injected companion planet, the off-axis PSF at the desired separation is scaled by the companion's flux ratio, and then added to the dark hole intensity.The separation of the injected companion is set to be 6.5 λ/D in the linear case (which is the middle of the dark hole) and 4.0 λ/D in the quadratic case (since the model is only valid within ∼ 5 λ/D).The position angles of both are set to be 0. Frames without the injected companion Singular Modes of Aq Figure 3. Top: The first 10 singular modes of A l as represented in the the detector plane intensity basis (linear scale).Bottom: The first 10 singular modes of Aq as represented in the the detector plane intensity basis (linear scale).The HLC design is nearly circularly symmetric, broken only by the six secondary mirror struts (which can also be seen in the Lyot stop).Because the quadratic transfer matrix depends only on the coronagraph operator C, its singular modes exhibit cosine and sine-like azimuthal behavior associated with circularly symmetric operators.However, the linear transfer matrix depends on both C as well as on the focal-plane electric field at the end of dark hole creation, which is random and not circularly symmetric.Thus, its singular modes show no such symmetry structures.These singular modes correspond to the intensity patterns most likely to be attributed to wavefront error.Meanwhile, the companion's intensity pattern (the PSF at its location) overlaps very little with these dominant modes, so its signal is mostly retained when the dominant modes are projected out.
are used for the control case.Figure 4 shows example data frames: the initial dark hole, example frames with the aberrations from both noise models applied, and the same frames with injected companion signals.The flux ratio of the companion is 2 × 10 −7 for frame with linear-regime errors and 5 × 10 −6 for the frame with quadratic-regime errors.These flux ratios correspond to particularly bright planets chosen to be visible by eye.
It is worthwhile to examine how well the response matrices calculated in Section 5.2 can reconstruct the intensity errors present in the synthetic data.Figure 5 compares the intensity error resulting from WFE as calculated from the optical model with the intensity error calculated by multiplying the WFE by the appropriate response matrix, for example frames in both the linear and the quadratic regimes.In both regimes, the response matrices largely reproduce the spatial structure of the intensity error from the optical model.

Processing Synthetic Data
The quantity |CE s0 | 2 is the initial dark hole intensity without any extra WFE applied (as determined from the data at the end of the dark-hole digging sequence on the reference star, for example).This nominal signal is first subtracted from each frame.Then, the pixels within the defined dark hole are gathered into the vector ∆I DH .The data is left-multiplied by the appropriate K matrix to obtain the observables O = K∆I DH .The data is also left-multiplied by the random matrix of the same size as K to obtain data whose dimension has been reduced randomly.For each case, the average of the data over the twenty frames is used as the final measurement, while the standard deviation of the frames is used as the measurement uncertainty.Note that the process outlined does not rely on reference stars or dithering by deformable mirrors, and can be used even on observations for which reference observations or wavefront diversity is unavailable.

Flux Ratio Detection Limits
Detection tests are applied to these measured intensities and observables in order to characterize the detectability of a companion with these measurements.Detection limits are determined using the Monte Carlo method.One thousand random datasets are generated for each noise model with a given flux ratio.Each dataset is processed as raw intensity data, and with each projection matrix with a different cutoff mode, and the ∆χ 2 r values are calculated for each case.Figure 6 shows example histograms of the resulting ∆χ 2 r values for a c = 5.4 × 10 −7 at 4.0 λ/D companion with the quadratic noise model, as well as the corresponding ROC curves, for the projection matrix with cutoff mode N m = 70 (which, as shown in Figure 7, is the optimal cutoff at this spatial separation).The ROC curve shows that while using the robust observables results in a FPR = 0.01 and TPR = 0.9 detection of the injected companion, both the raw intensity and the randomly dimensionally reduced data remain very far from detectability.
This process is repeated for a range of flux ratios (to a precision of two significant figures).The resulting FPR = 0.01 and TPR = 0.9 detection limits for both regimes, as a function of cutoff mode N m , are shown in Figure 7.Note that these flux ratio detection limits are not based on any statistical assumptions or extrapolations, but rather real FPRs and TPRs calculated by analyzing one thousand synthetically generated datasets, with injected companions of the given flux ratios and separations.The results show that with the linear-regime noise model, the robust observables increases the detectability of a companion at 6.5 λ/D by 28%.With the quadratic-regime noise model, using robust observables increases the detectability of a companion at 4.0 λ/D by over a factor of two, and the improvement is not particularly sensitive to N m beyond the first few modes.For the linear regime, this approach can also easily f.The difference between the response matrix prediction and the optical model prediction, plotted on the same scale as (d) and (e).Slight differences arise because the model includes both the linear and quadratic error terms while the matrix predictions only include one or the other, i.e. the linear matrix prediction neglects the contribution of the quadratic term and the quadratic matrix prediction neglects the contribution of the linear term (as well as the influence of any Zernikes past the first 100).While the linear matrix prediction is biased low near the peaks and the quadratic matrix prediction biased high overall, our method depends only on how well the spatial structure of the errors are reproduced.A relevant metric for characterizing the spatial overlap is the normalized inner product between the optical model prediction and the transfer matrix prediction, where a value of 1 indicates perfect spatial overlap and a value of 0 indicates perfect spatial orthogonality.In this case, the normalized inner product is 0.936 for the linear regime example and 0.985 for the quadratic regime example, sufficient for providing a quantifiable improvement in detection sensitivity.
be extended to companions throughout the entire dark hole, though significant computation would be required to optimize N m at all separations.For the quadratic regime, our model is only valid within ∼ 5 λ/D, though Appendix A discusses a method that can be used to extend the spatial coverage without incurring impractical computational costs.

Temporally Correlated WFE and Compatibility with Other Post-processing Techniques
This work aims to characterize the effect of using robust observables in isolation.Thus, only noise models in which the WFE is uncorrelated in time are examined, since additional post-processing techniques are typically used to handle time-correlated data.Robust observables are compatible with these other post-processing techniques, and can serve as an instrument-motivated prior in the overall post-processing strategy.For example, random errors can first be reduced by projecting the data into a subspace that is robust to wavefront error.Then, reference observations along with PCA-based methods such as KLIP (Soummer et al. 2012) can be used to calibrate static and quasi-static errors and de-correlate the frames in time.This is similar to the calibration approach used in non-redundant aperture masking (NRM) interferometry or kernel-phase interferometry, in which data is projected onto closure-phases or kernel-phases 0.004 0.002 0.000 0.002 0.004 0.006 0.008 and is located at 4.0 λ/D.Left: Histograms from using raw intensities compared to those from using quadratic robust observables with the optimal cutoff of Nm = 70.The histograms using raw intensity overlap significantly, making it difficult to distinguish between a model with a planet and a model without one, while the histograms using the robust observables are further separated and more distinguishable.Middle: Histograms for using raw intensities and a random projection matrix of the same size as the instrumentally-motivated projection.Both sets of histograms overlap significantly, and the random projection does not improve the distinguishability of the two models.Right: ROC curves corresponding to the histograms.Grey area indicates false positive rates which are not well sampled as they involve less than 3 datasets with false detections.The ROC curve shows that while using the robust observables results in a FPR = 0.01 and TPR = 0.9 detection of the injected planet, both the raw intensity and the randomly dimensionally reduced data remain very far from detectability.

Raw Intensity Observables
Figure 7. Flux ratio detection limits (FPR = 0.01, TPR = 0.90) for a binary companion (to two significant figures) as a function of cutoff mode.Upward triangles indicate that a projection matrix with the specified cutoff mode performs worse than using the raw intensity, which occurs when the modes the majority of the planet signal overlaps with have also been projected out.Left: Linear regime with a companion at 6.5 λ/D.The optimal cutoff mode is 2,727, which results in a detection limit of 2.8 × 10 −9 .Unshowable in log-log scale is the detection limit with Nm = 0, which, with observables, is 3.9 × 10 −9 .This is, as expected from the fact that no error modes are removed, equal to the raw intensity detection limit.Right: Left: Quadratic regime with a companion at 4.0 λ/D.The optimal cutoff mode is 70, which results in a detection limit of 5.4×10 −7 .Unshowable in log-log scale is the detection limit with Nm = 0, which, with observables, is 1.4 × 10 −6 .This is, as expected from the fact that no error modes are removed, equal to the raw intensity detection limit.
respectively, which are then calibrated based on reference observations (Martinache 2010;Ireland 2013;Pope et al. 2021).A more sophisticated approach would be to formulate post-processing as a statistical inference problem, where a least-squares fit with the reference frames makes one up term in the cost function, and a prior over the instrumental modes (e.g.weighted by the singular value spectrum) makes up another term.Ygouf et al. (2016) shows that for the time-varying wavefront error expected on the Roman Space Telescope HLC, classical PSF subtraction with a reference observation increases the contrast gain by a factor of a few to about ten, depending on the scenario.Future work includes investigating how much overall post-processing gain can be achieved when robust observables and calibration strategies are combined, and which hybrid strategies maximize the sensitivity that can be obtained with all available information.

PSD Engineering
The robust observables derived in this work are agnostic to the actual temporal or spatial PSD of the static and dynamical wavefront errors, and are intended to be applied when these PSDs are not well-known or imperfectly characterized.As of today, this is the case for all ground-based instruments (as predictions of the influence of the atmosphere are quite imperfect), and space-based missions (as HST and JWST observatory level key metrics for requirements are expressed in terms of encircled energy, not contrast).However, it has been proposed that for future space telescope coronagraphs, the telescope WFE PSD must comply with stringent requirements in order to facilitate exoplanet detection (Nemati et al. 2020).
For instance, the PASTIS approach (Leboulleux et al. 2018;Laginja et al. 2019) considers the effects of the quadratic response on the average intensity contrast over the entire dark hole (or region of interest), calculating which modes the coronagraph is most sensitive to in order to determine stability tolerances for the segments accordingly.Calculating robust observables for post-processing is akin to doing PASTIS backwards, where the modes the coronagraph is most sensitive to are calculated in order to project them out of the data.For such telescopes, that have PSDs engineered based on the instrument response, the additional gain from using robust observables will depend on how well the error modes are suppressed in hardware, as well as the timescales at which power in those modes leaks through.To some extent, robust observables will remain applicable to such future telescopes and instruments in the spatial and temporal sub-spaces in which they do not meet their requirements.

Model Accuracy
In this analysis, the model used to generate the instrument response matrices is exactly the same model that is used to generate the synthetic data.In a real observations, the instrument model will not exactly match the behavior of the actual instrument, and one future avenue to explore is how well a model must match the instrument in order for robust observables to work on real data.This technique's robustness can be investigated by first calculating the response matrices using one model, then changing the parameters of the model (e.g. the coronagraphic mask size and displacement, the DM alignment, the detector pixel scale) before generating synthetic data, and examining how well the robust observables work in the presence of model mismatch.
For instruments equipped with wavefront modulating devices such as deformable mirrors, however, the instrument response matrix may also be calculated experimentally.If a perturbation within the linear regime is applied, the difference in measured intensity can be directly registered into the appropriate column of the linear response matrix.The technique for experimentally building the quadratic response matrices is equivalent to the approach used for PASTIS (Laginja et al. 2019), with the difference that the measurements are not averaged over a dark hole, but rather maintained for every pixel.Additionally, some wavefront and control schemes such as implicit electric-field conjugation (Haffert et al. 2023) already involve an empirical measurement of the instrument response, which can be used to derive linear-regime robust observables without having to set aside additional calibration time.Experimentally building instrument response matrices circumvents the need to have a well-matched numerical model, and allows for the response matrices to capture effects in the real instrument.

CONCLUSIONS
A coronagraph model with linear and quadratic contributions of wavefront error to detector plane intensity is developed, and when either term is dominant, the coronagraph response can be approximated by a transfer matrix.A useful projection can be found from this transfer matrix that removes the dominant error modes, resulting in observables that are more robust to WFE in the regime of interest.These robust observables are extracted from synthetically generated data with the Hybrid Lyot Coronagraph of the Roman Space Telescope in both the linear and quadratic regimes.The performance of the robust observables is compared to that of the raw intensity data through calculations of their respective binary companion flux ratio detection limits.In these examples, using the robust observables significantly increases the sensitivity to the signal of a binary companion.A projection onto a robust subspace can in theory be combined with other families of post-processing algorithms.Hybrid post-processing approaches would incorporate information on the instrument response alongside the other available information (such as angular diversity, spectral diversity, reference observations, or WFC telemetry) to fully maximize the sensitivity to astrophysical signals in coronagraphic data; however, the approach outlined in this work can be applied to observational data and result in post-processing gains even if such additional information is unavailable.or spikes indicate that a projection matrix with the specified cutoff mode performs worse than using the raw intensity, which occurs when the modes the majority of the planet signal overlaps with have also been projected out.The full matrix achieves the best results with a cutoff mode of 70, leading to a detection limit of 5.4 × 10 −7 while the approximate matrix achieves the best results with a cutoff mode of 2,727, also leading to a detection limit of 5.4 × 10 −7 .Unshowable in log-log scale is the detection limit with Nm = 0, which, with observables, is 1.4 × 10 −6 .This is, as expected from the fact that no error modes are removed, the same as the raw intensity detection limit of 1.4 × 10 −6 .These results indicate that the approximation performs as well as the full model.
angles provide an indication of how much the subspace spanned by each of these projection matrices overlaps with the subspace spanned by the dominant modes identified by the full model.The results are shown in Figure 9.
The number of principle angles with value 0 is the dimension of overlap between the subspaces.As expected, the subspace angles between P and itself are all 0, meaning it overlaps completely with itself.Also as expected, the angles between K and P are all π/2, as K is orthogonal to P. Both of the random matrices have a random distribution of angles with P centered around π/4.Meanwhile, P ′ (the space of dominant modes derived from the approximate model) overlaps with P (the space of dominant modes derived from the full model) significantly more than random.Crucially, K ′ (the robust subspace from the approximate model) overlaps with P significantly less than random, and with the subspace orthogonal to P's significantly more than random, which is why data projected onto this subspace is still robust to wavefront error.This result shows why the approximate model, despite poorly predicting the detector intensity response, is nevertheless useful for identifying a subspace that overlaps significantly with the robust subspace of the full model.
We can thus use this approximation with all 528 Zernikes in our model to analyze spatial separations beyond the ∼ 5 λ/D spanned by the first 100 Zernikes.To demonstrate this, we build A ′ q528 according to Equation A1, and perform detection tests at a separation of 6.5 λ/D.The results are shown in Figure 10.
Our tests show that the approximation A ′ q528 can successfully increase signal-to-noise at spatial separations beyond the original regime of validity of A q .Thus, even though the input dimension of the quadratic model scales cumbersomely with the number of basis vectors, an approximation considering only norm-squared terms can still be used .The subspace angles between various projection matrices (onto dominant modes on the right, onto a robust subspace on the left) and P, the projection onto the dominant error modes determined from the full quadratic model.The number of principle angles with value 0 is the dimension of overlap between the subspaces.Angles with value π/2 indicate overlap with the subspace orthogonal to P. P ′ (the space of dominant modes derived from the approximate model) overlaps with P (the space of dominant modes derived from the full model) significantly more than random.Crucially, K ′ (the robust subspace from the approximate model) overlaps with the subspace orthogonal to P's significantly more than random, which is why data projected onto this subspace is still robust to wavefront error.
to find observables that are robust to quadratic wavefront error, and thus provide detection gains at farther spatial separations of interest.

Figure 1 .
Figure1.The CGI optical train and wavefront sensing and control architecture.The optical elements of the HLC mode of interest are depicted in the top panel.Before an observation, the high order wavefront sensing and control loop is performed on a bright reference star to generate a 'dark hole' (an area where starlight is suppressed).Then, the DM shapes are fixed, and the telescope slews to the target star for the observation.During the observation, wavefront errors accrue as a result of instrumental disturbances and drifts, the effects of which this work aims to mitigate in post-processing.Figure fromKasdin et al. (2020).

Figure 2 .
Figure 2. The singular values of A l (left) and Aq (right).Note that the transfer matrices are rectangular and have Npix = 5476 total singular modes, but the singular values beyond the size of the input dimension are all 0.

Figure 4
Figure 4. 1) Initial dark hole intensity achieved using electric field conjugation with the HLC.2) A single snapshot with linear-regime wavefront aberrations.3) The same snapshot with an injected companion with a flux ratio of 2 × 10 −7 at 6.5 λ/D (indicated with red circle).4) A single snapshot with quadratic-regime wavefront aberrations.5) The same snapshot with an injected companion with a flux ratio of 5 × 10 −6 at 4 λ/D (indicated with red circle).All intensities are shown in log 10 of raw contrast.

Figure 5
Figure 5. a. Example linear-regime intensity error from the optical model.b.Corresponding linear-regime intensity error reconstructed by response matrix A l , plotted on the same scale as (a).c.The difference between the response matrix prediction and the optical model prediction, plotted on the same scale as (a) and (b).d.Example quadratic-regime intensity error from the optical model.e. Corresponding quadratic-regime intensity error reconstructed by the response matrix Aq, plotted on the same scale as (d).f.The difference between the response matrix prediction and the optical model prediction, plotted on the same scale as (d) and (e).Slight differences arise because the model includes both the linear and quadratic error terms while the matrix predictions only include one or the other, i.e. the linear matrix prediction neglects the contribution of the quadratic term and the quadratic matrix prediction neglects the contribution of the linear term (as well as the influence of any Zernikes past the first 100).While the linear matrix prediction is biased low near the peaks and the quadratic matrix prediction biased high overall, our method depends only on how well the spatial structure of the errors are reproduced.A relevant metric for characterizing the spatial overlap is the normalized inner product between the optical model prediction and the transfer matrix prediction, where a value of 1 indicates perfect spatial overlap and a value of 0 indicates perfect spatial orthogonality.In this case, the normalized inner product is 0.936 for the linear regime example and 0.985 for the quadratic regime example, sufficient for providing a quantifiable improvement in detection sensitivity.

Figure 6 .
Figure6.Detection test results for the quadratic regime noise model.The companion planet considered has flux ratio of 5.4 × 10 −7 and is located at 4.0 λ/D.Left: Histograms from using raw intensities compared to those from using quadratic robust observables with the optimal cutoff of Nm = 70.The histograms using raw intensity overlap significantly, making it difficult to distinguish between a model with a planet and a model without one, while the histograms using the robust observables are further separated and more distinguishable.Middle: Histograms for using raw intensities and a random projection matrix of the same size as the instrumentally-motivated projection.Both sets of histograms overlap significantly, and the random projection does not improve the distinguishability of the two models.Right: ROC curves corresponding to the histograms.Grey area indicates false positive rates which are not well sampled as they involve less than 3 datasets with false detections.The ROC curve shows that while using the robust observables results in a FPR = 0.01 and TPR = 0.9 detection of the injected planet, both the raw intensity and the randomly dimensionally reduced data remain very far from detectability.

Figure 8 .
Figure8.Quadratic regime flux ratio detection limits (FPR = 0.01, TPR = 0.90), to two significant figures, as a function of cutoff mode, for a companion at 4.0 λ/D.Only the first 100 Zernikes are used in the model used to calculated the full and approximate quadratic transfer matrices, but WFE up to 538 Zernikes are included in the synthetic data.Upward triangles or spikes indicate that a projection matrix with the specified cutoff mode performs worse than using the raw intensity, which occurs when the modes the majority of the planet signal overlaps with have also been projected out.The full matrix achieves the best results with a cutoff mode of 70, leading to a detection limit of 5.4 × 10 −7 while the approximate matrix achieves the best results with a cutoff mode of 2,727, also leading to a detection limit of 5.4 × 10 −7 .Unshowable in log-log scale is the detection limit with Nm = 0, which, with observables, is 1.4 × 10 −6 .This is, as expected from the fact that no error modes are removed, the same as the raw intensity detection limit of 1.4 × 10 −6 .These results indicate that the approximation performs as well as the full model.
Figure9.The subspace angles between various projection matrices (onto dominant modes on the right, onto a robust subspace on the left) and P, the projection onto the dominant error modes determined from the full quadratic model.The number of principle angles with value 0 is the dimension of overlap between the subspaces.Angles with value π/2 indicate overlap with the subspace orthogonal to P. P ′ (the space of dominant modes derived from the approximate model) overlaps with P (the space of dominant modes derived from the full model) significantly more than random.Crucially, K ′ (the robust subspace from the approximate model) overlaps with the subspace orthogonal to P's significantly more than random, which is why data projected onto this subspace is still robust to wavefront error.

Table 1 .
Quantities and Dimensions for Analysis of the Roman Space Telescope HLC