text:
# Charge profile of surface doped C60

11institutetext: Theoretische Physik, Eidgenössische Technische Hochschule, CH-8093 Zürich, Switzerland 22institutetext: Laboratoire de Physique Quantique & UMR–CNRS 5626, Université Paul Sabatier, F-31062 Toulouse, France

## Abstract

We study the charge profile of a C60-FET (field effect transistor) as used in the experiments of Schön, Kloc and Batlogg. Using a tight-binding model, we calculate the charge profile treating the Coulomb interaction in a mean-field approximation. At low doping, the charge profile behaves similarly to the case of a continuous space-charge layer and becomes confined to a single interface layer for doping higher than $`0.3`$ electron (or hole) per C60 molecule. The morahedral disorder of the C60 molecules smoothens the structure in the density of states.

Surface doping of high quality organic crystals has recently been achieved by Schön, Kloc and Batlogg art:schoen0 ; art:schoen0b ; art:schoen0c . This breakthrough has led to novel devices, new superconductors and measurements of transport properties over unprecedented wide ranges of carrier concentration and temperature. A tunable field effect transistor (FET) is used to inject carriers near to the surface of the organic crystal. In this paper we report on calculations of the resulting charge profile. We chose to investigate the case of a C60-crystal which has a special interest in view of the recent observations of superconductivity with a high $`T_c`$ ($`52`$ K) in hole doped samples.

The FET used in the work of Schön et al. is illustrated in Fig. 1(a). In the present work, we consider a C60 crystal (fcc lattice) with a  plane parallel to the gate. Undoped C60 is a semiconductor with a 2 eV gap. However, when an electric potential is applied between the gate electrode and the source/drain electrodes, either electrons or holes accumulate on the interface between the C60 crystal and the gate dielectric. This leads to a doping of the C60 interface planes which allows current to flow parallel to the interface, between the source and the drain. The resistance of this channel can be measured as a function of temperature and doping. The FET can be represented as a planar capacitance with equal and opposite charges on the metal gate and the C60 planes near to the interface as shown in Fig. 1(b). We will calculate the charge profile as a function of the total charge induced at the interface (or equivalently on the metal gate) rather then in terms of the potential applied to the gate. We limit ourselves to the case of an ideal planar interface with no steps or imperfections.

We begin by introducing a tight-binding model to describe the electronic structure of the C60 crystal. The effects of Coulomb interactions are treated within a mean-field (or Hartree) approximation. The resulting density profile gives rise to highly structured density of states (DOS). However we show that the inclusion of morahedral (or orientational) disorder leads to a broadening and suppression of this structure in the DOS.

The C60 molecule, which has icosahedral symmetry, is almost spherical. Therefore, in the solid, it orders naturally in the close packed fcc structure. In the following we will assume the (hypothetical) unidirectional structure where all C60 molecules have the same orientation art:satpathy92 . The space group for this structure is $`Fm\overline{3}m`$ and its primitive cell contains one molecule. The actual low-temperature structure with space group $`Pa\overline{3}`$ is more complicated and has a unit cell with 4 molecules. However, in the present work, we are interested in the overall charge filling of the layers, for which details of the band structure are unimportant. The lattice constant of the cubic cell (containing two molecules) is $`a=14`$ Å for pure C60 art:liu91 . We introduce here the unit vectors $`\stackrel{}{e}_x`$, $`\stackrel{}{e}_y`$ and $`\stackrel{}{e}_z`$ which span the cubic cell. If one considers only a  layer, then the C60 molecules form a 2D square lattice. In this case the primitive cell has a side length of $`b=\frac{a}{\sqrt{2}}=10`$Å which is the distance between two neighboring molecules. The 2D primitive cell is spanned by the unit vectors $`\stackrel{}{e}_1=\frac{1}{\sqrt{2}}(\stackrel{}{e}_x+\stackrel{}{e}_y)`$ and $`\stackrel{}{e}_2=\frac{1}{\sqrt{2}}(\stackrel{}{e}_x+\stackrel{}{e}_y)`$.

The kinetic energy is well described by a tight-binding Hamiltonian with nearest neighbor (n.n.) hopping art:satpathy92 ; art:laouini95 :

$$H_{\text{TB}}=\underset{nj\alpha ,n^{}j^{}\alpha ^{}}{}t_{\alpha ^{}\alpha }(\stackrel{}{\delta })c_{n^{}j^{}\alpha ^{}}^{}c_{nj\alpha },$$
(1)
where $`t_{\alpha ^{}\alpha }(\stackrel{}{\delta })`$ are the hopping integrals depending on the relative position $`\stackrel{}{\delta }`$ of neighboring molecules. The label $`n`$ denotes the different layers starting with the $`n=1`$ layer closest to the gate. Within a layer $`n`$, the different sites are labelled by the index $`j`$. Finally, the index $`\alpha `$ denotes the different orbitals of the C60 molecule. The brackets $``$ indicate a summation over n.n. only. Summation over spin degrees of freedom is implicitly assumed. Hamiltonian (1) is invariant under translations parallel to the $`z=0`$ plane and can partly be diagonalized in the 2D k-space:

$$H_{\text{TB}}=\underset{\stackrel{}{k}}{}\underset{n\alpha ,n^{}\alpha ^{}}{}H_{n^{}\alpha ^{},n\alpha }(\stackrel{}{k})c_{n^{}\alpha ^{}}^{}(\stackrel{}{k})c_{n\alpha }(\stackrel{}{k}),$$
(2)
where $`|nn^{}|1`$ and $`\stackrel{}{k}`$ is the two-dimensional wave-vector parallel to the planes. In the case of the conduction band, the structure of the hopping integrals $`t_{\alpha ^{}\alpha }(\stackrel{}{\delta })`$ are given by the t1u symmetry of the LUMO (Lowest Unoccupied Molecular Orbital) which is threefold degenerate art:laouini95 . The corresponding wave-functions can be chosen such as to transform as x, y and z under the icosahedral symmetry group which reduces the number of independent hopping integrals to $`4`$. They are given in Table 1. The matrix elements $`H_{n^{}\alpha ^{},n\alpha }(\stackrel{}{k})`$ in (2) can be calculated explicitly which yields for intraplanar processes ($`n^{}=n`$),

$`H_{nx,nx}(\stackrel{}{k})`$ $`=`$ $`2t_{xx}[\mathrm{cos}(k_1b)+\mathrm{cos}(k_2b)],`$
$`H_{ny,ny}(\stackrel{}{k})`$ $`=`$ $`2t_{yy}[\mathrm{cos}(k_1b)+\mathrm{cos}(k_2b)],`$
$`H_{nx,ny}(\stackrel{}{k})`$ $`=`$ $`2t_{xy}[\mathrm{cos}(k_1b)\mathrm{cos}(k_2b)],`$ (3)
$`H_{nz,nz}(\stackrel{}{k})`$ $`=`$ $`2t_{zz}[\mathrm{cos}(k_1b)+\mathrm{cos}(k_2b)],`$
and for interplanar processes ($`n^{}=n+1`$),

$`H_{n+1x,nx}(\stackrel{}{k})`$ $`=`$ $`2t_{yy}\mathrm{cos}[(k_1k_2)b/2]+2t_{zz}\mathrm{cos}[(k_1+k_2)b/2],`$
$`H_{n+1y,ny}(\stackrel{}{k})`$ $`=`$ $`2t_{zz}\mathrm{cos}[(k_1k_2)b/2]+2t_{xx}\mathrm{cos}[(k_1+k_2)b/2],`$
$`H_{n+1z,nz}(\stackrel{}{k})`$ $`=`$ $`2t_{xx}\mathrm{cos}[(k_1k_2)b/2]+2t_{yy}\mathrm{cos}[(k_1+k_2)b/2],`$
$`H_{n+1x,nz}(\stackrel{}{k})`$ $`=`$ $`i\mathrm{\hspace{0.17em}2}t_{xy}\mathrm{sin}[(k_1k_2)b/2],`$ (4)
$`H_{n+1y,nz}(\stackrel{}{k})`$ $`=`$ $`i\mathrm{\hspace{0.17em}2}t_{xy}\mathrm{sin}[(k_1+k_2)b/2].`$
The band structure and DOS of a single layer (interplanar processes are turned off) are shown in Fig. 2. As a comparison, the 3D DOS is shown as well. The logarithmic van-Hove singularities are clearly apparent in the DOS. Moreover, one observes a symmetric DOS. This is explained by the substitution $`\stackrel{}{k}\stackrel{}{k}+(\frac{\pi }{b},\frac{\pi }{b})`$ which leads to $`HH`$ for the one-layer Hamiltonian (Charge profile of surface doped C60). It implies also that all bands at $`(\frac{\pi }{2b},\frac{\pi }{2b})`$ cross at zero energy. Coupling to neighboring layers breaks this symmetry. The structure of the hopping integrals for the valence band are somewhat more complicated due to the fivefold degeneracy of the HOMO (Highest Occupied Molecular Orbital) with hu symmetry. They are described in detail in Ref. art:laouini95  and similar matrix elements as in (Charge profile of surface doped C60) and (Charge profile of surface doped C60) were calculated with parameters taken from art:laouini95 .

The Coulomb interaction can be included by adding the following term to the tight-binding Hamiltonian (1),

$`H_{\text{Coul}}`$ $`=`$ $`{\displaystyle \underset{nj\alpha }{}}U_{\text{Gate}}(n)n_{nj\alpha }`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{nj\alpha ,n^{}j^{}\alpha ^{}}{}}V_{ee}(nj,n^{}j^{})n_{nj\alpha }n_{n^{}j^{}\alpha ^{}},`$
where $`n_{nj\alpha }`$ is the number operator. The first term is the potential due to the positive charge on the gate. The second term is the electron-electron Coulomb repulsion. Here we follow Antropov et al. art:antropov92  by using a screened Coulomb interaction for electrons on different sites and an on-site interaction $`U_0`$ for electrons on the same site,

$$V_{ee}(nj,n^{}j^{})=\{\begin{array}{cc}U_0& \text{if}\stackrel{}{R}_{nj}=\stackrel{}{R}_{n^{}j^{}}\hfill \\ \frac{e^2}{\epsilon \left|\stackrel{}{R}_{nj}\stackrel{}{R}_{n^{}j^{}}\right|}& \text{if}\stackrel{}{R}_{nj}\stackrel{}{R}_{n^{}j^{}}\hfill \end{array},$$
(6)
with $`\epsilon `$ being the dielectric constant of C60 and $`\stackrel{}{R}_{nj}`$ denoting the positions of the sites. Note that this interaction does not distinguish between different orbitals. The dielectric constant $`\epsilon `$ and the on-site interaction $`U_0`$ are easily accessible in the literature. We use the experimental estimation of $`\epsilon =4.4`$ by Hebard et al. art:hebard91b  and the theoretical value of $`U_0=1`$ eV by Antropov et al. art:antropov92 . Somewhat different values of $`\epsilon =3.66`$ and $`U_0=1.27`$ eV are proposed from theoretical calculations by Pederson et al. art:pederson92 .

The full Hamiltonian including the Coulomb term is a complicated many-body problem. We solve it within the mean-field approximation. From symmetry, solutions should be invariant under translations parallel to the planes. The mean-field charge distribution and the resulting potential are then independent of the in-plane site position $`j`$ and the wave-functions can be decomposed as

$$\mathrm{\Psi }_{\nu \stackrel{}{k}}(nj\alpha )=\frac{1}{\sqrt{N_{\mathrm{}}}}e^{i\stackrel{}{R}_{nj}^{\mathrm{}}\stackrel{}{k}}\psi _{\nu \stackrel{}{k}}(n\alpha ),$$
(7)
where $`\nu `$ is a band index and $`\stackrel{}{R}_{nj}^{\mathrm{}}`$ is the inplane component of $`\stackrel{}{R}_{nj}`$. $`N_{\mathrm{}}`$ is the number of sites per layer. For a fixed $`\stackrel{}{k}`$, the perpendicular part $`\psi _{\nu \stackrel{}{k}}(n\alpha )`$ of the wave-function is the solution of a Schrödinger equation with a partly discrete spectrum {$`E_{\nu \stackrel{}{k}}`$}. In the mean-field approximation, every occupied electronic $`\nu \stackrel{}{k}`$ state contributes $`\left|\psi _{\nu \stackrel{}{k}}(n\alpha )\right|^2/N_{\mathrm{}}`$ electrons to each orbital $`\alpha `$ in layer $`n`$. The sum over all these contributions leads, at zero temperature, to the electron distribution

$$\rho (n)=\frac{2}{N_{\mathrm{}}}\underset{E_{\nu \stackrel{}{k}}<E_F}{}\underset{\alpha }{}\left|\psi _{\nu \stackrel{}{k}}(n\alpha )\right|^2,$$
(8)
where $`\rho (n)`$ is the total number of electrons per site in layer $`n`$. The factor 2 comes from the spin degrees of freedom. Within the mean-field approximation, the Coulomb part (Charge profile of surface doped C60) can be written as $`H_{\text{Coul}}^{\text{MF}}=_{nj\alpha }U_{\text{MF}}(n)n_{nj\alpha }`$ where

$$U_{\text{MF}}(n)=U_{\text{Gate}}(n)+\underset{n^{}j^{}}{}\rho (n^{})V_{ee}(nj,n^{}j^{}).$$
(9)
The sum in the second term as well as the potential created by the gate are divergent. However, these divergences cancel when the overall system is neutral. In order to make this more explicit, we add and subtract a hypothetical potential created by uniformly charged layers whose total charge equals the actual charge of the C60 layers. The mean-field potential can then be rewritten as a sum of two terms

$$U_{\text{MF}}(n)=U_{\text{Cap}}(n)+U_{\text{Corr}}(n).$$
(10)
The first term, $`U_{\text{Cap}}`$, is given by the superposition of a series of planar capacitance formed by the negative charge $`e\rho (n^{})`$ on each layer $`n^{}`$ combined with the corresponding positive charge $`+e\rho (n^{})`$ on the gate. Choosing $`U_{\text{Cap}}(1)=0`$ yields

$$U_{\text{Cap}}(n)=\frac{4\pi e^2}{\epsilon \sqrt{2}b}\underset{n^{}}{}\rho (n^{})[\mathrm{min}(n,n^{})1],$$
(11)
where it has to be distinguished whether the given layer n lies inside or outside the capacitance with charges $`\pm e\rho (n^{})`$. Using $`\epsilon =4.4`$ and $`b=10`$ Å leads to a prefactor $`4\pi e^2/\epsilon \sqrt{2}b=2.9`$ eV. The second part of the mean-field potential, $`U_{\text{Corr}}`$, includes the corrections arising from the discreteness of the charge distribution and is given by

$`U_{\text{Corr}}(n)`$ $`=`$ $`{\displaystyle \underset{n^{}}{}}\rho (n^{})V_{\text{Corr}}(|nn^{}|)\text{with}`$ (12)
$`V_{\text{Corr}}(|nn^{}|)`$ $`=`$ $`{\displaystyle \underset{j^{}}{}}V_{ee}(nj,n^{}j^{}){\displaystyle \frac{d^2R_n^{}e^2}{b^2\epsilon \left|\stackrel{}{R}_n\stackrel{}{R}_n^{}\right|}},`$
where $`b^2`$ is the surface of the two-dimensional unit cell and $`\stackrel{}{R}_n`$ lies in the plane $`n`$. In equation (12) we expect the largest corrections for $`n^{}=n`$ since the difference of the potential created by a uniformly and a discretely charged plane is only sizeable at short distance. $`V_{\text{Corr}}`$ cannot be expressed analytically since it involves an infinite sum over a two-dimensional lattice. The result from numerical summation is given in Table 2. The function $`V_{\text{Corr}}(\mathrm{\Delta }n=|nn^{}|)`$ is rapidly decreasing and we will only retain the $`\mathrm{\Delta }n=0`$ term in the following. The correction term then becomes

$$U_{\text{Corr}}(n)=\rho (n)\left(3.9\frac{e^2}{\epsilon b}U_0\right).$$
(13)
In fact, the second part of this correction overestimates the on-site interaction $`U_0`$. The contribution of this term is proportional to $`\rho (n)`$ which reflects the fact that an electron feels its own mean-field. This is an artifact of the mean-field approximation. In order to improve our model, we require that the on-site interaction should only be effective when there is more than one electron on the molecule ($`\rho (n)>1`$), leading to the correction term

$$U_{\text{Corr}}(n)=3.9\frac{e^2}{\epsilon b}\rho (n)+U_0\mathrm{max}[0,\rho (n)1].$$
(14)
The mean-field Hamiltonian was solved self-consistently at zero temperature by numerical means. The mean-field potential (10) was used with the correction term (14). The resulting charge profile as a function of total charge is shown in Fig. 3. We have done the calculation for both, conduction and valence band, which led to similar results. In the following we focus on the conduction band. It can be seen that more than 98% of the total charge is confined to the first layer for doping higher than $`\rho _{\text{tot}}=0.3`$. Furthermore the confinement to the interface increases with the total charge. The mean distance of the charge distribution from the interface is defined as $`z_0=_{n1}nd\rho (n)/\rho _{\text{tot}}`$, where $`d=7`$Å is the width of one layer. For $`\rho _{\text{tot}}<0.1`$, this mean distance is found to follow a power law given by $`z_0[\text{Å}]3.8\rho _{\text{tot}}^{1/3}`$ with $`\rho _{\text{tot}}`$ in units of particles per area of the 2D unit cell. This power law behavior with respect to $`\rho _{\text{tot}}`$ is identical to the case of a standard space-charge layer in a continuous medium art:ando82 , which predicts $`z_0[\text{Å}]=3.02(\rho _{\text{tot}}m_z/\epsilon )^{1/3}=5.08\rho _{\text{tot}}^{1/3}`$, where the effective mass (perpendicular to the interface) of the bulk conduction band minima (X(1,0,0) points) is estimated to be $`m_z=0.92m_e`$ art:satpathy92 ; art:laouini95 . However, it is important to notice that our microscopic calculation gives a different prefactor. In addition, $`z_0`$ tends to saturate for $`\rho _{\text{tot}}>0.1`$, a regime that cannot be understood in the continuum model. In order to test the role of the correction term, we repeated the same calculation with $`U_{\text{MF}}=U_{\text{Cap}}`$, i.e. for uniformly charged planes. This led to a somewhat weaker confinement at high doping ($`10`$% less charge on the first layer for $`\rho _{\text{tot}}=1`$). At low doping the effect of $`U_{\text{Corr}}`$ vanishes because the charge is distributed over several layers. Knowing the charge profile and hence $`U_{\text{MF}}(n)`$, the band structure can be calculated. Fig. 4 shows the conduction band for $`\rho _{\text{tot}}=0.1`$. Few bands lie below the continuous spectrum and only the lowest of the 3 LUMO-bands is occupied in the region of the Q-point. In Fig. 5(a) the evolution of the discrete states at Q is shown as a function of doping. One observes that the Fermi level is always below the bottom of the continuous spectrum. This is true in general and reflects the fact that all electrons are bound to the interface since the overall system is neutral. Furthermore, only the lowest band is occupied at low doping which is again consistent with theory of space-charge layers art:ando82 . Of course, at doping higher than $`\rho _{\text{tot}}2`$ the second LUMO-band has to be filled as well. In addition, we note that the continuous spectrum shifts up sharply as the total charge increases. Therefore, at doping higher than $`\rho _{\text{tot}}0.3`$, the first layer becomes essentially decoupled from the continuous states of the subsequent layers. The energy levels at Q remain constant and the dispersion is given by the LUMO-bands of a single, isolated layer (Fig. 2(a)). The decoupling effect is due to the correction term (14) which drastically shifts the energy of the first layer with respect to the second one. If this term is omitted, as shown in Fig. 5(b), then the bottom of the continuous spectrum follows roughly the Fermi energy and the decoupling of the first layer is much less effective.

In A3C60, a scaling of the superconducting $`T_c`$ with increasing lattice constant and hence with the DOS is observed art:gunnarsson97 . However, the sharp peaks of the 2D DOS (Fig. 2(b)) are not reflected in the behavior of $`T_c`$ as measured by Schön et al. art:schoen0 . A possible explanation for this is the morahedral disorder which is considered in the following for the case of a single 2D layer.

The rotational motion of a C60 molecule in a crystal is described by a potential in the three Euler angle coordinates. This potential has, in addition to an absolute minimum, a local minimum which is only $`\mathrm{\Delta }=11`$ meV higher and corresponds to a non-equivalent orientation of the molecule art:david93 . Therefore, C60 molecules can flip into the second orientation causing the so-called morahedral disorder. The ratio between the number of flipped and unflipped molecules is given by a Boltzmann factor $`\mathrm{exp}(\mathrm{\Delta }/kT)`$. However, a freezing in is observed at 85 K. Below this temperature the crystal is in a glass phase with 15% of the molecules being flipped. If the problem is considered from a tight-binding point of view, then one expects the hopping integrals to be different for hopping between inequivalently oriented molecules. Hence, one gets a tight-binding model where different hopping integrals are distributed randomly. This will suppress structure in the DOS, as was calculated by Gelfand and Lu in the case of bulk C60 art:gelfand92b . Here, we consider the effect of morahedral disorder on the DOS of a single electron-doped  layer. The two types of inequivalent orientations lead to two sets of LUMO wave-functions. We shall denote the new ones on flipped molecules as $`|\stackrel{~}{x}`$, $`|\stackrel{~}{y}`$ and $`|\stackrel{~}{z}`$. These two sets can be related by a 90-degree rotation around the z-axes. The hopping integrals between inequivalently oriented molecules are given in Table 3. In order to calculate the DOS of the disordered system we set up a finite, two dimensional system which could be diagonalized exactly (for non-interacting electrons). The resulting DOS is shown in Fig. 6. The van-Hove singularities are washed out, but a rather sharp peak remains at the center of the band.

In conclusion, we have investigated the electronic properties of surface doped C60. The charge profile was calculated as a function of total charge in the system and it was found that the charge accumulation on the first layer increases with the total charge. At high doping, above $`0.3`$ electron per molecule, the first layer becomes essentially decoupled from the subsequent layers. This suggests that the electronic system should be well described by a single layer in this regime. In particular, the relevant DOS is then the one of a two dimensional system and hence substantially higher than in the bulk. Having chosen the unidirectional structure (one orientation per primitive cell), we found maximal values of the DOS near the center of the band. Introducing morahedral disorder led to a smoothened DOS which still shows a peak at half filling, i.e. in the region where Schön et al. measure the maximum $`T_c`$.

###### Acknowledgements.

We thank B. Batlogg for fruitful discussions throughout this work. D. P. also acknowledges support from the Center for Theoretical Studies and from the Physics Department at ETH-Zürich.


format: text

text:
# Morphology–transitions at vicinal Cu–surfaces based on entropic step–step interaction and diffusion along steps

## Abstract

An extension of the Burton–Cabrera–Frank model \[Phil. Trans. R. Soc. London, Ser. A 243, 299 (1951)\] including diffusion along steps and entropic step–step interaction is introduced. This extended model is successfully applied to simulate experiments done in the group of H.–J. Ernst at CEA Saclay, France, for vicinal Cu surfaces. In particular, the rise of two qualitatively different morphologies can be explained by the competition of growth directions implied by the four distinct driving and restoring forces of the model. Each of those forces sets a length scale for the dominant mode of a developing step instability. This explains that the wavelengthes of meandering–instabilities observed at vicinal surfaces experimentally are usually larger than the one predicted by Bales and Zangwill \[Phys. Rev. B 41, 5500 (1990)\] theoretically.

During molecular beam epitaxy (MBE) appropriate conditions for the controlled growth of vicinal surfaces can be fine–tuned. This way one tries to fabricate either atomistically flat or nanostructured surfaces. The ability to exert control on structuring along the growth direction and fabricate substrates with smallest–scale–built–in periodicies normal to the surface is well advanced. Efforts to still enhance abilities to determine functional properties of a grown substrate now focus on the lateral structuring within one layer of growth. One direction along this line is to make use of the inherent instabilities due to the dynamics of the growth process itself. A necessary first step is to understand the basic wavelengthes of those inherent instabilities just as well as the kind of morphologies which will develop.

More then ten years ago, Bales and Zangwill predicted that a growing vicinal surface should undergo a step meandering instability when kinetic step edge barriers suppress the attachments of atoms to decending steps. According to their analysis, a straight step is linearly unstable against perturbations with wavelengths larger then $`\lambda _c=2\pi (2\mathrm{\Gamma }L_\mathrm{\Delta })`$ and a fastest growing wavelength at $`\lambda _u=\sqrt{2}\lambda _c`$ ($`\mathrm{\Gamma }=\frac{\mathrm{\Omega }\gamma }{k_BT}`$, $`L_\mathrm{\Delta }=\frac{Dc_{eq}}{Fl^2}`$). Here $`F`$ denotes the deposition rate, $`D`$ the diffusion constant for diffusion on the terraces and $`c_{eq}`$ the equilibrium concentration, $`\gamma `$ the step stiffness, $`\mathrm{\Omega }`$ the atomic area and $`k_BT`$ the termal energy. Even though the meandering instability has meanwhile been observed in experiments and simulations , the quantitative prediction of the Bales–Zangwill analysis could not be recovered in many of those experimental findings. This point received much interest in view of the recent experiments by Maroutian et al. at CEA Saclay. Analytical efforts to resolve the disagreement between experimental measurements and theoretical prediction led to a precise study of the limit of desorptionless growth. For to derive a single evolution equation in the weakly nonlinear regime this case has the interesting feature of displaying a singularity in the spirit of multiscale expansion. As a result, meander–wavelengthes larger than the ones predicted by Bales and Zangwill can be explained, however, without reaching the order of experimental observations. Other efforts center around the investigation of an extra diffusion current along step edges. Such a current triggers an asymmetry in energy barriers for atoms attaching at kinks from different step edge directions, the so called Kink Ehrlich–Schwoebel effect (KESE). KESE can either destabilize or stabilize steps depending on whether the slope of the meander–instability is greater or less than one. In  calculations including a stabilizing KESE–current along steps were compared to the further thesis that islands might nucleate at step edges and thereby constitute a still different wavelength of instability. It appeared that the latter thesis fits best to the experimental wavelength–values.

New light was shed on these investigations by experiments which revealed that a morphology qualitatively different from the meandering one can develop on vicinal surface as well, with growth conditions being exactly the same as for the meandering–morphology except a different angle of miscut. It can be characterized by the absence of one global growth direction and thereby resembles the so called degenerated morphology. Thus the insufficiency of theoretical approaches to explain the difference between the wavelength of instability predicted by Bales and Zangwill and the one observed experimtentally was extended by their failure to explain this degenerated morphology. Note that the surface of the in–phase–meandering morphology displayed in the STM topograph in Fig. 1(b) was wrongly identified as Cu(1,1,17) and is in fact a Cu(0,2,24) surface. The true Cu(1,1,17) on the other hand is the one to display a degenerated morphology.

In this letter I will introduce a novel model which can explain the rise of two different morphologies as well as their basic wavelengthes including KESE currents plus entropic step–step repulsion in the classical Burton–Cabrera–Frank (BCF) model. It is the competition between the destabilizing and the stabilizing forces combined in this extended BCF model (EBCF model) which results into two different kinds of morphologies depending on the difference in magnitude of orthogonal driving forces. A careful analysis of this interaction also reveals a new basic wavelength of instability which is in the order of the experimental observations.

To get an impression of the basic ideas underlaying the EBCF let us start with the one sided BCF model, i.e. the model formulated by Burton, Cabrera and Frank in 1951 for the case of complete supression of attachment of adatoms to decending steps. This one–sided BCF model constitutes a moving boundary problem of a diffusion–relaxation equation for the dynamics of the adatoms on the terraces and two boundary conditions for the conservation of mass and the conservation of energy at the steps, respectively:

$`_tc=D^2c{\displaystyle \frac{1}{\tau }}c+F`$ (1)
$`c_{eq}=c_{eq}^0(1+\kappa \mathrm{\Omega }\gamma /k_BT)`$ (2)
$`v_n=D\mathrm{\Omega }{\displaystyle \frac{c}{n}}.`$ (3)
Here $`c`$ denotes the areal adatom density, $`c_{eq}^0`$ the equilibrium concentration for a straight step and $`\kappa `$ the step curvature. Equations (2) and (3) are to be evaluated at the front of each step. For vicinal Cu surfaces in the temperature range of the experiments under discussion a current involving the diffusion coefficient $`D_m`$ along the kinked steps is operative. Following the notation of Pierre–Louis et al. this diffusion along steps plus its anisotropy due to KESE can be included in the moving boundary problem above extending Eq.(3) by $`_xJ`$, with $`J=J_k+J_e`$ ($`J_k,J_e`$ as in Eq.(2)–(4) in the paper of Pierre–Louis et al.). Note that unlike the authors of  I consider $`\gamma `$ to be the step stiffness instead of the line tension as the precise measure of the step’s tendency to straighten and thus reduce its kink density. As a consequence it is taken into account that just as surface undergoing a faceting instability develops into the surface–tension–minimizing morphology, instabilities at step edges have ideal orientations which can minimize line tension as well. This is crucial for the model since this is exactly the point in which the Cu(1,1,17) and the Cu(0,2,24) vicinal surface differ. Both surfaces are made up of (001) terraces of monoatomic hight and mean terrace width $`l=21.7`$Å. However, for the Cu(1,1,17) surface straight step edges before the growth of instabilities are parallel to the $`<`$100$`>`$ direction, whereas they are parallel to $`<`$110$`>`$ in the case of Cu(0,2,24). As a consequence the direction of minimal step stiffness, which can be determined to be in $`<`$130$`>`$ via an embedded atom method (in analogy to  for surfaces), encloses an angle of appr. $`63.5^0`$ with the straight steps in case of Cu(0,2,24), whereas the angle is appr. $`71.5^0`$ in the case of Cu(1,1,17). Since growth proceeds in the direction of minimal step stiffness, these are the orientations for the initial formation of the instabilities. They trigger destabilizing $`J_k`$ currents of different magnitudes along the step edge forming in this direction. For to evaluate $`J_k`$ a quantitative value for the KESE–length $`L_s`$ (ref. , Eq.(3)) is required. The relevant kink energy barrier is the one directed in $`<`$$`\overline{1}10`$$`>`$, giving 0.518eV using an effective medium theory potential. The diffusion barrier for the jump of a single adatom along the step edge in this direction is 0.399eV, leading to a KESE–length $`L_s`$ of 201.2. The resulting $`j_k`$ current densities for Cu(0,2,24) versus Cu(1,1,17) differ by one order of magnitude: for Cu(0,2,24) $`j_k`$ takes a value of $`7.11810^7`$2s)-1, whereas at Cu(1,1,17) $`j_k=6.39610^8`$2s)-1. Despite their different magnitudes in both cases KESE currents are destabilizing, favoring unsaturating amplitude growth. Wavelengthes of instabilities turn out to be even less then in the BCF model. As a consequence taking into account KESE as only additional driving force in the BCF model is not well suited to explain the experimental work by T. Maroutian and coauthors. An obvious antagonist of the destabilizing KESE is the repulsion due to the succeeding step. Entropic just as well as elastic interactions have to be taken into account. Since the step interaction engergy $`A`$ considering an elastic or a dipole momentum is small compared to the entropic repulsion in the temperature ranges under discussion, it is sufficient to extend Eq.(2) by taking into account the surpression of step wandering:

$`c_{eq}=c_{eq}^0(1+{\displaystyle \frac{\kappa }{k_BT}}(\mathrm{\Omega }\gamma +{\displaystyle \frac{(\pi k_BT)^2}{6l^2\gamma }}))(\stackrel{~}{2})`$
$`v_n=D\mathrm{\Omega }{\displaystyle \frac{c}{n}}_xJ(\stackrel{~}{3})`$
The additional term $`(\pi k_BT)^2/6l^2\gamma `$ is the step interaction parameter with $`l`$ equal the width of the terraces. Equations $`(1)`$, $`(\stackrel{~}{2})`$ and $`(\stackrel{~}{3})`$ constitute the extended BCF model (EBCF). Its simulation with T=280K, F=3$`10^3`$ ML/s, l=21.7 Å, $`D_m=10^6`$cm2/s, $`c_{eq}=8.20810^6`$Å-1 and $`\gamma `$= 1.034$``$eV/Å yields the following two morphologies for the Cu(0,2,24) surface (Fig. 1) and the Cu(1,1,17) (Fig. 2), respectively.

The four components regulating growth of instabilities in the EBCF are

1. driving of growth via the gradient of the adatom diffusion field normal to the interface setting the length scale $`L_\mathrm{\Delta }`$ for a primary wavelength of instability
2. restore via step stiffness (corresponding length scale: $`\mathrm{\Gamma }`$)
3. driving of amplitude growth via KESE (length scale: $`L_s`$, ref.)
4. restore via entropic repulsion (length scale: $`L_S=\frac{6\gamma l^2}{k_BT\pi ^2}`$) .

In the situation depicted by Fig. 2 the magnitude of the KESE current derivative $`_xJ`$ is neglectable small. The remaining forces due entropic and step stiffness effects act in directions perpendicular to each other. The result of this competition of prependicular forces of equal strength is the degenerated morphology as observed in other contextes with an analogous competition of driving forces.

In contrary in Fig. 1 the KESE current does contribute to the evolution of surface morphology. The magnitude of the driving force due to this current is of an order 10 smaller then the magnitude of the adatom diffusion field gradient in direction of minimal step stiffness. The precise factor is evaluated via simulation and reads 11.129 for the simulation run done for Fig. 1 averaged over 500 timesteps. Due to the dominance of $`L_\mathrm{\Delta }`$ growth can proceed in the direction of minimal step stiffness.

For to understand the basic wavelength of the instability the two new length scales $`L_S`$ and $`L_s`$, set by entropic interaction as in Eq.($`\stackrel{~}{2}`$) and by an anisotrop diffusion current along the step as in Eq.($`\stackrel{~}{3}`$) respectively, have to be taken into account. Eq.(1),($`\stackrel{~}{2}`$) and ($`\stackrel{~}{3}`$) constitutes a system with a type I bifurcation. It’s dispersion relation can be evaluated in analogy to a dilute binary alloy undergoing directional solidification. For $`\mathrm{\Gamma }Ł_S`$ the resulting expression in leading order reads

$`\lambda _u=\sqrt{2}\pi ({\displaystyle \frac{L_\mathrm{\Delta }\mathrm{\Gamma }Ł_S}{4}})^{1/3}+2\sqrt{2}\pi {\displaystyle \frac{\mathrm{\Gamma }L_s}{3l_T}}.`$ (4)
It replaces Eq.(1) in . Expr.(4) is sufficient in the most interesting temperature range $`2.8`$ to $`3.8`$ x $`10^3K^1`$. It displays deviations from an Arrhenius–type behaviour, which increase with increasing temperature. Nevertheless it is in good agreement with the experimental data (Fig. 3). Simulations of the EBCF model support this expression as well.

The extended BCF explains two relevant new results for growth at vicinal surfaces, namely the experimentally observed deviation from the Bales–Zangwill instability and the rise of a degenerated morphology. These results are related to experimental findings in the group of H.–J. Ernst. The new basic wavelength predicted theoretically from the EBCF–model can be observed as initial perturbation of the step at Cu(0,2,24) just as well as Cu(1,1,17) surfaces. The qualitative difference between their morphologies, i.e. whether a surface displays an in–phase–meandering (Cu(0,2,24)) or a degenerated structure (Cu(1,1,17)) after deposition of a few monolayers, depends on the magnitude of the angle between the direction of destabilizing and restoring force. Time scales for the development from initial perturbation to the full degenerated morphology could not be compared with experiments due to the lack of experimental data. It seems an interesting question whether in situ transitions from one morphology to the other can be obtained via a change in the ratio of driving forces (e.g. by lowering temperature during the experiment). If this morphology transition displays features which resemble a true phase transition remains to be investigated.

The diffusion–relaxation Eq.(1) is solved on a quadratic grid. The interface is discretized separatly by curvilinear segments and respective interpolations from the interface to the grid. Details of the scheme and its parallelization will be given in a further paper.

I thank T. Maroutian and M. Rusanen for sending me a STM–picture of the true Cu (1,1,17) surface–morpology as well as a respective preprint prior to publication. Comments by H. Müller–Krumbhaar to focus my attention on a precise understanding of the basic wavelength of instability arising in the EBCF–model are gratefully acknowledged. Part of the computations were done on the T3E and the Origin 3800 of the URZ at Dresden Technical University.
Note–The two surfaces (Cu(1,1,17) and Cu(0,2,24)) have also been modeled in  on the basis of a kinetic Monte Carlo algorithm without taking into account entropic repulsion. The agreement with experimental findings remains restricted, nevertheless it is a very precise study of the different effects of KESE at the different surfaces.


text:
# Dynamical evolution of triplets of galaxies

## 1 Introduction

Galaxies in the universe tend to aggregate in different levels of structure, from binaries to large super-clusters. Systems of three galaxies, or triplets, constitute the smallest galaxy groups and are the boundary of $`N`$-body systems that cannot be modeled by analytical methods. On other hand, existing compact triplets of galaxies (Karachentsev 2000) are not very distinct from compact groups of galaxies (Hickson 1997) and their dynamical study may help in our understanding of the latter.

Karachentseva, Karachentsev & Shcherbanovskii \[Karachentseva et al. 1979\] were the first to compile a systematic catalog of triplets and to study some of their observational properties. At the present time, observational surveys in search of triplets and studies to determine their kinematical properties are being conducted by several groups (e.g., Karachentseva & Karachantsev 2000; Infante et al. 2000). On other hand, Heidt et al. \[Heidt et al. 1999\] have found that a strongly interacting triplet may be responsible for the production of a BL Lac object, and Iglesias-Páramo & Vílchez (1997) show an example of a probable triple interaction of galaxies in Hickson’s compact group HCG 95. All this observational research points toward the necessity of developing a better understanding of the dynamics of triplets.

Triplets maintain also a close relation with the 3-body problem of celestial mechanics and, hence, their study is also of theoretical interest (Valtonen & Mikkola 1991; Heinämäki et al. 1998; Muzzio, Wachlin & Carpintero 2000). A popular idea in the community (Valtonen & Flynn 2000) is that there are no stable states in a triplet, but as we will show later this is incorrect.

In general, the problem of triplets, or 3-galaxy problem, has not received much attention in contrast to the dynamics of compact groups of galaxies. Several dynamical studies of triplets have been carried out by different authors in order to explain some of their observational properties, but most of the up-to-date studies have used either a point-particle approach to model galaxy interactions (e.g., Chernin & Mikkola 1991) or included Chandrasekhar’s formula to model dynamical friction (Zheng et al. 1993).

In the 3-body simulations made by Chernin & Mikkola (1991) they found too many hierarchical structures (i.e. a close binary and a far away third body) that does not correspond to what is observed. Chernin et al. \[Chernin et al. 1994\] have suggested that dark matter may be the necessary ingredient to help 3-body simulations better reproduce these observations, based up on numerical experiments that show that the probability of binary formation is reduced if dark matter is present in the system \[Chernin et al. 1989\]. On other hand, Zheng et al. (1993) have shown that the excess of hierarchical structures may be reduced somewhat by the inclusion of Chandrasekhar’s dynamical friction formula in simulations of point-like galaxies; it is to note that Zheng et al. start their simulations with five ‘galaxies’.

The use of non self-consistent galaxies in previous studies lead to consider current results on the dynamics of triplets as a first approximation to the 3-galaxy problem. Effects like dynamical friction are not easy to model by an explicit-physics approach, which in any event requires fine tuning using self-consistent galaxies (Velázquez & White 1999). Therefore some doubts are cast on the results obtained when the self-gravity of galaxies is not taken into consideration. The natural following step in the modeling of triplets is to treat galaxies as self-gravitating systems, which will allow them to absorb orbital energy and angular momentum during an encounter.

An important quantity to determine in triplets, as well as in other systems of galaxies, is their gravitational mass. Although large amounts of dark matter are rather well established in groups and clusters of galaxies \[Bahcall 1999\], the situation in low-multiplicity systems is still a matter of some debate. For example, several authors have suggested that binary galaxies, even with the known uncertainties in their orbital parameters \[Binney & Tremaine 1987\], do not have much more dark matter than that probably associated with their individual haloes (e.g., Karachentsev 1990; Honma 1999).

On other hand, numerical simulations using point-particles indicate that the mass of a triplet obtained by the virial mass estimator is unreliable, with overestimates up to one or two orders of magnitude (e.g., Kiseleva & Chernin 1988; Chernin & Mikkola 1991). Dolgachev & Chernin \[Dolgachev & Chernin 1997\], carrying out simulations of wide-triplets (W-triplets, see $`\mathrm{\S }2`$) with free-fall initial conditions, conclude that these systems must hold a considerable amount of dark matter and that their masses are significantly underestimated if their non-stationary state is not taken into consideration. The assessment of typical methods used to determine mass, when using self-consistent galaxies and involving a non-stationary situation, may shed light towards understanding some of these discrepancies. It will be evident later that the self-gravity of galaxies introduces another degree of uncertainty in the determination of mass due to the transfer of orbital to internal energy during interactions.

In this contribution we perform a set of numerical experiments on the dynamics of triplets using a self-gravitating galaxy model. In particular, we will focus on global quantities that are directly comparable with observations in order to test our model. We divide this work as follows. In $`\mathrm{\S }2`$ we present an outline of some of the observational properties of compact and wide triplets, and our adopted set of definitions to estimate the dynamical quantities. In $`\mathrm{\S }3`$ the galaxy model and the initial conditions to be used in the simulations are presented. In $`\mathrm{\S }4`$ our results are presented and discussed. Here we consider the general aspects of the evolution of triplets, merging histories, geometrical properties, and the median values of their observationally related dynamical quantities. Finally, in $`\mathrm{\S }5`$ we summarize the main results of this work.

## 2 Observational setting

Karachentseva et al. \[Karachentseva et al. 1979\] have presented the first systematic study of the properties of triplets having a compact configuration in the sky (see also Karachentsev 2000; Karachentseva & Karachentsev 2000). According to the statistical criteria of Anosova \[Anosova 19871987\] to distinguish physical from optical triplets, which has become a standard in selecting physical triplets and is considered to be a strong selection criteria (Zheng et al. 1993), about 45 triplets out of 83 are probably physical; we will refer to the former as K-triplets or compact triplets.

Trofimov & Chernin \[Trofimov & Chernin 1995\] have analyzed the triplets in the catalogs of Maia, da Costa & Latham (1989) and Huchra & Geller (1982). For these triplets, they derive a median mean harmonic radius that is about an order of magnitude larger than that of K-triplets, thus calling them wide triplets (W-triplets). In Table 1 the median values of different dynamical quantities, both for K-triplets and W-triplets, are provided for future reference. The quantities written are: the mean harmonic radius, $`R_\mathrm{H}`$; the one-dimensional velocity dispersion, $`\sigma `$; the dimensionless crossing time, $`H_0\tau _\mathrm{c}`$, and the virial mass estimate $`M_\mathrm{v}`$. A Hubble constant of $`H_0=75`$ km s-1 Mpc-1 is assumed throughout this work.

The values for K-triplets in Table 1 were obtained using the information for $`\sigma `$ and $`R_\mathrm{H}`$ provided by Karachentsev et al. (1989), and for W-triplets’ those found in Trofimov & Chernin (1995) and Dolgachev & Chernin (1997). The differences of values in Table 1 with those found in the referred works arise from a difference in the definitions used for the dynamical parameters. The expressions used here, as well as that for the median mass estimator to be used later, are as follows (e.g., Nolthenius & White 1987; Heisler, Tremaine & Bahcall 1987):

$`{\displaystyle \frac{1}{R_\mathrm{H}}}`$ $`=`$ $`{\displaystyle \frac{2}{N(N1)}}{\displaystyle \frac{2}{\pi }}{\displaystyle \underset{i<j}{}}{\displaystyle \frac{1}{R_{ij}}},`$ (1)
$`\sigma ^2`$ $`=`$ $`{\displaystyle \frac{1}{N1}}{\displaystyle \underset{i}{}}(v_iv)^2,`$ (2)
$`\tau _\mathrm{c}`$ $`=`$ $`{\displaystyle \frac{2R_\mathrm{H}}{\sqrt{3}\sigma }},`$ (3)
where $`v=N^1v_i`$, $`N=3`$, and

$$M_\mathrm{v}=\frac{6}{G}R_\mathrm{H}\sigma ^2,M_{\mathrm{med}}=\frac{6.5}{G}\mathrm{MED}_{ij}[(v_iv_j)^2R_{ij}],$$
(4)
with $`v_i`$ being the line-of-sight velocity of a galaxy and $`R_{ij}`$ the projected pair-wise separation of two galaxies. In these formulae the bulk velocity of galaxies and their position vectors are assumed to be randomly oriented.

The expression for the virial mass yields the standard formula of the virial mass estimator, which leads to accurate results for the total mass of an equilibrium system made of point-like particles (e.g., Aceves & Perea 1999). The median mass estimator is used here due to its insensitivity to extreme values, and also to compare its performance against the virial mass estimator. In the expression for $`M_{\mathrm{med}}`$ a Michie’s equilibrium model was used to obtain the constant factor; although this factor is perhaps not true for triplets, we none the less use it as is normal practice and also to assess its accuracy in a triplet of self-consistent galaxies.

In Fig. 1, histograms for the different dynamical quantities are provided, normalized to the total number of triplets assumed physical in each of the above references. We note in passing the existence of overlapping regions between the K-triplets’ and W-triplets’ distributions.

In addition to the previous kinematical quantities, one may consider also the shape of the triangle formed by the galaxies in the sky to characterize a triplet; an idea introduced by Agekyan & Anosova (1968) in terms of an homologous map (AA-map). In this map one normalizes to unity the largest side of the triangle formed by the galaxies in the sky, and the remaining galaxy determines, when a new reference system’s axis is set along the longest side, a new set of coordinates which are to be plotted in such diagram. By symmetry arguments all triangles can be represented inside a particular region; see Fig. 2. In this AA-map hierarchical (H) configurations lie at the lower extreme right, Lagrangian (L) triangles at the upper corner, and alignments (A) or chain-like configurations lie at the bottom of the diagram.

Karachentsev’s triplets are scattered more or less randomly in the AA-map, while W-triplets show an absence of L configurations. An important result found by Chernin et al. \[Chernin et al. 1994\] was that projection effects are not responsible for the lack of excess of hierarchical triplets. This indicates that K-triplets are absent of many close pairs. These authors also indicate that W-triplets show only a weak excess of number density in the H-area of the AA-map.

## 3 Models for triplet evolution

### 3.1 Galaxy model

In this work the focus is on the overall dynamics of triplets. Hence we consider it convenient to use a spherically symmetric $`N`$-body system to model a galaxy. In particular, we have adopted a Plummer sphere. We consider that this galaxy model can represent an elliptical or, to a first approximation, the dark halo of a spiral, and can be used to obtain some information about the merging activity of triplets. We have not made any attempt to make an explicit difference between a dark and luminous component.

The mass enclosed within a radius $`r`$ and the phase-space density distribution function $`f()`$ for the Plummer model are given, respectively, by:

$$M(r)=\frac{M(r/R_0)^3}{[1+(r/R_0)^2]^{3/2}},f()=\frac{24\sqrt{2}R_0^2}{7\pi ^3G^5M^4}||^{7/2},$$
(5)
where $`M`$ is the total mass, $`R_0`$ its scale-length, and $``$ the energy per unit mass (Aarseth, Hénon & Wielen 1974). We constructed each numerical galaxy with $`N=3000`$ particles, with positions and velocities randomly sampled by an acceptance-rejection technique using the relations in (5). The $`N`$-body units used here are such that $`G=M=R_0=1`$, thus the corresponding unit of time that follows is $`t=(R_0^3/GM)^{1/2}`$; the total mass of a triplet is $`M_\mathrm{t}=3`$. We chose the dynamical time for the galaxy model to be $`t_\mathrm{d}=2\sqrt{2}(R_0^3/GM)^{1/2}`$; i.e., $`\sqrt{8}`$ in our units. The half-mass radius is $`R_\mathrm{h}1.3`$. Given that the mass distribution of a Plummer sphere falls very rapidly, $`\rho r^5`$, and about $`99`$ percent of the enclosed mass is contained within $`r10R_0`$ we take this value as the extent of our galaxy model.

In order to compare the results of the $`N`$-body simulations with those of astronomical systems we proceeded as follows. We denote by a subscript ‘$`n`$$`N`$-body values and with an ‘$`a`$’ the ones in astronomical units. We chose the astronomical unit of mass to be $`\mu `$ times the approximate mass of the disc of the Galaxy and $`l`$ times its disc scale length, respectively:

$$m_\mathrm{a}=[\mathrm{\hspace{0.17em}5.1}\times 10^{10}\mu \mathrm{M}_{}]m_\mathrm{n},r_\mathrm{a}=[\mathrm{\hspace{0.17em}4.5}l\mathrm{kpc}]r_\mathrm{n}.$$
(6)
These values imply the following time and velocity transformation formulae, respectively:

$`t_\mathrm{a}`$ $`=`$ $`[\mathrm{\hspace{0.33em}20}l^{3/2}\mu ^{1/2}\mathrm{Myr}]t_\mathrm{n},`$
$`v_\mathrm{a}`$ $`=`$ $`[\mathrm{\hspace{0.33em}221}l^{1/2}\mu ^{1/2}\mathrm{km}\mathrm{s}^1]v_\mathrm{n}.`$ (7)
If we assume each Plummer sphere represents the Galaxy, with a halo extent of $`135`$ kpc and a mass of $`5.5\times 10^{11}`$ M, we have $`\mu =10.78`$ and $`l=3.01`$ as the scaling factors. In particular, these values correspond to Model B of Kuijken & Dubinski \[Kuijken & Dubinski 1995\] for the Milky Way that has a mass ratio of dark to luminous matter $`10:1`$. This model will constitute our fiducial system to scale our $`N`$-body results to astronomical ones, unless otherwise stated. Therefore, for our fiducial system, we have: $`r_\mathrm{a}13.5r_\mathrm{n}\mathrm{kpc}`$, $`m_\mathrm{a}5.5\times 10^{11}m_\mathrm{n}\mathrm{M}_{}`$, $`t_\mathrm{a}32t_\mathrm{n}\mathrm{Myr}`$, and $`v_\mathrm{a}419v_\mathrm{n}\mathrm{km}\mathrm{s}^1`$.

We followed in isolation the numerical galaxy using a treecode (Barnes & Hernquist 1986) to test its stability. After experimenting with different values, and adopting a set of computational compromises, we decided on the following numerical parameters for our simulations, both for galaxies in isolation and during the triplet simulations: tolerance parameter $`\theta =0.9`$, softening parameter $`\epsilon =0.07R_\mathrm{h}/20`$, time step $`\mathrm{\Delta }t=0.1`$, and quadrupole corrections to the potential were included. With these parameters energy was conserved to within 0.3 percent during $`t/t_\mathrm{d}100`$ and to $`0.6`$ percent during $`t/t_\mathrm{d}200`$.

### 3.2 Initial conditions

Theoretically speaking, it will be highly desirable to obtain initial conditions (IC’s) from a cosmological scenario for the formation of triplets; however, this is out of the scope of the present work. None the less, we have chosen two typical IC’s for our triplets, namely: starting from (1) near ‘maximum expansion’ in a classical cosmology with $`\mathrm{\Lambda }=0`$, and (2) from virial equilibrium. We assume that galaxies are already formed at the moment of starting the simulations. We do not consider the possibility of secondary infall into the region where the triplets will evolve. This is a justifiable assumption in a low density universe similar the one where we appear to live (e.g., Bahcall 1999; Hradecky et al. 2000), where major infall of matter has ended at an earlier epoch ($`z1/\mathrm{\Omega }_0`$, with $`\mathrm{\Omega }_00.2`$).

For both types of IC’s, we do not consider the existence of a primordial common envelope of dark matter. This was motivated, on one hand, by the fact that some observations of galaxy clusters and groups tend to suggest that most of their dark matter is associated with the dark haloes of galaxies (e.g., Bahcall, Lubin & Droman 1995; Puche & Carignan 1991). On the other hand, there are indications that in binary galaxies their $`M/L`$-ratios might be rather small suggesting that large quantities of dark matter in a common halo are unnecessary \[Honma 1999\]. We make the hypothesis that these conditions prevailed once galaxies were already formed in a triplet, a condition that will have to be relaxed in future work where the effect of a self-gravitating common halo will have to be considered.

Results from simulations of groups of galaxies indicate that the merging rate is smaller when a common extended massive envelope of dark matter is included (Athanassoula, Makino & Bosma 1997). In the light of this result, we would expect that the number of mergers found in our numerical experiments to represent a sort of upper limit. The effect of a primordial concentrated dark halo is to increase the merging rate, so our results would tend to be a lower limit.

#### 3.2.1 Turnaround

For IC’s simulating the ‘turnaround’ moment for triplets, the centre of mass of each galaxy was sampled randomly from a uniform spherical density distribution of radius $`R_{\mathrm{max}}`$, which is assumed to be the size of the maximum radius of the density perturbation that led to the formation of a triplet. A small isotropic bulk velocity, derived from a Gaussian with mean zero and one-dimensional velocity dispersion $`\sigma `$, was imposed to each galaxy at $`t=0`$. In a sense, this isotropic $`\sigma `$ mimics a state of pre-virialization in our simulations (Davis & Peebles 1977; Łokas et al. 1996).

The numerical values of $`R_{\mathrm{max}}`$ and $`\sigma `$ were obtained as follows. In the spherical collapse model for the detachment of a density perturbation from the Hubble flow (Gunn & Gott 1972), the maximum expansion radius at a particular epoch $`t`$, for a standard cosmology with zero cosmological constant, $`\mathrm{\Lambda }=0`$, is:

$$R_{\mathrm{max}}(t)=\left(\frac{8GMt^2}{\pi ^2}\right)^{1/3},$$
(8)
where $`M`$ represents here the mass at a given epoch.

We make the hypothesis here that triplets are on the verge of complete collapse for the first time at the present epoch, $`t_0`$, which, in turn, leads to a particular value of $`R_{\mathrm{max}}`$. The present age of the universe $`t_0`$ depends, in the previous cosmology, on the current value of the Hubble constant $`H_0`$ and the dimensionless density parameter $`\mathrm{\Omega }_0`$. In particular, for an Einstein-de Sitter universe and a low-density universe we have, respectively (Gott et al. 1974):

$$H_0t_0=\{\begin{array}{cc}2/3\hfill & \hfill \mathrm{\Omega }_0=1\\ \mathrm{\hspace{0.33em}1}+\frac{1}{2}\mathrm{\Omega }_0\mathrm{ln}\mathrm{\Omega }_0\hfill & \hfill \mathrm{\Omega }_00\end{array}.$$

Therefore for $`\mathrm{\Omega }_00.2`$ we have $`H_0t_00.84`$ and, using our choice of the Hubble constant, $`t_00.84(13\mathrm{Gyr})10`$ Gyr. In Table 2, the turnaround radius (in Mpc) of a triplet of Galaxy-like objects ($`3\times 5.5\times 10^{11}`$ M), for different values of $`H_0`$ and $`\mathrm{\Omega }_0`$, are given; equation (8) was evaluated at a time $`t_0/2`$ to obtain $`R_{\mathrm{max}}`$ according to our stated hypothesis. Of the plausible values in Table 2, we have taken $`R_{\mathrm{max}}=500`$ kpc as our fiducial value for the simulations under this kind of initial conditions.

Our collapsing triplets are assumed to start with cold initial conditions satisfying a virial ratio of $`2T/|W|=1/4`$, as is common in simulations of groups of galaxies (e.g., Barnes 1985). A random velocity, sampled out of a Gaussian distribution with zero mean and 1-D velocity dispersion

$$\sigma =\frac{V_0}{2\sqrt{3}},\mathrm{with}V_0=\left(\frac{3GM_\mathrm{t}}{5R_{\mathrm{max}}}\right)^{1/2},$$
(9)
was assigned as bulk motion to each galaxy. Using our fiducial values for $`R_{\mathrm{max}}`$ and the triplet’s mass we obtain a value of $`\sigma 30`$ km s-1. This is consistent with the value found by Lake & Carlberg (1988) for the velocity dispersion of dark matter at the time of turnaround, $`35`$ km s-1.

#### 3.2.2 Virial

Adopting a realistic virial radius from which to sample the position of galaxies is not straightforward, because we are not considering the process of formation of triplets in a cosmological setting. Nevertheless, early point-like simulations made by Peebles (1970) suggest that by $`3t_\mathrm{c}/2`$, with $`t_\mathrm{c}`$ being the collapse time scale, a collapsing system has reached an approximate equilibrium state, after having rebounded, with a virial radius $`R_\mathrm{v}R_{\mathrm{max}}/2`$ (Gott & Rees 1975). According to Coles & Lucchin (1995), more recent $`N`$-body simulations of structure formation in an expanding universe tend to favour this time-scale.

Therefore, if we make the hypothesis that triplets at the present epoch are virialised, which is suggested by the crossing times of K-triplets, we have that $`R_\mathrm{v}=R_{\mathrm{max}}(t_0/3)/2`$; for our fiducial triplet we take $`R_\mathrm{v}=191`$ kpc. The bulk velocities of galaxies were chosen initially so as to satisfy the virial ratio $`2T/|W|=1`$ exactly.

### 3.3 Some computational aspects

To obtain the global dynamical properties of triplets, e.g., velocity dispersion and mean harmonic radius, we used the centre of each galaxy. We set the most bounded particle of the galaxy model to have 1 percent of its total mass, set its coordinates to be identically zero in the random realization, and identified this massive particle as the centre of the numerical galaxy. This numerical artifact allowed us to use this massive particle to trace the galaxy’s path and to compute its bulk velocity by time differentiation. A similar approach to identify the centre of a numerical galaxy was used by Aguilar & White (1985).

We performed 30 triplet simulations for each kind of initial conditions using the parameters and method discussed above. Each simulation differred in the initial seed used in the random realization to obtain the positions and velocities of galaxies at $`t=0`$. The collapse time for our fiducial triplet, in $`N`$-body units, is $`t_\mathrm{c}=\pi R_{\mathrm{max}}^{3/2}(2GM_\mathrm{t})^{1/2}290`$ (Gott 1975), and its half-mass free-fall time $`t_{\mathrm{ff}}=2^1\pi R_\mathrm{h}^{3/2}(GM_\mathrm{t})^{1/2}72`$ (Mamon 1990); here we have taken $`R_\mathrm{h}=R_{\mathrm{max}}/2`$. The dynamical time for initially virialised triplets is $`t_\mathrm{v}=(2R_\mathrm{v})^{3/2}(GM_\mathrm{t})^{1/2}130`$. Each triplet was evolved for $`t_\mathrm{n}=321.4`$ time units; i.e., $`10`$ Gyr. Energy was conserved to better than $`0.5`$ percent throughout this time for all simulations. Each run took $`4`$ hrs of cpu time in a Pentium-II 400MHz PC.

## 4 Results and Discussion

### 4.1 Global behaviour

In Figures 3 and 4 we display the XY-projection of the evolution of two particular triplet simulations starting from maximum expansion. These show some of the general qualitative characteristics found in the other simulations, including those starting in virial equilibrium. For example, a ‘binary’ galaxy is in general formed first while the third galaxy ‘orbits’ around it for some period of time. This binary instability is analogous to that found in the 3-body problem (e.g., Valtonen & Mikkola 1991) and, in this respect, supports a generalization of this behaviour to triplets of self-gravitating galaxies.

However, two important qualitative differences between the 3-galaxy and 3-body problem are found: (1) isolated triplets will eventually merge in their centre of mass provided enough time is given, and (2) no ‘sling-shot’ events were found for these initially bound systems. The physical reason for both of these characteristics stems from the fact that galaxies are able to absorb orbital kinetic energy, and angular momentum, and transfer it to their internal degrees of freedom. A clear manifestation of the previous is the large extent that matter in a merger can reach; e.g., in Fig. 4, the final system can occupy a region of $`1`$ Mpc in diameter. If we interpret our galaxy models as representing a spiral, this large region would be occupied primarily by dark matter and form a common halo by the present epoch. This has actually been observed in a simulation of a triplet of spirals (Aceves 2000b).

We emphasized above the issue of triplets being isolated because there are some indications that tidal interactions with large-scale structures might disrupt them over a Hubble time (Aceves 2000a). We believe that this may apply more properly to present-day W-triplets due to their large extent.

### 4.2 Merging histories

We have adopted a criterion of ‘loss-of-identity’ (Athanassoula et al. 1997) to determine when a pair of galaxies has merged. Although different merging criteria lead to differences in the number of mergers found at a particular time (e.g., García-Gómez & Athanassoula 1993), we consider the adopted criterion well-suited for our purposes, for it appears more closely related to an observational practice.

In particular, we consider that two galaxies to have lost their identity when

$$V_{ij}=|𝒗_j𝒗_i|<V_{\mathrm{rms}}/2,R_{ij}=|𝒓_j𝒓_i|<R_0/2;$$
(10)
where $`V_{ij}`$ and $`R_{ij}`$ are the relative three-dimensional velocity and distance between galaxies $`i`$ and $`j`$, and $`V_{\mathrm{rms}}`$ and $`R_0`$ are the root-mean-squared velocity and scale radius of a Plummer sphere, respectively. In the numerical implementation of criterion (10) a triple merger required a few time-steps to be detected, since Eq. 10 is a pair-wise recipe, but no unreasonable results were obtained.

In Figs. 5 and 6 we show the number of ‘galaxies’ remaining in the triplet at time $`t`$, both starting from maximum expansion and in virial equilibrium, respectively. A merged pair was counted as an individual ‘galaxy’. We note that the cosmic time, i.e., the time elapsed since the ‘big bang’, is different for each of the IC’s considered; this is because they were assumed to start a different epochs in the evolution of the universe. For IC’s starting at maximum expansion the cosmic time is $`t_{\mathrm{cos}}=t_0/2+t`$, where $`t`$ is the time span of the simulation (e.g., Ishizawa 1986); analogously, for initially virialised triplets we would have $`t_{\mathrm{cos}}=t_0+t`$.

It is interesting to note that by the present epoch, $`t_n160`$, only about 10 percent of triplets starting from maximum expansion have merged and about $`60`$ per cent have not suffered any mergers. At the end of the simulation, i.e. about 5 Gyr after $`t_0`$, about 30 percent have managed to survive without any merging.

Although there are not many observational studies regarding triplets’ characteristics, the previous results appear qualitatively consistent with the available observations. For instance, a rapid visual inspection of K-triplets in NASA’s database NED does not show strong triple mergers but their luminous components look fairly separated in most of the images, though this may be a consequence of the selection criteria. One example of a possible triple merger in process is in Hickson’s group 95 (Iglésias-Páramo & Vílchez 1997), but since this is a group one cannot even strictly compare such ‘triplet’ with our results, and another is that found by Heidt et al. (1999).

For initially virialised triplets, 4 ($`15`$%) merge completly in within a dynamical time, while 2 ($`7`$%) do not present any merger and 19 ($`60`$%) only present binary mergers in $`10`$ Gyr. This behaviour is probably related to the fact that virialised triplets tend to maintain their initial conditions longer, with higher velocity dispersion, hence being less disposed to satisfy the merging criterion (10). Obviously, these results would be altered if a primordial common dark halo had been introduced.

We do not find an overmerging problem for triplets under the kinds of initial conditions considered here, contrary to what is usually found in numerical simulations of compact groups of galaxies (Hickson 1997). The results obtained also show that there are stable states in the 3-galaxy problem, contrary to what is sometimes believed (Valtonen & Flynn 2000). We must recognize, however, that since we are not modeling the actual process of formation of a triplet, the above results need to be taken as a first approximation to the problem of merging activity in triplets.

### 4.3 Geometrical properties

Several authors have used the AA-map to characterize the geometrical properties of triplets and three-body systems (e.g., Chernin et al. 1994; Heinämäki et al. 1998). Although this approach does not provides much information on the dynamical state of triplets, since a triplet initially in the L-area can migrate toward other regions at later times (e.g., Ivanov, Filistov & Chernin 1995), can none the less provide us with another means to test our models.

In Figs. 7 and 8 we show the AA-map at different times for triplets starting from maximum expansion and virial equilibrium, respectively. Although it is somewhat difficult to assess quantitatively these results and to compare them with the observations of K and W-triplets, in part due to the small number of points for a statistical analysis, it is interesting to note that the AA-map indicates that no strong hierarchical structures are present at the current epoch ($`t6`$ Gyr) for collapsing triplets. On other hand, for initially virialised triplets we obtain a slight excess of triplets in the H-area at $`t=0`$ that is still present after a few giga-years of evolution.

We would like to point out also that it is possible that the lack of many observational triplets in the H-area of the AA-map may be due to an observational bias. Real galaxies are able to merge and disappear from the AA-map at some time. Once the haloes of galaxies start to overlap the merging process is rather rapid ($`1`$ Gyr), and this will reduce the probability of finding such extreme H-structures in a sky survey. If some W-triplets are truly physical hierarchical triplets one would expect them to disappear from the H-area of the AA-map in the order of a few Gyr.

The previous results lead us to consider that our models are not inconsistent with the AA-map of observed triplets, both compact and wide. However these results do not allow us to establish a preference for the current dynamical state of triplets. It appears, on other hand, that a large primordial common halo of dark matter is not required to avoid the the existence of hierarchical triplets, as suggested by other authors (e.g., Chernin et al. 1994).

### 4.4 Dynamical quantities

In order to compare our models with observations we have calculated, for both kinds of initial conditions, the one dimensional velocity dispersion, the mean harmonic radius, the dimensionless crossing time, and the virial and median mass along three orthogonal ‘lines-of-sight’ at different times; i.e., we mimic the observation of 90 triplets.

In Table 3, we present the median values of the above quantities in $`N`$-body units for triplets starting at maximum expansion, and in Table 4 the results for initially virialised triplets. The crossing time is expressed in terms of the $`N`$-body equivalent of our adopted Hubble time. The times in Tables 3 and 4 correspond to $`t_\mathrm{a}\{0.5,1,\mathrm{},9\}`$ Gyr in our fiducial system of units. In calculating the median of a given quantity at a certain time we have disregarded triplets that have a merger.

In Figs. 9 and 10, histograms for the dynamical quantities at times $`t_\mathrm{a}0.5,3,6`$ Gyr normalized to the initial number of ‘observed’ triplets are shown for collapsing and virialised triplets, respectively. This normalization was chosen to show the evolutionary trends in the ensemble of simulations. Downward and upward arrows indicate median observed values for K and W-triplets, respectively. The bar in the mass histograms denotes the total mass of the system in $`N`$-body units.

In general, as expected, the median properties of the virialised triplets tend to evolve much slower than those for the collapsing ones, which is perhaps most noticeable in the small change of the median velocity dispersion (Table 4). On the other hand, a significant evolution toward the establishment of a virial equilibrium among the translational bulk degrees of freedom of galaxies is observed for collapsing triplets. This is noticeable as their virial mass estimate approaches their true mass. However, because galaxies are self-gravitating systems, part of their mutual orbital kinetic energy is transformed into internal kinetic energy and, thus, no complete virial equilibrium in their translational bulk degrees of freedom can be achieved. A virial equilibrium, but in the complex mixture of the internal and bulk degrees of freedom of galaxies, tends to be established during their interactions.

From Tables 3 and 4, we observe that no overestimates of the median virial mass are obtained at any time, except at the end of the simulations for virialised triplets. Overestimates were found in few of our numerical experiments, but only for brief intervals of time when using the median mass estimator. At most a $`20`$ per cent overestimate was found for virialised IC’s and $`10`$ per cent for collapsing triplets. After examining each run, we found that several triplets were a factor at most of about $`4`$ higher than the true mass of the triplet, but this happened only along a particular line-of-sight while along the other two an underestimation always occurred. These overestimates happened before the galaxies overlapped of the galaxies; i.e., when galaxies still behaved somewhat as point-like particles and not as self-gravitating systems. These results are in contrast to others that indicate that overestimates of even a few orders of magnitudes may be obtained in triplets.

The results of Table 3 indicate that triplets arriving at the present epoch ($`t_\mathrm{n}160`$) from maximum expansion have their median virial masses underestimated by $`35`$ percent, with a somewhat better agreement when using the median mass estimator. The physical reason behind this underestimate, in these truly bounded systems, is the already-mentioned transfer of orbital to internal energy. This means, in observational practice, that the use of the bulk velocities of galaxies –even if we had three-dimensional information– has the intrinsic propensity to underestimate the total mass of the system. A similar conclusion was reached by Barnes (1985) in his study of small groups of galaxies. For initially virialised triplets (Table 4) the median mass estimator tends to provide a better estimate, though it may overestimate in some cases, by less than $`30`$ percent.

Overall, we consider the median mass estimator to perform better than the virial mass estimator in these low-multiplicity systems. The accuracy of the mass estimate depends in an important way on the dynamical status of the system. To derive a time-varying expression for the median or virial mass estimator is not the purpose here. See however Dolgachev & Chernin (1997) for a study developed to take this into account using point-like particles as a galaxy model.

At the present epoch, $`t_\mathrm{n}160`$, our models predict that triplets starting from maximum expansion should have a median $`R_\mathrm{H}290`$ kpc and $`H_0\tau _\mathrm{c}0.43`$. When comparing these values with the properties of K-triplets (Table 1), we find them higher by about a factor of 4 and 10, respectively. Therefore, it appears that our models yield results that are difficult to reconcile with the observations. This situation is worse when we note that a median of $`\sigma 52`$ km s-1 is predicted, a value which is about half of the observed one for K-triplets. A maximum of $`\sigma 64`$ km s-1 is obtained at the end of the simulations. For virialised IC’s, the maximum median is $`\sigma 80`$ km s-1, occurring at $`t=0`$. Afterwards, the median $`\sigma `$ tends to decrease somewhat due to the transfer of orbital to internal energy during the interactions of galaxies.

The previous results, taken at face value, tend to rule out the models presented for triplets with either kind of initial conditions considered. However, there are some observational matters and numerical aspects that indicate that our models might not be an unreasonable scenario for the global evolution of triplets. In $`\mathrm{\S }`$4.5 we will establish which kind of initial conditions we favour the most.

On other hand, when comparing our results to observations we are assuming implicitly that K-triplets form an homogeneous sample, both in luminosity ($``$ mass) and morphological type, but as can be seen from the properties of K-triplets \[Karachentsev et al. 1989\] this is not the case. In this respect our ensemble of numerical triplets forms, by construction, a homogeneous ‘catalog’. A mass-spectrum, for example, can have effects on the numerical values obtained for $`\sigma `$ and $`R_\mathrm{H}`$ that are not easy to estimate. This motivated us to estimate the velocity dispersion of K-triplets weighted by the luminosity of their member galaxies in order to compare it with our numerical results. Weighting by mass appears to be a more uncertain procedure. The median weighted velocity dispersion obtained was $`\sigma 90`$ km s-1 (about 25 per cent lower the value in Table 1), which in turn is closer to our results.

Different selection criteria can yield significantly different results for the dynamical quantities in galaxy systems. For example, in a recently-announced catalog of small galaxy groups by Makarov & Karachentsev (2000), using physically motivated selection criterion, the following median values, with quartiles, for triplets are reported: (a) velocity dispersion $`41_{18}^{+19}`$ km s-1, (b) mean harmonic radius $`191_{88}^{+157}`$ kpc, and dimensionless crossing time $`0.15_{0.09}^{+0.16}`$ (using $`H_0=70`$ km s-1 Mpc-1). These authors provide a harmonic radius value not corrected for projection effects (i.e. lower by a factor of $`\pi /2`$ than when using Eq. 1) and use the crossing time definition of Gott, Wrixon & Wannier (1973) (Makarov, private communication). The latter definition is $`t_\mathrm{c}=3R_\mathrm{H}(5^{3/2}\sigma )^1`$, that, when comparing with Eq. 3, yields $`\tau _\mathrm{c}4t_\mathrm{c}`$. Therefore, comparing our results at $`t_0`$ for collapsing triplets (Table 3) with the values found by Makarov & Karachentsev a much better agreement is found. Furthermore, we note that while all the previous comparisons were done among the median values, we find that $`10`$ per cent of collapsing triplets have $`\sigma 100`$ km s-1 at $`t_0`$ and small RH’s ($`\mathrm{\S }`$4.5).

Part of the disagreement of our median results with, for example, the observed velocity dispersion of K-triplets may be somewhat alleviated if we use a larger mass as our fiducial galaxy model, but it does not satisfactorily solve the problem. Zaritsky et al. (1997) after examining the distribution and kinematics of satellites around spiral galaxies conclude that the halo masses of spirals are in the excess of $`2\times 10^{12}`$ M within $`200`$ kpc. Using this type of galaxies, the scaling factors $`\{l,\mu \}`$ in Eq. 6 are $`\{4.4,39.2\}`$. From this we obtain that $`v_\mathrm{a}=660v_\mathrm{n}`$ km s-1, which leads to a median $`\sigma 80`$ km s-1 at $`t_0`$ for collapsing triplets. This new scaling provides a much better agreement with the velocity dispersion observed in K-triplets. However, the discrepancy with observations moves now to the $`R_\mathrm{H}`$ median value that results in $`425`$ kpc, i.e., about an order of magnitude larger than the observed value and marginally consistent with the results of Makarov & Karachentsev (1999). The time scale does not change by much. One may try to force a better agreement by selecting a galaxy with a mass of $`10^{12}`$ M within $`100`$ kpc, but this appears in conflict with the fact that not many galaxies satisfying these conditions, e.g., large ellipticals, exist in K-triplets. To properly address these scaling matters, more realistic simulations including a mass-spectrum would be required.

The models considered here do not properly apply to present day observed W-triplets. A better model for these would be, perhaps, to assume that they are in the early phases of their turnaround. With this assumption $`R_{\mathrm{max}}1`$ Mpc and $`\sigma 30`$ km s-1. The latter scenario has been considered briefly elsewhere (Aceves 1999a), where it was found that we will have to wait for about another Hubble time for W-triplets to reach a compact configuration similar to the one of K-triplets. See also Dolgachev & Chernin (1997). However, the current large $`\sigma 100`$ km s-1 for W-triplets does not seems to favour such possibility.

We note that W-triplets appear to be a less homogeneous sample than K-triplets, for they include triplets with a very low harmonic radius similar to those of K-triplets (Fig. 1). It is not the purpose here to make an analysis of the selection criteria used to compile W-triplets, but we point out that some of them would appear more properly classified as K-triplets. The overlapping in the RH-distribution is interesting in itself, however, since some triplets classified as W-triplets may constitute a link between those that are in the early phases of collapse and those that are in the more advanced phase.

Larger samples of triplets and more detail studies of their properties are required to build a more solid observational database for future comparisons with theoretical models. Efforts toward expanding the sample of triplets, for example, in the southern sky are presently being done by several groups (e.g., Karachentseva & Karachentsev 2000).

### 4.5 Evolutionary Trends

An important question to investigate is how K-triplets have attained their present configuration and dynamical properties. From the previous results, it appears that triplets starting from a diffuse configuration, followed by a collapse, is a probable scenario.

Perhaps one of the best methods to discern the evolution of a dynamical system is to look into its phase-space structure. However this is not possible in galaxy systems since we do not have yet enough information to construct a snapshot of their physical phase-space. None the less, we may construct a pseudo phase-plane for triplets by using their velocity dispersion and mean harmonic radius, compare them with the numerical results, and see if some evolutionary trends are discernible.

In Fig. 11, we show the pseudo phase-plane (hereafter, phase-plane for brevity) for our numerical triplets at different times spanned from $`t=0`$, both for collapsing and virialised initial conditions. In Fig. 12 the phase-planes of K and W-triplets are shown. Again it is difficult to assert some conclusions from a comparison of the phase-planes, but it is interesting to note that some tendencies exist. For example, K-triplets show, on one hand, the general trend followed by systems that had suffered a collapse out of a much wider configuration: a decrease in their $`\sigma `$ as $`R_\mathrm{H}`$ increases. This behaviour is similar to that exhibited by the numerical triplets that have started from maximum expansion in our simulations. On the other hand, W-triplets have a phase-plane that resembles qualitatively systems that are on the verge of suffering a gravitational collapse.

From the phase-plane of collapsing triplets shown in Fig. 11, we observe that about 10 percent of them have dynamical quantities ($`\sigma `$ and $`R_\mathrm{H}`$) similar to the median observed in K-triplets; c.f. Table 1. This is also observed in some triplets with virial IC’s at $`t_0`$.

In general, we consider that the results and data shown in these phase-planes, along with earlier discussions, suggest that: K-triples are probably the most extreme cases of triplets that started from a diffuse configuration, attained their current dynamical properties via a gravitational collapse process, have manage to survive without suffering major mergers, and are near virial equilibrium in their bulk degrees of freedom.

## 5 Summary and Conclusions

A simplified $`N`$-body model for the dynamical evolution of triplets, using self-consistent galaxies, has been proposed. Two types of initial conditions, assuming a classical cosmology with $`\mathrm{\Lambda }=0`$, have been investigated, namely ‘maximum expansion’ and virial equilibrium. Triplets starting from maximum expansion were imposed the condition that, assuming all their mass to be distributed homogeneously, would collapse at the present epoch, and those in virial equilibrium to be in such state now. No primordial common dark halo was included in our simulations.

The main results of this work are as follows:

1. The use of self-consistent galaxies to study the dynamics of triplets is of fundamental importance. No ‘sling-shot’ events occur during a triple interaction in an initially bound system, in contrast with what is a common event in the 3-body problem.
2. A rather low number of triple mergers ($`10`$ per cent) are expected to occur at the present epoch for systems starting from maximum expansion. This appears to be in concordance with observations. Initially virialised systems yield $`5`$ per cent of triple mergers in $`10`$ Gyr of evolution. Therefore, no overmerging problem is found for triplets. This contrasts to what it is usually found in numerical simulations of small compact groups.
3. Under both types of initial conditions considered, we find that the 3-galaxy problem has stable states; i.e., there are certain conditions that do not lead to a triple merger even after $`10`$ Gyr of evolution.
4. The geometrical properties of numerical triplets, as indicated by their homologous AA-map, are not inconsistent with those of observed triplets, both compact and wide. Physical triplets in the extreme hierarchical region of the AA-map would suffer a rapid dynamical evolution toward a binary merger leading to a depletion of triplets in such region, since they would no longer would be counted as triplets in an observational survey.
5. Median values of the dynamical parameters for triplets that have started as a diffuse system and are collapsing for the first time at the present epoch are not inconsistent with the data of Karachentsev’s compact triplets. A much better agreement is reached with recent data on triplets. We find that about $`10`$ percent of the simulated triplets reproduce well the K-triplets median dynamical quantities at the present epoch.
6. The median of the virial mass estimator does not in general overestimate mass, but, on the contrary, it underestimates it. For triplets in an advance stage of collapse at the present epoch an underestimate of $`35`$ percent is found. Therefore if K-triplets are truly bound physical systems, with no primordial common dark halo, their $`M/L`$-ratios are being underestimated by the same factor. Otherwise, the underestimate may be larger. The median mass estimator appears as a somewhat better mass estimator for triplets than the virial mass estimator.
7. The plane $`\sigma R_\mathrm{H}`$ of numerical triplets suggests that K-triplets have acquired their present configuration through a clustering process, and lead us to conclude that they are the triplets with the most extreme values in the literature, i.e., lowest harmonic radius and highest velocity dispersion. Some dynamical properties of K-triplets are also consistent with a virial equilibrium state in their bulk degrees of freedom at the present time.
8. Our results suggest that the presence of a primordial, large, massive, common dark halo might not be a strong requirement to account for the dynamical evolution of triplets, but they cannot exclude it.

In order to better compare with observations, future models should include, e.g., a mass-spectrum for galaxies, consider different morphological types, and investigate the role of a primordial common dark halo in a cosmological setting.

## Acknowledgments

The Spanish Ministry of Foreign Affairs (MUTIS program) and México’s CONACyT (Proyecto I35546-E) are thanked for financial support at different stages of this work. We thank greatly Gary Mamon for fruitful discussions on the issue of initial conditions. Jack Sulentic and Héctor Velázquez are also thanked for comments on this work, and Michael Richer for help with the english writing.


text:
# White Dwarfs: Contributors and Tracers of the Galactic Dark–Matter Halo

## 1 Introduction

Recently, Oppenheimer et al. (2001; O01 hereafter) found 38 white dwarfs (WD) in a sample of 99 high proper motion WDs, which were claimed to have kinematics inconsistent with that of both the old stellar disc and the thick disc, and therefore form a newly discovered halo population with an inferred density at least an order of magnitude higher than previously thought (e.g. Gould et al. 1998). Reid et al. (2001) have criticized this result and conclude that these WDs form the high-velocity tail of the thick disc, based on a comparison of their $`U`$$`V`$ velocity distribution with that of M-dwarf stars. Unfortunately, neither conclusion is statistically supported.

In this proceeding, we discuss a maximum-likelihood analysis of the complete sample of 99 WDs and derive the local phase-space density of thick-disc and halo WDs. This then allows us to draw more robust conclusions. For a discussion of the selection of the sample, we refer to O01. Details of the likelihood analysis and phase-space density estimate can be found in Koopmans & Blandford (2001; KB01 hereafter). For clarity, throughout the proceeding the low- and high-velocity dispersion components are referred to as the thick-disc and halo population, respectively. This does not imply that all WDs could not have originated from the thin or thick disc (e.g. Hansen 2001; KB01).

## 2 The White–Dwarf Velocity Distribution Function

We model the local ($`\stackrel{}{x}_0`$) velocity distribution function (VDF) as a superposition of two Schwarzschild (i.e. Maxwellian) VDFs, i.e. $`f_\mathrm{B}(\stackrel{}{x}_0,\stackrel{}{v})`$. We assume a constant space density of WDs throughout the surveyed volume. Each component has two free parameters, the radial velocity dispersion ($`\sigma _U`$) and the asymmetric drift ($`v_\mathrm{a}`$). In addition, the parameter $`r_n`$ is the ratio of the thick-disc to halo WD number densities ($`n`$). We assume that the vertex deviations are zero and that the ratios of the radial ($`U`$), azimuthal ($`V`$) and vertical ($`W`$) velocity dispersions of the ellipsoidal VDFs are $`\sigma _U`$:$`\sigma _V`$:$`\sigma _W`$=1:2/3:1/2, in agreement with observations (see KB01 for references). The probability that an observed WD with a velocity vector on the sky, $`\stackrel{}{p}=(v_l,v_b)`$, is drawn from this VDF becomes $`P=_{\mathrm{}}^{\mathrm{}}(v_l^2+v_b^2)^{3/2}f_\mathrm{B}(\stackrel{}{v}=\stackrel{}{p}+v_r\widehat{r})dv_r/_𝒱(v_l^2+v_b^2)^{3/2}f_\mathrm{B}(\stackrel{}{v}=\stackrel{}{p}+v_r\widehat{r})\mathrm{d}^3v`$, where $`v_r`$ is the radial velocity along the line-of-sight and $`v_l`$ and $`v_b`$ are the velocities projected on the sky in Galactic coordinates. This assumes that the WDs are all proper-motion limited ($``$90% are). $`𝒱`$ is then the velocity space in which the WDs could have been found, given the restriction that $`v_l^2+v_b^2>(\mu _0r)^2`$, where r is the distance to the WD and $`\mu _0`$ is the survey’s lower limit on the proper motion ($`\mu _0`$0.33 ′′ yr-1). For the $``$10% magnitude limited WDs, we find that their maximum detection volumes ($`V_{\mathrm{max}}`$) are nearly identical. For these, we modify the likelihood function, although the differences in parameter estimates are within the errors quoted below. By varying the five free parameters and optimizing the log-likelihood, $``$=$`_i\mathrm{log}(P_i)`$, of the sample, we solve for their most likely values and their error range. The results (excluding 3 thin-disc WDs) are: (a) $`\sigma _{U,\mathrm{d}}`$=62$`{}_{10}{}^{}{}_{}{}^{+8}`$ km s-1, (b) $`\sigma _{U,\mathrm{h}}`$=150$`{}_{40}{}^{}{}_{}{}^{+80}`$ km s-1, (c) $`v_{a,\mathrm{d}}`$=50$`{}_{11}{}^{}{}_{}{}^{+10}`$ km s-1, (d) $`v_{a,\mathrm{h}}`$=176$`{}_{80}{}^{}{}_{}{}^{+102}`$ km s-1and (e) $`r_n`$=16$`{}_{11}{}^{}{}_{}{}^{+30}`$ (KB01).

## 3 The Local Halo White–Dwarf Density

To normalize the local phase-space density of halo plus thick-disc WDs, i.e. $`n_{0,\mathrm{WD}}^{\mathrm{td}+\mathrm{h}}f_\mathrm{B}(\stackrel{}{x}_0,\stackrel{}{v})`$, it remains to estimate their local density, $`n_{0,\mathrm{WD}}^{\mathrm{td}+\mathrm{h}}`$. We do this in the conventional way, by summing the $`1/V_{\mathrm{max}}`$ values of the halo plus thick-disc WDs, where $`V_{\mathrm{max}}`$ is the smallest of the two volumes in which a WD could have been detected, when limited either by its proper motion or by its magnitude (e.g. O01; KB01). We find that $``$90% of the WDs are proper-motion limited and three WDs are very likely associated with the thin disc. The latter WDs are removed from the sample, as we are only interested in the halo plus thick-disc density, for which we determined the VDF. We find $`n_{0,\mathrm{WD}}^{\mathrm{td}+\mathrm{h}}=(1.9\pm 0.5)\times 10^3`$ pc-3, of which $`n_{0,\mathrm{WD}}^\mathrm{h}=1.1_{0.7}^{+2.1}\times 10^4`$ pc-3 belongs to the halo, given $`r_n`$ found from the maximum-likelihood analysis of the sample. The local halo plus thick-disc density estimate agrees within $``$$`\sigma `$ with the estimate from Reid et al. (2001). The halo density is 5.0$`{}_{3.2}{}^{}{}_{}{}^{+9.5}`$ times higher than the best previous estimate of $`\stackrel{~}{n}_{0,\mathrm{WD}}^\mathrm{h}=2.2\times 10^5`$ pc-3 from Gould et al. (1998), inferred from the mass function of halo stars. Locally, our density estimate amounts to 0.8$`{}_{0.5}{}^{}{}_{}{}^{+1.6}`$% of the nominal local dark-matter halo density of $`8\times 10^3`$ M pc-3 (Gates et al. 1995).

## 4 The Global Halo White–Dwarf Density

Given the local phase-space density of halo plus thick-disc WDs and a potential model for the Galaxy, we can estimate the global white-dwarf density in a large part of the Galactic halo, using Jean’s theorem (e.g. May & Binney 1986). For simplicity, we assume a spherical logarithmic potential, $`\mathrm{\Phi }(r)=v_\mathrm{c}^2\mathrm{ln}(r/r_\mathrm{c})`$, with rotation velocity $`v_\mathrm{c}`$=220 km s-1 and a solar radius of $`r_\mathrm{c}`$=8 kpc. For the halo WDs, which are less affected by the disc/bulge potentials this is probably a reasonable assumption, although this does not hold for the lower-velocity thick-disc WDs. A analytic expression for the phase-space density of halo WDs can then be derived for any point in halo that is in dynamic contact with the Solar neighborhood (KB01) . From this phase-space density model, we find an oblate ($`q=(c/a)_\rho 0.9`$) distribution of halo WDs with a total mass inside 50 kpc of $``$2.6$`\times `$$`10^9`$ M and a radial mass profile $`n(r)r^{3.0}`$. We expect $`q`$ to decrease further if a proper flattened Galactic potential is used. This halo WD mass amounts to $``$0.4% of the total Galactic mass inside 50 kpc, $``$4% of the Galactic stellar mass, or $`\mathrm{\Omega }_{\mathrm{WD}}10^4`$. Including disc WDs, the total mass of WDs in the Galaxy is $``$9% of the stellar mass, in agreement with standard population synthesis models (Hansen 2001; KB01).

## 5 The LMC Microlensing Optical Depth

Given our analytic expression of the halo WD phase-space density, we can integrate over the line-of-sight towards the Large Magellanic Cloud (LMC) to estimate the microlensing optical depth of this population. We find $`\tau _{\mathrm{WD}}^\mathrm{h}1.3\times 10^9`$, which is $`10^2`$ times lower than observed (Alcock et al. 2000). The thick-disc WDs have a $``$3 times higher optical depth. The integrated halo plus thick-disc WD optical depth is therefore still 1–2 orders of magnitude too small.

## 6 Conclusions, Discussion & Future

Modeling of the local phase-space density of halo plus thick-disc WDs indicates that there is a population of pressure-supported halo WDs with a local density 5.0$`{}_{3.2}{}^{}{}_{}{}^{+9.5}`$ times higher than previously estimated, although globally ($`r50`$ kpc) the mass contribution is negligible ($``$0.4%). Both the low- and high-velocity populations have similar color, magnitude and ages distributions (Hansen 2001; KB01). Two particularly interesting surveys to search for more high proper motion WDs are: (i) a wide survey with the Advanced Camera for Surveys (ACS) on the Hubble Space Telescope, which we estimate should find $``$5 WDs (KB01) and (ii) a similar survey as that by O01 towards the Galactic anti-center to constrain the WD velocities perpendicular to the Galactic plane.

The white-dwarf data set from O01 and results derived from it (i.e. our most likely model) have met a number of interal and external consistency checks (KB01) and agree with all observational constraints (e.g. direct observations, metal pollution of the ISM, microlensing, etc.) known to us. The color-color relation of the WDs shows them to be quite young (Hansen 2001) and have a birth rate roughly that expected from Galactic population synthesis models (KB01). Why then do a large fraction of these WDs have such high spatial velocities? To explain this we postulate the following; the majority of high-velocity WDs in the Galactic halo originate from multiple (i.e. $`N3`$) stellar systems in the Galactic disc. In case of a striple stellar system, stellar evolution (i.e. mass loss and/or transfer) of the (probably more massive) inner binary stars changes the ratios of orbital periods and eventually destabilizes the system. At that point, the lightest star, a WD in the inner binary, is ejected from the system into the halo – possibly through a slingshot of the outer star – leaving behind a recoiled compact binary (KB01). This postulate implies a direct relation between the radial surface density of the stellar disc, the star-formation rate and population synthesis models, which all seem to be in agreement (KB01), but still need considerable study.

## Acknowledgments

We thank Ben Oppenheimer for valuable discussions and sending their tabulated results. We thank David Graff and Andy Gould for pointing out an error in the normalisation of our likelihood function. LVEK thanks Priya Natarajan for organising a productive and well-organised meeting. This research has been supported by NSF AST–9900866 and STScI GO–06543.03–95A.

## References