The multi-wavelength phase curves of small bodies: Phase coloring

Context. Small bodies change their brightness due to different motives: Rotation along their axis or axes, combined with irregular shapes and/or changing surface properties, or changes in the geometry of observations. In this work, we tackle the problem of Phase curves, which show the change in brightness due to changes in the fraction of illuminated surface as seen by the observer. Aims. We aim to study the effect of the phase curves in the five wavelengths of the Sloan Digital Sky Survey in scores of objects (several tens of thousands), focusing particularly on the spectral slopes and the colors and their changes with phase angle. Methods. We used a Bayesian inference method and Monte Carlo techniques to retrieve the absolute magnitudes in five wavelengths, using the results to study the phase coloring effect in different bins of the semi-major axis. Results. We obtained absolute magnitudes in the five filters for over 40 000 objects. Although some outliers are identified, most of the usual color-color space is recovered by the data presented. We also detect a dual behavior in the spectral slopes, with a change at ${\alpha\approx}$ 5 deg.


Introduction
Large photometric surveys produce (and will produce) massive amounts of data, challenging us to develop new and better methods to make the most of their data products.Among the many surveys ongoing and to come, we would like to call attention to two: The Sloan Digital Sky Survey (SDSS, York et al. 2000) and the Legacy Survey of Space and Time (LSST, Ivezić et al. 2019) to be carried out by the Vera C. Rubin Observatory.The first one has been very productive for the small bodies community, while the second one will probably challenge all we know about small bodies due to the large number of objects to be observed and the depth of the survey (r ′ = 24.7 for a single visit, see Jones et al. 2009;Ivezić et al. 2019).Both surveys use similar photometric systems, the ugriz system (Fukugita et al. 1996) increasing the synergy between both datasets and with other surveys using similar photometric systems, for instance SkyMapper (Wolf et al. 2018), or different wavelength coverage, like the near-infrared (Popescu et al. 2016;Carry 2018).
Among the many science cases explored with the SDSS data, we will focus on small bodies' phase curves (PC) to study the phase coloring effect.A PC relates the change of apparent brightness, normalized to unit distances from the Sun and the observer, of a small body with the phase angle, α, and provides information regarding the micro and macro properties of the surfaces, as well as its absolute magnitude (H) (e.g., Nelson et al. 1998;Shkuratov et al. 2002;Grynko & Shkuratov 2008).The absolute magnitude in the V filter of a small body is defined as the object's apparent magnitude if illuminated by the Sun at 1 AU while being observed from a distance of 1 AU at opposition (i.e., α = 0 deg).H V is related to the size (D) and albedo (p V ) of the object through the well-known relation (Bowell & Lumme 1979): Taking advantage of the serendipitous observations of small bodies in the SDSS data at different epochs and orbital configurations, Alvarez-Candal et al. (2022a) developed a method to analyze the apparent magnitudes, together with the ephemeris at the time of observations, in particular topo-and heliocentric distances, ∆ and R, both in astronomical units (AU) and phase angle (α).The method produces multi-wavelength PCs and absolute magnitudes in all the filter system (H λ ).In this work, we further improve the algorithm and apply it to a different set of measurements derived from the SDSS survey, as described below.
The work we present here goes in line with a series of works analyzing massive databases, using different methods: Oszkiewicz et al. (2011) made a comprehensive study of the Minor Planets Center database, Vereš et al. (2015) analyzed the Pan-STARRS dataset producing about a quarter of a million absolute magnitudes.More recently, Mahlke et al. (2021) used ATLAS data to produce two-wavelength PCs, while Colazo et al. (2021) used GAIA (DR2) data to obtain H G for over 10 000 objects.They all intend to work with large databases, but they have yet to attempt a consistent analysis of the phase coloring along the Solar System.
Phase reddening, which we have been calling phase coloring until now, has been recognized as an important effect that may change the colors or the spectral slope of the objects (for example Taylor et al. 1971;Millis et al. 1976;Luu & Jewitt 1990).It has also been shown that it may impact the taxonomic classification of objects in the border between taxa (Sanchez et al. 2012).A clear understanding of the phase reddening effect is crucial for studying the space weathering acting on the small bodies' surfaces because it has a similar effect, and disentangling both is not straightforward (Sanchez et al. 2012).Finally, in this paper, we argue that the term phase reddening should be changed to the more descriptive phase coloring, terminology that will be used in the remaining text.
This paper is structured as follows: Description of the data in Sect.2, while the improvements made to the algorithm are described in Sect.3. The results are presented and discussed in Sect. 4 to draw the Conclusions in the last section.

Dataset
We use the data of small bodies extracted from the SDSS, but not the traditional and well-known Moving Objects Catalog (Ivezić et al. 2001;Jurić et al. 2002) frozen in its fourth release, with data obtained until March 2007.Instead, we use the re-analysis of the SDSS performed by Sergeyev & Carry (2021), hereafter S21, who released more than one million observations of almost 380 000 small bodies, multiplying by three the number of objects with respect to the last release of the Moving Objects Catalog.
The S21 catalog includes the point spread function magnitudes in the ugriz filter system and their respective uncertainties, σ m .We will use italics when discussing magnitudes in the text.S21 also includes the heliocentric and topocentric distances of the objects and the phase angle (along with other helpful information).As in Alvarez-Candal et al. (2022a), AC22 hereafter, we do not make any cut in the database other than removing any m with an associated error of σ m ≥ 1.We use data of objects with at least three different α and with a minimum span of ∆α ≥ 5.0 degrees (i.e., the distance between the maximum and the minimum α). Figure 1 shows the distribution in the orbital elements space of all objects with at least one H measured in any filter (see below).
To analyze the data, we built upon the algorithm developed in AC22, with a few modifications intended to improve accuracy, as discussed in the next section.

Analysis
We use the method described in AC22, updated as described below.The HG * 12 model (Penttilä et al. 2016) is used as the photometric model because it proved to be better suited to fit the data and because alternative models presented different problems when applied to the S21 dataset.The HG 1 G 2 model (Muinonen et al. 2010) needs well-sampled phase curve, which is usually not the case of the data analyzed in this work, the Shevchenko model (Belskaya & Shevchenko 2000) needs a good coverage of the phase curve in small phase angles; and, the HG model (Bowell et al. 1989) fails to represent phase curves of high or low-albedo asteroids.In this case, we did not use the Python adaptation of the HG * 12 model provided by the authors but a code based on the equations presented in Penttilä et al. (2016).
The adopted photometric model describes the phase curves as follows: (2) where Probability distribution function, PDF, of ∆m We followed a similar approach as in AC22, with one modification.Alvarez-Candal et al. (2022a) did not include any effect due to the possible rotational phase at which the objects could have been observed; this produced an artifact, a dip, in P A (∆m) as seen in Fig. 2 of AC22.The dip is alleviated when convolving with the errors in the apparent magnitudes, but it still assigned less probability than real to the region around ∆m = 0. We changed the code to account for this by multiplying a random sin(ϕ), where ϕ ∈ [0, 2π), simulating an observation at an unknown rotational phase.Equation 2 in AC22 is then re-written as where P A (∆m) = {P i A (∆m)} is the PDF of ∆m for object A in 150 bins of width 0.02 mag, meant to cover the whole range of reported ∆m ∈ [0.0, 3.0] (Warner et al. 2019).The denominator is a normalization factor.With this procedure, we created a PDF of possible values of ∆m that the object could have (see an example in Fig. 2 for asteroid 290301 ).
The rest of the process is the same as in AC22.Quoting: "To create the final PDF used to estimate the final range of possible solutions of H in all filters, we must also include the uncertainty in the measurement of m.In this case, we construct for each object and measure the distribution, P A (m), as a normal distribution with σ = σ m centered at m.The final PDF is the convolution of P A (m) and Eq.3: One vantage of including a factor sin (ϕ) in Eq. 3 is that it is not necessary to mirror and re-scale P A (m, ∆m) as AC22 needed because now the distribution naturally covers the whole space of possible ∆m.Computation of the phase curve As in AC22, we compute H and G for each object and filter using Eq. 2 generating different solutions of the PCs extracting 10 000 values randomly from P A (m, ∆m), the method is explained in full in AC22.
The outcome of the processing is the PDF of H λ and their corresponding Gs if there are at least 100 valid, that is, physically possible solutions.All PDFs are available for the readers upon request or in the Open Science Framework (OSF) project folder https://osf.io/43eny/.Table 1 shows a sample of the nominal values of the quantities and their uncertainties.The former is the median of the PDFs, while the uncertainty range indicates the 16th and 84th percentiles.The complete table is available online, but without format, in the OSF folder and will be available and formatted at the CDS upon acceptance.

Results and discussion
Following the above procedure, we computed H in all filters available in the S21 data (see numbers in Table 2) and also in the H V and H R using g, r, and i for each α when the magnitudes were available.The data was processed the same way as all the rest using V and R and their corresponding α.
Most objects have less than ten observations, with a median of four per filter.Noteworthy, over 100 objects have more than 20 observations, and one object was observed 70 times in the g, r, and i filters (Fig. 3, left panel).Most of the objects were observed with α min < 10 deg with a median span, ∆α, of 8.3 deg (Fig. 3, right panel).The vast majority of objects are main belt asteroids (Fig. 1), but there are objects with a > 100 AU at high eccentricities to fulfill the criteria set above to compute the PC.
The distributions of absolute magnitudes and G are shown in Fig. 4. As already seen in AC22, the H u distribution peaks at about two magnitudes fainter than those in the other filters (compatible with the solar (u−g) = 1.43).The G distributions all peak at about 0.6 with a FWHM of about 0.2.There are no apparent changes with wavelength other than a slight width increase towards bluer wavelengths.
The computed H g − H i , H g − H r , and H i − H z are displayed in Fig. 5 in color-color diagrams using density plots to simplify the figures.Overall, the distributions follow the same trend as in previous works (for example Alvarez-Candal et al. 2022a;Colazo et al. 2022;Sergeyev et al. 2023).It is worth mentioning that several objects seem to escape from the usual phase space because they have large uncertainties in their colors (seen as faint clouds in Fig. 5), and care must be taken when using this database.

Spectral slope
Alvarez-Candal et al. (2022a) established that the method was trustful, and it provided comparable results with other methods used to compute H (for example, Oszkiewicz et al. 2011or Vereš et al. 2015).Therefore, we jump straight into the central issue of this article: Phase coloring.
We compute first the spectral slope (slope for short) at different α for all objects with H u , H g , H r , and H i using a simple linear fit to the computed fluxes (see Eq. 5).H z was excluded from the fit because it falls on the silicate absorption feature in some objects, departing from the linear part of the spectrum.
The fluxes are computed using Notes.The first column shows the ID of the object and its Minor Planets Center ASCII designation.H u is the median of the probability distribution, while − and + designate the absolute magnitude at the 16th and 84th percentiles.The number of observations considered is shown in the fifth column.G u , G − u , and G + u from column six to eight.The last two columns show the minimum α and its total span.Flag -8 indicates objects with less than three observations in a given filter or ∆α < 5 deg, while flag -9 indicates that not enough valid solutions were obtained (that is, less than 100).The complete catalog is available at the CDS (https://cdsarc.u-strasbg.fr/) or upon request.where the solar colors were taken from the SDSS 2 .Note that in Eq. 5, it appears the object's color as an explicit function of the phase angle.The color was computed using the nominal values of Hs and Gs, feeding them into Eq.2. The computation was done in steps of 0.5 deg between α ∈ [0, 35] deg. Figure 6 shows the results for 42 168 objects.The y-axis shows a normalized spectral slope to ease comparison.This slope is computed as where ∆S = max[S (α) − S 0] − min[S (α) − S 0], and S 0 = S (α = 0 deg).It is immediately clear that there are two types of curves: (i) some slopes decrease with α until a critical value α c , and then revert the behavior (the BR group); (ii) other curves increase the slopes until a critical value, and then starts decreasing (the RB group).α c seems to lie between 4.5 and 5 degrees in all cases.
Behavior BR occurs for 22 858 objects and RB for 19 310.Note that there are no flat curves.Therefore, an object falls either on the BR or the RB group.In practice, if the normalization factor in Eq. 6 is removed, some curves will appear flat due to the small amplitude of ∆S .To highlight the maximum amount of change in S (α), we plot a histogram with the distribution of ∆S , as defined above, in Fig. 6 (left).It is clear that most objects do not display a radical change of S (α); the median value is 0.47 × 10 −4 nm −1 , within usual uncertainties in the measurements of spectral slopes.Nevertheless, it is important to note that there are almost 300 objects with ∆S > 5 × 10 −4 nm −1 , which is well above typical uncertainties.In any case, we would like to stress that the results displayed in Fig. 6 should be taken as a toy model representation of the results, as it does not include the full information in the probability distributions, although it is useful to glimpse the existence of different phase-coloring behaviors.The above results were obtained using the modeled data, which does not account for realistic observations.To add more "reality", we did a second experiment using the full probability .Change of spectral slope relative to the slope at α = 0 deg and normalized to its maximum amplitude to increase contrast using the full P(H λ ).In the continuous red line, we show the curve for 29030, whose colors make it possibly an S-complex asteroid, see below, while in the dot-dashed black line, we show the curve for 37478, a possible C-complex asteroid, according to its colors.
distribution of each absolute magnitude, extracting 1 000 different random values from each H and its corresponding G, obtaining a distribution of slopes per object and α.In this case, Fig. 6 cannot be reproduced because of the noise in the solutions.Nevertheless, inspecting some of the resulting curves visually, it is possible to detect that the behavior is discernible in high-quality data (Fig. 7).
Given the results shown in Fig. 6 and 7, we decided to separate the sample into two regimes and compute the change of α and S 0. The Spearman test shows that there is a weak negative anti-correlation with a r s = −0.01 and a p-value of 0.0026, in a way resembling (though much weaker) that detected for asteroids at small phase angles by Alvarez-Candal et al. (2022b), using a different photometric model.The right panel of the figure shows that above 4.5 deg, the general tendency is for objects to become redder when they are intrinsically redder in opposition.
In the next step, we checked for the coloring behavior in different bins of semi-major axis.The data was split into eight regions: (i) a ≤ 2.05 AU, (ii) 2.05 ≤ a ≤ 2.5 AU, (iii) 2.5 ≤ a ≤ 2.82 AU, (iv) 2.82 ≤ a ≤ 3.3 AU, (v) 3.3 ≤ a ≤ 3.7 AU, (vi) 3.7 ≤ a ≤ 4.5 AU, (vii) 4.5 ≤ a ≤ 5.5 AU, and (viii) a > 5.5 AU.We computed the average S (α) − S 0 value within each region per phase angle bin.The results are shown in Fig. 9.The figure clearly shows that the behavior of the average relative slope is not monotonic, neither with phase angle nor in the semimajor axis.Most curves show a reddening of the slopes up to 4.5 deg, then the relative slope decreases.The only "odd" behavior occurs in the semi-major axis occupied by the Jovian Trojans, 4.4 − 5.5 AU, where there is bluing and then reddening of the relative slopes.We checked that the distribution of α is similar in all semi-major axis bins to avoid introducing undesired biases.
The uncertainties of the quantities were not considered to avoid blurring the results too much.Note, however, that uncertainties were computed for the realistic slope as the standard deviation of the distributions and are available for the reader in the OSF folder.

Colors and phase angle
Not only the evolution of the slope with α is interesting.We also studied the change of the colors with the phase angle.Following a similar procedure as above to compute the more realistic slopes, we use the average relative color in the same semi-major axis bins as above and show the results for the H i − H z in Fig. 10 and for H g − H i in Fig. 11 The behavior of H i − H z indicates that, on average, the color decreases until about 5 deg in all semimajor axis bins, and it increases, getting redder with increasing α.Note that objects of all taxa are included in the averaging.
On the other hand, H g − H i seems to follow the curves of the average relative slope (Fig. 9), which is natural since this region is less affected by the 900 nm absorption band.Overall, the color seems to increase until the critical angle, and the color starts decreasing after α c .The same trends are seen in all semimajor axis bins, while the Hilda and Trojan regions have the less dramatic changes.
As mentioned, the uncertainties of the quantities were not included to avoid blurring the results.Nevertheless, the uncertainties were computed as the interval between the 16 th and the 84 th of the convolution of the PDFs of the colors, as explained in Sect.3.2 of Alvarez-Candal et al. (2022b).These uncertainties are available for the reader in the OSF folder.

The impact on the taxonomy classes
The previous sections established that spectral slopes and colors change with phase angle.The next step is to check whether these changes are sufficient to change the taxon of an object.To do so, a simple experiment was performed using the four taxa defined by Colazo et al. (2022), i.e., C, X, S, and V complexes.We crudely separated the H g − H i vs. H i − H z into the four regions each taxon covers (Fig. 12 provides a schematic view).0.6 0.8 1.0

Complex
N 0 The C-complex is located in We computed the number of objects in each box at α = 0 deg, N 0 , the number N c at α c , and the number of objects, N m , at the maximum phase angle (α max ).Notice that, according to Figs. 10 and 11, the maximum differences in color should happen between α c and α max .The results are shown in Table 3, where the differences are expressed in per mille.
As expected, the largest changes happen between N c and N m .Curiously, the V-complex is the only one that loses objects from opposition to α c to increment their numbers from α c to α max .In any case, this is a rough estimate that provides minimum numbers because the space covered by each box is larger than the actual taxon should occupy as well as we did not follow any particular object but just counted the number of objects at each step; therefore, it is likely that between steps more objects moved in or out of each region.(See also Appendix A.)

On the use of weighted averages
Weighted arithmetic averages are commonly used to compute colors (and other quantities) of Solar System objects (for example Hainaut et al. 2012;Sergeyev et al. 2023).The idea behind this approach is to prioritize data with high precision, the usual weight being the inverse of the uncertainty (or square uncertainty).This approach is valid if, and only if, all data are compatible.In this case, by compatible, we mean if all measurements were taken at the same phase angle, distances to the Sun and Earth, and rotational phase.Otherwise, some systematic may be introduced as objects tend to get brighter closer to the opposition, presumably increasing the precision of the measurements.We looked into S21's data and used the same criteria as mentioned above to select the data, but in this case, computing the weighted average in u−g, g−r, r−i, and i−z and then comparing these averaged value with the individual colors at each value of α and saved the phase angle (called α w ) where they were closer.Then, we compared the value of α w with the minimum phase angle for each object.The results displayed in Fig. 13 correspond to the g − r color measurements and are similar in the other colors.In the figure, it is clear that using a weighting scheme biases the colors towards low-phase angles (i.e., closer to opposition): 50% of the average colors coincide with the minimum α in the phase curve within 1 deg, while 3/4 of the sample are within 5 deg of the minimum phase angle.Therefore, we recommend care when using averaged values because these values are biased toward low-phase angle measurements carrying more weight.To further complicate the interpretation of results obtained with average colors, Alvarez-Candal et al. (2022a) showed that the average colors and the colors at opposition (i.e, the difference of absolute magnitudes that we call absolute colors) can have significant differences and that there is no one-to-one correspondence between them.

Conclusions
The terminology Phase reddening has been used for a while to describe an observational effect: objects observed at higher phase angles tended to have redder colors (or larger spectral slopes Luu & Jewitt 1990, for instance).Nevertheless, this was usually seen in a relatively low number of objects due to the inefficiency of traditional observation programs in scheduling repeated observations of the same object in different phase angles.For example, Nathues (2010), from a sample of over 90 Eunomia objects, only have six objects with two or more observations, including one with a clear bluing: (5785) Fulton.
The results we show in this paper point to different behaviors in different phase angles regimes: In terms of spectral slope, there seems to be a slight anti-correlation between S ′ α computed for α ≤ 4.5 deg and S 0, suggesting that intrinsically redder objects tend to get slightly bluer with increasing phase angle, while the opposite is seen for intrinsically bluer objects (see Fig. 8), these results are strikingly similar to these of Ayala-Loera et al. (2018) for TNOs and Alvarez-Candal et al. (2022b) for asteroids.On the other hand, S α for α > 5 deg shows, on average, that intrinsically redder objects tend to get redder in line with the traditional view of phase coloring.
There also seem to be two behaviors in terms of spectral slopes with respect to α (Figs 6, 7, and 9): The objects tend to have one behavior in the range [0, 4.5] deg, while the opposite in the range (4.5, 35] deg.There seems to be a slight predominance of objects that increase the spectral slope and then get bluer, the RB group.This behavior happens in almost all semimajor axis bins, but the effect is within error bars in typical spectral slopes measurements (the unit [10 4 nm −1 ] is equivalent to [% (1000 Å) −1 ]).We also cross-matched the data with the taxa determined by Colazo et al. (2022) to test if the two behaviors could be related to surface composition without finding any evidence.In the Appendix B we discriminate the behavior of the RB and BR groups to analyze their different demeanor.
Color-wise, there are similar behaviors: first, going in one direction to suddenly change in the opposite.The H g − H i correlates with the slope, which is expected because the i filter falls in a region relatively unaffected by absorption features.We find the behavior of the H i − H z more interesting because it targets the absorption band at 900 nm due to silicates.In this case, the color goes bluer for α ∈ [0, 4.5] deg, while the average color increases for α > 4.5 deg.Therefore, the band gets deeper or shallower in objects with an absorption feature just by phase coloring.Joining the results for both colors, H g − H i and H i − H z , there seem to be objects to get more concave for low-α and less concave for higher phase angle.Once again, we stress the orders of magnitude of the effect, which is at most about one-hundredth of magnitude.
One interesting point to consider is that α c coincides roughly with the onset of the non-linear part of the PC, where lies the opposition effect (OE).The OE happens due to a combination of shadow-hiding (Hapke 1963) and coherent back-scattering (Muinonen 1989) and produces a net increase in the brightness of the object when observed close to opposition (α ≈ 0 deg).The effect starts to appear below 5 to 6 deg (Belskaya & Shevchenko 2000), but it can start as close to the opposition as below 1 deg (for example, in trans-Neptunian objects, see, Verbiscer et al. 2022).The apparent coincidence is very suggestive, and it deserves further analysis beyond this work's scope.Perhaps is it here were the physical explanation for the phenomena described in this work lie: which scattering mechanism dominate and how does it behave with the wavelength.Unfortunately, we cannot provide a satisfactory answer yet.
The effects we report in this paper are unlikely to change how we see the distribution of small bodies' taxa in the Solar System (for example Mothé-Diniz et al. 2003;DeMeo & Carry 2013;Marsset et al. 2022;Sergeyev et al. 2023), at least not for the core of the groups, but it does add noise to the frontiers between taxa where some objects do change taxa.The numbers are a priory small, in the range of a few dozen per thousand.Nevertheless, with the arrival of massive multi-wavelength photometric surveys, this may change, as the precision in the measurements may reach the level of those shown in this work or close to it.On the other hand, with the large numbers, we can now statistically study the phase-related behavior of many objects, making possible analysis that could not have been done before.

Fig. 1 .
Fig. 1.Osculating orbital elements of all objects analyzed in this work with at least one H computed in any filter.

Fig. 3 .
Fig. 3. Observational circumstances of the data used in this work.Left panel: Cumulative distribution of the number of observations per filter for the objects in this work.Right panel: Minimum α vs. span in α covered by the objects in this work.

Fig. 4 .Fig. 5 .
Fig. 4. Distributions of computed quantities.Distributions of H (left panel) and G (right panel).u is shown in the blue line, g in green, r in red, i in purple, and z in black.

Fig. 6 .
Fig. 6.The two behaviors of the spectral slope.Left: Change of spectral slope relative to the slope at α = 0 deg and normalized to its maximum amplitude to increase contrast.The continuous line shows the average value of all objects with that behavior, while the dashed lines indicate the maximum and minimum values per α bin.Right: Distribution of ∆S .

Fig. 10 .Fig. 11 .
Fig. 10.Average relative H i − H z in bins of semi-major axis.The results for each bin are marked within the figures.

Fig. 12 .
Fig. 12. Schematic view of the boxes considered for each taxon.
is the reduced magnitude, H is the absolute magnitude, Φ i are known functions of α (see Penttilä et al. 2016 for the tables), and G i are the phase coefficients that provide the shape of the curve.In the HG * 12 model, G 1 and G 2 are parameterized by G * 12 according to G 1 = 0.84293649G *

Table 1 .
Sample of the catalog.

Table 2 .
Number of absolute magnitudes obtained

Table 3 .
Relative change of taxa