# Main-belt and Trojan Asteroid Phase Curves from the ATLAS Survey

James E. Robinson,1 Alan Fitzsimmons,2 David R. Young,2 Michele Bannister,3
Larry Denneau,4 Nicolas Erasmus, 5 Amanda Lawrence,4 Robert J. Siverd,4 John
Tonry,4
1Institute for Astronomy, The University of Edinburgh, Blackford Hill,
Edinburgh EH9 3HJ, UK
2Astrophysics Research Centre, Queen’s University Belfast, University Road,
BT7 1NN, Belfast, Northern Ireland
3School of Physical and Chemical Sciences – Te Kura Matū, University of
Canterbury, New Zealand
4 Institute for Astronomy, University of Hawaii, Honolulu, HI 96822, USA
5 South African Astronomical Observatory, Cape Town 7925, South Africa E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Sparse and serendipitous asteroid photometry obtained by wide field surveys
such as the Asteroid Terrestrial-impact Last Alert System (ATLAS) is a
valuable resource for studying the properties of large numbers of small Solar
System bodies. We have gathered a large database of ATLAS photometry in
wideband optical cyan and orange filters, consisting of 9.6$\times 10^{7}$
observations of 4.5$\times 10^{5}$ main belt asteroids and Jupiter Trojans. We
conduct a phase curve analysis of these asteroids considering each apparition
separately, allowing us to accurately reject outlying observations and to
remove apparitions and asteroids not suitable for phase curve determination.
We obtain a dataset of absolute magnitudes and phase parameters for over
100,000 selected asteroids observed by ATLAS, $\sim 66,000$ of which had
sufficient measurements to derive colours in the ATLAS filters. To demonstrate
the power of our dataset we consider the properties of the Nysa-Polana
complex, for which the ATLAS colours and phase parameters trace the S-like and
C-like compositions amongst family members. We also compare the properties of
the leading and trailing groups of Jupiter Trojans, finding no significant
differences in their phase parameters or colours as measured by ATLAS,
supporting the consensus that these groups were captured from a common source
population during planetary migration. Furthermore, we identify $\sim 9000$
asteroids that exhibit large shifts in derived absolute magnitude between
apparitions, indicating that these objects have both elongated shapes and spin
axes with obliquity $\sim 90\
$\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}$$.

###### keywords:

minor planets, asteroids: general – astronomical data bases: miscellaneous

††pubyear: 2023††pagerange: Main-belt and Trojan Asteroid Phase Curves from
the ATLAS Survey–C

##  1 Introduction

There are currently over $10^{6}$ main belt asteroids (MBAs) located between
the orbits of Mars and Jupiter that have been discovered and catalogued
through astrometric observations. Once catalogued, their orbits are known and
their location within the Solar System relative to the Sun and an Earth-based
observer can be calculated for any given time. Measuring the apparent
magnitude, $m_{\lambda}$, through photometric methods at a specific point in
time then allows for an estimate of the absolute magnitude, $H_{\lambda}$, by
fitting those observations to the formula

$H_{\lambda}=m_{\lambda}-5\log\left(R\Delta\right)-g_{\lambda}(\alpha)$ (1)

Here $R$ and $\Delta$ are the heliocentric and geocentric distances in au and
$\alpha$ is the Sun-asteroid-Earth phase angle, all calculated at the time
$m_{\lambda}$ was measured. The phase function, $g_{\lambda}(\alpha)$,
mathematically describes the macroscopic scattering law of sunlight from the
asteroidal surface as a function of $\alpha$. For all models of asteroid light
scattering, $g_{\lambda}(\alpha)$ behaves such that increasing phase angles
leads to decreasing brightness in a monotonic fashion. Furthermore, as phase
angle approaches zero the brightness of asteroids often increases rapidly in
an opposition surge. The main drivers of this effect are believed to be the
mechanisms of shadow hiding and coherent backscatter (e.g. Muinonen et al.,
2002; Muinonen et al., 2010; Penttilä et al., 2016). Therefore by definition
the absolute magnitude, $H_{\lambda}$, of a spherical asteroid is its
magnitude when observed at $1\,$\mathrm{a}\mathrm{u}$$ from both the Sun and
the Earth, whilst at opposition ($\alpha=0^{\circ}$). In reality most
asteroids are non-spherical, and so the measured $H_{\lambda}$ will also
depend on its orientation as viewed by the observer.

Due to the large number of known asteroids and the commensurate difficulty of
obtaining physical or spectral information for individual objects, the
measurement of $H_{\lambda}$ and $g_{\lambda}(\alpha)$ often provides the
first estimate of asteroid properties for a large fraction of the population.
Even then, these parameters may have significant uncertainties; for instance
if the photometric measurements that are used to derive $g_{\lambda}(\alpha)$
only cover a small range of phase angles the phase function may not be well
constrained. Most photometric measurements are currently reported by large
wide-field sky surveys such as Pan-STARRS (Panoramic Survey Telescope and
Rapid Response System; Kaiser et al., 2010), the Catalina and Mount Lemmon Sky
Survey (Christensen et al., 2012), and ATLAS (Asteroid Terrestrial-impact Last
Alert System; Tonry et al., 2018). All of these surveys use different filters,
which means the reported $m_{\lambda}$ requires transformation to a common
standard filter, normally Johnson $V$-band for historical reasons.
Additionally, the wide range of spectral types of asteroids (DeMeo et al.,
2009) means the correct colour-dependent photometric transformation varies
from asteroid to asteroid, which increases the photometric uncertainty if the
spectral type is unknown.

Finally, the brightness of an asteroid will change over time due to shape-
driven rotational brightness variations, and it may also experience orbital-
timescale changes in the observed $H$ due to the ‘apparition effect’. This
effect can occur for asteroids with non-spherical shapes as the aspect angle
(the angle between the spin axis and the Earth-asteroid vector) changes,
meaning the asteroid projects different cross-sectional areas to the observer.
These changes in projected area lead to variations in the absolute magnitude
and lightcurve amplitude measured between apparitions (Fernández-Valenzuela,
2022). However, the influence of rotational and apparition effects can be
reduced by deriving $g_{\lambda}(\alpha)$ from data restricted to a single
apparition, sampling the full range of rotational phases. This ensures that
any brightness variations are primarily due to rotational effects and results
in approximately homogeneous scatter around the true $g_{\lambda}(\alpha)$.

The current asteroid measurement catalogues held at the Minor Planet Center
(MPC) and the JPL Small Bodies Database apply general filter transformations
from the reported photometric bandpass to the $V$-band for all large surveys,
but they often do not correct for the rotational/apparition effects listed
above. Given the $\sim 4\times 10^{8}$ observations in these databases from
different facilities with a range of photometric accuracy and bandpasses, this
is not surprising. Instead, absolute magnitudes $H_{V}$ are calculated by
transforming to the $V$-band and fitting the parameterised $(H,G)$ system of
(Bowell et al., 1989), where the form of $g_{\lambda}(\alpha)$ depends only on
a phase parameter $G$. $H_{V}$ is calculated for most asteroids using an
assumed value of G = 0.15; this approach is used by the MPC and astorb
(Moskovitz et al., 2022). Often object-specific values of $G$ are only used
for the subset of asteroids with highly accurate known phase curves.

The combined uncertainty of photometric accuracy and conversion between
bandpasses is important as fitting the phase curve leads directly to the
measurement of $H_{\lambda}$, which in turn allows determination of the
diameter $D$

$D_{\mathrm{km}}=C_{\lambda}\frac{10^{-H_{\lambda}/5}}{\sqrt{p_{\lambda}}}$
(2)

where $p_{\lambda}$ is geometric albedo and $C_{\lambda}$ is a wavelength
dependent constant related to the absolute magnitude of the Sun; e.g. for
$V$-band observations $C_{V}=664.5\,$\mathrm{k}\mathrm{m}$$ (Pravec & Harris,
2007). In the past it has been shown that significant biases can exist for
$H_{\lambda}$ (Jurić et al., 2002; Pravec et al., 2012). While this has been
presumed to be due to historical photometric calibration or transformation
issues, it may also be partly due to ignoring rotational and/or apparition
effects on phase curve accuracy.

Aside from absolute magnitude, $g_{\lambda}(\alpha)$ is also important in the
study of asteroid properties. The form of $g_{\lambda}(\alpha)$ has been found
to be related to both surface composition and asteroid taxonomy (Belskaya &
Shevchenko, 2000; Oszkiewicz et al., 2012), and also the surface scattering
properties such as regolith particle size and overall roughness of the surface
(Muinonen et al., 2002). Therefore there are two benefits of measuring
$g_{\lambda}(\alpha)$ beyond establishing $H_{\lambda}$. First, this can
indicate approximate compositions for asteroids without multi-colour
photometry or spectroscopy. Second, as surface composition is related to
albedo (Usui et al., 2013), this can significantly decrease the uncertainty in
the estimated diameter of an asteroid when using an assumed albedo.

Recognising the importance of information obtained by measuring the phase
curves of asteroids, such measurements of asteroid ensembles have been made by
several authors. By obtaining high-quality photometry of a small number of
objects, Harris & Young (1989) and Lagerkvist & Magnusson (1990) showed that
phase function and taxonomy were related. Although the observed correlations
could be due to either surface composition or albedo as these two parameters
are also related, Belskaya & Shevchenko (2000) found that albedo was the
primary factor influencing the phase curve. Albedo is generally the primary
factor in determining whether single scattering (for more absorbant dark
surfaces) or multiple scattering (for more reflective bright surfaces) of
reflected sunlight dominates the observed photometry (Muinonen et al., 2002).
Moving to large scale analysis of much bigger databases of asteroid
photometry, Oszkiewicz et al. (2011, 2012) debiased the existing MPC catalogue
for inter-observatory offsets in reported magnitudes, deriving absolute
magnitudes and phase functions for approximately $5\times 10^{5}$ asteroids.
Their work clearly showed correlations between the phase parameters and
position in the main belt, similar to those found by colour surveys, although
uncertainties caused by the inhomogenity of the source photometry remained.
Vereš et al. (2015) reported phase function fits for over $2\times 10^{5}$
asteroids from photometry obtained using sparse photometry from the wide-field
Pan-STARRS survey. They used a Monte-Carlo technique to estimate parameter
uncertainties resulting from the unknown rotational lightcurve modulation.
Significant uncertainties in the fitted parameters were reported due to this
and the sparseness of the dataset.

More recently there have been a number of studies focused on extracting
asteroid phase curve parameters from additional wide-field surveys. In
particular Mahlke et al. (2021), hereafter M21, performed a phase curve
analysis for 95,000 asteroids observed by the ATLAS sky survey (Tonry et al.,
2018). While similar to Pan-STARRS where data is obtained with consistent
reduction methods by a single facility (consisting of multiple units), the use
of only 2 filters with ATLAS (compared to 6 with Pan-STARRS), together with a
wider field of view per exposure, allows for a significantly higher cadence of
observations. This makes the ATLAS dataset an excellent resource for the study
of asteroid phase curves, even though its sensitivity is lower than other sky
surveys. Performing Bayesian parameter inference Markov-Chain Monte Carlo
analyses of the asteroid phase curves, M21 obtained absolute magnitudes and
phase curves with robustly estimated uncertainties. In this paper, we describe
an independent phase curve analysis of $9.6\times 10^{7}$ photometric
measurements of 4.5$\times 10^{5}$ main belt and Trojan asteroids obtained by
the ATLAS survey, from December 2015 to January 2022. By doing so we aim to
provide a large dataset of reliable phase curve parameters and absolute
magnitudes for the use the asteroid community. As all observations were
obtained using the same facility and filter system, we avoid the potential
problems in combining data from multiple sources. Our work builds upon the
recent study by M21 as we have access to a longer baseline of observations.
Furthermore we approach the problem in a different manner whereby we consider
the phase curves during individual apparitions for all asteroids, rather than
fitting a single phase curve to all observations.

In section 2 we describe the ATLAS photometric observations and our algorithm
to fit phase curves to individual apparitions, including the filtering out of
outlying observations and inaccurate phase curves. In section 3 we present the
resulting Selected ATLAS Phase-curves (SAP) dataset and its distributions of
absolute magnitude, phase parameters and colours. Furthermore we compare our
results to those of previous studies. Section 4 demonstrates the power of our
dataset by considering the phase curve properties of asteroid families and
Jupiter Trojans observed by ATLAS. We also discuss the limiting factors in our
methodology and the implications for future surveys. Finally in section 5 we
present a summary of the work, highlighting the major results.

##  2 Methodology

###  2.1 ATLAS observations

Until recently, the ATLAS survey was composed of two $0.5\,$\mathrm{m}$$
robotic telescopes located on Mauna Loa and Haleakala, Hawaii, with a
separation of approximately 140 kilometres. Two additional Southern hemisphere
stations have recently been added in Chile and South Africa, however the
analysis in this work consists exclusively of data obtained with the original
Northern hemisphere units. The main goal of ATLAS is to detect potential
terrestrial impactors in the Near-Earth Object (NEO) population, including
those on their final impact trajectory (Tonry et al., 2018). It does this by
taking images of the night sky on an automated schedule and subtracting a pre-
determined background static sky wallpaper, thus leaving behind only
photometric sources that differ from the background sky. Moving objects can be
clearly identified in a series of these difference images. Apart from
discovering asteroids, the difference imaging technique has also proved
extremely valuable for a range of transient astronomy science; e.g. variable
stars (Heinze et al., 2018), supernovae (Smith et al., 2020), gravitational
wave source identification (Smartt et al., 2017) and also very slow motion
distant Solar System bodies (Dobson et al., 2023).

Each ATLAS unit consists of a custom Wright-Schmidt design with
$0.5\,$\mathrm{m}$$ aperture and $0.65\,$\mathrm{m}$$ spherical primary mirror
(Tonry et al., 2018). The detector is an STA-1600 $10.5\times
10.5\,$\mathrm{k}$$ CCD with a pixel scale of $1.86\arcsec$. ATLAS primarily
surveys in an “orange” $o$ and “cyan” $c$-filter, reaching down to $\sim
19.5\,$\mathrm{m}\mathrm{a}\mathrm{g}$$ for its typical $30\,$\mathrm{s}$$
exposures. These filters are approximately equivalent to PS1 filters $r+i$
($560-820\,$\mathrm{n}\mathrm{m}$$) and $g+r$
($420-650\,$\mathrm{n}\mathrm{m}$$) respectively. Generally the bluer
$c$-filter is reserved for new moon dark time during survey operations,
whereas the redder $o$-filter is more frequently used during bright time
observations (Smith et al., 2020), although this strategy has evolved over the
course of the survey.

ATLAS units scan across the available night sky (with an accessible
declination range of $-45^{\circ}<\delta<+90^{\circ}$ for the Hawaii units),
revisiting each field 4 times over timescales $\leq 1$ hr. An individual unit
can cover approximately 1/4 of the visible sky per night in this manner,
therefore the two Hawaii units together cover the accessible Northern sky
every 2 days. The addition of units in South Africa and Chile opens up the
Southern sky to ATLAS, increasing all-sky coverage and cadence of field
revisits and allows ATLAS to approach $24\,$\mathrm{h}\mathrm{r}$$ night sky
monitoring for Potentially Hazardous Asteroids (PHAs) and reducing
restrictions on observations due to local weather. Telescope operations at
each unit are fully automated (for monitoring weather, dome control and
scheduling) and each unit will proceed with observations when safe and able to
do so.

The ATLAS pipeline automatically reduces images and performs sky subtraction
to search for transient sources via difference imaging. Photometry of detected
sources is measured from the original reduced and calibrated image (using
dophot, Schechter et al., 1993; Alonso-García et al., 2012). Moving object
detection and linking of these sources is performed using the Moving Object
Processing System (MOPS, Denneau et al., 2013). Candidate moving object
tracklets, consisting of $>3$ (typically 4) detections, are matched against
known MPC sources (via MPChecker 111https://www.minorplanetcenter.net/cgi-
bin/checkmp.cgi). Candidate NEOs are assigned an NEO digest metric score (Keys
et al., 2019) and are filtered by a deep learning classifier (Chyba Rabeendran
& Denneau, 2021) prior to screening by a human observer and submission to the
MPC.

ATLAS acts in synergy with the other major sky surveys, such as Pan-STARRS,
enhancing sky coverage for NEO discovery but also obtaining serendipitous
observations of any minor planets in the survey footprint. In addition to
prompt NEO astrometry, ATLAS tracklets that have been matched against the
ephemerides of known objects are submitted en masse to the MPC at the end of
every night. As such, in addition to its primary mission of NEO discovery and
characterisation, ATLAS is a valuable source of astrometry and photometry for
a significant number of Solar System objects.

###  2.2 The rockAtlas database

For this study we generated a database of ATLAS detections matched to known
Solar System minor planets. The rockAtlas database (Young, 2023) has been
processing all ATLAS observations for a fixed list of known asteroids since
2017. The initial list of asteroids to be tracked was $\sim 445,000$ asteroids
from the astorb 222https://asteroid.lowell.edu/main/astorb/ dataset (Moskovitz
et al., 2022) with well defined orbits at that time. Accurate orbital
parameters were required as rockAtlas calculates the ephemerides of these
objects and matches them to ATLAS observations. When rockAtlas is running, a
night’s worth of observations is downloaded from the ATLAS automated detection
and photometry pipeline - dophot (Tonry et al., 2018). dophot provides details
on all ATLAS detections that night, such as time, position, magnitude etc.
These detections consist of moving Solar System bodies, static transients and
false detections. rockAtlas calculates precise ephemerides for all asteroids
in its watch list at the mid-time of each exposure using pyephem
333https://rhodesmill.org/pyephem/ and orbfit
444http://adams.dm.unipi.it/orbfit/. These ephemerides are then cross-matched
to the dophot detection positions using a HEALpix grid (Gorski et al., 2005;
Zonca et al., 2019). A dophot detection is associated with a moving object,
provided the separation from the predicted ephemeris at the exposure mid-time
is less than
$5\,$\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}$$, and it is
added to the rockAtlas database. Together with the ephemerides, pyephem and
orbfit provide the viewing geometry for each observation including
heliocentric and geocentric distances, the Sun-asteroid-Earth phase angle and
galactic latitude, which are appended to the database.

The rockAtlas database provides an ever growing dataset of ATLAS asteroid
photometry with all the additional information needed for studying asteroid
phase curves. It contains $\sim 9.6\times 10^{7}$ observations as of January
2022, comprised of $\sim 65\%$ in the $o$-filter and $\sim 35\%$ in the
$c$-filter. Asteroids in the rockAtlas database have a median of $\sim 130$
observations. Apparent magnitudes range from $\gtrsim
12\,$\mathrm{m}\mathrm{a}\mathrm{g}$$ down to a limiting magnitude of $\sim
19.5$, with a median observed magnitude of $\sim
17\,$\mathrm{m}\mathrm{a}\mathrm{g}$$. Overall the median photometric
uncertainty is $0.05$ magnitudes, however this uncertainty increases
substantially at fainter magnitudes; the median photometric uncertainty is 0.1
magnitudes for $m_{o}>17\,$\mathrm{m}\mathrm{a}\mathrm{g}$$.

In this work we focus on the phase curve properties of asteroids from the main
asteroid belt out to the Jupiter Trojans, which dominate the rockAtlas dataset
(figure 1). Other studies have been performed of Near Earth Objects
($q<1.3\,$\mathrm{A}\mathrm{U}$$) and Outer Solar System Objects
($a>Q_{\text{Jup}}=5.457\,$\mathrm{A}\mathrm{U}$$) observed with ATLAS (Heinze
et al., 2021; Dobson et al., 2021, 2023). The median maximum phase angle of
observations is $\sim 20^{\circ}$, which is the expected value for asteroids
on circular orbits within the outer main belt. As expected and shown in figure
1 the maximum observed phase angle decreases with semimajor axis. The more
numerous Main Belt asteroids typically have observations up to a phase angle
of $\sim 20^{\circ}$ whereas the most distant objects we consider in this
work, the Jupiter Trojans, are observed only up to $\alpha\sim 10^{\circ}$.
Observations of asteroid brightness at low phase angles ($<5^{\circ}$) are
essential to characterise any opposition brightness surge that may be present
in the phase curve as the phase angle decreases to zero. The median minimum
observed phase angle for objects in the rockAtlas database is $\alpha\sim
3^{\circ}$, making this a valuable dataset for accurately investigating
asteroid phase curves. One may refer to table 1 for details of the
observational circumstances of a sample of objects in the database.

Recently, direct public access to ATLAS asteroid photometry has been made
available via the online SSCAT
database555https://astroportal.ifa.hawaii.edu/atlas/sscat/. This database is
constructed in a similar manner to rockAtlas (by matching dophot photometry to
asteroid ephemerides) and therefore it is recommended that the same
considerations addressed in this work (e.g. dealing with outlying photometry,
see section 2.4) should be applied to future work making use of the SSCAT
database.

Alternatively, the ATLAS forced photometry server (Shingles et al., 2021)
provides access to ATLAS observations for static and moving astrophysical
objects. However, photometry retrieved by forced measurement at requested
ephemerides (such as the method used by Dobson et al., 2023) may differ from
the automatically generated dophot measurements. In addition, the
computational overhead of performing forced photometry queries would most
likely limit the number of objects we wish to consider, therefore we did not
make use of this tool in this work. For example, a forced photometry query for
a typical MBA takes approximately 20 minutes to match frames and determine
magnitudes, whereas the same query through rockAtlas takes only seconds to
retrieve magnitudes, distances and phase angles as many of the parameters have
already been calculated and stored in the database.

Figure 1:  Two dimensional histogram showing the number distribution of the
maximum observed phase angle of asteroids in the rockAtlas database
($o$-filter) as a function of the object’s semimajor axis (left-hand panel).
The maximum phase angle for a zero-inclination circular orbit,
$\text{max}(\alpha)=\sin^{-1}{1/a}$, is shown as a dotted red line. The total
distribution of maximum phase angle for the whole dataset is shown as a
vertical histogram in the right-hand panel. The Main Asteroid Belt dominates
the database; these asteroids are typically observed up to a maximum phase
angle of $\gtrsim 20^{\circ}$.

###  2.3 Phase curve functions

There have been many phase curve models devised in differing degrees of
complexity and physical basis. For this work we concentrate on two models. The
first is the $(H,G)$ system published by Bowell et al. (1989), hereafter
denoted as B89. This model has two primary advantages. It is based on the
observed surface scattering properties of well observed asteroids. Also,
having been in use for over 3 decades, it has acted as the de-facto standard
representation of asteroid photometric properties over that period. Phase
curves are fitted to the reduced magnitude, $V(\alpha)$, which is the apparent
magnitude corrected for distance effects such that only the phase function
$g_{\lambda}(\alpha)$ from equation 1 leads to deviation from the absolute
magnitude

$V(\alpha)=m-5\log(R\Delta)$ (3)

In the B89 system the change in reduced magnitude with respect to phase angle
is described by

$V(\alpha)=H-2.5\log\left(\left[1-G\right]\Phi_{\text{B89},1}(\alpha)+G\Phi_{\text{B89},2}(\alpha)\right)$
(4)

where $H$ is the absolute magnitude of the object (brightness at zero phase
scaled to a distance of $1\,$\mathrm{a}\mathrm{u}$$), $\Phi_{\text{B89},i}$
are the basis functions, and $G$ is the phase parameter that controls the form
of the model phase curve.

The second model we use is the $(H,G_{12}^{*})$ system from Penttilä et al.
(2016), hereafter referred to as P16. This was derived from the new IAU
standard $(H,G_{1},G_{2})$ system formulated by Muinonen et al. (2010), which
more accurately reflects the phase function at small phase angles and is
defined as

$\begin{split}V(\alpha)=&H-2.5\log(G_{1}\Phi_{\text{P16},1}(\alpha)+G_{2}\Phi_{\text{P16},2}(\alpha)\\\
&+\left[1-G_{1}-G_{2}\right]\Phi_{\text{P16},3}(\alpha))\end{split}$ (5)

where $\Phi_{\text{P16},i}$ are the basis functions in this system. In
particular $\Phi_{\text{P16},3}$ models the opposition effect. There are two
phase parameters $G_{1}$, $G_{2}$ and in the Penttilä et al. (2016) model
these parameters are combined via the single $G_{12}^{*}$ phase parameter;
such that $G_{1}=0.84293649G_{12}^{*}$ and $G_{2}=0.53513350(1-G_{12}^{*})$.
The P16 model is better suited for modelling the phase curves of objects with
lower quality, more sparse observations. The phase parameter $G_{12}^{*}$ must
not be confused with the similar $G_{12}$ parameter which was also defined by
Muinonen et al. (2010). In this work where we consider mainly the P16 model,
we use the notation $G_{12}^{*}\rightarrow G_{12}$ for simplicity. Likewise,
from here on we use $H_{12}$ to distinguish the P16 absolute magnitude from
the B89 absolute magnitude $H$.

###  2.4 Phase curve fitting algorithm

Using the photometric observations and associated viewing geometries in the
rockAtlas database we have measured the phase curve behaviour of asteroids in
both the $o$ and $c$ filters. For a given asteroid, all associated ATLAS
observations, ephemerides and related quantities are retrieved from rockAtlas.
We used the sbpy package (Mommert et al., 2019) to fit various phase curve
models to these observations, namely $HG$ (B89) and $H_{12}G_{12}$ (P16).

We build upon previous studies of large numbers of asteroid phase curves by
accounting for the apparition effect in asteroid brightness. We define an
apparition as the period of time when an asteroid is continuously visible in
the sky to the ATLAS survey, i.e. when it is not too close to the sun for
observation. We therefore group observations according to their solar
elongation angle and consider each apparition of an asteroid separately when
fitting phase curves (similar to Schemel & Brown, 2021, hereafter SB21). The
long baseline of the ATLAS survey from 2017 to the present day means many
asteroids have been observed for up to 6 different apparitions and there are
asteroids observed by ATLAS which clearly exhibit the apparition effect as
previously noted by M21. Furthermore, the presence of such an effect implies
the asteroid has a non-spherical shape and a tilted spin pole and so we would
expect to observe measureable lightcurve variations. We can see the change in
lightcurve amplitude with apparition as the projected cross-sectional area
changes (Fernández-Valenzuela, 2022), e.g. for asteroid (766) Moguntia in
figure 2. Identifying such objects for further study with targeted lightcurve
observations and more sophisticated lightcurve inversion methods would be
extremely useful to determine accurate physical properties. For example, work
by Ďurech et al. (2020) has already demonstrated the benefits of combining
sparse ATLAS data (which covers a range of different viewing geometries to
help with shape determination) with dense lightcurve data (which helps
constrain the asteroid rotation) using convex lightcurve inversion (CLI)
methods.

Due to the automated nature of ATLAS dophot processing and the matching of
observations by rockAtlas, anomalous magnitude measurements exist in the
database. Such outlying data points can be caused by false matching of non-
asteroid sources, background star contamination, cosmic rays, stellar
diffraction spikes and ghost reflections from bright sources etc. (Heinze et
al., 2021). In the case of false matching, dophot photometric detections may
have been erroneously matched to an asteroid by the rockAtlas code, either due
to inaccurate ephemerides or the matching of a nearby stellar/galactic source
when the asteroid is fainter than predicted and is not detected (e.g. at its
rotational lightcurve minimum). Therefore the first step was to reject any
data points where the separation of the orbfit ephemerides to the dophot
position was
$\theta>1\,$\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}$$ to
reduce the number of false matches. In addition we rejected any observations
with low galactic latitude, $b\leq 10^{\circ}$, as near the galactic plane the
higher density of background stars may either contaminate the measured
brightness or lead to rockAtlas-dophot matching errors. We also rejected any
photometry that was performed on frames with a recorded 5-sigma limiting
magnitude of $<17.5$ (e.g. observations made near full moon). Finally we
rejected data points with extremely low (or zero) uncertainties in the dophot
magnitude, $\sigma_{m}\leq 0.005\,$\mathrm{m}\mathrm{a}\mathrm{g}$$, as we do
not trust such measurements given the typical photometric accuracy of ATLAS.
Out of this series of filters, the ephemeris position check removed the bulk
of obviously erroneous observations (25% of all detections had $\theta>1$
arcsec), therefore this tighter constraint is worth possibly rejecting the
small number of true matches with separations
$1<\theta($\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}$)<5$.

We then performed additional filtering of outlying data to each individual
apparition, considering the $o$ and $c$-filters separately. To each apparition
we fitted the B89 $HG$ phase curve model using a Levenberg-Marquadt non-linear
least squares algorithm (astropy; Astropy Collaboration et al., 2022). For
this first stage of phase curve fitting the data points were not weighted by
their photometric uncertainties, in order to reduce the influence of extreme
outliers in brightness with very low uncertainties. The initial estimate for
$H$ was taken from astorb and colour corrected from $V$-band to the ATLAS
filter of the observations. This conversion depends on an asteroid’s colour
and taxonomic class but we use the values from Heinze et al. (2021) which are
derived from the mean colour of S and C type asteroids.

$\displaystyle\langle V-o\rangle$ $\displaystyle=0.332$ (6)
$\displaystyle\langle V-c\rangle$ $\displaystyle=0.054$ (7)

The slope parameter $G$ was set at $G=0.15$ which is the standard assumed
value for all asteroids (Marsden, 1985). Using this $HG$ fit we rejected any
observations with a residual greater than 3 magnitudes to remove the most
extreme outliers. This limit is higher than the maximum lightcurve amplitude
for asteroids in the Asteroid Lightcurve Database (LCDB; Warner et al., 2021)
and hence only rejects photometry that is unrealistic for rotating inert
asteroids. We note that M21 used a residual limit of 1 magnitude when
filtering their ATLAS data. We then calculated the standard deviation of the
remaining residuals and performed a 2-sigma clip to remove further outliers.
By considering apparitions separately and fitting each for $H$ we account for
possible apparition effects which would increase the residuals if all
observations were considered together. If an asteroid displays large shifts in
$H$ between apparitions a phase curve fit to the full set of observations
would not reflect the true residuals at each apparition, therefore outlier
rejection could erroneously remove part or all of a particular apparition.

After cleaning up the observations for an asteroid we consider the properties
of each apparition and selected a “primary” epoch to use for the main phase
curve fit. We scored each apparition by the total number of observations, the
number of observations at phase angles $\alpha<5^{\circ}$ and by the range of
phase angle coverage. The scores for each of these parameters were scaled by
subtracting the mean and scaling to unit variance. The total score of an
apparition was taken to be the sum of all three scores (with equal weighting).
The primary apparition was the one which had the highest total score and so
the best combination of all three constraints required to get a reliable phase
curve fit. Simply selecting the apparition with the most observations does not
guarantee good coverage of the opposition effect at low phase angles and will
not accurately constrain the slope parameter and absolute magnitude. We only
selected primary apparitions in the $o$-filter as the ATLAS dataset
preferentially contains observations in $o$ rather than $c$.

The primary apparition was then fitted for absolute magnitude and phase curve
parameters with the selected phase curve model from sbpy, where the initial
values for $H$ are taken from astorb and colour corrected (as described in the
data cuts above). Our code fits all available sbpy models but we restrict our
analysis in this work to B89 and P16 for simplicity. The Levenberg-Marquardt
least squares algorithm was used to fit the phase curve model to the data, now
weighting each data point by its photometric uncertainty. The uncertainties in
the fitted phase curve parameters are taken to be the square root of the
diagonal of the covariance matrix. All other apparitions in $o$ and $c$ are
fitted for $H$ only, with the slope parameters kept fixed at the value
obtained from the primary apparition fit. Therefore we obtain the phase curve
slope parameter from the primary apparition alone and implicitly assume that
the phase curve is the same in $o$ and $c$ within our measurement
uncertainties. In our investigation we originally attempted full phase curve
fits (absolute magnitude and slope parameter) to all observations, and also to
each apparition. For most objects this can be a valid approach as oftentimes
outlying photometry is easily caught with the steps above and there is little
variation between apparitions. However, for some objects missed outliers and
apparitions with poor observational coverage can skew the overall phase curve
fit. We found that obtaining the slope parameter from the primary apparition
alone did not lead to significant changes for most objects, however, this
approach did succeed in removing the smaller number of objects with clearly
unphysical slope parameters.

Unlike previous studies we do not constrain slope parameters when fitting the
B89 and P16 models. Normally the slope parameters are restricted to “physical”
values such that the phase curve is always monotonically increasing in
brightness as phase angle decreases. Imposing such constraints can lead to a
distribution in phase parameter with pile ups near the boundaries (e.g. at 0
and 1 for $G_{12}$, M21, see their figure 5). These edge cases may tend
towards the boundaries because they genuinely have surfaces described by
extremely low/high values of phase parameter. Otherwise they may have
insufficient observational coverage to accurately constrain the phase
parameter, or very extreme rotational variations affecting the phase curve
fit. Significant changes in viewing geometry and/or observed aspect of an
asteroid can also lead to deviations from expected phase curve behaviour even
after correction for rotational effects (Jackson et al., 2022). Therefore to
catch these possible interesting cases and to prevent an artificial build up
of phase parameter at the boundaries we allow the fitter to go to these
“unphysical” values.

Importantly, we do not consider rotational variations in brightness in our
phase curve model, similar to many previous studies of asteroid phase curves
(e.g. McNeill et al., 2021; Mahlke et al., 2021). We assume that the variation
in reduced magnitude due to lightcurve amplitude will be reflected in the
formal uncertainties of the phase curve fit and in most cases this amplitude
is low compared to variation caused by photometric uncertainty. The cadence of
ATLAS data alone can in some cases be used to find unique rotation periods.
For example Ďurech et al. (2020) were able to determine rotation periods from
ATLAS data alone for $\sim 5000$ asteroids out of an initial sample of $\sim
180,000$ objects. Furthermore Erasmus et al. (2021) searched ATLAS
observations of $\sim 250,000$ asteroids for extreme slow rotators
($\mathrm{period}>1000\ $\mathrm{h}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{s}$$)
and identified 39 objects with high-confidence period solutions. However,
additional dense in time photometry can assist greatly in determining unique
rotation states and rejecting alias effects, as shown by Ďurech et al. (2020)
who used periods from the LCDB to significantly reduce the parameter space to
be searched when determining the shape and spin state from ATLAS data.

###  2.5 Filtering of phase curves

We implemented our phase curve fitting algorithm in a relatively agnostic
manner to the observational data, which has varying coverage and quality. As
such, we know that not all fits in our full phase curve database will be
reliable and we need to identify the high quality fits with well defined
parameters. To aid this, we recorded a variety of metrics for each phase curve
fit, including the number of data points in each fitted apparition, the phase
angle coverage and the statistics of the residuals for the fitted phase curve
model. We used the values of these metrics to select and reject phase curves
and we tuned the acceptable values for each metric by considering the final
distributions of $G$, and its uncertainty $\sigma_{G}$ for the B89 fits.
Visual inspection of individual phase curves showed that the most extreme (and
probably unphysical) values of $G$ and $\sigma_{G}$ were caused by poor phase
angle coverage or large photometric variations which meant the observations
could not accurately constrain the phase curve behaviour of an asteroid. The
final acceptable values for each metric were therefore chosen such that the
vast majority of objects had $0<G<1$, without us having to perform a hard cut
on $G$ at these boundaries. The details of this process are given in appendix
A but the main restrictions were such that the primary apparition had to have
at least 50 observations, with at least one observation at low phase angle,
$\alpha<5^{\circ}$, after outlier rejection. All other apparitions required at
least 10 observations to be considered reliable. After developing our metric
cuts for the $HG$ model, we then applied them in the same manner to the
$H_{12}G_{12}$ data. Our intention was to not introduce biases by over-tuning
the metric parameters for a particular model, and to allow for a less biased
comparison between the two models.

By ensuring sensible distributions of phase curve parameters and uncertainties
for the overall dataset, when we take a large subsample for a particular
dynamical class we can assume the statistics of that distribution will be a
good representation of the overall properties of the class. Poorly constrained
phase curve parameters for individual objects may still exist, but with large
enough numbers we can compare the overall trends for different populations in
a statistical manner.

###  2.6 Apparition effects, absolute magnitudes and asteroid colours in ATLAS
filters

By considering phase curve fits of individual apparitions, we are able to
identify asteroids with significant variations in absolute magnitude between
apparitions. In figure 2 we show an example of the apparition effect for
asteroid (766) Moguntia. The presence of this effect in ATLAS observations was
noted by M21, and additionally by SB21 who also fitted separate apparitions in
their study of the Jupiter Trojans using Zwicky Transient Facility (ZTF, Bellm
et al., 2019) data. For each asteroid we used the statistics of the measured
$H$ and the phase curve residuals to identify significant variations across
several apparitions. We define a parameter

$F_{\textrm{app}}=\frac{\sigma(H)}{{\rm
median}\left[\sigma(m_{\textrm{obs}}-m_{\textrm{fit}})\right]}$ (8)

which is the ratio of the standard deviation of $H$ across the retained
apparitions, divided by the median of the standard deviation of the observed
minus fitted phase curve residuals in each apparition. This ratio compares the
variation of $H$ across several apparitions to the typical brightness
variation after phase curve correction, which is dominated by rotational
lightcurve effects and photometric uncertainty. We consider just the
$o$-filter data in this test due to the larger number of photometric
measurements compared to the $c$-filter. When the changes in $H$ are larger
than the phase curve corrected brightness variations, we assume that the
object is undergoing apparition effects. By inspecting the phase curves of
individual objects we found that a boundary of $F_{\textrm{app}}\geq 0.9$
succeeds in identifying the objects with the strongest apparition effects. The
remaining asteroids with $F_{\textrm{app}}<0.9$ do not present a significant
apparition effects compared to the rotational variation and ATLAS photometric
accuracy. Also, we note that 28% of the asteroids have only one suitable
apparition after the quality cuts and so cannot be searched for apparition
effects. Using the $F_{\textrm{app}}\geq 0.9$ metric, we found that 9% of our
asteroids displayed a significant apparition effect, where this is clearly a
lower limit to the true number that exhibit this phenomenon.

For objects without strong apparition effects the final value of absolute
magnitude is estimated as the median $H$ of all retained apparitions in either
the $o$ or $c$-filters; figure 3 shows this for typical MBA (4986) Osipovia.
We used the median to reduce the influence of possible outlying apparitions
that might skew the mean. For objects with strong apparition effects we should
consider the minimum value of $H$, i.e. the brightest estimate is closest to
the true size. However this value is only an upper limit on the true absolute
magnitude (i.e. a lower limit on equivalent size) and such objects require
further investigation to accurately determine their physical parameters. As
such we generally exclude the objects with strong apparition effects
($F_{\text{app}}\geq 0.9$) from our results when discussing absolute magnitude
and colour, although we make use of the phase parameters $G$, $G_{12}$ of
these objects when $H$ is not considered. For all objects in our dataset the
phase parameter is determined solely from the phase curve fit of the primary
apparition and so can be compared regardless of apparition effects.

Most of our analysis thus far has focused on the more frequently made
observations in the ATLAS $o$-filter. The majority of asteroids observed by
ATLAS also have a smaller number of measurements in the bluer $c$-filter. Thus
we can compare the estimates of $H$ in each filter to estimate a $H_{c}-H_{o}$
colour for a large sample of asteroids. Measurements in the $c$-filter suffer
from the same systematic errors discussed above for the $o$-filter. Therefore
we calculated colour only for objects with apparitions that satisfy the
quality control requirements in both $o$ and $c$-filters, excluding any
objects displaying strong apparition effects in the $o$-filter
($F_{\textrm{app}}\geq 0.9$) as we cannot define a unique $H$ for these
objects. The colour was then calculated as median($H_{c}$) - median($H_{o}$).
The ATLAS colour provides an additional diagnostic for comparing the surface
properties of different asteroid populations.

Figure 2:  ATLAS phase curves of asteroid (766) Moguntia, an object displaying
a large apparition effect. Observations in the $o$-filter are shown in the
left-hand panel. Each apparition has been colour coded and outlying
photometry/rejected apparitions are shown as small black markers (some of
which are obscured by the figure legend). The limits of the $y$-axis have been
set to focus on retained photometry, ignoring any more extreme outliers. The
corresponding $HG$ phase curve fit for each apparition is shown by the
matching coloured curves, with the primary apparition indicated by the
coloured star. The median $H$ of the retained apparitions is indicated by the
dashed black curve with parameters: $H_{o}=9.68$, $G=0.21$. Observations in
the $c$-filter are shown in the right-hand panel, where marker/line colour
have the same meanings as the left panel. The large shifts in $H$ between
apparitions are clear ($F_{\text{app}}=5.9$; equation 8).  Figure 3:  Phase
curves for (4986) Osipovia, a typical MBA observed by ATLAS. As in figure 2
the left and right-hand panels show the $o$ and $c$-filter observations
respectively. Marker/line colour indicates the apparition to which each
observation belongs, with the primary apparition marker by the coloured star
and the median value of $H$ indicated by the black dashed curve with
parameters: $H_{o}=13.21$, $G=0.30$

##  3 RESULTS

###  3.1 ATLAS phase curve database

Name | Number | $N_{o}$ | $N_{c}$ | max($t$) (mjd) | $N_{\textrm{app}}$ | max($\alpha_{o}$) | max($\alpha_{c}$) | $\Delta\alpha_{o}$ | $\Delta\alpha_{c}$ |  $a$ (au) | $e$ | $i(^{\circ})$
---|---|---|---|---|---|---|---|---|---|---|---|---
Moguntia | 766 | 1029 | 363 | 59432.3 | 5 | 21.24 | 21.22 | 19.68 | 19.73 | 3.02 | 0.10 | 10.090543
Priamus | 884 | 1297 | 420 | 59523.2 | 6 | 12.60 | 12.59 | 11.53 | 12.26 | 5.20 | 0.13 | 8.906102
Marietta | 2144 | 1309 | 433 | 59600.2 | 5 | 21.55 | 21.36 | 21.36 | 20.44 | 2.87 | 0.06 | 2.833143
Yatskiv | 2728 | 1192 | 335 | 59600.2 | 5 | 29.29 | 29.21 | 28.32 | 28.02 | 2.46 | 0.17 | 2.610173
Osipovia | 4986 | 652 | 246 | 59596.5 | 5 | 27.46 | 27.37 | 25.76 | 24.92 | 2.56 | 0.21 | 5.611378
1991 EN | 9857 | 1065 | 314 | 59598.3 | 6 | 11.57 | 11.45 | 10.68 | 9.28 | 5.16 | 0.03 | 19.612052
Table 1: Overview of the observations made by ATLAS and collected by rockAtlas, for a selected sample of 6 asteroids (this table will be available in its entirety online). We provide the name and MPC number of each asteroid, alongside the number of detections in the $o$ and $c$-filters ($N_{o}$, $N_{c}$) and the date of the last observation used in this work, $\text{max}(t)$. The maximum phase angle observed in each filter is given by $\text{max}(\alpha_{o})$ and $\text{max}(\alpha_{c})$, and the maximum minus minimum phase angle range for each filter is $\Delta\alpha_{o}$, $\Delta\alpha_{c}$. Finally, the orbital elements obtained from astorb are also included; semimajor axis $a$, eccentricity $e$ and inclination $i$.  Name |  Number |  $n_{\textrm{app}}$ |  $t_{\textrm{app},0}$ (mjd) |  max($\alpha_{\textrm{app},o}$) |  max($\alpha_{\textrm{app},c}$) |  $N_{\textrm{fit},o}$ |  $N_{\textrm{fit},c}$ |  $H_{o}$ |  $H_{c}$ |  $\sigma_{Ho}$ |  $\sigma_{Hc}$ |  $G$ |  $\sigma_{G}$ |  Primary Apparition |  Keepo |  Keepc-o
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Moguntia |  766 |  0 |  57380.7 |  17.90 |  18.00 |  47 |  59 |  9.28 |  9.66 |  6.93e-03 |  3.65e-03 |  0.21 |  NaN |  0 |  1 |  1
Moguntia |  766 |  1 |  57883.6 |  18.09 |  17.34 |  146 |  20 |  10.09 |  10.67 |  8.86e-03 |  2.34e-02 |  0.21 |  NaN |  0 |  1 |  0
Moguntia |  766 |  2 |  58298.6 |  20.52 |  20.53 |  199 |  49 |  9.25 |  9.61 |  2.00e-03 |  4.54e-03 |  0.21 |  NaN |  0 |  0 |  0
Moguntia |  766 |  3 |  58733.6 |  21.24 |  21.22 |  232 |  70 |  10.15 |  10.51 |  8.96e-03 |  1.61e-02 |  0.21 |  NaN |  0 |  1 |  1
Moguntia |  766 |  4 |  59190.6 |  18.29 |  18.46 |  244 |  138 |  9.27 |  9.63 |  2.42e-03 |  1.22e-03 |  0.21 |  2.87e-03 |  1 |  1 |  1
Priamus |  884 |  0 |  57384.4 |  9.92 |  9.94 |  41 |  74 |  8.59 |  8.99 |  1.93e-02 |  1.54e-02 |  0.23 |  NaN |  0 |  0 |  0
Priamus |  884 |  1 |  57673.6 |  10.30 |  10.30 |  145 |  88 |  8.43 |  8.86 |  5.70e-03 |  1.13e-02 |  0.23 |  NaN |  0 |  1 |  1
Priamus |  884 |  2 |  58037.6 |  10.94 |  10.78 |  179 |  22 |  8.39 |  8.72 |  3.90e-03 |  1.14e-02 |  0.23 |  NaN |  0 |  0 |  0
Priamus |  884 |  3 |  58460.6 |  11.82 |  11.82 |  292 |  35 |  8.39 |  8.75 |  9.36e-03 |  9.19e-03 |  0.23 |  1.48e-02 |  1 |  1 |  1
Priamus |  884 |  4 |  58867.6 |  12.60 |  12.59 |  200 |  80 |  8.50 |  8.86 |  7.11e-03 |  1.27e-02 |  0.23 |  NaN |  0 |  1 |  1
Priamus |  884 |  5 |  59255.6 |  12.05 |  12.37 |  20 |  14 |  8.55 |  8.96 |  2.54e-02 |  1.50e-02 |  0.23 |  NaN |  0 |  0 |  0
Marietta |  2144 |  0 |  57481.6 |  20.26 |  20.39 |  45 |  112 |  10.99 |  11.39 |  2.12e-02 |  1.33e-02 |  0.21 |  NaN |  0 |  1 |  1
Marietta |  2144 |  1 |  57943.6 |  21.55 |  21.36 |  175 |  27 |  11.06 |  11.39 |  1.14e-02 |  2.48e-02 |  0.21 |  NaN |  0 |  1 |  0
Marietta |  2144 |  2 |  58417.6 |  20.14 |  20.28 |  256 |  52 |  10.96 |  11.42 |  8.33e-03 |  1.77e-02 |  0.21 |  NaN |  0 |  1 |  1
Marietta |  2144 |  3 |  58846.6 |  19.45 |  19.45 |  309 |  126 |  11.08 |  11.41 |  8.92e-03 |  1.27e-02 |  0.21 |  NaN |  0 |  1 |  1
Marietta |  2144 |  4 |  59300.6 |  20.58 |  20.47 |  410 |  90 |  11.02 |  11.41 |  1.28e-02 |  1.39e-02 |  0.21 |  1.58e-02 |  1 |  1 |  0
Yatskiv |  2728 |  0 |  57400.5 |  22.19 |  21.77 |  49 |  75 |  12.42 |  12.56 |  1.66e-02 |  2.09e-02 |  0.11 |  NaN |  0 |  1 |  1
Yatskiv |  2728 |  1 |  57826.6 |  28.61 |  27.77 |  113 |  13 |  12.29 |  12.59 |  1.07e-02 |  1.65e-02 |  0.11 |  NaN |  0 |  1 |  0
Yatskiv |  2728 |  2 |  58385.6 |  20.64 |  19.43 |  131 |  28 |  12.30 |  12.65 |  1.23e-02 |  3.43e-02 |  0.11 |  NaN |  0 |  1 |  0
Yatskiv |  2728 |  3 |  58790.6 |  23.23 |  23.01 |  231 |  66 |  12.50 |  12.63 |  1.24e-02 |  1.82e-02 |  0.11 |  NaN |  0 |  1 |  1
Yatskiv |  2728 |  4 |  59276.6 |  29.29 |  29.21 |  404 |  85 |  12.31 |  12.56 |  1.55e-02 |  1.09e-02 |  0.11 |  1.35e-02 |  1 |  1 |  1
Osipovia |  4986 |  0 |  57591.6 |  27.46 |  27.37 |  48 |  46 |  13.21 |  13.63 |  1.53e-02 |  1.65e-02 |  0.30 |  NaN |  0 |  1 |  1
Osipovia |  4986 |  1 |  58037.6 |  24.71 |  23.86 |  158 |  37 |  13.23 |  13.63 |  1.01e-02 |  2.07e-02 |  0.30 |  NaN |  0 |  1 |  1
Osipovia |  4986 |  2 |  58515.7 |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  0 |  0
Osipovia |  4986 |  3 |  58974.6 |  25.95 |  25.79 |  191 |  96 |  13.19 |  13.58 |  1.82e-02 |  9.53e-03 |  0.30 |  2.31e-02 |  1 |  1 |  1
Osipovia |  4986 |  4 |  59493.6 |  26.30 |  24.83 |  110 |  35 |  13.22 |  13.60 |  1.10e-02 |  2.15e-02 |  0.30 |  NaN |  0 |  1 |  0
1991 EN |  9857 |  0 |  57461.6 |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  NaN |  0 |  0
1991 EN |  9857 |  1 |  57882.6 |  9.72 |  8.05 |  36 |  7 |  10.04 |  10.48 |  2.03e-02 |  6.54e-02 |  0.16 |  NaN |  0 |  0 |  0
1991 EN |  9857 |  2 |  58269.5 |  11.38 |  11.37 |  140 |  45 |  10.16 |  10.49 |  1.03e-02 |  1.41e-02 |  0.16 |  NaN |  0 |  1 |  1
1991 EN |  9857 |  3 |  58613.6 |  11.57 |  11.45 |  255 |  66 |  10.09 |  10.40 |  4.74e-03 |  1.36e-02 |  0.16 |  NaN |  0 |  1 |  1
1991 EN |  9857 |  4 |  58997.6 |  11.40 |  11.39 |  209 |  109 |  10.09 |  10.44 |  3.96e-03 |  7.08e-03 |  0.16 |  NaN |  0 |  1 |  1
1991 EN |  9857 |  5 |  59414.6 |  11.17 |  10.99 |  352 |  46 |  10.10 |  10.46 |  1.36e-02 |  1.34e-02 |  0.16 |  2.14e-02 |  1 |  1 |  1
Table 2: Results of the phase curve fitting algorithm to determine phase curves for each apparition of the selected objects from table 1. NaN entries indicate where a phase curve fit could not be determined (e.g. the apparition was within the galactic plane and all observations were rejected). We assign an index to each apparition for an object ($n_{\text{app}}$), where the first observation in an apparition was made at time $t_{\text{app},0}$. The maximum phase angle observed in an apparition is given by $\text{max}(\alpha_{\text{app},o})$ and $\text{max}(\alpha_{\text{app},c})$ for the $o$ and $c$-filters respectively. The number of detections used to fit the phase curve (after outlier rejection) is denoted by $N_{\text{fit},o}$ and $N_{\text{fit},c}$. We show the fitted phase curve parameters for the $HG$ B89 system, alongside their uncertainties ($\sigma_{H}$, $\sigma_{G}$) in both $o$ and $c$-filters. The “Primary Apparition” flag indicates which apparition was fitted for both absolute magnitude and phase parameter. We also flag whether an apparition passed the quality control checks (sections 2.4 & 2.5) in the $o$ filter (“Keepo”) or if it passed in both $o$ and $c$ (“Keepc-o”). Only selected columns are displayed here. The table of apparitions will be available in its entirety online, with additional columns describing in detail the phase curve fit to each apparition for both the B89 and P16 models.  Name | Number | $N_{\textrm{app},o}$ | $F_{\textrm{app}}$ | $H_{o}$ | $H_{c}$ | $\sigma_{Ho}$ | $\sigma_{Hc}$ | $G$ | $\sigma_{G}$ | $N_{\textrm{app,}c,o}$ | $H_{c}-H_{o}$
---|---|---|---|---|---|---|---|---|---|---|---
Moguntia | 766 | 4 | 5.90 | 9.68 | 10.08 | 2.42e-03 | 1.22e-03 | 0.21 | 2.87e-03 | 3 | 0.38
Priamus | 884 | 3 | 0.79 | 8.43 | 8.86 | 9.36e-03 | 9.19e-03 | 0.23 | 1.48e-02 | 3 | 0.43
Marietta | 2144 | 5 | 0.30 | 11.02 | 11.41 | 1.28e-02 | 1.39e-02 | 0.21 | 1.58e-02 | 3 | 0.43
Yatskiv | 2728 | 5 | 0.64 | 12.31 | 12.59 | 1.55e-02 | 1.09e-02 | 0.11 | 1.35e-02 | 3 | 0.15
Osipovia | 4986 | 4 | 0.14 | 13.21 | 13.62 | 1.82e-02 | 9.53e-03 | 0.30 | 2.31e-02 | 3 | 0.42
1991 EN | 9857 | 4 | 0.37 | 10.10 | 10.45 | 1.36e-02 | 1.34e-02 | 0.16 | 2.14e-02 | 4 | 0.35
Table 3: The final B89 phase curve parameters for the selected objects in tables 1, as determined from the individual apparitions in table 2. Column $N_{\text{app},o}$ provides the number of apparitions in the $o$ filter that were used and $F_{\text{app}}$ is the parameter measuring the strength of the apparition effect (equation 8). The overall absolute magnitude, phase parameter and the associated uncertainties for each object are presented as calculated using the methods in section 2.6. We also provide the number of apparitions that had high quality phase curves in both filters ($N_{\text{app},c,o}$), for which we could calculate an asteroid colour $H_{c}-H_{o}$. Available online.  Name | Number | $N_{\textrm{app},o}$ | $F_{\textrm{app}}$ | $H_{12o}$ | $H_{12c}$ | $\sigma_{H12o}$ | $\sigma_{H12c}$ | $G_{12}$ | $\sigma_{G12}$ | $N_{\textrm{app,}c,o}$ | $H_{12c}-H_{12o}$
---|---|---|---|---|---|---|---|---|---|---|---
Moguntia | 766 | 4 | 5.87 | 9.63 | 10.02 | 1.19e-03 | 1.28e-03 | 0.43 | 8.52e-03 | 2 | 0.37
Priamus | 884 | 3 | 0.74 | 8.40 | 8.83 | 1.01e-02 | 9.45e-03 | 0.69 | 6.22e-02 | 3 | 0.43
Marietta | 2144 | 4 | 0.35 | 10.96 | 11.37 | 7.70e-03 | 1.33e-02 | 0.63 | 6.01e-02 | 4 | 0.41
Yatskiv | 2728 | 5 | 0.62 | 12.36 | 12.63 | 6.24e-03 | 1.07e-02 | 0.73 | 4.17e-02 | 2 | 0.18
Osipovia | 4986 | 4 | 0.14 | 13.07 | 13.47 | 8.80e-03 | 9.46e-03 | 0.18 | 5.59e-02 | 3 | 0.42
1991 EN | 9857 | 4 | 0.37 | 10.12 | 10.47 | 1.49e-02 | 1.35e-02 | 0.74 | 8.47e-02 | 4 | 0.35
Table 4: The final $H_{12}G_{12}$ P16 phase curve parameters for the selected
objects from table 1, as determined from the individual apparitions in table
2. These columns have the same meanings as in table 3, but now for P16 phase
curves as opposed to B89. Available online.

The result of the analysis above are two databases. We have a database of
phase curve fits by apparition derived for 427,662 MBA and Jupiter Trojan
asteroids with ATLAS observations gathered by rockAtlas. A sample of these
results is presented in table 2. We also have the Selected ATLAS Phase-curves
(SAP) database containing the subset of filtered high quality fits for 106,055
unique asteroids; 93,373 and 91,103 in the B89 and P16 systems respectively).
Approximately 9% of these objects display strong apparition effects; 8450
(B89) and 8131 (P16). Samples of these datasets are provided in tables 3 and
4. In addition, we obtained $H_{c}-H_{o}$ colours for objects with suitable
absolute magnitudes in both filters, resulting in colours for 66143 asteroids
in the B89 system (64808 for P16). Please see the Data Availability statement
at the end of this manuscript.

---
Figure 4:  Log-scale distributions of $H$ (left), $G$ (middle) and
$H_{c}-H_{o}$ colour (right) for ATLAS phase curve fits that have passed the
quality control checks (section 2.5) and are included in the SAP dataset
(section 3.1). $H$ and $G$ are given in the $o$-filter and B89 model. In both
figures we show the underlying distributions of S and C type asteroids within
the ATLAS dataset (taxonomy from SsODNet; Berthier et al., 2023). The median
values of the S and C-type distributions are indicated by the vertical line of
matching colour and their values are as follows: $H_{o}(\mathrm{S})=14.83$,
$H_{o}(\mathrm{C})=14.52$, $G(\mathrm{S})=0.29$, $G(\mathrm{C})=0.21$,
$H_{c}-H_{o}(\mathrm{S})=0.39$, $H_{c}-H_{o}(\mathrm{C})=0.27$. In the middle
panel we indicate the boundaries on $G$ for a monotonically increasing B89
phase curve with vertical red dashed lines.

In figure 4 we present the distributions of $H_{o}$, $G$ and $H_{c}-H_{o}$
colour for asteroids observed by ATLAS in the SAP database. We also show the
distributions of asteroids that have either an S-type or C-type taxonomic
classification from the Solar system Open Database Network (SsODNet; Berthier
et al., 2023), as these are the most common asteroid compositional classes.
While the overall $H$ distribution is smooth with a single dominant peak, the
underlying S and C distributions peak at different values. This is an
observational bias; the S-types have higher albedos than the darker C-types
and are found more commonly in the inner asteroid belt (Usui et al., 2013;
DeMeo & Carry, 2014). Therefore the S-types are generally brighter than C
types and they are better represented in ATLAS observations down to fainter
absolute magnitudes.

The selection procedure made in section 2.5 has successfully constrained the
majority of the $G$ distribution to “physical” solutions between 0 and 1. Only
a small number of phase curves lie outside these limits ($0.3\%$ of the
dataset) however we did not wish to arbitrarily exclude them or other phase
curves near the boundaries from our results. Visual inspection of these
particular objects show that the fitted phase curves are consistent with the
observations, albeit generally for a more limited range of phase angles. For
$G$ outside the $(0,1)$ range it is usually not until higher phase angles than
those observed by ATLAS that the phase function becomes non-
monotonic/asymptotic.

As $G$ decreases the slope of the approximately linear component at phase
angles $\alpha\gtrsim 5^{\circ}$ gets steeper. In particular the outliers with
$G\lesssim 1$ have phase curves that are entirely compatible with the
observations. On the other hand, as $G$ increases the slope at larger phase
angles decreases and the model fails to describe the opposition effect at
phase angles $<5\,$\mathrm{d}\mathrm{e}\mathrm{g}$$. Our handful of objects
with $G\gtrsim 1$ have relatively flat phase curves over the observed phase
angle range and with rotational variation this allows for an unphysical phase
curve to be fitted, where brightness decreases at opposition. However, the
vast majority of objects in our dataset have sensible phase curve fits and we
can see a mildly bimodal distribution in $G$ caused by the underlying
distributions of S and C-type phase parameters. This bimodality primarily
arises due to the different albedos (and compositions) of the two taxonomic
classes; the brighter S-type asteroids have a mean geometric albedo of 0.21
whereas the dark C-types have typical albedos of 0.07 (Usui et al., 2013). As
shown by Belskaya & Shevchenko (2000) the gradient of the phase curve at large
phase angles is correlated with albedo. In the B89 model, the higher median
value of $G$ for the S-type asteroids leads to a flatter phase function, which
corresponds to a higher albedo as expected.

There is a correlation between $H$ and the range of $G$ we find in the data
(figure 5). This is partly due to the enhanced presence of S-type asteroids as
$H$ increases, as we are most sensitive to detecting small asteroids in the
inner main belt where S-types are most populous, and these objects typically
have higher values of $G$. However the main cause is a systematic effect; as
$H$ increases asteroid size is decreasing and we are generally considering
fainter photometric observations with higher uncertainty. This increased
uncertainty means there is more scatter in fitted phase parameters for fainter
objects with larger $H$. This $H$-$G$ correlation driven by photometric
uncertainty acts in addition to any uncertainty caused by rotational
lightcurve effects. For example, analysis of the $HG$ dataset published by
SB21 shows a similar trend.

Figure 5:  Two dimensional histogram which shows the variation of $G$ with
absolute magnitude $H$ for asteroids in the SAP dataset. There is a positive
correlation between $H$ and $G$, with an increasing dispersion in $G$ for
fainter objects (increasing $H$).

Figure 4 shows that there is a strong bimodality for the $H_{c}-H_{o}$ colour
distribution in the SAP database. As described in section 2.6 the data
presented is for the subset of objects that have observations of sufficient
quality in both ATLAS filters for one or more apparitions (with
$F_{\text{app}}<0.9$). The bimodality is again driven by the different surface
compositions of the dominant S and C type asteroids. We obtain median
$H_{c}-H_{o}$ colours of $0.39\pm 0.05$, $0.27\pm 0.06$ (uncertainties of one
standard deviation) for the objects in our dataset matched with S and C-type
classifications respectively (Berthier et al., 2023). These values are
compatible with the expected $H_{c}-H_{o}$ colours calculated by Erasmus et
al. (2020), who convolved the mean Bus-DeMeo spectra (DeMeo et al., 2009) with
the ATLAS filter response to obtain $H_{c}-H_{o}$ colours of 0.388 and 0.249
for S and C-type asteroids respectively. Figure 6 shows the 2-D bimodal
distribution of $H_{c}-H_{o}$ vs. $G$ in the SAP database. Taken together, the
phase parameter and colour can be extremely valuable in distinguishing between
the dominant S and C-type asteroid compositions (illustrated further in
section 4.2).

Figure 6:  Two dimensional histogram of $H_{c}-H_{o}$ colour vs. $G$. The
bimodal distribution of matched S and C-type asteroids (SsODNet; Berthier et
al., 2023) is highlighted by orange and green contours respectively, which
indicates the 25% level with respect to the largest histogram bin value of
their sub-distributions in $H_{c}-H_{o}$ vs. $G$ parameter space. The median
values of the S and C-types with ATLAS colours are $G=0.29,0.21$ and
$H_{c}-H_{o}=0.39,0.27\,$\mathrm{m}\mathrm{a}\mathrm{g}$$ respectively.

We have also derived the phase curve parameters of asteroids using the P16
$H_{12}G_{12}$ system, where we use the notation $H_{12}$ to denote absolute
magnitudes in this system. As described in section 2.5 the quality control
cuts were developed on the $HG$ dataset and then applied separately to the
residuals of the $H_{12}G_{12}$ phase curve fits. Therefore not all asteroids
are present in both datasets. Regardless, both phase curve models show broadly
similar results. As before the phase parameter and colour distributions are
both bimodal due to the dominant S and C-type asteroids, the former shown in
figure 7. The quality control metrics developed for $HG$ do not lead to as
strong restrictions on $G_{12}$ for “physical” values between 0 and 1. This
appears to be mostly due to inherent differences between the $HG$ and
$H_{12}G_{12}$ models. Figure 8 shows that there is a strong correlation
between $G_{12}$ and $G$ for objects in our dataset with a phase curve fit in
both models. We have verified that this is caused by the different
formulations of the two models and is not due to systematic effects in our
phase curve fitting procedure. As detailed in appendix B (figure 17), we
generated a range of ideal model $HG$ phase curves with $0\leq G\ \leq 1$ and
then fitted those simulated observations with the $H_{12}G_{12}$ model. The
locus of best fit $G_{12}$ values for phase curves with parameter $G$ is shown
in the top panel of figure 8. It can be seen that fitting P16 to observed
phase curves naturally gives rise to a range of $G_{12}$ larger than $(0,1)$ .

Figure 7:  Log-scale distribution of P16 $G_{12}$ phase parameters for the
quality controlled SAP dataset. Distributions of objects with S and C-type
taxonomic classifications from Berthier et al. (2023) are also shown. Similar
to figure 4 the median values of the S and C-type distributions are indicated
by coloured vertical lines (with values of $G_{12}(\mathrm{S})=0.22$,
$G_{12}(\mathrm{C})=0.51$) and the limits $0<G_{12}<1$ for a “physical” phase
curve are shown as vertical red dashed lines.

Additionally we find an offset between the absolute magnitudes in the two
systems as shown in figure 9. This is similar to the results of M21 (their
figure 4). There is a trend of increasing offset in $H$ for objects with
fainter absolute magnitudes. This means that applying the P16 model results in
lower values of $H$ (i.e. brighter, larger objects) than for the corresponding
B89 fit. We confirm that this effect also arises due to the inherent
differences between the two models. P16 is more sensitive at modelling
opposition effects and has a third basis function which accounts for non-
linear brightness surges at low phase angles (equation 5), which can lead to
brighter values of $H$. In the test described above the exact $HG$ phase curve
extends all the way to zero phase. The $H_{12}G_{12}$ model cannot match $HG$
for both the approximately linear behaviour at large phase angles and the non-
linear opposition brightening at low phase angles (figure 17). The $HG$ and
$H_{12}G_{12}$ models show the most significant divergence at low phase angles
which explains the differences in absolute magnitude between the two models.

We show the relation between $H_{12}-H$ and phase parameter $G$ for the phase
curves in our database in the lower panel of figure 8. We note that it is also
well described by a locus derived from fitting B89 phase curves with the P16
system (see again figure 17 in appendix B). Therefore the shift in absolute
magnitude is an inherent property of comparing two different phase curve
models rather than being caused by the phase angle coverage or photometric
uncertainty of ATLAS data. Clearly the P16 model more accurately captures the
opposition surge which is present for many asteroids, however, for historical
reasons the B89 model is often still in use (e.g. MPC, JPL, astorb). This $H$
difference is systematic across both $o$ and $c$-filter phase curve fits,
hence the resulting $H_{c}-H_{o}$ colours for both models are practically
indistinguishable. In the P16 system we find median $H_{c}-H_{o}$ colours that
are nearly identical to the B89 system: $0.39\pm 0.05$ and $0.27\pm 0.06$ for
S and C-types respectively (uncertainties of one standard deviation).

---
Figure 8:  Two dimensional histogram showing the relation between the $G_{12}$
and $G$ parameters for objects in our database (upper panel). A 3rd-order
polynomial is fitted to the data as shown by the black curve (equation 9a).
The red markers indicate the values of $G_{12}$ when a $H_{12}G_{12}$ model is
fitted to a $HG$ model (appendix B). This empirical test follows the trend in
the ATLAS data closely. The bottom panel shows the two dimensional histogram
of the difference in fitted absolute magnitude $H_{12}-H$ vs. $G$. Again, the
black curve is a 3rd-order polynomial fit (equation 9b) and the red markers
indicate the expected difference in absolute magnitude when $H_{12}G_{12}$ is
fitted to $HG$. The horizontal black line indicates the zero value of the
$H_{12}-H$ axis.  Figure 9:  Two dimensional histogram which shows the
variation of $\Delta H=H_{12}-H$ with $H$. The median difference in absolute
magnitude between the two phase curve models is $-0.11$ (indicated by the
horizontal red line). There is a negative correlation between $\Delta H$ and
$H$, as well as a trend for a larger dispersion in $\Delta H$ for fainter
objects (see also figure 8). The horizontal black line indicates the zero
value of the $H_{12}-H$ axis.

###  3.2 Comparison with previous studies

---
Figure 10:  Histograms showing the difference in $H$ (upper panel) and
$G_{12}$ (lower panel) of our phase curve fits in the P16 system, cross-
matched with the results of Mahlke et al. (2021) who also used ATLAS data. The
median of each distribution is indicated by the red vertical line with values
of -0.01 and -0.06 magnitudes for the upper and lower panels respectively. The
vertical black lines mark the zero value of the $x$-axis.

The most direct comparison of our results to published work is with the study
by M21, who also made use of ATLAS data for asteroid phase curve analysis.
They utilised the photometry for $\sim 180,000$ asteroids in the period of
2015-2019. From this they determined phase curves with the Muinonen et al.
(2010) $HG_{1}G_{2}$ and P16 $H_{12}G_{12}$ models for $\sim 94000$ asteroids
using MCMC Bayesian parameter inference, with fits in both $o$ and $c$
filters. Similar to our methods they did not directly account for rotational
effects in their phase curve model. However they did estimate the influence of
rotational and apparition driven magnitude variations on the uncertainties of
their phase curve fits. In this work we considered the phase curves of
separate apparitions and performed more rigorous rejection of outlying
photometry. Another important difference in methodology is that we did not
restrict phase curve parameters when fitting, whereas M21 restricted MCMC
walkers to bounds of $0<G_{12}<1$ to ensure monotonically increasing phase
curves. This leads to a pile up of $G_{12}$ values at the boundaries as seen
in figure 5 of M21. Despite these differences there is reasonable agreement
for absolute magnitudes of objects that are present in both datasets (for
49,558 matched objects); most objects also have similar values of $G_{12}$
(figure 10). The median difference in absolute magnitude
$H_{12,\text{Mahlke}}-H_{12}=-0.01\pm 0.08\,$\mathrm{m}\mathrm{a}\mathrm{g}$$
(one standard deviation) and the difference in $G_{12}$ is also relatively
near zero with a median $G_{12,\text{Mahlke}}-G_{12}=-0.06$ albeit with a
higher standard deviation of $0.31$.

In another study McNeill et al. (2021), hereafter McN21, published phase curve
properties and colours of 342 Jupiter Trojans asteroids using ATLAS
observations from 2015-2018. They fitted $HG$ B89 phase curves to the $o$ and
$c$-filter data, accounting for rotational brightness variation by repeatedly
fitting to a randomised 50% subsample of observations. In this work we
considered using a similar approach, however we found that it did not
significantly change the estimates of $HG$ (or their uncertainties). This is
because a sub-sample of a phase curve exhibiting strong rotational effects
will have similar residuals to the original data. The results of McN21 provide
the $G$ parameter and $c-o$ colour of their sample, which we were able to
match with 317 objects in our science dataset (253 with colours). Despite
using more data from the same survey, not all objects in the McN21 sample pass
our quality control checks. We applied a robust methodology to reject outlying
photometry and assessed the quality of phase curve fits to individual
apparitions. Furthermore we rejected any objects displaying strong apparition
effects from our colour analysis.

When we compare the matched objects we find a median phase parameter
difference of $G_{\text{McNeill}}-G=0.06$, with quite a large standard
deviation of 0.28 (figure 11). Similarly for the difference in $H_{c}-H_{o}$
colours, the median difference is $-0.008\pm 0.1$ (figure 12). These
difference distributions are centred relatively close to zero, however, the
spread in values between the two studies should not be surprising. These
differences primarily arise due to the differing methodology described above,
furthermore, we have made use of additional ATLAS observations that were taken
since their analysis was performed. In particular, the dataset of McN21
contains a small number of objects with $1<G<2$ for which we obtain more
realistic values of $G$.

Considering another study focused on the Jupiter Trojans, SB21 determined $HG$
B89 phase curves for a significant number of Jupiter Trojans asteroids ($1049$
objects) using ZTF data (Bellm et al., 2019). Importantly this study accounts
for rotational lightcurve variation by including the probability of an
observation being made at a particular phase for a given amplitude and
building a likelihood model for the apparent magnitude of each observation.
They determine posterior distributions for $H$, $G$ and $g-r$ colour, which
more accurately capture the true uncertainty in these parameters rather than
the formal fitted uncertainties quoted in this work. Furthermore they account
for offsets in $H$ and changing rotational amplitude for different apparitions
of an asteroid. We match Jupiter Trojans that appear in both our study and
that of SB21, for a total of 571 cross-matched objects.

Figure 11:  Histogram distributions showing the difference in $G$ (B89) for
phase curve fits in our SAP database cross-matched with the previous results
of Waszczak et al. (2015), McNeill et al. (2021) and Schemel & Brown (2021).
Each distribution is shown as a probability density to better compare datasets
with significantly different numbers of objects. The coloured $x$-axis ticks
indicate the median value of the corresponding distribution, whereas the
vertical black line indicates the zero value of the $x$-axis. The median
offset for each distribution is $G_{\mathrm{Schemel}}-G=-0.13$,
$G_{\mathrm{McNeill}}-G=0.06$, $G_{\mathrm{Waszczak}}-G=-0.07$.

There is reasonable agreement between $H$ given that we have not accounted for
rotational effects, the median difference in $H_{\text{Schemel}}-H$ is
$0.07\pm 0.13$. However, there is a significant offset in $G$ distributions:
$-0.13\pm 0.18$ (figure 11). This means that in general our work finds larger
values of $G$ for Jupiter Trojans, which corresponds to flatter phase curves
(at phase angles $>5\,$\mathrm{d}\mathrm{e}\mathrm{g}$$ in the approximately
linear regime). This arises due to the increased photometric uncertainty of
ATLAS observations for the dark and distant Jupiter Trojans, the objects with
the greatest semimajor axes in our study. This uncertainty, combined with
rotational brightness variation means that the ATLAS phase curves have a large
scatter around the expected phase curve behaviour, giving a flatter
appearance. Our simple methodology tends to find solutions with higher $G$
when such scatter is significant. This seems to be a common effect when
rotational effects are not accounted for. A recent study by Donaldson et al.
(2023) showed that the phase curve of the Jupiter family comet nucleus
162P/Siding-Spring (which has low albedo, similar to the Jupiter Trojans) got
steeper when the rotational brightness variation was subtracted from the phase
curve. Furthermore, it would appear that the scatter introduced by rotational
variation (and also photometric uncertainty) tends to “disguise” the
opposition brightness surge of an asteroid at low phase angles. However if
rotational effects were accounted for one would expect the opposition surge
(if present) to be more accurately resolved. In this case one should fit the
phase curve with a model such as P16 that better resolves this non-linear
surge when compared to B89 (figure 17).

We also consider the work of Waszczak et al. (2015) who simultaneously
determined rotation periods and phase curves of asteroids observed by the
Palomar Transient Factory (Law et al., 2009) in $R$ and $g$-filters, resulting
in a dataset of $\sim 9000$ objects. They determined phase curves in the two
parameter Shevchenko (1997) model, $HG$ B89 and the two parameter model of
Muinonen et al. (2010) (the precursor to P16). Therefore we compare their B89
phase curve parameters to ours for 5199 MBAs and Jupiter Trojans common to
both datasets (figure 11). The median difference in $G$ and its standard
deviation, $-0.07\pm 0.20$, are of similar magnitude to the difference with
McN21. However both Waszczak et al. (2015) and SB21 $\Delta G$ distributions
are shifted towards negative values. These works accounted for rotational
variation in their phase curve models, and for most objects obtain lower
values of $G$ compared to this work and that of McN21. This again implies that
magnitude dispersion due to rotation and photometric uncertainty tends to
flatten the fitted phase curve and push the distribution of $G$ to higher
values.

Figure 12:  Histogram comparing the ATLAS $H_{c}-H_{o}$ colour distribution of
objects in this work cross-matched with the previous colours ($c-o$) published
by McNeill et al. (2021), who also investigated the properties of the Jupiter
Trojans using ATLAS data. The median value of the distribution is -0.008
magnitudes, the vertical black line indicates the zero value of the $x$-axis.

##  4 DISCUSSION

###  4.1 The Selected ATLAS Phase-curves database (SAP)

The SAP database contains absolute magnitudes and phase curve parameters in
the B89 and/or P16 systems for 106055 main belt and Trojan asteroids. The
median measured value of $H_{o}$ changes from
$16.6\,$\mathrm{m}\mathrm{a}\mathrm{g}$$ at the inner edge of the main belt at
$a=1.9$ au to $13.0\,$\mathrm{m}\mathrm{a}\mathrm{g}$$ for the Hilda asteroids
at $a=4$ au, and to $11.6\,$\mathrm{m}\mathrm{a}\mathrm{g}$$ for the Jupiter
Trojan asteroids at $a=5.2$ au. The median uncertainties of the non-linear
least squares fits for all phase curves in the SAP database are
$\sigma_{G}=0.04$ and $\sigma_{G12}=0.11$. In section 4.2 we examine the phase
parameter and colour properties of asteroid families in the SAP database. Then
in section 4.3 we describe an initial analysis of Jupiter Trojan asteroids
using these data.

###  4.2 Asteroid families in the SAP database

To illustrate the usefulness of the SAP database, in figure 13 we show the B89
phase parameters as a function of proper orbital elements (obtained from the
astdys database666https://newton.spacedys.com/astdys/ accessed 3rd January
2023), highlighting the asteroid family members using the astdys
classifications. In this figure we have indicated the presence of non-family
MBAs using a greyscale two dimensional histogram and the dynamical family
members are overplotted as coloured markers, appearing as clusters in proper
element space. The marker colour in figure 13 provides the phase parameter $G$
of each asteroid and many of the dynamical families clearly share similar
phase parameters. This indicates a shared albedo/composition amongst family
members, which is to be expected when a family is created from the
catastrophic collisional fragmentation of a parent asteroid. There is a
general decrease in the $G$ of asteroid families as proper semimajor axis
increases, which is driven primarily by the increasing fraction of C-type
asteroids in the outer main belt.

We examine in more detail the Nysa-Polana asteroid family, which has
previously been identified as having a complex dynamical structure. It is
considered to be a large “clan” grouping also containing the low-numbered
asteroid Hertha, for which the group may also be named (Cellino et al., 2001;
Milani et al., 2014). In addition to the sub-structure of the group in proper
element space it has been shown that this large family also contains a mixture
of asteroid taxonomies. Erasmus et al. (2020), also using ATLAS data, derived
probabilities of an S-like vs. C-like composition for Nysa-Polana members
using their own measurements of the $H_{c}-H_{o}$ colour. They showed that
there is a clear two-part taxonomic mixture within this dynamical grouping.
Here we build upon this study by considering the combination of the
$H_{c}-H_{o}$ colour and phase parameter $G$ derived in this work for
identifying taxonomy. Figure 14 (left) shows the proper orbital elements of
the Nysa-Polana (Hertha) family members that are present in the SAP database.
The presence of two clusters in proper element space is clear and to indicate
the taxonomic diversity of these objects we have used a two dimensional
mapping of $H_{c}-H_{o}$ vs. $G$ to colour each marker. This colour map is
displayed in the right-hand panel of figure 14 and it illustrates that this
combination of $H_{c}-H_{o}$ and $G$ can provide a powerful tool for
identifying broad taxonomic information of asteroids observed by sparse, wide
field surveys. We have indicated the distribution of Nysa-Polana members in
this $H_{c}-H_{o}$ vs. $G$ space and indicated the location of the median S
and C-type parameters as measured from the full SAP database (figure 4). This
analysis shows that there is strong correlation between dynamical clustering
and taxonomic clustering in the Nysa-Polana group. These clusters have
distinct S-like and C-like features and supports the theory that this group is
made up of nested families with distinct parent bodies. Further analysis by
the community of the dynamical and collisional families present in the SAP
dataset is encouraged.

---
Figure 13:  Proper orbital element distributions of MBAs in the SAP database
(elements from astdys). The proper semimajor axis is plotted against the sine
of proper inclination (upper panel) and proper eccentricity (lower panel). The
63,197 non-family members in the dataset are displayed as a two dimensional
histogram in greyscale. The 25,479 asteroids that are a member of a family (as
classified by astdys) are colour-coded by $G$ with a diverging colour map,
where the colour transition is centred on the median value of $G$. The same
colour scheme is used for both plots and objects with phase parameters outside
the $1-99\%$ quantile of the family member $G$ distribution were excluded from
the plot in order to avoid stretching of the colour map by the small number of
more extreme values. As such we only display asteroid family members with
$0.03<G<0.60$ (where the median $G$ of this sample is 0.26). Some of the
larger asteroid families are labelled by black crosses at their mean proper
elements (Milani et al., 2014).  Figure 14: Here we show the astdys proper
orbital element distribution of Nysa-Polana (N-P) asteroid family members in
the SAP dataset, in proper eccentricity vs. proper semimajor axis space (left
panel). We indicate the proper elements of the main asteroids of the complex
(Nysa, Polana and Hertha) with shaped markers. Every other family member is
colour coded according to a two dimensional colour map, $H_{c}-H_{o}$ colour
against phase parameter $G$. This colour map is displayed in the right-hand
panel, where we have also marked the locations of the SAP Nysa-Polana family
members in this parameter space (black circles). The black contour indicates
the 10% level relative to the maximum count of the underlying two dimensional
histogram of the distribution. We have also indicated the median $H_{c}-H_{o}$
and $G$ of the S and C-type asteroids in the full SAP database (red cross and
plus respectively, values from figure 4).

###  4.3 Jupiter Trojan asteroids in the SAP database

Another important dynamical population to investigate with the SAP dataset are
the low-albedo, optically red Trojan asteroids in 1:1 mean-motion resonance
with Jupiter at orbital semimajor axes around $5.2\,$\mathrm{A}\mathrm{U}$$
(Marzari et al., 2002). The Jupiter Trojans are distributed in two groups that
are leading/trailing Jupiter by approximately
$60\,$\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}$$,
defined by the L4/L5 Lagrange points of the restricted three body problem. In
the Nice model of planetary migration (Gomes et al., 2005; Tsiganis et al.,
2005; Morbidelli et al., 2005) the Trojans formed exterior to the giant
planets and were scattered and captured onto their current orbits as the giant
planets dynamically evolved. Alongside $N$-body simulations supporting this
theory, observations of Jupiter Trojan surface composition has shown these
objects to be more similar to outer Solar System objects than the inner small
body population (Wong & Brown, 2016). There are over $10,000$ known Jupiter
Trojans (as listed by the MPC, complete down to around
$H=14\,$\mathrm{m}\mathrm{a}\mathrm{g}$$, Hendler & Malhotra, 2020), and
determination of the underlying size distribution indicates a number asymmetry
of $\sim 1.4$ between the L4/L5 groups (Grav et al., 2011). Within the Nice
model framework this could be caused by depletion of the trailing group by the
ejection of a massive planet from the early inner Solar System (Nesvorný et
al., 2013). Therefore there is interest in seeing how similar (or not) the
Jupiter Trojan groups are in terms of their physical properties as well. This
will help to determine if they indeed share the same parent population and
whether the number asymmetry has affected their properties through collisional
processing since depletion. Remote observations of Jupiter Trojans are
particularly relevant given the recent launch of the LUCY mission (Levison et
al., 2021). This mission will obtain in-situ observations of Jupiter Trojans
in both dynamical groups, which will provide a valuable anchor-point to the
studies made using disk-unresolved photometry.

The Jupiter Trojans have been well observed by the ATLAS survey, providing a
substantial and consistent photometric dataset with which to study these
objects. Notably, the Trojans have the lowest maximum phase angle
($\text{max}(\alpha)\leq 11^{\circ}$) in our database as they are the most
distant objects we have included in this study (figure 1). This naturally
limits the phase angle coverage, although we can still sample any opposition
surge for many of them, with observations of most Jupiter Trojans typically
reaching a phase angle of $1\,$\mathrm{d}\mathrm{e}\mathrm{g}$$. We note that
in 2017 the leading L4 cloud was at low galactic latitude and has since moved
out of the galactic plane, whereas the trailing L5 cloud has since been
approaching lower galactic latitudes. Hence both Jupiter Trojan clouds have
been sampled in ATLAS observations whilst away from the galactic plane, albeit
over different time periods, thanks to the long baseline of the survey.

Here we consider the phase curve properties of the Jupiter Trojans present in
our SAP dataset. This contains 740 Jupiter Trojans with high quality $HG$ fits
in the $o$-filter, of which 570 also have good fits in the $c$-filter suitable
for calculating a $H_{c}-H_{o}$ colour (section 2.6). 103/740 of these Jupiter
Trojans display strong apparition effects and our sample contains 1.98 more
Jupiter Trojans at L4 than L5. This is similar to the L4/L5 number ratio for
observed Jupiter Trojans (for counts of Jupiter Trojans recorded by the MPC,
Li et al., 2022) but is larger than the number ratio of $\sim 1.4$ when the
underlying size-frequency distributions are accounted for (Grav et al., 2011;
Uehata et al., 2022). We compare our distributions in absolute magnitude, $H$,
of these two groups (figure 15, left panels) using a two sample Kolmogorov-
Smirnov (KS) test. In a two sample KS test the parameter $D$ describes the
maximum distance (at a given value) between the cumulative distributions of
the two datasets being tested. The KS statistic is calculated with a $p$-value
which indicates the probability that the null hypothesis is true, i.e. that
both datasets were drawn from the same underlying distribution. Conventionally
a result is statistically significant when $p<0.05$ (see also equation 10).
For the L4/L5 $H$ distributions this test results in a KS statistic of $D=0.2$
(with a $p$-value of $3\times 10^{-6}$) from which we reject the null
hypothesis that the L4 and L5 absolute magnitudes are drawn from the same
distribution. There is no strong evidence in the literature that the L4, L5
groups should have different size distributions (e.g. Uehata et al., 2022). In
this case we determined that the differences in $H$ were introduced into our
study through the initial watchlist of asteroids monitored by rockAtlas. There
is no significant difference in the absolute magnitude distribution for the
L4/L5 Trojans when using $H_{V}$ from the full astorb dataset. However when we
compare the astorb $H_{V}$ of only the L4/L5 Trojans which were selected for
the rockAtlas watchlist, a KS test displays statistically significant
differences similar to those seen in figure 15 (left). Therefore we can rule
out observational bias between the L4 and L5 groups in the ATLAS survey as the
main cause of this effect, and attribute it instead to selection bias.

We performed similar tests on the remaining parameters retrieved from the SAP
dataset. For the B89 phase parameter $G$ the two sample KS test results are
$D=0.08$ and $p=0.18$. Likewise for $H_{c}-H_{o}$ colours (for the smaller
number of objects with good fits in both $o$ and $c$), the KS test results are
$D=0.08$ and $p=0.39$ (figure 15, middle and right panels respectively). These
$D$ values are near the boundary at which we would consider rejecting the null
hypothesis. To ensure that these results are not influenced by small number
statistics or dominated by specific objects in the dataset we performed
repeated KS tests on bootstrapped samples of each parameter. We assessed the
variation in $D$ for the resampled data, the results of which are detailed
further in appendix C (figure 18). Overall these tests confirm that there is
no statistical evidence for any differences in the distribution of phase
parameter $G$ or $H_{c}-H_{o}$ colour between the L4 and L5 Jupiter Trojans.
This is in contrast to the reported dissimilarity in $c-o$ colours of the
L4/L5 Jupiter Trojans by McN21. This difference in results is most likely due
to the shorter ATLAS survey baseline at the time of their study, which meant
that McN21 had fewer observations available for phase curve analysis.
Furthermore, although the statistically different $H$ distributions are due to
the rockAtlas watchlist, in this work we have not thoroughly accounted for
observational biases between L4 and L5 in the ATLAS dataset. Biases could have
been introduced by observations of the L4/L5 groups as they are moving away
from/towards the galactic plane over the course of the survey. We can see in
figure 15 that the ATLAS observations for the L4 group generally go down to
fainter absolute magnitudes. Therefore there could be size effects on $G$ and
$H_{c}-H_{o}$ in addition to the greater phase curve fitting uncertainty for
fainter objects (figure 5).

Figure 15:  Histogram distributions of $H$, $G$ and $H_{c}-H_{o}$ colour, in
the B89 system, for the L4 and L5 Jupiter Trojans in the SAP dataset (top row,
left, middle and right panels respectively). The distributions for the L4 and
L5 groups are shown in blue and orange respectively; the median of each
distribution is indicated by the vertical coloured line. The cumulative
distribution functions (CDFs) of each parameter are shown in the bottom row,
allowing us to compare the distributions via the two-sample Kolmogorov-Smirnov
test. The KS statistic $D$ is determined from the maximum difference between
the CDFs at a given value and when this is sufficiently large (at a level
determined by the $p$-value) we can reject the null hypothesis that the two
datasets are drawn from the same underlying distribution.

As a sanity check we repeated the above analysis but instead using the Jupiter
Trojan dataset published by SB21. At first we found that the KS tests for the
L4/L5 $H$ and $G$ distributions seemed to indicate statistically significant
differences (high $D$ values with low $p$). However these results are highly
dependent on the quality of the phase curve fit. SB21 report uncertainties
derived from the posterior distributions of their parameters, which are
generally larger than those reported in works like our own, but which are a
better representation of the true range of parameters compatible with the
observations. When we consider only their phase curves with better than median
uncertainties (using the same methods as SB21, in their work) then any
differences between the L4 and L5 groups become statistically insignificant,
similar to our results (figure 19).

Although the two dynamical clouds have similar distributions of phase curve
parameter, it has been shown by several authors that the Jupiter Trojans in
both clouds can be split into two compositional groups based on their colours
and/or spectra; the less-red and red sub-classes (Emery et al., 2011; Szabo et
al., 2007). Recently, it has been found by SB21 that these sub-classes show
distinct phase curves at optical wavelengths. As indicated in figure 16, this
bimodality in phase parameter $G$ is also present in our dataset. In this
figure we have plotted the phase parameter $G$ obtained from this work against
the $g-r$ colour derived from ZTF data by SB21, using only objects with better
than median uncertainties in both datasets (replicating their analysis, see
figure 6 of SB21, ). Of the 571 asteroids common to both datasets, 159 meet
this requirement. Similar to the previous work we find a strong positive
correlation between our $G$ values and their $g-r$ colours, with a Spearman
rank correlation coefficient of $\rho=0.51$ ($p=7\times 10^{-12}$). We then
divided objects into red and less-red groups by a ZTF colour boundary of
$g-r=0.565$ and analysed the corresponding distributions of $G$. We replicate
the results of SB21 whereby less-red Jupiter Trojans generally have lower $G$
than red Jupiter Trojans, albeit our phase parameters are uniformly shifted to
higher values as discussed previously in section 3.2. A KS test of the two
distributions prove that they are indeed distinct ($D=0.48$, $p=2\times
10^{-8}$).

We also search for correlation between $G$ and $H_{c}-H_{o}$ colour for the
570 Jupiter Trojans with colours in the ATLAS dataset. The Spearman rank
correlation coefficient of these parameters is $\rho=0.07$ ($p=0.08$), which
indicates that there is no significant correlation. This is due to the large
width of the ATLAS $o$ and $c$-filters which overlap around the spectral
“kink” which distinguishes the red and less-red types (Emery et al., 2011;
McNeill et al., 2021).

In summary, our database of ATLAS phase curves shows no significant
differences in phase parameter $G$ or $H_{c}-H_{o}$ colour between the L4 and
L5 Jupiter Trojans. This is in general agreement with previous studies and
supports the consensus that both groups of Jupiter Trojans were captured from
a common source population. Although there is evidence for a difference in the
average shapes of the L4/L5 groups (presumably due to differences in number
densities and collisional environment, McN21, ), ATLAS data shows no
corresponding difference in their surface properties caused by different
collision histories. Furthermore, although ATLAS $H_{c}-H_{o}$ colour can be a
useful diagnostic when distinguishing between the S and C-type asteroids
(figure 6 and also Erasmus et al., 2020) these broad filters are less
effective at differentiating more subtle spectral features, such as the red
and less-red Jupiter Trojans. Regardless, the phase parameter $G$ captures
information on the differences in surface properties of the two Jupiter Trojan
classes. When combined with more detailed colour information from other
surveys (figure 16) a stronger distinction between the red and less-red
Jupiter Trojans can be inferred, highlighting the benefits of utilising
information from multiple sources.

---
Figure 16:  Scatter plot of phase parameter $G$ from the SAP dataset against
the ZTF $g-r$ colours from Schemel & Brown (2021), for objects matched between
the two datasets (upper panel). Similar to their analysis we only include
objects from both datasets with lower than median uncertainties (compare to
their figure 6). The red vertical line corresponds to the colour boundary at
$g-r=0.565$ which divides the “red” and “less-red” Jupiter Trojan
compositional classes The lower panel shows the histogram distributions of
ATLAS $G$ for the red and less-red Jupiter Trojans from the panel above.

###  4.4 Limitations of the database

####  4.4.1 Outlying photometry

In order to reliably and automatically fit phase curves to a large number of
objects we developed various selection criteria to reject outlying
observations before and during the phase fitting algorithm (section 2.4). It
is necessary to remove any false matches caused by the automated nature of the
rockAtlas pipeline, however, we risk throwing away valid measurements, e.g.
increases in asteroid brightness due to activity/outburst or
rotation/apparition effects. Previous work such as M21 addressed this by
conducting a straight $1\,$\mathrm{m}\mathrm{a}\mathrm{g}$$ difference cut
between the predicted and observed magnitudes. We found that a combination of
data cuts based on magnitude difference and standard deviation of residuals
provided a more robust method of outlier rejection, which we confirmed through
a process of visual inspection of the extreme outliers in the resulting
dataset.

Furthermore, we considered additional properties of each observation. We
rejected all data points within $10^{\circ}$ of the galactic plane as these
measurements have a high probability of stellar contamination. In some cases
this led to entire apparitions being rejected as an object moved through the
galactic plane. On one hand such a rough cut could throw away valid
measurements, on the other when the Milky Way thickness exceeds our
$20^{\circ}$ span (e.g. near the galactic bulge) we may not be rejecting
contaminated data. However, we found that most galactic plane sources were
rejected anyway by our strict $\theta<1\arcsec$ constraint on rockAtlas
ephemeris matching.

Our methods for matching sources could be improved further by crossmatching
all detections with stationary source catalogues in order to flag data points
contaminated by nearby sources and to conduct a more precise rejection. This
is currently beyond the scope of this work, however, we have shown that our
results are generally consistent with existing literature which also used
ATLAS data (section 3.2). Therefore the rockAtlas observation matching is not
a significant source of uncertainty.

####  4.4.2 Catalogue contamination by poor quality fits

One of the complex aspects of this work was the production of a final high
quality phase curve dataset. We fitted all available apparitions and then
ideally one should be able to use the uncertainties of the fitted phase curve
parameters to assess the quality of the dataset. However as discussed
previously, we do not account for rotational brightness variation and
therefore the formal uncertainties produced by the Levenberg-Marquadt
algorithm do not reflect the true range of possible values for the phase curve
parameters. We therefore had to develop a series of cuts based on a variety of
metrics which we tried to choose in a relatively agnostic way so as to not
overly influence the final results (section 2.5 and appendix A). We did this
by considering especially the phase curves of objects with unusually high
values of $G$ and $\sigma_{G}$ and adjusting our metrics such that these truly
unphysical fits were rejected. We acknowledge that this approach relies on
somewhat subjective decisions. A more realistic estimate of the phase
parameter uncertainties could have been obtained by methods in which
rotational variation is accounted for (e.g. SB21, ). However, this would add
further computational cost to calculating phase parameters for each apparition
of a large number of objects. We did consider “quick” solutions to deal with
rotational effects, e.g. binning the observed magnitude over some timescale or
by fitting random samples of the observations, however we found that these
methods did not significantly improve the quality of the final dataset.

####  4.4.3 Phase curve wavelength dependence

We determined the phase parameter ($G$, $G_{12}$) for each object from its
primary apparition in the $o$-filter, and kept this phase parameter fixed when
fitting all other apparitions. There is evidence that phase curve properties
vary with wavelength in a “phase-coloring” effect (Alvarez-Candal et al.,
2022a; Ďurech et al., 2020). However in the ATLAS dataset the coverage in the
$c$-filter is much sparser compared to the $o$-filter. In most cases we lacked
sufficient data to obtain an independent phase curve for the $c$-filter data.
In order to obtain usable colour information we were forced to assume the
slope parameters to be wavelength independent in order to find $H_{o}$ and
$H_{c}$. This is a reasonable assumption given the accuracy of ATLAS
photometry and not accounting for rotational variation. These effects probably
dominate over any changes in phase parameter due to the wavelength of
observation.

####  4.4.4 Rotational lightcurve effects

While developing our methodology we attempted several simple methods to
incorporate rotational brightness variations into our phase curve fits. For
example we considered a Monte Carlo approach whereby each observation was re-
sampled based on its measurement uncertainty and phase curves were repeatedly
fit (similar to McN21, ). We also attempted to average out rotational
variation by binning reduced magnitudes nightly when fitting the phase curves.
However, we found that these methods did not adequately account for rotational
variations in the lightcurves. When resampling observations for objects
displaying clear lightcurve variations the photometric uncertainty is
generally small compared to the amplitude and so the fitted phase curve
parameters and their uncertainties did not differ significantly from our
original method. Similarly, binning the data did not have much effect on the
final parameters. Furthermore, binning observations of objects with a wide
possible range of rotational periods (minutes to days) compared to the ATLAS
observational cadence would likely introduce a number of additional biases to
our analysis. A more robust treatment would account for the unknown rotational
phase and amplitude of each observation and so the true uncertainties of the
phase curve parameters would be much wider but reflect better the true
distribution of possible parameters (e.g. Waszczak et al. (2015); Alvarez-
Candal et al. (2022b), SB21). However such methods are generally more
computationally intensive to implement, therefore in this work we stuck with a
simple phase curve fitting method, neglecting rotational effects as described
in section 2.4. This allowed us to consider a sizeable sample of asteroids,
and their individual apparitions, in a timely manner. As such, the formal fit
uncertainties produced in this work should be used to compare relative
goodness of fit within the dataset; they do not necessarily reflect the true
uncertainty in the fitted phase curve parameters if rotational effects are
strong.

###  4.5 Applicability to future surveys

In this work we have demonstrated that sparse wide field survey data can be
exploited to investigate and broadly compare the properties of different
asteroid populations. In particular, by considering the phase parameters and
surface colours derived in this work one may compare and contrast different
dynamical groups such as the asteroid families and Jupiter Trojans. However as
discussed above, in order to make more precise assessments one must take into
account the rotational variation introduced by the asteroid lightcurve, either
in a statistical sense (e.g. Waszczak et al. (2015); Alvarez-Candal et al.
(2022b), SB21) or with more detailed lightcurve shape modelling (e.g. Ďurech
et al., 2020). Despite our simple phase curve model, we have built upon
previous phase curve studies by considering the phase curve properties of
individual apparitions for a large sample of over $100,000$ asteroids in the
Main Asteroid Belt to Jupiter Trojan region, with data spanning a significant
temporal baseline since 2017. By doing so we have been able to identify
objects that are undergoing changes in observed absolute magnitude between
apparitions. By excluding these objects from our science analysis of absolute
magnitudes and colours we have improved the purity of our sample.

This variation in brightness known as the apparition effect is driven by the
shape and spin state of an object and means that closer study is required in
order to determine its true absolute magnitude. Studies such as Ďurech et al.
(2020) have shown that ATLAS sparse photometry can be analysed using
lightcurve inversion techniques, although in most cases additional dense
lightcurve observations are required to reduce the parameter space to be
searched. Of the 9168 objects we have identified with strong apparition
effects (in either B89 or P16) only 2828 have reported periods in SsODNet
(Berthier et al., 2023, accessed 15th March 2024). Therefore this list of
objects with demonstrated lightcurve effects makes for promising targets for
future dedicated observations in order to accurately define their physical
properties.

Assessing apparition effects informs us of the elongation and spin pole
orientation of an object. If we could disentangle the spin pole orientation
for a large database of objects showing apparition effects we could assess the
competing effects of collisional vs. YORP evolution on the spin pole (Pravec
et al., 2002). This could be achieved by searching for relations between
asteroid size and spin pole orientation, whereby the spin pole of asteroids is
expected to evolve away from the orbital plane at a rate which is inversely
proportional to size (e.g. Cibulková et al., 2016). Apparition effect objects
could also be relevant for comparing shapes and/or collisional evolution, in
particular for the L4/L5 Jupiter Trojans which show indications of different
shape distributions (McN21).

The presence of the aforementioned apparition effects must be accounted for
when analysing the results of wide field surveys expected to observe large
numbers of objects over a long baseline. In particular, LSST will discover
approximately $6\times 10^{6}$ MBAs over 10 years, with 100s of observations
for most objects. In addition to moving object discovery, LSST will provide
automatically derived properties of observed asteroids such as lightcurves,
phase curves and colours (LSST Science Collaboration et al., 2009). Given the
volume and depth of the survey it will be difficult to followup the expected
large numbers of faint asteroids, therefore LSST will likely be the only
source of data for many of these objects. Considering the apparition effect is
therefore essential in order to determine accurate absolute magnitudes, which
directly effects studies regarding the size distribution of small Solar System
bodies, such as the formation and evolution of planetesimals, properties of
asteroid families, planetary defence, etc.

Furthermore, our approach which has involved detailed analysis of phase curve
residuals offers a framework to monitor for unusual observations that might
indicate an object of interest. We considered a range of statistics for the
phase curve residuals, and in the vast majority of cases any rejected
apparitions were due to poor observational coverage or low quality photometry.
However, certain physical processes may also lead to deviations from the
expected phase curve behaviour, e.g. rapid changes in viewing geometry
(Jackson et al., 2022) or sudden increases in brightness due to the onset of
activity (Jewitt & Hsieh, 2022). For example, asteroid 2008 GO98 was found to
be active in July 2017 (Kokhirova et al., 2021). Our methodology produced a
most unusual phase curve for this object, with an extreme value of $G=1675$
which was determined from an apparition spanning from 8th June until 9th
December 2017. Future work in this area could involve identifying these
outliers in phase curve parameter space and using the ATLAS forced photometry
server to get cutouts to search for activity. Continuous updates of derived
physical properties (such as absolute magnitude, phase curve and colours) as
new observations are made will be an essential part of large surveys such as
LSST. Therefore learning how to identify truly unusual objects in this process
will greatly advance the study of relatively rare objects such as the active
asteroids.

##  5 Conclusions

We have compiled a database of ATLAS asteroid photometry by matching the
ephemerides of 4.5$\times 10^{5}$ asteroids with accurate orbital elements to
the detections produced by the ATLAS dophot pipeline, using rockAtlas (Young,
2023). We fitted phase curves in the $HG$, $H_{12}G_{12}$ models (Bowell et
al., 1989; Penttilä et al., 2016) to the apparition in the $o$-filter with the
best combination of number of observations and phase angle coverage. We then
fitted all other apparitions for $H$ with these derived phase parameters
fixed. This allowed us to assess any changes in $H$ during the survey and also
to determine the $H_{c}-H_{o}$ surface colour of many asteroids. We used
various metrics to select 106,055 asteroids with well defined phase curve fits
for the final Selected ATLAS Phase-curves (SAP) dataset, which is available
online.

Within the SAP dataset we identified around 9000 asteroids that show strong
apparition effects, i.e. large shifts in $H$ over the long baseline of ATLAS
observations. These objects most likely have elongated shapes and tilted spin
poles in order to display these effects. We suggest such asteroids should be
prioritised for followup lightcurve observations and analysis.

Our measured absolute magnitudes and phase curve parameters were compared to
previous studies, and we found that our analysis which does not account for
rotation is in general agreement with similar studies. However, our method
fails to accurately reproduce the Jupiter Trojan $G$ distribution compared to
more sophisticated methods that include rotational variation (Waszczak et al.,
2015; Schemel & Brown, 2021). We found that by not taking into account
rotational brightness variations when fitting the phase curves, the increased
dispersion and larger photometric uncertainties for the faint Jupiter Trojans
led to flatter phase curves in this work.

In spite of this, we found that the asteroid colours and phase parameters
provided in the SAP database can be a useful tool for comparing the properties
of various asteroid populations. We demonstrated this for main belt asteroid
families, whereby families can be distinguished by their differing
distributions of $G$. Our dataset illustrates the known trend towards lower
values of $G$ as semimajor axis increases, caused by the increasing fraction
of C-type asteroids in the outer main belt. Additionally, we examined in
closer detail the Nysa-Polana (Hertha) dynamical complex. The SAP database
allowed us to use both $H_{c}-H_{o}$ colour and $G$ to infer compositional
differences of the family members. We demonstrated that this combination of
physical parameters traces the differing compositions within the dynamical
sub-structure of the larger family group. Our results support theories that
the Nysa-Polana group is composed of nested families formed of at least two
parent bodies, which had different S-like and C-like compositions.

We also conducted a study of the Jupiter Trojans in the SAP database, the
dynamical and physical properties of which help to constrain the evolution of
the Solar System and the migration of the planets in particular. We found that
there were no statistically significant differences between the L4 and L5
groups for the distributions of $G$ and $H_{c}-H_{o}$. We did measure a
discrepancy in the distribution of $H$ between the L4 and L5 groups, however,
we determined that this resulted from an initial bias in asteroid selection.
Furthermore, we found that the broad $o$ and $c$-filters could not uniquely
distinguish between red and less-red Jupiter Trojan compositional classes,
however, when combined with colour measurements from other sources (e.g. ZTF
$g-r$ colours; Schemel & Brown, 2021) the difference in ATLAS phase parameter
$G$ between the two classes was emphasised. In general our results support the
consensus that the Jupiter Trojans were captured from a common source
population.

In this work we chose a phase curve model without rotational lightcurve
variations to better investigate individual apparitions of a large number of
asteroids. However, we recommend that future studies account for rotational
effects, e.g. via likelihood models including rotational variation (for sparse
data e.g. Schemel & Brown, 2021) and/or full CLI techniques (for objects with
additional lightcurve rotation periods e.g. Ďurech et al., 2020). Accurate
phase curve modelling with realistic uncertainties on absolute magnitude and
phase curve parameters are essential when comparing the optical properties of
different small body populations. This will be particularly relevant for
large, sparse surveys such as LSST which will discover and characterise
approximately 6$\times 10^{6}$ asteroids over its 10 year baseline (LSST
Science Collaboration et al., 2009). Many LSST objects will be too faint
and/or numerous for dedicated followup therefore extracting accurate
parameters from sparse data alone, and quickly highlighting interesting
outliers, is a pressing issue.

Author’s note: Since submission we have been made aware of development of the
sHG1G2 model (Carry et al., 2024); i.e. a “spinned” version of the
$HG_{1}G_{2}$ Muinonen et al. (2010) model. This enhanced phase curve model
also fits the oblateness and spin axis of an object in order to account for
shape driven rotational effects. Such a model will certainly help to address
the long-standing issues in determining accurate phase curve parameters from
sparse data.

## Acknowledgements

The authors wish to thank the reviewer, Josef Ďurech, for useful comments that
helped to improve this manuscript. This work has made use of data from the
Asteroid Terrestrial-impact Last Alert System (ATLAS) project. ATLAS is
primarily funded to search for near earth asteroids through NASA grants
NN12AR55G, 80NSSC18K0284, and 80NSSC18K1575; byproducts of the NEO search
include images and catalogs from the survey area. The ATLAS science products
have been made possible through the contributions of the University of Hawaii
Institute for Astronomy, the Queen’s University Belfast, the Space Telescope
Science Institute, the South African Astronomical Observatory (SAAO), and the
Millennium Institute of Astrophysics (MAS), Chile. This research has made use
of IMCCE’s SsODNet VO service (https://ssp.imcce.fr/webservices/ssodnet/).
J.E.R. acknowledges support from STFC award ST/T00021X/1 and the Royal
Society. A.F. acknowledges support from STFC award ST/T00021X/1. J.E.R. also
acknowledges the support and guidance of Cyrielle Opitom, whose discussions
helped improve the quality of the work. For the purpose of open access, the
author has applied a Creative Commons Attribution (CC BY) licence to any
Author Accepted Manuscript version arising from this submission. In addition
to software cited in the paper, the following packages were also used in this
work: `matplotlib` (Hunter, 2007), `numpy` (Harris et al., 2020), `scipy`
(Virtanen et al., 2020), `scikit-learn` (Pedregosa et al., 2011), `pandas`
(McKinney, 2010), `pds4_tools` (Nagdimunov, 2021), `pycolormap_2d` (Spinner,
2022).

## Data Availability

The databases constructed in this work will be available in their entirety via
the Queen’s University Belfast PURE open data repository
(https://pure.qub.ac.uk/en/datasets/), and will also be submitted to the NASA
PDS Small Bodies Node (https://pds-smallbodies.astro.umd.edu). They may also
be obtained via reasonable e-mail request to the lead authors.

## References

  * Alonso-García et al. (2012) Alonso-García J., Mateo M., Sen B., Banerjee M., Catelan M., Minniti D., von Braun K., 2012, The Astronomical Journal, 143, 70
  * Alvarez-Candal et al. (2022a) Alvarez-Candal A., Corral S. J., Colazo M., 2022a, Absolute Colors and Phase Coefficients of Asteroids (arxiv:2209.13246)
  * Alvarez-Candal et al. (2022b) Alvarez-Candal A., Benavidez P. G., Campo Bagatin A., Santana-Ros T., 2022b, Astronomy & Astrophysics, 657, A80
  * Astropy Collaboration et al. (2022) Astropy Collaboration et al., 2022, The Astrophysical Journal, 935, 167
  * Bellm et al. (2019) Bellm E. C., et al., 2019, Publications of the Astronomical Society of the Pacific, 131, 018002
  * Belskaya & Shevchenko (2000) Belskaya I. N., Shevchenko V. G., 2000, Icarus, 147, 94
  * Berthier et al. (2023) Berthier J., Carry B., Mahlke M., Normand J., 2023, Astronomy & Astrophysics, 671, A151
  * Bowell et al. (1989) Bowell E., Hapke B., Domingue D., Lumme K., Peltoniemi J., Harris A. W., 1989, Application of Photometric Models to Asteroids.
  * Carry et al. (2024) Carry B., Peloton J., Montagner R. L., Mahlke M., Berthier J., 2024, Combined Spin Orientation and Phase Function of Asteroids (arxiv:2403.20179)
  * Cellino et al. (2001) Cellino A., Zappalà V., Doressoundiram A., Di Martino M., Bendjoya Ph., Dotto E., Migliorini F., 2001, Icarus, 152, 225
  * Christensen et al. (2012) Christensen E., et al., 2012, in AAS/Division for Planetary Sciences Meeting Abstracts #44.
  * Chyba Rabeendran & Denneau (2021) Chyba Rabeendran A., Denneau L., 2021, Publications of the Astronomical Society of the Pacific, 133, 034501
  * Cibulková et al. (2016) Cibulková H., Ďurech J., Vokrouhlický D., Kaasalainen M., Oszkiewicz D. A., 2016, Astronomy & Astrophysics, 596, A57
  * DeMeo & Carry (2014) DeMeo F. E., Carry B., 2014, Nature, 505, 629
  * DeMeo et al. (2009) DeMeo F. E., Binzel R. P., Slivan S. M., Bus S. J., 2009, Icarus, 202, 160
  * Denneau et al. (2013) Denneau L., et al., 2013, Publications of the Astronomical Society of the Pacific, 125, 357
  * Dobson et al. (2021) Dobson M. M., et al., 2021, Research Notes of the AAS, 5, 211
  * Dobson et al. (2023) Dobson M. M., et al., 2023, The Planetary Science Journal, 4, 75
  * Donaldson et al. (2023) Donaldson A., Kokotanekova R., Rożek A., Snodgrass C., Gardener D., Green S. F., Masoumzadeh N., Robinson J., 2023, Monthly Notices of the Royal Astronomical Society, 521, 1518
  * Ďurech et al. (2020) Ďurech J., Tonry J., Erasmus N., Denneau L., Heinze A. N., Flewelling H., Vančo R., 2020, Astronomy & Astrophysics, 643, A59
  * Emery et al. (2011) Emery J. P., Burr D. M., Cruikshank D. P., 2011, The Astronomical Journal, 141, 25
  * Erasmus et al. (2020) Erasmus N., et al., 2020, The Astrophysical Journal Supplement Series, 247, 13
  * Erasmus et al. (2021) Erasmus N., et al., 2021, Monthly Notices of the Royal Astronomical Society, 506, 3872
  * Fernández-Valenzuela (2022) Fernández-Valenzuela E., 2022, Frontiers in Astronomy and Space Sciences, 9, 796004
  * Gomes et al. (2005) Gomes R., Levison H. F., Tsiganis K., Morbidelli A., 2005, Nature, 435, 466
  * Gorski et al. (2005) Gorski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, The Astrophysical Journal, 622, 759
  * Grav et al. (2011) Grav T., et al., 2011, The Astrophysical Journal, 742, 40
  * Harris & Young (1989) Harris A. W., Young J. W., 1989, Icarus, 81, 314
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Heinze et al. (2018) Heinze A. N., et al., 2018, The Astronomical Journal, 156, 241
  * Heinze et al. (2021) Heinze A. N., et al., 2021, The Planetary Science Journal, 2, 12
  * Hendler & Malhotra (2020) Hendler N. P., Malhotra R., 2020, The Planetary Science Journal, 1, 75
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Jackson et al. (2022) Jackson S. L., Rozitis B., Dover L. R., Green S. F., Kolb U. C., Andrews A. E., Lowry S. C., 2022, Monthly Notices of the Royal Astronomical Society, 513, 3076
  * Jewitt & Hsieh (2022) Jewitt D., Hsieh H. H., 2022, The Asteroid-Comet Continuum, doi:10.48550/arXiv.2203.01397
  * Jurić et al. (2002) Jurić M., et al., 2002, The Astronomical Journal, 124, 1776
  * Kaiser et al. (2010) Kaiser N., et al., 2010, in Stepp L. M., Gilmozzi R., Hall H. J., eds, SPIE Astronomical Telescopes + Instrumentation. San Diego, California, USA, p. 77330E, doi:10.1117/12.859188
  * Keys et al. (2019) Keys S., et al., 2019, Publications of the Astronomical Society of the Pacific, 131, 064501
  * Knuth (1998) Knuth D. E., 1998, The Art of Computer Programming, third edn. Computer Algorithms Vol. 2, Addison-Wesley
  * Kokhirova et al. (2021) Kokhirova G., Ivanova O., Rakhmatullaeva F.Dzh., Baransky A., Buriev A., 2021, Advances in Space Research, 67, 639
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, LSST Science Book, Version 2.0 (arxiv:0912.0201), doi:10.48550/arXiv.0912.0201
  * Lagerkvist & Magnusson (1990) Lagerkvist C.-I., Magnusson P., 1990, Astronomy and Astrophysics Supplement Series, 86, 119
  * Law et al. (2009) Law N. M., et al., 2009, Publications of the Astronomical Society of the Pacific, 121, 1395
  * Levison et al. (2021) Levison H. F., et al., 2021, The Planetary Science Journal, 2, 171
  * Li et al. (2022) Li J., Xia Z. J., Yoshida F., Georgakarakos N., Li X., 2022, Asymmetry in the Number of L4 and L5 Jupiter Trojans Driven by Jumping Jupiter (arxiv:2211.13877)
  * Mahlke et al. (2021) Mahlke M., Carry B., Denneau L., 2021, Icarus, 354, 114094
  * Marsden (1985) Marsden B. G., 1985, Minor Planet Circulars, 10193, 1
  * Marzari et al. (2002) Marzari F., Scholl H., Murray C., Lagerkvist C., 2002, Origin and Evolution of Trojan Asteroids
  * McKinney (2010) McKinney W., 2010, in Python in Science Conference. Austin, Texas, pp 56–61, doi:10.25080/Majora-92bf1922-00a
  * McNeill et al. (2021) McNeill A., et al., 2021, The Planetary Science Journal, 2, 6
  * Milani et al. (2014) Milani A., Cellino A., Knežević Z., Novaković B., Spoto F., Paolicchi P., 2014, Icarus, 239, 46
  * Mommert et al. (2019) Mommert M., et al., 2019, Journal of Open Source Software, 4, 1426
  * Morbidelli et al. (2005) Morbidelli A., Levison H. F., Tsiganis K., Gomes R., 2005, Nature, 435, 462
  * Moskovitz et al. (2022) Moskovitz N. A., Wasserman L., Burt B., Schottland R., Bowell E., Bailen M., Granvik M., 2022, The Astorb Database at Lowell Observatory (arxiv:2210.10217)
  * Muinonen et al. (2002) Muinonen K., Piironen J., Shkuratov Yu. G., Ovcharenko A., Clark B. E., 2002, Asteroid Photometric and Polarimetric Phase Effects
  * Muinonen et al. (2010) Muinonen K., Belskaya I. N., Cellino A., Delbò M., Levasseur-Regourd A.-C., Penttilä A., Tedesco E. F., 2010, Icarus, 209, 542
  * Nagdimunov (2021) Nagdimunov L., 2021, pds4-tools, https://pypi.org/project/pds4-tools/
  * Nesvorný et al. (2013) Nesvorný D., Vokrouhlický D., Morbidelli A., 2013, The Astrophysical Journal, 768, 45
  * Oszkiewicz et al. (2011) Oszkiewicz D., Muinonen K., Bowell E., Trilling D., Penttilä A., Pieniluoma T., Wasserman L., Enga M.-T., 2011, Journal of Quantitative Spectroscopy and Radiative Transfer, 112, 1919
  * Oszkiewicz et al. (2012) Oszkiewicz D. A., Bowell E., Wasserman L. H., Muinonen K., Penttilä A., Pieniluoma T., Trilling D. E., Thomas C. A., 2012, Icarus, 219, 283
  * Pedregosa et al. (2011) Pedregosa F., et al., 2011, Journal of Machine Learning Research, 12, 2825
  * Penttilä et al. (2016) Penttilä A., Shevchenko V., Wilkman O., Muinonen K., 2016, Planetary and Space Science, 123, 117
  * Pravec & Harris (2007) Pravec P., Harris A., 2007, Icarus, 190, 250
  * Pravec et al. (2002) Pravec P., Harris A. W., Michalowski T., 2002, Asteroid Rotations
  * Pravec et al. (2012) Pravec P., Harris A. W., Kušnirák P., Galád A., Hornoch K., 2012, Icarus, 221, 365
  * Schechter et al. (1993) Schechter P. L., Mateo M., Saha A., 1993, Publications of the Astronomical Society of the Pacific, 105, 1342
  * Schemel & Brown (2021) Schemel M., Brown M. E., 2021, The Planetary Science Journal, 2, 40
  * Shevchenko (1997) Shevchenko V. G., 1997, Solar System Research, 31, 219
  * Shingles et al. (2021) Shingles L., et al., 2021, Transient Name Server AstroNote, 7, 1
  * Smartt et al. (2017) Smartt S. J., et al., 2017, Nature, 551, 75
  * Smith et al. (2020) Smith K. W., et al., 2020, Publications of the Astronomical Society of the Pacific, 132, 085002
  * Spinner (2022) Spinner T., 2022, pycolormap-2d, https://pypi.org/project/pycolormap-2d/
  * Szabo et al. (2007) Szabo G. M., Ivezic Z., Juric M., Lupton R., 2007, Monthly Notices of the Royal Astronomical Society, 377, 1393
  * Tonry et al. (2018) Tonry J. L., et al., 2018, Publications of the Astronomical Society of the Pacific, 130, 064505
  * Tsiganis et al. (2005) Tsiganis K., Gomes R., Morbidelli A., Levison H. F., 2005, Nature, 435, 459
  * Uehata et al. (2022) Uehata K., Terai T., Ohtsuki K., Yoshida F., 2022, The Astronomical Journal, 163, 213
  * Usui et al. (2013) Usui F., Kasuga T., Hasegawa S., Ishiguro M., Kuroda D., Müller T. G., Ootsubo T., Matsuhara H., 2013, The Astrophysical Journal, 762, 56
  * Vereš et al. (2015) Vereš P., et al., 2015, Icarus, 261, 34
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Warner et al. (2021) Warner B. D., Harris A. W., Pravec P., 2021, NASA Planetary Data System, p. 10
  * Waszczak et al. (2015) Waszczak A., et al., 2015, The Astronomical Journal, 150, 75
  * Wong & Brown (2016) Wong I., Brown M. E., 2016, The Astronomical Journal, 152, 90
  * Young (2023) Young D. R., 2023, rockAtlas, doi:10.5281/zenodo.8037828, https://github.com/thespacedoctor/rockAtlas
  * Zonca et al. (2019) Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, Journal of Open Source Software, 4, 1298

##  Appendix A Metrics for high quality phase curves