© ESO 2020

1 Introduction

Outgassing causes nongravitational accelerations that affect both the trajectories and rotational properties of minor bodies (see for example Sekanina 1981). The spin-axis orientation of a comet is defined by two angles (see for example Fig. 1 in Yeomans et al. 2004), the equatorial obliquity (ɛ_c) and the cometocentric longitude or argument of the subsolar meridian at perihelion (ϕ_c). The subsolar meridian is the great circle intersecting the object’s rotation axis and the subsolar point, where the Sun is directly overhead as seen from the minor body.

The value of ɛ_c gives the angle between the orbital and equatorial planes of the comet. Retrograde comets have ɛ_c ∈ (90°,  180°); if the direction of rotation is consistent with that of its orbital motion, ɛ_c < 90° and the comet is prograde (Sekanina 1981). The value of ϕ_c is measured from the vernal equinox of the comet (projection of the ascending node of the orbital plane of the comet on its equator) to the subsolar meridian at perihelion in the direction of increasing true anomaly, f. Therefore, it gives the hour angle of the Sun as seen from the comet at perihelion. When ϕ_c ∈ (180°,  360°), the sunlit pole is the southern one; if ϕ_c ∈ (0°,  180°), the northern pole faces the Sun at perihelion (Sekanina 1981).

In addition to the spin-axis orientation defined by (ɛ_c,  ϕ_c), the thermal lag angle (η_c) measures how the region that dominates outgassing is shifted behind the subsolar meridian (Yeomans et al. 2004), and its value depends on the thermophysical properties of the sublimation phenomena produced in the comet and its rotational period (Sekanina 1981). Maximum insolation takes place at comet’s noon. If η_c ∈ (−90°, 0°), sublimation takes place mainly when the Sun is rising as seen from the surface of the comet; when the peak of activity is shifted toward the afternoon and evening hours, as seen from the comet, η_c ∈ (0°, 90°). The interval η_c ∈ (−90°, 90°) corresponds to the comet’s dayside; outside this range, we have the comet’s nightside. Sublimation of volatile materials, such as H₂O or CO, generates nongravitational accelerations that can be detected by analyzing astrometric data.

The nongravitational force is the result of directional mass loss and, therefore, its detectability decreases as the rotation rate increases; the orbital evolution of fast cometary nuclei with rotation periods of a few hours is nearly unaffected by outgassing (see for example Samarasinha et al. 2004). Marsden et al. (1973) discussed a mathematical formalism that developed the model originally proposed by Whipple (1950) to quantify the nongravitational effects on the orbital evolution of comets. Vaporization of volatiles produces a force with radial, transverse, and normal components of the form A g(r) with $g(r)=\alpha {(r/r_0)}^{-m}{(1+{(r/r_0)}^n)}^{-k}$, where r is the heliocentric distance of the comet (Marsden et al. 1973). For subsurface water ice, r₀ = 2.808 AU, k = 4.6142, m = 2.15, n = 5.093, and α = 0.1112620426 (Marsden et al. 1973); for carbon monoxide ice, r₀ = 5.0 AU, k = 2.6, m = 2.0, n = 3.0, and α = 0.0408373333128795 (Yabushita 1996; Micheli et al. 2018). The values of the nongravitational parameters – A₁ (radial), A₂ (transverse), and A₃ (normal) – depend on those of the angles ɛ_c, ϕ_c, and η_c.

Assuming a spherically symmetric nucleus and applying a linear precession model, the nongravitational parameters can be computed from astrometric observations of the comet, and ɛ_c, ϕ_c, and η_c can be obtained from the components of the nongravitational force by using an iterative least squares algorithm, as shown by Sitarski (1990). However, many observed comets and active asteroids are known to be nonspherical; in particular, interstellar object 1I/2017 U1 (‘Oumuamua) could be a very elongated (see for example Knight et al. 2017; Meech et al. 2017; Bolin et al. 2018; Drahus et al. 2018; Fraser et al. 2018) or flattened (Mashchenko 2019) body.

Sekanina (1984) developed the forced precession model of a nonspherical cometary nucleus to reproduce the effects of the torque that anisotropic outgassing exerts. This model includes additional parameters, and the direction of the spin pole (ɛ_c, ϕ_c) now depends on the oblateness of the nucleus: p = 1 − R_b∕R_a, where R_a is the equatorial radius of the object and R_b is its polar radius. When the nonspherical body rotates on its shorter axis, p > 0, we speak of an oblate nucleus; if the rotation of the nucleus is along its longer axis, p < 0, we have a prolate one.

In this work, we use the formulae of the forced precession model of a nonspherical nucleus described by Krolikowska et al. (1998). This model assumes that the nongravitational effects due to outgassing come from one dominant active area located on the surface of an ellipsoidal object of oblateness p that is spinning about a precessing axis at a given epoch with direction defined by (ɛ_c, ϕ_c); the outgassing area reaches its maximum activity at an angle η_c from the subsolar meridian as seen from the surface of the object.

The relevant part of the set of Eqs. (3) in Krolikowska et al. (1998) is: $$\begin{array}{lll}{C}_1\hfill & =\hfill & \frac{\cos {\eta }_{\text{c}}+(1-S-\cos {\eta }_{\text{c}}){\sin }^2{\varepsilon }_{\text{c}}{\sin }^2\Lambda }{\sqrt{1-S_1{\sin }^2\Psi }},\hfill \\ C_2\hfill & =\hfill & \frac{\sin {\eta }_{\text{c}}\cos {\varepsilon }_{\text{c}}+(1-S-\cos {\eta }_{\text{c}}){\sin }^2{\varepsilon }_{\text{c}}\sin \Lambda \cos \Lambda }{\sqrt{1-S_1{\sin }^2\Psi }},\hfill \\ C_3\hfill & =\hfill & \frac{(\sin {\eta }_{\text{c}}\cos \Lambda +(1-S-\cos {\eta }_{\text{c}})\cos {\varepsilon }_{\text{c}}\sin \Lambda )\sin {\varepsilon }_{\text{c}}}{\sqrt{1-S_1{\sin }^2\Psi }},\hfill \end{array}$$(1)

where S = p (2 − p), S₁ =S (2 − S), Λ = ϕ_c + f, and sin Ψ = sinɛ_c  sinΛ. For a spherical nucleus, p = 0, the set of Eqs. (1) collapses into the one shown on page 140 of Yeomans et al. (2004) and originally derived by Sekanina (1981). In both cases, the direction cosines are also given by the expressions C₁ = A₁∕A, C₂ = A₂∕A, and C₃ = A₃∕A, where A = $\sqrt{A_1^2+A_2^2+A_3^2}$. The nongravitational parameters (A₁, A₂, and A₃) can be computed from astrometric observations of the comet, and ɛ_c, ϕ_c, η_c, and p can be obtained from the components of the nongravitational force by using an iterative least squares algorithm, as shown in Krolikowska et al. (1998).

The orbit determinations of extrasolar minor bodies ‘Oumuamua and 2I/Borisov (see Tables 1–3, and 5) include the usual orbital elements – eccentricity, e, perihelion distance, q, inclination, i, longitude of the ascending node, Ω, argument of perihelion, ω, and time of perihelionpassage, τ – and the nongravitational parameters (see above). If an orbit determination is meant to reproduce actual observations accurately and to make reliable ephemeris predictions into the past and the future, one has to consider how the orbital parameters (nongravitational ones included) affect one another using the covariance matrix whose elements indicate the level to which two given parameters vary together. In this context, the mean values of the parameters and the covariance matrix define a hyperellipsoid in multidimensional space. Such data can be used to derive distributions of the angles (ɛ_c, ϕ_c, and η_c), and the oblateness parameter, p, compatible with the observations via Monte Carlo random search starting with input data generated as described by de la Fuente Marcos & de la Fuente Marcos (2015).

Here, we outline a procedure to derive the spin-axis orientation and thermal lag angle of a comet from its nongravitational orbit determination that includes the mutual uncertainties between the various parameters. This approach is applied to study the orientation of the spin axes of extrasolar minor bodies ‘Oumuamua and 2I/Borisov. This paper is organized as follows. In Sect. 2, we present the data used in our analyses, discussing and validating our techniques. In Sect. 3, we applyour methodology to ‘Oumuamua’s data. The spin-axis orientation of 2I/Borisov is studied in Sect. 4. Our results are placed within the context of previous work in Sect. 5 and our conclusions are summarized in Sect. 6.

2 Data description and methods

In this work, we find two types of assumptions that affect the results obtained. On the one hand, we have those linked to the model used to constrain the three angles (ɛ_c, ϕ_c, and η_c), and the value of the oblateness; on the other, we find those associated with the public orbit determination, that has not been computed by us, such as the one species dominating the outgassing.

The assumptions behind the forced precession model of a nonspherical cometary nucleus are described in detail in Sekanina (1984) and Krolikowska et al. (1998). The ones affecting the orbit determination are validated numerically by their ability to reproduce the available observations and predict the future positions of the object. If the orbit determination is capable of providing the actual position of the object within a reasonable level of uncertainty (often less than a few arcseconds), then one has to assume that the starting hypotheses (used to compute the orbit determination) are sufficiently justified. However, it is certainly possible that some starting hypotheses could be formally inconsistent with evidence from other sources and still produce numerically correct results, in the sense of being consistent with the observations.

2.1 Data sources

In this paper, we use publicly available data from Jet Propulsion Laboratory’s (JPL) Small-Body Database (SBDB, Giorgini 2015)¹ to investigate the spin-axis orientations of 1I/2017 U1 (‘Oumuamua) and 2I/Borisov. The data include the orbital elements and nongravitational parameters shown in Tables 1–3, and their associated covariance matrices. The data were used to construct the distributions of the angles ɛ_c, ϕ_c, and η_c, and in some cases that of the oblateness. Additional data have been retrieved from JPL’s SBDB or the Minor Planet Center (MPC, Rudenko 2016) using the tools provided by the Python package Astroquery (Ginsburg et al. 2019) and the Python module sbpy (Mommert et al. 2019).

JPL’s SBDB explicitly states that the orbit determination of ‘Oumuamua in Table 1 is identical to solution 7c in Micheli et al. (2018) that assumes a nongravitational acceleration that uses normalization factors appropriate for CO driven sublimation (see above). However, we have to acknowledge that outgassing of carbon based molecules, specifically CO and CO₂ has been definitively ruled out as accelerants by NASA’s Spitzer Space Telescope non-detection of ‘Oumuamua (Trilling et al. 2018). The authors of this analysis suggest that outgassing of another gas species, perhaps H₂O, may be responsible for the observed nongravitational acceleration; however, the producer of the orbit determination in Table 1, D. Farnocchia, is listed as a coauthor of both Micheli et al. (2018) and Trilling et al. (2018). In their Table 1, ‘Oumuamua ISSI Team et al. (2019) show that the upper limits for CO from different sources are clearly orders of magnitude lower than those of other species. In any case, this inconsistency may not have any major impact on our results.

The primary objective of this work is not in contesting the orbit determinations computed by others, but in exploring their associated effects on the possible rotational properties of the objects. As we use the numerical mean values and their associated covariance matrices, and they are able to reproduce observational results (the astrometry) within reasonable accuracy levels, we consider their numerical values as essentially valid even if the original hypotheses used to compute them could be invalid. The subject of the actual driver of the detected outgassing remains controversial in the case of ‘Oumuamua(see for example Flekkøy et al. 2019; Hui & Knight 2019; Seligman & Laughlin 2020; Hoang & Loeb 2020).

As for the orbit determinations of 2I/Borisov in Tables 2 and 3, JPL’s SBDB states that the standard model – in other words, assuming a nongravitational acceleration that uses normalization factors appropriate for H₂O driven sublimation (see above) – has been applied to compute the nongravitational parameters. A water ice driver is supported by Sekanina (2019). However, it has already been confirmed that this comet exhibits an unusually high CO abundance (Cordiner et al. 2020). Although the orbital elements in Tables 2 and 3 are, in statistical terms, mutually consistent, we cannot say the same about the nongravitational parameters. The orbit determination in Table 2 – for which the last observation used was acquired on January 6, 2020 – shows a statistically significant (above the 6.8-σ level) value of the nongravitational normal acceleration parameter (but the others are statistically compatible with zero). This is often the case when cometary outgassing comes from a single active region located at arotation pole and when the obliquity of the orbit plane to the equatorial plane is near 50° or 130° (see Sect. 3.1 in Yeomans et al. 2004).

In sharp contrast, the orbit determination in Table 3 has no value of the normal nongravitational parameter; therefore, we have to assume that it is compatible with zero. In addition, the value of A₂ is consistent with that in Table 2 but now the dominant nongravitational acceleration is radial (statistically significant above the 17-σ level), as it is in the case of ‘Oumuamua (see Table 1). The orbit determination in Table 3 includes observations acquired in January and February. Assuming that both orbit determinations are correct, they may represent two snapshots in the dynamical life of 2I/Borisov that may reflect the rapidly evolving nature of the outgassing processes taken place at the surface of the comet. In fact, there is robust observational evidence to support this interpretation: a sequence of outbursts was observed by multiple observers in February-March (Bolin et al. 2020a; Drahus et al. 2020; Jewitt et al. 2020b,c; Zhang et al. 2020).

2.2 Methodology

The first step of the approach used in this paper is an implementation of the Monte Carlo using the Covariance Matrix (MCCM, Bordovitsyna et al. 2001; Avdyushev & Banshchikova 2007) methodology in which a Monte Carlo process generates control or clone orbits based on the nominal orbit, but adding random noise on each orbital parameter by making use of the covariance matrix. Considering a covariance matrix as computed by JPL’s Solar System Dynamics Group, Horizons On-Line Ephemeris System, the vector including the mean values of the orbital parameters at a given epoch t₀ is of the form v = (e, q, τ, Ω, ω, i, A₁, A₂, A₃). If C is the covariance matrix at the same epoch associated with the nominal orbital solution that is symmetric and positive-semidefinite, then C = A A^T where A is a lower triangular matrix with real and positive diagonal elements, A^T is the transpose of A. In the case studied here, these matrices were 9 × 9 and the expressions of the individual elements a_ij of matrix A are shown in Appendix A.

If r is a vector made of univariate Gaussian random numbers – components r_i with i = 1, 9 generated using the Box-Muller method (Box & Muller 1958) – the sets of initial orbital elements of the control orbits can be computed using the expressions (assuming the structure provided by JPL’s Horizons On-Line Ephemeris System pointed out above), v_c =v + A r: $$\begin{array}{ll}{e}_{\text{c}}\hfill & =e+a_{11}{r}_1\hfill \\ q_{\text{c}}\hfill & =q+a_{22}{r}_2+a_{21}{r}_1\hfill \\ {\tau }_{\text{c}}\hfill & =\tau +a_{33}{r}_3+a_{32}{r}_2+a_{31}{r}_1\hfill \\ {\Omega }_{\text{c}}\hfill & =\Omega +a_{44}{r}_4+a_{43}{r}_3+a_{42}{r}_2+a_{41}{r}_1\hfill \\ {\omega }_{\text{c}}\hfill & =\omega +a_{55}{r}_5+a_{54}{r}_4+a_{53}{r}_3+a_{52}{r}_2+a_{51}{r}_1\hfill \\ i_{\text{c}}\hfill & =i+a_{66}{r}_6+a_{65}{r}_5+a_{64}{r}_4+a_{63}{r}_3+a_{62}{r}_2+a_{61}{r}_1\hfill \\ A_{1\text{c}}\hfill & =A_1+a_{77}{r}_7+a_{76}{r}_6+a_{75}{r}_5+a_{74}{r}_4+a_{73}{r}_3+a_{72}{r}_2\hfill \\ \hfill & +a_{71}{r}_1\hfill \\ A_{2\text{c}}\hfill & =A_2+a_{88}{r}_8+a_{87}{r}_7+a_{86}{r}_6+a_{85}{r}_5+a_{84}{r}_4+a_{83}{r}_3\hfill \\ \hfill & +a_{82}{r}_2+a_{81}{r}_1\hfill \\ A_{3\text{c}}\hfill & =A_3+a_{99}{r}_9+a_{98}{r}_8+a_{97}{r}_7+a_{96}{r}_6+a_{95}{r}_5+a_{94}{r}_4\hfill \\ \hfill & +a_{93}{r}_3+a_{92}{r}_2+a_{91}{r}_1.\hfill \end{array}$$(2)

This approach and set of equations have already been used to study the dynamics of 2I/Borisov as presented in de León et al. (2020).

In order to reconstruct the distributions of ɛ_c, ϕ_c, and η_c, and the oblateness, we generated 10⁵ instances of a given orbit (unless stated otherwise) using the set of Eqs. (2). For each one of them, we estimated the associated values of the angles (and p, if not imposed by additional evidence) by applying a Monte Carlo random search algorithm. The random search uses sets of values of the three angles (and p in relevant cases) that are generated randomly (uniformly) within the relevant intervals (see Sect. 1) until a stopping criterion is met. For example, if we consider the direction cosine for the radial component of the force (see equations in Sect. 1): $$C_1=\frac{\cos {\eta }_{\text{c}}+(1-S-\cos {\eta }_{\text{c}}){\sin }^2{\varepsilon }_{\text{c}}{\sin }^2\Lambda }{\sqrt{1-S_1{\sin }^2\Psi }},$$(3)

that is also given by the expression C₁ = A₁∕A, where A = $\sqrt{A_1^2+A_2^2+A_3^2}$. We can construct C_1c using data from MCCM and then try values of ɛ_c, ϕ_c, η_c, and p (to compute S and S₁) on Eq. (3) until |C_1c − C₁| < Δ, where Δ is a certain tolerance value that in our calculations and for practical reasons was set to 0.0001 (our results do not depend significantly on the tolerance value). A set of ɛ_c, ϕ_c, η_c, and p is regarded as a valid solution if |C_ic − C_i| < Δ, with i = 1, 3 (the radial, transverse, and normal direction cosines). We applied this technique and stopping criterion to obtain 10⁵ sets of angles and p (one per orbital instance) that we used to generate the distributions of angles and oblateness. This approach is robust when working with multimodal distributions that may signal the existence of multiple solutions. An iterative procedure like the one discussed by Sitarski (1990) may not be able to obtain all the statistically viable solutions.

Our technique delivers the distribution of spin-axis orientations defined by (ɛ_c, ϕ_c), but other alternative approaches discussed in the literature provide the ecliptic coordinates, (λ, β), of the spin-axis orientation. We can explore the associated distribution of pole directions by taking into account the dot product of the unit vector in the direction of the rotation pole – u_r = (cosβ_p cosλ_p, cosβ_p sinλ_p, sinβ_p) – with that of the orbital pole – u_o = (sini_c sinΩ_c, −sini_c cosΩ_c, cosi_c) – that gives u_r ⋅u_o = cosɛ_c. Therefore, the expression that can be used to apply a random search technique with a stopping criterion to find (λ_p,  β_p) is: cos β_p cosλ_p sini_c sinΩ_c − cosβ_p sinλ_p sini_c cosΩ_c + sinβ_p cosi_c = cosɛ_c. From here, we can compute the equivalent distribution in equatorial coordinates, (α, δ), using the conversion tools provided by Astropy (Astropy Collaboration 2013, 2018).

In order to analyze the results, we produced histograms using the Matplotlib library (Hunter 2007) with two different sets of bins computed using Astropy (Astropy Collaboration 2013, 2018): one by applying the Freedman and Diaconis rule (Freedman & Diaconis 1981) and a second one using the Bayesian Blocks technique (Scargle et al. 2013). While the Freedman and Diaconis rule produces a constant bin size, the Bayesian Blocks technique allows varying bin widths that work well for complicated distributions like the ones studied here. We consider two different statistically motivated histograms to show that our conclusions do not depend on the choice of bin width. On the other hand and instead of using frequency-based histograms, we considered counts to form a probability density so the area under the histogram will sum to one. In addition, we construct diagrams to visualize better the statistical significance of the spin-axis orientations – as (ɛ_c, ϕ_c), (λ_p, β_p), or (α_p, δ_p) maps – representing kernel-density estimates using Gaussian kernels (see for example Scott 1992) implemented by SciPy (Virtanen et al. 2020).

The spin-axis orientations obtained using the methodology described above are referred to the same epoch as the associated orbit determination. We have not studied the time evolution of the rotational properties of minor bodies as described by, for example, Krolikowska et al. (1998). However, for an object following a hyperbolic path and in absence of outbursts, fragmentation events, or close encounters with massive planets, the overall spin-axis orientation is expected to remain relatively stable as it is often the case for each apparition (not for multiple consecutive apparitions)of well-studied comets, see Sects. 3.5–3.8 in Yeomans et al. (2004). In this context, we may interpret the results obtained for a given epoch as the most probable values of the rotational properties around that epoch.

Fig. 1 Distributions of equatorial obliquity and cometocentric longitude of the Sun at perihelion for 51P/Harrington. Top left panel: histogram of equatorial obliquity, ɛ51P/Harrington, with bins computed using the Freedman and Diaconis rule, while top right panel: uses the Bayesian Blocks technique. Similarly, the cometocentric longitude, ϕ51P/Harrington, histograms displayed in the lower panels also use the Freedman and Diaconis rule (left) and Bayesian Blocks (right). Histograms are based on data from Table 4 (see the text for details).

2.3 Validation: comet 51P/Harrington

The nongravitational motion of comet 51P/Harrington has been studied by Sitarski (1996). This comet has experienced multiple fragmentation events since its discovery on August 14, 1953 (Sekanina 1997; Kidger & Manteca 2002). The work by Sitarski (1996) is referred to the epoch February 12, 1995 and the data-arc spans from discovery time to January 28, 1995. The values of the relevant properties are ɛ_{51P/Harrington} = 77. °6 ± 1. °1, ϕ_{51P/Harrington} = 92° ± 4°, η_{51P/Harrington} = 18. °4 ± 1. °4, and p_{51P/Harrington} = −0.32 ± 0.08 (Sitarski 1996), and they have been computed using Sekanina’s forced precession model of a nonspherical cometary nucleus. Sitarski (1996) considered asymmetric cometary activity with respect to perihelion in his calculations and A = 7.9192 × 10⁻⁹ ± 1.987 × 10⁻¹⁰ AU d⁻², but the values of the individual nongravitational acceleration parameters were not provided.

Table 4 shows the latest orbit determination available for comet 51P/Harrington. JPL’s SBDB states that the orbit determination has been computed using the standard model (see above). When comparing the data in Table 4 with those in Tables 1–3, we realize that the statistical relevance of the nongravitational parameters is substantially higher in the case of comet 51P/Harrington. It is also important to emphasize that 51P/Harrington is a periodic comet that has been studied for multiple decades (its orbital period is 7.16 yr); interstellar objects can only be studied once and they have comparatively short visibility windows.

Applying the procedure outlined in Sect. 2 using as input data those in Table 4 and the associated covariance matrix, we generated 10⁴ instances of a given orbit using the set of Eqs. (2). The results for the spin-axis orientation (ɛ_{51P/Harrington}, ϕ_{51P/Harrington}) are shown in Fig. 1. The top panels display the distribution of the equatorial obliquity that is unimodal with median and 16th and 84th percentiles of ${82}^{°}{}_{-2^{°}}^{+6^{°}}$; the bottom panels show a bimodal distribution for the cometocentric longitude with average maxima at about 55° and 235° (therefore, separated by 180°).

These results indicate that the effective equatorial plane of 51P/Harrington during its most recent perihelion was nearly perpendicular to its orbital plane; on the other hand, if ϕ_{51P/Harrington} ~ 55° the northern pole of the comet was facing the Sun at perihelion, or the southern one if ϕ_{51P/Harrington} ~ 235°. Our result for the value of the equatorial obliquity is fully consistent with that in Sitarski (1996), in particular with the trend shown in his Fig. 1. The value ofϕ_{51P/Harrington} in Fig. 1, bottom panels, is however inconsistent with the one in Sitarski (1996) and the trend in his Fig. 2; the median and 16th and 84th percentiles of ϕ_{51P/Harrington} are ${93}^{°}{}_{-{38}^{°}}^{+{145}^{°}}$. It can be argued that the differences in ϕ_{51P/Harrington} could be theresult of the fragmentation events pointed out above.

Figure 2 shows an asymmetric distribution for η_{51P/Harrington}, with no clear maximum. The median and 16th and 84th percentiles of the distribution are ${-38}^{°}{}_{-{89}^{°}}^{+{135}^{°}}$. This impliesthat the maximum outgassing was taking place when the Sun was rising as seen from the surface of the comet. Again, this result is inconsistent with the one in Sitarski (1996), but it might be a by-product of the fragmentation events pointed out above if the fresh surface has a different value of the thermal inertia.

Figure 3 shows a bimodal distribution for the oblateness of the nucleus of 51P/Harrington. The median and 16th and 84th percentiles of the distribution are $-{0.50}_{-1.5}^{+0.7}$, with maxima at about 0.06 and −0.6. Our results favor an oblate nucleus although a prolate one cannot be discarded. Considering the large dispersion, our median value is statistically consistent with the results in Sitarski (1996) and our secondary maximum is marginally consistent with his result. Again, the fact that the surface of the comet has experienced multiple alterations since the data used by Sitarski (1996) were acquired and that some of these episodes may have even altered the shape of 51P/Harrington make it difficult to validate our results against those in Sitarski (1996). This is particularly obvious when we consider that the value of A in Table 4 is nearly five times larger than the one computed by Sitarski (1996). In addition, the model used to derive the values of the parameters in Table 4 considered symmetric – not asymmetric as in Sitarski (1996) – cometary activity with respect to perihelion.

Taking into account the dot product of the unit vector in the direction of the rotation pole with that of the orbital pole and its relationship with ɛ_{51P/Harrington} as pointed out above, we have computed the distributions of (λ_p, β_p) and (α_p, δ_p). Figure 4 shows the resulting distributions in ecliptic coordinates and Fig. 5 shows the maps of statistical significance computed as described in the second to last paragraph of Sect. 2.2. The most probable spin-axis direction of 51P/Harrington in equatorial coordinates could be (260°, −13°) or (120°, +20°). On the other hand, the equatorial coordinates of the Sun when the comet reached perihelion on August 12, 2015 were approximately (141. °8, +15. °0), somewhat consistent with the orientation of one of the poles.

Fig. 2 Distribution of thermal lag angle for 51P/Harrington. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 4.

Fig. 3 Distribution of values of the oblateness of the nucleus for 51P/Harrington. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 4.

Fig. 4 Distributions of spin-axis orientations, (λp, βp), for 51P/Harrington. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 4.

3 Interstellar object 1I/2017 U1 (‘Oumuamua)

We applied the procedure outlined in Sect. 2 to produce the distributions of ɛ_c, ϕ_c, η_c, and p_c for interstellar object 1I/2017 U1 (‘Oumuamua) considering the orbit determination in Table 1 that is discussed by Micheli et al. (2018). The model favored by these authors assumes that ‘Oumuamua’s comet-like outgassing is driven by CO, not H₂O (see Sect. 1 and the discussion in Sect. 2). The results for the spin-axis orientation (ɛ_‘Oumuamua, ϕ_‘Oumuamua) are displayed in Fig. 6, which summarizes the results of 10⁵ Monte Carlo experiments.

Figure 6, top panels, shows the distribution of the equatorial obliquity. The median and 16th and 84th percentiles of the distribution are ${93}^{°}{}_{-{13}^{°}}^{+8^{°}}$. This result indicates that the effective equatorial plane of ‘Oumuamua was nearly perpendicular to its orbital plane. Figure 6, bottom panels, shows a bimodal distribution with average maxima at about 140° and 320° (therefore, separated by 180°). If ϕ_‘Oumuamua ~ 140° the northern pole of ‘Oumuamua was facing the Sun at perihelion, if ϕ_‘Oumuamua ~ 320° was the southern one. In general, our methodology cannot always distinguish which pole was facing the Sun when ‘Oumuamua reached perihelion.

Figure 7 also shows a unimodal distribution for η_‘Oumuamua. The median and 16th and 84th percentiles of the distribution are ${0.08}^{°}{}_{-{55}^{°}}^{+{52}^{°}}$. This impliesthat the maximum outgassing was taking place mainly at ‘Oumuamua’s noon, when the Sun was at its highest as seen from the surface of the interstellar visitor.

Figure 8 shows a unimodal distribution for the oblateness of ‘Oumuamua with median and 16th and 84th percentiles of $-{0.13}_{-1.56}^{+0.22}$. This result slightly favors a fusiform shape for ‘Oumuamua over a disk-like one. The actual shape of ‘Oumuamua remains a controversial subject (see for example the discussion in Vazan & Sari 2020). Although multiple studies have pointed out that the available observational data strongly favor that ‘Oumuamua has a very elongated shape (see for example Knight et al. 2017; Bolin et al. 2018; Drahus et al. 2018; Fraser et al. 2018), Belton et al. (2018) stated that its shape might be fusiform, but a disk-like appearance could not be ruled out.

Mashchenko (2019) presented a very detailed modeling of the light curve of ‘Oumuamua to show that his best-fitting model, with a probability of 91%, was a thin disk with p = 0.836; his second choice, with a probability of 16%, was a thin spindle with p = −6.7. The thin-disk model was significantly more successful at reproducing the observed variations in apparent magnitude over time. However, the results presented by Mashchenko (2019) have been contested by Zhang & Lin (2020) by arguing that multiple lines of evidence favor the very prolate or cigar-like shape against the very oblate or disk-like shape. In particular, they pointed out that the fusiform case is energetically more stable. In the following, we test the impact of both models on the values of the rotational parameters of ‘Oumuamua. We first assume that p = −6.7 together with the data in Table 1, then we consider p = 0.836.

Figures 9 and 10, show the distributions of ɛ_‘Oumuamua, ϕ_‘Oumuamua and η_‘Oumuamua under the assumption of p = −6.7; we generated 10⁴ instances of a given orbit. The results for the orientation of the pole are similar to those in Fig. 6: the equatorial obliquity becomes 93° ± 3° and ϕ_‘Oumuamua ~ 140° or ~ 320°. However, Fig. 10 shows a relatively flat distribution for η_‘Oumuamua; in other words, outgassing takes place evenly throughout a rotation period. The actual location of the pole is studied in Figs. 11 and 12: under the thin-spindle scenario, the most probable spin-axis direction of ‘Oumuamua in equatorial coordinates could be (280°, +46°), with approximate Galacticcoordinates l = 75°, b = +21°, that point slightly above the Galactic disk, toward the constellation of Lyra.

Figures 13 and 14, show the distributions of ɛ_‘Oumuamua, ϕ_‘Oumuamua, and η_‘Oumuamua for 10⁴ instances of the orbit, under the assumption of p = 0.836. The results for the orientation of the pole are very different from those in Fig. 6: the equatorial obliquity is now ${16}^{°}{}_{-7^{°}}^{+{80}^{°}}$ and ϕ_‘Oumuamua ~ 67° or ~247°. However, the distribution of η_‘Oumuamua in Fig. 14 resembles that in Fig. 7; maximum outgassing would take place at ‘Oumuamua’s noon with η_‘Oumuamua = 3° ± 2°. The location of the pole under the thin-disk scenario is shown in Figs. 15 and 16; the most probable spin-axis direction of ‘Oumuamua in equatorial coordinates could be (312°, −50°). This is consistent with the estimate of the orbital pole of ‘Oumuamua in de la Fuente Marcos & de la Fuente Marcos (2017a,b): compare Fig. 16, middle panel, with the middle panel in Fig. 1 of de la Fuente Marcos & de la Fuente Marcos (2017b), both in ecliptic coordinates. The approximate Galactic coordinates of the rotational pole are l = 349°, b = −39°, that point slightly below the Galactic bulge, toward the constellation Indus.

The rotational properties provided by our modeling are very different depending on how extreme the shape of ‘Oumuamua is; this suggests that the spin-axis orientation of ‘Oumuamua is tightly controlled by its shape. In the Solar System, most cometary nuclei are closer to a prolate configuration, not oblate (Samarasinha et al. 2004). Although the work by Mashchenko (2019) is widely regarded as the best available at the moment, the arguments put forward by Zhang & Lin (2020) may indicate that a final answer on ‘Oumuamua’s shape still remains elusive.

Fig. 5 Gaussian kernel density estimation of the spin-axis orientations computed using the SciPy library (Virtanen et al. 2020). Top panel: (ɛc, ϕc) map. The results in ecliptic coordinates, (λp, βp), are shown in the middle panel and those in equatorial coordinates, (αp, δp), are displayed in the bottom panel that includes the actual observations (in red) of comet 51P/Harrington available from the MPC.

Fig. 6 Distributions of equatorial obliquity and cometocentric longitude of the Sun at perihelion for 1I/2017 U1 (‘Oumuamua). Top left panel: histogram of equatorial obliquity, ɛ‘Oumuamua, with bins computed using the Freedman and Diaconis rule, while top right panel: uses the Bayesian Blocks technique. Similarly, the cometocentric longitude, ϕ‘Oumuamua, histograms displayed in the lower panels also use the Freedman and Diaconis rule (left) and Bayesian Blocks (right). Histograms are based on data from Table 1 (see the text for details).

Fig. 7 Distribution of thermal lag angle for 1I/2017 U1 (‘Oumuamua). Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 1.

Fig. 8 Distribution of values of the oblateness of interstellar object 1I/2017 U1 (‘Oumuamua). Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel uses the Bayesian Blocks technique. Histograms are based on data from Table 1.

Fig. 9 Same as Fig. 6 but assuming a value for the oblateness of interstellar object 1I/2017 U1 (‘Oumuamua) equal to −6.7 as suggestedby Mashchenko (2019).

Fig. 10 Same as Fig. 7 but assuming a value for the oblateness of interstellar object 1I/2017 U1 (‘Oumuamua) equal to −6.7 as suggestedby Mashchenko (2019).

Fig. 11 Distributions of spin-axis orientations, (λp, βp), for 1I/2017 U1 (‘Oumuamua) assuming a very prolate shape, with p = −6.7. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 1.

Fig. 12 Gaussian kernel density estimation of the spin-axis orientations computed using the SciPy library (Virtanen et al. 2020) for a very prolate case (p = −6.7). Top panel: (ɛc, ϕc) map. The results in ecliptic coordinates, (λp, βp), are shown in the middle panel and those in equatorial coordinates, (αp, δp), are displayed in the bottom panel that includes the actual observations (in red) of interstellar object 1I/2017 U1 (‘Oumuamua) available from the MPC.

4 Interstellar comet 2I/Borisov

It is unlikelythat additional data on 1I/2017 U1 (‘Oumuamua) may eventually emerge; in sharp contrast, interstellar comet 2I/Borisov is still being actively observed. The two most recent orbit determinations of 2I/Borisov are shown in Tables 2 and 3 and both consider nongravitational effects driven by sublimation of water ice (see Sect. 1). However, Ye et al. (2020) argued that a model based on sublimation of CO could be more consistent with the observations compared to a H₂O model. In addition, Bodewits et al. (2020) have pointed out that the coma of 2I/Borisov contains significantly more CO than H₂O gas, with abundances of at least 173%, more than three times higher than previously measured for any comet in the inner (<2.5 AU) Solar System. However, the water production rate prior to perihelion was found to increase faster than for most known dynamicallynew comets (Xing et al. 2020); a hyperactive nucleus was also favored by McKay et al. (2020). Again and as considered above, an orbit determination may produce numerically correct results (in the sense of correct ephemeris predictions) even if the hypotheses used to carry out the calculations are physically questionable or even invalid.

The orbit determination in Table 2 corresponds to an epoch nearly a month past perihelion and shows that only the normal acceleration may be playing a significant role on the dynamical evolution of 2I/Borisov as the other components have values statistically compatible with zero. Using data from Table 2 as input, weobtained the distributions in Figs. 17 and 18. The distribution of ɛ_2I/Borisov is unimodal but wide (see Fig. 17, top panels); there is one statistically significant maximum and the median and 16th and 84th percentiles are ${59}^{°}{}_{-{18}^{°}}^{+{15}^{°}}$. Figure 17, bottom panels, shows a multimodal distribution for ϕ_2I/Borisov; the absolute maxima are found at about 90° and 270°, with secondary maxima at 40° and 220° (therefore, each pair separated by 180°), and the associated dispersions are relatively small. Figure 18 displays a multimodal distribution for η_2I/Borisov, with a dominant maximum at about −130°, which is on the nightside of the comet.

The most probable value of the equatorial obliquity of 2I/Borisov is sometimes found for periodic comets in the Solar System, when outgassing comes from a single active region located at a rotation pole (see Sect. 3.1 in Yeomans et al. 2004). Surprisingly, outgassing appears to be restricted to the nightside of the comet. Most comets tend to exhibit afternoon activity (Sekanina 1981). The lack of dayside activity seems problematic when considering that ‘Oumuamua’s activity was restricted to the dayside for the very oblate case; however, outbursts from the nightside of comets 9P/Tempel 1 (Farnham et al. 2007, 2013), 103P/Hartley 2 (A’Hearn et al. 2011; Bruck Syal et al. 2013), and 67P/Churyumov-Gerasimenko (Knollenberg et al. 2016; Rinaldi et al. 2019) have been observed.

Figure 19 shows that the most probable value for the oblateness parameter for 2I/Borisov corresponds to that of a very oblate body with median and 16th and 84th percentiles of ${0.6}_{-0.7}^{+0.3}$; the peak of the distribution is found at p = 0.95. This is even more extreme than the value of p = 0.836 favored by Mashchenko (2019) for 1I/2017 U1 (‘Oumuamua).

The location of the pole is shown in Figs. 20 and 21. The most probable spin-axis direction of 2I/Borisov in equatorial coordinates could be (275°, +65°), see Fig. 21, bottom panel, with approximate Galactic coordinates l = 95°, b = +28° that point slightly above the Galactic disk, toward the constellation Draco.

Our results are quite different when considering the most recent orbit determination of 2I/Borisov in Table 3. The distribution of ɛ_2I/Borisov in Fig. 22, top panels, is still unimodal but now the median and 16th and 84th percentiles are 90° ± 52°. The distribution of ϕ_2I/Borisov in Fig. 22, bottom panels, is multimodal as the one in Fig. 17, bottom panels, but now the dominant maxima are found at about 100° and 280°. The distribution of η_2I/Borisov in Fig. 23 is very different from that in Fig. 18 and shows peaks at −20° and 20°, with median and 16th and 84th percentiles of 0° ± 29°; the maximum of activity takes place around noon as seen from the comet.

Figure 24 shows that the distribution of the values for the oblateness parameter of 2I/Borisov is bimodal with an absolute maximum of about 0.34 and a secondary peak at −0.36. Although our analysis favors an oblate shape for 2I/Borisov, a prolate one cannot be ruled out.

The location of the pole is shown in Figs. 25 and 26. The most probable spin-axis direction of 2I/Borisov in equatorial coordinates could be (231°, +30°), see Fig. 26, bottom panel, with approximate Galactic coordinates l = 47°, b = +56° that point well above the Galactic disk, toward the constellation Corona Borealis.

Fig. 13 Same as Fig. 6 but assuming a value for the oblateness of interstellar object 1I/2017 U1 (‘Oumuamua) equal to 0.836 as suggested by Mashchenko (2019).

Fig. 14 Same as Fig. 7 but assuming a value for the oblateness of interstellar object 1I/2017 U1 (‘Oumuamua) equal to 0.836 as suggested by Mashchenko (2019).

Fig. 15 Distributions of spin-axis orientations, (λp, βp), for 1I/2017 U1 (‘Oumuamua) assuming a very oblate shape, with p = 0.836. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 1.

Fig. 16 Gaussian kernel density estimation of the spin-axis orientations computed using the SciPy library (Virtanen et al. 2020) for a very oblate case (p = 0.836). Top panel: (ɛc, ϕc) map. The results in ecliptic coordinates, (λp, βp), are shown in the middle panel and those in equatorial coordinates, (αp, δp), are displayed in the bottom panel that includes the actual observations (in red) of interstellar object 1I/2017 U1 (‘Oumuamua) available from the MPC.

Fig. 17 Distributions of equatorial obliquity and cometocentric longitude of the Sun at perihelion for 2I/Borisov. These results correspond to the orbit determination of comet 2I/Borisov in Table 2.

Fig. 18 Distribution of thermal lag angle for 2I/Borisov. These results correspond to the orbit determination of comet 2I/Borisov in Table 2.

Fig. 19 Distribution of the oblateness parameter for 2I/Borisov. These results correspond to the orbit determination of comet 2I/Borisov in Table 2.

Fig. 20 Distributions of spin-axis orientations, (λp, βp), for 2I/Borisov. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 2.

5 Discussion

The main objective of the research presented here is not to argue against or in favor of the various orbit determinations (and their underlying outgassing models) used to generate input data for our analysis, but to point out the consequences and implications that a given nongravitational solution may have on certain comet properties such as the spin-axis orientation, the thermal lag angle, or the value of the oblateness. On the other hand, both 1I/2017 U1 (‘Oumuamua) and 2I/Borisov are single-apparition minor bodies and it is unclear how many of the trends observed forwell-studied periodic comets and discussed in the literature can be expected to be present for these extrasolar visitors. Sekanina (1981) argues that most periodic comets tend to have values of the obliquity close to 90° and the spinaxis tends to point toward the Sun at perihelion, which results in the most extreme insolation possible. In addition, the forced precession model of a nonspherical cometary nucleus provides better results when applied to data from several consecutive apparitions of the body under study (Krolikowska et al. 1998), but this is clearly not possible when considering interstellar objects.

Fig. 21 Gaussian kernel density estimation of the spin-axis orientations of 2I/Borisov computed using the SciPylibrary (Virtanen et al. 2020) and based on data from Table 2. Top panel: (ɛc, ϕc) map. The results in ecliptic coordinates, (λp, βp), are shown in the middle panel and those in equatorial coordinates, (αp, δp), are displayed in the bottom panel that includes the actual observations (in red) of 2I/Borisov available from the MPC. Already published spin-axis directions are plotted as filled orange symbols with white error bars (when available): triangles (Ye et al. 2020), square (Manzini et al. 2020), star (Bolin et al. 2020b), and circle (Kim et al. 2020).

5.1 1I/2017 U1 (‘Oumuamua)

Although noestimates of the orientation of ‘Oumuamua’s rotation pole have been made public yet, Fitzsimmons et al. (2018) assumed an obliquity of 0° in their analyses and Micheli et al. (2018) considered 45°, but additional modeling by Zhou (2019) indicated that ‘Oumuamua’s effective obliquity when it entered the Solar System may have been 93. °5 (this result has been removed from the final, published version of this work, Zhou 2020). The value assumed by Fitzsimmons et al. (2018) is marginally consistent with the one computed here for a very oblate shape (see Fig. 13, top panel); the value obtained by Zhou (2019) coincides with our estimate in the case of a very prolate shape (see Fig. 9, top panel).

Seligman et al. (2019) argued that ‘Oumuamua’s outgassing is made of H₂O molecules and that the venting of this material is expected to be collimated toward the Sun. We found a small lag angle between the subsolar meridian and the direction of maximum mass ejection for the case of a very oblate shape, but the very prolate case shows no statistically significant value of the lag angle, which can be interpreted as lack of correlation between level of activity and elevation of the Sun over the local horizon. On the other hand, the equatorial coordinates of the Sun when ‘Oumuamua reached perihelion on September 9, 2017 were approximately (167. °8, +5. °2), therefore none of our predictions regarding the spin-axis direction are close to the position favored (see above) by Sekanina (1981).

Rafikov (2018) showed that under standard conditions, if ‘Oumuamua is an elongated body experiencing outgassing, it should have followed a rapid and violent rotational evolution leading to rotational fission and probably catastrophic disruption prior to its discovery. The fact is that ‘Oumuamua was discovered and Seligman et al. (2019) explained that being long-term structurally stable, outgassing, and having an extremely oblong appearance could be simultaneously possible if the active region tracks the subsolar point (maximal venting occurs when both are aligned). The modeling presented by Mashchenko (2019) is consistent with the discussion in Seligman et al. (2019).

Our analysis of the rotational properties of ‘Oumuamua in Sect. 3 shows that, for the same orbit determination in Table 1, a very prolate or cigar-like shape with p = −6.7 leads to outgassing taking place evenly throughout the rotation period (see Fig. 10); in sharp contrast, a very oblate or disk-like shape with p = 0.836 leads to a strong correlation between time of maximum activity or outgassing and alignment between active region and subsolar point (see Fig. 14) as demanded by Seligman et al. (2019). On the other hand, the Solar System is no strange to very oblate or flattened bodies, 486958 Arrokoth (2014 MU₆₉), a trans-Neptunian contact binary, has been found to have an average oblateness of p = 0.695, with 0.808 for the large lobe and 0.582 for the small lobe (Amarante & Winter 2020).

It can be argued that the values of the parameters obtained when no value of p is imposed (${\varepsilon }_{\text{`Oumuamua}}={93}^{°}{}_{-{13}^{°}}^{+8^{°}}$, ${\eta }_{\text{`Oumuamua}}={0.08}^{°}{}_{-{55}^{°}}^{+{52}^{°}}$) lead to a trivial solution to the problem under study here. The effective equatorial plane of ‘Oumuamua was nearly perpendicular to its orbital plane and the maximum outgassing was taking place at ‘Oumuamua’s noon, when the Sun was at its highest as seen from the surface of the interstellar visitor. In this case and if ‘Oumuamua is fusiform, its body would be parallel to its orbital plane and the thermal lag could be effectively zero, then the outgassing is confined only to the radial direction and the two poles are indistinct. This is caused by the fact that both the transverse and normal components of the force are statistically compatible with zero. The inclusion of information on p leads to more complex solutions.

The orbit determination in Table 1 corresponds to 74 d past perihelion and shows that only the radial acceleration may have played a significant role on the dynamical evolution of ‘Oumuamua as it visited the inner Solar System because the other components have rather uncertain values, statistically compatible with zero (Micheli et al. 2018). ‘Oumuamuais in a tumbling rotational state (Belton et al. 2018; Drahus et al. 2018; Fraser et al. 2018). Belton et al. (2018) argued for two periods – a short one of 8.7 h and a long one (slow precession) of 54.5 h that are the result of tumbling – and a long, narrow rod-like shape. Drahus et al. (2018) found a variable light curve with a periodicity of 7.56 h consistent with a very elongated shape resulting from a catastrophic collision in the distant past. Extensive reviews on ‘Oumuamua’s properties can be found in Hainaut et al. (2018) and ‘Oumuamua ISSI Team et al. (2019).

Fig. 22 Distributions of equatorial obliquity and cometocentric longitude of the Sun at perihelion for 2I/Borisov. Same as Fig. 17 but for the orbit determination of comet 2I/Borisov in Table 3.

Fig. 23 Distribution of thermal lag angle for 2I/Borisov. Same as Fig. 18 but for the orbit determination of comet 2I/Borisov in Table 3.

Fig. 24 Distribution of the oblateness parameter for 2I/Borisov. Same as Fig. 19 but for the orbit determination of comet 2I/Borisov in Table 3.

Fig. 25 Distributions of spin-axis orientations, (λp,  βp), for 2I/Borisov. Same as Fig. 20 but based on data from Table 3.

Fig. 26 Same as Fig. 21 but based on data from Table 3.

5.2 2I/Borisov

Although no estimates on the spin-axis orientation of ‘Oumuamua have been presented in the literature, there are a number of released results for 2I/Borisov although they are all mutually exclusive. Ye et al. (2020) found two equally preferred locations (see Fig. 7 in their paper) for the pole’s right ascension and declination: (340°, +30°) and (205°, −55°). Kim et al. (2020) discussed a persistent asymmetry in the dust coma that could be best explained by a thermal lag on the rotating nucleus causing peak mass loss to take place during the comet nucleus afternoon; these authors calculated a value for the obliquity of 2I/Borisov of 30° and estimated a pole direction of (α, δ) = (205°, +52°). Manzini et al. (2020) computed a spin-axis direction (α, δ) = (260° ± 15°, −35° ± 10°). Bolin & Lisse (2020) estimated a value of the oblateness ≤−0.5 and argued that a jet is close to the rotation pole that points toward (α, δ) = (322° ± 10°, +37° ± 10°), mentioning that the rotation of the nucleus of comet 2I/Borisov seems to occur on a single axis and that it appears to be not chaotic.

Together with our own estimates in the form of maps of statistical significance, the various estimates for the spin-axis orientation of 2I/Borisov are plotted in Figs. 21 and 26, bottom panels, as filled orange symbols with white error bars (when available): triangles (Ye et al. 2020), square (Manzini et al. 2020), star (Bolin & Lisse 2020), and circle (Kim et al. 2020). The published solutions for the spin-axis orientation of 2I/Borisov are quite different and they may signal actual short-term changes in the values of the rotational properties of this object or problems with the interpretation of some of the observed features. Our estimate based on Table 2 (see Fig. 21, bottom panel) is inconsistent with any of the other four, but the one based on the orbit determination in Table 3 (see Fig. 26, bottom panel) is marginally consistent with the result obtained by Kim et al. (2020) although the associated values of the obliquity are quite different, 30° versus 90°. Our value of the oblateness for the solution in Table 3 is inconsistent with the one computed by Bolin & Lisse (2020). On the other hand, the equatorial coordinates of the Sun when 2I/Borisov reached perihelion on December 8, 2019 were approximately (254. °6, -22. °7), so none of our predictions are close to the position favored by Sekanina (1981). This location is also inconsistent with the estimates provided by Ye et al. (2020), Bolin & Lisse (2020), and Kim et al. (2020); it is, however, somewhat consistent with the spin-axis direction computed by Manzini et al. (2020).

Most of the observed properties of interstellar comet 2I/Borisov have been found to be remarkably similar to those of known Solar System comets (see for example de León et al. 2019, 2020; Fitzsimmons et al. 2019; Opitom et al. 2019; Cremonese et al. 2020; Guzik et al. 2020; Jewitt et al. 2020a; Kareta et al. 2020; Yang et al. 2020). Recent estimates, placed the rotation period of 2I/Borisov at about 5.3 h or perhaps 10.6 h (Bolin & Lisse 2020). Gladman et al. (2019) had found a very slow rotation rate with a period of about 13 d. If confirmed, this would be the longest rotation period ever observed in a comet. Slowly rotating Solar System comets, such as 1P/Halley and 109P/Swift-Tuttle, have rotation periods of about 2.8 d (see for example Table 1 in Samarasinha et al. 2004); comet 29P/Schwassmann-Wachmann 1 – well known for its enormous and random outbursts – was reported to have a rotation period of about 5 d (see for example Sekanina 1981), but recent determinations place its spin rate close to 1.4 d (see for example Table 1 in Samarasinha et al. 2004). JPL’s SBDB indicates that, among known comets, P/2006 HR30 (Siding Spring) has the longest rotation period with 2.9 d.

The thermal lag angle is relatively small for our results based on Table 3 and this suggests that outgassing for 2I/Borisov appears to have been active only when the Sun was rather high over the local horizon. As for the large lag angle obtained in the case of the results based on Table 2, it suggests that an insulating crust of nonvolatile material might have been present over most of the surface of this comet prior to the sequence of outbursts pointed out above; no significant change in the slightly reddish color of the comet was observed from late September 2019 through late January 2020 (Hui et al. 2020). Alternatively, this could be the result of its putative long period as the effects of the outgassing from the comet’s dayside may cancel out and only discrete nightside outbursts may result in a net contribution to the value of nongravitational normal acceleration parameter. Yet another scenario may involve an active vent close to the south pole of the comet that, due to the viewing geometry, remained on the nightside of the comet during the time interval covered by the observations. Both Manzini et al. (2020) and Bolin & Lisse (2020) argued that such a scenario may have been possible.

An improved orbit determination of 2I/Borisov was published on August 21, 2020 (see Table 5), when this work was under review, that also considers nongravitational effects driven by sublimation of water ice (see Sect. 1). This solution is close to that in Table 3 although it assummes asymmetric nongravitational forces including an additional parameter, dτ, to account for an asymmetry of the comet’s outgassing (see for example Yeomans et al. 2004). We have repeated the calculations discussed in Sect. 4 using the data in Table 5 as input and found results that are similar (see Fig. 27) to those obtained for the orbit determination in Table 3. The inclusion of the nongravitational perihelion offset has no effects on the mathematical expressions used to perform our calculations (see Sect. 1).

Fig. 27 Same as Fig. 26, bottom panel, but based on data from Table 5.

6 Conclusions

In this paper, we have shown that the most probable values of the equatorial obliquity and the cometocentric longitude of the Sun at perihelion as well as the value of the thermal lag angle and the oblateness of a comet can be computed from an orbit determination that includes nongravitational parameters using a Monte Carlo random search approach within the framework of the forced precession model of a nonspherical cometary nucleus. The algorithm receives input from a MCCM process thatincludes the uncertainties of the orbit determination in a self-consistent way. The resulting combined algorithm is applied to investigate the spin-axis orientations of extrasolar minor bodies 1I/2017 U1 (‘Oumuamua) and 2I/Borisov.

For ‘Oumuamua and using data from Table 1, we found that its equatorial obliquity could be about 93°, if it has a very prolate (fusiform) shape, or close to 16°, if it is very oblate (disk-like). Outgassing for ‘Oumuamua appears to have been active only when the Sun was highest over the local horizon, as the thermal lag angle is small, if the body has a thin-disk appearance; for a cigar-like shape, outgassing may have been evenly distributed in time. Our calculations suggest that the most probable spin-axis direction of ‘Oumuamua in equatorial coordinates could be (280°, +46°) if very prolate or (312°, −50°) if very oblate. Our analysis favors a prolate shape. No previous estimates of the direction of the polar axis of ‘Oumuamua are available.

For 2I/Borisov, three different orbit determinations gave values of 59° (using data from Table 2) and 90° (using data from Tables 3 and 5) for its obliquity, with most probable poles pointing near (275°, +65°) and (231°, +30°), respectively. Although our analysis favors an oblate shape (p = 0.34) for 2I/Borisov, a prolate one cannot be ruled out. Our estimates for the spin-axis direction of 2I/Borisov are inconsistent with published determinations (except perhaps with the one in Kim et al. 2020), which already were mutually exclusive.

No further data may be collected from ‘Oumuamua, but new data on 2I/Borisov will be made public and they may lead to better, morereliable estimates of its rotational properties. The new data may improve what is already known and/or perhaps provide new snapshots of the evolution of its rotation.

Acknowledgements

We thank the referee for her/his constructive, detailed, and insightful report that included very helpful suggestions regarding the presentation of this paper and the interpretation of our results, J. de León, J. Licandro, and M. Serra-Ricart for comments on 2I/Borisov’s ongoing cometary activity, and A. I. Gómez de Castro for providing access to computing facilities. Part of the calculations and the data analysis were completed on the Brigit HPC server of the “Universidad Complutense de Madrid” (UCM), and we thank S. Cano Alsúa for his help during this stage. This research was partially supported by the Spanish “Ministerio de Economía y Competitividad” (MINECO) under grant ESP2017-87813-R. In preparation of this paper, we made use of the NASA Astrophysics Data System, the ASTRO-PH e-print server, and the MPC data server. This research made use of Astropy², a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

Appendix A Elements of the matrix A

If the elements of the covariance matrix C are written as c_ij and those of A as a_ij (so C = A A^T), where those are the entries in the ith row and jth column, they are related by the following expressions: $$a_{11}=\sqrt{c_{11}}{a}_{21}=c_{12}/a_{11}{a}_{31}=c_{13}/a_{11}{a}_{41}=c_{14}/a_{11}{a}_{51}=c_{15}/a_{11}{a}_{61}=c_{16}/a_{11}{a}_{71}=c_{17}/a_{11}{a}_{81}=c_{18}/a_{11}{a}_{91}=c_{19}/a_{11}$$ $$a_{22}=\sqrt{c_{22}-a_{21}^2}{a}_{32}=(c_{23}-a_{21}{a}_{31})/a_{22}{a}_{42}=(c_{24}-a_{21}{a}_{41})/a_{22}{a}_{52}=(c_{25}-a_{21}{a}_{51})/a_{22}{a}_{62}=(c_{26}-a_{21}{a}_{61})/a_{22}{a}_{72}=(c_{27}-a_{21}{a}_{71})/a_{22}{a}_{82}=(c_{28}-a_{21}{a}_{81})/a_{22}{a}_{92}=(c_{29}-a_{21}{a}_{91})/a_{22}$$ $$a_{33}=\sqrt{c_{33}-a_{31}^2-a_{32}^2}{a}_{43}=(c_{34}-a_{31}{a}_{41}-a_{32}{a}_{42})/a_{33}{a}_{53}=(c_{35}-a_{31}{a}_{51}-a_{32}{a}_{52})/a_{33}{a}_{63}=(c_{36}-a_{31}{a}_{61}-a_{32}{a}_{62})/a_{33}{a}_{73}=(c_{37}-a_{31}{a}_{71}-a_{32}{a}_{72})/a_{33}{a}_{83}=(c_{38}-a_{31}{a}_{81}-a_{32}{a}_{82})/a_{33}{a}_{93}=(c_{39}-a_{31}{a}_{91}-a_{32}{a}_{92})/a_{33}$$ $$a_{44}=\sqrt{c_{44}-a_{41}^2-a_{42}^2-a_{43}^2}{a}_{54}=(c_{45}-a_{41}{a}_{51}-a_{42}{a}_{52}-a_{43}{a}_{53})/a_{44}{a}_{64}=(c_{46}-a_{41}{a}_{61}-a_{42}{a}_{62}-a_{43}{a}_{63})/a_{44}{a}_{74}=(c_{47}-a_{41}{a}_{71}-a_{42}{a}_{72}-a_{43}{a}_{73})/a_{44}{a}_{84}=(c_{48}-a_{41}{a}_{81}-a_{42}{a}_{82}-a_{43}{a}_{83})/a_{44}{a}_{94}=(c_{49}-a_{41}{a}_{91}-a_{42}{a}_{92}-a_{43}{a}_{93})/a_{44}$$ $$a_{55}=\sqrt{c_{55}-a_{51}^2-a_{52}^2-a_{53}^2-a_{54}^2}{a}_{65}=(c_{56}-a_{51}{a}_{61}-a_{52}{a}_{62}-a_{53}{a}_{63}-a_{54}{a}_{64})/a_{55}{a}_{75}=(c_{57}-a_{51}{a}_{71}-a_{52}{a}_{72}-a_{53}{a}_{73}-a_{54}{a}_{74})/a_{55}{a}_{85}=(c_{58}-a_{51}{a}_{81}-a_{52}{a}_{82}-a_{53}{a}_{83}-a_{54}{a}_{84})/a_{55}{a}_{95}=(c_{59}-a_{51}{a}_{91}-a_{52}{a}_{92}-a_{53}{a}_{93}-a_{54}{a}_{94})/a_{55}$$ $$a_{66}=\sqrt{c_{66}-a_{61}^2-a_{62}^2-a_{63}^2-a_{64}^2-a_{65}^2}{a}_{76}=(c_{67}-a_{61}{a}_{71}-a_{62}{a}_{72}-a_{63}{a}_{73}-a_{64}{a}_{74}-a_{65}{a}_{75})/a_{66}{a}_{86}=(c_{68}-a_{61}{a}_{81}-a_{62}{a}_{82}-a_{63}{a}_{83}-a_{64}{a}_{84}-a_{65}{a}_{85})/a_{66}{a}_{96}=(c_{69}-a_{61}{a}_{91}-a_{62}{a}_{92}-a_{63}{a}_{93}-a_{64}{a}_{94}-a_{65}{a}_{95})/a_{66}$$ $$a_{77}=\sqrt{c_{77}-a_{71}^2-a_{72}^2-a_{73}^2-a_{74}^2-a_{75}^2-a_{76}^2}{a}_{87}=(c_{78}-a_{71}{a}_{81}-a_{72}{a}_{82}-a_{73}{a}_{83}-a_{74}{a}_{84}-a_{75}{a}_{85}-a_{76}{a}_{86})/a_{77}{a}_{97}=(c_{79}-a_{71}{a}_{91}-a_{72}{a}_{92}-a_{73}{a}_{93}-a_{74}{a}_{94}-a_{75}{a}_{95}-a_{76}{a}_{96})/a_{77}$$ $$a_{88}=\sqrt{c_{88}-a_{81}^2-a_{82}^2-a_{83}^2-a_{84}^2-a_{85}^2-a_{86}^2-a_{87}^2}{a}_{98}=(c_{89}-a_{81}{a}_{91}-a_{82}{a}_{92}-a_{83}{a}_{93}-a_{84}{a}_{94}-a_{85}{a}_{95}-a_{86}{a}_{96}-a_{87}{a}_{97})/a_{88}$$ $$a_{99}=\sqrt{c_{99}-a_{91}^2-a_{92}^2-a_{93}^2-a_{94}^2-a_{95}^2-a_{96}^2-a_{97}^2-a_{98}^2}.$$

References

All Tables

All Figures

Fig. 1 Distributions of equatorial obliquity and cometocentric longitude of the Sun at perihelion for 51P/Harrington. Top left panel: histogram of equatorial obliquity, ɛ51P/Harrington, with bins computed using the Freedman and Diaconis rule, while top right panel: uses the Bayesian Blocks technique. Similarly, the cometocentric longitude, ϕ51P/Harrington, histograms displayed in the lower panels also use the Freedman and Diaconis rule (left) and Bayesian Blocks (right). Histograms are based on data from Table 4 (see the text for details).

In the text

	Fig. 2 Distribution of thermal lag angle for 51P/Harrington. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 4.
In the text

	Fig. 3 Distribution of values of the oblateness of the nucleus for 51P/Harrington. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 4.
In the text

	Fig. 4 Distributions of spin-axis orientations, (λp, βp), for 51P/Harrington. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 4.
In the text

Fig. 5 Gaussian kernel density estimation of the spin-axis orientations computed using the SciPy library (Virtanen et al. 2020). Top panel: (ɛc, ϕc) map. The results in ecliptic coordinates, (λp, βp), are shown in the middle panel and those in equatorial coordinates, (αp, δp), are displayed in the bottom panel that includes the actual observations (in red) of comet 51P/Harrington available from the MPC.

In the text

Fig. 6 Distributions of equatorial obliquity and cometocentric longitude of the Sun at perihelion for 1I/2017 U1 (‘Oumuamua). Top left panel: histogram of equatorial obliquity, ɛ‘Oumuamua, with bins computed using the Freedman and Diaconis rule, while top right panel: uses the Bayesian Blocks technique. Similarly, the cometocentric longitude, ϕ‘Oumuamua, histograms displayed in the lower panels also use the Freedman and Diaconis rule (left) and Bayesian Blocks (right). Histograms are based on data from Table 1 (see the text for details).

In the text

	Fig. 7 Distribution of thermal lag angle for 1I/2017 U1 (‘Oumuamua). Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 1.
In the text

	Fig. 8 Distribution of values of the oblateness of interstellar object 1I/2017 U1 (‘Oumuamua). Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel uses the Bayesian Blocks technique. Histograms are based on data from Table 1.
In the text

	Fig. 9 Same as Fig. 6 but assuming a value for the oblateness of interstellar object 1I/2017 U1 (‘Oumuamua) equal to −6.7 as suggestedby Mashchenko (2019).
In the text

	Fig. 10 Same as Fig. 7 but assuming a value for the oblateness of interstellar object 1I/2017 U1 (‘Oumuamua) equal to −6.7 as suggestedby Mashchenko (2019).
In the text

	Fig. 11 Distributions of spin-axis orientations, (λp, βp), for 1I/2017 U1 (‘Oumuamua) assuming a very prolate shape, with p = −6.7. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 1.
In the text

Fig. 12 Gaussian kernel density estimation of the spin-axis orientations computed using the SciPy library (Virtanen et al. 2020) for a very prolate case (p = −6.7). Top panel: (ɛc, ϕc) map. The results in ecliptic coordinates, (λp, βp), are shown in the middle panel and those in equatorial coordinates, (αp, δp), are displayed in the bottom panel that includes the actual observations (in red) of interstellar object 1I/2017 U1 (‘Oumuamua) available from the MPC.

In the text

	Fig. 13 Same as Fig. 6 but assuming a value for the oblateness of interstellar object 1I/2017 U1 (‘Oumuamua) equal to 0.836 as suggested by Mashchenko (2019).
In the text

	Fig. 14 Same as Fig. 7 but assuming a value for the oblateness of interstellar object 1I/2017 U1 (‘Oumuamua) equal to 0.836 as suggested by Mashchenko (2019).
In the text

	Fig. 15 Distributions of spin-axis orientations, (λp, βp), for 1I/2017 U1 (‘Oumuamua) assuming a very oblate shape, with p = 0.836. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 1.
In the text

Fig. 16 Gaussian kernel density estimation of the spin-axis orientations computed using the SciPy library (Virtanen et al. 2020) for a very oblate case (p = 0.836). Top panel: (ɛc, ϕc) map. The results in ecliptic coordinates, (λp, βp), are shown in the middle panel and those in equatorial coordinates, (αp, δp), are displayed in the bottom panel that includes the actual observations (in red) of interstellar object 1I/2017 U1 (‘Oumuamua) available from the MPC.

In the text

	Fig. 17 Distributions of equatorial obliquity and cometocentric longitude of the Sun at perihelion for 2I/Borisov. These results correspond to the orbit determination of comet 2I/Borisov in Table 2.
In the text

	Fig. 18 Distribution of thermal lag angle for 2I/Borisov. These results correspond to the orbit determination of comet 2I/Borisov in Table 2.
In the text

	Fig. 19 Distribution of the oblateness parameter for 2I/Borisov. These results correspond to the orbit determination of comet 2I/Borisov in Table 2.
In the text

	Fig. 20 Distributions of spin-axis orientations, (λp, βp), for 2I/Borisov. Left panel: histogram with bins computed using the Freedman and Diaconis rule, while right panel: uses the Bayesian Blocks technique. Histograms are based on data from Table 2.
In the text

Fig. 21 Gaussian kernel density estimation of the spin-axis orientations of 2I/Borisov computed using the SciPylibrary (Virtanen et al. 2020) and based on data from Table 2. Top panel: (ɛc, ϕc) map. The results in ecliptic coordinates, (λp, βp), are shown in the middle panel and those in equatorial coordinates, (αp, δp), are displayed in the bottom panel that includes the actual observations (in red) of 2I/Borisov available from the MPC. Already published spin-axis directions are plotted as filled orange symbols with white error bars (when available): triangles (Ye et al. 2020), square (Manzini et al. 2020), star (Bolin et al. 2020b), and circle (Kim et al. 2020).

In the text

	Fig. 22 Distributions of equatorial obliquity and cometocentric longitude of the Sun at perihelion for 2I/Borisov. Same as Fig. 17 but for the orbit determination of comet 2I/Borisov in Table 3.
In the text

	Fig. 23 Distribution of thermal lag angle for 2I/Borisov. Same as Fig. 18 but for the orbit determination of comet 2I/Borisov in Table 3.
In the text

	Fig. 24 Distribution of the oblateness parameter for 2I/Borisov. Same as Fig. 19 but for the orbit determination of comet 2I/Borisov in Table 3.
In the text

	Fig. 25 Distributions of spin-axis orientations, (λp,  βp), for 2I/Borisov. Same as Fig. 20 but based on data from Table 3.
In the text

	Fig. 26 Same as Fig. 21 but based on data from Table 3.
In the text

	Fig. 27 Same as Fig. 26, bottom panel, but based on data from Table 5.
In the text