A&A 366, 418-427 (2001)
DOI: 10.1051/0004-6361:20000256

The orbital evolution of binary galaxies

R. Chan - S. Junqueira

1 - Observatório Nacional, Rua Ganeral José Cristino 77, São Cristóvão, CEP 20921-400, Rio de Janeiro, Brazil
2 - University of Massachusetts, Dept. of Astronomy, GRC, Amherst, MA 01003-4525, USA

Received 1 August 2000 / Accepted 9 November 2000

We present the results of self-consistent numerical simulations performed to study the orbital circularization of binary galaxies. We have generalized a previous model (Junqueira & de Freitas Pacheco 1994) and confirmed partially their results. The orbital evolution of pairs of galaxies is faster when we consider interacting pairs with contacting "live'' galaxy halos but the circularization time remains larger than the Hubble time. Besides, the time behavior of the orbits has changed in comparison with previous work because of tidal forces and dynamical friction acting on the halos.

Key words: galaxies: interactions - galaxies: kinematics and dynamics - galaxies: halos

1 Introduction

The physical interaction between pairs of galaxies plays a fundamental role in their dynamical evolution, since the present distribution of binary galaxies is a synthesis of pairs at various stages of orbits that are consequence of unknown initial conditions. Such evolution probably begins from the initial separation of the galaxy from the Hubble flow, reaches the turnabout point and begins the collision with another galaxy. The subsequent steps of this process depend on the initial orbital parameters and on the effect of dynamical friction due to the overlap of the dark halos. However, this scenario is very difficult to be verified by the observations, because the orbital periods of binary galaxies are very long and because we have only projected quantities.

Despite these difficulties, observationally the study of binary galaxies remains a powerful method for determining galaxy masses. The observations of these systems are limited to the radial velocity difference and to the projected pair separation thus, in general, a statistical analysis of the data for a given sample is required to take into account the projection effects. As a consequence, the results refer to the sample as a whole rather than any specific pair of galaxies. Assuming the two galaxies to be point masses in Keplerian orbits, we can write the total binary mass as

 \begin{displaymath}M_1+M_2={{V_z^2 R_{\rm p} } \over {G f(\alpha,e,i,\Omega)}},
\end{displaymath} (1)

where M1 and M2 are the masses of the galaxies, G is the gravitational constant, Vz is the radial velocity difference, $R_{\rm p}$ is the linear projected separation, $f(\alpha,e,i,\Omega)$ is a function that characterizes the properties of the orbit, $\alpha$ is the orbital phase angle, e is the eccentricity of the orbit, i is the inclination angle and $\Omega$ is the longitude of the pericenter. As stressed by Zaritsky & White (1994), this estimation method can be statistically robust if the binary sample has the following properties: (1) the binaries have random inclination angles, orbital phases and longitudes of the pericenter; and (2) the distribution of orbital eccentricities is known. Many works have used this technique in order to calculate the total binary mass (e.g., Schweizer 1987; Soares 1990; Honma 1999). However, de Freitas Pacheco & Junqueira (1988) have shown that the phase angle is not a statistically independent variable, a fact confirmed by Zaritsky & White (1994). Besides, Schweizer (1987) has shown that the eccentricity probability is given by $p(e){\rm d}e=2e{\rm d}e$ and Zaritsky & White (1994) have found a near-isotropic distribution of velocities, with eccentricities between 0.50 and 0.88. Sharp (1990), analyzing several statistical studies of binary galaxies, has emphasized that masses derived from binary galaxy studies depend critically on the orbital eccentricity. White et al. (1983) presented dynamical studies of binary galaxies and developed suitable models and analysis techniques, but as noticed by Zaritsky & White (1994), they were unable to constrain the orbital eccentricity in their sample, and this indeterminacy led to a large uncertainty in mass of their binary galaxy ensemble. Whereas statistical analyses of observational data indicate an average orbital eccentricity of 0.25 (Karachentsev 1990; Wiren et al. 1996), some authors (Chengalur et al. 1995; Chatterjee & Magalinsky 1997) suggest that the morphology of their studied pairs implies that the galaxies are on high eccentricity.

Many numerical works have demonstrated the importance of the dynamical friction for the physics of interacting pairs of galaxies (White 1983; Aguilar & White 1986; Borne et al. 1988; Hernquist & Weinberg 1989; McGlynn 1990; Barnes 1990; Prugniel & Combes 1992; Bartlett & Charlton 1995). The theoretical formula derived by Chandrasekhar (Binney & Tremaine 1987, p. 425) for the dynamical friction, corresponds to a systematic tendency of deceleration of a body in the direction of its motion. Such friction is a consequence of the fluctuating field acting on the body due to the gravitational interaction with near neighbors. Chandrasekar's formula describes a local phenomenon and when applied to a satellite orbiting around a massive galaxy and moving inside a homogeneous background, not always gives results compatible with the numerical experiments (Combes 1992). In previous work, Junqueira & de Freitas Pacheco (1994) have presented numerical experiments showing the dynamics of noncontacting pairs of galaxies in closed orbits. In their simulations, the orbital decay was studied considering a pair of two spherical elliptical galaxies. With static isothermal dark halos the dynamical evolution of the galaxies was obtained using a numerical code very similar to the "multiple three-body algorithm'' (Borne 1984), which was not a self-consistent method. Estimations of the circularization time were obtained determining the variation of the orbital eccentricity for two subsequent positions in the trajectory of the galaxies along the orbits. The results obtained with this method suggest that the circularization times are extremely long, certainly longer than the Hubble time. However, these results were not conclusive because it was not considered the self-gravity of the systems and the dynamical friction of the halos in contact.

The uncertainties of the orbital parameters claim for a more accurate study of the orbital evolution of binary galaxies, since the derived masses depend sensitively on the eccentricity and on the galaxy orbits themselves. One way to investigate this problem is to estimate, for a sample of simulated pairs of galaxies on initial elliptical orbits of different eccentricities, the time of circularization and the geometry of the orbits. The galaxy models should have "live'' dark halos in order to abandon the hypothesis of point masses that most of the works use in the statistical methods for determining the mass of binary galaxies.

In this work we present full self-consistent N-body simulations of binary galaxies in closed orbits performed in order to study the induced effects of dynamical friction on orbiting galaxy pairs. For this purpose we study the influence of the distinct dynamical parameters on the circularization time of the orbits, when two deformable galaxies with "live'' halo are considered. In this case, we expect that tidal forces and dynamical friction will play the major role in dissipating energy. The elliptical galaxies are modeled assuming an isothermal dark halo constituted of "dark particles'' and a luminous component constituted of "luminous particles''. With these new experiments we considered pairs in contact and we verified if these conditions can significantly affect the estimation of the circularization time and the orbits of the galaxies.

This work is organized as follows. In Sect. 2 we present the models of the elliptical galaxies and in Sect. 3 we describe the computational method. In Sect. 4 we describe the initial conditions of the pairs of galaxies. In Sect. 5 we present the eccentricity calculation methods. In Sect. 6 we obtain new estimations of the circularization time and in Sect. 7 we discuss the results and summarize the conclusions.

2 The models

Our simulations assume two possible different primary and secondary galaxies with N=6000 and $N=30\,000$ equal-mass particles, total mass of 2 1011 $M_{\odot}$ and 1012 $M_{\odot}$, whose parameters are shown in Table 1. Hereinafter the value of the gravitational constant is G = 1 and the Hubble constant is H0 = 100 kms-1/Mpc. Our units are [L]=1 Mpc, [V]=100 kms-1, [T]=h-1H0-1 = 9.778 109 years and [M]= 2.325 1012 $M_{\odot}$, where h=H0/100. In these units, the galaxy masses $2~10^{11}~M_{\odot}$ and $10^{12}~M_{\odot}$ are 0.086 and 0.430 mass units, respectively.

Table 1: Elliptical galaxies
Type $m_{\rm l}$ $n_{\rm l}$ $r_{\rm e}$ $r_{\rm l}$ $m_{\rm h}$ $n_{\rm h}$ $a_{\rm h}$ $r_{\rm h}$   $m_{\rm t}$ $n_{\rm t}$
g1 0.0286 2000 0.0013 0.0094 0.0574 4000 0.0026 0.0188   0.0860 6000
g2 0.1433 10000 0.0022 0.0160 0.2867 20000 0.0044 0.0320   0.4300 30000

$m_{\rm l}$ is the luminous mass, $n_{\rm l}$ is the number of luminous particles, $r_{\rm e}$ is the effective radius as defined in de Vaucouleurs' law, $r_{\rm l}$ is the luminous exterior cutoff, $m_{\rm h}$ is the halo mass, $n_{\rm h}$ is the number of particles in the halo, $a_{\rm h}$ is the characteristic radius of the halo, $r_{\rm h}$ is the cutoff radius of the halo, $m_{\rm t}$ is the total mass of the galaxy and $n_{\rm t}$ is the total number of particles.

To generate the initial conditions, the spatial distribution of the luminous particles of the galaxies follows the de Vaucouleurs' law (Hernquist 1990) and the dark halo particles follow an isothermal density distribution (van Albada et al. 1985). Although the density of the Hernquist model falls rapidly in its outer regions, there are a few particles at extremely large distances. This causes numerical artifacts such as the inner core of the galaxy being far removed from its center of mass. Hence, we truncated the Hernquist model at $r_{\rm l} = 160$ kpc, for the most massive galaxy. We have also truncated its dark halo at $r_{\rm h} = 320$ kpc.

We have used a 1/3 power scaling with mass ensuring that the galaxies within the cutoff radius have similar density. This is motivated from the cosmological consideration that only objects that collapse at the same time have the same mean density within their turnaround radius (Gunn & Gott 1972). Although the currently favored cosmological paradigm of hierarchical clustering suggest that smaller objects will tend to be denser (Navarro et al. 1996, 1997). Thus,

 \begin{displaymath}r_{\rm h}(g_1) = r_{\rm h}(g_2) \left[ m_{\rm t}(g_1) \over m_{\rm t}(g_2) \right]^{1/3},
\end{displaymath} (2)

where g2 denotes the most massive galaxy.

We have also assumed that the mass-luminosity ratio of the galaxies is constant and equal to $m_{\rm h}/m_{\rm l} = 2$. This value has been chosen in such a way that the number of the luminous particles of the less massive galaxy not be very low. A very low number of luminous particles could affect the statistics of the results by numerical noise.

The density of luminous particles is given by (Hernquist 1990)

 \begin{displaymath}\rho_{\rm l}(r) = {{m_{\rm l} \over 2\pi a_{\rm l}^3} {1 \ove...
...over a_{\rm l} \right)\left(1+{r \over a_{\rm l}}\right)^3}}},
\end{displaymath} (3)

where $a_{\rm l}=r_{\rm e}/1.8153$, $r_{\rm e}$ is the effective radius as defined in de Vaucouleurs' law and $m_{\rm l}$ is the total luminous mass.

The density of dark halo particles is assumed to be isothermal (van Albada et al. 1985) and is given by

 \begin{displaymath}\rho_{\rm h}(r) = {{\rho_0} \over {1+\left({r \over a_{\rm h}}\right)^2}},
\end{displaymath} (4)

where $\rho_0$ is the central density. The central density of the halo can be written as

 \begin{displaymath}\rho_0 = {{m_{\rm h}} \over {4\pi a_{\rm h}^3 \left[ {r_{\rm ...
...\arctan \left(
{r_{\rm h} \over a_{\rm h}} \right) \right]}},
\end{displaymath} (5)

where $m_{\rm h}$ is the total halo mass, $r_{\rm h}$ is the cutoff radius and $a_{\rm h}$ is the core radius.

The Jeans' equation for a spherical system with isotropic velocity distribution (Binney & Tremaine 1987, p. 204) is given by

 \begin{displaymath}{1 \over \rho(r)} {{\rm d} \over {\rm d}r} \left[ \rho(r) \bar v^2(r) \right] = -{G M(r) \over r^2},
\end{displaymath} (6)

where $M(r)=4\pi \int^r_0 \rho(s) s^2 {\rm d}s$. Since the galaxy has two different component, we have to substitute $\rho(r)$ and M(r) by

 \begin{displaymath}\rho(r)= \rho_{\rm l}(r) + \rho_{\rm h}(r),
\end{displaymath} (7)


 \begin{displaymath}M(r)= M_{\rm l}(r) + M_{\rm h}(r),
\end{displaymath} (8)


 \begin{displaymath}M_{\rm l}(r)=4\pi \int^r_0 \rho_{\rm l}(s) s^2 {\rm d}s,
\end{displaymath} (9)


 \begin{displaymath}M_{\rm h}(r)=4\pi \int^r_0 \rho_{\rm h}(s) s^2 {\rm d}s.
\end{displaymath} (10)

Using Eqs. (6-10) we obtain that the velocity dispersion is given by
$\displaystyle \bar v^2(r)$ = $\displaystyle {G \over {\rho_{\rm l}(r) + \rho_{\rm h}(r)}}$  
    $\displaystyle \times \int^\infty_r {[\rho_{\rm l}(s)+\rho_{\rm h}(s)]
[M_{\rm l}(s)+M_{\rm h}(s)] \over s^2}{\rm d}s.$ (11)

Once integrated the Eq. (11) for each radius r, we constructed ten galaxies following the density distribution Eq. (7) and following this velocity dispersion. After that, we calculated the virial parameter $\kappa=2K/\vert W\vert$, where K is the total kinetic energy and W is total potential energy. We have considered the galaxy distribution with the parameter $\kappa$ as close as possible to 1. In general, the best value of $\kappa$ is between 0.99 and 1.01. The distinction between luminous and dark particles is made because we only want to study the orbital evolution of the luminous part of the galaxy, since the dark halo is not observable.

3 The computational method

The time evolution was carried out using a version of the fully vectorized Tree-Code (Hernquist 1987). We assume, in this work, the tolerance parameter $\theta$ the value 0.577, in order to avoid pathological situations when using the Tree-Code as described by Salmon & Warren (1994). We have also used the quadrupole corrections to the gravitational force.

The Tree-Code does not exactly conserve energy or momentum, but the energy is conserved to better than 0.3% over an entire evolution, and the center of the mass moves a distance of at most $2 \epsilon$ (the softening parameter) from the initial system center, with a time step size $\Delta t = 5~10^{-5}$ and $\epsilon=8~10^{-4}$. This softening scale is motivated by a tradeoff between accuracy and computational speed. We can note from Table 1 that this softening scale is smaller than the effective radius $r_{\rm e}$, thus the luminous cores of the galaxies have being followed accurately.

The Tree-Code, as parameterized requires 1000 Cray J90 CPU hours for $N = 60\,000$ particles of two massive galaxies, evolving the system to one Hubble time. Because of this large necessary amount of computing time we had to trade off between the total possible number of particles, experiments and the total CPU time available.

We remark that we have evolved the galaxies alone during one Hubble time in order to confirm that they are stable.

4 The initial conditions

The initial conditions were chosen in such a way that we start with a primary galaxy model and changed all the possible quantities, such as the initial pericentric distance, the initial eccentricity, in order to form the initial pair. There are two model of galaxies whose physical properties are summarized in Table 1. The initial pericentric distance is assumed to be a multiple factor (two or four times) of the sum of the primary galaxy halo radius and the secondary galaxy halo radius. The chosen initial eccentricities were 0.2, 0.4 and 0.6. We have performed some experiments with 0.9 but the initial apocentric distance was so large that the orbits of the galaxies were very short, during a Hubble time.

The initial position of the secondary is given by

 \begin{displaymath}{r_{\rm a}}(t=0) = a(1+e),
\end{displaymath} (12)

where the major semi-axis is $a=r_{\rm p}/(1-e)$, e is the eccentricity and $r_{\rm p}$ is the pericenter distance.

The initial tangential velocity is given by

 \begin{displaymath}{v_{\rm a}}(t=0) =\sqrt {{{(M_1+M_2)} \over a} \left( {1-e} \over {1+e} \right)}.
\end{displaymath} (13)

In Table 2 we can see all the initial conditions of the simulations - the initial apocentric distances $r_{\rm a}$ and the initial apocentric velocities $v_{\rm a}$ of the pairs.

Table 2: Initial conditions
No. prim. gal. sec. gal. $n_{\rm p}$ $r_{\rm p}$ e $r_{\rm a}$ $v_{\rm a}$
1 g1 g1 2 0.0752 0.6 0.3008 0.4783
2 g1 g1 2 0.0752 0.4 0.1755 0.7669
3 g1 g1 2 0.0752 0.2 0.1128 1.1045
4 g1 g1 4 0.1504 0.6 0.6016 0.3382
5 g1 g1 4 0.1504 0.4 0.3509 0.5423
6 g1 g1 4 0.1504 0.2 0.2256 0.7810
7 g1 g2 2 0.1016 0.6 0.4064 0.7127
8 g1 g2 2 0.1016 0.4 0.2371 1.1428
9 g1 g2 2 0.1016 0.2 0.1524 1.6458
10 g1 g2 4 0.2032 0.6 0.8128 0.5039
11 g1 g2 4 0.2032 0.4 0.4741 0.8081
12 g1 g2 4 0.2032 0.2 0.3048 1.1638
13 g2 g2 2 0.1280 0.6 0.5120 0.8197
14 g2 g2 2 0.1280 0.4 0.2987 1.3144
15 g2 g2 2 0.1280 0.2 0.1920 1.8930
16 g2 g2 4 0.2560 0.6 1.0240 0.5796
17 g2 g2 4 0.2560 0.4 0.5973 0.9294
18 g2 g2 4 0.2560 0.2 0.3840 1.3385

$n_{\rm p}$ is the number which multiplies the sum of the halo radius of the primary ( $r_{{\rm hp}}$) and the secondary galaxy ( $r_{{\rm hs}}$), giving the pericenter distance. Thus, $r_{\rm p}=n_{\rm p}(r_{{\rm hp}}+r_{{\rm hs}})$.

The initial state of the halos is unrealistic because there is a sharp boundary between the galaxy and empty space, but the initial separation was always larger than the combined halo sizes to allow them to relax to a smoother more realistic profile before they begin to overlap.

\end{figure} Figure 1: Snapshots at the last pericenter time for differentexperiments
Open with DEXTER

\end{figure} Figure 2: Time evolution of the distances between the centers of the primary and secondary galaxies for the experiments of Fig. 1
Open with DEXTER

The overlapping of the halos can be seen in Fig. 1 where we show the snapshots of four different experiments at their last pericentric position. The last pericenter time is determined using Fig. 2 and estimating the last minimum in the distance between the primary and secondary galaxies. We call attention to the fact that the orbits are different in comparison with the ones obtained by Junqueira & de Freitas Pacheco (1994). In the latter, the apocenter and pericenter positions were located always along the x-axis, with phase angles $\alpha=0$ or $\alpha=\pi$, depending on the orientation of the orbit. However, in our case the apocenter and pericenter times are not at these angles anymore, fact that can be easily seen in Fig. 1, which shows snapshots at the pericenter time (Exp. 3 and 15). Thus, the orbits of the galaxies are not well behaved as the ones in Junqueira & de Freitas Pacheco (1994), where the foci of the "instantaneous'' elliptical orbits remained at rest at the center of mass of the systems. In this work, as can be seem in Sect. 5.2, we had to generalize the time evolution of the semi-axis a and the coordinates of the focus of the orbit, not only the eccentricity, in order to fit adequately the orbits. This orbital behavior can be understood as a precession of the pericenter and of the apocenter position due to an extend spherical mass distribution (Binney & Tremaine 1987, p. 106). Besides the tidal forces and the dynamical friction may play an additional roles in the temporal evolution of the orbits.

5 The eccentricities

In order to calculate the time evolution of the eccentricity, we have compared two different methods. The first method is the same presented in Junqueira & de Freitas Pacheco (1994). Self-consistent simulations generate new features that did not appear in the experiments with "multiple three-body method'', because the limitations of this kind of code constraint us to non-contact pairs. In order to take into account these features, we have developed the second method that is a generalization of the previous one. Both methods are described below.

Due to the tidal interaction between the galaxies, the center of mass of the luminous particles has being far removed from the highest density center. Thus, in all measurements we take hereinafter as the center of the system the "density center'' (Farouki & Salpeter 1994) defined by

 \begin{displaymath}\vec r_{\rm c} = {{{\sum_{i=1}^{n_{\rm l}}} w_{i} \vec r_{i}} \over {{\sum_{i=1}^{n_{\rm l}}} w_{i}}},
\end{displaymath} (14)

where the weight wi for each luminous particle i is taken to be the third-largest value of |ri - rj|-3 over all other luminous particles $j \ne i$.

\end{figure} Figure 3: The orbits of the binary galaxies, for the experiments  1-18. The solid line represents the orbit of the primary galaxy and the dashed line represents the secondary galaxy, where $r_{\rm a}$ is the apocentric distance
Open with DEXTER

In the Fig. 3 we show the orbital evolution of the 18 pairs of galaxies during a Hubble time, whose initial conditions are shown in Table 2. The experiments in Table 2 (numbers 1-6) are very similar to the experiments numbers 13-18 but they are not simply repeats. The experiments numbers 1-6 have the less massive galaxy as primary and secondary, while in the experiments numbers 13-18 have the more massive galaxy as primary and secondary galaxies. The effects of this difference appear more clearly comparing the experiments number 3 and number 15, where we have more orbital periods. We can notice that the orbit in solid line in experiment 3 goes further along than the solid line orbit in experiment 15. This effect is due to the action of dynamical friction, which is greater in experiment 15 (more massive galaxy) than in experiment 3 (less massive galaxy). This can be seen clearly from the circularization time that it will be discussed in Sect. 6.

5.1 First method

The first method (Junqueira & de Freitas Pacheco 1994) is based in the equation of an ellipse in polar coordinates, given by

 \begin{displaymath}r={{a(1-e^2)} \over {1+e \cos \alpha}}\cdot
\end{displaymath} (15)

\end{figure} Figure 4: The time evolution of the eccentricities of the orbits up to one Hubble time, for the experiments 1-18. Dots represent the results obtained using the first method as described in the text. The straight solid lines represent the fitted curves to these points
Open with DEXTER

For two very near subsequent positions in the orbit and using Eq. (15) we can calculate the time evolution of the eccentricity, assuming that the semi-axis a and the eccentricity e are constant for these two points. Thus, we obtain

 \begin{displaymath}e={{r_{k+1} - r_k} \over {r_k \cos \alpha_{k} - r_{k+1} \cos \alpha_{k+1}}},
\end{displaymath} (16)

where the subscripts k and k+1 represent any two subsequent points in the orbit. Before utilizing the Eq. (16), in order to obtain the eccentricities, the orbits have been smoothed using a cubic spline interpolation method (Press et al. 1992). The first method has been applied in Junqueira & de Freitas Pacheco (1994) successfully, first, because the simulated systems were not in physical contact (condition imposed by the non-self-consistent code) and, second, because their halos were static, not evolving in time. Due to these facts the orbital decay rate was not very high and it did not appear the strong discontinuities in the time evolution of the eccentricities around the apocenter and pericenter time. However, as can be seen from Fig. 4, this is not an efficient method to obtain the time evolution of the eccentricity for all experiments when we consider contact pairs in self-consistent simulations. In these figures dots represent the eccentricities calculated using the Eq. (16) and straight lines are the fitted curves, excluding the points with 0>e>1. Considering orbits shown in Fig. 3 and the eccentricities in Fig. 4, we verified for some experiments (for example in the experiment 3), although the orbit at the end of the evolution is visibly almost circular, the eccentricity increases with the time. In order to avoid this problem we have developed another more general method to obtain the eccentricities of the orbits.

5.2 Second method

The second method consists in assuming that the major semi-axis a and the eccentricity e vary linearly with the angle $\alpha=\alpha(t)$ in Eq. (15). Another hypothesis is that the focus of the ellipse change also linearly on the orbit plane with $\alpha (t)$. Thus,

 \begin{displaymath}a(t)=a_0+a_1 \alpha(t),
\end{displaymath} (17)

 \begin{displaymath}e(t)=e_0+e_1 \alpha(t),
\end{displaymath} (18)

 \begin{displaymath}x_{\rm f}(t)=x_{0}+x_{1} \alpha(t),
\end{displaymath} (19)

 \begin{displaymath}y_{\rm f}(t)=y_{0}+y_{1} \alpha(t),
\end{displaymath} (20)

where a0 is the initial major semi-axis, e0, is the initial eccentricity, x0=0 and y0=0 are the initial position of the focus of the orbit, the quantities a1, e1, x1, y1 are the position angle rate of these quantities.

In order to verify how each of these parameters affect the orbits, and consequently the eccentricity calculation, we have constructed four synthetic orbits changing one of the parameters and fixing the other ones. The synthetic orbits are shown in Fig. 5 for an ellipse with major semi-axis a=1 and eccentricity e=0.5. In Fig. 5a we show the orbit where the major semi-axis evolves in time with a(t)=1-2 10-4 t and the eccentricity and the position of the focus are remained static. In Fig. 5b the eccentricity of the orbit changes with e(t)=0.5-2 10-4 t and the major semi-axis and the focus position do not evolve in time. In Fig. 5c the focus of the ellipse with a=1 and e=0.5 change in time with $x_{\rm f}=-2~10^{-4}~t$ and $y_{\rm f}=0$. Finally, in Fig. 5d we present the orbit where all the three quantities, i.e., a, e and $x_{\rm f}$ evolve in time with the same orbital parameter rate.

\end{figure} Figure 5: Theoretical curves as described in the text. The synthetic orbits begin always at the position $x(0)/r_{\rm a}=-1.5$ and $y(0)/r_{\rm a}=0$
Open with DEXTER

\end{figure} Figure 6: The eccentricities calculated for the synthetic curves shown in Fig. 5, considering a time rate of the orbital parameters, a, e and $x_{\rm f}$, multiplied by a factor of 0.2
Open with DEXTER

In order to distinguish between the synthetic orbits and the results of the binary galaxy simulations, hereinafter, we will call the latter one by real data. For the synthetic orbits we calculated the time evolution of the eccentricities utilizing the first method, as we have done for the real eccentricities. The synthetic eccentricities are shown in Figs. 6a-6d, where we have decreased the rates by a factor of 0.2 for the sake of comparison with the real eccentricities. Several new features are shown, for example, the discontinuity around the pericenter and apocenter time. The larger are the time rates of the orbital parameters, a, e and $x\rm _f$, the larger are the discontinuities. Comparing the time evolution of the eccentricities, the real ones with the synthetic ones, we can notice that the former eccentricities do not present exactly the same behavior as the synthetic eccentricities do. This is due to the fact that the real orbits have a combination of features of the synthetic orbits, where several orbital parameters vary simultaneously. We point out that the Fig. 5d shows a synthetic orbit where the three orbital parameters change simultaneously with the same time rate. However, we notice that this situation is an ideal one. In the real orbits all the orbital parameters may vary simultaneously but with different time rates. The second method consists in fitting the orbits allowing that each orbital parameter varies independently.

\end{figure} Figure 7: The solid line represents the orbit of the primary galaxy of the experiment 9 and the dashed line represents the fitted curve using the second method as described in the text
Open with DEXTER

As an example of a fitted orbit, we show the orbit of the primary galaxy of the experiment 9 in Fig. 7. We call attention that, as can be seen in Table 3, other experiments have a better fit with a smaller $\chi^2$ but these cases are the simulations whose orbits have not completed a period.


Table 3: Fitted parameters
No. $a_1~({\rm rad}^{-1})$ $e_1~({\rm rad}^{-1})$ $x_{1}~({\rm rad}^{-1})$ $y_{1}~({\rm rad}^{-1})$ $\chi^2$ $a_{\rm e}$ $b_{\rm e}$ $t_{\rm c}$
1 -2.7441 10-3 -1.7887 10-2 -7.1571 10-4   3.5709 10-4   1.8385 10-7 0.6000 -0.1041 5.763
2 -2.1289 10-3 -8.2113 10-3 -5.0706 10-5   2.4677 10-4   7.6612 10-7 0.3484 -0.09325 3.736
3 -1.4601 10-3 -1.4381 10-2 -4.0039 10-4   1.2020 10-4   1.0542 10-6 0.1097 -0.2824 0.4883
4   2.1685 10-3   1.1322 10-3 -2.8417 10-3   1.9706 10-4   3.1501 10-8 0.6071 0.007169  
5   2.7160 10-3   5.4489 10-4 -3.9955 10-4   3.9715 10-4   2.7562 10-8 0.4034 0.001665  
6   1.8211 10-4 -2.1144 10-2 -8.5703 10-4 -1.3919 10-4   4.1349 10-8 0.2000 -0.1050 1.905
7 -6.0345 10-3 -7.3769 10-3 -6.8749 10-5   2.1217 10-4   7.2144 10-7 0.5537 -0.04433 12.49
8 -3.0634 10-3 -8.1268 10-3   1.5556 10-5 -3.2051 10-4   5.9883 10-7 0.3489 -0.09391 3.716
9 -2.2244 10-3 -9.5825 10-3 -2.5612 10-5 -4.8902 10-4   1.6476 10-6 0.1398 -0.1796 0.7786
10 -4.0244 10-2   2.7424 10-3 -1.4765 10-2 -9.1042 10-3   3.4555 10-8 0.6172 0.002074  
11 -4.4068 10-3 -6.9471 10-3   1.2804 10-4 -1.3335 10-4   3.4657 10-8 0.3564 -0.02496 14.28
12 -2.0087 10-3 -6.5533 10-3   3.9831 10-4 -3.1209 10-4   1.2262 10-7 0.1588 -0.03486 4.557
13 -5.3235 10-3 -3.1723 10-2 -9.7516 10-4 -1.9609 10-3   3.0551 10-6 0.5999 -0.1863 3.220
14 -3.1380 10-3 -2.4780 10-2 -3.9949 10-4 -2.5221 10-4   3.7029 10-6 0.3999 -0.2703 1.479
15 -2.0407 10-3 -2.6919 10-2 -1.6638 10-4 -2.3280 10-4   3.0496 10-6 0.1999 -0.4942 0.4044
16 -3.4327 10-2 -4.7343 10-2 -6.6893 10-2   1.2627 10-2   2.2690 10-8 0.6000 -0.03219 18.64
17   2.5961 10-3 -5.9489 10-2 -8.4285 10-3   6.1018 10-4   1.2778 10-7 0.4000 -0.1731 2.311
18   1.0951 10-3 -4.1757  10-2 -3.0408 10-3 -2.7699 10-3   9.7108 10-8 0.1999 -0.2049 0.9759

\end{figure} Figure 8: Comparison between the two methods for obtaining the time evolution of the eccentricities, for the experiments 1-18. The solid lines represent the results obtained using the first method while the dashed lines represent the fitted curves obtained from the second method
Open with DEXTER

The eccentricity fittings using the two methods are shown in Fig. 8 as straight lines. We can observe that the second method (dashed lines) agrees with the first method (solid lines), for example for the experiments 5, 8, but they disagree completely for the experiment 3 as shown in Fig. 8. The reason is that in self-consistent experiments, for closed contact pairs, all orbital parameters change faster and a more general procedure (second method) is necessary in order to consider the new features that appear.

As a final remark of this section we should mention that we have tested the hypothesis of the time evolution of the precession but we have not presented it. If we assume that the longitude of the pericenter also changes linearly with the time, we can show that the fitting is worst than fixing null the longitude of the pericenter. For example, for the experiment 9 we obtained $\chi^2=5.991~10^{-6}$ that is bigger than the previous value (Table 3).

6 The circularization time

Once we have fitted the curves $r=r(\alpha)$, we have gotten for each value of $\alpha$, a value of e and t, for each experiment. With these values of e and t we have fitted a linear function of t and obtaining directly $e(t)=a_{\rm e}+b_{\rm e} t$. The circularization time is calculated assuming $a_{\rm e}+b_{\rm e} t_{\rm c} = 0$, that gives us $t_{\rm c} = -a_{\rm e} / b_{\rm e}$. The fitted parameters are shown in Table 3, where the circularization times for the experiments 4, 5 and 10 are not shown because the slopes of the fitted eccentricities are positive. This is due to the very short orbit of these simulations. We notice that the circularization time is larger than the Hubble time, except some simulations with initial eccentricity e=0.2. Comparing the circularization times of the experiments 1-6 with 13-18, we can notice that the latter ones are approximately twice faster than the former ones (except the experiment 3). However, using the Chandrasekhar's formula we could expect that the latter models would evolve, not twice but, five times faster, since there has been a factor five increase in mass. The long circularization times obtained in this work might result from the fact that the galaxies are equal mass. Besides, the Chandrasekhar's formula is valid for a homogeneous medium, which in our case is very far from that condition. Using the pericentric distance in Table 2 and the fitted values of the eccentricities and the circularization time of the Table 3, we can fit an equation given by

$\displaystyle t_{\rm c}$ = $\displaystyle 10^{(-0.03 \pm 0.92)} r_{\rm p}^{(0.53 \pm 0.17)}$ (21)
    $\displaystyle \times e^{(1.56 \pm 0.10)}
{\left( {M_1 \over M_2} \right)}^{(0.64 \pm 0.05)},$ (22)

where $\chi^2=0.09$, the correlation coefficient is 0.72, $t_{\rm c}$ is in units of Hubble time and $r_{\rm p}$ in units of kpc. Comparing this equation to the Eq. (9) in Junqueira & de Freitas Pacheco (1994), we can see that the major difference comes from the pericenter distance exponent of $r_{\rm p}$, where in our case it is smaller than their one. Thus, we can conclude that our binary galaxies circularize faster, due to the fact that our halos are in contact and they are not static in time. In recent work, Colpi (1998), studying the accretion of a satellite onto a spherical galaxy using a linear response theory in order to modify and solve the equations of motion, has shown that the time of binary orbital decay has a power-law dependence on the mass ratio of the satellite ($M_{\rm s}$) and the primary ($M_{\rm p}$) given by $(M_{\rm p}/M_{\rm s})^{0.4}$. Although the methods be different and the satellite be a point mass, it is interesting to remark that our result is close to her one.

7 Conclusions

We have performed several self-consistent numerical simulations of galaxy pairs in order to study the orbital evolution of contact systems in closed orbits. We had not included deep interpenetrating encounters because in this case the galaxies would merge quickly (timescale shorter than a Hubble time). This is corroborated by the fact that we do not see in the universe so many close interacting binary galaxies in comparison with open pairs. The aim of this work was to study the open binary galaxies and to clarify the controversy of the orbital circularization and the eccentricity probability distribution of the population of the observed pairs of galaxies. In particular we investigated the timescale of the circularization when dark halos are overlapping. Our simulations have shown that the circularization time is larger than the Hubble time, except some experiments with initial eccentricity equal to 0.2. This result confirms partially the result of Junqueira & de Freitas Pacheco (1994) but the circularization rate of the orbits in the present work is faster, due to the dynamical friction. Because the time required to circularize many of the orbits we have discussed is larger than a Hubble time, it seems quite possible that some real pairs are still at the beginning stages of their orbital evolution. Bartlett & Charlton (1995) performing numerical simulations of close pairs of galaxies in initial parabolic orbits found that: (1) the galaxies trace out quasi-elliptical paths after the first closest approach; (2) there is little evidence for significant orbital circularization during the subsequent bounces. Our simulations confirm these results, although they have started their experiments with parabolic orbits that evolve to quasi-elliptical orbits. Chatterjee & Magalinsky (1997) pointed out that the observational data indicate the existence of very wide physical galaxy pairs. These pairs with the close ones form a single population and they are characterized by highly radial orbits just in the final stages. According to these authors the circularization occurs very late, and our results are in agreement with this conclusion, however we have analyzed systems in initially elliptical orbits that evolve to a more complex geometry.

In a recent work, van den Bosch et al. (1999) studied the effect of the dynamical friction in the distributions of orbital eccentricities in a variety of spherical potentials. They found a weak dependence of the decay times on eccentricities with absence of any significant amount of circularization and concluded that the dynamical friction does not lead to strong changes in the distribution of orbital eccentricities. They show that the timescale for dynamical friction is shorter for more eccentric orbits but that dependence is considerably weaker than claimed previously by Lacey & Cole (1993), based on an analytical integration of Chandrasekhar's formula. They also showed that, contrary to common belief, dynamical friction does not lead to circularization. Dynamical friction leads to only a very moderate change in the distribution of orbital eccentricities over time. These results agree with ours in the sense of having very long timescale of the circularization, although they have used a point mass satellite with a small fraction of the mass of the primary galaxy. We call attention that the van den Bosch's study was aimed at satellite rather than binary galaxies. In contrast, we have developed a formalism for characterizing the eccentricity of orbits within extended mass distributions.

As we have seen, our orbits are not well behaved elliptical ones as in Junqueira & de Freitas Pacheco (1994). The time evolution of the orbits of realistic galaxies is very important to determine the mass of binary galaxies. In order to calculate the total binary mass statistically we need to know the probability of finding the galaxies at a given phase angle, which is not so obvious to find because of the kind of orbits we have found in this work. Considering that our simulations do not present a single maximum central luminosity density at the end of the experiments, a sign that the merger of the galaxies has occurred (Bartlett & Charlton 1995), we can say that the timescale of this phenomenon is greater than the studied interval of time.

A final remark that we can do is about the internal structure of the galaxies. We have obtained that the de Vaucouleurs' density power law, at least until the radius $r < 5 r_{\rm e}$, does not change with the interaction with the other galaxy in a Hubble time interval.

One of the authors (S. J.) acknowledges the financial support from the Brazilian FAPERJ (No. E-26/150.245/99). We are also grateful to the generous amounts of computer time provided by NACAD-COPPE/UFRJ (CRAY J90) and by LNCC/CNPq (SP2 IBM RISC6000). The authors would like to thank Dr. Lars Hernquist for several suggestions that improved this work and Dr. Shashi Kanbur for reading the manuscript.


Copyright ESO 2001