© ESO, 2015

1. Introduction

The existence of analog comet clouds in different planetary systems spread out across the Galaxy is suggested by recent evidences of dusty excess on debris disks in exoplanetary systems observed by Spitzer at 70 μm (Greaves & Wyatt 2010). The Galactic tide may re-inject the Oort Cloud comets towards the inner part of the planetary system producing a comet flux with possible impacts on it. The effects of the Galactic tide on the Oort Cloud comets were studied by several authors (e.g., Heisler & Tremaine 1986; Matese & Whitman 1992; Fouchard et al. 2006). Usually they are separated into radial, transverse and orthogonal (to the disk plane) components. At solar distance from the Galactic center the orthogonal tide has been estimated to be about ten times more effective than the radial one, in perturbing the Oort Cloud comets into the solar system (Heisler & Tremaine 1986; Morbidelli 2005). At this distance, it appears to be the dominant perturber today and also over the future long term (Heisler 1990). In spite of that, the contribution of planar tides (radial and transverse) becomes non negligible when the distance of the parent star with respect to the Galactic center changes, as shown by Masi et al. (2009) for the Oort cloud comets and by Veras & Evans (2013) in their analysis on dynamics for extrasolar planets. In particular the planar component of the tide may cause a great change if the central star decreases its distance and/or its mass with respect to the present solar one. Moreover they will certainly be non negligible in presence of spiral arm perturbations, as also Kaib et al. (2011) have already pointed out and we confirm in our preliminary results (De Biasi 2014). Before the inclusion of orthogonal tides, only in this preliminary Part I, we will consider the sole planar tide effects produced by the main dynamical Galactic components: bulge, disk and dark matter halo. To probe how much the tidal outputs are system independent, we check the consistency among the descriptions led in three different reference systems: the Galactic one and two heliocentric systems with and without Hill’s approximation developed in 3D axisymmetric potential. A relevant fall-out of this analysis is related to the comparison of the tidal effects from the Galaxy axisymmetric components, separately considered: bulge, disk and dark matter halo. The Hill’s approximation shows the capability to guess a priori the importance of each of them simply looking at their different contributions to the total rotational velocity of the Sun on the Galaxy disk and at their corresponding contributions to the logarithmic velocity gradients.

Inside each reference system we consider different values of parameters: position inside Galaxy of the cloud central star (the analog of the Sun), rotational direction and initial Galactic longitude of the cometary orbit. The spiral arms introduction will severely affect comet orbits and the movement of the Sun causing a possible migration of it (Kaib et al. 2011) from the inner limit (≃4 kpc) to the present position on the Galaxy disk (8 kpc). To test preliminarily the consequences that a migration may entail on cometary dynamics, we will consider here the action of the Galactic tide at the extreme values of the previous range of distance: 4 and 8 kpc.

The discovery of extra solar planetary systems led to questioning about the special conditions in order to define a Galactic habitable zone (hereafter GHZ). The birth and the development of Life, as we know it, are based on a very fine balance between many factors and conditions that involve different areas of study, from chemistry to dynamics ones as we will address in the discussion.

Next sections will describe the model that has been built to reproduce the comet orbits under the Galactic tidal effects, in particular: Sect. 2 describes the way in which the Galactic potential has been shaped using its most important dynamical components: bulge, disk and dark matter halo; in Sect. 3 the three reference systems in use are explained, with particular attention to the advantages and the disadvantages of each of them; Sect. 4 is reserved to the numerical integration (initial conditions and integrators); Sect. 5 contains the outputs of the study, the discussion about the influence of the different parameters and the comparison among the reference systems; finally Sect. 6 is devoted to conclusions with some suggestions for future work.

2. Galaxy tide on Oort Cloud

Comets in the Oort Cloud are only weakly bound to the Sun’s gravity, and for this reason they are easily perturbed. The inner Oort cloud comets have aphelion Q< 20 000 AU, while the outer ones have Q> 20 000 AU. At large distances from the Sun, the Galactic tidal force (Masi et al. 2009; Matese & Lissauer 2004; Heisler & Tremaine 1986) is one of the dominant perturbers of comet orbits with stellar perturbations (Hills 1981; Matese & Lissauer 2002) and giant molecular clouds (GMCs; Bailey 1983). We will focus our attention to the first one neglecting, in this context, the other two.

Recent numerical simulations have demonstrated that the Sun may have migrated through the Galactic disk from an inner zone to the current one (Kaib et al. 2011). This means that the solar system could have been subjected to tidal perturbations that were significantly different from those which characterize its current Galactic environment (Sellwood & Binney 2002). This new perspective and the discovery of many extra solar systems, makes it necessary to study the influence that a Galactic environment, in positions closer to the Galactic center than the solar one. In order to simulate how the Galaxy tides may change their dynamical effects on comet orbits during migration we will consider two different limit distances from the Galactic center: 4 kpc and 8 kpc, assuming the last one as the current solar position.

2.1. The Galactic potentials

The starting point for the study of the dynamical behavior of a comet cloud in the Galactic potential is to build a realistic matter distribution for the most important dynamical components of the Milky Way. In the present study we adopt a simple azimutally symmetric Galactic model, and decompose the total mass distribution in its different Galactic components, i.e. the bulge, the disk and the dark matter (DM) halo. This approach does not take into account of irregularities and asymmetries in the mass distribution (like the spiral structure, the bulge triaxial form) since for the present aims, their effects can be neglected. This makes the mass model easier to integrate.

An axisymmetric model has been used for the total Galactic potential Φ(R,z), in which R is the Galactocentric radius and z the distance above the plane of the disk. The complete dynamical effect of the Galaxy could be expressed by the sum of the potentials due to bulge (Φ_BG), disk (Φ_D) and DM halo (Φ_DH): $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{Gal}}\mathrm{=}{\mathrm{\Phi }}_{\mathrm{BG}}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{D}}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{DH}}\mathit{.}& & \end{array}$(1)The bulge is one of the less extended components of our Galaxy. We may consider as a good approximation the spherical model of the Plummer sphere (Binney & Tremaine 2008), that has the following analytical form for its gravitational potential: $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{BG}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{BG}}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{r}}_{\mathrm{c}}^{\mathrm{2}}}}& & \end{array}$(2)where M_BG is the bulge total mass and r_c is the core radius (Flynn et al. 1996).

The disk presents a mixed composition in age, metallicity, velocity and velocity dispersion of its stars that complicates the disk modeling based on radial and vertical mass-luminosity distributions. A wide sequence of opportunities is available in the literature (Binney & Tremaine 2008) to model a three-dimensional axisymmetric disk distinguishing the two main contributions due to the thin and the thick disk as discovered by Burstein (1979). Because its existence as a distinct component appears to be not enough sure (Bovy et al. 2013) and owing to its contribution to the total surface density is only of about 5% at solar distance (Dehnen & Binney 1998; see De Biasi 2010, Chap. 2), as first approximation we limit ourselves only to consider the thin disk. We start from the main reference which is the infinitely thin Freeman model (Freeman 1970), having an exponential radial surface mass density distribution: $\begin{array}{ccc}\mathrm{\Sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}{\mathrm{e}}^{\mathrm{-}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{d}}}\mathit{,}& & \end{array}$(3)where Σ₀ ≃ 359 M_⊙ pc^-2 is the central disk surface mass density and R_d = 4.1 kpc is the disk scale length in agreement with the observation (Lewis & Freeman 1989) and with the surface density at the solar distance Σ₀(R_⊙) ≃ 50 M_⊙ pc^-2 (see Kuijken & Gilmore 1991).

The corresponding potential for the mass distribution (3) is: $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{D}}\mathrm{\left(}\mathit{R}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\mathit{\pi G}{\mathrm{\Sigma }}_{\mathrm{0}}\mathit{R}[{\mathit{I}}_{\mathrm{0}}\mathrm{(}\mathit{R}\mathit{/}\mathrm{2}{\mathit{R}}_{\mathrm{d}}\mathrm{)}{\mathit{K}}_{\mathrm{1}}\mathrm{(}\mathit{R}\mathit{/}\mathrm{2}{\mathit{R}}_{\mathrm{d}}\mathrm{)}\\ & & \mathrm{-}{\mathit{I}}_{\mathrm{1}}\mathrm{(}\mathit{R}\mathit{/}\mathrm{2}{\mathit{R}}_{\mathrm{d}}\mathrm{)}{\mathit{K}}_{\mathrm{0}}\mathrm{(}\mathit{R}\mathit{/}\mathrm{2}{\mathit{R}}_{\mathrm{d}}\mathrm{)}]\mathit{,}\end{array}$(4)where I_n and K_n are modified Bessel functions of different order n. Even if using Hankel transform and the cylindrical Bessel function of order zero, J_o, we may obtain a two-dimensional gravitational potential Φ_D(R,z), this kind of modeling is affected by two difficulties: the typical vertical extension of our comet orbits turns to be less than the scale height observed for the Galaxy disk, which is about 200 pc. The consequence is that only a fraction of the whole disk mass is involved in the tide; moreover this parametrization of the disk potential is slow to integrate. We may overcome both the problems combining three Miyamoto-Nagai disks (Miyamoto & Nagai 1975) of differing scale lengths and masses: $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{MN}}\mathrm{\left(}\mathit{R,z}\mathrm{\right)}\mathrm{=}\mathrm{-}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathrm{3}}\frac{\mathit{G}{\mathit{M}}_{\mathit{n}}}{\sqrt{{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathrm{[}{\mathit{a}}_{\mathit{n}}\mathrm{+}{\sqrt{{\mathit{b}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}}^{\mathrm{]}}}^{\mathrm{2}}}}\mathrm{·}& & \end{array}$(5)The parameter b is related to the disk scale height, a_n to the disk scale lengths and M_n are the masses of the three disk combined components. We list the values of the constants a_n, b, M_n, according with Flynn et al. (1996), in the Table 1. It needs to stress that the potential trend given by Eq. (4) results in fair agreement with that of Eq. (5) at the mid-plane (z = 0).

The Galactic mass distribution of the DM halo is still uncertain and there are several models in the literature. The general representation of the DM spherical distribution is given by the Zhao’s radial density profiles (Zhao 1996): $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\hspace{0.17em}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{S}}}{{\mathrm{(}{\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{H}}}}^{\mathrm{)}}}^{\mathit{\gamma }}{\left[\mathrm{1}\mathrm{+}{\mathrm{(}{\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{H}}}}^{\mathrm{)}}}^{\mathit{\alpha }}\right]}^{\frac{\mathit{\beta }\mathrm{-}\mathit{\gamma }}{\mathit{\alpha }}}}& & \end{array}$(6)where α, β and γ are three slope parameters, while r_H and ρ_S are two scale parameters¹. This pattern models the density as a double power law with a limit slope − γ toward the center of the halo (r → 0) and − β to the infinity (r → ∞).

Inside this general class of matter distribution we consider, instead of classical Navarro-Frenk-White (NFW; Navarro et al. 1997) density profile, the modified pseudo-isothermal (MPI) profile (Spano et al. 2008). That has been introduced to obtain the best fits of rotation curves in spirals and for LSB (low surface brightness galaxies), probably DM dominated. The choice between the two previous models for the DM halo lies inside the wide cusp/core chapter (Bindoni 2008). The MPI may be expressed as the general formulation of Zhao (Eq. (6)) with (α,β,γ) = (2,3,0).

The corresponding gravitational potential is given by: $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{DH}}^{\mathrm{MPI}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{0}}{\mathit{r}}_{\mathrm{H}}^{\mathrm{2}}\left(\frac{\ln \left(\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{H}}\mathrm{+}\sqrt{\mathrm{1}\mathrm{+}{\mathrm{(}{\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{H}}}}^{\mathrm{)}}}^{\mathrm{2}}}\right)}{\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{H}}}\mathrm{-}\frac{\mathrm{1}}{\sqrt{\mathrm{1}\mathrm{+}{\mathrm{(}{\frac{{\mathit{R}}_{\mathrm{vir}}}{{\mathit{r}}_{\mathrm{H}}}}^{\mathrm{)}}}^{\mathrm{2}}}}\right)\mathit{,}& & \\ & & \end{array}$(7)with ρ₀ the central dark halo density.

Fig. 1 Contributions to the circular velocity plotted vs. r from the different Galaxy components. DM halo has a MPI profile (see text).

In order to choose suitable values for the parameters that appear within the DM halo model, we follow the pattern proposed by Klypin et al. (2002), in particular the model A1 (with no exchange of angular momentum), bounded to the cosmological constraint of the collapse (or formation) redshift, z_F. In a CDM cosmology, with H₀ = 75 km s^-1 Mpc^-1 and Ω₀ = 1, the virial mass of the DM halo has been set at: M_vir = 10¹²M_⊙. Using the subroutine of Navarro et al. (1997), we obtain: z_F = 2.68, with the corresponding value for the overdensity δ_c = 1.179 × 10⁵ (in units of critical density at z = 0) for the NFW profile. Then also the MPI one is determined assuming both profiles reach the same density value at their corresponding scale radius r_H.

To conclude our Galaxy modeling we have to verify if the trend for the total circular velocity obtained turns to be in agreement with that observed in the solar environment. The values provided for the total circular velocity next to the solar system is below the mean observed value of V_tot(R_⊙) = 220 km s^-1 (see Table 1), however it is inside the measured error bar for these distances (Klypin et al. 2002).

The contributions to the total circular velocity trend vs. r for the different Galaxy components are plotted in Fig. 1.

It is to be noticed that while the tidal effect due to the DM halo appears to become relevant only at distance larger than 15 kpc (see Sect. 5.1) the heavy bulge here considered, plays a dominant role, in the plane, already at 4 kpc, as shown also by the mass distribution for each Galactic components in (Fig. 2). Modeling the Galactic potential at this shorter distance, turns then to be critical for our aims more than at largest ones. Undoubtedly the dynamical model of Bovy & Rix (2013), based on large sets of phase-space data on individual stars, has to be considered a very substantial contribution to build up the gravitational potential of the Milky Way in the range 4 kpc ≤ R ≤ 9 kpc. They modeled an Henquist bulge with lower mass than typically considered but their rotation curve at R ≥ 4 kpc is similar to our one obtained according to Binney & Tremaine (2008).

Fig. 2 Contributions to the mass distribution plotted vs. r from the different Galaxy components. DM halo has a MPI profile (see text).

3. Different reference systems

To describe the tidal effects of each Galactic component on the comet orbits, we consider three different reference systems for the perturbed motion equations:

• the inertial system with the origin in the Galaxy center (Masi et al. 2009);

• the pseudo-inertial system with the Sun as the origin, having the Galactic components as perturbation terms of the Sun-comet system (Danby 1988);

• the synodic system which has the Sun as origin and which is solar co-rotating (Binney & Tremaine 2008).

The co-rotating system has been considered in Hill’s approximation developed for an axisymmetric potential also in 3D-dimensions. The integration in different reference frames allow us to focus on distinct aspects of the Galaxy tidal perturbation, changing the point of view from which we can analyze the problem of the perturbed comet orbits. We can check how much the results are system-independent, the numerical accuracy of each of them and to test which is the best theoretical framework to describe the topic. In following sections we try to explain briefly the integration of motion equations in each different reference frame listed above. We have limited the integration to the planar tide imposing that the comet orbit remains on the Galactic equatorial plane as the solar one. Moreover no planetary effects are included and we do not treat encounters with GMCs or stellar passages.

3.1. Inertial system with origin on Galactic center

Introducing cartesian galactocentric coordinates (x,y,z), according to Eq. (1), the total contribution of the Galactic gravitational potential at the comet position (x_c,y_c) is: $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{Gal}}\mathrm{\left(}{\mathit{x}}_{\mathrm{c}}\mathit{,}{\mathit{y}}_{\mathrm{c}}\mathrm{\right)}\mathrm{=}{\mathrm{\Phi }}_{\mathrm{B}}\mathrm{\left(}{\mathit{x}}_{\mathrm{c}}\mathit{,}{\mathit{y}}_{\mathrm{c}}\mathrm{\right)}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{DH}}\mathrm{\left(}{\mathit{x}}_{\mathrm{c}}\mathit{,}{\mathit{y}}_{\mathrm{c}}\mathrm{\right)}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{D}}\mathrm{\left(}{\mathit{x}}_{\mathrm{c}}\mathit{,}{\mathit{y}}_{\mathrm{c}}\mathrm{\right)}\mathit{.}& & \end{array}$(8)The gravitational potential, due to Sun (or in general to the host star) at the galactocentric position (x_∗,y_∗) on the comet, is expressed by: $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{\ast }}\mathrm{(}{\mathit{x}}_{\mathrm{c}}\mathrm{-}{\mathit{x}}_{\mathrm{\ast }}\mathit{,}{\mathit{y}}_{\mathrm{c}}\mathrm{-}{\mathit{y}}_{\mathrm{\ast }}\mathrm{)}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{\sqrt{\mathrm{(}{\mathit{x}}_{\mathrm{c}}\mathrm{-}{\mathit{x}}_{\mathrm{\ast }}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{y}}_{\mathrm{c}}\mathrm{-}{\mathit{y}}_{\mathrm{\ast }}{\mathrm{)}}^{\mathrm{2}}}}\mathrm{·}& & \end{array}$(9)Then the total effect of the Galactic and stellar components on the cometary orbit is: $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{Tot}}\mathrm{\left(}{\mathit{x}}_{\mathrm{c}}\mathit{,}{\mathit{y}}_{\mathrm{c}}\mathit{,}{\mathit{x}}_{\mathrm{\ast }}\mathit{,}{\mathit{y}}_{\mathrm{\ast }}\mathrm{\right)}\mathrm{=}{\mathrm{\Phi }}_{\mathrm{Gal}}\mathrm{\left(}{\mathit{x}}_{\mathrm{c}}\mathit{,}{\mathit{y}}_{\mathrm{c}}\mathrm{\right)}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{\ast }}\mathrm{(}{\mathit{x}}_{\mathrm{c}}\mathrm{-}{\mathit{x}}_{\mathrm{\ast }}\mathit{,}{\mathit{y}}_{\mathrm{c}}\mathrm{-}{\mathit{y}}_{\mathrm{\ast }}\mathrm{)}\mathit{.}& & \end{array}$(10)Applying the usual equations of motion: a_c = −∇Φ_Tot, a_∗ = −∇Φ_Gal, for the comet and the star accelerations, respectively with the corresponding initial conditions, a differential system for each of them is built up describing their motion around the Galactic center. After the integration the galactocentric cometary orbit is transformed into an heliocentric one of coordinates (x_cH(t),y_cH(t)), making the following difference:$\begin{array}{ccc}{\mathit{x}}_{\mathrm{cH}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{x}}_{\mathrm{c}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}{\mathit{x}}_{\mathrm{\ast }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{;}{\mathit{y}}_{\mathrm{cH}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{y}}_{\mathrm{c}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}{\mathit{y}}_{\mathrm{\ast }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{.}& & \end{array}$(11)This treatment avoids the calculations of the non-inertial components due to the rotation of the Sun around the Galaxy, with the advantage to obtain the exact magnitude of Galactic tide effect on the comet orbit without assumption of any approximation. On the other side it suffers of a higher loss of significant numerical digits, because the heliocentric cometary motion is obtained by subtraction of two very close numbers, namely the Sun’s and the comet’s galactocentric position.

3.2. Pseudo-inertial system with origin on the Sun

The natural choice to study the motion of bodies internal to our solar system is an heliocentric framework. In this case the reference system is of course not inertial and the non-inertial forces need to be added.

The first approach we present is the pseudo-inertial one. Following Danby (1988; see De Biasi 2010) the dominant body becomes the Sun around which the comet describes a Keplerian motion at the presence of perturbing terms due to the direct and non-direct accelerations. The last ones include the fictitious forces justifying the name of the reference system.

To describe the cometary orbit we have to integrate two different equations: the first which yields the heliocentric cometary motion around the Sun perturbed by the presence of the Galaxy, $\begin{array}{ccc}{\begin{array}{c}{¨}\\ {r}\end{array}}_{\mathrm{C}}\mathrm{+}\mathit{G}\mathrm{(}{\mathit{m}}_{\mathrm{S}}\mathrm{+}{\mathit{m}}_{\mathrm{C}}\mathrm{)}\frac{{{r}}_{\mathrm{C}}}{{\mathit{r}}_{\mathrm{C}}^{\mathrm{3}}}\mathrm{=}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{G}}\mathrm{\left(}{{r}}_{\mathrm{CG}}\mathrm{\right)}\mathrm{-}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{G}}\mathrm{\left(}{{r}}_{\mathrm{G}}\mathrm{\right)}\mathit{,}& & \end{array}$(12)and the second one that expresses the motion of the Galactic components around the Sun perturbed by the comet (with absolutely negligible effects, being typically, m_C = 10^-17 ÷ 10^-14M_⊙): $\begin{array}{ccc}{\begin{array}{c}{¨}\\ {r}\end{array}}_{\mathrm{G}}\mathrm{+}\mathit{G}{\mathit{m}}_{\mathrm{S}}\frac{{{r}}_{\mathrm{G}}}{{\mathit{r}}_{\mathrm{G}}^{\mathrm{3}}}\mathrm{+}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{G}}\mathrm{\left(}{{r}}_{\mathrm{G}}\mathrm{\right)}\mathrm{=}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{G}}\mathrm{\left(}{{r}}_{\mathrm{GC}}\mathrm{\right)}\mathrm{-}\mathit{G}{\mathit{m}}_{\mathrm{C}}\frac{{{r}}_{\mathrm{C}}}{{\mathit{r}}_{\mathrm{C}}^{\mathrm{3}}}\mathit{,}& & \end{array}$(13)where the subscript C, G and S are referred to cometary, Galactic and solar positions respectively. Note that in Eq. (12) the Galactic perturbation term is still a difference between very close values. Therefore, this approach is subject to a numerical truncation as the previous one. To avoid this type of numerical problem we have applied a regularization, in particular we used the Potter’s reformulation of Encke’s method (Battin 1964). Through the use of regularization we can integrate the pseudo-inertial motion equations with the advantage of not introducing any approximations in the description of the comet orbit.

3.3. Co-rotating system in Hill’s approximation

Fig. 3 Picture of the initial location for a comet orbit with Galactic longitude $\frac{\mathrm{3}\mathit{\pi }}{\mathrm{2}}$ in the reference system with the origin on the Galactic center (X,Y). The heliocentric system (x,y) in Hill’s approximation is also shown.

The Hill’s approximation (Binney & Tremaine 2008) adapts the formalism of the restricted three-body problem for the case in which the dimension of satellite system is much smaller than the distance to the center of the host system. The previous condition suits our problem in which the system Sun-comet achieves the maximum size of 1 pc, while the distances from the Galactic center are about three orders of magnitude larger. In this type of situation we may also assume that the variation of the gravitational potential along the cometary orbits is very smooth, making possible the application of the distant-tide approximation (Binney & Tremaine 2008, Chap. 8).

According to their treatment a spherically symmetric host system has been considered, with a gravitational potential Φ(R) at the distance R from its own center. Of course our problem does not enjoy this kind of symmetry owing to the presence of the disk. So we have re-considered the Hill’s approximation in 3D-dimensions changing the symmetry of the host system potential from a spherical to an axial one (see Appendix A). We follow the analysis of the comet motion in a co-rotating Sun centered system in which the x-y plane coincides to the stellar orbital plane, $\hat{{e}}{\mathit{x}}_{}$ points directly away from the center of the host system and $\hat{{e}}{\mathit{y}}_{}$ points in direction of the orbital motion of the satellite (see Fig. 3).

This reference system, called synodic, rotates with the frequency ${{\Omega }}_{\mathrm{0}}\mathrm{\equiv }{\mathrm{\Omega }}_{\mathrm{0}}\hat{{e}}{}_{\mathit{z}}$ and in this system the acceleration of a particle is described by the following expression: $\begin{array}{ccc}\frac{\mathrm{d}\mathrm{2}{x}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\mathrm{=}\mathrm{-}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{T}}\mathrm{-}\mathrm{2}{{\Omega }}_{\mathrm{0}}\mathrm{\times }\frac{\mathrm{d}{x}}{\mathrm{d}\mathit{t}}\mathrm{-}{{\Omega }}_{\mathrm{0}}\mathrm{\times }\mathrm{(}{{\Omega }}_{\mathrm{0}}\mathrm{\times }{x}\mathrm{)}\mathit{,}& & \end{array}$(14)in which the effect of the Coriolis force (second term on the right side) is separated from the centrifugal one (third term). Φ_T is the total potential acting on the comet. It may be decomposed in two different parts as follows: $\begin{array}{ccc}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{T}}\mathrm{=}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{s}}\mathrm{+}\sum _{\mathit{j,k}\mathrm{=}\mathrm{1}}^{\mathrm{3}}{\mathrm{\Phi }}_{\mathit{jk}}{\mathit{x}}_{\mathit{k}}\mathit{x̂}\mathit{j}& & \end{array}$(15)where ∇Φ_s, equal to GM_⊙/r² in the solar case, represents the contribution of gravitational potential due to the satellite system, and the second term comes from the distant-tide approximation once expanded the gravitational potential of the host system Φ(R) in Taylor series, starting from the center-of-mass of satellite system. This term, developed in an axisymmetric potential, expresses the Galactic tide along the three directions (x,y,z) as variation (at the second order) of the gravitational force exerted by the host system, from the solar to the comet distance. By comparison with the previous formulation (Heisler & Tremaine 1986; Levison et al. 2001; Brasser et al. 2010), Eq. (14) turns out to include the contribution of the variation of centrifugal force to the Galactic tide, and moreover the Coriolis force in order to obtain the complete system of equations for the comet motion, according to Binney & Tremaine (2008) Hill’s approximation.

In our coordinate system the center of the host system is in X = (−R₀,0,0) and Ω₀ = (0,0,Ω₀). Introducing the satellite system of reference x = (x,y,z) (Levison et al. 2001) it turns that: ${\mathrm{\Phi }}_{\mathit{xx}}\mathrm{=}{\mathrm{\Phi }}^{\mathrm{\prime \prime }}\mathrm{\left(}{\mathit{R}}_{\mathrm{0}}\mathrm{\right)}\mathrm{;}{\mathrm{\Phi }}_{\mathit{yy}}\mathrm{=}\frac{{\mathrm{\Phi }}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{R}}_{\mathrm{0}}\mathrm{\right)}}{{\mathit{R}}_{\mathrm{0}}}\ne {\mathrm{\Phi }}_{\mathit{zz}}$²; Φ_xy = Φ_xz = Φ_yz = 0. Then, the motion equations may be expressed as follows (see Appendix A): $\begin{array}{ccc}\{\begin{array}{c}\\ & \\ \\ & \\ \\ & \end{array}& & \end{array}$(16)where, according to (A.16), , is the density in the Sun’s neighborhood at z = 0 (Brasser et al. 2010). Then we could use the relation ${\mathrm{\Phi }}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{R}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}{\mathit{R}}_{\mathrm{0}}{\mathrm{\Omega }}_{\mathrm{0}}^{\mathrm{2}}$, due to the assumption of a circular motion of the Sun around the Galactic center and introduce the Oort’s constants: $\begin{array}{ccc}\mathit{A}\mathrm{\left(}\mathit{R}\mathrm{\right)}& \mathrm{\equiv }& \frac{\mathrm{1}}{\mathrm{2}}\left(\frac{{\mathit{v}}_{\mathrm{c}}}{\mathit{R}}\mathrm{-}\frac{\mathrm{d}{\mathit{v}}_{\mathrm{c}}}{\mathrm{d}\mathit{R}}\right)\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{R}\frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}\mathit{R}}\mathit{,}\\ \mathit{B}\mathrm{\left(}\mathit{R}\mathrm{\right)}& \mathrm{\equiv }& \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{{\mathit{v}}_{\mathrm{c}}}{\mathit{R}}\mathrm{+}\frac{\mathrm{d}{\mathit{v}}_{\mathrm{c}}}{\mathrm{d}\mathit{R}}\right)\mathrm{=}\mathrm{-}\left(\mathrm{\Omega }\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{R}\frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}\mathit{R}}\right)\mathit{,}\end{array}$(17)so that: Ω_o = A − B is the angular speed of Galactic rotation measured at the distance R = R_o from the Galactic center.

3.3.1. On tide contributions in Hill’s approximation

Due to the additive contribution of each Galaxy component to the total gravitational potential (see Eq. (1)) and to the assumption of the Sun’s circular motion, it follows that the corresponding circular and angular velocity contribution of each dynamical component, v_ci and Ω_oi, respectively, is: $\begin{array}{ccc}\frac{{\mathit{v}}_{\mathrm{ci}}^{\mathrm{2}}}{{\mathit{R}}_{\mathrm{o}}}\mathrm{=}{\left(\mathrm{-}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathit{i}}\right)}_{{\mathit{R}}_{\mathrm{o}}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{oi}}^{\mathrm{2}}{\mathit{R}}_{\mathrm{o}}\mathrm{;}\mathit{i}\mathrm{=}\mathrm{BG}\mathit{,}\mathrm{D}\mathit{,}\mathrm{DH}\mathit{.}& & \end{array}$(18)Then, in the first equation of the system (16), the following transformation holds for the tidal term in square brackets as soon as it refers to a specific dynamical component: $\begin{array}{ccc}{\mathrm{\Omega }}_{\mathrm{oi}}^{\mathrm{2}}\mathrm{-}{\mathrm{\Phi }}_{\mathit{i}}^{\mathrm{\prime \prime }}\mathrm{=}{\left(\mathrm{-}\mathrm{2}{\mathrm{\Omega }}_{\mathit{i}}{\mathrm{\Omega }}_{\mathit{i}}^{\mathrm{\prime }}\mathit{R}\right)}_{{\mathit{R}}_{\mathrm{o}}}\mathrm{=}\mathrm{2}\frac{{\mathit{v}}_{\mathrm{ci}}^{\mathrm{2}}}{{\mathit{R}}_{\mathrm{o}}^{\mathrm{2}}}{\left(\mathrm{1}\mathrm{-}\frac{\mathrm{d}\ln \hspace{0.17em}{\mathit{v}}_{\mathrm{ci}}}{\mathrm{d}\ln \mathit{R}}\right)}_{{\mathit{R}}_{\mathrm{o}}}\mathrm{·}& & \end{array}$(19)The x-component of Galaxy tide turns to be depending on the circular velocity contribution of each Galactic component and on its logarithmic gradient. From the other side the y-component Galaxy tide (see second equation of system (16)), due to the corresponding variation of centrifugal force, turns out to be zero so that the term in brackets:$\begin{array}{ccc}{\mathrm{\Omega }}_{\mathrm{oi}}^{\mathrm{2}}\mathrm{-}\frac{{\mathrm{\Phi }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{R}}_{\mathrm{0}}\mathrm{\right)}}{{\mathit{R}}_{\mathrm{0}}}& & \end{array}$(20)vanishes no matter what Galactic component considered. From system (16) it is clear that, for a particle at rest in the synodic system and in a cylindrical potential, the planar tide is only due to the difference between the variation, along the x-axis, of the centrifugal force and of the gravitational one due to the host system (the term in square brackets in the first equation of system (16)). Conversely the analogous component along the y-axis is completely compensated.

This last formulation has many advantages: first of all it avoids the problem of the loss of significant numerical digits present in both of the previous approaches. Further, if the aim is the total Galactic contribution to the tide, we have only one system of differential equations to integrate. Indeed the contribution of all the Galaxy components is inside the total angular velocity Ω₀ thanks to the additivity of the potentials and to the assumption of circular motion of the Sun. For the same reasons an analogous set of equations holds for each Galactic component, allowing us to underline their specific dynamical contribution to the total tide. The equations of system (16) have only to be transformed taking into account Eqs. (18) and (19).

On the other side this treatment is an approximated one and the assumptions on which it is based could not be fulfilled in every cases of our analysis.

Finally, in order to compare the results obtained in the synodic system with those in the two previous cases, the classic rotation matrix has to be used: $\begin{array}{ccc}\mathrm{𝒜}\mathrm{=}\left(\begin{array}{c}\\ & \\ & \end{array}\right)\mathit{.}& & \end{array}$(21)

4. Integration

4.1. Initial conditions

We have considered a comet body as a test particle for the Galactic tide and chosen a comet belonging to the outer shell of the Oort Cloud, where the solar gravitational force is lower and the Galactic perturbations are then more evident. The comet has initial aphelion of 140 000 AU, initial perihelion of 2000 AU, inclination equal to zero and Galactic longitude which starts from $\frac{\mathrm{3}\mathit{\pi }}{\mathrm{2}}$ (see Fig. 3) covering 360° with a step of 15°. The motion direction may be direct or retrograde. For each reference system the orbit has been integrated for 100 Myr.

4.2. Integration strategy

The implementation of the comet orbit is obtained by using two different integrators. The first one based on a conventional, explicit Runge-Kutta with variable order and variable integration step size, from the Wolfram Mathematica library. The method selects by default the Runge-Kutta coefficient order according to the initial values and the local error tolerances, in balance against the work (number of function evaluations and integration steps) of the method for each coefficient set. The second integrator is based on the Taylor series, TIDES software (Barrio 2005), in particular on the evaluation for the time Taylor series of the variables obtained by an iterative way, using the decomposition of the derivatives by automatic differentiation method. It presents a variable order and a variable step size as the first method. An analysis of the accuracy of the integrators has been made on a total integration time of 1.8 Gyr, about half the Earth’s age, an appropriate time scale to study the evolution of the comet motion. The accuracy has been checked by comparing the perihelion positions of the integrated unperturbed Keplerian orbit and the analytic one for both integrators. The result of this analysis is shown in Fig. 4. This shows that the percent error order of magnitude on the integrated position is 10^-6−10^-7, increasing with the integration time, for the explicit Runge-Kutta, while the accuracy for TIDES is six order of magnitude higher. This underlines TIDES as a superb integrator for comet orbits compared to standard integration schemes. We met some disadvantages of this integrator in dealing with a more complicate mathematical expression for the potentials as it will be in the case of a 3D spiral arm potential perturbation.

Fig. 4 Accuracy for the integrator based on an explicit Runge-Kutta (top panel) and the Taylor series development (bottom panel). Both evaluations have been done by comparing the perihelion positions of analytic unperturbed Keplerian orbit and the integrated ones (see text).

We absolutely need this numerical integrator when we are dealing with time spans of the order of Gyr (Parts II and III) whereas in this preliminary Part I we will use both the integrators also in the case in which the typical time span is of the order of 100 Myr only to check the consistency of some relevant results (e.g. Sect. 5.2).

5. Discussion

Growing evidences about the possibility of a migration of our Sun and the already strengthen discovery of different exoplanets in the Galaxy, address many studies on the consequences that a different Galactic environment may produce on a planetary system as the solar one and an eventually cometary cloud around it. Brasser et al. (2010) investigate the formation processes and the internal structure of the Oort cloud by numerical simulations for different collocations of Sun (2, 8, 14 and 20 kpc). The crucial importance of the collocation with respect to the center of the Galaxy are also stressed in Veras & Evans (2013), that pointed out the not negligible effect of the planar components of the Galactic tide on the planetary dynamics for system collocated at distance less then 3.5 kpc. We focus our attention on the cometary evolution under the effects of a Galactic tide produced by a environment different from the solar current one, in order to understand how a different collocation of a planetary system may influence the cometary injection process due to the Galaxy. Our aim is to investigate the effects of the planar Galactic tide on the comet orbits in three different reference systems underlining the weight of each single Galactic component to the total perturbation effects. The conditions that could amplify or reduce the tidal perturbations turn out to be depending on the three main parameters: the star distance from the Galactic center, the Galactic longitude and the direction of motion of the comet orbit. Numerical integration of highly eccentric orbits of individual comets around a star, which is assumed to move on a circular Galactic orbit, are performed.

5.1. The role of different components and parameters

The major effect of the Galactic tide is the perihelion variation, which manifests itself through both the precession of the major axis and the variation of the perihelion distance q. To test the influence of the distance from the Galactic center on q, we have integrated the orbit with the same initial values quoted in Sect. 4.1 (longitude = 3π/ 2; direct motion) for different distances, decomposing the total tidal perturbation into the contributions of the single Galactic components. To take into account that the Sun (or the central star in the analogs of Oort Cloud) may have undergone a possible migration owing to the spiral arms effect (which will be introduced in Parts II and III) we investigate the tides at two limit distances r = 4, and 8 kpc, without considering shorter distances in order to preserve the spherical structure for the bulge.

Fig. 5 Zoom of the perihelion zone for the comet orbit at 8 kpc from the Galactic center (Hill’s approximation). The different contributions from the Galaxy components to the perihelion variation are shown. The order of time-sequence is marked by numbers in the “Total” case.

Fig. 6 Zoom of the perihelion zone for the comet orbit at 4 kpc from the Galactic center (Hill’s approximation). The different contributions from the Galaxy components to the perihelion variation are shown. The order of time-sequence is marked by numbers in the “Total” case.

The strongest perturbation effect is caused by the Galactic disk at 8 kpc (Fig. 5). The total Galactic tidal perturbation increases moving toward the center of the Galaxy (Fig. 6). Moreover for a comet orbit closer to the Galactic center (4 kpc) the strongest perturbation is due to the bulge, while the disk has only a secondary influence. Even if this result is shown for the Hill’s approximation it turns to be completely independent of the reference system we are considering (see also Sect. 5.2).

To understand how the role of the different Galactic components changes with the distance, we refer to the powerful framework of Hill’s approximation. In it we are able to loock clearly how are the tidal trends. The tidal term in square brackets in the first equation of the system (16) gives us the x-component of Galaxy tide (the y-component is always compensated) separated in each Galactic components by Eq. (19). In Fig. 7 the contribution to the total Galaxy tide due to the bulge, disk and DM halo is shown. The tide from the bulge dominates until about 7 kpc from the center, that of the disk is prevailing in the range of about: 7–15 kpc, whereas for larger Galactic distance the DM halo plays its dominant tidal effect over the previous ones.

It is remarkable that one may predict what will be the Galactic component producing the most tidal effect on the Oort Cloud analogs around new extra solar planetary systems, simply looking at the position of its central star on the plot of Fig. 1 translated by Eq. (19) into that of Fig. 7. That may lead to a redefinition of the boundaries for the Galactic habitable zone. Several Galactic-scale factors can influence indeed the belonging to GHZ for a star (Gonzalez et al. 2001; Gonzalez 2005; Lineweaver et al. 2004; Gowanlock et al. 2010): they include those which are relevant to the formation of planets, such as radial disk metallicity gradient, and moreover the events that can threaten Life such as nearby supernovae and gamma ray bursts. There are at least two ways in which they might be involved into the process of Life formation and development. On one side they could upset the planetary conditions necessary for Life, but on the other side they also might deliver water and prebiotic organic compounds (Parts II–III of the series).

The dependence on the Galactic longitude and on the direction of revolution has been tested. The tidal effects turn out to be different for direct or retrograde comet orbits (see Fig. 9). Indeed, for orbits with the same initial conditions, the Galactic tide can cause an increase or a reduction of the perihelion distance, with an associated change of the orbital shape due to the perturbation of the eccentricity, depending on the direction of revolution.

Fig. 7 Tides of Galaxy components (Eq. (19)): from bulge (red), disk (blue), DM halo (dark), vs. distance from the Galactic center (in kpc). Bulge dominates until about 7 kpc, disk in about the range: 7–15, DM halo over 15 kpc. An amplification factor of 103 has been used on the y-axis.

In order to quantify this variation as function of Galactic longitude and direction of motion, we have calculated the mean perihelion distance ${\overline{\mathit{q}}}_{\mathrm{G}}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ by averaging over an integration time of 1 Gyr (Kómar et al. 2009), for different longitudes and both direct and retrograde motion as $\begin{array}{ccc}{\overline{\mathit{q}}}_{\mathrm{G}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{t}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}{\mathit{q}}_{\mathrm{G}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{d}\mathit{t,}& & \end{array}$(22)where q_G(t) is the instantaneous perihelion distance. The analysis shows a strong dependence of the comet perihelion distance variation on the combination of Galactic longitude and direction of motion. As shown in Fig. 9, comet orbits with the same Galactic longitude, but opposite direction of motion, have opposite signs of mean perihelion variation. Note that this dependence is stronger as the Sun-comet system approaches the Galactic center.

Fig. 8 Comparison between a comet orbit (zoom of the perihelion zone) in the three different reference systems, for an integration time of 100 Myr, at 8 kpc from the Galactic center.

These features are not easy to explain analytically due the intricate form of the several tidal actions at work. However, if we stipulate to reduce the complexity of the present system to that of three point masses, we can try to explain qualitatively the behavior on the basis of a simple application of the Lagrange planetary equations (LPE). Following the approach developed by Kozai (1973) to study the lunisolar perturbations on an Earth satellite, consider the third-body disturbing function $\begin{array}{ccc}\mathrm{ℛ}\mathrm{=}\frac{\mathrm{15}}{\mathrm{8}}{{\mathit{n}}^{\mathrm{\prime }}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}{\mathit{e}}^{\mathrm{2}}\left(\cos \mathrm{2}{\mathit{\lambda }}^{\mathrm{\prime }}\cos \mathrm{2}\mathit{\omega }\mathrm{+}\sin \mathrm{2}{\mathit{\lambda }}^{\mathrm{\prime }}\sin \mathrm{2}\mathit{\omega }\right)\mathit{,}& & \end{array}$(23)which has been restricted to the leading long-period terms for motion on a plane. Here, λ′ in the longitude of the disturbing, or third body, in its circular orbit and its mean motion.

Fig. 9 Mean perihelion trend for different distances (R0 = 4, 8 kpc) from the Galactic center as function of longitude for the two different directions of motion.

Since it is well known that the semimajor axis is only affected by short-period perturbations, the time rate of change of the perihelion distance is . If we now use the LPE for the eccentricity (Murray & Dermott 2000) we find that the regulating equation for $\overline{\mathit{q}}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ is $\begin{array}{ccc}\frac{\mathrm{d}\hspace{0.17em}\overline{\mathit{q}}}{\mathrm{d}\hspace{0.17em}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathrm{15}}{\mathrm{4}}\frac{{{\mathit{n}}^{\mathrm{\prime }}}^{\mathrm{2}}}{\mathit{n}}\mathit{ae}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\sin \mathrm{2}\mathrm{(}\mathit{\omega }\mathrm{-}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{)}\mathit{,}& & \end{array}$(24)which integrates to $\begin{array}{ccc}\overline{\mathit{q}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\overline{\mathit{q}}}_{\mathrm{0}}\left[\mathrm{1}\mathrm{+}\frac{\mathrm{15}}{\mathrm{8}}\frac{{{\mathit{n}}^{\mathrm{\prime }}}^{\mathrm{2}}}{\mathit{n}\mathrm{(}\mathrm{\pm }\mathit{\omega ̇}\mathrm{-}\mathit{\lambda ̇}\mathrm{\prime }\mathrm{)}}\frac{\mathit{e}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\mathrm{)}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}}\cos \mathrm{2}\mathrm{(}\mathit{\omega }\mathrm{-}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{)}\right]\mathit{,}& & \end{array}$(25)where we have used ${\overline{\mathit{q}}}_{\mathrm{0}}\mathrm{=}\mathit{a}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{e}\mathrm{)}$ as the integration constant. The ± sign appearing before the rate of the argument of perihelion at the denominator of this equation accounts for the opposite directions of precession associated with direct (upper sign) and retrograde (lower sign) motion. The functional form of the mean perihelion distance exhibited by Eq. (25) is consistent with the behavior of ${\overline{\mathit{q}}}_{\mathrm{G}}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ shown in Fig. 9. In fact, it predicts two different mean levels of the average perihelion distance, depending on the direction of motion. At the same time it also shows that the mean perihelion distance $\overline{\mathit{q}}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ varies periodically with twice the longitude of the disturbing body. Finally, the amplitude of the perturbation in the perihelion distance is seen to depend directly on the square of the mean motion of the perturber, which is consistent with the increase in the perturbations as the system gets closer to the Galactic center. We conclude that the essential features of the dependence of the mean perihelion distance on the longitude and the distance of the disturbing Galactic center is qualitatively similar to the behavior of a particle on an eccentric orbit under the gravitational action of a third particle at different longitudinal positions.

5.2. Comparison among reference systems: numerical accuracy and systematic deviations

We test the agreement between the different reference systems to see how the outputs are system independent. To do this we have to evaluate: i) the numerical accuracy due to the integration methods; ii) the systematic deviations due to the difference among the reference systems. We start with the integration of a comet orbit at the actual distance of the solar system from the Galactic center (8 kpc) with a Galactic longitude $\frac{\mathrm{3}\mathit{\pi }}{\mathrm{2}}$ and with direct motion under the effects of the total Galactic tide. Looking at Fig. 8, at first sight it appears a quite satisfactory agreement among the outputs of the three reference systems on the perihelion zone, without any detectable difference between the numerical accuracy of Runge-Kutta and TIDES methods (integration time is short with respect to the time scale considered in Fig. 4).

But to quantify the possible systematic deviations among the references we go further as follows: the two first reference systems encompass no approximations, so from this point of view, they may be considered two exact equivalent descriptions of the comet movements. Conversely the Hill’s approximation, by which we get the best meaningful description, is intrinsically affected by the distant-tide approximation. In order to investigate more deeply at which extent the Hill’s approximation starts to be seriously affected, we consider only the comparison between the two heliocentric reference systems (pseudo-inertial and synodic in Hill’s approximation). We integrate now cometary orbits with initial conditions for aphelion and perihelion as the previous ones (see Sect. 4.1) but changing the Galactic longitude for each orbit by a step of 15° and considering two different sets of integrations at the star distance from the Galactic center: 8 kpc and 4 kpc.

The results for the distance of 8 kpc, are shown in Fig. 9, where we see relatively good agreement between the two considered systems. The difference between two series of mean perihelion distances are also evaluated and quantified in Fig. 10: at 8 kpc the relative deviation reaches the maximum of 5%. The disagreement begins to increase when we approach the Galactic center at 4 kpc. The relative discrepancy between the two sets reaches indeed the maximum value of 33%. We performed the integration using both Runge-Kutta and TIDES without significant change in the conclusion. It proofs that the gap between the two different sets is not due to the numerical accuracy but to the main assumptions present in the Hill’s approximation that becomes increasingly bad approaching the bulge. We have to underline that it does not mean the general tide trends of Fig. 7, performed into Hill’s approximation, fail. We have indeed checked that the dominance of one tide from a Galaxy component over the others is the same in all our reference systems. Indeed also at the critical value of 4 kpc the tidal bulge effect dominates with respect of that from the disk with a relative ratio of about 50. Then the Hill’s approximation reveals its limit in the orbit integration without affecting the general tidal description.

Fig. 10 Relative discrepancy in the mean perihelion distances between the two heliocentric reference systems at 8 kpc and at 4 kpc from the Galactic center. The black points are relative to the direct orbits, while the grey ones are relative to the retrograde orbits.

6. Conclusions and future work

The Oort Cloud, because of its peripheral position inside our planetary system, is particularly sensitive to the tidal action of the Galaxy especially when the Sun (or the central star) is located closer to the Galactic center. Then we examine if the real comet distribution, so strongly concentrated around a semimajor axis a = 10⁴ AU, could be the result of a strong depletion of the comet population already occurred due to the past injection effect produced by the migration of Sun from inside toward the present location as we will consider in the next papers of the present series.

We have simulated the orbits of comet clouds using three different reference systems and two different integrators, around a star with the same mass of the Sun but with different Galactic locations. In conclusion we can affirm that: i) the disk has the strongest weight in the tidal perturbation for the actual solar distance; ii) the Galactic planar tides become much more effective in perturbing comets as soon as the central star approaches the Galactic center; iii) under about 7 kpc the tide of bulge dominates; iv) the cometary behavior for different Galactic longitudes and rotation directions indicates that the perihelion distance and the shape of orbits are strongly influenced by these two factors. The mean perihelion trend exhibits a surprising high regularity as function of the orbit longitude with a sinusoidal period of π at fixed direction of motion. An opposite direction of precession causes a phase change of π/ 2 in respect to the previous one. All the three analysed parameters: distance, longitude and rotation direction appear to play an important role in order to make a comprehensive discussion about the Galactic tidal perturbation on comet orbits.

An interesting future prospective may be opened by the incoming GAIA data, that should provide a more complete and accurate census for motion and distribution of stars in the solar neighborhood, allowing also a better determination for the Galactic field.

Acknowledgments

We like to thank Roberto Caimmi for fruitful discussions, constructive comments and very helpful suggestions. Our special gratitude is to Chris Flynn who gave us the opportunity to improve and clarify better the relevance of some parts of the paper.

References

Appendix A: Hill’s approximation in 3D axisymmetric potential

First of all, referring to Fig. 3, we consider the general gravitational potential Φ due to the host system on the point ${X}\mathrm{\equiv }\mathrm{(}\hat{{X}}\mathrm{+}{x}\mathrm{)}$ of the satellite system which has the center of mass at $\hat{{X}}$ (see Fig. 8.7 of Binney & Tremaine 2008, Chap. 8). By expansion in Taylor series we get: $\begin{array}{ccc}\mathrm{\Phi }\mathrm{(}\hat{{X}}\mathrm{+}{x}\mathrm{)}\mathrm{\simeq }\mathrm{\Phi }\mathrm{\left(}\hat{{X}}\mathrm{\right)}\mathrm{+}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{3}}\mathrm{(}{x}{\mathrm{)}}_{\mathit{k}}{\left(\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial }{\mathit{X}}_{\mathit{k}}}\right)}_{\hat{{X}}}\mathrm{·}& & \end{array}$(A.1)In the case of the solar system inside the Galaxy it turns to be , in the coordinate system (X,Y,Z) having the origin into the Galaxy center. Taking into account the X_j (j = 1,2,3) component of the gradient ∇Φ, Eq. (A.1) yields: $\begin{array}{ccc}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial }{\mathit{X}}_{\mathit{j}}}\mathrm{\simeq }{\left(\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial }{\mathit{X}}_{\mathit{j}}}\right)}_{{{R}}_{\mathrm{o}}}\mathrm{+}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{3}}\mathrm{(}{x}{\mathrm{)}}_{\mathit{k}}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial }{\mathit{X}}_{\mathit{j}}\mathit{\partial }{\mathit{X}}_{\mathit{k}}}\right)}_{{{R}}_{\mathrm{o}}}\mathrm{=}& & \\ {\mathrm{(}{\mathrm{\Phi }}_{\mathit{j}}\mathrm{)}}_{{{R}}_{\mathrm{o}}}\mathrm{+}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{3}}\mathrm{(}{x}{\mathrm{)}}_{\mathit{k}}{\mathrm{(}{\mathrm{\Phi }}_{\mathit{jk}}\mathrm{)}}_{{{R}}_{\mathrm{o}}}\mathrm{·}& & \end{array}$(A.2)The first term of Taylor’s expansion disappears as soon as we refer the movement of an object in the satellite system to the center of mass of the same solar system, choosing it as origin of reference system x ≡ (x,y,z) (Fig. 3). Indeed for a system of β particles and total mass M_S, the contribution to the barycentre acceleration from the term in the second order derivatives, vanishes by the definition of center of mass. It remains only the first order term which gives: $\begin{array}{ccc}\mathrm{-}\frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{S}}}\sum _{\mathit{\beta }}\mathrm{(}{\mathit{m}}_{\mathit{\beta }}{\mathrm{\Phi }}_{\mathit{j}}\mathrm{)}\mathrm{=}\mathrm{-}{\mathrm{\Phi }}_{\mathit{j}}\mathit{.}& & \end{array}$(A.3)So that the acceleration on one of these particle once referred to the barycentre yields: $\begin{array}{ccc}{a}\mathrm{=}\mathrm{-}{\mathrm{\Phi }}_{\mathit{j}}\mathrm{-}\sum _{\mathit{j,k}}{\mathrm{\Phi }}_{\mathit{jk}}\mathrm{(}{x}{\mathrm{)}}_{\mathit{k}}\mathrm{+}{\mathrm{\Phi }}_{\mathit{j}}\mathrm{=}\mathrm{-}\sum _{\mathit{j,k}}{\mathrm{\Phi }}_{\mathit{jk}}\mathrm{(}{x}{\mathrm{)}}_{\mathit{k}}\mathit{.}& & \end{array}$(A.4)Considering then two new reference systems with the center on Sun: X^′,x, which differ from the old reference one simply by a translation along X, so that: $\begin{array}{ccc}{\mathit{X}}^{\mathrm{\prime }}\mathrm{=}\mathit{X}\mathrm{-}{\mathit{R}}_{\mathrm{o}}\mathrm{=}\mathit{x}\mathrm{;}{\mathit{Y}}^{\mathrm{\prime }}\mathrm{=}\mathit{Y}\mathrm{=}\mathit{y}\mathrm{;}{\mathit{Z}}^{\mathrm{\prime }}\mathrm{=}\mathit{Z}\mathrm{=}\mathit{z}\mathit{.}& & \end{array}$Due to the same differentials in the old and new systems we may calculate the force components of the host system, with the center located now at ${{X}}_{{c}}^{{\prime }}\mathrm{\equiv }\mathrm{(}\mathrm{-}{\mathit{R}}_{\mathrm{o}}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{)}$, over the solar system as follows: $\begin{array}{ccc}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial X}}\mathrm{\simeq }\mathit{x}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial }{\mathit{X}}^{\mathrm{2}}}\right)}_{{{R}}_{\mathrm{o}}}\mathrm{+}\mathit{y}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial X\partial Y}}\right)}_{{{R}}_{\mathrm{o}}}\mathrm{+}\mathit{z}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial X\partial Z}}\right)}_{{{R}}_{\mathrm{o}}}& & \\ \frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial Y}}\mathrm{\simeq }\mathit{y}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial }{\mathit{Y}}^{\mathrm{2}}}\right)}_{{{R}}_{\mathrm{o}}}\mathrm{+}\mathit{z}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial Y\partial Z}}\right)}_{{{R}}_{\mathrm{o}}}\mathrm{+}\mathit{x}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial Y\partial X}}\right)}_{{{R}}_{\mathrm{o}}}& & \\ \frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial Z}}\mathrm{\simeq }\mathit{z}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial }{\mathit{Z}}^{\mathrm{2}}}\right)}_{{{R}}_{\mathrm{o}}}\mathrm{+}\mathit{x}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial Z\partial X}}\right)}_{{{R}}_{\mathrm{o}}}\mathrm{+}\mathit{y}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial Z\partial Y}}\right)}_{{{R}}_{\mathrm{o}}}& & \end{array}$without distinguishing the derivatives respect to the system X from those respect to X^′.

We consider now the case of an axisymmetric potential: $\begin{array}{ccc}\mathrm{\Phi }& \mathrm{=}& \mathrm{\Phi }\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}\mathrm{;}\mathit{R}\mathrm{=}\sqrt{\mathrm{(}{\mathit{X}}^{\mathrm{\prime }}\mathrm{+}{\mathit{R}}_{\mathrm{o}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{Y}}^{\mathrm{\prime }\mathrm{2}}}\mathrm{;}{\mathit{R}}_{\mathrm{\odot }}\mathrm{=}{\mathit{R}}_{\mathrm{o}}\mathrm{;}\\ \frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial R}}& \mathrm{=}& {\mathrm{\Phi }}^{\mathrm{\prime }}\mathrm{;}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial }{\mathit{R}}^{\mathrm{2}}}\mathrm{=}{\mathrm{\Phi }}^{\mathrm{\prime \prime }}\end{array}$then: $\begin{array}{ccc}& & \frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial }{\mathit{X}}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial R}}\frac{\mathit{\partial R}}{\mathit{\partial }{\mathit{X}}^{\mathrm{\prime }}}\mathrm{=}{\mathrm{\Phi }}^{\mathrm{\prime }}\frac{{\mathit{X}}^{\mathrm{\prime }}\mathrm{+}{\mathit{R}}_{\mathrm{o}}}{\mathit{R}}\\ & & \frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial }{\mathit{Y}}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial R}}\frac{\mathit{\partial R}}{\mathit{\partial }{\mathit{Y}}^{\mathrm{\prime }}}\mathrm{=}{\mathrm{\Phi }}^{\mathrm{\prime }}\frac{{\mathit{Y}}^{\mathrm{\prime }}}{\mathit{R}}\\ & & \frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial }{\mathit{X}}_{\mathit{i}}^{\mathrm{\prime }}\mathit{\partial }{\mathit{X}}_{\mathit{j}}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{X}}_{\mathit{i}}^{\mathrm{\prime }}}\left({\mathrm{\Phi }}^{\mathrm{\prime }}\frac{{\mathit{X}}_{\mathit{j}}^{\mathrm{\prime }}}{\mathit{R}}\right)\mathrm{=}\left({\mathrm{\Phi }}^{\mathrm{\prime \prime }}\mathrm{-}\frac{{\mathrm{\Phi }}^{\mathrm{\prime }}}{\mathit{R}}\right)\frac{{\mathit{X}}_{\mathit{i}}^{\mathrm{\prime }}{\mathit{X}}_{\mathit{j}}^{\mathrm{\prime }}}{{\mathit{R}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathrm{\Phi }}^{\mathrm{\prime }}}{\mathit{R}}{\mathit{\delta }}_{\mathit{ij}}\mathrm{;}\\ & & \mathrm{(}\mathit{i,j}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{;}{\mathit{X}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{=}{\mathit{X}}^{\mathrm{\prime }}\mathrm{+}{\mathit{R}}_{\mathrm{o}}\mathit{,}{\mathit{X}}_{\mathrm{2}}^{\mathrm{\prime }}\mathrm{=}{\mathit{Y}}^{\mathrm{\prime }}\mathrm{)}\mathit{.} \end{array}$From the last relationships we obtain that on the equatorial plane: $\begin{array}{ccc}\mathrm{(}\mathrm{\Phi }{\mathrm{)}}_{\mathit{xx}}& \mathrm{=}& {\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial }{\mathit{X}}^{\mathrm{\prime }\mathrm{2}}}\right)}_{{\mathit{R}}_{\mathrm{\odot }}}\mathrm{=}{\mathrm{\Phi }}^{\mathrm{\prime \prime }}\mathrm{\left(}{\mathit{R}}_{\mathrm{o}}\mathrm{\right)}\\ \mathrm{(}\mathrm{\Phi }{\mathrm{)}}_{\mathit{yy}}& \mathrm{=}& {\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial }{\mathit{Y}}^{\mathrm{\prime }\mathrm{2}}}\right)}_{{\mathit{R}}_{\mathrm{\odot }}}\mathrm{=}\frac{{\mathrm{\Phi }}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{R}}_{\mathrm{o}}\mathrm{\right)}}{{\mathit{R}}_{\mathrm{o}}}\mathrm{;}\mathrm{(}\mathrm{\Phi }{\mathrm{)}}_{\mathit{xy}}\mathrm{=}\mathrm{(}\mathrm{\Phi }{\mathrm{)}}_{\mathit{yx}}\mathrm{=}\mathrm{0.}\end{array}$For the vertical z-component of gravitational force: $\begin{array}{ccc}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial Z\partial }{\mathit{X}}_{\mathit{j}}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathit{\partial }}{\mathit{\partial Z}}\left({\mathrm{\Phi }}^{\mathrm{\prime }}\frac{{\mathit{X}}_{\mathit{j}}^{\mathrm{\prime }}}{\mathit{R}}\right)\mathrm{=}\frac{{\mathit{X}}_{\mathit{j}}^{\mathrm{\prime }}}{\mathit{R}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial Z\partial R}}\mathrm{·}& & \end{array}$(A.13)It means: $\begin{array}{ccc}\mathrm{(}\mathrm{\Phi }{\mathrm{)}}_{\mathit{xz}}& \mathrm{=}& {\left(\frac{{\mathit{X}}^{\mathrm{\prime }}\mathrm{+}{\mathit{R}}_{\mathrm{o}}}{\mathit{R}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial Z\partial R}}\right)}_{{\mathit{R}}_{\mathrm{\odot }}}\mathrm{=}{\left(\frac{{\mathit{X}}^{\mathrm{\prime }}\mathrm{+}{\mathit{R}}_{\mathrm{o}}}{\mathit{R}}\frac{\mathit{\partial }{\mathrm{\Phi }}^{\mathrm{\prime }}}{\mathit{\partial Z}}\right)}_{{\mathit{R}}_{\mathrm{\odot }}}\\ \mathrm{(}\mathrm{\Phi }{\mathrm{)}}_{\mathit{yz}}& \mathrm{=}& {\left(\frac{{\mathit{Y}}^{\mathrm{\prime }}}{\mathit{R}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial Z\partial R}}\right)}_{{\mathit{R}}_{\mathrm{\odot }}}\mathrm{=}{\left(\frac{{\mathit{Y}}^{\mathrm{\prime }}}{\mathit{R}}\frac{\mathit{\partial }{\mathrm{\Phi }}^{\mathrm{\prime }}}{\mathit{\partial Z}}\right)}_{{\mathit{R}}_{\mathrm{\odot }}}\end{array}$which are both equal zero due to the assumption of a circular motion of Sun on the Galactic plane around the Galactic center: Φ′ = RΩ² and to the symmetry of the system respect to the equatorial plane according to which the force component (A.7) on the plane (z = 0) has to be zero.

Remembering Poisson’s equation it follows: $\begin{array}{ccc}\mathrm{(}\mathrm{\Phi }{\mathrm{)}}_{\mathit{zz}}\mathrm{=}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial }{\mathit{Z}}^{\mathrm{2}}}\right)}_{{\mathit{R}}_{\mathrm{\odot }}}\mathrm{\simeq }\mathrm{4}\mathit{\pi G}\mathit{\rho ̅}& & \end{array}$(A.16)where is the density in the Sun’s neighborhood at z = 0. From the Eqs. ((A.12), (A.16)) the breaking of spherical symmetry is manifest.

All Tables

All Figures

	Fig. 1 Contributions to the circular velocity plotted vs. r from the different Galaxy components. DM halo has a MPI profile (see text).
In the text

	Fig. 2 Contributions to the mass distribution plotted vs. r from the different Galaxy components. DM halo has a MPI profile (see text).
In the text

	Fig. 3 Picture of the initial location for a comet orbit with Galactic longitude $\frac{\mathrm{3}\mathit{\pi }}{\mathrm{2}}$ in the reference system with the origin on the Galactic center (X,Y). The heliocentric system (x,y) in Hill’s approximation is also shown.
In the text

	Fig. 4 Accuracy for the integrator based on an explicit Runge-Kutta (top panel) and the Taylor series development (bottom panel). Both evaluations have been done by comparing the perihelion positions of analytic unperturbed Keplerian orbit and the integrated ones (see text).
In the text

	Fig. 5 Zoom of the perihelion zone for the comet orbit at 8 kpc from the Galactic center (Hill’s approximation). The different contributions from the Galaxy components to the perihelion variation are shown. The order of time-sequence is marked by numbers in the “Total” case.
In the text

	Fig. 6 Zoom of the perihelion zone for the comet orbit at 4 kpc from the Galactic center (Hill’s approximation). The different contributions from the Galaxy components to the perihelion variation are shown. The order of time-sequence is marked by numbers in the “Total” case.
In the text

	Fig. 7 Tides of Galaxy components (Eq. (19)): from bulge (red), disk (blue), DM halo (dark), vs. distance from the Galactic center (in kpc). Bulge dominates until about 7 kpc, disk in about the range: 7–15, DM halo over 15 kpc. An amplification factor of 103 has been used on the y-axis.
In the text

	Fig. 8 Comparison between a comet orbit (zoom of the perihelion zone) in the three different reference systems, for an integration time of 100 Myr, at 8 kpc from the Galactic center.
In the text

	Fig. 9 Mean perihelion trend for different distances (R0 = 4, 8 kpc) from the Galactic center as function of longitude for the two different directions of motion.
In the text

	Fig. 10 Relative discrepancy in the mean perihelion distances between the two heliocentric reference systems at 8 kpc and at 4 kpc from the Galactic center. The black points are relative to the direct orbits, while the grey ones are relative to the retrograde orbits.
In the text