Phase Curves of Kuiper Belt Objects, Centaurs, and Jupiter-family Comets from the ATLAS Survey

All ATLAS apparent magnitude measurements were transformed to reduced magnitudes to remove distance effects. The reduced magnitude, M, is defined as

$$\begin{eqnarray}&&M=m-5{\mathrm{log}}_{10}\left(r{\rm{\Delta }}\right),\end{eqnarray} \tag{ 2 }$$

where m is the apparent magnitude, and r and Δ are the heliocentric and geocentric distances of the object at the time of each observation, respectively. Distances and solar phase angles for each observation of a given object were obtained from JPL Horizons. The uncertainty in reduced magnitude was taken to be that of the measured apparent magnitude, due to the comparatively negligible uncertainties of the geocentric and heliocentric distances. Data points with uncertainties larger than 85% of all data points were removed. These observations were likely taken during poor observing conditions, causing their large uncertainties. Each object's data set was then clipped to remove outlying points beyond 3σ from the median-reduced magnitude to exclude measurements from contaminating flux, poor image subtraction, or poor calibration. All data in a given filter were treated as a single bin of phase angle unless observations extended to phase angles >20° or <0.1°, where phase-curve profiles become nonlinear, in which cases data were split into two or four bins of phase angle, respectively. All data in each bin were sigma-clipped separately until convergence.

Beyond very small phase angles (α ∼ 0.1°; Belskaya et al. 2008) where the opposition effect is significant, KBO, Centaur, and JFC phase curves exhibit approximately linear profiles when observed from Earth. We apply a weighted linear fit to the phase curve of each object in each ATLAS filter, using the polyfit function from the NumPy Python package (van der Walt et al. 2011; Harris et al. 2020) to fit a first-degree polynomial to the reduced magnitudes, M(α), according to

$$\begin{eqnarray}&&M(\alpha )=H+\alpha \beta ,\end{eqnarray} \tag{ 3 }$$

where α is the solar phase angle of the object at the time of measurement, β is the linear phase coefficient, and H is the absolute magnitude of the object. We consider only data at phase angles α ≥ 0.1° to remove any influence of an opposition surge on the resulting β value. Each data point was weighted by the inverse of its magnitude uncertainty. The best-fit values of both β and H and their associated uncertainties were computed via Monte Carlo simulation. We generate 10⁵ synthetic phase curves with synthetic reduced-magnitude measurements, M_i,synthetic, calculated according to

$$\begin{eqnarray}&&{M}_{i,\mathrm{synthetic}}={M}_{i}+\mathrm{rand}[0,\delta {M}_{i}],\end{eqnarray} \tag{ 4 }$$

where rand[0, δ M_i ] is a random number generated from a normal distribution centered at zero of standard deviation equal to δ M_i , the uncertainty of each reduced-magnitude measurement. The best-fit β and H values were taken to be the medians of the distributions of all 10⁵ results for each parameter. These distributions tended to be approximately Gaussian in profile, thus we took the associated uncertainties of these parameters to be half of the range that includes the central 95% of their corresponding distributions, yielding 2σ uncertainty values. Figure 4 shows an example for the KBO Quaoar of an ATLAS phase curve of our sample and the associated uncertainty distributions. For objects for which changes in brightness due to their periodic light curve or outbursts of cometary activity can be identified (applicable to five objects in total), we correct their data sets for these effects (see Sections 5.3 and 5.4 for light curve and cometary activity correction, respectively) and reapply the phase-curve fitting algorithm to the corrected data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (18 panels)

Tarball of figure set images

We aim to correct for the effect of rotational modulation on the phase curves wherever possible by searching for periodicities corresponding to the period of the rotational light curve of an object for all ATLAS data in our sample. For most objects in our sample the uncertainty in the reduced-magnitude measurements exceeds the literature peak-to-peak amplitude of the rotational light curve. This obscures any rotational modulation, thus precluding any correction of this effect on the phase curve of the object. However, for some objects in our sample we observe substructure in the residuals to the best-fitting linear phase-curve function, which extends beyond the associated uncertainties of the brightness measurements, potentially due to rotational modulation.

To search for periodicities in the light curve, we use the Lomb–Scargle periodogram algorithm (Lomb 1976; Scargle 1982) due to the irregular time spacing of the ATLAS observations. For each filter the best-fit linear phase curve is subtracted from each data set. The resulting data, corrected also for light travel time, were then passed into the Lomb–Scargle periodogram function (LombScargle) from the astropy Python package (Astropy Collaboration et al. 2013, 2018), weighted by the reduced-magnitude uncertainty. Frequency grids were chosen with minimum and maximum rotation frequencies of 0.1 day⁻¹ (240 hr period) and 24 day⁻¹ (1 hr period), with 10⁷ periodogram evaluations between these frequency extremes. Potential rotation periods for each object in the sample were searched for in the c and o data sets separately. For a given filter the maximum peak in the power spectrum was taken to be the object's best-fitting period value, with its associated uncertainty equal to half the difference between the period values corresponding to the closest two stationary points either side of the maximum peak.

We highlight that the best-fit period value is subject to uncertainty due to aliasing, background noise, and uncertainties in the measured magnitudes (VanderPlas 2018). A method of quantifying the validity of a given rotation period is to calculate the false-alarm probability of the maximum peak of the power spectrum: the probability that Gaussian noise could produce a signal of equal height to a given peak in the power spectrum (VanderPlas 2018). For the rotation period value yielded for each filter, we calculated false-alarm probabilities, choosing the Baluev method (Baluev 2008) for its suitable accuracy and computational efficiency (VanderPlas 2018). Only objects whose periods in both filters came from peaks higher than the 95% false-alarm probability level (>2σ results) were considered to be plausible.

Large photometric uncertainties in the data plus the generally fewer c-filter observations compared to the o filter could mean that significantly different best-fit period values are yielded per filter. We set a criterion that the period values from both filters must match to <0.5 hr to be considered valid. For any object which fits this criterion, we check if their rotation periods match values aliased to Earth's sidereal period (12, 24, 48, etc., hr) to <0.1 hr, to ensure we are not detecting periodicity due to sampling rate. If so, we check the rotation phase coverage of the data when phase-folded to the rotation period in each filter. If the phase-folded data of an object are grouped into clusters with gaps in between, or cover <50% of the rotation phase span, we reject that object's rotation period value.

Three objects in our sample satisfy the above criteria for valid rotation periods: Pluto, Haumea, and Hidalgo. For these objects, we hence consider only the rotation period yielded by the filter with more data points (o). The ATLAS rotation period for Pluto is consistent with values from the literature (Tholen & Tedesco 1994; Tholen & Buie 1997; Sheppard & Jewitt 2002; Schaefer et al. 2008; Buie et al. 2010; Benecchi et al. 2018) and, to within uncertainties, the ATLAS rotation periods equal the half-values of the literature rotation periods of Haumea (Rabinowitz et al. 2006; Lacerda et al. 2008; Lellouch et al. 2010; Thirouin et al. 2010; Lockwood et al. 2014) and Hidalgo (Hanuš et al. 2016; Santos-Sanz et al. 2017). We note that uncertainties in our reduced-magnitude measurements could obscure double-peaked light curves of elongated or irregular-shaped objects, causing them to appear as single peaked and the rotation period value yielded by the Lomb–Scargle algorithm to equate to half the true rotation period. Both Haumea (Rabinowitz et al. 2006; Lacerda et al. 2008; Lellouch et al. 2010; Thirouin et al. 2010; Lockwood et al. 2014) and Hidalgo (Ďurech et al. 2007; Hanuš et al. 2016) exhibit shape elongation. Our measured period values for these objects are consistent with the half the reported rotation periods, thus we consider the true rotation period of these two objects to be twice the value output from the Lomb–Scargle periodogram. The ATLAS rotation period values for Pluto, Haumea, and Hidalgo, both those yielded by the maximum periodogram peak and the final chosen rotation period values, are listed in Table 5, in addition to those from the literature.

Some rotation periods may not appear in individual filter data sets, but may be discerned when the data sets are combined to augment sampling of any rotational light curve. We repeat our Lomb–Scargle periodogram algorithm on both the ATLAS c-filter and o-filter data combined. We assume that the objects in our sample exhibit no variability in color over their rotation periods, and correct the data set in each filter for color according to the difference in the median reduced magnitudes of each filter per object. We consider only Lomb–Scargle periodogram spectra maxima higher than the 95% false-alarm probability level. For objects which satisfy this criterion, we consider only those where the data in both filters span the full rotation phase space (ensuring that the rotation period value is not being driven by one filter only) and where the standard deviation of the data points exceeds twice the median magnitude uncertainty. If the object's phase-folded light curve exhibits gaps in the phase-folded light curve, and/or the rotation period value lies close to those aliased to Earth's sidereal period (12, 24, 48 etc. hr) to <0.1 hr, we reject the object's rotation period value. We recover rotation periods for Pluto, Haumea, and Hidalgo, all consistent to those found by separate analysis of each filter. We discover a rotation period value of 5.4749 ± 0.0015 hr from the Lomb–Scargle periodogram for the JFC 2012 TW236. The Lomb–Scargle periodogram power spectra and the data sets of 2012 TW236 phase-folded to the obtained rotation periods are shown in Figure 5. Due to the small size of this object, it is likely irregularly shaped with a double-peaked rotational light curve, so we consider the rotation period of 2012 TW236 to be 10.950 ± 0.003 hr, to our knowledge the first measured value for this object.

Zoom In Zoom Out Reset image size

We phase-fold the data sets of Pluto, Haumea, Hidalgo, and 2012 TW236 in both filters to their selected rotation period value. We fit a sinusoid of period equal to the chosen rotation period value to the data using the lmfit Python package (Newville et al. 2014), weighted by the uncertainty of each magnitude. This fitted function is subtracted off the reduced-magnitude data to eliminate the rotational modulation, and remove any outlying data in the corrected data set to 3σ. The linear phase-curve model for the object is refitted again on the corrected data to measure the phase coefficient and absolute magnitude, corrected for the effect of the object's rotation on its axis. The effect of correcting for an object's rotational light curve is illustrated in Figure 6, showing the significant decrease in spread of the data about the best-fitting linear phase-curve function. The rms values of the residuals to the best-fitting linear function decrease from 0.096 to 0.066 mag in the c filter and from 0.093 to 0.057 mag in the o filter when the effects of Haumea's rotational light curve is removed. Table 4 lists the phase coefficient and absolute magnitude values for each object for which rotational modulation was corrected, before and after correction. Table 5 lists the amplitudes of the best-fitting sinusoid light-curve function in both c and o, the uncertainty being calculated from the covariance matrix of the resulting fit. We note that our amplitudes of rotational modulation differ significantly both between filters and from those in the literature, though the large uncertainties in ATLAS reduced magnitudes are a possible cause of this.

Zoom In Zoom Out Reset image size

The Centaur and JFC populations contain objects known to exhibit cometary activity (Lowry et al. 2008; Peixinho et al. 2020, and references therein). It is thus plausible that some objects in our sample may have exhibited activity, forming a coma which may partly or fully obscure an object's surface from direct observation. This would brighten the object's apparent magnitude, and could alter an object's phase curve for a period of time. Given the number of data points per object, we can use the ATLAS phase curves to search for deviations that may be indicative of possible cometary outgassing or dust production. To search for cometary outbursts, we select for observations that were significantly brighter than predicted from the measured phase curve. We assume that most of an object's data are taken from epochs of (relative) quiescence potentially punctuated by shorter epochs of activity or outbursts. We consider this assumption to be valid as most of our objects exhibit reduced magnitudes which are consistent in value across separate apparitions. We subtract the best-fitting phase-curve function from the data set (corrected for rotational light curve if possible), and select observations brighter than 2σ from the object's median brightness across the ATLAS baseline.

Active epochs of Centaurs and JFCs exhibit timescales ranging from weeks to years (Peixinho et al. 2020, and references therein). Such activity would manifest as a series of bright observations across consecutive nights. To detect such cometary activity, we utilize the ability of ATLAS to observe a given target multiple times per night, and select nights containing at least three observations, of which ≥75% have magnitudes brighter than 2σ from the object's median brightness. We term such nights of observations "bright nights" henceforth. We further analyze these bright nights to determine if cometary activity is the source of the brightening. We note that objects may exhibit activity which does not produce short increases in brightness, and we would not be able to detect this activity using the described algorithm. To detect any such lower-level cometary activity, we consider only observations whose phase-corrected reduced magnitudes lie <2σ from the object's sigma-clipped median brightness. We consider each data point in turn, and bin all data taken 30 days after each observation. If at least 14 observations (the approximate average number of ATLAS exposures per week) are taken in any given interval, we consider the median brightness of all data in that interval. If the binned median brightness exceeds the object's sigma-clipped median brightness by 0.2 mag and 1 standard deviation of the binned data, we consider all data in this interval to be a long outburst. We remove any data which correspond to epochs of cometary activity from the data sets and reapply the phase-curve fitting algorithm to the remaining data. If any instance of cometary activity is found to extend over time, we remove all data in both filters that falls within the time range of the activity.

Several objects in our sample were observed across phase angle ranges comparable to those of main belt asteroids. This allowed us to fit their data from ATLAS with phase-curve functions used for these asteroids. The resulting best-fit parameters for our sample could then be compared to those of the asteroid population. We utilize the three-parameter HG₁ G₂ phase-curve model developed by Muinonen et al. (2010), where H is the absolute magnitude of the phase curve, and G₁ and G₂ both quantify the slopes of the phase angle functions that govern the overall phase-curve shape, analogous to the linear phase coefficient β. We utilize this function as per the findings of Mahlke et al. (2021) that taxonomic complexes could be distinguished using the HG₁ G₂ function applied to ATLAS data. Mahlke et al. (2021) also found the related ${{HG}}_{12}^{* }$ asteroid phase-curve function (Penttilä et al. 2016) failed to distinguish between complexes for ATLAS data, thus we do not apply the ${{HG}}_{12}^{* }$ function to our data. We select only phase curves in a given filter to fit with the HG₁ G₂ model which satisfy the following criteria:

• 1.  
  Observations of α ≥ 12°, to ensure we are sampling the section of the phase curve where curvature becomes significant and the accuracy of linear fits is reduced.
• 2.  
  Phase curves must contain ≥35 data points at phase angles 0 ≤ α ≤ 5°, to ensure we have sufficient sampling of the region where the opposition surge can become significant, which may affect the profile of the fitted phase-curve model.

The o-filter phase curves of Narcissus, Hidalgo, and 2014 QA43, and the c-filter phase curve of Hidalgo satisfy these criteria. We apply the HG₁ G₂ model to the reduced-magnitude values of these data sets using the HG1G2 function from the Python photometry module sbpy.photometry (Mommert et al. 2019). We use a Monte Carlo algorithm to measure the absolute magnitude and G₁ and G₂ parameters, generating 10⁵ synthetic phase curves by varying each reduced-magnitude value by a random fraction of its associated uncertainty (assuming a normal uncertainty distribution) and fitting the HG₁ G₂ model to each synthetic data set. The final values for H, G₁, and G₂ are the median of the resulting distribution of all 10⁵ values of each parameter, with their associated uncertainties being half of the range between the values which enclose 95% of all values in the distribution, yielding 2σ uncertainties. If an object's data satisfy the requirements for fitting with the HG₁ G₂, we retroactively use the resulting fit instead of the linear model for eliminating phase angle dependence when searching for rotation periods and cometary activity in the full data set. For the object Narcissus, due to near-identical resulting chi-square values, we use both the linear and HG₁ G₂ fits when searching for rotation period and cometary activity, finding that the choice of fitted function has no effect on the results yielded for Narcissus by either algorithm.

We measured the linear phase coefficient β and absolute magnitude H in both ATLAS c and o filters for all 18 objects in our sample across the full range of phase angles spanned by each data set (excluding opposition surges at α < 0.1°). Table 6 lists our measured β and H values and 2σ errors for all objects in our sample across the full phase angle range. Figures for all c and o ATLAS phase curves and the two-dimensional histograms of the fitted phase-curve parameters are shown in Figure 4 and its online figure set. Linear fits applied to the full range of phase angles spanned by the data were used to correct for variability due to object rotation and cometary activity for all objects, apart from the o-filter data sets of Narcissus and 2014 QA43 and the c- and o-filter data sets of Hidalgo, which were corrected using the HG₁ G₂ function. Four such objects—2012 TW236, 2014 QA43, 2016 ND21, and Narcissus—have data sets which extend to phase angles beyond α ∼ 10° where nonlinearity could become significant, but the small number of objects this applies to means we consider a linear fit suitable for variability correction. However, for the remainder of our analysis, including a comparison to the Kuiper Belt and is associated populations, we restrict ourselves to measuring phase coefficients and absolute magnitudes for a fixed phase angle range of 0.1° ≤ α ≤ 6°, as detailed in Section 6.2.

The objects in our sample span semimajor axes of 5.7 au ≤ a ≤ 68.1 au, resulting in a wide range of maximum phase angles observable from Earth. ATLAS observations thus sample different sections of the phase curves depending on an object's heliocentric distance. To compare phase curves over a fixed range of phase angle, we restrict our data sets (corrected for rotational modulation and cometary activity where possible and excluding data at α ≤ 0.1°) to include only observations taken at phase angles α ≤ 6°, and consider only objects with >50 data points in this range in both c and o. We apply our methods to the resulting phase curves of 13 objects. The resulting absolute magnitudes and phase coefficients in c and o are listed in Table 7. Comparing the resulting values to those measured over the maximum phase angle range, we find that the phase coefficient and absolute magnitude values change beyond 2σ error bars only for the Transition Objects Echeclus and 2009 YF7 in the c filter. Figures for all c and o ATLAS phase curves for data in the phase angle range 0.1 ≤ α ≤ 6 are shown in Figure 7 and its online figure set.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (13 panels)

Tarball of figure set images

We note that our β values differ by >2σ between the ATLAS filters for many objects. Significant variation of β with wavelength has been previously reported by several multifilter phase-curve studies of Centaurs and KBOs (Buratti et al. 2003; Rabinowitz et al. 2007; Schaefer et al. 2009; Alvarez-Candal et al. 2016; Ayala-Loera et al. 2018; Alvarez-Candal et al. 2019). Figure 8 shows the ATLAS β values compared to those from the literature in the Johnson–Cousins BVRI photometric system, for which the ATLAS filters c and o exhibit significant overlap in wavelength range with B and V filters (c filter) and R and I filters (o filter), respectively. Though the ATLAS β values often differ in value from the literature, for most objects they overlap within the 2σ error bars or lie between the extremes of the range of the literature values for a given object. Our ATLAS linear phase coefficients are consistent with those from the literature, as expected, and we attribute any disparity to be due to the different wavelength ranges of the filters of the ATLAS and Johnson–Cousins photometric systems. We also note that our measured phase coefficient and absolute magnitudes for Pluto are those of the Pluto-Charon system, the components of which cannot be resolved at the ATLAS pixel scale.

Zoom In Zoom Out Reset image size

Figure 9 shows the best-fitting linear phase-curve models for all objects in our sample, divided by dynamical class and photometric filter. The β values indicate relatively flat phase curves across our entire sample, with values ranging from −0.01 to 0.18 mag/deg in c, and 0.02 to 0.106 mag/deg in o. Such flat phase curves have been previously observed for several JFCs, Centaurs, and Pluto-sized KBOs (Rabinowitz et al. 2006, 2007; Kokotanekova et al. 2017, 2018; Verbiscer et al. 2022) and are thought to be indicative of low-albedo, organic-rich Centaur and JFC surfaces (Nelson et al. 2000; Shkuratov et al. 2002; Cruikshank & Dalle Ore 2003; Rabinowitz et al. 2007) and high-albedo KBO surfaces covered in granular volatile ices (Brown et al. 2005; Licandro et al. 2006; Rabinowitz et al. 2006, 2007; Trujillo et al. 2007). Furthermore, we find that for both the c and o filters, the phase-coefficient values of the Centaurs and Transition Objects all fall in the same range as those of the KBO population when measured across the same phase angle range 0.1° ≤ α ≤ 6°. All but one of our samples exhibit positive phase-coefficient slopes in both ATLAS filters. Huya exhibits a marginally negative phase coefficient in the c filter, as shown in Figure 10, but is consistent within 2σ of having a flat or positive slope. Alvarez-Candal et al. (2016), Ayala-Loera et al. (2018), and Alvarez-Candal et al. (2019) reported negative phase coefficients for several objects in their sample, including Chiron, Eris, Haumea, Huya, and Makemake. β < 0 would imply an unphysical dimming near opposition. Alvarez-Candal et al. (2016) and Ayala-Loera et al. (2018) speculated this effect may be due to unidentified ring systems, epochs of cometary activity, or poorly corrected rotational modulation. We instead measure phase-coefficient values for these objects which are either positive slopes or consistent with positive slopes to within 2σ. Furthermore, Huya's negative c-filter phase coefficient is derived from the second sparsest data set in our sample of N_c = 68 observations, comparable in size to the V-filter data set from Ayala-Loera et al. (2018; N_V = 45), which also yielded a negative phase-coefficient value. This tendency toward positive β values with increasing data set size suggests negative phase-coefficient values are most likely due to small data sets with large scatter that are unable to sufficiently sample the phase-curve profile, and not a real property of these objects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We explore potential correlations between phase coefficients β_c and β_o , absolute magnitudes H_c and H_o , color index H_c –H_o , and relative phase coefficients β_c –β_o derived from the ATLAS phase curves. We analyze parameter values measured across the phase angle range 0.1 ≤ α ≤ 6 deg to ensure all parameters are measured over the same region of the phase curve. We make use of the Spearman rank correlation test to search for monotonic relationships between pairs of variables, using the spearmanr function from the Python module scipy (Virtanen et al. 2020). The spearmanr function calculates two parameters. The first is the Spearman rank correlation coefficient, r_s , which quantifies the strength and direction of any detected monotonic relationship between two variables; ∣r_s ∣ → 1 defines a strong correlation, whereas r_s ≈ 0 implies no correlation. The second parameter is the null hypothesis probability value of that correlation, or p-value, ${P}_{{r}_{s}}$, which quantifies the probability that any detected correlation is derived from an inherently uncorrelated distribution of data. We note, however, that the Spearman rank correlation test itself does not consider uncertainties in the analyzed variables. To account for this, for every parameter pair tested we perform a Monte Carlo simulation on the parameter values. For each data point we vary the nominal value within associated uncertainties, assuming every uncertainty follows a normal distribution with a standard deviation equal to the nominal value of the uncertainty itself, perform the rank statistic on the resulting synthetic values, and repeat 10⁵ times. We accumulate all the associated p-values and calculate the percentage of simulations for which ${P}_{{r}_{s}}\gt 0.05$, i.e., we could not rule out an uncorrelated population to 95% probability. We consider a correlation to exist if ∣r_s ∣ ≥ 0.5 and <5% of simulations return ${P}_{{r}_{s}}\gt 0.05$. A correlation that satisfies one of these criteria is considered tentative, and we reject any correlation if neither is satisfied or ${P}_{{r}_{s}}\gt 95 \%$. We apply the Spearman rank correlation test to our full sample and subsamples: Centaurs/JFCs/Transition Objects; KBOs; dwarf planets (H ≤ 3); and non-dwarf planet objects from the KBO, Centaur, JFC, and Transition Object populations (H > 3). Figure 11 shows plots of the tested ATLAS parameters for the full object sample, the Centaurs, JFCs, Transition Objects, the KBOs, the dwarf planets, and the non-dwarf planet objects in our sample.

Zoom In Zoom Out Reset image size

Table 8 lists the parameters from our Spearman rank correlation tests. We find no strong correlations between β_c and H, β_o and H, and β_o and H_c –H_o for the full sample and all subpopulations. However, we find a tentative negative correlation between β_c and H_c –H_o for the full sample and all subpopulations as we cannot reject the null hypothesis to >95% confidence. A tentative negative correlation between β_o and H_c –H_o is also found for the Centaur/JFC/Transition Object and non-dwarf planet populations. We also find a tentative negative correlation between β_c –β_o and H_c –H_o for the full sample, the KBO population, and the non-dwarf planet population, yet no significant correlation between these parameters for the dwarf planet population, and a tentative correlation that can almost be ruled out to 95% confidence for the Centaur/JFC/Transition Objects subpopulation. No strong correlation is identified for any tested parameter pair. We also find no significant change in correlation results when applying a correlation search algorithm to parameter values derived from the full phase angle range.

Table 8. Results of Correlation Investigation Using ATLAS Data in Phase Angle Range 0.1° ≤ α ≤ 6°

Parameters	rs	${P}_{{r}_{s}}$	% Simulations with ${P}_{{r}_{s}}\gt 0.05$	N	Correlation?
Full Sample
βc versus Hc	0.214	0.482	99.985	13	No
βo versus Ho	0.159	0.603	99.99	13	No
βc versus Hc –Ho	−0.747	0.003	13.931	13	Tentative
βo versus Hc –Ho	−0.297	0.325	95.694	13	No
βc –βo versus Hc –Ho	−0.648	0.017	52.726	13	Tentative
Centaurs/JFCs/Transition Objects
βc versus Hc	0.257	0.623	100.0	6	No
βo versus Ho	−0.314	0.544	99.792	6	No
βc versus Hc –Ho	−0.771	0.072	63.171	6	Tentative
βo versus Hc –Ho	−0.657	0.156	73.967	6	Tentative
βc –βo versus Hc –Ho	−0.771	0.072	94.894	6	Tentative
KBOs
βc versus Hc	−0.286	0.535	99.94	7	No
βo versus Ho	0.643	0.119	99.538	7	No
βc versus Hc –Ho	−0.75	0.052	67.321	7	Tentative
βo versus Hc –Ho	−0.036	0.939	99.975	7	No
βc –βo versus Hc –Ho	−0.786	0.036	72.76	7	Tentative
Dwarf Planets
βc versus Hc	0.143	0.787	99.3	6	No
βo versus Ho	0.771	0.072	98.216	6	No
βc versus Hc –Ho	−0.714	0.111	68.087	6	Tentative
βo versus Hc –Ho	−0.086	0.872	99.969	6	No
βc –βo versus Hc –Ho	−0.714	0.111	97.497	6	No
Non-dwarf Planets
βc versus Hc	0.536	0.215	100.0	7	No
βo versus Ho	−0.071	0.879	97.804	7	No
βc versus Hc –Ho	−0.857	0.014	43.9	7	Tentative
βo versus Hc –Ho	−0.679	0.094	79.644	7	Tentative
βc –βo versus Hc –Ho	−0.857	0.014	82.876	7	Tentative

Note. Spearman rank correlation coefficients r_s , associated null hypothesis probability values ${P}_{{r}_{s}}$, percentage of Monte Carlo simulations where ${P}_{{r}_{s}}\gt 0.05$ (i.e., we cannot reject a random distribution to 95% confidence), and number of data points analyzed for correlation tests between phase-curve parameters using ATLAS data.

Download table as:  ASCIITypeset image

We also apply our algorithm to the data of Ayala-Loera et al. (2018), who reported phase coefficients and absolute magnitudes for 114 objects in Johnson–Cousins V and R filters, which exhibit overlap in wavelength range with the ATLAS c and o filters. We select only objects from their study with β values consistent with positive to within 2σ and divide their object sample according to our definitions. Table 9 lists the results from applying our algorithm to their data. For β and H, we find tentative negative correlations for the Centaurs/JFCs/Transition Objects subpopulation in both the V and R filters. For β and H_V −H_R , we find no correlation in V for any subpopulation, yet find a tentative positive correlation for the Centaurs/JFCs/Transition Objects and the dwarf planets subpopulation. For β_V −β_R and H_V −H_R , we find highly significant negative correlations for the full sample, the KBOs, and the non-dwarf planet subpopulations.

Table 9. Results of Correlation Investigation Using Data from Ayala-Loera et al. (2018)

Parameters	rs	${P}_{{r}_{s}}$	% Simulations with ${P}_{{r}_{s}}\gt 0.05$	N	Correlation?
Full Sample
βV versus HV	−0.278	0.005	18.208	100	No
βR versus HR	−0.176	0.081	77.217	99	No
βV versus HV −HR	−0.219	0.037	81.36	91	No
βR versus HV −HR	0.408	0.0	10.796	91	No
βV − βR versus HV −HR	−0.747	0.0	0.133	91	Yes
Centaurs/JFCs/Transition Objects
βV versus HV	−0.508	0.026	14.087	19	Tentative
βR versus HR	−0.521	0.022	13.153	19	Tentative
βV versus HV −HR	−0.449	0.054	90.363	19	No
βR versus HV −HR	−0.573	0.010	75.758	19	Tentative
βV −βR versus HV −HR	−0.307	0.201	84.616	19	No
KBOs
βV versus HV	−0.152	0.177	99.297	81	No
βR versus HR	−0.110	0.330	99.959	80	No
βV versus HV −HR	−0.271	0.022	75.219	72	No
βR versus HV −HR	0.418	0.0	19.681	72	No
βV −βR versus HV −HR	−0.802	0.0	0.248	72	Yes
Dwarf Planets
βV versus HV	0.143	0.787	100.0	6	No
βR versus HR	0.6	0.208	99.623	6	No
βV versus HV −HR	0.429	0.397	86.778	6	No
βR versus HV −HR	0.771	0.072	86.046	6	Tentative
βV −βR versus HV −HR	−0.714	0.111	96.087	6	No
Non-Dwarf Planets
βV versus HV	−0.241	0.019	56.606	94	No
βR versus HR	−0.225	0.031	35.844	93	No
βV versus HV −HR	−0.230	0.034	81.014	85	No
βR versus HV −HR	0.377	0.0	19.221	85	No
βV −βR versus HV −HR	−0.745	0.0	0.28	85	Yes

Note. Spearman rank correlation coefficients r_s , associated null hypothesis probability values ${P}_{{r}_{s}}$, percentage of Monte Carlo simulations where ${P}_{{r}_{s}}\gt 0.05$ (i.e., we cannot reject a random distribution to 95% confidence), number of data points N analyzed for correlation tests between phase-curve parameters from Ayala-Loera et al. (2018), and if the parameters exhibit a correlation.

Download table as:  ASCIITypeset image

We find a highly significant negative correlation between relative phase coefficient and color index using data from Ayala-Loera et al. (2018) for objects with phase-coefficient values consistent with positive to within 2σ. We also find similar negative correlations with ATLAS data for the full sample, the KBOs, and the non-dwarf planets subpopulations. However, these correlations are tentative, and an inherently uncorrelated distribution cannot be ruled out to 95% confidence. Jackknife resampling of the Ayala-Loera et al. (2018) sample, using subsample sizes of N = 13 (for the full ATLAS sample) and N = 7 (for both the KBO and non-dwarf planet subpopulations), yields no correlation for 16.6% and 53.4% of simulations, respectively. Thus, the difference in size between our sample and that of Ayala-Loera et al. (2018) is likely responsible for this discrepancy in correlation certainty. Furthermore, our ATLAS data sets usually have many more observations than those of Ayala-Loera et al. (2018), and are taken in filters of different wavelength ranges. We therefore consider our findings to be consistent with those of Ayala-Loera et al. (2018) and attribute any difference in correlation certainty to be due to differing sample sizes, data set sizes, and filter wavelength ranges. An increase in the number of objects with phase coefficient and absolute magnitude measurements from data sets comparable to those afforded by ATLAS would allow us to potentially recover this negative correlation between relative phase coefficient—the difference in phase-coefficient value between broadband filters of differing wavelength range—and color index.

Ayala-Loera et al. (2018) reported the strong negative correlation between relative phase coefficient and color index as persisting even when their sample was divided by absolute magnitude and semimajor axis. However, their sample was divided by semimajor axis and absolute magnitude at a = 40 au and H = 4.5, in contrast to our thresholds of a = 30 au and H = 3. Dividing our sample by the subpopulation definitions of Ayala-Loera et al. (2018) using H values calculated from ATLAS data, we find tentative correlations for the a < 40 au and H > 4.5 subpopulations, and no correlation for the H < 4.5 au subpopulation, with the a > 40 au population being too small for an accurate rank statistic. Applying our algorithm to their sample divided by their definitions, we also find strong negative correlations between relative phase coefficient and color index for the H > 4.5 population and a > 40 au population, and tentative negative correlations for the a < 40 au population and H < 4.5 au populations. Dividing their sample by our definitions of subpopulations, we find no correlation for Centaurs/JFCs/Transition Objects or the dwarf planets subpopulations, yet a strong correlation for the KBOs and the non-dwarf planet subpopulation, which is dominated in number by KBOs. This difference in observed correlation between populations may point to a difference in substructure between the surfaces of these objects, likely caused by an evolutionary effect due to differing insolation. An increased sample size similar to those of Alvarez-Candal et al. (2016), Ayala-Loera et al. (2018), and Alvarez-Candal et al. (2019), with precisions of parameter values comparable to those afforded by the ATLAS data sets, would help to better recover and study this negative correlation between relative phase coefficient and color index, which may in fact be a universal property across different populations of small solar system objects (Alvarez-Candal et al. 2022), and confirm the existence of this possible distinction.

We searched for instances of cometary activity in the ATLAS data sets for all 18 objects in our sample. Most of our samples show no sign of cometary activity. A small number of objects, namely Chiron, Echeclus, Pluto, and Quaoar, exhibit "bright nights" as defined in Section 5.4, suggesting potential cometary outbursts or long-term activity. Table 10 lists the bright observations per object identified by our algorithm as potentially indicative of cometary activity, either in the form of short outbursts or longer-term sustained brightening. We discuss each object below.

Echeclus. We detect three separate instances of bright nights across the ATLAS baseline, as shown in Figure 12. The first, seen in the c filter on 2018 January 17 (MJD 58135) and the o filter from 2017 December 24 (MJD 58111) to 2018 January 2 (MJD 58120), coincides in time with a previously reported 2017 December outburst (James 2018; Kareta et al. 2019). During this outburst, Echeclus exhibited a fan-shaped coma that remained clearly visible in the ATLAS images for approximately 1 month, as seen in Figure 13. This activity seemed to build to a second outburst observed as two consecutive bright nights, 2018 March 6 (MJD 58183) and 2018 March 8 (MJD 58185). Analysis of the FWHM of the PSF of Echeclus during the night of 2018 March 6 (the night during which the PSF of Echeclus was more clearly visible) reveals significant spread compared to those of background stars, as shown in Figure 14. The PSF FWHM of Echeclus compared to that of the background stars for this date are listed in Table 11, indicative of a visible coma. The third instance appears as a single bright night in the o filter on 2018 April 9; however, analysis of the ATLAS images reveals Echeclus to lie close to a poorly subtracted background star during this time, as seen in Figure 15, making this a result of background stellar flux contamination and not a comet-like outburst.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Chiron. We detect two separate instances of bright nights, highlighted in Figure 16. The first is a single bright night at 2016 December 22 (MJD 57744) observed in c; visual inspection reveals Chiron lies close to a poorly subtracted background star, as seen in Figure 15, thus this brightening is due to stellar contamination. The second instance is an increase in brightness constituting three consecutive bright nights starting at 2021 June 18 (MJD 59383), the first three nights of observation of a new apparition of Chiron. No coma or tail is visible in the ATLAS images, and Chiron's flux signal appears stellar (see Figures 17 and 18 and Table 11), yet this brightness increase exceeds Chiron's amplitude of rotational modulation. This new epoch of increased activity is reported in Dobson et al. (2021a) and Dobson et al. (2021b). Further analysis of Chiron is beyond the scope of this study, and will be the subject of a later paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Pluto. We detect two bright nights at 2019 July 5 (MJD 58669) and 2020 October 7 (MJD 59129), highlighted in Figure 19, yet neither are caused by cometary activity. Visual inspection of ATLAS images of each night shows that Pluto is contaminated by background stars, as seen in Figure 15.

Zoom In Zoom Out Reset image size

Quaoar. We detect a single bright night at 2021 March 22 (MJD 59295), highlighted in Figure 20. We find Quaoar lies close to a poorly subtracted background star in all observations from this night, as seen in Figure 15, contaminating its photometry.

Zoom In Zoom Out Reset image size

We also detect two further instances of brightening for both Echeclus and 2014 QA43 which may correspond to long-term activity, described below.

Echeclus. We detect an instance of sustained brightening of 0.2 mag and >1 standard deviation of the data from the sigma-clipped median ranging from 2017 November 10 (MJD 58067) to 2018 February 12 (MJD 58161), as highlighted in Figure 12. This coincides in time with Echeclus's >2σ outburst in 2017 December, and a fan-shaped coma is clearly visible in the ATLAS images during this time, as seen in Figure 13.

2014 QA43. We detect an instance of brightening ranging in time from 2020 September 16 (MJD 59108) to 2020 December 7 (MJD 59190), spanning a large fraction of the observations from that apparition, as shown in Figure 21. The median phase-corrected reduced magnitude of this apparition is brighter than that of previous years of data by approximately 0.3 mag. Measurements of the FWHM of 2014 QA43 reveal no significant broadening compared to background stars, as seen in Figure 22. It remains unclear if this brightening is due to cometary activity not visible in the images or due to the changing aspect angle of the object.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In conclusion, only Chiron and Echeclus are found to exhibit definite comet-like activity. It is not surprising that we find these objects have been active during the ATLAS survey, as both Chiron (Bus et al. 1988; Luu & Jewitt 1990; Dahlgren et al. 1991; Marcialis & Buratti 1993; Elliot et al. 1995; Lazzaro et al. 1997; Silva & Cellone 2001; Duffard et al. 2002) and Echeclus (Choi et al. 2006a, 2006b; Jaeger et al. 2011; Rousselot et al. 2016; James 2018; Seccull et al. 2019; Kareta et al. 2019) have well-recorded histories of comet-like activity in the literature. Both these objects have been classified by previous studies as belonging to the Centaur population, members of which have been known to exhibit comet-like outbursts. However, unlike comets, Centaur activity is decoupled from perihelion. The mechanism responsible for cometary activity on Centaurs currently remains elusive at present. Jewitt (2009) found that from a population of 23 Centaurs, the median perihelion of active Centaurs was smaller than that of inactive ones, and proposed a transition from amorphous to crystalline water ice to be the cause of Centaur activity. Considering only objects in our sample which satisfy the Centaur definition of Jewitt (2009), we find that Echeclus and Chiron have perihelia between those of Jupiter and Saturn, with those of our inactive Centaurs extending to beyond that of Saturn. Though our sample is too small to accurately quantify median perihelia, our observations of Centaur activity appear to be consistent with the findings of Jewitt (2009).

All of the variable (or comet-like activity) observations in the ATLAS data set that correspond to instances of potential cometary activity were removed from the data set before phase-curve analyses were conducted. If any found instance of possible cometary activity was extended in time, we remove all data points in both filters between the start date and end date of the activity. Furthermore, any data identified by our algorithm as a bright night yet showing background star contamination were removed from the data set. We note that the number of observations taken during an outburst remains a limiting factor in detecting such activity. For example, despite occurring during the ATLAS observation baseline, a previously recorded outburst of Echeclus in 2016 August (Rousselot et al. 2016) did not cause a >2σ increase in the ATLAS data, and observations were too sparse for it to be detected as lower-level cometary activity. Furthermore, any lower-level cometary activity that persisted for several consecutive apparitions would likely remain undetected, thereby reducing the accuracy of any applied phase-curve fit. However, our algorithm has allowed us to detect an instance of cometary activity on Chiron despite displaying no visible coma or comet-like tail. Such activity would remain undetected from analysis of ATLAS postage stamps, yet is readily detectable in our algorithm. Therefore, we present our phase-curve fits as the best analysis that can be carried out within the limits of the data set.

We apply the HG₁ G₂ phase-curve model to the c-filter data set of Hidalgo and the o-filter phase curve of Hidalgo, Narcissus, and 2014 QA43. The resulting fits to the data are shown in Figure 23 and its online figure set. Table 12 lists the absolute magnitudes H and slope parameters G₁ and G₂ of these phase curves. Our best-fit slope parameters are consistent within their uncertainties with monotonic functions of positive slope, as predicted from surface reflectance models (Muinonen et al. 2010; Mahlke et al. 2021). Mahlke et al. (2021) utilized serendipitous ATLAS observations of main belt asteroids to fit their phase curves with the HG₁ G₂ function to determine possible correlations of slope parameters with main belt asteroid (MBA) taxonomic class. We aim to compare our best-fit slope parameters with those from the taxonomic classes to infer the surface types of Hidalgo, Narcissus, and 2014 QA43. Figure 24 displays the kernel density estimation (KDE) plots of G₁ and G₂ values of 19,708 MBAs as measured by Mahlke et al. (2021) sorted by taxonomic complex, with our derived slope parameters for our three JFCs overplotted with their associated error bars. We include the c-filter slope parameter of Hidalgo for completeness, but for analysis we consider only the slope parameters for our objects in the o filter as they were derived from more observations. As seen in Figure 24, the o-filter slope parameters of all three fitted JFCs are consistent within uncertainties with the KDE distributions of all asteroid taxonomic complexes to 3σ, with greater correspondences to the carbonaceous C-type complexes. The surfaces of Narcissus and 2014 QA43 remain unclassified, while Hidalgo has been previously classified as a D-type asteroid (Tholen 1984; Belskaya et al. 2017; Geem et al. 2022). Our G₁ and G₂ parameters for Hidalgo are consistent with its previous classification, lying within the 95% confidence interval, albeit just outside the core of the distribution of the D-type complex.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (4 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size