A&A 366, 418-427 (2001)
DOI: 10.1051/0004-6361:20000256
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
Abstract
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
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
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.
Our simulations assume two possible different primary and secondary galaxies with N=6000 and equal-mass particles, total mass of 2 10^{11} and 10^{12} , whose parameters are shown in Table 1. Hereinafter the value of the gravitational constant is G = 1 and the Hubble constant is H_{0} = 100 kms^{-1}/Mpc. Our units are [L]=1 Mpc, [V]=100 kms^{-1}, [T]=h^{-1}H_{0}^{-1} = 9.778 10^{9} years and [M]= 2.325 10^{12} , where h=H_{0}/100. In these units, the galaxy masses and are 0.086 and 0.430 mass units, respectively.
Type | |||||||||||
g_{1} | 0.0286 | 2000 | 0.0013 | 0.0094 | 0.0574 | 4000 | 0.0026 | 0.0188 | 0.0860 | 6000 | |
g_{2} | 0.1433 | 10000 | 0.0022 | 0.0160 | 0.2867 | 20000 | 0.0044 | 0.0320 | 0.4300 | 30000 |
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 kpc, for the most massive galaxy. We have also truncated its dark halo at 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,
We have also assumed that the mass-luminosity ratio of the galaxies is constant and equal to . 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)
The density of dark halo particles is assumed to be isothermal (van Albada et al. 1985)
and is given by
The Jeans' equation for a spherical system with isotropic velocity distribution
(Binney & Tremaine 1987, p. 204) is given by
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 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 (the softening parameter) from the initial system center, with a time step size and . 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 , thus the luminous cores of the galaxies have being followed accurately.
The Tree-Code, as parameterized requires 1000 Cray J90 CPU hours for 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.
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
The initial tangential velocity is given by
No. | prim. gal. | sec. gal. | e | ||||
1 | g_{1} | g_{1} | 2 | 0.0752 | 0.6 | 0.3008 | 0.4783 |
2 | g_{1} | g_{1} | 2 | 0.0752 | 0.4 | 0.1755 | 0.7669 |
3 | g_{1} | g_{1} | 2 | 0.0752 | 0.2 | 0.1128 | 1.1045 |
4 | g_{1} | g_{1} | 4 | 0.1504 | 0.6 | 0.6016 | 0.3382 |
5 | g_{1} | g_{1} | 4 | 0.1504 | 0.4 | 0.3509 | 0.5423 |
6 | g_{1} | g_{1} | 4 | 0.1504 | 0.2 | 0.2256 | 0.7810 |
7 | g_{1} | g_{2} | 2 | 0.1016 | 0.6 | 0.4064 | 0.7127 |
8 | g_{1} | g_{2} | 2 | 0.1016 | 0.4 | 0.2371 | 1.1428 |
9 | g_{1} | g_{2} | 2 | 0.1016 | 0.2 | 0.1524 | 1.6458 |
10 | g_{1} | g_{2} | 4 | 0.2032 | 0.6 | 0.8128 | 0.5039 |
11 | g_{1} | g_{2} | 4 | 0.2032 | 0.4 | 0.4741 | 0.8081 |
12 | g_{1} | g_{2} | 4 | 0.2032 | 0.2 | 0.3048 | 1.1638 |
13 | g_{2} | g_{2} | 2 | 0.1280 | 0.6 | 0.5120 | 0.8197 |
14 | g_{2} | g_{2} | 2 | 0.1280 | 0.4 | 0.2987 | 1.3144 |
15 | g_{2} | g_{2} | 2 | 0.1280 | 0.2 | 0.1920 | 1.8930 |
16 | g_{2} | g_{2} | 4 | 0.2560 | 0.6 | 1.0240 | 0.5796 |
17 | g_{2} | g_{2} | 4 | 0.2560 | 0.4 | 0.5973 | 0.9294 |
18 | g_{2} | g_{2} | 4 | 0.2560 | 0.2 | 0.3840 | 1.3385 |
is the number which multiplies the sum of the halo radius
of the primary (
)
and the secondary galaxy (
), giving the
pericenter distance. Thus,
.
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.
Figure 1: Snapshots at the last pericenter time for differentexperiments | |
Open with DEXTER |
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 or , 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.
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
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 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.
The first method (Junqueira & de Freitas Pacheco 1994) is based in the
equation of an ellipse in polar coordinates, given by
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
The second method consists in assuming that the
major semi-axis a and the eccentricity e vary linearly with the angle
in Eq. (15). Another hypothesis is that
the focus of the ellipse change also linearly on the orbit plane with
.
Thus,
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
and
.
Finally, in Fig. 5d we present the orbit where all
the three quantities, i.e., a, e and
evolve in time with the same orbital
parameter rate.
Figure 5: Theoretical curves as described in the text. The synthetic orbits begin always at the position and | |
Open with DEXTER |
Figure 6: The eccentricities calculated for the synthetic curves shown in Fig. 5, considering a time rate of the orbital parameters, a, e and , 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 , 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.
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 but these cases are the simulations whose orbits have not completed a period.
No. | ||||||||
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 |
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 that is bigger than the previous value (Table 3).
Once we have fitted the curves
,
we have gotten for each value of
,
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
.
The circularization time is calculated assuming
,
that gives
us
.
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
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 , does not change with the interaction with the other galaxy in a Hubble time interval.
Acknowledgements
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.