Issue 
A&A
Volume 642, October 2020



Article Number  A138  
Number of page(s)  19  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202038036  
Published online  13 October 2020 
Asteroid lightcurve inversion with Bayesian inference
^{1}
Department of Physics, University of Helsinki, Gustaf Hällströmin katu 2a, PO Box 64, 00014 U.
Helsinki, Finland
email: Karri.Muinonen@Helsinki.Fi
^{2}
Finnish Geospatial Research Institute FGI, Geodeetinrinne 2,
02430 Masala, Finland
^{3}
Space Systems Finland, Kappelitie 6, 02200 Espoo, Finland
^{4}
Yunnan Observatories, CAS, PO Box 110, Kunming 650216, PR China
^{5}
School of Astronomy and Space science, University of Chinese Academy of Sciences, Beijing 100049, PR China
^{6}
INAF, Osservatorio Astrofisico di Torino, Strada Osservatorio 20,
10025 Pino Torinese (TO), Italy
Received:
27
March
2020
Accepted:
9
August
2020
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 groundbased 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 Markovchain Monte Carlo sampler (MCMC) with a novel proposal probability density function based on the simulation of virtual observations giving rise to virtual leastsquares solutions. We utilized these leastsquares 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 groundbased 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 Sunobjectobserver 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.
Key words: minor planets, asteroids: general / radiative transfer / scattering / methods: numerical / methods: statistical
© K. Muinonen et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 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 diskintegrated 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 Sunasteroidobserver 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 Kaasalainen & Torppa (2001) and Kaasalainen et al. (2001). 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 inversion result, while the facet approach was shown to describe sharp features in lightcurves better. Convex shape models of highly nonconvex test objects were found to agree well with the convex hulls of the original shapes (Kaasalainen & Torppa 2001). 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 nonconvex, 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 (Viikinkoski et al. 2015), 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, diskresolved 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 photometricobservations. 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 analysisof real data (Kaasalainen & Ďurech 2007; Cellino et al. 2009; Carbognani et al. 2012; SantanaRos et al. 2015). The paper by SantanaRos 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 (SantanaRos et al. 2015). The 180degree ambiguity affects the inversion, in particular, for asteroids in lowinclination 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 groundbased 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 LommelSeeliger 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 LommelSeeliger scattering law. Again, the obtained results of photometric inversion have been found to be quite accurate. Muinonen et al. (2015) complemented these methods with a full Markovchain Monte Carlo treatment (MCMC) of the lightcurve inversion problem using LommelSeeliger ellipsoids. The genetic algorithm with the LommelSeeliger scattering ellipsoid (Cellino et al. 2015) has been implemented in the CU4 SolarSystemObject LongTerm pipeline (SSOLT) 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 groundbased 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 LommelSeeliger ellipsoids and showing also the first encouraging results using the socalled 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 leastsquares 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 phasecurve parameters. We formulate the statistical inversion problem with a detailed assessment of both random and systematic uncertainties. In Sect. 3, we describe theMCMC 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.
2 Statistical inversion
2.1 Asteroid model
We consideran asteroid in principalaxis rotation about its axis of maximum inertia, and denote the rotation period by P, the pole orientation in ecliptic longitude and latitude by (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 (1)
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 LommelSeeliger reflection coefficient (subscript LS) is (e.g., Lumme & Bowell 1981; Wilkman et al. 2015), (2)
where is the singlescattering albedo (proportion of incident light scattered, ), P_{11} is the singlescattering phase function (from the 4 × 4 scattering phase matrix), and α is the phase angle. We note that α derives unequivocally from ι, ɛ, and ϕ. R_{LS}, the firstorder multiplescattering approximation of the radiative transfer equation for a semiinfinite planeparallel medium (e.g., Chandrasekhar 1960), is applicable to weakly scattering media: this means that the intensity terms , 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 (3)
The recent works by, for example, Muinonen et al. (2018) and Väisänen et al. (2019) provide advanced multiplescattering models for planetary regoliths.
The diskintegrated brightness L equals the surface integral (4)
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, (5)
where Ω_{+} stands for the part of the unit sphere both illuminated and visible. With the LommelSeeliger reflection coefficient, the calculation of L in Eq. (5) can be carried out analytically, (6)
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 diskintegrated brightness of the asteroid and the diskintegrated brightness of a normally illuminated Lambertian disk, both with the same effective diameter D, in the exact backscattering direction α = 0°: (7)
The geometric albedo can obtain values larger than unity (p ≥ 0). For example, fora plane mirror oriented towards the light source, the geometric albedo approaches infinity.
Consider the singlescattering albedo and phase function P_{11}(α) in the LommelSeeliger reflection coefficient R_{LS} (Eq. (2)). On one hand, the diskintegrated brightness phase curve of an asteroid in a given apparition can be described by the H, G_{1}, G_{2} phase function (Muinonen et al. 2010; Penttilä et al. 2016). On the other hand, the reflection coefficient R_{LS} in Eq. (2) includes the 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 LommelSeeliger sphere is, by using L(0) from Eq. (6) in Eq. (7), (8)
Thus, it is reasonable to assume that the model (9)
can serve well in asteroid phasecurve analyses. A normalization to the LommelSeeliger 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 (Muinonen et al. 2010).
The V band magnitude of a spherical asteroid is (11)
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 10based 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): (12)
The linear dependence can be ascertained, for example, from the phase curves in Fig. 1 (see also Table 1). The phaseangle 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 lighttimecorrected directions of the Sun and the observer, respectively, as seen in the principal axes reference frame of the ellipsoid. The diskintegrated brightness of the LommelSeeliger ellipsoid is given by (13)
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 (15)
where the Y_{lm}’s are the orthogonal spherical harmonics and the complexvalued coefficients s_{lm} ensure a realvalued surface density (the superscript * refers to complex conjugation): (16)
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 independent real or imaginary parts of the coefficients s_{lm}. When utilizing realvalued spherical harmonics, we utilize wellbehaving normalized associated Legendre functions. This entails multiplication of the normal associated Legendre functions by .
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 LommelSeeliger reflection coefficient (Eq. (2)) instead of the frequentlyused combinationof the LommelSeeliger and Lambert coefficients (e.g., Kaasalainen & Torppa 2001; Kaasalainen et al. 2001). 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 singlescattering albedo and phase function in Eq. (2) do not vary across the surface.
Fig. 1 Photometric phase curves in relative magnitudes pertaining to the five taxonomical classes of asteroids according to Penttilä et al. (2016): E (solid blue line, shallow), S/M (dashed red line), C (dotted yellow line), P (dashdotted purple line), and D (solid green line, steep). Normalization at the phase angles of α = 0° (left) and α = 20° (right). See Table 1. 
Photometric G_{1} and G_{2} parameters for different asteroid taxonomic classes in the singleparameter phase function by Penttilä et al. (2016) based on the photometric observations reported in Shevchenko et al. (2016).
2.2 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}, …, , where the parameters are, respectively, the rotation period, ecliptic pole longitude, ecliptic pole latitude, rotational phase at t_{0}, shape parameters , 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 setof parameters. The total number of parameters equals (17)
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 representthe physical diameter of an equalvolume 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 twoparameter H, G_{12} phase function (Muinonen et al. 2010), including G_{12} among the parameters. Furthermore, we set l_{min} = 0. In this case,for general convex shapes, (18)
and our parameters are . For triaxial ellipsoid shapes, we have . We may also utilize the singleparameter 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 SantanaRos et al. (2015).
Let N and K be the number of observations and number of lightcurves, respectively. Then, (19)
where N_{k} specifies the number of observations in lightcurve k (k = 1, …, K). Let (20)
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: (21)
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 phasecurve 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 (22)
The covariance matrices Λ_{ɛ} and Λ_{υ} are blockdiagonal. As for Λ_{ɛ}, it may or may not include offdiagonal 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 offdiagonal 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 “ObservedComputed” (OC) residual magnitudes ΔM(P), (23)
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 (24)
where χ^{2} measures the OC 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 (25)
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): (26)
Furthermore, define the magnitude differences m_{obs,kj} and m_{kj}(P) with the help of Eq. (26): (27)
Introduce then the relative brightnesses ℓ_{obs,kj} and ℓ_{kj}(P): (28)
In the case of large uncertainties υ, that is, in the case of what can be called relative photometry, (29)
Recalling that Λ_{υ} includes high correlations among the uncertainties for the observations of a given lightcurve, (30)
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): (31)
In the present case of large systematic uncertainties, the absolute magnitude information is lost. Subsequently, in the last two steps of Eq. (31), we rescale 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 leastsquares 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 leastsquares 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, (32)
For a diagonal matrix Λ_{ɛ}, (33)
If the uncertainties ɛ are small, the OC magnitude difference can be linearized into a difference in the corresponding brightnesses L_{obs,kj} and L_{kj} (P): (34)
where 2.5lge ≈ 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 (35)
where σ_{ɛ,k} now describes the uncertainty of the observations in lightcurve k and is taken to account for the factor of 2.5lge. 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 Muinonen et al. (2015) and Torppa et al. (2018), and derives from defining the inverse problem for magnitudes (Eq. (21)) instead of brightnesses.
3 Numerical inverse methods
3.1 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 (36)
where σ_{0,k} denotes the rmsvalue 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 (37)
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 lightcurve contributes a χ^{2} value of unity in the inverse problem, whereas each other lightcurve k is considered to split into approximately (38)
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 (39)
where, again, σ_{0,k} denotes the rmsvalue 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 lightcurve corresponds to the one with the minimum number of observations and contributes a χ^{2} value equal to (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 (40)
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 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 that depends on the shape parameters included in P: (41)
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 value of Eq. (41) should vanish. Here the triangular discretization needs to have a resolution high enough to warrant an accurate computation of the diskintegrated 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.
3.2 Virtualobservation MCMC sampler
Markovchain 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 MetropolisHastings algorithm that is based on the computation of the ratio a_{r} : (42)
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]: (43)
that is, the proposed parameters are accepted with the probability of min (1, a_{r}). After a number of transitions in the socalled 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 virtualobservation MCMC methods, earlier studied for asteroid orbital inversion (Muinonen et al. 2012) and lightcurve inversion with LommelSeeliger ellipsoids (Muinonen et al. 2015), 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}, (44)
Thereafter, virtual leastsquares parameters P_{v} are obtained by minimizing (45)
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 Ndimensional integral (46)
where p(M_{v}) is the Gaussian p.d.f. for the virtual observations and δ is Dirac’s Delta function.
A randomwalk MCMCsampler is obtained with the help of the symmetric proposal p.d.f., defined as the 2(N+ N_{p})dimensional integral (47)
where the integration over N_{p} parameters results in a (2N + N_{p})dimensional convolution integral (48)
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 MetropolisHastings 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 (j = 1, 2, 3…, N_{v}). Parameter differences are then readily available from (49)
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 differences available, that is, a number that can become extremely large for reasonable values of N_{v}. Any pair of virtual leastsquares 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 (50)
The proposals can now be computed during the MCMC sampling using the N_{v} virtual leastsquares 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 leastsquares 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 virtualobservation and Gaussian proposals results in a hybrid proposal probability density that will be used in the MCMC sampling below.
3.3 Virtualobservation importance sampler
In what follows, we describe a virtualobservation randomwalk importance sampler, utilizing the probability density for the virtual leastsquares 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), (51)
where P_{0} denotes the leastsquares parameters and the cutoff is set on the basis of the number of parameters. For the present inverse problem, we may define with the help of probability thresholds for multivariate normal statistics: (52)
where P_{c} denotes the cutoff 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 (53)
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 (54)
that is, transitions are always accepted within the given regime but never across the boundary.
The proposal procedure utilized in the randomwalk 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 (55)
For highdimensional inverse problems, with the virtualobservation importance sampler, it is challenging to generate sufficient numbers of samples in the highprobabilitymass 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.
3.4 Leastsquares optimization
In the case of general convex shapes, the LevenbergMarquardt leastsquares optimization method (Press et al. 1992) is utilized in solving for the leastsquares parameter sets described in the preceding subsections. In order to guarantee stable convergence towards the leastsquares solution, an additional regularization term is incorporated in the total χ^{2} describing the goodness of fit. The numerical integration of the diskintegrated brightness is carried out using a triangulation of the unit sphere. The regularizing 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 leastsquares optimization is carried out using the NelderMead downhill simplex method (Press et al. 1992). These methods are described, in detail, by Muinonen et al. (2015).
3.5 Solution of the Minkowski problem
Generating the threedimensional 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 socalled tangent plane (the plane that tangents the surface at (θ, φ)): (56)
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 (57)
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 (58)
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 (59)
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 quicklook analysis, we have developed an alternative method for the derivation of the polyhedral shape from the Gaussian surface density. We start from a triangulated model of a spherical asteroid and treat the coordinates of its N_{n} nodes i = 1, …, N_{n} as the 3N_{n} free parameters. The derivation of the leastsquares 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.
Fig. 2 Example photometric lightcurves for asteroid (21) Lutetia. We show the groundbased photometric lightcurves#41 and #3 and the Gaia lightcurve (blue circles; left, middle, and right, respectively; see Table A.1) together with the bestfit 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. 
4 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 highresolution 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 moderatealbedo, possibly Mclass asteroid (Xcomplex asteroid in more modern classification) with 7 GDR2 observations, Bilkis is a lowalbedo C/D/Pclass asteroid with a substantial number of 32 GDR2 observations, and Proserpina is a moderatetohighalbedo Sclass asteroid with a larger number of 39 GDR2 observations. To summarize, the selected asteroids represent different taxonomical classes with differing numbers of GDR2 observations.
4.1 Initial data analysis
Table 2 includes the Tholen taxonomical classification of the asteroids and summarizes the groundbased and spacebased observations. For Lutetia, Proserpina, and Bilkis, there are dense groundbased lightcurves and sparse GDR2 lightcurves. For Proserpina, there is additionally sparse groundbased data (Ďurech et al. 2010; Bowell et al. 2014). Example lightcurves, including the Gaia lightcurve, are shown in Figs. 2–4 for the three asteroids.
For Lutetia, Proserpina, and Bilkis, we have utilized altogether 50, 29, and 17 dense groundbased lightcurves and one sparse GDR2 lightcurve, respectively (Tables A.1–A.3). For the first two asteroids, we have extracted the groundbased observations from the DAMIT database (Ďurech et al. 2010). Additionally, the data setfor Proserpina includes a sparse groundbased lightcurve also extracted from the DAMIT database (Ďurech et al. 2010; Bowell et al. 2014). For Bilkis, the 17 groundbased lightcurves are from Wang et al. (2016) and the Minor Planet Center. The number and total timespan of the observations are, respectively for Lutetia, Proserpina, and Bilkis, 4019 observations in about 53 yr, 7266 observations in over 51 yr, 1876 observations in about 15 yr.
Tables A.1–A.3 include a considerable amount of information about the characteristics of the photometric observations, including the references. For each individual lightcurve, the tables include the time span, mean sampling time interval, number of observations, and the effective number of observations N_{k,eff} as computed based on the uncertainty model of Sect. 3.1. Furthermore, the tables include the number of nodes for the cubic spline fits based on the Bayesian information criterion, the corresponding rmsvalues of the spline fit in relative brightness and relative magnitude, and the resulting initial uncertainties of the observations. For comparison, the rmsvalues and final uncertainties for the bestfit model from convex inversion are also given (see below for additional information). The spline fits reveal the varying level of intrinsic uncertainties. These uncertainties correlate with the final postfit uncertainties.
Consider first the case of Lutetia in Table A.1. For the 50 dense groundbased 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.
For (26) Proserpina (Table A.2), the number of observations in dense groundbased 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 groundbased lightcurve #30 receives the same weight as the sparse spacebased 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 groundbased lightcurves varies from 25 to 505. The effective number of observationsvaries from 1 to approximately 14.
Lightcurve characteristics for asteroids (21) Lutetia, (26) Proserpina, and (585) Bilkis incorporated in the present study.
Fig. 3 Example photometric lightcurves for asteroid (26) Proserpina. We depict the groundbased lightcurves #12 and #5 and the Gaia lightcurve. For more information, see Fig. 2 and Table A.2. 
Fig. 4 Example photometric lightcurves for asteroid (585) Bilkis. We depict the groundbased lightcurves #6 and #4 and the Gaia lightcurve. For more information, see Fig. 2 and Table A.3. 
Fig. 5 Markovchain Monte Carlo sampling of the rotational, shape, and photometric phasecurve parameters for asteroid (21) Lutetia. The marginal probabilitity densities are characterized by an unbiased subsample of 1000 parameter sets selected from the complete sample of 30 000 parameter 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}) vs. rotationperiod (down to the right). The red, dotted line marks the least squares rotation period from convex inversion. 
Rotation period P (h), pole longitude λ (°) and latitude β (°), and phasecurve slope parameter β_{S} (mag rad^{−1}) for (21) Lutetia.
Fig. 6 Quicklook leastsquares convex shape model for asteroid (21) Lutetia (top, two rows) accompanied by an example shape model from MCMC sampling (bottom, two rows). 
Fig. 7 Leastsquares convex shape model for asteroid (21) Lutetia reconstructed using the conventional Minkowski problem solver (top two rows, see Sect. 3.5). For comparison, we show the 6kresolution shape model (bottom two rows) by Farnham (2013). 
4.2 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 Proserpinaand 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 observationswith the help of rmsvalues 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 of the groundbased lightcurves. For the Gaia data, we set the uncertainty at σ_{ɛ} = 0.01.
The results are based on the utilization of 2500 virtualleastsquares convex inversion solutions for Lutetia and Bilkis and 2000 solutions for Proserpina. Random Gaussian errors are introduced for each lightcurve with the standard deviation equal to the rmsvalue of the lightcurve. A hybrid proposal probability density with equal weights for the virtualobservation and Gaussian proposals is utilized as the proposal p.d.f. in MCMC sampling (see Sect. 3.2). The hybrid modeling removes the need to derive large numbers of virtual leastsquares solutions. Altogether 30 000 MCMC sample solutions have been generated for each asteroid and unbiased subsamples of 1000 solutions randomly selected from the complete set of solutions are used to depict marginal p.d.f.s. For all three cases, we have verified, with the virtualobservation importance sampler of Sect. 3.3, 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 obser vation: 1 − σ uncertainty as computed from the MCMC samples and the complete minimummaximum envelope following from the samples. Overall, the results are consistent.
4.2.1 (21) Lutetia
Figure 5 illustrates, for asteroid Lutetia, MCMC randomwalk sampling results for the rotation period, pole longitude, pole latitude, example lowdegree realvalued spherical harmonics coefficients a_{20} and b_{21} (derived from the complexvalued 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 leastsquares solutions computed for Lutetia show similar but slightly more constrained variegation.
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. (2010, 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 leastsquares shape modelto the nonconvex 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 zaxis 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; Kaasalainen & Torppa 2001) 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.
4.2.2 (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 groundbased lightcurve originating from the Lowell Observatory photometricdatabase (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 fullH, 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, structures show up in the twodimensional plots of sample solutions so that even more than 30 000 sample solutions should be incorporated. But, as for Lutetia, the statistical moments and correlations are converging with the current number of samples. Overall, the 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 leastsquares 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 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.
4.2.3 (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 (Wang et al. 2016) 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.
5 Conclusions
We have assessed the statistical inversion (Bayesian inference) of asteroid lightcurves for the rotation periods, pole orientations, convex shapes, and photometric phasecurve parameters. We have developed randomwalk 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 diskintegrated 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 possible to further refine the model for the observational uncertainties and weights if deemed necessary.
The Markovchain Monte Carlo and importance samplers, using a simplified model for the observational uncertainties, have been implemented in the Gaia AddedValue Interface Platform (GAVIP) by Torppa et al. (2018). The GAVIP portal complements and improves on the CPUtimerestricted fast computational tools based on triaxial LommelSeeliger 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 highprecision, milliarcsecondtomicro arcsecond stellar and asteroid astrometry by the Gaia mission is likely to result in accelerated groundbased occultation observationsfor 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 taxonomicalconsiderations based on geometric albedos and phase curves. Last but not least, together with the spectroscopy from the Gaia mission and groundbased polarimetric observations, they will allow for physicsbased interpretation of asteroid surface properties.
Acknowledgments
Research supported, in part, by the Academy of Finland Grants No. 1298137 and 1325805, by the ERC Advanced Grant No 320773 (SAEMPL), and by the National Natural Science Foundation of China Grant No. 11673063. We would like to thank Dr. Olli Wilkman, Mr. Tuomo Pieniluoma, and Dr. Maria Gritsevich for valuable help during the course of the research. A.C. worked on this project using funds from European COST Action CA18104 and from the University of Helsinki. We are thankful to the anonymous referees for their insightful and constructive reviews. This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia.
Appendix A Tables of lightcurve characteristics
Lightcurve characteristics for (21) Lutetia.
References
 Bartczak, P., & Dudziński, G. 2018, MNRAS, 473, 5050 [NASA ADS] [CrossRef] [Google Scholar]
 Barucci, M. A., Cellino, A., De Sanctis, C., et al. 1992, A&A, 266, 385 [Google Scholar]
 Bowell, E., Oszkiewicz, D., Wasserman, L., et al. 2014, MAPS, 49, 95 [Google Scholar]
 Carbognani, A., Tanga, P., Cellino, A., et al. 2012, Planet. Space. Sci., 73, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Carry, B., Kaasalainen, M., Leyrat, C., et al. 2010, A&A, 523, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carry, B., Kaasalainen, M., Merline, W., et al. 2012, Planet. Space. Sci., 66, 200 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cellino, A., Zappalà, V., & Farinella, P. 1989, Icarus, 78, 298 [NASA ADS] [CrossRef] [Google Scholar]
 Cellino, A., Hestroffer, D., Tanga, P., Mottola, S., & Dell’Oro, A. 2009, A&A, 506, 935 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cellino, A., Muinonen, K., Hestroffer, D., & Carbognani, A. 2015, Planet. Space. Sci., 118, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Cellino, A., Hestroffer, D., Lu, X.P., Muinonen, K., & Tanga, P. 2019, A&A, 631, A67 [CrossRef] [EDP Sciences] [Google Scholar]
 Chandrasekhar, S. 1960, Radiative Transfer (New York: Dover) [Google Scholar]
 Chang, Y., & Chang, C. 1963, Acta Astron. Sin., 11, 139 [Google Scholar]
 Chang, Y., Zhou, X.h., Yang, X.y., et al. 1981, Acta Astron. Sin., 22, 169 (Chinese. English translation – Chin. Astron. Astrophys. 5, 434 [Google Scholar]
 Denchev, P., Magnusson, P., & Donchev, Z. 1998, Planet. Space. Sci., 46, 673 [NASA ADS] [CrossRef] [Google Scholar]
 Dotto, E., Barucci, M., Fulchignoni, M., et al. 1992, A&AS, 95, 195 [Google Scholar]
 Ďurech, J., & Hanuš, J. 2018, A&A, 620, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ďurech, J., Carry, B., Delbo, D., Kaasalainen, M., & Viikinkoski, M. 2015, in Asteroids IV, eds. P. Michel, F. DeMeo, & W. Bottke Jr. (Tucson, AZ: University of Arizona Press), 183 [Google Scholar]
 Ďurech, J., Sidorin, V., & Kaasalainen, M. 2010, A&A, 513, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Farnham, T. 2013, Shape model of asteroid 21 Lutetia, NASA Planetary Data System, rOAOSINAC/OSIWAC5LUTETIASHAPEV1.0 [Google Scholar]
 Gaia Collaboration (Spoto, F., et al.) 2018a, A&A, 616, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018b, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hanuš, J., Durech, J., Oszkiewicz, D. A., et al. 2016, A&A, 586, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kaasalainen, M., & Torppa, J. 2001, Icarus, 153, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Kaasalainen, M., & Ďurech, J. 2007, Proc. IAU Symp. 236, 151 [Google Scholar]
 Kaasalainen, M., Lamberg, L., & Lumme, K. 1992a, A&A, 259, 333 [NASA ADS] [Google Scholar]
 Kaasalainen, M., Lamberg, L., Lumme, K., & Bowell, E. 1992b, A&A, 259, 318 [Google Scholar]
 Kaasalainen, M., Torppa, J., & Muinonen, K. 2001, Icarus, 153, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Kaasalainen, M., Mottola, S., & Fulchignoni, M. 2002, in Asteroids III, eds. W. Bottke, A. Cellino, P. Paolicchi, & R. Binzel (Tucson, AZ: University of Arizona Press), 139 [Google Scholar]
 Lagerkvist, C.I., Erikson, A., Debehogne, H., et al. 1995, A&AS, 113, 115 [NASA ADS] [Google Scholar]
 Lamberg, L. 1993, PhD thesis, University of Helsinki, Finland [Google Scholar]
 Lumme, K., & Bowell, E. 1981, AJ, 86, 1694 [NASA ADS] [CrossRef] [Google Scholar]
 Lupishko, D. F., Belskaya, I. N., & Tupieva, F. A. 1983, Sov. Astron. Lett., 9, 691 [Google Scholar]
 Lupishko, D., Velichko, F., Belskaya, I., & Shevchenko, V. 1987, Kinem. Fiz. Nebesn. Tel., 3, 36 [Google Scholar]
 Minkowski, H. 1903, Math. Ann., 57, 447 [CrossRef] [Google Scholar]
 Muinonen, K., & Bowell, E. 1993, Icarus, 104, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Muinonen, K., & Lumme, K. 2015, A&A, 584, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Muinonen, K., Belskaya, I., Cellino, A., et al. 2010, Icarus, 209, 542 [NASA ADS] [CrossRef] [Google Scholar]
 Muinonen, K., Granvik, M., Oszkiewicz, D., Pieniluoma, T., & Pentikäinen, H. 2012, Planet. Space. Sci., 73, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Muinonen, K., Wilkman, O., Cellino, A., Wang, X., & Wang, Y. 2015, Planet. Space. Sci., 118, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Muinonen, K., Markkanen, J., Väisänen, T., Peltoniemi, J., & Penttilä, A. 2018, Opt. Lett., 43, 683 [NASA ADS] [CrossRef] [Google Scholar]
 O’Hagan, A., & Forster, J. J. 2004, Kendall’s Advanced Theory of Statistics: Bayesian Inference, 2nd edn. (New Jersey: John Wiley & Sons), 2B [Google Scholar]
 Oszkiewicz, D., Muinonen, K., Virtanen, J., & Granvik, M. 2010, Metroid. Planet. Sci., 44, 1897 [NASA ADS] [CrossRef] [Google Scholar]
 Penttilä, A., Shevchenko, V., Wilkman, O., & Muinonen, K. 2016, Planet. Space. Sci., 123, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in C: The Art of Scientific Computing, 2nd edn. (Cambridge, USA: Cambridge University Press) [Google Scholar]
 Russell, H. N. 1906, ApJ, 24, 1 [NASA ADS] [CrossRef] [Google Scholar]
 SantanaRos, T., Bartczak, P., Michałowski, T., Tanga, P., & Cellino, A. 2015, MNRAS, 450, 333 [NASA ADS] [CrossRef] [Google Scholar]
 Scaltriti, F., & Zappalà, V. 1979, Icarus, 39, 124 [CrossRef] [Google Scholar]
 Shevchenko, V., Belskaya, I., Muinonen, K., et al. 2016, Planet. Space. Sci., 123, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Torppa, J., Kaasalainen, M., Michałowski, T., et al. 2003, Icarus, 164, 346 [NASA ADS] [CrossRef] [Google Scholar]
 Torppa, J., Hentunen, V.P., Pääkkönen, P., Kehusmaa, P., & Muinonen, K. 2008, Icarus, 198, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Torppa, J., Granvik, M., Penttilä, A., et al. 2018, Adv. Space Res., 62, 464 [NASA ADS] [CrossRef] [Google Scholar]
 Väisänen, T., Markkanen, J., Penttilä, A., & Muinonen, K. 2019, PLoS One, 14, 1 [CrossRef] [Google Scholar]
 Vernazza, P., Jorda, L., Sevecek, P., et al. 2020, Nat. Astron., 4, 136 [CrossRef] [Google Scholar]
 Viikinkoski, M., Kaasalainen, M., & Durech, J. 2015, A&A, 576, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, X., Muinonen, K., & Wang, Y. 2015a, Planet. Space. Sci., 118, 242 [CrossRef] [Google Scholar]
 Wang, X., Muinonen, K., Wang, Y., et al. 2015b, A&A, 581, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, A., Wang, X.B., Muinonen, K., Han, X., & Wang, Y.B. 2016, Res. A&A, 16 [Google Scholar]
 Wilkman, O., Muinonen, K., & Peltoniemi, J. 2015, Planet. Space. Sci., 118, 250 [CrossRef] [Google Scholar]
 Zappala, V., di Martino, M., Knezevic, Z., & Djurasevic, G. 1984, A&A, 130, 208 [NASA ADS] [Google Scholar]
All Tables
Photometric G_{1} and G_{2} parameters for different asteroid taxonomic classes in the singleparameter phase function by Penttilä et al. (2016) based on the photometric observations reported in Shevchenko et al. (2016).
Lightcurve characteristics for asteroids (21) Lutetia, (26) Proserpina, and (585) Bilkis incorporated in the present study.
Rotation period P (h), pole longitude λ (°) and latitude β (°), and phasecurve slope parameter β_{S} (mag rad^{−1}) for (21) Lutetia.
All Figures
Fig. 1 Photometric phase curves in relative magnitudes pertaining to the five taxonomical classes of asteroids according to Penttilä et al. (2016): E (solid blue line, shallow), S/M (dashed red line), C (dotted yellow line), P (dashdotted purple line), and D (solid green line, steep). Normalization at the phase angles of α = 0° (left) and α = 20° (right). See Table 1. 

In the text 
Fig. 2 Example photometric lightcurves for asteroid (21) Lutetia. We show the groundbased photometric lightcurves#41 and #3 and the Gaia lightcurve (blue circles; left, middle, and right, respectively; see Table A.1) together with the bestfit 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. 

In the text 
Fig. 3 Example photometric lightcurves for asteroid (26) Proserpina. We depict the groundbased lightcurves #12 and #5 and the Gaia lightcurve. For more information, see Fig. 2 and Table A.2. 

In the text 
Fig. 4 Example photometric lightcurves for asteroid (585) Bilkis. We depict the groundbased lightcurves #6 and #4 and the Gaia lightcurve. For more information, see Fig. 2 and Table A.3. 

In the text 
Fig. 5 Markovchain Monte Carlo sampling of the rotational, shape, and photometric phasecurve parameters for asteroid (21) Lutetia. The marginal probabilitity densities are characterized by an unbiased subsample of 1000 parameter sets selected from the complete sample of 30 000 parameter 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}) vs. rotationperiod (down to the right). The red, dotted line marks the least squares rotation period from convex inversion. 

In the text 
Fig. 6 Quicklook leastsquares convex shape model for asteroid (21) Lutetia (top, two rows) accompanied by an example shape model from MCMC sampling (bottom, two rows). 

In the text 
Fig. 7 Leastsquares convex shape model for asteroid (21) Lutetia reconstructed using the conventional Minkowski problem solver (top two rows, see Sect. 3.5). For comparison, we show the 6kresolution shape model (bottom two rows) by Farnham (2013). 

In the text 
Fig. 8 As in Fig. 5 for asteroid (26) Proserpina. 

In the text 
Fig. 9 As in Fig. 5 for asteroid (26) Proserpina using the GDR2 observations only. 

In the text 
Fig. 10 As in Fig. 6 for asteroid (26) Proserpina. 

In the text 
Fig. 11 As in Fig. 5 for asteroid (585) Bilkis. 

In the text 
Fig. 12 As in Fig. 9 for asteroid (585) Bilkis. 

In the text 
Fig. 13 As in Fig. 6 for asteroid (585) Bilkis. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.