© ESO 2022

1 Introduction

Long-period comets (LPCs) travel on near-parabolic orbits with huge heliocentric distances ahead and behind them. When they are far from the Sun, they are gravitationally influenced by passing stars and by the tidal action of the overall Galactic potential. While the importance of stellar perturbations in their motion were recognized long ago (Öpik 1932; Oort 1950), the differential influence of the Galaxy as a whole started to be investigated at the end of the 20th century. Among first papers presenting the practical way to include Galactic perturbations in LPCs dynamics we mention that by Heisler & Tremaine (1986). They derived rather simple formulae accounting for the Galactic disk tidal action, which they recognized to be dominating. The problem of LPC motion under the Galactic disk influence can be approximated as an integrable one if the mean anomaly is removed by a perturbation method. Such an averaged solution was studied by a number of authors seeking its qualitative properties and approximate long-term solution. For more details, see, for example, Matese & Whitman (1992), Breiter et al. (1996), or a review paper by Fouchard et al. (2005) and references therein.

Adding the Galactic central bulge influence, which was first proposed by Levison et al. (2001), leads to a qualitatively different problem that remains non-integrable even after the removal of the short-period terms. Thus, even the secular dynamics of LPCs requires numerical integration and only basic facts about the equilibria of this model have been reported (Breiter et al. 2008).

In a series of papers, Królikowska and Dybczyński (e.g., Królikowska & Dybczyński 2020 and references therein), the authors investigate past and future motion of a large sample of the observed LPCs, taking into account both the Galactic and stellar perturbations based upon the available stellar data and the Galaxy potential models. In the Hipparcos era, there was no star identified that would significantly perturbing LPC motion and, additionally, a simplified Galactic potential that only accounted for the tidal effect of the Galactic disk was widely used.

Our knowledge of the Galactic and stellar neighborhood of our planetary system has expanded significantly in recent years. Several models of the Galactic gravitational potential have been proposed; for instance, Bajkova & Bobylev (2017) present a review. The number of recognized stars potentially perturbing LPC motion has also significantly increased, particularly thanks to the Gaia mission results (Gaia Collaboration 2016a). Based on the Gaia DR2 catalog (Gaia Collaboration 2018), the StePPeD database of potential perturbers was presented by Wysoczańska et al. (2020b). It consists of about 650 stars or stellar systems. Following the publication of the Gaia EDR3 data, StePPeD database is undergoing a substantial revision and the results are available in part.

Using the modern Galactic potential model proposed by Irrgang et al. (2013), it was shown by Dybczyński & Berski (2015) how a star trajectory could be calculated by means of the direct numerical integration in the galactocentric rectangular frame. This method was recently used to find the first stars that significantly perturb the LPC motion (Wysoczańska et al. 2020a). Over the course of that work, it appeared that calculating LPC motion together with stars and under the overall Galactic potential (as the N-body problem with N of the order of 400) encounters two serious problems, making the task very difficult.

The first problem comes from the fact that in the galactocentric frame, the positions of a comet and the Sun can be almost identical up to the first ten significant digits. This loss in precision leads to erroneous results in many cases, with regard to the long-term numerical integrations of the LPC motion. The only remedy is to switch to a quadruple precision, but this drastically worsens the second problem, namely, unacceptably time-consuming calculations. The N-body code with hundreds of interacting bodies requires hours to complete calculations for a single LPC, whereas we want to estimate the uncertainty of our results by means of the Monte Carlo sampling, that is, replacing a comet or a star with thousands of its clones drawn from the respective covariance matrices. We found this task nearly impossible, and this was the main impulse to develop alternative methods and algorithms described in this manuscript.

In Sect. 2, we describe a new way of formulating the equations of motion of a small body in the Solar System in a heliocentric frame, while maintaining the Galactic tides force intact. Section 3 shows the method for including stellar perturbations in this dynamical model. The next two sections describe the overall picture of the Galactic and stellar perturbations on LPC motion obtained with the proposed methods and present a detailed example of a comet dynamics investigation available with our approach. Section 6 shows how effectively we can estimate the uncertainty of the above-mentioned results and describes the difficulties arising from continuously imprecise stellar data in some cases. Our summary and conclusions are given in Sect. 7.

2 Exact differential formulae for Galactic tides in the heliocentric frame

In this paper, we use the term “galactocentric frame” for the reference frame whose origin is placed in the center of the Galaxy, whereas the axes are parallel to those of the standard Galactic heliocentric frame (withthe same directions). This means that in the such defined galactocentric frame, the OX axis is directed opposite to the Sun, the OY is directed in accordance to the Sun’s rotation around the Galactic center, and the OZ axis is directed to the north Galactic pole, which makes this system a right-handed one. For the details and numerical parameters of the frame orientation and translation, we refer to Dybczyński & Berski (2015).

We describe the galactocentric coordinates of the Sun as: $$R_{\odot }={\left(X_{\odot },Y_{\odot },Z_{\odot }\right)}^{\text{T}},$$(1)

and the galactocentric coordinates of a comet (or another point mass, e.g., a star) as $$R={\left(X,Y,Z\right)}^{\text{T}}.$$(2)

If we denote the heliocentric coordinates of a comet as: $$r={\left(x,y,z\right)}^{\text{T}},$$(3)

we have: $$R=R_{\odot }+r.$$(4)

The gravitational potential of the Galaxy Φ(ρ, Z) considered here is based on Model I from Irrgang et al. (2013), with parameters listed in Table 1, see also Dybczyński & Berski (2015). It is a revised version of the Allen & Santillan (1991) potential and is expressed as the sum of three components: a central bulge component Φ_b (R) (spherically symmetric), an axisymmetric disk Φ_d(ρ, Z), and a massive spherical Galactic halo Φ_h(R) (dark matter included). Thus: $$\Phi (\rho ,Z)={\Phi }_{\text{b}}(R)+{\Phi }_{\text{d}}(\rho ,Z)+{\Phi }_{\text{h}}(R),$$(5)

where (ρ, φ, Z) are the galactocentric cylindrical coordinates, $\rho =\sqrt{X^2+Y^2}$, and $R=\sqrt{{\rho }^2+Z^2}$ is a galactocentric spherical radius introduced by Eq. (2). In a very recent paper by Bovy (2020) this model of a Galactic potential is the only one which is shown to be in agreement with the Solar System’s galactocentric acceleration deduced from the Gaia EDR3 data (Gaia Collaboration 2021b). It should be noted that the formulae derived below can be applied to other particular Galactic potential models, with different numerical constants replacing those from Table 1.

For the Galactic bulge and disk components, the formulae are given in the form proposed by Miyamoto & Nagai (1975): $${\Phi }_{\text{b}}(R)=\frac{-G{M}_{\text{b}}}{\sqrt{R^2+b_{\text{b}}^2}}={\Phi }_b(X,Y,Z)=\frac{-G{M}_{\text{b}}}{\sqrt{X^2+Y^2+Z^2+b_{\text{b}}^2}},$$(6) $$\begin{array}{ll}{\Phi }_{\text{d}}(\rho ,Z)\hfill & =\frac{-G{M}_{\text{d}}}{\sqrt{{\rho }^2+{\left(a_{\text{d}}+\sqrt{Z^2+b_{\text{d}}^2}\right)}^2}}\hfill \\ \hfill & =\frac{-G{M}_{\text{d}}}{\sqrt{X^2+Y^2+{\left(a_{\text{d}}+\sqrt{Z^2+b_{\text{d}}^2}\right)}^2}}.\hfill \end{array}$$(7)

For the Galactic halo, Irrgang et al. (2013) proposed the following form: $${\Phi }_{\text{h}}(R)=\{\begin{array}{ll}\frac{G{M}_{\text{h}}}{a_{\text{h}}}\left(\frac{1}{\gamma -1}\ln \left(\frac{1+{\left(\frac{R}{a_h}\right)}^{(\gamma -1)}}{1+{\left(\frac{\Lambda }{a_{\text{h}}}\right)}^{(\gamma -1)}}\right)-\frac{{\left(\frac{\Lambda }{a_{\text{h}}}\right)}^{(\gamma -1)}}{1+{\left(\frac{\Lambda }{a_{\text{h}}}\right)}^{(\gamma -1)}}\right),\hfill & \text{ifR}<\Lambda ,\hfill \\ -\frac{G{M}_{\text{h}}}{R}\frac{{\left(\frac{\Lambda }{a_{\text{h}}}\right)}^{(\gamma -1)}}{1+{\left(\frac{\Lambda }{a_{\text{h}}}\right)}^{(\gamma -1)}},\hfill & \text{elsewhere},\hfill \end{array}$$(8)

which, after adopting γ = 2, and choosing the first equation in (8) because we do not allow a comet to go as far as R = Λ = 200 kpc, reduces to: $$\begin{array}{ll}{\Phi }_{\text{h}}(X,Y,Z)\hfill & =\frac{G{M}_{\text{h}}}{a_{\text{h}}}\left(\ln \left(\frac{a_{\text{h}}+\sqrt{X^2+Y^2+Z^2}}{a_{\text{h}}+\Lambda }\right)-\frac{\Lambda }{a_{\text{h}}+\Lambda }\right).\hfill \end{array}$$(9)

In the absence of other bodies, the equations of a comet and the Sun motion in the Galactic frame can be written as: $$\begin{array}{ll}\stackrel{\mathrm{..}}{R}\hfill & =-\nabla \Phi (R)-\frac{G{M}_{\odot }}{r^3}r\hfill \\ \hfill & =F_{\text{b}}(R)+F_{\text{d}}(R)+F_{\text{h}}(R)-\frac{G{M}_{\odot }}{r^3}r,\hfill \end{array}$$(10) $$\begin{array}{ll}{\stackrel{\mathrm{..}}{R}}_{\odot }\hfill & =-\nabla \Phi (R_{\odot })=F_{\text{b}}(R_{\odot })+F_{\text{d}}(R_{\odot })+F_{\text{h}}(R_{\odot }),\hfill \end{array}$$(11)

where Φ is the overall potential of the Galaxy (5), and the subscripts of F, similarly to those of Φ, refer to their division into three commonly used parts: central bulge, disk, and halo.

The heliocentric motion of a comet results from $$\stackrel{\mathrm{..}}{r}=\stackrel{\mathrm{..}}{R}-{\stackrel{\mathrm{..}}{R}}_{\odot }=-\frac{G{M}_{\odot }}{r^3}r+f(R_{\odot },r),$$(12)

where $$\begin{array}{ll}f(R_{\odot },r)=\hfill & \left[F_{\text{b}}(R_{\odot }+r)-F_{\text{b}}(R_{\odot })\right]\hfill \\ \hfill & +\left[F_{\text{d}}(R_{\odot }+r)-F_{\text{d}}(R_{\odot })\right]\hfill \\ \hfill & +\left[F_{\text{h}}(R_{\odot }+r)-F_{\text{h}}(R_{\odot })\right]=f_{\text{b}}+f_{\text{d}}+f_{\text{h}}.\hfill \end{array}$$(13)

To avoid the subtraction of almost equal terms and the resulting loss of significant digits, the expressions for all parts of f need to be appropriately transformed.

Let us introduce the following quantities: $$\begin{array}{ll}{p}_1\hfill & =\sqrt{R_{\odot }^2+b_{\text{b}}^2},\hfill \end{array}$$(14) $$\begin{array}{ll}{q}_1\hfill & =\sqrt{R^2+b_{\text{b}}^2},\hfill \end{array}$$(15) $$\begin{array}{ll}{p}_2\hfill & =\sqrt{X_{\odot }^2+Y_{\odot }^2+{(a_{\text{d}}+p_3)}^2},\hfill \end{array}$$(16) $$\begin{array}{ll}{q}_2\hfill & =\sqrt{X^2+Y^2+{(a_{\text{d}}+q_3)}^2},\hfill \end{array}$$(17) $$\begin{array}{ll}{p}_3\hfill & =\sqrt{Z_{\odot }^2+b_{\text{d}}^2},\hfill \end{array}$$(18) $$\begin{array}{ll}{q}_3\hfill & =\sqrt{Z^2+b_{\text{d}}^2},\hfill \end{array}$$(19) $$\begin{array}{ll}\sigma \hfill & =r^2+2r\cdot R_{\odot }.\hfill \end{array}$$(20)

Then, the components of a comet’s heliocentric acceleration per unit mass due to the Galactic tide can be expressed as: $$f_{\text{b}}=\frac{G{M}_{\text{b}}}{q_1^3}\left[\frac{\sigma \left(p_1^4+q_1^4+p_1^2{q}_1^2\right)}{p_1^3\left(p_1^3+q_1^3\right)}{R}_{\odot }-r\right],$$(21) $$f_{\text{h}}=\frac{G{M}_{\text{h}}}{a_{\text{h}}R\left(a_{\text{h}}+R\right)}\left[\left(1+\frac{a_{\text{h}}}{R+R_{\odot }}\right)\frac{\sigma }{R_{\odot }\left(a_{\text{h}}+R_{\odot }\right)}{R}_{\odot }-r\right],$$(22) $$f_{\text{d}}=\frac{G{M}_{\text{d}}}{q_2^3}\left[\frac{p_2^4+q_2^4+p_2^2{q}_2^2}{p_2^3\left(p_2^3+q_2^3\right)}\left(\sigma +2a_{\text{d}}\frac{z\left(z+2Z_{\odot }\right)}{p_3+q_3}\right)R_{\odot }-r+P\right],$$(23)

where $$P=\frac{a_{\text{d}}}{q_3}\left(\begin{array}{c}0\\ 0\\ \frac{A}{p_2^3{p}_3\left(p_2^3{p}_3+q_2^3{q}_3\right)}{Z}_{\odot }-z\end{array}\right),$$(24)

and $$A=p_3^2\left(\sigma +2a_d\frac{z\left(z+2Z_{\odot }\right)}{p_3+q_3}\right)\left(p_2^4+q_2^4+p_2^2{q}_2^2\right)+q_2^{6}z\left(z+2Z_{\odot }\right).$$(25)

The strict differential formulae can also be derived for an alternative model of the halo component of the Galactic potential, Model III in Irrgang et al. (2013), which is deduced from the halo density profile of Navarro et al. (1997), hereafter NFW, that is, $${\Phi }_{\text{H}}(R)=-\frac{G{M}_{\text{H}}}{R}\ln \left(1+\frac{R}{a_{\text{H}}}\right).$$(26)

For the NFW halo potential model, from Eq. (26) we have $$F_{\text{H}}=-\nabla {\Phi }_{\text{H}}(R)=-\frac{G{M}_{\text{H}}}{R^3}\left[\ln \left(1+\frac{R}{a_{\text{H}}}\right)-\frac{R}{a_{\text{H}}+R}\right]R,$$(27)

and f_H =F_H(R) −F_H(R_⊙), should be transformed to $$f_{\text{H}}=\frac{G{M}_{\text{H}}}{R^3}\left[W_1{R}_{\odot }+W_{2}r\right],$$(28) $$W_2=\frac{R}{a_{\text{H}}+R}-\ln \left(1+\frac{R}{a_{\text{H}}}\right),$$(29) $$\begin{array}{l}{W}_1=\ln \left(1-\frac{\sigma }{\left(a_{\text{H}}+R\right)\left(R_{\odot }+R\right)}\right)-\frac{a_{\text{H}}\sigma R}{\left(a_{\text{H}}+R\right)\left(a_{\text{H}}+R_{\odot }\right)R_{\odot }^2}\\ +\frac{\sigma \left(R_{\odot }^2+R^2+R_{\odot }R\right)}{\left(R_{\odot }+R\right)R_{\odot }^3}\left[\ln \left(1+\frac{R_{\odot }}{a_{\text{H}}}\right)-\frac{R{R}_{\odot }}{\left(a_{\text{H}}+R\right)\left(a_{\text{H}}+R_{\odot }\right)}\right].\end{array}$$(30)

If it does not degrade the accuracy (e.g. in quadruple arithmetic), we can use a simpler form of the first logarithm in W₁ : $$\ln \left(1-\frac{\sigma }{\left(a_{\text{H}}+R\right)\left(R_{\odot }+R\right)}\right)=\ln \left(\frac{a_{\text{H}}+R_{\odot }}{a_{\text{H}}+R}\right).$$(31)

The derivation of the above equations is elementary. The outline can be found in Appendix A.

If one aim is to study only the Galactic tide influence on Solar System bodies, the above proposed formulae offer a fast and accurate solution. However, they require galactocentric positions of the Sun R_⊙ (t). It is possible to construct more or less approximate analytical formulae to achieve this, for example, by obtaining a polynomial approximation of the results of direct numerical integration of the Sun’s galactocentric motion. Our tests showed that such ephemeris can be quite compact. For example, the backward motion of the Sun in the Galaxy spanning 30 Myr can be successfully described by six subsequent intervals covered by approximating polynomials of degree 5, individually for each coordinate. In this case. the whole backward ephemeris of the Sun takes only 120 double precision numbers. We offer an example of such an ephemeris in Appendix B.

Since we are interested in both Galactic and stellar perturbations, we decided to prepare a compound source for the Sun and star positions, described in Sect. 3.1.

3 Simultaneous Galactic and stellar perturbations

When searching for the source region or regions of the observed LPCs, we have to study their past dynamics, starting from their “original” orbit and taking into account both Galactic and stellar perturbations. There is an ongoing project lead by Królikowska and Dybczyński, aimed at revealing the dynamical history of the observed LPCs; for more, see Królikowska & Dybczyński (2020) and references therein.

Adding stellar perturbations to the heliocentric comet equations of motion is much simpler than carrying out this process for Galactic tides. It is reduced to including a sum of point-mass interactions with all stars described by their heliocentric positions. This requires an extension of (12) into $$\stackrel{\mathrm{..}}{r}=-\frac{G{M}_{\odot }}{r^3}r+f(R_{\odot },r)+{\sum }_{s=1}^N\left[\frac{G{M}_s(r_s-r)}{|r_s-r{|}^3}\right],$$(32)

where r_s =R_s −R_⊙, and the summation in the last term is over all N stars included in the model. The use of this equation requires the knowledge of the positions of the Sun and of all the stars in the galactocentric frame at each moment covered by the numerical integration. Obtaining these Sun and star positions requires a numerical integration of their galactocentric motion that takes into account their mutual interactions and an overall Galactic potential. We note that this should not be repeated each time we investigate a subsequent comet. For this purpose, we decided to build a numerical ephemeris of the Sun and stars.

3.1 Numerical ephemeris of the Sun and the stars galactocentric motion

For the numerical integration, we use the well-known and widely recommended Gauss-Radau integrator published by Everhart (1985). Although we have tested various orders of the integrator, the practice shows that the classical 15th-order RA15 works best for the numerical integration of a comet motion in the heliocentric frame. The integrator implements a collocation method, which means that at every step (referred to as a “sequence” by Everhart), a polynomial of time is constructed to approximate the integrated body motion. Although the polynomial is evaluated only at the end of each sequence, its coefficients may also serve to create a dense output: similarly to the Taylor series methods, the collocation methods are capable of providing the position and velocity at any time with no additional cost (even though the accuracy between the endpoints is formally lower than in the Taylor polynomial case). We decided to use this property of the Gauss-Radau integrator to construct the numerical ephemeris by simply recording all polynomial coefficients at the end of each sequence when integrating the Sun and stars in the galactocentric frame.

The list of potential stellar perturbers used to construct an ephemeris is based on the publicly available stellar database¹ (version 2.3) announced by Wysoczańska et al. (2020b). To make our proposition more up to date, we refined this list using the latest Gaia EDR3 catalog (Gaia Collaboration 2021a) and making this partial update publicly available as the StePPeD release 3.0. It is worth mentioning that thanks to the more accurate and complete data in the Gaia EDR3 catalog, we identified 43 additional double stars. As a result, the number of double or multiple systems included in our list of potential perturbers increased to 81. In all available cases, the center of mass of the system is used. The work on the complete update of the whole StePPeD database by incorporating Gaia EDR3 results which also include some new stars is in progress. As a result of this refinement, the whole number of stars or stellar systems taken into account is 407, but in order to optimize the ephemeris file size and the speed of its use, we decided to limit the number of perturbers by the following rules.

First of all, we produce two separate ephemerides: one for integrating a comet motion backward and one for the forward integration. This is fully justified because studying cometary past dynamics (starting from the “original” orbit) and its future dynamics (starting from the “future” orbit) are always separate tasks.

Moreover, most of the perturbers are now much farther than the 4 pc limit assumed when constructing the database. As a result, in the ephemeris dedicated for forward integrations, we include 25 stars or stellar system that are currently closer than 4 pc, and 175 perturbers that will pass the Sun closer than 4 pc in the future, which makes the total number of perturbers in this ephemeris equal 200 plus the Sun. In the same manner, but looking to the past, we restricted the number of perturbers in the backward ephemeris to 25 + 207 = 232 objects plus the Sun.

For the purpose of storing ephemeris data in binary files, we used integer numbers (4 bytes), double precision floating point numbers (8 bytes) and 8 byte character strings for perturber names. Both ephemerides have the same format of the binary file. It starts with the header: two integer numbers and a table of perturber names and masses, structured as described in Fig.1.

Next, the main body of the ephemeris follows. It consists of a large number of “one-step data blocks”. Each of them contains two epochs, expressed as Julian Days and stored as double precision floating point numbers. These are the ends of the validity time interval of polynomial coefficients stored in the rest of one-step data block, see Fig. 2.

Polynomial coefficients are stored in body data blocks, each containing 10 double precision floating point numbers for each of the three coordinates of the body in the galactocentric frame, expressed in parsecs. As a result, each body data block is 240 bytes long. In the backward ephemeris there are 233 body data blocks in each one-step data block. In the forward ephemeris the number of bodies is only 201.

We decided for both ephemerides to cover 30 Myr, which is comparable with the orbital period of a comet having the semi-major axis of 100 000 au, but only for a handful of observed LPCs we have to perform a numerical integration spanning more than 10 Myr. As a result we have 2323 one-step data blocks in the backward ephemeris, and 2382 in the forward one. The total size of the ephemerides files is 128 MB and 115 MB, respectively. The averagestep size is on the level of 10 000 yr but it goes down as low as to ten years during the closest encounters.

It should be noted here that these ephemerides were obtained from the precise numerical integration (80 bits arithmetic, i.e., long double in the C language) of the stars and the Sun in the galactocentric frame under the overall Galactic potential and taking into account all mutual gravitational interactions. The internal precision was kept on the level of 10⁻¹¹ in terms of thegalactocentric positions in parsecs, which gives 15 significant digits. Starting data for all bodies are based on nominal stellar data taken in most cases from the Gaia EDR3 catalog and are referred to the epoch 2016.0. For such details as the Sun initial conditions, numerical parameters, frame definition, or dynamical model, we refer to Dybczyński & Berski (2015).

The stellar ephemeris proposed here is designed for screening tests when searching for detectable star-comet interactions, using long lists of comets and potential stellar perturbers. It can also be used when analyzing the influence of the cometary data uncertainties. However, finally, when investigating a certain case of the close encounter, it is necessary to also take into accountthe stellar data uncertainty that cannot be performed with this ephemeris, based on nominal stellar data.

Both ephemerides binary files as well as the documentation and examples of a code using them are publicly available on the Internet² as a part of our stellar database site (see the download section). The ephemerides (and the code if necessary) will be updated each time we update our database of potential stellar perturbers.

Fig. 1 Structure of the header part of the binary ephemeris file; “Serial” is the (4 byte integer) code of the ephemeris version (1255 at present), “number” is the (4 byte integer) number of bodies included. They are followed with the table of “number” rows, containing 8 byte character string with the perturber name and double precision floating point number containing the mass of the perturber.

Fig. 2 Structure of the “one-step data blocks” of the binary ephemeris file; t1 and t2 are the ends of the validity time interval of polynomial coefficients stored in “body data blocks”. The number of these blocks equals the number of bodies included in the ephemeris. Each “body data block” is 240 bytes long (3 × 10 double precision floating point numbers).

3.2 Computer implementation for cometary dynamics

As described in Królikowska & Dybczyński (2020) for each studied LPCs we would like to follow their motion to the past, starting from their “original” orbits, and to the future, starting from “future” orbits. These integrations are typically stopped at the previous or next perihelion, or when a comet leaves a heliocentric sphere of the 120 000 au radius. The previous perihelion distance value can be used to discriminate between dynamically old and new comets.

While preparing the paper cited above, we encountered problems with these calculations resulting from losing significant digits due to a very small difference between the Sun and comet positions expressed in the galactocentric frame. This enforced us to perform a high accuracy numerical integration which appeared to be extremely time consuming. As a result, we decided to switch off all mutual interactions between stellar perturbers in the dynamical model. This greatly speeded up our calculations and in the great majority of cases the resulting comet “previous” and “future” orbits are correct. But there are several cases where the mutual inter-stellar gravitational action changes a particular perturber trajectory in a manner, which changes the resulting “previous” or “next” orbits slightly. Those problems stimulated us to develop the alternative methods described in this paper.

During the ephemeris production, when integrating backwards, we recorded nine encounters between stars closer than 0.2 pc, and only two of them are the encounters with the Sun. Three of the encounters were closer than 0.1 pc: stars Gaia EDR3 2946037094762449664 and Gaia EDR3 5261593808165974784 approached each other at a distance of 0.05 pc 0.89 Myr ago, the stars TYC 160-136-1 and Gaia EDR3 5949463366461232896 miss distance was 0.09 pc 1.7 Myr ago, and the star HIP 31626 passed near the double star p Eridani at a distance of 0.04 pc 0.42 Myr ago. It is worth noting that the same star HIP 31626 met another one, Gaia EDR3 2730508416002618752 a little bit later. When analyzing the future motion we recorded seven mutual star encounters below 0.2 pc. In the rare cases of such close encounters, even small modifications of stellar trajectories (on the level of 0.1 pc) can induce substantial changes in cometary motion, as shown in Wysoczańska et al. (2020a). It should be stressed that the above data for the closest encounters are obtained from the nominal stellar data, which have a different level of accuracy.

With the exact differential formulae described in Sect. 2, and having numerical ephemerides containing the Sun and thepotential stellar perturbers’ positional galactocentric data, the construction of a computer program for studying the comet’s long-term motion is not difficult. Preparing the code for Galactic tides is quite straightforward, save for one exception: if we uses the NFW variant of the Galactic halo model, it is necessary to use the special numerical function for calculating ln (1 + ξ) when the value of ξ is small (e.g., function log1p in the C/C++ language mathematical library) to evaluate Eqs. (30) and (29).

It is also worth observing that the step size of the stellar ephemeris is typically several orders larger than the step size in the comet heliocentric motion integration. For the optimal use of the ephemeris, the calling function should use static variables to monitor whether a new set of coefficients should be read from a file.

In our code, we use internally dedicated units for keeping the optimal numerical precision, namely: a solar mass, 10³ yr, and 10⁴ au. When reading the ephemeris file, the rescaling is necessary as well as the transformation of stellar positions from galactocentric to heliocentric, with the Sun position included.

In all the above, we use the term “heliocentric” for simplicity but, in fact, we work within the frame, with the origin at the barycenter of the Solar System. It is a widely accepted practice that the “original” and “future” comet orbits are presented in a barycentric frame (see Królikowska & Dybczyński 2020 for details) and we keep our calculations of the “previous” and “next” orbits in the same reference frame. As a result, we used a central body mass equal to 1.001341838222754 (value taken from the DE430 JPL planetary ephemeris constants, see Folkner et al. 2014).

The comparison with our previous, direct galactocentric integrations of a comet motion showed the improvement of an internal precision by at least six orders of magnitude. Moreover, the present code works over 100 times faster. It is worth mentioning that obtaining in the galactocentric integration the accuracy comparable to the present code requires quadruple precision, which is approximately 200 times slower.

4 Overall stellar impact on the cometary orbit evolution

Using the proposed methods we performed an overall test of the stellar perturbations effect on the past motion of LPCs. For all 270 preferred original comet orbits taken from the version 3.0 of the CODE catalog (Królikowska & Dybczyński 2020), we performed a series of numerical integrations. Each comet was first integrated under the influence of the Galactic tide and all 232 stellar perturbers taken from the release 3.0 of the StePPeD database. Then we repeated this calculation 232 times with each individual stellar perturber removed. We compared the comet previous perihelion distance obtained from the full model and the one with a particular star omitted. In this way, we reasonably measured the influence of each individual perturber on each comet’s past motion.

In Table 2, we list five stellar perturbers which showed the strongest influence on the analyzed sample of LPC orbits. The first three columns contain the star identifications. In the fourth column, we present the number of comets for which we observed a change in the previous perihelion distance greater than 10% in the numerical experiment described above. The same, but for changes greater than 50%, is showed in the fifth column. Next three columns show the parameters of each particular star closest passage near the Sun, and the estimate of its mass used in our calculations. In the last column we show the relative velocity at the moment of the star closest approach.

Solar proximity parameters together with their uncertainties were obtained from a series of the extensive numerical Monte Carlo simulations. In the precise numerical integration performed in the galactocentric frame, we use all 232 perturbers and take into account all their mutual interactions. For the uncertainty estimation we repeated this calculation replacing a particular star with at least of 10 000 of its clones (using nominal data for the remaining objects). These clones were drawn according the six-dimensional covariance matrix from Gaia EDR3³. For each perturber we generated more than 10000 of its clones, numerically integrated their past motion, and stopped them at their closest approach to the Sun. This method of calculating the uncertainty estimation is obviously simplified since we take into account only one star’s uncertainties, while the rest move along their nominal tracks. Fortunately, the close star-star encounters are very rare (and frequently with a high relative velocity), so this simplification still offers the correct order of the scatter of the results. As it is described in detail in Dybczyński & Berski (2015), we calculated the coordinates of the clone cloud centroid (calculating a simple mean in each heliocentric coordinate, i.e., $\overline{x}$, ȳ, $\overline{z}$), and then obtaining its distance from the Sun: $$\text{mindist}=\sqrt{{\overline{x}}^2+{\overline{y}}^2+{\overline{y}}^2}.$$

As the uncertainty of this expected minimal distance, we measure the radius of a sphere around the centroid ($\overline{x}$,ȳ,$\overline{z}$) containing 90% of individual clone proximity points. The epochs of the closest approach are presented in Table 2 as their mean and standard deviation since the distribution of these results is close to Gaussian in all five cases.

If the presented uncertainty in the minimal distance is greater than the value of the expected proximity distance (as in the two cases included in Table 2), we should interpret this fact by expecting some fraction of clones to pass on the other side of the Sun than the centroid. In other words, the position of the Sun is inside the sphere of 90% of clones, and as a result, an arbitrarily small proximity distance is statistically possible. We show an example of such a situation in a case of the first and the most important star from Table 2, namely, P0230 (HD 7977); see Fig. 3.

The possibility of a close approach of this star to the Sun was first pointed out by Bailer-Jones et al. (2018). In earlier papers, this star was omitted because its parallax was first measured by Gaia mission and published in the Gaia DR1 catalog Gaia Collaboration (2016b) and its radial velocity appeared only in Gaia DR2 (Gaia Collaboration 2018). In Bailer-Jones et al. (2018), based on the Gaia DR2 data, the authors estimated the nominal closest distance of P0230 to the Sun as 0.43 pc.

That trajectory of P0230 was recognized by Wysoczańska et al. (2020a) as passing very close to the long period comet C/2002 A3 LINEAR, possibly imparting a strong gravitational influence on the comet motion. Following updates to the Gaia EDR3 data (Gaia Collaboration 2021a), the situation has changed substantially. Now the nominal minimal Sun-star distance is only 0.03 pc (see the first row in Table 2), and based on Fig. 3, we clearly see that due to the still significant uncertainties in the stellar data for P0230, this distance might happen to be arbitrarily small or significantly larger, and the geometry of this passage is highly uncertain.

According to the data in Table 2, 145 of LPCs listed in the CODE (Królikowska & Dybczyński 2020) database might have had their previous perihelion distance changed by more than 10% due to the gravitational action of this star, and 134 of them by more than 50%. Taking into account that only 129 comets from our list have their previous orbital period longer than 2.8 Myr, it is evident that all comets that exhibited a chance of an encounter with P0230 via our calculations were affected.

The importance of this perturber comes from its relatively large mass, estimated to be equal to 1.08 M_⊙ (Stassun et al. 2019), but mainly from its very close flyby near the Sun. This is a well-known effect that the change of a comet orbit is the result of the difference between a gravitational pull induced by a passing star on a comet, and that gained by the Sun. In the case of such a close P0230 flyby, effectively, all comets can be indirectly affected. This important stellar passage near the Sun and its possible consequences will be studied in detail in a separate paper which is in preparation.

In the next section, we present a detailed analysis of an individual example: the past dynamical evolution of the comet C/2015 XY1 Lemmon. It properly illustrates the cumulative influence of two different stars listed in Table 2.

Fig. 3 Distribution of the cloud of 16 545 clones of P0230 stopped at the closest proximity to the Sun about 2.80 Myr ago. All points are projected onto the maximum scatter plane X′ Y′. The green dot marks the nominal approach position, while the red one depicts the position of the Sun. The blue circle shows the widely adopted approximate boundary of the Oort cometary cloud at 0.5 pc from the Sun. Orange dot and arrows showthe position of C/2015 XY1 at the moment of the nominal P0230 approach; see Sect. 5 for details.

5 Example result in detail: past dynamics of comet C/2015 XY1 Lemmon

When reviewing the effect of stellar perturbations on all LPCs included in the CODE database, we found a very interesting case of the past motion of C/2015 XY1 Lemmon. This comet has an orbit of the greatest quality 1a+ (see Królikowska & Dybczyński 2018 for an explanation). As the starting point we used the preferred original orbit of this comet taken directly from the CODE database (Królikowska & Dybczyński 2020) For convenience, we quote the elements of this orbit in Table 3. We found that during the orbital period preceding the observed perihelion passage, this comet significantly changed its orbit twice, as a result of perturbations caused by stars P0506 and P0230. Based on nominal data, it was 1.09 Myr ago that P0506 passed 0.005 pc (~1100 au) from this comet. However, much earlier, that is, 2.8 Myr ago, when the comet was near its aphelion, it received another change of its heliocentric orbit as a result of the strong interaction between the Sun and P0230.

If we exclude these two events from our simulation, the past dynamical evolution of C/2015 XY1 would look like in Fig. 4. All perturbations from the other 230 stars were taken into account here, but they appeared to be small. In Figs 4–6, we keep the same scales for a simpler comparison. The left vertical axis is expressed in au and corresponds to the osculating perihelion distance plot (q, green line) as well as the heliocentric distance plot (r, thin, vertical blue line, visible only at the previous perihelion). The right vertical axis is expressed in degrees and describes the evolution of the osculating inclination (i, fuchsia line) and of the argument of perihelion (ω, red line). Both these angular elements are expressed in the Galactic reference frame. Thick lines show the dynamical evolution under the simultaneous stellar and Galactic actions, while thin lines depict the evolution when all stellar perturbations are excluded. In Fig. 4, it is clearly shown that after excluding the two most important perturbers, P0506 and P0230, the comet orbital elements evolve very similarly to the case when only Galactic tide is considered.

However, as it was mentioned above, our calculations show that 1.09 Myr ago this comet motion was strongly perturbed by the star P0506. The data in the Gaia EDR3 catalog for this star are much more precise than before, so its trajectory near the Sun is much better defined, as it is shown in Fig. 7. Additionally, a new and more precise radial velocity of this star was obtained by Errmann et al. (2020). In Fig. 7, we show 25 000 clones of this star (black dots) stopped at their closest approach to the Sun, and projected onto the X′ Y′ maximum scatter plane obtained by a principal component analysis (see Dybczyński & Berski 2015 for details). A green dotmarks the nominal position of the star at its closest approach, a red dot shows the position of the Sun. An inset at the right upper corner shows the cloud of the star clones enlarged ten times and the position of C/2015 XY1 (the orange cross) at the moment of the closest star-comet encounter, projected onto the same X′ Y′ plane. The uncertainty of a comet position is below 0.001 pc (~200 au). This uncertainty was obtained by the numerical integration of 5000 clones (hereafter virtual comets, VCs) of C/2015 XY1. The VCs were generated by the method proposed by Sitarski (1998), based on the osculation orbit of C/2015 XY1 and then propagated back in time to obtain a set of original orbits reflecting the observational uncertainties in comet positional observations. The set of these VCs can be downloaded from the CODE catalog web page⁴. Such an uncertainty estimation can be done quite simply and fast, using the stellar ephemeris described in Sect. 3.1.

The comet in Fig. 7 is slightly above the X′ Y′ maximum scatter plane and its nominal minimal distance from this star is only about 1120 au. Ninety percent of star clones are situated in the sphere with a radius of 0.013 pc (~2000 au) around the P0506 nominal proximity point. The cometary and stellar uncertainties taken into account simultaneously result in the minimal comet-star distance in the range of 13–5840 au obtained from a set of 12 000 random pairs of star and comet clones.

The nominal orbit evolution with perturbations from P0506 taken into account (i.e., now only P0230 is excluded) is shown in Fig. 5. One can notice a strong star-comet interaction 1.09 Myr ago. The angular elements do not change significantly,but an orbital period increases to almost 7 Myr, and the previous perihelion distance increases from 14 to ~230 au. This last effect might have been interpreted as a proof that this comet is dynamically new, but this is not true, since another significant orbit change occurred about 2.8 Myr ago due to P0230.

This nominal orbit evolution, compatible with our current knowledge, is depicted in Fig. 6. There we see the significant C/2015 XY1 orbit change 2.8 Myr ago. The orbital period is again a little longer, there is a big change in the inclination, and the perihelion distance is significantly cut down to ~50 au. However, the orbital evolution, depicted in Fig. 6, should be interpreted only as one of the possible scenarios due to the still great uncertainties on the P0230 passage geometry near the Sun.

Fig. 4 Dynamical past evolution of the nominal orbit of C/2015 XY1. Here P0506 and P0230 are excluded from the model. 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. Angular elements are expressed in the Galactic frame.

Fig. 5 Past evolution of C/2015 XY1 with only the star P0230 excluded. The strong interaction with P0506 is evident.

Fig. 6 Past changes in nominal orbital elements of C/2015 XY1 with all potential stellar perturbers included. The strong interaction of C/2015 XY1 with P0506 and P0230 is shown.

Fig. 7 Distribution of the cloud of 25 000 clones of P0506 stopped at the closest proximity to the Sun 1.09 Myr ago. All points are projected onto the maximum scatter plane X′ Y′. The green dot marks the nominal approach position, while the red one depicts the position of the Sun. The blue circle shows the widely adopted approximate boundary of the Oort cometary cloud at 0.5 pc from the Sun. An inset at the upper right shows the star clones cloud enlarged ten times to show the position of C/2015 XY1 (an orange cross) at the closest approach.

6 Influence of the stellar and cometary data uncertainties

Taking into account the relatively large uncertainties of stellar data for P0230, let us focus on the influence of the data uncertainty in the case of the close encounter between C/2015 XY1 and P0506. This encounter occurred 1.09 Myr ago, so we decided to analyze osculating elements of this comet orbit exactly at 2.0 Myr ago and to show the individual and combined effect of the uncertainties.

Before showing the distributions of orbital elements, we present the comparison between the influence of only cometary and only stellar data uncertainties on the minimum star-comet distance in the studied case. In Fig. 8, we show two distributions of the minimum star-comet distances. The blue one is based on 5000 VCs of C/2015 XY1 obtained as described in the previous section and the nominal stellar data. The green one shows the distribution of the same parameter, but based on random 12 000 star clones and the comet nominal data (before plotting, the histograms were normalized to the same vertical scale). The blue histogram extends from q = 792 to 1797 au with the maximum around the nominal value of 1100 au. The green histogram extends from q = 12 to 5573 au, and 146 clones from the right tail are out of the border of the plot. The maximum value for the green distribution is significantly larger and lays somewhere between 1200 and 1600 au. It is worth mentioning that the distribution of the minimum star-comet distance, when both a comet and a star were randomly drawn from their respective sources, is almost identical to the green one, and extends from q = 13 to 5837 au (with 152 clones above 4000 au). We do not present this distribution, since both histograms would overlap.

To study the spread of osculating elements at 2.0 Myr ago, we also started by analyzing the effect of thecometary data uncertainty. In Fig. 9, we present the plot of osculating 1∕a versus q for all 5000 VCs of C/2015 XY1 augmented with the marginal distribution of the perihelion distance, which is the most perturbed element during the P0506 passage (see Fig. 5). For each VC, we measured the minimum comet-star distance, and we use a black dot when this distance was smaller than 1000 au, and a red dot otherwise. In this plot, there are 1079 black dots and 3921 red dots. The blue cross marks the nominal result. It is clear from the q histogram that the most probable value is a little smaller than the nominal one.

However, to observe the overall uncertainty of C/2015 XY1 orbit at 2.0 Myr ago, we had to add the effect of the stellar data uncertainty. The result of a calculation performed for 12 000 pairs of a star clone and a comet clone acted by the former, is presented inFig. 10. Apart from black and red dots, we also use green dots in this plot, which corresponds to all pairs with the proximity distance greater than 2000 au. There are 1725 black, 5111 red, and 4324 green dots in this plot. For an easier comparison, we keep the same scales as in Fig. 9 and as a result, 841 black dots from the right tail of the distribution are not visible. These absent dots are extremely spread in both elements. For 1∕a, the covered range is from −1431.32 to 38.78, so some small fraction of hyperbolic orbits were also obtained. The largest perihelion distance obtained in this simulation is 46 160 au.

In Fig. 10, we observe that the most probable value of a perihelion distance is significantly smaller (~50 au) than the nominal q = 221 au (blue cross).This effect results from a large number of cases when a star-comet minimum distance is larger than the nominal (allgreen dots and most of the red dots). This is well illustrated in the inset of Fig. 7.

In light of the above results, it should be obvious that tracing the motion of C/2015 XY1 further to the past, and calculating the results of the strong perturbation caused by P0230 2.8 Myr ago, seems useless since almost all physically acceptable results will be possible after taking into account P0230 uncertainties. For example, according to our current knowledge, we cannot say anything about the previous perihelion distance of this comet – even a hyperbolic past orbit is possible due to the strong interaction with P0506.

Fig. 8 Comparison of the cometary and stellar data uncertainty influence on the minimum star-comet distance distribution.

Fig. 9 Effect of the cometary data uncertainty on C/2015 XY1 orbital elements recorded 2.0 Myr ago. See text for a detailed explanation.

Fig. 10 Simultaneous effect of both cometary and stellar data uncertainties. The scales are the same as in Fig. 9.

7 Summary and conclusions

Here, we propose a set of algorithms and methods for practical use in the Solar System small body dynamics investigation when Galactic and stellar perturbations have to be taken into account. The proposed approach allows for fast and precise numerical integration of the equations of motion formulated in the Solar System barycentric reference frame. The past and future motion of LPCs for millions years can be studied in detail and the reliability of the results can be effectively controlled using the Monte Carlo sampling of the input data.

The equations of motion derived in Sect. 2 were carefully tested by a comparison with the result of a quadruple precision calculation performed in the galactocentric frame. These equations give the full available precision. The numerical stellar ephemerides described in Sect. 3.1 are optimized for the standard double precision calculations. Their application increases the precision by 6 significant digits, while reducing the computation time by a factor 100. There is no trade-off between the precision and speed – both the quantities are improved. These calculations are over 20 000 times faster than their previous versions performed in the quadruple precision, which is necessary to obtain a comparable precision in the galactocentric frame. The numerical ephemerides of a set of potential stellar perturbers, as well as the example C code for their application and the approximate polynomial galactocentric Sun ephemerides are publicly available.

To present possible applications of the proposed methods, we performed the overall tests of the stellar perturbations effect on all LPCs included in the CODE database. We show which stars from our current list are the most effective comet motion perturbers. We also present a detailed analysis of one particular example – the past motion of comet C/2015 XY1. This comet orbit was strongly perturbed as a consequence of a very close passage of the star P0506 about 1.8 Myr ago, but it was also perturbed at an earlier time by P0230. We described the influence of both comet and stellar parameter uncertainties on the possibility to obtain this comet past dynamic state.

Acknowledgements

We are very grateful to the reviewer, Marc Fouchard, for his insightful comments that have helped us improve this manuscript. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research has also made use of the VizieR catalog access tool, CDS, Strasbourg, France (DOI: 10.26093/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 partially 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 Heliocentric accelerations

The transformation of equations of motion to a form that does not require subtraction of almost equal terms dates back to the Encke’s method of numerical integration, where the relative acceleration with respect to the nominal Keplerian orbit had to be computed (Brouwer & Clemence 1961). Here, the context is different and the potentials are more complicated, but the spirit and some of the details remain similar.

For example, transforming f_b, initially defined through Eq. (13), we proceed as follows. First, the coefficients of r and R_⊙ are separated: $$\frac{R}{{\left(R^2+b_b^2\right)}^{\frac{3}{2}}}-\frac{R_{\odot }}{{\left(R_{\odot }^2+b_b^2\right)}^{\frac{3}{2}}}=\frac{1}{q_1^3}\left(\frac{p_1^3-q_1^3}{p_1^3}{R}_{\odot }+r\right),$$(A.1)

where p₁ and q₁, given by Eqs. (14,15), are almost equal. Then, the problematic term is converted into $$\frac{p_1^3-q_1^3}{p_1^3}=\frac{p_1^6-q_1^6}{p_1^3\left(p_1^3+q_1^3\right)}.$$(A.2)

Further, $$p_1^6-q_1^6=\left(p_1^2-q_1^2\right)\left(p_1^4+p_1^2{q}_1^2+q_1^4\right),$$(A.3)

and the first factor is neatly transformed into $$\begin{array}{lll}{p}_1^2-q_1^2\hfill & =\hfill & \left(R^2+b_b^2\right)-\left(R_{\odot }^2+b_b^2\right)\hfill \\ \hfill & =\hfill & \left(R_{\odot }^2+2r\cdot R_{\odot }+r^2\right)-R_{\odot }^2\hfill \\ \hfill & =\hfill & 2r\cdot R_{\odot }+r^2.\hfill \end{array}$$(A.4)

Even if the two components in the last line have different signs, they are far from having the same order of magnitude. The same device has been applied for the f_d.

The treatment of f_h is also simple. The coefficient of R_⊙ involves $$1-\frac{\left(a_h+R\right)R}{\left(a_h+R_{\odot }\right)R_{\odot }}=\frac{a_h\left(R_{\odot }-R\right)+R_{\odot }^2-R^2}{\left(a_h+R_{\odot }\right)R_{\odot }},$$(A.5)

where $$R_{\odot }-R=\frac{R_{\odot }^2-R^2}{R_{\odot }+R},$$(A.6)

and $$R_{\odot }^2-R^2=-2r\cdot R_{\odot }-r^2.$$(A.7)

In the NFW model of the halo potential, one meets a new problem of subtracting logarithms with almost equal arguments or, equivalently,of a logarithm whose argument is a fraction with almost equal numerator and denominator. This obstacle is bypassed by the conversion: $$\begin{array}{lll}\ln \frac{a_H+R_{\odot }}{a_H+R}\hfill & =\hfill & \ln \frac{\left[\left(a_H+R\right)+\left(R_{\odot }-R\right)\right]\left(R+R_{\odot }\right)}{\left(a_H+R\right)\left(R+R_{\odot }\right)}\hfill \\ \hfill & =\hfill & \ln \left[1+\frac{R_{\odot }^2-R^2}{\left(a_H+R\right)\left(R+R_{\odot }\right)}\right],\hfill \end{array}$$(A.8)

and recalling Eq. (A.7), we find the form ln (1 + ξ), where |ξ| ≪ 1 does not involve a problematic subtraction. We note, however, that direct application of the natural logarithm function, implemented in standard compilers, may also degrade the accuracy. For small values of ξ, it is recommended to use specially designed functions like log1p in C/C++, or to resort to the power series expansion.

Appendix B Polynomial ephemeris of the Sun

For the cases where only the Galactic tides are to be taken into account in a small body motion in the Solar System, it is sufficient to have only the Sun’s galactocentric positions calculated in advance while using the formulae proposed in Sect. 2. The simplest and fasted way is to use polynomial approximations of these positions obtained on the basis of the precise numerical integration. In this appendix, we propose such ephemerides, both for backward and forward small body dynamics studies. The coefficients of the fitted polynomials (all of them of the fifth degree) are presented in Tables B.1 and B.2. In these tables, ephemerides for the past and future 30 Myr of the Sun’s galactocentric motion are presented. In both cases the whole interval is divided into six parts to ensure the similar accuracy of the whole ephemeris. The maximum difference between numerical integration results and positions obtained with these ephemerides is below 1 × 10⁻⁴ parsecs. It is worth recalling here that the distance of the Sun from the Galactic center is known with the uncertainty of hundreds of parsecs and for the vertical Sun position only the sign (positive) is certain. This makes the proposed ephemerides accuracy quite satisfactory.

The structure of both tables is identical: we present six polynomial segments covering the whole 30 Myr intervals. Each segment is described in the appropriate table with four rows: the start and the end of the validity time interval expressed in Myr, and three rows with the polynomial coefficients for X_⊙, Y_⊙, and Z_⊙ galactocentric positions of the Sun.

If we denote the coefficients in any row as a₀, a₁, a₂, …, a₅, then each coordinate in each interval can be obtained from the following formula: $$a(t_0+t)=(((((a_{5}t)+a_4)t+a_3)t+a_2)t+a_1)t+a_0,$$(B.1)

where for greater speed and better precision, we used the Horner scheme. In this formula, t₀ is the initial epoch of the time interval of interest, t is the time from the beginning of the interval to the point of interest, and a represent any of the X_⊙, Y_⊙ or Z_⊙ Galactocentric coordinate of the Sun at the epoch (t₀ + t).

It should be noted that both these ephemerides are obtained by polynomial least square fitting to the results from a precise numerical integration of the Sun and all 407 stars, based on a list from the StePPeD database but updated to the Gaia EDR3 results. All mutual interactions between stars were taken into account. As a result, the division of the whole intervals of 30 Myr is different and unequal because it depends on these mutual interactions during star-star encounters or star-sun interactions. We note that the computer readable files with these ephemerides are publicly available at the StePPeD database internet site.

References

All Tables

All Figures

	Fig. 1 Structure of the header part of the binary ephemeris file; “Serial” is the (4 byte integer) code of the ephemeris version (1255 at present), “number” is the (4 byte integer) number of bodies included. They are followed with the table of “number” rows, containing 8 byte character string with the perturber name and double precision floating point number containing the mass of the perturber.
In the text

	Fig. 2 Structure of the “one-step data blocks” of the binary ephemeris file; t1 and t2 are the ends of the validity time interval of polynomial coefficients stored in “body data blocks”. The number of these blocks equals the number of bodies included in the ephemeris. Each “body data block” is 240 bytes long (3 × 10 double precision floating point numbers).
In the text

Fig. 3 Distribution of the cloud of 16 545 clones of P0230 stopped at the closest proximity to the Sun about 2.80 Myr ago. All points are projected onto the maximum scatter plane X′ Y′. The green dot marks the nominal approach position, while the red one depicts the position of the Sun. The blue circle shows the widely adopted approximate boundary of the Oort cometary cloud at 0.5 pc from the Sun. Orange dot and arrows showthe position of C/2015 XY1 at the moment of the nominal P0230 approach; see Sect. 5 for details.

In the text

Fig. 4 Dynamical past evolution of the nominal orbit of C/2015 XY1. Here P0506 and P0230 are excluded from the model. 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. Angular elements are expressed in the Galactic frame.

In the text

	Fig. 5 Past evolution of C/2015 XY1 with only the star P0230 excluded. The strong interaction with P0506 is evident.
In the text

	Fig. 6 Past changes in nominal orbital elements of C/2015 XY1 with all potential stellar perturbers included. The strong interaction of C/2015 XY1 with P0506 and P0230 is shown.
In the text

Fig. 7 Distribution of the cloud of 25 000 clones of P0506 stopped at the closest proximity to the Sun 1.09 Myr ago. All points are projected onto the maximum scatter plane X′ Y′. The green dot marks the nominal approach position, while the red one depicts the position of the Sun. The blue circle shows the widely adopted approximate boundary of the Oort cometary cloud at 0.5 pc from the Sun. An inset at the upper right shows the star clones cloud enlarged ten times to show the position of C/2015 XY1 (an orange cross) at the closest approach.

In the text

	Fig. 8 Comparison of the cometary and stellar data uncertainty influence on the minimum star-comet distance distribution.
In the text

	Fig. 9 Effect of the cometary data uncertainty on C/2015 XY1 orbital elements recorded 2.0 Myr ago. See text for a detailed explanation.
In the text

	Fig. 10 Simultaneous effect of both cometary and stellar data uncertainties. The scales are the same as in Fig. 9.
In the text