# Properties of a Previously Unidentified Instrumental Signature in Kepler/K2
that was Confused for AGN Variability

Jackeline Moreno Department of Physics, Drexel University, 32 S. 32nd Street,
Philadelphia, PA 19104, USA  Rachel Buttry  Department of Physics, Drexel
University, 32 S. 32nd Street, Philadelphia, PA 19104, USA  Department of
Physics, Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, PA 15213,
USA  John O’Brien  Department of Physics, Drexel University, 32 S. 32nd
Street, Philadelphia, PA 19104, USA  Department of Physics and Astronomy,
Michigan State University, East Lansing, MI 48824, USA  Michael S. Vogeley
Department of Physics, Drexel University, 32 S. 32nd Street, Philadelphia, PA
19104, USA  Gordon T. Richards Department of Physics, Drexel University, 32 S.
32nd Street, Philadelphia, PA 19104, USA  Krista Lynne Smith Department of
Physics, Southern Methodist University, Dallas, TX 75205, USA

###### Abstract

The Kepler satellite potentially provides the highest precision photometry of
active galactic nuclei (AGN) available to investigate short-timescale optical
variability. We targeted quasars from the Sloan Digital Sky Survey that lie in
the fields of view of the Kepler/K2 campaigns. Based on those observations, we
report the discovery and properties of a previously unidentified instrumental
signature in K2. Systematic errors in K2, beyond those due to the motion of
the detector, plague our AGN and other faint-target, guest-observer science
proposals. Weakly illuminated pixels are dominated by low frequency trends
that are both non-astrophysical and correlated from object to object. A
critical clue to understanding this instrumental noise is that different
targets observed in the same channels of Campaign 8 (rear facing) and Campaign
16 (forward facing) had nearly identical light curves after time reversal of
one of the campaigns. This observation strongly suggests that the underlying
problem relates to the relative Sun-spacecraft-field orientation, which was
approximately the same on day 1 of Campaign 8 as the last day of Campaign 16.
Furthermore, we measure that the instrumental signature lags in time as a
function of radius from the center of the detector, crossing channel
boundaries. Systematics documented in this investigation are unlikely to be
due to Moiré noise, rolling band, or pointing jitter. Instead this work
strongly suggests temperature-dependent focus changes that are further subject
to channel variations. Further characterization of this signature is crucial
for rehabilitating K2 data for use in investigations of AGN light curves.

††journal:

##  1 Introduction

The Kepler satellite was originally conceived as an exoplanet finding mission.
Its rich archival data (observing 2009-2018) has revolutionized planet
discovery and characterization. Kepler’s second mission, K2, has also impacted
many other research fields including astrosiemsology, stellar activity,
microlensing, white dwarfs, supernovae, active galactic nuclei (AGN), and
more, together comprising $\sim
60\%$111https://keplerscience.arc.nasa.gov/publications.html#breakdown-by-
subject of all publications using K2 data thus far. However, poorly understood
instrumental systematics posed a significant impediment to an efficient
exploitation of K2’s thousands of light curves. Part of the problem is that
noise removal is optimized differently depending on the science objective.

Noise removal takes on one definition for planet transit detection and another
for characterization of a faint “unknown” signal, as is the case for AGN in
the short timescale regime ($<100$ days; optical wavelengths). A planet
transit manifests as a strictly periodic signal, characterized by a peak at
the orbital frequency in the power spectral density (PSD) of a light curve. In
contrast, AGN exhibit intrinsically noisy fluctuations that are best
characterized by the entire shape of the empirical power spectrum. The problem
of cleaning data of instrumental noise at nearly all observed frequencies is
critical to enabling science that relies on the entire PSD, which is the case
for aperiodic or stochastic variability.

K2 data provides the most evenly sampled high-cadence (30 minutes)
observations for a large sample of AGN in the optical (over 4000 targets). If
these data are properly mitigated for systematics, the observed sample may be
studied jointly with data from the Zwicky Transit Factory (Bellm et al., 2019)
and the Deep Drilling Fields of the Vera Rubin Observatory (Ivezić et al.,
2008) to probe supermassive black hole (SMBH) physics on timescales of days to
months. The physics in this regime is extremely valuable to constraining the
corona-accretion disk relationship (e.g., Giustini & Proga, 2019) across the
broad spectral taxonomy of AGN (blazars, BL LAC, Seyfert I, Seyfert II, radio
quiet quasars, radio-loud quasars, weak-lined quasars, broad-absorption line
quasars, and more).

Characterizing the short-timescale variability properties of AGN offers the
promise of yielding deeper insights into the accretion physics of SMBHs and
can potentially constrain physical scales. Most importantly, Kepler/K2’s
extremely regular 30-minute cadence and sensitivity (to quasars) will remain
unmatched for the foreseeable decade. Between Mushotzky et al. 2011,Wehrle et
al. 2013, Revalski et al. 2014, Edelson et al. 2014, Shaya et al. 2015, Chen &
Wang 2015, Kasliwal et al. 2015, Dobrotka et al. 2017, and Smith et al. 2018a,
twenty-one (21) Kepler AGN have already revolutionized our understanding of
AGN variability, while $\sim$ 4000 K2 AGN remain un- or under-studied.

Because of the failure of the second of four reaction wheels in 2013, causing
Kepler to be re-purposed as the “K2” mission, the telescope was constrained to
observe along the orbital plane while steadying pointings with solar radiation
pressure and periodic thruster firings. Several re-processing pipelines were
developed to correct pointing errors known as arc-drift (e.g., Vanderburg &
Johnson 2014; Armstrong et al. 2015; Crossfield et al. 2015; Foreman-Mackey et
al. 2015; Huang et al. 2015; Lund et al. 2015; Aigrain et al. 2015, Libralato
et al. 2015). Custom reprocessing software and recommendations were also
published in the literature (Kinemuchi et al., 2012). Although many of these
pipelines were optimized specifically for exoplanet transit detection, the
majority were highly successful at producing arc-drift-free light curves and
enabling science for many types of highly variable astrophysical signals. Many
of these re-extracted light curves are publicly available on
K2-MAST222https://archive.stsci.edu/missions-and-
data/k2333https://keplerscience.arc.nasa.gov/new-keplerk2-high-level-science-
products-available.html.

Large-sample science that may be possible with data in multiple K2 campaigns
requires a streamlined approach that manual custom re-processing cannot scale
to. However, no single prescription to date rehabilitates K2 data for the
study of all types of astrophysical variability. Correlated Kepler/K2 noise
systematics may arise from a combination of error sources such as rolling
band, Moiré noise, faint blended background sources (Van Cleve & Caldwell,
2016), and point spread function (PSF) changes (Libralato et al., 2015). To
characterize systematics, we primarily take advantage of correlations across
multiple sources on the same channel. Using all pixels downlinked by K2 (per
campaign), we detect features that exhibit dependencies on campaign, channel,
magnitude, and distance from the boresight that cause hundreds of targets to
be correlated with each other. Low-frequency and red-noise systematics, which
are not mitigated by arc-drift corrections (Vanderburg & Johnson 2014, Luger
et al. 2018), especially plague AGN and other faint targets (O’Brien et al.,
2018).

Unfortunately, there is no pipeline that sufficiently cleans K2 data for
large-sample AGN variability studies. In this work, we discuss possible ways
to identify systematics in K2 light curves that have been confused for
astrophysical variability in published literature (e.g., Aranzana et al. 2018
for K2 AGN and Dobrotka et al. 2017 for Kepler AGN followed by a retraction in
Dobrotka et al. 2019). On particular channels, systematics may dominate even
at the optimal Kepler magnitude (12th magnitude stars). Finally, we review
available light-curve-processing software that may provide suitable starting
points for rehabilitating K2 data and we provide recommendations for continued
development of reprocessing software that may potentially produce a fully
usable light curve PSDs for scientific investigation.

This paper is organized as follows. Section 2 describes K2 AGN light curves
that are available in Campaigns 0-19 with an emphasis on targets that overlap
with Sloan Digital Sky Survey (SDSS; York et al. 2000) fields. Section 3
discusses known systematics in the first phase of Kepler that were corrected
with cotrending basis vectors (CBVs); Section 4 presents our investigation of
remaining systematics in K2 and discusses possible approaches for mitigating
errors more effectively than CBVs. Section 5 demonstrates the performance of
popular reprocessing methods such as the self-flat-fielding prescription of
Vanderburg & Johnson (2014) and the Everest 2.0 pipeline of Luger et al.
(2018). Finally Sections 6 and 7 summarize our AGN science motivations to
rehabilitate K2 data and the most important findings in this work.

##  2 Data

The Kepler telescope focused on a fixed 116 square degree field of view (FOV)
between Cygnus and Lyra for 4 years prior to K2. The focal plane array
consisted of forty-two $2200\times 1024$ pixel CCDs444For focal plane
labelling and signal map per CCD, see Figures 24 and 25 in
https://archive.stsci.edu/kepler/manuals/KSCI-19033-001.pdf. Each CCD had a
half-maximum bandpass of 435 nm to 845 nm and was read out through two
channels (84 CCD channels also referred to as “module outputs”). The long-
cadence data consisted of 270 frame (6.02 second exposures) integrations for
an epoch every $\sim$ 30 minutes (Van Cleve & Caldwell, 2016). The first
Kepler mission observed continuously for a period of $\sim$90 days (a time
frame referred to as a quarter), then downlinked to Earth, followed by a 90
degree spacecraft rotation before returning to observing mode. A typical
quarter contained $\approx 4000$ epochs known as cadences. However, not all 96
million pixels per cadence were downlinked due to limited bandwidth. A
specific image cutout, referred to as a postage stamp or mask, of fixed size
for the duration of the quarter or campaign was saved for storage. The longest
AGN light curve available, Zwicky 229-015, contains 14 quarters of data (3.9
yrs; e.g., Kasliwal et al. 2015).

A systematic effort was undertaken by several groups to identify bright AGN in
the Kepler Field of View (FOV). By cross-matching observations from ROSAT,
2MASS, and WISE, Edelson & Malkan (2012) identified AGN in the Kepler FOV.
Radio-loud AGN were identified in the Kepler FOV using the VLBA (Wehrle et
al., 2013). Further efforts to identify AGN candidates brought the total
number up to 87 (Kepler GO40041 Edelson). Twenty of these AGN are
spectroscopically confirmed (7 cataloged in Véron-Cetty & Véron (2010) and 13
by Edelson & Malkan 2012). The KSwAGS X-ray survey of the Kepler field
conducted with SWIFT discovered several candidate AGN, 13 of which have been
spectroscopically verified, and many X-ray bright stellar sources (Smith et
al., 2015). Data for all observations acquired during quarters 0 through 17
are publicly available on the Kepler MAST archive (Data Release 25).

Figure 1: Left: Apparent magnitude distribution of AGNs from our K2 proposal that were observed. Right: Redshift distribution of AGNs from our K2 proposal that were observed. Quasar candidates with no available spectroscopic redshifts have a value of -1. The confusion limit of Kepler is reported at about 20th magnitude. Table 1: AGN Targets by Campaign Campaign | Obs. Targs | zmax | imagmax
---|---|---|---
8 | 444 |  $\sim$4 |  $\sim$19.5
10 | 426 |  $\sim$4 |  $\sim$19.5
12 | 337 |  $\sim$4 |  $\sim$20.5
14 | 622 |  $\sim$4 |  $\sim$20.0
16 | 1047 |  $\sim$4 |  $\sim$20.0
18 | 123 |  $\sim$3 |  $\sim$19.0
19 | 1169 |  $\sim$6 |  $\sim$22.0

Pointing jitter was a minor source of error during the first Kepler mission,
causing a maximal 0$\farcs$55 drift of the target centroids across pixels
(3$\farcs$92 in length). After the failure of the second of four reaction
wheels in May 2013, engineers (Ball Aerospace; McCalmont et al. 2015) devised
a way to stabilize pointings. Spacecraft roll error resulted from residual
radiation pressure torques, known in the literature as arc-drift (Cleve et
al., 2016). The Kepler pipeline detrending module (PDC; Smith et al. 2012) was
not updated to correct for arc-drift and the first K2-handbook detailing data
quality was published in 2017 (and updated in 2018; Van Cleve & Bryson 2017).

Kepler’s second mission, K2, deployed an observing strategy consisting of 19
subsequent observing fields along the ecliptic with an $\approx 80$ day
duration. Each of these observation periods is referred to as a campaign.
Campaign 0 provided a commissioning data set used to evaluate the pointing
jitter and necessary adjustments. It pointed at the Galactic ‘Anti-Center’ and
observed 113 AGN targets of various types ranging from relatively ordinary
Seyferts 1 to highly variable blazars. The Campaign 1 FOV pointed at a region
of the North Galactic Cap. The total number of extra-Galactic targets observed
was $\sim$ 550 with $\sim$ 300 thought to harbor AGN. Campaigns 2 $\&$ 3
observed another 266 AGN of various sub-types including bright quasars and
more blazars.

Our group proposed additional observations of known quasars beginning with
Campaign 8. These quasar targets were selected from the SDSS-I/II (Schneider
et al., 2010) and SDSS-III quasar catalogs (Pâris et al., 2017). Robust
photometric quasar candidates (with positions accurate to $\sim$50 mas) are
also drawn from Richards et al. (2004), Bovy et al. (2011), Richards et al.
(2015), and Peters et al. (2015), where the Peters et al. catalog specifically
targeted quasar candidates due to their variability. Thus, our targets are a
heterogeneous K2 sample of known quasars and high-confidence photometric
quasar candidates (most with robust photometric redshifts). In short, we
proposed to target all SDSS quasars (and many quasar candidates) that
overlapped with the K2 field of view in each campaign starting with Campaign
8. The focus was on objects in SDSS Stripe 82 where there already exist $\sim$
100 epochs of photometry (Annis et al., 2014) that can be coupled with K2
photometry to increase the time baseline for investigations of quasar
variability. Figure 1 illustrates the magnitude and redshift distribution of
SDSS quasars targeted by our group in Campaigns 8, 10, 12, 14, 16, 18, and 19.
These targets include objects that are fainter than the typical range of
magnitudes included in the standard pipeline processing (e.g., Luger et al.,
2018). Table 1 gives the number of observed targets and approximate maximum
redshift and faint limit ($i$ band) for each campaign. Note that many more
targets were proposed than were accepted and the process for sub-selecting the
observed targets was not uniform (sometimes favoring bright targets and known
quasars, but other times to reduce high spatial density such as in the Stripe
82 area).

##  3 Kepler/K2 Cotrending Basis Vectors

Low frequency instrumental systematics in the first Kepler mission could be
removed (to an extent) with cotrending basis vectors (CBVs) generated by the
data calibration module, PDC (Stumpe et al. 2012; Smith et al. 2012). CBVs
were calculated for each observing quarter and per channel with a principal
component analysis (PCA) of an ensemble of “quiet” stars. PDC-MAP (maximum
a-posterior) supplied CBV-corrected light curves with weights derived from the
most correlated neighboring sources to individual targets for a predetermined
aperture (using a minimum of 5 vectors).

CBV subtraction was, however, problematic for preserving intrinsic
astrophysical variability. Higher-order vectors (3rd or higher) added variance
and high frequency features (Smith et al., 2012). Any coincidental
correlations between the intrinsic signal and CBVs was projected onto the
model used for correction. CBVs do not mitigate spatially dependent
systematics such as (boresight-dependent) focus changes, rolling band, or
Moiré noise.

The challenge of error mitigation in the first phase of Kepler was typically
approached with custom manual re-processing or calibration with ground-based
data when available (Smith et al. 2018a and Kasliwal et al. 2015,
respectively). For AGN, custom-aperture or PSF photometry was particularly
important in order to minimize noise from background sources in crowded
regions and from extended host galaxies. Custom-aperture optimization and CBV
subtraction was possible with PyKE (Still & Barclay, 2012) or similar software
but without the advantage of simultaneously fitting neighboring sources as in
PDC-MAP that could mitigate issues of over-fitting CBV components.

The compromise of minimizing known and new instrumental noise while preserving
the astrophysical signal of interest carried over into Kepler’s second
mission, K2. While correcting complex arc-drift systematics became the
calibration priority, few publications focused on methodologies for removing
other types of instrumental noise or documenting evidence of Kepler’s natural
degradation.

CBVs are available for K2 in the EVEREST 2.0 package (Luger et al., 2018)
calculated with the SysRem algorithm (Tamuz et al., 2005) from an arc-drift-
corrected light curve ensemble. While the first- and second-order CBVs do an
excellent job at removing the lowest frequency trends, they cannot account for
some of the more troubling and recognizable instrumental signatures that we
identify in our AGN and other variable star samples. For both Kepler and K2
data, CBVs removed the lowest frequency systematics very well, but residual
instrumental trends are at best mixed with variance by the PCA-detrending
approach. Theoretically, PCA cannot account for any local spatial dependence
within a CCD (channel) or for time-dependent systematics across the focal
plane such as temperature-dependent effects. There is a critical need for
other calibration algorithms in order to fully exploit the rich K2 archive.

##  4 Correlated K2-SSF AGN

Figure 2: Expected AGN Variability (simulated). Columns (left to right) show
example light curves with increasing characteristic intrinsic timescales
$\tau_{decay}=10,25,{\rm and}50$ days. These simulated lightcurves include a
constant photometric noise of $5\%$ and constant variability amplitude. The
black curve outlines a smoothed trend through each lightcurve to show that
relatively low freqency trends do not correlate even within columns for
“objects” having the same characteristic timescales. We do not expect AGN to
vary synchronously, except if that trend is due to correlated noise in the
CCD.  Figure 3: Clustering of full postage-stamp AGN light curves from K2
Campaign 1, referenced in Aranzana et al. (2018). The Vanderburg and Johnson
light curves retain similar patterns that can be detected by the same
clustering algorithm. The three different colors highlight just a few of many
“types” of systematic trends that cause independent AGN to be correlated with
each other.

The first evidence that we encountered for the presence of residual
systematics in K2 data after arc-drift removal came from visually identifying
strong correlations between independent target light curves. Inspection of AGN
light curves from K2 indicated that certain patterns were commonly seen,
despite the fact that all AGN light curves are expected to be unique. Further
investigation shows that those common patterns reveal strong systematic
effects in the data.

To illustrate the problem, here we compare some simulated light curves that,
by design, lack features correlated in time, with some K2 light curves that
clearly exhibit correlated features. In Figure 2, we illustrate the AGN
variability that we expected to observe with K2. Following Moreno et al.
(2019), each column is generated with a simple stochastic model (damped random
walk, hereafter DRW) corresponding to characteristic timescales $\tau_{decay}$
$=10,25,$ and $50$ days, respectively. Although the light curves in each
column have the same underlying parameters, each is unique and there are not
obvious correlations between the different light curves in terms of the
locations of the peaks and valleys.

In contrast, in Figure 3 we plot light curves identified using an
agglomerative clustering analysis that detects high correlation in the K2 AGN
sample used for scientific analysis by Aranzana et al. (2018). A visualization
of the raw and smooth AGN light curves reveals trends that exhibit multiple
peaks and toughs in phase across several sources. Such patterns are not
expected between unique AGNs and indicate that systematics are completely
dominating these light curves. Simply put, any matching variability patterns
have nothing to do with the astrophysical targets.

The groups represented in Figure 3 show light curves from clusters with the
strongest covariance. In the right-most panel, we plot cluster members that
exhibit a mixture of astrophysical signal and instrumental trends. Although
this clustering analysis was executed on full postage-stamp raw light curves,
rerunning the analysis with K2-SFF (Vanderburg & Johnson 2014) from MAST also
reveals the same correlated behavior that we attribute to instrumental
systematics. Thus, we find strong evidence that faint-object K2 light curves
show patterns that are instrumental, not astrophysical, in origin. These
systematics are not removed by any currently used analysis pipeline. The goal
of this work is to identify degrees of freedom that may prove useful for
characterizing and removing this time dependent correlated noise signature
from faint targets like AGN.

###  4.1 Channel Dependence and Spacecraft Orientation

As a next step in our analysis, we test for channel dependence of the
instrumental systematics in K2 light curves. In Figure 4, we show median light
curves (normalized flux) per channel for Campaign 8 (orange) and time-reversed
for Campaign 16 (purple). Features resembling an “M” or “W” (or other complex
patterns) that are synchronized in time reveal channel-dependent systematics
that dominate over the astrophysical variability of the sources.

Figure 4: Median channel light curves for K2 Campaigns 8 and 16. The median
light curves are calculated by computing the total-postage-stamp light curves
for all objects between 13th and 20th magnitude that are observed on a given
channel, then taking the median of those light curves. Using all of the pixels
in the postage stamp for each object results in individual light curves that
are dominated by the systematic noise found in the background pixels. Taking
the median of these light curves will result in a light curve that maintains
any pattern that is similar across the individual light curves. The median
light curves for campaign 8 (orange) and campaign 16 (purple) are plotted on
their respective channels. The channels outlined in red are labeled as being
affected the worst by rolling band according to the Campaign 0 release notes.
The campaign 16 median light curves are time-reversed to illustrate the
mirroring of the systematics seen in the two campaigns.

For a given channel in Campaign 8 or 16, we calculate the curves plotted in
Figure 4 as follows:

  * 1

We calculate all the light curves for objects within the 13–20 magnitude range
on that channel, using all of the pixels in each object’s postage stamp (also
interpolating any missing/spurious cadences).

  * 2

We normalize the individual light curves by subtracting the mean flux of the
object and dividing by the standard deviation. This normalization was done to
ensure that the median light curves would not be disproportionately affected
by brighter objects.

  * 3

For each channel, we take the median value across all the light curve
ensembles, evaluated at each cadence.

Using the total postage stamp for a large range of source magnitudes ensures
that the systematics from the background pixels will dominate at each cadence.

Note the strong correlations after time-reversal between targets that are
unique to the fields corresponding to Campaign 8 (rear facing) and Campaign 16
(forward facing). This observation is perhaps the most crucial clue to
understanding the correlated noise shown in Figure 3. It strongly suggests
that the underlying problem relates to the relative Sun-spacecraft-field
orientation which was approximately the same on day 1 of Campaign 8 and the
last day of Campaign 16 and vice versa. It further means that these
systematics are unlikely to be due to Moiré, rolling band, dust particles,
etc.—unless those are somehow dependent on Sun-spacecraft-field orientation
and are largely repeatable. Instead it is strongly suggestive of temperature-
dependent focus changes (that are further subject to differences between the
channels). This pattern may explain the observation by Armstrong et al. (2015)
and Aigrain et al. (2015) regarding the need to split campaigns into separate
time segments on either side of when the spacecraft was perpendicular to the
sun (see Howell et al. 2014 Figure 2 for a visualization of the geometry).

###  4.2 Rolling band

According to the Kepler Handbook, rolling band can be a significant source of
additive alias noise that depends on channel, local detector electronics (LDE)
board temperature, amplifier temperature, and excitation by bright sources.
The rolling band may add as much as 20 flux counts per
pixel555https://keplerscience.arc.nasa.gov/new-in-lightkurve-identifying-time-
variable-background-noise.html and has been measured to propagate across
pixels on the timescale of hours (Cleve et al., 2016) which was critically
important to Kepler’s original mission. It is also possible that increases in
flux due to rolling band may occur on timescales of days to weeks which are
much more detrimental to the K2 mission, heavily obscuring the signal of faint
or extragalactic targets. In the case of AGN, 20 counts of noise per pixel
would completely obscure the astrophysical signal if all of the source
variability is in the range of one hundred flux counts. Examples of such
ranges exist among the 21 confirmed AGN of the original Kepler mission (see
Figure 12 in Smith et al. 2018b). Any additive instrumental signature that
modulates the flux of distant AGN over time by a few tens of counts per pixel
could entirely dominate the shape of the K2 light curve.

The channels outlined in red in Figure 4 are flagged for worst rolling band
(Campaign 0 release notes; Van Cleve & Bryson 2017); however, these channels
are not necessarily the most dominated by the type of systematics that we aim
to characterize. Since rolling band flags were generated at commissioning, it
is possible that noise characterization was only exploited at timescales most
likely to compromise the original mission (detection of planet transits) and
that low frequency modes of the rolling band are poorly characterized across
each CCD. Additionally, channels may not have degraded uniformly and
temperature changes may follow very different trends in the K2 mission than
they did during pre-launch commissioning. Thus, channels that are unflagged
may still be strongly affected by rolling band.

Regardless of timescale, hours or weeks, much of the rolling band should in
theory (and as per handbook recommendations) be mitigated by careful
background subtraction. Rolling band is characteristically responsible for a
flux increase that propagates across CCD rows, but the systematics seen in
Figures 3 and 4 remain even after background subtraction (Everest light
curves, see Luger et al. 2016 section 3.1 for details of photometry). The low
frequency patterns seen in Figures 3 and 4 are either overlooked modes of
rolling band that are not sufficiently mitigated by a regular background
subtraction or they are an altogether different type of noise that turned up
post-launch, especially if they cross channel boundaries.

###  4.3 Spatial Dependence

There are a number of arguments against rolling band being the main cause of
systematic noise in our AGN sample. The main argument is that we find a
spatial dependence relative to the boresight rather than along the read-out
direction of the CCDs. We also note a time dependence of that spatial
dependence, which we focus on in this section. In Figure 5, we gather a set of
noise-dominated light curves from channel 55. The light curves are ordered by
decreasing pixel row index (Y). The sequence of dip features (starting with
three such dips in the top panel) shift by a few days with decreasing pixel
row index (ending with only two dip features in the bottom panel example).
Such a dependence may appear to be consistent with “rolling band” despite
channel 55 not being identified as having high rolling band noise. However,
the timescale of the shift is days not hours (although there exist examples of
what is thought to be rolling band on such timescales), and this shifting
sequence of dips in the light curves is seen in reverse in Campaign 16.
Furthermore, we find evidence (below) that the pattern crosses channel
boundaries and is related to boresight distance rather than local CCD
direction.

Figure 5: Campaign 8 Channel 55 objects ordered by their Y-position on the
CCD. ‘X’ is the pixel column and ‘Y’ is the pixel row. The red line marks the
position of the central dip in the light curve of the 4th object. There is a
position- (but not magnitude-) dependent phase shift in the systematics.

In Figure 6 we illustrate a slightly different pattern affecting another
channel (33). We measure the time when a peak and dip pair occurred in
Campaign 8 for all point sources greater than 17th Kepler magnitude. In the
top row, the images at the light curve peak and dip rule out neighbors,
asteroids, and extreme pointing errors as the origin of these features. In the
middle panel, we show that the peak and dip are measured as the maximum and
subsequent minimum of the light curve as marked by the blue and orange
vertical lines corresponding to the epochs of each image in the top row
respectively. In the bottom panel, we plot the dip epochs, $t_{min}$, as a
function of pixel row and column for the full sample of light curves. The dip
feature indeed shifts as a function of pixel row in this particular channel.
We test for global spatial dependence across the focal plane next by
transforming row and column pixel coordinates to radius. The diversity of
patterns, for example the number, width and amplitude of peaks and dips, makes
it more difficult to measure the time dependence and spatial dependence of
such features.

Figure 6: Pixel row and column dependence of a time-shifted dip feature
appearing in multiple light curves in channel 33. Top: Image of the PSF at the
maximum and minimum of the light curve. Middle: An example light curve, here
from channel 33 (module 11), showing the location of the maximum and minimum
features as illustrated in the images. The vertical lines indicate where we
detect the those features and where each of those images were taken. Bottom:
Timing of the minimum with pixel row (purple) and column (in grey) for channel
33. Rolling band is characteristically responsible for a flux increase that
propagates across CCD rows. The features indeed shift as a function of pixel
row and column, but this trend is a broader function of radius from the center
of the FOV as we will illustrate in the following figures. Figure 7: Timing a
sequence of two dip features in the light curves of point sources with fainter
than 16th Kepler magnitude for all channels in modules 11, 12, 13, and 14. The
scatter plots depict the epochs of two dip features per target as a function
of radius from the FOV center. Channels are colored such that the gradient
from light to dark indicates increasing radius. For example, navy points on
the left are channels 33 and 36 and on the right are channels 51 and 52. The
lack of randomness indicates that the detected minima/dips are systematic and
form trends that cross channel boundaries. The channels tested correspond to
horizontal (x-axis) modules across the FOV diagram in Figure 4. Figure 8: Same
analysis as in Figure 7 but for vertical (y-axis) modules 8, 13, 18, and 23
shown in Figure 4. Failed module 3 (channels 5–8) is shaded in grey. Time
dependence and radial dependence of these systematics show the least variance
at the outer-most radius.

Figure 4 provides a big-picture view of systematics that vary in pattern and
phase across the CCD array with no obvious positional dependence. Instrumental
trends are most likely superposed with astrophysical signal and are harder to
characterize from anything but very faint targets, especially if the noise is
additive. To measure the time-shift of peak and dip sequences such as “N”,
“M”, or “W” shaped features, we use a wavelet algorithm to time the occurrence
of several local minima per light curve, similar to our approach in Figure 6
for channel 33. Flux minima or dips due to astrophysical signal are expected
to occur randomly in time, whereas any one or more minima that occur at about
the same time for several light curves reveal an instrumental origin.

To test if the timing of sequential dips (for example a “W” feature) exhibit a
spatial dependence from channel to channel, we plot the times, $t_{min}$, for
the last two dip features per light curve as a function of radius from the FOV
center. Figures 7 and 8 show timings of the dip features for light curves
observed in modules 11 (which includes channels 33–36 as seen in Figure 4),
12, 13, 14, and 15 (channels 49–52), which span the full diameter of the focal
plane, and perpendicular modules 8 (channels 21–24), 13, 18, and 23 (channels
77–80; short of a full diameter due to failed module 3)666See
https://keplerscience.arc.nasa.gov/the-kepler-space-telescope.html. We account
for the CCD rotations and readout directions of the four channels in each
module according to rotations noted in Figure 24 of Van Cleve & Caldwell
(2016) and we track radius in units of pixel coordinates.

Figures 7 and 8 reveal a time dependence and radial dependence for systematics
that in fact cross channel boundaries. Transitioning colors indicate these
channel boundaries. On some channels this approach reveals more than two
clusters. Since clusters point to correlated dips in the light curves of
multiple targets, it is possible that the depths of these dips vary across the
light curve sample, such that for a subset of light curves the first two dips
were the most prominent “W” shape and for other light curves the last two dips
were the most prominent “W” shape. Thus, some channels appear to have up to
three or four clusters which reflects diversity in the depth width of the
instrumental signatures on that channel. However uncertain the source of the
instrumental signature, it is empirically clear that the signature crosses
channel boundaries and cannot be rolling band noise.

###  4.4 Magnitude Dependence

Figure 9: Pixel magnitudes in Campaign 8, channel 55. In the left panel, we
show all epoch-1 postage stamps. An inset that was positioned on a randomly
selected group of targets shows that very few pixels are truly empty. In the
right panel, all channel 55 pixels for all epochs are plotted as a magnitude
histogram. The distribution of pixels corresponding to the targets shown in
the inset is represented by the purple bins (the absolute counts of this
distribution are scaled for the benefit of making them visible). There are
very few truly empty pixels that can be used for computing local or global
backgrounds.

Since it is possible that the systematics are additive and thereby are more
dominant in faint sources with low counts, we explored a magnitude dependence
at the pixel level and the target level. The majority of bright K2 stellar
sources reside in the magnitude range $12<K_{p}<15$, AGNs reside in the
magnitude range $15<K_{p}<21$, and background pixels generally have magnitudes
$K_{p}>18$ depending on the size of the postage stamp and the extent of the
PSF. Faint background sources $\sim 19$th magnitude may also reside in every
pixel (Van Cleve & Caldwell, 2016).

In Figure 9, we selected a worst-case example, channel 55 (see Figure 4), from
Campaign 8 to investigate the typical distribution of pixel illumination,
measured in magnitudes. We show the distribution of Kepler magnitudes for
individual pixels in order to quantify the percentage of “empty” pixels that
may be used to characterize instrumental trends in each CCD background. The
left panel includes an inset with example postage stamps at a snapshot in
time, demonstrating that the postage stamps contain very few truly empty
pixels that can be used for computing local backgrounds. The raw magnitude
distribution of pixels in the inset is also shown in the right panel by the
purple bins. Pixels in these four postage stamps range from $21-22.2$ Kepler
magntiudes. Although nearly $40-50\%$ of pixels are fainter than the faintest
pixels in these four postage stamps, these faint sources likely have no
background pixels and their PSFs may even jitter out of the cutout boundaries.

To examine how systematic effects reveal themselves differently as a function
of target magnitude, we compute and stack full-postage-stamp flux light curves
and select the 5th, 10th, 15th, 25th, 50th, 75th, and 85th percentile at each
epoch. In Figure 10, we plot these percentile light curves in units of
fractional change and label each with its corresponding mean (target)
magnitude for channel 55 (left panel; edge of focal plane) and channel 62
(right panel; center of focal plane). Figure 10 reveals that fainter magnitude
targets have larger fractional changes than brighter targets, especially on
bad channels. Channel $55$ shows that targets with Kepler magnitudes fainter
than $18th$ exhibit amplitudes greater than $5\%$ while on better channels,
such as $62$ located at the center of the focal plane array (FPA), faint
targets exhibit fractional changes at $\sim 2\%$ or better.

Figure 10: Target magnitude dependence on channel 55 (edge of FPA) and channel
62 (center of FPA) in Campaign 8. Each panel shows flux-percentile target
light curves in units of fractional change. Percentile light curves are
labelled by their corresponding mean target magnitude. A 13th magnitude target
on channel $62$ has noise fluctuations with amplitudes less than $1\%$. On the
other hand, in channel 55, targets fainter than magnitude 18 regularly have
multiple dips at the fractional amplitude of $\sim 5\%$ or greater

.

We display the results of a related analysis for all of the channels (from
Campaigns 8 and 16) in Figure 11, where correlated systematics can be seen
across the K2 FOV. In this presentation we have stacked 500 (magnitude-binned)
median light curves generated from “empty” pixels (the faintest 80% shown in
Figure 9) to produce intensity maps for the entire CCD array for Campaigns 8
and 16 using the method developed by O’Brien et al. (2018). Non-random
structure in these intensity maps show pixel correlations (in time; x-axis in
each channel) due to the instrument as a function of magnitude (y-axis in each
channel). We applied flux smoothing and whitening (transformed to having a
mean of zero and a variance equal to one) in order to uniformly compare the
strength of systematics across magnitudes and channels. We use a divergent
purple-to-orange colormap to represent trends above and below the median flux
value, respectively. Saturated colors from purple to orange indicate a
decreasing trend in each magnitude-binned light curve as seen in Campaign 8.
For example, shifting dips show up as the orange Florida shape (channel 15) or
striping feature (channel 55) in the left-most upper block in Figure 11.
Trends in Campaign 16 are seen to be reversed in time as compared to Campaign
8 for corresponding channels.

In all campaigns (not shown here), modules 14, 15, 19, 20 are least affected
by systematics (i.e., intensity maps contain the most white), whereas modules
6, 11, 16 (on the left edge of Figure 11) exhibit the strongest intensity
systematics. Scanning each intensity map from the top-down, the direction of
increasing magnitude (dimmer pixel illumination), we observe structure
corresponding to peaks (purple) and dips (orange) in median light curves.
Strikingly, all of this correlation or pattern structure comes from pixels
with magnitudes from $19-30$ which are used for local background subtraction.
We attempted to use these magnitude-binned background light curves to correct
various stellar and AGN light curves, but it was unclear how to evaluate an
optimal fit without suffering the same failings of a PCA-approach. The
magnitude and spatial systematic dependencies discovered in these analyses
must be considered during aperture section and arc-drift removal steps.

Figure 11:  Pixel-correlation-intensity maps as a function of magnitude. Left:
Campaign 8; right: Campaign 16. Magnitude-dependent systematics in C8 are seen
to be reversed in time for C16. The $y$-axis of each channel corresponds to
magnitude bins. Each row, or magnitude-binned background light curve is
generated from an averaging of background pixels. Each row on the intensity
map is a magnitude-specific global background, stacked together to visualize
systematics across the CCD array as a function of magnitude. Purple
corresponds to normalized flux above the mean. Orange to purple indicates an
increasing trend in each magnitude-binned light curve. Fainter magnitudes
(towards the bottom of each channel) exhibit more complicated systematics.

Three inner modules exhibit a trend in Figure 11 that is opposite to the rest
of the modules in the array. Modules 12, 13, and 18 exhibit a negative slope,
purple to orange in Campaign 8, while all other modules have a positive trend,
orange to purple. This result provides evidence that focus is a dominant
source of error responsible for long-term tends. Given the spacecraft
orientation in Campaign 16, for example, outer modules may have started the
campaign in better focus and ended out of focus with a broader PSF, while
inner modules may have had the opposite trend with a narrowing PSF or vice-
versa.

Figure 12: Custom SFF reprocessing examples. Each column is a unique target
from Campaign 8 (1st and 3rd columns) or 16 (2nd and 4th columns). Each row is
a step in a typical extraction procedure. _top row:_ Raw light curves taken
from the complete postage stamp. _middle row:_ Raw light curves taken from the
4 brightest pixels in the postage stamp at any given time. _bottom row:_
(lightkurve) SFF light curve using 4 brightest pixels as aperture. The
decrease in aperture size significantly increases arc drift which is
“corrected” by the application of SFF. The smaller aperture removes the long-
term systematics we see, but the light curve is still dominated by the effects
of arc drift, even after SFF (as is clearly seen by the “mirrored” nature of
the light curves illustrated here and in Figure 4).

###  4.5 Optimizing the Aperture

Aperture photometry is a technique that assumes a useful shape such as a
square or circle rather than a PSF model to collect and sum observed counts.
After estimating the total counts from all pixels within some defined area, an
estimate of the background can be subtracted to compute the flux of the
source. An optimal aperture captures as much of the signal (as close to 100%
of the true PSF) as possible to maximize S/N. Fainter sources typically
require larger apertures. However, as the number of pixels increases, the PSF
wings of a faint source have increasingly larger errors. These considerations
are particularly relevant for AGN in crowded fields and local (low luminosity)
AGN that do not significantly outshine their host galaxy. In such cases, flux
extraction of the faint AGN becomes prohibitively more difficult since the
outer wings of the PSF are convolved with the extended galaxy profile. The K2
systematics discussed in the previous subsection present a tertiary motivation
for PSF modelling or developing a careful aperture optimization scheme.

Careful frame-by-frame aperture photometry may be a way to mitigate time-
dependent focus changes in K2 light curves. Empirically, we find that the
aperture size highly influences the strength of the instrumental signature
under investigation. Typically for AGN, we would want very small apertures to
reduce blended background sources in dense fields. We demonstrate the issue
with this approach in Figure 12 using the lightkurve package (Lightkurve
Collaboration et al., 2018)777lightkurve also contains functions for centroid
position information and PSF engineering data..

Figure 13: PLD and SFF corrections from the python lightkurve package (without
background subtraction) applied to three objects from campaign 8, channel 55.
Each column is a unique target. _top row:_ Raw light curves extracted from the
complete postage stamp. _second row:_ VJ14 SFF light curve using the full
postage stamp. _third row:_ (lightkurve) SFF light curve using the full
postage stamp. _fourth row:_ (EVEREST 2.0) PLD + GP process corrected light
curve using the full postage stamp. _bottom row:_ (lightkurve) PLD extracted
light curve without GP kernel. Instrumental systematic effects persist after
arc-drift removal in all rows.

We observe a trade-off between arc-drift correction performance and the
inclusion of complex systematics that are difficult to detrend. In Figure 12,
we compare 4-pixel moving aperture light curves for four unique targets
(columns) to full postage-stamp light curves (little arc drift but maximal
error contributed by more pixels). In the first row of Figure 12, we plot full
postage-stamp light curves for pairs of objects with mirrored systematics in
Campaigns 8 and 16. In the second row, we show the 4-pixel moving aperture for
each object. Finally, in the third row, we show the SFF-corrected moving
aperture light curves. The smaller aperture minimizes the inclusion of
instrumental trend, as seen in the second row, but the light curve is still
dominated by the effects of arc drift, even after the application of the SFF
correction in the third row.

Arc-drift requires larger apertures to be properly corrected, but because
Kepler pixels are large this ultimately includes more sources of error. An
additional difficulty for aperture photometry and PRF modelling using
commissioning data (Bryson et al., 2010) is that the true PSF may change with
time due to temperature dependent focus changes throughout the course of a
campaign as we have been suggesting. A fixed aperture size, although moving to
mitigate arc-drift, would fail to preserve a constant percentage of the true
PSF as it breathes in response to temperature changes experienced as a
consequence of the Sun-spacecraft orientation.

PSF modelling may be a crucial part of the error mitigation process for faint
targets at the largest distance from the center of the FOV. However, is
unclear how dramatic PSF changes may be throughout a campaign and if they are
radially dependent given initial focus non-uniformities. It is also possible
that the PSF may go out of focus to the point that the signal of interest
spreads beyond the edges of a postage stamp. There is no easy way, to our
knowledge, to deal with these cases.

##  5 K2 Re-processing Pipelines

Two of the most successful and popular arc-drift-corrected K2 data sets
available through MAST are the self-flat fielding (SFF) method of Vanderburg &
Johnson (2014) (henceforth VJ14) and the pixel level decorrelation (PLD)
method implemented in EVEREST by Luger et al. (2018). Both of these methods
have also been implemented in the lightkurve package (Lightkurve Collaboration
et al., 2018). In this section, we present and discuss SSF and EVEREST 2.0
public light curves, showing that these data products are successfully free of
arc-drift, but still require instrumental noise removal at low to mid-range
frequencies prior to employing time series analyses such as estimating power
spectra.

The SFF method corrects for systematics correlated with centroid arclength
(drift) by fitting a third order B-spline to position offsets over windows of
$\sim 1.5$ days across the light curve. VJ14 light curves are plotted in the
second row of Figure 13 (for multiple Kepler magnitudes). For comparison, raw
light curves are shown in the first row. The VJ14 pipeline first defines an
optimal stationary aperture by defining a circle of pixels around the target
centroid where radius is a function of magnitude. A median background is
estimated from excluded pixels and subtracted from the optimal aperture light
curve. Multiple circular and pixel response function based apertures are
available in MAST. K2 SFF light curves are continuum normalized and no flux
errors are available. A flexible implementation of the SSF method is available
in the lightkurve package.

Pixel level decorrelation (PLD; Deming et al. 2015) is another highly
successful method for arc-drift removal (Luger et al., 2018). PLD is a Taylor-
series expansion of fractional flux changes in individual pixels and pixel
products intended to model covariances in flux variations (across pixels)
resulting from instrumental motion. In EVEREST 1.0, PLD is calculated from
pixels within each individual target postage stamp. In EVEREST 2.0, PLD
weights are computed from an ensemble of neighboring stars on the same channel
(limited to bright stars in the magnitude range of $11<K_{p}<13$) and thus is
referred to as the neighboring-PLD (nPLD). EVEREST 2.0 also makes use of a
Matern 3/2 Gaussian process (GP) kernel to account for systematic covariance
between cadences when fitting the noise model. The PSD of a Matern 3/2 kernel
has an exact closed-form representation, with power law slope $-4$ in the
limit of high frequencies and constant for low frequencies, with a knee set by
the natural frequency tuning parameter $\omega_{0}$ (Foreman-Mackey et al.,
2017). The choice of a Matern 3/2 in EVEREST may therefore imprint an
artificial power law slope into the observed PSD. We experimentally verified
that these power law slopes appear in finite sample sizes by drawing random
sequences from Gaussian Process kernels with known input $\omega_{0}$ and time
sampling comparable to the K2
lightcurves888https://github.com/BrownDwarf/probabilisticAGN. EVEREST 2.0
light curves are plotted in the fourth row of Figure 13.

Figure 14: EVEREST 2.0 applied to a quiet star. The top left panel is the raw
light curve, using all of the pixels on the postage stamp as the aperture. The
bottom left panel is the light curve resulting from PLD method and a Matern
$3/2$ Gaussian process applied using the same aperture. The right panels are
the corresponding PSDs. The black dotted line represents the PSD slope of a
damped random walk (DRW or red noise). The Matern $3/2$ kernel adds in
variability that results in a PSD slope that artificially mimics a DRW. Figure
15: PSD logarithmic slopes as a function of magnitude for non-AGN objects in
Campaign 8 before (left panel; raw full-postage stamp light curve PSDs) and
after (right panel) EVEREST 2.0 reprocessing using the Matern $3/2$ kernel.
Here we see that a naive use of the Matern 3/2 kernel artificially steepens
the PSD slopes of $9483$ non-AGN targets (from various Guest Observer
proposals) to the value of $-2$ which corresponds to red noise spectra. An
incorrect GP model selected for AGN will similarly produce an artificial PSD
slope.

Naive use of the Matern 3/2 kernel results in a light curve with added
variability at low frequencies (red noise). In Figure 14, we show an example
raw light curve with little variability for which use of the Matern 3/2 kernel
results in added power at low frequencies. In the left panels, we see the raw
light curve with arc drift (top) and the EVEREST 2.0 light curve (bottom)
available through MAST. In the right panels, we plot the power spectral
density (PSD) estimated with the Lomb-Scargle periodogram and a power law
model fit above and below the noise floor. The PSD slope of the raw light
curve resembles white noise with a slope approximating zero, while the PSD of
the EVEREST 2.0 reprocessed light curve exhibits a slope equal to that
expected for AGN variability (PSD $\propto f^{-2}$), but this variability is
not real.

In the case of AGNs and low frequency stellar variables such as M dwarfs,
artificial power is added by the use of the Matern 3/2 kernel. In Figure 15,
we show the slope of the power spectrum in EVEREST 2.0 lightcurves as a
function of Kepler magnitude. Nearly all of the targets, represented by
density contours, in Figure 15, need to be reprocessed with the EVEREST 2.0
package or lightkurve PLD routine with the Matern 3/2 amplitude parameter set
to zero (or GP modeling toggled off) in order to benefit from the n-PLD
method. See also Saunders et al. (2019) who performed an assessment of the
impact of different detrending algorithms on the PSD and determined that PLD
with a careful choice of kernel is a promising method for recovering signals
that can be modeled by a DRW.

For AGN short timescale variability (and K2 data in particular), we are most
interested in estimating the high frequency slope of the empirical PSD. GP
processes require prior knowledge of an appropriate model (kernel) for the
astrophysical signal. However, Kepler/K2 observes a unique timescale regime
that has never been studied or surveyed ever before for a complete AGN sample.
Additionally, the Kepler/K2 30-minute cadence makes the empirical PSD the best
tool for informing our stochastic model selection for long-timescale ($>100$
day) variability studies. Thus there is a strong desire to rehabilitate the
light curves of AGNs observed during the K2 campaigns as discussed in the next
section.

##  6 Why Rehabilitate K2 Data?

Figure 16: The approximate locations of AGN relative to the K2 field of view.
The translations of different campaign fields of views to overlap added in
slight rotation error in the AGN positions. The coloring of the channels is
based on the median absolute deviation (MAD) of the median channel light
curves (after subtracting a linear trend). The MAD values are calculated from
the smoothed median channel light curves shown in Figure 4. The red outline
indicates a bad module. Module 4 failed after Campaign 10. The cluster of
objects in Campaign 19 are in SDSS “stripe 82” (York et al., 2000).

Kepler/K2 AGNs may offer significant insights to the overall understanding of
AGN variability acquired from investigations of ground-based light curves. PSD
(and structure function; SF) features at timescales less than 50 days have
been debated in multiple investigations of ground-based light curves. Short-
timescale AGN variability is difficult to disentangle from the effects of
gappy, irregular cadences using methods such as the Lomb-Scargle periodogram
and structure functions. K2 will be able to reveal if the AGN PSD slope in the
high frequency regime differs significantly from the low frequency regime
constrained with ground-based surveys.

The simplest stochastic model for simulating optical AGN variability is the
damped random walk (DRW). The DRW is, however, an inflexible model with a
fixed PSD power law shape (PSD $\propto\nu^{-2}$). Empirical PSD and SF
analyses of Kepler AGN from the mission’s first phase showed a range of PSD
shapes (curvature and slope transitions) (Mushotzky et al. 2011; Kasliwal et
al. 2015, Kasliwal et al. 2017; Smith et al. 2018a). Rehabilitating K2 light
curves would provide the opportunity to characterize the high-freqency PSD
slope for a large statistical sample of AGNs. In Figure 16, we show 4000 AGN
observed in Campaigns 8, 10, 12, 14, 16, 18, and 19 over-plotted on the Kepler
CCD array.

The (relatively) high frequency regime of AGN-PSDs can be investigated as a
function of luminosity, black hole mass, and accretion rate. The
characteristic frequency or timescale for a slope transition may also be
related to physical scales or black hole mass as is the case for AGN X-ray
variability (Uttley et al., 2002). Recently Simm et al. (2016), Caplar et al.
(2017), Zinn et al. (2017), Smith et al. (2018a) and Moreno et al (2021, in
prep.) find evidence of a range of behaviors in the high frequency regime of
AGN PSDs estimated with multiple methods and in multiple data sets. These
findings imply that the DRW is not a universal model of AGN variability, but
we do not necessarily understand what model is a more useful characterization.

Should a procedure for rehabilitating light curves for faint K2 targets become
feasible, Figure 16 provides key guidance with regard to the location and
potential light curve quality of all the SDSS quasars that were targeted by
our team. Clearly the targets are not uniformly distributed just as the
channels are not all of uniform photometric quality. AGN appearing in channels
outlined in cyan (light) blue may be the best choices for further
investigation, while those those AGN appearing in channels outlined in purple
and magenta, may be beyond repair.

##  7 Summary

Kepler/K2 has the potential to enable large-sample science programs for AGN,
stars, and extra-galactic time domain phenomena. These large-sample programs
are currently devastated by systematics that persist even after pointing
errors known as arc-drift are corrected. The systematics may be similar to
trends observed in the first phase of Kepler such as rolling band, Moiré
noise, faint blended background sources, and changes in point spread function
(PSF). For variability studies requiring the full light curve PSD, there is no
way to measure the performance trade-off of removing instrumental trends and
scattering the astrophysical signal. Cotrending basis vectors (or PCA
components), which were popular for detrending Kepler light curves, suffer
this particular failing because they do not orthogonally account for magnitude
or spatial information which we have shown to contribute significantly to the
shape of the systematics. Careful error mitigation is even more critical to
enable science for K2 programs at magnitudes and light curve cadence that are
unmatched by other telescopes.

We note two particularly useful characteristics of the systematics for readers
interested in rehabilitating these data. First, there is a trade-off between
arc-drift and the instrumental trends that we wish to remove. Arc-drift is
easier to remove and in some cases not present in larger aperture definitions.
However, the inclusion of more pixels in the aperture leads to stronger, more
complex instrumental trends. Second, there appear to be some modules that are
significantly less affected by systematics according to empty pixels as shown
in Figure 11. Modules 14, 15, 19, and 20 may provide the best data quality,
whereas modules 6, 11, and 16 exemplify the worst. This module quality, and by
extension data quality, is independent of campaign or telescope orientation.

The mirroring between Campaigns 8 and 16, in Figure 4, suggests that the
telescope-Sun orientation may be responsible for an increasing or decreasing
trend in some modules. The direction of the trend depends on whether the
spacecraft was forward-facing or rear-facing, but the slope and amplitude of
systematic features are consistently larger in poorer quality channels
regardless of the spacecraft orientation.

An ideal reprocessing pipeline would consider campaign, module/channel number,
magnitude, similar magnitude neighbors, aperture optimization, local
background subtraction, arc-drift correction, and distance from the boresight.
The lightkurve package provides most of the functionality to get started.

The Kepler/K2 observations are a powerful legacy dataset for investigations of
variable objects. The statistical criteria for doing AGN variability studies
are rather different from the requirements for exoplanet discovery, for which
the experiment was designed, thus some of the problems we identify above are
particularly challenging. More precise characterization and amelioration of
these systematic effects are beyond the scope of this work and, we believe,
requires in-depth investigation by teams experienced in optical
instrumentation. Nevertheless, we anticipate that the effort required to clean
Kepler/K2 light curves of the systematics we describe above will pay
worthwhile dividends not only for AGN science, but also expand the parameter
space of exoplanet detections that lie at even lower signal to noise than has
been achieved to date.

We acknowledge support from NASA grant NNX17AF18G and NASA K2 guest observer
grant 80NSSC19K0249. This paper includes data collected by the K2 mission.
Funding for the K2 mission is provided by the NASA Science Mission
directorate. M.S.V. acknowledges support from the Ambrose Mondell Foundation
during sabbatical leave at the Institute for Advanced Study. We thank Michael
Gully-Santiago, Nicholas Saunders, Gert Barentsen, Patricia Boyd, and Thomas
Barclay for their contributions to the discussion about this work.

## References

  * Aigrain et al. (2015) Aigrain, S., Hodgkin, S. T., Irwin, M. J., Lewis, J. R., & Roberts, S. J. 2015, MNRAS, 447, 2880, doi: 10.1093/mnras/stu2638
  * Annis et al. (2014) Annis, J., Soares-Santos, M., Strauss, M. A., et al. 2014, ApJ, 794, 120, doi: 10.1088/0004-637X/794/2/120
  * Aranzana et al. (2018) Aranzana, E., Körding, E., Uttley, P., Scaringi, S., & Bloemen, S. 2018, MNRAS, 476, 2501, doi: 10.1093/mnras/sty413
  * Armstrong et al. (2015) Armstrong, D. J., Kirk, J., Lam, K. W. F., et al. 2015, A&A, 579, A19, doi: 10.1051/0004-6361/201525889
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002, doi: 10.1088/1538-3873/aaecbe
  * Bovy et al. (2011) Bovy, J., Hennawi, J. F., Hogg, D. W., et al. 2011, ApJ, 729, 141, doi: 10.1088/0004-637X/729/2/141
  * Bryson et al. (2010) Bryson, S. T., Tenenbaum, P., Jenkins, J. M., et al. 2010, The Astrophysical Journal, 713, L97, doi: 10.1088/2041-8205/713/2/l97
  * Caplar et al. (2017) Caplar, N., Lilly, S. J., & Trakhtenbrot, B. 2017, ApJ, 834, 111, doi: 10.3847/1538-4357/834/2/111
  * Chen & Wang (2015) Chen, X.-Y., & Wang, J.-X. 2015, ApJ, 805, 80, doi: 10.1088/0004-637X/805/1/80
  * Cleve et al. (2016) Cleve, J. E. V., Howell, S. B., Smith, J. C., et al. 2016, Publications of the Astronomical Society of the Pacific, 128, 075002, doi: 10.1088/1538-3873/128/965/075002
  * Crossfield et al. (2015) Crossfield, I. J. M., Petigura, E., Schlieder, J. E., et al. 2015, ApJ, 804, 10, doi: 10.1088/0004-637X/804/1/10
  * Deming et al. (2015) Deming, D., Knutson, H., Kammer, J., et al. 2015, The Astrophysical Journal, 805, 132, doi: 10.1088/0004-637x/805/2/132
  * Dobrotka et al. (2017) Dobrotka, A., Antonuccio-Delogu, V., & Bajčičáková, I. 2017, MNRAS, 470, 2439, doi: 10.1093/mnras/stx961
  * Dobrotka et al. (2019) Dobrotka, A., Bezák, P., Revalski, M., & Strémy, M. 2019, MNRAS, 483, 38, doi: 10.1093/mnras/sty3074
  * Edelson & Malkan (2012) Edelson, R., & Malkan, M. 2012, ApJ, 751, 52, doi: 10.1088/0004-637X/751/1/52
  * Edelson et al. (2014) Edelson, R., Vaughan, S., Malkan, M., et al. 2014, ApJ, 795, 2, doi: 10.1088/0004-637X/795/1/2
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220, doi: 10.3847/1538-3881/aa9332
  * Foreman-Mackey et al. (2015) Foreman-Mackey, D., Montet, B. T., Hogg, D. W., et al. 2015, ApJ, 806, 215, doi: 10.1088/0004-637X/806/2/215
  * Giustini & Proga (2019) Giustini, M., & Proga, D. 2019, A&A, 630, A94, doi: 10.1051/0004-6361/201833810
  * Howell et al. (2014) Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398, doi: 10.1086/676406
  * Huang et al. (2015) Huang, C. X., Penev, K., Hartman, J. D., et al. 2015, MNRAS, 454, 4159, doi: 10.1093/mnras/stv2257
  * Ivezić et al. (2008) Ivezić, Ž., Axelrod, T., Becker, A. C., et al. 2008, in American Institute of Physics Conference Series, ed. C. A. L. Bailer-Jones, Vol. 1082, 359–365, doi: 10.1063/1.3059076
  * Kasliwal et al. (2017) Kasliwal, V. P., Vogeley, M. S., & Richards, G. T. 2017, MNRAS, 470, 3027, doi: 10.1093/mnras/stx1420
  * Kasliwal et al. (2015) Kasliwal, V. P., Vogeley, M. S., Richards, G. T., Williams, J., & Carini, M. T. 2015, MNRAS, 453, 2075, doi: 10.1093/mnras/stv1797
  * Kinemuchi et al. (2012) Kinemuchi, K., Barclay, T., Fanelli, M., et al. 2012, PASP, 124, 963, doi: 10.1086/667603
  * Libralato et al. (2015) Libralato, M., Bedin, L. R., Nardiello, D., & Piotto, G. 2015, Monthly Notices of the Royal Astronomical Society, 456, 1137, doi: 10.1093/mnras/stv2628
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python, Astrophysics Source Code Library.  http://ascl.net/1812.013
  * Luger et al. (2016) Luger, R., Agol, E., Kruse, E., et al. 2016, The Astronomical Journal, 152, 100, doi: 10.3847/0004-6256/152/4/100
  * Luger et al. (2018) Luger, R., Kruse, E., Foreman-Mackey, D., Agol, E., & Saunders, N. 2018, AJ, 156, 99, doi: 10.3847/1538-3881/aad230
  * Lund et al. (2015) Lund, M. N., Handberg, R., Davies, G. R., Chaplin, W. J., & Jones, C. D. 2015, ApJ, 806, 30, doi: 10.1088/0004-637X/806/1/30
  * McCalmont et al. (2015) McCalmont, K. M., Larson, K. A., Peterson, C. A., & Ross, S. E. 2015
  * Moreno et al. (2019) Moreno, J., Vogeley, M. S., Richards, G. T., & Yu, W. 2019, Publications of the Astronomical Society of the Pacific, 131, 063001, doi: 10.1088/1538-3873/ab1597
  * Mushotzky et al. (2011) Mushotzky, R. F., Edelson, R., Baumgartner, W., & Gandhi, P. 2011, ApJ, 743, L12, doi: 10.1088/2041-8205/743/1/L12
  * O’Brien et al. (2018) O’Brien, J. T., Moreno, J., Richards, G. T., & Vogeley, M. S. 2018, Research Notes of the American Astronomical Society, 2, 127, doi: 10.3847/2515-5172/aad331
  * Pâris et al. (2017) Pâris, I., Petitjean, P., Ross, N. P., et al. 2017, A&A, 597, A79, doi: 10.1051/0004-6361/201527999
  * Peters et al. (2015) Peters, C. M., Richards, G. T., Myers, A. D., et al. 2015, ApJ, 811, 95, doi: 10.1088/0004-637X/811/2/95
  * Revalski et al. (2014) Revalski, M., Nowak, D., Wiita, P. J., Wehrle, A. E., & Unwin, S. C. 2014, ApJ, 785, 60, doi: 10.1088/0004-637X/785/1/60
  * Richards et al. (2004) Richards, G. T., Nichol, R. C., Gray, A. G., et al. 2004, ApJS, 155, 257, doi: 10.1086/425356
  * Richards et al. (2015) Richards, G. T., Myers, A. D., Peters, C. M., et al. 2015, ApJS, 219, 39, doi: 10.1088/0067-0049/219/2/39
  * Saunders et al. (2019) Saunders, N., Luger, R., & Barnes, R. 2019, AJ, 157, 197, doi: 10.3847/1538-3881/ab12e4
  * Schneider et al. (2010) Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, AJ, 139, 2360, doi: 10.1088/0004-6256/139/6/2360
  * Shaya et al. (2015) Shaya, E. J., Olling, R., & Mushotzky, R. 2015, AJ, 150, 188, doi: 10.1088/0004-6256/150/6/188
  * Simm et al. (2016) Simm, T., Salvato, M., Saglia, R., et al. 2016, A&A, 585, A129, doi: 10.1051/0004-6361/201527353
  * Smith et al. (2012) Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000, doi: 10.1086/667697
  * Smith et al. (2018a) Smith, K. L., Mushotzky, R. F., Boyd, P. T., et al. 2018a, ApJ, 857, 141, doi: 10.3847/1538-4357/aab88d
  * Smith et al. (2018b) Smith, K. L., Mushotzky, R. F., Boyd, P. T., & Wagoner, R. V. 2018b, ApJ, 860, L10, doi: 10.3847/2041-8213/aac88c
  * Smith et al. (2015) Smith, K. L., Boyd, P. T., Mushotzky, R. F., et al. 2015, AJ, 150, 126, doi: 10.1088/0004-6256/150/4/126
  * Still & Barclay (2012) Still, M., & Barclay, T. 2012, PyKE: Reduction and analysis of Kepler Simple Aperture Photometry data, Astrophysics Source Code Library.  http://ascl.net/1208.004
  * Stumpe et al. (2012) Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985, doi: 10.1086/667698
  * Tamuz et al. (2005) Tamuz, O., Mazeh, T., & Zucker, S. 2005, MNRAS, 356, 1466, doi: 10.1111/j.1365-2966.2004.08585.x
  * Uttley et al. (2002) Uttley, P., McHardy, I. M., & Papadakis, I. E. 2002, MNRAS, 332, 231, doi: 10.1046/j.1365-8711.2002.05298.x
  * Van Cleve & Bryson (2017) Van Cleve, J. E., & Bryson, S. T. 2017, K2 Handbook, Tech. rep.
  * Van Cleve & Caldwell (2016) Van Cleve, J. E., & Caldwell, D. A. 2016, Kepler Instrument Handbook, Tech. rep.
  * Vanderburg & Johnson (2014) Vanderburg, A., & Johnson, J. A. 2014, PASP, 126, 948, doi: 10.1086/678764
  * Véron-Cetty & Véron (2010) Véron-Cetty, M.-P., & Véron, P. 2010, A&A, 518, A10, doi: 10.1051/0004-6361/201014188
  * Wehrle et al. (2013) Wehrle, A. E., Wiita, P. J., Unwin, S. C., et al. 2013, ApJ, 773, 89, doi: 10.1088/0004-637X/773/2/89
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513
  * Zinn et al. (2017) Zinn, J. C., Kochanek, C. S., Kozłowski, S., et al. 2017, MNRAS, 468, 2189, doi: 10.1093/mnras/stx586


# Joint nonstationary blind source separation and spectral analysis

Adrien Meynard1
1Aix Marseille Univ, CNRS, Centrale Marseille, I2M, Marseille, France.

###### Abstract

We address a nonstationary blind source separation (BSS) problem. The model
includes both nonstationary sources and mixing. Therefore, we introduce an
algorithm for joint BSS and estimation of stationarity-breaking deformations
and spectra. Finally, its performances are evaluated on a synthetic example.

##  1 Introduction: model and background

The BSS problem, originally introduced in a stationary context, has also been
discussed in nonstationary situations. Extensions to nonstationary signals
have been proposed, based on time-frequency analysis (see [1], chap. 9 in [2]
and references therein), or based on mutual information [3]. The BSS of a
nonstationary mixtures of stationary signals have also been studied. For
instance, in [4], the authors explore the convolutive BSS problem. In the
following, we tackle a doubly nonstationary BSS problem, and propose a
demultiplexing algorithm adapted to a specific class of nonstationary signals
mixed by a instantaneous nonstationary mixing matrix.

###  1.1 Nonstationarity

The nonstationary signals of interest here are deformed versions of stationary
signals.

Let $x$ denote a stationary signal, modeled as a realization of a stationary
random process with power spectrum denoted by $\mathscr{S}_{X}$. Acting on $x$
with a stationarity-breaking operator yields a nonstationary signal denoted by
$y$. Various classes of stationarity-breaking operators are relevant to model
physical phenomena (e.g. frequency modulation [5], amplitude modulation [6]).
We focus here on the time warping operator denoted by $\mathcal{D}_{\gamma}$
and defined by:

$y(t)=\mathcal{D}_{\gamma}x(t)=\sqrt{\gamma^{\prime}(t)}\,x(\gamma(t))\ ,$ (1)

where $\gamma\in C^{2}$ is a strictly increasing smooth function. Such
deformations can model nonstationary physical phenomena as diverse as Doppler
effect, speed variations of an engine, animal vocalization or speech [7, 6].

The wavelet transform is a natural tool to analyze such signals. Hence, the
wavelet transform $\mathcal{W}_{x}$ of the signal $x$ is defined by:

$\mathcal{W}_{x}(s,\tau)=\int_{\mathbb{R}}x(t)q^{-s/2}\overline{\psi}\left(\dfrac{t-\tau}{q^{s}}\right)dt\quad\text{with}\quad
q>1\ .$ (2)

In that framework, it can be shown that the respective wavelet transforms
$\mathcal{W}_{y}$ and $\mathcal{W}_{x}$ of $y$ and $x$ are approximately
related by

$\mathcal{W}_{y}(s,\tau)\approx\mathcal{W}_{x}(s+\log_{q}(\gamma^{\prime}(\tau)),\gamma(\tau))\
.$ (3)

In the following, we make the assumption that $x$ is a realization of a
stationary random process $X$. In such a setting, the approximation error can
be controlled thanks to the decay properties of the wavelet $\psi$, and the
variations of $\gamma^{\prime}$. In [5, 6], corresponding quantitative error
bounds are given.

###  1.2 Blind source separation

The problem we consider is the BSS of nonstationary signals modeled by
equation (1).

We investigate the case where the number of sources and the number of
observations are equal and denoted by $N$. The sources are additionally
assumed to be independent. Let
${\mathbf{y}}(t),{\mathbf{z}}(t)\in\mathbb{R}^{N}$ denote the column vectors
containing respectively all the sources and observations at time $t$. Then,
the mixture is written as

${\mathbf{z}}(t)={\mathbf{A}}(t){\mathbf{y}}(t)\ ,$ (4)

where ${\mathbf{A}}(t)\in\mathbb{R}^{N\times N}$ denotes the time varying
mixing matrix, assumed to be invertible. This model generalizes the amplitude
modulation model in the case $N=1$ detailed in [6]. For example, this model
can be appropriate in bioacoustics to describe the BSS of a howling wolf pack
[8, 9].

Our goal is to determine jointly the mixing matrix ${\mathbf{A}}(t)$, the time
warping functions $\gamma_{i}(t)$, and the spectra of the stationary sources
$\mathscr{S}_{X_{i}}$ for $i=1,\ldots,N$ from the observations
${\mathbf{z}}(t)$.

Let us consider a fixed time $\tau$, then for each observation $z_{i}$, we
denote by ${\mathbf{w}}_{z_{i},\tau}=\mathcal{W}_{z_{i}}({\mathbf{s}},\tau)$
the row vector containing the values of the wavelet transform for a vector of
scales ${\mathbf{s}}$ (of size denoted by $M_{s}$). Then, all these vectors
are gathered into a $N\times M_{s}$ matrix ${\mathbf{w}}_{{\mathbf{z}},\tau}$
such that
${\mathbf{w}}_{{\mathbf{z}},\tau}=\left({\mathbf{w}}_{z_{1},\tau}^{T}\cdots{\mathbf{w}}_{z_{N},\tau}^{T}\right)^{T}$.
The same operation is applied to the wavelet transform of the sources. The
matrix ${\mathbf{A}}(t)$ is assumed to vary slowly with respect to the
oscillations of the signals. It can be shown that the linear relation (4)
becomes in this new setting a relationship between the wavelet transforms of
${\mathbf{y}}$ and ${\mathbf{z}}$ of the form

${\mathbf{w}}_{{\mathbf{z}},\tau}\approx{\mathbf{A}}(\tau){\mathbf{w}}_{{\mathbf{y}},\tau}\
.$ (5)

Aside from the terms controlling the error bound in (3), the error bound in
(5) is also controlled by the variations of the mixing matrix coefficients.

##  2 Estimation procedure

Approximation equations (3) and (5) allow us to write an approximate
likelihood in the Gaussian case (see [10] for more details on this approach).

The estimation procedure is based upon discrete wavelet transforms, time-
varying parameters are therefore estimated on a discrete time grid $D$. In the
following, the estimation procedure is described for a given $\tau\in D$. For
the sake of simplicity, we introduce the following notations:
${\mathbf{B}}_{\tau}={\mathbf{A}}(\tau)^{-1}$,
$\theta_{i,\tau}=\log_{q}\left(\gamma^{\prime}_{i}(\tau)\right)$ and
$\boldsymbol{\theta}_{\tau}=(\theta_{1,\tau}\cdots\theta_{N,\tau})^{T}$.

###  2.1 Probabilistic setting

It follows from the Gaussianity assumption on $X$ that
${\mathbf{w}}_{y_{i},\tau}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma}_{i}(\theta_{i,\tau}))$,
where

$\left[\boldsymbol{\Sigma}_{i}(\theta_{i,\tau})\right]_{kk^{\prime}}=q^{\frac{s_{k}+s_{k^{\prime}}}{2}}\\!\\!\int_{\mathbb{R}}\\!\mathscr{S}_{X_{i}}(q^{-\theta_{i,\tau}}\xi)\overline{\hat{\psi}}\left(q^{s_{k}}\xi\right)\hat{\psi}\left(q^{s_{k^{\prime}}}\xi\right)d\xi.$

Let $p_{V}$ denote generically the probability density function of a random
vector $V$. Then, the source independence hypothesis gives the following
opposite of the log-likelihood:

$\displaystyle\ell_{\tau}({\mathbf{B}}_{\tau},\boldsymbol{\theta}_{\tau})\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle-\log(p_{{\mathbf{w}}_{z,\tau}|({\mathbf{B}}_{\tau},\boldsymbol{\theta}_{\tau})}({\mathbf{w}}_{{\mathbf{z}},\tau};{\mathbf{B}}_{\tau},\boldsymbol{\theta}_{\tau}))+c$
$\displaystyle=$ $\displaystyle-
M_{s}\log|\det({\mathbf{B}}_{\tau})|\\!+\\!\dfrac{1}{2}\sum_{i=1}^{N}\log\left|\det\boldsymbol{\Sigma}_{i}(\theta_{i,\tau}))\right|$
$\displaystyle+\dfrac{1}{2}\sum_{i=1}^{N}[{\mathbf{B}}_{\tau}{\mathbf{w}}_{{\mathbf{z}},\tau}]_{i\cdot}\boldsymbol{\Sigma}_{i}(\theta_{i,\tau})^{-1}[{\mathbf{B}}_{\tau}{\mathbf{w}}_{{\mathbf{z}},\tau}]_{i\cdot}^{H}\
,$

where $[{\mathbf{M}}]_{i\cdot}$ denotes the $i$-th line of the matrix
${\mathbf{M}}$, and ${\mathbf{M}}^{H}$ is its conjugate transpose. Maximum
likelihood (ML) estimates, i.e. minimizers of
$\ell_{\tau}({\mathbf{B}}_{\tau},\boldsymbol{\theta}_{\tau})$, can be
evaluated numerically.

However, in order to take into account the smoothness assumption on the mixing
matrix with respect to time, we switch to the Bayesian framework and introduce
a prior $p_{{\mathbf{B}}_{\tau}}$ on the unmixing matrix ${\mathbf{B}}_{\tau}$
(assuming i.i.d. matrix coefficients). We choose for $p_{{\mathbf{B}}_{\tau}}$
a uniform distribution centered on ${\mathbf{B}}_{\tau-\Delta_{\tau}}$, and
with support $2\epsilon_{B}\Delta_{\tau}$. Then, the maximum a posteriori
(MAP) estimate $\tilde{\mathbf{B}}_{\tau}$ can be written as the solution of
the problem

$\tilde{\mathbf{B}}_{\tau}=\arg\min_{{\mathbf{B}}_{\tau}}\
\ell_{\tau}({\mathbf{B}}_{\tau},\boldsymbol{\theta}_{\tau})\quad\text{s.t.}\quad\|{\mathbf{B}}_{\tau}-{\mathbf{B}}_{\tau-\Delta_{\tau}}\|_{\infty}\leq\epsilon_{B}\Delta_{\tau}.$
(6)

This problem is consistent with the smoothness hypothesis on
${\mathbf{B}}_{\tau}$. Indeed, assuming $\Delta_{\tau}$ is small, the
constraint in equation (6) is almost equivalent to
$\|{\mathbf{B}}_{\tau}^{\prime}\|_{\infty}\leq\epsilon_{B}$.

Concerning the time warping estimation, we choose not to give a prior on
$\boldsymbol{\theta}_{\tau}$. Thus, $\tilde{\boldsymbol{\theta}}_{\tau}$ is
the ML estimation of $\boldsymbol{\theta}_{\tau}$.

###  2.2 Estimation algorithm

The estimation strategy is to alternate the estimations of
${\mathbf{B}}_{\tau}$, $\boldsymbol{\theta}_{\tau}$ and the spectra. The
algorithm 1 (named JEFAS-BSS) synthesizes all the estimation steps which are
described below.

  * •

Mixing matrix estimation. In practice, we numerically solve the problem (6).
Besides, because of the assumption of slow variations of the matrix
coefficients, we make the approximation that ${\mathbf{B}}_{\tau}$ is constant
on the interval $I_{\tau}=[\tau-\Delta_{\tau}/2,\ \tau+\Delta_{\tau}/2[$.
Finally, the estimated sources $\tilde{\mathbf{y}}_{\tau}$ are obtained via
$\tilde{\mathbf{y}}_{\tau}(t)=\tilde{\mathbf{B}}_{\tau}{\mathbf{z}}(t)$ where
$t\in I_{\tau}$. Notice that for each interval $I_{\tau}$, a new matrix
$\tilde{\mathbf{B}}_{\tau}$ is applied to the observations. Due to the source
ordering indeterminacy, a reordering method has to be introduced to connect
consecutive segments of each source signal. We use for that the Gale-Shapley
stable marriage algorithm [11] which constructs stable matchings between
consecutive time slices source estimations. The ranking criterion is based on
the comparison of the dot products between normalized Fourier spectra of these
slices.

  * •

Deformations and spectra estimations. For each source, the joint estimation of
$\left\\{\theta_{i,\tau}\right\\}_{\tau\in D}$ and $\mathscr{S}_{X_{i}}$ is
obtained via the JEFAS algorithm (which is detailed in [6]). For this purpose,
the input wavelet transform ${\mathbf{w}}_{{\mathbf{y}}}$ of the source
$y_{i}$ is replaced with its estimate
$\left\\{{\mathbf{B}}_{\tau}{\mathbf{w}}_{{\mathbf{z}},\tau}\right\\}_{\tau\in
D}$.

Regarding initialization, a basic method is to use a stationary BSS method on
observations to obtain a first unmixing matrix estimate. For instance, SOBI
[12] is a stationary BSS algorithm which can give an initial unmixing matrix.
A better initial matrix can be obtained by piecewise SOBI estimates on non
overlapping segments (called p-SOBI), where the stationarity assumption makes
more sense.

Algorithm 1 JEFAS-BSS

Initialization: Obtain $\tilde{\mathbf{B}}^{(0)}_{\tau}$ by means of the
p-SOBI algorithm. Compute the estimated source
$\tilde{\mathbf{y}}^{(0)}(\tau)=\tilde{\mathbf{B}}^{(0)}_{\tau}{\mathbf{z}}(\tau)$.

$\bullet$ $k\leftarrow 1$

while stopping criterion is false and $k\leq k_{max}$ do

$\bullet$ For $i=1,\ldots,N$, estimate parameters
$\tilde{\theta}^{(k)}_{i,\tau},\ \forall\tau\in D$ and spectrum
$\tilde{\mathscr{S}}_{X_{i}}^{(k)}$ applying JEFAS algorithm to
$\tilde{y}_{i}^{(k-1)}$.

for $\tau=0,\Delta_{\tau},\ldots,T$ do

$\bullet$ Estimate $\tilde{\mathbf{B}}^{(k)}_{\tau}$: solve (6) replacing
$\boldsymbol{\theta}_{\tau}$ and $\mathscr{S}_{X}$ with their current
estimations $\tilde{\boldsymbol{\theta}}_{\tau}^{(k)}$ and
$\left\\{\tilde{\mathscr{S}}_{X_{i}}^{(k)}\right\\}_{i=1,\ldots,N}$.

end for

$\bullet$ Estimate the sources $\tilde{\mathbf{y}}^{(k)}$.

$\bullet$ $k\leftarrow k+1$

end while

The convergence is monitored using the Source to Interference Ratio (SIR)
introduced in [13]. For a given estimated source, SIR quantifies the presence
of interferences from the other true sources. As we do not have access to the
ground truth sources, we use as stopping criterion the SIR between
$\tilde{\mathbf{y}}^{(k-1)}$ and $\tilde{\mathbf{y}}^{(k)}$ (instead of
${\mathbf{y}}$) which gives an evaluation of the BSS update, and is therefore
a relevant convergence assessment.

##  3 Results

We construct a synthetic example to evaluate the performances of the
algorithm. The two sources are band-pass filtered white noise, with time-
varying bandwidth. The mixing matrix coefficients are sinusoidally varying
over time. On the left of figure 1, the wavelet transforms of both
observations are displayed.

The evolution of the convergence criterion through iterations of JEFAS-BSS is
displayed in figure 1 (top-right). We can empirically note that our algorithm
converges in a small number of iterations. Indeed, after 15 iterations the
convergence criterion is around 100 dB meaning the BSS update is negligible.

Finally, we evaluate the performances of the BSS algorithms (we refer to [6]
for the evaluation of the performances of the deformations and spectra
estimations). The Amari index [14] is a measure of divergence between the
matrix $\tilde{\mathbf{B}}_{\tau}{\mathbf{A}}_{\tau}$ and the identity matrix.
The closer to zero the Amari index the better. On the bottom-right of figure
1, we display the evolution of the Amari index through time for each BSS
algorithm. In table 1, we also compare the SIR, the SDR (Source to Distortion
Ratio [13]), and the time-averaged Amari index of the BSS algorithms. Those
different criteria show that BSS-JEFAS performances are higher than those of
SOBI and p-SOBI. Besides, in average, p-SOBI gives a better Amari index than
SOBI, which is understandable because it takes into account the
nonstationarity of the mixing matrix. Nonetheless, the SIR and SDR of p-SOBI
are worse than those of SOBI. Indeed, because this method does not take into
account the regularity of ${\mathbf{B}}_{\tau}$, the connections between
slices are sensitive to discontinuities and create distortion in the estimated
sources.

Figure 1: Synthetic example. Left: scalograms of the observations. Top-right: Convergence criterion evolution. Bottom-right: Amari index evolution. Criterion | SOBI | p-SOBI | JEFAS-BSS
---|---|---|---
SIR (dB) | $28.55$ | $15.04$ | $46.55$
SDR (dB) | $16.60$ | $-4.53$ | $37.69$
Amari index | $4.63\times 10^{-2}$ | $1.74\times\\!10^{-2}$ | $1.40\times 10^{-4}$
Table 1: Comparison of the performances between BSS algorithms: standard SOBI,
piecewise SOBI and the proposed algorithm.

## References

  * [1] A. Belouchrani, M. G. Amin, N. Thirion-Moreau, and Y. D. Zhang, “Source separation and localization using time-frequency distributions: An overview,” IEEE Signal Processing Magazine, vol. 30, pp. 97–107, Nov 2013\.
  * [2] C. Jutten and P. Comon, Séparation de sources 2. Au-delà de l’aveugle et applications. Hermes, 2007.
  * [3] D.-T. Pham and J.-F. Cardoso, “Blind separation of instantaneous mixtures of nonstationary sources,” IEEE Transactions on Signal Processing, vol. 49, pp. 1837–1848, Sep 2001.
  * [4] L. Parra and C. Spence, “Convolutive blind separation of non-stationary sources,” IEEE Transactions on Speech and Audio Processing, vol. 8, pp. 320–327, May 2000.
  * [5] A. Meynard and B. Torrésani, “Spectral estimation for non-stationary signal classes,” in Sampling Theory and Applications, Proceedings of SampTA17, (Tallinn, Estonia), July 2017.
  * [6] A. Meynard and B. Torresani, “Spectral analysis for nonstationary audio,” IEEE/ACM Transactions on Audio, Speech, and Language Processing, 2018.
  * [7] D. Stowell, Computational Analysis of Sound Scenes and Events, ch. Computational Bioacoustic Scene Analysis, pp. 303–333.  Springer, 2018.
  * [8] D. Passilongo, L. Mattioli, E. Bassi, L. Szabó, and M. Apollonio, “Visualizing sound: counting wolves by using a spectral view of the chorus howling,” Frontiers in Zoology, vol. 12, p. 22, Sep 2015.
  * [9] M. Papin, J. Pichenot, F. Guérold, and E. Germain, “Acoustic localization at large scales: a promising method for grey wolf monitoring,” Frontiers in Zoology, vol. 15, p. 11, Apr 2018.
  * [10] J.-F. Cardoso, “Blind signal separation: statistical principles,” Proceedings of the IEEE, vol. 86, pp. 2009–2025, Oct 1998.
  * [11] D. Gale and L. S. Shapley, “College admissions and the stability of marriage,” The American Mathematical Monthly, vol. 69, pp. 9–15, January 1962.
  * [12] A. Belouchrani, K. Abed-Meraim, J.-F. Cardoso, and E. Moulines, “A blind source separation technique using second-order statistics,” IEEE Transactions on Signal Processing, vol. 45, pp. 434–444, Feb 1997.
  * [13] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” IEEE Transactions on Audio, Speech, and Language Processing, vol. 14, pp. 1462–1469, July 2006.
  * [14] S.-i. Amari, A. Cichocki, and H. H. Yang, “A new learning algorithm for blind signal separation,” Advances in Neural Information Processing Systems, vol. 8, pp. 757–763, 1996.