Asteroid lightcurve inversion with Bayesian inference

Context. We assess statistical inversion of asteroid rotation periods, pole orientations, shapes, and phase curve parameters from photometric lightcurve observations, here sparse data from the ESA Gaia space mission (Data Release 2) or dense and sparse data from ground-based observing programs. Aims. Assuming general convex shapes, we develop inverse methods for characterizing the Bayesian a posteriori probability density of the parameters (unknowns). We consider both random and systematic uncertainties (errors) in the observations, and assign weights to the observations with the help of Bayesian a priori probability densities. Methods. For general convex shapes comprising large numbers of parameters, we developed a Markov-chain Monte Carlo sampler (MCMC) with a novel proposal probability density function based on the simulation of virtual observations giving rise to virtual least-squares solutions. We utilized these least-squares solutions to construct a proposal probability density for MCMC sampling. For inverse methods involving triaxial ellipsoids, we update the uncertainty model for the observations. Results. We demonstrate the utilization of the inverse methods for three asteroids with Gaia photometry from Data Release 2: (21) Lutetia, (26) Proserpina, and (585) Bilkis. First, we validated the convex inverse methods using the combined ground-based and Gaia data for Lutetia, arriving at rotation and shape models in agreement with those derived with the help of Rosetta space mission data. Second, we applied the convex inverse methods to Proserpina and Bilkis, illustrating the potential of the Gaia photometry for setting constraints on asteroid light scattering as a function of the phase angle (the Sun-object-observer angle). Third, with the help of triaxial ellipsoid inversion as applied to Gaia photometry only, we provide additional proof that the absolute Gaia photometry alone can yield meaningful photometric slope parameters. Fourth, for (585) Bilkis, we report, with 1-σ uncertainties, a refined rotation period of (8.5750559± 0.0000026) h, pole longitude of 320.6◦ ± 1.2◦, pole latitude of −25.6◦ ± 1.7◦, and the first shape model and its uncertainties from convex inversion. Conclusions. We conclude that the inverse methods provide realistic uncertainty estimators for the lightcurve inversion problem and that the Gaia photometry can provide an asteroid taxonomy based on the phase curves.


Introduction
Asteroids are irregularly shaped Solar System bodies rotating, typically, about their axes of maximum inertia. Their surfaces are assumed to be covered by regolith, that is, particulate material in the size scales of microns to meters. Asteroids offer information about the origin and evolution of the Solar System, may provide valuable space resources, and represent an impact hazard for life on the Earth.
The most plentiful source of data for the physical characterization of asteroids is photometry: the measurement of the diskintegrated brightness of the asteroid. An asteroid's lightcurve, that is, its observed brightness as a function of time, depends on the shape and spin state of the asteroid, as well as the scattering properties of its surface. It follows that these properties can be estimated from observations, to the extent allowed by a given data set. The phase curve of the asteroid results from multiple observations during an apparition. The phase curve describes the variation of the disk-integrated brightness with the solar phase angle, the angle between the Sun and the observer as seen from the object.
The theory of asteroid lightcurve inversion was presented by Kaasalainen et al. (1992a,b) and Lamberg (1993), demonstrating the estimation of the Gaussian surface density (inverse of total curvature) from lightcurve data using lightcurves observed in varying Sun-asteroid-observer geometries. Preliminary numerical inverse methods for retrieving convex shapes were initiated simultaneously, accompanied by a successful application to (951) Gaspra by Barucci et al. (1992). Robust inverse methods were established a decade later by  and . Russell (1906), considering lightcurves of the opposition geometry only, had underscored the illposedness of the inverse problem. The progress almost a century later resulted in no less than a paradigm change.
The convex methods concern shape inversion both through the estimation of the Gaussian surface density and by directly solving for the facet areas of a triangulated, polyhedral shape. The surface density approach was found to lead to a more stable A138, page 1 of 19 inversion result, while the facet approach was shown to describe sharp features in lightcurves better. Convex shape models of highly non-convex test objects were found to agree well with the convex hulls of the original shapes . Limiting the shape modeling to convex shapes leads to a much more stable inverse problem, compared to the general nonconvex case. Most of the shape models produced so far through lightcurve inversion are convex. Ďurech et al. (2010) present a database of shape models, both convex and non-convex, derived primarily from lightcurves. Some of the shape models obtained are refined by other observational data, such as radar range and Doppler data, infrared interferometry, and stellar occultations , which provide additional constraints, especially on the size of the asteroid. The use of these additional data sources together with photometry is described inĎurech et al. (2015). Recently, disk-resolved observations with adaptive optics (Vernazza et al. 2020) have become a new data source for the largest asteroids. More information on recent advances in nonconvex inversion is available in works by Viikinkoski et al. (2015) and Bartczak & Dudziński (2018).
Due to the fact that the illuminated areas of ideal triaxial ellipsoids observed in any illumination condition can be computed using analytical formulas, this particular shape has been assumed to be the best choice to quickly process large data sets of sparse photometric observations. The most important current application is the inversion of sparse photometric observations of asteroids carried out by the Gaia mission. On one hand, the main limitation of this method is clearly the fact of using a shape model that is, a priori, not completely realistic. On the other hand, a triaxial ellipsoid shape is a reasonably flexible choice, and is adequate to fit a large variety of complex shapes well, as also confirmed by Torppa et al. (2008) and Carbognani et al. (2012).
The effectiveness of an approach to the inversion problem based on a triaxial ellipsoid shape model and the development of a genetic algorithm to the simultaneous determination of the spin period, the pole coordinates, the b/a and c/a axial ratios for the triaxial shape (ellipsoid semiaxes a ≥ b ≥ c), and also a supposedly linear phase angle -mean magnitude relation, has been proven in a wide variety of numerical experiments and analysis of real data (Kaasalainen &Ďurech 2007;Cellino et al. 2009;Carbognani et al. 2012;Santana-Ros et al. 2015). The paper by Santana-Ros et al. (2015) includes analysis of the inversion of a large set of simulations (more than 10 000 simulated asteroids) of sparse data corresponding to the sky survey by Gaia. The simulations were carried out considering both ideal triaxial ellipsoid (the "easy" case) and complex shapes, combined with random photometric uncertainties larger than those predicted for most real asteroids observed by Gaia. The results showed that the inversion is reliable even in the most difficult cases (irregular shapes and large photometric uncertainties), although a 180degree ambiguity in the determination of the ecliptic longitude of the pole is inherent to the chosen triaxial ellipsoid shape (Santana-Ros et al. 2015). The 180-degree ambiguity affects the inversion, in particular, for asteroids in low-inclination orbits.
The paper by Cellino et al. (2009) included the first application of the same inversion algorithm to the data set of photometric measurements of asteroids obtained by the HIPPARCOS satellite. In spite of the HIPPARCOS data being of low quality in terms of both the number of measurements per object and the photometric accuracy of the data, correct determinations of the spin period and decent determination of the pole coordinates (by a comparison with available ground-based data based on full lightcurves) could be obtained for little less than 50% of the cases.
A limitation of the abovementioned results has been the assumption of a purely geometric relation between incident and scattered light. More realistic simulations using a Lommel-Seeliger surface reflectance model computed in the case of ideal triaxial ellipsoid shape have been developed (Muinonen & Lumme 2015). Cellino et al. (2015) applied the triaxial ellipsoid shape model and the genetic inversion algorithm to a number of simulated data of complex Gaussian shapes and Lommel-Seeliger scattering law. Again, the obtained results of photometric inversion have been found to be quite accurate.  complemented these methods with a full Markov-chain Monte Carlo treatment (MCMC) of the lightcurve inversion problem using Lommel-Seeliger ellipsoids. The genetic algorithm with the Lommel-Seeliger scattering ellipsoid  has been implemented in the CU4 Solar-System-Object Long-Term pipeline (SSO-LT) of the Gaia Data Processing and Analysis Consortium (DPAC).
We note that initial assessment of the Gaia Data Release 2 (GDR2) photometric data on asteroids (Gaia Collaboration 2018a,b) has been offered byĎurech & Hanuš (2018), who provided 173 asteroid models, of which 129 were new, by combining GDR2 photometry and ground-based lightcurves. Considering GDR2 data alone, they concluded that rotation periods could be derived even in the case of small numbers of observations (<20). However, only for more than 30 observations, the periods were likely to be correct. Similar results were obtained more recently by Cellino et al. (2019), who completed their first study of the GDR2 data, proving the satisfactory performance of the Lommel-Seeliger ellipsoids and showing also the first encouraging results using the so-called cellinoid shape models (Cellino et al. 1989).
The objective of the classical statistical inference is to gain knowledge of the unknown parameter values of a model on the basis of the observed data. The virtual observation method put forward by Wang et al. (2015a,b) for asteroid lightcurve inversion is based on the classical inference. Random samples of parameters are derived from virtual observations by least-squares minimization. The virtual observations in each night are generated randomly by a Monte Carlo method from the original data by adding Gaussian random noise with an assumed variance, which is the uncertainty in the photometric observations on a single night. Based on the distribution of the parameters, statistical descriptors such as the mean and variance can be estimated.
If the means of the distributions are close to the parameters estimated from the original data, unbiased estimates for the parameters are obtained. In fact, the parameter uncertainties derive from two parts: the uncertainties in the data and the uncertainties in the model. As the virtual observations are generated according to the quality of the observations, the effect of the uncertainties in the data on the estimate of the parameters can be investigated by the variances of the distributions, resulting from the propagation of the uncertainties through a given model. Here, in Bayesian statistical inference, we take the distribution of the virtual solutions for the parameters as the proposal distribution in MCMC.
Our current study is organized as follows. In Sect. 2, we present forward modeling for the asteroid, consisting of the rotational parameters, the triaxial ellipsoid shape, the convex shape expressed in terms of the Gaussian surface density, and the photometric phase-curve parameters. We formulate the statistical inversion problem with a detailed assessment of both random and systematic uncertainties. In Sect. 3, we describe the MCMC sampler for the characterization of the a posteriori probability densities. We offer a concise review of the conventional convex inversion method, including the solution to the Minkowski problem of retrieving the shape from the Gaussian surface density. In Sect. 4, we apply the methods to the asteroids (21) Lutetia, (26) Proserpina, and (585) Bilkis. We offer our conclusions and future prospects in Sect. 5.

Asteroid model
We consider an asteroid in principal-axis rotation about its axis of maximum inertia, and denote the 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 shape, we incorporate either triaxial ellipsoidal shapes or general convex shapes.
The diffuse reflection coefficient R of an asteroid's surface element relates the incident flux density πF 0 (including π in accordance with a common definition) and the emergent intensity I as where ι and are the angles of incidence and emergence as measured from the outward normal vector of the element, and φ 0 and φ are the corresponding azimuthal angles. It is customary to measure φ so that the backscattering direction (the direction for the source of light) is with φ = 0 • . With the common assumption of a geometrically isotropic surface, specifying φ 0 is unnecessary. The Lommel-Seeliger reflection coefficient (subscript LS) is (e.g., Lumme & Bowell 1981;Wilkman et al. 2015), whereω is the single-scattering albedo (proportion of incident light scattered, 0 ≤ω ≤ 1), P 11 is the single-scattering phase function (from the 4 × 4 scattering phase matrix), and α is the phase angle. We note that α derives unequivocally from ι, , and φ. R LS , the first-order multiple-scattering approximation of the radiative transfer equation for a semi-infinite plane-parallel medium (e.g., Chandrasekhar 1960), is applicable to weakly scattering media: this means that the intensity terms [ω k ], k ≥ 2 are assumed negligible. When omitting polarization effects, the scattering phase function P 11 provides the angular distribution of scattered light in an individual interaction and is normalized so that The recent works by, for example, Muinonen et al. (2018) and Väisänen et al. (2019) provide advanced multiple-scattering models for planetary regoliths. The disk-integrated brightness L equals the surface integral where A + stands for the part of the surface that is both illuminated by the light source and visible to the observer. For a nonspherical asteroid, L typically depends strongly on the orientation of the asteroid. As a consequence, the phase curve is generally not a function. For a spherical asteroid with diameter D, where Ω + stands for the part of the unit sphere both illuminated and visible. With the Lommel-Seeliger reflection coefficient, the calculation of L in Eq. (5) can be carried out analytically, where we have defined a phase function Φ LS (α) normalized to unity at zero phase angle (as are all functions Φ from now on).
The geometric albedo p is, by definition, the ratio of the disk-integrated brightness of the asteroid and the disk-integrated brightness of a normally illuminated Lambertian disk, both with the same effective diameter D, in the exact backscattering direction α = 0 • : The geometric albedo can obtain values larger than unity (p ≥ 0). For example, for a plane mirror oriented towards the light source, the geometric albedo approaches infinity. Consider the single-scattering albedoω and phase function P 11 (α) in the Lommel-Seeliger reflection coefficient R LS (Eq. (2)). On one hand, the disk-integrated brightness phase curve of an asteroid in a given apparition can be described by the H, G 1 , G 2 phase function Φ HG 1 G 2 Penttilä et al. 2016). On the other hand, the reflection coefficient R LS in Eq. (2) includes the (µ + µ 0 ) −1 term that produces the phase function Φ LS given in Eq. (6). Thus, modeling P 11 (α) so that the phase curve from Eq. (2) for a spherical asteroid is precisely that of the H, G 1 , G 2 phase function is an attractive alternative. Furthermore, the geometric albedo p of a Lommel-Seeliger sphere is, by using L(0) from Eq. (6) in Eq. (7), Thus, it is reasonable to assume that the model resulting in can serve well in asteroid phase-curve analyses. A normalization to the Lommel-Seeliger phase curve of a spherical asteroid is convenient, especially for inverse methods, as it is independent of the actual shape being derived. Furthermore, the H,G 1 ,G 2 phase function allows one to obtain estimates for the phase integral and the Bond albedo ). The V-band magnitude of a spherical asteroid is where H is the absolute magnitude, r and ∆ are the heliocentric and topocentric distances (in au), D is the diameter (in km), and lg denotes the 10-based logarithm. It is known from the observations of asteroids that the phase dependence in the magnitude scale is approximately linear for 20 • ≤ α ≤ 50 • . In order to characterize the steepness of the phase curve at moderate phase angles, we thus introduce the slope parameter β S (in mag/• or mag/rad): The linear dependence can be ascertained, for example, from the phase curves in Fig. 1 (see also Table 1). The phase-angle range of validity is conservative: for certain G 1 , G 2 parameters, the linear dependence can start from smaller phase angles or can extend to larger phase angles. Consider the triaxial ellipsoidal shape model for an asteroid with semiaxes a, b, and c, and denote C = diag(a −2 , b −2 , c −2 ) (C is a diagonal 3 × 3 matrix). Let e and e ⊕ be the unit vectors in the light-time-corrected directions of the Sun and the observer, respectively, as seen in the principal axes reference frame of the ellipsoid. The disk-integrated brightness of the Lommel-Seeliger ellipsoid is given by where Notes. Also shown are the slope parameters β S at α = 20 • . See Fig. 1.
and ln is the natural logarithm. Equation (13) is an alternative expression for the one in Muinonen & Lumme (2015). The general convex shape of an asteroid is modeled with the Gaussian surface density G, requiring it to be positive everywhere. In the spherical coordinate system (θ, ϕ) describing the directions of the outward unit normal vectors, G is expressed as where the Y lm 's are the orthogonal spherical harmonics and the complex-valued coefficients s lm ensure a real-valued surface density (the superscript * refers to complex conjugation): s l,−m = (−1) m s * lm , Im(s l0 ) = 0, l = 0, . . . , l max , m = − l, . . . , 0, . . . , l.
The series starts from a minimum degree l min (usually l min = 0) and is truncated at a maximum degree l max (typically l max ≤ 8, see Sect. 4.2 below), and there are altogether (l max + 1) 2 − l 2 min independent real or imaginary parts of the coefficients s lm .
A138, page 4 of 19 When utilizing real-valued spherical harmonics, we utilize wellbehaving normalized associated Legendre functions. This entails multiplication of the normal associated Legendre functions by √ (2l + 1)/2 √ (l − m)!/(l + m)!. The Gaussian surface density gives the relative proportion of surface normal vectors in the given directions specified by the normal coordinates. The convex shape can be derived from the Gaussian surface density unambiguously via the solution of the Minkowski problem (Kaasalainen et al. 1992a,b;Lamberg 1993), see Sect. 3.1.
We have chosen to use the pure Lommel-Seeliger reflection coefficient (Eq. (2)) instead of the frequently-used combination of the Lommel-Seeliger and Lambert coefficients (e.g., . First, we have wanted to maintain consistency with the methods based on the ellipsoids and general convex shapes. Second, in terms of physics, it continues to be unclear whether the Lambert component would significantly improve the reflection coefficient. The choice affects the shape retrieved from the lightcurves (Sect. 4.2.1). Third, for simplicity, we assume that the surfaces of individual asteroids are homogeneous in terms of the scattering properties, implying that the single-scattering albedo and phase function in Eq. (2) do not vary across the surface.

Inverse problem
In the case of general convex shape models, we denote the rotation, size, shape, and scattering parameters of an asteroid at a given epoch t 0 by the vector P = (P, λ, β, ϕ 0 , s 00 , . . ., s l max l max , p, G 1 , G 2 , D) T , where the parameters are, respectively, the rotation period, ecliptic pole longitude, ecliptic pole latitude, rotational phase at t 0 , shape parameters s 00 , . . . , s l max l max , geometric albedo, two parameters of the H, G 1 , G 2 phase function, and the physical scale parameter (e.g., diameter for a spherical asteroid). Currently, we omit the geometric albedo p and the diameter D from the set of parameters. The total number of parameters equals N P = (l max + 1) 2 − l 2 min + 6.
Here the rotational phase ϕ 0 is strictly redundant, as the convex shape model will accommodate for any change in ϕ 0 . We will nevertheless keep it in the list of parameters to maintain consistency with the triaxial ellipsoid model (see below). In the case of triaxial ellipsoid shape models, the parameters of an asteroid reduce to the vector P = (P, λ, β, ϕ 0 , a, b, c, p, G 1 , G 2 , D) T . Furthermore, we set a = 1, consider b and c to be ratioed against a, and take D to represent the physical diameter of an equal-volume spherical asteroid. Again, p and D are presently excluded from the parameter set so there are maximum 8 parameters remaining in the inverse problem involving triaxial ellipsoid shapes.
For sparse photometry spanning varying illumination and observation geometries, we may utilize the two-parameter H, G 12 phase function , including G 12 among the parameters. Furthermore, we set l min = 0. In this case, for general convex shapes, and our parameters are P = (P, λ, β, ϕ 0 , s 00 , . . . , s l max l max , G 12 ) T . For triaxial ellipsoid shapes, we have P = (P, λ, β, ϕ 0 , b, c, G 12 ) T . We may also utilize the single-parameter phase function provided by Penttilä et al. (2016): in this case, the G 1 , G 2 parameters are fixed based on five predefined sets of values (Table 1).
Alternatively, we may use the linear phase dependence of Eq. (12) and replace G 12 with β S . The latter is the choice adopted in the inversion of sparse photometric data by Cellino et al. (2009) and Santana-Ros et al. (2015). Let N and K be the number of observations and number of lightcurves, respectively. Then, where N k specifies the number of observations in lightcurve k (k = 1, . . . , K). Let denote, respectively, the observed and computed magnitudes. It is assumed that there are only vanishingly small uncertainties in the times of the observations and that all times have been corrected for the light time.
In our standard approach of statistical inversion, the observation equation links together the observed and computed magnitudes for the given parameters P: Here and υ stand for two kinds of uncertainties. First, represents the uncertainty that can be assumed random from one observation to another within a lightcurve. The uncertainty can have correlations among the observations of a given lightcurve but cannot have correlations among the observations from two different lightcurves. Second, υ represents the uncertainty that can be assumed random from one photometric lightcurve to another but that is equal for all observations within a given lightcurve. The uncertainty υ is thus systematic within each lightcurve. We assume that the probability densities for and υ, that is, p and p υ , are Gaussian with zero means and covariance matrices Λ and Λ υ . First, the Gaussian hypothesis for p is supported by the phase-curve studies in Penttilä et al. (2016). Second, for p υ , the actual choice plays no role in the final analysis, and we may adopt the Gaussian hypothesis for convenience. It then follows that the probability density for the combined uncertainty p +υ is also Gaussian with zero means and covariance matrix The covariance matrices Λ and Λ υ are block-diagonal. As for Λ , it may or may not include off-diagonal elements within a given block, describing correlations among the uncertainties of observations of a given lightcurve. As for Λ υ , it incorporates correlations approaching unity for all the off-diagonal matrix elements within a given block. Let p p be the a posteriori probability density function (p.d.f.) for the parameters. Within the Bayesian framework (cf. Muinonen & Bowell 1993), p p is proportional to the a priori and observational uncertainty p.d.f.s p pr and p +υ , the latter being evaluated for the "Observed-Computed" (O-C) residual magnitudes ∆M(P), p p (P) ∝ p pr (P)p +υ (∆M(P)), A138, page 5 of 19 A&A 642, A138 (2020) Even though p +υ is here related to the observations, it can also describe the uncertainties deriving from the possible 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.
It is enlightening to study the limits of either large or negligible uncertainties υ, the uncertainties that are systematic within a given lightcurve. Let us introduce, for each lightcurve, the mean observed and computed magnitudes M obs,k0 and M k0 (P) and their difference ∆M k0 (P): Furthermore, define the magnitude differences m obs,k j and m k j (P) with the help of Eq. (26): Introduce then the relative brightnesses obs,k j and k j (P): In the case of large uncertainties υ, that is, in the case of what can be called relative photometry, Recalling that Λ υ includes high correlations among the uncertainties for the observations of a given lightcurve, We do not know the actual value of σ υ,k , the standard deviation for the uncertainties υ in the lightcurve k. We note that the values of the residuals in the brackets are approximately equal and large in the absolute sense. In order to make use of Eq. (30), we utilize the relative and mean magnitudes in Eq. (27): In the present case of large systematic uncertainties, the absolute magnitude information is lost. Subsequently, in the last two steps of Eq. (31), we re-scale the observed and theoretical mean magnitudes to coincide with one another, resulting in ∆M k0 (P) → 0.
In other words, we can introduce a systematic correction (or, formally, an additional parameter) for each lightcurve to exactly result in ∆M k0 (P) = 0. Furthermore, we have assumed an equal standard deviation σ υ for all lightcurves. Equation (31) shows that the least-squares optimization for relative photometry can conveniently be carried out by giving an equal weight for each lightcurve (assuming σ υ,k ≡ σ υ for all k). However, the equation also shows that the probabilistic treatment runs into difficulties, since, with the missing information about σ υ , there is no meaningful weighting scheme for χ 2 (P) to be utilized in establishing the a posteriori p.d.f. in Eq. (24). To summarize the case of predominating uncertainties υ, what has been described above provides the justification for the known thumb rule that, for large numbers of lightcurves taken relative to one another, the convex least-squares optimization converges when each lightcurve obtains an equal weight.
In the case of negligible uncertainties υ, that is, in the case of what can be termed absolute photometry, For a diagonal matrix Λ , If the uncertainties are small, the O-C magnitude difference can be linearized into a difference in the corresponding brightnesses L obs,k j and L k j (P): where 2.5 lg e ≈ 1.0857. Thus, the inverse problem is linearized and the uniqueness properties of convex inversion are conserved (Kaasalainen et al. 1992b;Lamberg 1993).
With the help of Eqs. (26)-(28), and (34), Eq. (33) leads to the expressions where σ ,k now describes the uncertainty of the observations in lightcurve k and is taken to account for the factor of 2.5 lg e. In Eq. (35), the χ 2 -value is computed using the differences in the observed and computed relative brightnesses, relative to the observed relative brightnesses. This differs from the earlier models of the uncertainties, for example, in  and Torppa et al. (2018), and derives from defining the inverse problem for magnitudes (Eq. (21)) instead of brightnesses.

Model for uncertainties
For the uncertainty modeling, Eq. (35) provides a starting point in the cases of both predominating and negligible systematic uncertainties. The modeling is required to account for the following points. First, there are two different classes of lightcurves entering the inverse problem, that is, lightcurves densely and sparsely sampled in time. Second, there are different numbers of observations in the lightcurves, the difference sometimes amounting to orders of magnitude. Third, the time intervals spanned by the lightcurves vary for both classes of lightcurves. For dense lightcurves, the intervals can vary from a fraction of an asteroid's rotational period to several periods. For sparse lightcurves, the intervals may vary from days, weeks, and months to several decades. Fourth, the sampling rates of the lightcurves in time are different for each lightcurve. Fifth, due to the presence of ∆M k0 (P) in Eq. (35), we may account for it in the absolute photometry or set it to zero in relative photometry, leading to a meaningful probabilistic analysis in both cases. Sixth, the uncertainty modeling must result in consistent treatments of the inverse problem in cases, where, for one reason or another, different sets of lightcurve data are utilized. In Bayesian inference, the regularization of the inverse problem and the weighting of the observations must be carried out with the help of the a priori p.d.f.s. Considering Eqs. (35) and (24), one may arrive at different σ ,k -values and different a posteriori p.d.f.s with an a priori p.d.f. that utilizes the χ 2 of Eq. (35) in its exponent with an opposite sign and an overall scaling factor. Additionally, the a priori p.d.f. must then utilize a part that enforces it to vanish far away from the regime with meaningful probability mass. For brevity, we choose to omit the detailed mathematical formulation of such an a priori p.d.f. here: it is sufficient to state that the p.d.f. exists and can be formulated if necessary.
For dense lightcurves, we devise the model for the observational uncertainties and, consequently, the observational weights based on the idea that the lightcurve with the sparsest sampling obtains a basic weight of unity and that the lightcurves with denser sampling are considered as consisting of a number of independent lightcurves corresponding to the sparsest sampling. In other words, the lightcurve with the sparsest sampling is the basic individual observation in the inverse problem. Consequently, we set where σ 0,k denotes the rms-value as computed using the square form in Eq. (35) 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 and k marking the lightcurve with the lowest sampling rate in time and T k standing for the time span of the lightcurve k. Equation (36) can be interpreted in the following way. The lightcurvẽ k contributes a χ 2 -value of unity in the inverse problem, whereas each other lightcurve k is considered to split into approximately independent lightcurves. For sparse lightcurves, we devise the model for the uncertainties and weights by identifying the lightcurve with the smallest number of observations and assigning the same maximum weight to the remaining lightcurves with larger numbers of observations. Thus, we set where, again, σ 0,k denotes the rms-value as computed using the square form in Eq. (35) and σ pr,k denotes the a priori threshold value for the uncertainty in the lightcurve k. Now the lightcurvẽ k corresponds to the one with the minimum number of observations and contributes a χ 2 -value equal to N˜k (when σ 0,k ≥ σ pr,k ), whereas each other lightcurve k is scaled to result in the same χ 2 -value at maximum. The final a posteriori p.d.f. is the product of the p.d.f.s for the two classes. The two p.d.f.s can be considered as a priori p.d.f.s for each other. The final χ 2 takes the form where the primes and double primes denote the dense and sparse lightcurves, respectively. Based on the principles above, it is possible to develop a consistent uncertainty model also for the case of dense, absolute lightcurves. After finding the sparsest sampling in the A138, page 7 of 19 A&A 642, A138 (2020) lightcurves, one can introduce revised weights for the individual observations in each absolute lightcurve with a procedure analogous to the one introduced for the individual relative lightcurves above.
We regularize the inversion for convex shapes with χ 2 R (P) that depends on the shape parameters s 00 , . . . , s l max l max included in P: where A j (P) and n(P) denote the area and unit outward normal vector of the triangular element j = 1, . . . , N . In order for A j (P) and n(P) to describe a convex shape, the χ 2 R (P)-value of Eq. (41) should vanish. Here the triangular discretization needs to have a resolution high enough to warrant an accurate computation of the disk-integrated brightnesses. In the spirit of the Bayesian inference, σ reg describes the standard deviation for the regularization term. We return to specifying the actual value of σ reg in Sect. 4.2.
Finally, when considering sparse observational data sets, like the GDR2 photometric data sets, the convex inversion is additionally regularized by constraining the spectral power in the spherical harmonics expansion for the Gaussian surface density. After experimenting with the present asteroid examples, we require that the spectral power in each degree of spherical harmonics must not exceed a value of 0.5. In the present examples, this does not substantially limit the resulting spectra of shapes, but prevents the MCMC sampling from wandering into unrealistic parameter regimes.

Virtual-observation MCMC sampler
Markov-chain Monte Carlo methods (cf., O'Hagan & Forster 2004) provide the means for sampling unnormalized p.d.f.s such as that in Eq. (24). We utilize the Metropolis-Hastings algorithm that is based on the computation of the ratio a r : a r = p p (P )p t (P j ; P ) p p (P j )p t (P ; P j ) .
Here P j and P denote the current and proposed parameters in a Markov chain, respectively, and p t (P ; P j ) is the proposal p.d.f. from P j to P (subscript "t" stands for transition). The proposed parameters P are accepted or rejected with the help of a uniform random deviate y ∈ [0, 1]: that is, the proposed parameters are accepted with the probability of min(1, a r ). After a number of transitions in the so-called burnin phase, the Markov chain, in the case of success, converges to sampling the target p.d.f. p p . For convergence monitoring, there are various diagnostics tools available (see, e.g., Oszkiewicz et al. 2010).
In the virtual-observation MCMC methods, earlier studied for asteroid orbital inversion (Muinonen et al. 2012) and lightcurve inversion with Lommel-Seeliger ellipsoids , the potentially complex solution phase space is first characterized statistically to arrive at a relevant proposal p.d.f. A set of virtual observations M v is generated from the original observations M obs through the addition of Gaussian random uncertainties v with zero means and covariance matrix Λ v , Thereafter, virtual least-squares parameters P v are obtained by minimizing where the covariance matrix Λ + Λ v covers both real and virtual uncertainties, Λ representing the model of uncertainties in Sect. 2.2.
As v is a vector of Gaussian random variables, the parameters P v are random variables, too. The p.d.f. for P v is formally the N-dimensional integral where p(M v ) is the Gaussian p.d.f. for the virtual observations and δ is Dirac's Delta function.
A random-walk MCMC-sampler is obtained with the help of the symmetric proposal p.d.f., defined as the 2(N + N p )dimensional integral where the integration over N p parameters results in a (2N + N p )dimensional convolution integral It is noteworthy that, due to the symmetry of the proposal p.d.f., these integrals need not be numerically evaluated. The procedure above meets the mathematical requirements for a symmetric proposal p.d.f. in a Metropolis-Hastings algorithm (O'Hagan & Forster 2004). However, it may be inconvenient to generate such proposals during the MCMC sampling. A straightforward means for utilizing the symmetric proposal probability density in Eq. (48) is to compute large numbers of parameters P v before the MCMC sampling (cf. Wang et al. 2015a,b). The generation of virtual observation sets L v and the derivation of P v are repeated until there are N v 1 sets of parameters P ( j) v ( j = 1, 2, 3 . . . , N v ). Parameter differences are then readily available from These differences offer a finite number of relevant transitions from one parameter set to another in the phase space. There are of the order of N 2 v differences available, that is, a number that can become extremely large for reasonable values of N v . Any pair of virtual least-squares parameters consumed for a difference should be replaced by a new pair of parameters computed as described above, but, in practice, if N v is large enough, the replacement can be unnecessary.
The ratio a r for the decision criterion is expressed in the concise form The proposals can now be computed during the MCMC sampling using the N v virtual least-squares element sets P v, j . The covariance matrix Λ v can be scaled from Λ . For dense lightcurves, it is also possible to generate the virtual observations from cubic spline fits to the lightcurves, in which case we may set Λ v = Λ .
The virtual least-squares solutions can be utilized to compute a covariance matrix for the parameters. It is straightforward to incorporate the covariance matrix in a multivariate Gaussian proposal probability density. Combining the virtual-observation and Gaussian proposals results in a hybrid proposal probability density that will be used in the MCMC sampling below.

Virtual-observation importance sampler
In what follows, we describe a virtual-observation random-walk importance sampler, utilizing the probability density for the virtual least-squares parameters in Eq. (46). Instead of drawing random samples from the a posteriori probability density for the parameters (Eq. (24)), we can draw random samples of parameters from a uniform distribution within a given regime in χ 2 (P), where P 0 denotes the least-squares parameters and the cut-off ∆χ 2 c is set on the basis of the number of parameters. For the present inverse problem, we may define ∆χ 2 c with the help of probability thresholds for multivariate normal statistics: where P c denotes the cut-off probability value. It is recognized that Eq. (52) describes the probability mass for normal statistics only. Nevertheless, we consider it useful in approximately evaluating the statistical distance in our probability densities. For N P ≈ 50, we may choose allowing for more than 99.9999% of the probability mass to be covered if the probability density were Gaussian. For the acceptance and rejection of parameters, we have that is, transitions are always accepted within the given ∆χ 2 c regime but never across the ∆χ 2 c boundary. The proposal procedure utilized in the random-walk importance sampler is analogous to the MCMC sampler described in Sect. 3.1: only the decision criterion in Eq. (50) is replaced by the one in Eq. (54). Since the parameter phase space is uniformly sampled, the weights of the individual solutions are w j ∝ p p (P j ), j = 1, 2, 3, . . . , N.
For high-dimensional inverse problems, with the virtualobservation importance sampler, it is challenging to generate sufficient numbers of samples in the high-probability-mass regime of the phase space. In the present work, we utilize the importance sampler in the mapping of the regime of solutions. This allows us to evaluate the convergence of the MCMC sampling described in Sect. 3.2.

Least-squares optimization
In the case of general convex shapes, the Levenberg-Marquardt least-squares optimization method (Press et al. 1992) is utilized in solving for the least-squares parameter sets described in the preceding subsections. In order to guarantee stable convergence towards the least-squares solution, an additional regularization term is incorporated in the total χ 2 describing the goodness of fit. The numerical integration of the disk-integrated brightness is carried out using a triangulation of the unit sphere. The regularizing χ 2 R in Eq. (41) measures how well the Gaussian surface density with the given triangulation describes a real convex shape. There is extensive literature available for additional information on the optimization and regularization (e.g., Kaasalainen et al. 2002). In the case of triaxial ellipsoid shapes, the least-squares optimization is carried out using the Nelder-Mead downhill simplex method (Press et al. 1992). These methods are described, in detail, by .

Solution of the Minkowski problem
Generating the three-dimensional shape of an object from its Gaussian surface density G was published by Minkowski (1903) and first applied in asteroid lightcurve inversion by Kaasalainen et al. (1992a,b). Two quantities introduced in the context are the mixed volume and the support function. The support function is defined as the distance, from the origin, of the so-called tangent plane (the plane that tangents the surface at (θ, ϕ)): where n is the surface normal vector and r the radius vector. Here θ and ϕ are the spherical polar coordinates of the surface normal vector. The mixed volume of two bodies E and F is defined as According to Minkowski, for a constant volume of body E, V(E, F) reaches its minimum, when E and F refer to the same convex shape. Thus, when the Gaussian surface density and the corresponding surface normal vectors are known for F, the corresponding surface distance from the origin can be computed by minimizing the integral in Eq. (57). An analogous and computationally easier way is to maximize the volume of E, while keeping the mixed volume constant. For a discretized surface, the volume to be maximized is then where A j (l) is the area of the facet j as computed from l, the distances of the facet planes from the origin. The mixed volume for a discretized surface is where a j is the area of facet j on surface F. The constraint of keeping the mixed volume constant is fulfilled by projecting the gradients onto the constraint plane. For the present study, for a quick-look analysis, we have developed an alternative method for the derivation of the polyhedral shape from the Gaussian surface density. We start from A138, page 9 of 19 A&A 642, A138 (2020)  (21) Lutetia. We show the ground-based photometric lightcurves #41 and #3 and the Gaia lightcurve (blue circles; left, middle, and right, respectively; see Table A.1) together with the best-fit convex model lightcurves (red crosses) with 1-σ standard deviations (red bars) and uncertainty envelopes using all 30 000 MCMC sample solutions (red points). m, t, and k stand for relative magnitude, time from the first lightcurve observation, and lightcurve observation counter, respectively. a triangulated model of a spherical asteroid and treat the coordinates (x i , y i , z i ) T of its N n nodes i = 1, . . . , N n as the 3N n free parameters. The derivation of the least-squares solution for the coordinates is carried out via Monte Carlo optimization. At each iteration, the convexity of the proposed shape is verified before computing the χ 2 -deviation for the Gaussian surface density, that is, for the amounts of area in given discrete directions (see above). The optimization is carried out in a localized way for a single node at a time.

Results and discussion
The inverse methods presented in Sect. 3 are here applied to the photometric data of (21) Lutetia, (585) Bilkis, and (26) Proserpina. The example asteroids have been selected on the basis of careful consideration. Lutetia represents an asteroid for which high-resolution shape models are available from a space mission flyby, here from the flyby of the ESA Rosetta mission. In terms of the Tholen taxonomical classification, Lutetia represents a moderate-albedo, possibly M-class asteroid (X-complex asteroid in more modern classification) with 7 GDR2 observations, Bilkis is a low-albedo C/D/P-class asteroid with a substantial number of 32 GDR2 observations, and Proserpina is a moderateto-high-albedo S-class asteroid with a larger number of 39 GDR2 observations. To summarize, the selected asteroids represent different taxonomical classes with differing numbers of GDR2 observations.  .1-A.3). For the first two asteroids, we have extracted the ground-based observations from the DAMIT database (Ďurech et al. 2010). Additionally, the data set for Proserpina includes a sparse ground-based lightcurve also extracted from the DAMIT database (Ďurech Notes. "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. For Lutetia and Bilkis, the uppermost numbers refer to the ground-based observations, the ones in the middle to the GDR2 observations, and the lowermost ones refer to the combined observations. For Proserpina, the second row describes an additional sparse ground-based lightcurve. For the references to the observations, see Tables A.1  fits reveal the varying level of intrinsic uncertainties. These uncertainties correlate with the final post-fit uncertainties. Consider first the case of Lutetia in Table A.1. For the 50 dense ground-based lightcurves, the effective number of observations varies from 1 to 46, whereas the original number of observations varies from 6 to 647. On one hand, lightcurve #41 obtains the unit weight in terms of the effective number of observations, that is, the lowest weight among the lightcurves. On the other hand, lightcurve #42 obtains the highest effective number of observations, being valued at the level of approximately 46 times more than lightcurve #41. However, in terms of the cubic spline fits, the numbers of nodes dictated by the Bayesian information criterion are essentially the same, 11 nodes for lightcurve #41 and 12 nodes for lightcurve #42. Both lightcurves are accurately represented by the spline fit and turn out to be well fitted by the model, independently of the substantially differing final uncertainty standard deviations.

Initial data analysis
For (26) Proserpina (Table A.2), the number of observations in dense ground-based lightcurves varies from 25 to 403, whereas the effective number of observations varies from 1 to approximately 7. As dictated by the uncertainty model, the sparse ground-based lightcurve #30 receives the same weight as the sparse space-based lightcurve #31. However, the former receives a considerably lower weight due to the large uncertainty related to the data set. In the case of (585) Bilkis (Table A.3), the number of observations in dense ground-based lightcurves varies from 25 to 505. The effective number of observations varies from 1 to approximately 14.

Application of inverse methods
We apply the convex inversion methods to the complete data sets of all three example asteroids and to the GDR2 data sets of Proserpina and Bilkis. The resulting shape models are available from the authors. Additionally, we apply the triaxial ellipsoid inverse methods to the GDR2 data sets of all three asteroids.
For all three asteroids, in the construction of the proposal probability densities for convex inversion, we have generated virtual observations with the help of rms-values for each individual lightcurve starting from the cubic spline fits to the lightcurves. The final uncertainty models in MCMC sampling utilize individual uncertainty values for each lightcurve (see Tables A.1-A.3).
For the spherical harmonics expansion of the Gaussian surface density, we select l min = 0. The value of l max can be selected, for example, on the basis of the Bayesian information criterion. However, for securing numerical stability in the inversion, it is usually required that l max ≤ 8. We utilize l max = 6 throughout the work. As to the regularization of the inversion based on Eq. (41), the value of σ reg = 1 has worked well in convex inversion. For the present example asteroids, we are not utilizing any a priori threshold value σ pr,k (Eq. (36)) in the case A138, page 12 of 19  A&A 642, A138 (2020) that the extents of the clouds of MCMC sample solutions are realistic. Nevertheless, we would like to stress that, for Proserpina and Bilkis, the present results correspond to one potential solution regime, with the other regime differing from the present one about 180 • in pole longitude. For 30 000 MCMC samples with the current numbers of photometric observations, the simulations for each asteroid took some tens of hours of sequential computing time.
In addition to the photometric observations, Figs. 2-4 indicate uncertainty envelopes as derived from the MCMC simulations. Two kinds of envelopes are shown for each observation: 1 − σ uncertainty as computed from the MCMC samples and the complete minimum-maximum envelope following from the samples. Overall, the results are consistent. Figure 5 illustrates, for asteroid Lutetia, MCMC random-walk sampling results for the rotation period, pole longitude, pole latitude, example low-degree real-valued spherical harmonics coefficients a 20 and b 21 (derived from the complex-valued coefficients s 20 and s 21 ), and the photometric slope using the complete set of observations. The observations are characterized in Tables 2 and A.1 and in Fig. 2. It is clear that both a 20 and b 21 show substantial variegation that is nevertheless realistic: the 2500 virtual least-squares solutions computed for Lutetia show similar but slightly more constrained variegation.

(21) Lutetia
Overall, Fig. 5 shows that, even for 30 000 sample solutions, certain structure still remains in the illustrations of the marginal p.d.f.s. However, moments and correlations computed for the parameters are already converging. Table 3 collects, for Lutetia, the mean rotational and slope parameters with their the 1-σ uncertainties as computed from the full sample of 30 000 solutions. First of all, the rotation period and pole orientation are in fair agreement with the original results from convex inversion by Torppa et al. (2003): they obtained P = 8.165455 h, λ = 39 • , and β = 3 • based on 32 lightcurves. The rotation period and pole orientation are in excellent agreement with those by Carry et al. (2010Carry et al. ( , 2012, who give P = 8.168270 h, λ = 52 • , and β = −6 • . Our results are realistic for the rotational parameters: they indicate uncertainties of about a degree in the pole orientation and an uncertainty in the rotation period of some 9.6 × 10 −7 h ≈ 8.3 × 10 −2 s. The period uncertainty is small, comparable to uncertainties derived from space mission flybys.
The spherical harmonics coefficients a 20 and b 21 show substantial variegation. However, as illustrated in Fig. 6 for the least squares solution from convex inversion and one sample MCMC solution, the shapes are similar overall but show scaling in certain directions allowed for by the observations and, in particular, the illumination and observation geometries within the data sets. In Fig. 7, we compare our least-squares shape model to the non-convex shape model by Farnham (2013): the agreement validates our convex inversion method. As discussed in Sect. 2.1, the choice of the reflection coefficient affects the retrieved shape model. In our shape model, the extent along the z-axis is larger than that in the true shape. Our shape model in Fig. 7 has been produced from the spherical harmonics coefficients of the Gaussian surface density with the Minkowski solver (Kaasalainen et al. 1992a,b; reviewed in Sect. 3.5. Consider next the photometric slope parameter β S in Fig. 5 and in Table 3. Recalling that there are only 7 GDR2 observations included, the MCMC result is already accurate at β S = (2.050 ± 0.050) mag rad −1 . In order to provide an alternative  estimate for β S , we have carried out MCMC sampling with the triaxial ellipsoid model for solely the Gaia observations by assuming the rotation period and rotational pole orientation known from the solution with the complete set of observations. We obtain β S = (1.75 ± 0.25) mag rad −1 , which is in agreement with the result with all the observations included.

(26) Proserpina
For (26) Proserpina, we have carried out the MCMC parameter sampling in ways largely similar to that for (21) Lutetia, but with one important difference: since the sparse ground-based lightcurve originating from the Lowell Observatory photometric database (Bowell et al. 2014) includes observations also at small phase angles, we have utilized the full H, G 1 , G 2 phase function in the modeling of the photometric phase curve. The value for the slope β S then derives from the full H, G 1 , G 2 phase function. Figure 8 illustrates the marginal p.d.f.s for the rotational parameters, example spherical harmonics coefficients, and the photometric slope for the complete set of observations characterized in Tables 2 and A.2 and in Fig. 3. In Fig. 8 Fig. 12. As in Fig. 9 for asteroid (585) Bilkis.
uncertainties in the angular parameters are larger than those for Lutetia, which is in agreement with the smaller number of dense lightcurves for Proserpina. Figure 9 shows the convex inversion results using the GDR2 data only. There is clearly more extent and structure in the distributions, and there is a tendency for the rotation period to drift away from the value obtained with the complete data set. Figure 9 suggests multimodal distributions. We stress that the results for the latter case have been constrained, in terms of the plausible shape solutions (see Sect. 3.1), into the neighborhood of the least-squares solution using all the data. Using the Gaia data only can give rise to larger extents of the distributions when the constraints are relaxed. Table 4 illustrates the rotation period, pole orientation, and photometric slope and their 1-σ uncertainties for Proserpina. The results are again realistic for the rotational parameters: the uncertainty in the period is 6.4 × 10 −6 h ≈ 5.5 × 10 −1 s. The A138, page 15 of 19 A&A 642, A138 (2020) spherical harmonics coefficients a 20 and b 21 show similar substantial variegation as in Fig. 5 for Lutetia. The least squares and MCMC example shapes in Fig. 10 are similar, deviations are nevertheless visible. The rotational parameters agree well with those derived by Hanuš et al. (2016): they obtained P = 13.10977 h, λ = 88 • , and β = −52 • . The photometric slope parameter β S in Fig. 8 and in Table 4 shows values clearly smaller than those for Lutetia in Fig. 5 and in Table 3. In terms of the 1-σ uncertainties, the difference is significant. The analysis of the complete and GDR2 data sets with convex inversion gives β S = (1.669 ± 0.050) mag rad −1 and β S = (1.720 ± 0.109) mag rad −1 , respectively. The analysis of the Gaia data only with ellipsoid inversion gives β S = (1.62 ± 0.15) mag rad −1 . The analysis of the Gaia data only with ellipsoid inversion and assuming the rotation parameters fixed gives β S = (1.53 ± 0.15) mag rad −1 . Given their uncertainties, the four photometric slope values are in mutual agreement.

(585) Bilkis
We have carried out the MCMC sampling of the parameters in a way largely similar to that for (21) Lutetia. Figure 11 illustrates the marginal p.d.f.s for the rotational parameters, example spherical harmonics coefficients, and the photometric slope for the complete set of observations characterized in Tables 2 and A.3 and in Fig. 4. Our work provides the first analysis of Bilkis using convex inversion. Of the three example asteroids, the case of Bilkis exhibits the largest uncertainties in the pole orientation. Figure 12 shows the convex inversion results for Bilkis using the GDR2 data only. As compared to Fig. 11, there is clearly more extent and structure in the distributions. Multimodal distributions are suggested.
For (585) Bilkis, Table 5 illustrates the 1-σ uncertainties for the rotational and slope parameters. The results are again realistic for the rotational parameters. In particular, we have refined the original rotation period of Bilkis based on ellipsoid inversion  to P = (8.5750559 ± 0.0000026) h with a pole orientation of λ = 320.6 • ± 1.2 • and β = −25.6 • ± 1.7 • . Whereas the mirror pole orientation of λ + 180 • and −β cannot be ruled out, the present rotation period provides a substantially improved fit to the Gaia data. The uncertainty in the rotation period is thus about 2.2 × 10 −1 s. Again, the spherical harmonics coefficients a 20 and b 21 show substantial variegation, but the differences in the least squares and MCMC example shapes in Fig. 13 are realistic.
The photometric slope parameter β S in Fig. 11 and in Table 5 shows values clearly larger than those for Lutetia and Proserpina. This is not unexpected, because the three objects belong to different taxonomic classes. In terms of the 1-σ uncertainties, the deviation from the Lutetia value is statistically significant, and even more so for the deviation from the Proserpina value. The analysis of the complete data set with convex inversion gives β S = (2.292 ± 0.090) mag rad −1 , whereas the analysis of the Gaia data only gives β S = (2.152 ± 0.119) mag rad −1 . The analysis of the Gaia data only with ellipsoid inversion gives β S = (2.52 ± 0.10) mag rad −1 . The analysis of the Gaia data only with ellipsoid inversion and assuming the rotation parameters fixed gives β S = (2.43 ± 0.20) mag rad −1 . The photometric slope values are in satisfactory mutual agreement and agree with the taxonomical class of the asteroid (Table 2). Last but not least, for Bilkis, Lutetia, and Proserpina, the correlations among all the parameters remain only low or moderate: there are no extremely high correlations.

Conclusions
We have assessed the statistical inversion (Bayesian inference) of asteroid lightcurves for the rotation periods, pole orientations, convex shapes, and photometric phase-curve parameters. We have developed random-walk MCMC and importance samplers for the inverse problem. The samplers have been validated via case studies for asteroids (21) Lutetia, (26) Proserpina, and (585) Bilkis, and they are applicable to photometric data sparse or dense in time.
The future prospects for an extensive retrieval of asteroid photometric slopes from disk-integrated phase curves using the GDR2 and forecoming Gaia Data Release 3 (GDR3) data are outstanding. Whereas it is not feasible to derive pole orientations solely from the GDR2 photometry due to insufficient coverage in ecliptic longitudes, the GDR3 photometry can be extensive enough to facilitate pole retrievals. As to the computational demands, a massive application to tens of thousands of asteroids can be carried out using parallel computing, for example, by assigning a specific computer processor core to each asteroid.
It is possible to develop the MCMC and importance inverse methods further based on, for example, independence sampling and kernel estimation. In the potential future developments, the present methods can provide benchmarks. Likewise, it is A138, page 16 of 19 possible to further refine the model for the observational uncertainties and weights if deemed necessary.
The Markov-chain Monte Carlo and importance samplers, using a simplified model for the observational uncertainties, have been implemented in the Gaia Added-Value Interface Platform (GAVIP) by Torppa et al. (2018). The GAVIP portal complements and improves on the CPU-time-restricted fast computational tools based on triaxial Lommel-Seeliger ellipsoids implemented by Gaia DPAC. The user community can utilize the GAVIP methods in their own analyses of Gaia photometry. Currently, the GDR2 data needs to be extracted from the ESA GDR2 portal and then formatted and copied to the GAVIP portal.
Finally, the high-precision, milliarcsecond-to-microarcsecond stellar and asteroid astrometry by the Gaia mission is likely to result in accelerated ground-based occultation observations for the retrieval of asteroid sizes. These sizes promise to remove the ambiguity in resolving asteroid sizes and geometric albedos. Combining the geometric albedos with the phase curves will constitute no less than a major leap in asteroid taxonomical considerations based on geometric albedos and phase curves. Last but not least, together with the spectroscopy from the Gaia mission and ground-based polarimetric observations, they will allow for physics-based interpretation of asteroid surface properties. Notes. First, we give the lightcurve identifier (k), time span (T k ), mean sampling time interval (∆T k ), number of observations (N k ), and effective number of observations (N k,eff , Eq. (38)). Second, we give the number of nodes for the cubic spline fit based on the Bayesian information criterion (N BIC )), the rms-value of the spline fit in relative brightness (rms( )) and relative magnitude (rms(m)), and the resulting initial uncertainty (σ "ini"). Third, we give the resulting rms-values and uncertainty (σ "fin") for the best-fit model from convex inversion. In the third column, for example, 0.4658(−2) stands for 0.4658 × 10 −2 .