Robust Detrending of Spatially Correlated Systematics in Kepler Light Curves Using Low-rank Methods Jamila S. Taaki1 , Athol J. Kemball2,3 , and Farzad Kamalabadi1 1 Department of Electrical and Computer Engineering, University of Illinois at Urbana-Champaign 306 N. Wright Street MC 702, Urbana, IL 61801-2918, USA; xiaziyna@gmail.com 2 Department of Astronomy, University of Illinois at Urbana-Champaign 1002 W. Green Street, Urbana, IL 61801-3074, USA 3 School of Physics, University of the Witwatersrand, P.O. Box Wits, Johannesburg, South Africa Received 2023 July 22; revised 2023 November 13; accepted 2023 November 28; published 2024 January 15 Abstract Light curves produced by wide-field exoplanet transit surveys such as CoRoT, Kepler, and the Transiting Exoplanet Survey Satellite are affected by sensor-wide systematic noise, which is correlated both spatiotemporally and with other instrumental parameters such as the photometric magnitude. Robust and effective systematics mitigation is necessary to achieve the level of photometric accuracy required to detect exoplanet transits and to faithfully recover other forms of intrinsic astrophysical variability. We demonstrate the feasibility of a new exploratory algorithm to remove spatially correlated systematic noise and detrend light curves obtained from wide- field transit surveys. This spatial systematics algorithm is data-driven and fits a low-rank linear model for the systematics conditioned on a total-variation spatial constraint. The total-variation constraint models spatial systematic structure across the sensor on a foundational level. The fit is performed using gradient descent applied to, a variable reduced least-squares penalty and a modified form of total-variation prior; both the systematics basis vectors and their weighting coefficients are iteratively varied. The algorithm was numerically evaluated against a reference principal component analysis, using both signal injection on a selected Kepler dataset, as well as full simulations within the same Kepler coordinate framework. We develop our algorithm to reduce the overfitting of astrophysical variability over longer signal timescales (days) while performing comparably relative to the reference method for exoplanet transit timescales. The algorithm performance and application are assessed, and future development is outlined. Unified Astronomy Thesaurus concepts: Exoplanet detection methods (489); Astronomy data analysis (1858); Wide-field telescopes (1800); Surveys (1671); Transit photometry (1709) 1. Introduction Transit surveys have contributed significantly to exoplanet science over the past two decades, uncovering exoplanet population statistics at current survey completeness limits (Batalha 2014; Bryson et al. 2020). Space-based surveys, including CoRoT (Auvergne et al. 2009), Kepler (Borucki et al. 2010), and the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) have led to the discovery of over 54004 exoplanets to date. These wide-field surveys include a large number of candidate systems and require a photometric precision after standard calibration and processing (residual photometric precision) adequate to detect exoplanet transits (Deeg & Alonso 2018). Achieving low residual photometric noise requires the optimization of both the instrument design and the associated algorithms for instrumental noise removal. In this paper, we describe a method for addressing spatially correlated systematic noise across a wide-field imaging sensor used for exoplanet transit detection. If an exoplanet transit induces a fractional flux density variation !F in observations, with a differential photometric precision σp (after processing) then the signal-to-noise ratio is ~ s / S N F p (Deeg & Alonso 2018). For example, an Earth–Sun transit requires σp∼ 80 parts-per-million (ppm) over the transit duration T ∼ hours (Caldwell et al. 2010). The Kepler telescope achieved σp∼ 30 ppm over T= 6.5hr for stars with Kepler magnitude Kp∼ 12 (Koch et al. 2010; Gilliland et al. 2011; Christiansen et al. 2012). TESS achieved σp∼ 230 ppm over T ∼ 1 hr for a star with TESS magnitude Tp∼ 10, which is sufficient to detect super-Earths around bright stars (Fausnaugh 2018). For reference, the Hubble Space Telescopes Space Telescope Imaging Spectrograph (HST/STIS) can attain σp∼ 120 ppm over T= 45 minutes (Demory et al. 2015), and non-wide-field ground-based telescopes can reach a precision sufficient to detect large Jovian exoplanets (Tregloan-Reed & Southworth 2013; Stefansson et al. 2017). The instrumental response will vary over a range of timescales and across various spatial scales producing systematic noise effects, as expected from general principles (McLean 2008). The primary data products for the Kepler and TESS telescopes include light curves derived using optimized simple aperture photometry (SAP) or image data providing pixel-level light curves (Jenkins 2017; Tenenbaum & Jenkins 2018). Algorithmic innovations to suppress the residual differential photometric precision σp, arising from unmodeled systematic noise effects, are therefore critical to detect weak exoplanet signals. These systematic effects have diverse physical origins (Jenkins 2017; Tenenbaum & Jenkins 2018). A single pixel exhibits a nonuni- form sensitivity across its surface (Toyozumi & Ashley 2005; Hedges et al. 2021). In addition, there are variations in pixel-to- pixel sensitivity and instrumental pixel response (Van Cleve & Caldwell 2016) including CCD pattern, read, photon, and The Astronomical Journal, 167:60 (29pp), 2024 February https://doi.org/10.3847/1538-3881/ad1110 © 2024. The Author(s). Published by the American Astronomical Society. Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 4 https://exoplanetarchive.ipac.caltech.edu/ 1 https://orcid.org/0000-0001-5475-1975 https://orcid.org/0000-0001-5475-1975 https://orcid.org/0000-0001-5475-1975 https://orcid.org/0000-0001-6233-8347 https://orcid.org/0000-0001-6233-8347 https://orcid.org/0000-0001-6233-8347 mailto:xiaziyna@gmail.com http://astrothesaurus.org/uat/489 http://astrothesaurus.org/uat/1858 http://astrothesaurus.org/uat/1800 http://astrothesaurus.org/uat/1671 http://astrothesaurus.org/uat/1709 https://doi.org/10.3847/1538-3881/ad1110 https://crossmark.crossref.org/dialog/?doi=10.3847/1538-3881/ad1110&domain=pdf&date_stamp=2024-01-15 https://crossmark.crossref.org/dialog/?doi=10.3847/1538-3881/ad1110&domain=pdf&date_stamp=2024-01-15 http://creativecommons.org/licenses/by/4.0/ https://exoplanetarchive.ipac.caltech.edu/ https://exoplanetarchive.ipac.caltech.edu/ https://exoplanetarchive.ipac.caltech.edu/ quantization noise (McLean 2008; Caldwell et al. 2010; Gilliland et al. 2011); instrumental CCD terms are expected to vary by CCD output channel or module. A net pixel response function (PRF; Bryson et al. 2010) captures the spatially variant response across the detector including the spatial variation of the optical PSF and instrumental detector sensitivity, among related factors (Jenkins 2017). The spatial response of the detector may, however, include unmodeled error due to pointing jitter, focus changes, uncorrected differential velocity aberration terms, and their interaction (Van Cleve & Caldwell 2016). Thermal effects as well as discrete spacecraft guidance and downlink control events, including momentum dumps, may produce temporal and spatial variation in the systematics, some abrupt (Van Cleve & Caldwell 2016; Vanderspek 2018). Further, bright astrophysical sources may produce pixel saturation or broadening of the PRF (Van Cleve & Caldwell 2016; Vanderspek 2018). The first correction for noise effects is typically applied in early robust calibration pipelines (Jenkins 2017; Tenenbaum & Jenkins 2018). However residual spatially variant and time- variant systematics unavoidably remain in the calibrated data products given the complex instrumental response at the current level of photometric accuracy. We denote these residual systematics contributions to each light curve i at discrete time sample n as li[n] (Section 2.2) where the data may be either pixel or SAP light curves. The systematics li[n] has the same discrete time sampling as the measured light curve. The residual systematics vary spatially over the sensor, reflected in equivalent notation lx,y[n], where light curve i maps to sensor position (x, y). Petigura & Marcy (2012) show an example of spatially varying residual systematic noise in Kepler data by characterizing correlations across the sensor. Moreno et al. (2021) similarly characterize spatial systematics correlation in Kepler/K2 data. An a priori analytic instrumental model for lx,y[n] is intractable in practice, and data-driven approaches, informed by physically realistic instrumental assumptions, provide the most effective approach to mitigate the residual systematic noise in the target light curves. There is a rich history in the literature on this subject, which we review briefly only to place our method in context. A class of methods (hereinafter external parameter decorr- elation methods; Bakos et al. 2007) model the functional form of the residual systematics lx,y[n] as correlates of other explanatory variables such as pointing error estimates or ancillary engineering data (PDC-LS, Twicken et al. 2010). Charbonneau et al. (2005) demonstrated Spitzer5 aperture photometry correction using decorrelation with pointing errors estimated using target centroiding. There has been substantial development in this latter area, especially inspired by the pointing challenges of the Kepler K2 mission and pixel-level detrending, including the work described by Vanderburg & Johnson (2014), Huang et al. (2015), Aigrain et al. (2015), Lund et al. (2015) and Crossfield et al. (2015). External parameter decorrelation methods implicitly include spatiotem- poral variability in the residual systematics via the encoding of their functional form in the proxy variables. A second class of methods (hereinafter cotrending methods) makes the foundational assumption that a set of measured light curves, or their computed vector basis, comprise an efficient basis of regressors over which to expand the unknown functional form of lx,y[n] as a low-rank linear model. This implicitly assumes that the residual systematics are temporally and spatially correlated across the sensor; this is similarly a physically reasonable assumption. The trend-filtering algorithm (TFA; Kovács et al. 2005) directly uses a uniform selection of external light curves as regressors in a least-squares fit to detrend an individual target light curve. The cotrending approach is analogous to the joint coupled least-squares solution for the product of stellar extinction coefficients and airmasses as posed in a global solution for wide-field multiobject photometry (Sysrem; Tamuz et al. 2005); this solution is equivalent to a principal component analysis (PCA) of target observations. A refinement including additive external parameters is provided by SARS (Ofir et al. 2010). Cotrending methods have seen significant algorithmic refinement. The PDC-MAP algorithm (Smith et al. 2012; Stumpe et al. 2012) constructs a set of cotrending basis vectors (CBV) using singular value decomposition (SVD) on a template set of light curves selected for their high degree of correlation (therefore capturing foundational instrumental trends) and quiescence. Each CBV is fit against the target light curve using Bayesian maximum a posteriori (MAP) methods to obtain relative coefficient weightings for each basis term, the sum of which is then subtracted from the target time series to remove systematic effects. An empirical Bayesian prior is used to constrain overfitting, constructed on the variation of coefficient value over stellar magnitude and spatial position on the sensor. The CBV are orthogonal by mathematical construction and, therefore, do not map directly to constituent instrumental effects; this dilutes the variable dependence in the empirical coefficient prior. Instrumental effects are often separated by characteristic timescale and this mapping can be improved (MS-MAP; Stumpe et al. 2014) by deriving separate CBV for different timescale wavelet sub- bands; this provides better temporal separation in the spatiotemporal dependence of lx,y[n]. A net composite systematic correction is applied as the multiscale combination across each subband. The Astrophysically Robust Correction algorithm (ARC; Roberts et al. 2013; Aigrain et al. 2017) is related conceptually to PDC-MAP but uses automatic rele- vance coefficient priors to maximize model evidence for instrumental trends. Additional cotrending methods include the Causal Pixel Model (CPM; Wang et al. 2016), which uses a template (or training) set of regressor pixel-data light curves that are separated in time from the target light-curve sample being corrected for systematics. The training set is spatially separated from the target but the target light curve itself is included to include autoregressive information. Both measures suppress residual variability except on short transit timescales. A successor method is described by Hattori et al. (2022). The method described by Foreman-mackey et al. (2015) uses a large training set to derive an expanded CBV but uses a joint estimation of systematics and the transit signal in order to better constrain the problem. A related Bayesian formulation is described by Taaki et al. (2020). The pixel-level decorrelation algorithm (PLD; Deming et al. 2015) was developed to implicitly correct photometric decorrelation with pointing errors in Spitzer data. PLD uses pixel light-curve regressors within a summed aperture to model SAP light curves. The regressors encode pointing errors without the need for explicit centroiding as in external parameter decorrelation approaches. The PLD method was extended by Luger et al. (2016) for application to Kepler K2 data by including higher-order terms in the SAP flux dependence on pointing error, by modeling astrophysical variability using a5 https://www.spitzer.caltech.edu/ 2 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi https://www.spitzer.caltech.edu/ Gaussian process, and by constructing a PCA basis from the regressors to better constrain the model. The method was further improved (Luger et al. 2018) by including L2 coefficient regularization, which was also used in CPM, and by incorporating PLD vectors from nearby stars (nPLD) to improve sensitivity to pointing error for faint stars. The cotrending methods described above address challenges in estimating intrinsic astrophysical variability. Foundationally it is difficult to perfectly separate systematics and astrophysical variability a priori in a light-curve regressor or derived basis set. Therefore the linear systematics model may include and absorb true astrophysical variability (overfitting) or inject spurious variability into detrended light curves. There may also be an incidental correlation between a relatively clean systematics basis and astrophysical variability (Smith et al. 2018). Overfitting can be reduced by judicious use of prior constraints that represent a physically realistic instrumental response; these constraints may also significantly improve the condition of the problem. In the cotrending methods described above, spatial structure is incorporated into the systematic noise model using several different approaches. In the simplest form, the spatial variation constraint is imposed implicitly by cotrending targets across a discrete sensor region as a whole (e.g., a CCD output channel; Smith et al. 2012; Roberts et al. 2013; Stumpe et al. 2014; Luger et al. 2018; Lund et al. 2021). Spatial variation has also been incorporated by restricting regressors by proximity to the target light curve (Luger et al. 2018; Lund et al. 2021), or by directly introducing a parametric dependence of fitted systema- tics on sensor position in an empirical prior (Smith et al. 2012; Stumpe et al. 2014). In this paper, we present an exploratory algorithm to solve for a cotrending low-rank linear systematics model while incorporating a spatial constraint across the sensor at a foundational level. We refer to this as the spatial systematics method in what follows. The spatial constraint is of generalized total variation (TV) from Rudin et al. (1992). This TV constraint on basis vector coefficients promotes the correlation of adjacent neighbors on the sensor and is weighted spatially across the sensor while permitting discontinuities as found in wide-field imaging sensors at module or channel edges. As such, the spatial constraint has a realistic physical motivation and therefore a strong potential to reduce overfitting. The fit the overall SAP light curves is performed by minimizing the sum of the least-squares residual between the light curves and the linear systematics model, and, the total-variation spatial constraint. We use variable elimination (Golub & Pereyra 1973) to reformulate the optimization problem as a function of the weighted coefficients only, introducing stability to the mini- mization (Golub & Pereyra 2003; Shearer & Gilbert 2013). An approximate closed-form gradient of the objective function is derived and minimized via gradient descent. In this iterative solution both the weighted coefficients and the basis vectors, which functionally are defined by the coefficients, are varied. The algorithm as implemented is available as a public Python package on the GitHub repository6 or via PyPI under the package name spatial-detrend.7 We numerically evaluate the performance of our method in the context of the Kepler sensor and against PCA as a reference method. We use both injected signals in real long-cadence (LC) Kepler data and fully simulated data in this evaluation. As noted above (Smith et al. 2012, PDC-MAP) Kepler systematics depend on stellar magnitude. This is also found for CoRoT light curves (Mazeh et al. 2009; Ofir et al. 2010). This effect likely arises due to pixel saturation effects. Kepler operated in the magnitude band Kp ä 9 to 15 (Koch et al. 2010). In this work for simplicity we do not include a prior on stellar magnitude dependence but instead select light curves from a narrow fixed magnitude band Kp ä 12 to 13 within which magnitude-dependent systematics such as saturation are not expected to substantially vary (Smith et al. 2012); this therefore allows for the spatial variability systematics to be isolated. Our numerical evaluation demonstrates that the spatial systematics algorithm has reduced overfitting of astrophysical variability, particularly on day-long timescales, compared to the reference PCA method. The paper is organized as follows. Section 2 introduces the light-curve signal model, our method for spatial systematics inference, and an analysis of spatial correlation across the Kepler sensor, and describes our experimental tests to numerically evaluate the algorithm. In Section 3 we report the results of our numerical evaluation. These results are discussed in Section 4 and conclusions are presented in Section 5. A table of commonly used symbols is provided in Appendix A for reference. 2. Methods As described above, the data products for wide-field exoplanet transit surveys may comprise: (i) pixel-level image data either within target apertures or for full frames; and, (ii) SAP target light curves (Jenkins 2017; Tenenbaum & Jenkins 2018). Here we use Kepler SAP target light curves with a sampling cadence of 30 minutes produced by the Kepler science data processing pipeline but not including post- processing by the Kepler Pre-Search Data Conditioning (PDC) module (Twicken et al. 2010; Smith et al. 2012; Stumpe et al. 2012); the latter light curves with full Kepler post- processing are distinguished as PDC-SAP light curves. We present the details of the spatial systematics method that is the subject of this paper. 2.1. Matrix Notation We define common notation here that is used throughout the paper. The ith row of matrix M is denoted as [ ]M i and the jth column as [ ]·M j, . The matrix element corresponding to the ith row and the jth column is denoted as [ ]M i j, or Mi,j. The transpose of a matrixM is denoted asMT. An identity matrix of size N×N is denoted as 1N. The Lp norm of a vector x of length I is defined as: = å =( ∣ ∣ ) x xp i I i p p 1 1 . If Î ´M I J , we use ∥M∥p,q to denote the Lp norm applied to each column [ ]·M j, , followed by the Lq norm applied to this vector, such that = å å= =( ( ∣ ∣ ) ) M Mp q j J i I i j p q p q , 1 1 , 1 . The Frobenius norm is defined as =   M MF 2,2 1 , such that = å å= = ∣ ∣ M MF j J i I i j 2 1 1 , 2 (Golub & Van Loan 2013). 2.2. Light-curve Decomposition A collection of target light curves {yi: iä I} are obtained on a sensor, each light curve yi (over target index set I) is a length-N time series. Each light curve is represented as the sum of a 6 https://github.com/xiaziyna/spatial-detrend 7 https://pypi.org/project/spatial-detrend 3 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi https://github.com/xiaziyna/spatial-detrend https://pypi.org/project/spatial-detrend systematics term li and a statistical noise term ni: = + ( )y l n . 1i i i In matrix form the data model takes the form Y= L+N, where yi, li, and ni form the columns of matrices Î ´Y N I , Î ´L N I , and Î ´N N I representing the light curves, systematic noise, and statistical noise, respectively. The goal of this work is to form an accurate estimate of the systematic noise li. Sources of systematic errors are described in Section 1; see also Auvergne et al. (2009), Van Cleve & Caldwell (2016), and Vanderspek (2018). Statistical noise originates due to a mixture of instrumental error and astrophysical variability. We approximate statistical noise ni as white Gaussian noise  s~ ( )n 0,i i 2 as additive statistical noise sources may reasonably approach  s( )0, i 2 under the central limit theorem (Grinstead & Snell 2012). This is a common assumption in this domain and implicit in least-squares minimization approaches (Smith et al. 2012; Stumpe et al. 2012, 2014; Aigrain et al. 2017). The noise level σi is unknown a priori, but may reasonably be estimated from filtered and coarsely detrended light curves. In what follows, light curves are normalized by an estimate of σi so that ~ ( )n 0, 1i for mathematical convenience. Light curves may include intermittent sources of systematic noise that produce outliers and are difficult to model. Their influence may be minimized by prior constraints (a further motivation for our approach here), but data-driven outlier filtering is generally necessary in this domain (Section 2.5). 2.3. Systematics Model Low-rank model: We adopt a cotrending basis model: å= = [ ] [ ] ( )l vn c n , 2i k K i k k 1 where (K= N) describes the rank of the systematic noise model, {vk: k ä K} are a set of basis vectors shared by light curves, and each ci k is a coefficient weighting of vk for light curve i. The coefficient vector for light curve i is = [ ]c c c c, ,...i i i i k T1 2 . In matrix form: = ( )L VC, 3 where Î ´L N I is defined above, the columns of Î ´V N K are the basis vectors {vk: k äK}, and the columns of Î ´C K I are the coefficients {ci: iä I}. Graphically, ¼ = ¼ ¼ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ( )l l l v v v c c c . 4I K I1 2 1 2 1 2 We define a column-normalized coefficient matrix C̄ of the form =¯   c c ci i i 2 : = ¼ ⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ¯ ∣ ∣ ∣ ∣ ∣ ∣ ( )     C . 5 c c c c c c I I 1 1 2 2 2 2 2 The rank of L may be less than the number of independent systematic noise sources. As ( ) VCrank { ( ) ( )}V Cmin rank , rank , it is possible to have a more expansive basis set V representing many individual noise effects, but for the relative weightings C between light curves to have a degenerate and therefore low-rank structure and consequently for L to be low rank. SVD/PCA systematics estimation: SVD and PCA are integral to many cotrending approaches (Tamuz et al. 2005; Thatte et al. 2010; Petigura & Marcy 2012; Smith et al. 2012; Stumpe et al. 2012) and are summarized in Appendix B. Under a white noise model for N the maximum likelihood (ML) estimate of a matrix L (not necessarily low-rank) given light curves Y=N+ L is equivalent to minimizing the least-squares residual = -( ) ∣∣ ∣∣V C Y Lf , F 2 : º -( ∣ ) ∣∣ ∣∣ ( )Y L Y Lparg max arg min , 6 L L F 2 where p(Y|L) is the likelihood of Y given L (Srebro 2004). For the data model in Equation (1), without any further constraints, a rank-K optimal least-squares solution of L can be found by rank-thresholding the SVD or PCA decomposition of Y per the Eckart–Young–Mirsky theorem (Appendix B) and in this case takes the form: L= VKCK. In what follows we use the SVD/PCA method primarily as a comparative method (Section 2.5). Spatial systematics model: As described in Section 1, the residual systematics are spatially correlated across the sensor. In addition, the white noise model for N is an idealization and the systematics cannot then be perfectly separated by rank thresholding alone. In our spatial systematics algorithm, we retain a low-rank systematics data model L while recognizing the incomplete model for N. We introduce a spatial side constraint of total-variation measure (Appendix C) on the normalized coefficient matrix C̄ to facilitate a more robust separation of the unknown astrophysical signals. The total- variation constraint admits a small number of bounded discontinuities within generally smooth functions (Rudin et al. 1992; Vogel 2002) and promotes correlation between neighboring light-curve systematics. This choice of spatial constraint is experimentally motivated in Section 2.5. In this formulation, the estimate for L= VC is obtained by minimizing an objective function comprising the least-square residual = -( ) ∣∣ ∣∣V C Y VCf , F 2 and the total-variation penalty con- straint =( ) ¯ C D Cg W p p 2, with p ä [1, 2] (Appendix C): + = - + { ( ) ( )} {∣∣ ∣∣ ∣∣ ¯ ∣∣ } ( ) ( ) ( )   V C C Y VC D C f garg min , arg min . 7 V C VC V C VC W K K F p p , : rank , : rank 2 2, Here DW=WD where D is a difference operator, as defined in Appendix C, and W is a diagonal matrix of weights wx, y between [0, 1] to model nonuniform spatial correlation across the sensor. We assume each light curve i has a pixel position i→ (x, y) ä (X, Y), where Î +X Y, . Further, we assume uniform spatial cells of size (Δx, Δy) across the sensor and that there is a one-to-one mapping between light curve i and spatial cell. The weighted difference operator DW applied to C̄ computes for each k äK and cell (x, y), the weighted difference of neighboring coefficients - - Î+ +[ ¯ ¯ ¯ ¯ ] c c c c,x y k x y k x y k x y k , , 1 , 1, 2 as the columns of Î ´¯ · ·D CW K X Y2 and where = ¯ ¯ ( )c cx y k i x y k , , . This is depicted graphically in Figure 1. In expanded form ¯ D CW p p 2, = ∑käK∑(x,y) - + -+ +(∣ ¯ ¯ ∣ ∣ ¯ ¯ ∣ )w c c c cx y x y k x y k x y k x y k, , , 1 2 , 1, 2 p 2 . 4 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi When p= 2 minimizing the total-variation spatial constraint is equivalent to maximizing the correlation between neighbor- ing coefficient vectors; this is shown in Appendix C.1. Placing the spatial prior on the normalized coefficients C̄ instead of the full systematics model L= VC has mathematical advantages as described in Section 2.4 and is shown in Appendix E to promote correlation between neighboring systematics estimates by proxy. The prior is calculated using normalized coefficients C̄ to avoid the objective being minimized to a trivial minimum. If instead the prior was of the form  D CW p p 2, any solution C can be replaced with a solution a¢ =C C for α< 1 that produces a lower value of the spatial prior a¢ =   D C D CW Wp p p p 2, 2, . Taking ¢ = a V V , there is no change to the least-squares penalty - = - ¢ ¢   Y VC Y V C ; therefore a trivial minimum can be found as α→ 0. Basis vectors obtained with SVD and PCA are constrained to be orthogonal; however, orthogonality is not a necessary constraint as a property of a systematics model, nor to achieve a minimal optimization cost. Our model does not restrict V or C to be orthogonal. 2.4. Systematics Inference The systematics L=VC are inferred by minimizing the objective function defined in Equation (7). This function is defined over variables Πδ ´ V C,N K K I , and com- prises the least-squares residual = -( )  V C Y VCf , F 2 , which depends on both variables V, C, and the total-variation penalty constraint =( ) ¯ C D Cg W p p 2, , which only depends on C. The least-squares residual f (V, C) is not separable over V, C. Coupled least-squares problems of this type can be solved by iterative alternating solutions for V and C (Wold 1973; Tamuz et al. 2005) or by using variable elimination to remove dependence on one of the variables (Golub & Pereyra 2003). There is no evidence that either is preferable and we adopt the latter variable-projection method here. For fixed V or fixed C, the least-squares residual f (V, C) has an analytic solution for the complementary free variable. Therefore we can eliminate variable dependence on V in the least-squares residual so that the objective can be minimized over C alone. For a value of C, the minimizing value V of the least-squares residual is denoted as a function h(C): =( ) ( ) ( )C V Ch farg min , . 8 V This function h(C) can be inserted as V into the least-squares residual so that the original objective depends only on C: +  +{ ( ) ( )} { ( ( ) ) ( )} ( )V C C C C Cf g f h gmin , min , . 9 V C C, Using the theory described in Golub & Pereyra (2003), as elaborated in Appendix F.1, the variable-reduced objective takes the form:  - +{∣∣( ) ∣∣ ∣∣ ¯ ∣∣ } ( ) ( ) †  C C Y D Carg min . 10 C C W K I T T F p p : rank 2 2, The minimum of the derived objective does not have an analytic solution but can be obtained numerically using gradient-descent optimization methods. However, the Lp norm in the total-variation constraint g(C) is not directly differenti- able as it is discontinuous for zero-valued entries of ¯D CW when p= 1. Proximal methods are standard for obtaining the minima of objective functions with nondifferentiable penalties such as L1 norms, as described in Parikh & Boyd (2013). In particular, the Split–Bregman method (Goldstein & Osher 2009) is well suited to total-variation-regularized problems. A preferred approach is to replace the nondifferenti- able Lp norm with a continuously differentiable Huber loss function as an approximation (Vogel 2002). Our general gradient-descent algorithm to minimize the variable-reduced objective is described in Algorithm 1. We refer to the variable- reduced objective function in Equation (10) as w(·), the gradient of w(·) with respect to C is denoted ∇w(·). The gradient of the variable-reduced least-squares term ¢ = = -( ) ( ( ) ) ∣∣( ) ∣∣†C C C C C Yf f h , I T T F 2 in the objective is denoted  ¢(·)f and the gradient of the total-variation constraint term =( ) ∣∣ ¯ ∣∣C D Cg W p p 2, as ∇g(·). These gradients are used to iteratively update Ct at each timestep t. We note that this variable-projection algorithm also updates V implicitly, denoted V|C. The details of the Huber loss function and the derived gradients are described in Appendix F.2. Algorithm 1. Variable-projection Gradient Descent for the Spatial Systematics Algorithm Initialize C for rank K and weight matrix W . The step size at may be fixed or vary at each step t. while -+∣ ( ) ( )∣ C Cw wt t1 do 1)  =  ¢ + ( ) ( ) ( )C C Cw f gt t t 2) a= - + ( )C C Cwt t t t1 end = ∣L V CC with = -( )∣V YC CCC T T 1. 2.5. Algorithm Implementation and Evaluation In this section, we describe further details concerning the implementation and evaluation of the spatial systematics algorithm presented in this paper. As noted above, this method has been developed to explore new algorithmic approaches to exoplanet transit detection. Our prime focus is to demonstrate Figure 1. Visual representation of the coefficient matrix C̄. Each image layer C̄k represents the coefficient weights for a basis term k, ([ ¯ ]C k) organized by 2D sensor cell position. The mapping of light curve i to cell position i → (x, y) is implicit in this notation, as described in the main text. 5 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi the feasibility of the algorithm and to provide an initial assessment of its performance. Accordingly, this work is not designed to support or claim the optimality of this method over existing specialized detrending approaches. The spatial systematics algorithm was evaluated numerically on a selected subset of Kepler test data using both injection tests with time-variable signals (Section 2.5.4) and fully simulated light curves in the same Kepler coordinates (Section 2.5.5). Further, a high-level comparison between the spatial systematics algorithm and several standard detrending methods was performed over the same Kepler test data (Section 2.5.6). The Kepler test data are described in Section 2.5.1. The data analysis framework used for the numerical evaluation tests is depicted in Figure 2. The free parameters to be chosen, shown in our core Algorithm 1 include the initial value of C0, the model rank K, the total-variation norm parameter p, the gradient step size α, and the spatial weighting matrix W. The analysis framework and initial values are discussed in further detail below. The data analysis was performed in Python and using a single CPU core. The computational time complexity of the spatial algorithm acting on a collection of light curves I= X× Y is quadratic in |I| and linear in the number of gradient step iterations, and the model rank K. The leading order time complexity of a gradient step in Algorithm 1 defined in Appendix F.2 is O(|I|2K ). As for difference matrices Î - ´( )Dx Y X YX1 , Î - ´( )Dy Y X YX1 , per gradient step, a matrix multiplication of Î ´D Dx T x I I and Î ´D Dy T y I I with Î ´CT I K is performed. The difference matrices Dx and Dy are sparse in form and for either matrix approximately ∣ ∣I 2 fraction of elements are nonzero. In our implementation, matrices Dx and Dy are generated and stored in compressed- sparse-row (CSR) format with Scipy sparse. 2.5.1. Kepler Test Data The test data were selected from the Kepler long-cadence SAP light curves; these have an integration time of 29.4 minutes (Jenkins et al. 2010). The mapping of a stellar target to approximate sensor position repeats every four Kepler quarters (Aigrain et al. 2017) and accordingly to investigate the consistency of spatial trends over several quarters, we selected quarters (Q hereinafter) Q6, Q10, and Q14 separated at this cadence. Quarter Q6 was chosen to match Petigura & Marcy (2012) and allow for intercomparison with their prior work on spatial correlation; the decision to include successor quarters (Q10, Q14) rather than preceding quarters (Q2) was also made due to calibration issues during early Kepler operations (Van Cleve 2010). Stellar targets were selected within a range of Kepler magnitude Kpä [12, 13] and with a combined differential photometric precision over 12 hr CDPP12h� 40 (Christiansen et al. 2012),8 resulting in 6179, 6286, and 6049 light curves for Q6, Q10, and Q14, respectively. These target selection criteria were informed by Petigura & Marcy (2012); however, we do not believe that our numerical algorithm evaluation is highly sensitive to these exact parameter choices. The narrow range in photometric magnitude was used to isolate the known dependence of systematics on magnitude (Stumpe et al. 2014) and to allow us a focus on spatial systematics in our algorithm. Preprocessing: As noted in Section 2.2 the assumption of Gaussian statistical noise N is an idealization and preprocessing is necessary to address several data conditioning issues, as described here. First, any missing data values, and four visually identified outlier cadences, were substituted using linear interpolation between the preceding and following data samples spanning each gap. Median normalization was then applied to each light curve as = -ˆ ( )y y ymed 1, where med (y) is the median of y, and the subscript on light-curve index i (Equation (1)) is dropped for clarity. Here, each preprocessing step is sequential with generic input y and output ŷ. We then removed a first-order linear trend from each light curve as systematics are known to depend on timescale (Stumpe et al. 2014). This is evident in Figure 3, where there is no clear correlation between short- and long-timescale effects. Remov- ing the first-order trend allows us a focus on spatial systematics in our algorithm. While it is possible to apply our algorithm to each timescale bandwidth separately in a decoupled manner that generalization is beyond the scope of the current paper. After linear detrending, the first-order difference variance s = - "-({ })z y y nvar :z n n 2 1 and variance σ2= var(y) were computed for each light curve, where n ä N here denotes the time sample index in an individual light curve. The light curves were then filtered to retain only those that are in the lowest 90% in both sz 2 and σ2, ranked separately. In a second filtering step, light curves are removed if they do not satisfy a minimum correlation requirement. Each light curve yi must have a correlation coefficient ρ greater than 0.6 with at least 10 other light curves yj: j≠ i, where the correlation coefficient is calculated as r =( )    y y, y y y yi j i T j i j . To identify outlier samples the light curves were then coarsely detrended using exploratory PCA9 (Appendix B) and all samples exceeding a 3σ point-to- point scatter were then flagged. This exploratory PCA detrending was not retained however; it was only used to flag and linearly interpolate outlier samples in the input data at this point in the pre-preprocessing. These initial preprocessing and filtering steps resulted in selecting 4749/6179 (76%), 4850/ 6286 (77%), and 4741/6049 (78%) of light curves for Q6, Q10, and Q14, respectively. The final step in preprocessing concerned mapping the data to a regular spatial grid on the sensor. The Kepler sensor comprises 25 modules, of which 21 contain two CCDs across four output channels (Van Cleve & Caldwell 2016). Each CCD is 2200× 1024 pixels in size and each module contains a total Figure 2. Data analysis framework for numerical evaluation of the spatial systematics algorithm. The data flow includes input data (purple box), preprocessing (gray box; Section 2.5.1), and application of the spatial systematics algorithm (blue box; Sections 2.4 and 2.5). 8 Light curve data were obtained from the MAST archive:archive.stsci.edu/ kepler/data_search/search.php. 9 https://scikit-learn.org/stable/ 6 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi http://archive.stsci.edu/kepler/data_search/search.php http://archive.stsci.edu/kepler/data_search/search.php https://scikit-learn.org/stable/ of 2200× 2048 CCD pixels (Van Cleve & Caldwell 2016). The pixel position of each target light curve on a CCD was mapped to a rescaled square global pixel coordinate grid across the sensor of size 11,000× 11,000 for ease of analysis. In this process each CCD row coordinate was rescaled by a factor 1100/1024 and gaps between CCDs were removed; this resulted in convenient module coordinates of size 2200× 2200. The total-variation constraint depends only on light-curve adjacency, so this rescaling has no impact on our algorithm. The positions of the preprocessed test data light curves on the Kepler sensor are shown in blue in Figure 4. This Figure uses global pixel coordinates and is further labeled by Kepler module number. The difference operators Dx and Dy are easier to implement on a rectangular data grid and light curves on Figure 3. Representative sample of 15 light curves from Q10 shown before (a) and after (b) linear detrending and outlier removal. The x-axis index is in the unit of long-cadence time samples (Δt = 29.4 minutes). The dominant systematics visible are due to the quarterly roll and earth point recoveries (Stumpe et al. 2012). Figure 4. Positions of all (4850/6286) preprocessed Kepler test data light curves for Q10 plotted on the Kepler sensor in global pixel coordinates (11,000×11,000). Kepler modules are demarcated by bold dashed lines and labeled by module number. The (30 × 50) spatial grid is drawn in dashed lines with individual cells of size (220 × 220) pixels. All light curves with coordinates included in the rectangular gridded region are shown in blue. The subset of light curves that were gridded, and therefore processed by the spatial systematics algorithm, are shown in red at their original positions on the CCD. These light curves were either gridded to the (30 × 50) spatial cell in which they fall or to a nearby empty spatial cell, as described in the main text. After gridding, no spatial cell contained multiple light curves and there were no empty spatial cells. Light curves on modules excluded from the data grid are shown in gray. There are no targets in module 3 due to module failure. 7 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi modules (2, 3, 4, 22, 23, 24) were discarded, with light curves on a total of 15 modules retained. However, a rectangular grid is not an algorithmic requirement. As described in Section 2.3 our implementation of the total-variation constraint requires a one-to-one mapping between light curve i and individual cells in a regular spatial grid. We adopted a spatial cell size of 220× 220 pixels for this gridding, comprising (30× 50) cells across the remaining modules with (10× 10) cells per module; this grid is shown in Figure 4. The light curve mapped to each cell was selected randomly from the preprocessed light curves falling within each cell. Where a cell contained no light curves, the light curve closest to the center of the cell was selected; this occurred for 12% of 30× 50 targets. The 1500 gridded light curves are shown in red in Figure 4. Spatial structure on Kepler sensor: The pairwise neighbor correlation of the gridded light-curve sample across the sensor was measured to inform our choice of algorithm parameter W (Equation 7). The correlation between vectors p and q is denoted r =( )    p q, p q p q T . We denote the set of gridded cell positions within each module as x ä XM and y ä YM, respec- tively. The set of gridded light curves within a module is then the Cartesian set product M= XM× YM. The mean pairwise correlation r̄M between all gridded light curves (ym, yn) in a module was computed as: år r= - Î ´ ¹ ¯ ∣ ∣ ∣ ∣ ( ) ( ) ( ) ( ) y y M M 1 , . 11M m n M M m n m n2 , , The mean pairwise neighboring correlation within a module r̄MXY was calculated as: å r r r = - ´ + Î - ´ - + + ¯ ( ) ( ) ( ) ( ) ( ) ( ) y y y y X Y 1 2 2 , , . 12 M M M x y X Y x y x y x y x y , 1 1 , , 1 , 1, M M XY The spatial correlation structure across the gridded (30× 50) region of the Kepler sensor is shown in Figure 5. The mean pairwise correlation r̄M for each module is tabulated in black in this Figure, and the difference r r-¯ ¯M MXY is shown in red. The mean pairwise neighbor correlation r̄MXY is generally higher than the mean pairwise correlation r̄M across the module as a whole. The mean pairwise correlation for all gridded light curves over all modules was computed to be r =¯ 0.26 and is notably lower than the per-module correlations. Moreno et al. (2021) describe how time-dependent systematics, which lags with the radius, may plausibly give rise to spatial correlations as observed in Figure 5. We further explored the spatial structure of the Kepler sensor by performing an exploratory PCA decomposition of the gridded light curves for Q6, Q10, and Q14. The resulting leading PCA basis vector v1 for each quarter is shown in Figure 6 along with a color map of the leading coefficient value ci 1 in each cell. We note the implicit mapping of light-curve index i to cell position i→ (x, y) ä (X, Y) defined in Section 2.3. The leading systematics basis vector is informative of general systematics as it is the strongest term. As shown in Figure 6 the leading coefficient values are highly spatially correlated, with a blocked structure per module, with anti-correlations present between modules. This effect is explained by Petigura & Marcy (2012) as due to PSF variation from the momentum cycle. The overall spatial structure is also persistent across these quarters. As discussed in Section 1 and in Moreno et al. (2021), Petigura & Marcy (2012), and Lund et al. (2021), the systematics exhibit spatial dependence that is key to our choice of a total-variation constraint. 2.5.2. Algorithm Initial Parameter Values As described above, there are several free parameters to be chosen in our method, including the choice of initialization C0, Figure 5. Spatial correlation structure within modules for Q6, Q10, and Q14 across the (30 × 50) gridded region of the Kepler sensor. The module number is shown in the top left corner of each module. The mean pairwise correlation r̄M within each module is tabulated in black. The difference r r-¯ ¯M MXY is tabulated in red, where r̄MXY is the mean pairwise correlation of neighboring light curves within each module. The shaded color intensity is proportional to r̄MXY . In most modules, the correlation between neighboring light curves is greater than the mean pairwise correlation between all pairs of light curves on the module. 8 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi the model rank K, the choice of prior norm p, the gradient step size α, and the prior weighting matrix W. The weighting matrixW was formed as neighboring pairwise correlations of gridded light curves yi scaled by a factor 0.2 as r r= ++ +[ ( ) ( )]y y y yw , ,x y x y x y x y x y , 0.2 2 , 1, , , 1 . A comprehensive parameter optimization was not performed; however, heuristic evaluation for these data (Q6, Q10, and Q14) showed slightly improved convergence for this choice of W compared to uniform weighting. The exact model rank is unknown and for these Kepler test data, a model rank of K= 20 was chosen based on the singular value distribution of the light-curve sample. This singular value distribution was consistent with that shown in Smith et al. (2012) where a Kepler CBV rank of K= 16 was adopted. The initial C0 was obtained from PCA applied to Y. The choice of p in the total-variation constraint determines the degree of smoothness in the spatial correlation and the tolerance of discontinuities (Appendix C). In Appendix D the choice of p is interpreted in a Bayesian framework, where for p= 1 the spatial constraint is equivalent to a Laplacian prior on the difference of coefficients, and for p= 2 the spatial constraint is equivalent to a Gaussian on the difference of coefficients. Referring to Figure 6, PCA fitted coefficients representative of the underlying systematic structure show a mostly spatially uniform structure with discontinuities at module edges. We adopted a value p= 1.1, chosen empirically as the observed coefficients are closer to Laplacian in form. Gradient step size α t at iteration t was set using a backtracing line search (Armijo 1966) starting from α t= 0.1. This was found to allow for the minimization to progress adequately during initial iterations but with fine-tuning in later iterations, thereby saving computation time. Our stopping condition was based on the difference in our cost function in Equation 10 between successive iterates falling to a negligible level ∥w(Ct+1)−w(Ct)∥� 10−5. Figure 6. Leading PCA basis vector v1 for quarters Q6, Q10, and Q14 (top) and an associated color map of the leading PCA coefficient ¬( )ci x y, 1 at each spatial cell position (bottom). The spatial structure of the leading coefficient term is persistent between quarters. 9 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi 2.5.3. Simulated Astrophysical Signals We describe here the functional form of the simulated transient and variable signals used in the injection tests and full simulations described below. We denote an individual simulated astrophysical signal in light curve form as vector a defined on n ä {1,.., N} and relative to normalized light curve y: ∥y∥= 1. The functional form of the simulated signals was chosen to be either a sinusoid (as), a periodic exoplanet transit (at), or a transient flare (af). Figure 7 shows a simulated example signal of each type. These cover a representative range of timescales. The simulated signals were randomly generated relative to light curve y as follows: 1. Sine wave: pb=[ ] ( )a n A nsin 2s . The amplitude A was drawn from the uniform distribution U(0.005∥y∥, 0.015∥y∥). The angular frequency β was drawn from ( )U , N N 4 8 . 2. Simulated exoplanet transit signal: at[n]. A transit signal is simulated as a periodic repeating transit profile b[n], with randomized parameters: orbital period P, epoch n0, duration d, and transit depth δ. d = - -⎧ ⎨⎩ [ ] · [( ) ] ( ) ( )  a n b n n P n n P dmod , if mod 0, otherwise . 13 t 0 0 The transit profile b[n] is an empirical tapered symmetric profile, defined on the first half of the transit Î [ ]n 1, d 2 as: = - -[ ] ( )b n n0.4 exp 0.3 1 and on the second half Î +[ ]n d1,d 2 as = - - -[ ] ( ( ))b n d n0.4 exp 0.3 1. The depth of this profile generally ranges from −0.6 at transit ingress and egress to −1 at the midpoint. Transits were simulated to allow for at least three transits to occur in the light curve. The period P was drawn from the uniform distribution U(96, 360), which corresponds to a range of 4–15 days. The transit epoch n0 is drawn from the uniform distribution U(0, P). The duration d of the transit profile was drawn from the uniform distribution U (8, 24), which corresponds to a range of 4–12 hr. The depth δ is selected between U(0.002∥y∥, 0.016∥y∥). 3. Simulated flare: af[n]. The time sample index np of the flare peak was drawn from the uniform distribution U(1, N). The flare profile is simulated as a rising exponential up to np followed by a slower decaying exponential. The overall signal is multiplied by a random walk stochastic process x[n], generated as + = +[ ] [ ]x n x n1 = =[ ] [ ] ( ) [ ]w n w n x: 0, 1 , 0 0. The simulated flare is defined on the interval n ä [1, np] as = - -( )[ ] [ ] ·( ) a n x n expf n n N 20 p and on the interval n ä [np, N] as = - -( )[ ] [ ] ·( ) a n x n expf n n N 10 p . The simulated signals vary in strength but are generally not weaker than the sample variance of the first-order difference across each light curve (Section 2.5.1), and not greater than the magnitude of systematics. Each injected signal has on average 30% of the energy of the systematics term =    0.3a y s . 2.5.4. Experiment A: Injected Transient and Variable Signals In this experiment (A) we numerically evaluated the performance of the spatial systematics algorithm in the recovery of simulated signals (Section 2.5.3) injected into the gridded Kepler test data (Figure 2), here selected for quarter Q10. We denote the preprocessed gridded Kepler test data here as matrix Y (Section 2.3). We simulated nine astrophysical signals, with random parameters drawn as described in Section 2.5.3, divided equally into three types: sinusoid, flare, and exoplanet transit. Each of the nine individual signals was randomly assigned to a separate randomly selected light curve in the test data in a one-to-one mapping. A priori the Kepler test data contain unknown astrophysical variability. Therefore choosing a small number of injected signals matches this realistic incidence and avoids significantly biasing the sample. We represent the injected signals as matrix A, with the same shape as Y, and where A is null except for nine randomly selected columns containing the simulated injection signals. The post-injection data are denoted: Y→ Y+ A. A least-squares PCA solution was obtained using the data model Y= VPCACPCA+N (Section 2.3), yielding detrended light curves using PCA: ¢ = -Y Y V CPCA PCA PCA. The post- injection Kepler test data Y were then detrended using the spatial systematics algorithm (Algorithm 1) using CPCA to initialize C and yielding detrended light curves ¢Y . For both methods, we adopt a model rank K= 20. This model rank is empirically sufficient to represent the systematics and is smaller Figure 7. Example simulated sine as, transit at, and flare af signals as defined in Section 2.5.3. The inset panel (blue outline) shows the tapered transit profile at higher temporal resolution. 10 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi than both the number of light curves 1500 and the length N∼ 4000 of each light curve. This experiment was performed over an ensemble of ten instances of the random parameters defining the injected signals (Section 2.5.3) and random light-curve assignment from Y. The performance of both the spatial systematics and PCA algorithm was evaluated in terms of the estimated level of residual systematics after detrending and in terms of their recovery of injected signals, as described in further detail below. 2.5.5. Experiment B: Fully Simulated Light Curves In this experiment (B) we fully simulated the light curves using only the coordinate framework of the Kepler test data. No simple generative model exists to fully simulate light curves due to the complex origin of residual systematics across the sensor, as discussed in Section 1. Accordingly, we adopted a simple low-rank linear model for the systematics L=VC with rank K= 4. The basis vectors V were selected at approximately equal spacing as index set (1, 5, 11, 16) from the 16 Kepler cotrending basis vectors (CBV) published for Q10 and channel 30 (Stumpe et al. 2012). The coefficients C were simulated with a random mix of discontinuities and smooth features, broadly comparable to the spatial structure found for the Kepler sensor test data but not constrained to be an exact match, including module gridding (Figure 6). We chose this spatial structure to be representative but distinct from the Kepler test data. The basis vectors V and coefficients C used in the full simulations are shown in Figure 8. Median normalization was applied to each simulated systematic term as = -ˆ ( )l l lmed 1i i i , where med(li) is the median of li for iä I. We note that a PCA decomposition of the simulated systematics will not be identical. The simulated basis vectors and coefficients are rank K= 4; however, the simulated systematics are rank Ks= 5. An extra basis term is introduced in a small number of light curves due to the median normalization of the simulated systematics. The corresponding basis vector is a constant offset. We simulated 1500 light curves in the 30× 50 cell region, chosen to mirror the Kepler sensor region and discretization used in Section 2.5.1. A total of 300 astrophysical signals were simulated and are represented as matrix A as introduced in Section 2.5.4. Equal proportions of Figure 8. Basis vectors V (upper) and coefficients C (lower) used for the full simulation of light curves in experiment B using a rank K = 4 model. The basis vectors were selected at approximately equal spacing as index sets (1, 5, 11, 16) from the 16 Kepler CBV published for Q10 and channel 30 (Stumpe et al. 2012). The coefficient spatial dependence was simulated as a mixture of random smooth features and sharp discontinuities. Both components are arbitrarily chosen and the discontinuities are not constrained to match the Kepler module gridding. 11 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi sinusoid, flare, and exoplanet transit signals were simulated with random signal parameters drawn as described in Section 2.5.3. Each simulated signal was injected into a single randomly selected light curve in a one-to-one mapping. The light curves Y (Section 2.3) were simulated as Y=A+VC+N, where N is iid Gaussian noise  s~[ ] ( )N 0,n i, 2 , where s = å Î ( ) ∣ ∣ lvar I i I i 2 1 4 . The simula- tions were performed over one run as no ensemble was necessary given the sufficient sample size of simulated and injected astrophysical signals. As in Experiment A (Section 2.5.4), both the PCA method and the spatial systematics algorithm were used to detrend the data Y using model rank K= 8, yielding ¢YPCA and ¢Y , respectively; as before the PCA coefficients were used to initialize the spatial systematics algorithm (Algorithm 1). This experiment deliberately simulates a high level of overfitting as the fitted model rank exceeds the rank of the simulated systematics. The performance of the algorithms was evaluated in terms of the accuracy with which the known systematics and astrophysical signals were recovered and the level of residual noise in the detrended light curves, as described in further detail below. 2.5.6. Experiment C: Comparison with Standard Detrending Methods We performed a high-level comparison against several standard detrending methods including the Kepler cotrending basis vector method (CBV; Smith et al. 2012; Stumpe et al. 2012), self-flat-fielding (SFF; Vanderburg & Johnson 2014), and pixel-level-decorrelation (PLD; Deming et al. 2015; Luger et al. 2016, 2018), all as implemented in the Lightkurve package (Lightkurve Collaboration et al. 2018). We caution that this is a high-level comparison only and is not intended to assess the optimality of any method. For this comparison, the gridded Kepler targets (Figure 2) selected for Q10 were used. The spatial systematics method was applied to these data with rank K= 20 after preprocessing (Section 2.5.1; Figure 2) and with no injection of astrophysical signals. However, for the Lightkurve methods, the data for this target list were down- loaded and processed inside that package as is customary. In general, for all Lightkurve methods, we used the recommended parameters in the Lightkurve tutorial.10 Speci- fically for the CBV and SFF detrending methods SAP preprocessing was performed using the remove_nans and remove_outliers functions. The PLD detrending method uses pixel-level data and no additional preprocessing was applied. The Lightkurve module masks bad quality-flagged cadences and removes ∼100 cadences per light curve for these test data. We did not remove these cadences in our processing using the spatial systematics method (Figure 2). However, in our calculation of performance metrics, we use the intersection of nonmasked cadences among Lightkurve light curves. This excludes a further ∼100 cadences per light curve. For the Lightkurve CBV detrending we used all single-scale CBVs available K= 16. For Lightkurve SFF detrending we specified a window of width 401 cadences for preprocessing; this removes long-term variability and flattens each light curve. The Lightkurve SFF detrending was performed using 20 such windows. No other Lightkurve input parameters were specified. The broad relative performance of the spatial systematics algorithm against the comparison of Lightkurve methods was evaluated in terms of the estimated level of residual systematics after detrending (Section 3.1). 3. Results In this section, we describe the results of experiments A, B, and C. 3.1. Performance Metrics The metrics used to numerically evaluate algorithm performance fall into the following two broad categories. Residual Noise in Detrended Light Curves: This section concerns metrics used to assess the level of residual systematic noise in detrended light curves. The combined differential photometric precision (CDPP) over 6 hr is used in the Kepler pipeline and serves as an estimate of the level of white noise remaining after detrending on a transit timescale (Christiansen et al. 2012). Following Aigrain et al. (2016, 2017) we used the procedure described in Gilliland et al. (2011) to compute an approximate measure of the CDPP in our detrended light curves. A 2D quadratic Savitsky–Golay filter was applied to the detrended light curves to remove low-frequency variability and the data were then averaged over 6 hr bins. The approximate CDPP was then computed as the standard deviation of the binned time series and scaled by a value of 1.168. The median CDPP across the set of detrended light curves is reported. The second metric for residual systematic noise was adopted from the goodness metric defined by Stumpe et al. (2012) and is denoted here as ¢Gy . This metric is formed as the average absolute cubed pairwise correlation between detrended light curves as r= å ¢ ¢¢ - ¹ ∣ ( )∣ ( ) y yG ,y N N i j i j 2 1 3. Lower values of ¢Gy suggest a lower residual systematic noise in detrended light curves. The cube downweights correlations, which may be spuriously low if all residual astrophysical variability is removed and higher correlations make a more significant contribution to this metric. Cross-Correlation against Known Signals and Systematics: In this section, we consider metrics used to measure the degree to which known astrophysical signals or systematics are recovered in the injection tests or full simulations. The injection matrix A defined above is a null matrix with a subset of columns containing simulated vector light curves that are either sinusoids (as), flares (af), or exoplanet transit signals (at), as defined in Section 2.5.3. Each simulated light curve is mapped to a random light curve in a one-to-one mapping. We denote the mean correlation of all simulated astrophysical signals of type a[m],mä {s, t, f} against their associated detrended light curves ¢[ ]y m and averaged over all simulation runs as r ¢¯ ( )[ ] [ ]a y,m m . This is a measure of the degree to which the simulated or injected astrophysical signal has been recovered from the detrended light curve. When averaged across all signal types this is denoted r ¢¯ ( )a y, a . We denote the mean absolute correlation of all simulated astrophysical signals of type a[m] against the set of all detrended light curves where no astrophysical signal was injected ¢yac and averaged over all simulation runs, as r ¢ Î∣ ¯ ( )∣ { }[ ]a y m s t f, , , ,am c . The absolute value is used as the simulated signal may be anti- correlated with the estimated systematics. This correlation measures whether a detrending algorithm is corrupting light curves that do not contain an injected astrophysical signal ¢( )yac . 10 https://docs.lightkurve.org/tutorials/index.html#removing-instrumental- noise 12 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi https://docs.lightkurve.org/tutorials/index.html#removing-instrumental-noise https://docs.lightkurve.org/tutorials/index.html#removing-instrumental-noise The unprocessed light curves before detrending ( )yac may, however, be incidentally correlated with an injected astrophysical signal and therefore we normalize by the mean absolute correlation of the nondetrended light curves r r ¢∣ ¯ ( ) ∣ ∣ ¯ ( ) ∣ a y a y , , a a c c . Values below unity indicate that the detrending algorithm is not introducing spurious correlation with the injected signals beyond that present incidentally in the nondetrended light curves. In Experiment B, the constituent systematic light-curve vectors li, (i ä I) comprising matrix L are also known, in addition to the injected astrophysical signals a[m], m ä {s, t, f}. We denote the estimated systematics vector for light curve i as l̂i. The mean correlation of the simulated li and estimated systematics vectors l̂i averaged over all light curves is denoted r̄ ( ˆ)l l, . If this metric is restricted to light curves where an astrophysical signal a[m], m ä {s, f, t} was injected, it is denoted r̄ ( ˆ )l l,a a . 3.2. Experiment A Experiment A, as described in Section 2.5.4, was performed to numerically evaluate the spatial systematics algorithm in terms of its ability to recover transient and variable signals injected into Kepler test data. A standard PCA decomposition method was used as a comparison. The convergence of the spatial systematics algorithm for five representative runs of this experiment is shown in Figure 9. This figure plots the mean correlation r ¢¯ ( )[ ] [ ]a y,m m (over the light curve) of the injected astrophysical signal types mä {s, t, f} as a function of the spatial systematics iteration number. This mean correlation is a proxy for the degree of recovery of injected astrophysical signals, as described in Section 3.1. Figure 9 shows the general convergence of the spatial systematics algorithm within ∼30 iterations. The figure shows signal recovery performance comparable to or exceeding the reference PCA method in the majority of cases but not without exception for individual flare and transit runs. The performance of the spatial systematics algorithm is shown in more detail for a single representative run in Figures 10–12. Figure 10 shows the estimated basis vectors vk obtained with the spatial systematics and reference PCA algorithms. The estimated basis vectors vk are generally similar but differences increase at higher-order k. Figure 11 shows the Figure 9. Scaled mean correlation of injected signals and associated detrended light curves r ¢¯ ( )[ ] [ ]a y,m m for Experiment A for sine, flare, and transit signal types m ä {s, t, f} (across columns) for each of five separate simulation runs (across rows). This mean correlation, defined in the main text, is plotted for the spatial systematics algorithm (blue) and the PCA decomposition method (red). The mean correlation is scaled by the maximum value achieved by either method; this maximum is shown in the legend of each subplot. The PCA algorithm is noniterative and is depicted as a straight line accordingly. The x-axis label is the iteration number of the spatial systematics algorithm. The spatial systematics algorithm is generally convergent and typically has comparable or improved performance relative to PCA but this is not universal for all runs. 13 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi Figure 10. Estimated basis vectors vk, k ä {1, 20}(K = 20) obtained for run #0, a single representative run of Experiment A, shown here for the spatial systematics algorithm (blue) and the reference PCA method (red). The spatial basis vectors and PCA basis vectors closely overlap for the first two terms. Figure 11. Basis vector coefficients Î ={ }( )c k K, 1 ,.., 5 20i k for light curve i → (x, y) ä (X, Y) obtained using the spatial systematics (top row) and PCA algorithms (bottom row) for run #0, a representative run of Experiment A. The coefficient index k is shown in the top left of each subfigure. The x- and y-axes are in units of gridded spatial cells. 14 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi fitted coefficients Î ={ } ( )c k K, 1 ,.., 5 20i k obtained by the spatial systematics and PCA algorithms across the gridded light curves i→ (x, y) ä (X, Y). This figure shows that the total- variation constraint in the spatial systematics algorithm enforces smoothness in the coefficients while preserving discontinuities. In Figure 12, examples of detrended light curves are shown over injected astrophysical signal type. In this example, the spatial systematics algorithm has preserved the injected astrophysical variability to a greater degree than the reference PCA method. However, the spatial systematics algorithm has a higher residual scatter in the detrended light curves than PCA, as is visible near cadence 3000 for the light curve injected with a sine signal. The summary noise metrics (Section 3.1) for Experiment A for the final iteration of the spatial systematics algorithm and averaged over all runs in the ensemble are listed in Table 1. This table shows that for these simulations the spatial systematics algorithm generally outperforms the reference PCA method with caveats discussed in further detail in Section 4. 3.3. Experiment B Experiment B (Section 2.5.5) was designed to numerically evaluate the spatial systematics algorithm using fully simulated data generated using only the Kepler test data coordinates. The algorithm was evaluated in terms of its ability to recover the known astrophysical signals and systematics used in simulating the light curves. The simulated systematics are rank Ks= 5 but were estimated with rank K= 8 as described above. The estimated and simulated basis vectors vk, k ä {1,..,8} for Experiment B are shown in Figure 13. The first four estimated basis vectors are identical between the spatial systematics and PCA algorithms. Higher-order basis vectors deviate between the two methods and are also poorly constrained. However, fitted coefficients for these basis terms (K> 4) are small. For the spatial systematics algorithm, these fitted coefficients contri- bute ∼4% in total magnitude (å å= =   c ck K k k K k5 1 ), while for PCA the contribution is ∼11%. Figure 14 shows the true simulated coefficients and those estimated using the PCA and spatial systematics algorithms for the first 5/8 basis terms. The coefficients estimated using PCA have a more speckled spatial substructure than the smoother distribution obtained by the spatial systematics algorithm. Example individual detrended light curves for sine, transit, and flare simulated signal types are shown in Figure 15. Summary results for Experiment B, in terms of the metrics described in Section 3.1, are provided in Table 2. 3.4. Experiment C Experiment C was designed to provide a high-level comparison of the spatial systematics algorithm against the standard CBV, SFF, and PLD methods (Section 2.5.6). The Figure 12. Three sample light curves across injected signal type a[m] m ä {s, t, f} (sine, transit, and flare) from run #0, a representative run of Experiment A. The figure shows the preprocessed light curve before detrending y (upper row), the detrended light curve using PCA ¢yPCA (middle row), and the detrended light curve using the spatial systematics algorithm ¢y (lower row). Simulated injected signals are overlaid in black. The x-axis is Kepler's long-cadence time sample index. This sample shows an instance where PCA detrending removed the injected astrophysical sine and flare signals, while spatial detrending preserved those signals. 15 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi Table 1 Experiment A Results Sine Recovery Sine Corruption Transit Recovery Transit Corruption Flare Recovery Flare Corruption Signal Recovery Residual Systematics r ¢¯ ( )a y,s s r r¢∣ ¯ ( )∣ ∣ ¯ ( )∣a y a y, ,a as sc c r ¢¯ ( )a y,t t r r¢∣ ¯ ( )∣ ∣ ¯ ( )∣a y a y, ,a at tc c r ¢¯ ( )a y,f f r r¢∣ ¯ ( )∣ ∣ ¯ ( )∣a y a y, ,a af fc c r ¢¯ ( )a y, a ¢Gy CDPP6h (ppm) PCA 0.14 ± 0.04 0.12 ± 0.02 0.45 ± 0.08 0.60 ± 0.11 PCA 0.35 ± 0.13 0.09 ± 0.03 0.31 ± 0.08 (6.61 ± 0.12) ×10−4 38.91 ± 0.03 Spatial 0.57 ± 0.11 0.85 ± 0.09 0.44 ± 0.06 0.80 ± 0.12 Spatial 0.64 ± 0.15 0.54 ± 0.09 0.55 ± 0.11 (1.73 ± 0.04) ×10−2 40.33 ± 0.14 Ratio(PCA/Spatial) 0.25 L 1.01 L Ratio(PCA/Spatial) 0.54 L 0.56 0.04 0.97 Note. Performance metrics for detrended light curves, as defined in Section 3.1. As described there, recovery indicates the degree to which injected signals are preserved correctly in the detrended light curves. Corruption indicates the degree to which injected signals incorrectly appear in light curves that contain no injected signals. The uncertainties listed are the standard deviation of the value over the ensemble of ten runs. For reference, the correlation of nondetrended, preprocessed, and noninjected light curves yac with simulated astrophysical signals a[m] have values r =∣ ¯ ( )∣a y, 0.07as c , r =∣ ¯ ( )∣a y, 0.04at c , r =∣ ¯ ( )∣a y, 0.18af c for sine, transit, and flare signals, respectively. For nondetrended preprocessed light curves, the goodness metric is Gy is 0.26 ± 0.01 and the median CDPP6h is 51.21 ± 0.14 ppm. 16 T h e A stro n o m ica l Jo u rn a l, 167:60 (29pp), 2024 F ebruary T aaki, K em ball, & K am alabadi performance metrics ¢Gy and CDPP6h for the residual noise in the detrended light curves for this experiment are shown in Table 3. A composite plot of CDPP6h versus ¢Gy across all methods in Experiment C is shown in Figure 16. These results show that the CDPP6h metric is broadly comparable across all methods while the ¢Gy metric has a greater dependence on the detrending method. As noted in Section 2.5.6, the data used in this comparison have undergone different preprocessing, and furthermore detrending methods have not been optimized in each specified case. The purpose of this comparison is to provide a high-level reference comparison and not to show the general optimality of any single method. 4. Discussion The numerical evaluation of the spatial systematics algo- rithm presented here, using injected signals in Kepler data (Experiment A) and full simulations (Experiment B), shows that in these cases the algorithm matches or outperforms the reference PCA method for astrophysical signal recovery and achieves comparable performance for systematics removal. This is indicated by higher values of r ¢¯ ( )a y, a and r̄ ( ˆ )l l,a a , respectively, in Tables 1 and 2. The improved recovery of injected astrophysical signals by the spatial systematics algorithm is direct evidence that the algorithm is less prone to erroneously absorbing true astrophysical variability in the systematics (overfitting) than the reference PCA method. Astrophysical variability and systematics are not separable a priori and overfitting is thus a foundational concern in the detrending of exoplanet transit light curves (Stumpe et al. 2012; Smith et al. 2018). Overfitting may remove true astrophysical variability from detrended light curves and also introduce spurious variability into other detrended light curves in the sample as these corrupted systematics are applied. A number of detrending approaches infer a basis representative of systematic effects from a correlated set of light curves Petigura & Marcy (2012), Foreman-Mackey et al. (2015) and Kepler PDC-MAP Smith et al. (2012), Stumpe et al. (2012). The set of light curves from which a set of basis vectors is constructed can be robustly filtered to exclude outliers and those with clear astrophysical variability as a means to mitigate overfitting, as in PDC-MAP (Stumpe et al. 2012). A basis consisting only of systematics (which is difficult to realize) can however still be overfitted to astrophysical features (Smith et al. 2018), as it is unlikely that a basis is completely orthogonal to all possible astrophysical signals. Spatial dependence among light curves has been identified by Petigura & Marcy (2012) and Moreno et al. (2021), whereby the latter work explains spatial correlations as time-delayed systematic effects traversing the sensor. The modified total variation prior applies a spatial correlation constraint across the sensor, which is well suited for inference of local effects of this form. Functionally, the spatial systematics model depends on systematic coefficients alone, this model is differentiable, and derived gradients are computationally tractable. The algorithm iteratively refines both the fitted coefficients and the overall systematics model. These algorithm properties help to separate the systematics and astrophysical variability and thereby mitigate overfitting. A position-based prior is used in PDC-MAP to reduce the Figure 13. Estimated (left) and known simulated (right) basis vectors vk, k ä {1,..,8} for Experiment B. The estimated basis vectors (left) are shown for the spatial systematics algorithm (blue) and the reference PCA method (red). The x-axis is in units of Kepler long-cadence time sample. The estimated spatial and PCA basis vectors (left) for the first four terms overlap one another in these plots. 17 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi overfitting of cotrending basis vectors to a light curve. The spatial systematics algorithm is novel in that it uses a spatial constraint computed directly from systematic estimates, in total-variation form, and further iteratively refines the overall systematics solution for a collection of light curves. We emphasize the comparative performance to standard PCA and while we perform a high-level comparison to specialized detrending methods (Experiment C), we do not claim optimality for this exploratory algorithm development. The reduced overfitting is most evident in Tables 1 and 2 for sine r ¢¯ ( )a y,s s and flare r ¢¯ ( )a y,f f signals and comparable for transit signals r ¢¯ ( )a y,t t . This is also visible in the sample detrended light curves shown in Figures 12 and 15. In summary, the longer-duration astrophysical signals (sine, flare) have reduced overfitting relative to the transit signals, which are of shorter duration. The shorter the duration of an astrophysical signal relative to the timescale of systematics present in Kepler light curves (Figure 10), the less linearly dependent on an inferred systematics basis and therefore the smaller the expected improvement over the reference PCA method. Notably, the spatial systematics algorithm performs well over all astrophysical signal characteristic timescales considered here. The spatial systematics algorithm contains no explicit model for transient or variable astrophysical signals; however, it performs well for the signal types considered here (Section 2.5.3). The consistent relative performance between signal types in Experiments A and B further supports the preceding proposed causal link between the signal and characteristic residual systematic timescales, rather than the exact functional form of the astrophysical signal. Further evidence for the improved mitigation of overfitting by the spatial systematics algorithm over the reference PCA method is provided by Experiment B. As a full simulation, this experiment includes no unknown astrophysical variability beyond the simulated astrophysical signals. In Figure 14 the true simulated coefficients shown in the top row are spatially smooth but with realistic module discontinuities. In contrast, the coefficients in the second row estimated using the PCA method show speckling in smooth regions. These discrepant coefficient values correspond to cells containing light curves with simulated astrophysical signals that are overfitted by the PCA method. The coefficients estimated by the spatial systematics method shown in the bottom row of this figure Figure 14. Coefficients Î ={ }( )c k K, 1 ,.., 5 8i k for light curve i→ (x, y) ä (X, Y) for Experiment B. This figure shows the true simulated coefficients (top row), the coefficients obtained using the PCA method (middle row), and the coefficients obtained using the spatial systematics algorithm (bottom row). The coefficient index k is shown in the top left of each subfigure. The x- and y-axes are in units of gridded spatial cells. Coefficients fitted with PCA show speckling, indicative of the overfitting of injected astrophysical signals. 18 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi generally have reduced speckling and reflect the true coefficients with greater fidelity. Figure 17 shows the cells where an astrophysical signal was injected into a simulated light curve, alongside fitted coefficients from PCA and the spatial systematics method. This figure supports the argument that PCA has overfitted light curves containing astrophysical variability. However, we note that not all individual cells containing injected astrophysical signals are overfitted by PCA. A light curve and the associated detrending results corresp- onding to an overfit cell in Figure 17 are shown in Figure 18. The spatial systematics algorithm varies both the coefficients and basis vectors (Figure 13) in the iterative fit. Empirically, in Experiment B the spatial systematics algorithm estimates systematics l̂ with a correlation r̄ ( ˆ)l l, , against true simulated systematics l, which is comparable to the reference PCA method. The correlation r̄ ( ˆ )l l,a a against systematics for light curves containing a simulated astrophysical signal, and there- fore a higher risk of overfitting, is slightly higher for the spatial systematics algorithm relative to the reference PCA method. This implies the spatial systematics algorithm shows reduced overfitting relative to the reference PCA method. We believe this is due to the spatial total-variation constraint appropriately estimating spatially correlated systematics. As a definitional consequence of the reduction in overfitting, the detrended light curves obtained using the spatial systema- tics algorithm will retain a higher degree of true astrophysical variability relative to those obtained using the reference PCA method in addition to a component for uncorrected systematics, in keeping with all detrending algorithms. We note that PCA is guaranteed to maximally flatten a collection of light curves as the PCA solution is equivalent to a minimal least-squares residual between fitted systematics and light curves (Equation (6)). PCA implicitly assumes a white noise model and is susceptible to overfitting of astrophysical signals (Candès et al. 2011) because astrophysical signals are not appropriately modeled as white noise Pont et al. (2006). In both experiments, the two algorithms obtained low goodness metric values <¢G 2%y indicating largely successful systematics removal in the detrended light curves. The residual systematics metrics ¢Gy and CDPP6h are plotted for run #0 of Experiment A per light curve and over-detrending algorithm in Figures 19 and 20. Figure 19 shows that the distribution of the goodness metric ¢Gy is broadly larger than that for PCA. This metric (Section 3.1) is nonlinear and pairwise, and therefore sensitive to those light curves with poorly modeled systematics (Figure 20). For example, individual light curves may be poorly modeled by the spatial systematics algorithm if their systematics exhibit strong temporal variation insufficiently captured by the fixed spatial model. It is also possible that the spatial systematics algorithm retains a greater degree of astrophysical variability that biases ¢Gy through incidental correlation. We lack sufficient evidence to state either claim strongly here. Future generalizations of the constraint terms are described below as potential future work. Figure 15. Three sample light curves across injected signal type a[m] m ä {s, t, f} (sine, transit, and flare) for Experiment B. The figure shows the simulated light curve before detrending y (upper row), the detrended light curve using PCA ¢yPCA (middle row), and the detrended light curve using the spatial systematics algorithm ¢y (lower row). Simulated injected signals are overlaid in black. The x-axis is Kepler's long-cadence time sample index. This sample shows an instance where PCA detrending moderately removed the injected astrophysical sine and flare signals, while spatial detrending more clearly preserved those signals. 19 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi In Experiment B, the goodness metric value ¢Gy was lower for the spatial method than the reference PCA method. For both experiments (Tables 1 and 2) the two algorithms obtained very similar CDPP6h metrics and a reduction from the nondetrended preprocessed light curves, indicating that both methods were equivalently successful at suppressing transit timescale systematics. In Experiment A the measures of residual systematics ¢Gy and CDPP6h are higher for the spatial systematics method compared to the PCA method. As noted above, these values for the detrended light curves may be elevated primarily due to improved preservation of true astrophysical variability or alternatively due to unmodeled residual systematics. In Figure 12 the example detrended light curve from Experiment A obtained by the spatial systematics algorithm in the first column shows an example of a large residual systematic scatter not present in the PCA detrended light curve; however, the overall variability is dominated in this case by the injected astrophysical signal. We argue, that although the spatial method may retain a greater level of residual systematics, comparably the major component of retained light-curve variability is astrophysical. As discussed above, Experiment B demonstrated a more accurate estimation both of systematics l̂ and of astrophysical signals a in light curves containing astrophysical signals. Experiment B, although simulated with a simplified data model, therefore aids the interpretation of metric performance. We note also that the improved preservation of true astrophysical variability in the light curves detrended by the spatial systematics algorithm may reduce the correlations r ¢ ί ( ) { }[ ] [ ]a y m s f t, , , ,m m in Tables 1 and 2 discussed above and that these may be underestimated as a result. In Experiment A, both the spatial systematics algorithm and the reference PCA method produce values r r¢ <∣ ¯ ( )∣ ∣ ¯ ( )∣[ ] [ ]a y a y, , 1a am mc c , for mä {s, t, f}. Therefore, as discussed in Section 3.1, this suggests neither algorithm has increased the average incidental correlation between the simulated astrophysical signals and the light curves in which no astrophysical signals were injected (the noninjected subsample). Equivalently, neither method has introduced spurious astrophysical content into this noninjected subsample of light curves. However, the ratio r r¢∣ ¯ ( )∣ ∣ ¯ ( )∣[ ] [ ]a y a y, ,a am mc c is higher for the spatial systematics algorithm compared to the PCA method (Tables 1 and 2), particularly for sine and flare mä {s, f} astrophysical signals, which have a longer duration. In Figure 21 we plot the light curves in the noninjected subsample for which the spatial systematics algorithm produces Table 2 Experiment B Results Systematics Recovery Signal Recovery Residual Systematics r̄ ( ˆ )l l,a a r̄ ( ˆ)l l, r ¢¯ ( )a y,s s r ¢¯ ( )a y,t t r ¢¯ ( )a y,f f r ¢¯ ( )a y, a ¢Gy CDPP6h PCA 0.96 ± 0.08 0.99 ± 0.03 0.52 ± 0.30 0.28 ± 0.15 0.39 ± 0.17 0.39 ± 0.24 1.71 × 10−4 0.39 ± 0.01 Spatial 0.98 ± 0.08 0.99 ± 0.03 0.63 ± 0.22 0.28 ± 0.15 0.48 ± 0.19 0.47 ± 0.24 1.55 × 10−4 0.40 ± 0.01 Ratio (PCA/Spatial) 0.98 0.99 0.82 0.98 0.79 0.82 1.10 0.98 Note. Performance metrics for detrended light curves, as defined in Section 3.1 and used in Table 1. As described there, recovery indicates the degree to which injected signals are preserved correctly in the detrended light curves. The uncertainties listed are the standard deviation of the value over the ensemble of simulated light curves. No uncertainty is listed for ¢Gy as it is computed over the ensemble of light curves. The goodness metric value computed over the nondetrended simulated light curves Gy is 0.23. The goodness metric value computed between all simulated astrophysical signals Ga is 2.2 × 10−2. Table 3 Experiment C Results Residual Systematics ¢Gy CDPP6h μ σ Spatial 1.13 × 10−2 38.37 60.17 CBV 8.95 × 10−4 33.16 97.60 SFF 1.45 × 10−3 30.78 24.16 PLD 3.92 × 10−1 36.49 41.44 Note. Residual systematics metrics for detrended light curves produced in Experiment C. The CDPP6h mean over the ensemble of simulated light curves is denoted μ with sample standard deviation σ. No uncertainty is listed for ¢Gy as it is computed over the ensemble of light curves. Figure 16. CDPP6h (ppm) vs. goodness metric ¢Gy per target light curve shown for the spatial systematics method and comparison standard detrending methods CBV, SFF, and PLD. Empirical distribution functions for each light-curve type are plotted over ¢Gy (top) and CDPP6h (right). The CDPP6h metric is broadly comparable across methods while the ¢Gy metric is more strongly dependent on the detrending method. 20 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi the highest incidental correlation r∣ ( )∣[ ]a ymax ,y amac c , mä {s, t, f} between the noninjected spatial detrended light curve and any injected astrophysical signal used in the simulation run. This figure shows the worst cases of incidental correlation; in general, the incidental correlation average is low. As described above, the spatial systematics algorithm is more successful at preserving astrophysical variability in detrended light curves and will therefore retain rare incidentally correlated signals, particularly if their timescale matches the simulated astro- physical signal. We believe this explains the higher values of this ratio in the case of sine and flare signal types. Section 3.4 describes the results of Experiment C, a high- level comparison between the spatial systematics algorithm and the standard detrending methods CBV, SFF, and PLD (Light- kurve Collaboration et al. 2018). As noted above, the default recommended parameters were used for each method as opposed to a customized optimization. In addition, individual methods apply different default preprocessing steps. Accord- ingly, this high-level comparison cannot address the optimality of different detrending methods. In this experiment, the CDPP6h residual systematics metric, which is computed per light curve, shows broadly comparable (�8 ppm between mean values) detrending performance between the spatial systematics algorithm and the reference detrending methods CBV, SFF, and PLD (Table 3; Figure 16). The spatial systematics method has a slightly increased mean Figure 17. Gridded light-curve cells i→ (x, y) ä (X, Y) where a simulated astrophysical signal a[m] was injected for Experiment B. The x- and y-axes are in units of gridded spatial cells. The column at left depicts the magnitude of the injected astrophysical signal ∥a[m]∥ as the color intensity of each cell on a log scale (blue color wedge at left). The middle and right columns depict coefficients ci 5 for light curve i→ (x, y) ä (X, Y) at each spatial cell (from Figure 14) for the reference PCA method and the spatial systematics method, respectively. For these columns, the magnitude of the coefficient is shown as color intensity (color wedge at right). A black outline around a cell indicates that the corresponding light curve contained an injected astrophysical signal. The rank 5 coefficient was chosen for clarity of interpretation as this coefficient should be closer to zero. Coefficients fitted with PCA exhibit some overfitting of injected astrophysical signals. Figure 18. Simulated light curve y before detrending (left), the detrended light curve using PCA ¢yPCA (middle), and the detrended light curve using the spatial systematics algorithm ¢y (right), all at the position of the maximum value of [ ]CPCA 5 in Figure 17. Simulated injected signals are overlaid in black. The x-axis is Kepler's long-cadence time sample index. This shows an instance where PCA detrending overfits the injected astrophysical sine signal relative to the spatial systematics algorithm. 21 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi CDPP6h compared to other methods possibly due to instances of residual scatter where the spatial model imperfectly corrects strong temporal systematic variation. The ¢Gy residual systema- tics metric differs in distribution across the detrending method (Figure 16). By definition (Section 3.1) and as discussed earlier, this metric is nonlinear and pairwise; accordingly, light curves with unmodeled systematics are more heavily weighted. The PLD method has the highest mean value of ¢Gy (Figure 16) but this method retains long-term variability, which is fit by a spline term in the data model, and therefore is expected to have an inherently higher ¢Gy distribution. In contrast, in preproces- sing, SFF light curves are flattened using a Savitsky–Golay filter (Lightkurve Collaboration et al. 2018), and SFF is therefore broadly comparable to the CBV method in this metric. The spatial systematics algorithm has a higher mean value of ¢Gy and broader distribution than the CBV and SFF methods (Figure 16). As discussed earlier for Experiments A and B, this may be due to the spatial systematics method retaining a higher level of astrophysical variability or due to the nonlinear contribution to ¢Gy of light curves that do not fit the spatial systematics data model assumptions. Generalizations to make the constraint terms more complete are discussed below as potential future work. Our current work has several limitations that would benefit from exploration in future work. As described in Section 2.5.2, we chose tunable algorithm parameters, including model rank, convergence criteria, optimization step size, weighting matrix, and initial coefficient values, among other parameters, based on feasibility and empirical evaluation alone. It would be beneficial to sample this algorithm parameter space over a broader range, for completeness. In conjunction, our numerical evaluations relied on an associated set of injected astrophysical signal types (Section 2.5.3), which could be expanded. Generalizing the sensor spatial discretization to nonuniform layouts is left for future work. It has been noted that light curves at small spatial separations may suffer from the blending of astrophysical signals (Kovács et al. 2005; Hattori et al. 2022) requiring future consideration when using the spatial systema- tics algorithm on densely populated fields. Although some light curves in our data sample have small spatial separation, this is not a major concern in the current work, as the cell separation on the spatial grid used significantly exceeds the PSF width. The sparse population of the field is likely due to the narrow (and bright) magnitude range selected. There is scope to generalize the constraint term to incorporate nonneighboring pixel spatial correlations or weightings of the spatial prior to excluding targets within a pixel separation where astrophysical blending could occur. In addition, the constraint terms may be generalized to include the known magnitude correlation (Smith et al. 2012; Stumpe et al. 2012) and other parametric dependencies (Moreno et al. 2021). As noted earlier, the Kepler test data were selected in a narrow magnitude range to isolate the effect of magnitude dependence in the current work. 5. Conclusions In this paper, we present an exploratory algorithm for detrending light curves obtained in wide-field exoplanet transit surveys. This spatial systematics algorithm fits a low-rank linear systematics model to a collection of light curves while also including a total-variation spatial constraint across the sensor at a foundational level. The resulting objective function is reduced using variable elimination (Golub & Pereyra 1973), which also stabilizes the solution (Golub & Pereyra 2003; Shearer & Gilbert 2013). An approximate closed-form gradient was developed for minimization by gradient descent; this particular formulation is a modification of total-variation optimization approaches (Vogel 2002). The spatial systematics algorithm was numerically evaluated relative to a reference PCA method using both injection tests with Kepler data (Experiment A) and full simulations including simulated systematics, astrophysical signals, and statistical noise within the same Kepler coordinate framework (Experiment B). The principal conclusions of the paper are: 1. The spatial systematics algorithm showed reduced over- fitting of intrinsic astrophysical variability relative to the PCA method in the two numerical evaluation studies. The reduced overfitting was demonstrated by a comparable or significantly higher correlation between the known astrophysical signals and the detrended light curves in both experiments. In Experiment B the known simulated systematics were comparably or more accurately esti- mated relative to the PCA method, particularly for light curves containing injected signals. We argue that the reduced overfitting likely arises from the physically realistic total-variation constraint. In addition, the algo- rithm simultaneously varies both the basis vectors and their coefficient weights. Both factors likely contribute to the improved separation of systematics and astrophysical signals. 2. A marked reduction in overfitting relative to PCA was found for slowly varying astrophysical signals while comparable performance was achieved for shorter exoplanet transit signals. We argue that the longer- Figure 19. CDPP6h vs. goodness metric ¢Gy per target light curve for a single run #0 of Experiment A. Three light-curve types are shown, including nondetrended preprocessed light curves (black), light curves detrended using PCA (red), and light curves detrended using the spatial systematics algorithm (blue). Empirical distribution functions for each light-curve type are plotted over ¢Gy (top) and CDPP6h (right). 22 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi Figure 20. Metrics ¢Gy (left) and CDPP6h (right) for the reference PCA method vs. the spatial systematics method. Both plots are per target light curve for a single run #0 of Experiment A. Empirical distribution functions for each light-curve metric are plotted over the spatial systematics method (horizontal axis) and the PCA method (vertical axis). Although correlated, there is significant scatter in the left plot due to the intrinsic nonlinear and pairwise nature of the metric ¢Gy . It is accordingly more sensitive to light curves with poorly modeled systematics or greater retained astrophysical variability. Figure 21. Noninjected light curves yac for which the spatial systematics algorithm produces the highest incidental correlation r ( )[ ]a ymax ,y amac c , m ä {s, t, f} between the nondetrended light curve yac and any injected astrophysical signal a[m] used in the simulation run. The top row shows the nondetrended light curve and the astrophysical signal (injected into another light curve) with which the incidental correlation is highest. The detrended light curves ¢yac are shown as obtained by the reference PCA method (middle row) and the spatial systematics algorithm (bottom row). Simulated injected signals are overlaid in black. In these worst cases, it can be seen that the high level of correlation of the spatial detrended light curve ¢y with an astrophysical signal a is incidental, as the original light curve y itself is incidentally correlated with a. This suggests that high values for the metric of corruption r r¢∣ ¯ ( )∣ ∣ ¯ ( )∣a y a y, ,a ac c may result incidentally for the spatial systematics algorithm because it retains a higher level of astrophysical variability. 23 The Astronomical Journal, 167:60 (29pp), 2024 February Taaki, Kemball, & Kamalabadi duration