Open Access
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/0004-6361/202038036
Published online 13 October 2020

© K. Muinonen et al. 2020

Licence Creative CommonsOpen 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 disk-integrated 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 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 non-convex 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 non-convex 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 (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, 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 non-convex 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 ba and ca axial ratios for the triaxial shape (ellipsoid semiaxes abc), 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; 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 180-degree 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. Muinonen et al. (2015) 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 (Cellino et al. 2015) 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 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 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 (J2000.0, T stands for transpose), and the rotational phase at a given epoch t0 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 πF0 (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 Lommel-Seeliger reflection coefficient (subscript LS) is (e.g., Lumme & Bowell 1981; Wilkman et al. 2015), (2)

where is the single-scattering albedo (proportion of incident light scattered, ), P11 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 ϕ. RLS, 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 ≥ 2 are assumed negligible. When omitting polarization effects, the scattering phase function P11 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 multiple-scattering models for planetary regoliths.

The disk-integrated 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 Lommel-Seeliger 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 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°: (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 single-scattering albedo and phase function P11(α) in the Lommel-Seeliger reflection coefficient RLS (Eq. (2)). On one hand, the disk-integrated brightness phase curve of an asteroid in a given apparition can be described by the H, G1, G2 phase function (Muinonen et al. 2010; Penttilä et al. 2016). On the other hand, the reflection coefficient RLS in Eq. (2) includes the term that produces the phase function ΦLS given in Eq. (6). Thus, modeling P11(α) so that the phase curve from Eq. (2) for a spherical asteroid is precisely that of the H, G1, G2 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), (8)

Thus, it is reasonable to assume that the model (9)

resulting in (10)

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, G1, G2 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 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): (12)

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 G1, G2 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 (13)

where (14)

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 Ylm’s are the orthogonal spherical harmonics and the complex-valued coefficients slm ensure a real-valued surface density (the superscript * refers to complex conjugation): (16)

The series starts from a minimum degree lmin (usually lmin = 0) and is truncated at a maximum degree lmax (typically lmax ≤ 8, see Sect. 4.2 below), and there are altogether independent real or imaginary parts of the coefficients slm. When utilizing real-valued spherical harmonics, we utilize well-behaving 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 Lommel-Seeliger reflection coefficient (Eq. (2)) instead of the frequently-used combinationof the Lommel-Seeliger 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 single-scattering albedo and phase function in Eq. (2) do not vary across the surface.

thumbnail 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 (dash-dotted purple line), and D (solid green line, steep). Normalization at the phase angles of α = 0° (left) and α = 20° (right). See Table 1.

Table 1

Photometric G1 and G2 parameters for different asteroid taxonomic classes in the single-parameter 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 t0 by the vector P = (P, λ, β, φ0, s00, …, , where the parameters are, respectively, the rotation period, ecliptic pole longitude, ecliptic pole latitude, rotational phase at t0, shape parameters , geometric albedo, two parameters of the H, G1, G2 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, G1, G2, 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 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, G12 phase function (Muinonen et al. 2010), including G12 among the parameters. Furthermore, we set lmin = 0. In this case,for general convex shapes, (18)

and our parameters are . For triaxial ellipsoid shapes, we have . We may also utilize the single-parameter phase function provided by Penttilä et al. (2016): in this case, the G1, G2 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 G12 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, (19)

where Nk 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 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 (22)

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 pp be the a posteriori probability density function (p.d.f.) for the parameters. Within the Bayesian framework (cf. Muinonen & Bowell 1993), pp is proportional to the a priori and observational uncertainty p.d.f.s ppr and pɛ+υ, the latter being evaluated for the “Observed-Computed” (O-C) 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 ppr will describe, for example, the regularization needed in convex inversion. The final a posteriori p.d.f. is thus (24)

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 (25)

where Mobs,k, Mk (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 Mobs,k0 and Mk0(P) and their difference ΔMk0(P): (26)

Furthermore, define the magnitude differences mobs,kj and mkj(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 re-scale the observed and theoretical mean magnitudes to coincide with one another, resulting in ΔMk0(P) → 0. In other words, we can introduce a systematic correction (or, formally, an additional parameter) for each lightcurve to exactly result in ΔMk0(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, (32)

For a diagonal matrix Λɛ, (33)

If the uncertainties ɛ are small, the O-C magnitude difference can be linearized into a difference in the corresponding brightnesses Lobs,kj and Lkj (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 ΔMk0(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 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, ΔTk is the mean sampling time interval of the lightcurve k and (37)

marking the lightcurve with the lowest sampling rate in time and Tk 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 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 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 Aj(P) and n(P) denote the area and unit outward normal vector of the triangular element j = 1, …, N. In order for Aj(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 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.

3.2 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 ar : (42)

Here Pj and P′ denote the current and proposed parameters in a Markov chain, respectively, and pt (P′;Pj) is the proposal p.d.f. from Pj 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, ar). After a number of transitions in the so-called burn-in phase, the Markov chain, in the case of success, converges to sampling the target p.d.f. pp. 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 (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 Mv is generated from the original observations Mobs through the addition of Gaussian random uncertainties ɛv with zero means and covariance matrix Λv, (44)

Thereafter, virtual least-squares parameters Pv 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 Pv are random variables, too. The p.d.f. for Pv is formally the N-dimensional integral (46)

where p(Mv) 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+ Np)-dimensional integral (47)

where the integration over Np parameters results in a (2N + Np)-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 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 Pv before the MCMC sampling (cf. Wang et al. 2015a,b). The generation of virtual observation sets Lv and the derivation of Pv are repeated until there are Nv ≫ 1 sets of parameters (j = 1, 2, 3…, Nv). 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 Nv. 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 Nv is large enough, the replacement can be unnecessary.

The ratio ar for the decision criterion is expressed in the concise form (50)

The proposals can now be computed during the MCMC sampling using the Nv virtual least-squares element sets Pv,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.

3.3 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), (51)

where P0 denotes the least-squares parameters and the cut-off 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 Pc 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 NP ≈ 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 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 (55)

For high-dimensional inverse problems, with the virtual-observation 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.

3.4 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 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 Muinonen et al. (2015).

3.5 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 (θ, φ)): (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 Aj(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 aj 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 a triangulated model of a spherical asteroid and treat the coordinates of its Nn nodes i = 1, …, Nn as the 3Nn 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.

thumbnail Fig. 2

Example photometric lightcurves for asteroid (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.

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 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 moderate-to-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.

4.1 Initial data analysis

Table 2 includes the Tholen taxonomical classification of the asteroids and summarizes the ground-based and space-based observations. For Lutetia, Proserpina, and Bilkis, there are dense ground-based lightcurves and sparse GDR2 lightcurves. For Proserpina, there is additionally sparse ground-based data (Ďurech et al. 2010; Bowell et al. 2014). Example lightcurves, including the Gaia lightcurve, are shown in Figs. 24 for the three asteroids.

For Lutetia, Proserpina, and Bilkis, we have utilized altogether 50, 29, and 17 dense ground-based lightcurves and one sparse GDR2 lightcurve, respectively (Tables A.1A.3). For the first two asteroids, we have extracted the ground-based observations from the DAMIT database (Ďurech et al. 2010). Additionally, the data setfor Proserpina includes a sparse ground-based lightcurve also extracted from the DAMIT database (Ďurech et al. 2010; Bowell et al. 2014). For Bilkis, the 17 ground-based 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.1A.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 Nk,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 rms-values of the spline fit in relative brightness and relative magnitude, and the resulting initial uncertainties of the observations. For comparison, the rms-values and final uncertainties for the best-fit 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 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.

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 observationsvaries from 1 to approximately 14.

Table 2

Lightcurve characteristics for asteroids (21) Lutetia, (26) Proserpina, and (585) Bilkis incorporated in the present study.

thumbnail Fig. 3

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

thumbnail Fig. 4

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

thumbnail Fig. 5

Markov-chain Monte Carlo sampling of the rotational, shape, and photometric phase-curve parameters for asteroid (21) Lutetia. The marginal probabilitity densities are characterized by an unbiased sub-sample 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 a20 and b21 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.

Table 3

Rotation period P (h), pole longitude λ (°) and latitude β (°), and phase-curve slope parameter βS (mag rad−1) for (21) Lutetia.

thumbnail Fig. 6

Quick-look least-squares convex shape model for asteroid (21) Lutetia (top, two rows) accompanied by an example shape model from MCMC sampling (bottom, two rows).

thumbnail Fig. 7

Least-squares 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 6k-resolution 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 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.1A.3).

For the spherical harmonics expansion of the Gaussian surface density, we select lmin = 0. The value of lmax 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 lmax ≤ 8. We utilize lmax = 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 ground-based lightcurves. For the Gaia data, we set the uncertainty at σɛ = 0.01.

The results are based on the utilization of 2500 virtual-least-squares 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 rms-value of the lightcurve. A hybrid proposal probability density with equal weights for the virtual-observation 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 least-squares solutions. Altogether 30 000 MCMC sample solutions have been generated for each asteroid and unbiased sub-samples 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 virtual-observation 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. 24 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 minimum-maximum envelope following from the samples. Overall, the results are consistent.

thumbnail Fig. 8

As in Fig. 5 for asteroid (26) Proserpina.

thumbnail Fig. 9

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

4.2.1 (21) Lutetia

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 a20 and b21 (derived from the complex-valued coefficients s20 and s21), 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 a20 and b21 show substantial variegation that is nevertheless realistic: the 2500 virtual least-squares 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 a20 and b21 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 modelto 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; 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.

thumbnail Fig. 10

As in Fig. 6 for asteroid (26) Proserpina.

Table 4

Asin Table 3 for (26) Proserpina.

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 ground-based lightcurve originating from the Lowell Observatory photometricdatabase (Bowell et al. 2014) includes observations also at small phase angles, we have utilized the full H, G1, G2 phase function in the modeling of the photometric phase curve. The value for the slope βS then derives from the fullH, G1, G2 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 two-dimensional 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 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 spherical harmonics coefficients a20 and b21 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.

thumbnail Fig. 11

As in Fig. 5 for asteroid (585) Bilkis.

thumbnail Fig. 12

As in Fig. 9 for asteroid (585) Bilkis.

Table 5

Asin Table 3 for (585) Bilkis.

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 a20 and b21 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.

thumbnail Fig. 13

As in Fig. 6 for asteroid (585) Bilkis.

5 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 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-micro- arcsecond stellar and asteroid astrometry by the Gaia mission is likely to result in accelerated ground-based 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 ground-based polarimetric observations, they will allow for physics-based 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

Table A.1

Lightcurve characteristics for (21) Lutetia.

Table A.2

As in Table A.1 for (26) Proserpina.

Table A.3

As in Table A.1 for (585) Bilkis.

References

  1. Bartczak, P., & Dudziński, G. 2018, MNRAS, 473, 5050 [NASA ADS] [CrossRef] [Google Scholar]
  2. Barucci, M. A., Cellino, A., De Sanctis, C., et al. 1992, A&A, 266, 385 [Google Scholar]
  3. Bowell, E., Oszkiewicz, D., Wasserman, L., et al. 2014, MAPS, 49, 95 [Google Scholar]
  4. Carbognani, A., Tanga, P., Cellino, A., et al. 2012, Planet. Space. Sci., 73, 80 [NASA ADS] [CrossRef] [Google Scholar]
  5. Carry, B., Kaasalainen, M., Leyrat, C., et al. 2010, A&A, 523, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Carry, B., Kaasalainen, M., Merline, W., et al. 2012, Planet. Space. Sci., 66, 200 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Cellino, A., Zappalà, V., & Farinella, P. 1989, Icarus, 78, 298 [NASA ADS] [CrossRef] [Google Scholar]
  8. Cellino, A., Hestroffer, D., Tanga, P., Mottola, S., & Dell’Oro, A. 2009, A&A, 506, 935 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Cellino, A., Muinonen, K., Hestroffer, D., & Carbognani, A. 2015, Planet. Space. Sci., 118, 221 [NASA ADS] [CrossRef] [Google Scholar]
  10. Cellino, A., Hestroffer, D., Lu, X.-P., Muinonen, K., & Tanga, P. 2019, A&A, 631, A67 [CrossRef] [EDP Sciences] [Google Scholar]
  11. Chandrasekhar, S. 1960, Radiative Transfer (New York: Dover) [Google Scholar]
  12. Chang, Y., & Chang, C. 1963, Acta Astron. Sin., 11, 139 [Google Scholar]
  13. 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]
  14. Denchev, P., Magnusson, P., & Donchev, Z. 1998, Planet. Space. Sci., 46, 673 [NASA ADS] [CrossRef] [Google Scholar]
  15. Dotto, E., Barucci, M., Fulchignoni, M., et al. 1992, A&AS, 95, 195 [Google Scholar]
  16. Ďurech, J., & Hanuš, J. 2018, A&A, 620, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Ď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]
  18. Ďurech, J., Sidorin, V., & Kaasalainen, M. 2010, A&A, 513, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Farnham, T. 2013, Shape model of asteroid 21 Lutetia, NASA Planetary Data System, rO-A-OSINAC/OSIWAC-5-LUTETIA-SHAPE-V1.0 [Google Scholar]
  20. Gaia Collaboration (Spoto, F., et al.) 2018a, A&A, 616, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Gaia Collaboration (Brown, A. G. A., et al.) 2018b, A&A, 616, A1 [Google Scholar]
  22. Hanuš, J., Durech, J., Oszkiewicz, D. A., et al. 2016, A&A, 586, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Kaasalainen, M., & Torppa, J. 2001, Icarus, 153, 24 [Google Scholar]
  24. Kaasalainen, M., & Ďurech, J. 2007, Proc. IAU Symp. 236, 151 [Google Scholar]
  25. Kaasalainen, M., Lamberg, L., & Lumme, K. 1992a, A&A, 259, 333 [NASA ADS] [Google Scholar]
  26. Kaasalainen, M., Lamberg, L., Lumme, K., & Bowell, E. 1992b, A&A, 259, 318 [Google Scholar]
  27. Kaasalainen, M., Torppa, J., & Muinonen, K. 2001, Icarus, 153, 37 [NASA ADS] [CrossRef] [Google Scholar]
  28. 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]
  29. Lagerkvist, C.-I., Erikson, A., Debehogne, H., et al. 1995, A&AS, 113, 115 [NASA ADS] [Google Scholar]
  30. Lamberg, L. 1993, PhD thesis, University of Helsinki, Finland [Google Scholar]
  31. Lumme, K., & Bowell, E. 1981, AJ, 86, 1694 [NASA ADS] [CrossRef] [Google Scholar]
  32. Lupishko, D. F., Belskaya, I. N., & Tupieva, F. A. 1983, Sov. Astron. Lett., 9, 691 [Google Scholar]
  33. Lupishko, D., Velichko, F., Belskaya, I., & Shevchenko, V. 1987, Kinem. Fiz. Nebesn. Tel., 3, 36 [Google Scholar]
  34. Minkowski, H. 1903, Math. Ann., 57, 447 [CrossRef] [Google Scholar]
  35. Muinonen, K., & Bowell, E. 1993, Icarus, 104, 255 [NASA ADS] [CrossRef] [Google Scholar]
  36. Muinonen, K., & Lumme, K. 2015, A&A, 584, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Muinonen, K., Belskaya, I., Cellino, A., et al. 2010, Icarus, 209, 542 [NASA ADS] [CrossRef] [Google Scholar]
  38. Muinonen, K., Granvik, M., Oszkiewicz, D., Pieniluoma, T., & Pentikäinen, H. 2012, Planet. Space. Sci., 73, 15 [NASA ADS] [CrossRef] [Google Scholar]
  39. Muinonen, K., Wilkman, O., Cellino, A., Wang, X., & Wang, Y. 2015, Planet. Space. Sci., 118, 227 [NASA ADS] [CrossRef] [Google Scholar]
  40. Muinonen, K., Markkanen, J., Väisänen, T., Peltoniemi, J., & Penttilä, A. 2018, Opt. Lett., 43, 683 [NASA ADS] [CrossRef] [Google Scholar]
  41. 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]
  42. Oszkiewicz, D., Muinonen, K., Virtanen, J., & Granvik, M. 2010, Metroid. Planet. Sci., 44, 1897 [NASA ADS] [CrossRef] [Google Scholar]
  43. Penttilä, A., Shevchenko, V., Wilkman, O., & Muinonen, K. 2016, Planet. Space. Sci., 123, 117 [NASA ADS] [CrossRef] [Google Scholar]
  44. 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]
  45. Russell, H. N. 1906, ApJ, 24, 1 [NASA ADS] [CrossRef] [Google Scholar]
  46. Santana-Ros, T., Bartczak, P., Michałowski, T., Tanga, P., & Cellino, A. 2015, MNRAS, 450, 333 [NASA ADS] [CrossRef] [Google Scholar]
  47. Scaltriti, F., & Zappalà, V. 1979, Icarus, 39, 124 [CrossRef] [Google Scholar]
  48. Shevchenko, V., Belskaya, I., Muinonen, K., et al. 2016, Planet. Space. Sci., 123, 101 [NASA ADS] [CrossRef] [Google Scholar]
  49. Torppa, J., Kaasalainen, M., Michałowski, T., et al. 2003, Icarus, 164, 346 [NASA ADS] [CrossRef] [Google Scholar]
  50. Torppa, J., Hentunen, V.-P., Pääkkönen, P., Kehusmaa, P., & Muinonen, K. 2008, Icarus, 198, 91 [NASA ADS] [CrossRef] [Google Scholar]
  51. Torppa, J., Granvik, M., Penttilä, A., et al. 2018, Adv. Space Res., 62, 464 [NASA ADS] [CrossRef] [Google Scholar]
  52. Väisänen, T., Markkanen, J., Penttilä, A., & Muinonen, K. 2019, PLoS One, 14, 1 [CrossRef] [Google Scholar]
  53. Vernazza, P., Jorda, L., Sevecek, P., et al. 2020, Nat. Astron., 4, 136 [CrossRef] [Google Scholar]
  54. Viikinkoski, M., Kaasalainen, M., & Durech, J. 2015, A&A, 576, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Wang, X., Muinonen, K., & Wang, Y. 2015a, Planet. Space. Sci., 118, 242 [CrossRef] [Google Scholar]
  56. Wang, X., Muinonen, K., Wang, Y., et al. 2015b, A&A, 581, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Wang, A., Wang, X.-B., Muinonen, K., Han, X., & Wang, Y.-B. 2016, Res. A&A, 16 [Google Scholar]
  58. Wilkman, O., Muinonen, K., & Peltoniemi, J. 2015, Planet. Space. Sci., 118, 250 [CrossRef] [Google Scholar]
  59. Zappala, V., di Martino, M., Knezevic, Z., & Djurasevic, G. 1984, A&A, 130, 208 [NASA ADS] [Google Scholar]

All Tables

Table 1

Photometric G1 and G2 parameters for different asteroid taxonomic classes in the single-parameter phase function by Penttilä et al. (2016) based on the photometric observations reported in Shevchenko et al. (2016).

Table 2

Lightcurve characteristics for asteroids (21) Lutetia, (26) Proserpina, and (585) Bilkis incorporated in the present study.

Table 3

Rotation period P (h), pole longitude λ (°) and latitude β (°), and phase-curve slope parameter βS (mag rad−1) for (21) Lutetia.

Table 4

Asin Table 3 for (26) Proserpina.

Table 5

Asin Table 3 for (585) Bilkis.

Table A.1

Lightcurve characteristics for (21) Lutetia.

Table A.2

As in Table A.1 for (26) Proserpina.

Table A.3

As in Table A.1 for (585) Bilkis.

All Figures

thumbnail 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 (dash-dotted 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
thumbnail Fig. 2

Example photometric lightcurves for asteroid (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.

In the text
thumbnail Fig. 3

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

In the text
thumbnail Fig. 4

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

In the text
thumbnail Fig. 5

Markov-chain Monte Carlo sampling of the rotational, shape, and photometric phase-curve parameters for asteroid (21) Lutetia. The marginal probabilitity densities are characterized by an unbiased sub-sample 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 a20 and b21 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
thumbnail Fig. 6

Quick-look least-squares 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
thumbnail Fig. 7

Least-squares 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 6k-resolution shape model (bottom two rows) by Farnham (2013).

In the text
thumbnail Fig. 8

As in Fig. 5 for asteroid (26) Proserpina.

In the text
thumbnail Fig. 9

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

In the text
thumbnail Fig. 10

As in Fig. 6 for asteroid (26) Proserpina.

In the text
thumbnail Fig. 11

As in Fig. 5 for asteroid (585) Bilkis.

In the text
thumbnail Fig. 12

As in Fig. 9 for asteroid (585) Bilkis.

In the text
thumbnail Fig. 13

As in Fig. 6 for asteroid (585) Bilkis.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.