format: text

text:
# Tidal Effects in Clusters of Galaxies

## 1 Introduction

This work is motivated by the wealth of recent observational discoveries of the strong morphological, chemical, and dynamical evolution of galaxies in clusters. Transformations of luminous galaxies, the existence of separate massive dark objects, and the evolution of cluster substructure intensify the need for theoretical understanding. Several different phenomena have been proposed as the “likely” explanation (see §2). While all of them take place in one form or another, it has eluded our understanding which process is primarily responsible for the evolution.

I argue that tidal interactions with neighboring galaxies and with the dark cluster halo are the favored scenario. Tides operate even before the cluster dynamically relaxes and affect both stars and gas in the infalling galaxies. Tidal effects add enough random motion to transform thin galactic disks into the kinematically hot thick-disk configurations characteristic of S0 galaxies. Tidal disruption of low density galaxies may provide a large fraction of the diffuse intracluster light.

The relative importance of tides depends on the cosmological context. In an $`\mathrm{\Omega }_0=1`$ universe, gravitating structures assemble late and the tidal effects are only now starting to play a role in shaping the galaxies, whereas in an open universe, $`\mathrm{\Omega }_00.4`$, galaxies have formed early and experienced a substantial metamorphosis. The currently fashionable model, a flat universe with $`\mathrm{\Omega }_0=0.4`$ and a cosmological constant, provides the intermediate case. In order to incorporate these effects into the galactic dynamics self-consistently, I perform numerical simulations of the evolution of a medium-size cluster of galaxies in the three cosmologies.

The choice of a numerical scheme is very important here, as the required dynamical range of simulation is huge: from hundred parsecs to model the galactic dynamics to tens of Mpc to derive the cluster tidal field correctly. In order to achieve such a range with the existing codes, I introduce a new scheme of separating the cluster tidal field from the galactic model. First, I run a lower resolution cosmological simulation that captures the details of cluster formation and its time-varying potential. At some early redshift, $`z_g=5`$, I identify the probable galactic halos and follow their trajectories until the present. Along each trajectory, I calculate and store the components of the tidal shear tensor. Then, I resimulate individual galaxies with a very high resolution applying the cluster tidal field as an external perturbation. The results of the stellar dynamical galaxy simulations are presented in a companion paper (Gnedin 2003, hereafter Paper II). No single simulation of a cluster of galaxies (§2.5) has been able to achieve comparable resolution on the galactic scale.

I describe a careful procedure to extract the tidal field from the grid PM simulation. Because of the particle noise, it is necessary to apply a smoothing filter extending over a range of several grid cells. Also, large halos contribute significantly to the tidal field around their centers. Since we are interested only in the external tidal force on the galaxy, this self-contribution needs to be removed. Special care is taken to evaluate the contribution of the galactic particles and subtract it from the total tidal tensor.

Also, the resulting tidal field depends on the correct calculation of the galactic trajectories. I use two independent methods to follow the trajectories: in one, the galactic centers of mass are assumed to be test particles following the overall cluster potential; in the other, the center of mass is calculated using the densest subset of particles constituting an actual halo. The sampling of halos at each time step provides an excellent time resolution of the tidal field, of order $`10^7`$ yr. Finally, the distribution of halos and their relative velocities enable us to estimate the contribution of close encounters to the external tidal force.

In §2, I review the present observational status of galaxies in clusters at various redshifts and the theoretical mechanisms connecting pieces of the puzzle. In §3, I describe the numerical method used to simulate the self-consistent evolution of the three clusters starting with the primordial fluctuations. In §4, I identify the galactic halos and calculate their trajectories. In §5, I study the possible mergers and close encounters between galaxies, including the effects of halo clustering. Finally in §6, I calculate the tidal field along the galaxy trajectories and use it in §7 to estimate the amount of tidal heating of galaxies.

## 2 Scenarios of Galactic Evolution

The idea that the environment plays a role in the evolution of galaxies has been recognized long ago. Spitzer & Baade (1951) suggested that collisions between galaxies might be responsible for transforming spirals into the gas-free S0s. Toomre (1977) advocated a merger of two spirals to produce an elliptical. Richstone (1976) argued for the truncation of dark matter halos of galaxies by the tidal field of the cluster. Gunn & Gott (1972) proposed the stripping of gas by the ram pressure of the intracluster medium. Finally, Merritt (1983) and most recently Moore at al. (1998,1999) considered tidal interactions between the galaxies and the dense cluster core and between the galaxies themselves. I review the observational picture and then discuss briefly the pros and cons of each of these scenarios.

### 2.1 “Spirals vs Ellipticals” or “Spirals vs S0s”?

The Hubble Space Telescope has been instrumental in shaping our current understanding of distant clusters of galaxies. The unique resolving power of the HST brings us exquisite details of the distorted profiles of the otherwise normal disk galaxies. Strongly depleted in nearby clusters, spirals make up to 50% of the galaxies in clusters at redshift $`z0.5`$ (Oemler et al., 1997; Ellis et al., 1997; Smail et al., 1997). These results are similar to the Butcher–Oemler (1978) effect, an abundance of blue irregular galaxies at an intermediate redshift. Spectroscopic evidence (Dressler et al., 1997) suggests bursts of recent, or even current, star formation. This picture is drastically different from the rich present-day clusters, populated mostly by featureless elliptical or S0 galaxies which follow the local morphology-density relation (Dressler, 1980). It is well established now that a strong chemical and morphological evolution has taken place in clusters of galaxies.

What are the options? Did the $`z=0.5`$ spirals become the $`z=0`$ ellipticals? On the contrary, the HST data show that the fraction of elliptical galaxies increases with redshift, although not as strongly as that of spirals. Both of these effects occur at the expense of S0s. Their population is depleted by a factor of two or three at $`z0.5`$ (Dressler et al., 1997). Observations suggest that star formation has ceased in the ellipticals by $`z0.6`$, but continued in the S0s at least until $`z0.3`$ (van Dokkum et al., 1998). The ellipticals must have exhausted their supply of fresh gas at an earlier time, as indicated also by their very high central densities (Kormendy, 1989).

Even at low redshift the situation is far from being smooth (Conselice & Gallagher, 1999). New ground-based imaging of the Coma and five other nearby clusters reveals a variety of distorted, tailed, “line”, and dwarf galaxies. These are indicators of the ongoing dynamical evolution. Many galaxies are disturbed and some show blue colors. Also, individual groups of dwarf galaxies have been found in Coma. These results were not expected for the smooth and regular clusters.

A new observational puzzle has been presented by the discovery of massive dark objects in the cluster CL 0024+1654, through the technique of image reconstruction of gravitational lensing (Tyson, Kochanski & Dell’Antonio, 1998; Kochanski, Dell’Antonio & Tyson, 1998). The dark objects are not associated with any of the luminous galaxies in the cluster and may contain as much mass as $`10^{10}M_{\mathrm{}}`$. One of the possible explanations is that these objects are just “failed” galaxies, in which gaseous dissipation was insufficient to form stars. Alternatively, they can be stripped dark matter halos of other galaxies. Analytical arguments and previous studies (§2.3) suggest that extended galactic halos may be effectively swept away by the tidal field of the cluster.

### 2.2 Mergers of Galaxies

Collisions of galaxies are assisted by the dense environment of the cluster. However, in a virialized cluster, galaxies move with the relative velocities of the order 1000 km s-1 and even on a head-on trajectory will pass through each other without doing much harm to the stellar disks. The gas can be shocked and stripped away in such a collision (Spitzer & Baade, 1951), which would leave two gas-free galaxies. But a simple estimate argues that close encounters are very infrequent (however see §5 and Appendix A).

Consider a virialized cluster of galaxies with a one-dimensional velocity dispersion $`\sigma _{cl}=1000`$ km s-1 and a virial radius $`R_{cl}=1`$ Mpc. Assuming that the cluster contains $`N_g=1000`$ galaxies uniformly distributed inside this radius, the number of close encounters a given galaxy of size $`R_g=10`$ kpc will experience in a Hubble time ($`t_H=10^{10}`$ yr) is

$$N_{\mathrm{enc},1}\frac{N_g}{\frac{4\pi }{3}R_{cl}^3}\pi R_g^2\sqrt{2}\sigma _{cl}t_H1.$$
(1)
Here I have neglected gravitational focusing and assumed the relative velocity to be $`\sqrt{2}\sigma _{cl}`$. Thus a galaxy is expected to encounter one other galaxy over the course of its evolution. However, the probability of a merger of two galaxies is much smaller:

$$P_{\mathrm{mer}}N_{\mathrm{enc},1}\left(\frac{\sigma _g}{\sigma _{cl}}\right)^410^3$$
(2)
for the galactic velocity dispersion $`\sigma _g=200\mathrm{km}\mathrm{s}^1`$. Here, three powers of $`\sigma _g/\sigma _{cl}`$ come from the phase space density and another power of $`\sigma _g/\sigma _{cl}`$ comes from the limitation of the velocity that leads to a merger (Binney & Tremaine, 1987).

Also, ellipticals usually have much higher central densities than spiral galaxies. So if two spirals merge, it seems hard to close the gap in density without strong dissipation (Carlberg, 1986). As Jim Gunn once said (Gunn, 1987), “I do not think you can make rocks by merging clouds.” Thus the merger scenario appears unlikely to dominate the dynamics of galaxies, although mergers could be important for the formation of the central giant cD galaxy (Ostriker & Hausman, 1977; Merritt, 1985; Lauer, 1988) and the ultraluminous infrared galaxies and AGNs (e.g., Canalizo & Stockton 2001).

### 2.3 Truncation of Dark Matter Halos

Richstone (1976), White (1976), and Merritt (1983), among others, have investigated the truncation of dark galactic halos by the tidal field of the cluster. In the simplest case, where both the cluster and the galaxy halos have singular isothermal distributions, the tidal radius is approximately

$$R_t\frac{\sigma _g}{\sigma _{cl}}R_p,$$
(3)
where $`R_p`$ is the distance of closest approach to the center of the cluster.

Real clusters are not perfectly smooth. Any substructure left after the hierarchical formation of the cluster will amplify the truncation effect. It is expected to control the amount of dark matter, and therefore, the mass of the galaxies. On the other hand, tidal truncation should not significantly affect the dynamics of the galactic disks which are much more compact. Morphological evolution of luminous galaxies must be due to something else.

### 2.4 Ram-Pressure Stripping of Gas

Gunn & Gott (1972) argued that the ram pressure of the intracluster medium (ICM) might sweep away gas from the galactic disks. X-ray observations show the presence of a significant amount of hot gas at the virial temperature, contributing as much as 10% of the total cluster mass. While the stars move freely through this media, the gas is subject to the ram pressure $`P_{\mathrm{ram}}\rho _{\mathrm{icm}}v^2`$, where $`\rho _{\mathrm{icm}}`$ is the density of the ICM and $`v\sigma _{cl}`$ is the velocity of the galaxy. In order to strip the gas, the ram pressure should exceed the restoring gravitational force per unit area of the disk:

$$\rho _{\mathrm{icm}}v^2>2\pi G\mathrm{\Sigma }_s\mathrm{\Sigma }_g,$$
(4)
where $`\mathrm{\Sigma }_s`$ is the surface density of the stars, and $`\mathrm{\Sigma }_g`$ is the surface density of the gas. For example, if the gas with an average density 1 cm-3 is distributed over a scale height of 100 pc and the stellar density is typical of the solar neighborhood, $`\mathrm{\Sigma }_s75M_{\mathrm{}}`$ pc-2, the ICM needs to be as dense as $`n_{\mathrm{icm}}>4\times 10^4`$ cm-3. Such densities are observed in the cores of rich clusters, and numerical experiments (Kundić et al., 1993; Abadi, Moore & Bower, 1999; Quilis, Moore & Bower, 2000) show that a typical spiral there may lose most of its atomic gas.

Spiral galaxies in clusters have long been known to be deficient of atomic hydrogen (Krumm & Salpeter, 1976; Haynes, Giovanelli & Chincarini, 1984). These “anemic” galaxies (van den Bergh, 1991) have normal molecular hydrogen abundances but lack the atomic component in the outer parts of their disks. Since the H2 gas resides in dense molecular clouds, it is more able to withstand the external pressure; on the other hand, the i gas is usually distributed much more smoothly at lower densities and is susceptible to the ram pressure stripping. But Valluri & Jog (1991) noticed an interesting feature in three rich Abell clusters and the Virgo cluster: the i deficiency is stronger for larger galaxies. Since bigger spirals (excluding the LSBs) usually have a higher surface density, they should be able to protect their gaseous component better. The observed trend is the opposite. Having considered various possible mechanisms of gas removal, Valluri and Jog conclude that the i deficiency can only be explained by tidal interactions.

### 2.5 Tidal Interactions

Tides are a critical factor of the cluster environment. Tidal effects depend on relative sizes of the system and the perturber. Consider stellar systems from smallest to largest scales: stars are $`10^8`$ times smaller than star clusters and experience almost no tidal effects (except in close binaries). Star clusters are $`10^4`$ times smaller than galaxies and experience significant tidal effects, which however operate on a long time scale (Gnedin & Ostriker, 1997). In contrast, galaxies are just 30 to 100 times smaller than clusters of galaxies. There tidal effects ought to be important!

Semi-analytic models of galactic encounters (Merritt, 1983) predict the segregation of galaxies by mass and the formation of a central cD galaxy. The giant cD grows by swallowing other massive galaxies in a process that widens the gap between the first and the second ranking galaxies (Ostriker & Hausman, 1977). Mass segregation in clusters is also supported by numerical simulations (Frenk et al., 1996). Recently, Dubinski (1998) confirmed that the central galaxy, rising by mergers of its neighbors, follows a de Vaucouleurs profile and would be classified as a giant elliptical.

Restricted numerical simulations of Byrd & Valtonen (1990) and Valluri (1993) show that tidal compression by the cluster halo produces spiral arms and tidal tails in some of the representative models of disk galaxies.

More recently, Moore at al. (1996b,1998) simulated the dynamics of galaxies in a fixed, isothermal cluster. They found that fast close encounters with massive galaxies tend to destroy many dwarfs, and even large galaxies lose their thin disks. So the spiral galaxies are essentially transformed into the ellipticals or dwarf spheroidals. However, the potential of a real cluster grows in time, which leads to two competing effects: (1) galaxies enter the cluster at a later epoch and have less time to interact (depending on cosmological model), and (2) the time-varying cluster potential gives rise to large-scale perturbations in the tidal field and causes tidal shocks comparable to those from massive galaxies. New simulations by Moore at al. (1999) include the hierarchical cluster formation and show weaker evolution of galaxies. Whether the galaxy–galaxy or the cluster–galaxy interactions dominate has not yet been understood.

Similar tidal effects operate on smaller scale in the Local Group and may lead to the morphological evolution of the dwarf satellite galaxies (Mayer et al., 2001).

## 3 Numerical Method

I use a Particle-Mesh code (Cen, 1992) to simulate the formation and evolution of three clusters of galaxies in the critical ($`\mathrm{\Omega }_0=1`$), open ($`\mathrm{\Omega }_0=0.4`$), and flat ($`\mathrm{\Omega }_0+\mathrm{\Omega }_\mathrm{\Lambda }=1`$) cosmological models described in Table 1. The code solves Poisson’s equation using FFT on a rectangular grid and evaluates the density field using the Cloud-in-Cell (CIC) method. The code runs fast in parallel, allowing a large simulation with $`256^3`$ particles distributed on a $`512^3`$ grid. Some of the PM results have been tested with the Adaptive Mesh Refinement code (Bryan et al., 1995). Specifically, I checked galaxy trajectories and the calculation of the tidal force.

The comoving size of the simulation box is $`32h^1`$ Mpc and the grid cell size is $`\mathrm{\Delta }=62.5h^1`$ kpc. This is a compromise between the desire to reach a reasonably high resolution in the cluster center and the need to have enough computational volume to include late infall. The latter is important because a very strong initial density perturbation may lead to more mass being bound to the cluster than is available in the whole box, $`9.1\times 10^{15}\mathrm{\Omega }_0h^1M_{\mathrm{}}`$, causing the cluster to “under-form”.

Tests with different amplitudes of the power spectrum, $`\sigma _8`$, show that the dark matter velocity dispersion artificially saturates for $`\sigma _81`$ in the low $`\mathrm{\Omega }_0`$ models. Therefore, I adopt $`\sigma _8=1`$ for these models, in rough agreement with the relation derived from the comparison of large scale N-body simulations with the all-sky survey of X-ray clusters, $`\sigma _80.5\mathrm{\Omega }_0^{0.5}`$ (Eke, Cole & Frenk, 1996; Pen, 1998). Then the virial mass of the simulated clusters is about $`4\times 10^{14}M_{\mathrm{}}`$, or 7% of the total mass within the computational box. The actual amount of mass bound to the clusters may be a factor of 2–3 higher, up to 20% of the total mass. In the critical density model, the normalization $`\sigma _8=0.5`$ is chosen such that the cluster velocity dispersion is the same as in the other two models.

The simulations start at redshift $`z_i=20`$. The time step is determined by the requirement that no particle crosses more than a half of the grid cell, $`0.5\mathrm{\Delta }`$. Also, it limits the expansion parameter $`a(1+z)^1`$ to grow by no more than $`0.05(1+z_i)^1`$. I have checked that increasing the time step by a factor of two does not change noticeably the structure of the resulting cluster.

The look-back time as a function of redshift (Peebles 1993, p. 314) can be calculated analytically in each of the three models. For $`\mathrm{\Omega }_0=1`$, the standard result is $`H_0t(z)=\frac{2}{3}(1+z)^{3/2}`$. For $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$,

$$H_0t(z)=\frac{(1+\mathrm{\Omega }_0z)^{1/2}}{(1\mathrm{\Omega }_0)(1+z)}+\frac{\mathrm{\Omega }_0}{2(1\mathrm{\Omega }_0)^{3/2}}\mathrm{ln}\frac{f^{1/2}1}{f^{1/2}+1},f\frac{1+\mathrm{\Omega }_0z}{1\mathrm{\Omega }_0}.$$
(5)
For $`\mathrm{\Omega }_0+\mathrm{\Omega }_\mathrm{\Lambda }=1`$,

$$H_0t(z)=\frac{2}{3(1\mathrm{\Omega }_0)^{1/2}}\mathrm{ln}\frac{1+(1+f)^{1/2}}{f^{1/2}},f\frac{\mathrm{\Omega }_0}{1\mathrm{\Omega }_0}(1+z)^3.$$
(6)
The Hubble constant is assumed to be $`H_0=65`$ km s-1 Mpc-1 ($`h=0.65`$), according to the time delay measurements in a gravitational lens system 0957+561 (Kundić et al., 1997) and other independent methods (Freedman et al., 1994).

On 16 processors of SGI Origin 2000, the full simulation with a $`512^3`$ grid runs for about 36 hours. A large fraction of the time is spent on the calculation of the tidal field due to and around the identified halos. All simulations presented in this paper were done using the 64-processor SGI Origin 2000 supercomputer at the Princeton University Observatory.

### 3.1 Constrained Initial Conditions

The cluster simulations use the constrained Gaussian random initial conditions, generated using the method of Hoffman & Ribak (1991). The actual realization is taken from the cluster comparison project of Frenk et al. (1999). In that work the initial density perturbations were chosen to produce a $`3\sigma `$ peak with a Gaussian filter of radius $`5h^1`$ Mpc. Retaining the same random phases from the original realization, I have renormalized the perturbations according to the power spectrum of the chosen cosmological models.

For each of the three models, I evaluate the power spectrum $`P(k)_{\mathrm{new}}`$, including the effects of baryons, using the fitting expressions of Bardeen et al. (1986) and Hu & Sugiyama (1996). On the cluster and sub-cluster scales, the power spectra follow approximately a power law $`P(k)k^{2.5}`$. Then I transform the initial displacements $`d`$ into the Fourier space and renormalize their amplitudes according to the new power spectrum:

$$d_k^{\mathrm{new}}=d_k^{\mathrm{old}}\left(\frac{P(k)_{\mathrm{new}}}{P(k)_{\mathrm{old}}}\right)^{1/2},$$
(7)
where the old power spectrum, $`P(k)_{\mathrm{old}}`$, is from Frenk et al. (1999). Both power spectra are normalized to the corresponding value of $`\sigma _8^2(z_i)`$ at the initial redshift $`z_i`$. Finally, I transform the perturbation field back to real space and obtain the initial displacements for the dark matter particles.

### 3.2 Structure of the Clusters

The parameters of the resulting clusters are given in Table 2. The virialized region is taken to be a sphere originating at the cluster center, with a mean overdensity $`\delta =200`$. The average mass of the clusters is $`4\times 10^{14}M_{\mathrm{}}`$, and the one-dimensional velocity dispersion is 660 km s-1.

Figure 1 shows the projected surface density of the three clusters at the present epoch. The $`\mathrm{\Omega }_0=1`$ cluster has not yet relaxed and still continues to accrete matter. Note a large infalling group of galaxies to the upper right from the center on the enlarged panel. In contrast, the two low $`\mathrm{\Omega }_0`$ clusters appear more spherical and regular, and also very similar to each other; the cosmological constant does not seem to affect the cluster properties.

Figure 2 shows that the density profiles of the three clusters are very similar, with a slight excursion of the $`\mathrm{\Omega }_0=1`$ cluster at $`r2h^1`$ Mpc due to the infalling group. This group is clearly visible on the velocity dispersion plot, Figure 3.

The decrease of the velocity dispersion at small distances from the center is partly real and partly numerical. The high resolution dark matter simulations in Frenk et al. (1999) show a similar decrease; therefore our models should exhibit the trend. Also, such a decrease exists in the popular NFW model (Navarro, Frenk & White, 1997). On the other hand, present simulations lack resolution in the cluster core and cannot be trusted at small radii. To illustrate the resolution effects, I have repeated the simulations of Cluster III on smaller grids, $`256^3`$ and $`128^3`$ cells (Table 3), with the correspondingly fewer particles, $`N_{\mathrm{grid}}=2^3N_\mathrm{p}`$. Figure 4 shows that the density profile is reproduced well to the resolution limit of each simulation, but the velocity dispersion starts to decline much earlier in the lower resolution runs. Nevertheless, the difference between the $`512^3`$ and the $`256^3`$ runs is small and confirms the accuracy of the larger simulation.

The simulated clusters are similar to the Virgo cluster in size and velocity dispersion. ROSAT results (Böhringer, 1994) indicate that Virgo’s X-ray mass is $`(1.55.5)\times 10^{14}M_{\mathrm{}}`$, while the optical velocity dispersion of galaxies is between 581 and 721 km s-1 (Huchra, 1985). These values are quite common for the nearby clusters, as shown by the optical and X-ray survey of Girardi et al. (1998). Often cited in theoretical work, the Coma cluster is much larger, with the one-dimensional velocity dispersion about 1000 km s-1 (Merritt & Gebhardt, 1994) and the virial mass of order $`2\times 10^{15}M_{\mathrm{}}`$ (Böhringer, 1994). Our models, therefore, represent modest but more common clusters of galaxies. Their velocity dispersion is characteristic of the Abell richness class 1 (see Postman, Huchra & Geller 1992).

## 4 Galaxies in the Simulation

The definition of a galaxy in the cluster simulation is not trivial. While at high redshift dark matter halos are significantly denser than the surrounding media, at low redshift the density contrast effectively disappears. This effect is known as overmerging of halos (e.g., Moore, Katz & Lake 1996a; Klypin et al. 1999) due to the insufficient mass resolution. Therefore, it is difficult and, to some extent, uncertain to try to identify the galaxy halos from the particle distribution at all times. Instead, I have implemented the following simple procedure, which can be extended to a higher degree of sophistication.

I identify the dark matter halos at a time early enough that the density contrast is large, but late enough that galaxies are known to exist. A fiducial epoch at which the central parts of large galaxies are assumed to have formed is taken to be $`z_g=5`$.

In order to obtain a statistically significant sample, I choose 100 most massive halos identified by the group finding algorithm. Then, I add 100 massless particles to the PM simulation with the coordinates and velocities of the halo centers of mass. These new “galactic” particles do not affect either the density or the potential calculation but rather act as test particles. I follow the motion of these particles and output the tidal field along their trajectories at each time step.

By not using real particles to calculate galaxy trajectories, we effectively avoid the overmerging problem: the test particles will not merge even if they come very close to each other. The drawback is that due to a finite time step the test particles may gradually depart from the bulk of the galactic particles. To guard against this effect, I have also implemented another scheme.

The second scheme involves real particles constituting the halo. I calculate the halo density at each particle’s position using the same CIC kernel that is used in the main PM calculation. Then I take the densest subset of particles, 1/8 of the total, to calculate the center of mass directly. These dense particles presumably lie deep in the potential well of the halo and serve as robust indicators of the center. In case of the two galaxies merging, their centers defined this way would clearly indicate it. However, I do not attempt to resolve the internal structure of galaxies in these simulations and only trace their centers of mass.

### 4.1 Group Finding Algorithm

Dark matter halos are identified using the coordinate-free algorithm HOP proposed by Eisenstein & Hut (1998). It is similar in spirit to DENMAX (Gelb & Bertschinger, 1994) but uses densities at the particle positions instead of the regular grid. I have extended HOP by including a procedure for the removal of particles which are not gravitationally bound to the group. This procedure, similar to the SKID algorithm (Weinberg, Hernquist & Katz, 1997), is outlined below.

The gravitational energy of the particles in a group is calculated by direct summation of the Plummer potentials with a smoothing length corresponding to one grid cell of the cluster simulation. The kinetic energy includes the Hubble expansion term. The total energy per unit mass is

$$E_i=\frac{1}{2}\left[\text{v}_i\text{v}_{\mathrm{com}}+H_g(\text{r}_i\text{r}_{\mathrm{com}})\right]^2\underset{j}{}\frac{Gm_p}{((\text{r}_i\text{r}_j)^2+ϵ^2)^{1/2}},$$
(8)
where $`m_p`$, $`\text{r}_i`$, $`\text{v}_i`$ are the particle mass, physical coordinate, and velocity, respectively, and $`H_g`$ is the Hubble constant at redshift $`z_g=5`$. The index “com” refers to the center of mass of the group. The particles are unbound as follows. I find particle $`m`$ with the maximum energy and discard it, if $`E_m>0`$. I recalculate the center of mass using the remaining particles and update the potential energy. Then, I find a new highest energy particle and repeat the procedure. In the end, all remaining particles are bound or else the group is dispersed.

A straightforward implementation of the unbinding procedure works well, although somewhat slowly for groups containing more than $`10^4`$ particles. In the longest case, the unbinding algorithm took about 40 minutes on a single SGI R10000 processor.

The second part of the HOP algorithm, REGROUP, uses three different overdensity parameters to define the outer boundary of the group ($`\delta _{outer}`$), the saddle point for merging two groups together ($`\delta _{saddle}`$), and the minimum density of a viable group ($`\delta _{peak}`$). After some experimenting, I have chosen the following parameters to produce the reasonable-size halos: \[$`\delta _{outer},\delta _{saddle},\delta _{peak}`$\] =  for the $`\mathrm{\Omega }_0=1`$ model, and  for the low $`\mathrm{\Omega }_0`$ models. The different sets of parameters are necessary because of the different masses of particles ($`m_p=5.4\times 10^8\mathrm{\Omega }_0h^1M_{\mathrm{}}`$).

Figures 57 illustrate the location of particles in the identified halos. Notice the filamentary structure of bound particles indicating the directions along which the halos continue to assemble. Most galaxies in the central part of the simulation volume (shown on the Figures) fall eventually into the forming cluster. Because of the difference in particle masses, the $`\mathrm{\Omega }_0=1`$ halos contain significantly fewer particles than the $`\mathrm{\Omega }_0=0.4`$ halos.

### 4.2 Halo Mass Function

The mass function of the identified halos can be compared with the prediction of the Press-Schechter theory (Press & Schechter, 1974). I have calculated the expected number density of galaxies at $`z_g`$ using equations (A4-A9) of Mo, Mao & White (1998).

Figure 8 shows the halo mass functions in the three clusters. Note the relative depletion of low mass halos in the simulations. This is partly due to the truncation of the input power spectrum at the Nyquist frequency, which effectively prevents the formation of small-scale objects. For comparison, the lower right panel shows the mass function obtained in the lower resolution, $`256^3`$ simulation of Cluster III. As expected, the smaller halos are depleted more strongly than the larger ones.

### 4.3 Galaxy Trajectories

Figures 911 show the trajectories of two selected halos in each cluster. These galaxies are used for the high resolution simulations in Paper II, representing a large spiral and a dwarf spheroidal, respectively. Their trajectories are quite irregular. In the $`\mathrm{\Omega }_0=1`$ case, the galaxies are just entering the cluster at the present epoch, although not on a free-fall path. The orbits are similar to those shown in Frenk et al. (1996). On the other hand, in the low $`\mathrm{\Omega }_0`$ cases the galaxies have enough time to move around the cluster center several times, more so in the open model. These differences will have major effects on the dynamical evolution of the galaxies.

In general, there is a good agreement between the trajectories calculated using two different methods, especially for the less massive galaxies. For the very massive halos, the “test particle” motion occasionally tends to oscillate around the trajectory determined by the densest particle subset. I have therefore decided to use the latter as the more reliable trajectories.

## 5 Mergers and Close Encounters

The hierarchical formation of clusters and the survival of substructure lead to more frequent collisions of galaxies and lower relative velocities at early epochs. Preferentially radial infall implies that the final distribution of galaxies is very dense. Most of the galaxies lie close to the cluster center, while the rest form a scarcely populated halo.

The number density of massive galaxies ($`M_{\mathrm{min}}10^{11}M_{\mathrm{}}`$) can be fitted by a power law with a central core:

$$n_g(r)=\frac{n_0}{(r+r_0)^\gamma },$$
(9)
normalized to the total number of galaxies $`N_g`$ within the virial radius of the cluster, $`R_{\mathrm{vir}}`$. Not all initially identified halos end up in the clusters, so with $`N_g<100`$ trying to estimate the core radius $`r_0`$ and the slope $`\gamma `$ by fitting a histogram of galaxy positions yields highly uncertain results. Instead, I use a maximum likelihood (ML) method analogous to that in Gnedin (1997). The parameters $`r_0`$ and $`\gamma `$ are to be determined by searching for the maximum of the ML function, $`=e^{N_g}_{i=1}^{N_g}4\pi r_i^2n(r_i)`$. Requiring the two partial derivatives of $`\mathrm{log}(r_0,\gamma )`$ to vanish gives a system of two integral equations, which have been solved numerically using the software package Maple. The resulting parameters are given in Table 4.

The slopes of the density profile are steep, $`\gamma =2.44.2`$. The massive galaxies are more concentrated in the inner parts of the clusters than the dark matter which follows an approximately isothermal profile, $`\rho _{dm}r^2`$. Also, the core radii in all three clusters are very small, below the force resolution scale. At such small distances from the center the distribution of galaxies can be affected by the giant cD galaxy which may swallow its neighbors. To check this, I have repeated the ML analysis excluding all galaxies within one resolution element, $`\mathrm{\Delta }`$, from the center. This forces the core radii to be at least one $`\mathrm{\Delta }`$. The resulting values (Table 4) are very close to this limit, suggesting that the true core radii can indeed be quite small. Thus the effect of the central giant on the overall galaxy distribution seems to be small.

### 5.1 Number of Close Encounters

The detailed trajectories of galaxies can be used to calculate the number of close encounters during the epoch of cluster formation. In order to detect all encounters, the trajectories are linearly interpolated between the successive time steps. The comoving separations are converted into physical kiloparsecs using the appropriate value of the Hubble constant at each time step. If a galaxy is found to have merged (see Appendix A), its subsequent trajectory is ignored.

Table 5 shows the number of encounters, $`N_{\mathrm{enc}}`$, as a function of the impact parameter, $`R_{\mathrm{enc}}`$. It is defined as a minimum distance between the centers of mass of the galaxies. The last column gives the “total” number of encounters for $`R_{\mathrm{enc}}=40`$ kpc (it is assumed that no mergers are possible for the separations larger than 40 kpc). Table 5 also shows the average number of close encounters per galaxy, $`N_{\mathrm{enc},1}`$. Since each encounter involves two galaxies, this number is $`N_{\mathrm{enc},1}2N_{\mathrm{enc}}/N_g`$.

Figure 12 shows the number of collisions per galaxy within $`R_{\mathrm{enc}}=20`$ kpc. Note that in the low $`\mathrm{\Omega }_0`$ clusters, most collisions take place near the end of the simulation, when the galaxies have assembled inside the cluster. In Cluster I, the galaxies are closer to each other initially and have close encounters uniformly distributed in time. Several distinct peaks, notably in Clusters I and III, can be due to the new groups of galaxies entering the clusters.

The number of close encounters does not scale very well with $`R_{\mathrm{enc}}^2`$, as predicted by equation (1). Since most encounters take place in the already formed cluster, any mergers prior to that remove prospective galaxies and reduce the number of subsequent collisions. To test that equation (1) would have been valid in the absence of mergers, I count the number of encounters disallowing mergers. The new numbers, $`\stackrel{~}{N}_{\mathrm{enc}}`$, scale as $`R_{\mathrm{enc}}^2`$ in agreement with the prediction (see Table 5).

However, the total number of encounters is much higher than expected. The effect is mainly due to the concentration of galaxies at the center of the cluster, implying a smaller effective volume. Equation (1) can be rewritten in a general form:

$$N_{\mathrm{enc}}=\pi R_g^2V_{\mathrm{enc}}t_Hn_g^2(r)𝑑V,$$
(10)
where the integration extends over the virial volume of the cluster, $`V`$. The mean value of the number density squared, $`n_g^2`$, is substantially higher than the square of the mean density, $`n_g^2`$, because of the galaxy clustering. Their ratio is the clumping factor

$$𝒞\frac{n_g^2}{n_g^2}=\frac{Vn_g^2(r)𝑑V}{\left(n_g(r)𝑑V\right)^2}.$$
(11)
Using the azimuthally-averaged density profile (9) and the auxiliary functions

$$f_n(r_0,\gamma )_0^1\frac{x^2dx}{(x+r_0/R_{\mathrm{vir}})^{n\gamma }},$$
(12)
I find $`𝒞=f_0f_2/f_1^2`$ (see Table 4). Now, the number of close encounters is

$$N_{\mathrm{enc}}=320\left(\frac{𝒞}{1000}\right)\left(\frac{N_g}{65}\right)^2\left(\frac{R_{cl}}{2\mathrm{Mpc}}\right)^3\left(\frac{R_{\mathrm{enc}}}{10\mathrm{kpc}}\right)^2\left(\frac{V_{\mathrm{enc}}}{800\mathrm{km}\mathrm{s}^1}\right).$$
(13)
This is much higher than the naive estimate (1) and agrees well with the results for Cluster III.

Because of the higher degree of clustering in the low $`\mathrm{\Omega }_0`$ models, more than half of all encounters occur within the inner 0.4 Mpc from the center. In contrast, in the $`\mathrm{\Omega }_0=1`$ model collisions take place throughout the cluster, leading to a smaller value of the clumping factor.

### 5.2 Relative Velocities

Figure 13 shows the relative velocities of galaxies at the time of collision. The distributions are skewed, peaking at low values and extending large tails into the high velocities. A characteristic value, the median, is around 800 km s-1. Again, there is a strong difference between Cluster I and the low $`\mathrm{\Omega }_0`$ clusters. Table 6 demonstrates that while the peak velocity is higher in Cluster I, the median and mean values are substantially lower than those in Clusters II and III. The very high velocities are almost unattainable in the $`\mathrm{\Omega }_0=1`$ model.

Table 7 shows that the velocity dispersion of galaxies is smaller than that of dark matter and varies strongly with the position in the cluster. Due to a small number of identified galaxies I have split them in three equally-populated radial bins. In Clusters II and III, the velocity gradient is strong and inverted, with the lowest dispersion within the inner cluster core. In Cluster I the effect is less conspicuous. Since the samples are small, Table 7 serves only to indicate the trend. In contrast, “test particle” galaxies tracing the dark matter potential show a normal gradient weakly rising towards the central bin (from 550 km s-1 to 750 km s-1 in Cluster II). Thus the inverted gradient is a property of only real galaxies.

A review of the literature reveals that numerical simulations often produce “galaxies” with the pairwise velocity dispersion lower than that of dark matter (Colín, Klypin & Kravtsov 2000; Summers, Davis & Evrard 1995; Carlberg 1994; but see also Okamoto & Habe 1999; Diaferio et al. 1999; Ghigna et al. 2000). In a specific study of the relative velocities at the time of galaxy collisions, Tormen, Diaferio & Syer (1998) find that the peak velocity is only $`V_{\mathrm{rel}}2.5\sigma _g500`$ km s-1 (their Fig. 7). There are several reason for this: (i) A statistical bias of hierarchical structure formation, whereas galaxies form around highly-clustered density peaks. Therefore, large galaxies are always around the cluster core and do not gain velocity from gravitational infall. (ii) A physical bias relating to gas cooling and star formation (Cen & Ostriker, 2000) makes younger galaxies move more slowly than the older ones. (iii) Dynamical friction reduces the angular momentum of large galaxies and draws them towards the center. In the present simulations, the final distances from the center are anti-correlated with galaxy masses (Spearman’s correlation coefficient is $`r_s0.3`$ for Clusters II and III, but only $`r_s0.1`$ for Cluster I). Thus there is a moderate amount of mass segregation in the clusters formed at earlier times.

Finally, numerical effects may lead to low galactic velocities. The PM code loses force resolution in the two inner grid cells, within which the galaxies do not “feel” each other’s gravity. Also, if the halo particles departed from its center, the center-of-mass velocity would be significantly reduced. However, if they were scattered completely randomly in the cluster the average velocities would have fallen below 100 km s-1, which is clearly ruled out.

At present, observations neither confirm nor reject the velocity bias. Carlberg, Yee & Ellingson (1997), in a survey of 1150 galaxies from 14 clusters, find an almost constant projected velocity dispersion up to the virial radius. The line-of-sight projection may, however, wash out a decrease of the true velocity dispersion near the center. Another large survey of 15 rich clusters (Dressler & Shectman, 1988) shows a significant amount of substructure with the lower pairwise velocities of galaxies, as well as a modest decline (of the order 20%) of the azimuthally averaged dispersion in the central regions. So the current observational status of the velocity bias is inconclusive.

### 5.3 Merger History

The low relative velocities of galaxies make mergers more likely. Without having enough resolution to study mergers in detail, I adopt a simple prescription based on the impact characteristics of the galaxies. It consists of the three merger criteria described in Appendix A.

Figure 12 and Table 5 show the expected number of mergers per galaxy, $`N_{\mathrm{mer},1}`$. It is calculated as $`N_{\mathrm{mer}}/N_g`$, since each merger removes only one galaxy. The calculated merger probability is high. At impact parameters $`R_{\mathrm{enc}}<10`$ kpc it is already $`N_{\mathrm{mer},1}10\%`$, and the total value reaches $`20\%`$. Most of the mergers take place shortly after galaxies enter the clusters (even earlier in Cluster I). Once the clusters virialize, galactic velocities become too high and prohibit mergers. Yet, the number of mergers is only a very small fraction, about 0.1%, of the number of close encounters (cf eq. ).

Ghigna et al. (1998) find very few mergers within the virial radius of their high-resolution cluster. This can be reconciled with our results because Ghigna et al. (1998) count only the recent encounters in the virialized cluster ($`z<0.5`$) while our numbers include all encounters starting from $`z_g=5`$. Note that most mergers occur at relatively early times (Figure 12). In the three cases considered, the last merger occurred at 2.3, 1.0, and 3.3 Gyr before the present, respectively. The merger probability of 22–29% may reflect the instances of the formation of large elliptical galaxies at high redshifts. They are expected to take place naturally in the centers of future clusters.

## 6 Tidal Field Around Galaxies

The tidal field of the cluster is calculated along the trajectories of the identified galaxies. The tidal force $`F_{\mathrm{tid},\alpha }`$ and the shear tensor $`F_{\alpha \beta }`$ are

$$F_{\mathrm{tid},\alpha }\left(\frac{d^2\mathrm{\Phi }}{dr_\alpha dr_\beta }\right)_0r_\beta F_{\alpha \beta }r_\beta ,$$
(14)
where $`\mathrm{\Phi }`$ is the external potential and r is the radius-vector in the galactic reference frame. By definition tidal force is zero at the center of mass of the galaxy. I evaluate the tidal tensor at each time step for each identified halo, using the tri-linear interpolation scheme with a smoothing filter described in Appendix B.

Being coordinate-independent, the trace of the tidal tensor serves as a good measure of the tidal field. Figures 1416 show the trace $`F_{\alpha \alpha }`$ around the pairs of large and dwarf galaxies chosen in §4.3 from each cluster simulation. The time variation of the tidal force is extremely irregular. There are many strong peaks, reaching the amplitude of 200 Gyr-2 and higher. The largest peaks are usually narrow, with the characteristic duration of $`10^8`$ yr, while the smaller peaks last longer. The occurrence of the two types of peaks is clearly differentiated in time: the larger happen early, when massive galaxies are distinctly denser than the background; the smaller happen at later stages, when the halos merge and the gravitational potential becomes smoother.

The amplitude of the tidal field is the largest in the $`\mathrm{\Omega }_0=0.4`$ simulation, although strong perturbations continue for a longer time in the $`\mathrm{\Lambda }`$ model. In the $`\mathrm{\Omega }_0=1`$ model, the overall amplitude is smaller but there is more temporal variation. The main distinction between the galaxies in the low and high $`\mathrm{\Omega }_0`$ models is the time spent within the cluster. In the latter case, the galaxies are just arriving to the center, while in the low $`\mathrm{\Omega }_0`$ models the galaxies have entered earlier and experienced more tidal interaction.

The overall amplitude of the tidal force is determined by the epoch of cluster formation. For the same virial overdensity, the cluster is denser on average when the mean density of the Universe is high. In the low $`\mathrm{\Omega }_0`$ models clusters form at a higher redshift and therefore exert stronger tidal forces than in the $`\mathrm{\Omega }_0=1`$ model. Since the subsequent growth of the clusters is balanced by the drop of the mean density, the tidal force is expected to be of the same order of magnitude throughout the evolution.

Note that the peaks of the tidal force are not always correlated with closest approach to the cluster center. The galaxies are being tidally “shocked” at 1 Mpc from the center as often as at 400 kpc. Even though the center of the cluster has been moving slowly in the simulation box, it is still likely that close encounters with massive galaxies or groups of galaxies happen relatively far from the center. As the cluster grows, such tidal encounters occur within a fairly large (comoving) volume.

The amplitude of the tidal force tells us about the effective scale of perturbation. For a smooth cluster, the tidal tensor is of the order

$$F_{\alpha \beta }\frac{GM(r)}{r^3}\frac{\sigma _{cl}^2}{r^2},$$
(15)
where $`M(r)`$ is the enclosed mass within radius $`r`$. In the simulated clusters the three-dimensional velocity dispersion is $`\sigma _{cl}=\sqrt{3}\times 660`$ km s-1. If a characteristic scale, the cluster core radius, is 200 kpc then the expected tidal perturbation is

$$F_{\alpha \beta }32\left(\frac{\sigma _{cl}}{1100\mathrm{km}\mathrm{s}^1}\right)^2\left(\frac{r}{200\mathrm{kpc}}\right)^2\mathrm{Gyr}^2.$$
(16)
The measured amplitude of the tidal force (cf. Figures 1416) is about 150 Gyr-2. Turning the argument around, such tidal force corresponds to the effective scale of perturbation

$$r_p=92\left(\frac{\sigma _{cl}}{1100\mathrm{km}\mathrm{s}^1}\right)\left(\frac{F_{\alpha \beta }}{150\mathrm{Gyr}^2}\right)^{1/2}\mathrm{kpc}.$$
(17)
Smaller than the core radius but larger than individual galaxies, this scale is appropriate for small groups of galaxies. As discussed in §2.1, groups maintain their identity in many clusters. Even such a relaxed cluster as Coma shows a detectable amount of substructure (Colless & Dunn, 1996; Conselice & Gallagher, 1999). Since large peaks of the tidal force dominate the total amount of heating (see §7), I conclude that the tidal effects are produced primarily by the interaction of galaxies with the remaining substructure in clusters.

## 7 Tidal Heating

The main effect of external tidal forces on a galaxy is a positive energy change of stars and dark matter, or tidal heating. The amount of tidal heating can be estimated analytically, at least for the spherically-symmetric component. For a single tidal shock, the energy input is the ensemble average of the velocity change squared, $`\frac{1}{2}(\mathrm{\Delta }v)^2`$, where $`\mathrm{\Delta }v`$ is a time integral of the tidal force times the reduction factor that takes into account the conservation of adiabatic invariants of stellar orbits (Gnedin, Hernquist & Ostriker, 1999). In the simulations, the tidal force has many peaks of variable duration. In addition, different components of the tidal tensor are not perfectly correlated so that they experience maxima and minima at different times. Clearly, this complex behavior cannot be fully described as simply as a single tidal perturbation. Instead, I design a cumulative parameter, based on the semi-analytic theory of tidal shocks, which provides a useful scaling relation for the amount of tidal heating.

Assume that each peak of the tidal force can be considered a separate “tidal shock”. The total amount of heating is then the sum of the contributions from all peaks. Moreover, since the velocity changes add in quadrature, all components of the tidal tensor can be treated separately. Then it follows from the semi-analytic theory (cf. Gnedin & Ostriker 1999) that the total energy per unit mass of the system changes by

$$\mathrm{\Delta }E=\frac{1}{6}I_{\mathrm{tid}}r^2,$$
(18)
where $`r`$ is a characteristic size of the galaxy and $`I_{\mathrm{tid}}`$ is the tidal heating parameter

$$I_{\mathrm{tid}}\underset{n}{}\underset{\alpha ,\beta }{}\left(F_{\alpha \beta }𝑑t\right)_n^2\left(1+\frac{\tau _n^2}{t_{\mathrm{dyn}}^2}\right)^{3/2},$$
(19)
where the sum extends over all peaks $`n`$ and all components of the tidal tensor, $`\alpha ,\beta =\{x,y,z\}`$. Here $`\tau _n`$ is the effective duration of peak $`n`$ for each value of $`\alpha `$ and $`\beta `$, and $`t_{\mathrm{dyn}}`$ is the half-mass dynamical time. For the galaxy model from Appendix A, the dynamical time scales with the halo mass as $`t_{\mathrm{dyn}}M_g^{1/4}`$. The inner integration is limited to the time intervals between the successive minima of the tidal force. In carrying out the integration, I have linearly interpolated $`F_{\alpha \beta }`$ between the output points to match the time step of the galactic simulations presented in Paper II. This provides the most direct way to compare the analytical estimate $`I_{\mathrm{tid}}`$ with the results of self-consistent N-body simulations. The agreement is generally good.

Figure 17 shows the distribution of parameters $`I_{\mathrm{tid}}`$ for the three clusters. There are more galaxies with the low value of $`I_{\mathrm{tid}}`$ than with the high – tidal effects are different for different galaxies. The histograms display a distinct trend: the low $`\mathrm{\Omega }_0`$ clusters have stronger tidal heating than the $`\mathrm{\Omega }_0=1`$ cluster. A least-squares fit to the histograms, using a convenient functional form $`N(I_{\mathrm{tid}})=N_0\mathrm{exp}(I_{\mathrm{tid}}/I_0)`$, yields the following characteristic parameters:

$$\begin{array}{cc}I_0=8\times 10^2\text{Gyr}^2\hfill & \text{(Cluster I),}\hfill \\ I_0=5\times 10^3\text{Gyr}^2\hfill & \text{(Cluster II),}\hfill \\ I_0=2\times 10^3\text{Gyr}^2\hfill & \text{(Cluster III).}\hfill \end{array}$$
(20)
The timescale of perturbation, $`\tau _n`$, is very important in determining the total amount of heating. Without the adiabatic correction factor, the values of $`I_0`$ would have been at least an order of magnitude higher. Overall, the open $`\mathrm{\Omega }_0=0.4`$ model has the highest heating rates and the $`\mathrm{\Omega }_0=1`$ model the lowest. Thus galaxies should experience considerably more dynamical evolution in a low density universe than in the $`\mathrm{\Omega }_0=1`$ case.

A quick estimate whether tidal heating can turn a disk galaxy into an elliptical is the increase of the velocity dispersion induced by the energy change (eq. ): $`\frac{1}{2}\sigma _{\mathrm{tid}}^2\mathrm{\Delta }E`$. For a typical tidal parameter $`I_{\mathrm{tid}}=10^3`$ Gyr-2, the inner region at 3 kpc will be heated to $`\sigma _{\mathrm{tid}}=54`$ km s-1. This is not enough to completely randomize stellar velocities in a large galaxy, but enough to turn a thin disk into a thick one (cf. §7.2). In order to reach the velocity dispersion of a large elliptical galaxy, 250 km s-1, the tidal heating parameter would need to be $`2.2\times 10^4`$ Gyr-2. This is above the maximum for our clusters, and therefore the ellipticals are unlikely to be created by tidal interactions.

### 7.1 Tidal Effects from Close Encounters

The smoothing filter used in the calculation of the tidal field (Appendix B) automatically removes the effect of close galactic encounters. Their additional contribution to the tidal heating rate can be estimated as follows. The number of close encounters per galaxy in the simulation as a function of impact parameter $`R=R_{\mathrm{enc}}`$ is (cf. eq. )

$$N_{\mathrm{enc},1}10\left(\frac{R}{10\text{kpc}}\right)^2\left(\frac{V_{\mathrm{enc}}}{800\mathrm{km}\mathrm{s}^1}\right).$$
(21)
The amplitude of the tidal force due to a galaxy of mass $`M_g`$ is of the order $`F_{\alpha \beta }GM_g(R)/R^3`$. The duration of the encounter is approximately $`R/V_{\mathrm{enc}}`$, so that the contribution to the heating parameter at a distance $`R`$ is

$$I_{\mathrm{enc}}(R)=\left(F_{\alpha \beta }𝑑t\right)^2\left(\frac{GM_g(R)}{R^2V_{\mathrm{enc}}}\right)^2.$$
(22)
The total contribution of all encounters is the integral over all impact parameters weighted by the number of encounters, $`N_{\mathrm{enc},1}`$. The perturbing halo can be represented by a truncated isothermal sphere with the density profile given by equation (A1). The lower limit of integration is set by the condition $`N_{\mathrm{enc},1}>1`$, i.e. that at least one encounter takes place; this gives $`R_{\mathrm{min}}3`$ kpc. The upper limit can be taken as the tidal radius, $`R_{\mathrm{max}}30`$ kpc, but in the end the heating rate depends only logarithmically on the integration limits.

$`I_{\mathrm{enc}}`$ $`=`$ $`{\displaystyle _{R_{\mathrm{min}}}^{R_{\mathrm{max}}}}I_{\mathrm{enc}}(R)𝑑N_{\mathrm{enc},1}(R)`$ (23)
$``$ $`4.9\times 10^2\left({\displaystyle \frac{\sigma _g}{200\text{km s}^1}}\right)^4\left({\displaystyle \frac{V_{\mathrm{enc}}}{800\text{km s}^1}}\right)^1\text{Gyr}^2,`$
where $`\sigma _g`$ is the galactic velocity dispersion. This value of $`I_{\mathrm{enc}}`$ is comparable to, but smaller than, the characteristic parameter $`I_0`$. Thus close encounters contribute an important part but do not dominate the tidal heating.

### 7.2 Self-consistent N-body simulations

These analytical estimates will be substantiated in Paper II (Gnedin, 2003) which presents self-consistent N-body simulations of several galaxies selected from each cluster sample. Each galaxy is resimulated with $`2\times 10^6`$ particles, supplementing the internal dynamics with the external tidal forces of the cluster. Disks of large spiral galaxies survive tidal heating but thicken in the vertical direction by a factor of two, while their dark matter halos are truncated at roughly 30 kpc for $`\sigma _g=250`$ km s-1. At the end of the simulation they resemble featureless S0 galaxies. In contrast, low surface brightness galaxies are almost completely disrupted by tidal forces.

Most of the results can be explained qualitatively by the critical tidal density, corresponding to the trace of the tidal tensor via Poisson’s equation: $`4\pi G\rho _{\mathrm{tid}}F_{\alpha \alpha }`$. It is similar to the instantaneous density of the Roche lobe, but reflects the highest tidal peaks along the galactic orbit. This density $`\rho _{\mathrm{tid}}`$ determines the truncation radius of the halos and the amount of luminous and dark mass stripped by tidal heating.

Galactic N-body simulations confirm that the amount of tidal heating is set by $`I_{\mathrm{tid}}`$ (or $`\rho _{\mathrm{tid}}`$). The evolution of galaxies is stronger on average in Cluster II. However, reflecting the wide distribution of tidal parameters in each cluster (Fig. 17), the heating on different orbits within the same cluster varies as much as between the different cluster models.

## 8 Conclusions

Tidal interactions may explain the observed dynamical and morphological evolution of galaxies in clusters. Using the simulations of hierarchical cluster formation, I calculate the tidal forces along the galactic trajectories starting at redshift $`z_g=5`$. The tidal field shows a lot of temporal variation. The peaks of the tidal force do not always correlate with the closest distance to the cluster center and are likely to be due to the infalling groups of galaxies or individual massive galaxies. The tidal peaks are stronger and more impulsive in the low $`\mathrm{\Omega }_0`$ models. Also, in those models the cluster forms at earlier epochs and the galaxies experience more interaction. Thus overall tidal effects are stronger in the low $`\mathrm{\Omega }_0`$ clusters.

Massive galaxies are strongly clustered at the end of the simulations. Together with the remaining substructure, this significantly increases the probability of galactic collisions. The number of encounters within 10 kpc is about 10 per galaxy in the Hubble time. About 0.1% of such encounters is expected to result in a merger, defined by the analytical merger criteria in §A. This gives a high merger probability of the identified halos, 22–29% per halo, although most mergers take place before the clusters virialize.

The amount of dynamical evolution of galaxies induced by tidal interactions can be estimated using the tidal heating parameter, $`I_{\mathrm{tid}}`$. The typical values of $`I_{\mathrm{tid}}10^3`$ Gyr-2 result in an increase of the stellar velocity dispersion by $`50`$ km s-1. At present this is insufficient to turn a large disk galaxy into an elliptical galaxy. However, in the future a further secular evolution and transformation into the ellipticals seems inevitable. Close galactic encounters cannot be resolved in present simulations, but they are estimated to contribute a fraction of 10% to 50% of the total heating rate.

I would like to thank my thesis advisor J. P. Ostriker for his support and guidance, Nick Gnedin for expertise in cosmology, Renuye Cen for his PM code, Greg Bryan for sharing his AMR code, Mike Blanton for his excellent Points visualization software, and Jeremy Goodman, David Spergel, Neta Bahcall, and Scott Tremaine for valuable comments. I acknowledge the support from NSF grants AST 93-18185 and AST 94-24416. This work was submitted in partial fulfillment of the Ph.D. requirements at Princeton University.

## Appendix A Merger Criteria and Merging Time

This appendix describes the merger criteria used in §5.3. The physical processes leading to merger are the gravitational capture of two galaxies and the subsequent dynamical friction of their halos. The energy exchange during the first encounter determines whether the galaxies form a bound pair, while dynamical friction determines the timescale for final coalescence. Also, the tidal forces of the cluster can disrupt even a bound system at its point of maximum separation.

For clarity, I model the halos as truncated isothermal spheres with the parameters scaled to the initial galaxy mass, $`M_i`$. Assume that the velocity dispersion scales as $`\sigma _g=\sigma _{}(M_i/M_{})^{1/4}`$, in agreement with the Tully-Fisher and Faber-Jackson relations. Here $`\sigma _{}=250`$ km s-1 and $`M_{}=2\times 10^{12}M_{\mathrm{}}`$ are the fiducial circular velocity and mass of an $`L_{}`$ galaxy. Furthermore, assume that the velocity dispersion does not change in the course of evolution but the density distribution is truncated at a radius $`R_t=R_{t,}\sigma _g/\sigma _{}`$. The latter assumption is confirmed by the self-consistent N-body simulations presented in Paper II, and $`R_{t,}`$ is determined from the simulations. The density profile of the halos is

$$\rho (r)=\frac{\sigma _g^2}{4\pi G(r^2+R_c^2)},\text{for}r<R_t$$
(A1)
where $`R_c`$ is the core radius. The enclosed mass scales almost linearly with $`R_t`$. Let galaxy 1 be more massive initially than galaxy 2; it will still be more massive after the truncation, $`M_1M_2`$. As usual, the two-body problem can be written in terms of the reduced particle with the mass $`\mu M_1M_2/(M_1+M_2)`$. The dynamics is characterized by the radius of encounter, $`\text{R}_{\mathrm{enc}}=\text{r}_1\text{r}_2`$, and the relative velocity, $`\text{V}_{\mathrm{enc}}=\text{v}_1\text{v}_2`$.

### A.1 Orbital Energy

Dynamical friction works when the halos overlap. For that the distance of closest approach during the first encounter, $`R_{\mathrm{enc}}`$, must be smaller than the tidal radius of the larger halo, $`R_{t,1}`$. Then the potential energy of the binary system is given by

$$W=\mathrm{\Phi }_1(\text{r}+\text{R}_{\mathrm{enc}})\rho _2(\text{r})𝑑\text{r}=\frac{\sigma _1^2M_2}{1\chi _2},$$
(A2)
where $`\chi _2\gamma _2\mathrm{arctan}\gamma _2^1`$, $`\gamma _2R_{c,2}/R_{t,2}`$ and

$$=\frac{1}{2}_0^1\frac{x^2dx}{x^2+\gamma _2^2}_0^\pi \mathrm{sin}\theta d\theta \left[1\frac{R_{c,1}}{|\text{r}+\text{R}_{\mathrm{enc}}|}\mathrm{arctan}\frac{|\text{r}+\text{R}_{\mathrm{enc}}|}{R_{c,1}}+\frac{1}{2}\mathrm{ln}\frac{R_{t,1}^2+R_{c,1}^2}{(\text{r}+\text{R}_{\mathrm{enc}})^2+R_{c,1}^2}\right],$$
(A3)
where $`x=r/R_{t,2}`$. The integral $``$ varies slowly with $`R_{\mathrm{enc}}`$ and is always of order unity. In the limit $`R_{c,1},R_{t,2}R_{t,1}`$ and also in the case $`R_{t,1}R_{t,2}`$, the integral evaluates to $`(1\chi _1)(1\chi _2)`$. Thus the orbital energy of the reduced particle is

$$E_{\mathrm{orb}}\frac{1}{2}\mu V_{\mathrm{enc}}^2M_2\sigma _1^2(1\chi _1).$$
(A4)
Note that it does not depend on the radius of encounter as long as the two halos overlap. Also, being enclosed within the larger halo, the smaller halo is further truncated according to $`R_{t,2}R_{\mathrm{enc}}\sigma _2/\sigma _1`$.

The dissipation of orbital energy into the internal stellar motion can produce a bound orbit even if initially $`E_{\mathrm{orb}}>0`$. This is a non-linear interaction and cannot be calculated analytically. Instead, an approximate functional dependence can be derived and then fitted to the results of N-body experiments. The velocity drag during the encounter follows from the dynamical friction analysis (Binney & Tremaine 1987; also S. Tremaine, 1998 private communication):

$$\frac{dV_{\mathrm{enc}}}{dt}\frac{4\pi G^2M_2\rho _1}{V_{\mathrm{enc}}^2},$$
(A5)
where the slow-varying Coulomb logarithm is neglected. The encounter lasts of the order $`\mathrm{\Delta }tR_{\mathrm{enc}}/V_{\mathrm{enc}}`$, so that the energy change is

$$\mathrm{\Delta }E\mu V_{\mathrm{rel}}\frac{dV_{\mathrm{enc}}}{dt}\mathrm{\Delta }t\mu \frac{M_2}{M_1}\frac{\sigma _1^4}{V_{\mathrm{enc}}^2+\sigma _1^2}.$$
(A6)
The last factor assures that equation (A6) remains valid even in the limit where the relative velocity is lower than the velocity dispersion of the larger galaxy, $`V_{\mathrm{enc}}\sigma _1`$.

The normalization is provided by the simulations of Makino & Hut (1997) who considered a hyperbolic encounter of two equal-mass galaxies. I derive $`\mathrm{\Delta }E`$ as the initial orbital energy leading to a bound pair of galaxies after the first passage. For the encounters with the impact parameter close to the galactic tidal radius, the correction factor for equation (A6) is about 2.4. Thus the galaxies become bound if they satisfy the capture criterion:

$$E_{\mathrm{orb}}<\mathrm{\Delta }E=2.4\mu \frac{M_2}{M_1}\frac{\sigma _1^4}{V_{\mathrm{enc}}^2+\sigma _1^2}.$$
(A7)
### A.2 Dynamical Friction Time

Even if two galaxies satisfy the capture criterion, they may not have enough time to merge until the present. To estimate the merger timescale via dynamical friction, assume that the velocity of the reduced particle is $`v(\sigma _1^2+\sigma _2^2)^{1/2}`$. For circular orbits of radius $`r`$, the rate of loss of the specific angular momentum is given by the Chandrasekhar formula (eq. \[A5\]):

$$\frac{dJ}{dt}r\frac{dv}{dt}\frac{M_2}{M_1}\frac{\sigma _1^4}{v^2}.$$
(A8)
For eccentric orbits, $`r`$ should be taken as the distance of closest approach. Recent work extends Chandrasekhar’s calculation and includes a time-dependent response to the perturbation (Séguin & Dupraz, 1994; Colpi & Pallavicini, 1998; Domínguez-Tenreiro & Gómez-Flechoso, 1998). However, the complicated dynamics of the interaction cannot be easily parametrized. Instead, I employ the simple functional form given by equation (A8) and use numerical simulations again to obtain correct normalization. Taking now $`dJ/dtvdr/dt`$, yields the inspiraling time

$$t_{DF}\frac{r}{\sigma _1}\left(\frac{v}{\sigma _2}\right)^3.$$
(A9)
The normalization follows from the simulations of Dubinski, Mihos & Hernquist (1999) who studied the formation of tidal tails by the interacting equal-mass galaxies. From the merging time of their models (their Table 3), I deduce $`t_{DF}r`$ and $`t_{DF}\sigma _1^{0.5}`$, roughly in agreement with equation (A9). The resulting timescale is

$$t_{DF}=0.6\left(\frac{r}{10\text{kpc}}\right)\left(\frac{\sigma _1}{250\text{km s}^1}\right)^1\frac{1}{\eta ^3}\left(\frac{1+\eta ^2}{2}\right)^{3/2}\text{Gyr},$$
(A10)
where $`\eta \sigma _2/\sigma _1`$. If the time of collision $`t_{\mathrm{coll}}`$ does not satisfy the dynamical friction criterion $`t_{\mathrm{coll}}+t_{DF}<t_H`$, the merger is disallowed.

### A.3 Tidal Disruption

According to numerical simulations of isolated mergers (Mihos, 1999), after a first encounter the galaxies separate on a highly elongated orbit before coming back for a final unification. In the cluster, the external tidal force at the distance of maximum separation can prevail the gravitational attraction of the two galaxies and disrupt the pair.

Tidal disruption in clusters is important for the following reason. Halos are truncated at the radius $`R_{t,1}`$ where the tidal force balances the gravitational pull of the galaxy. If after the encounter the smaller galaxy moves away more than $`R_{t,1}`$ from the center of the larger galaxy, the external tidal field would be stronger than the mutual attraction.

There is still a possibility, however, that the pair will merge. The tidal field varies along the galaxy trajectory and the truncation is determined by the peak values. If the collision takes place when the external tidal field is weak, the galaxies may remain bound.

Quantitatively, the distance of maximum separation $`r_m`$ can be determined from the orbital energy, which must now be negative after including the dissipated term (eq. \[A7\]). If the halos detach completely after the encounter,

$$E_{\mathrm{orb}}=\frac{\mu J^2}{2r_m^2}\frac{GM_1M_2}{r_m}<0.$$
(A11)
The angular momentum term can usually be ignored, so that the maximum distance is

$$\frac{r_m}{R_{t,1}}(1\chi _1)(1+M_2/M_1)\frac{\mu \sigma _1^2}{|E_{\mathrm{orb}}|}.$$
(A12)
If the value of $`r_m`$ determined from equation (A12) is smaller than $`R_{t,1}`$, equation (A11) does not apply and the galaxies never separate but undergo a quick merger.

The external tidal force on the larger galaxy can be written as $`F_{\mathrm{tid}}fr_m\sigma _1^2/R_{t,1}^2`$, where $`f=1`$ corresponds to the peak value responsible for the truncation. The galaxy pair remains bound at the distance of maximum separation if it satisfies the binding criterion:

$$\frac{F_{\mathrm{grav}}}{F_{\mathrm{tid}}}=\frac{(1\chi _1)(1+M_2/M_1)}{f}\left(\frac{R_{t,1}}{r_m}\right)^3>1.$$
(A13)
## Appendix B Calculating the Tidal Field

This appendix describes the procedure for calculating the tidal field in the simulations, as well as some numerical tests. The second derivative of the potential is evaluated using tri-linear interpolation from the eight nearest mesh points enclosing the center of mass of the galaxy. In order to remove the effects of small-scale noise, I apply the Savitzky-Golay smoothing filters of the 4th order (Press et al. 1992, §14.8, routine savgol). These filters replace each data point $`\mathrm{\Phi }_i`$ with a linear combination of itself and several neighbors:

$$\mathrm{\Phi }_i^{}=\underset{n=N_L}{\overset{N_R}{}}c_n^{(0)}\mathrm{\Phi }_{i+n},$$
(B1)
where $`N_L`$ and $`N_R`$ are the number of points used “to the left” and “to the right”, respectively. The first and second derivatives of the potential are obtained through similar expressions with different coefficients $`c_n^{(1)}`$ and $`c_n^{(2)}`$. The double derivatives, $`F_{xx}`$, $`F_{yy}`$, $`F_{zz}`$, are calculated at the mesh points $`(i,j,k)`$ as

$$F_{xx}(i,j,k)=\frac{1}{\mathrm{\Delta }^2}\underset{n=N_L}{\overset{N_R}{}}c_n^{(2)}\mathrm{\Phi }(i+n,j,k),$$
(B2)
where $`\mathrm{\Delta }`$ is the mesh size. The mixed derivatives, $`F_{xy}`$, $`F_{xz}`$, and $`F_{yz}`$, are calculated using the smoothing in both directions:

$$F_{xy}(i,j,k)=\frac{1}{\mathrm{\Delta }^2}\underset{n_1=N_L}{\overset{N_R}{}}\underset{n_2=N_L}{\overset{N_R}{}}c_{n_1}^{(1)}c_{n_2}^{(1)}\mathrm{\Phi }(i+n_1,j+n_2,k).$$
(B3)
The coefficients are chosen such that the 4th derivative (the smoothing order) of the field is preserved. I take $`N_L=N_R=4`$, which corresponds to the effective smoothing of the tidal field on the scale of $`250500h^1`$ kpc.

The accuracy of the above smoothing has been tested against the analytic NFW model (Navarro, Frenk & White, 1997) with the parameters of the Cluster II simulation:

$$\rho (r)=\frac{\rho _0}{(r/r_s)(1+r/r_s)^2}.$$
(B4)
For the virial radius $`R_{\mathrm{vir}}1.5h^1`$ Mpc and the “concentration” relation $`R_{\mathrm{vir}}=10r_s`$, the scale radius is $`r_s=0.15h^1`$ Mpc. In the $`512^3`$ simulation the grid cell size is $`\mathrm{\Delta }=62.5h^1`$ kpc, so that $`r_s2.5\mathrm{\Delta }`$. The exact tidal force in the NFW model is readily calculated using equation (5) of Gnedin, Hernquist & Ostriker (1999).

Figure 18 shows that our method reproduces well the mixed derivative $`F_{xy}`$ down to 3 cell sizes from the center, and somewhat worse the double derivative $`F_{xx}`$. In the case of the mixed derivative, the filter may seem to be less accurate than the straightforward two-point finite difference. But in real simulations with numerical noise the smoothing filter becomes more effective and removes small-scale artifacts.

### B.1 Test: The Pancake Simulation

As another check of accuracy of the calculation of the tidal field, I have run a test simulation of a one-dimensional Zeldovich’s pancake. The initial conditions are set up as described in Bryan et al. (1995). The wavelength of the sine wave perturbation equals the size of the simulation box; the maximum displacement is in the center. The simulation is run until a caustic of infinite density forms. From Poisson’s equation it follows that the only non-zero component of the tidal tensor is $`F_{xx}=4\pi G\mathrm{\Delta }\rho `$, where $`\mathrm{\Delta }\rho `$ is the density perturbation.

I have run two simulations, with the high ($`512^3`$ grid) and low ($`128^3`$ grid) resolution. The density profile in both simulations is reproduced accurately. The calculated tidal force matches the exact solution in the beginning while the perturbation is small, but deviates at later times. The smoothing filter cannot resolve a very narrow caustic. As expected, the higher-resolution simulation departs from the exact solution later than the lower-resolution one, at 70% vs 20% of the total simulation time. A similar effect of the resolution loss in high density regions in real simulations is emphasized again in the next section.

### B.2 Subtracting Galactic Self-Contribution

Large halos may contribute significantly to the amplitude of the local tidal field. If the halo extends over more than a couple of cells, it will contaminate the calculation of the external tidal force. It is therefore necessary to correct for the galactic self-contribution.

For this purpose, I have implemented a potential solver for an isolated galaxy. On a rectangular grid of $`16^3`$ cells, the density of the galactic particles is assigned using the CIC method. Then, the potential is obtained by the FFT method with isolated boundary conditions. Applying the same procedure for calculating the tidal tensor as above, I evaluate and subtract the halo self-contribution.

The FFT solver for Poisson’s equation with isolated boundary conditions is implemented following Hockney & Eastwood (1981). The density grid is doubled in all three directions and padded with zeros in the added space. Green’s function is translated correctly to the new sections in real space and then Fourier transformed and stored for the potential calculation. The trick with extending the grid is necessary in order to avoid the multiple images which appear in the usual FFT with periodic boundary conditions. The resulting Green’s function $`G_{\mathrm{isol}}(k)`$ differs from the usual form, $`G(k)=k^2`$, by cutting off the high frequency tail. Afterwards the algorithm proceeds as usual: Fourier transforming the density field, multiplying by $`G_{\mathrm{isol}}(k)`$, and transforming back to real space.

The isolated FFT procedure has been tested on the same $`16^3`$ grid, using $`2^{15}`$ particles distributed in the central cell according to the Hernquist (1990) model:

$$\rho (r)=\frac{M}{2\pi }\frac{a}{r}\frac{1}{(r+a)^3},$$
(B5)
with a tiny core radius, $`a=0.1\mathrm{\Delta }`$. Essentially, this is a point mass represented by particles. The algorithm reproduces the exact potential with very high accuracy, all the way to the edge of the simulation box. (For comparison, the periodic FFT produces large errors on most of the grid.) The force calculation is limited in the center at $`2\mathrm{\Delta }`$, as in any grid-based method, but behaves perfectly well at larger radii. The force calculated with the periodic FFT shows much more noise.

Another test has been applied to make sure the procedure works in a real simulation. The $`16^3`$ density field was calculated using the particles from a large galaxy, which were all located within a single grid cell and resembled a smoothed point source. Although with increased scatter, the new method is still in a good agreement with the exact potential and force.

Figure 19 shows the self-contribution of a large galaxy in the Cluster III simulation. For the ease of comparison, I plot the trace of the tidal tensor ($`F_{xx}+F_{yy}+F_{zz}`$) which is generally proportional to the density of galactic particles, in accord with Poisson’s equation. The smoothing filter misses the extreme variations of the density, but otherwise the amplitude of the tidal force is calculated correctly.

Also, Figure 19 shows the scaling $`(1+z)^3`$, which corresponds to a constant comoving density. The halo density follows this scaling quite closely indicating that the galaxy expands in real space. This numerical artifact is due to the loss of force resolution within a grid cell. The particles cannot collapse below 1 $`\mathrm{\Delta }`$, and therefore in comoving coordinates the galaxy continues to occupy the same volume.


text:
# Crack fronts and damage in glass at the nanometer scale.

## 1 Introduction

Despite many works and recent progresses, materials fracture still exhibits many puzzling aspects. Two different observations seem particularly difficult to reconcile:

* Brittle materials - whose most common example is glass - break abruptly, without first deforming, while in ductile materials - like metallic alloys - fracture is preceded by large plastic deformations. In this latter case, the crack is observed  to progress through the coalescence of micrometric damage cavities nucleated from microstructural defects (second phase precipitates, grain boundaries…).
* However, quantitative studies reveal that both types of fracture surfaces have a very similar morphology. Fracture surfaces have been shown to be self-affine objects  for both brittle and ductile materials. Two self affine regimes coexist : at small length scales (up to a length $`\xi _c`$), the roughness exponent $`\zeta `$ is close to 0.5, whereas at large length scale (from $`\xi _c`$ to a larger length $`\xi `$), $`\zeta `$ is around 0.8. These values of $`\zeta `$ are universal, but $`\xi _c`$ as well as $`\xi `$ depend on the kind of material considered, ranging from a few nanometres for glass  to a few centimetres for concrete and rocks . It is worthy to note that, for ductile materials, $`\xi _c`$ is of the order of the typical size of the damage cavities when they coalesce .

This strongly suggests that the observed transition between the two self-affine regimes reflects a change from an intra-cavity structure to an inter-cavity one. Such a scenario was first proposed in , based on the observation of the roughness of growing cavities before coalescence in an aluminium alloy and further developed in a recent paper .

How to reconciliate these two observations? A natural supposition is that a ductile-like - up to now undetected - fracture process exists in brittle materials like glass. In this case, the size of the associated damage cavities at coalescence is expected to be equal to the crossover length $`\xi _c`$, i.e. of the order of ten nanometers. The need for an experimental observation of such effect was at the genesis of the experiment described in the following. Our new experimental set-up, based on Atomic Force Microscopy (AFM), is described in section 2. Huge velocity fluctuations in the progression of the tip were detected (section 3) and explained by the observation of nanometric cavities developing ahead of the crack tip (section 4). These cavities were proved to be damage cavities by using two independent methods (section 5). Finally, the possible origins of these cavities are conjectured and the consequences of such ductile fracture mode in glass are discussed (section 6).

## 2 Experimental set-up

The design of our experiments should meet several successive requirements to achieve the goal described above. Those are listed below:

1. The apparatus should have a very high spatial resolution (a few tenths of nanometers) and work in real space and real time (no post-mortem studies). Moreover, it should deal with electrically non conductive materials. Surface treatments - which limit the spatial resolution (by metal deposition for instance) - are also prohibited. The AFM technology appears then to be the ideal tool that meets all these requirements. Due to the nanometric resolution of AFM we measure crack speeds lower than $`10^{12}\mathrm{m}.\mathrm{s}^1`$ with a reasonable acquisition time (15 minutes between each acquired image).
2. The main drawback of AFM resides in its recording time, that cannot be set lower than a few minutes without significant loss in the spatial resolution. It forces us to work in the so-called ’stress corrosion’ or ’sub-critical growth’ regime . The crack behaviour is then controlled by the chemically activated processes with the surrounding environment. A rigorous control of atmosphere and temperature is then necessary.
3. The crack progression should be stable and well monitored. This leads us to design a mechanical DCDC (Double Cleavage Drilled Compression) system .

The experimental set-up is illustrated in Fig. 1. All the experiments are performed at a constant temperature of $`22.0\pm 0.5^{}\mathrm{C}`$ in a leak-proof chamber under an atmosphere composed of pure nitrogen and water vapour at a relative humidity level of $`42\pm 1\%`$ after preliminary out-gassing.

Lithium alumino-silicate glass-ceramics were studied. The main advantage of this kind of materials is that its structure can be controlled through an applied thermal treatment . Two types of structures are investigated below. The first one (sample A) corresponds to a pure glassy state as identified by X-ray diffraction (Fig. 2a). A single thermal treatment at $`660^{}\mathrm{C}`$ is performed in order to remove residual stresses. The second one (sample B) is classified as a slightly unglassy state in which small crystals of $`\beta `$-quartz phase have formed (Fig. 2b) in the glassy phase. The structure of sample B is obtained after a two-step thermal treatment at plateau temperatures of $`T_1`$=$`750^{}\mathrm{C}`$ and $`T_2`$=$`900^{}\mathrm{C}`$. Both the size and the volume fraction of crystalline grains are evaluated by imaging the samples surface by AFM after a stay in fluorhydric acid (concentration ranging from $`0.4\%`$ to $`2\%`$) for $`30\mathrm{s}`$, which dissolves the amorphous phase faster than the crystalline one .

Fractures are then performed on the DCDC set-up : Parallelepipedic ($`4\times 4\times 40\mathrm{mm}^3`$) samples (Fig. 3) are designed with a cylindrical hole (radius $`a=0.5\mathrm{mm}`$) drilled in the centre and perpendicularly to the $`4\times 40\mathrm{mm}^2`$ surface (Fig. 3). The hole axis defines the $`z`$-direction. The $`x`$-axis (resp. $`y`$-axis) is parallel to the $`40\mathrm{mm}`$ (resp. $`4\mathrm{mm}`$) side of the $`4\times 40\mathrm{mm}^2`$ surfaces.

In all cases, the $`4\times 40\mathrm{mm}^2`$ surfaces are optically polished (the measured RMS roughness is around $`0.25\mathrm{nm}`$ for a $`10\times 10\mu \mathrm{m}^2`$ scan size). A compressive load is applied perpendicularly to the $`4\times 4\mathrm{mm}^2`$ surfaces. The external stress $`\sigma `$ is gradually increased by the slow constant displacement ($`0.02\mathrm{mm}/\mathrm{min}`$) of the jaws of the compressive machine (Fig. 1). Once the two cracks are initiated - symmetrically to the hole axis -, the jaws displacement is stopped. The crack then propagates along the $`x`$-axis in the symmetry plane of the sample parallel to the ($`x`$,$`z`$) plane. In this geometry, the stress intensity factor $`K_I`$ is given by : $`K_I=\sigma \sqrt{a}/(0.375c/a+2)`$, where $`c`$ is the crack length (Fig. 3).

The crack motion within the ($`x`$,$`y`$) surface is monitored by our experimental system combining optical microscopy and AFM. Optical image processing gives the position of the crack tip (with a spatial resolution of 5 micrometers) and consequently the velocity for $`v`$ ranging from $`10^6`$ to $`10^9\mathrm{m}.\mathrm{s}^1`$. Let us note that the measured ”instantaneous” velocity is actually averaged over a temporal window the width of which is given by time taken by the crack tip to advance on a distance equal to the spatial resolution (typically 5000s for a mean speed of $`10^9\mathrm{m}.\mathrm{s}^1`$). By AFM measurements  \- performed in a high amplitude resonant mode (”tapping” mode) -, one probes the crack tip neighbourhood at magnifications ranging from $`75\times 75\mathrm{nm}^2`$ to $`5\times 5\mu \mathrm{m}^2`$. The crack tip can be tracked at velocities ranging from $`10^9`$ to $`10^{12}\mathrm{m}.\mathrm{s}^1`$. In this latter case, the minimum width of the temporal window is reduced by two orders of magnitude thanks to the increased spatial resolution.

## 3 Crack tip velocity

At the very first moments, the crack propagates very quickly. In this regime, the crack velocity $`v`$ is independent of the chemical composition of the surrounding environment . As the crack length $`c`$ increases, $`K_I`$ decreases, and $`v`$ decreases quickly. Under vacuum, the crack stops for $`K_I`$ smaller than a critical value $`K_{Ic}`$ referred to as the toughness of the material. But in a humid atmosphere, the corrosive action of water on glass allows for slow crack propagation at much lower values of the stress . When $`K_I`$ becomes smaller than the fracture toughness $`K_{Ic}`$ (i.e. in the stress corrosion regime), the crack motion is slow enough to be monitored by our experimental system.

Figure 4 shows the variation of the velocity $`v`$ as a function of the stress intensity factor $`K_I`$ for both the amorphous and the partially devitrified specimens. The observed exponential behaviour is compatible with a stress enhanced chemical activated process..

The velocities measured for the partially devitrified sample (Fig. 4b) are shifted toward positive $`K_I`$ when compared to the data relative to the amorphous sample. To understand this shift, we probe  the surface crack path at sub-micrometric scales (Fig. 5). In the partially devitrified specimen, the crystalline germs deflect the crack (Fig. 5b). As the crack tip keeps going through the amorphous phase, the stress enhanced chemical activated process is identical in both samples (A and B) and consequently the slope of the semilogarithmic curves $`v(K_I)`$ are the same in both the amorphous and the partially devitrified specimens. However, the deflections of the crack by the crystals induce local mode II and mode III components in the local stress intensity factor, which toughens the material . Consequently, $`K_{Ic}`$ is larger in the devitrified sample, which shifts the $`v(K_I)`$ curve. A similar effect of the size of heterogeneities has been observed for the fatigue of metallic alloys .

The lowest velocities - only reachable thanks to very high spatial resolution of AFM- shown in Fig. 4 exhibit important fluctuations - of the order of the average velocity - for both specimens . Consequently, they cannot be related to the crystalline heterogeneities but are more likely inherent to the amorphous phase. The velocity fluctuations shown in figure 4 are of the order of $`5.10^{11}\mathrm{m}/\mathrm{s}`$, i.e. of the same order of the magnitude as the lowest measured averaged velocities : this fluctuations are the same in both specimens.

## 4 Evidence of nanometric cavities ahead of the crack tip

To understand the origin of the velocity fluctuation, we probed the neighbourhood of the crack tip at the nanometer scale in the amorphous specimen (Fig. 6). This clearly reveals the presence of cavities of typically $`20\mathrm{nm}`$ in length and $`5\mathrm{nm}`$ in width ahead of the crack tip (Fig. 6a) .

These cavities grow with time (Fig. 6b) until they coalesce (Fig. 6c). At these nanometric scales, the crack front does not propagate regularly, but intermittently through the merging of the nano-scale cavities, which explains the large fluctuations observed for the lowest velocitites (Fig. 4).

## 5 Nanometric cavities and nano-ductility.

To ensure that the spots observed ahead of the crack tip are actually damage cavities which grow further and coalesce with the main crack leading to failure, we use the Fracture Surface Topography Analysis (FRASTA) technique first introduced to study damage in metallic alloys . These authors showed that, in a ductile scenario, each cavity initiation is accompanied by local irreversible plastic deformations printed in relief on the developing fracture surfaces (the crack lines when the method is applied in two dimensions as in the present case) that should remain visible after the cavities have coalesced and the crack has crossed. The FRASTA method is designed to analyse the mismatch between the lower fracture surface (the lower crack line when the method is applied in two dimensions as in the present case) and the upper one (the upper crack line), mismatch due to the presence of cavities. For that purpose, the crack lines are first determined by binarising the image of the sample after fracture (Fig. 7a), and the unbroken material is reconstituted virtually by placing numerically the lower line over the upper one (Fig. 7b). By translating the lower crack line gradually in the direction of decreasing $`y`$ (Fig. 7b), - this displacement is what actually occurs during fracture -, one can see the cavities appearing and growing in the chronological order. The structure obtained for a given displacement, i.e. at a given time, is superimposed on images recorded prior to failure and shown to correspond actually to cavities observed at this given time (Fig. 7c). This clearly indicates that the spots shape is determined well ahead of the crack tip, before the first ones can be actually observed, which provides a rather strong argument to relate these nano-scale spots to damage cavities.

This nano-scale ductility is also confirmed by the study of the displacement field around the crack tip (see also references  for related discussion): for a slit-like plane crack in an ideal Hookean continuum solid, all the components of the stress tensor are predicted to decrease as $`r^{1/2}`$ with the distance $`r`$ from the crack tip . If the macroscopic stress-strain relation work on the nanometer scale the depth $`\delta `$ (see the caption of Fig. 8 for a precise definition) should thus decrease as $`r^{1/2}`$ in the vicinity of the crack tip during crack propagation. To be more precise, $`\delta =Ar^{1/2}\mathrm{cos}(\theta /2)`$ for mode I fractures  where $`A`$ is a pre-factor depending only on the applied load and on the specimen geometry, and $`\theta `$ is the angle between the direction of crack propagation and the $`\stackrel{}{r}`$ direction. Measurements of depth profiles have been performed on $`1\times 1\mu \mathrm{m}^2`$ AFM topographical frames (Fig. 8a) along the direction of crack propagation (Fig. 8b) and perpendicularly to it (Fig. 8c). For both profiles, $`\delta `$ departs from the linear elastic $`r^{1/2}`$ scaling for $`r`$ smaller than a threshold $`r_c`$ highly dependent on $`\theta `$: for $`\theta =0^{}`$, $`r_c=100\mathrm{nm}`$ while for $`\theta =90^{}`$, $`r_c=20\mathrm{nm}`$. These short range departures from the linear elastic behaviour may be related to the presence of cavities although other phenomena could be responsible for this discrepancy (as for instance viscous effect, chemical effect, …). However, the fact that the order of magnitude of the ratio $`r_c(\theta =0^{})/r_c(\theta =90^{})`$ \- much higher than the ratio of the cosine terms in the linear elastic expression of $`\delta `$ \- is close to the aspect ratio of the observed damage cavities, strongly suggests a correlation between damage and non linear elasticity.

## 6 Discussion

Same fracture experiments performed in amorphous silica specimens reveal similar damage cavities. This suggests that their existence does not depend on the precise chemical composition of the studied glass. The origin of the nucleation of cavities should be found more likely in the amorphous structure, which contains inherent atomic density fluctuations at the nanometer scale. Such a scenario was indeed predicted by Molecular Dynamics Simulation  that evidenced atomic density fluctuations in the structure of simulated amorphous silica : The Si and O atoms are shown to form silica tetrahedra connected together to build rings of different sizes ranging from 3 to 9 tetrahedra. At larger length scales, ranging from $`1.5\mathrm{nm}`$ to $`6\mathrm{nm}`$, the density of these rings is found to fluctuate with high density areas surrounded by low density areas. Moreover, the Molecular Dynamics simulations of van Brutzel  show that, at this length scale, crack propagates by growth and coalescence of small cavities which appear in areas with low density of rings, ahead of the crack tip. They behave as stress concentrators and grow under the stress imposed by the presence of the main crack to give birth to the cavities actually observed in the AFM frames.

The comparison of these theoretical results and the experimental evidences reported here strongly suggests that the crack advance in glass occurs by a ductile-like process - at a temperature far below $`T_g`$ \- bound to nucleation, growth and coalescence of damage cavities. The nucleation of these damage cavities is initiated very likely within the low density nanometric areas. The origin of the universality in the existence of two roughness regimes observed on post-mortem glass fracture surfaces  can now be explained: The fracture advances through the nucleation, growth and coalescence of damage cavities in both brittle and ductile material. The only difference resides in the length scales of the damage cavities : those related to ’brittle’ materials are approximately three orders of magnitude smaller than those relative to metallic alloys.

These observations strongly support the scenario proposed by Bouchaud et al. : The coalescence of the cavities is likely responsible for the observed roughness exponent of the $`\zeta =0.8`$ regime observed at larger length scales (from $`\xi _c`$ to $`100\mathrm{nm}`$). The order of magnitude of $`\xi _c`$ is indeed in good agreement with our present observations, where cavities at coalescence are a few tens of nanometers wide. On the other side, the $`\zeta =0.5`$ regime observed for smaller length scales is likely related to an irregular start-stop motion characteristic of the non equilibrium dynamic critical-point at which the crack starts to advance. This may generate diffusely damped crack front corrugations naturally leading to a steady state roughening of fracture surfaces. This effect is predicted to occur even though the macroscopic velocity of the crack front is small provided that the instantaneous velocity during an ‘avalanche’ in a strongly heterogeneous medium may reach much higher values. This may have been at the origin of the important fluctuations in crack velocities reported in the present paper.

Here, let us note that AFM observations are performed on the sample surface, where the mechanical state is different from that of the bulk. Hence, the observed sizes and growth rates of cavities at the surface may well differ from those in the bulk. New experiments using the FRASTA method in three dimensions applied to the post-mortem study of the fracture surfaces are currently being performed, in order to have access to the three dimensional structure of bulk damage and its evolution. Through this new set of experiments, one should be able to correlate also the damage structure to the fracture surface morphology.

The topographic study of post-mortem surfaces will also be a good experimental test of the validity of Bouchaud et al. . This scenario suggests that small length scale exponent $`\zeta =0.5`$ is due to the existence of diffusively broadened crack front waves in localised depinning events. Other important predictions are that roughness amplitude should be strongly anisotropic and this effect should be very pronounced at low crack growth velocities. That model also predicts that the physics of the formation and coalescence of the cavities is responsible for the observed roughness exponent of $`\zeta =0.8`$. Such a scenario was already proved to be realistic thanks to observations of the roughness of growing cavities before coalescence in an aluminium alloy . It will be of course crucial -even if much more experimentally difficult due to much lower $`\xi _c`$ value- to check these hypotheses in the case of nano-ductile materials.

Moreover, the structure of damage, which influences macroscopic mechanical properties such as fracture toughness and lifetime, should then be linked to the glass composition and nanostructure . Complementary analyses addressing the question of the chemical bonds on the fracture surface will also be performed.

## 7 Conclusion

Quasi-static fracture in glassy materials has been studied in real time at sub-micrometric length scales in a wide velocity range. The crack velocities are measured as a function of the stress intensity factor in both amorphous materials and a partially devitrified one. Important fluctuations are reported for the lowest velocities. They are conjectured to be related to the presence - inside of the glassy phase - of nano-scale damage cavities observed in real time ahead of the crack tip. The influence of these damage cavities on mechanical properties in the non purely elastic zone in the vicinity of the crack tip have been presented. These results confirm the scenario proposed by Bouchaud et al.  to explain the origin of the two self-affine regimes observed on fracture surfaces. The implications of such a ductile fracture mode in glass on the morphology of the fracture surfaces have been discussed. However it should be noted that our investigations are performed on the sample surface while the fracture surface morphology is related to the fracture of bulk. Consequently, it would be interesting to investigate the 3D distribution of damage cavities. Work in this direction is currently in progress.

Finally, the fact that glass, at temperatures far below the glass transition temperature $`T_g`$, joins the class of damageable materials should have important consequences for its mechanical properties. In applications, the design of structures using glass might be modified to take this behaviour into account, especially for slow crack propagation processes.

The similarity between the damage modes of materials as different as glass and metallic alloys is an important clue to understand the origin of puzzling universal behaviors, hence shedding new light into the basic physical mechanisms of fracture.

## References