# The Survival of The Core Fundamental Plane Against Galactic Mergers

## 1 Introduction

In the framework of hierarchical galaxy formation, massive elliptical galaxy centers are merger archives, in the sense that the history of each component is imprinted upon it.

The stellar Roche limit suggests that the accretion of dense galaxies by larger, more tenuous ones would leave dense galaxies nearly intact. Hence, the central density of any massive galaxy should vary widely as a function of the extent of prior merging. Instead, the global parameters of elliptical galaxies vary simply with the luminosity. Elliptical galaxies exist in a well-defined slice of global parameter space called the Fundamental Plane (Dressler et al. 1987, Djorgovski and Davis 1987). Similarly, the central parameters of galaxies are also arranged on a plane called the core Fundamental Plane (Faber 1997).

On the core Fundamental Plane, small galaxies can be $`10^6`$ times more dense than large galaxies (Faber, 1997). This enormous dynamic range in density presents a formidable challenge for hierarchical merging because in this theory, the dense center that results from the nearly intact accretion of a small galaxy by a bigger one will cause the remnant to evolve off the Fundamental Plane. Consequently, the persistence of the Fundamental Plane is an indication either that mergers are uncommon or that additional forces are at work to destroy the secondary.

There is considerable evidence that mergers do occur. Observations of brightest cluster galaxies indicate that massive ellipticals evolve by accretion of small satellites at a rate of up to $`0.2L_{}`$ per Gyr (Lauer 1988, Blakesee and Tonry 1992). Simulations of galaxy cluster evolution also demonstrate that mergers occur, although the importance of this process in the creation of cD galaxies has been debated (Malumuth and Richstone, 1983; Merritt, 1984). So, while it is clear from the persistence of the Fundamental Plane that mergers must destroy the smaller satellites, it is not clear how such dense systems are destroyed.

Recently, Weinberg (1997) has argued that internal heating of the secondary can cause it to evaporate well inside the tidal radius. This heating occurs when some of the orbits of stars in the secondary resonate with the time-dependent tidal forces imposed by the primary. To demonstrate this, Weinberg perturbs a solution of the Boltzmann equation and recomputes the gravitational field of the secondary for King models on the global Fundamental Plane. He states that a satellite with a mass ratio as high as 100:1 will dissipate entirely by the time it falls to the primary center (Weinberg, 1997). As we will show later, although we agree with Weinberg’s result for the specific problem he did, his models lack the density contrast observed for Fundamental Plane galaxy cores, so his simulations have limited application to the problem of core Fundamental Plane persistence.

Mergers of galaxies with very different densities are expensive to simulate with a straightforward N-body approach due to the disparity in two important timescales. First, we must consider the orbital decay of the secondary into the primary, which is of order of the dynamical friction timescale, and varies as $`M_1/M_2\rho _1^{1/2}`$. We must also consider the orbital periods of stars within the secondary, which vary as $`\rho _2^{1/2}`$ where $`\rho _1`$ and $`\rho _2`$ are the densities of the primary and secondary. So, even for a modest mass ratio, tracking the stars in each galaxy over the entire merger requires on the the order of $`100M_1/M_2(\rho _1/\rho _2)^{1/2}`$ timesteps. There are few published direct simulations of this type of merger for this reason. In one important case, however, Balcells and Quinn (1990) did simulate an elliptical galaxy merger with a mass ratio of 10:1 (and density ratio of approximately 1:14). Contrary to Weinberg’s result, the secondary remained intact inside its tidal radius as it sank to the center of the primary, although the orbital decay was perhaps too fast to produce resonant heating.

We have designed an approach to determine the conditions under which core Fundamental Plane secondaries survive mergers. This approach can efficiently merge systems with mass ratios as high as 1000, a range unexplored by standard simulations. In this paper, we describe this method and apply it to the problem of secondary survival. We review the core Fundamental Plane in $`\mathrm{\S }`$ 2, our galaxy models in $`\mathrm{\S }`$ 3.1, the approximation method in $`\mathrm{\S }`$ 3.3, and our code in $`\mathrm{\S }`$ 3.4. The tests of the method and results of our simulations can be found in $`\mathrm{\S }`$ 4. Briefly, it appears that core Fundamental Plane secondaries survive mergers in which the mass ratio is 10:1 or larger. Section 5 discusses the implications of our results on the persistence of the core Fundamental Plane and previews future work.

## 2 The Core Fundamental Plane

In the parameter space composed of the effective radius, $`r_e`$, the effective surface brightness, $`\mu _e`$, and the central velocity dispersion, $`\sigma _0`$, elliptical galaxies are confined to a plane — the “Fundamental Plane” or FP (Dressler et al. (1987), Djorgovski & Davis, 1987). Recently Jorgensen et al. (1996) characterized the FP as $`\mathrm{log}r_e=1.24\mathrm{log}\sigma _00.82I_e+\gamma `$, where $`\gamma `$ may depend on environment. A popular interpretation of the Fundamental Plane is that it is a manifestation of virial equilibrium together with an orderly and slow dependence of galaxy mass–to–light ratio on luminosity. Faber et al. (1997) report the existence of an analogous plane for elliptical galaxy centers based on core parameters, though this core Fundamental Plane has more scatter than the global Fundamental Plane. The existence of such a plane indicates that the processes which act upon galaxies are primarily scale-independent. More importantly, this plane provides an important constraint on the merger history of elliptical galaxies, because if a galaxy merges, the remnant must also lie upon the plane. In particular, this implies that when a small dense galaxy collides with a massive galaxy, the remnant must have the low density center that is characteristic of a massive primary.

When simulating the mergers of real galaxies, then, it is imperative that the initial models lie on both planes. Although the Fundamental plane is 2-dimensional, we adopt the convention of Faber (1997) and represent the core galaxies as a 1-dimensional family dependent on total luminosity. This is a convenient approximation which is made possible by the tight correlation between global and break luminosity and by the nearly edge-on projection of the plane onto these axes. The core data can be fit as follows:

$`\mathrm{log}r_b`$ $`=`$ $`0.46M_v7.7`$ (1)   
$`\mathrm{log}L_b`$ $`=`$ $`1.3\mathrm{log}L_e5.097`$ (2)   
$`\mathrm{log}r_{\mathrm{eff}}`$ $`=`$ $`0.75\mathrm{log}r_b+2.06,`$ (3)   
where $`r_b`$ is the break radius, $`M_v`$ is the absolute visual magnitude, $`L_b`$ is the visual luminosity at the break radius, $`L_e`$ is the visual luminosity at the effective radius, and $`r_{eff}`$ is the effective radius. The following additional relations were obtained (Faber et al. 1997) by assuming a slightly varying mass-to-light ratio along the Fundamental Plane of $`M/LL_e^{0.25}`$,

$`\sigma _0`$ $``$ $`L_e^{0.2}`$ (4)   
$`\rho _b`$ $``$ $`\sigma _0^2/r_b^2L_e^{1.9}.`$ (5)   
This last equation is the crux of the difference between our simulations and Weinberg’s. Weinberg’s expression of the Fundamental plane is $`\rho _bM^{0.5}`$, where M is the mass of the galaxy, so for reasonable choices of $`M/L`$, his Fundamental Plane has a much shallower dependence on density with radius. See figure 1 for best-fit projections of the core Fundamental Plane to the Lauer et al. (1995) core data.

## 3 Methods

### 3.1 Modeling the Core Fundamental Plane

The disruption of the secondary galaxy centers will depend on the density ratio between the primary and the secondary, and on their density profiles. We therefore construct models designed to reflect the range in profiles observed on the core Fundamental Plane. We use $`\eta `$ models for this purpose (Tremaine et al. 1994, Dehnen, 1994). These are simple, spherical systems defined by their density distributions:

$$\rho _\eta (r)=\frac{M_{tot}}{r_{b}^{}{}_{}{}^{3}}\stackrel{}{\rho }(\eta ,r/r_b),$$
(6)   
where

$$\stackrel{}{\rho }(\eta ,x)\frac{1}{4\pi }\frac{\eta }{x^{3\eta }(1+x)^{1+\eta }},0<\eta 3,$$
(7)   
where $`\eta `$ is the single dimensionless parameter that defines the slope of the inner density profile, $`\rho `$ is the stellar density, $`M_{tot}`$ is the total mass, $`r_b`$ is the break radius, and x is the dimensionless radius $`x=r/r_b`$.

We fit the the deprojected stellar density profiles provided by Gebhardt et al. (1996) to determine $`\eta `$ as a function of the absolute magnitude $`M_v`$. The data (see Figure 1) show a correlation of the structural parameter $`\eta `$ with $`M_v`$. The ridgeline of the correlation was fit by eye as follows:

$$\eta =\{\begin{array}{cc}3.00.375M_v9.1,\hfill & \text{if }M_v19.5\text{;}\hfill \\ 3.00.15M_v4.7,\hfill & \text{if }M_v19.5\hfill \end{array}.$$
(8)   
While eta models are well-suited to cover the range of luminous matter densities in elliptical galaxy centers, the model’s density profile drops off more rapidly ($`\rho r^4`$) than a real galaxy’s at large radii. Although we anticipate that the envelope of a secondary would be completely stripped in an encounter with a primary, and is therefore irrelevant, we were concerned that truncating the envelope might have an unforseen effect on our calculations. So, we mimicked an $`r^{1/4}`$ law outer profile by constructing a basic set of galaxies with eta model envelopes. We call these galaxies double $`\eta `$ models. Their density is defined as follows:

$$\rho (r)=\frac{M_c}{r_{c}^{}{}_{}{}^{3}}\stackrel{}{\rho }(\eta _c,r/r_c)+\frac{M_e}{r_{e}^{}{}_{}{}^{3}}\stackrel{}{\rho }(\eta _e,r/r_e),$$
(9)   
where the subscript c and e corresponds to eta model parameters chosen at the core and effective radius, respectively. The core radius in equation 9 is chosen to be mimic the break radius on the core Fundamental Plane. The envelope term parameters are chosen to be consistent with global Fundamental Plane parameters and with $`r_er_{c}^{}{}_{}{}^{1.33}`$ (see equation 3) (Faber, 1997). For $`r_e`$ on the order of $`100r_c`$, $`\rho (r)`$ provides a reasonable approximation to a deVaucoleur profile.

### 3.2 Populating the Galaxies

To populate our galaxies, we require the phase space distribution function. From our spherically symmetric initial density profile, we can derive the initial distribution function, $`f(\epsilon )`$, via Eddington’s formula:

$$f(\epsilon )=\frac{1}{\sqrt{8}\pi ^2}\left[_0^\epsilon \frac{d^2\rho }{d\mathrm{\Psi }^2}\frac{d\mathrm{\Psi }}{\sqrt{(\epsilon \mathrm{\Psi })}}+\frac{1}{\sqrt{\epsilon }}\left(\frac{d\rho }{d\mathrm{\Psi }}\right)_{\mathrm{\Psi }=0}\right],$$
(10)   
where $`\epsilon =E+\mathrm{\Phi }_0`$ is the relative energy, and $`\mathrm{\Psi }=\mathrm{\Phi }+\mathrm{\Phi }_0`$ is the relative potential (Binney and Tremaine, 1987). $`\mathrm{\Phi }_0`$ is chosen so that the distribution function is positive for $`\epsilon >0`$ and zero for $`\epsilon 0`$. The second term in equation 12 is zero, because $`\mathrm{\Psi }=0`$ when $`r\mathrm{}`$, and as $`r\mathrm{}`$, $`\rho r^4\mathrm{\Psi }^4`$. Determining $`f(\epsilon )`$, then, simply requires knowledge of $`d^2\rho /d\mathrm{\Psi }^2`$.

Defining $`u=1/r`$,

$$\frac{d^2\rho }{d\mathrm{\Psi }^2}=\frac{d}{du}\left(\frac{d\rho }{d\mathrm{\Psi }}\right)/\frac{d\mathrm{\Psi }}{du}.$$
(11)   
For unscaled $`\eta `$ models,

$$\frac{d}{du}\left[\frac{d\rho }{d\mathrm{\Psi }}\right]=\left[\frac{12u^2+4(3\eta )u^3}{z}\frac{4u^3+(3\eta )u^4z^{}}{z^2}\right],$$
(12)   
where $`z=(1+u)^2`$, and $`z^{}=dz/du`$.

$$\frac{d\mathrm{\Psi }}{du}=\{\begin{array}{cc}(1+u)^\eta ,\hfill & \text{if}\text{ }\eta 1.0,\hfill \\ 1/(1+u),\hfill & \text{if}\text{ }\eta =1.0\hfill \end{array}.$$
(13)   
See Tremaine et al, 1994 for an equivalent expression of the second derivative. The second derivative for the double $`\eta `$ model case has many more terms, but follows directly from the above derivation. Given the distribution function, we know the number of particles within the phase space $`d^3\mathrm{v}d^3\mathrm{x}`$. To populate the secondary, then, we choose a random scalar radius from the mass profile, and determine the speed by choosing a random scalar velocity from:

$$f(v)dv=4\pi v^2f(\mathrm{\Psi }\frac{1}{2}v^2)dv.$$
(14)   
Finally, cartesian position and velocity components are derived from uncorrelated random orientations of r and v.

These two component galaxies were simulated with 5000 particles. In these simulations, the outer envelope was stripped away on the first encounter with the center of the primary, and we took advantage of this result by conducting additional simulations of the inner portions of the galaxies with an inner eta model of 2000 particles. This markedly improved our spatial resolution. These simulations contain a reasonable particle number, but in principle, the simulations can and will be repeated with more particles as computational speed increases. Model parameters are given in Table 1. Figure 1 illustrates the placement of the models on the core Fundamental Plane.

### 3.3 The Tidal Force on the Secondary

A straightforward technique to simulate mergers is to follow the interaction of the primary and secondary with an N-body code. However, when the galaxies involved have very different masses, the interaction takes an unreasonably long time to simulate in this straightforward manner, as described in the introduction. We have developed an approximation method which dramatically decreases the computational time required to simulate high mass contrast mergers. Our approximation hinges on the assumption that the potential of a high mass galaxy is not significantly altered by the impact of a low mass object. This approximation allows us to concentrate on the behavior of the secondary only. The merger, then, can be characterized as the secondary’s response to the external force it experiences as its orbit decays in the field of the primary.

Clearly, the external force on the secondary depends upon its position in the primary galaxy. We determine the orbital decay trajectory by solving the equations of motion for a point particle in the field of a galaxy of mass $`M_1`$:

$$\stackrel{}{F}_{12}(\stackrel{}{S})=\frac{GM_1(\stackrel{}{S})m_2}{\stackrel{}{S}^2}F_{\mathrm{fric}}(\stackrel{}{S}),$$
(15)   
where $`\stackrel{}{S}`$ is the vector from the center of the primary and to the secondary center, and $`F_{\mathrm{fric}}`$ is a modified Chandrasekhar dynamical friction formula. This also presupposes an unchanging primary, and neglects the effect of mass loss in the secondary.

$$F_{\mathrm{fric}}(\stackrel{}{S})=f_{\mathrm{drag}}\frac{4\pi \mathrm{ln}\mathrm{\Lambda }G^2\rho _1m_{2}^{}{}_{}{}^{2}}{\stackrel{}{v_2}^3}\left[\mathrm{erf}(X)\frac{2X}{\sqrt{\pi }}e^{X^2}\right]\stackrel{}{v}_2,$$
(16)   
where $`\stackrel{}{v}_2`$ is the velocity of the secondary’s center of mass, $`Xv_2/(\sqrt{2\sigma })`$, $`\sigma =\sqrt{0.4GM_1/r_{1,\mathrm{eff}}}`$, $`\mathrm{\Lambda }`$ is the Coulomb logarithm which was set to $`M_1/M_2`$, and $`f_{\mathrm{drag}}`$ is a drag coefficient. This drag coefficient is necessary, because the decay times generated by numerical treecode experiments were longer and exhibited more pericenter passes than the unmodified Chandrasekhar formula. Since the number of pericenter passes is correlated with the amount of destruction in a particular secondary, we chose to match the number of pericenter passes and the overall decay time by adjusting the analytical dynamical friction force by drag coefficient $`f_{\mathrm{drag}}`$. Figure 3 shows the treecode-generated trajectory and several analytically-derived trajectories with different drag coefficients.

Once the orbital decay trajectory is computed, the external force on particles within the secondary can be written as:

$$\stackrel{}{F}_\mathrm{S}(\stackrel{}{r})=\frac{GM_1(\stackrel{}{S}+\stackrel{}{r})m_2}{(|\stackrel{}{S}+\stackrel{}{r}|+ϵ)^3},$$
(17)   
where $`\stackrel{}{r}`$ is the vector which points from the secondary center to a secondary particle, and $`ϵ`$ is a softening parameter which we chose to be close to the core radius of the secondary. We advance the system in the secondary frame, so the appropriate force is the sum of the tidal force and the secondary’s self-gravity:

$$\stackrel{}{F}_{tot}(\stackrel{}{r})=\stackrel{}{F}_{in}(\stackrel{}{r})+\stackrel{}{F}_\mathrm{S}(\stackrel{}{r})\stackrel{}{F}_\mathrm{S}(0).$$
(18)   
We tested the assumption that the primary is unchanged by the secondary for each mass ratio by dropping a secondary point mass into a 10000 particle primary using either Hernquist’s treecode (1987) or self-consistent field code (1992). The density profile of the primary did not change for the 100:1 and 10:1 cases (see Figure 2), so the approximation is a valid one for those mass ratios. This method is more computationally efficient as the mass ratio increases, and is over 50 times faster than the Hernquist treecode in the 10:1 mass ratio.

### 3.4 Approximating the Self-Gravity of the Secondary

We constructed a self-consistent particle-field code which treats the mass distribution as a sphere and uses a logarithmic grid. This force is estimated by first convolving the mass distribution with an adaptive kernel function (Silverman, 1986) to get a density estimate. An adaptive kernel provides a non-parametric density estimate on a grid. We use an Epanechnikov kernel (an inverted parabola) for the density estimate:

$$K_e=\{\begin{array}{cc}\frac{3}{4(1r^2)}\hfill & \text{if }r<1\text{ }\text{,}\hfill \\ 0\hfill & \text{otherwise}\text{,}\hfill \end{array}$$
(19)   
where $`r=|xx_i|/h`$, x is a particle position, $`x_i`$ is a grid point, and h is the window width. An initial density estimate is obtained, then the window width of each grid element is adjusted according to the initial density such that low density regions have a large window width, and high density regions have greater mass resolution. This gives us a smoothed $`M(r)`$ from which force and potential are easily computed (Merritt, 1996). Although it is true that this code, which forces spherical symmetry, is not an appropriate choice if one wishes to model the shape of the highly non-spherical process of satellite disruption, we were concerned about whether a particle was bound to the secondary, not with obtaining an accurate shape for the stripped debris. Furthermore, since the tidal force is a very strong function of r, the core of the secondary, where the calculation matters most, is spherically symmetric.

Particle trajectories are advanced in the secondary frame with a single timestep leapfrog method.

We tested this code by following the total energy of each secondary without the influence of an external force. For double $`\eta `$ secondaries, energy was conserved to within 1 $`\%`$ over 1700 core crossing times, the time in which the 100:1 mass ratio case would have decayed had an external force been present. The mean energy and root mean square energy of the individual stars were conserved to within 2 $`\%`$ over this interval. Inner $`\eta `$ models exhibited more heating, initially, because the kernel inaccurately determined the density very near the center of the system. Since inner $`\eta `$ models have a better resolution, there were a small number of particles that were close enough to the center of the secondary for this inaccuracy to matter. We restored energy conservation by setting the mass distribution to that of the inner $`\eta `$ model inside $`0.01r_b`$, and reflecting any particles entering this sphere back outward with the same energy and angular momentum. These modifications to the code affected approximately 10 particles, but markedly improved the energy efficiency. See figure 4 for $`\delta E_{\mathrm{tot}}`$ vs. crossing time.

### 3.5 Initial Conditions

Each secondary was followed as an isolated system for several dynamical times to ensure its virial equilibrium. The merger simulation begins as the secondary is launched from $`3r_{\mathrm{half}}`$ of the primary with an initial velocity that depends upon the orbit desired; for a plunging orbit, we launch the secondary with zero velocity. We define an angular momentum parameter $`\kappa L/L_{\mathrm{circ}}`$, where L is the initial orbital angular momentum, and $`L_{\mathrm{circ}}`$ is the orbital angular momentum that the secondary would have if it were on a circular orbit about the primary at the launch point. Each simulation is followed until the secondary’s apocenter is less than $`0.1r_b`$ of the secondary.

## 4 Results

Since our goal is to investigate the breakup or survival of dense secondaries orbiting in primaries, we initially focused on the encounters we thought would be the most destructive. The destruction of the secondary, however, must occur without leaving a lasting imprint on the remnant, and this depends on where in the primary a secondary is stripped. For example, in a decaying circular orbit, the tidal force from the primary would gain strength as the secondary’s orbit decays. In this case, even if substantial damage occurred, the debris would continue to orbit near the center of the primary, and the density of the remnant would be large near the center. The orbits most effective for dispersing the secondary over a large volume must be the deeply plunging ones, since these encounters would carry the secondary past the center of the primary on its first orbit with considerable orbital energy. Our basic set of experiments, therefore, were parabolic encounters with no angular momentum. For an overview of the different experiments conducted, see table 2.

### 4.1 Basic Set of Experiments

Anticipating that the most damaging effect would be due to tidal forces experienced as the secondary passes through the center of the primary, we launched a series of nearly parabolic, plunging orbits. In this basic set of experiments, we investigated three mass ratios: 100:1, 10:1, and 2.5:1. These mass ratios correspond to density ratios of approximately 1:51, 1:20, and 1:1 at a radius of 1 pc, and the density ratio increases at smaller radii. The secondary was modeled as a double eta model in order to represent the envelope. We focus first on the 100:1 simulation. Figure 5 shows the time evolution of the secondary in the force field of the primary. Here, we learn that most of the mass is stripped at pericenter, and comes off impulsively in a cloud which continues to expand as the secondary crosses the primary center again. In figure 6, we illustrate the change in the secondary center from beginning to end. The secondary remained intact inside the tidal radius, $`r_{\mathrm{tide}}^{}{}_{}{}^{3}M_1/M_2S^3`$. Figure 7 illustrates the change in the secondary for the 10:1 mass ratio. In the 2.5:1 experiment, the secondary and primary have similar central densities, and the secondary is destroyed (figure 8).

We note that the primary and secondary in the 2.5:1 mass ratio are close enough in mass for the approximation scheme to fail, because the primary is not expected to remain unaffected by a merger of similar mass. We address this fact by conducting a treecode experiment, which we discuss in $`\mathrm{\S }`$ 4.5.

### 4.2 Comparison of Drag Coefficients

Recall the drag coefficient, $`f_{\mathrm{drag}}`$ in $`\mathrm{\S }`$ 3.3 which adjusts the strength of Chandrasekhar’s dynamical friction force to be comparable to treecode generated orbital decay times. Changing the strength of the frictional force, though, also changes the velocity of the secondary, and therefore both the magnitude of the shock on the secondary as it passes through pericenter and the number of pericenter passes. We tested the importance of this effect by simply varying the drag coefficient on the 100:1 fiducial experiment. In 2 simulations with $`f_{\mathrm{drag}}=0.1`$ and $`f_{\mathrm{drag}}=0.2`$, we saw no significant difference in the secondary’s inner density profile, although it is apparent that the secondary suffers more tidal stripping as it goes through more pericenter passes (see figure 9). Consequently, we set $`f_{\mathrm{drag}}`$ to the value determined by a treecode simulation of a point secondary in a responsive primary for each mass ratio in all further experiments. This provided a more accurate timescale for the decay, but did not affect the survival of the secondary.

### 4.3 Comparison between Inner $`\eta `$ and Double $`\eta `$ Models

To achieve better sampling of the central regions of the secondary, we populated only the center of the galaxy with an inner $`\eta `$ model. While this method excludes the secondary envelope, it is a reasonable plan because our basic set of experiments determined that the envelope is stripped away quickly without markedly changing the density profile at the center (see $`\mathrm{\S }`$ 4.1 and figures 6-8). Note that the orbital decay due to dynamical friction is the same as the previous experiment. We only followed the stars in the inner part of the galaxy but used the total mass of the secondary to determine the orbital decay.

As in the fiducial set of experiments, these simulations show that the secondary remains intact in the 100:1 and 10:1 mass ratio, and that the secondary breaks up in the 2.5:1 mass ratio (figures 10-12).

### 4.4 Non-radial Secondary Orbits

To explore the effect that different orbits have upon the disruption of the secondary, we launched secondaries in the 10:1 mass ratios on $`\kappa L/L_{\mathrm{circ}}=0.2`$ and $`\kappa =0.6`$ orbits, and the 100:1 mass ratio on a $`\kappa =0.2`$ orbit. The ratio of second apocenter to first pericenter is approximately 66 in the 100:1, $`\kappa =0.2`$ orbit, and 34 and 16 for the 10:1, $`\kappa =0.2`$ and $`\kappa =0.6`$ orbits, respectively. The computation time for these simulations required us to use inner $`\eta `$ models and to follow only the bound particles. Still, the $`\kappa =0.6`$ run took approximately 2.5 months. In every case, the secondary’s innermost density profile remained essentially intact; the secondary center was steeper and more dense than the primary, so the center of the remnant is dominated by the secondary (Figures 13-15). For a comparison of the effect that orbit has upon the disruption of the secondary, we overplot the final secondary profile for each orbit in the 100:1 and 10:1 inner $`\eta `$ experiments (figure 16).

### 4.5 Treecode Comparison

Our method works best at high mass contrasts, because for a given primary, the lower the secondary mass, the less the primary is affected by a merger. In addition, the gravitational force on the secondary is better approximated as a tide when it is small compared to the self-gravity of the secondary, as it is for a small, dense, secondary. On the other hand, the treecode is superior at low mass-contrasts because it can follow changes in the primary. We compare the results from these two problems at a mass ratio of 2.5:1 We launched a 2000 particle double $`\eta `$ secondary into a 5000 particle double $`\eta `$ primary along the same orbit as described in the fiducial experiments and followed the behavior with the Hernquist treecode.

In the treecode simulation the primary center is highly disrupted by the secondary. Our method cannot reproduce this effect. In spite of the marked difference between the behavior of the primary in each method, the qualitative behavior of the secondary is the same: the secondary central density decreases to that of the primary (figure 17). The fact that we obtain the same qualitative results provides some confidence in the method. In addition, the fact that our method fails to address the change in the potential of the primary in a clear warning that the method is quantitatively unreliable when the mass ratios are very close to 1.

### 4.6 Comparison to Weinberg’s Results

A key test of our method is to duplicate Weinberg’s experiment using our approximation method. Weinberg chose to represent the Fundamental Plane with King model galaxies (King, 1966).

King models are a single parameter set of solutions of the Vlasov equation, which can be labeled by the dimensionless potential $`W_0=V_0/(2\sigma ^2)`$ (where $`V_0`$ and $`\sigma `$ are the central potential and dispersion parameter), or equivalently by the ratio of “tidal” to core radii $`r_t/r_c`$. The mass distribution in a King model is given by $`M(r)=\mu (r/r_c)\rho _0r_c^3`$, where $`\rho _0`$ is the central density and $`\mu `$ is a dimensionless function. Weinberg adopted a single choice of $`W_0`$ in his experiments, so his models all had the same dimensionless mass distribution but scaled up or down in both density and length scale.

We chose to duplicate one of Weinberg’s 10:1 mass ratio simulations. In this experiment, he populated two $`W_0=9.5`$ King models with densities defined by the global Fundamental Plane: $`\rho M^{0.5}`$. He determined the orbital decay with Chandrasekhar’s analytical dynamical friction formula, using the Coulomb logarithm $`\mathrm{ln}\mathrm{\Lambda }=max[\mathrm{ln}(R_{\mathrm{circ}}/r_{\mathrm{half}}),0.1]`$, which caused the orbit to decay in 10 orbits for a choice of $`\kappa =0.1`$. We prefer our handling of the orbital dynamical friction because it is calibrated by the treecode. However, we wish to critically compare our calculation of the effect of the fluctuating tide to Weinberg’s problem, so we selected a choice of $`f_{\mathrm{drag}}`$ such that our orbit also decayed in 10 passes for $`\kappa =0.1`$. In all dimensionless characteristics, our experimental set up is the same. See table 3 for the parameters used in this experiment. <sup>1</sup><sup>1</sup>1 There is a typographical error in the Weinberg preprint that causes in an inconsistency between the King model concentration and the tidal radius; this prevented us from taking Weinberg’s precise radius in parsecs and his Coulomb logarithm. We circumvented this problem by creating two $`W_0=9.5`$ King models that differed in mass by 10 and density by a factor of 3, consistent with Weinberg’s model.

Weinberg contends that the secondary is disrupted by the time it reaches the center, with most of the disruption due to stripping near the primary center. In this particular case, we obtain the same result (figure 18); that is, the secondary is stripped by nearly a factor of 10, and $`\rho _2=\rho _1`$. This shows is that when we do exactly the same experiment as Weinberg, we get the same result. Therefore, the differences between our results and Weinberg’s are not due to differences in calculation methods. Instead, we believe the differences are due to different choices of galaxy models and the use of these models to mimic the Fundamental Plane. In our paper, we carefully selected each galaxy to fall on the global and core Fundamental Planes, assuming light follows mass within the galaxies.

## 5 Discussion

We have developed an approximation scheme that allows us to investigate the mergers of systems with high density ratios. We can use this method to simulate merging systems with central density contrasts as high as 1500 with less than a month of computer time on an UltraSparc.

Our simulations show that disruption will not occur for secondaries that are much denser than their primaries. While resonant heating may enhance tidal shocking in the disruption of secondary galaxies, this only appears to matter when the density contrast between primary and secondary is not too large. Weinberg’s King model secondaries evaporated because they do not possess the high densities that real secondaries exhibit on the core Fundamental Plane.

This result presents a quandary. We suspect, from observations of merger rates, that dense secondaries do merge with more diffuse primaries. In the process, the secondary must be destroyed, since observations of the centers of massive merger remnants rarely show dense, small secondaries. Instead, nearly every large remnant has a low-density center with a density profile that flattens out inside a break radius; this is a key aspect of the core Fundamental Plane. The inability of purely stellar mergers to destroy the secondaries suggests that an important component of the problem is missing.

One possible solution to this quandary may be massive black holes. The widespread presence of massive black holes in galaxy centers has been suggested, for example, by Kormendy & Richstone (1995), Magorrian (in press), Richstone (in press). Black holes act as excellent scatterers, and their presence in the centers of massive primaries may provide tidal forces sufficient to disrupt a dense secondary. We will present extensive simulations of the effect of black holes in a second paper. However, alternative possibilities include radial variations of M/L within individual galaxies, with an amplitude strongly dependent on galaxy luminosity, in the sense that luminous galaxies possess sharp concentrations of unseen mass. In order to be effective, the mass concentration would have to be more dramatic than the M/L variation with luminosity implied by the core Fundamental Plane.

Support for this work was provided by the Space Science Telescope Institute, through general observer grant GO-06099.05-94A, and by NASA through a theory grant G-NAG5-2758. DR thanks the Guggenheim Foundation and the Ambrose Monell Foundation for support at the IAS. We thank the members of the NUKER collaboration for helpful conversations. We also thank the referee, Martin Weinberg, for his insightful suggestions.