# Dust and Volatiles in the Disintegrating Comet C/2019 Y4 (ATLAS)

Ruining Zhao CAS Key Laboratory of Optical Astronomy, National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China;
<EMAIL_ADDRESS><EMAIL_ADDRESS>School of Astronomy and Space Sciences,
University of Chinese Academy of Sciences, Beijing 100049, China  Aigen Li
Department of Physics and Astronomy, University of Missouri, Columbia, MO
65211, USA<EMAIL_ADDRESS>Bin Yang Instituto de Estudios Astrofísicos,
Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Santiago,
Chile<EMAIL_ADDRESS>Liang Wang Nanjing Institute of Astronomical
Optics & Technology, Chinese Academy of Sciences, Nanjing 210042, China  CAS
Key Laboratory of Astronomical Optics & Technology, Nanjing Institute of
Astronomical Optics & Technology, Chinese Academy of Sciences, Nanjing 210042,
China  Huijuan Wang CAS Key Laboratory of Optical Astronomy, National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China
School of Astronomy and Space Sciences, University of Chinese Academy of
Sciences, Beijing 100049, China  Yu-Juan Liu  CAS Key Laboratory of Optical
Astronomy, National Astronomical Observatories, Chinese Academy of Sciences,
Beijing 100101, China  Jifeng Liu CAS Key Laboratory of Optical Astronomy,
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, China  New Cornerstone Science Laboratory, National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100012, China  School of
Astronomy and Space Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China  Institute for Frontiers in Astronomy and Astrophysics,
Beijing Normal University, Beijing, 102206, China

###### Abstract

C/2019 Y4 (ATLAS) is an Oort cloud comet with an orbital period of
$\sim\,$5895$\,{\rm yr}$. Starting in March 2020, its nucleus underwent
disintegration. In order to investigate the gas and dust properties of C/2019
Y4 (ATLAS) during its disintegration, we obtained long-slit spectra at
3600–8700$\,{\rm\AA}$ and $BVRI$ multi-band images with the Xinglong
2.16-Meter Telescope in April 2020. Our observations revealed that C/2019 Y4
(ATLAS) exhibited strong emission bands of CN, C2, C3, and NH2 which are
superimposed on a dust scattering continuum, typical of cometary spectra in
the optical. The production rates of CN, C2, and C3 derived using the Haser
model and the corresponding C2/CN and C3/CN ratios suggest that C/2019 Y4
(ATLAS) is a “typical” Oort cloud comet under the A’Hearn classification,
although it appears less dusty as revealed by the $Af\rho$ quantities. Its
dust-scattering reflectivity is slightly red, with a gradient of $\sim\,$5%
per $10^{3}\,{\rm\AA}$. We model the reflectivity gradient in terms of porous
dust and find that the red color is accounted for by porous dust.

Long period comets (933); Comets (280); Small solar system bodies (1469)

††journal: ApJ††facilities: Beijing:2.16m (BFOSC)††software: astropy (Astropy
Collaboration et al., 2013), ccdproc (Bardyn et al., 2017), photutils (Bradley
et al., 2019), astroquery (Ginsburg et al., 2019), specutils (Astropy-
Specutils Development Team, 2019), sbpy (Mommert et al., 2019), lmfit
(Newville et al., 2014), emcee (Foreman-Mackey et al., 2013),

##  1 Introduction

Cometary nuclei are the most primitive objects in the solar system and have
preserved pristine materials from the presolar molecular cloud and from the
early stages of the protosolar nebula. Therefore, the origin of cometary
nuclei is closely linked to the origin of the solar system. The chemical and
physical properties of dust and volatile materials in cometary nuclei and
comae shed light on the nature of the primordial interstellar materials
present during the solar system formation, and provide clues to the processes
of incorporation of these materials into cometary nuclei, as well as the
chemical and thermodynamic conditions in the outer solar nebula in which
comets formed.

In the absence of a direct analysis of cometary nucleus materials, one often
resorts to cometary comae formed by volatiles outgassed from the nucleus and
dust particles dragged out by the expanding gas. Once formed and stored in the
Oort cloud or the Kuiper belt, however, comets are exposed for $\sim\,$4.5
billion years to the flux of Galactic cosmic rays and may develop a
substantial “crust” of non-volatile materials (e.g., see Strazzulla, 1999).
For short-period comets, such a crust could also result from the solar
radiation during their numerous close approaches to the Sun (e.g., see Li &
Greenberg, 1998a). In this scenario, only the sub-surface materials still
remain pristine. Therefore, disintegrating comets provide us an unique
opportunity to gain access to the pristine, sub-surface materials.

Comet C/2019 Y4 (ATLAS) (hereafter, 19Y4), discovered by the Asteroid
Terrestrial-impact Last Alert System (ATLAS) at the Mauna Loa Observatory on
2019 December 28, is a long-period comet. According to the JPL orbital
solution #17, it takes 19Y4 $5895\pm 23\,{\rm yr}$ to orbit the Sun in an
elliptical trajectory of eccentricity $e=0.999$, inclination $i=45.4^{\circ}$,
and perihelion distance $q=0.253\,{\rm
AU}$.111https://ssd.jpl.nasa.gov/tools/sbdb_lookup.html#/?sstr=2019Y4 Its
initial light curve exhibited such a steep brightening that, if continued,
19Y4 would have been visible to naked eyes when approaching its perihelion in
2020 late May.222http://astro.vanbuitenen.nl/comet/2019Y4 However, the steep
brightening terminated in 2020 mid-March, when 19Y4 started to disintegrate,
as revealed by the nongravitational effect and a blue color likely caused by
the release of a large amount of volatiles (Hui & Ye, 2020). The polarization
increased dramatically in 2020 late-March as well, indicating a large increase
in the amount of carbonaceous material (Zubko et al., 2020). On 2020 April 6,
Ye & Zhang (2020) and Steele et al. (2020) successively reported an elongated
nucleus. On 2020 April 13, the nucleus of 19Y4 had already split into at least
five condensations (Sekanina et al., 2020). Optical spectroscopy on 2020 April
14 and 16 revealed emission bands of CN, C2, C3, and NH2. The C2/CN and C3/CN
production-rate ratios suggested a “typical” Oort cloud comet under the
A’Hearn classification (Ivanova et al., 2021). Continuing disintegration was
reported by Ye & Hui (2020), and the evolution of the fragment clusters was
investigated in detail through high angular resolution images obtained with
the Hubble Space Telescope (HST) on 2020 April 20 and 23 (Ye et al., 2021).

To explore the properties of the pristine volatiles and dust beneath the
refractory crust, we have performed long-slit spectroscopic and broad-band
photometric observations of 19Y4 during its disintegration. This paper reports
and analyzes the spectroscopic and photometric data and is organized as
follows. The observation and data reduction are described in §2 and the
results are reported in §3. The gas and dust productions are discussed in §4.1
and §4.2, respectively. The dust-scattered light is modeled in §4.3 to infer
the properties of the dust. We summarize our major conclusions in §5.

##  2 Observation and Data Reduction

###  2.1 Long-slit Spectroscopy

We obtained long-slit optical spectra of 19Y4 with the Beijing Faint Object
Spectrograph and Camera (BFOSC) on board the 2.16-Meter Telescope at Xinglong
Station. BFOSC has a $2048\times 2048$ pixel2 CCD installed, subtending a
$9.36\arcmin\times 9.36\arcmin$ field of view (FOV). The pixel scale is thus
$\sim\,$0.274$\arcsec$ and no binning was performed. The slit length of
9.4$\arcmin$ was fixed along the north-south orientation. We note that the
north-south orientation is parallactic only when the target is on the
meridian. So our spectra are somewhat subjected to atmospheric dispersion (see
§3). A slit width of 1.8$\arcsec$ or 2.3$\arcsec$ was chosen, depending on the
seeing at the time of observation. The G4 grism, covering a wavelength range
of 3600–8700$\,{\rm\AA}$, was used. It has a resolving power of $\sim$265 at
[O iii] $\lambda$4959$\,{\rm\AA}$ and $\lambda$5007$\,{\rm\AA}$, if combined
with a slit width of 2.3$\arcsec$ (Fan et al., 2016).

We performed spectroscopic observations of 19Y4 in four nights: 2020 April
6.51, 13.53, 20.53, and 23.58 (UT). Although the disintegration of the nucleus
had already started in mid-March 2020 (Hui & Ye, 2020), no fragment was
resolved in our observation on April 6.51 (UT), and the slit was centered on
the brightest pixels of the elongated nucleus. On April 13.53 (UT), the slit
was centered on the fragment 19Y4-A (see Figure 1a). On April 20.53 (UT),
fragments 19Y4-A and 19Y4-B were resolved and they roughly aligned along the
north-south orientation (see Figure 1b), therefore, the slit was over both of
them. On April 23.58 (UT), the slit was centered on the brightest fragment
19Y4-B. Non-sidereal tracking was used. We alternated between short exposures
of 600$\,{\rm s}$ and slit view inspections in order to ensure the guiding
accuracy. The iron-argon arc was used as the wavelength calibrator. In
addition, we have also obtained the long-slit spectra of several flux
standards: Feige 34, He 3, HR 4554, HD 109995, HZ 44 and Feige 98. The
spectroscopy log is shown in Table 1. Also tabulated in Table 1 are the
kinematic parameters of 19Y4, including its phase angle ($\phi$), heliocentric
($r_{\rm h}$), and geocentric distances ($\Delta$) as well as their temporal
changing rates $\dot{r_{\rm h}}$ and $\dot{\Delta}$, obtained from the JPL
Horizons On-Line Ephemeris System web-
interface.333https://ssd.jpl.nasa.gov/horizons.cgi

Following the general data reduction procedures, we reduced the long-slit
spectroscopic data with Python. First, we took bias subtraction, corrected for
flat fields, and removed cosmic rays. Then, we corrected for the geometric
curvature by fitting the curved arc lines with the B-spline functions.
Finally, we combined multiple short exposures in each night, and calibrated
their wavelengths and fluxes to derive fully reduced 2-D spectra.

###  2.2 Broadband Photometry

BFOSC was also used to perform broadband photometry over 19Y4 on 2020 April
23.55 (UT). The CCD configuration is the same as that used in the long-slit
spectroscopy (see §2.1). We made three exposures of 300$\,{\rm s}$ each, with
the BVRI filters in the Johnson-Cousins system. Again, non-sidereal tracking
was used. The photometry log is also shown in Table 1.

We developed a Python pipeline to reduce the photometric data. First, we
corrected all the science frames for bias, flat field and cosmic ray pixels.
Then, we smoothed the corrected frames to detect and center the peaks of the
trailed field stars, which were flattened under non-sidereal tracking. We
performed photometry at each centroid on the corrected frames with rectangular
aperture and annulus to obtain instrumental magnitudes. Next, we used the
centroids of the star trails to derive plate solutions via Astrometry.net
(Lang et al., 2010), and the instrumental magnitudes to derive flux zeropoints
by cross-match with the Sloan Digital Sky Survey Data Release 16 catalog
(Ahumada et al., 2022). The BVRI fluxes and errors of the field stars were
constructed from the Sloan gri fluxes and errors according to the
transformation relations given in Chonis & Gaskell (2008). The flux calibrated
frames were finally combined according to their filters.

##  3 Results

We adopt two apertures, different in width, to extract 1-D spectra: the narrow
one, with a single aperture, has a width of $\sim\,$27.5$\arcsec$
(corresponding to a physical scale of $\sim\,$2.4 $\times$ $10^{4}\,{\rm
km}$); the broad one, divided into 40 sub-apertures, has a width
$\sim\,$165$\arcsec$ in total (corresponding to a physical scale of
$\sim\,$1.2 $\times$ $10^{5}\,{\rm km}$; see Figure 1). The narrow aperture is
chosen to maximize the signal-to-noise ratio in the dust continuum bands, and
the extracted spectra are used to determine the reflectivity gradients (see
below). We show in Figure 2 the narrow aperture spectra of 19Y4. The spectra
are composed of emission bands from gases and sunlight scattered by dust
grains. With the help of an emission line catalog (Brown et al., 1996), we
identify the $\Delta\nu=0$ bands of CN, C2, and C3, the $\Delta\nu=+1$ band of
C2, and a series of bands of NH2. The sub-aperture spectra are used to both
estimate the influence of atmospheric dispersion on the reflectivity gradients
(see below) and model the gas and dust production (see §4.1 and §4.2).

Figure 3 shows the BVRI images of 19Y4 obtained with BFOSC on 2020 April 23.55
(UT). As the band goes from B to I, the coma of 19Y4 shows a less-extended
morphology, indicating the broadband colors vary from inner to outer coma. To
characterize the color variation, we derive $\left(B-V\right)$ and
$\left(V-R\right)$ profiles of the coma by performing photometry with a slit-
like aperture along the north-south orientation. The aperture is divided into
a series of sub-apertures, $4.1\arcsec$ in height and $2.3\arcsec$ in width
(see Figure 3), which are the same as that used to extract sub-aperture
spectra. The color profiles are shown in Figure 4. In general, the
$\left(B-V\right)$ profile is slightly redder and the $\left(V-R\right)$
profile is significantly bluer than the solar colors of
$\left(B-V\right)_{\odot}=0.653\pm 0.005$ and
$\left(V-R\right)_{\odot}=0.352\pm 0.007$ (Ramírez et al., 2012). More
specifically, the $\left(B-V\right)$ profile has a local minimum of 0.68 near
the nucleus and becomes redder outwards. After reaching $\sim\,$0.9 at
$\sim\,$30$\arcsec$ from the nucleus, it starts to turn blue. The
$\left(V-R\right)$ profile has a maximum of $-0.07$ near the nucleus and
becomes bluer outwards. As the broad BVR bands cover multiple gas emission
lines and dust continuum, such color variation should not be simply attributed
to the intrinsic changes in the nature of gas or dust. Instead, it may merely
reflects the spread of gases and dust in the coma. To be more quantitative, we
estimate from the spectra that the gas emission bands contribute $\sim\,$50%
and $\sim\,$25% of the VR fluxes, respectively. This suggests that the blue
color of $\left(V-R\right)$ does not necessarily indicate a blue dust
continuum of 19Y4.

To characterize the dust color, we define seven dust continuum bands free from
gas emission and telluric absorption: 5225–5255$\,{\rm\AA}$,
5800–5850$\,{\rm\AA}$, 6410–6450$\,{\rm\AA}$, 6800–6850$\,{\rm\AA}$,
7420–7450$\,{\rm\AA}$, 7540–7580$\,{\rm\AA}$, 7795–7821$\,{\rm\AA}$, and
8000–8050$\,{\rm\AA}$, and derive the reflectivity gradient, following Jewitt
& Meech (1986). Observationally, the reflectivity is defined as

$S^{\rm obs}_{\lambda}(r_{\rm h})=\frac{F^{\rm sca}_{\lambda}(r_{\rm
h})}{F^{\sun}_{\lambda}(1\,{\rm AU})}\times\left(\frac{r_{\rm h}}{\rm
AU}\right)^{2}~{}~{},$ (1)

where $F^{\rm sca}_{\lambda}(r_{\rm h})$ is the dust scattered continuum of
19Y4 at a heliocentric distance of $r_{\rm h}$; $F^{\sun}_{\lambda}(1\,{\rm
AU})$ is the solar flux at $r_{\rm h}=1\,{\rm AU}$. Here we use the solar
reference spectrum from CALSPEC (Bohlin et al., 2014). Since the spectra
obtained at different epochs have different dust continuum levels, the
reflectivity is further normalized to $\langle S^{\rm obs}\rangle$, the mean
reflectivity in the observed wavelength range. We fit $S^{\rm
obs}_{\lambda}/\langle S^{\rm obs}\rangle$ with linear functions and determine
the reflectivity gradient, defined as the slope of the best-fit.

As shown in Figure 5, the reflectivity gradients are negative at wavelengths
5000–8200$\,{\rm\AA}$ for all the four epochs, indicating that 19Y4 is
significantly bluer than the Sun for which, by definition, the reflectivity
gradient is 0$\%$ per $10^{3}\,{\rm\AA}$. It is also bluer than typical comets
for which the reflectivity gradient is $\sim\,$5–18% per $10^{3}\,{\rm\AA}$ at
wavelengths 3500–6500$\,{\rm\AA}$ (see Jewitt & Meech, 1986). We note that Hui
& Ye (2020) also found that 19Y4 at one point had a color that was much bluer
than the Sun. However, as we mentioned in §2.1, the spectra of 19Y4 were
affected by atmospheric dispersion and so was the color. Therefore, correction
is needed to recover the true reflectivity gradient.

To estimate the influence of atmospheric dispersion on the reflectivity
gradient, we further derive spectroscopic $\left(B-V\right)$ and
$\left(V-R\right)$ profiles on April 23.58 (UT) by convolving the sub-aperture
spectra with the the transmission functions of the BVR filters, and compare
them in Figure 4 with the photometric color profiles on the same date. If the
spectra were free from atmospheric dispersion, there should not be any
difference between the spectroscopic and photometric color profiles, as the
apertures used to derive the color profiles are the same. Nevertheless, as
shown in Figure 4, the spectroscopic $\left(B-V\right)$ and $\left(V-R\right)$
colors are systematically bluer, indicating the influence of the atmospheric
dispersion on the colors. To the first-order approximation, we assume that the
color deviation due to the atmospheric dispersion varies linearly with the
wavelength, and determine a linear correction function by minimizing the
systematic differences in the color profiles. We note that such a correction
is only applied to the spectrum obtained on 2020 April 23.58 (UT) and, as
shown in Figure 5, the recovered reflectivity gradient is $\sim\,$5% per
$10^{3}\,{\rm\AA}$ at wavelengths 5000–8200$\,{\rm\AA}$. The results are shown
in Figure 5 and will be used to constrain the dust properties in §4.3.

##  4 Discussion

###  4.1 Gas Production

The gas production rate and scale length of a given volatile species are
derived by modeling the “observed” column density profile, which, for species
$j$ at a projected distance $\rho$, is

$N^{\rm obs}_{j}=\frac{4\pi I_{j}(\rho)}{g_{j}}~{}~{},$ (2)

where $I_{j}(\rho)$ is the continuum-subtracted radiance of the species
measured at $\rho$ and integrated over a bandpass, and $g_{j}$ is the
fluorescence efficiency of the species averaged over the band (Langland-Shula
& Smith, 2011). We average $N^{\rm obs}_{j}$ over the two pairs of the sub-
apertures symmetric about the nucleus. In this work, we derive the “observed”
column density profiles of CN, C3 and C2 through their $\Delta\nu=0$ bands.
The bandpasses and fluorescence efficiencies are taken from Langland-Shula &
Smith (2011) except for CN ($\Delta\nu=0$), the fluorescence efficiency of
Schleicher (2010) is adopted for which the Swings effect (Swings, 1941) is
taken into account. We show in Figure 6 the “observed” column density profiles
of CN, C3, and C2.

In the literature, cometary gas production rates, scale lengths and their
dependencies on $r_{\rm h}$ have been carefully investigated in numerous
statistical works (e.g., see A’Hearn et al., 1995; Fink & Hicks, 1996;
Langland-Shula & Smith, 2011). We apply those results as prior knowledge for
Bayesian inference studies. To this end, we assume a log-uniform prior for the
gas production rate $Q$ and normal priors for both primary ($l_{0}$) and
product ($l_{1}$) scale lengths, i.e., $\log
Q\sim\mathcal{U}\left(22,28\right)$,
$l_{0}\sim\mathcal{N}(\hat{l}_{0},\hat{\sigma}^{2}_{l_{0}})$ and
$l_{1}\sim\mathcal{N}(\hat{l}_{1},\hat{\sigma}^{2}_{l_{1}})$, where
$\hat{l}_{0},\hat{\sigma}_{l_{0}}$ and $\hat{l}_{1},\hat{\sigma}_{l_{1}}$ are
taken from Langland-Shula & Smith (2011). Here $\mathcal{U}$ and $\mathcal{N}$
are probability functions. With $\log Q\sim\mathcal{U}\left(22,28\right)$, we
assume that the production rate $Q$ is uniformly distributed between
$10^{22}\,{\rm s}^{-1}$ and $10^{28}\,{\rm s}^{-1}$ in the log space. As we do
not have any knowledge about $Q$, we should assume such a non-informative
prior. $l_{0}$ and $l_{1}$ are scale lengths for both of which we assume
Gaussian priors, as their values have been well investigated in previous
statistical studies. The Gaussian parameters
$\hat{l}_{0},\hat{\sigma}_{l_{0}}$ and $\hat{l}_{1},\hat{\sigma}_{l_{1}}$
(i.e., mean and standard deviation) are taken from Langland-Shula & Smith
(2011). The use of Bayesian inference is justified later in this section. The
corresponding prior probability density function is denoted as
$p\left(Q,l_{0},l_{1}\right)$. From Bayes’ theorem, the posterior density is
proportional to the prior multiplied by a likelihood function, i.e.,

$p\left(Q,l_{0},l_{1}\middle|N^{\rm obs};\sigma_{N^{\rm
obs}},\rho\right)\propto p\left(N^{\rm
obs}\middle|Q,l_{0},l_{1};\sigma_{N^{\rm obs}},\rho\right)\times
p\left(Q,l_{0},l_{1}\right)~{}~{},$ (3)

where $\sigma_{N^{\rm obs}}$ is the uncertainty of $N^{\rm obs}$. If
$\sigma_{N^{\rm obs}}$ is assumed to be Gaussian, the likelihood function can
be defined as

$p\left(N^{\rm obs}\middle|Q,l_{0},l_{1};\sigma_{N^{\rm
obs}},\rho\right)=\prod_{i}\frac{1}{\sqrt{2\pi}\sigma_{N^{\rm
obs}_{i}}}\exp\left[-\frac{\left(N^{\rm mod}_{i}-N^{\rm
obs}_{i}\right)^{2}}{2\sigma^{2}_{N^{\rm obs}_{i}}}\right]~{}~{},$ (4)

where $N^{\rm mod}_{i}$ is the modeled column density in the $i$-th sub-
aperture. With a set of $(Q,l_{0},l_{1})$ given, $N^{\rm mod}_{i}$ can be
derived from eqs. (11)–(13) in Langland-Shula & Smith (2011). The joint
posterior distributions are then sampled using the Markov chain Monte Carlo
(MCMC) method. The $50\%$ (median) and $50\pm 34\%$ percentiles in the
marginalized distributions of $Q$, $l_{0}$ and $l_{1}$ are taken as the best-
fit values and the uncertainties, respectively. Note that, in the fitting,
$N^{\rm obs}_{i}$ with $\rho_{i}<5\times 10^{3}\,{\rm km}$ are excluded,
considering that the gas distribution in the inner coma may not be spherically
symmetric due to the disintegration and also that “holes” may be present in
the profiles of C2 (see Langland-Shula & Smith, 2011).

The best-fits to the “observed” column density profiles are shown in Figure 6,
with the best-fit model parameters and uncertainties listed in Table 2. The
fits to C3 ($\Delta\nu=0$) on 2020 April 6.51 and 23.58 (UT) fail due to the
low signal-to-noise ratio of the observational data, and are thus not shown.
Typically, to simultaneously determine both $l_{0}$ and $l_{1}$, the
“observed” profiles are required to extend further than $l_{1}$ (Fray et al.,
2005). However, this is not the case here. If the widely adopted
$\chi^{2}$-minimizing method is used, several $\left(l_{0},l_{1}\right)$ pairs
could lead to equally good fits (Fray et al., 2005; Langland-Shula & Smith,
2011). Therefore, it is the Bayesian model or, more specifically, the prior
that helps to find unique best-fits.

Since G4 does not cover any strong ${\rm OH}$ bands, we take $Q({\rm CN})$ as
an indicator of the gas production rate of 19Y4. We show in Figure 7a the
dependency of $Q({\rm CN})$ on $r_{\rm h}$. In general, the gas production
rate in 19Y4 is “typical” in the sense that both the production rate and its
dependency on $r_{\rm h}$ agree with that of the 26 comets reported in
Langland-Shula & Smith (2011).

We estimate a mean ratio of $\log\left[Q({\rm C_{2}})/Q({\rm
CN})\right]\approx 0.2$ and $\log\left[Q({\rm C_{3}})/Q({\rm
CN})\right]\approx-1.2$ at $r_{\rm h}=1\,{\rm AU}$, if we assume
$\log\left[Q({\rm C_{2}})/Q({\rm CN})\right]$ and $\log\left[Q({\rm
C_{3}})/Q({\rm CN})\right]$ follow a power-law dependence on $r_{\rm h}$ and
take the power-law slopes of Langland-Shula & Smith (2011). The C2/CN and
C3/CN production rate ratios estimated for 19Y4 suggest a “typical” comet
instead of a “depleted” one under the classification of A’Hearn et al. (1995),
consistent with that found by Ivanova et al. (2021).

###  4.2 Dust Production

Initially introduced in A’Hearn et al. (1984), the so-called $Af\rho$
quantities are often used to describe cometary dust production, where $A$ is
the dust reflectivity (i.e., albedo), $f$ is the dust filling factor and
$\rho$ is the projected linear radius of the FOV at the comet as discussed in
§4.1. The $Af\rho$ quantity within a circular aperture with a projected radius
of $\rho$ is determined directly from observation as

$Af\rho=\frac{(2\Delta)^{2}}{\rho}\left(\frac{r_{\rm h}}{\rm
AU}\right)^{2}\frac{F_{\rm c}}{F_{\sun}(1\,{\rm AU})}\frac{1}{\Phi_{\rm
HM}(\phi)}~{}~{},$ (5)

where $F_{\rm c}$ is the cometary flux in a chosen dust continuum band,
$F_{\sun}(1\,{\rm AU})$ is the solar flux at $r_{\rm h}=1\,{\rm AU}$ in the
same band, and $\Phi_{\rm HM}(\phi)$ is the normalized Halley-Marcus phase
function (Schleicher & Bair, 2011). In most cases, $F_{\rm c}$ is obtained
from narrow-band photometry. But it can also be constructed from long-slit
spectroscopy by considering a geometric correction function

$F_{\rm c}=\sum_{i}G(\rho_{i})F(\rho_{i})~{}~{},$ (6)

where $\rho_{i}$ is the projection distance from the $i$-th sub-aperture to
the comet nucleus, and $G(\rho_{i})$ is the geometric correction function
(Langland-Shula & Smith, 2011). The summation is over all the sub-apertures
within $\rho$ (i.e., $\rho_{i}<\rho$). In this work, we calculate $Af\rho$
from the constructed $F_{\rm c}$ in the $\lambda 4850\,{\rm\AA}$ and $\lambda
5240\,{\rm\AA}$ continuum bands refined in Langland-Shula & Smith (2011) with
an aperture radius of 10000$\,{\rm km}$. The results are tabulated in Table 3.
We should admit that our method of reconstructing $Af\rho$ from long-slit
spectra which only sample a strip along certain part of the coma has a
drawback of not accounting for the coma’s asymmetry, even more so for a
disintegrating comet. However, detailed analysis of dust production during
disintegration requires high resolution, spatially resolved photometry and
detailed modeling, which is beyond the scope of this work. Nevertheless, we
note that our $Af\rho$ values closely agree with that from narrow-band
photometry measured with TRAPPIST (Y. Moulane, private conversation).

We obtain a mean $Af\rho$ for 19Y4 by averaging over that in the $\lambda
4850\,{\rm\AA}$ and $\lambda 5240\,{\rm\AA}$ continuum bands and examine the
dependency of $Af\rho$ on $r_{\rm h}$ in Figure 7b. For comparison, we also
show in Figure 7b the $Af\rho$ quantities of the 26 comets studied in
Langland-Shula & Smith (2011). While 19Y4 appears to be less dusty, it
exhibits a $r_{\rm h}$-dependence more or less resembling that of the
Langland-Shula & Smith (2011) samples for which a power-law with an index
$5.3\pm 0.6$ is suggested, even though most comets in the sample of Langland-
Shula & Smith (2011) did not disintegrate when observed.

###  4.3 Scattered Sunlight

The dust-scattered light observationally determined in §3 can be used to infer
the dust size, compositonal, and structural properties. To reproduce the
normalized reflectivity $S^{\rm obs}_{\lambda}/\langle S^{\rm obs}\rangle$
(see §3 and Figure 8), we adopt the cometary dust model of Li & Greenberg
(1998b), which assumes cometary dust grains to be porous aggregates of small
astronomical silicates, amorphous carbon, and vacuum. We consider two cases
for the porosity (i.e., the fractional volumn of vacuum): $f_{\rm{vac}}=0.9$
and $f_{\rm{vac}}=0.5$, which are roughly the upper and lower limits for
cometary dust grains (see Greenberg & Hage, 1990; Greenberg, 1998). For porous
dust resulting from coagulational growth, one would expect a porosity higher
than $f_{\rm vac}=0.5$ (see Blum & Wurm, 2008). On the other hand,
polarimetric measurements suggest that, for most of the time, amorphous
silicates in the coma of 19Y4 account for a fractional volumn ($f_{\rm sil}$)
of 0.17–0.28 (Zubko et al., 2020). Therefore, we assume a bulk volume mixing
ratio of $f_{\rm carb}/f_{\rm sil}=3$ for our dust model. The refractive
indices of amorphous silicates and amorphous carbon are taken, respectively,
from Draine & Lee (1984) and Rouleau & Martin (1991).

We use Mie theory combined with the Bruggman effective medium theory (Bohren &
Huffman, 1983) to calculate $C_{\rm sca}(\lambda,a)$ and $g(\lambda,a)$, the
scattering cross sections and asymmetry factors of spherical porous aggregates
of radii $a$ at wavelength $\lambda$. Here the aggregate size $a$ refers to
the radius of the sphere encompassing the entire aggregate (we assume that all
grains are spherical in shape). We assume a power-law dust size distribution,
$dn/da\propto a^{-\alpha}$, over the size range $a_{\rm min}\leq a\leq a_{\rm
max}$. The model reflectivity, $S^{\rm mod}(\lambda)$, is calculated as

$S^{\rm mod}(\lambda)\propto\int^{a_{\rm max}}_{a_{\rm min}}C_{\rm
sca}(\lambda,a)\Phi_{\rm
HG}\left[g(\lambda,a);\theta\right]\frac{dn}{da}da~{}~{},$ (7)

where $\theta$, defined as $\pi-\phi$, is the scattering angle, and $\Phi_{\rm
HG}\left[g(\lambda,a);\theta\right]$ is the Henyey-Greenstein phase function
(Henyey & Greenstein, 1941). Note that the right-hand side differs from the
absolute reflectivity by some scaling factor. However, similar to $S^{\rm
obs}(\lambda)$, $S^{\rm mod}(\lambda)$ is normalized with respect to the
averaged value over the defined continuum wavelengths to produce $S^{\rm
mod}_{\lambda}/\langle S^{\rm mod}\rangle$ during which the scaling factor is
cancelled out. It is apparent that the model reflectivity $S^{\rm
mod}(\lambda)$ contains information about the dust size, composition and
morphology since the scattering cross section $C_{\rm sca}(\lambda,a)$ and the
asymmetry factor $g(\lambda,a)$ are dependent on the dust size $a$,
composition (through the refractive index), and morphology (e.g., porosity).

We calculate normalized reflectivities with different size distributions. Four
size ranges are considered, different in the lower cutoff, i.e., $a_{\rm
min}=0.1$, 1, 10, and 100$\,{\rm\mu m}$. The upper cutoff ($a_{\rm max}$) is
fixed to 1000$\,{\rm\mu m}$. The resulting reflectivity gradient is not
sensitive to the choice of $a_{\rm max}$ since the scattering in the optical
by such mm-sized or larger grains is gray (e.g., see Li, 2008). Given a size
range, we vary $\alpha$ from 2.5 to 4.1, the flattest and the steepest power-
law for cometary dust grains from space measurements (McDonnell et al., 1991;
Lasue et al., 2009; Price et al., 2010). The parameters are tabulated in Table
4.

We compare the model reflectivities in Figure 8 with the reflectivity derived
from observations in §3. Models with same $a_{\rm min}$ but different $\alpha$
are grouped and shown as shadows of different colors. Figure 8 shows that the
“observed” reflectivity gradient (i.e., $5\%$ per $10^{3}\,{\rm\AA}$) can be
explained by different combinations of dust parameters. For highly porous dust
with $f_{\rm vac}=0.9$, models with $a_{\rm min}=1\,{\rm\mu m}$ and $a_{\rm
min}=10\,{\rm\mu m}$ are able to reproduce the “observed” reflectivity
gradient. For dust with a lower porosity of $f_{\rm vac}=0.5$, while models
with $a_{\rm min}=10\,{\rm\mu m}$ result in reflectivity gradients too flat to
be comparable with observed, the “observed” reflectivity gradient can be
reproduced by dust with $a_{\rm min}=1\,{\rm\mu m}$. Apparently, although the
limited observational constraint does not allow us to uniquely determine the
dust parameters, the porous dust model is capable of accounting for the
“observed” reflectivity gradient with appropriate combination of parameters.
In general, the reflectivity gradient $S^{\rm mod}(\lambda)$ decreases as the
dust becomes more porous or $a_{\rm min}$ becomes smaller. The variation of
$S^{\rm mod}(\lambda)$ with $\alpha$, the power index of the dust size
distribution, depends on the choice of $a_{\rm min}$. For $a_{\rm
min}<1\,{\rm\mu m}$, $S^{\rm mod}(\lambda)$ becomes bluer as $\alpha$
increases, due to the increasing role played by submicron-sized grains which
scatter more at shorter wavelengths. By increasing $a_{\rm min}$ to several
microns (and larger), $S^{\rm mod}(\lambda)$ becomes redder as $\alpha$
increases since these micron-sized grains scatter more effectively at longer
wavelengths. With $a_{\rm min}$ reaching several tens of microns, $S^{\rm
mod}(\lambda)$ becomes gray and does not vary with $\alpha$ since these grains
are in the geometrical optics limit and their scattering does not vary much
with wavelength. Apparently, $a_{\rm min}$ plays a more important role than
$\alpha$ in determining $S^{\rm mod}(\lambda)$.

##  5 Summary

We presented long-slit optical spectra at $\sim\,$3600–8700$\,{\rm\AA}$ and
multi-band $BVRI$ images of the disintegrating comet 19Y4, from which we
derived the gas and dust production rates and modeled the dust-scattering
reflectivity. Our principal results are as follows:

  1. 1.

All of the spectra taken during the disintegration showed strong gas emission.
The production rates of CN, C2, and C3 and the corresponding C2/CN and C3/CN
production-rate ratios were derived. Both the gas production rates and their
dependencies on $r_{\rm h}$ are consistent with that of “typical” Oort cloud
comets, not that of carbon-chain “depleted” Kuiper Belt comets.

  2. 2.

The dust-scattering reflectivity of 19Y4 was somewhat red, with a gradient of
$\sim\,$5% per $10^{3}\,{\rm\AA}$. The reflectivity gradient was modeled in
terms of porous dust. The red color is accounted for by porous dust.

We thank the anonymous referees for their valuable comments and suggestions
which improved the quality and presentation of this work. RNZ thanks Dan Yu
for her love and companionship over the past eight years, and dedicates this
article to their marriage. We acknowledge the support of the staff of the
Xinglong 2.16m telescope. This work was partially Supported by the Open
Project Program of the CAS Key Laboratory of Optical Astronomy, National
Astronomical Observatories, Chinese Academy of Sciences. RNZ and JFL are
supported by NSFC through grant numbers 11988101 and 11933004. LW is supported
by NSFC through grant number U2031144. JFL acknowledges the support from the
New Cornerstone Science Foundation through the New Cornerstone Investigator
Program and the XPLORER PRIZE.

Table 1: Log of BFOSC observations. Date of 2020 | Object | Grism/Filter | Slit Width | Exposure | Airmass |  $r_{\rm h}$aaHeliocentric distance | $\dot{r_{\rm h}}$ |  $\Delta$bbGeocentric distance | $\dot{\Delta}$ |  $\phi$ccPhase angle
---|---|---|---|---|---|---|---|---|---|---
(UT) |  |  | ($\arcsec$) | (s) |  | (AU) | (${\rm{km\ s^{-1}}}$) | (AU) | (${\rm{km\ s^{-1}}}$) | (∘)
Apr. 6.51 | 19Y4 | G4 | 1.8 | $6\times 600$ | 1.215 | 1.37 | -32.4 | 1.03 | -5.55 | 46.6
Apr. 6.54 | Feige 34 | 600 | 1.012 |  |  |  |  |
Apr. 13.53 | 19Y4 | G4 | 2.3 | $6\times 600$ | 1.279 | 1.24 | -33.7 | 1.00 | -6.03 | 51.8
Apr. 13.56 | He 3 | 600 | 1.389 |  |  |  |  |
Apr. 20.53 | 19Y4 | G4 | 2.3 | $4\times 600$ | 1.404 | 1.10 | -35.4 | 0.98 | -7.33 | 57.2
Apr. 20.56 | HR 4554 | 5 | 1.037 |  |  |  |  |
Apr. 20.59 | HD 109995 | 120 | 1.014 |  |  |  |  |
Apr. 23.55 | 19Y4 | BVRI |  | $3\times 300$ | 1.526 | 1.04 | -36.1 | 0.97 | -8.10 | 59.9
Apr. 23.58 | 19Y4 | G4 | 1.8 | $1\times 600$ | 1.763 | 1.04 | -36.1 | 0.97 | -8.10 | 59.9
Apr. 23.59 | HD 109995 | 30 | 1.005 |  |  |  |  |
Apr. 23.67 | HZ 44 | 360 | 1.012 |  |  |  |  |
Apr. 23.84 | Feige 98 | 600 | 1.412 |  |  |  |  |
Table 2: Gas production rates and scale lengths.aaThe production rates of Ivanova et al. (2021) were merely upper limits. Also, when Ivanova et al. (2021) derived the production rates from the Haser model, they adopted a set of scale lengths which are different from ours. Note that our rates are consistent with that derived from the TRAPPIST observations (Y. Moulane, private conversation). Date of 2020 | $r_{\rm h}$ | CN | C3 | C2 | Ref.bbReference: (1) This work; (2) Ivanova et al. (2021)
---|---|---|---|---|---
|  | $Q$ | $l_{0}$ | $l_{1}$ | $Q$ | $l_{0}$ | $l_{1}$ | $Q$ | $l_{0}$ | $l_{1}$ |
(UT) | (AU) | ($10^{25}\,{\rm s^{-1}}$) | ($10^{4}\,{\rm km}$) | ($10^{4}\,{\rm km}$) | ($10^{25}\,{\rm s^{-1}}$) | ($10^{4}\,{\rm km}$) | ($10^{4}\,{\rm km}$) | ($10^{25}\,{\rm s^{-1}}$) | ($10^{4}\,{\rm km}$) | ($10^{4}\,{\rm km}$) |
Apr. 6.51 | 1.37 | $4.4^{+0.3}_{-0.3}$ | $6.5^{+0.5}_{-0.5}$ | $38^{+6}_{-6}$ | $-$ | $-$ | $-$ | $3.5^{+0.4}_{-0.4}$ | $5.7^{+0.7}_{-0.6}$ | $13^{+3}_{-2}$ | 1
Apr. 13.53 | 1.24 | $3.7^{+0.2}_{-0.2}$ | $3.5^{+0.2}_{-0.2}$ | $35^{+8}_{-6}$ | $0.21^{+0.01}_{-0.01}$ | $1.0^{+0.1}_{-0.1}$ | $27^{+5}_{-4}$ | $4.1^{+0.1}_{-0.2}$ | $7.3^{+0.2}_{-0.3}$ | $8.0^{+0.3}_{-0.2}$ | 1
Apr. 14.74 | 1.21 |  $<$1.05 | 3.2 | 43 |  $<$0.015 | 0.15 | 8.8 |  $<$0.83 | 2.4 | 16 | 2
Apr. 16.77 | 1.17 |  $<$1.12 | 3.0 | 41 |  $<$0.015 | 0.14 | 8.2 |  $<$0.79 | 2.2 | 15 | 2
Apr. 20.53 | 1.10 | $1.8^{+0.3}_{-0.3}$ | $3.2^{+0.6}_{-0.5}$ | $8^{+2}_{-2}$ | $0.11^{+0.01}_{-0.01}$ | $1.3^{+0.1}_{-0.1}$ | $16^{+4}_{-3}$ | $3.4^{+0.1}_{-0.2}$ | $7.2^{+0.3}_{-0.3}$ | $8.7^{+0.5}_{-0.4}$ | 1
Apr. 23.58 | 1.04 | $1.3^{+0.1}_{-0.1}$ | $4.1^{+0.3}_{-0.3}$ | $30^{+6}_{-5}$ | $-$ | $-$ | $-$ | $1.5^{+0.1}_{-0.1}$ | $5.6^{+0.5}_{-0.4}$ | $10^{+1}_{-1}$ | 1
Table 3: $Af\rho$ within $\rho=10000\,{\rm km}$ in $\lambda 4850\,{\rm\AA}$ and $\lambda 5240\,{\rm\AA}$ bands.  Date of 2020 | $r_{\rm h}$ | $Af\rho$
---|---|---
|  | $\lambda 4850\,{\rm\AA}$ | $\lambda 5240\,{\rm\AA}$
(UT) | (AU) | (cm) | (cm)
Apr. 6.51 | 1.37 | 972$\pm$10 | 753$\pm$9
Apr. 13.53 | 1.24 | 288$\pm$1 | 199$\pm$1
Apr. 20.53 | 1.10 | 155$\pm$2 | 97$\pm$1
Apr. 23.58 | 1.04 | 145$\pm$3 | 73$\pm$3
Table 4: Dust model parameters. $f_{\rm vac}$ | $f_{\rm carb}\left/f_{\rm sil}\right.$ | $a_{\rm min}$ | $a_{\rm max}$ | $\alpha$
---|---|---|---|---
|  | ($\,{\rm\mu m}$) | ($\,{\rm\mu m}$) |
|  | 0.1 |  |
0.5 | 3 | 1 | 1000 | 2.5-4.1
|  | 10 |  |
0.9 | 3 | 0.1 | 1000 | 2.5-4.1
1
10
100
Figure 1: Slit views of C/2019 Y4 (ATLAS) on 2020 April 13.53 and 20.53 (UT).
The ephemeris positions of the two major fragments A and B (from MPEC 2020-L06
and 2023-J29, respectively) are marked. On April 13.53 (UT), the slit was
centered on A. On April 20.53 (UT), the slit was over both A and B. The
rectangles, 4.1$\arcsec$ in height and 2.3$\arcsec$ in width, along the slits
are sub-apertures. They were used to extract sub-aperture spectra (see §3).
The arrows in the upper right corner mark the north and the east directions,
and those in the lower right corners mark the velocity and comet-Sun vectors.
Figure 2: Narrow-aperture BFOSC spectra of C/2019 Y4 (ATLAS) obtained on 2020
April 6.51 (blue line), April 13.53 (yellow line), April 20.53 (green line),
and April 23.58 (UT; red line). For comparison, the spectra are vertically
shifted by some constants. Strong emission bands are marked.  Figure 3: BVRI
images of C/2019 Y4 (ATLAS) on 2020 April 23.55 (UT). The coma shows a less-
extended morphology from B to I, indicating the broadband colors vary from
inner to outer coma. The white rectangle is a slit-like aperture. It is
divided into a series of sub-apertures, 4.1″in height and 2.3″in width, to
derive $\left(B-V\right)$ and $\left(V-R\right)$ profiles (see §3). The arrows
in the upper right corner mark the north and the east directions, and those in
the lower right corners mark the velocity and comet-Sun vectors.

Figure 4: Comparison of the $\left(B-V\right)$ (left panel) and
$\left(V-R\right)$ (right panel) profiles of C/2019 Y4 (ATLAS) derived from
photometry (blue open squares) and spectroscopy (yellow open squares), both on
2020 April 23 (UT). Small sections of the spectroscopic profiles are
contaminated by field stars (grey open squares). Also shown are the solar
colors (red lines). The spectroscopic profiles systematically deviate from the
photometric profiles by $-0.19$ and $-0.20$ in $\left(B-V\right)$ and
$\left(V-R\right)$, respectively, which are shown as black open squares with
fitted black dashed lines. Note that the error bars only account for the
propagated Poisson errors. The omitted systematic errors can be approximated
by the standard deviations of $\Delta\left(B-V\right)$ and
$\Delta\left(V-R\right)$, which are $\sim\,$0.15 and 0.08 mag, respectively.
Figure 5: The “observed” reflectivities ($S^{\rm obs}_{\lambda}/\langle S^{\rm
obs}\rangle$; gray lines) of C/2019 Y4 (ATLAS) on 2020 April 6.51 (a), 13.53
(b), 20.53 (c), and 23.58 (UT; d). The gradients are determined by fitting
those in the dust continuum bands (black open squares with error bars) with
linear functions (red lines with light red shadows representing the fitting
uncertainties). As the spectra of C/2019 Y4 (ATLAS) are somewhat subjected to
atmospheric dispersion, the gradients need correction. To the first-order
approximation, the reflectivity gradient on April 23.58 (UT) is corrected to
$\sim\,$5% per $10^{3}\,{\rm\AA}$ (blue line; see §3 for details). All
quantities are normalized to their mean values in the observed wavelength
range. Figure 6: Comparison of the “observed” column density profiles (filled
black circles with error bars) of CN ($\Delta\nu=0$; upper panels), C3
($\Delta\nu=0$; middle panels), and C2($\Delta\nu=0$; bottom panels) with
model fits (solid red lines with light red shadows representing the model
uncertainties). In each panel the model parameters are labelled. The open
black circles are not considered for modeling (see §4.1 for details).  Figure
7: Top panel (a): Variation of the CN production rate Q(CN) of C/2019 Y4
(ATLAS) (yellow circles) with heliocentric distance $r_{\rm h}$. Also shown
are the CN production rates of the Langland-Shula & Smith (2011) sample of 26
comets (black circles). Middle panel (b): Same as (a) but for $Af\rho$. Also
shown is a linear fit to the $r_{\rm h}$-dependence of $Af\rho$ for the
Langland-Shula & Smith (2011) sample: $d(\log Af\rho)/d(r_{\rm h})=5.3\pm 0.6$
(black dashed line). Bottom panel: Same as (a) but for $Af\rho/Q({\rm CN})$.
The power-law has an index of $4.6\pm 0.8$.  Figure 8: Comparison of the dust
models with the recovered reflectivity gradient (black line). For the dust,
two cases are considered for the porosity (i.e., the fractional volumn of
vacuum): $f_{\rm vac}=0.9$ (left panel) and $f_{\rm vac}=0.5$ (right panel).
The volume mixing ratio of amorphous silicate and amorphous carbon is fixed to
$f_{\rm carb}/f_{\rm sil}=3$. For the size distribution, four size ranges are
considered, different in the lower cutoff, i.e., $a_{\rm min}=0.1$ (gray
shadow), 1 (blue shadow), 10 (red shadow), and 100$\,{\rm\mu m}$ (yellow
shadow). The upper cutoff ($a_{\rm max}$) is fixed to 1000$\,{\rm\mu m}$.
Given a size range, the power-law index $\alpha$ is varied from 2.5 to 4.1.
The recovered reflectivity gradient ($\sim\,$5% per $10^{3}\,{\rm\AA}$) is
barely explained by both $1\leq a\leq 1000\,{\rm\mu m}$ and $10\leq a\leq
1000\,{\rm\mu m}$ when $f_{\rm vac}=0.9$, while well explained by the size
range $1\leq a\leq 1000\,{\rm\mu m}$ when $f_{\rm vac}=0.5$. For both cases,
larger or smaller $a_{\rm min}$ leads to flatter gradient (or “gray” in
color).

## References

  * A’Hearn et al. (1995) A’Hearn, M. F., Millis, R. C., Schleicher, D. O., Osip, D. J., & Birch, P. V. 1995, Icarus, 118, 223, doi: 10.1006/icar.1995.1190
  * A’Hearn et al. (1984) A’Hearn, M. F., Schleicher, D. G., Millis, R. L., Feldman, P. D., & Thompson, D. T. 1984, AJ, 89, 579, doi: 10.1086/113552
  * Ahumada et al. (2022) Ahumada, R., Allende Prieto, C., Almeida, A., & et al. 2022, VizieR Online Data Catalog, V/154
  * 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-Specutils Development Team (2019) Astropy-Specutils Development Team. 2019, Specutils: Spectroscopic analysis and reduction, Astrophysics Source Code Library, record ascl:1902.012.  http://ascl.net/1902.012
  * Bardyn et al. (2017) Bardyn, A., Baklouti, D., Cottin, H., et al. 2017, MNRAS, 469, S712, doi: 10.1093/mnras/stx2640
  * Blum & Wurm (2008) Blum, J., & Wurm, G. 2008, ARA&A, 46, 21, doi: 10.1146/annurev.astro.46.060407.145152
  * Bohlin et al. (2014) Bohlin, R. C., Gordon, K. D., & Tremblay, P. E. 2014, PASP, 126, 711, doi: 10.1086/677655
  * Bohren & Huffman (1983) Bohren, C. F., & Huffman, D. R. 1983, Absorption and scattering of light by small particles
  * Bradley et al. (2019) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2019, astropy/photutils: v0.6, doi: 10.5281/zenodo.2533376
  * Brown et al. (1996) Brown, M. E., Bouchez, A. H., Spinrad, A. H., & Johns-Krull, C. M. 1996, AJ, 112, 1197, doi: 10.1086/118090
  * Chonis & Gaskell (2008) Chonis, T. S., & Gaskell, C. M. 2008, AJ, 135, 264, doi: 10.1088/0004-6256/135/1/264
  * Draine & Lee (1984) Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89, doi: 10.1086/162480
  * Fan et al. (2016) Fan, Z., Wang, H., Jiang, X., et al. 2016, PASP, 128, 115005, doi: 10.1088/1538-3873/128/969/115005
  * Fink & Hicks (1996) Fink, U., & Hicks, M. D. 1996, ApJ, 459, 729, doi: 10.1086/176938
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fray et al. (2005) Fray, N., Bénilan, Y., Cottin, H., Gazeau, M. C., & Crovisier, J. 2005, Planet. Space Sci., 53, 1243, doi: 10.1016/j.pss.2005.06.005
  * Ginsburg et al. (2019) Ginsburg, A., Sipőcz, B. M., Brasseur, C. E., et al. 2019, AJ, 157, 98, doi: 10.3847/1538-3881/aafc33
  * Greenberg (1998) Greenberg, J. M. 1998, A&A, 330, 375
  * Greenberg & Hage (1990) Greenberg, J. M., & Hage, J. I. 1990, ApJ, 361, 260, doi: 10.1086/169191
  * Henyey & Greenstein (1941) Henyey, L. G., & Greenstein, J. L. 1941, ApJ, 93, 70, doi: 10.1086/144246
  * Hui & Ye (2020) Hui, M.-T., & Ye, Q.-Z. 2020, AJ, 160, 91, doi: 10.3847/1538-3881/ab9d81
  * Ivanova et al. (2021) Ivanova, O., Luk’yanyk, I., Tomko, D., & Moiseev, A. 2021, MNRAS, 507, 5376, doi: 10.1093/mnras/stab2306
  * Jewitt & Meech (1986) Jewitt, D., & Meech, K. J. 1986, ApJ, 310, 937, doi: 10.1086/164745
  * Lang et al. (2010) Lang, D., Hogg, D. W., Mierle, K., Blanton, M., & Roweis, S. 2010, AJ, 139, 1782, doi: 10.1088/0004-6256/139/5/1782
  * Langland-Shula & Smith (2011) Langland-Shula, L. E., & Smith, G. H. 2011, Icarus, 213, 280, doi: 10.1016/j.icarus.2011.02.007
  * Lasue et al. (2009) Lasue, J., Levasseur-Regourd, A. C., Hadamcik, E., & Alcouffe, G. 2009, Icarus, 199, 129, doi: 10.1016/j.icarus.2008.09.008
  * Li (2008) Li, A. 2008, in Small bodies in planetary systems (Berlin: Springer), 167–188
  * Li & Greenberg (1998a) Li, A., & Greenberg, J. M. 1998a, A&A, 338, 364
  * Li & Greenberg (1998b) —. 1998b, ApJ, 498, L83, doi: 10.1086/311308
  * McDonnell et al. (1991) McDonnell, J. A. M., Lamy, P. L., & Pankiewicz, G. S. 1991, Physical Properties of Cometary Dust, ed. J. Newburn, R. L., M. Neugebauer, & J. Rahe, Vol. 167, 1043
  * Mommert et al. (2019) Mommert, M., Kelley, M., de Val-Borro, M., et al. 2019, The Journal of Open Source Software, 4, 1426, doi: 10.21105/joss.01426
  * 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
  * Price et al. (2010) Price, M. C., Kearsley, A. T., Burchell, M. J., et al. 2010, Meteoritics and Planetary Science, 45, 1409, doi: 10.1111/j.1945-5100.2010.01104.x
  * Ramírez et al. (2012) Ramírez, I., Michel, R., Sefako, R., et al. 2012, ApJ, 752, 5, doi: 10.1088/0004-637X/752/1/5
  * Rouleau & Martin (1991) Rouleau, F., & Martin, P. G. 1991, ApJ, 377, 526, doi: 10.1086/170382
  * Schleicher (2010) Schleicher, D. G. 2010, AJ, 140, 973, doi: 10.1088/0004-6256/140/4/973
  * Schleicher & Bair (2011) Schleicher, D. G., & Bair, A. N. 2011, AJ, 141, 177, doi: 10.1088/0004-6256/141/6/177
  * Sekanina et al. (2020) Sekanina, Z., Guido, E., & Bryssinck, E. 2020, Central Bureau Electronic Telegrams, 4751, 1
  * Steele et al. (2020) Steele, I. A., Smith, R. J., & Marchant, J. 2020, The Astronomer’s Telegram, 13622, 1
  * Strazzulla (1999) Strazzulla, G. 1999, Space Sci. Rev., 90, 269, doi: 10.1023/A:1005270703284
  * Swings (1941) Swings, P. 1941, Lick Observatory Bulletin, 508, 131, doi: 10.5479/ADS/bib/1941LicOB.19.131S
  * Ye & Hui (2020) Ye, Q., & Hui, M.-T. 2020, The Astronomer’s Telegram, 13651, 1
  * Ye & Zhang (2020) Ye, Q., & Zhang, Q. 2020, The Astronomer’s Telegram, 13620, 1
  * Ye et al. (2021) Ye, Q., Jewitt, D., Hui, M.-T., et al. 2021, AJ, 162, 70, doi: 10.3847/1538-3881/abfec3
  * Zubko et al. (2020) Zubko, E., Zheltobryukhov, M., Chornaya, E., et al. 2020, MNRAS, 497, 1536, doi: 10.1093/mnras/staa1725


$v_{s}v_{i,1}v_{i,2}\ldots v_{i,k-3}v_{i,k-2}v_{i,k-1}v^{c}_{j}$ for some $i$
(we call this a $\top$-path for $x_{i}$), or of the form
$v_{s}v_{i,1}v_{i,2}\ldots v_{i,k-4}v_{i,k-2}v_{i,k-1}v_{i,k}v^{c}_{j}$ for
some variable $x_{i}$ (we will call this a $\bot$-path for $x_{i}$). If we
have a $\top$-path for $x_{i}$ in $(G,\lambda^{\prime})$ then
$\lambda^{\prime}(v_{i,k-2}v_{i,k-1})\leq\delta+k-1$, and if we have a
$\bot$-path for $x_{i}$ in $(G,\lambda^{\prime})$ then
$\lambda^{\prime}(v_{i,k-2}v_{i,k-1}\geq 2\delta+k-1$.

Next we show that $((G,\lambda),v_{s},k,1,m)$ is a yes-instance for TSEP if
and only if $\bigwedge_{j}C_{j}$ is satisfiable. First, assume that
$\bigwedge_{j}C_{j}$ is satisfiable and that
$\phi:\\{x_{i}\\}\mapsto\\{\top,\bot\\}$ is a truth assignment that satisfies
$\bigwedge_{j}C_{j}$, and construct $\lambda^{\prime}$ as follows.

$\displaystyle\forall i\quad$ $\displaystyle\lambda^{\prime}(v_{s}v_{i,1})=1$
$\displaystyle\forall i\,\forall\ell\in[k-4]\quad$
$\displaystyle\lambda^{\prime}(v_{i,\ell}v_{i,\ell+1})=\ell+1$
$\displaystyle\forall i\quad$
$\displaystyle\lambda^{\prime}(v_{i,k-4}v_{i,k-2})=2\delta+k-2$
$\displaystyle\forall i\quad$
$\displaystyle\lambda^{\prime}(v_{i,k-1}v_{i,k})=2\delta+k$
$\displaystyle\forall i,j\quad$
$\displaystyle\lambda^{\prime}(v_{i,k-1}v^{c}_{j})=\delta+k$
$\displaystyle\forall i,j\quad$
$\displaystyle\lambda^{\prime}(v_{i,k}v^{c}_{j})=2\delta+k+1$
$\displaystyle\forall i\,\text{ where }\phi(x_{i})=\top\quad$
$\displaystyle\lambda^{\prime}(v_{i,k-2}v_{i,k-1})=k-1$ $\displaystyle\forall
i\,\text{ where }\phi(x_{i})=\bot\quad$
$\displaystyle\lambda^{\prime}(v_{i,k-2}v_{i,k-1})=2\delta+k-1$

This satisfies $\lambda^{\prime}$ is a $(\delta,\zeta)$-perturbation of
$\lambda$, and note that in particular, if $\phi(x_{i})=\top$ then we have a
$\top$-path for $x_{i}$ in $(G,\lambda^{\prime})$, and if $\phi(x_{i})=\bot$
then we have a $\bot$-path for $x_{i}$ in $(G,\lambda^{\prime})$.

We now consider shortest paths to various groups of vertices. First, for any
variable $x_{i}$, if $\phi(x_{i})=\top$ then $v_{s}v_{i,1}v_{i,2}\ldots
v_{i,k-3}v_{i,k-2}v_{i,k-1}v_{i,k}$ is a temporal path of length $k$, and if
$\phi(x_{i})=\bot$ then $v_{s}v_{i,1}v_{i,2}\ldots
v_{i,k-4}v_{i,k-2}v_{i,k-1}v_{i,k}$ is a temporal path of length $k-1$ and
$v_{s}v_{i,1}v_{i,2}\ldots v_{i,k-4}v_{i,k-3}$ is a temporal path of length
$k-3$. It thus remains to see whether vertices of the form $v^{c}_{j}$ are
reachable in paths of length $k$. Take an arbitrary such vertex $v^{c}_{j}$;
as $\phi$ satisfies $\bigwedge_{j}C_{j}$, at least one literal in $C_{j}$ must
be true. If this literal is of the form $x_{i}$ (i.e., $x_{i}$ appears non-
negated in $C_{j}$), then $\phi(x_{i})=\top$ and so we have a $\top$-path for
$x_{i}$ in $(G,\lambda^{\prime})$ and thus the vertex $v^{c}_{j}$ is reachable
by a path of length $k$. Next, if this literal is of the form
$\overline{x_{i}}$ (i.e., $x_{i}$ appears negated in $C_{j}$), then
$\phi(x_{i})=\bot$ and so we have a $\bot$-path for $x_{i}$ in
$(G,\lambda^{\prime})$ and thus the vertex $v^{c}_{j}$ is reachable by a path
of length $k$. Thus as $\lambda^{\prime}$ is a $(\delta,m)$-perturbation,
$((G,\lambda),v_{s},k,\delta,m)$ is a yes-instance for TSEP.

Next, assume that $((G,\lambda),v_{s},k,1,m)$ is a yes-instance for TSEP. Then
there must be a $(\delta,m)$-perturbation $\lambda^{\prime}$ of $\lambda$,
such that for any clause vertex $v^{c}_{j}$, there is a path of length $k$
from $v_{s}$ to $v^{c}_{j}$. We now define
$\phi:\\{x_{i}\\}\mapsto\\{\top,\bot\\}$. For each edge of the form
$v_{i,k-2}v_{i,k-1}$, if $\lambda^{\prime}(v_{i,k-2}v_{i,k-1})\leq\delta+k-1$
then let $\phi(x_{i})=\top$, else let $\phi(x_{i})=\bot$. Now consider an
arbitrary clause $C_{j}$, and its associated clause vertex $v^{c}_{j}$. There
must exist a temporal path of length $k$ in $(G,\lambda^{\prime})$ from
$v_{s}$ to $v^{c}_{j}$, and by the construction of $(G,\lambda)$ and the fact
that $\lambda^{\prime}$ is a $\delta$-perturbation, it must be that this path
is either either a $\top$-path for some variable $x_{i}$, or a $\bot$-path for
some variable $x_{i}$. If this path is a $\top$-path for some variable $x_{i}$
then $\lambda^{\prime}(v_{i,k-2}v_{i,k-1}\leq\delta+k-1$ and so
$\phi(x_{i})=\top$, and by construction, the literal $x_{i}$ appears non-
negated in $C_{j}$, so $C_{j}$ is satisfied. Alternatively, if this path is a
$\bot$-path for some variable $x_{i}$ then
$\lambda^{\prime}(v_{i,k-2}v_{i,k-1}\geq 2\delta+k-1>\delta+k-1$ and so
$\phi(x_{i})=\bot$, and by construction, the literal $x_{i}$ appears negated
in $C_{j}$, so $C_{j}$ is satisfied. Since we took an arbitrary clause vertex
$v^{c}_{j}$, it follows that $\phi$ satisfies $\bigwedge_{j}C_{j}$.

$v_{s}$$v_{1,1}$$v_{1,2}$$v_{1,k-4}$$v_{1,k-3}$$v_{1,k-2}$$v_{1,k-1}$$v_{1,k}$$v^{c}_{1}$$v^{c}_{2}$$1$$2$$k-3$$3\delta+k-2$$k-2$$\delta+k-1$$3\delta+k$$k$$3\delta+k+1$
Figure 2: Diagram of a partial construction for Theorem 5.1 showing the
relevant vertices for variable $x_{1}$ and clauses $C_{1}$ and $C_{2}$. In
this construction, clause $C_{1}$ contains the literal $x_{1}$ and clause
$C_{2}$ contains the literal $\overline{x_{1}}$. This diagram is labelled with
time labels from $\lambda$, and the dashed line indicates a path on the
vertices $v_{i,2}\ldots v_{i,k-4}$, where
$\lambda(v_{i,\ell}v_{i,\ell+1})=\ell+2$. Note that, as $\delta=1$, if
$\lambda^{\prime}$ is a $1$-perturbation of $\lambda$ and contains a path of
either length or duration $k$ from $v_{s}$ to $v^{c}_{1}$ through $v_{i,k-1}$
it must be that $\lambda^{\prime}(v_{i,k-2}v_{i,k-1})=k-1$, while if
$(G,\lambda^{\prime})$ is to contain a path of length or duration $k$ from
$v_{s}$ to $v^{c}_{2}$ through vertex $v_{i,k-1}$ it must that
$\lambda^{\prime}(v_{i,k-2}v_{i,k-1})=k$. This dichotomy translates to setting
the variable $x_{i}$ to either true (if
$\lambda^{\prime}(v_{i,k-2}v_{i,k-1})=k-1$) or false (if
$\lambda^{\prime}(v_{i,k-2}v_{i,k-1})=k$).

We give an analogous result for TFaEP, which again involves a reduction from
SAT.

######  Theorem 5.3.

For any integers $k\geq 2$ and any $\delta\geq 1$, TFaEP is -hard even if
$\zeta=m$.

######  Proof 5.4.

We show how, given integers $k\geq 2$ and $\delta$, an instance
$((G,\lambda),k,\delta,m)$ (where $m=|E(G)|$) of TFaEP can be constructed from
an instance $\bigwedge C_{j}$ of SAT such that $((G,\lambda),k,\delta,m))$ is
a yes-instance of TFaEP if and only if $\bigwedge C_{j}$ is a yes-instance of
SAT. First, we take our instance $\bigwedge C_{j}$ of SAT and construct our
temporal graph $(G,\lambda)$, beginning by creating the vertex $v_{s}$. Then
for each variable $x_{i}$, define $v_{i},0=v_{s}$, and create the vertices
$v_{i,\ell}$ for $\ell\in[k-1]$. Lastly, for each clause $C_{j}$ in the SAT
instance create the vertex $v^{c}_{j}$. We next create the temporal edges.
First, for each variable $x_{i}$ add the edge $(v_{s}v_{i,1},1)$. Then for
each variable $x_{i}$, and for each $\ell\in[k-2]$, add the edge
$(v_{i,\ell}v_{i,\ell+1},\ell+1)$. Lastly, for each variable $x_{i}$, and each
clause $C_{j}$, if $x_{i}$ appears as a literal (i.e., non-negated) in
$C_{j}$, add the edge $(v_{i,k-1}v^{c}_{j},k+2\delta)$, and if
$\overline{x_{i}}$ appears as a literal in $C_{j}$, add the edge
$(v_{i,k-1}v^{c}_{j},k-1)$. A diagram of this gadget is shown in Figure 3.

We now note some properties of this construction for an arbitrary
$\delta$-perturbation $\lambda^{\prime}$ of $\lambda$. First, consider some
variable $x_{i}$, and some clause $C_{j}$. There can be no temporal path from
$v_{s}$ to $v^{c}_{j}$ of length strictly less than $k$, and so for any
temporal path from $v_{s}$ to $v^{c}_{j}$ in $(G,\lambda^{\prime})$ of
duration at most $k$, if edge $e^{\prime}$ immediately follows edge $e$, then
$\lambda^{\prime}(e^{\prime})=\lambda^{\prime}(e)+1$. Now consider some
variable $x_{i}$, and some clause $C_{j}$ such that $x_{i}$ only appears non-
negated in $C_{j}$. To have a valid temporal path of duration $k$ in
$(G,\lambda^{\prime})$ from $v_{s}$ to $v^{c}_{j}$ through $v_{i,k-1}$, we
must have $\lambda^{\prime}(v_{i,k-2}v_{i,k-1})=k-1+\delta$ and
$\lambda^{\prime}(v_{i,k-1}v^{c}_{j})=k+\delta$. On the other hand, if $x_{i}$
only appears negated in $C_{j}$, then for a valid temporal path of duration
$k$ in $(G,\lambda^{\prime})$ from $v_{s}$ to $v^{c}_{j}$ through $v_{i,k-1}$
to exist, we must have $\lambda^{\prime}(v_{i,k-1}v^{c}_{j})\leq k-1+\delta$,
and so we must have $\lambda^{\prime}(v_{i,k-2}v_{i,k-1})\leq k-2+\delta$. We
later use this dichotomy to encode whether the variable $x_{i}$ is true or
false.

We now show that $((G,\lambda),v_{s},k,\delta,m)$ is a yes-instance to TFaEP
if and only if $\bigwedge_{j}C_{j}$ is a yes-instance to SAT. First, assume
that $\bigwedge_{j}C_{j}$ is satisfiable, and let $\phi$ be a function mapping
variables to $\\{\top,\bot\\}$ that satisfies $\bigwedge_{j}C_{j}$. We will
use $\phi$ to construct a a $\delta$-perturbation $\lambda^{\prime}$ of
$\lambda$. For each variable $x_{i}$, if $\phi(x_{i})=\top$ then let
$\lambda^{\prime}(v_{s}v_{i,1},1+\delta)$, for each $\ell\in[k-2]$ let
$\lambda^{\prime}(v_{i,\ell}v_{i,\ell+1})=\ell+1+\delta$, and for each clause
$C_{j}$ that contains $x_{i}$ as a literal (i.e. non-negated) let
$\lambda^{\prime}(v_{i,k-1}v^{c}_{j})=k+\delta$. For each variable $x_{i}$, if
$\phi(x_{i})=\bot$ then let $\lambda(v_{s}v_{i,1},1)$, for each $\ell\in[k-2]$
let $\lambda^{\prime}(v_{i,\ell}v_{i,\ell+1})=\ell+1$, and for each clause
$C_{j}$ that contains $x_{i}$ as a literal (i.e. non-negated) let
$\lambda^{\prime}(v_{i,k-1}v^{c}_{j})=k$. We have not yet defined
$\lambda^{\prime}$ for some edges $e$ yet; for such edges, let
$\lambda^{\prime}(e)=\lambda(e)$. By this construction, $\lambda^{\prime}$ is
a $\delta$-perturbation of $\lambda$, so we now consider
$\operatorname{ecc^{Fa}}((G,\lambda^{\prime}),v_{s})$. First, for any variable
$x_{i}$, the path $v_{s}v_{i,1}\ldots v_{i,k-2}v_{i,k-1}$ is a path of
duration $k-1$ in $(G,\lambda^{\prime})$, so we need only consider vertices of
the form $v^{c}_{j}$. Consider an arbitrary such vertex $v^{c}_{j}$, and let
$C_{j}$ be its associated clause. As $\phi$ satisfies $C_{j}$, there must be
some variable $x_{i}$ such that either $\phi(x_{i})=\top$ and $x_{i}$ appears
as a literal in $C_{j}$, or $\phi(x_{i})=\bot$ and $\overline{x_{i}}$ appears
as a literal in $C_{j}$. In either case, the path $v_{s}v_{i,1}\ldots
v_{i,k-2}v_{i,k-1}v^{c}_{j}$ is a temporal path of duration $k$ in
$(G,\lambda^{\prime})$, and so
$\operatorname{ecc^{Fa}}((G,\lambda^{\prime}),v_{s})=k$ and so
$((G,\lambda),v_{s},k,\delta,m)$ is a yes-instance to TFaEP.

Next, assume that $((G,\lambda),v_{s},k,\delta,m)$ is a yes-instance to TFaEP,
and let $\lambda^{\prime}$ be the associating perturbation that achieves
$\operatorname{ecc^{Fa}}((G,\lambda^{\prime}),v_{s})\leq k$. Construct a
$\phi$ mapping variables to $\\{\top,\bot\\}$ as follows. If
$\lambda^{\prime}(v_{i,k-2}v_{i,k-1})=k-1+\delta$ let $\phi(x_{i})=\top$, else
let $\phi(x_{i})=\top$. Now take an arbitrary clause $C_{j}$, and associated
vertex $v^{c}_{j}$. There is a temporal path of duration $k$ in
$(G,\lambda^{\prime})$ from $v_{s}$ to $v^{c}_{j}$, as we have a yes-instance
to TFaEP. By construction, such a path must use an edge of the form
$v_{i,k-2}v_{i,k-1}$ for some $i$ corresponding to a variable $i$. If
$\lambda^{\prime}(v_{i,k-2}v_{i,k-1})=k-1+\delta$ then $\phi(x_{i})=\top$ and
the edge $v_{i,k-1}v^{c}_{j}$ must be used at time $k+\delta$, and so by
construction (and as $\lambda^{\prime}$ is a $\delta$-perturbation) it must be
that $x_{i}$ appears non-negated in $C_{j}$, and so $C_{j}$ is satisfied.
Alternatively, $\lambda^{\prime}(v_{i,k-2}v_{i,k-1})<k-1+\delta$ and
$\phi(x_{i})=\bot$, meaning we must have
$\lambda^{\prime}(v_{i,k-1}v^{c}_{j})<k+\delta$, and so by construction (and
as $\lambda^{\prime}$ is a $\delta$-perturbation) it must be that $x_{i}$
appears negated in $C_{j}$, and so $C_{j}$ is satisfied. As we chose $C_{j}$
arbitrarily, it follows that $\bigwedge_{j}C_{j}$ is satisfied by $\phi$, and
so is a yes-instance to SAT.

$v_{s}$$v_{1,1}$$v_{1,2}$$v_{1,k-2}$$v_{1,k-1}$$v^{c}_{1}$$v^{c}_{2}$$1$$2$$k-1$$k+2\delta$$k-1$
Figure 3: Diagram of a partial construction for Theorem 5.3 showing the
relevant vertices for variable $x_{1}$ and clauses $C_{1}$ and $C_{2}$. In
this construction, clause $C_{1}$ contains the literal $x_{1}$ and clause
$C_{2}$ contains the literal $\overline{x_{1}}$. This diagram is labelled with
time labels from $\lambda$, and the dashed line indicates a path on the
vertices $v_{1,2}\ldots v_{1,k-4}$, where
$\lambda(v_{1,\ell}v_{i,\ell+1})=\ell+1$. Note that, as $\delta=1$, if
$\lambda^{\prime}$ is a $\delta$-perturbation of $\lambda$ and contains a path
of duration at most $k$ from $v_{s}$ to $v^{c}_{1}$ through $v_{1,k-1}$ it
must be that $\lambda^{\prime}(v_{1,k-2}v_{1,k-1})=k+\delta-1$, while if
$(G,\lambda^{\prime})$ is to contain a path of duration at most $k$ from
$v_{s}$ to $v^{c}_{2}$ through vertex $v_{1,k-1}$ it must that
$\lambda^{\prime}(v_{1,k-2}v_{1,k-1})<k+\delta-1$. This dichotomy translates
to setting the variable $x_{1}$ to either true (if
$\lambda^{\prime}(v_{1,k-2}v_{1,k-1})=k-1+\delta$) or false (if
$\lambda^{\prime}(v_{1,k-2}v_{1,k-1})\leq k-2+\delta$).

We conclude this section with one positive result for TSEP and TFaEP. Recall
that TRLP becomes trivial when both $\zeta$ and $\delta$ are sufficiently
large (every instance is a yes-instance). While neither TSEP nor TFaEP becomes
trivial under the same conditions, we now show that both problems can be
solved in polynomial time when both $\zeta$ and $\delta$ are large.

Our algorithms make use of ubiquitous foremost temporal path trees. Recall
that a ubiquitous foremost temporal path is a foremost temporal path $P$ such
that all prefixes of $P$ are themselves foremost temporal paths, and that a
ubiquitous foremost temporal path tree from a source vertex $v_{s}$ is a
temporal graph $(G_{T},\lambda_{T})$ such that $G_{T}$ is a tree,
$\tau(G_{T},\lambda_{T})=1$ and for any $v_{t}\in V(G)$, a foremost temporal
path in $(G_{T},\lambda_{T})$ from $v_{s}$ to $v_{t}$ arrives at $v_{t}$ no
later than any temporal path. Our algorithms further rely on the following
observation: for any edge $e$, the possible set of times it can be perturbed
to includes the range $\\{1,\ldots,r_{e}\\}$ for some suitably large integer
$r_{e}$, which means that along any sequence of (non-temporal) edges
$(e_{1},e_{2},\ldots)$ that can form a temporal path under some suitable
perturbation, we can assume that $\lambda(e_{i})=i$. We write
$d_{G}(v_{1},v_{2})$ for the number of edges in the shortest path between
$v_{1}$ and $v_{2}$ in $G$.

######  Theorem 5.5.

Each of TSEP and TFaEP can be solved in time $O(m(\log
n+\log\delta\tau(G,\lambda)))$ on instances
$((G,\lambda),v_{s},k,\delta,\zeta)$ if $\delta\geq T(G,\lambda)$ and
$\zeta\geq n-1$.

######  Proof 5.6.

Define $\lambda^{\prime}$ by $\lambda^{\prime}(e)=\\{\lambda(e)+p\mid
p\in[-\delta,\delta]\\}$ for $e\in E(G)$, noting that as $\delta\geq
T(G,\lambda)$, for each $e\in E(G)$ there is some positive integer $r_{e}$
such that $\lambda^{\prime}(e)=\\{1,2,\ldots,r_{e}\\}$. Then use Theorem 1.1
to calculate a ubiquitous foremost temporal path tree $(G_{T},\lambda_{T})$ on
$v_{s}$ in the temporal graph $(G,\lambda^{\prime})$. in $O(m(\log
n+\log\tau(G,\lambda^{\prime})))=O(m(\log n+\log\delta\tau(G,\lambda)))$ time.
By setting $\lambda^{\prime\prime}(e)=\lambda_{T}(e)$ for $e\in E(G_{T})$ and
$\lambda^{\prime\prime}(e)=\lambda(e)$ for $e\not\in E(G_{T})$, we know that
$\lambda^{\prime\prime}$ is a $(\delta,\zeta)$-perturbation of $\lambda$ with
the property that if a temporal path exists in $(G_{T},\lambda_{T})$, then the
same path must also exist in $(G,\lambda^{\prime\prime})$. Note that while
$V(G_{T})\subseteq V(G)$, it is possible that $V(G)\neq V(G_{T})$: in this
case the instance of either TSEP or TFaEP is a no-instance, so we assume from
now on that $V(G_{T})=V(G)$. We claim that for each goal vertex $v_{g}\in
V(G)$, the unique temporal path from $v_{s}$ to $v_{g}$ in
$(G_{T},\lambda_{T})$ arrives at time $d_{G_{T}}(v_{s},v_{g})$, where
$d_{G_{T}}(v_{s},v_{g})$ is the distance from $v_{s}$ to $v_{g}$ in the graph
$G_{T}$.

We show this by induction on $d_{G_{T}}(v_{s},v_{g})$. If
$d_{G_{T}}(v_{s},v_{g})=1$, then $v_{s}v_{g}\in E(G_{T})$, and as
$1\in\lambda^{\prime}(v_{s}v_{g})$, there is a foremost path from $v_{s}$ to
$v_{g}$ that arrives at time 1. By definition of a ubiquitous foremost
temporal path tree, the unique temporal path from $v_{s}$ to $v_{g}$ must
therefore arrive at time 1, completing the base case for the induction.

We next proceed by induction. Let $d_{G_{T}}(v_{s},v_{g})=d>1$, and let
$v_{1}$ be a vertex such that $v_{1}v_{g}\in E(G_{T})$ and
$d_{G_{T}}(v_{s},v_{1})=d-1$. Such a vertex $v_{1}$ must exist, else there can
be no path of length $d$ from $v_{s}$ to $v_{g}$ in $G_{T}$. By induction, the
unique temporal path from $v_{s}$ to $v_{1}$ in $(G_{T},\lambda_{T})$ arrives
at time $d-1$. Additionally, there is a temporal path from $v_{s}$ to $v_{g}$
in $(G_{T},\lambda_{T})$ (else we would have $v_{g}\not\in V(G_{T})$. Let
$d^{\prime}$ be the arrival time of the unique foremost temporal path in
$(G_{T},\lambda_{T})$ from $v_{s}$ to $v_{g}$. As
$\lambda^{\prime}(v_{1}v_{g})=\\{1,2,\ldots,r_{e}\\}$ for some positive
integer $r_{e}$, if $d^{\prime}>d$ then $d^{\prime}\in\\{1,2,\ldots,r_{e}\\}$
and so we must have $d\in\\{1,2,\ldots,r_{e}\\}$ as both are positive
integers. But then the edge $v_{1}v_{g}$ could be used at time $d$ instead,
breaking our assumption. Therefore the unique temporal path from $v_{s}$ to
$v_{g}$ in $(G_{T},\lambda_{T})$ arrives at time $d$.

Noting that, as $(G_{T},\lambda_{T})$ is a ubiquitous foremost temporal path
tree, no temporal path in any $(\delta,\zeta)$-perturbation of $\lambda$ can
arrive earlier than $d_{G_{T}}(v_{s},v_{g})$, we can also deduce that no
temporal path in any $(\delta,\zeta)$-perturbation could be shorter than
$d_{G_{T}}(v_{s},v_{g})$ (else it could arrive earlier than
$d_{G_{T}}(v_{s},v_{g})$) and no temporal path in any
$(\delta,\zeta)$-perturbation can have a lower duration (by the same
argument). It therefore remains to check that either
$\operatorname{ecc^{S}}{((G_{T},\lambda_{T}),v_{s})}\leq k$ (if we are
considering TSEP) or $\operatorname{ecc^{Fa}}{((G_{T},\lambda_{T}),v_{s})}\leq
k$ (if we are considering TFaEP). As $(G_{T},\lambda_{T})$ is a ubiquitous
foremost temporal path tree, each of those can be checked in $O(n)$ time.

##  6 Conclusion and Open Questions

Motivated by the uncertainty and errors intrinsic in recording real-world
temporal graphs, we have investigated temporal graph reachability problems
under perturbation, in particular asking whether there exists a perturbation
of time-edges of a temporal graph that results in a minimum required
reachability. We find that when many perturbations are allowed, this can be
resolved efficiently, but that the problem is, in general, -hard and
$W[2]$-hard with respect to $\zeta$. Furthermore, under the assumption of ETH,
we show that there is no $f(\zeta)n^{o(\zeta)}$-time algorithm for our
problem. We have tractability results when the underlying graph is a tree or
when it has bounded treewidth and the lifetime and size of permitted
perturbations are also bounded by constants. For contrast, we show that
allowing many perturbations is not sufficient to give tractability of related
eccentricity problems, where the question is not just whether there is some
source vertex that reaches some specified number of vertices, but whether it
reaches vertices within some travel duration (fastest) or in few steps
(shortest).

An obvious avenue for further research is on other temporal graph problems
with the same style of perturbations. However, inspired by the motivation to
make temporal algorithms more suitable for real-world temporal contact data,
we suggest exploring temporal graph modification problems that are robust to
perturbations: for example, removing (by, e.g. vaccination) vertices from a
temporal network to limit reachability in a way that has guarantees under
possible bounds of perturbation.

## References

  * [1] Brenda Baker and Robert Shostak.  Gossips and telephones.  Discrete Mathematics, 2(3):191–193, June 1972.  URL: https://www.sciencedirect.com/science/article/pii/0012365X72900015, doi:10.1016/0012-365X(72)90001-5.
  * [2] Surender Baswana, Keerti Choudhary, and Liam Roditty.  Fault tolerant subgraph for single source reachability: generic and optimal.  In Proceedings of the forty-eighth annual ACM symposium on Theory of Computing, STOC ’16, pages 509–518, New York, NY, USA, June 2016\. Association for Computing Machinery.  URL: https://dl.acm.org/doi/10.1145/2897518.2897648, doi:10.1145/2897518.2897648.
  * [3] Davide Bilò, Gianlorenzo D’Angelo, Luciano Gualà, Stefano Leucci, and Mirko Rossi.  Blackout-Tolerant Temporal Spanners.  In Thomas Erlebach and Michael Segal, editors, Algorithmics of Wireless Networks, Lecture Notes in Computer Science, pages 31–44, Cham, 2022. Springer International Publishing.  doi:10.1007/978-3-031-22050-0_3.
  * [4] Hans L. Bodlaender.  A Linear-Time Algorithm for Finding Tree-Decompositions of Small Treewidth.  SIAM Journal on Computing, 25(6):1305–1317, December 1996.  Publisher: Society for Industrial and Applied Mathematics.  URL: https://epubs.siam.org/doi/10.1137/S0097539793251219, doi:10.1137/S0097539793251219.
  * [5] Binh-Minh Bui-Xuan, Afonso Ferreira, and Aubin Jarry.  Computing shortest, fastest, and foremost journeys in dynamic networks.  International Journal of Foundations of Computer Science, 14(02):267–285, April 2003.  doi:10.1142/s0129054103001728.
  * [6] Arnaud Casteigts.  A Journey through Dynamic Networks (with Excursions).  thesis, Université de Bordeaux, June 2018.  URL: https://hal.science/tel-01883384.
  * [7] Arnaud Casteigts, Swan Dubois, Franck Petit, and John M. Robson.  Robustness: A new form of heredity motivated by dynamic networks.  Theoretical Computer Science, 806:429–445, February 2020.  URL: https://www.sciencedirect.com/science/article/pii/S0304397519304979, doi:10.1016/j.tcs.2019.08.008.
  * [8] GL Chaters, PCD Johnson, Sarah Cleaveland, Joseph Crispell, WA De Glanville, Thomas Doherty, Louise Matthews, Sibylle Mohr, OM Nyasebwa, Gianluigi Rossi, et al.  Analysing livestock network data for infectious disease control: an argument for routine data collection in emerging economies.  Philosophical Transactions of the Royal Society B, 374(1776):20180264, 2019.
  * [9] Marek Cygan, Fedor V. Fomin, Łukasz Kowalik, Daniel Lokshtanov, Dániel Marx, Marcin Pilipczuk, Michał Pilipczuk, and Saket Saurabh.  Parameterized Algorithms.  Springer International Publishing, Cham, 2015.  URL: http://link.springer.com/10.1007/978-3-319-21275-3, doi:10.1007/978-3-319-21275-3.
  * [10] Argyrios Deligkas, Michelle Döring, Eduard Eiben, Tiger-Lily Goldsmith, and George Skretas.  Being an Influencer is Hard: The Complexity of Influence Maximization in Temporal Graphs with a Fixed Source, March 2023.  arXiv:2303.11703 [cs].  URL: http://arxiv.org/abs/2303.11703.
  * [11] Argyrios Deligkas, Eduard Eiben, and George Skretas.  Minimizing reachability times on temporal graphs via shifting labels.  In Proceedings of the Thirty-Second International Joint Conference on Artificial Intelligence, IJCAI ’23, pages 5333–5340, August 2023.  doi:10.24963/ijcai.2023/592.
  * [12] Argyrios Deligkas and Igor Potapov.  Optimizing reachability sets in temporal graphs by delaying.  Information and Computation, 285:104890, May 2022.  URL: https://www.sciencedirect.com/science/article/pii/S0890540122000323, doi:10.1016/j.ic.2022.104890.
  * [13] Michael Dinitz and Robert Krauthgamer.  Fault-tolerant spanners: better and simpler.  In Proceedings of the 30th annual ACM SIGACT-SIGOPS symposium on Principles of distributed computing, PODC ’11, pages 169–178, New York, NY, USA, June 2011. Association for Computing Machinery.  URL: https://dl.acm.org/doi/10.1145/1993806.1993830, doi:10.1145/1993806.1993830.
  * [14] Swan Dubois, Laurent Feuilloley, Franck Petit, and Mikaël Rabie.  When Should You Wait Before Updating? - Toward a Robustness Refinement.  In David Doty and Paul Spirakis, editors, 2nd Symposium on Algorithmic Foundations of Dynamic Networks (SAND 2023), volume 257 of Leibniz International Proceedings in Informatics (LIPIcs), pages 7:1–7:15, Dagstuhl, Germany, 2023. Schloss Dagstuhl – Leibniz-Zentrum für Informatik.  ISSN: 1868-8969.  URL: https://drops.dagstuhl.de/opus/volltexte/2023/17943, doi:10.4230/LIPIcs.SAND.2023.7.
  * [15] Krzysztof Dudziński and Stanisław Walukiewicz.  Exact methods for the knapsack problem and its generalizations.  European Journal of Operational Research, 28(1):3–21, January 1987\.  URL: https://www.sciencedirect.com/science/article/pii/0377221787901652, doi:10.1016/0377-2217(87)90165-2.
  * [16] Jessica Enright, Kitty Meeks, George B. Mertzios, and Viktor Zamaraev.  Deleting edges to restrict the size of an epidemic in temporal networks.  Journal of Computer and System Sciences, 119:60–77, August 2021\.  URL: https://www.sciencedirect.com/science/article/pii/S0022000021000155, doi:10.1016/j.jcss.2021.01.007.
  * [17] Jessica Enright, Kitty Meeks, and Fiona Skerman.  Assigning times to minimise reachability in temporal graphs.  Journal of Computer and System Sciences, 115:169–186, February 2021\.  URL: https://www.sciencedirect.com/science/article/pii/S0022000020300799, doi:10.1016/j.jcss.2020.08.001.
  * [18] Thomas Erlebach, Murilo Santos de Lima, Nicole Megow, and Jens Schlöter.  Learning-Augmented Query Policies for Minimum Spanning Tree with Uncertainty.  In Shiri Chechik, Gonzalo Navarro, Eva Rotenberg, and Grzegorz Herman, editors, 30th Annual European Symposium on Algorithms (ESA 2022), volume 244 of Leibniz International Proceedings in Informatics (LIPIcs), pages 49:1–49:18, Dagstuhl, Germany, 2022. Schloss Dagstuhl – Leibniz-Zentrum für Informatik.  ISSN: 1868-8969.  URL: https://drops.dagstuhl.de/opus/volltexte/2022/16987, doi:10.4230/LIPIcs.ESA.2022.49.
  * [19] Thomas Erlebach and Michael Hoffmann.  Minimum Spanning Tree Verification Under Uncertainty.  In Dieter Kratsch and Ioan Todinca, editors, Graph-Theoretic Concepts in Computer Science, Lecture Notes in Computer Science, pages 164–175, Cham, 2014. Springer International Publishing.  doi:10.1007/978-3-319-12340-0_14.
  * [20] Thomas Erlebach, Michael Hoffmann, and Murilo Santos de Lima.  Round-Competitive Algorithms for Uncertainty Problems with Parallel Queries.  Algorithmica, September 2022.  doi:10.1007/s00453-022-01035-6.
  * [21] Thomas Erlebach, Michael Hoffmann, Danny Krizanc, Matús Mihal’ák, and Rajeev Raman.  Computing Minimum Spanning Trees with Uncertainty, February 2008\.  arXiv:0802.2855 [cs].  URL: http://arxiv.org/abs/0802.2855, doi:10.48550/arXiv.0802.2855.
  * [22] Helen R Fielding, Trevelyan J McKinley, Matthew J Silk, Richard J Delahay, and Robbie A McDonald.  Contact chains of cattle farms in great britain.  Royal Society Open Science, 6(2):180719, 2019.
  * [23] Eugen Füchsle, Hendrik Molter, Rolf Niedermeier, and Malte Renken.  Delay-Robust Routes in Temporal Graphs, January 2022.  arXiv:2201.05390 [cs].  URL: http://arxiv.org/abs/2201.05390.
  * [24] Samuel D. Hand, Jessica Enright, and Kitty Meeks.  Making Life More Confusing for Firefighters, May 2022.  Conference Name: 11th International Conference on Fun with Algorithms (FUN 2022) ISBN: 9783959772327 ISSN: 1868-8969 Meeting Name: 11th International Conference on Fun with Algorithms (FUN 2022) Place: Sicily, Italy.  URL: https://eprints.gla.ac.uk/269008/.
  * [25] Petter Holme.  Network reachability of real-world contact sequences.  Physical Review E, 71(4):046119, 2005.
  * [26] David Kempe, Jon Kleinberg, and Amit Kumar.  Connectivity and inference problems for temporal networks.  Journal of Computer and System Sciences, 64(4):820–842, June 2002\.  doi:10.1006/jcss.2002.1829.
  * [27] David C. Kutner and Laura Larios-Jones.  Temporal Reachability Dominating Sets: Contagion in Temporal Graphs.  In Konstantinos Georgiou and Evangelos Kranakis, editors, Algorithmics of Wireless Networks, Lecture Notes in Computer Science, pages 101–116, Cham, 2023. Springer Nature Switzerland.  doi:10.1007/978-3-031-48882-5_8.
  * [28] William Lochet, Daniel Lokshtanov, Pranabendu Misra, Saket Saurabh, Roohani Sharma, and Meirav Zehavi.  Fault Tolerant Subgraphs with Applications in Kernelization.  In Thomas Vidick, editor, 11th Innovations in Theoretical Computer Science Conference (ITCS 2020), volume 151 of Leibniz International Proceedings in Informatics (LIPIcs), pages 47:1–47:22, Dagstuhl, Germany, 2020. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.  ISSN: 1868-8969.  URL: https://drops.dagstuhl.de/opus/volltexte/2020/11732, doi:10.4230/LIPIcs.ITCS.2020.47.
  * [29] George B. Mertzios, Othon Michail, and Paul G. Spirakis.  Temporal Network Optimization Subject to Connectivity Constraints.  Algorithmica, 81(4):1416–1449, April 2019.  doi:10.1007/s00453-018-0478-6.
  * [30] Hendrik Molter, Malte Renken, and Philipp Zschoche.  Temporal Reachability Minimization: Delaying vs. Deleting.  In Filippo Bonchi and Simon J. Puglisi, editors, 46th International Symposium on Mathematical Foundations of Computer Science (MFCS 2021), volume 202 of Leibniz International Proceedings in Informatics (LIPIcs), pages 76:1–76:15, Dagstuhl, Germany, 2021. Schloss Dagstuhl – Leibniz-Zentrum für Informatik.  ISSN: 1868-8969.  URL: https://drops.dagstuhl.de/opus/volltexte/2021/14516, doi:10.4230/LIPIcs.MFCS.2021.76.
  * [31] Ashleigh Myall, Robert L Peach, Yu Wan, Siddharth Mookerjee, Elita Jauneikaite, Frankie Bolt, James Price, Frances Davies, Andrea Y Weiße, Alison Holmes, et al.  Characterising contact in disease outbreaks via a network model of spatial-temporal proximity.  medRxiv, pages 2021–04, 2021.
  * [32] Ran Raz and Shmuel Safra.  A sub-constant error-probability low-degree test, and a sub-constant error-probability pcp characterization of np.  In Proceedings of the Twenty-Ninth Annual ACM Symposium on Theory of Computing, STOC ’97, page 475–484, New York, NY, USA, 1997. Association for Computing Machinery.  doi:10.1145/258533.258641.
  * [33] Steven S. Skiena.  The Algorithm Design Manual.  Springer Publishing Company, Incorporated, 2nd edition, 2008.
  * [34] Huanhuan Wu, James Cheng, Silu Huang, Yiping Ke, Yi Lu, and Yanyan Xu.  Path problems in temporal graphs.  Proc. VLDB Endow., 7(9):721–732, may 2014.  doi:10.14778/2732939.2732945.