Asteroid Photometric Phase Functions From Bayesian Lightcurve Inversion

Photometry is an important tool for characterizing the physical properties of asteroids. An asteroid’s photometric lightcurve and phase curve refer to the variation of the asteroid’s disk-integrated brightness in time and in phase angle (the Sun-asteroid-observer angle), respectively. They depend on the asteroid’s shape, rotation, and surface scattering properties, and the geometry of illumination and observation. We present Bayesian lightcurve inversion methods for the retrieval of the asteroid’s phase function, the unambiguous phase curve of a spherical object with surface scattering properties equal to those of the asteroid. A collection of such phase functions can give rise to a photometric taxonomy for asteroids. In the inverse problem, first, there are four classes of lightcurves that require individual error models. The photometric observations can be absolute or relative and they can have dense or sparse cadence in comparison to the rotation period of the asteroid. Second, the observations extend over varying phase angle ranges, requiring different phase function models. Asteroid photometry from the European Space Agency Gaia space mission extends, typically, over a range of phase angles, where the phase curve tends to be linear on the magnitude scale. Photometry from ground-based observing programs can reach small phase angles, where the asteroids show an opposition effect, a nonlinear increase of brightness on the magnitude scale towards zero phase angle. We provide error models for all four classes of lightcurves and make use of linear or linear-exponential phase functions for phase angles below 50°. We apply the inverse methods to sparse absolute Gaia and dense relative ground-based lightcurves and obtain absolute magnitudes and phase functions, with uncertainties, for ∼500 asteroids. Finally, we assess the lightcurve inversion problem for dense absolute photometry with the help of a numerical simulation for a Gaussian-random-sphere asteroid.


INTRODUCTION
Asteroids are small Solar System objects that originate from times preceding planet formation. Typically, when rotating about its principal axis of inertia, an asteroid exhibits a periodic change of brightness caused by the varying part of its surface being both illuminated by the Sun and visible to the observer. On one hand, a photometric lightcurve is the result of photometric observations extending over a time span covering a substantial part of the rotation period or more. On the other hand, a photometric phase curve is the result of observations of the change in an asteroid's apparent brightness obtained at different epochs during a single apparition, corresponding to a slowly changing Sun-observerobject geometry.
Consequently, photometric lightcurves provide immediate information about the rotation period. Historically, asteroid lightcurve observations have been the first tool intensively applied to derive physical properties of these objects. In addition to the retrieval of the rotation period, lightcurve amplitudes have been used to derive rough characteristics of the shape (elongated shapes generally producing lightcurves with larger amplitudes) and, when lightcurves obtained at different apparitions were available, to derive estimates of the direction of the asteroid's pole orientation (Taylor, 1979;Magnusson et al., 1989;Kaasalainen et al., 2002;Durech et al., 2015). Some early statistical analyses of lightcurve amplitudes and periods led also to the discovery of the existence of equilibrium shapes among large asteroids (Farinella et al., 1981), and opened the way to the discovery of the importance of the general phenomenon of asteroid collisional evolution (Farinella et al., 1982).
Photometric phase curves provide information about the intrinsic light-scattering properties of the surface that are intimately related to the regolith composition and structure. Since the composition and structure determine also the reflectance properties observed at different wavelengths, phase curves are strictly related to the taxonomic classification of asteroids. Phase curves are often represented with the help of a few parameters. The classical two-parameter H, G magnitude system has been used for a long time (Bowell et al., 1989). H is the asteroid's absolute magnitude, namely the magnitude (corresponding to unit distance from the Sun and the observer) measured at zero phase angle, and G is a parameter describing the overall variation of magnitude at different phase angles, including a nonlinear brightness surge at phase angles smaller than 10°(opposition effect). In recent years, it has been replaced by the so-called H, G 1 , G 2 magnitude system, a threeparameter model developed by Muinonen et al. (2010) to remove the caveats of the H, G system in the case of low-albedo and highalbedo asteroids. Further refinements and applications of the H, G 1 , G 2 system have been published, among others, by Penttilä et al. (2016) and Shevchenko et al. (2016).
Gaia Data Release 3 (DR3) is imminent (13 June 2022), including extensive astrometric, photometric, and spectroscopic observations of small Solar System objects, primarily asteroids. It is important to prepare to analyze such a large amount of new data. In the field of asteroid photometry, in particular, significant progress has been made recently in what concerns the capability of obtaining efficient processing and reliable interpretation of sparse photometric data. Martikainen et al. (2021) carried out an analysis of photometric data combined from ground-based lightcurves and observations published in Gaia Data Release 2 (DR2) in order to invert the data for reliable estimates of shape, rotational properties, and phase curves for a large number of objects belonging to a variety of asteroid taxonomic classes. They derived photometric phase curve slopes, rates of brightness change on the magnitude scale, for more than 300 asteroids in the so-called reference geometry of equatorial illumination and observation (Kaasalainen et al., 2001). Expanding on the study in Muinonen et al. (2020), Martikainen et al. (2021) provided unequivocal proof that the projection to similar illumination and observation conditions was needed to enable unbiased comparative studies of asteroid phase curves. It is possible to strive towards minimizing the biases by incorporating all practical geometries of illumination and observation (Oszkiewicz et al., 2011). However, for asteroids with their individual pole orientations and orbits, different illumination and observation geometries are sampled, and biases remain. In spite of these biases, Oszkiewicz et al. (2011) and Mahlke et al. (2021) have successfully related phase curve parameters from massive observing programs to the taxonomical classes.
We extend the work by Martikainen et al. (2021) by providing computational tools for the derivation of asteroid phase functions using fictitious spherical asteroids with equal surface properties. In earlier works, we have described these phase functions as being proper phase functions that describe the intrinsic properties of the surfaces. We apply the methods to the asteroids studied by Martikainen et al. (2021), by starting from their results of Markov-chain Monte Carlo (MCMC) lightcurve inversion for some 500 asteroids. Furthermore, we generalize the studies by Martikainen et al. (2021) and Muinonen et al. (2020) by incorporating a combined linear-exponential model of phase functions on the magnitude scale.
In earlier work (Muinonen et al., 2020), models of observational uncertainties were developed for dense relative photometry and sparse relative photometry. The former entailed ground-based lightcurves that were treated, in lightcurve inversion, on a relative magnitude scale. The latter comprised lightcurves of sparse Gaia photometry that were incorporated on a relative magnitude scale, too. Martikainen et al. (2021) then treated the Gaia photometry in the absolute sense, deriving absolute magnitudes for a large number of asteroids. In the present work, we provide a complete set of four models for observational uncertainties, including models for dense relative, sparse relative, dense absolute, and sparse absolute lightcurves.
The paper is organized as follows. Section 2 describes the theoretical framework in asteroid lightcurve inversion. Particular attention is paid to the different surface scattering models utilized in the forward and inverse problems. In Sect. 3, error models are presented for dense and sparse lightcurves (relative or absolute), and the retrieval of absolute magnitudes and phase functions is outlined. Section 4 first provides the application of the methods tõ 500 asteroids with both sparse absolute Gaia photometry and dense relative ground-based photometry. Thereafter, a numerical forward simulation of dense absolute photometry with a Gaussian-random-sphere asteroid model (GS-asteroid) is described, followed by inversion results and their discussion. Conclusions and future prospects are offered in Sect. 5.

Asteroid Modeling
Consider an asteroid in principal-axis rotation. We denote the asteroid's rotation period by P, the pole orientation in ecliptic longitude and latitude by (λ,β) T (J2000.0, T stands for transpose), and the rotational phase at a given epoch t 0 by φ 0 . As to the asteroid's shape, in the inverse problem, we consider both ellipsoidal shapes and general convex shapes. In the forward problem, we incorporate star-like shapes described by the Gaussian-random-sphere geometry. The apparent V-band magnitude (in mag) of a fictitious spherical asteroid (diameter D in km and geometric albedo pV) at a given phase angle α is where H is the absolute magnitude (in mag), r and Δ are the asteroid's heliocentric and topocentric distances (in au), Φ is the phase function (Φ(0°) = 1), and lg denotes the 10-based logarithm. The corresponding reduced V-band magnitude is projected to r = Δ = 1 au: Scattering of light from a surface element on an asteroid is described by the diffuse reflection coefficient R that relates the incident solar flux density πF 0 to the emergent, scattered intensity I: Here ι and ϵ are the angles of incidence and emergence as measured from the outward normal vector of the element, and ϕ 0 and ϕ denote the respective azimuthal angles. For a geometrically isotropic surface, it is unnecessary to specify ϕ 0 and we set the coordinate system so that ϕ 0 = 0°. Consequently, the backscattering direction, the direction for the source of light, is with ϕ = 0°. The Lommel-Seeliger surface reflection coefficient (subscript LS) derives from radiative transfer (e.g., Lumme and Bowell, 1981): where p is the geometric albedo for the wavelength band considered.
The diffuse reflection coefficient for dark particulate media (subscript PM)-such as planetary regoliths of low-albedo asteroids-can be expressed in the form where p is the geometric albedo, α is the phase angle, andω and P 11 are the single-scattering albedo and single-scattering phase function (different from the phase function of the asteroid). The function Φ S represents the corrections to radiative transfer due to the dense packing of the particulate medium, for example, the corrections due to shadowing among the particles typically much larger than wavelength. Φ 11 and Φ S , as all functions Φ of the present study, are normalized to unity at backscattering: The reflection coefficient in Eq. 5 belongs to a class of photometric models consisting of a Lommel-Seeliger-type volume-element part and a part describing scattering among volume elements in a particulate medium (e.g., Lumme and Bowell, 1981;Muinonen and Lumme, 1991).  provide an extensive set of numerical computations of the dense-packing correction Φ S using particulate media of opaque spherical particles. They complete the modeling with the help of a fractional-Brownian-motion model for the roughness of interface between the medium and free space. There are altogether three parameters: the packing density of the particles v, the fractal Hurst exponent H fBm , and the amplitude σ fBm . Smaller H fBm and higher σ fBm imply rougher interfaces with stronger effects of interface roughness. The function Φ S is an unknown function for asteroid surfaces. In what follows, as in Muinonen et al. (2015), we utilize Φ S estimated for the lunar mare regolith by Wilkman et al. (2014).
Consider next the function Φ 11 in Eqs 4, 5 and the large numbers of existing phase curve observations of asteroids. Indeed, we may consider that the principal form for the phase function of a fictitious spherical asteroid is already known to be described by, for example, the H, G 1 , G 2 phase function. Consequently, we introduce denominators in Eqs 4, 5 that cancel the inherent phase function that would result from the assumption of Φ 11 = 1 (isotropic scattering): where Φ represents an unknown function that should nevertheless be close to the empirically known phase functions of asteroids. In other words, the functions Φ LS and Φ PM are the phase functions for a spherical asteroid with reflection coefficients R LS and R PM assuming Φ 11 = 1, respectively. . Their ratio is also given. By approximating that Φ PM can be valid for the higher-albedo lunar highland regions, too, the ratio points to the part of the lunar opposition effect at small phase angles unexplained by shadowing and rather explained by the coherent backscattering mechanism (Muinonen et al., 2002). Significantly different single-scattering phase functions Φ/Φ LS or Φ/Φ PM result from utilizing R LS or R PM : as R PM assigns a significant part of phase curve steepness to shadowing, the resulting single-scattering phase function is shallower. In order to model the function Φ in Eq. 7 for lightcurve inversion, we make use of a linear-exponential model on the magnitude scale: where m 0 and α 0 are the amplitude and angular width of the opposition effect, respectively, and β 0 is a slope parameter. In the linear-exponential model of Eq. 8, the photometric slope at α = 20°(the phase angle chosen for comparative studies) equals For phase angles outside the angular regime of the opposition effect, 10°≤ α ≤ 50°, a linear model can be utilized: The H, G 12 phase function is a two-parameter phase function developed for scarce photometric data (Muinonen et al., 2010). As in Martikainen et al. (2021), we rule out increasing brightness with increasing phase angle and enforce the absence of an opposition effect for G 12 values that would result in negative weights 1 − G 1 − G 2 : The extension does not affect the H, G 12 phase function within its approximate nominal range of 0 ≤ G 12 ≤ 1. Frontiers in Astronomy and Space Sciences | www.frontiersin.org June 2022 | Volume 9 | Article 821125

Inverse Problem
Consider the free parameters (unknowns) of the asteroid forward model. In the case of the ellipsoid shape model, the parameters are described by the vector where the N P = 9 free parameters are, respectively, the rotation period, ecliptic pole longitude, ecliptic pole latitude, rotational phase at a given epoch t 0 , two ellipsoid axial ratios (with a > b > c denoting the semiaxes), and three parameters of the linearexponential phase function. In the case of the convex shape models, the parameters are where the shape is described by the (l max + 1) 2 − l 2 min spherical harmonics coefficients s 00 , . . . , s lmaxlmax for the Gaussian surface density. Here φ 0 and the real-valued s 00 can be fixed as the rotational phase is immersed in the spherical harmonics coefficients and the asteroid size information is omitted. Thus, the total number of parameters equals Let p p be the a posteriori probability density function (p.d.f.) for the parameters. Within the Bayesian framework (cf. Muinonen et al., 2020), p p is proportional to the a priori and observational uncertainty p.d.f.s p pr and p ϵ+υ , ϵ and υ referring to random and systematic uncertainties. p ϵ+υ is evaluated for the "Observed-Computed" (O-C) residual magnitudes ΔM(P), Even though p ϵ+υ is here related to the observations, it can also describe the uncertainties deriving from the shortcomings in the physical model. It is currently assumed that p ϵ+υ is Gaussian and that p pr will describe, for example, the regularization needed in convex inversion. The final a posteriori p.d.f. is thus where χ 2 measures the O-C distance in terms of the model for the uncertainties. The observation vector is composed of a number of lightcurves with their varying numbers of magnitudes, and the uncertainties are assumed to be uncorrelated between the lightcurves. We may thus rephrase χ 2 (P) as where M obs,k , M k (P), and Λ ϵ+υ,k pertain to the observations, computations, and the covariance matrix for the uncertainties in lightcurve k, the total number of lightcurves being K.
In detail, we simplify the χ 2 -value in Eq. 17 to the form where the σ ϵ,k values describe the uncertainty (and weight) of the N k observations in lightcurve k, and M obs,kj and M kj (P) are the observed and computed magnitudes. For small relative uncertainties in brightness, the χ 2 -value can be approximated by (Muinonen et al., 2020) where ΔM k0 (P) denote the O-C difference of the mean magnitudes in lightcurve k, and ℓ obs,kj and ℓ kj (P) are the observed and computed brightnesses relative to the brightnesses corresponding to the lightcurve mean magnitude lg~e 0.43429. Thus, in Eq. 19, the χ 2 -value is computed using the differences in the observed and computed relative brightnesses, relative to the observed relative brightnesses.

NUMERICAL METHODS
The MCMC inverse methods are based on proposal probability densities characterized with the help of the so-called virtual leastsquares solutions and are described in detail by Muinonen et al. (2020). In what follows, we introduce the error models for the FIGURE 2 | Example photometric lightcurves for asteroid (167) Urda (blue circles) on the magnitude scale. We show the ground-based dense lightcurve #8 (left) and the Gaia sparse lightcurve (middle and right). Red crosses stand for the best-fit convex model lightcurves with the 1-σ uncertainties (red bars), and the red points are the uncertainty envelopes using 5,000 MCMC sample solutions. The symbols t, α, and k denote the time, phase angle, and lightcurve observation counter.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org June 2022 | Volume 9 | Article 821125 observations and the algorithm for the retrieval of absolute magnitudes and phase functions from DR2 photometry.

Error Models for Observations
We define dense lightcurves as having photometric points timed, on average, less than the asteroid's estimated rotation period apart. For dense relative lightcurves, the points are calibrated against each other but not against absolute standards. We set (Muinonen et al., 2020) where σ 0,k denotes the rms-value and σ pr,k denotes the a priori threshold value for the uncertainty in the lightcurve k, obtainable, for example, from cubic spline fits to the lightcurves. Furthermore, ΔT k is the mean sampling time interval of the lightcurve k, wherek marks the lightcurve with the lowest sampling rate in time and T k stands for the time span of the lightcurve k. Equation 20 can be interpreted in the following way. The lightcurvek contributes a χ 2 -value of unity in the inverse problem, whereas each other lightcurve k is considered to split into approximately N k,eff independent lightcurves. For dense absolute lightcurves composed of K lightcurves, the points are calibrated both against each other and absolute standards. We start by setting where σ 0,k and σ pr,k are defined as for Eq. 20. The definition of N k,eff is the key point of the error model. For simplicity, consider the case σ 0,k ≥ σ pr,k . First, if we set N k,eff = 1, each lightcurve obtains an equal weight of unity. This is the most conservative model with systematic errors dominating over the random errors. Second, we may set N k,eff = 1 for the lightcurvek that has the minimum number of observations and consider the other lightcurves as composed of an effective number of lightcurves like lightcurvek, balancing the weights of lightcurves with drastically different numbers of observations. We recall that the systematic errors include, in addition to the observational errors, the effects of simplified forward modeling, in particular, for the phase function in Eqs 8-10. We define sparse lightcurves as having photometric points timed, on average, more than the estimated rotation period apart. For sparse relative lightcurves, as for dense relative lightcurves, the points are calibrated against each other but not against absolute standards. In accordance with Muinonen et al. (2020), we set N k,eff Nk in an error model coinciding with the form in Eq. 22, the indexk now denoting the sparse lightcurve with the smallest number of observations. Finally, for sparse absolute lightcurves, we introduce a model that coincides with the  one for the sparse relative lightcurves. Now the observations are, however, calibrated against absolute standards. The complete a posteriori p.d.f. is the product of the p.d.f.s for the four classes so that the corresponding χ 2 -values are added together: or FIGURE 5 | The topmost graphs depict σ β S (20°) as a function of β S (20°) using 358 asteroids (left) and highlighting 13 new asteroids presently processed (right). The middle graphs show σ β S (20°) as a function of β S (20°) for 97 asteroids with known Tholen classes (left), highlighting six new asteroids (right). The bottom graphs provide a zoom-in of the middle graphs. Frontiers in Astronomy and Space Sciences | www.frontiersin.org June 2022 | Volume 9 | Article 821125 7 where the index i describes the dense relative (i = 1), sparse relative (i = 2), dense absolute (i = 3), and sparse absolute photometry (i = 4). For dense and sparse relative lightcurves, the relative brightnesses ℓ obs,ikj and ℓ ikj (P) are computed using the mean-magnitude brightnesses of each lightcurve. However, for dense and sparse absolute lightcurves, they are computed using the mean magnitude of the entire absolute lightcurve data. Note that the sparse relative lightcurves are weighted equally to the absolute photometric data but they enter the inverse problem separately with no regard to the absolute level of brightness.

Absolute Magnitudes and Phase Functions
We refine the algorithm in Martikainen et al. (2021) for the derivation of absolute magnitudes and phase functions from Gaia photometry. As in their study, first, we start from the phase function slope parameter β S retrieved from MCMC inversion, recalling that β S describes the intrinsic surface-element properties of an asteroid. Second, using the full asteroid model available from the inversion, we move to the reference     (1,0) and G EQ (1,0) for the reference geometry and equal-sphere case, respectively, the derived absolute magnitudes H G using the H, G 1 , G 2 fit, and the derived G 1 and G 2 parameters for the equal-sphere photometric phase functions. The uncertainties are given in units of the last digit shown.

Asteroid
Tholen  (74) FIGURE 8 | The distributions of the different Tholen classes (black dots) and their range (black crosses) in the G 1 and G 2 parameter space based on the equal-sphere phase functions. The grey dots represent single asteroids and the green line shows how the G 12 parameter maps into G 1 and G 2 in the H, G 12 magnitude system.  4 | Lightcurve characteristics for the simulated Gaussian-sphere asteroid (GS-asteroid). "Class" denotes the Tholen taxonomical class, N and K denote the numbers of observations and lightcurves, respectively, and T obs is the time span of the observations. The uppermost numbers refer to the simulated dense ground-based observations, the ones in the middle to the simulated sparse observations, and the lowermost ones refer to the combined observations. improvement, we fit, directly, the H, G 12 phase function to the magnitudes of the brightness maxima. The resulting H and G 12 parameters then allow us to predict the lightcurve brightness maxima at arbitrary phase angles within 0°-120°. Fourth, a reference phase curve is computed for selected phase angles by averaging the magnitudes over one full rotation. Fifth, magnitudes of lightcurve brightness maxima are computed for the selected phase angles, together with the values for the integrated disk function (∝ μ 0 /(μ + μ 0 ) in Eqs 3, 4) and singlescattering phase function, all in the magnitude scale for the geometries corresponding to the brightness maxima. Sixth, a phase function is computed for an equal-projected-area spherical asteroid with the help of the single-scattering phase function, made possible by the fact that the mean projected area of a convex object in random orientation is one fourth of its total surface area (van de Hulst, 1957). An absolute magnitude then follows from the prediction to zero phase angle. Seventh, we fit the equal-sphere phase function using the full H, G 1 , G 2 phase function. Finally, we compute the slopes of the mean-magnitude reference phase curve (β ref ) and the phase function (β S ) at the 20°phase angle. Repeating the computations above for all the MCMC solutions allows us to obtain uncertainty estimates for the absolute magnitudes and phase functions throughout the entire phase curve analysis.

RESULTS AND DISCUSSION
The inverse methods are here applied to the~500 asteroids with ground-based and Gaia DR2 lightcurve data and to the simulated lightcurve data of the GS-asteroid.

Asteroids of Gaia Data Release 2
Martikainen et al. (2021) carried out MCMC convex lightcurve inversion for some 500 asteroids with photometric data deriving from Gaia DR2. In the present study, we make use their 5,000 MCMC samples for each asteroid and revisit the photometric phase function retrieval with improvements described in Sections 2 and 3. The number of successful phase function retrievals increased and 13 new asteroids were incorporated into the analysis, bringing the total number of successful retrievals to 358. The improvements further allowed us to derive uncertainties for the reference phase curves and final equal-sphere phase functions and absolute magnitudes. Whereas Martikainen et al. (2021) documented photometric slopes at 20°phase angle for the reference geometry of equatorial illumination and observation, we proceed to provide the full photometric 6 | Lightcurve characteristics for the GS-asteroid in the case of relative photometry. First, we give the lightcurve identifier (k), effective number of observations (N k,eff , Eq. 20), and the resulting initial uncertainties (σ ϵ,k , "ini"). Second, we give the final uncertainties (σ ϵ,k , "fin") and rms-values (rms(m), "fin") for the best-fit model from convex inversion. phase function of a fictitious spherical object with surface scattering properties equal to those of the asteroid. Note that the earlier work addressed the estimation of the equal-sphere photometric slope across the phase-angle range of the Gaia observations, by incorporating the corresponding slope parameter in MCMC lightcurve inversion. As the phase-angle ranges differed for each asteroid, the slopes also differed a priori, sometimes substantially. The remedy utilized by Martikainen et al. (2021) was to move to the reference geometry and project the photometric slopes to the same phase angle of 20°. Due to the fact that, in the reference geometry, we replace the linear model by the H, G 12 phase function and repeat the analysis separately for each MCMC sample solution, we are capable of carrying out the analysis of uncertainties all the way to the equalsphere phase functions. Complete tables including all the present results for the Gaia DR2 phase functions are available from the authors. Figure 2 shows example lightcurves for asteroid (167) Urda together with best-fit convex model predictions and MCMC uncertainty envelopes. The sparse Gaia DR2 lightcurve is accurately modeled with realistic model uncertainties. Modeling for the dense ground-based lightcurve is satisfactory. Figure 3 presents histograms of the β S (20°) and β ref values that were obtained for 358 asteroids. To reiterate, β S (20°) refers to the photometric slope from equal-sphere phase function obtained at the 20°phase angle, and β ref refers to the slope of the mean-magnitude reference phase curve at the same phase angle. The β ref values show no significant differences compared to the values derived in the previous study by Martikainen et al. (2021). The values of β S and β ref range from 1.3 to 5.3 mag/rad and are almost identical which can be seen in Figures 3, 4 with the exception of two noticeable outliers. The asteroid (15172) 3086 P-L shows a deviation of 0.19 mag/rad. It has only 3 lightcurves with a total of 398 observations. The phase-angle range of the Gaia lightcurve is 8°. Asteroid (5902) Talima shows a slightly larger deviation of 0.23 mag/rad. It has only 3 lightcurves totaling 29 observations. Again, the phase-angle range of the Gaia lightcurve is 8°. For these objects, the small number of ground-based lightcurves gives rise to large uncertainties.
For the vast majority of the asteroids presently studied, the absolute magnitude referring to random orientation is brighter than the absolute magnitude referring to the reference geometry. The explanation derives from a bias in the reference geometry. For the sake of clarity, let us characterize the shape of the asteroid by the semiaxes a ≥ b ≥ c as in the ellipsoidal model. On one hand, in the reference geometry, an asteroid rotating about its axis c of maximum inertia will never be viewed with its maximum projected area pointing in equatorial directions, whereas, in random orientation, maximum-projected-area viewing is enabled. Consequently, the asteroid exhibits a larger mean projected area in random orientation than in the reference geometry and the resulting absolute magnitude is brighter. On the other hand, in the reference geometry, in the unlikely case of an asteroid rotating about its long axis a of minimum inertia, FIGURE 9 | Example photometric lightcurves for the Gaussian-random-sphere asteroid (GS-asteroid). We show the dense photometric lightcurves #12 and #5 and the sparse lightcurve (blue circles; left, middle, and right, respectively; see Tables 5-7) together with the best-fit convex model lightcurves (red crosses) with 1-σ standard deviations (red bars) and uncertainty envelopes using all 30,000 Markov-chain Monte Carlo (MCMC) sample solutions (red points). We depict the results for relative (top) and absolute lightcurves (bottom). m, t, and k stand for relative magnitude, time from the first lightcurve observation, and lightcurve observation counter, respectively.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org June 2022 | Volume 9 | Article 821125 there is a bias towards the asteroid appearing larger in the reference geometry than in random orientation. As most of the asteroids are rotating about their axis of maximum inertia, at the population level, the absolute magnitude is predominantly brighter in the case of random orientation.
In Figure 4, we compare the β S (20°) and β ref values. For every asteroid, the value of the β ref is marginally greater than or equal to the value of β S (20°). The explanation can again be related to a bias in the reference geometry. There, the minimum projected area of an ellipsoidal asteroid appears in equatorial viewing. Recalling the increase of lightcurve amplitude with phase angle assessed by Zappalà et al. (1990), the reference geometry provides the steepest phase curve dependence available for lightcurve brightness minima of a given asteroid. As to the lightcurve brightness maxima, phase curve dependences shallower than those in the reference geometry are available in random orientation, as geometries involving the maximum projected area of the asteroid are exposed. It is reasonable that the mean-magnitude reference phase curve turns out to be steeper than the meanmagnitude random-orientation phase curve. Conversely, for asteroids rotating about their axis of minimum inertia, the random-orientation phase curves are steeper than the reference phase curves.
The uncertainties for β S (20°) as a function of β S (20°) itself are depicted in Figure 5. Martikainen et al. (2021) showed that, despite of all the asteroid classes overlapping with each other to an extent, the C-complex and the S-complex asteroids form their own distinctive groups. This can be seen more clearly in our results: the C-complex asteroids have larger β S (20°) values and are located on the right, whereas the S-complex asteroids have smaller β S (20°) and are located on the left. The X-complex and end member asteroids are spread more widely as both groups are heterogeneous, containing asteroids with different albedos and compositions. In the middle and bottom panel of Figure 5, we focus for simplicity on the Tholen taxonomic classification, and do not check the possibility that some of the objects could belong to a different class according to more modern taxonomic classifications.
The relation between the derived G-band absolute magnitudes G(1,0) and the V-band absolute magnitudes H for the 358 asteroids is displayed in Figure 6. The V-band absolute magnitudes were retrieved from the Jet Propulsion Laboratory Small-Body Database. The figure shows that our derived G-band magnitudes are nearly perfectly correlated save for a few outliers. These are most likely caused by insufficient lightcurve data. Asteroid (1368) Numidia shows a deviation of 0.81 mag from its respective H value, and has only 4 lightcurves that totaled in 40 observations. The phase-angle range of the Gaia lightcurve is only 4°, which is bound to result in uncertainties. Asteroid (446) Aeternitas shows the largest deviation of 1.37 mag. Despite having 128 lightcurves totaling 1794 observations, the phase-angle range of the Gaia lightcurve is again just 4°. Table 1 describes the taxonomical classes, numbers of observations, and phase-angle ranges for 12 example asteroids TABLE 8 | Rotation period P (h), pole longitude λ (°) and latitude β (°), rotational phase ϕ (°), and phase-curve parameters m 0 (mag) and β 0 (mag/rad) for the simulated Gaussian-sphere asteroid. The parameter values on the topmost line (GS-model) represent those utilized to simulate the observations. For the simulation of the scattering model, see text. On the lines that follow, the parameter values represent least-squares solutions (LS) and means from MCMC sampling (MCMC) as obtained from ellipsoid inversion (EI) and convex inversion (CXI) in the cases of relative (Relative) and absolute simulated observations (Absolute). For the MCMC solutions, the parameter 1-σ standard deviations are given in the units of the last digit shown.  (47) TABLE 9 | Correlation coefficients for the rotation period P (h), pole longitude λ (°) and latitude β (°), rotational phase ϕ (°), and phase-curve parameters m 0 (mag) and β 0 (mag/rad) for MCMC lightcurve inversion using ellipsoids (EI for ellipsoid inversion) and convex shapes (CXI for convex inversion) in the case of the simulated GS asteroid. The upper (lower) triangles depict the coefficients in the case of relative (absolute) simulated observations.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org June 2022 | Volume 9 | Article 821125 FIGURE 10 | MCMC sampling of the rotational, shape, and photometric phase-curve parameters for the GS-asteroid simulation in the case of relative lightcurves. The marginal probabilitity densities are characterized by an unbiased sub-sample of 1,000 parameter sets selected from the complete sample of 30,000 sets. We depict rotational pole longitude (λ) and latitude (β) vs. rotation period (P; up to the left and middle, respectively), pole longitude vs. pole latitude (up to the right), example spherical harmonics coefficients a 20 and b 21 vs. rotation period (down to the left and middle, respectively), and photometric slope (β S (20°)) vs. rotation period (down to the right). The red, dotted line marks the true rotation period. Figure 10 in the case of absolute photometry.

FIGURE 11 | As in
Frontiers in Astronomy and Space Sciences | www.frontiersin.org June 2022 | Volume 9 | Article 821125 13 chosen for a more detailed study as in Martikainen et al. (2021). Thereafter, Figure 7 showcases the random-orientation phase functions along with the H, G 1 , G 2 fits for asteroids (55) Pandora, (95) Arethusa, (245) Vera, (596) Scheila, (731) Sorga, and (1251) Hedera. These asteroids were chosen to demonstrate different Tholen classes and to compare our results to the phase functions shown in Martikainen et al. (2021). The H, G 1 , G 2 phase functions were fitted using the phase angles of 0.0°, 0.3°, 1.0°, 4.0°, 8.0°, 12.0°, 20.0°, 30.0°, and 50.0°. A noticeable change to results by Martikainen et al. (2021) comes up in the order of the phase functions in terms of their slope at 20°. In particular, the asteroids (245) Vera and (1251) Hedera have exchanged places. We are still involved in analyzing scarce Gaia DR2 data and expect to be able to offer more solid conclusions based on the forthcoming Gaia DR3 data. Table 2 provides a detailed comparison, including uncertainties, of the phase function slopes β S , β ref , and β S (20°). It is clear that the latter two are closer to one another than to β S that has been computed over a different phase-angle range for each asteroid. When comparing to phase-angle ranges documented in Table 1, a systematic behavior is seen: for phaseangle ranges located towards phase angles smaller than 20°, those local slopes are higher than the slope β S (20°). This derives from the smoothly onsetting opposition effect prominent at small phase angles. For phase-angle ranges centered near 20°, the two slopes are similar. Finally, asteroid (245) Vera shows a phase function strikingly similar to that of the Moon. In Table 3, absolute magnitudes with uncertainties are collected for all 12 example asteroids. For most of the asteroids, the reference absolute magnitude is fainter than the random-orientation one, but there are two exceptions, (246) Asporina and (1251) Hedera. Whether these asteroid rotate about their axes of minimum inertia remains as an open question that may well be resolved by the forthcoming Gaia DR3 data.
We present the derived G 1 and G 2 parameters for the different Tholen classes in Figure 8 together with the G 12 parameter mapping to G 1 and G 2 in the H, G 12 magnitude system. It is possible to distinguish three groups in the G 1 , G 2 parameter space: the first group contains the M and S classes, the second group contains the D class, and the third group consists of the C, G, and F classes. The figure also includes A, P, X, and B classes that each contains a single asteroid. Our results differ from those of Martikainen et al. (2021) as they pointed out only two distinctive groups. The updated computations suggest that, for the M and S classes, the G 1 parameter is smaller and the G 2 parameter is larger, moving the group toward the upper left. The opposite happens to the third group as the G 1 parameter is larger and the G 2 parameter smaller than in the previous computations moving the group to the lower right. As above, we must emphasize the preliminary nature of the studies by us and Martikainen et al. (2021), due to the scarce data available. We also recall that, in producing the plot in Figure 8, we did not check the possibility that, according to more modern taxonomic classifications, some of the objects we have analyzed could belong to a different class.

Gaussian-Sphere Asteroid
The Gaussian-random-sphere asteroid simulation is carried out using the observational cadence of asteroid (26) Proserpina both for the Gaia DR2 and ground-based observations (Durech et al., 2010;Gaia Collaboration et al., 2018;Muinonen et al., 2020). It involves the shape and rotation model in Torppa and Muinonen (2005) and the so-called lunar mare surface scattering model  of the form in Eq. 5. The asteroid is assumed to have a photometric phase function coinciding with the H, G 1 , G 2 fit to the lunar phase curve (Muinonen et al., 2010). Basic information about the simulated observations is offered in Table 4. Zero-mean Gaussian random errors with a standard deviation of 0.01 mag were added to the numerically integrated magnitudes.
There are 29 dense lightcurves with altogether 6,664 data points and a single sparse lightcurve with 39 data points simulated for the GS-asteroid ( Table 4). The total number of points is thus 6,703, spanning more than 51 years. The data set covers a phase-angle range from 1.77°to 21.62°, altogether 19.85°. Note that the single sparse lightcurve covers the phase angle range of~14-22°and does not constrain the opposition effect. Table 5 describes the properties of the individual lightcurves in terms of their time spans, mean sampling time intervals, and number of data points. Additionally, cubic spline fits dictated by the Bayesian Information Criterion (Liddle, 2007) for the numbers of evenly spaced spline nodes are presented for each lightcurve. The rms-values of the fits are in agreement with the 0.01-mag standard deviation used to simulate the errors.
Lightcurve inversion for the GS-asteroid data was carried out using Lommel-Seeliger ellipsoids and general convex shapes for two principal cases: one where all of the lightcurves were assumed to represent relative photometry and the other where all were assumed to represent absolute photometry. In addition, two error models were studied for each of the principal cases: the first model assigned a unit weight for each ground-based lightcurve, whereas the second model differed for the relative and absolute photometry as described in Section 3. Consequently, the analysis consisted of 8 different case studies.
In the case of convex inversion, after in-depth studies using least squares and MCMC, the error model for relative photometry was chosen to be the conventional one presented in Muinonen et al. (2020), whereas the error model for absolute photometry was taken to be the one with a unit weight for each dense ground-based lightcurve. The latter, conservative error model was required to account for, primarily, the systematic errors arising from the forward modeling of the inverse problem. Likewise, in ellipsoid inversion concerning both relative and absolute photometry, a unit weight was assigned for each dense ground-based lightcurve.
Least-squares fitting and MCMC sampling were carried out for all cases. In final MCMC inversion, 30,000 samples were computed using Lommel-Seeliger ellipsoids and 30,000 samples were computed using general convex shapes. Altogether, of the order of a month of sequential CPU time was required. Tables 6, 7 describe, for convex inversion, the details of the initial and final models of uncertainties in the cases of relative photometry and absolute photometry, respectively. The tables include the effective number of observations for each lightcurve, to be compared to the actual number of observations ( Table 5). The final rms-values on the magnitude scale indicate that, even with convex inversion, the goodness of fit cannot reach the standard deviation used to simulate the errors. This is due to the forward modeling incorporated in the inversion, in particular, due to the restrictive linear-exponential phase function and the pure Lommel-Seeliger scattering model. Indeed, we can see that the fits are better for the case of relative photometry, due to the less stringent requirements on the phase function. Figure 9 shows examples of simulated lightcurves and the modeling results from MCMC sampling. The overall match between the simulated observations and modeling is quite satisfactory, but there are some systematics present, in particular, for the dense lightcurve #12 in the case of relative photometry. The conservative error model for absolute photometry shows up as large uncertainty envelopes using all of the MCMC samples, the standard deviations being significantly smaller. Table 8 gives the true rotational parameters used in the generation of the simulated data. Thereafter, it includes the results from least-squares and MCMC inversion, using ellipsoids and general convex shapes, for the rotational and phase function parameters in the cases of relative and absolute simulated data. Note that the angular width of the opposition effect has been assumed to be α 0 = 3°in the cases of absolute data. The reason is that the present data set has been insufficient to reliably predict the angular width. With the fixed value of α 0 = 3°, we have β S (20°) ≈ β 0 + m 0 × 0.02431/ rad from Eq. 9. For the absolute photometry, the mean values from MCMC sampling are β S (20°) = 1.725 mag/rad using ellipsoid inversion and β S (20°) = 1.670 mag/rad using convex inversion. Both values remain reasonably close to the corresponding values of β 0 . Based on the H, G 1 , G 2 fit to the lunar phase curve (Muinonen et al., 2010), the lunar photometric slope at 20°phase angle is 1.659 mag/rad. We obtain β S (20°) = 2.06 ± 0.23, 1.825 ± 0.063, 1.725 ± 0.085, and 1.670 ± 0.041 in the cases of ellipsoid and convex inversion with relative data and ellipsoid and convex inversion with absolute data, respectively. With absolute photometry, both convex and ellipsoid inversion yield accurate results. For relative photometry, the retrieved slopes differ more substantially from the lunar slope. Based on the uncertainties, the lunar slope is, however, encompassed with reasonable probability. Table 9 gives the correlation coefficient for the rotational and phase function parameters, as computed from the MCMC samples. Considerable and realistic correlations are found for pole longitude and latitude and for m 0 and β 0 in the case of absolute photometry. In ellipsoid inversion, the rotational phase correlates with the pole longitude.
It is clear that, with ellipsoids, there can be statistically significant differences from the true values of the parameters. However, with the current error modeling, ellipsoid inversion works in an acceptable way both for relative and absolute photometry. For more solid conclusions, additional examples must be studied. Convex inversion works well for both relative and absolute data.
Figures 10, 11 depict MCMC samples from convex inversion for the rotational and photometric parameters, and two spherical harmonics coefficients of the Gaussian surface density. Whereas the true rotation period is close to the centers of the distributions for relative photometry, there are apparent biases for absolute data. Since we have simulated the lightcurve data by using an advanced particulate medium scattering model, there are strictly no true values for the photometric parameters to compare with.

CONCLUSION
We have provided models of observational uncertainties for four classes of lightcurves, that is, for dense relative, sparse relative, dense absolute, and sparse absolute lightcurves. We have studied the models using an asteroid simulation based on the Gaussianrandom-sphere shape model, making use of linear or linearexponential phase function models for phase angles below 50°on the magnitude scale.
We have extended the Bayesian lightcurve inversion methods to allow for full retrieval of phase functions with their uncertainties and applied the methods to a set of some 500 asteroids with both ground-based and Gaia photometry. In comparison to earlier work, certain differences are seen in the mapping from apparent phase curves to intrinsic phase functions. These differences may derive, to an extent, from the scarce amount of Gaia DR2 photometry, and we expect to resolve the issues with the forthcoming Gaia DR3 photometry.
Finally, an asteroid's spectrum from the ultraviolet (UV) through the visible (Vis) to the near-infrared (NIR) regime is known to vary with phase angle, evidently due to the fact that the phase function depends on the intrinsic brightness of the asteroid at a given wavelength. In the present work, we have confirmed that the shape and rotational parameters of the asteroid and the geometry of illumination and observation have an effect on the asteroid's photometric phase function. The effect is different for different geometric albedos. Thus, the asteroid's UV-Vis-NIR spectrum can depend on the shape and rotation parameters and vary as a function of the geometry of illumination and observation. It can be particularly enlightening to consider laboratory spectrometric measurements, such as those by Cloutis et al. (2012), extended to analogs of complete asteroids. The phase curve effects on the UV-Vis-NIR spectrometry and spectropolarimetry of asteroids remain as an open problem for future theoretical, observational, and experimental studies.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
KM has been responsible for the article concept, programming, simulations, analysis, and writing the article. EU, together with JM, have contributed many of the figures. EU, JM, and AC have participated in writing the article and JM, AP, and XW have Frontiers in Astronomy and Space Sciences | www.frontiersin.org helped in the simulations and analysis. All have offered comments on the article.