Orbit decay in encounters between anisotropic spherical galaxies of equal mass

We use the Hernquist (1990) density profile to model a simple, spherically symmetric galaxy. This profile is fully characterized by two parameters, the mass $M$ and scale radius $a$. The density is

The corresponding gravitational potential is

where $G$ is the gravitational constant.

For simplicity, we use units in which $G=M=a=1$. This entails no loss of generality, as the resulting models can be rescaled to any desired system of units. In these units, the half-mass radius is $r_{\mathrm{h}}=1+\sqrt{2}$. A circular orbit at this radius has period $P_{\mathrm{h}}\simeq 33.332$; this is a characteristic timescale for dynamical evolution.

At large radii, $\rho_{\mathrm{H}}(r)$ falls off as $\rho\propto r^{-4}$. When this density profile is sampled with $N$ particles, the outermost one will have radius $r\sim aN$. Such extreme outliers are numerically inconvenient, so we use a general tapering function (Barnes, 2012, eq. 16) to smoothly taper the density profile, starting at a radius $b=100a$. This effectively limits the outermost particle to a radius of roughly $500a$. To keep the total mass of the model fixed, matter at large radii is redistributed within $r<b$; this increases the density within radius $b$ by just under $1$ per cent. In what follows, $\rho(r)$ is the tapered density profile.

Gravitational interactions between particles are softened at short range using the usual Plummer formula: at a distance $R$, the potential of a particle is proportional to $(R^{2}+\epsilon^{2})^{-1/2}$, where $\epsilon$ is a small ‘softening length’. This biases the overall gravitational potential, so $\Phi_{\mathrm{H}}(r)$ is no longer accurate. Since we need accurate potentials to compute velocity distributions (§ 2.2), the tapered density profile $\rho(r)$ is numerically convolved with a Plummer kernel (Barnes, 2012, eq. 6), and the result is used to integrate the potential as a function of radius. In what follows, $\Phi(r)$ is the resulting potential.

We use the Osipkov–Merritt method to construct radially anisotropic DFs (distribution functions) for our models (Osipkov, 1979; Merritt, 1985). For a spherical isotropic system, the DF is $f(\mathbf{r},\mathbf{v})=f(E)$, where $E\equiv\frac{1}{2}|\mathbf{v}|^{2}+\Phi(|\mathbf{r}|)$ is the specific binding energy. The Osipkov–Merritt DF is $f_{\mathrm{a}}(\mathbf{r},\mathbf{v})=f_{\mathrm{a}}(Q)$, where $Q\equiv E+\frac{1}{2}J^{2}/r_{a}^{2}$; here, $J$ is the specific angular momentum and $r_{\mathrm{a}}$ is the anisotropy radius, a parameter which controls the degree of anisotropy. Resolving the velocity into a radial component $v_{\mathrm{r}}$ and a tangential component $v_{\mathrm{t}}$, we obtain

The velocity distribution at any radius is stratified on prolate spheroids. For $r\ll r_{\mathrm{a}}$, the quantity $Q\simeq E$ and the velocity distribution is approximately isotropic. Conversely, for $r\gg r_{\mathrm{a}}$, transverse velocities are weighted more than radial velocities, so the velocity distribution is preferentially radial. The anisotropy at radius $r$ is quantified by

where the final equality is specific to Osipkov–Merritt DFs. Isotropic distributions have $\beta=0$, while radially anisotropic distributions have $0<\beta\leq 1$.

To compute the DF, given a density $\rho(r)$, potential $\Phi(r)$, and anisotropy radius $r_{\mathrm{a}}$, we first set

Since the potential is a monotonic function of radius, $\tilde{\rho}$ can be considered a function of $\Phi$. The DF is then given by

In practice, these steps must be performed numerically, since the potential $\Phi(r)$ is not available analytically. We use spline interpolation to invert $\Phi(r)$ and compute $\mathrm{d}\tilde{\rho}/\mathrm{d}\Phi$. An adaptive routine is used to compute the integral in equation (6); the result is tabulated as a function of $Q$, and spline interpolation is used for the final differentiation. When given $\rho_{\mathrm{H}}(r)$ and $\Phi_{\mathrm{H}}(r)$ as inputs, our numerical procedure reproduces the corresponding Osipkov–Merritt DF (Hernquist, 1990) with relative errors $<0.1$ per cent.

We use self-consistent N-body simulations to follow the interactions and mergers of our anisotropic models. This approach represents the DF as a collection of $N$ particles; particle $i=1,\dots,N$ has mass $m_{i}=M/N$, position $\mathbf{r}_{i}$, and velocity $\mathbf{v}_{i}$. Each particle obeys Newton’s laws of motion in the gravitational field generated by all other particles.

Particles are initialized so that the probability of finding one with coordinates $(\mathbf{r}_{i},\mathbf{v}_{i})$ is proportional to the value of $f_{\mathrm{a}}(\mathbf{r}_{i},\mathbf{v}_{i})$. Each particle is initialized independently as follows. First, the particle’s radius $r_{i}$ is chosen by generating a random number $M_{i}$ uniformly distributed in the range $[0,M)$, and then solving $M_{i}=M(r_{i})$, where $M(r^{\prime})=\int_{0}^{r^{\prime}}4\pi r^{2}\rho(r)dr$ is the cumulative mass profile. We then generate a random unit vector $\hat{\mathbf{u}}_{i}$ and set the particle position $\mathbf{r}_{i}=r_{i}\,\hat{\mathbf{u}}_{i}$. Second, at radius $r_{i}$, we use rejection sampling to pick a speed $v^{\prime}_{i}$ from the speed distribution $v^{2}f_{\mathrm{a}}(\Phi(r_{i})+\frac{1}{2}v^{2})$. Another random unit vector $\hat{\mathbf{u}}^{\prime}_{i}$ is generated, and we set $\mathbf{v}^{\prime}_{i}=v^{\prime}_{i}\,\hat{\mathbf{u}}^{\prime}_{i}$. Finally, the tangential components of $\mathbf{v}^{\prime}_{i}$ are rescaled by a factor of $(1+r_{i}^{2}/r_{\mathrm{a}}^{2})^{-1/2}$ to yield the particle velocity $\mathbf{v}_{i}$.

The first and second steps of this procedure, which suffice to generate isotropic N-body realizations, were extensively tested by Barnes (2012). The final step, which transforms an isotropic (spherical) velocity distribution into the spheroidal distribution required by the Osipkov–Merritt formalism, is new. Stability tests, to be described in the next section, confirm that our initial conditions are very close to equilibrium.

In our simulations we use $N=65536=2^{16}$ particles for each spherical model. Power-of-two values for $N$ ensure that the particle mass $m_{i}$, and integer multiples thereof, are exactly represented as floating point numbers. Our specific choice for $N$ represents a compromise between sampling accuracy and computing resources. Since particles are initialized independently, integrated quantities will have Monte-Carlo uncertainties of order $N^{-1/2}$, or $\sim 0.4$ per cent for our value of $N$. ‘Quiet-start’ procedures which suppress Monte-Carlo fluctuations are available for disc simulations (e.g., Debattista & Sellwood, 2000), but not for a general $3$-D problem like the encounters we study here. We therefore ran multiple realizations of each experiment with different random seeds to initialize the particle coordinates, and used run-to-run variation to estimate the uncertainty due to sampling.

Gravitational forces are computed using a tree algorithm (Barnes & Hut, 1986). In this algorithm, the force on each particle is approximated by a sum of $O(\log N)$ terms, which represent nearby particles and a hierarchy of cubical cells at larger distances. To be included in this hierarchy, a cell at a distance $R$ from a given particle must satisfy $d<\theta R$, where $\theta$ is a fixed parameter which controls the trade-off between speed and accuracy. We set $\theta=0.75$ and expand the gravitational field of each cell to quadrupole order. With these choices, the median acceleration error is $\sim 0.03$ per cent.

As noted above, we adopt Plummer softening to limit gravitational accelerations due to short-range interactions. We use a softening length of $\epsilon=1/64$, which is small enough to resolve the inner $\rho\propto r^{-1}$ slope of the Hernquist model (Barnes, 2012). Compared to the original potential (equation 2), this softening biases the central potential upward by $\sim 3$ per cent. However, at radius $r=a$, the potential bias is only $\sim 0.1$ per cent, and at larger radii it’s even smaller. Since the models in our simulated encounters are typically separated by distances greater that $a$, softening is unlikely to have any significant impact on our results.

We use a time-centred leapfrog with a uniform time-step to compute particle trajectories. This simple scheme is adequate for our purposes since softening limits the rate of change of particle accelerations. Our time-step is $\delta t=1/64=2^{-6}$ (again, a power-of-two value guarantees that times have exact representations). With this $\delta t$, our simulations conserve energy and angular momentum to better than $0.1$ per cent.

All our calculations were run on $4$ GHz Intel processors, using single-precision arithmetic. A typical encounter simulation, following $131072$ particles for $500$ time units, took $\sim 58$ hours of CPU time.

Spherical systems with predominantly radial orbits are often unstable, spontaneously evolving into roughly prolate configurations (Antonov 1973; Merritt & Aguilar 1985; Barnes et al. 1986). While it’s certainly possible to simulate encounters of unstable systems, it’s not easy to motivate such experiments as analogs of encounters between real galaxies; to the extent that galaxies can be described as equilibrium systems, we expect them be stable or nearly so. Thus our initial simulations were designed to determine the maximum amount of anisotropy our models can bear without becoming noticeably unstable.

In addition to visually inspecting individual simulations, we quantified the effects of the radial orbit instability as follows. After generating each model realization (§ 2.3), all $N$ particles are sorted by binding energy $E$; this ordering will be largely preserved as the system evolves. Sectioning the particle array into $8$ equal segments thus defines $8$ Lagrangian volumes or ‘shells’ ordered by binding energy. The shape of shell $j$ is measured by computing the minor-to-major axis ratio $(c/a)_{j}$ from eigenvalues of the shell’s moment of inertia tensor. The most tightly-bound shells are relatively poor indicators of instability, remaining nearly spherical even when more outlying shells are clearly distorted. We therefore used the $7^{\mathrm{th}}$ shell (i.e., the $75^{\mathrm{th}}$ to $87.5^{\mathrm{th}}$ percentile), which primarily samples the mass distribution at radii between $\sim 5a$ and $\sim 14a$, to track the development of the radial-orbit instability. (The $8^{\mathrm{th}}$ shell was not useful, as it includes a few very distant particles which individually bias the moment of inertia tensor; these particles have orbital periods much longer than the duration of our simulations.)

A trial run with $r_{\mathrm{a}}=1$ indicated that this model is unstable. We therefore tested anisotropy radii $r_{\mathrm{a}}=2^{k/8}$ where $k>0$ is an integer. This scheme samples $r_{\mathrm{a}}$ finely for small $k$ and more coarsely for larger $k$. Single realizations with $k=1,\dots,8$ were run using the N-body parameters given in § 2.3. Models with $k\geq 4$ (i.e., $r_{\mathrm{a}}\geq 1.414$) appeared stable, with $(c/a)_{7}$ values fluctuating between $\sim 0.95$ and unity. In contrast, models with $k\leq 2$ (i.e., $r_{\mathrm{a}}\leq 1.189$) clearly evolved towards more non-spherical shapes, typically falling below $(c/a)_{7}<0.85$ after $500$ time units ($\sim 15$ orbit periods at the half-mass radius).

To better constrain the boundary of stability, we ran another three realizations of each model with $k=1,\dots,4$. Fig. 1 presents the results. The additional runs confirm that models with $k\leq 2$ are unstable, while those with $k\geq 4$ appear to be stable. As a group, the models with $k=3$ (i.e., $r_{\mathrm{a}}=1.297$) typically maintain axial ratios $(c/a)_{7}\simeq 0.95$. Individual realizations fluctuate about this value, but none of the four realizations displays the systematically decreasing $(c/a)_{7}$ ratios characteristic of the unstable models.

After completing our tests, we learned that Meza & Zamorano (1997) and Buyle et al. (2007) also studied the stability of anisotropic Hernquist models with Osipkov–Merrit DFs. Meza & Zamorano estimated the stability boundary at $r_{\mathrm{a}}\simeq 1.0$, while Buyle et al. placed it at $r_{\mathrm{a}}\simeq 1.1$. Our results indicate that both of these models are unstable. The reason for this discrepancy may lie in the method used to detect the instability. Both of the earlier studies quantified the shapes of their simulations using the inner $\sim 70$ per cent of the mass distribution. Measuring shapes at smaller radii may have led these studies to overlook weak instabilities, which manifest most strongly at large radii.

We adopt the $r_{\mathrm{a}}=1.297$ model as the most anisotropic which is dynamically stable enough to be useful in this work. We cannot preclude the possibility that this model is unstable on timescales longer than the duration of our experiments; nor can we preclude an instability which manifests only at extremely large radii. Rather, we argue that such weak instabilities would have little practical consequence for our experiments. In support of this, we note that the merging behavior of our ‘maximally anisotropic’ model is continuous with the behavior of less anisotropic models.

Our experiments mimic encounters of galaxies falling together on asymptotically parabolic trajectories. Initial positions and velocities are obtained from an idealized parabolic orbit, with pericentric separation $r_{\mathrm{p}}$, of two point-like bodies each of mass $M=1$. We adopt a coordinate system in which these bodies move clockwise in the $(x,y)$ plane. At pericentre, which defines time $t=0$, they are symmetrically placed on the $x$-axis, with positions $x=\pm\textstyle\frac{1}{2}r_{\mathrm{p}}$ and velocities $v_{y}=\mp\sqrt{GM/r_{\mathrm{p}}}$. We then track the bodies back to a chosen starting time $t_{0}<0$, yielding their initial separation vector $\mathbf{r}_{0}$ and relative velocity vector $\mathbf{v}_{0}$.

We next replace these point-like bodies with N-body realizations of the anisotropic DF $f_{\mathrm{a}}(\mathbf{r},\mathbf{v})$ for a given anisotropy radius $r_{\mathrm{a}}$. A different random number seed is used for each realization. These realizations are then translated by $\pm\frac{1}{2}\mathbf{r}_{0}$ in position, and $\pm\frac{1}{2}\mathbf{v}_{0}$ in velocity, and superimposed. Thus, the initial conditions realize a system with the DF

Note that in constructing our initial conditions, we have not subtracted centre-of-mass positions and velocities from the individual galaxy realizations. As noted in § 2.3, bulk quantities have uncertainties due to Monte-Carlo statistics; as a result, each galaxy realization has an initial offset in position and velocity, of order $O(N^{-1/2})$, with respect to its assigned coordinates. While these offsets are small, they do have measurable consequences, as we show in Appendix A.

The choice of starting time requires some attention. To be entirely consistent with the chosen two-body orbit, the galaxies should ideally start so far apart that they initially interact like disjoint spherical systems. This isn’t possible for models with extended outer envelopes; some overlap has to be tolerated. We quantified the degree of overlap using

where $U_{\mathrm{tot}}$ is the total potential energy of the initial configuration, and $U_{1}$ and $U_{2}$ are the internal potential energies of the two models. Separate spherical systems yield $\delta U=0$; overlap correlates with $\delta U>0$.

Our experiments start $t_{0}=-200$ time units before pericentre, unless another time is explicitly stated. This places the models $r_{0}\simeq 69$ length units apart. For this configuration, the resulting $\delta U$ was $\sim 1.3$ per cent of $GM^{2}/r_{0}$, which seemed a bit large but acceptable inasmuch as the same bias applies to all models with the same $t_{0}$ and $r_{\mathrm{p}}$. To check, we re-ran some key experiments starting at $t_{0}=-400$. This initially placed the models $\sim 111$ length units apart; the resulting $\delta U$ was $\sim 0.42$ per cent of $GM^{2}/r_{0}$.

Table 1 lists parameters for our encounter simulations. We ran at least four realizations of each encounter, and an additional four for a subset. The extra realizations were run to test both choices for $t_{0}$ at the extremes for $r_{\mathrm{p}}$ and $r_{\mathrm{a}}$, and to anchor our results by accurately characterizing the most extreme cases.

Fig. 4 compares trajectories from encounters between anisotropic models with their isotropic counterparts. Here, the colored curves represent models with maximal anisotropy $r_{\mathrm{a}}=1.297$, while the light grey curves represent isotropic models. Each panel shows results for a different choice of idealized pericentric separation $r_{\mathrm{p}}$. Across the entire set of encounters, it’s evident that anisotropy has a strong effect on orbit decay; this effect is much larger than the run-to-run variation between models in each ensemble. The anisotropic models all attain much smaller separations after first passage, and their subsequent trajectories are more nearly radial than those of their isotropic counterparts. Indeed, the anisotropic encounters with $r_{\mathrm{p}}=1$ come to a near-standstill at apocentre, and fall back together on nearly linear trajectories. This is an example of ‘radialization’ – the tendency for orbits to become more eccentric as they decay. Radialization in unequal-mass encounters is discussed by Amorisco (2017) and Vasiliev et al. (2022); the latter study showed that radial anisotropy in the host accelerates the radialization of satellite orbits. Barnes (2016) illustrates an extreme example of radialization, in which two equal-mass galaxies with isotropic halos reverse their sense of orbital motion after first passage; passages closer than those considered here would likely exhibit this reversal as well.

Careful inspection of Fig. 4 reveals a curious trend; during the initial approach, the trajectories of the anisotropic models lie inside the idealized parabolic orbits, while those of the isotropic models generally fall on the outside. This is not the result of run-to-run variations or a problem with the initial conditions; orbit decay is modifying the trajectories of these models even before first passage.

Fig. 5 shows the same simulations, with separations plotted as functions of time. This figure again shows that the anisotropic simulations have much smaller orbits after the first passage. As one might expect, the encounters between highly anisotropic systems also reach apocentre and fall back to a second pericentre much faster than their isotropic counterparts.

We focus on first pericentre passages in Fig. 6, which plots relative velocities $v_{\mathrm{p1}}$ vs. separations $r_{\mathrm{p1}}$ for all three values of $r_{\mathrm{p}}$. This figure shows results for isotropic models (black; $r_{\mathrm{a}}=\infty$), moderately anisotropic models (blue; $r_{\mathrm{a}}=2$), and maximally anisotropic models (red; $r_{\mathrm{a}}=1.297$). Between $N_{\mathrm{real}}=4$ and $8$ realizations of each encounter are plotted as points to show run-to-run variations. For each ensemble, we plot a box centered on the sample means $\overline{r_{\mathrm{p1}}}$ and $\overline{v_{\mathrm{p1}}}$. Each box extends $\pm 2N_{\mathrm{real}}^{-1/2}s_{r}$ in the horizontal direction, and $\pm 2N_{\mathrm{real}}^{-1/2}s_{v}$ in the vertical direction, where $s_{r}$ and $s_{v}$ are sample standard deviations in $r_{\mathrm{p1}}$ and $v_{\mathrm{p1}}$, respectively (estimated assuming $1$ degree of freedom). Thus, to the extent that the central limit theorem applies for small samples, the boxes indicate $\sim 2\sigma$ uncertainties in the sample means. It is evident that pericentric separation and velocity are both influenced by anisotropy, with smaller $r_{\mathrm{a}}$ values yielding closer and faster passages.

The systematic trends with $r_{\mathrm{a}}$ seen here imply that a significant amount of orbital evolution is taking place during the approach to first passage. For equal-mass encounters, the specific orbital angular momentum at first passage is $J_{\mathrm{p1}}=\textstyle\frac{1}{2}r_{\mathrm{p1}}v_{\mathrm{p1}}$; the curves plotted in Fig. 6 are contours of constant angular momentum, passing through the mean $\overline{r_{\mathrm{p1}}}$ and $\overline{v_{\mathrm{p1}}}$ of each sample. Anisotropic models clearly lose more orbital angular momentum than their isotropic counterparts. This difference is most dramatic for close passages; among the $r_{\mathrm{p}}=1$ encounters, the most anisotropic models have lost $\sim 42$ per cent of their initial orbital angular momentum, while their isotropic counterparts have lost only $\sim 19$ per cent. The latter is consistent with a previous study of encounters between disc galaxy models dominated by isotropic dark haloes, which found average losses of $\sim 16$ per cent (Barnes, 2016).

Orbit decay in galaxy encounters is driven by the transfer of energy and angular momentum. If galaxies could not deform during tidal encounters, this transfer would be impossible. We therefore studied the shape of our galaxies prior to first pericentre, when they have not yet intermingled. Tidal deformation is already visible in Fig. 3, but the response of each galaxy is hard to follow, as they partially overlap. By contouring the surface density for just one of the two galaxies, we can accurately visualize its response to the other galaxy’s tidal field; Fig. 9 shows examples. The anisotropic simulations in the top row are conspicuously elongated, roughly towards the other galaxy, from $t=-60$ onward. The isotropic simulations in the bottom row show little sign of distortion until $t=-20$.

Fig. 10 compares shape responses of anisotropic and isotropic encounters for all $r_{\mathrm{p}}$ values. As in Sec. 2.4, we section each realization by binding energy into $8$ shells; here we plot $(c/a)_{7}$, the minor-to-major axial ratio for the $7^{\mathrm{th}}$ shell. This probes the tidal response at the relatively large radii which drive orbit decay. As Fig. 9 already suggests, the anisotropic systems show a consistently earlier and stronger response to the other galaxy’s tidal field. We observe a similar pattern of behavior for more tightly-bound shells, although the response is weaker since tides have less effect at smaller radii.

Starting time has a clear effect on the orbital evolution of anisotropic models, and in Fig. 11 we show how starting time influences shape response. For isotropic models, the choice of starting time makes little difference; these systems remain nearly spherical until they draw close to each other. In contrast, starting time has a large effect on the shapes of the anisotropic models; the encounters started at $t_{0}=-400$ are already exhibiting significant distortions (e.g., $(c/a)_{7}\simeq 0.9$) by $t_{0}=-200$.

Tidal distortion is necessary for orbit decay, but in order to actually transfer angular momentum, it’s also necessary for this distortion to lag behind the current tidal forcing. This lag arises because a stellar system doesn’t respond instantaneously to tidal perturbations; rather, a velocity perturbation imposed at time $t$ produces a density response at times $t^{\prime}>t$. In the case of a pair of galaxies making their initial approach, the major axis of each galaxy’s tidal distortion does not point directly toward its companion, but rather toward where the companion was at some earlier time. This misalignment produces a tidal torque which transfers orbital angular momentum to internal rotation.

Fig. 12 shows how orientations evolve during the approach to first pericentre. We use the simulations starting at $t_{0}=-400$, which follow pre-encounter evolution more completely. Each galaxy is represented by a thin curve plotting the angle $\theta_{7}$ of its $7^{\mathrm{th}}$ shell, measured clockwise from the $+x$ axis; this curve is plotted from the time when the shell’s axial ratio $(c/a)_{7}$ falls below $0.9$ ($0.95$) for the anisotropic (isotropic) models, respectively. The heavier solid curves are means computed over the entire sample. Finally the filled circles show angles to companion galaxies, again measured from the $+x$ axis.

This figure illustrates the delayed response necessary for orbit decay. It also points up another factor which plausibly accelerates orbit decay in encounters of anisotropic galaxies: the size of the lag. At $t\ll t_{0}=-400$ the companion galaxy’s bearing would have been $\theta\simeq 0$ degrees (i.e., along the $+x$ axis), but by $t=t_{0}$ it has swung clockwise, reaching $\theta\simeq 11$ to $22$ degrees for $r_{\mathrm{p}}=1$ to $4$, respectively. Until shortly before pericenter ($t=0$), the distance between the galaxies scales approximately as $r\propto|t|^{2/3}$, and their mutual angular velocity scales as $\dot{\theta}\propto Jr^{-2}\propto|t|^{-4/3}$. Anisotropic galaxies respond earlier to the tidal forcing (growing as $r^{-3}$), and their outer regions align roughly with the bearing of the other galaxy at early times. In contrast, the isotropic galaxies don’t respond until later, by which time the companion galaxy has swung further away from the $x$ axis. As a result, anisotropic galaxies lag more than their isotropic counterparts. This increased misalignment, along with the stronger tidal distortions already described, enable the anisotropic encounters to lose orbital angular momentum more rapidly.

To examine the effects of anisotropy in a less complicated setting, we simulated the response of single galaxy disturbed by a tide-like field. This allowed us to focus on the response without trying to disentangle the detailed interaction of two galaxies. We implemented the tide-like field by impulsively altering particle velocities. Inasmuch as any perturbation may be approximated by a series of impulses, we expect that the response to a single impulse to be governed by the same physics which is relevant during the initial stages of a tidal encounter.

Our first experiments adopted a linear tidal field aligned with the $x$-axis, mimicking the effect of a brief tidal perturbation due to a large mass located far away. However we found that this simple linear prescription had an oversized effect on particles at large radii. In particular, loosely bound particles far from the centre were given such large velocities that they became unbound.

To reduce the magnitude of the velocity perturbation on outlying particles, we next tried the following:

The parameter $A$ controls the strength of the perturbation. We rather arbitrarily set $r_{\mathrm{0}}=1+\sqrt{2}$, which is the half mass radius of our original Hernquist model. In effect, this scheme follows the linear relationship at small radii but reduces the perturbation strength for $r>r_{\mathrm{0}}$.

After perturbing the velocities, we simulated the model evolution and quantified the shape using the same methodology as in § 2.4. Experimenting with a range of values, we found that $A=0.01$ produced a clear response which relaxed in a few hundred time units. The $7^{\mathrm{th}}$ shell gave the most unambiguous results, with more tightly bound shells behaving in a similar but less dramatic fashion.

Fig. 13 compares results for three different levels of anisotropy. Clearly, anisotropic models are more easily deformed than their isotropic counterparts. The $(b/a)_{7}$ ratios closely track the $(c/a)_{7}$ ratios, indicating these distortions are prolate, as expected given the symmetry of the imposed perturbation. These results may be compared with the results shown in Fig. 10. Here the models are responding to a brief initial perturbation rather than the increasing tidal fields generated during initial approach, but in both cases the anisotropic models exhibit the strongest response. We argue that this stronger response accounts for the difference in orbit decay.

Our investigation of radial anisotropy in encounters of spherical galaxy models has shown that anisotropic models merge faster than their isotropic counterparts. In comparison to isotropic models, anisotropic models are more susceptible to tidal deformation. This accelerates the process of momentum transfer and orbit decay. Measurements of deformation confirm that anisotropic models become more elongated than their isotropic counterparts. During encounters, anisotropic models arrive at first passage having already lost significant amounts of angular momentum.