© K. Muinonen et al. 2020

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 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; 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 ${(\lambda ,\beta )}^T$ (J2000.0, T stands for transpose), and the rotational phase at a given epoch t₀ by φ₀. 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₀ (including π in accordance with a common definition) and the emergent intensity I as $$\begin{array}{l}I(\mu ,\varphi ;{\mu }_0,{\varphi }_0)={\mu }_{0}R(\mu ,\varphi ;{\mu }_0,{\varphi }_0)F_0,\\ {\mu }_0=\cos \iota ,\mu =\cos \varepsilon ,\end{array}$$(1)

where ι and ɛ are the angles of incidence and emergence as measured from the outward normal vector of the element, and ϕ₀ 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 ϕ₀ is unnecessary.

The Lommel-Seeliger reflection coefficient (subscript LS) is (e.g., Lumme & Bowell 1981; Wilkman et al. 2015), $$\begin{array}{lll}{R}_{\text{LS}}(\mu ,{\mu }_0,\varphi )\hfill & =\hfill & \frac{1}{4}\tilde{\omega }{P}_{11}(\alpha )\frac{1}{\mu +{\mu }_0},\hfill \end{array}$$(2)

where $\tilde{\omega }$ is the single-scattering albedo (proportion of incident light scattered, $0\le \tilde{\omega }\le 1$), P₁₁ is the single-scattering phase function (from the 4 × 4 scattering phase matrix), and α is the phase angle. We note that α derives unequivocally from ι, ɛ, and ϕ. R_LS, the first-order multiple-scattering approximation of the radiative transfer equation for a semi-infinite plane-parallel medium (e.g., Chandrasekhar 1960), is applicable to weakly scattering media: this means that the intensity terms $[{\tilde{\omega }}^k]$, k ≥ 2 are assumed negligible. When omitting polarization effects, the scattering phase function P₁₁ provides the angular distribution of scattered light in an individual interaction and is normalized so that $$\begin{array}{lll}{\displaystyle {\int }_{(4\pi )}\frac{\text{d}\Omega }{4\pi }}{P}_{11}(\alpha )\hfill & =\hfill & 1.\hfill \end{array}$$(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 $$\begin{array}{lll}L(\alpha )\hfill & =\hfill & {\displaystyle {\int }_{A_{+}}\text{d}}A\mu I(\mu ,{\mu }_0,\alpha )\hfill \\ \hfill & =\hfill & {\displaystyle {\int }_{A_{+}}\text{d}}A\mu {\mu }_{0}R(\mu ,{\mu }_0,\alpha )F_0,\hfill \end{array}$$(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, $$\begin{array}{lll}L(\alpha )\hfill & =\hfill & \frac{1}{4}{D}^2{\displaystyle {\int }_{{\Omega }_{+}}\text{d}}\Omega \mu {\mu }_{0}R(\mu ,{\mu }_0,\alpha )F_0,\hfill \end{array}$$(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, $$\begin{array}{l}L(\alpha )=\frac{1}{32}\pi F_0{D}^2\tilde{\omega }{P}_{11}(\alpha ){\Phi }_{\text{LS}}(\alpha ),\\ {\Phi }_{\text{LS}}(\alpha )=1-\sin \frac{1}{2}\alpha \tan \frac{1}{2}\alpha \ln \left(\cot \frac{1}{4}\alpha \right),\end{array}$$(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°: $$\begin{array}{lll}p\hfill & =\hfill & 4\frac{L(0)}{D^2\pi F_0}.\hfill \end{array}$$(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 $\tilde{\omega }$ and phase function P₁₁(α) in the Lommel-Seeliger reflection coefficient R_LS (Eq. (2)). On one hand, the disk-integrated brightness phase curve of an asteroid in a given apparition can be described by the H, G₁, G₂ phase function ${\Phi }_{H{G}_1{G}_2}$ (Muinonen et al. 2010; Penttilä et al. 2016). On the other hand, the reflection coefficient R_LS in Eq. (2) includes the ${(\mu +{\mu }_0)}^{-1}$ term that produces the phase function Φ_LS given in Eq. (6). Thus, modeling P₁₁(α) so that the phase curve from Eq. (2) for a spherical asteroid is precisely that of the H, G₁, G₂ 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), $$\begin{array}{lll}p\hfill & =\hfill & \frac{1}{8}\tilde{\omega }{P}_{11}(0^{°}).\hfill \end{array}$$(8)

Thus, it is reasonable to assume that the model $$\begin{array}{lll}\frac{1}{8}\tilde{\omega }{P}_{11}(\alpha )\hfill & =\hfill & p\frac{{\Phi }_{H{G}_1{G}_2}(\alpha )}{{\Phi }_{\text{LS}}(\alpha )},\hfill \end{array}$$(9)

resulting in $$\begin{array}{lll}{R}_{\text{LS}}(\mu ,{\mu }_0,\varphi )\hfill & =\hfill & 2p\frac{{\Phi }_{H{G}_1{G}_2}(\alpha )}{{\Phi }_{\text{LS}}(\alpha )}\frac{1}{\mu +{\mu }_0},\hfill \end{array}$$(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, G₁, G₂ 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 $$\begin{array}{l}V(\alpha )=H-[2.5\lg {\Phi }_{H{G}_1{G}_2}(\alpha )-5\lg (r\Delta )]\text{mag},\\ H=2.5(6.259-2\lg D-\lg p)\text{mag},\end{array}$$(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): $$\begin{array}{l}{\Phi }_{H{G}_1{G}_2}(\alpha )\approx {10}^{-0.4(m_{20}+{\beta }_{\text{S}}\Delta \alpha )},{20}^{°}\le \alpha \le {50}^{°},\\ m_{20}=-2.5\lg {\Phi }_{H{G}_1{G}_2}({20}^{°}),\\ \Delta \alpha =\alpha -{20}^{°}.\end{array}$$(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 G₁, G₂ 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⁻², b⁻², c⁻²) (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 $$\begin{array}{lll}L(\alpha )\hfill & =\hfill & \frac{1}{8}\pi F_0\tilde{\omega }{P}_{11}(\alpha )abc\{\frac{(S_{\odot }+S_{\oplus })(S_{\odot }{S}_{\oplus }+s_{\odot \oplus })}{S}\hfill \\ \hfill & \hfill & -\frac{S_{\odot }^2{S}_{\oplus }^2-s_{\odot \oplus }^2}{S^3}\ln \left|\frac{S+S_{\odot }+S_{\oplus }}{S-S_{\odot }-S_{\oplus }}\right|\},\hfill \end{array}$$(13)

where $$\begin{array}{l}S=\sqrt{S_{\odot }^2+S_{\oplus }^2+2s_{\odot \oplus }},\\ S_{\odot }=\sqrt{e_{\odot }^{\text{T}}C{e}_{\odot }},\\ S_{\oplus }=\sqrt{e_{\oplus }^{\text{T}}C{e}_{\oplus }},\\ s_{\odot \oplus }=e_{\odot }^{\text{T}}C{e}_{\oplus },\end{array}$$(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 $$\begin{array}{lll}G(\theta ,\phi )\hfill & =\hfill & \exp \left[{\sum }_{l=l_{\text{min}}}^{l_{\text{max}}}{\sum }_{m=-l}^l{s}_{lm}{Y}_{lm}(\theta ,\phi )\right],\hfill \end{array}$$(15)

where the Y_lm’s are the orthogonal spherical harmonics and the complex-valued coefficients s_lm ensure a real-valued surface density (the superscript * refers to complex conjugation): $$\begin{array}{l}{s}_{l,-m}={(-1)}^m{s}_{lm}^{\ast },\text{Im}(s_{l0})=0,\\ l=0,\dots ,l_{\text{max}},m=-l,\dots ,0,\dots ,l.\end{array}$$(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 ${(l_{\text{max}}+1)}^2-l_{\text{min}}^2$ independent real or imaginary parts of the coefficients s_lm. When utilizing real-valued spherical harmonics, we utilize well-behaving normalized associated Legendre functions. This entails multiplication of the normal associated Legendre functions by $\sqrt{(2l+1)/2}\sqrt{(l-m)!/(l+m)!}$.

The Gaussian surface density gives the relative proportion of surface normal vectors in the given directions specified by the normal coordinates. The convex shape can be derived from the Gaussian surface density unambiguously via the solution of the Minkowski problem (Kaasalainen et al. 1992a,b; Lamberg 1993), see Sect. 3.1.

We have chosen to use the pure Lommel-Seeliger reflection coefficient (Eq. (2)) instead of the frequently-used 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.

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.

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₀ by the vector P = (P, λ, β, φ₀, s₀₀, …, $s_{l_{\text{max}}{l}_{\text{max}}},p,G_1,G_2,D{)}^T$, where the parameters are, respectively, the rotation period, ecliptic pole longitude, ecliptic pole latitude, rotational phase at t₀, shape parameters $s_{00},\dots ,s_{l_{\text{max}}{l}_{\text{max}}}$, geometric albedo, two parameters of the H, G₁, G₂ 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 $$\begin{array}{lll}{N}_{\text{P}}\hfill & =\hfill & {(l_{\text{max}}+1)}^2-l_{\text{min}}^2+6.\hfill \end{array}$$(17)

Here the rotational phase φ₀ is strictly redundant, as the convex shape model will accommodate for any change in φ₀. 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, λ, β, φ₀, a, b, c, p, G₁, G₂, 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, G₁₂ phase function (Muinonen et al. 2010), including G₁₂ among the parameters. Furthermore, we set l_min = 0. In this case,for general convex shapes, $$\begin{array}{lll}{N}_{\text{P}}\hfill & =\hfill & {(l_{\text{max}}+1)}^2+5,\hfill \end{array}$$(18)

and our parameters are $P={(P,\lambda ,\beta ,{\phi }_0,s_{00},\dots ,s_{l_{\text{max}}{l}_{\text{max}}},G_{12})}^T$. For triaxial ellipsoid shapes, we have $P={(P,\lambda ,\beta ,{\phi }_0,b,c,G_{12})}^T$. We may also utilize the single-parameter phase function provided by Penttilä et al. (2016): in this case, the G₁, G₂ 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₁₂ 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, $$\begin{array}{lll}N\hfill & =\hfill & {\sum }_{k=1}^K{N}_k,\hfill \end{array}$$(19)

where N_k specifies the number of observations in lightcurve k (k = 1, …, K). Let $$\begin{array}{l}{M}_{\text{obs}}=(M_{\text{obs},11},\dots ,M_{\text{obs},1N_1};\dots ;\\ M_{\text{obs},K1},\dots ,M_{\text{obs},K{N}_K}{)}^{\text{T}},\\ M(P)=(M_{11}(P),\dots ,M_{1N_1}(P);\dots ;\\ M_{K1}(P),\dots ,M_{K{N}_K}(P){)}^{\text{T}}\end{array}$$(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: $$\begin{array}{lll}{M}_{\text{obs}}\hfill & =\hfill & M(P)+\varepsilon +\upsilon .\hfill \end{array}$$(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 $$\begin{array}{lll}{\Lambda }_{\varepsilon +\upsilon }\hfill & =\hfill & {\Lambda }_{\varepsilon }+{\Lambda }_{\upsilon }.\hfill \end{array}$$(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 p_p be the a posteriori probability density function (p.d.f.) for the parameters. Within the Bayesian framework (cf. Muinonen & Bowell 1993), p_p is proportional to the a priori and observational uncertainty p.d.f.s p_pr and p_ɛ+υ, the latter being evaluated for the “Observed-Computed” (O-C) residual magnitudes ΔM(P), $$\begin{array}{l}{p}_{\text{p}}(P)\propto p_{\text{pr}}(P)p_{\varepsilon +\upsilon }(\Delta M(P)),\\ \Delta M(P)=M_{\text{obs}}-M(P).\end{array}$$(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 $$\begin{array}{lll}{p}_{\text{p}}(P)\hfill & \propto \hfill & p_{\text{pr}}(P)\exp \left[-\frac{1}{2}{\chi }^2(P)\right],\hfill \\ {\chi }^2(P)\hfill & =\hfill & \Delta M^T(P){\Lambda }_{\varepsilon +\upsilon }^{-1}\Delta M(P).\hfill \end{array}$$(24)

where χ² 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 χ²(P) as $$\begin{array}{l}{\chi }^2(P)={\sum }_{k=1}^K\Delta M_k^T(P){\Lambda }_{\varepsilon +\upsilon ,k}^{-1}\Delta M_k(P),\\ \Delta M_k(P)=M_{\text{obs},k}-M_k(P),\end{array}$$(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): $$\begin{array}{l}{M}_{\text{obs},k0}=\frac{1}{N_k}{\sum }_{j=1}^{N_k}{M}_{\text{obs},kj},\\ M_{k0}(P)=\frac{1}{N_k}{\sum }_{j=1}^{N_k}{M}_{kj}(P),\\ \Delta M_{k0}(P)=M_{\text{obs},k0}-M_{k0}(P),\\ k=1,2,\dots ,K.\end{array}$$(26)

Furthermore, define the magnitude differences m_obs,kj and m_kj(P) with the help of Eq. (26): $$\begin{array}{lll}{m}_{\text{obs},kj}\hfill & =\hfill & M_{\text{obs},kj}-M_{\text{obs},k0},\hfill \\ m_{kj}(P)\hfill & =\hfill & M_{kj}(P)-M_{k0}(P).\hfill \end{array}$$(27)

Introduce then the relative brightnesses ℓ_obs,kj and ℓ_kj(P): $$\begin{array}{lll}{\mathcal{l}}_{\text{obs},kj}\hfill & =\hfill & {10}^{-0.4m_{\text{obs},kj}},\hfill \\ {\mathcal{l}}_{kj}(P)\hfill & =\hfill & {10}^{-0.4m_{kj}(P)}.\hfill \end{array}$$(28)

In the case of large uncertainties υ, that is, in the case of what can be called relative photometry, $$\begin{array}{lll}{\Lambda }_{\varepsilon +\upsilon }\hfill & \approx \hfill & {\Lambda }_{\upsilon }.\hfill \end{array}$$(29)

Recalling that Λ_υ includes high correlations among the uncertainties for the observations of a given lightcurve, $$\begin{array}{lll}{\chi }^2(P)\hfill & =\hfill & {\sum }_{k=1}^K\Delta M_k^T(P){\Lambda }_{\upsilon ,k}^{-1}\Delta M_k(P)\hfill \\ \hfill & \approx \hfill & {\sum }_{k=1}^K\frac{1}{N_k{\sigma }_{\upsilon ,k}^2}{\sum }_{j=1}^{N_k}{\left[M_{\text{obs},kj}-M_{kj}(P)\right]}^2.\hfill \end{array}$$(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): $$\begin{array}{lll}{\chi }^2(P)\hfill & =\hfill & {\sum }_{k=1}^K\frac{1}{{\sigma }_{\upsilon ,k}^2}\{\Delta M_{k0}^2(P)+\frac{1}{N_k}{\sum }_{j=1}^{N_k}{\left[m_{\text{obs},kj}-m_{kj}(P)\right]}^2\hfill \\ \hfill & \hfill & -2\Delta M_{k0}(P){\sum }_{j=1}^{N_k}\left[m_{\text{obs},kj}-m_{kj}(P)\right]\}\hfill \\ \hfill & \to \hfill & {\sum }_{k=1}^K\frac{1}{N_k{\sigma }_{\upsilon ,k}^2}{\sum }_{j=1}^{N_k}{\left[m_{\text{obs},kj}-m_{kj}(P)\right]}^2\hfill \\ \hfill & \to \hfill & \frac{1}{{\sigma }_{\upsilon }^2}{\sum }_{k=1}^K\frac{1}{N_k}{\sum }_{j=1}^{N_k}{\left[m_{\text{obs},kj}-m_{kj}(P)\right]}^2.\hfill \end{array}$$(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 ΔM_k0(P) → 0. In other words, we can introduce a systematic correction (or, formally, an additional parameter) for each lightcurve to exactly result in ΔM_k0(P) = 0. Furthermore, we have assumed an equal standard deviation σ_υ for all lightcurves.

Equation (31) shows that the least-squares optimization for relative photometry can conveniently be carried out by giving an equal weight for each lightcurve (assuming σ_υ,k ≡ σ_υ for all k). However, the equation also shows that the probabilistic treatment runs into difficulties, since, with the missing information about σ_υ, there is no meaningful weighting scheme for χ²(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, $$\begin{array}{lll}{\Lambda }_{\varepsilon +\upsilon }\hfill & \approx \hfill & {\Lambda }_{\varepsilon }.\hfill \end{array}$$(32)

For a diagonal matrix Λ_ɛ, $$\begin{array}{lll}{\chi }^2(P)\hfill & =\hfill & {\sum }_{k=1}^K{\sum }_{j=1}^{N_k}\frac{1}{{\sigma }_{\varepsilon ,kj}^2}{\left[M_{\text{obs},kj}-M_{kj}(P)\right]}^2.\hfill \end{array}$$(33)

If the uncertainties ɛ are small, the O-C magnitude difference can be linearized into a difference in the corresponding brightnesses L_obs,kj and L_kj (P): $$\begin{array}{lll}{M}_{\text{obs},kj}-M_{kj}(P)\hfill & =\hfill & 2.5\lg \frac{L_{kj}(P)}{L_{\text{obs},kj}}\hfill \\ \hfill & \approx \hfill & -(2.5\lg \text{e})\left[\frac{L_{\text{obs},kj}-L_{kj}(P)}{L_{\text{obs},kj}}\right],\hfill \end{array}$$(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 $$\begin{array}{lll}{\chi }^2(P)\hfill & =\hfill & {(2.5\lg \text{e})}^2{\sum }_{k=1}^K{\sum }_{j=1}^{N_k}\frac{1}{{\sigma }_{\varepsilon ,kj}^2}{\left[\frac{L_{\text{obs},kj}-L_{kj}(P)}{L_{\text{obs},kj}}\right]}^2\hfill \\ \hfill & =\hfill & {(2.5\lg \text{e})}^2{\sum }_{k=1}^K{\sum }_{j=1}^{N_k}\frac{1}{{\sigma }_{\varepsilon ,kj}^2}{\left[\frac{{\mathcal{l}}_{\text{obs},kj}-{\mathcal{l}}_{kj}(P){10}^{0.4\Delta M_{k0}(P)}}{{\mathcal{l}}_{\text{obs},kj}}\right]}^2\hfill \\ \hfill & =\hfill & {\sum }_{k=1}^K\frac{1}{{\sigma }_{\varepsilon ,k}^2}{\sum }_{j=1}^{N_k}{\left[\frac{{\mathcal{l}}_{\text{obs},kj}-{\mathcal{l}}_{kj}(P){10}^{0.4\Delta M_{k0}(P)}}{{\mathcal{l}}_{\text{obs},kj}}\right]}^2,\hfill \end{array}$$(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 χ²-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 χ² 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 $$\begin{array}{l}{\sigma }_{\varepsilon ,k}=\sqrt{\frac{N_k\Delta T_k}{\Delta T_{\tilde{k}}}}\max ({\sigma }_{0,k},{\sigma }_{\text{pr},k}),\\ k=1,\dots ,K,\end{array}$$(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, ΔT_k is the mean sampling time interval of the lightcurve k and $$\begin{array}{lll}\Delta T_{\tilde{k}}\hfill & =\hfill & {\max }_{k=1,\dots ,K}\Delta T_k,\hfill \\ \Delta T_k\hfill & =\hfill & \frac{T_k}{N_k-1},\hfill \end{array}$$(37)

$\tilde{k}$ marking the lightcurve with the lowest sampling rate in time and T_k standing for the time span of the lightcurve k. Equation (36) can be interpreted in the following way. The lightcurve $\tilde{k}$ contributes a χ² -value of unity in the inverse problem, whereas each other lightcurve k is considered to split into approximately $$\begin{array}{lll}{N}_{k,\text{eff}}\hfill & =\hfill & \frac{\Delta T_{\tilde{k}}}{\Delta T_k}\hfill \end{array}$$(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 $$\begin{array}{l}{\sigma }_{\varepsilon ,k}=\sqrt{\frac{N_k}{N_{\tilde{k}}}}\max ({\sigma }_{0,k},{\sigma }_{\text{pr},k}),\\ k=1,\dots ,K,\end{array}$$(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 $\tilde{k}$ corresponds to the one with the minimum number of observations and contributes a χ² -value equal to $N_{\tilde{k}}$ (when σ_0,k≥ σ_pr,k), whereas each other lightcurve k is scaled to result in the same χ²-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 χ² takes the form $$\begin{array}{lll}{\chi }^2(P)\hfill & =\hfill & {\sum }_{k=1}^{K^{\prime }}\frac{1}{{{\sigma }^{\prime }}_{\varepsilon ,k}^2}{\sum }_{j=1}^{{N^{\prime }}_k}{\left[\frac{{{\mathcal{l}}^{\prime }}_{\text{obs},kj}-{{\mathcal{l}}^{\prime }}_{kj}(P){10}^{0.4\Delta {M^{\prime }}_{k0}(P)}}{{{\mathcal{l}}^{\prime }}_{\text{obs},kj}}\right]}^2\hfill \\ \hfill & \hfill & +{\sum }_{k=1}^{K^{″}}\frac{1}{{{\sigma }^{″}}_{\varepsilon ,k}^2}{\sum }_{j=1}^{{N^{″}}_k}{\left[\frac{{{\mathcal{l}}^{″}}_{\text{obs},kj}-{{\mathcal{l}}^{″}}_{kj}(P){10}^{0.4\Delta {M^{″}}_{k0}(P)}}{{{\mathcal{l}}^{″}}_{\text{obs},kj}}\right]}^2,\hfill \end{array}$$(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 ${\chi }_{\text{R}}^2(P)$ that depends on the shape parameters $s_{00},\dots ,s_{l_{\text{max}}{l}_{\text{max}}}$ included in P: $$\begin{array}{lll}{\chi }_{\text{R}}^2(P)\hfill & =\hfill & \frac{1}{{\sigma }_{\text{reg}}^2}{\left[{\sum }_{j=1}^{N_{△}}{A}_j(P)n_j(P)\right]}^2,\hfill \end{array}$$(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 ${\chi }_{\text{R}}^2(P)$-value of Eq. (41) should vanish. Here the triangular discretization needs to have a resolution high enough to warrant an accurate computation of the disk-integrated brightnesses. In the spirit of the Bayesian inference, σ_reg describes the standard deviation for the regularization term. We return to specifying the actual value of σ_reg in Sect. 4.2.

Finally, when considering sparse observational data sets, like the GDR2 photometric data sets, the convex inversion is additionally regularized by constraining the spectral power in the spherical harmonics expansion for the Gaussian surface density. After experimenting with the present asteroid examples, we require that the spectral power in each degree of spherical harmonics must not exceed a value of 0.5. In the present examples, this does not substantially limit the resulting spectra of shapes, but prevents the MCMC sampling from wandering into unrealistic parameter regimes.

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 a_r : $$\begin{array}{lll}{a}_r\hfill & =\hfill & \frac{p_{\text{p}}(P^{\prime })p_{\text{t}}(P_j;P^{\prime })}{p_{\text{p}}(P_j)p_{\text{t}}(P^{\prime };P_j)}.\hfill \end{array}$$(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]: $$\begin{array}{lll}{P}_{j+1}\hfill & =\hfill & \{\begin{array}{ll}{P}^{\prime },\hfill & y\le a_r,\hfill \\ P_j,\hfill & y>a_r,\hfill \end{array}\hfill \end{array}$$(43)

that is, the proposed parameters are accepted with the probability of min (1, a_r). 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. p_p. For convergence monitoring, there are various diagnostics tools available (see, e.g., Oszkiewicz et al. 2010).

In the virtual-observation MCMC methods, earlier studied for asteroid orbital inversion (Muinonen et al. 2012) and lightcurve inversion with Lommel-Seeliger ellipsoids (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, $$\begin{array}{lll}{M}_{\text{v}}\hfill & =\hfill & M_{\text{obs}}+{\varepsilon }_{\text{v}}.\hfill \end{array}$$(44)

Thereafter, virtual least-squares parameters P_v are obtained by minimizing $$\begin{array}{lll}{\chi }_{\text{v}}^2(P)\hfill & =\hfill & {\left[M_{\text{v}}-M(P)\right]}^{\text{T}}{({\Lambda }_{\varepsilon }+{\Lambda }_{\text{v}})}^{-1}\left[M_{\text{v}}-M(P)\right],\hfill \end{array}$$(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 N-dimensional integral $$\begin{array}{lll}{p}_{\text{V}}(P_{\text{v}})\hfill & =\hfill & {\displaystyle \int \text{d}}{M}_{\text{v}}\delta (P_{\text{v}}-P_{\text{v}}(M_{\text{v}}))p(M_{\text{v}}),\hfill \end{array}$$(46)

where p(M_v) is the Gaussian p.d.f. for the virtual observations and δ is Dirac’s Delta function.

A random-walk MCMC-sampler is obtained with the help of the symmetric proposal p.d.f., defined as the 2(N+ N_p)-dimensional integral $$\begin{array}{lll}{p}_{\text{t}}(\Delta P)\hfill & =\hfill & {\displaystyle \int {\displaystyle \int \text{d}}}{P}_{\text{v}}\text{d}{P}^{\prime }{}_{\text{v}}\delta \left(\Delta P-(P_{\text{v}}-P^{\prime }{}_{\text{v}})\right)p_{\text{V}}(P_{\text{v}})p_{\text{V}}(P^{\prime }{}_{\text{v}}),\hfill \\ \hfill & \hfill & \hfill \end{array}$$(47)

where the integration over N_p parameters results in a (2N + N_p)-dimensional convolution integral $$\begin{array}{lll}{p}_{\text{t}}(\Delta P)\hfill & =\hfill & {\displaystyle \int \text{d}}{P}_{\text{v}}{p}_{\text{V}}(P_{\text{v}})p_{\text{V}}(P_{\text{v}}-\Delta P).\hfill \end{array}$$(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 P_v before the MCMC sampling (cf. Wang et al. 2015a,b). The generation of virtual observation sets L_v and the derivation of P_v are repeated until there are N_v ≫ 1 sets of parameters $P_{\text{v}}^{(j)}$ (j = 1, 2, 3…, N_v). Parameter differences are then readily available from $$\begin{array}{lll}\Delta P_{jk}=P_{\text{v},j}-P_{\text{v},k},\hfill & \hfill & j,k=1,2,3\dots ,N_{\text{v}};j\ne k.\hfill \end{array}$$(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 $N_{\text{v}}^2$ differences available, that is, a number that can become extremely large for reasonable values of N_v. Any pair of virtual least-squares parameters consumed for a difference should be replaced by a new pair of parameters computed as described above, but, in practice, if N_v is large enough, the replacement can be unnecessary.

The ratio a_r for the decision criterion is expressed in the concise form $$\begin{array}{lll}{a}_r\hfill & =\hfill & \frac{p_{\text{p}}(P^{\prime })}{p_{\text{p}}(P_j)}.\hfill \end{array}$$(50)

The proposals can now be computed during the MCMC sampling using the N_v virtual least-squares element sets P_v,j.

The covariance matrix Λ_v can be scaled from Λ_ɛ. For dense lightcurves, it is also possible to generate the virtual observations from cubic spline fits to the lightcurves, in which case we may set Λ_v = Λ_ɛ.

The virtual least-squares solutions can be utilized to compute a covariance matrix for the parameters. It is straightforward to incorporate the covariance matrix in a multivariate Gaussian proposal probability density. Combining the virtual-observation and Gaussian proposals results in a hybrid proposal probability density that will be used in the MCMC sampling below.

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 χ² (P), $$\begin{array}{lll}\Delta {\chi }^2(P)={\chi }^2(P)-{\chi }^2(P_0)\hfill & \le \hfill & \Delta {\chi }_{\text{c}}^2,\hfill \end{array}$$(51)

where P₀ denotes the least-squares parameters and the cut-off $\Delta {\chi }_{\text{c}}^2$ is set on the basis of the number of parameters. For the present inverse problem, we may define $\Delta {\chi }_{\text{c}}^2$ with the help of probability thresholds for multivariate normal statistics: $$\begin{array}{lll}{\left[\text{Prob}(-\Delta {\chi }_{\text{c}}\le \Delta \chi \le \Delta {\chi }_{\text{c}})\right]}^{N_{\text{P}}}\hfill & =\hfill & P_{\text{c}},\hfill \end{array}$$(52)

where P_c denotes the cut-off probability value. It is recognized that Eq. (52) describes the probability mass for normal statistics only. Nevertheless, we consider it useful in approximately evaluating the statistical distance in our probability densities. For N_P ≈ 50, we may choose $$\begin{array}{lll}\Delta {\chi }_{\text{c}}^2\hfill & \approx \hfill & 100,\hfill \end{array}$$(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 $$\begin{array}{lll}{a}_r\hfill & =\hfill & \{\begin{array}{ll}0,\hfill & \Delta {\chi }^2\ge \Delta {\chi }_{\text{c}}^2,\hfill \\ \text{const}.\ne 0,\hfill & \Delta {\chi }^2<\Delta {\chi }_{\text{c}}^2,\hfill \end{array}\hfill \end{array}$$(54)

that is, transitions are always accepted within the given $\Delta {\chi }_{\text{c}}^2$ regime but never across the $\Delta {\chi }_{\text{c}}^2$ 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 $$\begin{array}{lll}{w}_j\hfill & \propto \hfill & p_{\text{p}}(P_j),j=1,2,3,\dots ,N.\hfill \end{array}$$(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 χ² 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 ${\chi }_{\text{R}}^2$ 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 (θ, φ)): $$\begin{array}{lll}\varrho (\theta ,\phi )\hfill & =\hfill & n(\theta ,\phi )\cdot r(\theta ,\phi ),\hfill \end{array}$$(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 $$\begin{array}{lll}V(\text{E},\text{F})\hfill & =\hfill & \frac{1}{3}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }d}}\theta d\phi \sin \theta {\varrho }_{\text{E}}(\theta ,\phi )G_{\text{F}}(\theta ,\phi ).\hfill \end{array}$$(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 $$\begin{array}{lll}V(E)\hfill & =\hfill & \frac{1}{3}{\sum }_{j=1}^{N_f}{l}_j{A}_j(l),\hfill \end{array}$$(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 $$\begin{array}{lll}V(E,F)\hfill & =\hfill & \frac{1}{3}{\sum }_{j=1}^{N_f}{l}_j(E)a_j(F),\hfill \end{array}$$(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 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 ${(x_i,y_i,z_i)}^T$ of its N_n nodes i = 1, …, N_n as the 3N_n free parameters. The derivation of the least-squares solution for the coordinates is carried out via Monte Carlo optimization. At each iteration, the convexity of the proposed shape is verified before computing the χ² -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 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. 2–4 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.1–A.3). For the first two asteroids, we have extracted the ground-based observations from the DAMIT database (Ďurech et al. 2010). Additionally, the data 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.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 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.

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.

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.

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.

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).

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.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 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. 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 minimum-maximum envelope following from the samples. Overall, the results are consistent.

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

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 a₂₀ and b₂₁ (derived from the complex-valued coefficients s₂₀ and s₂₁), 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₂₀ and b₂₁ 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⁻⁷ h ≈ 8.3 × 10⁻² s. The period uncertainty is small, comparable to uncertainties derived from space mission flybys.

The spherical harmonics coefficients a₂₀ and b₂₁ 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⁻¹. 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⁻¹, which is in agreement with the result with all the observations included.

Fig. 10 As in Fig. 6 for asteroid (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, G₁, G₂ phase function in the modeling of the photometric phase curve. The value for the slope β_S then derives from the fullH, G₁, G₂ 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⁻⁶ h ≈ 5.5 × 10⁻¹ s. The spherical harmonics coefficients a₂₀ and b₂₁ 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⁻¹ and β_S = (1.720 ± 0.109) mag rad⁻¹, respectively. The analysis of the Gaia data only with ellipsoid inversion gives β_S = (1.62 ± 0.15) mag rad⁻¹. The analysis of the Gaia data only with ellipsoid inversion and assuming the rotation parameters fixed gives β_S = (1.53 ± 0.15) mag rad⁻¹. Given their uncertainties, the four photometric slope values are in mutual agreement.

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

Fig. 12 As in Fig. 9 for asteroid (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⁻¹ s. Again, the spherical harmonics coefficients a₂₀ and b₂₁ 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⁻¹, whereas the analysis of the Gaia data only gives β_S = (2.152 ± 0.119) mag rad⁻¹. The analysis of the Gaia data only with ellipsoid inversion gives β_S = (2.52 ± 0.10) mag rad⁻¹. The analysis of the Gaia data only with ellipsoid inversion and assuming the rotation parameters fixed gives β_S = (2.43 ± 0.20) mag rad⁻¹. 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.

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

References

All Tables

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

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

	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

	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

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

	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

	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

	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