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

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 × 10 7 observations of 4.5 × 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, ∼ 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 ∼ 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 ∼ 90 degrees.


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,   , through photometric methods at a specific point in time then allows for an estimate of the absolute magnitude,   , by fitting those observations to the formula   =   − 5 log (Δ) −   () (1) Here  and Δ are the heliocentric and geocentric distances in au and  is the Sun-asteroid-Earth phase angle, all calculated at the time   was measured.The phase function,   (), mathematically describes the macroscopic scattering law of sunlight from the asteroidal surface as a function of .For all models of asteroid light scattering,   () 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 ★ E-mail: james.robinson@ed.ac.uk are believed to be the mechanisms of shadow hiding and coherent backscatter (e.g.Muinonen et al. 2002Muinonen et al. , 2010;;Penttilä et al. 2016).
Therefore by definition the absolute magnitude,   , of a spherical asteroid is its magnitude when observed at 1 au from both the Sun and the Earth, whilst at opposition ( = 0 • ).In reality most asteroids are non-spherical, and so the measured   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   and   () 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   () 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 Terrestrialimpact Last Alert System; Tonry et al. 2018).All of these surveys use different filters, which means the reported   requires transformation to a common standard filter, normally Johnson -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  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   () 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   ().
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 -band for all large surveys, but they often do not correct for the rotational/apparition effects listed above.Given the ∼ 4 × 10 8 observations in these databases from different facilities with a range of photometric accuracy and bandpasses, this is not surprising.Instead, absolute magnitudes   are calculated by transforming to the -band and fitting the parameterised (, ) system of (Bowell et al. 1989), where the form of   () depends only on a phase parameter .  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  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   , which in turn allows determination of the diameter   km =   10 −  /5 √ (2) where   is geometric albedo and   is a wavelength dependent constant related to the absolute magnitude of the Sun; e.g. for band observations   = 664.5 km (Pravec & Harris 2007).In the past it has been shown that significant biases can exist for   (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,   () is also important in the study of asteroid properties.The form of   () 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   () beyond establishing   .First, this can indicate approximate compositions for asteroids without multicolour 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. (2011Oszkiewicz et al. ( , 2012) ) debiased the existing MPC catalogue for interobservatory offsets in reported magnitudes, deriving absolute magnitudes and phase functions for approximately 5 × 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 × 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×10 7 photometric measurements of 4.5×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.

ATLAS observations
Until recently, the ATLAS survey was composed of two 0.5 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 predetermined 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 m aperture and 0.65 m spherical primary mirror (Tonry et al. 2018).The detector is an STA-1600 10.5 × 10.5 k CCD with a pixel scale of 1.86 ′′ .ATLAS primarily surveys in an "orange"  and "cyan" -filter, reaching down to ∼ 19.5 mag for its typical 30 s exposures.These filters are approximately equivalent to PS1 filters  +  (560 − 820 nm) and  +  (420 − 650 nm) respectively.Generally the bluer -filter is reserved for new moon dark time during survey operations, whereas the redder -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 • <  < +90 • for the Hawaii units), revisiting each field 4 times over timescales ≤ 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 hr 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 1 ).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 1 https://www.minorplanetcenter.net/cgi-bin/checkmp.cgi 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.

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 ∼ 445, 000 asteroids from the astorb2 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 pipelinedophot (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 pyephem3 and orbfit4 .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 arcsec, 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 ∼ 9.6 × 10 7 observations as of January 2022, comprised of ∼ 65% in the -filter and ∼ 35% in the -filter.Asteroids in the rockAtlas database have a median of ∼ 130 observations.Apparent magnitudes range from ≳ 12 mag down to a limiting magnitude of ∼ 19.5, with a median observed magnitude of ∼ 17 mag.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   > 17 mag.
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 ( < 1.3 AU) and Outer Solar System Objects ( >  Jup = 5.457 AU) observed with ATLAS (Heinze et al. 2021;Dobson et al. 2021Dobson et al. , 2023)).The median maximum phase angle of observations is ∼ 20 • , 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 ∼ 20 • whereas the most distant objects we consider in this work, the Jupiter Trojans, are observed only up to  ∼ 10 • .Observations of asteroid brightness at low phase angles (< 5 • ) 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  ∼ 3 • , 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 database 5 .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.

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 (, ) 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,  (), which is the apparent magnitude corrected for distance effects such that only the phase function   () from equation 1 leads to deviation from the absolute magnitude In the B89 system the change in reduced magnitude with respect to phase angle is described by where  is the absolute magnitude of the object (brightness at zero phase scaled to a distance of 1 au), Φ B89, are the basis functions, and  is the phase parameter that controls the form of the model phase curve.
The second model we use is the (,  * 12 ) system from Penttilä et al. (2016), hereafter referred to as P16.This was derived from the new IAU standard (,  1 ,  2 ) system formulated by Muinonen et al. (2010), which more accurately reflects the phase function at small phase angles and is defined as where Φ P16, are the basis functions in this system.In particular Φ P16,3 models the opposition effect.There are two phase parameters  1 ,  2 and in the Penttilä et al. (2016)  ).The P16 model is better suited for modelling the phase curves of objects with lower quality, more sparse observations.The phase parameter  * 12 must not be confused with the similar  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  * 12 →  12 for simplicity.Likewise, from here on we use  12 to distinguish the P16 absolute magnitude from the B89 absolute magnitude .

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  and  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  (B89) and  12  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  > 1 arcsec to reduce the number of false matches.In addition we rejected any observations with low galactic latitude,  ≤ 10 • , 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,   ≤ 0.005 mag, 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  > 1 arcsec), therefore this tighter constraint is worth possibly rejecting the small number of true matches with separations 1 <  (arcsec) < 5.
We then performed additional filtering of outlying data to each individual apparition, considering the  and -filters separately.To each apparition we fitted the B89  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  was taken from astorb and colour corrected from -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.
The slope parameter  was set at  = 0.15 which is the standard assumed value for all asteroids (Marsden 1985).Using this  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  we account for possible apparition effects which would increase the residuals if all observations were considered together.If an asteroid displays large shifts in  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  < 5 • 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 -filter as the ATLAS dataset preferentially contains observations in  rather than .
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  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  and  are fitted for  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  and  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  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 ∼ 5000 asteroids out of an initial sample of ∼ 180, 000 objects.Furthermore Erasmus et al. ( 2021) searched ATLAS observations of ∼ 250, 000 asteroids for extreme slow rotators (period > 1000 hours) 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.

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 , and its uncertainty   for the B89 fits.Visual inspection of individual phase curves showed that the most extreme (and probably unphysical) values of  and   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 <  < 1, without us having to perform a hard cut on  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,  < 5 • , after outlier rejection.All other apparitions required at least 10 observations to be considered reliable.After developing our metric cuts for the  model, we then applied them in the same manner to the  12  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.

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  and the phase curve residuals to identify significant variations across several apparitions.We define a parameter which is the ratio of the standard deviation of  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  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 -filter data in this test due to the larger number of photometric measurements compared to the -filter.When the changes in  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  app ≥ 0.9 succeeds in identifying the objects with the strongest apparition effects.The remaining asteroids with  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  app ≥ 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  of all retained apparitions in either the  or -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 , 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 ( app ≥ 0.9) from our results when discussing absolute magnitude and colour, although we make use of the phase parameters ,  12 of these objects when  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 -filter.The majority of asteroids observed by ATLAS also have a smaller number of measurements in the bluer -filter.Thus we can compare the estimates of  in each filter to estimate a   −   colour for a large sample of asteroids.Measurements in the -filter suffer from the same systematic errors discussed above for the -filter.Therefore we calculated colour only for objects with apparitions that satisfy the quality control requirements in both  and -filters, excluding any objects displaying strong apparition effects in the -filter ( app ≥ 0.9) as we cannot define a unique  for these objects.The colour was then calculated as median(  ) -median(  ).The ATLAS colour provides an additional diagnostic for comparing the surface properties of different asteroid populations.

ATLAS phase curve database
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   −   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.
In figure 4 we present the distributions of   ,  and   −   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  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  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  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  decreases the slope of the approximately linear component at phase angles  ≳ 5 • gets steeper.In particular the outliers with  ≲ 1 have phase curves that are entirely compatible with the observations.On the other hand, as  increases the slope at larger phase angles decreases and the model fails to describe the opposition effect at phase angles < 5 deg.Our handful of objects with  ≳ 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  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  for the S-type asteroids leads to a flatter phase function, which corresponds to a higher albedo as expected.
There is a correlation between  and the range of  we find in the data (figure 5).This is partly due to the enhanced presence of S-type asteroids as  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 .However the main cause is a systematic effect; as  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 .This - correlation driven by photometric uncertainty acts in addition to any uncertainty caused by rotational lightcurve effects.For example, analysis of the  dataset published by SB21 shows a similar trend.
Figure 4 shows that there is a strong bimodality for the   −   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  app < 0.9).The bimodality is again driven by the different surface compositions of the dominant S and C type asteroids.We obtain median   −   colours of 0.39 ± 0.05, 0.27 ± 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   −   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   −  colours of 0.388 and 0.249 for S and C-type asteroids respectively.Figure 6 shows the 2-D bimodal distribution of   −   vs.  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).
We have also derived the phase curve parameters of asteroids using the P16  12  12 system, where we use the notation  12 to denote absolute magnitudes in this system.As described in section 2.5 the quality control cuts were developed on the  dataset and then applied separately to the residuals of the  12  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  do not lead to as strong restrictions on  12 for "physical" values between 0 and 1.This appears to be mostly due to inherent differences between the  and  12  12 models.Figure 8 shows that there is a strong correlation between  12 and  for objects in our dataset with a phase curve fit    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 ( app ), where the first observation in an apparition was made at time  app,0 .The maximum phase angle observed in an apparition is given by max(  app, ) and max(  app, ) for the  and -filters respectively.The number of detections used to fit the phase curve (after outlier rejection) is denoted by  fit, and  fit, .We show the fitted phase curve parameters for the  B89 system, alongside their uncertainties (  ,   ) in both  and -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  filter ("Keep  ") or if it passed in both  and  ("Keep − ").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.provides the number of apparitions in the  filter that were used and  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 ( app,, ), for which we could calculate an asteroid colour   −   .Available online.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  for objects with fainter absolute magnitudes.This means that applying the P16 model results in lower values of  (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 .In the test described above the exact  phase curve extends all the way to zero phase.The  12  12 model cannot match  for both the approximately linear behaviour at large phase angles and the nonlinear opposition brightening at low phase angles (figure B1).The  and  12  12 models show the most significant divergence at low phase angles which explains the differences in absolute magnitude between the two models.

Name
We show the relation between  12 −  and phase parameter  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 B1 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  difference is systematic across both  and -filter phase curve fits, hence the resulting   −   colours for both models are practically indistinguishable.In the P16 system we find median   −   colours that are nearly identical to the B89 system: 0.39 ± 0.05 and 0.27 ± 0.06 for S and C-types respectively (uncertainties of one standard deviation).

Comparison with previous studies
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 ∼ 180, 000 asteroids in the period of 2015-2019.From this they determined phase curves with the Muinonen et al. ( 2010)  1  2 and P16  12  12 models for ∼ 94000 asteroids using MCMC Bayesian parameter inference, with fits in both  and  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 <  12 < 1 to ensure monotonically increasing phase curves.This leads to a pile up of  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  12 (figure 10).The median difference in absolute magnitude  12,Mahlke −  12 = −0.01 ± 0.08 mag (one standard deviation) and the difference in  12 is also relatively near zero with a median  12,Mahlke −  12 = −0.06albeit with a higher standard deviation of 0.31.
In tational 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  (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  parameter and  −  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  McNeill −  = 0.06, with quite a large standard deviation of 0.28 (figure 11).Similarly for the difference in   −   colours, the median difference is −0.008 ± 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 <  < 2 for which we obtain more realistic values of .
Considering another study focused on the Jupiter Trojans, SB21 determined  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 ,  and  −  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  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.
There is reasonable agreement between  given that we have not accounted for rotational effects, the median difference in  Schemel − is 0.07±0.13.However, there is a significant offset in  distributions: −0.13 ± 0.18 (figure 11).This means that in general our work finds larger values of  for Jupiter Trojans, which corresponds to flatter phase curves (at phase angles > 5 deg 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  when such scatter is  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 -axis.
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 B1).
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  and -filters, resulting in a dataset of ∼ 9000 objects.They determined phase curves in the two parameter Shevchenko (1997)  parameters to ours for 5199 MBAs and Jupiter Trojans common to both datasets (figure 11).The median difference in  and its standard deviation, −0.07 ± 0.20, are of similar magnitude to the difference with McN21.However both Waszczak et al. (2015) and SB21 Δ 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  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  to higher values.

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   changes from 16.6 mag at the inner edge of the main belt at  = 1.9 au to 13.0 mag for the Hilda asteroids at  = 4 au, and to 11.6 mag for the Jupiter Trojan asteroids at  = 5.2 au.The median uncertainties of the non-linear least squares fits for all phase curves in the SAP database are   = 0.04 and  12 = 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.

Asteroid families in the SAP database
To illustrate the usefulness of the SAP database, in figure 13 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 -axis.
(obtained from the astdys database6 ), 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  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  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   −   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   −   colour and phase parameter  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   −   vs.  to colour each marker.This colour map is displayed in the right-hand panel of figure 14 and it illustrates that this combination of   −   and  can provide a powerful tool for identifying broad taxonomic information of asteroids observed by sparse, wide field where the colour transition is centred on the median value of .The same colour scheme is used for both plots and objects with phase parameters outside the 1 − 99% quantile of the family member  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 <  < 0.60 (where the median  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).
surveys.We have indicated the distribution of Nysa-Polana members in this   −  vs.  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.

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 AU (Marzari et al. 2002).The Jupiter Trojans are distributed in two groups that are leading/trailing Jupiter by approximately 60 degrees, 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 -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  = 14 mag, Hendler & Malhotra 2020), and determination of the underlying size distribution indicates a number asymmetry of ∼ 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 (max() ≤ 11 • ) 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 deg.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  fits in the -filter, of which 570 also have good fits in the -filter suitable for calculating a   −   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 ∼ 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, , of these two groups (figure 15, left panels) using a two sample Kolmogorov-Smirnov (KS) test.In a two sample KS test the parameter  describes the maximum distance (at a given value) between the cumulative distributions of the two datasets being 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,   −   colour against phase parameter .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   −   and  of the S and C-type asteroids in the full SAP database (red cross and plus respectively, values from figure 4).
tested.The KS statistic is calculated with a -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  < 0.05 (see also equation C1).For the L4/L5  distributions this test results in a KS statistic of  = 0.2 (with a -value of 3 × 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  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   from the full astorb dataset.However when we compare the astorb   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  the two sample KS test results are  = 0.08 and  = 0.18.Likewise for   −   colours (for the smaller number of objects with good fits in both  and ), the KS test results are  = 0.08 and  = 0.39 (figure 15, middle and right panels respectively).These  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  for the resampled data, the results of which are detailed further in appendix C (figure C1).Overall these tests confirm that there is no statistical evidence for any differences in the distribution of phase parameter  or   −   colour between the L4 and L5 Jupiter Trojans.This is in contrast to the reported dissimilarity in − 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  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  and   −   in addition to the greater phase curve fitting uncertainty for fainter objects (figure 5).
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  and  distributions seemed to indicate statistically significant differences (high  values with low ).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 C2).
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  is also present in our dataset.In this figure we have plotted the phase parameter  obtained from this work against the  −  colour derived from ZTF data by SB21, using only objects with better than median uncertainties in both datasets (replicating their analysis, see figure 6  The KS statistic  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 -value) we can reject the null hypothesis that the two datasets are drawn from the same underlying distribution.
to both datasets, 159 meet this requirement.Similar to the previous work we find a strong positive correlation between our  values and their  −  colours, with a Spearman rank correlation coefficient of  = 0.51 ( = 7 × 10 −12 ).We then divided objects into red and lessred groups by a ZTF colour boundary of  −  = 0.565 and analysed the corresponding distributions of .We replicate the results of SB21 whereby less-red Jupiter Trojans generally have lower  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 ( = 0.48,  = 2 × 10 −8 ).
We also search for correlation between  and   −   colour for the 570 Jupiter Trojans with colours in the ATLAS dataset.The Spearman rank correlation coefficient of these parameters is  = 0.07 ( = 0.08), which indicates that there is no significant correlation.This is due to the large width of the ATLAS  and -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 sig-nificant differences in phase parameter  or   −   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   −   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  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.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  −  = 0.565 which divides the "red" and "less-red" Jupiter Trojan compositional classes The lower panel shows the histogram distributions of ATLAS  for the red and less-red Jupiter Trojans from the panel above.

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 mag 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 • 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 • 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  < 1 ′′ 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.

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  and   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.

Phase curve wavelength dependence
We determined the phase parameter (,  12 ) for each object from its primary apparition in the -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 -filter is much sparser compared to the -filter.In most cases we lacked sufficient data to obtain an independent phase curve for the -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   and   .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.

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.

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 × 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  = 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.

CONCLUSIONS
We have compiled a database of ATLAS asteroid photometry by matching the ephemerides of 4.5×10 5 asteroids with accurate or-bital elements to the detections produced by the ATLAS dophot pipeline, using rockAtlas (Young 2023).We fitted phase curves in the ,  12  12 models (Bowell et al. 1989;Penttilä et al. 2016) to the apparition in the -filter with the best combination of number of observations and phase angle coverage.We then fitted all other apparitions for  with these derived phase parameters fixed.This allowed us to assess any changes in  during the survey and also to determine the   −   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  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  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 .Our dataset illustrates the known trend towards lower values of  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   −   colour and  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  and   −   .We did measure a discrepancy in the distribution of  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  and -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  −  colours; Schemel & Brown 2021) the difference in ATLAS phase parameter  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×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 sHG 1 G 2 model (Carry et al. 2024); i.e. a "spinned" version of the  1  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.

APPENDIX A: METRICS FOR HIGH QUALITY PHASE CURVES
The temporal and phase angle coverage of an asteroid imaged by ATLAS will vary depending on many factors, e.g.accessible viewing geometry relative to the Earth, sky location relative to bright objects (such as the Moon, Jupiter and Saturn) and/or missed observations due to weather/technical issues.Likewise the quality of observations obtained by dophot will depend on factors such as object brightness, photometric seeing and/or sky brightness.Furthermore, the matching of dophot detections to known asteroid ephemerides may lead to false matches with stellar sources in the high density galactic plane or glints near bright objects and other such artifacts.Therefore, poor quality phase curves can occur when there is not enough observational coverage to adequately constrain the phase curve, photometric uncertainty is high and/or there is significant outlying photometric observations.
We performed data cuts on the observations of each asteroid to account for obvious photometric outliers (section 2.4), although our methods required sufficient data in each apparition such that the initial  ( = 0.15) fit provided a reasonable description of the general phase curve behaviour.Considering the individual apparitions of each asteroid reduces the number of observations used for each phase curve fit, as such after performing our phase curve fitting algorithm there will many apparitions which do not produce accurate phase curve fits.Therefore we used a series of cuts based on various metrics describing the observational coverage and residuals of each phase curve fit to reject poor quality apparitions from our final dataset.By analysing the residuals of the phase curve fit at each apparition we were able to identify and exclude individual apparitions for which an accurate phase curve could not be obtained.Therefore one bad apparition did not necessarily lead to complete rejection of an object.The primary apparition provides us with the fitted phase parameter ( or  12 ) for each asteroid, therefore these apparitions had somewhat stricter requirements on the number of data points and phase angle coverage.The phase parameter was kept fixed when fitting the non-primary apparitions.For these fits the most common phase curve error was unphysical shifts in  caused by outlying photometry missed by our data selection, poor phase angle coverage, and/or large rotational variations driven by asteroid shape.We found that these problematic apparitions were best characterised by the statistical properties of their residuals.For a high quality phase curve fit of an asteroid with low rotational amplitude we would expect the residuals to be approximately normally distributed around a residual of 0 with low standard deviation 7 .As such the standard deviation and mean value of the residuals proved to be valuable metrics.We also performed a Kolmogorov-Smirnov test to identify distributions of 7 In reality the likelihood function of brightness for a sinusoidal lightcurve peaks at the maximum and minimum amplitude (see figure 1 of SB21) although this would not be clearly seen for low amplitude lightcurves given the typical accuracy of ATLAS photometry.residuals that differed significantly from a normal distribution.Furthermore, we used the Pearson correlation coefficient of the residuals as a function of phase angle as a metric.This metric should ideally be near zero (indicating no correlation) but will approach 1 or −1 when there is strong correlation in the residuals.Correlated residuals implies that the phase curve model was not consistent with the observations.
The values of each metric cut were chosen through a process of trial and error alongside detailed visual inspection of individual B89 phase curve fits.As far as possible we tried to avoid choosing values in a way that directly influenced the final distributions of phase curve parameters, i.e. we did not simply cut phase curves with  < 0 or  > 1.We did however choose values that led to a distribution such that the vast majority (97.7%) of  fits had 0 <  < 1.This also removed the objects with extremely high uncertainties in , while not removing an excessive number of objects from the final database.We did have to perform a hard cut in   ,   as a final step, as there was no other combination of metric cuts that could consistently remove obvious poor quality phase curves.This approach was justified by detailed visual inspection of many phase curves and their residuals from which it was clear that phase curves with extreme  and   were unphysical, and as expected were caused by a combination of poor coverage, bad photometry, and/or unaccounted rotational variation.
The cuts on the primary apparition (which was always in the filter) for each asteroid were as follows: • In order to compare the properties of the  (B89) and  12  12 (P16) phase curve systems we generated model  phase curves across a phase angle range of 0 → 25 deg, which is the typical phase angle range for asteroids observed by ATLAS (figure 1).We generated a number of phase curves with slope parameters 0 <  < 1.We then performed a linear least squares fit of the  12  12 model to each  phase curve, the results of which are shown in figure B1.Both models are in general agreement at larger phase angles where brightening is approximately linear, however, the differing treatment of the opposition brightness surge (through their basis functions, equations 4 and 5) by the two models is clear as phase angles decreases.From this simple analysis we obtain an empirical mapping of  →  12 which provides a good description of the relation between  and  12 for the SAP dataset (figure 8).We .Plots of B89 model  phase curves (blue curves), which were generated with a range of 0 ≤  ≤ 1.Each  phase curve was then fitted with the P16  12  12 model (orange curves).The two models have best agreement at large phase angles, however, they vary strongly at lower phase angles leading to differences in absolute magnitude.performed a 3rd-order polynomial fit for each relation, the parameters of which are as follows  12 = −0.7451 3 + 2.577 2 − 3.716 + 1.124 (B1a) This behaviour also explains why our phase curve fits for  12  12 more frequently stray into "unphysical" solutions, i.e. solutions outside of 0 ≤  12 ≤ 1 (figure 7).Furthermore, these model fits clearly illustrate the origin of the different absolute magnitudes obtained with the two phase curve systems (figure 8), where the relation between Δ =  12 −  and  is reflected well in the ATLAS data.This idealised test shows that even if there are observations all the way down to zero phase the fitted  12  12 model will differ significantly from the  fit.Therefore Δ in our dataset is not primarily caused by a lack of information at low phase angles, rather an inherent difference between the two models.

APPENDIX C: JUPITER TROJAN KS BOOTSTRAP TESTS
Here we present the results of performing KS tests on bootstrapped data (resampled with replacement) in order to analyse the distributions of Jupiter Trojan parameters ,  and   −   .For each parameter, 1000 bootstrapped KS tests were performed and the distribution of the resulting  values is shown in figure C1.In these figures we compare the median  value for each set of KS tests and compare this to the critical value of  required for rejecting the null hypothesis at a significance of 0.05.where  and  are the number of data points in the two samples being considered, and the constant of 1.36 refers to a significance of 5% (Knuth 1998).As described in section 4.3 we found that this method was is more accurate than relying on the results of a single KS test of what is a relatively small set of objects.We confirm differences in the distribution of  between the L4 and L5 groups as observed by ATLAS; the distribution of  for resampling tests is consistently higher than  crit (figure C1, left panel).In contrast, although some tests for slope parameter and colour exceed  crit , the distributions of  are more often below or near  crit (figure C1, middle and right panels).As such this indicates that any differences in   −   colour and phase parameter  are not statistically significant enough to reject the null hypothesis.Furthermore, we repeated this analysis for the values of ,  and  −  colour published by SB21 (figure C2).These results show that differences in the distributions of  for the L4/L5 groups is not present in their analysis of ZTF data.This is as expected given that the initial rockAtlas watchlist used in this study started with different  distributions of the L4/L5 Trojans (section 4.3).Similar to our results there are no statistically significant differences between L4/L5 for  or  −  in ZTF data (SB21)

Figure 1 .
Figure 1.Two dimensional histogram showing the number distribution of the maximum observed phase angle of asteroids in the rockAtlas database (filter) as a function of the object's semimajor axis (left-hand panel).The maximum phase angle for a zero-inclination circular orbit, max( ) = sin −1 1/, 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 ≳ 20 • .

Figure 2 .
Figure2.ATLAS phase curves of asteroid (766) Moguntia, an object displaying a large apparition effect.Observations in the -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 -axis have been set to focus on retained photometry, ignoring any more extreme outliers.The corresponding  phase curve fit for each apparition is shown by the matching coloured curves, with the primary apparition indicated by the coloured star.The median  of the retained apparitions is indicated by the dashed black curve with parameters:   = 9.68,  = 0.21.Observations in the -filter are shown in the right-hand panel, where marker/line colour have the same meanings as the left panel.The large shifts in  between apparitions are clear ( app = 5.9; equation 8).

Figure 3 .
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  and -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  indicated by the black dashed curve with parameters:   = 13.21, = 0.30

Figure 4 .Figure 5 .
Figure 4. Log-scale distributions of  (left),  (middle) and   −   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). and  are given in the -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:   (S) = 14.83,   (C) = 14.52,  (S) = 0.29,  (C) = 0.21,   −   (S) = 0.39,   −   (C) = 0.27.In the middle panel we indicate the boundaries on  for a monotonically increasing B89 phase curve with vertical red dashed lines.

Figure 6 .
Figure 6.Two dimensional histogram of   −   colour vs. .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 subdistributions in   −   vs.  parameter space.The median values of the S and C-types with ATLAS colours are  = 0.29, 0.21 and   −   = 0.39, 0.27 mag respectively.

Figure 7 .
Figure 7. Log-scale distribution of P16  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  12 (S) = 0.22,  12 (C) = 0.51) and the limits 0 <  12 < 1 for a "physical" phase curve are shown as vertical red dashed lines.

Figure 8 .
Figure 8. Two dimensional histogram showing the relation between the  12 and  parameters for objects in our database (upper panel).A 3rd-order polynomial is fitted to the data as shown by the black curve (equation B1a).The red markers indicate the values of  12 when a  12  12 model is fitted to a  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  12 −  vs. .Again, the black curve is a 3rd-order polynomial fit (equation B1b) and the red markers indicate the expected difference in absolute magnitude when  12  12 is fitted to .The horizontal black line indicates the zero value of the  12 −  axis.

Figure 9 .
Figure 9. Two dimensional histogram which shows the variation of Δ = 12 −  with .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 Δ and , as well as a trend for a larger dispersion in Δ for fainter objects (see also figure8).The horizontal black line indicates the zero value of the  12 −  axis.

Figure 10 .
Figure10.Histograms showing the difference in  (upper panel) and  12 (lower panel) of our phase curve fits in the P16 system, cross-matched with the results ofMahlke 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 -axis.

Figure 11 .
Figure 11.Histogram distributions showing the difference in  (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 -axis ticks indicate the median value of the corresponding distribution, whereas the vertical black line indicates the zero value of the -axis.The median offset for each distribution is  Schemel −  = −0.13, McNeill −  = 0.06,  Waszczak −  = −0.07.

Figure 12 .
Figure 12.Histogram comparing the ATLAS   −   colour distribution of objects in this work cross-matched with the previous colours ( − ) 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 -axis.

Figure 13 .
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  with a diverging colour map,where the colour transition is centred on the median value of .The same colour scheme is used for both plots and objects with phase parameters outside the 1 − 99% quantile of the family member  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 <  < 0.60 (where the median  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 .
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,   −   colour against phase parameter .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   −   and  of the S and C-type asteroids in the full SAP database (red cross and plus respectively, values from figure4).

Figure 15 .
Figure 15.Histogram distributions of ,  and   −   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  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 -value) we can reject the null hypothesis that the two datasets are drawn from the same underlying distribution.

Figure 16 .
Figure16.Scatter plot of phase parameter  from the SAP dataset against the ZTF  −  colours fromSchemel & 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 figure6).The red vertical line corresponds to the colour boundary at  −  = 0.565 which divides the "red" and "less-red" Jupiter Trojan compositional classes The lower panel shows the histogram distributions of ATLAS  for the red and less-red Jupiter Trojans from the panel above.
Number of data points (after outlier rejection) ≥ 50 • ≥ 1 observation at phase angle < 5 deg • Phase angle range; max( app ) − min( app ) ≥ 5 deg • Absolute magnitude uncertainty   (or  12 )< 0.05 • Phase parameter uncertainty   (or  12 ) < 0.2 All apparitions (primary and non primary,  and -filters) had to satisfy the following requirements: • Number of observations (after outlier rejection) ≥ 10 • Phase angle ratio (the ratio of phase angle range of this apparition to the phase angle range of all observations in the same filter for the asteroid) ≥ 0.3 • Standard deviation of residuals < 0.4 • Mean value of residuals < 0.2 • Absolute value of the Pearson correlation coefficient of residuals < 0.25 • Kolmogorov-Smirnov statistic (for normal distribution) < 0.3 APPENDIX B: COMPARING  AND  12  12 Figure B1.Plots of B89 model  phase curves (blue curves), which were generated with a range of 0 ≤  ≤ 1.Each  phase curve was then fitted with the P16  12  12 model (orange curves).The two models have best agreement at large phase angles, however, they vary strongly at lower phase angles leading to differences in absolute magnitude.

Figure C1 .Figure C2 .
Figure C1.Distributions of  values for KS tests comparing ,  and   −   colour of the L4 and L5 Jupiter Trojans in the SAP dataset (left, middle and right panels respectively; compare to figure 15 CDFs).The vertical blue line indicates the median  of each distribution and the red line indicates the critical value  crit for rejection of the null hypothesis that the two distributions are the same (equation C1).

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  and -filters (  ,   ) and the date of the last observation used in this work, max( ).The maximum phase angle observed in each filter is given by max(   ) and max(   ), and the maximum minus minimum phase angle range for each filter is Δ  , Δ  .Finally, the orbital elements obtained from astorb are also included; semimajor axis , eccentricity  and inclination .

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

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  app,

Table 4 .
The final  12  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.
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 (figureB1), we generated a range of ideal model  phase curves with 0 ≤  ≤ 1 and then fitted those simulated observations with the  12  12 model.The locus of best fit  12 values for phase curves with parameter  is shown in the top panel of figure8.It can be seen that fitting P16 to observed phase curves naturally gives rise to a range of  12 larger than (0, 1) .