© ESO 2022

1 Introduction

We present the long-term dynamical evolution of two unusual long-period comets (LPCs), C/2014 UN₂₇₁ (Bernardinelli-Bernstein) and C/2017 K2 (PANSTARRS). We follow the abbreviation convention proposed by Bernardinelli et al. (2021) and for the brevity refer to these comets as BB and K2, respectively. In both cases, the first observations were taken at extremely large heliocentric distances. In addition, both comets attracted scientific attention because of their pronounced cometary activity well beyond the orbit of Saturn.

The discovery of BB was announced in June 2021 by Bernardinelli & Bernstein (2021) on the basis of astrometric observations spanning 20 nights from 2014 October 20 to 2018 November 8 (29.0–23.7 au from the Sun). A little later, at a distance of over 20 au from the Sun, observers reported the cometary activity of BB (Farnham et al. 2021; Kokotanekova et al. 2021; Bernardinelli et al. 2021). In all these papers, the authors estimate that BB might be one of the largest comets ever discovered. Backward integration of the BB orbit by Bernardinelli et al. (2021) under their assumed Galactic tidal model, and including perturbations from eight stars taken from a list of Sun–star encounters published by Bailer-Jones et al. (2018), yields a perihelion distance of ~18 au during its previous perihelion passage 3.5 Myr ago.

The comet K2 observational campaigns have been conducted from the moment of its discovery in May 2017. Undoubtedly, its activity, most of all related to the existence of CO ice, was found at heliocentric distances well beyond 20 au (Meech et al. 2017; Jewitt et al. 2021; Yang et al. 2021).

BB will pass through perihelion (10.95 au from the Sun) on 2031 January 22, whereas K2 will reach its perihelion on 2022 December 19, at 1.80 au from the Sun. Therefore, both have been observed for many years on their pre-perihelion orbital leg, and, in addition, at great distances from the Sun. As a result, they are favorable candidates for studying their past motion and their origin. Moreover, such a large perihelion distance of BB allows us to constrain the future orbital evolution of this comet on the basis of pre-perihelion data because we can assume that nongravita-tional (hereafter NG) forces will be negligible in the context of the orbital motion of this comet.

Both comets have an original barycentric semimajor axis of greater than 10 000 au, and so they belong to the so-called Oort spike. It is widely accepted that it is necessary to take into account both Galactic and stellar perturbations when investigating the long-term dynamical evolution of such elongated orbits, far outside our planetary system. A new, reliable, precise, and fast method of calculating the effect of these perturbations on LPC motion was recently proposed by Dybczyński & Breiter (2022). We use these methods in what follows.

After applying the methods mentioned above, we realized that increasing knowledge on stars visiting the solar neighborhood becomes more important in long-term dynamical studies of LPCs. In recent years, we have benefited from the Gaia astro-metric mission (Gaia Collaboration 2016). Using its latest third data release (Gaia Collaboration 2021), we show the effect of some particularly important stars passing near to the Sun on the motions of BB and K2. We also present the impact of the stellar data uncertainty on our results.

In the abstract of their paper on K2, Jewitt et al. (2021) wrote: Nongravitational acceleration in C/2017 K2 and similarly distant comets, while presently unmeasured, may limit the accuracy with which we can infer the properties of the Oort cloud from the orbits of long-period comets. However, as we discussed in Królikowska et al. (2012), uncertainties related to the NG acceleration in the motion of the comet can be eliminated or substantially minimized by using the pre-perihelion part of the data to calculate the osculating orbit. Uncertainties will be even further reduced when the pre-perihelion part of data is limited to large heliocentric distances. At this point, the original orbits will also not be burdened by the uncertainty related to the subsequent orbital change under the influence of increasing NG forces as the comet approaches the Sun.

Comets BB and K2 are ideal for such investigations based on distant data before perihelion. Astrometric measurements starting from distances exceeding 20 au are available for both objects. This guarantees that the distant parts of their data arcs are long enough to determine orbits of the highest quality class (Królikowska & Dybczyński 2018c). One of the goals of this study is to show how we can deal with good-quality orbits by limiting the data arcs for orbit determination to the properly selected part of pre-perihelion orbit. Another goal is to verify how the past motion of a particular comet fits the general picture of cometary reservoir evolution; see Vokrouhlický et al. (2019) and Dones et al. (2015) for extensive reviews.

The structure of this paper is as follows. In Sect. 2, we discuss the choice of appropriate data arcs to obtain the starting orbits to study the origin of both comets. Our current state of knowledge on the potential stellar perturbers of LPC motion is briefly described in Sect. 3 while in Sect. 4 we present our methods of dealing with the orbital parameter uncertainties for comets and stars in our long-term dynamical studies. We also present data on all the particular stars mentioned in this paper. Section 5 describes BB past orbit evolution and its parameters at the previous perihelion, while Sect. 6 offers the same for the future BB dynamics. We discuss changes in past orbit of K2 in detail in Sect. 7. Section 8 contains a summary and our conclusions.

2 Orbit determination for both comets

We present a detailed description of the positional data used in this paper for both comets in Table 1. This table also includes values of an original 1/a. The model of motion and the methods used for osculating orbit determination are described in Królikowska & Dybczyński (2017) and references therein.

2.1 The orbit of BB

At the end of September 2021, comet BB was about 20 au from the Sun, and we took the entire set of astrometric data available at the IAU Minor Planet Center¹. We supplement this data with the single precovery positional detection of Bernardinelli et al. (2021). Thanks to this single precovery measurement, the data arc of BB is extended by about 3.7 yr, going back to 2010 November 15. Additionally, Bernardinelli et al. (2021) suggested some corrections of positional data, and we applied them according to Table 2 in their paper. We find that the precovery measurement and the corrections to other positions result in only a slight change in the orbital elements; for example, the value of the original 1/a changed by less than 1%. We conclude that the orbit of BB is now quite well constrained.

As expected, we find no evidence of NG effects in the motion of BB, which is still beyond the orbit of Saturn. Table 1 also shows the solution “i1” based on a data arc that is 3 months shorter and has almost two times fewer positional measurements than the data arc used in the case of the “b8” solution. The two 1/a_ori values do not differ in statistical terms (at one sigma level). Furthermore, based on a data arc that is almost 1 year shorter than that used for solution “i1” Bernardinelli et al. (2021) obtained the value of 49.50 × 10⁻⁶ au⁻¹ for 1/a_ori, which indicates the high compatibility of the 1/a_ori as a function of the increasing data arc taken for the orbit determination.

2.2 Orbit of K2 for the backward dynamical evolution

For K2, the situation is much more interesting in the context of detecting NG effects in the motion of this comet. Therefore, it is more complex from the perspective of the dynamical study of the origin of K2.

Our first published attempt to study the past dynamics of K2 was based on the 4.6-yr data arc (Królikowska & Dybczyński 2018a, hereafter KD18a). Almost a year later, we updated our calculations based on a longer interval of observations and published our results as an addendum to our preprint available at arXiv (Królikowska & Dybczyński 2018b, hereafter KD18b). These two orbital solutions were presented as Solution A1 in KD18a and Solution A1-new in the addendum of KD18b. In the present paper, we refer to these as “a8” and “a9,” respectively; see Table 1. Here, we analyze whether or not a longer observational arc will improve the reliability of the original K2 orbit. We decided to include three additional solutions in Table 1, “a6,” “c5,” and “b5,” based on longer data arcs.

We find that some traces of NG effects can be seen in the longer data arcs considered here; see solution “c5.” For this orbit, we obtained the following NG parameters: A₁ = 28.414 ± 4.968, A₂ = 11.610 ± 1.089, and A₃ = 4.5795 ± 0.7221, in units of 10⁻⁸ au day⁻², where the g(r)-like function reflects CO-driven sublimation and was applied as described in Królikowska (2020). Fortunately, we found that this NG orbit gives a very similar original 1/a_ori to that in the purely gravitational case “a6” based on the same interval. In addition, the NG solution for the longest data arc (the same as in the case of the GR solution marked “b5”) results in a small and uncertain radial component of NG acceleration (using the CO-driven formula) and allowing for NG effects does not eliminate some trends visible in [O-C]; and so it is not considered here. By the analysis of many NG and GR orbits, we finally concluded that the best solution for studying long-term past dynamics of K2 is the solution “a9”. Thus, solutions based on longer data arcs (“a6”, “c5”, and “b5”) are only used here to identify the extent to which longer data arcs can change our inference about the dynamical status of K2. Choosing the “a9” solution is a rather cautious approach, because the uncertainty on 1/a_ori is two times greater than for the “c5” solution.

2.3 Future orbital evolution predictions

Today we do not know how the NG acceleration will change the current orbit of K2 when the comet passes the perihelion at a distance of 1.8 au from the Sun in December 2022. Therefore, its future dynamics can only be discussed when the comet is on its post-perihelion trajectory. The situation is different in the case of BB, because its perihelion distance is equal to almost 11 au. Thus, it seems that predictions made now will be correct also in the aspect of a future dynamical evolution for this comet.

3 Stellar perturbers

Our knowledge of stars that can perturb the motion of LPCs has expanded in recent years. Wysoczańska et al. (2020) have introduced a publicly available StePPeD database² containing data on potential stellar perturbers of such LPC motion. This database was revised several times in 2020, which finally led to StePPeD release 2.3, which is mostly based on the Gaia DR2 catalog (Gaia Collaboration 2018). However, after the Gaia Early Data Release 3 was made available (Gaia Collaboration 2021) it became clear that the next substantial update of the StePPeD database was necessary.

The first step of this upgrade was finished in September 2021, and a new StePPeD release 3.0 was announced. This first step was limited to the same list of stars as in version 2.3, but for almost all of them, the new data from the Gaia EDR3 catalog have been incorporated. This work is in progress to add new stellar perturber candidates.

We note here that a comparison of the data in DR2 and EDR3 sometimes revealed substantial changes. Almost 25% of the stars in our list have a new parallax that differs by more than 100% from the value previously provided (i.e., that in DR2); nine of them have a negative parallax in EDR3. Over one-third of the stars in our database have a RUWE parameter of greater than 1.4, which according to the Gaia documentation indicates possibly less accurate data (Fabricius et al. 2021). For 5% of our stars, there is no parallax or proper motion in EDR3 and we still have to rely on the Gaia DR2 data.

In November 2021, radial velocities of a dozen stars were updated, resulting in the release 3.1 of the StePPeD database (see footnote 1); for details see the “Changelog” available at the StePPeD Internet page. All calculations presented in this study are based on this latest version of the StePPeD database.

During our study of the long-term dynamical evolution of BB and K2, we recognized very strong perturbations caused by stars HD 7977 and Gliese 710 (named P0230 and P0107 in the StePPeD database, respectively) and weaker but potentially important perturbations from five other stars. The details for all seven stars are shown in Table 2, including their identifiers, parameters of the closest approach to the Sun, and mass estimates. We describe the minimal Sun–star distance in two different ways. The value presented in the fourth column, “mindist,” is calculated in a special way: this is the distance from the Sun to the centroid of a cloud of 10 000 star clones drawn according to the data uncertainties and using the respective covariance matrix. As the uncertainty on the mindist parameter, here we present a radius of a sphere around that centroid, which includes 90% of the star clones. Uncertainties greater than the mindist directly indicate that the clone cloud surrounds the Sun (see Col. 4 for P0230). In the fifth column, we present a formal statistical description of the Sun–star distance set for all stellar clones using three percentiles: 5% (p05), the median, and 95% (p95). These are also expressed in parsecs.

We present a graphic comparison of the stellar data uncertainties and their minimal distances to the Sun in Fig. 1. For each star listed in Table 2, we plot the positions of its 10 000 clones at their closest approach to the Sun, projected on the plane of the maximum scatter; for details of such calculations; see Dybczyński & Berski (2015). The seven different plots of the spread of star clones were then merged into a single image. Cen-troids of all clouds are aligned along the horizontal line, keeping the correct distance from the Sun and maintaining the same scale in the spread of clone clouds. The whole clone cloud of P0107 (C in Fig. 1) is hidden under its centroid black dot because of the extremely precise data for this star. It should be stressed that the seven stellar close approaches to the Sun presented in Fig. 1 happened or will happen in different epochs spread over an interval of 6 Myr; see the sixth column of Table 2.

There is an important qualitative difference between the per-turbative action of P0230 and P0107 and the remaining five stars mentioned above. According to the most recent data, approaches of these two stars to the Sun are so close that they noticeably perturb its galactic trajectory. As a result, these perturbations affect the motions of all Solar System bodies. The resulting heliocentric orbit changes depend mostly on the heliocentric velocity of a given body. LPCs, when far from their perihelion and therefore moving very slowly, are the best candidates for the largest change in their perihelion distance. We show this effect in the following sections.

Fig. 1 Comparison of nominal minimal distance from the Sun and spread of star clones at their closest approach for seven stars mentioned in this paper. Black dots in the center of each clone cloud mark the nominal star position at the closest Sun–star approach; see text on how this composite figure was obtained. All 10 000 clones of star C are hidden under its black dot because the stellar data are extremely precise in this case. Label meanings: A: P0509, B: P0230 (HD 7977), C: P0107 (Gliese 710), D: P0506, E: P0508, F: P0417 (Ton 214), G: P0111 (HIP 94512); see Table 2 for details of each star. The blue circle marks the habitually adopted outer limit of the Oort cloud.

4 Dealing with the uncertainties

In the next three sections, we describe the past and future motions of BB and the past motion of K2 in detail. For both comets, we also estimated the influence of their orbital uncertainty on a previous or subsequent orbit for all solutions presented in Table 1. In addition, for the preferred solutions, we performed simulations to observe the effects of stellar uncertainties in various cases. To this purpose, we extensively use the methods proposed recently by Dybczyński & Breiter (2022) and the stellar data from the latest release 3.1 of the StePPeD database. In determining an osculating orbit, we obtain the covariance matrix that allows us to construct comet orbit clones satisfying the observational constraints.

Dealing with the cometary orbit uncertainty is relatively straightforward. During an osculating orbit determination, we obtain the covariance matrix, which allows us to draw comet orbit clones satisfying the observational constraints; details of this procedure can be found in Sitarski (1998). In all simulations presented in this paper, we use 5000 comet clones plus a nominal orbit. We numerically propagated each clone in this set back and forth to obtain original and future swarms of barycentric orbit clones. Investigating the previous (or next) orbit, we repeat the backward (or forward) numerical integration for each comet clone exactly as described in Dybczyński & Breiter (2022), using the latest stellar ephemerides from the StePPeD site. In all cases, we take into account both the overall Galactic gravitational potential and stellar perturbations from a set of 407 potential perturbers. Due to the use of advanced algorithms, this calculation is precise and very fast; for 5001 comet orbits, we obtain the results in 15–20 min on a standard PC. For the sake of comparison, we also present the results of calculations for comet clones with all stellar perturbations excluded.

Analyzing the effect of stellar data uncertainties requires a more complicated and much more time-consuming approach. Each individual calculation consists of two stages. First, we generate a tailored stellar ephemeris and then use a standard code (described above) for a comet motion propagation using this individualized stellar ephemeris. To obtain this special ephemeris, we replace the nominal data for the studied star with its clone drawn from the respective six-dimensional covari-ance matrix from Gaia EDR3³. We then generate the stellar ephemeris using nominal data for the remaining stars. In this way, we still take into account all possible star–star interactions, which are very rare but happen. In all simulations presented in this paper, we use 10 000 star clones for the star in question plus the nominal star tracks for the rest of them. The tailored stellar ephemeris can be produced for a shorter time than the publicly available one (30 Myr back and forth), but it still requires considerable CPU time, typically about 20 s. At the second stage, the subsequent comet motion integration takes only milliseconds, but a whole investigation consisting of 10 000 cases takes over 50 h on a single CPU core.

We need the mass of a star in order to calculate the effect of its perturbation on a comet motion. In the StePPeD database, there are mass estimates included together with their uncertainties and we used these values. These data are taken from numerous different sources. However, the uncertainties on the mass estimate were obtained by various considerably different methods. Their nature is highly nonlinear, often asymmetric, and with unknown error distribution. Therefore, the statistical interpretation of the mass estimate uncertainties is difficult; we decided not to draw random mass values from these data, and in all calculations we utilize the nominal mass estimate.

However, we should keep in mind that the mass uncertainty is an additional source of the approximate nature of our results. For example, the mass estimate for the most important star, P0230, is taken from Stassun et al. (2019) and is described as 1.080 M_⊙ with a symmetric uncertainty of ± 0.136 M_⊙. The nominal value of the previous perihelion distance of K2 is 441 au (for the “a9” solution) but if we use a lower limit mass for P0230, that is, 0.944 M_⊙, we obtain a lower value of q_prev = 331 au.

Using the above procedure, we can theoretically investigate the simultaneous effect of two or more star uncertainties without additional computational cost, because the process of drawing a star clone is very fast. However, the spread of the results grows considerably. Therefore, it is necessary to calculate a much larger number of cases to achieve a statistically valuable result. Fortunately, in all cases described in the following sections, only one star strongly perturbs a comet motion during the previous or next orbital period.

Fig. 2 Past and future dynamical evolution of C/2014 UN271 nominal orbit (“b8” solution). changes in a perihelion distance (green), an inclination (fuchsia), and an argument of perihelion (red) are shown. the thick lines depict the result of the full dynamical model while thin lines show the evolution of elements in the absence of any stellar perturbations, i.e., only the galactic perturbations are taken into account. angular elements are expressed in a galactic frame. we also show names of stars that have a significant effect on this dynamical evolution.

5 C/2014 UN₂₇₁ past orbit evolution

Starting from the up-to-date osculating orbit of BB (solution “b8” in Table 1), we calculated the original and future orbits at a distance well outside the planetary perturbation zone (as usual, we used a distance of 250 au from the Sun). The elements of these orbits are presented in Tables A.1 and A.2. Using these original and future orbital elements, we calculated the dynamical evolution of the nominal BB orbit for one orbital period to the past and to the future. This is presented in Fig. 2. In the past motion, we see several small perturbations from passing stars and one very strong orbit change. The two small interactions were with P0509 at –0.67 Myr and P0506 at –1.08 Myr; see Table 2 for more information on these stars. However, the most prominent perturbation was caused 2.47 Myr ago by the star P0230. As shown in Fig. 2, the nominal previous perihelion distance is equal to 238 au and BB was at the previous perihelion nominally 2.71 Myr ago. If we exclude all stellar perturbations (the thin lines in Fig. 2), then the previous perihelion distance is found to be much smaller, only 15 au.

It should be stressed that the dynamical evolution depicted in Fig. 2 is based on the nominal BB data and nominal data for all 232 stars included in the dynamical model of the past comet motion. To obtain a more realistic picture, we should estimate the influence of the cometary and stellar data uncertainties and add the results to the above discussion. In the first step, we calculated the effect of the BB orbit uncertainty using the methods described in Sect. 4. The result of this numerical experiment is summarized in Fig. 3. The previous perihelion distance appears to be very close to the nominal value and can be described as 237.6 ± 32.9 au, because its distribution can be quite well approximated with the Gaussian one; see also the BB-prev-B row in Table 3.

As is clearly shown in Fig. 2, apart from a series of weak stellar perturbations, the strongest one is caused by P0230. This perturbation is also important because it is a perturbation of the motion of the Sun, and so it acts on the comet indirectly and does not depend on the star–comet distance. Instead, it strongly depends on the minimal distance of the Sun–star encounter and its geometry. In Fig. 1 we present the effect of P0230 uncertainties (Star B) on these parameters, drawing 10 000 clones of this star and stopping their motion at the closest approach to the Sun. From this picture, we can see that the minimum distance can be arbitrarily small and that the direction of the impulse that acts on the Sun is unknown.

To observe the possible effect of P0230 uncertainties on the past BB motion, we performed a calculation in which we used 10 000 pairs: a P0230 clone and a BB clone as described in Sect. 4. The result is shown in Fig. 4. The spread of the previous orbital parameters is considerable, and to make this plot readable, we have to omit a very long tail of points placed to the right. This omitted set of points consists of 27 cases of negative 1/a_prev and 87 large values for elliptic orbits. For the sake of comparison, we superimposed the swarm of blue dots over the main plot, which show the distribution of the previous BB orbit when only the nominal star P0230 acts on 5001 comet clones. This blue set of points is the same as that shown in black in Fig. 3.

The smallest 1/a_prev value obtained in this numerical experiment is equal to –558.7 in the same units as in the plot. The interval of q_prev spreads from 0.021 to 20908.3 au. The distribution of q_prev is strongly non-Gaussian and we can describe it with three deciles, 10th, 50th (median), and 90th: 12.5–49.6–252.4 au.

In Fig. 4, it can be seen that the nominal q_prev = 237.74 au is far from the maximum of the q_prev distribution. The reason for this comes from the geometry of the cloud of P0230 clones with respect to the Sun (star B in Fig. 1). The nominal minimum distance between P0230 and the Sun is 0.014 pc. However, if we ignore the geometry and analyze the one-dimensional minimal distance distribution, the resulting median is equal to 0.032 au (see the fifth column in Table 2), and over 87% of the values are greater than the nominal one. As a result, the much weaker perturbation producing a smaller q_prev is much more probable.

The conclusion from this section is clear: Until much more precise data for P0230 are available, we will be unable to describe the BB orbit at its previous perihelion in a definitive manner. For this reason, we look forward to the next Gaia data release, which is announced for June 2022.

Fig. 3 Influence of BB orbit uncertainty on parameters of the previous orbit. black dots show the results of the backward numerical integration of 5001 clones of BB stopped at the previous perihelion with all Galactic and stellar perturbations included. The center of the green circle marks the nominal result. The red histogram presents a marginal distribution of qprev along with the best Gaussian approximation (green dots in the middle of each bar). The blue histogram is a marginal distribution of 1/aprev.

Fig. 4 Influence of P0230 data uncertainty on parameters of the previous BB orbit – in this simulation we used 10 000 pairs: a comet clone and a P0230 clone. The result of each pair is represented by a black dot. Omitted are 114 points from a long right tail, among which 27 have negative 1/aprev. The marginal distribution of qprev is performed for all black points and shown as the red histogram below the main plot. For the sake of comparison, the main result is superimposed, with the blue swarm of dots showing the result of an additional experiment: 5001 BB orbits and nominal P0230 track. The nominal result is in the center of a green circle.

6 C/2014 UN₂₇₁ future orbit evolution

In the future BB motion part of Fig. 2, we can observe a strong perturbation caused by P0107 and a barely visible, weaker one by P0111. Similarly to the past BB motion, the strongest perturbation is an indirect one, because it results from a very close passage of P0107 near the Sun in 1.29 Myr. However, in this favorable case, the prediction of the future is much easier than the prediction of the past. This comes from the very precise data available for P0107. As mentioned above, the whole P0107 cloud of 10 000 clones in Fig. 1 is hidden under the black dot that marks the star’s nominal position at the closest encounter with the Sun.

Using the methods described in Sect. 4 we performed three numerical simulations to observe the effect of cometary and stellar data uncertainties on the parameters of the next BB orbit: BB-next-B: all stars on their nominal tracks acting on each of the 5001 BB clones; BB-next-C: the same star set but in each run P0107 is replaced by one of 10000 clones, acting on the nominal BB orbit; and BB-next-D: 10000 pairs of P0107 and BB clones. The numerical results are presented in Table 3. As the influence of cometary orbit uncertainties appeared comparable to that of stellar data, we summarized all three simulations in one composite plot shown in Fig. 5. We plotted the results of simulation BB-next-D as black dots superimposed with orange dots from BB-next-C, and the results of simulation BB-next-B as blue dots on top of the previous two sets. As is clearly shown in Fig. 5, stellar perturbations only slightly enlarge the interval of q_next values resulting from the BB orbital uncertainty. As was stated at the end of Sect. 3, the perturbation of P0107 on BB is only indirect, resulting from an impulse gained by the Sun. This is well illustrated by the fact that in our BB-next-D simulation, the smallest distance between a BB clone and a P0107 clone was over 55 000 au.

For the future BB motion, our conclusion is that its next perihelion could be described by deciles (9.13–11.06–13.50) au, that is, with a median value almost identical to the nominal one (11.11 au). The semimajor axis of the next BB orbit is larger than for the current apparition, mainly due to the planetary perturbations. Our simulation BB-next-D gives 1/a_next = (35.40–36.16–36.90) ×10⁻⁶ au⁻¹ which corresponds to a semimajor axis of about 27 700 au and an orbital period of 4.6 Myr.

If one compares the original and future BB orbits (see Tables A.1 and A.2), it can be noticed that while planetary perturbations will change the semimajor axis of this comet orbit by almost 30%, its perihelion distance will remain almost unchanged. After passing the planetary zone, the orbital period of BB will increase to 4.6 Myr, and consequently this comet will reappear among planets.

Fig. 5 Composite picture of the influence of cometary and stellar uncertainties on the BB orbit recorded at the next perihelion passage. Blue dots mark 5001 clones of BB, each perturbed by all stars on their nominal tracks. Orange dots represent the nominal BB orbit perturbed by 10 000 clones of P0107. Black dots present 10 000 pairs of the interacting BB and P0107 clones. The red marginal distribution histogram of qnext corresponds to the black cloud of points. The result for nominal stellar and cometary data is in the center of the green circle.

7 Past dynamical evolution of the orbit of C/2017 K2

As explained in detail in Sect. 2, we study only the past longterm dynamical evolution of K2 and base our calculations on the original orbit obtained using the solution “a9.” We present the nominal past K2 orbit evolution in Fig. 6. Here we notice a series of moderate stellar perturbations and a very strong perturbation caused by P0230, 2.47 Myr ago. This solution gives 1/a_ori = 35.48 ± 2.33 × 10⁻⁶ au⁻¹, which corresponds to the previous orbital period of 4.73 Myr. As a result, the indirect perturbation from P0230 happened almost exactly when the comet was at its aphelion and was moving very slowly.

In Fig. 6, it is possible to see local orbit changes due to the action of P0509 at –0.55 Myr and P0506 at –1.08 Myr, as in the case of BB. In addition, small orbit changes resulting from a passage of P0508, 0.98 Myr ago, and P0417, 1.47 Myr ago, can also be observed. More information about these stellar perturbers can be found in Table 2. The accumulated stellar perturbations increased the nominal K2 previous perihelion distance from 12 au obtained when only Galactic potential is taken into account (thin lines in Fig. 6) up to 441 au. We verified numerically that, in the absence of P0230, all other stellar perturbations cancel each other and the resulting previous perihelion distance would be almost unaffected by stars.

As shown above, the nominal previous perihelion is only a highly approximate qualitative result. We must elucidate the possible influence of cometary and stellar data uncertainties on the final result. The influence of the uncertainty related to the K2 orbit (for the solution “a9”) on our results is presented in Fig. 7. We used the methods described in Sect. 4. The numerical results of this calculation are included in Table 4. In Fig. 7, we also present marginal distributions of 1/a_prev (blue) and q_prev (red). As the best estimate, we obtained from these calculations: 1/a_prev = (34.88 ± 2.61) × 10⁻⁶ au⁻¹ and q_prev = (373.5–441.6–490.0) au. The latter is in the form of three deciles because its distribution is not Gaussian. It is easy to note that the range of the previous perihelion distance is well defined here, and its nominal and median values are close to each other. We performed similar numerical experiments for the remaining orbital solutions and the results can be reviewed in Table 4.

From Fig. 6 one can read that during the K2 backward numerical integration, from the listed stars, we first meet P0509 (star A in Fig. 1). The data uncertainties for this star are rather large. It has no parallax or proper motion data in Gaia EDR3, which probably causes problems in the observational data reduction, and so we use all the astrometry for this star from DR2. We decided to check the importance of this star in our calculations. We draw 10 000 clones of P0509 to observe the spread of the result when calculating past K2 orbit parameters at the previous perihelion. The smallest star–comet distance in this experiment was 765 au. Results of this simulation are presented in Fig. 8 and Table 5 as K2-prev-C. As this stellar perturbation is direct, that is, P0509 has a strong effect on a comet because it passes relatively close to it, we also analyzed the star–comet distance distribution. We use two different colors for dots representing the individual clone results: orange dots for clones with a minimal distance from K2 smaller than 20 000 au and black for more distant ones. Figure 8 shows that the dependence of the strength of the perturbation on the star–comet distance is not obvious because it also strongly depends on the star–comet–Sun geometry. Moreover, some P0509 clones passed much closer to the Sun than the nominal one, and therefore an indirect perturbation is also possible. We omitted 15 cases as the extreme outliers. The smallest q_prev is equal to 3.6 au while the omitted largest one is equal to 9510 au. Despite such a large spread, the q_prev distribution shown as a red histogram in Fig. 8 is quite compact. It can be described by three deciles as: (407.1–447.6–489.1) au. Still, we are close to the nominal value of q_prev.

The situation is quite different when we investigate the impact of the uncertainties of P0230 on our results. We performed two additional simulations using the techniques presented in Sect. 4, labeled K2-prev-D: a comet in nominal orbit and 10 000 clones of P0230, and K2-prev-E: 10 000 pairs of star and comet clones. Numerical results from all simulations performed for K2, and those for orbital solutions other than “a9,” are presented in Tables 4 and 5. The results based on the “a9” solution are summarized in a composite plot shown in Fig. 9. We plotted the results of K2-prev-E as black dots overprinted with orange dots from K2-prev-D. For the sake of comparison, we added the results of K2-prev-C (the same as presented in Fig. 7) as blue dots on top of the previous two sets.

Several interesting features can be seen in Fig. 9. The K2 orbit uncertainties are responsible mainly for the spread in 1/a_prev (blue dots). The swarm of P0230 clones, when acting on a nominal K2 orbit, causes mainly the spread in q_prev (orange dots) across a wide range of values from 0.0005 au to 19 350.4 au. The third simulation, K2-prev-E (black dots), spread the results in both elements. Here, we observe a similar effect to that in Fig. 4: the maximum of the q_prev distribution occurs for a substantially smaller value of the previous perihelion distance than that for the nominal orbit. The reason is also the same: due to the geometry of the P0230 swarm of random clones, the nominal star–Sun distance is much smaller than the most probable value. As a result, the weaker perturbation is more frequent than that for the nominal orbit of K2 and a star track.

Despite the general indirect nature of the influence of P0230 on K2, there are several star clones in K2-prev-E that pass close to K2. The minimum star–comet distance in this numerical experiment is 7707 au. Such close passages produce outlier results with negative 1/a_prev values; the smallest value is –621 in the units used in the plot. When preparing Fig. 9 we omitted 226 points, 27 of them with negative 1/a_prev. The largest value of the right tail of the black swarm is q_prev = 59 312.6 au.

Based on the “a9” solution, the K2-prev-E simulation presents the best of our knowledge on the previous K2 orbit at the time of writing. We obtained q_prev described by deciles: (5.163–144.0–1065) au, which, contrary to our previous opinion, suggests that K2 is probably a dynamically new comet. However, taking into account the sensitivity of this result to the P0230 data change, we should wait for much more precise measurements for this star before coming to any definitive conclusions.

Fig. 6 Dynamical past evolution of the nominal orbit of C/2017 K2 based on the “a9” orbital solution. Several individual stellar perturbations are marked. The meanings of the colors and line thicknesses are the same as in Fig. 2.

Fig. 7 Influence of K2 orbit uncertainty (for solution “a9”) on the parameters of its previous orbit. The number of points, color, and symbol meanings are the same as in Fig. 3.

Fig. 8 Influence of P0509 data uncertainty on the parameters of the previous K2 nominal orbit (solution “a9”). Black dots mark 8827 cases where the distance between K2 and P0509 is greater than 20 000 au while the 1149 orange dots show the opposite cases. Two black and 23 orange points are omitted as the extreme outliers: the maximum qprev in this simulation is equal to 9510 au. The red marginal qprev distribution was made for all 9976 plotted points.

Fig. 9 Composite picture of the influence of cometary and stellar uncertainties on the K2 orbit recorded at the previous perihelion passage. Blue dots mark 5001 K2 clones perturbed by all stars on their nominal tracks. Orange dots represent the nominal K2 orbit perturbed by 10 000 P0230 clones. Black dots present 10 000 pairs of K2 and P0230 clones. The red marginal distribution histogram of qprev corresponds to the black cloud of points. The result for nominal stellar and cometary data is in the center of the green circle. Omitted are 226 outlier points from the far right tail of the black swarm and 219 from the orange one.

8 Summary and conclusions

We performed an extensive study of the past and future motion of comet C/2014 UN₂₇₁ (Bernardinelli-Bernstein) and the past motion of C/2017 K2 (PANSTARRS). For each comet, we obtain a series of osculating orbits based on data arcs of different lengths. In the case of BB, we did not find detectable effects related to the NG acceleration caused by the sublimation of ices from the comet surface. In K2 motion, measurable effects of NG acceleration have now been noted for the data arc covering large heliocentric distances down to 7.8 au (solution c5 in Table 1). From all orbits listed in Table 1, we selected the most appropriate ones for a dynamical study of each comet motion outside the planetary zone. In the case of BB, the most appropriate solution is “b8,” based on all observations available at the time of calculation (October 2021), and for K2, the best solution is “a9,” based on the longest data arc without the visible NG effect on the orbit. However, we also performed simulations for all remaining solutions given in Table 1 for comparison purposes.

The analysis of a series of K2 orbital solutions based on different data arcs clearly shows that in the case of LPCs discovered at large heliocentric distances, the starting orbits for the study of their origin should be determined from the pre-perihelion leg of the orbit and limited to large distances from the Sun. Moreover, in cases such as K2, it is necessary to individually balance the benefits of limiting the action of NG acceleration (the farther from the Sun the better) and the quality of the obtained orbit (the longer the data arc the better). This approach, extended to LPCs discovered at very large heliocentric distances, confirms our previous published conclusions (Królikowska et al. 2012; Królikowska 2020).

Based on the selected preferred orbital solutions, we obtained a nominal previous perihelion distance value. For K2, nominal q_prev = 441 au and for BB nominal q_prev = 237 au, taking into account the perturbing effect of the full Galactic tidal field and all currently known potential stellar perturbers. The influence of orbital uncertainties on the preferred comet orbital solutions is quite small in both cases and resulted in q_prev = (373.5–441.6–490.0) au for K2 (the 10th, median, and 90th percentiles are used here) and q_prev = 237.6 ± 32.9 au for BB. At this stage, we were able to conclude that both K2 and BB are dynamically new comets. However, as shown above, the influence of the stellar data uncertainties on our results must be considered. For the past motion of K2, we separately analyzed the effects caused by stars P0509 and P0230. For the past motion of BB, we also experimented with P0230, and for its future motion, we analyzed the effect of P0107 uncertainties. Details of the above-mentioned results and several numerical experiments with the stellar perturbations can be found in Tables 3, 4, and 5 and the corresponding figures.

The list of stars that noticeably change the orbit evolution of BB and K2 is presented in Table 2 together with their parameters of the close Sun–star approaches and their mass estimates, based on release 3.1 of the StePPeD database. The past motion of both studied comets is dominated by a perturbation caused by the star P0230 (HD 7977). Our most advanced simulations, K2-prev-E and BB-prev-D, resulted in significant changes in the previous perihelion distance values obtained: q_prev = (5.163–144.0–1065.) au for K2 and q_prev = (12.50–49.62–252.4) au for BB. Looking at the median values, we can still classify both comets as dynamically new. However, the observed high sensitivity of these results to the uncertainty on the P0230 data means that this conclusion should be taken as preliminary and uncertain. We hope to receive more precise data for P0230.

It is worth mentioning that due to the specific nature of perturbations caused by P0107 and P0230 they can be the significant perturbers of many LPCs. This effect comes from their very close passage near the Sun, which generates a strong velocity impulse that indirectly affects all Solar System bodies. As a result, many heliocentric small body orbits have been or will be changed. The strength of this perturbation depends mainly on the small body heliocentric velocity. Therefore, LPCs at their aphe-lia are the best candidates to track important changes in their orbits.

Our numerical Monte Carlo simulations of the impact of stellar data uncertainties on the long-term evolution of comet orbits clearly show that the accuracy of the stellar data is still insufficient in many cases. On the contrary, the accuracy of the orbits of the contemporary observed comets seems to be satisfactory. There are seven stars that significantly perturb the motion of the two comets studied here, but we only have satisfying 6D data for two of them. We hope that future observational attempts —for example, the next Gaia data release – will improve this situation.

Acknowledgements

This research has made use of positional data of analyzed comets provided by the International Astronomical Union’s Minor Planet Center. This research has also made use of the SIMBAD database, operated at CDS, Strasbourg, France, and the VizieR catalog access tool, CDS, Strasbourg, France (DOI: 18.26893/cds/vizier). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The calculations which led to this work were in some part performed with the support from the project “GAVIP-GC: processing resources for Gaia data analysis” funded by European Space Agency (4000120180/17/NL/CBi).

Appendix A Orbital data for C/2014 UN₂₇₁ (Bernardinelli-Bernstein)

Appendix B Orbital data for C/2017 K2 (PANSTARRS)

References

All Tables

All Figures

Fig. 1 Comparison of nominal minimal distance from the Sun and spread of star clones at their closest approach for seven stars mentioned in this paper. Black dots in the center of each clone cloud mark the nominal star position at the closest Sun–star approach; see text on how this composite figure was obtained. All 10 000 clones of star C are hidden under its black dot because the stellar data are extremely precise in this case. Label meanings: A: P0509, B: P0230 (HD 7977), C: P0107 (Gliese 710), D: P0506, E: P0508, F: P0417 (Ton 214), G: P0111 (HIP 94512); see Table 2 for details of each star. The blue circle marks the habitually adopted outer limit of the Oort cloud.

In the text

Fig. 2 Past and future dynamical evolution of C/2014 UN271 nominal orbit (“b8” solution). changes in a perihelion distance (green), an inclination (fuchsia), and an argument of perihelion (red) are shown. the thick lines depict the result of the full dynamical model while thin lines show the evolution of elements in the absence of any stellar perturbations, i.e., only the galactic perturbations are taken into account. angular elements are expressed in a galactic frame. we also show names of stars that have a significant effect on this dynamical evolution.

In the text

Fig. 3 Influence of BB orbit uncertainty on parameters of the previous orbit. black dots show the results of the backward numerical integration of 5001 clones of BB stopped at the previous perihelion with all Galactic and stellar perturbations included. The center of the green circle marks the nominal result. The red histogram presents a marginal distribution of qprev along with the best Gaussian approximation (green dots in the middle of each bar). The blue histogram is a marginal distribution of 1/aprev.

In the text

Fig. 4 Influence of P0230 data uncertainty on parameters of the previous BB orbit – in this simulation we used 10 000 pairs: a comet clone and a P0230 clone. The result of each pair is represented by a black dot. Omitted are 114 points from a long right tail, among which 27 have negative 1/aprev. The marginal distribution of qprev is performed for all black points and shown as the red histogram below the main plot. For the sake of comparison, the main result is superimposed, with the blue swarm of dots showing the result of an additional experiment: 5001 BB orbits and nominal P0230 track. The nominal result is in the center of a green circle.

In the text

Fig. 5 Composite picture of the influence of cometary and stellar uncertainties on the BB orbit recorded at the next perihelion passage. Blue dots mark 5001 clones of BB, each perturbed by all stars on their nominal tracks. Orange dots represent the nominal BB orbit perturbed by 10 000 clones of P0107. Black dots present 10 000 pairs of the interacting BB and P0107 clones. The red marginal distribution histogram of qnext corresponds to the black cloud of points. The result for nominal stellar and cometary data is in the center of the green circle.

In the text

	Fig. 6 Dynamical past evolution of the nominal orbit of C/2017 K2 based on the “a9” orbital solution. Several individual stellar perturbations are marked. The meanings of the colors and line thicknesses are the same as in Fig. 2.
In the text

	Fig. 7 Influence of K2 orbit uncertainty (for solution “a9”) on the parameters of its previous orbit. The number of points, color, and symbol meanings are the same as in Fig. 3.
In the text

Fig. 8 Influence of P0509 data uncertainty on the parameters of the previous K2 nominal orbit (solution “a9”). Black dots mark 8827 cases where the distance between K2 and P0509 is greater than 20 000 au while the 1149 orange dots show the opposite cases. Two black and 23 orange points are omitted as the extreme outliers: the maximum qprev in this simulation is equal to 9510 au. The red marginal qprev distribution was made for all 9976 plotted points.

In the text

Fig. 9 Composite picture of the influence of cometary and stellar uncertainties on the K2 orbit recorded at the previous perihelion passage. Blue dots mark 5001 K2 clones perturbed by all stars on their nominal tracks. Orange dots represent the nominal K2 orbit perturbed by 10 000 P0230 clones. Black dots present 10 000 pairs of K2 and P0230 clones. The red marginal distribution histogram of qprev corresponds to the black cloud of points. The result for nominal stellar and cometary data is in the center of the green circle. Omitted are 226 outlier points from the far right tail of the black swarm and 219 from the orange one.

In the text