# Recurring Planetary Debris Transits and Circumstellar Gas around White Dwarf
ZTF J0328–1219

Zachary P. Vanderbosch Department of Astronomy, University of Texas at Austin,
Austin, TX-78712, USA  McDonald Observatory, Fort Davis, TX-79734, USA  Saul
Rappaport Department of Physics and Kavli Institute for Astrophysics and Space
Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Joseph A. Guidry Department of Astronomy, University of Texas at Austin,
Austin, TX-78712, USA  Bruce L. Gary Hereford Arizona Observatory, Hereford,
AZ 85615, USA  Simon Blouin Los Alamos National Laboratory, PO Box 1663, Los
Alamos, NM 87545, USA  Thomas G. Kaye Foundation for Scientific Advancement,
Sierra Vista, AZ, USA  Alycia J. Weinberger Earth and Planets Laboratory,
Carnegie Institution for Science, 5241 Broad Branch Rd NW, Washington, DC
20015, USA  Carl Melis Center for Astrophysics and Space Sciences, University
of California, San Diego, CA 92093-0424, USA  Beth L. Klein Department of
Physics and Astronomy, University of California, Los Angeles, CA 90095-1562,
USA  B. Zuckerman Department of Physics and Astronomy, University of
California, Los Angeles, CA 90095-1562, USA  Andrew Vanderburg Department of
Astronomy, University of Wisconsin-Madison, Madison, WI 53706, USA  J. J.
Hermes Department of Astronomy, Boston University, 725 Commonwealth Ave.,
Boston, MA 02215, USA  Ryan J. Hegedus University of North Carolina at Chapel
Hill, Department of Physics and Astronomy, Chapel Hill, NC 27599, USA
Matthew. R. Burleigh School of Physics and Astronomy, University of Leicester,
Leicester, LE1 7RH, UK  Ramotholo Sefako South African Astronomical
Observatory, P.O Box 9, Observatory 7935, Cape Town, South Africa  Hannah L.
Worters South African Astronomical Observatory, P.O Box 9, Observatory 7935,
Cape Town, South Africa  Tyler M. Heintz Department of Astronomy, Boston
University, 725 Commonwealth Ave., Boston, MA 02215, USA

(Received March 20, 2021; Revised May 25, 2021; Accepted June 2, 2021)

###### Abstract

We present follow-up photometry and spectroscopy of ZTF J0328$-$1219
strengthening its status as a white dwarf exhibiting transiting planetary
debris. Using TESS and Zwicky Transient Facility photometry, along with
follow-up high speed photometry from various observatories, we find evidence
for two significant periods of variability at 9.937 and 11.2 hr. We interpret
these as most likely the orbital periods of different debris clumps. Changes
in the detailed dip structures within the light curves are observed on
nightly, weekly, and monthly timescales, reminiscent of the dynamic behavior
observed in the first white dwarf discovered to harbor a disintegrating
asteroid, WD 1145+017. We fit previously published spectroscopy along with
broadband photometry to obtain new atmospheric parameters for the white dwarf,
with $M_{\star}=0.731\pm 0.023\,M_{\odot}$, $T_{\mathrm{eff}}=7630\pm 140\,$K,
and [Ca/He] $=-9.55\pm 0.12$. With new high-resolution spectroscopy, we detect
prominent and narrow Na D absorption features likely of circumstellar origin,
with velocities $21.4\pm 1.0$ km s-1 blue-shifted relative to atmospheric
lines. We attribute the periodically modulated photometric signal to dusty
effluents from small orbiting bodies such as asteroids or comets, but are
unable to identify the most likely material that is being sublimated, or
otherwise ejected, as the environmental temperatures range from roughly 400 K
to 600 K.

White dwarf stars — Transits — Planetesimals — Circumstellar dust —
Circumstellar gas — Roche limit — Tidal disruption

††journal: ApJ††facilities: TESS, PO:1.2m (ZTF), Gaia, McD:Struve (ProEM), DCT
(DeVeny), Magellan:Clay (MIKE), SOAR (Goodman), SAAO, ADS, CDS. ††software:
Astropy (Astropy Collaboration et al., 2013, 2018), iraf (National Optical
Astronomy Observatories), LMFIT (Newville et al., 2014), phot2lc (Vanderbosch
et al. in prep.), Photutils (Bradley et al., 2020), Starlink (Currie et al.,
2014), CDS’s (Strasbourg, France) SIMBAD and VizieR online pages and tables,
and the NASA Astrophysics Data System (ADS) repositories.

##  1 Introduction

White dwarfs, the most common endpoint of stellar evolution, are likely to
host planetary systems which have survived post-main sequence evolution.
Indirect evidence for the presence of these planets includes atmospheric metal
pollution in 30% of white dwarfs cooler than 20,000 K (Zuckerman et al., 2003,
2010; Koester et al., 2014), infrared excess caused by a warm circumstellar
debris disc in more than 50 white dwarfs (Zuckerman & Becklin, 1987; Mullally
et al., 2007; Dennihy et al., 2020a; Xu et al., 2020), and in a small number
of these objects, metallic emission lines due to a hot gaseous component of
the debris disk (Gänsicke et al., 2006; Dennihy et al., 2020b; Melis et al.,
2020; Gentile Fusillo et al., 2020). These observational signatures are
generally agreed to be the result of the tidal disruption of a planetary body
after being gravitationally perturbed onto a highly eccentric orbit by one or
more of the surviving outer planets (Debes & Sigurdsson, 2002; Jura, 2003).

Another way to observe this planetary debris is via the transit method, when
dust-emitting debris clumps passing in front of the white dwarf momentarily
attenuate some of the stellar flux along our line of sight. To date, only two
objects have been confirmed to exhibit recurring transit events due to
planetary debris, WD 1145+017 with debris in 4.5–4.9 hr orbits (Vanderburg et
al., 2015), and ZTF J0139+5245 with debris in a ${\approx}100$ day orbit
(Vanderbosch et al., 2020). The more than 2-order of magnitude difference in
orbital periods may be indicative of different evolutionary stages during the
standard tidal disruption process, though other processes such as rotational
fission (Veras et al., 2020) may be at play in ZTF J0139+5245. Regardless of
the mechanism that is disrupting planetary bodies, it is likely that a broad
distribution of orbital periods could be observable among this class of
objects, and confirming more such objects is needed to empirically fill out
and understand this distribution.

Now in the age of large scale time-domain surveys, the likelihood of finding
new transiting planetary debris candidates has increased dramatically.
Publicly available light curves from both current and future surveys such as
the Zwicky Transient Facility (ZTF, Masci et al., 2019; Bellm et al., 2019),
the Catalina Real-Time Transient Survey (CRTS, Drake et al., 2009), the
Transiting Exoplanet Survey Satellite (TESS, Ricker et al., 2014), and the
Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST, Ivezić et
al., 2019), can be assessed for variability and cross-matched with catalogs of
known and candidate white dwarfs such as Gaia (e.g. Gentile Fusillo et al.,
2019) and the Montreal White Dwarf Database (MWDD, Dufour et al., 2017) to
identify good candidates for transiting debris systems.

In a recent study by Guidry et al. (2021), five new candidates for white
dwarfs exhibiting variability caused by transiting planetary debris were
identified using variability metrics based on Gaia DR2 and ZTF DR3 photometry.
Guidry et al. (2021) obtained follow up spectroscopy and time series
photometry to identify the spectral types of each object and confirm the
presence of variability. One of these five objects, ZTF J032833.52$-$121945.27
(hereafter ZTF J0328$-$1219), was identified as a metal-polluted, DZ white
dwarf111DZ follows the white dwarf classification system of Sion et al. (1983)
indicating a degenerate star whose dominant optical spectral features are from
elements heavier than He. with short-timescale, $\approx$10 % amplitude
variability in its 3-hr light curve from McDonald Observatory. Given the
promising evidence for transiting debris, we targeted this object for detailed
follow-up observations.

In this work, we present new and archival observations and identify periodic
optical variability and atomic absorption features likely of circumstellar
origin in ZTF J0328$-$1219, strengthening its classification as a white dwarf
with transiting planetary debris. In Section 2 we describe the new
observations obtained along with archival photometry, in Section 3 we discuss
the periodic variability, new white dwarf atmospheric parameters derived from
photometric and spectroscopic observations, and the detection of circumstellar
gas, in Section 4 we discuss the possible nature of ZTF J0328$-$1219 and
potential causes of the periodic variability, and in Section 5 we summarize
our results and draw some conclusions.

Table 1: Journal of Time-Series Photometric Observations of ZTF J0328$-$1219 Date | Start | Telescope | Filter | $t_{\mathrm{exp}}$ | Length |  | Date | Start | Telescope | Filter | $t_{\mathrm{exp}}$ | Length
---|---|---|---|---|---|---|---|---|---|---|---|---
| (UTC) |  |  | (s) | (hr) |  |  | (UTC) |  |  | (s) | (hr)
2020 Oct 19 | 06:30:09 | McD | BG40 | 10 | 3.19 |  | 2021 Jan 03 | 02:19:14 | JBO | clear | 60 | 4.88
2020 Dec 02 | 06:16:58 | HAO | clear | 90 | 4.29 |  | 2021 Jan 05 | 01:58:16 | JBO | clear | 60 | 5.82
2020 Dec 03 | 01:32:32 | HAO | clear | 90 | 9.01 |  | 2021 Jan 06 | 02:09:04 | JBO | clear | 60 | 6.23
2020 Dec 04 | 04:38:48 | HAO | clear | 90 | 5.88 |  | 2021 Jan 08 | 01:38:02 | JBO | clear | 60 | 5.90
2020 Dec 05 | 01:33:34 | HAO | clear | 90 | 8.92 |  | 2021 Jan 12 | 01:39:36 | HAO | clear | 90 | 5.80
2020 Dec 06 | 01:28:18 | HAO | clear | 90 | 8.89 |  | 2021 Jan 13 | 01:38:10 | HAO | clear | 90 | 5.69
2020 Dec 07 | 02:42:08 | HAO | clear | 90 | 7.35 |  | 2021 Jan 14 | 01:38:43 | HAO | clear | 90 | 5.69
2020 Dec 12 | 01:25:19 | HAO | clear | 90 | 3.49 |  | 2021 Jan 15 | 01:48:57 | HAO | clear | 90 | 5.60
| 03:22:49 | McD | BG40 | 10 | 5.59 |  | 2021 Jan 16 | 01:38:55 | HAO | clear | 90 | 5.75
2020 Dec 13 | 01:30:24 | HAO | clear | 90 | 8.40 |  | 2021 Jan 17 | 01:39:15 | HAO | clear | 90 | 4.91
| 02:45:03 | McD | BG40 | 10 | 5.90 |  | 2021 Jan 18 | 01:42:10 | HAO | clear | 90 | 5.47
| 19:08:02 | SAAO | V | 60 | 3.05 |  | 2021 Jan 27 | 01:47:34 | HAO | clear | 90 | 4.79
2020 Dec 14 | 02:36:19 | McD | BG40 | 20 | 6.56 |  | 2021 Jan 31 | 01:43:18 | HAO | clear | 90 | 4.56
2020 Dec 15 | 01:33:39 | HAO | clear | 90 | 8.29 |  | 2021 Feb 01 | 01:43:19 | HAO | clear | 90 | 2.93
| 02:32:49 | McD | BG40 | 20 | 6.63 |  | 2021 Feb 05 | 01:52:07 | HAO | clear | 90 | 4.12
2020 Dec 16 | 01:58:31 | HAO | clear | 90 | 7.87 |  | 2021 Feb 06 | 01:54:34 | HAO | clear | 90 | 3.88
2020 Dec 17 | 01:27:23 | HAO | clear | 90 | 8.43 |  | 2021 Feb 07 | 01:52:29 | HAO | clear | 90 | 4.04
| 03:53:41 | JBO | clear | 60 | 5.40 |  | 2021 Feb 09 | 01:31:56 | McD | BG40 | 5 | 3.76
2020 Dec 20 | 02:59:04 | JBO | clear | 60 | 6.24 |  |  | 01:51:44 | HAO | clear | 90 | 4.03
2020 Dec 21 | 04:50:08 | JBO | clear | 60 | 3.99 |  | 2021 Feb 10 | 01:26:58 | McD | BG40 | 5 | 3.82
2020 Dec 22 | 04:24:24 | JBO | clear | 60 | 4.02 |  | 2021 Feb 11 | 01:35:51 | McD | BG40 | 10 | 3.63
2020 Dec 31 | 04:40:12 | JBO | clear | 60 | 3.59 |  | 2021 Feb 12 | 01:28:04 | McD | BG40 | 10 | 3.80
2021 Jan 02 | 05:16:32 | JBO | clear | 60 | 3.02 |  |  | 01:55:02 | HAO | clear | 90 | 3.90

##  2 Observations

###  2.1 Time-Series Photometry

McDonald/ProEM: We acquired high-speed time-series photometry of ZTF
J0328$-$1219 on nine nights, 2020 October 19, 2020 December 12–15, and 2021
February 9–12, with the McDonald Observatory (McD) 2.1-m Otto Struve telescope
using the Princeton Instruments ProEM frame-transfer CCD. Located at
Cassegrain focus, the camera has a $1.6\times 1.6$ arcmin field of view. We
used $4\times 4$ binning resulting in 0.38″ per binned pixel. We exclusively
used the blue-bandpass BG40 filter, with exposure times of 5, 10 or 20-s.
Using standard calibration images acquired each night, we bias, dark, and flat
field corrected each science exposure with standard iraf routines. We
performed circular aperture photometry on ZTF J0328$-$1219 and the one
available comparison star in our field of view using the iraf package ccd_hsp
(Kanaan et al., 2002), with aperture radii ranging from 1.5 to 10 pixels in
half-pixel steps.

We extracted light curves generated from each aperture radius using the
phot2lc222phot2lc.readthedocs.io python package, and performed differential
photometry to remove atmospheric effects. With phot2lc we clipped $\geq
5\sigma$ outliers within a sliding window of width 30-min, and divided each
light curve by a low-order polynomial of degree two or less to remove airmass
trends. We selected the optimal aperture size based on the light curve with
the lowest average point-to-point scatter. Lastly, we used Astropy (Astropy
Collaboration et al., 2018, 2013) within phot2lc to apply barycentric
corrections to the mid-exposure time stamps of each image. A summary of
McDonald observations can be found in Table 1.

HAO: We acquired time-series photometry on 26 nights between 2020 December
2–17 (10 nights), 2021 January 12–18 (7 nights), and 2021 January 27 –
February 12 (8 nights) using the Astro-Tech 16-in Ritchey-Chretien telescope
located at Hereford Arizona Observatory (HAO). We used an SBIG ST-10XME CCD
camera with a KAF 3200E main chip to acquire a series of 90-s exposures with
9-s readout times each night using a clear filter. We performed bias, dark,
and flat-field calibrations, along with circular aperture photometry, using
the Maxim DL program. We tested aperture radii of 4, 5, and 6 pixels on the
target and comparison stars, choosing a weighted-average of the light curves
corresponding to apertures which maximized the signal-to-noise of each night’s
observations. Six comparison stars were used to correct for atmospheric
effects and calibrate the brightness of ZTF J0328-1219 using magnitudes and
B–V colors from the AAVSO Photometric All-Sky Survey (APASS) archive. We
applied extinction corrections to our target as a function of airmass, using
an empirically determined correction function based on our selection of
comparison stars. Lastly, we applied barycentric timestamp corrections using
Jason Eastman’s Barycentric Julian Date web
page333http://astroutils.astronomy.ohio-state.edu/time/utc2bjd.html. A summary
of HAO observations can be found in Table 1.

Table 2: Journal of Spectroscopic Observations of ZTF J0328$-$1219

Date | Facility/Instrument | $t_{\mathrm{exp}}$ | Seeing | Airmass | Slit Width | Wavelengths | $R$
---|---|---|---|---|---|---|---
|  | (s) | (″) |  | (″) | (Å) | ($\lambda/\Delta\lambda$)
2020 Nov 16 | LDT/DeVeny | 6$\times$180 | 1.7 | 1.50 | 1.0 | 3700$-$6800 | 1200
2020 Dec 07 | Magellan/MIKE | 2$\times$1800 | 0.5 | 1.06 | 1.0 | 3500$-$5060B† | $28{,}000$
|  |  |  |  |  | 5000$-$9400R† | $22{,}000$
2020 Dec 21 | SOAR/Goodman | 7$\times$600 | 1.2 | 1.16 | 1.0 | 7560$-$8750 | 3000

††footnotetext: B and R correspond to the blue and red arms of the MIKE
spectrograph.

JBO: We acquired time-series photometry on 10 nights between 2020 December 16
and 2021 January 7, using a 32-in Ritchey-Chretien telescope located at Junk
Bond Observatory (JBO). On each night we acquired a series of 60-s exposures
without any filter. We performed bias, dark, and flat-field corrections, along
with circular aperture photometry, using the AstroImageJ program (Collins et
al., 2017). We performed differential photometry using the same six comparison
stars used for the HAO observations, but no detrending procedures were
applied. We performed barycentric corrections to the mid-exposure timestamps
of each image using Astropy. A summary of JBO observations can be found in
Table 1.

SAAO: We acquired time-series photometry on one night, 2020 December 13, with
the South African Astronomical Observatory (SAAO)’s 1-m telescope called
Lesedi, and the Sutherland High-speed Optical Camera (SHOC; Coppejans et al.,
2013). The target was observed continuously for three hours with 60-s
exposures in V band. $2\times 2$ binning was used, giving $0.67$ arcsecs per
binned pixel. The field of view of the SHOC camera on Lesedi is $5.7\times
5.7$ arcmin squared, allowing several comparison stars of similar brightness
to be monitored for differential photometry. Conditions were very good, with
seeing $\approx 1.2$″. The images were bias and flat-field corrected using the
local SAAO SHOC pipeline, which is driven by python scripts running iraf tasks
(pyfits and pyraf). Aperture photometry was performed using the Starlink
package autophotom. We used a 3 pixel radius aperture that maximized the
signal to noise, and the background was measured in an annulus surrounding
this aperture. Differential photometry was then performed using three
comparison stars. The resultant light curve was detrended for air mass effects
with a second-order polynomial, and a barycentric correction was applied to
the time stamps using Astropy routines. A summary of SAAO observations can be
found in Table 1.

TESS: ZTF J0328$-$1219 (TIC 93031595) was observed by TESS in Sector 4 (2018
October 18 – November 14) and in Sector 31 (2020 October 22 – November 18).
During Sector 4, ZTF J0328$-$1219 was not specifically targeted for
observations, so the pixels surrounding ZTF J0328$-$1219 were only downloaded
after being co-added to an exposure time of 30 minutes. By the time target
lists were assembled for Sector 31, however, ZTF J0328$-$1219 had been
identified as a white dwarf candidate by Gaia DR2 Gentile Fusillo et al.
(2019), so it was targeted with TESS’s two-minute cadence mode. The Sector 4
and 31 light curves have durations of 25.8 and 25.4 days, and data gaps
lasting 6.2 and 2.2 days, respectively.

We extracted the Sector 4 TESS light curves with a custom aperture photometry
pipeline. We downloaded the pixels surrounding ZTF J0328$-$1219 using the
TESSCut interface at the Mikulski Archive for Space Telescopes and extracted a
light curve from a single-pixel stationary aperture (the closest pixel to the
R.A. and Decl. of ZTF J0328$-$1219). We corrected for spacecraft systematics
and scattered light by decorrelating the light curve against the mean and
standard deviation of the spacecraft quaternion time series within each
exposure (see Vanderburg et al., 2019) and the time series of background flux
outside the photometric aperture. Since ZTF J0328$-$1219 was observed with
two-minute cadence observations in Sector 31, its light curve was extracted by
the official Science Processing Operation Center (SPOC) pipeline (Jenkins,
2015; Jenkins et al., 2016) based at NASA Ames Research Center. The SPOC
pipeline works by performing simple (stationary) aperture photometry and
removing common-mode systematic errors with its Pre-search Data Conditioning
module (Smith et al., 2012; Stumpe et al., 2014).

ZTF: We queried the public ZTF survey (Masci et al., 2019; Bellm et al., 2019)
for $g$ and $r$-band DR5 photometry using the NASA/IPAC Infared Science
Archive (IRSA). Following recommendations in the ZTF Science Data System
Explanatory
Supplement444http://web.ipac.caltech.edu/staff/fmasci/ztf/ztf_pipelines_deliverables.pdf,
we removed poor quality detections by requiring $\textit{catflags}=0$ and
$|\textit{sharp}|<0.5$. The resulting ZTF light curve has 348 data points
($g=166$, $r=182$), all with exposure times of 30-s and spanning 898 days
between 2018 August 4 and 2021 January 18. Following light curve filtering, we
applied barycentric corrections to the time stamps using Astropy.

###  2.2 Spectroscopy

Table 3: MIKE Comparison Star Observations

Target | Source ID | Type | $G$ | Distance | Sky Sep.† | Date Obs. | $t_{\mathrm{exp}}$ | Airmass
---|---|---|---|---|---|---|---|---
| (Gaia eDR3) |  | (mag) | (pc) | (′) |  | (s) |
HD 22243 | 5163917352681876736 | A0 | 6.23 | 118.3 | 172.5 | 2020 Dec 08 | 36  | 1.06 |  |  |
HD 21917 | 5161837970035375488 | G6V | 8.75 | 73.5 | 45.0 | 2020 Dec 08 | 96  | 1.05 |  |  |
HD 21634 | 5161919059017811072 | F3V | 8.71 | 132.6 | 17.3 | 2020 Dec 08 | 120  | 1.05 |  |  |
Gaia J0328$-$1216 | 5161809451452454784 | — | 12.15 | 122.2 | 3.7 | 2020 Dec 08 | 210  | 1.05 |  |  |

Note. — Distances are the photogeometric values from Bailer-Jones et al.
(2021), and Sky Sep. is the on-sky separation between each comparison star and
ZTF J0328$-$1219.

We detail here all of the spectroscopic observations obtained and used in this
work. See Tables 2 and 3 for a summary of the instrument and observational
details.

SOAR/Goodman: We carried out spectroscopic observations of ZTF J0328$-$1219 on
2020 December 21 using the Goodman Spectrograph (Clemens et al., 2004) mounted
on the 4.1-m Southern Astrophysical Research (SOAR) telescope. We used the Red
Camera, the 1200l-R grating, a 1.0″ long slit, and the GG-495 order-blocking
filter to obtain a resolving power of $R{\approx}3000$ over the
$7560-8750\,\mathrm{\AA}$ wavelength region. We obtained seven consecutive
600-s exposures, for a total exposure time of 1.17 hr. We bias-subtracted,
flat-fielded, wavelength-calibrated, and optimally extracted (Marsh, 1989)
each individual spectrum using a set of custom Python tools based on the
zzceti_pipeline555https://github.com/joshfuchs/ZZCeti_pipeline reduction
package described in Fuchs (2017), before combining them into a final
spectrum. We performed wavelength calibrations using the abundant night sky
emission lines in this wavelength region, as identified by Osterbrock et al.
(1996). Lastly, we applied a heliocentric velocity correction to the combined
spectrum.

LDT/DeVeny: To obtain improved atmospheric parameters for ZTF J0328$-$1219, we
used the same spectrum presented in Guidry et al. (2021), which was obtained
using the DeVeny Spectrograph (Bida et al., 2014) mounted on the 4.3-m Lowell
Discovery Telescope (LDT). On 2020 November 16 with $1.7\arcsec$ seeing, we
acquired $6\times 180$ s exposures using the 300 line mm-1 DV2 grating and a
1.0″ slit to achieve a resolving power of $R{\approx}1200$ over 3700$-$6800 Å.
We debiased and flat-fielded the exposures using standard starlink routines
(Currie et al., 2014), and optimally extracted the spectrum (Horne, 1986)
using the software pamela. Using molly (Marsh, 1989) we applied heliocentric
velocity corrections to the wavelengths.

Magellan/MIKE: We observed ZTF J0328$-$1219 on 2020 December 7 using the
Magellan Inamori Kyocera Echelle (MIKE) spectrograph (Bernstein et al., 2003)
on the 6.5-m Magellan Clay Telescope at Las Campanas Observatory. Observations
were done at an airmass of 1.06. Two consecutive integrations of 1800-s each
were acquired with the 1′′ wide slit for a resolving power of about
$R{\approx}28{,}000$ in the blue channel (3500$-$5060Å) and
$R{\approx}22{,}000$ in the red channel (5000$-$9400Å). The signal-to-noise
(S/N) achieved was about 12 at 3960Å, 27 at 5890Å, and 30 at 6563Å. Wavelength
calibration was done with ThAr lamps taken just before and after the science
exposures. The data were extracted, flat-fielded, and wavelength calibrated
using the Carnegie Python MIKE pipeline which incorporates methods described
in Kelson et al. (2000) and Kelson (2003). We applied Heliocentric velocity
corrections calculated by the MIKE pipeline to each spectrum. On 2020 December
8, the same instrument setup was used to observe four stars near the line of
sight to ZTF J0328$-$1219 to probe the interstellar medium (ISM). Details of
the comparison stars and their observations can be found in Table 3.

Table 4: Astrometric and Photometric Properties of ZTF J0328$-$1219

Parameter | Value | Uncertainty | Ref.
---|---|---|---
Gaia eDR3 | 5161807767825277184 | – | 1
R.A. (J2016) | $03^{\mathrm{h}}28^{\mathrm{m}}33^{\mathrm{s}}.63$ | – | 1
decl. (J2016) | ${-}12^{\circ}19^{\prime}45^{\prime\prime}.49$ | – | 1
$\mu_{\mathrm{RA}}$ (mas yr-1) | 110.739 | 0.055 | 1
$\mu_{\mathrm{dec}}$ (mas yr-1) |  ${-}$14.39 | 0.046 | 1
parallax (mas) | 23.038 | 0.056 | 1
$d$ (pc) | 43.32 | 0.12 | 2
$G$ | 16.6186 | 0.0035 | 1
$G_{\rm{BP}}$ | 16.732 | 0.015 | 1
$G_{\rm{RP}}$ | 16.446 | 0.011 | 1
galex nuv | 18.91 | 0.05 | 3
$g$ | 16.722 | 0.021 | 4
$r$ | 16.6750 | 0.0057 | 4
$i$ | 16.7754 | 0.0062 | 4
$z$ | 16.958 | 0.016 | 4
$y$ | 16.951 | 0.011 | 4
$J$ | 16.2969 | 0.0074 | 5
$K_{s}$ | 16.188 | 0.035 | 5
$W1$ | 16.05 | 0.05 | 6
$W2$ | 15.88 | 0.10 | 6

Note. — [1] Gaia Collaboration et al. (2020), [2] Bailer-Jones et al. (2021),
[3] Bianchi et al. (2017), [4] Chambers et al. (2016), [5] McMahon et al.
(2013), [6] Schlafly et al. (2019)

Figure 1: Left: Periodograms for different selections of time-series
photometry, each vertically shifted for clarity, with the combined periodogram
shown in black at bottom. MHJS stands for the combination of McDonald, HAO,
JBO, and SAAO photometry, and horizontal red dash-dotted lines denote the 0.1%
false-alarm probability levels for each selection of data. The $P_{A}$ (9.937
hr) period and its first seven harmonics are denoted with vertical dashed
lines, while the vertical dotted lines denote the $P_{B}$ (11.2 hr) period and
its first four harmonics. Vertical red-hatched lines are placed at 11.6 and
23.2 $\mu$Hz to mark the 1-day and 0.5-day periods, respectively. Peaks at
these frequencies only show up in the ground-based data, and we consider them
to be artifacts caused by the diurnal sampling of those light curves. The
over-plotted blue periodogram at bottom represents the combined data with the
best-fit $P_{\mathrm{A}}$ and $P_{\mathrm{B}}$ values determined with Period04
prewhitened, along with 30 harmonics related to $P_{\mathrm{A}}$. Some
significant residual power is still present after prewhitening, especially at
harmonics of $P_{\mathrm{A}}$. For the MHJS periodogram, we excluded the
single McDonald light curve from 2020 October 19 to reduce the amount of
aliasing in the spectral window. Right: The spectral window of each data set,
with the same x-axis scaling as the left panel, and amplitudes scaled to the
highest peak in the corresponding periodogram.

###  2.3 SED Photometry

To construct a spectral energy distribution (SED) for ZTF J0328$-$1219, we
retrieved photometry from the Galaxy Evolution Explorer (GALEX, Martin et al.
2005; Morrissey et al. 2005), Gaia EDR3 (Gaia Collaboration et al., 2016,
2020), the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS1
or PS1, Chambers et al., 2016; Magnier et al., 2013), the VISTA Hemisphere
Survey (VHS, McMahon et al., 2013), and NASA’s Wide-field Infrared Survey
Explorer (WISE, Wright et al., 2010). A summary of all the catalogued
astrometry and photometry can be found in Table 4. For the WISE W1 and W2
photometry, we used the values reported in the unWISE catalogue (Schlafly et
al., 2019) which uses deeper imaging and improved modeling of crowded regions
compared to the ALLWISE catalogue (Cutri et al., 2021). For the distance, we
used the photogeometric value reported in Bailer-Jones et al. (2021).

##  3 Results

Table 5: Period Identifications

| Period04 | Gaussian Centroid
---|---|---
Parameter |  | ZTFg+r | TESS-04 | TESS-31 | MHJS† | Combined | LS‡ | BLS‡ | SH‡ | Plavchan
$P_{\mathrm{A}}$ | (hr) | $9.9374$ | $9.9314$ | $9.9243$ | $9.9324$ | $9.9370^{*}$ | $9.9246$ | $9.9225$ | $9.9161$ | $9.9191$
$\sigma_{\mathrm{A}}^{**}$ | (hr) | $\pm 0.0003$ | $\pm 0.0074$ | $\pm 0.0041$ | $\pm 0.0012$ | $\pm 0.0001$ | $\pm 0.0002$ | $\pm 0.0002$ | $\pm 0.0003$ | $\pm 0.0001$
$P_{\mathrm{B}}$ | (hr) | — | $11.549^{*}$ | $11.162^{*}$ | — | $11.1113$ | $11.1703$ | $11.1823$ | $11.1881$ | $11.1494$
$\sigma_{\mathrm{B}}^{**}$ | (hr) | — | $\pm 0.029$ | $\pm 0.011$ | — | $\pm 0.0002$ | $\pm 0.0002$ | $\pm 0.0006$ | $\pm 0.0017$ | $\pm 0.0001$
HWHM${}_{\mathrm{A}}^{\dagger\dagger}$ | (hr) | — | — | — | — | — | $0.052$ | $0.065$ | $0.057$ | $0.053$
HWHM${}_{\mathrm{B}}^{\dagger\dagger}$ | (hr) | — | — | — | — | — | $0.097$ | $0.111$ | $0.110$ | $0.093$

††footnotetext: The combined McDonald, HAO, JBO, and SAAO data.

‡‡footnotetext: LS = Lomb-Scargle Periodogram, BLS = Box Least-Squares
Periodogram, SH = Summed Harmonics FFT.

$*$$*$footnotetext: Values used to generate phase-folded light curves seen in
Figures 2, 3, and 4

$**$$**$footnotetext:  Formal least squares uncertainties for the Period04
periods and for the centroids of the periodogram peaks.

††††footnotetext:  The measured half width at half maximum of the periodogram
peaks which we take as a measure of the variations in the period and/or the
photometric modulation profile with time.

###  3.1 Periodic Variability

After the initial identification of ZTF J0328-1219 as a transiting planetary
debris candidate due to its metal pollution and short-timescale variability
(Guidry et al., 2021), we set out to investigate whether this object, like the
two known dust cloud transiting systems WD 1145+017 and ZTF J0139+5245,
exhibits periodic variability related to the orbital period of planetary
debris. The discovery of variability in ZTF J0328$-$1219 originally came from
just 3 hr of time-series photometry during a single night, during which no
significant periodicity was observed.

To identify possible periods on longer timescales, we first analyzed the light
curves from TESS Sectors 4 and 31. We used Astropy routines to calculate
unnormalized666See “PSD normalization” at the Astropy Website Lomb-Scargle
periodograms (LS, Lomb, 1976; Scargle, 1982) of the TESS light curves (see
Figure 1), which we then converted from power ($p$) to fractional amplitude
units ($A$) with the relation $A=\sqrt{4p/N}$, where $N$ is the total number
of light curve data points. We estimated the false-alarm probability level of
each periodogram using a bootstrap method (VanderPlas, 2018; Bell et al.,
2019), where we randomly resampled each light curve 10,000 times with
replacement, and with identical time samplings as the original light curves.
In 99.9% of cases, the maximum peak of the resulting periodograms did not
exceed 1.59% amplitude for Sector 4, and 0.85% for Sector 31.

In both TESS sectors, a significant peak is present at around 27.9 $\mu$Hz, or
9.93 hr (hereafter labeled the A-period, $P_{A}$), along with several higher
frequency peaks corresponding to harmonics of this period. In the TESS Sector
31 data, an additional significant peak at around 24.8 $\mu$Hz, or 11.2 hr
(hereafter labeled the B-period, $P_{B}$) is present, and also appears in the
Sector 4 data at lower significance. Some power appears to be present at the
first four harmonics of the B-period in TESS Sector 4, but do not show up in
TESS Sector 31.

Using the Period04 application (Lenz & Breger, 2005), we took the frequencies
and amplitudes of the highest peaks near 9.93 and 11.2 hr, and optimized them
with a non-linear least squares fit of two independent sine curves to the
light curves. We present the resulting periods and their formal least squares
uncertainties ($\sigma_{A}$ and $\sigma_{B}$, Montgomery & Odonoghue, 1999) in
Table 5. The resulting A-periods for the two sectors agree within 1-$\sigma$,
while the B-periods exhibit a small but significant difference of 0.39 hr, or
about 9.7$\sigma$. We also examined the stability of the periods and
amplitudes within each Sector by breaking the light curves up into five
roughly 5-day segments and performing a least squares fit to each segment
holding the A and B periods fixed at their best fit values for the whole
sector. Both periods exhibit phase stability consistent with no change in
period within each sector, though the B-period amplitude dropped significantly
in Sector 31 between the first and last segment, from 2.1$\pm$0.3% to
0.3$\pm$0.3% respectively.

In Figure 1, we also show the periodograms for the combined McD, HAO, JBO, and
SAAO data (labeled MHJS), for the combined ZTF $g+r$ data, and for the
combination of all data sets presented in this work. For the MHJS periodogram,
we only use data acquired from 2020 December through 2021 February to help
reduce cycle count errors, which meant excluding the single McDonald run
acquired in 2020 October. We use the same bootstrap method described above to
define a 0.1% false-alarm probability level for the ZTF data resulting in a
2.25% threshold. For the multi-site MHJS data, however, we adopt four times
the average periodogram amplitude between 0 and 300 $\mu$Hz as an estimate of
the 0.1 % false-alarm probability level (Breger et al., 1993; Kuschnig et al.,
1997), resulting in a 0.93 % threshold.

Both the ZTF and MHJS periodograms exhibit peaks aligned with the A-period and
some of its harmonics. Several significant peaks at daily aliases of the
A-period are also observed, which result from the spectral windows of the ZTF
and MHJS light curves. Normally these aliases would make it difficult to
determine the true period of variation, but the unambiguous detection of
variability at 9.93 hr by TESS helps identify the correct peaks in these
ground-based observations.

Figure 2: Individual light curves from McDonald Observatory (solid circles)
and SAAO (open circles) phased to the best-fit $P_{\mathrm{A}}$ value of
9.9370 hr. Each successive light curve is vertically shifted down for easier
comparison between nights and labeled by their corresponding dates. The number
of $P_{A}$ cycles that have occurred since $\mathrm{BJD_{TDB}}=2459141.504938$
at the start of each night is given in brackets next to each date.

We also searched for significant peaks corresponding to the B-period in the
ZTF and MHJS periodograms. After prewhitening $P_{A}$ and 30 of its harmonics
from the MHJS light curves, a small peak does show up near 23.5 $\mu$Hz (11.8
hr) just below the significance threshold. Given that $P_{B}$ changes by about
0.39 hr between TESS Sectors 4 and 31, it is possible that this peak is also
related to variability in ZTF J0328$-$1219\. However, this peak is also
consistent with being an alias of a ${\approx}$1-day period, that may result
from the daytime gaps in our ground-based observations coupled with slight
differences in the mean flux level between nights and facilities. Peaks
corresponding to a 1-day period (11.6 $\mu$Hz) and its nearby daily alias
(23.2 $\mu$Hz) also show up quite prominently in the ground-based ZTF data
(see peaks marked by vertical red-hatched lines in Figure 1). These peaks do
not show up, however, in either of the TESS light curves, which do not suffer
from daytime observational gaps. Thus, we consider any peaks that show up at
the 1-day period, or aliases thereof, systematic artifacts caused by the
diurnal sampling of our ground based data, and not related to variability
within ZTF J0328$-$1219.

We again used Period04 to optimize the periods for the ZTF, MHJS, and combined
light curves. For ZTF and MHJS, we fit only a single sine curve for $P_{A}$
since the $P_{B}$ period was not significant in their periodograms. The
resulting periods and formal least squares uncertainties are shown in Table 5.
These uncertainties are likely underestimated, as only the ZTF and TESS Sector
4 $P_{A}$ values agree within 3-sigma with the $P_{A}$ values determined from
the combined data set.

We also attempt period determinations with a separate method, using a Gaussian
function to fit the peaks in the power spectra. This method was motivated by
noticing that the observed periodogram peaks for the combined data set are
significantly broader than would be expected for stable periodic variability,
and that significant residual power remains near the A-period and its
harmonics even after prewhitening the best-fit period and many of its
harmonics from the data (see Figure 1). This suggests period or amplitude
modulations are present in this system, or that multiple closely spaced
periods exist near $P_{A}$ and $P_{B}$, any of which would make precise period
determinations difficult for this system.

We only perform Gaussian fits for the combined data set, and use a non-linear
least-squares routine to optimize the fit. We use the Gaussian center as a
measure of the average period, and the Gaussian half-width-half-max (HWHM) as
a measure of the effects that period modulations, amplitude modulations, or
multiple closely space periods may have on the peak widths. We perform this
fit on four different types of power spectra, a Lomb-Scargle periodogram, a
box least-squares periodogram, a summed harmonics fast FT777For the summed-
harmonics FFT, we assume that each frequency bin in a standard FT is the 10th
harmonic of a corresponding base frequency, and we add the amplitudes of all 9
lower harmonics to what is already in that bin. This summation of the first 10
harmonics can be generalized to any choice of N harmonics., and a Plavchan
periodogram (Plavchan et al., 2008), which may each be sensitive to different
features in the light curves. The results of the fits are summarized in Table
5, where we report the center and HWHM values for the Gaussian fits.

We also fit a Gaussian to the spectral window of the combined data set, and
find that on average the spectral window is about three times narrower than
the HWHM values for each period and for each type of power spectrum. The
A-periods from the Gaussian fits appear shifted to lower values relative to
the Period04 results, while the B-periods are closer to the Period04 results
for TESS Sector 31 and the combined data set.

With the above analysis, we attempt to convey that given the evolution of the
light curve modulation profile (see Section 3.2), the possible changes in the
periods over time, and the potential presence of multiple, closely-spaced
periods, a precise period determination is difficult to obtain for ZTF
J0328$-$1219\. Even if both $P_{\mathrm{A}}$ and $P_{\mathrm{B}}$ are highly
stable over time, changes in the light curve modulation profile will
themselves make precise determinations difficult. For simplicity, however, in
the remainder of this work we use the best-fit $P_{\mathrm{A}}$ value
determined using Period04 for the combined light curve, 9.9370 hr, to
construct phase-folded light curves and determine orbital parameters. Since
$P_{\mathrm{B}}$ shows a significant change between the two TESS sectors, we
use the best-fit periods determined with Period04 from each individual sector
to generate their phase-folded light curves (11.549 and 11.162 hr for Sectors
4 and 31, respectively, and opt to use the Sector 31 B-period to determine
orbital parameters.

###  3.2 Phase-Folded Light Curves

Figure 3: Top: Light curves from TESS (red and green triangles), ZTF (orange
pentagons), HAO (blue open circles), JBO (purple stars), and McDonald
Observatory (dark blue circles) showing the total time baseline covered by
each data set. The TESS Sector 31 data shown in this panel are binned into
10-minute exposures to reduce scatter. Bottom: Light curves phase folded on
the best-fit A-period, 9.9370 hr, shown in rough chronological order from top
to bottom, and labeled at right by their source. The TESS Sector 4 and 31
light curves are divided into first- and second-half segments to show changes
in light curve structure over their durations, while the ZTF data are divided
into three groups by observing season. The HAO and McDonald data are grouped
by month. For the ZTF light curves, we display the $g$-band (filled symbols)
and $r$-band (open symbols) separately. Each light curve has been vertically
shifted for clarity, and two full phases have been plotted.

In Figure 2 we present the detailed phased light curves obtained from McDonald
and SAAO, and in Figure 3, we present the binned, phase-folded light curves
obtained from all of our photometry sources. The light curves have been folded
on the best-fit A-period, $P_{\mathrm{A}}=9.9370\,$hr, identified in Section
3.1 and marked with an asterisk in Table 5. The phases are given with respect
to an arbitrary date, $\mathrm{BJD_{TDB}}=2459141.504938$, chosen close to the
start of the first McDonald run on 2020 October 19. We only generate binned,
phase-folded light curves averaged over several days and nights for TESS, ZTF,
HAO, and JBO due to having lower S/N or sparse sampling. TESS Sector 31 would
be an ideal data set to observe cycle-to-cycle variations in dip structures
due to its 2-minute cadence and near-continuous coverage over about 26 days,
but ZTF J0328$-$1219 is too dim to observe individual dips within a single A
or B-period cycle. For example, the largest variations observed in ZTF
J0328$-$1219 have amplitudes of about 10%, but the average point-to-point
scatter among the TESS Sector 31 data points in about 18%, so noise clearly
dominates unless averaged out over several cycles.

McDonald light curves from the four consecutive nights in 2020 December and
the four consecutive nights in 2021 February each span nearly 8 full cycles of
$P_{\mathrm{A}}$, and show both broad and narrow dip-like features with depths
of 5$-$10% that line up in phase but also vary significantly in depth and
shape between nights. This behavior is reminiscent of the dynamic orbit-to-
orbit variations in dip structure seen in WD 1145+017 (Gänsicke et al., 2016;
Rappaport et al., 2016; Gary et al., 2017). The narrowest of the McDonald dips
have durations of about 20 minutes, while the broadest features extend more
than an hour. In contrast to the observed behavior of WD 1145+017 which often
shows intervals of inactivity between successive dip events, ZTF J0328$-$1219
appears to exhibit constant photometric variability when observed at both high
speed and high signal-to-noise, likely the result of being occulted by dust
clouds populating a wide range of phases within the orbit.

In Figure 3, the phase-folded and binned light curves from individual data
sets, which span a total of about 2.5 years, also show variations in light
curve structure over a broad range of timescales. We split each TESS light
curve into first and second half components, each between 7.4 and 12.2 days
long, to search for changes in light curve structure on these timescales.
Between the first and second halves of TESS Sector 31, a large dip appeared at
phase 0.2, which then disappeared again in the HAO observations which began
just 16 days after TESS Sector 31 completed. At the same time, a sharp dip
feature around phase 0.45 seems to persist from Sector 31 through our McDonald
data (about 54 days total), only recently becoming less pronounced in JBO
observations which began on 2020 December 16.

Since high-speed observations have revealed that dips can occur on timescales
as short as 20 minutes, and changes in the modulation profile can occur on day
to week-long timescales, the 30-min long exposure times of the TESS Sector 4
data will inevitably tend to smooth the light curve structure, while the
$\approx$7-month long time baselines covered by each season of ZTF data will
average out changes in the modulation profile. These effects are noticeable in
their folded light curves which show considerably less detail, but still show
general trends which are in agreement with the more resolved structures seen
in the TESS Sector 31, McDonald, HAO, and JBO light curves. For example, the
regions of lowest flux in TESS Sector 4 and ZTF are aligned in phase with the
regions of highest dip activity in the other light curves.

Figure 4: Left: The TESS Sector 4 (top) and Sector 31 (bottom) light curves
phase-folded on the best-fit B-periods determined from each sector, and
vertically shifted for clarity. We again split each TESS light curve into
first- and second-half components to display changes that occur throughout
their ${\approx}$26-day durations. The amount of data contained within each
phase-folded segment is shown to the left of each light curve in parentheses.
Right: The corresponding periodograms of each TESS sector, where solid lines
indicate the first-half light curves, and dashed lines the second-half.
Frequencies corresponding to the A and B periods ($f_{A}$ and $f_{B}$) are
marked with vertical dotted lines. In TESS Sector 31, the B-period underwent a
significant reduction in amplitude from 2.3$\pm$0.2% in the first half to
0.8$\pm$0.2% in the second half.

Despite our extensive follow-up observations, we lack a conclusive $P_{B}$
detection in the periodograms of the MHJS photometry, nor do we readily
identify any persistent features in our highest S/N light curves from McDonald
and SAAO when folding on the B-period. This is possibly due to the relatively
short durations of our individual observing runs, which fail to capture the
full 11.2 hr B-period in one night. Thus, if the 11.2 hr period results from
one or more localized dip features, these dips may not always be observed each
night. On the other hand, if the 11.2 hr period results from broad features,
these may be diminished in amplitude by the routines we use to remove airmass
trends in some of our nightly runs. Additionally, given the changes in period
and harmonic structure of the $P_{B}$ peaks between the two TESS sectors, it
is also possible that the activity level at this period changes with time, and
our follow-up observations occurred during an interval of low activity at this
period.

Since the TESS observations cover a large number of both $P_{A}$ and $P_{B}$
cycles without any daytime interruptions, they offer a better chance at
averaging out variations related to the A-period and revealing structure when
phase-folded on the B-period. In Figure 4, we show the TESS light curves
phase-folded on the best-fit $P_{B}$ values determined using Period04 for each
sector (see Table 5). We again divide each sector into first and second halves
to reveal changes in structure over the ${\approx}$26-day duration of each
sector, and also show the periodograms associated with each TESS light curve
segment.

In both halves of the Sector 4 light curves, the B-period variations are only
weakly detected, but they are clearly detected in the first half of Sector 31
where they appear nearly sinusoidal with an amplitude of 2.3$\pm$0.2%. In the
second half of Sector 31, however, the B-period amplitude decreased to
0.8$\pm$0.2%, potentially revealing a rapid reduction in the activity level
associated with this period just 16 days prior to the start of our HAO
observations. It is possible the activity level remained low throughout our
follow-up observations with HAO, McDonald, JBO, and SAAO, though this is
difficult to confirm given the previously mentioned issues surrounding our
ability to detect $P_{B}$ variability with ground based observations. Future
ground-based observations with continuous coverage, such as with the Whole
Earth Telescope (Nather et al., 1990), would greatly aid in the detection and
detailed characterization of the $P_{B}$ modulations.

###  3.3 White Dwarf Atmospheric Properties

Previous estimates of the effective temperature and mass of ZTF J0328$-$1219
have been attempted using Gaia DR2 photometry and parallax while assuming a
pure H or He atmosphere (Gentile Fusillo et al., 2019). Since being identified
as a metal-polluted He-atmosphere white dwarf whose dominant spectral features
come from elements heavier than He (DZ), these estimates are no longer
reliable for ZTF J0328$-$1219\. Using the same LDT spectrum presented in
Guidry et al. (2021), along with broadband optical and near-infrared
photometry from Gaia eDR3 and PS1, we determine new atmospheric parameters for
ZTF J0328$-$1219.

Table 6: Atmospheric Properties of ZTF J0328$-$1219  Parameter | Value | Uncertainty
---|---|---
$T_{\mathrm{eff}}$ (K) | 7630 | 140
$\log(g)$ (cgs) | 8.245 | 0.035
$M_{\star}$ ($M_{\odot}$) | 0.731 | 0.023
$R_{\star}$ ($R_{\odot}$) | 0.0107 | 0.0002
[Ca/He] |  $-$9.55 | 0.12
[H/He] |  ${<}{-}$3.5 | –
$\tau_{\mathrm{cool}}$ (Gyr) | 1.84 | 0.17
Figure 5: The LDT spectrum (black) along with the best-fit model (red), with
an inset zoomed in on the Ca H and K features.

We follow the fitting procedure outlined in Coutu et al. (2019) and employ the
DZ model atmospheres of Blouin et al. (2018). We first derive an effective
temperature, $T_{\mathrm{eff}}$, and solid angle, $\pi R^{2}/D^{2}$, using the
Gaia eDR3 and PS1 photometry. From this solid angle and the Gaia parallax, we
obtain the white dwarf radius. We then evaluate the corresponding white dwarf
mass and $\log(g)$ using theoretical white dwarf structure models (Fontaine et
al., 2001), assuming a C/O core, a helium envelope of $\log(M_{\rm
He}/M_{\star})=-2$, and a thin hydrogen envelope of $\log(M_{\rm
H}/M_{\star})=-10$. Keeping $T_{\mathrm{eff}}$ and $\log(g)$ fixed to those
values, we then adjust the calcium abundance, [Ca/He]
$=\log(n_{\mathrm{Ca}}/n_{\mathrm{He}})$ where $n$ is the number density, to
the LDT spectrum. H is included in our model atmospheres assuming an abundance
[H/He] $=\log(n_{\mathrm{H}}/n_{\mathrm{He}})$ that corresponds to the
visibility limit of the H$\alpha$ line (Coutu et al., 2019), which is not
detected in the LDT spectrum. In addition, all elements from C to Cu are
included assuming chondritic abundance ratios with respect to Ca. An iterative
procedure which alternates between these two steps (the photometric fit and
the adjustment of the abundances to the spectrum) is then performed to
converge on a best-fit solution. The final values are summarized in Table 6,
while Figure 5 shows the DZ model fit to the LDT spectrum, while Figure 6
shows the DZ model fit the to SED.

Figure 6: Top: The observed SED for ZTF J0328$-$1219 including GALEX (purple),
Pan-STARRS1 (blue), Gaia eDR3 (green), VHS (orange), and unWISE (red)
photometry. The best-fit white dwarf model is plotted in grey, with calculated
photometric fluxes in each passband shown with dark grey squares. The dark and
light-colored error bars show the 1-sigma and 3-sigma uncertainties,
respectively, for each measurement. We only used the PS1 and Gaia photometry
to constrain the white dwarf atmospheric fit, which produces a slightly poorer
fit to the VHS and WISE photometry. Bottom: Residuals between the observed and
model flux values in each passband. We find no evidence for an infrared excess
in the WISE $W1$ and $W2$ bands within their 3-sigma uncertainty limits,
though a small but significant excess is present in the VHS $J$-band.

For ZTF J0328$-$1219 we derive $T_{\mathrm{eff}}=7630\pm 140\,$K,
$M_{\star}=0.731\pm 0.023\,M_{\odot}$, and $R_{\star}=0.0107\pm
0.0002\,R_{\odot}=1.17\pm 0.02\,R_{\earth}$. Due to the inclusion of
atmospheric metals and H, this is significantly cooler and less massive than
the pure He-atmosphere parameters derived by Gentile Fusillo et al. (2019)
using Gaia DR2 photometry and parallax ($T_{\mathrm{eff}}=8551\pm 161\,$K and
$M_{\star}=0.855\pm 0.029\,M_{\odot}$). This trend is in line with previous
results (e.g., Coutu et al., 2019) and is due to the increased He- free–free
opacity resulting from the additional free electrons provided by the metals
and H. We chose not to use the VHS and unWISE photometry in our fitting
process so we could better assess whether an infrared excess is present. If we
include these data in the fit, however, we find a 2-sigma drop in
$T_{\mathrm{eff}}$ by about 300 K, producing a better fit to the VHS and
unWISE photometry while slightly degrading the fit to the Gaia and PS1
photometry.

We estimate a Ca abundance of [Ca/He] $=-9.55\pm 0.12$, consistent with
abundances of known metal-polluted He-atmosphere white dwarfs in this
temperature range (Blouin, 2020; Coutu et al., 2019; Hollands et al., 2017;
Dufour et al., 2007), and put an upper limit on the H abundance of [H/He]
$<-3.5$ based on the non-detection of Balmer lines in the LDT spectrum.
Assuming the white dwarf has followed a standard single-star evolutionary
pathway and has a thin H envelope ($\log M_{\rm H}/M_{\star}=-10$), we find a
white dwarf cooling age of $\tau_{\mathrm{cool}}=1.84\pm 0.17\,$Gyr using the
evolution models described in Blouin et al. (2020).

Using the He convection zone masses ($M_{\mathrm{CVZ}}$) and Ca diffusion
timescales ($\tau_{\mathrm{Ca}}$) for helium-atmosphere white dwarfs in
Koester (2009), we find $\log(M_{\mathrm{CVZ}}/M_{\star})=-4.75$ and
$\log(\tau_{\mathrm{Ca}}[\mathrm{yr}])=6.46$ at the effective temperature of
ZTF J0328$-$1219\. Assuming heavy metal abundance ratios outlined in Zuckerman
et al. (2010), this implies that $4.4\times 10^{21}\,$g of heavy metals
currently reside within the outer He convection zone, or about 0.5% the mass
of Ceres.

Steady-state accretion episodes around white dwarfs are expected to last up to
$10^{4}-10^{6}\,$yr (Girven et al., 2012), and often times are expected to be
much shorter (Veras & Heng, 2020), making it highly unlikely that the current
[Ca/He] is being sustained with steady-state accretion. More likely, the
[Ca/He] represents some combination of residual Ca from previous accretion
episodes, and new Ca being presently accreted. Steady-state accretion would
imply a total mass accretion rate of $\log(\mathrm{\dot{M}[{\rm
g\,s}^{-1}]})=7.7$, consistent with accretion rates for white dwarfs with
similar temperatures and atmospheric compositions (Dufour et al., 2007;
Farihi, 2016; Hollands et al., 2017; Xu et al., 2019), but the current
accretion rate might actually be orders of magnitude higher or lower. Compared
to the steady-state accretion rate of WD 1145+017 at
$\log(\mathrm{\dot{M}[{\rm g\,s}^{-1}]})=10.6$ (Xu et al., 2016) whose
diffusion timescales are of order $10^{5}$ yr, however, ZTF J0328$-$1219 is
much lower. This might be the result of less overall accretion in ZTF
J0328$-$1219 due to having lower-mass bodies in orbit, or perhaps due to the
bodies orbiting further away from the white dwarf and undergoing less
disruption, but it is also possible that past accretion episodes have
contributed to the currently observed abundances in these objects given their
long diffusion timescales.

We note that even though the DZ model provides a good fit to the observed Ca H
and K lines, some caution must still be taken when interpreting the calculated
Ca abundance. WD 1145+017 exhibits prominent circumstellar absorption
components (Xu et al., 2016) which evolve significantly over time (Redfield et
al., 2017; Cauley et al., 2018). We do detect some narrow atomic absorption
features consistent with circumstellar absorption in our high-resolution MIKE
spectrum (see more detailed analysis in Section 3.4), the strongest of which
have equivalent widths of about 0.2Å. This is small compared to the
${\approx}45$Å combined equivalent width of the photospheric Ca H and K lines,
so we suspect the impact of circumstellar absorption features on our measured
[Ca/He] is negligible.

We also note a weak H$\alpha$ feature is present in our high-resolution MIKE
spectrum (see Section 3.4 and Figure 9) and we attempt to include it in our
atmospheric modeling of the white dwarf for better constraints on [H/He], but
the narrow line is poorly fit with current atmospheric models. To match the
depth of the narrow H$\alpha$ feature, our models suggest that a broad
H$\alpha$ component should also be present, which is not observed. Such
mismatches have been seen before in metal polluted white dwarfs, such as the
narrow He-I lines observed in GD 362 (Zuckerman et al., 2007; Tremblay et al.,
2010) and in GD 16 and PG 1225-079 (Klein et al., 2011). Therefore, we do not
consider the mismatch a sign of a non-atmospheric origin for H$\alpha$, but
rather a difficulty with atmospheric modeling. We still believe that our
calculated upper limit for [H/He] of $-$3.5 at the photosphere is valid since
it is based on the lack of a broad H$\alpha$ absorption feature in the low-
resolution LDT spectrum.

###  3.4 Circumstellar Gas and Dust

In addition to atmospheric metal pollution, other observational signatures of
planetary debris around white dwarfs include infrared excesses due to the
presence of warm dust disks (e.g. Zuckerman & Becklin, 1987; Dennihy et al.,
2020a), metallic emission lines from hot gaseous components of the dust disks
(e.g. Gänsicke et al., 2006; Melis et al., 2020), and metallic absorption
lines from circumstellar gas (e.g. Debes et al., 2012; Xu et al., 2016). WD
1145+017, the prototypical transiting debris white dwarf, exhibits an infrared
excess (Vanderburg et al., 2015) and circumstellar gas absorption (Xu et al.,
2016), but has not been observed to exhibit metallic emission lines.

We find no evidence for gaseous emission in any of our new spectroscopic
observations from MIKE or SOAR (see Figure 7). Also, in the WISE $W1$ and $W2$
bands, the observed flux agrees within 3-sigma with the best-fit DZ model (see
Figure 6). Using unWISE $W1$ and $W2$ photometry, Xu et al. (2020) also do not
find a significant magnitude or color excess for ZTF J0328$-$1219, although
they assumed a pure hydrogen-atmosphere white dwarf model. We do find a small
but significant excess of 0.04 mJy (11.1-sigma) in the VHS $J$-band. Larger
$J$-band flux excesses relative to $W1$ and $W2$ excesses are sometimes seen
for white dwarfs with cool sub-stellar companions (see Figure 1 of Xu et al.,
2020), but this may also be the result of underestimated $J$-band
uncertainties coupled with a white dwarf atmospheric fit that is slightly too
hot. If we include the VHS and unWISE photometry in the atmospheric fitting
process described in Section 3.3, the best-fit white dwarf model becomes about
300K cooler, and the VHS $J$-band then agrees within 4.3-sigma. To further
investigate whether an infrared excess is present in ZTF J0328$-$1219,
higher-S/N observations in the near and mid-infrared, such as with JWST, will
be needed.

Figure 7: The SOAR (top) and MIKE (bottom) continuum-normalized spectra
centered on the Ca infrared triplet region, and vertically shifted for
clarity. The best-fit white dwarf model spectrum convolved with a Gaussian
kernel to the respective resolutions of SOAR and MIKE is overplotted in blue
and velocity-shifted to match the measured velocity of the H$\alpha$ line
(58.8 km s-1). Vertical black lines denote the rest-wavelengths of each Ca
infrared triplet component. CCD fringing is prevalent in the MIKE spectrum at
these wavelengths, producing many of the broad features not observed in the
SOAR spectrum, though detections of absorption at 8542Å, and possibly at
8662Å, are seen in both. We find no evidence for Ca triplet emission in these
observations.

We detect prominent and narrow Na-I D lines in the MIKE spectrum of ZTF
J0328$-$1219, with a full-width-half-maximum (FWHM) of $0.24$Å and equivalent
widths of $0.18\pm 0.01\,$Å and $0.24\pm 0.01\,$Å for D1 (5895Å) and D2
(5889Å), respectively (see Figure 8). To determine whether these lines are of
interstellar or circumstellar origin, we compare their line shapes and
velocities to the Na D lines observed in our nearby comparison stars. Using
the 100 $\mu$m IRAS dust maps (Schlegel et al., 1998) and extinction values
calculated by Schlafly & Finkbeiner (2011), we find that the line-of-sight
Galactic extinction at the coordinates of ZTF J0328$-$1219 is $A_{V}=0.221$.
At a distance of just 43 pc, however, it is unlikely that a significant amount
of ISM material lies between us and ZTF J0328$-$1219, and the 3D dust maps of
(Green et al., 2019) predict no extinction at the coordinates and distance of
ZTF J0328$-$1219\. The two comparison stars nearest to ZTF J0328$-$1219 on
sky, Gaia J0328$-$1216 and HD 21634, are the most suitable for comparing Na
line shapes and velocities since they would likely sample the same ISM
material and have the most similar line-of-sight Galactic extinctions compared
to ZTF J0328$-$1219 with $A_{V}=0.234$ and $A_{V}=0.253$, respectively. Both
comparison stars are about three times further away than ZTF J0328$-$1219 so
they might experience more ISM absorption, but the 3D dust maps of Green et
al. (2019) again predict no extinction at the coordinates and distances of
these objects. The same is true for the two comparison stars farthest from ZTF
J0328$-$1219 on sky, HD 22243 and HD 21917, but which have lower line-of-sight
Galactic extinctions ($A_{V}$ of 0.091 and 0.156, respectively).

Figure 8: The continuum-normalized MIKE spectra centered on the Na D doublet
for ZTF J0328$-$1219 (solid black line) and the two nearest comparison stars,
HD 21634 (dotted blue line) and Gaia J0328$-$1216 (dash-dot red line). The
vertical dashed lines denote the rest wavelengths of the Na D doublet. Neither
comparison star shows any Na D absorption — in addition to their photospheric
components — at the same velocity as those in ZTF J0328$-$1219, ruling out an
interstellar origin.

To establish our comparison stars as reliable radial velocity (RV) references,
we first compared their measured velocities to catalogued Gaia DR2 values. One
RV standard, GJ 908 (spectral type M1 at $-$71.13 km s-1), was observed on
2020 December 1. We measured the velocities of our comparison stars from 22
orders of our blue spectra by cross-correlation with GJ 908. Our measured
radial velocities are close to those in DR2: Gaia J0328$-$1216: $1.08\pm 0.14$
versus $2.3\pm 0.93$ km s-1, HD 21634: $-6.35\pm 0.42$ versus $-7.37\pm 0.49$
km s-1, and HD 21917: $-58.0\pm 0.2$ versus $-57.9\pm 0.18$ km s-1. We exclude
HD 22243 from our RV analysis due to its exceptionally broad lines, though we
still use it to search for narrow Na features indicative of ISM absorption. We
also measured the velocities of several individual lines in our comparison
stars by fitting them with Voigt profiles. We did this for the Na-I D doublet,
the Ca-I 6717Å line, and H$\alpha$, and found velocities in each case that are
consistent with the cross-correlation velocities, confirming photospheric
origins for these lines. Hereafter, we use our measured, cross-correlation RVs
to compare with ZTF J0328$-$1219.

Table 7: Radial Velocities of ZTF J0328$-$1219 from MIKE

Line | $\lambda_{\mathrm{rest}}$ | RV | Origin
---|---|---|---
| (Å) | (km s-1) |
H$\alpha$ | 6562 | 58.8$\pm$0.8 | atmospheric
Na-I | 5889, 5895 | 37.4$\pm$0.5 | circumstellar
Ca-II† | 8542 | 40.1$\pm$1.5 | circumstellar
Ca-II‡ | 3933, 3968 | 51.5$\pm$3.5 | blended?

††footnotetext: Velocity shown is for the narrow line component.

‡‡footnotetext: Velocity shown is for the narrow H and K line cores.

Figure 9: The continuum-normalized MIKE spectrum (solid black) of ZTF
J0328$-$1219 near various lines with measured velocities (see Table 7). For
the 8542Å Ca triplet component, we also show the continuum-normalized, lower-
resolution SOAR spectrum (dashed red) which detects the same feature. The 3933
and 3968Å panels (Ca-II K and H respectively) highlight the narrow line
components, but also exhibit a small portion of the broad components of these
lines which cause the upward slopes away from line center. Vertical lines
denote the measured velocities of the Na and H$\alpha$ features from MIKE
spectra, with the dashed line showing the 37.4 km s-1 velocity of the Na
lines, and the dash-dotted line showing the 58.8 km s-1 velocity of the
H$\alpha$ feature.

We measured the absolute radial velocity of the Na D lines in ZTF J0328$-$1219
by cross-correlation against the stellar Na lines of four comparison stars: HD
21917, HD 21634, Gaia J0328$-$1216, and GJ 908. We used only the wavelength
region between 5886 and 5899 Å for cross correlation, and we do not perform
any corrections to the RVs for the gravitational redshift caused by the white
dwarf. The average absolute radial velocity of the Na D lines in ZTF
J0328$-$1219 is $37.4\pm 0.5$ km s-1.

We also detect a weak H$\alpha$ line in ZTF J0328$-$1219 with an equivalent
width of only $0.08\pm 0.01\,$Å, but a FWHM of 0.6 Å (see Figure 9). We
measure its velocity to be $58.8\pm 0.8$ km s-1, a difference of $21.4\pm 1.0$
km s-1 with respect to the Na D features in the MIKE spectrum. In the blue-
side MIKE spectrum, we also detect Ca H and K lines with broad wings and
narrow components at their cores, consistent with the LDT observations. The
S/N in the Ca H and K order is not good enough, however, to measure the
velocities precisely, and cross-correlation with the reference stars does not
work well because of the vastly different shapes of the lines. Still, we
attempt to measure the H and K velocity by fitting a Voigt profile to the line
core and obtain a velocity of $51.5\pm 3.5$ km s-1. We caution that in
addition to the noisy detection, the cores of these lines may be blended with
circumstellar absorption features if present, further complicating this
velocity measurement.

Two of the Ca infrared triplet lines, 8542 and 8662Å, are also detected in the
MIKE spectrum. The spectrum is considerably noisier at these wavelengths and
significant CCD fringing is present, but we detect similar features in our
lower-resolution SOAR spectrum, which has only 0.3% peak-to-peak fringing at
these wavelengths, leading us to believe these features are real. The stronger
of the two lines at 8542Å has a narrow component with FWHM ${\sim}0.25$Å,
similar to that of the Na D lines, superimposed on a broader, slightly redder
component (see Figure 9). By simultaneously fitting the narrow and broad line
components with two Voigt profiles of different widths, and comparing them to
the position of the stellar lines in the comparison stars, we measure the
velocity of the narrow and broad components of the Ca 8542Å line as $40.1\pm
1.5$ and $59.5\pm 3.5$ km s-1, respectively. In our SOAR spectrum, the 8542Å
feature is similar in strength to the MIKE feature, though the narrow
component is not detected at the lower resolution (see Figure 9). We once
again fit a Voigt profile to the 8542Å component in the SOAR spectrum, and
measure a velocity of $59.6\pm 8.9$ km s-1, in agreement with the broad
component observed by MIKE.

The Na D lines provide the most compelling evidence for a detection of
circumstellar gas absorption. They are much stronger than predicted by our DZ
model, and much narrower than would be expected if photospheric. We can also
rule out interstellar absorption since none of the comparison stars show
similarly narrow Na D absorption — in addition to their photospheric
absorption — at the same velocity as ZTF J0328$-$1219 (see Figure 8). This
includes the two nearest comparison stars, HD 21634 and Gaia J0328-1216, which
are most likely to sample the same ISM material as ZTF J0328$-$1219, and also
the two comparison stars further away, HD 22243 and HD 21917.

The H$\alpha$ feature is likely to be in the white dwarf atmosphere given its
velocity is 21 km s-1 red-shifted with respect to the Na D and calcium triplet
lines. This difference is smaller than expected since the gravitational
redshift for a 0.73 $M_{\odot}$ white dwarf would be about 44 km s-1, but this
also assumes that the Na D lines come from circumstellar gas traveling
perpendicular to our line of sight towards the white dwarf, which may not be
the case if the gas has an eccentric orbital configuration. Smaller velocity
differences may also be observed if the gas is located sufficiently deep in
the gravitational well of the white dwarf, but this effect has only been
observed for high ionization species (Gänsicke et al., 2012), and likely does
not apply here. Long-term spectroscopic monitoring and modeling of the Na D
line velocities, such as the studies that have been done for WD 1145+017
(Redfield et al., 2017; Cauley et al., 2018; Fortin-Archambault et al., 2020)
and WD 1124-293 (Debes et al., 2012; Steele et al., 2020), can be used to
probe the geometry and column densities of the circumstellar gas in ZTF
J0328$-$1219.

The narrow component of the Ca-II 8542Å feature is likely of circumstellar
origin as well, since its velocity is most in agreement with the Na D lines
and has a similar FWHM, but the broader component detected by both MIKE and
SOAR is more in agreement with the H$\alpha$ velocity. It is currently unclear
if that suggests some of the Ca infrared triplet absorption has an atmospheric
origin, since the DZ model predicts much weaker absorption than what is
observed. It is possible that the circumstellar gas exhibits an asymmetric
Doppler-broadened line profile that can account for both of the 8542Å line
components (e.g. Xu et al., 2016), but it would be puzzling why a similar
velocity distribution is not observed for the much stronger Na D lines. With
continued spectroscopic monitoring at high resolution, it will be interesting
to see if any velocity variations in the Na D lines also occur to one or both
components of the Ca-II 8542Å line.

##  4 Discussion

###  4.1 Orbiting Debris

The evidence for planetary debris in orbit around ZTF J0328$-$1219 consists of
atmospheric metal pollution, periodic optical variability at 9.93 and 11.2 hr,
and the presence of circumstellar gas. As seen in Figures 2 and 3, the
characteristics of the optical variability closely resemble those in WD
1145+017 (see Gänsicke et al., 2016; Rappaport et al., 2016, 2018), where
small variations in the detailed transit shapes are likely occurring on orbit-
to-orbit timescales, while more dramatic changes can occur on week, month, and
year-long timescales. In spite of these similarities, ZTF J0328$-$1219 also
exhibits unique characteristics which bring into question the source of its
variability. Notably, ZTF J0328$-$1219 appears almost continuously variable
(see Figure 2), whereas WD 1145+017 in most cases exhibits intervals of
quiescence between successive dip events. This may suggest that ZTF
J0328$-$1219 is continuously being occulted by clouds of planetary debris
which populate a wide range of orbital phases, but here we first consider
alternate explanations before exploring the planetary debris scenario further.

###  4.2 Other Possible Explanations for Variability

####  4.2.1 Rotational Modulation

White dwarfs are known to exhibit periodic photometric variability due to
stellar rotation combined with magnetism or some form of surface inhomogeneity
(see e.g. Maoz et al., 2015; Kilic et al., 2015; Hermes et al., 2017a; Reding
et al., 2020; Gänsicke et al., 2020). Typical white dwarf rotation periods
range from 0.5–2.0 days (Hermes et al., 2017b), so both the 9.93 and 11.2 hr
periods are feasible options for a rotation period in ZTF J0328$-$1219\.
Already, however, the presence of two periods eliminates rotational modulation
as the sole source of variability, since it can only account for one of the
observed periods. Additionally, the vast majority of known and suspected
rotationally modulated white dwarfs exhibit persistent and near-sinusoidal
variations in their light curves, but ZTF J0328$-$1219 exhibits both highly
non-sinusoidal variations and rapid evolution in light curve structure on
timescales of 1-day or less.

Due to the complexity of the light curve structure, however, we have yet to
identify individual dip-like features in our ground-based followup
observations associated with the 11.2 hr B-period. The TESS light curves show
a significant detection at this period, but are too noisy to analyze
individual B-period cycles in detail. The phase-folded light curves show some
details of the B-period modulation profile, which appears to change
significantly in both shape and amplitude between TESS sectors and between the
first and second halves of individual sectors (see Figure 4), behavior which
is not expected for rotationally modulated white dwarfs. These changes are
also evidenced by the presence of harmonics in the TESS Sector 4 periodogram,
indicating sharper features in the light curve, which are not seen in the TESS
Sector 31 periodogram. Lastly, the broader-than-expected peaks seen in the
periodogram (see Section 3.1) also suggest changes to the modulation profile,
changes in the period, or the presence of multiple closely spaced periods, all
of which rule against rotational modulation. Therefore, we strongly suspect
the 11.2 hr period to most likely be associated with transiting dust-emitting
debris, but we still consider rotational modulation a possibility for $P_{B}$
until more detailed observations are carried out which can better characterize
the 11.2 hr modulations.

####  4.2.2 Cataclysmic Variable

The near-constant photometric variability in ZTF J0328$-$1219 looks
qualitatively similar to the photometric variability observed in some
cataclysmic variables (CV). Unlike CVs, however, ZTF J0328$-$1219 does not
exhibit any spectral emission lines that would indicate mass accretion from a
companion stellar or sub-stellar object (e.g. Szkody et al., 2009), or
emission from the chromosphere of an irradiated companion (e.g. Longstaff et
al., 2017, 2019).

In addition, the presence of a sub-stellar companion is strongly disfavored
due to the lack of a significant infrared excess (Figure 6). We generated
composite SEDs for a white dwarf and L/T-dwarfs using the best-fit DZ model
from Section 3.3 and L and T-dwarf templates from the SpeX Prism library
(Burgasser, 2014). We combined the spectra following the procedure described
in Casewell et al. (2017), normalizing the white dwarf model to absolute
$g$-band magnitudes from Holberg & Bergeron (2006), and the L and T-dwarf
templates to absolute $J$-band magnitudes from Dupuy & Liu (2012).

The SpeX templates extend out to just 2.6 $\mu$m, so we exclude the WISE
photometry from our comparison with the composite SEDs. Within the 3-sigma
uncertainty limits of the $J$-band and $K_{\mathrm{s}}$-band photometry, we
can exclude the presence of a companion with spectral type earlier than T7. We
have not carried out any time series spectroscopy to look for radial velocity
variations, but we note that the radial velocity of the Ca 8542Å line from the
SOAR spectrum and the broad component of the Ca 8542Å line from the MIKE
spectrum are consistent with no change in velocity.

###  4.3 Source of the Dusty Effluents

Ultimately, we find that the periodic variability at both the 9.93 and 11.2 hr
periods is best explained by planetary debris transiting the white dwarf. We
further postulate that the photometric variability itself is caused by dust,
either entrained in the orbiting debris or emitted continuously by one or more
discrete bodies. This general scenario is supported by the presence of
atmospheric metal pollution, the highly non-sinusoidal and volatile nature of
the periodic photometric variability, and the detection of circumstellar gas.

The orbital semi-major axes corresponding to periods of 9.937 and 11.162 hr
are $2.11\,R_{\odot}$ and $2.28\,R_{\odot}$, respectively888If we assume that
the debris is being fed by the disruption of two planetesimals orbiting at
$2.11\,R_{\odot}$ and $2.28\,R_{\odot}$ from the host star, respectively, we
can set a rough upper limit on the masses of the these bodies if we require
that the two orbits be separated by more than $\sim$10 Hill’s sphere radii of
the bodies. This is a minimum requirement for long-term dynamical stability.
From this, we can conclude that the masses of these bodies must be less than
Jupiter’s mass, which is not interestingly restrictive.. At these distances,
and in the case where no shielding of the white dwarf flux by intervening dust
or gas occurs, the equilibrium temperatures of absorbing bodies could range
from 385 K to 700 K for orbital periods near 9.93 hr, and somewhat lower for
orbiting bodies with periods near 11.2 hr. These temperatures are summarized
in Table 8, and raise the issue of what kinds of materials could survive and
continue to emit dusty effluents in this environment. We have also included in
Table 8 the Chiang & Goldreich (1997) formalism for equilibrium temperatures
within a flat, optically thick disk that can partially shield an orbiting body
from direct radiation coming from the host star. In order for this to be
applicable, there must be a dusty disk, and the body should be physically
smaller than the disk height. We currently have no direct evidence for the
former and little information about the latter.

In order to guide our thinking about what types of bodies might emit dust at
these temperatures, we consider the properties, including composition, of
comets and asteroids. In Table 9 we present a broad-brush and simplified, but
hopefully useful, comparison between these two types of bodies that are seen
in the Solar System. Comets are mostly fragile material whose main volatiles
are H2O, CO, and CO2 ices. The sublimation of these bodies occurs well below
the environmental temperature of ZTF J0328$-$1219, and hence would not last
very long before completely evaporating. By contrast, many of the minerals
expected to be found in asteroids, some examples of which are given in Table
9, have sublimation temperatures above about 1200 K. These temperatures are
too high to expect substantial emission of heavy molecules that might condense
into dust in $\sim$10-hr obits around ZTF J0328$-$1219.

We know of no obvious candidate materials with sublimation temperatures in the
vicinity of $400-600$ K. Solutions to this apparent paradox include the
following possibilities. (1) If the material responsible for the dust is in a
ring of debris, then collisions can both eject pulverized fine particles into
the surrounding space, or expose pristine volatiles lying beneath the surface
of the body which had previously shielded this material from the stellar
radiation. (2) Volatile materials beneath the surface layers, and therefore
largely shielded from the stellar radiation, may nonetheless be slowly heated
to the point where the pressure buildup blows away and ejects the overlying
material. (3) In a related scenario, if orbit migration of a comet is rapid
compared to the timescale for heat to penetrate to the comet center, an exotic
model invoking deep pits can be invoked to preserve volatile-based dust
release at the present time. (4) If the dust-emitting body is indeed orbiting
completely inside a dusty disk, then the equilibrium temperature of the body
given by the Chiang & Goldreich (1997) expression may be sufficiently low that
dust-ladened ices of CO2 and H2O may survive the radiation of the white dwarf.
Lastly, (5) if a body is nearly filling its critical potential lobe, and the
orbit is slightly eccentric, then fragments of it may peel off at periastron
passages.

Table 8: Equilibrium Temperature on Orbit

Dependencea |  $T_{\rm eq}$(9.93 hr) |  $T_{\rm eq}$(11.2 hr)
---|---|---
$T_{\rm BB}^{b}$ | $\simeq T_{\star}(R_{\star}/2d)^{1/2}$ | 385 K | 370 K
$T_{\rm hemi}^{c}$ | $\simeq T_{\star}(R_{\star}/\sqrt{2}d)^{1/2}$ | 458 K | 440 K
$T_{\rm ssp}^{d}$ | $\simeq T_{\star}(R_{\star}/d)^{1/2}$ | 545 K | 525 K
$T_{\rm ssp}$ crater bottomse | 600 K | 575 K
$T_{\rm Rayl}^{f}$ | $\simeq T_{\star}(R_{\star}/2d)^{2/5}$ | 700 K | 677 K
$T_{\rm disk}^{g}$ | $\simeq\xi T_{\star}(R_{\star}/d)^{3/4}$ | 99 K | 93 K

aafootnotetext: $R_{\star}$ and $T_{\star}$ are the white dwarf radius and
temperature, and $d$ is the orbital distance of the heated body from the white
dwarf.

bbfootnotetext: $T_{\rm BB}$ is the blackbody temperature of a body that is
small enough to be heated to a uniform temperature, but large compared to the
wavelength of the incident radiation (see Xu et al., 2018).

ccfootnotetext: Temperature of a body heated over a hemisphere to a uniform
temperature over that hemisphere.

ddfootnotetext: Temperature of a large body at its substellar point.

eefootnotetext: Crater bottoms near the substellar point may be heated to
temperatures $\sim$10% hotter than $T_{\rm ssp}$ (see Hansen 1977; Keihm, S.,
personal communication).

fffootnotetext: Temperature of a particle much smaller than the wavelength of
the incident radiation (Rayleigh limit; see Xu et al., 2018).

ggfootnotetext: Temperature in the interior of a gas-dust disk (Chiang &
Goldreich, 1997) where the leading coefficient is $\xi=(2/3\pi)^{1/4}$. This
entry below the horizontal line is the only one where shielding of the
orbiting body by a disk is invoked.

Table 9: Asteroids vs. Comets

Parameter | Asteroida | Cometa
---|---|---
Mass | up to $10^{24}$ g | up to $3\times 10^{17}$ g
Radius | up to 500 km | up to 30 km
Density |  $1.0-3.5\,$g cm-3 |  $0.3-0.6\,$g cm-3
Composition | rocky | gas ices
| minerals | dust
Lifetime $@$ $10^{9}$ g/s |  $0-10$ Myr | 10 days $-$ 10 yr
Illustrative material | Fe (1546 K) | H2O (252 K)
(Sublimation temp.)b: | SiO (1586 K) | CO2 (194 K)
| Fe2SiO4 (1585 K) | CO (70 K)
| MgSiO3 (1822 K) | …
| Mg2SiO4 (1904 K) | …
| SiO2 (1977 K) | …
| Al2O3 (1969 K) | …
| SiC (2098 K) | …
| C (2585 K) | …

aafootnotetext: Generic properties of asteroids and comets taken from Carry
(2012).

bbfootnotetext: Minerals and characteristic sublimation temperatures taken
from van Lieshout et al. (2014).

As it regards such close-in orbiting bodies, the original Roche limit (Roche,
1849) can be recast as an expression for the minimum allowed orbital period as
a function of the mean density of the orbiting body

$P_{\rm min}\simeq\sqrt{\frac{3\pi(2.44)^{3}}{G\bar{\rho}}}\simeq
12.6\left(\frac{{\rm g~{}cm}^{-3}}{\bar{\rho}}\right)^{1/2}~{}{\rm hr}$ (1)

(see e.g. Rappaport et al., 2013, 2021). A period of 9.93 hrs therefore
requires a minimum mean density of only 1.6 $\mathrm{g\,cm}^{-3}$. By
contrast, in the case of WD 1145+017 with orbital periods near 4.5 hr, the
minimum mean density required is 7.8 $\mathrm{g\,cm}^{-3}$. This all assumes
that the bodies are not held together by material forces. In the case of
rubble piles (Veras et al. 2014, 2017), the individual chunks of material are
relatively free to peel off the host body when the size of the body exceeds
its critical potential lobe. Equation (1) is plotted in Figure 10 to help
visualize where ZTF J0328$-$1219 and WD 1145+017 fall with respect to the
curve.

Figure 10: Minimum required mean density of the orbiting body as a function of
its orbital period—described by Equation (1). The typical density ranges for
comets ($0.3{-}0.6\;\mathrm{g\,cm^{-3}}$), asteroids
($1.0{-}3.5\;\mathrm{g\,cm^{-3}}$), and terrestrial planets
($4.0{-}5.5\;\mathrm{g\,cm^{-3}}$) are indicated on the left-hand side. The
dot-dashed and dotted curves are approximate modifications to Equation (1)
multiplied by a factor of $(1-e)^{-3/2}$ for eccentricities of 0.1 and 0.2,
respectively, which illustrate how higher eccentricities at fixed periods push
the required density higher. Vertical red lines mark the two periods
identified in ZTF J0328$-$1219 while the vertical blue line is a
characteristic period of WD 1145+017.

In terms of the material released by Roche breakup, we envision that smaller
bodies might be produced in the following way. The orbiting bodies are likely
naturally in slightly eccentric orbits. In that case the size of the critical
potential lobe expands and contracts by a fractional amount $\pm e$ around the
orbit, where $e$ is its orbital eccentricity. In that case, the fractional
volume of the rubble pile that is exposed to ejection of material at
periastron passage is $\lesssim$ $3e$. Thus, for example, if $e=10^{-4}$, then
some $\lesssim 0.0003$ of the mass of the rubble pile is susceptible to
breaking off each orbit. In that case, the rubble pile itself could have a
lifetime of only a few thousand orbits, or only a few years.

Another promising mechanism for producing dust even in quite cool environments
might be a set of enhanced cascading collisions. This scenario starts with a
catastrophic collision between two substantial bodies which, after a number of
subsequent orbits, produces an asymmetric debris disk with a region of high
density that remains fixed near the original collision point (Jackson et al.,
2014). Subsequent collisions within this high density region may be capable of
producing fresh clouds of dust that last several orbits before spreading out
in azimuth due to differential orbital velocities. While collisions provide a
viable mechanism for producing dust, however, they have thus far only been
examined theoretically in the contexts of debris disk morphology (Jackson et
al., 2014) and dust and gas disk production and accretion onto white dwarfs
(Kenyon & Bromley, 2017a, b), and do not attempt to produce semi-persistent
dust clouds that could yield the types of transits observed in WD 1145+017 or
ZTF J0328$-$1219.

####  4.3.1 Eccentric Debris Orbits

Here we briefly consider some of the ramifications if the dust-emitting bodies
are in substantially eccentric orbits. For a fixed observed orbital period,
the semi-major axis of the orbit is fixed, regardless of the eccentricity. The
distance of closest approach, at which point the body is most vulnerable to
tidal breakup, is $a(1-e)$, where $a$ is the semi-major axis of the orbit.

One might guess that for eccentric orbits the equilibrium temperatures of a
body would be higher than for a circular orbit of the same semi-major axis,
due to the fact that the orbiting body spends time closer to the irradiating
source. However, while the time-averaged stellar flux experienced by a body
throughout an orbit increases with increasing eccentricity, the average
distance also increases. The net result is an overall decrease in the time-
averaged equilibrium temperature with increasing eccentricity (Méndez &
Rivera-Valentín, 2017). As shown by Méndez & Rivera-Valentín (2017), for
eccentricities of $e<0.5$, and assuming constant albedo throughout the orbit,
the equilibrium temperature drops by only about 1% or less compared to a
circular orbit, and is lowered by no more than 10% even out to eccentricities
approaching unity. Thus, the equilibrium temperatures presented in Table 8 are
upper limits in the case of circular orbits, but the effect of eccentric
orbits on equilibrium temperatures is relatively small.

Another effect of eccentric orbits is that they change the $P_{\rm
min}(\bar{\rho})$ relation give by Equation (1). For an eccentric orbit we can
roughly estimate999This estimate is based on the assumption that the linear
size of the critical potential surface at periastron is proportional to the
distance of closest approach. that the minimum period in Equation (1) is
increased by a factor of $(1-e)^{-3/2}$; conversely, this increases the
minimum mean density required to attain an orbit as short as is observed. In
particular, for $e=0.2$, the mean density for an object in a 9.93-hr orbit
would become 3.1 g cm-3 as opposed to 1.6 g cm-3 for a circular orbit. The
effect of eccentricity on the $P_{\rm min}(\bar{\rho})$ relation is
illustrated graphically in Figure 10.

A question may also arise about the stability of two nearby orbits where an
eccentricity of one or both objects will lead to crossing orbits. To address
this issue we have carried out a few numerical simulations with up to a half
dozen low-mass bodies orbiting a white dwarf. Two of the orbits replicated
those of the 9.93 and 11.2 hr orbits observed in ZTF J0328$-$1219\. We tested
a set of orbits where all bodies had different semi-major axes, but one body
had an eccentricity of 0.2, leading to crossing orbits. We also tested a set
of orbits where two bodies have the same semi-major axis, but one is
eccentric. For each simulation, all orbiting bodies were given the same mass,
and we ran simulations with masses of $10^{20}$, $10^{23}$, $10^{26}$, and
$10^{28}$ g. The first two of these, where the orbiting bodies have masses
comparable to Solar System asteroids, led to stable orbits for more than 3
years (i.e., longer than ZTF J0328$-$1219 has been observed), while the third
of these, with masses comparable to the Moon, shows signs of perturbed orbits.
For the last case, with masses roughly equal to Earth, the orbits were grossly
unstable.

###  4.4 Dust Required to Produce the Photometric Variations

The minimum amount of dust that would be required to produce the photometric
fluctuations that we see can be estimated as follows. The minimum total
effective area that is blocked around the orbit by dust is given by

$\mathcal{A}\gtrsim 4\pi dR_{\star}f\tau$ (2)

where $d$ is the orbital radius of the dust emitting body, $R_{\star}$ is the
white dwarf radius, $f$ is the fraction of the orbit where photometric
modulation is observed, and $\tau$ is the optical depth of the dust in the
visible band. In the optically thin limit, the total cross section of all the
dust particles, $\sigma_{\rm tot}$, for a given total mass in dust, $M_{\rm
dust}$ is

$\sigma_{\rm tot}\simeq N_{\rm grain}\sigma_{\rm grain}\simeq\frac{3M_{\rm
dust}\sigma_{\rm grain}}{4\pi\rho_{d}s^{3}}$ (3)

where $N_{\rm grain}$ is the total number of dust grains blocking light from
the WD, $\sigma_{\rm grain}$ is the cross section of a single dust grain,
$\rho_{d}$ is the mean bulk density of the dust particles, and $s$ is the mean
effective size of a dust grain. We can write a general expression for the
minimum required mass in dust by equating Equations (2) and (3). We then find:

$M_{\rm dust}\gtrsim\frac{16\pi dsR_{\star}f\tau\rho_{d}}{3(\sigma_{\rm
grain}/\sigma_{\rm geom})}$ (4)

where the grain cross section is expressed in terms of its geometric cross
section. The most efficient grain scattering cross sections occur near $\sim$1
micron size particles, where $\sigma_{\rm grain}$ is still $\approx\sigma_{\rm
geom}$. Thus, we have

$M_{\rm dust}\gtrsim 16dsR_{\star}\rho_{d}f\tau~{}\simeq
10^{16}\left(\frac{d}{2\,R_{\odot}}\right)\left(\frac{f\tau}{0.025}\right)~{}\mathrm{g}$
(5)

and we have taken $s\simeq 1\,\mu$m and $\rho_{d}\simeq 3$ g/cc. Finally, if
we somewhat arbitrarily assign a lifetime for the dust-grains of about 1 month
(e.g., the timescale for the photometric modulations to undergo substantial
change, then this leads to an inferred grain production rate of $\dot{M}_{\rm
dust}\simeq 5\times 10^{9}$ g s-1. Similar approaches to estimating dust mass
and loss rates have been discussed in, for example, Chen & Jura (2001),
Rappaport et al. (2012), Vanderburg et al. (2015), and numerous other papers.

We have limited information about the optical depths of the ZTF J0328$-$1219
dust clouds. A dip with depth of 10% can be accounted for by either a cloud
larger than the WD disk with an optical depth of 10%, or an optically thick
cloud that obstructs 10% of the WD disk – or some combination of these two
extreme cases. The optically thick case will exhibit dip depths that are the
same at all wavelengths. The optically thin case can exhibit anything between
no depth dependence upon wavelength and a significant decrease of depth with
increasing wavelength. Any observed wavelength dependence requires that the
particle size distribution (PSD) include a significant component that is
smaller than optical wavelengths (in addition to requiring optical depth
$\lesssim 1$). For WD 1145+017 it has been predicted (Xu et al., 2018) that
the PSD should not include such small particles. Such particles are expected
to become hotter than the volatile temperatures for typical materials because
their emissivity is small at the particles’ black body wavelength while they
absorb starlight as well as large particles.

The orbiting debris in ZTF J0328$-$1219 is far enough from its white dwarf
that no such evaporation of small particles occurs, so we can expect to find
that the ZTF J0328$-$1219 dust clouds consist of the full spectrum of the PSD.
The most relevant observations we have for constraining wavelength dependence
of the dip depth is the ZTF set of $g$ and $r$ band measurements (shown in
Figure 3). Using 30-bin phased light curves in each filter, we find that the
ZTF $g$ and $r$ band fluxes, when plotted against each other, have a slope of
$(g-\langle g\rangle)\simeq(0.38\pm 0.10)(r-\langle r\rangle).$ This is
contrary to what would be expected if the dust clouds included a significant
component of small particles ($<$0.3 micron radius). We note, however, that
the ZTF $g$ and $r$ band measurements are not made simultaneously and cover a
broad time baseline, which may impact the observed amplitude correlation.
Simultaneous multi-band observations would be ideal for addressing this
important question.

##  5 Summary & Conclusions

In this work we have presented new and archival observations which confirm ZTF
J0328$-$1219 as the third white dwarf to exhibit recurring transits caused by
orbiting debris clouds. Photometry spanning more than 900 days was taken from
TESS Sectors 4 and 31, as well as ground-based observations from ZTF, McDonald
Observatory, SAAO, HAO, and JBO. The periods underlying the dusty transit
behavior are 9.93 hr (A-period) and 11.2 hr (B-period), with the former
periodicity being dominant. The peaks in a coherent Lomb-Scargle transform of
the photometric data are about three times wider than would have been found
for constant periods and stationary orbital modulation profiles. At this point
we cannot determine whether pure amplitude and shape modulations are
responsible for the broad peaks, or if slightly varying periods or the
presence of multiple closely spaced periods for several transiting dust clouds
are also responsible.

These periods are roughly twice as long as those observed in WD 1145+017, the
first of the white dwarfs found to be transited by multiple dusty bodies
(Vanderburg et al., 2015), but still about two orders-of-magnitude shorter
compared to the ${\approx}100$-day period observed in ZTF J0139+5245
(Vanderbosch et al., 2020). The nearly periodic dips in flux seen in WD
1145+017 differ in two somewhat subtle ways from those observed thus far in
ZTF J0328$-$1219\. First, the dips seen in WD 1145+017 tend to be more
localized with substantial portions of the orbital cycle where the flux level
is constant between dips. In ZTF J0328$-$1219 there is variability throughout
essentially the entire periodic cycle. Second, the modulation depths in WD
1145+017 can reach 60% of the flux. Thus far, the modulation amplitudes in ZTF
J0328$-$1219 are about 10% or less.

The existence of two periods argues strongly in favor of multiple dust-
emitting orbiting bodies. In this regard, ZTF J0328$-$1219 is similar to WD
1145+017, which also has multiple periodicities. The two periods in ZTF
J0328$-$1219, plus the changing modulation profiles exhibiting highly non-
sinusoidal variations, argue against the periodicity being due to a white
dwarf rotational modulation. If rotational modulation is involved, it would
more likely be at the 11.2 hr B-period since we have yet to identify any dip-
like features with high-S/N follow-up observations. The B-period detections by
TESS, however, suggest changes to both the period and modulation profile over
time, characteristics which are inconsistent with rotational modulation.

We have used both an LDT spectrum and SED photometry to carefully characterize
the white dwarf atmospheric parameters. We find $M_{\star}=0.731\pm
0.023\,M_{\odot}$, $T_{\mathrm{eff}}=7630\pm 140\,$K, [Ca/He] $=-9.55\pm
0.12$, and [H/He] $<-3.5$. With new spectroscopic observations from MIKE and
SOAR, we have also detected narrow Na D doublet features with velocities
21.4$\pm$1.0 km s-1 blue-shifted with respect to atmospheric features, which
we determine to be of circumstellar origin. These lines are weak relative to
the photospheric Ca H and K lines, so we suspect the impact of circumstellar
absorption on our calculated white dwarf atmospheric parameters to be
insignificant. A weak H$\alpha$ line which we believe to be of atmospheric
origin is also detected, though it is poorly fit by our DZ model.

Finally, we would like to advocate for further long-term monitoring of this
source. Due to its brightness, ZTF J0328$-$1219 is accessible to good amateur
photometry with 16-in telescopes, as demonstrated by our HAO and JBO
photometry. This will help resolve the issue of whether the two periods we
have found are unique, whether they are long-term stable, and whether the dips
have many slightly different periods, as is the case for WD 1145+017.
Additionally, continued high-resolution spectroscopic monitoring will be
useful for monitoring the velocities of the detected circumstellar features.
Velocity variations over time may help constrain the eccentricity of the
orbiting gaseous debris. Lastly, mid-infrared observations with JWST will be
useful to better constrain the presence of a dusty debris disk in this system.

We thank the anonymous referee for a detailed review which led to significant
improvements of this paper. We also thank Rob Robinson (Univ. Texas, Austin)
for helpful comments related to the properties of cataclysmic variables, and
we thank Mathieu Choukroun (Caltech/JPL), Stephen J. Keihm (Caltech/JPL),
Vishnu Reddy (Univ. Arizona) and Richard Binzel (MIT) for helpful discussions
about solar system asteroids. We also thank Erik Dennihy (Gemini South) for an
independent look at both the SED and SOAR spectra, and thank John Kuehne,
Brian Roman, Coyne Gibson, Dave Doss, Anita Cochran, and all other McDonald
observing support staff for making remote observing on the McDonald 2.1-m
telescope possible during COVID-related travel restrictions. Z.P.V.
acknowledges support from the Wootton Center for Astrophysical Plasma
Properties under U.S. Department of Energy cooperative agreement number DE-
NA0003843. S.B. acknowledges support from the Laboratory Directed Research and
Development program of Los Alamos National Laboratory under project number
20190624PRD2. C.M. and B.Z. acknowledge support from NSF grants SPG-1826583
and SPG-1826550. B.L.K. was supported by the M. Hildred Blewett Fellowship of
the American Physical Society. J.J.H and T.M.H. acknowledge support from the
National Science Foundation under Grant No. AST-1908119. This paper employs
data collected by the Transiting Exoplanet Survey Satellite (TESS) mission.
Funding for the TESS mission is provided by the NASA’s Science Mission
Directorate. This research was made possible through the use of the AAVSO
Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences
Fund and NSF AST-1412587. This work has made 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). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. This work has also made use
of observations obtained with the Samuel Oschin 48-inch Telescope at the
Palomar Observatory as part of the Zwicky Transient Facility (ZTF) project.
ZTF is supported by the NSF under Grant No. AST-1440341 and collaborating
institutions. Observations obtained using the 2.1 m Otto Struve Telescope at
the McDonald Observatory operated by The University of Texas at Austin are
included in this work. Observations made using the Goodman Spectrograph
mounted on the Southern Astrophysical Research (SOAR) telescope are also used
in this work. We have also made use of observations obtained using the DeVeny
Spectrograph mounted on the Lowell Discovery Telescope (LDT) at Lowell
Observatory. This paper uses observations made at the South African
Astronomical Observatory (SAAO).

## References

  * 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 et al. (2019) Bell, K. J., Córsico, A. H., Bischoff-Kim, A., et al. 2019, A&A, 632, A42, doi: 10.1051/0004-6361/201936340
  * 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
  * Bernstein et al. (2003) Bernstein, R., Shectman, S. A., Gunnels, S. M., Mochnacki, S., & Athey, A. E. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye & A. F. M. Moorwood, 1694–1704, doi: 10.1117/12.461502
  * Bianchi et al. (2017) Bianchi, L., Shiao, B., & Thilker, D. 2017, ApJS, 230, 24, doi: 10.3847/1538-4365/aa7053
  * Bida et al. (2014) Bida, T. A., Dunham, E. W., Massey, P., & Roe, H. G. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 91472N, doi: 10.1117/12.2056872
  * Blouin (2020) Blouin, S. 2020, MNRAS, 496, 1881, doi: 10.1093/mnras/staa1689
  * Blouin et al. (2020) Blouin, S., Daligault, J., Saumon, D., Bédard, A., & Brassard, P. 2020, A&A, 640, L11, doi: 10.1051/0004-6361/202038879
  * Blouin et al. (2018) Blouin, S., Dufour, P., & Allard, N. F. 2018, ApJ, 863, 184, doi: 10.3847/1538-4357/aad4a9
  * Bradley et al. (2020) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, astropy/photutils: 1.0.0, 1.0.0, Zenodo, doi: 10.5281/zenodo.4044744
  * Breger et al. (1993) Breger, M., Stich, J., Garrido, R., et al. 1993, A&A, 271, 482
  * Burgasser (2014) Burgasser, A. J. 2014, in Astronomical Society of India Conference Series, Vol. 11, Astronomical Society of India Conference Series, 7–16.  https://arxiv.org/abs/1406.4887
  * Carry (2012) Carry, B. 2012, Planet. Space Sci., 73, 98, doi: 10.1016/j.pss.2012.03.009
  * Casewell et al. (2017) Casewell, S. L., Geier, S., & Lodieu, N. 2017, in Astronomical Society of the Pacific Conference Series, Vol. 509, 20th European White Dwarf Workshop, ed. P. E. Tremblay, B. Gaensicke, & T. Marsh, 133
  * Cauley et al. (2018) Cauley, P. W., Farihi, J., Redfield, S., et al. 2018, ApJ, 852, L22, doi: 10.3847/2041-8213/aaa3d9
  * 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 & Jura (2001) Chen, C. H., & Jura, M. 2001, ApJ, 560, L171, doi: 10.1086/324057
  * Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368, doi: 10.1086/304869
  * Clemens et al. (2004) Clemens, J. C., Crain, J. A., & Anderson, R. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5492, Ground-based Instrumentation for Astronomy, ed. A. F. M. Moorwood & M. Iye, 331–340, doi: 10.1117/12.550069
  * Collins et al. (2017) Collins, K. A., Kielkopf, J. F., Stassun, K. G., & Hessman, F. V. 2017, AJ, 153, 77, doi: 10.3847/1538-3881/153/2/77
  * Coppejans et al. (2013) Coppejans, R., Gulbis, A. A. S., Kotze, M. M., et al. 2013, PASP, 125, 976, doi: 10.1086/672156
  * Coutu et al. (2019) Coutu, S., Dufour, P., Bergeron, P., et al. 2019, ApJ, 885, 74, doi: 10.3847/1538-4357/ab46b9
  * Currie et al. (2014) Currie, M. J., Berry, D. S., Jenness, T., et al. 2014, in Astronomical Society of the Pacific Conference Series, Vol. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay, 391
  * Cutri et al. (2021) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR Online Data Catalog, II/328
  * Debes et al. (2012) Debes, J. H., Kilic, M., Faedi, F., et al. 2012, ApJ, 754, 59, doi: 10.1088/0004-637X/754/1/59
  * Debes & Sigurdsson (2002) Debes, J. H., & Sigurdsson, S. 2002, ApJ, 572, 556, doi: 10.1086/340291
  * Dennihy et al. (2020a) Dennihy, E., Farihi, J., Fusillo, N. P. G., & Debes, J. H. 2020a, ApJ, 891, 97, doi: 10.3847/1538-4357/ab7249
  * Dennihy et al. (2020b) Dennihy, E., Xu, S., Lai, S., et al. 2020b, ApJ, 905, 5, doi: 10.3847/1538-4357/abc339
  * Drake et al. (2009) Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2009, ApJ, 696, 870, doi: 10.1088/0004-637X/696/1/870
  * Dufour et al. (2017) Dufour, P., Blouin, S., Coutu, S., et al. 2017, in Astronomical Society of the Pacific Conference Series, Vol. 509, 20th European White Dwarf Workshop, ed. P. E. Tremblay, B. Gaensicke, & T. Marsh, 3.  https://arxiv.org/abs/1610.00986
  * Dufour et al. (2007) Dufour, P., Bergeron, P., Liebert, J., et al. 2007, ApJ, 663, 1291, doi: 10.1086/518468
  * Dupuy & Liu (2012) Dupuy, T. J., & Liu, M. C. 2012, ApJS, 201, 19, doi: 10.1088/0067-0049/201/2/19
  * Farihi (2016) Farihi, J. 2016, New A Rev., 71, 9, doi: 10.1016/j.newar.2016.03.001
  * Fontaine et al. (2001) Fontaine, G., Brassard, P., & Bergeron, P. 2001, PASP, 113, 409, doi: 10.1086/319535
  * Fortin-Archambault et al. (2020) Fortin-Archambault, M., Dufour, P., & Xu, S. 2020, ApJ, 888, 47, doi: 10.3847/1538-4357/ab585a
  * Fuchs (2017) Fuchs, J. T. 2017, PhD thesis, The University of North Carolina at Chapel Hill
  * Gaia Collaboration et al. (2020) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2020, arXiv e-prints, arXiv:2012.01533.  https://arxiv.org/abs/2012.01533
  * 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
  * Gänsicke et al. (2012) Gänsicke, B. T., Koester, D., Farihi, J., et al. 2012, MNRAS, 424, 333, doi: 10.1111/j.1365-2966.2012.21201.x
  * Gänsicke et al. (2006) Gänsicke, B. T., Marsh, T. R., Southworth, J., & Rebassa-Mansergas, A. 2006, Science, 314, 1908, doi: 10.1126/science.1135033
  * Gänsicke et al. (2020) Gänsicke, B. T., Rodríguez-Gil, P., Gentile Fusillo, N. P., et al. 2020, MNRAS, 499, 2564, doi: 10.1093/mnras/staa2969
  * Gänsicke et al. (2016) Gänsicke, B. T., Aungwerojwit, A., Marsh, T. R., et al. 2016, ApJ, 818, L7, doi: 10.3847/2041-8205/818/1/L7
  * Gary et al. (2017) Gary, B. L., Rappaport, S., Kaye, T. G., Alonso, R., & Hambschs, F. J. 2017, MNRAS, 465, 3267, doi: 10.1093/mnras/stw2921
  * Gentile Fusillo et al. (2019) Gentile Fusillo, N. P., Tremblay, P.-E., Gänsicke, B. T., et al. 2019, MNRAS, 482, 4570, doi: 10.1093/mnras/sty3016
  * Gentile Fusillo et al. (2020) Gentile Fusillo, N. P., Manser, C. J., Gänsicke, B. T., et al. 2020, arXiv e-prints, arXiv:2010.13807.  https://arxiv.org/abs/2010.13807
  * Girven et al. (2012) Girven, J., Brinkworth, C. S., Farihi, J., et al. 2012, ApJ, 749, 154, doi: 10.1088/0004-637X/749/2/154
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93, doi: 10.3847/1538-4357/ab5362
  * Guidry et al. (2021) Guidry, J. A., Vanderbosch, Z. P., Hermes, J. J., et al. 2021, ApJ, 912, 125, doi: 10.3847/1538-4357/abee68
  * Hansen (1977) Hansen, O. L. 1977, Icarus, 31, 456, doi: 10.1016/0019-1035(77)90149-X
  * Hermes et al. (2017a) Hermes, J. J., Kawaler, S. D., Bischoff-Kim, A., et al. 2017a, ApJ, 835, 277, doi: 10.3847/1538-4357/835/2/277
  * Hermes et al. (2017b) Hermes, J. J., Gänsicke, B. T., Kawaler, S. D., et al. 2017b, The Astrophysical Journal Supplement Series, 232, 23, doi: 10.3847/1538-4365/aa8bb5
  * Holberg & Bergeron (2006) Holberg, J. B., & Bergeron, P. 2006, AJ, 132, 1221, doi: 10.1086/505938
  * Hollands et al. (2017) Hollands, M. A., Koester, D., Alekseev, V., Herbert, E. L., & Gänsicke, B. T. 2017, MNRAS, 467, 4970, doi: 10.1093/mnras/stx250
  * Horne (1986) Horne, K. 1986, PASP, 98, 609, doi: 10.1086/131801
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jackson et al. (2014) Jackson, A. P., Wyatt, M. C., Bonsor, A., & Veras, D. 2014, MNRAS, 440, 3757, doi: 10.1093/mnras/stu476
  * Jenkins (2015) Jenkins, J. M. 2015, in AAS/Division for Extreme Solar Systems Abstracts, Vol. 3, AAS/Division for Extreme Solar Systems Abstracts, 106.05
  * Jenkins et al. (2016) Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in Proc. SPIE, Vol. 9913, Software and Cyberinfrastructure for Astronomy IV, 99133E, doi: 10.1117/12.2233418
  * Jura (2003) Jura, M. 2003, ApJ, 584, L91, doi: 10.1086/374036
  * Kanaan et al. (2002) Kanaan, A., Kepler, S. O., & Winget, D. E. 2002, A&A, 389, 896, doi: 10.1051/0004-6361:20020485
  * Kelson (2003) Kelson, D. D. 2003, PASP, 115, 688, doi: 10.1086/375502
  * Kelson et al. (2000) Kelson, D. D., Illingworth, G. D., van Dokkum, P. G., & Franx, M. 2000, ApJ, 531, 137, doi: 10.1086/308460
  * Kenyon & Bromley (2017a) Kenyon, S. J., & Bromley, B. C. 2017a, ApJ, 844, 116, doi: 10.3847/1538-4357/aa7b85
  * Kenyon & Bromley (2017b) —. 2017b, ApJ, 850, 50, doi: 10.3847/1538-4357/aa9570
  * Kilic et al. (2015) Kilic, M., Gianninas, A., Bell, K. J., et al. 2015, ApJ, 814, L31, doi: 10.1088/2041-8205/814/2/L31
  * Klein et al. (2011) Klein, B., Jura, M., Koester, D., & Zuckerman, B. 2011, ApJ, 741, 64, doi: 10.1088/0004-637X/741/1/64
  * Koester (2009) Koester, D. 2009, A&A, 498, 517, doi: 10.1051/0004-6361/200811468
  * Koester et al. (2014) Koester, D., Gänsicke, B. T., & Farihi, J. 2014, A&A, 566, A34, doi: 10.1051/0004-6361/201423691
  * Kuschnig et al. (1997) Kuschnig, R., Weiss, W. W., Gruber, R., Bely, P. Y., & Jenkner, H. 1997, A&A, 328, 544
  * Lenz & Breger (2005) Lenz, P., & Breger, M. 2005, Communications in Asteroseismology, 146, 53, doi: 10.1553/cia146s53
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447, doi: 10.1007/BF00648343
  * Longstaff et al. (2017) Longstaff, E. S., Casewell, S. L., Wynn, G. A., Maxted, P. F. L., & Helling, C. 2017, MNRAS, 471, 1728, doi: 10.1093/mnras/stx1786
  * Longstaff et al. (2019) Longstaff, E. S., Casewell, S. L., Wynn, G. A., et al. 2019, MNRAS, 484, 2566, doi: 10.1093/mnras/stz127
  * Magnier et al. (2013) Magnier, E. A., Schlafly, E., Finkbeiner, D., et al. 2013, ApJS, 205, 20, doi: 10.1088/0067-0049/205/2/20
  * Maoz et al. (2015) Maoz, D., Mazeh, T., & McQuillan, A. 2015, MNRAS, 447, 1749, doi: 10.1093/mnras/stu2577
  * Marsh (1989) Marsh, T. R. 1989, PASP, 101, 1032, doi: 10.1086/132570
  * Martin et al. (2005) Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJ, 619, L1, doi: 10.1086/426387
  * Masci et al. (2019) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003, doi: 10.1088/1538-3873/aae8ac
  * McMahon et al. (2013) McMahon, R. G., Banerji, M., Gonzalez, E., et al. 2013, The Messenger, 154, 35
  * Melis et al. (2020) Melis, C., Klein, B., Doyle, A. E., et al. 2020, ApJ, 905, 56, doi: 10.3847/1538-4357/abbdfa
  * Méndez & Rivera-Valentín (2017) Méndez, A., & Rivera-Valentín, E. G. 2017, ApJ, 837, L1, doi: 10.3847/2041-8213/aa5f13
  * Montgomery & Odonoghue (1999) Montgomery, M. H., & Odonoghue, D. 1999, Delta Scuti Star Newsletter, 13, 28
  * Morrissey et al. (2005) Morrissey, P., Schiminovich, D., Barlow, T. A., et al. 2005, ApJ, 619, L7, doi: 10.1086/424734
  * Mullally et al. (2007) Mullally, F., Kilic, M., Reach, W. T., et al. 2007, ApJS, 171, 206, doi: 10.1086/511858
  * Nather et al. (1990) Nather, R. E., Winget, D. E., Clemens, J. C., Hansen, C. J., & Hine, B. P. 1990, ApJ, 361, 309, doi: 10.1086/169196
  * Newville et al. (2014) Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, 0.8.0, Zenodo, doi: 10.5281/zenodo.11813
  * Osterbrock et al. (1996) Osterbrock, D. E., Fulbright, J. P., Martel, A. R., et al. 1996, PASP, 108, 277, doi: 10.1086/133722
  * Plavchan et al. (2008) Plavchan, P., Jura, M., Kirkpatrick, J. D., Cutri, R. M., & Gallagher, S. C. 2008, ApJS, 175, 191, doi: 10.1086/523644
  * Rappaport et al. (2016) Rappaport, S., Gary, B. L., Kaye, T., et al. 2016, MNRAS, 458, 3904, doi: 10.1093/mnras/stw612
  * Rappaport et al. (2018) Rappaport, S., Gary, B. L., Vanderburg, A., et al. 2018, MNRAS, 474, 933, doi: 10.1093/mnras/stx2663
  * Rappaport et al. (2021) Rappaport, S., Sanchis-Ojeda, R., Rogers, L. A., Levine, A., & J, W. 2021, submitted to ApJ
  * Rappaport et al. (2013) Rappaport, S., Sanchis-Ojeda, R., Rogers, L. A., Levine, A., & Winn, J. N. 2013, ApJ, 773, L15, doi: 10.1088/2041-8205/773/1/L15
  * Rappaport et al. (2012) Rappaport, S., Levine, A., Chiang, E., et al. 2012, ApJ, 752, 1, doi: 10.1088/0004-637X/752/1/1
  * Redfield et al. (2017) Redfield, S., Farihi, J., Cauley, P. W., et al. 2017, ApJ, 839, 42, doi: 10.3847/1538-4357/aa68a0
  * Reding et al. (2020) Reding, J. S., Hermes, J. J., Vanderbosch, Z., et al. 2020, ApJ, 894, 19, doi: 10.3847/1538-4357/ab8239
  * Ricker et al. (2014) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 1 , doi: 10.1117/1.JATIS.1.1.014003
  * Roche (1849) Roche, E. 1849, Académie des Sciences de Montpellier: Mémoires de la Section des Sciences, 1, 243
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835, doi: 10.1086/160554
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schlafly et al. (2019) Schlafly, E. F., Meisner, A. M., & Green, G. M. 2019, ApJS, 240, 30, doi: 10.3847/1538-4365/aafbea
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Sion et al. (1983) Sion, E. M., Greenstein, J. L., Landstreet, J. D., et al. 1983, ApJ, 269, 253, doi: 10.1086/161036
  * Smith et al. (2012) Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000, doi: 10.1086/667697
  * Steele et al. (2020) Steele, A., Debes, J., Xu, S., Yeh, S., & Dufour, P. 2020, arXiv e-prints, arXiv:2010.12667.  https://arxiv.org/abs/2010.12667
  * Stumpe et al. (2014) Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100, doi: 10.1086/674989
  * Szkody et al. (2009) Szkody, P., Anderson, S. F., Hayden, M., et al. 2009, AJ, 137, 4011, doi: 10.1088/0004-6256/137/4/4011
  * Tremblay et al. (2010) Tremblay, P. E., Bergeron, P., Kalirai, J. S., & Gianninas, A. 2010, ApJ, 712, 1345, doi: 10.1088/0004-637X/712/2/1345
  * van Lieshout et al. (2014) van Lieshout, R., Min, M., & Dominik, C. 2014, A&A, 572, A76, doi: 10.1051/0004-6361/201424876
  * Vanderbosch et al. (2020) Vanderbosch, Z., Hermes, J. J., Dennihy, E., et al. 2020, ApJ, 897, 171, doi: 10.3847/1538-4357/ab9649
  * Vanderburg et al. (2015) Vanderburg, A., Johnson, J. A., Rappaport, S., et al. 2015, Nature, 526, 546, doi: 10.1038/nature15527
  * Vanderburg et al. (2019) Vanderburg, A., Huang, C. X., Rodriguez, J. E., et al. 2019, ApJ, 881, L19, doi: 10.3847/2041-8213/ab322d
  * VanderPlas (2018) VanderPlas, J. T. 2018, ApJS, 236, 16, doi: 10.3847/1538-4365/aab766
  * Veras et al. (2017) Veras, D., Carter, P. J., Leinhardt, Z. M., & Gänsicke, B. T. 2017, MNRAS, 465, 1008, doi: 10.1093/mnras/stw2748
  * Veras & Heng (2020) Veras, D., & Heng, K. 2020, MNRAS, 496, 2292, doi: 10.1093/mnras/staa1632
  * Veras et al. (2014) Veras, D., Leinhardt, Z. M., Bonsor, A., & Gänsicke, B. T. 2014, MNRAS, 445, 2244, doi: 10.1093/mnras/stu1871
  * Veras et al. (2020) Veras, D., McDonald, C. H., & Makarov, V. V. 2020, MNRAS, 492, 5291, doi: 10.1093/mnras/staa243
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Xu et al. (2019) Xu, S., Dufour, P., Klein, B., et al. 2019, AJ, 158, 242, doi: 10.3847/1538-3881/ab4cee
  * Xu et al. (2016) Xu, S., Jura, M., Dufour, P., & Zuckerman, B. 2016, ApJ, 816, L22, doi: 10.3847/2041-8205/816/2/L22
  * Xu et al. (2020) Xu, S., Lai, S., & Dennihy, E. 2020, ApJ, 902, 127, doi: 10.3847/1538-4357/abb3fc
  * Xu et al. (2018) Xu, S., Rappaport, S., van Lieshout, R., et al. 2018, MNRAS, 474, 4795, doi: 10.1093/mnras/stx3023
  * Zuckerman & Becklin (1987) Zuckerman, B., & Becklin, E. E. 1987, Nature, 330, 138, doi: 10.1038/330138a0
  * Zuckerman et al. (2007) Zuckerman, B., Koester, D., Melis, C., Hansen, B. M., & Jura, M. 2007, ApJ, 671, 872, doi: 10.1086/522223
  * Zuckerman et al. (2003) Zuckerman, B., Koester, D., Reid, I. N., & Hünsch, M. 2003, ApJ, 596, 477, doi: 10.1086/377492
  * Zuckerman et al. (2010) Zuckerman, B., Melis, C., Klein, B., Koester, D., & Jura, M. 2010, ApJ, 722, 725, doi: 10.1088/0004-637X/722/1/725