# Identifying Periodic Variable Stars and Eclipsing Binary Systems with Long-
Term Las Cumbres Observatory Photometric Monitoring of ZTF J0139+5245

Aniket Sanghi The University of Texas at Austin, Department of Astronomy, 2515
Speedway, C1400, Austin, TX 78712, USA  Zachary P. Vanderbosch The University
of Texas at Austin, Department of Astronomy, 2515 Speedway, C1400, Austin, TX
78712, USA  McDonald Observatory, Fort Davis, TX 79734, USA  Michael H.
Montgomery The University of Texas at Austin, Department of Astronomy, 2515
Speedway, C1400, Austin, TX 78712, USA  McDonald Observatory, Fort Davis, TX
79734, USA

(Received 2021 May 07; Revised 2021 July 08; Accepted 2021 July 27)

###### Abstract

We present the results of our search for variable stars using the long-term
Las Cumbres Observatory (LCO) monitoring of white dwarf ZTF J0139+5245 with
the two 1.0-m telescope nodes located at McDonald Observatory using the
Sinistro imaging instrument. In this search, we find 38 variable sources, of
which 27 are newly discovered or newly classified (71%) based on comparisons
with previously published catalogs, thereby increasing the number of
detections in the field-of-view under consideration by a factor of $\approx$
2.5. We find that the improved photometric precision per-exposure due to
longer exposure time for LCO images combined with the greater time-sampling of
LCO photometry enables us to increase the total number of detections in this
field-of-view. Each LCO image covers a field-of-view of $26^{\prime}\times
26^{\prime}$ and observes a region close to the Galactic plane
($b=-9.4^{\circ}$) abundant in stars with an average stellar density of
$\approx 8$ arcmin-2. We perform aperture photometry and Fourier analysis on
over 2000 stars across 1560 LCO images spanning 537 days to find 28 candidate
BY Draconis variables, 3 candidate eclipsing binaries of type EA, and 7
candidate eclipsing binaries of type EW. In assigning preliminary
classifications to our detections, we demonstrate the applicability of the
Gaia color-magnitude diagram (CMD) as a powerful classification tool for
variable star studies.

Time domain astronomy (2109) — Periodic variable stars (1213) — Surveys (1671)
— Eclipsing binary stars (444) — BY Draconis stars (190)

##  1 Introduction

Time-domain astronomy exploits the photometric variability of astronomical
sources to probe their underlying physical mechanisms as well as details of
interactions with other objects. Variable stars have profound implications in
astronomy and are a key driver of research beyond the realm of time-domain
astronomy. Pulsating stars such as Cepheids and RR Lyrae serve as “standard
candles” and can be used for distance determinations on the cosmic scale
(Leavitt & Pickering, 1912; Pietrzyński et al., 2013, 2019; Riess et al.,
2018). Stellar surface inhomogeneities on rotational variables, such as BY
Draconis variable stars, can be used to infer their rotational periods to
study stellar angular momentum evolution in large samples (Lanzafame et al.,
2018). The orbital kinematics of binary star systems allow the determination
of companion masses, which in turn can be utilized to indirectly estimate
other fundamental stellar parameters (Torres et al., 2010). Variable stars
have also been extensively used as tracers of the structure, kinematics,
chemical composition, and evolution of the Milky Way and other nearby galaxies
(e.g., Skowron et al., 2019; Chen et al., 2019; Jacyszyn-Dobrzeniecka et al.,
2016, 2017).

The above applications are enabled by a revolutionary transformation underway
in the field of astronomy — a collective move towards the era of open-access
data. Combining this openness with an improvement in the depth and efficiency
of time-domain surveys, specifically, the use of wide-field CCD imagers, has
led to significant advancements in the field, with the number of detected
variables increasing by a factor of $\approx 2$ every year. Past surveys have
generated large amounts of data used to fuel variable star discoveries and
analyses. Over a period of nearly 20 years, the Optical Gravitational Lensing
Experiment (OGLE: Udalski et al., 1994) has detected more than 900,000
variables in the Magellanic Clouds, the Galactic bulge, and the Galactic
plane, revolutionizing the study of periodic stars and eclipsing systems. More
recently, Chen et al. (2020) have leveraged the Zwicky Transient Facility’s
(ZTF: Masci et al., 2018) large field-of-view and faint limiting magnitude to
detect more than 700,000 variable sources in the northern sky. Further, a
total of nearly 200,000 variable stars have been discovered by surveys such as
the All-Sky Automated Survey for Supernovae (ASAS-SN: Shappee et al., 2014),
the Catalina Sky Survey (CSS: Larson et al., 2003), the Massive Compact Halo
Object survey (MACHO: Alcock et al., 1993), Pan-STARRS1 (Chambers et al.,
2016; Flewelling et al., 2020; Magnier et al., 2020a, b, c; Waters et al.,
2020), and the Gaia mission (Gaia Collaboration et al., 2016)

Figure 1: SDSS color image of the sky representing the coverage and field-of-
view (FOV) of our LCO image frame (red) and the ZTF survey (dark blue). We
also overplot the region covered by our query of the Pan-STARRS1 DR2 catalogue
(cyan). The center of the LCO image frame is marked with a yellow star. Our
images cover a small patch of sky that falls within a ZTF chip gap and is thus
previously unexplored from the time-domain perspective beyond very sparsely
sampled surveys like Pan-STARRS1.

Uniquely positioned in this data-driven revolution in time-domain astrophysics
is the Las Cumbres Observatory (LCO). LCO utilizes a global distribution of
observational facilities purpose-built to study transient and periodically
variable objects at optical and near-IR wavelengths (Brown et al., 2013) by
implementing uniform instrumentation across its network and operating fully-
robotically around-the-clock. It provides near continuous coverage of the
night sky in both the northern and southern hemispheres, and utilizes a
dynamic observation scheduling system (Saunders et al., 2014) to enable
extensive photometric and spectroscopic studies of transient and periodically
variable sources. Consequently, LCO has had a strong impact on several fields
in astronomy. It has enabled follow-up observations for exoplanet transit
surveys such as the Kilodegree Extremely Little Telescope survey (e.g. Ahlers
et al., 2020; Stevens et al., 2020; Rodríguez Martínez et al., 2020), which
has discovered a number of planets around bright stars. LCO is leading the
detection and tracking of near Earth objects (NEOs) and asteroids (e.g. Lister
et al., 2017, 2021) and has also made the first optical follow-up observations
to a kilonova observed by the LIGO gravitational-wave observatory (Arcavi et
al., 2017a, b; McCully et al., 2017).

Since June 2019, the LCO network has been in use to acquire near-nightly
images of ZTF J013906.17+524536.89, a white dwarf which exhibits transits
likely caused by planetary debris (Vanderbosch et al., 2020). Each LCO image
covers a field-of-view of $26^{\prime}\times 26^{\prime}$ and observes a
region close to the Galactic plane ($b=-9.4^{\circ}$) abundant in stars with
an average stellar density of $\approx 8$ arcmin-2. Thus, besides their
primary purpose, these observations provide a unique opportunity to leverage
the high stellar density of the field to detect and classify variable sources
in a relatively unexplored patch of sky from the time-domain perspective.
Crucially, this observational dataset enables us to demonstrate the potential
for discovery in long-term imaging observations as the astronomy community
prepares for the Legacy Survey of Space and Time (LSST) to be undertaken by
the Vera C. Rubin Observatory (Ivezić et al., 2019).

Figure 2: Plot of a full frame and zoomed in LCO image with sources (black)
overlaid with the locations of Pan-STARRS1 DR2 sources (red circles). The
cross-matched sources remaining after application of source cuts, on which
aperture photometry is performed, are marked with blue crosses. The x and y
axes are in units of pixels. The variable sources identified in this work are
marked with yellow circles and the exclusion regions around bright stars are
shown as white boxes.

In this paper, we present the results of our search for periodic variable
stars and eclipsing systems in LCO photometric data. We organize this paper
into the following sections. Section 2 provides an overview of our
observations, source selection and aperture photometry methods, and period
analysis techniques that enable us to uncover variable stars in our dataset.
Section 3 details the preliminary classifications of the variable stars based
on considerations of their periods, light curve morphology, and position on
the Gaia EDR3 color-magnitude diagram. Section 4 discusses and compares our
results with two prominent surveys that observed in our field-of-view.
Finally, we conclude in section 5.

##  2 Methods

###  2.1 Observations

We utilize observations from the LCO 1.0-m telescope network obtained as part
of an extended monitoring program for the white dwarf ZTF J013906.17+524536.89
(hereafter ZTF J0139+5245, Vanderbosch et al., 2020) which, due to its
northern location and relatively faint $g=18.5$ magnitude, could only be
observed by the two 1.0-m telescope nodes located at McDonald Observatory
(telescope codes ELP 06 and 08). Images were acquired near-nightly using the
Sinistro imaging instrument (Brown et al., 2013) whenever ZTF J0139+5245 was
seasonably available. For each visit, a sequence of 3$-$6 consecutive images
was acquired in both the LCO $gp$ and $rp$ bands, with occasional images in
the $ip$ band. Exposure times were typically between 2 and 3 minutes to
achieve a signal-to-noise ratio $>30$ for the intended target, ZTF J0139+5245.
We have a total of 1560 images spanning 537 days, with a median separation
between visits of 0.98 days. Completed observations were bias, dark, and flat-
field corrected via the LCO BANZAI pipeline.

###  2.2 Source Selection and Aperture Photometry

To identify good candidates for aperture photometry within our images, we
first queried the second data release of the Pan-STARRS1 (PS1) survey
(Chambers et al., 2016) for sources within the field of view covered by an
average LCO exposure. Since our images are all centered roughly on the
coordinates of ZTF J0139+5245 ($\alpha=24.77633493\,$deg,
$\delta=+52.7602692\,$deg) with a field of view of $26^{\prime}\times
26^{\prime}$, we performed a cone search within the PS1 catalogue centered on
these coordinates with a search radius of $20^{\prime}$ to ensure full
coverage of each LCO exposure. This query resulted in a total of 5181 sources
in our image frame’s field-of-view.

First, we cut this PS1 query down by excluding extremely bright stars and the
objects surrounding them in boxes of width 100–220 pixels (see Figure 2) to
avoid spurious detections of variability arising from scattering, diffraction,
or bleeding of light from bright stars into nearby pixels in the CCD. Second,
we applied magnitude limits of 12.0–19.25 to the PS1 $r$-band to exclude
saturated bright objects and nearly undetectable dim objects to ensure
accurate photometry. We used the $r$-band for magnitude cuts since the
majority of objects in our field-of-view are red and have much fainter
$g$-band magnitudes, typically between 20 and 20.5 mag for the faintest
objects. This way we ensure that we can use both the $g$ and $r$-band LCO
photometry for all of our objects. Third, we applied a source separation limit
of $12$ arcseconds to exclude crowded objects that can cause errors in
brightness measurements and centroid finding while performing aperture
photometry.

Figure 3: A $G$ vs $BP-RP$ color-magnitude diagram (CMD) depicting locations
of variable sources found, for which Gaia Early Data Release 3 (EDR3) data are
available. The background CMD is also constructed using the Gaia EDR3 dataset.
Rough regions representing the location of common variable stars and systems
in the CMD are drawn as ellipses for reference (Gaia Collaboration et al.,
2019). Marker colors do not correspond to the color scheme adopted for the
representation of different variability types.

We reduced our sample size to 2227 sources in the photometry target list after
the application of the above cuts. A total of 2954 PS1 sources were excluded
from consideration of which $>50\%$ were removed on application of the
magnitude cut. We cross-matched this sample with each individual image, each
time excluding any objects which came within 50 pixels of the CCD edge.
Sources were located by converting each source RA and Dec coordinate to pixel
coordinates using the LCO image header WCS information, followed by a 2D
Gaussian fit to identify the centroid. A 1D Gaussian was then fit to the
radial profile of each PSF to determine the median FWHM of the image. Using
the Astropy-affiliated photutils package (Astropy Collaboration et al., 2013;
Caswell et al., 2019), we then performed circular aperture photometry on each
source with an aperture radius equal to the median FWHM. Local sky subtraction
was performed using a circular annulus centered on each source with inner and
outer radii of 32 and 48 pixels respectively. The median sky counts per pixel
within each annulus were subtracted from the respective apertures.

To calibrate our LCO photometry onto the PS1 magnitude scale, we followed the
methods described in Vanderbosch et al. (2020). The difference between PS1 and
instrumental magnitudes ($m_{d}=m_{\mathrm{PS1}}-m_{\mathrm{LCO}}$) was
measured while filtering for outliers and likely galaxy candidates. We then
solved for a zero-point offset ($z$) and color term ($c$) with a least squares
fit to $m_{d}=z+c\,(g_{\mathrm{PS1}}-r_{\mathrm{PS1}})$.

We also obtained additional time-series photometric data for each source (if
available) by querying ZTF DR4 to improve the accuracy of our period analysis.
Following recommendations in the ZTF Science Data System Explanatory
Supplement111http://web.ipac.caltech.edu/staff/fmasci/ztf/ztf_pipelines_deliverables.pdf,
for objects with ZTF photometry we removed poor quality detections by
requiring $\textit{catflags}=0$ and $|\textit{sharp}|<0.5$.

###  2.3 Period Analysis of Candidate Variables

To assess variability in our sample of objects covered by the LCO images, we
first generated Lomb-Scargle periodograms (Lomb, 1976; Scargle, 1982). We
applied barycentric time corrections to all images from both LCO and ZTF using
the Astropy Python package (Astropy Collaboration et al., 2018, 2013), to
ensure that any stable periodicities present in our objects are accurately
recovered. For each object, we used the frequency of the highest periodogram
peak to also generate a phase-folded light curve. We then visually inspected
both the light curves and periodograms for each object to identify promising
candidates for variability. Given our small sample size ($\approx\,$2000), we
visually selected objects based on qualitative features in order to maximize
the true positive fraction of candidate variables in our FOV — objects
possessing well-defined, continuously variable light curves as well as objects
with a clear set of peaks in their periodograms indicating an underlying
period. We refer the reader to Debosscher et al. (2007); Kim et al. (2011,
2014); Woźniak et al. (2004); Masci et al. (2014); Armstrong et al. (2016);
Holl et al. (2018); Jayasinghe et al. (2019, 2020); Heinze et al. (2018) for
more detailed discussions on statistical selection techniques, more suitable
for survey-size samples of sources larger than what is considered in our
study.

For 216 objects in our sample, we observed significant periodogram peaks
recurring at near 1-day periodicities. These are usually considered to be
artifacts of the window function and hence were excluded from our list of
candidates. We find two possible exceptions to the above (sources 23 and 28)
where the light curves indicated potential astrophysical variability. However,
we note that there is an ambiguity associated with these detections due to the
difficulty in distinguishing real variable sources from artifacts at near
1-day periods. In general, we identified candidate variables as objects with
significant periodogram peaks which did not occur at aliases of 1-day. For
these sources, we then identified the underlying period(s) using Lomb-Scargle
periodograms generated with the python package Pyriod (Bell, 2020). Finally,
with the refined periods from Pyriod, we re-generated the phase-folded light
curves (Figures 4, 5) for each candidate variable star.

We make a couple of cautionary comments. Using the above period determination
technique, the most significant peaks in the periodograms of eclipsing
binaries will sometimes correspond to half the orbital period of the system if
ellipsoidal variations or two equally-sized primary and secondary transits are
present. Accurate period determination for binaries may thus involve setting
the folding frequency to half the peak frequency obtained from the
periodograms. For sources whose variability may be caused by pulsations,
rotation, chromospheric activity, etc., it is unclear whether the folding
frequency should be modified, as above, to determine the period. This
highlights the fundamental limitation in determining the period of objects
without a classification. For such objects in our dataset, wherever necessary,
we present the light curve generated with no changes to the folding frequency,
i.e., by setting the folding frequency equal to the peak frequency in the
periodogram. We also apply a $5\sigma$ limit to the data points in each light
curve using the sigma_clip function from astropy.stats to exclude outliers
that could increase the noise levels and potentially mask true peaks in the
Lomb-Scargle periodogram. We then note that the periodograms in Figures 4, 5
are generated after the application of sigma clipping but the corresponding
light curves still present the clipped data points for reference.

###  2.4 Gaia Early Data Release 3 (EDR3) Color-Magnitude Diagram

We obtain photometry for our candidate variable sources from the publicly
accessible Gaia EDR3 dataset (Gaia Collaboration et al., 2021). We note that
the large fractional parallax uncertainties precluded a simple inversion of
the parallax measurement provided by Gaia to calculate the distance to each of
the variable sources. Instead, we determined the absolute $G$ band mean
magnitude using the geometric distance estimates provided by Bailer-Jones et
al. (2021). Errors on the relevant quantities quoted in Table 1 are then
derived using standard propagation of errors.

Figure 4: Light curves and periodograms for BY Draconis Variables identified
in our set of variable sources. The left panel in each sub-figure represents
the phase-folded light-curve data in $gp$, $rp$, and $ip$ LCO bands as well as
$g$ and $r$ ZTF bands, if available, with respective error measurements. The
right panel in each sub-figure represents the periodogram results with a red
diamond marking the peak frequency and a blue triangle marking the folding
frequency used. If the folding frequency is equal to the peak frequency, the
two markers are presented one on top of another. If the folding frequency is
half the peak frequency, the markers are presented in the same line, at their
respective frequencies. A spectral window is displayed as an inset plot
calculated based on the time-sampling of our observations. A common legend for
the light curves is presented. In each light curve, the survey data and the
final phase folded period (in days) used is noted.

##  3 Results

We present the 38 variable sources found by visual inspection of light curves
and periodograms obtained from aperture photometry of sources cross-matched
with the PS1 Survey, after application of source data cuts, in Table 1 along
with their preliminary classifications. The phase-folded light curves and
periodogram results are shown in Figures 4 and 5. The primary source of
variability observed is due to rotation and eclipses. We describe certain
general characteristics of the sources found, grouped by their preliminary
classification, below.

Figure 5: The same as Figure 4 but for Eclipsing Binaries.

###  3.1 Variability due to eclipses and rotation

  1. 1.

BY Draconis Variables: These are main-sequence stars with late spectral types
(K and M) that exhibit quasi-periodic light curves due to spots and
chromospheric activity. The fluctuations are low amplitude, generally less
than 0.5 mag. Their periods range from a fraction of a day to nearly 120 days.
We classify sources (1)–(28) by this definition. The periods of our candidates
are in the range 0.1–20d. This classification is supported by the location of
these sources on the Gaia EDR3 color-magnitude diagram (G3CMD) along the main
sequence. Since nearly every K and M star has a light curve characteristic of
a BY Draconis, this category of variables is expected to have the largest
number of detections in our sample, as we do indeed find. The light curves and
corresponding periodograms of all BY Draconis variables are presented in
Figure 4. We note that sources (23) and (28) have a period very close to a 1
day period characteristic of artifacts in the window function. Thus, it is
possible that while the light curves of these sources indicate potential
variability, these might not be true BY Draconis variables. This highlights
our limitation of confirming true variables with periods near the 1 day
period.

  2. 2.

Eclipsing Binaries of type EA: These are binaries with spherical or slightly
ellipsoidal components with well-separated, nearly constant light curves in
between minima. Secondary minima can be absent in these types. The prototype
for this class is Algol. We classify sources (29)–(31) by this definition.
This classification is supported by the locations of these sources on the
G3CMD along the main sequence. The periods of our candidates are in the range
0.4–0.8d. The light curves and corresponding periodograms of all eclipsing
binaries are presented in Figure 5.

  3. 3.

Eclipsing Binaries of type EW: These are binaries where the components are
nearly or actually in contact and both minima are virtually equally strong.
Onsets and ends of minima are not well defined in the light curves of these
objects. Their periods are generally less than one day. The prototype for this
class is W Ursae Majoris. We classify sources (32)–(38) using this definition.
This classification is supported by the location of these sources on the Gaia
EDR3 color-magnitude diagram (G3CMD) along the main sequence, similar to type
EA eclipsing binaries. The periods of our candidates are in the range
0.25–98d.

###  3.2 Variability due to pulsations

We note that the light curves of BY Draconis variables, as presented in Figure
4, are often similar to first overtone Cepheids, and in the case of low period
variables, similar to $\delta$ Scuti stars, leading to an ambiguity in
assigning a preliminary classification. Gaia Collaboration et al. (2019) find
constrained regions in a $G$ vs $BP-RP$ Gaia color-magnitude diagram for each
type of variable star. We can then eliminate ambiguity in our classifications
using the location of our objects in the CMD presented in Figure 3. We observe
that none of our sources lie in the regions dominated by classical cepheids
(Ceph), type II cephieds (CephII), or long period variables (LPV). There are a
few sources that fall within the $\delta$ Scuti ($\delta$Sct) region, but
these either possess transits characteristic of eclipsing binary systems or
periods outside the range of periods expected for $\delta$ Scuti variables.
Moreover, since K and M dwarfs are magnetically active, a certain amount of
scatter and irregularity in the shape of the light curve is expected compared
to the smoother light curves of Cepheids and $\delta$ Scutis. This allows us
to further strengthen our preliminary classification in certain cases. Lastly,
we note that source (16) is not presented in the Gaia CMD due to the lack of a
parallax measurement and hence holds a larger uncertainty with respect to its
preliminary classification.

Table 1: Variable Sources Found in LCO Images

Index | Pan-STARRS ID | R.A. (J2000) | Decl. (J2000) | Period | $\langle g\rangle$ | $\langle r\rangle$ | $\langle i\rangle$ | Classification | $G$ | $BP-RP$
---|---|---|---|---|---|---|---|---|---|---
|  | (deg) | (deg) | (days) | (mag) | (mag) | (mag) |  | (mag) | (mag)
(1) | 171470250272081779 | 25.02717 | 52.89257 | 2.288 | 19.15 | 18.35 | 17.96 | BYDrac | $7.83_{-0.37}^{+0.66}$ | $1.50\pm 0.04$
(2) | 171180244883024042 | 24.48828 | 52.65282 | 0.583 | 19.51 | 18.33 | 17.65 | BYDrac | $7.23_{-0.39}^{+0.49}$ | $2.05\pm 0.04$
(3) | 171340249027671181 | 24.90277 | 52.78374 | 11.526 | 18.82 | 18.07 | 17.74 | BYDrac | $6.19_{-0.48}^{+0.55}$ | $1.30\pm 0.02$
(4) | 171110250391044651 | 25.03911 | 52.59500 | 11.765 | 19.10 | 17.82 | 17.21 | BYDrac | $8.04_{-0.21}^{+0.20}$ | $1.86\pm 0.03$
(5) | 171140247106853549 | 24.71066 | 52.61904 | 0.853 | 17.08 | 16.11 | 15.64 | BYDrac | $6.38_{-0.10}^{+0.09}$ | $1.59\pm 0.01$
(6) | 171120249053056178 | 24.90530 | 52.60457 | 1.903 | 19.45 | 18.04 | 17.05 | BYDrac | $8.68_{-0.16}^{+0.17}$ | $2.46\pm 0.04$
(7) | 171490245651510736 | 24.56514 | 52.90840 | 0.694 | 17.29 | 16.09 | 15.45 | BYDrac | $7.99_{-0.04}^{+0.03}$ | $1.90\pm 0.01$
(8) | 171230246907620573 | 24.69076 | 52.69156 | 10.581 | 17.24 | 16.27 | 15.86 | BYDrac | $6.98_{-0.07}^{+0.09}$ | $1.53\pm 0.01$
(9) | 171360247742360866 | 24.77424 | 52.80015 | 12.007 | 16.53 | 15.78 | 15.46 | BYDrac | $6.14_{-0.09}^{+0.08}$ | $1.27\pm 0.01$
(10) | 171290245026655286 | 24.50266 | 52.74551 | 0.133 | 15.39 | 14.92 | 14.73 | BYDrac | $2.85_{-0.23}^{+0.24}$ | $0.97\pm 0.01$
(11) | 171200248226272288 | 24.82264 | 52.66799 | 6.909 | 15.67 | 14.90 | 14.56 | BYDrac | $5.68_{-0.09}^{+0.09}$ | $1.33\pm 0.01$
(12) | 171480246728186487 | 24.67282 | 52.90483 | 0.162 | 19.32 | 18.74 | 18.44 | BYDrac | $5.77_{-0.89}^{+0.97}$ | $1.17\pm 0.06$
(13) | 171490245557812506 | 24.55579 | 52.90988 | 1.307 | 19.39 | 18.53 | 18.14 | BYDrac | $6.69_{-0.56}^{+0.82}$ | $1.52\pm 0.04$
(14) | 171310250983433394 | 25.09835 | 52.76059 | 6.133 | 15.31 | 14.61 | 14.30 | BYDrac | $2.12_{-0.77}^{+1.06}$ | $1.23\pm 0.01$
(15) | 171320246965463740 | 24.69654 | 52.76920 | 16.500 | 17.48 | 16.39 | 15.90 | BYDrac | $7.58_{-0.06}^{+0.08}$ | $1.67\pm 0.01$
(16) | 171240249993444675aaGaia EDR3 does not quote a parallax measurement. | 24.99945 | 52.70332 | 19.021 | 16.93 | 16.13 | 15.79 | BYDrac | $\cdots$ | $\cdots$
(17) | 171360244167352479 | 24.41674 | 52.80152 | 6.774 | 18.59 | 17.45 | 16.93 | BYDrac | $7.58_{-0.18}^{+0.19}$ | $1.71\pm 0.01$
(18) | 171470245496025935 | 24.54960 | 52.89607 | 4.291 | 20.25 | 18.87 | 17.57 | BYDrac | $8.63_{-0.24}^{+0.22}$ | $2.68\pm 0.04$
(19) | 171190247395279443 | 24.73952 | 52.66563 | 2.097 | 18.93 | 18.13 | 17.73 | BYDrac | $6.54_{-0.44}^{+0.71}$ | $1.41\pm 0.03$
(20) | 171160250745668048 | 25.07457 | 52.63950 | 19.201 | 17.56 | 16.22 | 15.35 | BYDrac | $8.56_{-0.03}^{+0.03}$ | $2.17\pm 0.01$
(21) | 171380246037492681 | 24.60373 | 52.81834 | 14.151 | 17.02 | 16.22 | 15.88 | BYDrac | $5.99_{-0.12}^{+0.16}$ | $1.35\pm 0.01$
(22) | 171430244750309507 | 24.47503 | 52.86571 | 9.377 | 18.81 | 17.46 | 16.72 | BYDrac | $8.27_{-0.11}^{+0.14}$ | $2.08\pm 0.02$
(23) | 171050249931526875bbPossible artifacts of the window function with a 1 day period. | 24.99314 | 52.54684 | 0.993 | 15.75 | 15.10 | 14.82 | BYDrac | $4.64_{-0.08}^{+0.08}$ | $1.17\pm 0.01$
(24) | 171470246482432949 | 24.64823 | 52.89356 | 20.345 | 16.67 | 15.78 | 15.41 | BYDrac | $6.55_{-0.06}^{+0.05}$ | $1.43\pm 0.01$
(25) | 171540244285547663 | 24.42854 | 52.95582 | 0.635 | 17.08 | 16.19 | 15.76 | BYDrac | $5.23_{-0.14}^{+0.17}$ | $1.50\pm 0.01$
(26) | 171200245752030093 | 24.57520 | 52.66616 | 6.751 | 16.15 | 15.37 | 15.04 | BYDrac | $6.21_{-0.04}^{+0.04}$ | $1.33\pm 0.01$
(27) | 171380245939520760 | 24.59393 | 52.81674 | 4.541 | 18.43 | 17.65 | 17.30 | BYDrac | $6.04_{-0.32}^{+0.40}$ | $1.33\pm 0.02$
(28) | 171360251049928675bbPossible artifacts of the window function with a 1 day period. | 25.10500 | 52.80664 | 1.088 | 16.33 | 15.97 | 15.84 | BYDrac | $3.22_{-0.30}^{+0.49}$ | $0.81\pm 0.01$
(29) | 171320246599753528 | 24.65998 | 52.76903 | 0.747 | 17.46 | 16.73 | 16.34 | EA | $5.41_{-0.23}^{+0.28}$ | $1.34\pm 0.02$
(30) | 171370244878453108 | 24.48785 | 52.81034 | 0.804 | 19.41 | 18.86 | 18.58 | EA | $5.59_{-0.55}^{+0.76}$ | $1.05\pm 0.08$
(31) | 171560244345567723 | 24.43454 | 52.97254 | 0.422 | 17.06 | 16.54 | 16.32 | EA | $4.36_{-0.34}^{+0.41}$ | $1.07\pm 0.04$
(32) | 171060248174965790 | 24.81747 | 52.55424 | 0.257 | 17.72 | 16.82 | 16.41 | EW | $5.34_{-0.26}^{+0.47}$ | $1.45\pm 0.02$
(33) | 171480244845204434 | 24.48451 | 52.90315 | 0.383 | 17.74 | 16.98 | 16.62 | EW | $5.95_{-0.32}^{+0.42}$ | $1.35\pm 0.04$
(34) | 171240250341981190 | 25.03427 | 52.70046 | 0.665 | 16.55 | 16.10 | 15.99 | EW | $3.49_{-0.56}^{+0.69}$ | $0.90\pm 0.03$
(35) | 171480250400887540 | 25.04008 | 52.90570 | 0.295 | 19.43 | 18.76 | 18.46 | EW | $6.42_{-0.61}^{+0.93}$ | $1.22\pm 0.05$
(36) | 171510247091430730 | 24.70913 | 52.92505 | 0.628 | 15.07 | 14.80 | 14.71 | EW | $2.96_{-0.16}^{+0.13}$ | $0.71\pm 0.01$
(37) | 171070245379218746 | 24.53793 | 52.56505 | 97.387 | 17.65 | 16.71 | 16.25 | EW | $2.89_{-0.52}^{+0.96}$ | $1.60\pm 0.02$
(38) | 171440251042963672 | 25.10429 | 52.86914 | 3.298 | 14.90 | 14.43 | 14.23 | EW | $3.60_{-0.08}^{+0.08}$ | $0.95\pm 0.01$

Note. — EA and EW represent sub-types of Eclipsing Binaries and BYDrac
represents a BY Draconis Variable.

Figure 6: A comparison of the periodograms generated for each of our variable
stars using three sets of survey data — ZTF DR2 (red), ZTF DR5 (yellow), and
LCO (black). The diamond markers represent the peak frequency obtained with
the periodogram generated for each dataset. The first column of the panel
comprising sources (15), (17), (29) represents examples of sources for which
the peak frequency obtained agrees across all three datasets. The second
column of the panel comprising sources (2), (6), (8) represents examples of
sources for which the peak frequency obtained with LCO agrees only with that
obtained with ZTF DR5. The third column of the panel comprising sources (5),
(10), (19) represents examples of sources for which the peak frequency
obtained with LCO does not agree with ZTF DR2 and DR5.

##  4 Discussion

###  4.1 ZTF Catalog of Variable Stars

We performed a box search of size 25 arcminutes on our FOV centered on the
object J013906.17+524536.8 in the catalog of ZTF variable stars provided by
Chen et al. (2020, hereafter C20). We then cross-matched our list of candidate
variable stars with the results of the box search. In our image’s FOV, C20
have identified 7 variable sources, i.e., approximately a factor of five less
than our identification of 38 variable sources. Out of these seven sources,
our cross-match yielded four matching identifications: (5), (30), (32), (37).
This suggests that both sets of data are missing variable stars in this FOV.
Further investigation revealed that our crowding-limitation selection
criterion eliminated 3/7 sources in C20’s coverage of our FOV.

C20 classify (5) as a BY Draconis variable with period 0.8556985d, (30) as an
EA-type eclipsing binary system with period 0.8038496d, (32) as an EW-type
eclipsing binary system with period 0.257176d, and (37) as a Semiregular
variable with period 50.2286076d. Our classifications and derived periods for
these objects are in agreement with the above as presented in Table 1, with
the exception of (37) which we identify as an EW type eclipsing binary system
with a folded period of 97.387d.

We generated periodograms for all our candidate variable objects with three
sets of data — ZTF DR2, ZTF DR5, and LCO — to investigate the lack of variable
objects identified with ZTF. We find that there is either insufficient or no
ZTF DR2 data available for 19 sources. C20 utilize ZTF DR2 and hence would not
have found these sources. DR5 improves the availability of data over DR2, with
only 6 sources having either insufficient or no data. We present a subset of
the generated periodograms in Figure 6 under three categories — example
sources where the peak frequency is in agreement across all three datasets,
example sources where the peak frequency for LCO agrees with that of ZTF DR5
but not with ZTF DR2, and example sources where the peak frequency across all
three datasets are in disagreement. We note that with the exception of (21),
the peak frequency obtained with DR2 and DR5 simultaneously agree only when it
also agrees with the peak frequency obtained with LCO data. Agreement of our
peaks with DR5 but not DR2 as well as the observed lower noise level with DR5
data compared with DR2 data, in all cases presented, highlights an improvement
in the data quality in the newer ZTF data release.

We highlight two possible reasons for the higher number of detections obtained
with our LCO observations compared to the ZTF Catalog of Variable Stars
compiled by C20 using ZTF DR2 observations as well as compared to ZTF DR5.

First, we quantify and compare the photometric precision per exposure for LCO
and ZTF photometry using the Median Absolute Deviation (MAD) indicator. Figure
7 compares the MAD of $g$ and $r$-band light curves obtained with LCO with
those obtained with ZTF DR2 and DR5. We exclude all variable sources
identified in this work along with the white dwarf with transiting planetary
debris — ZTF J0139+5245. We observe that our LCO photometry has a lower MAD on
average at essentially all magnitudes, and performs significantly better in
the $g$-band at faint magnitudes. In the $r$-band, the LCO photometric
precision is comparable to ZTF for most of our objects, though we do not
sample the same range of magnitudes as in the $g$-band due to our imposed
magnitude limits (see Section 2.2). Aside from the main groups of objects in
Figure 7, we also find several outliers with much higher MAD values in LCO
than ZTF. We investigated the light curves and LCO images of these objects and
found that they all fall into one of two categories which make their LCO
photometry unreliable: (1) sources which are frequently found near to or on
the edge of the CCD, and (2) faint sources contaminated by nearby bright
stars. These outliers represent a small subset of our total sample, but
indicate that the source selection process we used in Section 2.2 is not
perfect at removing potentially problematic sources. While several factors
contribute to the photometric precision for a given set of observations, such
as atmospheric transmission, CCD quantum efficiency, telescope aperture size,
and many others, we suspect that a major factor improving the photometric
precision of our LCO observations compared to ZTF observations is the
difference in exposure times. Our LCO exposure times range from 2$-$3 minutes
while the ZTF exposure times are always 30 seconds.

Figure 7: Comparison of ZTF and LCO photometry. Each panel compares the median
absolute deviations (MAD) of the $g$-band (top panels) and $r$-band (bottom
panels) light curves of LCO with those of ZTF DR2 (left panels) and DR5 (right
panels). A diagonal line in each panel represents a 1:1 relationship. The data
points are color coded by their median ZTF DR5 magnitudes. All variable
sources identified in this work and the white dwarf with transiting planetary
debris – ZTF J0139+5245 – are excluded from these figures. In general our LCO
photometry exhibits improved photometric precision per-exposure compared to
ZTF, especially for faint sources in the $g$-band, likely due to our longer
exposure times of 2$-$3 minutes compared to 30 seconds for ZTF. Since we used
the PS1 $r$-band to perform a magnitude cut at 19.25 mag, the $r$-band
magnitudes for our sample do not go nearly as faint as the $g$-band since our
sources are mostly red. The outliers in these panels that have much higher MAD
in LCO than ZTF are primarily sources found near to CCD edges and have
unreliable photometry, or are faint sources contaminated by nearby bright
stars.

Second, we find that the photometric time-sampling of our observations with
LCO is greater than both ZTF DR5 and DR2 by about a factor of 2.5 — we have an
average of 1141 photometric datapoints across a time baseline of 537 days with
LCO ($\approx$ 2.12 per day), 553 photometric datapoints across a time
baseline of 719 days with ZTF DR5 ($\approx$ 0.77 per day), and 293
photometric datapoints across a time baseline of 368 days with ZTF DR2
($\approx$ 0.79 per day).

The result of a longer exposure time and denser photometric time sampling with
the LCO dataset is lower noise levels in the periodograms compared to ZTF. We
quantify the noise levels using the statistical mean of the periodogram data.
We find that, on average, the ZTF noise level is $\approx$ 3.6 times greater
than the LCO noise level. Consequently, this implies that a 3$\sigma$ peak
detection in the LCO periodograms would correspond to a 0.8$\sigma$ non-
detection in the ZTF periodograms, assuming similar peak amplitudes.
Regardless of the exact reason, however, the demonstrated improvement in
photometric precision and noise levels enables us to detect lower-amplitude
variability in the objects covered by our LCO images.

###  4.2 ATLAS Catalog of Variable Stars

We performed a box search of size 25 arcminutes on our FOV centered on the
object J013906.17+524536.8 in the ATLAS catalog of variable stars provided by
Heinze et al. (2018, hereafter H18). We then cross-matched our list of
candidate variable stars with the results of the box search. In our image’s
FOV, H18 have identified 84 variable sources. We note that they classify 77/84
of the detections as dubious sources, defined as stars that might not be real
variables. Our candidate variable star list confirms 6/7 of the remaining
sources that were assigned a classifications other than dubious. Our
investigation revealed that our crowding limitation selection criterion
eliminated one of these 7 sources from consideration. Moreover, we have
confirmed variability in two additional matches that were classified as
dubious by H18. Out of the remaining 75 dubious sources, 44 sources made it
through our selection criterion for visual inspection. We found no variability
in these objects. The remaining 31 dubious sources have been excluded from
consideration by our selection criterion and thus their variability can
neither be confirmed nor denied.

The 6/7 sources in H18 with classifications other than dubious matched with
(31), (32), (33), (34), (36), and (38) in our candidate variable star list. We
confirmed variability in source (7), classified as a BY Draconis variable, and
(29), classified as an EA type eclipsing binary system, which were classified
as dubious sources by H18.

H18 classify (32), (34), and (36) as SINE variables with periods of 0.257171d,
0.332667d, and 0.627621d respectively. SINE variables are defined as
sinusoidal variables, likely dominated by ellipsoidal variables, in which
stars exhibit simple sine-wave variability with little residual noise. Our
classifications and period identifications for the above sources are in
agreement, with the exception of source (34) for which we identified secondary
eclipses and hence a period of 0.665d. While H18 broadly classify the above
sources as SINE variables based on features identified using machine-learning
algorithms, a technique very well-suited to large sample sizes, given our
small sample size, we are able to classify these objects more specifically as
Eclipsing Binaries of type EW.

H18 classify (31) as a CBH with a period of 0.422072d. CBH stands for close
binary, half period. These stars are contact or near-contact eclipsing
binaries for which the Fourier fit has settled on half the correct period and
hence has overlapped the primary and secondary eclipses. They also classify
(33) as an IRR with a period of 0.382738d, where IRR stands for “irregular”
variables. This class serves as a “catch-all” bin for objects that do not seem
to fit into any of their more specific categories. Most of the stars
classified as IRR do not show coherent variations that can be folded cleanly
with a single period.

Our identified periods for the above sources are in agreement but our
classification for (33) is in disagreement. Based on the light curve presented
in Figure 5, we classify (33) as an Eclipsing Binary of type EW.

##  5 Conclusion

In this paper, we have presented the results of aperture photometry and
Fourier analysis performed on 2227 sources within 1560 LCO images spanning 537
days and cross-matched with the PS1 catalogue. Our main findings are as
follows:

1\. We have identified 38 variable sources and have classified them on the
basis of characteristic light curve properties, periods, as well as their
locations within the Gaia CMD: 28 candidate BY Draconis variables, 3 candidate
eclipsing binaries of type EA, and 7 candidate eclipsing binaries of type EW.
The complete results are presented in Table 1 and Figures 4 and 5. In
assigning preliminary classifications to the detected sources, we have
demonstrated the applicability of the Gaia CMD as a powerful classification
tool for future surveys to utilize.

2\. We determined that 27 of the 38 detected variable sources are newly
discovered or newly classified (71%) based on comparisons with previously
published catalogs, thereby increasing the number of detections in the field-
of-view under consideration by a factor of $\approx$ 2.5. In addition, we
confirmed variability in two additional sources previously observed but
classified as dubious sources.

3\. Using the median absolute standard deviations (MAD) of the $g$-band and
$r$-band light curves obtained with LCO and ZTF photometric data, we
demonstrate that, in general, our LCO photometry exhibits improved photometric
precision per-exposure compared to ZTF, especially for faint sources in
$g$-band. We attribute this improvement to LCO’s longer exposure times of
2$-$3 minutes compared to the 30 seconds of ZTF. The greater time-sampling of
LCO photometry combined with the longer exposure time reduces the noise levels
in our periodograms and improves our sensitivity to potential periodicities by
a factor of $\approx$ 3.6. This enables us to increase the total number of
detections in this field-of-view.

4\. Finally, since the majority of sources we detect do not have spectral
classifications, we note that it may be beneficial to conduct follow-up
spectroscopic observations to further solidify the variable star
classifications we provide. Further, additional time-series photometry for
objects with the least significant periodogram peaks may enable better
characterization of these sources.

## Acknowledgements

This work was conducted as part of the White Dwarf Stars Research Stream
(Montgomery et al., 2010; Winget & Montgomery, 2013) in the Freshman Research
Initiative (FRI) program (Clark et al., 2016; Rodenbusch et al., 2017) in the
College of Natural Sciences at the University of Texas at Austin. A.S.
acknowledges financial support from the Freshman Research Initiative Summer
Fellowship awarded by the College of Natural Sciences at the University of
Texas at Austin.

Z.P.V. and M.H.M. acknowledge support from the United States Department of
Energy under grant DE-SC0010623, the National Science Foundation under grant
AST-1707419, and the Wootton Center for Astrophysical Plasma Properties under
the United States Department of Energy collaborative agreement DE-NA0003843.
M.H.M. acknowledges support from the NASA ADAP program under grant
80NSSC20K0455.

This work is based on observations obtained with the Samuel Oschin 48-inch
Telescope at the Palomar Observatory as part of the Zwicky Transient Facility
project, which is supported by NSF grant No. AST-1440341 and the participating
institutions of the ZTF collaboration. We make use of observations and
reduction software from the LCO network, and include data from the Pan-STARRS1
Surveys (PS1) and the PS1 public science archive, which have been made
possible through contributions by participating institutions of the PS1
collaboration (https://panstarrs.stsci.edu/). We also make use of data from
the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia)
processed by the Gaia Data Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium).

_Additional Software/Resources_ : Astropy (Astropy Collaboration et al., 2013,
2018), Photoutils (Bradley et al., 2019), Pyriod maintained by Dr. Keaton Bell
(Bell, 2020), the NASA Astrophysics Data System (ADS), and SIMBAD and VizieR
(operated at CDS, Strasbourg, France).

## References

  * Ahlers et al. (2020) Ahlers, J. P., Johnson, M. C., Stassun, K. G., et al. 2020, AJ, 160, 4, doi: 10.3847/1538-3881/ab8fa3
  * Alcock et al. (1993) Alcock, C., Akerlof, C. W., Allsman, R. A., et al. 1993, Nature, 365, 621, doi: 10.1038/365621a0
  * Arcavi et al. (2017a) Arcavi, I., Hosseinzadeh, G., Howell, D. A., et al. 2017a, Nature, 551, 64, doi: 10.1038/nature24291
  * Arcavi et al. (2017b) Arcavi, I., McCully, C., Hosseinzadeh, G., et al. 2017b, ApJ, 848, L33, doi: 10.3847/2041-8213/aa910f
  * Armstrong et al. (2016) Armstrong, D. J., Kirk, J., Lam, K. W. F., et al. 2016, MNRAS, 456, 2260, doi: 10.1093/mnras/stv2836
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147, doi: 10.3847/1538-3881/abd806
  * Bell (2020) Bell, K. J. 2020, in American Astronomical Society Meeting Abstracts, Vol. 235, American Astronomical Society Meeting Abstracts #235, 106.06
  * Bradley et al. (2019) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2019, astropy/photutils: v0.7.2, v0.7.2, Zenodo, doi: 10.5281/zenodo.3568287
  * Brown et al. (2013) Brown, T. M., Baliber, N., Bianco, F. B., et al. 2013, PASP, 125, 1031, doi: 10.1086/673168
  * Caswell et al. (2019) Caswell, T. A., Droettboom, M., Hunter, J., et al. 2019, matplotlib/matplotlib v3.1.0, v3.1.0, Zenodo, doi: 10.5281/zenodo.2893252
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560.  https://arxiv.org/abs/1612.05560
  * Chen et al. (2019) Chen, X., Wang, S., Deng, L., et al. 2019, Nature Astronomy, 3, 320, doi: 10.1038/s41550-018-0686-7
  * Chen et al. (2020) —. 2020, ApJS, 249, 18, doi: 10.3847/1538-4365/ab9cae
  * Clark et al. (2016) Clark, G., Russell, J., Enyeart, P., et al. 2016, PLoS Biology, 14, 1, doi: 10.1371/journal.pbio.1002368
  * Debosscher et al. (2007) Debosscher, J., Sarro, L. M., Aerts, C., et al. 2007, A&A, 475, 1159, doi: 10.1051/0004-6361:20077638
  * Flewelling et al. (2020) Flewelling, H. A., Magnier, E. A., Chambers, K. C., et al. 2020, ApJS, 251, 7, doi: 10.3847/1538-4365/abb82d
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2019) Gaia Collaboration, Eyer, L., Rimoldini, L., et al. 2019, A&A, 623, A110, doi: 10.1051/0004-6361/201833304
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Heinze et al. (2018) Heinze, A. N., Tonry, J. L., Denneau, L., et al. 2018, AJ, 156, 241, doi: 10.3847/1538-3881/aae47f
  * Holl et al. (2018) Holl, B., Audard, M., Nienartowicz, K., et al. 2018, A&A, 618, A30, doi: 10.1051/0004-6361/201832892
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jacyszyn-Dobrzeniecka et al. (2016) Jacyszyn-Dobrzeniecka, A. M., Skowron, D. M., Mróz, P., et al. 2016, Acta Astron., 66, 149.  https://arxiv.org/abs/1602.09141
  * Jacyszyn-Dobrzeniecka et al. (2017) —. 2017, Acta Astron., 67, 1, doi: 10.32023/0001-5237/67.1.1
  * Jayasinghe et al. (2019) Jayasinghe, T., Stanek, K. Z., Kochanek, C. S., et al. 2019, MNRAS, 486, 1907, doi: 10.1093/mnras/stz844
  * Jayasinghe et al. (2020) —. 2020, MNRAS, 491, 13, doi: 10.1093/mnras/stz2711
  * Kim et al. (2014) Kim, D.-W., Protopapas, P., Bailer-Jones, C. A. L., et al. 2014, A&A, 566, A43, doi: 10.1051/0004-6361/201323252
  * Kim et al. (2011) Kim, D.-W., Protopapas, P., Byun, Y.-I., et al. 2011, ApJ, 735, 68, doi: 10.1088/0004-637X/735/2/68
  * Lanzafame et al. (2018) Lanzafame, A. C., Distefano, E., Messina, S., et al. 2018, A&A, 616, A16, doi: 10.1051/0004-6361/201833334
  * Larson et al. (2003) Larson, S., Beshore, E., Hill, R., et al. 2003, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 35, AAS/Division for Planetary Sciences Meeting Abstracts #35, 36.04
  * Leavitt & Pickering (1912) Leavitt, H. S., & Pickering, E. C. 1912, Harvard College Observatory Circular, 173, 1
  * Lister et al. (2017) Lister, T., Greenstreet, S., Gomez, E., Christensen, E. J., & Larson, S. M. 2017, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 49, AAS/Division for Planetary Sciences Meeting Abstracts #49, 110.10
  * Lister et al. (2021) Lister, T. A., Gomez, E., Chatelain, J., et al. 2021, arXiv e-prints, arXiv:2102.10144.  https://arxiv.org/abs/2102.10144
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447, doi: 10.1007/BF00648343
  * Magnier et al. (2020a) Magnier, E. A., Chambers, K. C., Flewelling, H. A., et al. 2020a, ApJS, 251, 3, doi: 10.3847/1538-4365/abb829
  * Magnier et al. (2020b) Magnier, E. A., Schlafly, E. F., Finkbeiner, D. P., et al. 2020b, ApJS, 251, 6, doi: 10.3847/1538-4365/abb82a
  * Magnier et al. (2020c) Magnier, E. A., Sweeney, W. E., Chambers, K. C., et al. 2020c, ApJS, 251, 5, doi: 10.3847/1538-4365/abb82c
  * Masci et al. (2014) Masci, F. J., Hoffman, D. I., Grillmair, C. J., & Cutri, R. M. 2014, AJ, 148, 21, doi: 10.1088/0004-6256/148/1/21
  * Masci et al. (2018) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2018, Publications of the Astronomical Society of the Pacific, 131, 018003, doi: 10.1088/1538-3873/aae8ac
  * McCully et al. (2017) McCully, C., Hiramatsu, D., Howell, D. A., et al. 2017, ApJ, 848, L32, doi: 10.3847/2041-8213/aa9111
  * Montgomery et al. (2010) Montgomery, M. H., Winget, D. E., Allen, A., et al. 2010, in Astronomical Society of the Pacific Conference Series, Vol. 432, New Horizons in Astronomy: Frank N. Bash Symposium 2009, ed. L. M. Stanford, J. D. Green, L. Hao, & Y. Mao, 225
  * Pietrzyński et al. (2013) Pietrzyński, G., Graczyk, D., Gieren, W., et al. 2013, Nature, 495, 76, doi: 10.1038/nature11878
  * Pietrzyński et al. (2019) Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Nature, 567, 200, doi: 10.1038/s41586-019-0999-4
  * Riess et al. (2018) Riess, A. G., Casertano, S., Yuan, W., et al. 2018, ApJ, 861, 126, doi: 10.3847/1538-4357/aac82e
  * Rodenbusch et al. (2017) Rodenbusch, S. E., Hernandez, P. R., Simmons, S. L., & Dolan, E. L. 2017, CBE Life Sci Educ., 15, 1, doi: 10.1187/cbe.16-03-0117
  * Rodríguez Martínez et al. (2020) Rodríguez Martínez, R., Gaudi, B. S., Rodriguez, J. E., et al. 2020, AJ, 160, 111, doi: 10.3847/1538-3881/ab9f2d
  * Saunders et al. (2014) Saunders, E. S., Lampoudi, S., Lister, T. A., Norbury, M., & Walker, Z. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9149, Observatory Operations: Strategies, Processes, and Systems V, ed. A. B. Peck, C. R. Benn, & R. L. Seaman, 91490E
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835, doi: 10.1086/160554
  * Shappee et al. (2014) Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48, doi: 10.1088/0004-637X/788/1/48
  * Skowron et al. (2019) Skowron, D. M., Skowron, J., Mróz, P., et al. 2019, Science, 365, 478, doi: 10.1126/science.aau3181
  * Stevens et al. (2020) Stevens, D. J., Zhou, G., Johnson, M. C., et al. 2020, MNRAS, 499, 3775, doi: 10.1093/mnras/staa3142
  * Torres et al. (2010) Torres, G., Andersen, J., & Giménez, A. 2010, A&A Rev., 18, 67, doi: 10.1007/s00159-009-0025-1
  * Udalski et al. (1994) Udalski, A., Szymanski, M., Kaluzny, J., et al. 1994, ApJ, 426, L69, doi: 10.1086/187342
  * Vanderbosch et al. (2020) Vanderbosch, Z., Hermes, J. J., Dennihy, E., et al. 2020, ApJ, 897, 171, doi: 10.3847/1538-4357/ab9649
  * Waters et al. (2020) Waters, C. Z., Magnier, E. A., Price, P. A., et al. 2020, ApJS, 251, 4, doi: 10.3847/1538-4365/abb82b
  * Winget & Montgomery (2013) Winget, D. E., & Montgomery, M. H. 2013, in American Astronomical Society Meeting Abstracts, Vol. 221, American Astronomical Society Meeting Abstracts, #115.01
  * Woźniak et al. (2004) Woźniak, P. R., Williams, S. J., Vestrand, W. T., & Gupta, V. 2004, AJ, 128, 2965, doi: 10.1086/425526


# Feature Space Hijacking Attacks against Differentially Private Split
Learning

Grzegorz Gawron
LiveRamp, Privacy Tech R&D
<EMAIL_ADDRESS>Philip Stubbings
LiveRamp, Privacy Tech R&D
<EMAIL_ADDRESS>
###### Abstract

Split learning and differential privacy are technologies with growing
potential to help with privacy-compliant advanced analytics on distributed
datasets. Attacks against split learning are an important evaluation tool and
have been receiving increased research attention recently. This work’s
contribution is applying a recent feature space hijacking attack (FSHA) to the
learning process of a split neural network enhanced with differential privacy
(DP), using a client-side off-the-shelf DP optimizer. The FSHA attack obtains
client’s private data reconstruction with low error rates at arbitrarily set
DP $\epsilon$ levels. We also experiment with dimensionality reduction as a
potential attack risk mitigation and show that it might help to some extent.
We discuss the reasons why differential privacy is not an effective protection
in this setting and mention potential other risk mitigation methods.

##  1 Introduction

Privacy compliant data processing has been receiving an increasing amount of
attention recently. Take a consortium of hospitals as a motivational example.
Each owns a data set describing treatments given to patients, together with
the treatments’ outcomes. There is a potential great value in learning from
all of the data sets taken together. However, the hospitals are bound by
privacy regulations that prevent them from sharing the data in raw form. For a
thorough discussion of distributed learning from patient data see Vepakomma et
al. [2018a].

Split learning (see Gupta and Raskar [2018]), a type of federated learning,
has been proposed to address scenarios like this. It is a neural network
training mechanism, which splits the network into modules corresponding to the
private data silos, so that each module is trained within the silo and no raw
data is exposed outside. Another reason for growing popularity of split
learning is it’s performance. Since only the processed or smashed data of
lower dimensions is being exchanged among silos and their coordinator, the
process is more efficient than some of the more established mechanisms of
distributed learning, such as federated averaging. Yet, relying on smashing
the data as the only mechanism of protecting privacy has been shown by Erdogan
et al. [2021b] to be misguided and that even an honest-but-curious actor
having access to smashed data can reconstruct the private inputs. Vepakomma et
al. [2018b] have surveyed various methods to deal with split learning privacy
leaks having, what they call a no peek learning, as their target.

Figure 1: Plots comparing the private input (top row) and its corresponding
reconstruction (bottom row) in a successful FSHA attack on MNIST.  Figure 2:
FSHA split network architecture with the DP logic implemented by a DPOptimizer
from Tensorflow-Privacy package

In this paper we focus on another technique which can be used to enhance split
learning’s privacy compliance: differential privacy (DP). Its promise is to
provide a rigorous mathematical privacy guarantee (see Dwork et al. [2014]).
Abadi et al. [2016] showed that training with DP is feasible, while having
control over model performance degradation and privacy cost. There has also
been work by Papernot et al. [2020] on fine tuning deep learning with DP to
minimise the performance cost with some good results. Geyer et al. [2017] is
another example of the work on optimizing the DP deep learning process.

Theoretical privacy guarantees, when applied in practice, need to be backed by
automatic evaluation, which helps to guard against plain bugs or other
implementation issues. There are multiple, already well established libraries
available, which provide the necessary evaluation components, e.g.
Adversarial-Robustness-Toolbox111https://github.com/Trusted-AI/adversarial-
robustness-toolbox or Tensorflow-
Privacy222https://github.com/tensorflow/privacy. For a comprehensive overview
of attack and defences literature and evaluation best practices we refer the
reader to Carlini et al. [2019].

Attack based evaluation is particularly relevant for newly emerging methods
like split learning. There have already been first papers published on this
subject. One recent attack method is feature space hijacking attack FSHA by
Pasquini et al. [2021], where a malicious attacker taking over the server
module is able to reconstruct all of the client private input data with
arbitrarily low reconstruction error on MNIST and other visual datasets. See
Figure 1 for an illustration of a successful FSHA attack on a MNIST dataset,
based on an experiment reproduction from the original paper - without DP
applied.

As a testament to the rapid pace of research progress in the field, methods
guarding against such attacks have already been presented, see Erdogan et al.
[2021a].

Our contributions. This work contributes by applying differential privacy to
the client module of a split neural network, running a suite of FSHA attacks
against it and drawing conclusions for differentially private split learning
system implementations. To our knowledge this is the first FSHA style attack
attempt at such a setup. Hitaj et al. [2017] have shown the limits of applying
DP to general collaborative deep learning and we validate the vulnerability of
split learning in that respect.

We reuse the software prototype from Pasquini et al. [2021] found at
https://github.com/pasquini-dario/SplitNN_FSHA and aim to introduce minimal
code changes to reach our goal. We also follow the authors’ choice of using
Tensorflow and use its sub-package Tensorflow-Privacy for the application of
DP clipping and noise. See McMahan et al. [2019] for the description of this
package.

##  2 Architecture

In what follows we briefly describe the original FSHA attack split learning
architecture and refer the reader to the original paper of Pasquini et al.
[2021] for a complete discussion of the components (we use the same notation
as in the original paper). We restrict our investigation to a situation where
the attacker takes over only the server. The clients are not taken over and
believe they are part of a legitimate cooperation. However, as highlighted in
the original paper, the same attack would work successfully by taking over one
of the clients only. The example neural network in Figure 2 is split into two
main components: client and server (in general there could be multiple
clients). The client owns their private data stored in a silo $X_{priv}$ and
is responsible for training the $f$ model component of the network. The server
runs the attacker’s code, which instead of running its part of the split
learning protocol, trains two components: discriminator $D$ and an auto-
encoder $\tilde{f}$ and $\tilde{f}^{-1}$. The auto-encoder, when given inputs
from $X_{pub}$, encodes them using $\tilde{f}$ to an internal representation
(feature space) and decodes using $\tilde{f}^{-1}$ with the loss set, so as to
learn to reconstruct the original input data. Finally, discriminator $D$ is
responsible for classifying the input as coming from either $f$ model or the
$\tilde{f}$ model. The ultimate attacker’s target is for the discriminator $D$
to force the private model $f$’s output to come from the same feature space as
the outputs of $\tilde{f}$ and the $\tilde{f}^{-1}$ to be able to decode it
with minimal error, just as it learned to decode the $\tilde{f}$.

The standard way of applying differential privacy during neural network
training is on the gradients during back propagation. In our scenario the most
natural place to apply the DP clipping and noise to the $f$ gradients is just
before applying the gradient updates sent back by server’s discriminator $D$
as illustrated by the green fragments in Figure 2. This allows the client data
owner to control the noise application process and configure it according to
their privacy preferences.

##  3 Experiments

###  3.1 Experimental Setup

We enhanced the FSHA code with DP by introducing only minimal changes, with
the aim of being able to easily compare the results. We used TensorFlow-
privacy python package to add a differentially private optimizer
DPAdamGaussianOptimizer. We parametrised DP with $\delta=0$ and various levels
of $\epsilon$ from $1$ to $100$.

We deployed several GCP Notebooks instances with the following hardware: (a)
a2-highgpu-1g (Accelerator Optimized: 1 NVIDIA Tesla A100 GPU, 12 vCPUs, 85GB
RAM); and (b) n1-standard-4 (4 vCPUs, 15 GB RAM, 1 NVIDIA Tesla T4). GPU load
was $30$ percent on (a) and $80$ percent on (b) with similar run times, so we
settled on the latter to optimize the experiment cost.

(a) 10 000 iterations, DP with $\epsilon=10$.

(b) 10 000 iterations, DP with $\epsilon=100$.

Figure 3: Private input (top row) and its corresponding reconstruction (bottom
row) for various epsilon settings and 10 000 iterations.

###  3.2 MNIST Results

We start with the MNIST dataset because it is widely used for benchmarking and
its simplicity helps to clearly show the impact of the application of DP.
Note, that we focus on showing the reconstruction visualisations, as they
nicely communicate the results. The diagrams of the actual reconstruction
errors can be found in the Appendix.

(a) 100 000 iterations, DP with $\epsilon=0.5$.

(b) 100 000 iterations, DP with $\epsilon=10$.

Figure 4: Plots comparing the private input (top row) and its corresponding
reconstruction (bottom row) for various epsilon settings and 100 000
iterations.

We reproduced the results from the original paper, without any differential
privacy applied. Figure 1 shows a successful FSHA reconstruction attack after
10 000 iterations of $64$ element mini-batches. This corresponds to $10$
epochs of training, since there are $N=\numprint{60000}$ samples in the
training set. The next step is to run the same number of iterations using DP
with $\epsilon=10$ and $\delta=1/N$. Figure 3(a) shows that the reconstruction
was far less successful, however the reconstruction error was still declining
after the training stopped. This suggested that more iterations might make the
attack successful. As expected the attack against a split learning with
$\epsilon=100$ is able to reconstruct the private data with much less error
(Figure 3(b)). We continued with the training process running for 100 000
iterations ($100$ epochs) and at $\epsilon=10$ it was enough to produce very
low-error reconstructions, as illustrated in Figure 4(b). Finally, we set
$\epsilon$ to a much lower value $0.5$ and obtained the reconstruction results
shown in Figure 4(a). Although the error is much higher than with the previous
experiment (where $\epsilon=10$) the reconstruction error was still steadily
declining, so by increasing the number of iterations we could expect the error
to go down as per previous experiments.

See the Appendix for a more complete record of reconstruction plots.

###  3.3 Fashion-MNIST Results

Beyond the Fashion-MNIST attacks reproduced from Pasquini et al. [2021] we
have experimented with removing some of the image classes from the public data
set and applying differential privacy in both settings.

(a) no DP, iterations $=\numprint{100000}$

(b) $\epsilon=10$, iterations $=\numprint{100000}$

(c) no DP, iterations $=\numprint{100000}$

(d) $\epsilon=10$, iterations $=\numprint{100000}$

Figure 5: FSHA against Fashion-MNIST with category $0$ (tops) missing from the
public data set in Figures 5(a) and 5(b) and category $8$ (bags) missing in
Figures 5(c) and 5(d).

We first removed the category $0$ representing t-shirts/tops and ran the
attacks with and without the DP. We present the relevant results in Figure
5(a) and Figure 5(b). In both cases the shapes were mostly visible, but the DP
version had less logo details visible (see the two pictures to the right).

There are, however, other categories in the public data set (see Figure 10)
that are similar to the t-shirts, so the process had more chance to learn the
relevant features to also reconstruct t-shirts. In an attempt to minimise this
effect, we remove the category which is least similar to others: bags
(category = $8$). The result is illustrated in Figures 5(c) and 5(d).
Interestingly, the version without DP reconstructs some of the bags as shoes,
perhaps because they activate similar latent features learned by the attack
model. However we can recognise certain features coming from the private data
e.g., the circle on the right t-shirt. With DP added we see no ’bag as shoe’
effect. The images are more blurred, but the general shapes are reconstructed.
This might be due to regularisation effect of DP. It is important to note that
it should be possible to improve this results by running more training
iterations (we stopped at 100 000).

To summarise, Fashion-MNIST data set with some of the classes removed from the
public data set achieves lower reconstruction precision, especially with DP
applied. This helps demonstrate the point that the reconstructions are really
generated out of the feature space driven by the access to similar public
data. The closer the public and private data distributions are, the more
successful the attacks.

###  3.4 Mitigating the attack with dimensionality reduction

We explored the application of DP during the model training process as a means
to mitigate the attack and demonstrated that DP can at most delay FSHA
convergence.

In an un-trusted environment in which it can be assumed a client or server
will attempt to exploit the learning protocol by means of the FSHA method,
additional steps should be taken to protect the underlying data. One
possibility is to apply dimensionality reduction techniques directly to the
data prior to training to a level which still yields acceptable model
accuracy.

(a) Image compression

(b) FSHA recovered images

Figure 6: PCA image compression and recovery. $1^{st}$ row - no compression,
$2^{nd}$ row - principal components = 4, $3^{rd}$ row - principal components =
2.

To test this hypothesis, we re-run the FSHA attack up to 100k iterations on
image data which has been compressed using various levels of PCA principal
components as shown in Figure 6(a) (also in Appendix Figure 11(a)). We make
use of the CelebA dataset and define the classification task to predict
whether the celebrity is likely male or female.

The resulting recovered images shown in Figure 6(b) (also in Appendix Figure
11(b)) demonstrate that applying dimensionality reduction to the data can
mitigate the attack to a degree in which the celebrities are no longer
identifiable. One should note though, that some reconstructed image features
might still be obvious and unique to certain individuals.

##  4 Conclusions

The FSHA mechanism used is very robust. It only requires that the client part
of the model can be successfully trained; and for any useful usage that must
be the case. Once the attacker obtains a sufficiently good client model they
can reconstruct the private data which are passed from the client to server
during the training process. As Pasquini et al. [2021] notice, they could also
infer certain features out of the private data, if that’s what they choose to
do.

It might seem surprising that the DP doesn’t give enough protection against
FSHA. However, the DP noise applied during the training is designed to prevent
the model from memorizing data so that the private data cannot be retrieved
from the trained model itself. It must still allow for training useful models.
After a successful training, a good DP model must be able to achieve high
inference accuracy. The inferred values might well be sensitive and constitute
a privacy leak. Hence, having access to the model being trained allows an
attacker to either reconstruct or infer properties of the private data. We
refer the reader to Hitaj et al. [2017] for a more thorough discussion of this
topic.

Employing some dimensionality reduction techniques might make the potential
reconstructions useless, but we must remember that it might not always produce
a model with good enough utility and even more importantly the attacker might
still infer sensitive properties from the private data.

Detection methods such as ones proposed by Erdogan et al. [2021a] might also
be used, but they are not bullet proof and require constant attention to new
detection method workarounds used by potential attackers.

As for now, to ensure that the split learning process is safe we should
mitigate the risk of someone taking over the server (or indeed any other split
network component). The risk mitigation might rely on enforcing a secure
compute environment with e.g. only a certified version of the training code
running on all the modules of the split network. It also follows that the
application of the DP trained split neural network model should happen in an
environment of trust. That is, the owners of private data should trust that
inferred outputs are appropriately taken care of. It is safe, for instance, if
the model is applied directly by the data owner on their own data and they own
the inferred results.

## Acknowledgements

This work has benefited from valuable input during the 2021 OpenDP Fellows
Program (see https://opendp.org/fellows-program).

## References

  * Abadi et al. [2016] M. Abadi, A. Chu, I. Goodfellow, H. B. McMahan, I. Mironov, K. Talwar, and L. Zhang.  Deep Learning with Differential Privacy.  _Proceedings of the ACM Conference on Computer and Communications Security_ , 24-28-Octo:308–318, jul 2016.  doi:10.1145/2976749.2978318.  URL http://arxiv.org/abs/1607.00133http://dx.doi.org/10.1145/2976749.2978318.
  * Carlini et al. [2019] N. Carlini, A. Athalye, N. Papernot, W. Brendel, J. Rauber, D. Tsipras, I. Goodfellow, A. Madry, and A. Kurakin.  On evaluating adversarial robustness, 2019.
  * Dwork et al. [2014] C. Dwork, A. Roth, C. Dwork, and A. Roth.  The Algorithmic Foundations of Differential Privacy.  _Foundations and Trends in Theoretical Computer Science_ , 9:211–407, 2014.  doi:10.1561/0400000042.
  * Erdogan et al. [2021a] E. Erdogan, A. Kupcu, and A. E. Cicek.  Splitguard: Detecting and mitigating training-hijacking attacks in split learning.  _IACR Cryptol. ePrint Arch._ , 2021:1080, 2021a.
  * Erdogan et al. [2021b] E. Erdogan, A. Kupcu, and A. E. Cicek.  Unsplit: Data-oblivious model inversion, model stealing, and label inference attacks against split learning.  _IACR Cryptol. ePrint Arch._ , 2021:1074, 2021b.
  * Geyer et al. [2017] R. C. Geyer, T. Klein, and M. Nabi.  Differentially private federated learning: A client level perspective, 2017.
  * Gupta and Raskar [2018] O. Gupta and R. Raskar.  Distributed learning of deep neural network over multiple agents, 2018\.
  * Hitaj et al. [2017] B. Hitaj, G. Ateniese, and F. Perez-Cruz.  Deep Models Under the GAN: Information Leakage from Collaborative Deep Learning.  _Proceedings of the ACM Conference on Computer and Communications Security_ , pages 603–618, feb 2017.  URL https://arxiv.org/abs/1702.07464v3.
  * McMahan et al. [2019] H. B. McMahan, G. Andrew, U. Erlingsson, S. Chien, I. Mironov, N. Papernot, and P. Kairouz.  A general approach to adding differential privacy to iterative training procedures, 2019.
  * Papernot et al. [2020] N. Papernot, S. Chien, S. Song, A. Thakurta, and U. Erlingsson.  Making the shoe fit: Architectures, initializations, and tuning for learning with privacy.  https://openreview.net/forum?id=rJg851rYwH, 2020.  Accessed: 2021-11-22.
  * Pasquini et al. [2021] D. Pasquini, G. Ateniese, and M. Bernaschi.  Unleashing the tiger: Inference attacks on split learning, 2021.
  * Vepakomma et al. [2018a] P. Vepakomma, O. Gupta, T. Swedish, and R. Raskar.  Split learning for health: Distributed deep learning without sharing raw patient data, 2018a.
  * Vepakomma et al. [2018b] P. Vepakomma, T. Swedish, R. Raskar, O. Gupta, and A. Dubey.  No peek: A survey of private distributed deep learning, 2018b.

##  Appendix A Appendix

Below we include the complete figures from our experiments, where they didn’t
fit in the main body of the paper.

(a) 10 000 iterations, no DP applied (original reconstruction).

(b) 100 000 iterations, DP with $\epsilon=0.5$.

(c) 10 000 iterations, DP with $\epsilon=100$.

(d) 100 000 iterations, DP with $\epsilon=10$.

(e) 10 000 iterations, DP with $\epsilon=10$.

Figure 7: Reconstruction errors (vertical axis) against iterations (horizontal
axis) for various iteration levels and DP application parameters.  Figure 8:
Reconstruction attempts for $\epsilon=10$ and 100 000 iterations. Odd rows
show the private data fed to the private local $f$ model and the even rows
below show the corresponding images reconstructed using attacker’s $f^{-1}$
model. The lower the line, the more training iterations there were. Figure 9:
Reconstruction attempts for $\epsilon=0.5$ and 100 000 iterations. Odd rows
(with little titles) show the private data fed to the private local $f$ model
and the even rows below show the corresponding images reconstructed using
attacker’s $f^{-1}$ model. The lower the line, the more training iterations
there were. Figure 10: A few examples per each of the ten Fashion-MNIST
categories

(a) Image compression

(b) FSHA recovered images

Figure 11: PCA image compression and recovery. $1^{st}$ row - no compression,
$2^{nd}$ row - principal components = 4, $3^{rd}$ row - principal components =
2.