# On the response of a star cluster to a tidal perturbation

Luis A. Martinez-Medina1, Mark Gieles2,3, Oleg Y. Gnedin4, and Hui Li5
1Instituto de Astronomía, Universidad Nacional Autónoma de México, A.P.
70–264, 04510, México, CDMX, México.
2ICREA, Pg. Lluís Companys 23, E08010 Barcelona, Spain
3Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB),
Martí Franquès 1, E08028 Barcelona, Spain
4Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA
5Department of Physics, Kavli Institute for Astrophysics and Space Research,
MIT, Cambridge, MA 02139, USA Contact<EMAIL_ADDRESS>Hubble Fellow

(Released )

###### Abstract

We study the response of star clusters to individual tidal perturbations using
controlled $N$-body simulations. We consider perturbations by a moving point
mass and by a disc, and vary the duration of the perturbation as well as the
cluster density profile. For fast perturbations (i.e. ‘shocks’), the cluster
gains energy in agreement with theoretical predictions in the impulsive limit.
For slow disc perturbations, the energy gain is lower, and this has previously
been attributed to adiabatic damping. However, the energy gain due to slow
perturbations by a point-mass is similar to, or larger than that due to fast
shocks, which is not expected because adiabatic damping should be almost
independent of the nature of the tides. We show that the geometric distortion
of the cluster during slow perturbations is of comparable importance for the
energy gain as adiabatic damping, and that the combined effect can
qualitatively explain the results. The half-mass radius of the bound stars
after a shock increases up to $\sim$7% for low-concentration clusters, and
decreases $\sim$3% for the most concentrated ones. The fractional mass loss is
a non-linear function of the energy gain, and depends on the nature of the
tides and most strongly on the cluster density profile, making semi-analytic
model predictions for cluster lifetimes extremely sensitive to the adopted
density profile.

###### keywords:

galaxies: star clusters: general — globular clusters: general — open clusters
and associations: general

††pubyear: 2020††pagerange: On the response of a star cluster to a tidal
perturbation–On the response of a star cluster to a tidal perturbation

##  1 Introduction

Theories of globular cluster (GC) formation throughout cosmic time are
becoming sophisticated enough (e.g., Muratov & Gnedin, 2010; Li & Gnedin,
2014; Kruijssen, 2015; Li et al., 2017; Li et al., 2018; Li et al., 2019;
Pfeffer et al., 2018; Choksi et al., 2018; Choksi & Gnedin, 2019) that
accurate modeling of cluster dynamics in the early evolution has become more
important. GCs are expected to form in environments where the gas density is
high (Elmegreen & Efremov, 1997), which leads to tidal perturbations by
passing molecular gas clouds soon after GC formation.

Much of our understanding of the response of a star cluster to a tidal
perturbation is based on results of fast variation of the tidal field, i.e.
‘tidal shocks’ (Spitzer, 1958; Ostriker et al., 1972). However, when the
duration of the shock is longer than the crossing timescale of stars within
the cluster – which is short for dense clusters – the response of the cluster
is no longer described by the impulse approximation. During a slow
perturbation, the stars conserve orbital actions, reducing the energy gain.
The relevant parameter is the adiabatic parameter $x$, which is the ratio of
the duration of the perturbation $\tau$ and the cluster dynamical time at the
half-mass radius $t_{\rm dyn,h}$. Spitzer (1958) showed that for increasing
$x$ the effect of the perturbation is adiabatically damped, leading to an
exponential decrease of the energy gain with increasing $x$. In a series of
important studies (Weinberg, 1994a, b, c), Weinberg showed that Spitzer’s
result underestimates the energy gain of slow perturbations because of
resonances, leading to a power-law decrease of the energy gain with $x$.
Understanding how the total energy gain of the cluster depends on $x$ makes it
possible to generalize results from the impulse approximation to the adiabatic
regime by introducing the concept of adiabatic correction.

Analytically, resonances in the adiabatic regime have been explored only for a
star cluster passing through a one-dimensional slab that represents a galaxy
disc (Weinberg, 1994b, c). Similarly, adiabatic corrections have been
quantified by $N$-body simulations only for disc perturbations by Gnedin &
Ostriker (1999, hereafter GO99) and for extended spherical perturbers by
Gnedin et al. (1999), who found that the cluster energy gain decreases as
$x^{-3}$. The resulting fitting expressions are often applied to other types
of perturbers (Gnedin, 2003; Prieto & Gnedin, 2008; Pfeffer et al., 2018),
however the adiabatic corrections have not been quantified for other
perturbers. There is still a large unexplored region of the parameter space,
such as the nature of the perturbers and the density profile of the cluster.
In addition to disc-crossing, GCs also experience significant tidal forces
from compact galactic structures such as galactic bulge or giant molecular
clouds. These structures can be modeled as point mass (PM) since most
encounters are distant ones such that their sizes are typically much smaller
than the impact parameter to the cluster orbit. Gieles et al. (2006) showed
that for a population of giant molecular clouds, most of the energy comes from
penetrating encounters, for which the PM approximation does not hold. The
results of such encounters can be understood from the combination of our 2
assumed models. The distinction of different perturbers is meaningful in this
context. The tidal forces of a disc perturbation are compressive in the $z$
direction, while half of the tidal forces of a passing PM are extensive,
possibly leading to different forms of adiabatic corrections.

The aim of this work is to go beyond the impulse approximation and explore the
adiabatic regime of tidal perturbation through a set of controlled $N$-body
experiments. We investigate the dependence of the energy gain on the nature of
the perturber and properties of the cluster. This paper is organized as
follows. In Section 2 we describe the setup of the $N$-body models and the
parameterization of the tidal perturbations. Section 3 contains our results of
a study of the energy gain, mass loss, radius and density change. In Section 4
we present a discussion on the intrinsic differences between a PM perturbation
and a disc perturbation, and their impact on the adiabatic corrections.
Finally, we summarize our results in Section 5.

##  2 Description of the $N$-body experiments

###  2.1 $N$-body code

We use a state-of-the-art direct $N$-body code nbody6tt, a modified version of
the direct $N$-body integrator nbody6 (Aarseth, 2003) optimised for use with
Graphics Processing Units (GPUs, Nitadori & Aarseth, 2012). It solves pairwise
gravitational interactions between stars in the cluster. In order to apply
single tidal perturbations to each cluster we use ‘Mode A’ in nbody6tt, which
applies tidal forces in the tidal approximation via user-defined tidal tensors
as a function of time (Renaud, Gieles & Boily, 2011). We adopt the canonical
Hénon (1971) $N$-body units: $G=M_{0}=-4E_{0}=1$, where $G$ is the
gravitational constant, $M_{0}$ and $E_{0}$ are the total initial mass and
energy of the cluster, respectively. For clusters in virial equilibrium this
results in an initial virial radius of $r_{\rm v}\equiv-
GM_{0}^{2}/(2W_{0})=1$, where $W_{0}=2E_{0}$ is the initial gravitational
energy.

Shock $\Delta E_{\mathrm{imp}}/|E_{0}|$ | Shock duration $\tau$ | Cluster $W_{0}$
---|---|---
0.1 | 0.05, 0.1, 0.25, 0.5, 1, 2, 3, 4, 5 | 2, 4, 6, 8
0.03 | 0.05, 0.1, 0.25, 0.5, 1, 2, 3, 4, 5 | 4
0.01 | 0.05, 0.1, 0.25, 0.5, 1, 2, 3, 4, 5 | 4
1.0 | 0.05, 0.5 | 2, 4, 6, 8
0.5 | 0.05, 0.5 | 2, 4, 6, 8
0.1 | 0.05, 0.5 | 2, 4, 6, 8
0.05 | 0.05, 0.5 | 2, 4, 6, 8
Table 1: Parameter space explored in our $N$-body models. It encompasses a
total of 172 simulations, 86 for each perturber: PM and disc. Figure 1: The
fractional energy gain as a function of $x=\tau/t_{\rm dyn,h}$, computed after
the end of the disc (left) and PM tidal perturbation (right), with the
expected $\Delta E_{\mathrm{imp}}/|E_{0}|=0.1$ in the impulse approximation.
Solid orange line in the left panel shows the adiabatic correction fit from
GO99.

###  2.2 Setup of the perturbers and tidal tensors

We perform experiments with two types of tidal perturbers: (i) a PM with mass
$M_{\rm p}$, and (ii) an infinite disc with vertical density profile
$\rho(Z)=\rho_{\rm p}\exp(-Z^{2}/H^{2})$, where $H$ is the scale height and
$\rho_{\rm p}$ is the mid-plane density.

The duration of the perturbation is defined as $\tau\equiv b/V$ for the PM
case, where $b$ is the impact parameter and distance of closest approach
(which in the impulsive regime is the same as the impact parameter at
infinity) and as $\tau\equiv H/V$ for the disc shock. To achieve a certain
energy gain in the impulse approximation $\Delta E_{\mathrm{imp}}$, the
parameter $M_{\rm p}$ in $N$-body units is found from equation (9) of Spitzer
(1958):

$M_{\rm p}=\frac{b^{2}V}{2G}\sqrt{\frac{3\Delta E_{\mathrm{imp}}}{\langle
r^{2}\rangle}}.$ (1)

For the disc shock we find that $\Delta E_{\mathrm{imp}}=(2/3)\langle
r^{2}\rangle g_{\rm m}^{2}/V^{2}$ (Ostriker et al., 1972), where $g_{\rm
m}=2\pi^{3/2}G\rho_{\rm p}H$ is the maximum acceleration, such that

$\rho_{\rm p}=\frac{V}{GH}\sqrt{\frac{3\Delta
E_{\mathrm{imp}}}{8\pi^{3}\langle r^{2}\rangle}}.$ (2)

In all simulations we adopt $b=H=20$, but we note that in the tidal
approximation these values are not important, because the strength of the
tides are set by the density (see equations 3 and 4 below).

Each perturbation is applied to the cluster by computing the corresponding
tidal tensor as a function of time. As the tidal tensor depends on the mass
$M_{\rm p}$ for a PM shock, and on the density $\rho_{\rm p}$ for a disc
shock, the inputs to compute the tensor are: the desired fractional energy
gain in the impulsive approximation $\Delta E_{\mathrm{imp}}/|E_{0}|$, the
duration $\tau$, the mean-square cluster radius $\langle r^{2}\rangle$ (which
depends on the cluster model), and the type of the perturber (i.e. PM or
disc). With this input, $M_{\rm p}$ and $\rho_{\rm p}$ are computed from
equations (1) and (2), respectively.

Then we compute the tidal tensor of a PM shock as:

${\rm
T}_{ij}\equiv-\frac{\partial^{2}\Phi}{\partial{X_{i}}\partial{X_{j}}}=\frac{GM_{\rm
p}}{R^{3}}\left(\frac{3{X_{i}}{X_{j}}}{R^{2}}-\delta_{ij}\right),$ (3)

where $X_{i}$ are the coordinates of the PM with respect to the cluster, and
$R^{2}=\sum_{i}X_{i}^{2}$. The time dependent values of $X_{i}$ were computed
for an orbit with $b$ and $V$ at infinity, meaning that as the orbit evolves
the minimum distance (i.e. closest approach) is smaller than $b$ and the
relative velocity at closest approach is larger. For large $x$ this effect is
important.

The tidal tensor of the disc is computed as:

${\rm T}_{ij}=-4\pi G\rho_{\rm
p}\exp\left(-\frac{Z^{2}}{H^{2}}\right)\,\delta_{i3}\delta_{3j}.$ (4)

For the disc crossing we did not evolve the orbit but instead assumed a
constant $V$. This was done because the infinity disc has a constant
acceleration for $|Z|/H>>1$, which implies a strong, nonphysical correlation
between the starting $Z$ and the ratio of the initial velocity and the maximum
velocity.

In terms of the tidal tensor, the expected energy change in the impulse
approximation can be written as

$\Delta E_{\rm imp}=\frac{1}{6}\,I_{\mathrm{tid}}\,\langle r^{2}\rangle,$ (5)

where the tidal heating parameter (Gnedin, 2003)

$I_{\mathrm{tid}}\equiv\sum_{i,j}\left(\int T_{ij}\,{\rm d}t\right)^{2}\ ,$
(6)

evaluates to

$I_{\mathrm{tid}}=\frac{4g_{\rm m}^{2}}{V^{2}}=\frac{16\pi^{3}G^{2}\rho_{\rm
p}^{2}H^{2}}{V^{2}}\ ,$ (7)

for the disc shock, and to

$I_{\mathrm{tid}}=\frac{8G^{2}M_{\rm p}^{2}}{b^{4}V^{2}}\ ,$ (8)

for the PM shock.

We define the time of impact (or closest approach to the perturber) as $t=0$,
such that most of the perturbation occurs between the times $-\tau$ and
$\tau$. The values of $\Delta E_{\mathrm{imp}}/|E_{0}|$ and $\tau$ that we
consider are given in Table 1.

It is worth mentioning that the temporal resolution for the sampling of the
tensors scales with the duration of the encounter; it is 0.5 for the slowest
perturbations, $\tau=5$, and improves up to 0.01 for the fastest shocks,
$\tau=0.05$.

###  2.3 Cluster setup

We explore different density profiles for the cluster, given by King (1966)
models with dimensionless central concentration parameter $W_{0}=[2,4,6,8]$.
For these models, we find that the half-mass radius is $r_{\rm
h}=[0.849,0.827,0.804,0.871]$ and ${\langle
r^{2}\rangle}=[1.027,1.184,1.693,3.242]$ in $N$-body units. Each cluster model
was generated using limepy (Gieles & Zocchi, 2015), and discretized into
100,000 equal-mass particles.

For each cluster model and type of perturber we explore nine values of the
shock duration, $\tau$, from fast to slow perturbations. The regime at which
the shock occurs is given by the adiabatic parameter

$x\equiv\frac{\tau}{t_{\rm dyn,h}},$ (9)

where $t_{\rm dyn,h}$ is defined as in GO99:

$t_{\rm dyn,h}\equiv\left(\frac{\pi^{2}r_{\rm h}^{3}}{2GM}\right)^{1/2}.$ (10)

We find $t_{\rm dyn,h}\simeq[1.74,1.67,1.60,1.81]$ in $N$-body units, for
$W_{0}=[2,4,6,8]$. Note that despite different density profiles, all cluster
models we consider have similar values of $t_{\rm dyn,h}$ as well as of a
related quantity – the average density at the half-mass radius.

Fast perturbations (i.e. shocks) have $x\ll 1$, while slow perturbations have
$x\gtrsim 1$. We constructed models with different $x$ by changing the value
of $\tau$. This combination of parameters gives us a total of 172 $N$-body
models, summarized in Table 1. All simulations are run for
$2t_{\mathrm{imp}}$+40 $N$-body times, where $t_{\mathrm{imp}}={\rm
int}(5\tau+0.5)$ and the tensors are setup such that the maximum of the
perturbation occurs at $t_{\mathrm{imp}}$. When plotting results as a function
of time, we subtract $t_{\mathrm{imp}}$ from the $N$-body time, such that
$t=0$ corresponds to the middle of the perturbation. By starting the
perturbation at $t=-5\tau$, we miss the contribution from the tidal forces in
the time interval $-\infty<t/\tau<-5$. By integrating the analytic expressions
for the tidal acceleration of Spitzer (1958) over this time interval, we
estimate that it contributes approximately 1% to the total velocity change,
and can therefore be safely ignored.

In order to test this estimation we have therefore repeated our slowest PM
encounter ($x=3$) with a $W_{0}=4$ cluster, and added 40 time units to the
“pre-shock” part of the original simulation. Hence, this new simulation runs
for 40 + 10$\tau$ \+ 40 time units, i.e., now the encounter is symmetric. With
this new simulation we found that the originally truncated part of the
encounter accounts for 0.02% of the total energy change. This result confirms
that the truncated part of the encounters has a negligible effect.

##  3 Results

In this section we describe the results of our $N$-body models and focus on
the change in $E$, $M$ and $r_{\rm h}$ as a result of the perturbation.

###  3.1 Energy gain for different cluster models

First, we quantify the total energy gain of all stars (bound and unbound) in a
cluster model, $\Delta E$, due to a single PM or disc tidal perturbation.

Fig. 1 shows the fractional energy gain relative to the initial energy,
$\Delta E/|E_{0}|$, computed at the end of the simulation for the 72 models
with the expected $\Delta E_{\mathrm{imp}}/|E_{0}|=0.1$. For fast shocks
($x\lesssim 0.1$) there is a good agreement with the expected values from the
impulse approximation for all models. For the disc perturbations, the energy
gain due to slow perturbations ($x\gtrsim 1$) is 1-2 orders of magnitudes
lower than what is expected from the impulsive approximation. This reduced
energy gain has been reported before and attributed to adiabatic damping.

Here it is worth mentioning that for the PM case we are using the values of
$x$ computed at infinity; it is because as the orbit evolves, the relative
velocity at closest approach is larger that the initial $V$, causing the
actual $x$ value to change.

In the left panel of Fig. 1 we overplot the result of GO99 for disc perturbers
and clusters with $W_{0}=4$

$\frac{\Delta E}{|E_{0}|}=\frac{\Delta
E_{\mathrm{imp}}}{|E_{0}|}\left(1+x^{2}\right)^{-3/2}.$ (11)

The asymptotic behaviour $\Delta E\propto x^{-3}$ was derived by GO99 in the
limit of slow perturbations, $x\gg 1$. One power of $\tau$ comes from the
expectation of the linear perturbation theory, and two powers of $\tau$ come
from the normalization of the shock amplitude. Together this results in the
expected scaling $\Delta E\propto\tau^{-3}\propto x^{-3}$. Equation (11)
matches the results of the GO99 $N$-body simulations for a cluster with
$W_{0}=4$.

We find that equation (missing) 11 gives a fair description of the reduction
in energy change for our disc $W_{0}=4$ models, however, it does not explain
the results for the clusters with different $W_{0}$. The reduction is less for
more concentrated clusters, which could be because these clusters have more
stars with large periods which are in the impulsive regime. Alternatively, the
larger envelopes of more concentrated models could lead to more geometrical
distortion during slow perturbations. We discuss this further in Section 4.1.

However, we notice that by slightly changing this equation, we can still
describe our disc perturbation models. To capture the dependence of our
results on the density profile, we include an extra parameter $\epsilon$ that
depends on $W_{0}$:

$\frac{\Delta E}{|E_{0}|}=\frac{\Delta
E_{\mathrm{imp}}}{|E_{0}|}\left(1+(\epsilon x)^{2}\right)^{-3/2}.$ (12)

Fig. 2 shows the fits to our results using equation (missing) 12, where the
different values of $\epsilon$ illustrate the importance of including an
additional parameter to take into account the density profile of the cluster.
We find $\epsilon=[1.46,1.14,0.74,0.55]$ for $W_{0}=[2,4,6,8]$.

The results for slow PM perturbations ($x\gtrsim 1$) are very different than
those for the disc: for all cluster concentrations the energy gain is larger
than what is expected from the impulsive approximation. It is also clear that
none of the PM models can be described by equation (missing) 11.

Figure 2: Energy gain of clusters of different $W_{0}$ for a disc tidal
perturbation. Dashed lines are our fits given by equation (missing) 12.

Therefore we find that the difference with the impulsive prediction depends on
the type of perturber and on the density profile of the cluster. We discuss
possible causes of this behaviour in Section 4.

###  3.2 Energy gain for different perturbation strengths

Here we investigate a possible dependence of the cluster response on the
perturbation amplitude. We apply perturbations of different strength to the
$W_{0}=4$ model. For both perturbers, PM and disc, we take three values for
the perturbation strength, $\Delta E_{\mathrm{imp}}/|E_{0}|=[0.01,0.03,0.1]$.

To directly compare these cases, Fig. 3 shows the ratio $\Delta E/\Delta
E_{\mathrm{imp}}$ as a function of $x$. For the disc case, the cluster
response is almost insensitive to the perturbation strength and is well
described by equation (missing) 11. In contrast, for the PM case, the cluster
response depends on the amplitude of the perturbation in the regime of $x>1$.

Figure 3: Fractional energy gain relative to the impulse approximation. Each
panel shows three sets of $N$-body models that correspond to the shock
strengths 0.1, 0.03 and 0.01, for cluster with $W_{0}=4$. The orange line
shows the adiabatic correction fit from GO99.

It is worth mentioning that GO99 tested equation (missing) 11 in the range of
disc shock strength from $\sim 10^{-3}$ to $\sim 1$. Here we reproduce that
result for the $W_{0}=4$ models perturbed by a disc.

###  3.3 Mass loss

The induced energy gain on the cluster causes some stars to escape the system,
either by increasing their kinetic energy or reducing the binding energy. As a
consequence, an imprint of energy to the cluster translates into a mass lost.

In our simulations $\Delta M$ comes from the particles that become unbound
after the tidal perturbation. We define a particle as unbound if its specific
energy $\mathcal{E}=0.5v^{2}+\phi$ is positive at the end of the
simulation,111Note that we use a different symbol $\mathcal{E}$ to distinguish
the specific energy of an individual particle from the total energy of the
whole cluster $E$. where $v$ is the velocity (measured in the non-rotating
frame of the cluster’s center of mass) and $\phi$ the specific potential due
to the other stars (bound and unbound).

Figure 4: Fractional change in mass as a function of $x$, for a perturbation
strength $\Delta E_{\mathrm{imp}}/|E_{0}|=0.1$

In Fig. 4 we show the fractional mass loss as a function of $x$, for the disc
(left) and PM (right) perturbations for $\Delta E_{\mathrm{imp}}/|E_{0}|=0.1$.
Notice that the value of the mass loss depends on the cluster concentration.
The fractional mass loss is a factor of 20(400) higher for the disc(PM) shocks
for $W_{0}=8$ clusters than for $W_{0}=2$ clusters. For stronger shocks of
$\Delta E_{\mathrm{imp}}/|E_{0}|\simeq 0.5$ this difference reduces to a
factor of $\sim 2$, and for $\Delta E_{\mathrm{imp}}/|E_{0}|\simeq 1$ there is
no $W_{0}$ dependence (see Fig. 6).

To understand the $W_{0}$-dependence for weak shocks, we consider the fact
that the maximum energy of particles in models with different $W_{0}$ is
different: $\mathcal{E}_{\rm max}=-GM/r_{\rm t}$ and $r_{\rm t}$ for $W_{0}=8$
is about 3 times larger than for a $W_{0}=2$ model. One therefore expects that
a $W_{0}=2$ model requires a larger $\Delta\mathcal{E}$ to unbind at least one
particle. We need to keep in mind that the experiments were setup such that
clusters have the same $\Delta E_{\mathrm{imp}}/|E_{0}|$, and because $\Delta
E_{\mathrm{imp}}=\frac{1}{6}I_{\mathrm{tid}}\langle r^{2}\rangle$ and $\langle
r^{2}\rangle$ is approximately 3 times smaller for $W_{0}=2$ than for
$W_{0}=8$, the value of $I_{\mathrm{tid}}$ is 3 times larger for $W_{0}=2$
models. This means that at $r_{\rm t}$ the expected energy gain
($\Delta\mathcal{E}_{\rm t}\propto I_{\mathrm{tid}}\,r_{\rm t}^{2}$) is 3
times larger for $W_{0}=8$ models, thereby unbinding more stars that are near
$r_{\rm t}$.

Meanwhile, in the adiabatic regime and for disc shocks, the cluster mass loss
decreases as the perturbation becomes slower, following a trend similar to the
energy gain shown in Fig. 1. For the PM perturbation $|\Delta M|/M_{0}$
behaves in the same way as for the disc in the impulsive regime, i.e., it is
nearly constant, and its value strongly depends on $W_{0}$. However, for slow
tidal perturbations, the cluster mass loss increases, opposite to what happens
with the disc.

It is worth pointing out that although fast PM and disc shocks converge to the
same $\Delta E$ value (see Fig. 1), they do not converge to the same $\Delta
M$ value. Fig. 4 shows that, in the impulsive regime and for the same $\Delta
E$, the cluster systematically losses more mass under a disc shock than under
a PM shock; this difference is negligible for high concentration clusters
($W_{0}=8$), but becomes significant as the concentration decreases. For low
concentrations ($W_{0}=2$) the cluster mass loss under a fast disc shock can
be more than an order of magnitude higher than under a PM fast shock, for the
same energy gain.

This discrepancy is likely the result of the different nature of the tidal
forces: for the disc, the tidal forces are compressive along the $z$-axis,
whereas for the PM case, half of the tides act along the $x$-axis (extensive)
and the other half acts along the $z$-axis (compressive). This means that for
the disc $\Delta v=\Delta v_{z}$, while for the PM $\Delta v\simeq\sqrt{\Delta
v_{x}^{2}+\Delta v_{z}^{2}}$ and $\langle\Delta v_{x}^{2}\rangle=\langle\Delta
v_{z}^{2}\rangle=0.5\langle\Delta v^{2}\rangle$. Therefore the maximum energy
gain of individual particles, $\Delta\mathcal{E}$, in the PM case is half of
the maximum $\Delta\mathcal{E}$ in the disc case. We illustrate this in Fig. 5
by showing the $\Delta\mathcal{E}$ as a function of their initial energy
$\mathcal{E}_{0}$ at various moments in time for both a PM and a disc case.
The first column shows $\Delta\mathcal{E}$ in the middle of the shock, and the
evolution proceeds to the right. The dashed lines indicate the boundary
between bound and unbound. Soon after the shock ($t=1$), the loosely bound
stars in the disc case have increased more in energy than in the PM case,
resulting in more unbound stars (red dots).

Figure 5: Energy gain of particles at different times (in $N$-body units) for
an impulsive shock ($x=0.03$) due to a point-mass (top row) and a disc (bottom
row). In both cases the cluster is described by a $W_{0}=2$ and the strength
of the shock is $\Delta E_{\mathrm{imp}}/|E_{0}|=0.1$. Stars that have become
unbound are above the dashed line and are shown in red. More stars become
unbound for a disc shock because all tidal forces are exerted along a single
axis (see text for details).

Figure 6: Fractional mass loss vs. energy change on a linear (top panel) and
log-log scale (bottom panel). $N$-body results shown by symbols are
differentiated by shock type and cluster density profile. Lines show linear
regression fits, described in the text.

Fig. 6 shows the fractional mass loss as a function of the energy gain for the
64 simulations described in the second part of Table 1, i.e., four values of
the perturbation strength, two values of the duration, four values of $W_{0}$,
and two perturber types. We show this relation first on a linear scale because
this emphasizes the large values of mass loss which are of primary interest in
studies of cluster evolution. Then we show it on a log-log scale to emphasize
the weak shock regime.

There is large dispersion in mass loss at low $\Delta E/|E_{0}|$, which makes
it difficult to find a simple relation between mass loss and energy gain. At
large energy gains ($\Delta E/|E_{0}|>0.1$), the fractional mass loss is well
approximated by a linear relation

$\frac{|\Delta M|}{M_{0}}\simeq f\,\frac{\Delta E}{|E_{0}|},\hskip
22.76219ptf=0.217\,\pm 0.005,$ (13)

with rms of the fit $\sigma=0.022$. This result is consistent with the
relation found by Gieles et al. (2006) ($f=0.22$) for Plummer (1911) models.
It is worth mentioning that, although the error in the fit is relatively
small, at low values of $\Delta E/|E_{0}|$ the scatter is not random: instead,
the models with $W_{0}\leq 4$ systematically fall bellow the linear fit, and
those with $W_{0}\geq 6$ fall on or above the fit.

We fit the linear mass–energy change relation individually for different
$W_{0}$ models and find the proportionality factor $f$ to vary between 0.18
and 0.22 for the disc shocks, and between 0.20 and 0.26 for the PM shocks.
Interestingly, the dependence on $W_{0}$ is not monotonic: the $W_{0}=4$ model
has the smallest slope. Most concentrated clusters have the largest mass loss,
as already seen in Fig. 4. The adiabatic parameter $x$ plays a smaller role,
for the range of values we considered $(x<1)$. These results also show the
dependence on the perturber type. The mean values are $f=0.195\pm 0.006$ for
disc, $0.233\pm 0.008$ for PM, but they are still relatively close.

Fig. 6 also indicates that the relation between mass and energy change can be
alternatively described as non-linear. Indeed, we fit a power-law relation

$\frac{|\Delta M|}{M_{0}}=f_{\beta}\,\left(\frac{\Delta
E}{|E_{0}|}\right)^{\beta},\hskip 17.07164ptf_{\beta}\approx 0.25$ (14)

and find the slope systematically decreasing from $\beta=2.01\pm 0.13$ for
$W_{0}=2$ to $\beta=0.88\pm 0.06$ for $W_{0}=8$. The constant of
proportionality of this relation is tightly constrained to be
$f_{\beta}=0.23-0.26$. These relations are shown in the bottom panel of Fig.
6. We can see again that the more concentrated cluster models show larger
relative mass loss and the less steep relation.

Despite the present scatter and non-linearity for weak shocks, the relation
(13) for strong shocks ($\Delta E/|E_{0}|>0.1$) can be used in analytical or
sub-grid models of cluster disruption, when the energy change from tidal
shocks can be calculated. In addition to specifying the information about the
shocks ($I_{\mathrm{tid}}$), this requires assumptions about the cluster
profile determining $\langle r^{2}\rangle$. To eliminate the latter assumption
we can rephrase the mass loss of a cluster of a given density profile due to a
perturbation with strength $I_{\mathrm{tid}}$.

In Fig. 7 we plot the fractional mass loss as a function $I_{\mathrm{tid}}$.
The scatter is much larger than in Fig. 4, because high(low) concentration
clusters have moved to the left(right). This highlights the importance of the
dependence of $\Delta E$ (and therefore $\Delta M$) on $\langle r^{2}\rangle$
of the cluster, which can vary by a factor of 3 for clusters with the same
$r_{\rm v}$ (and almost the same $r_{\rm h}$).

We find that the relation between $\Delta E$ and $I_{\mathrm{tid}}$ is almost
exactly linear for a given cluster model $W_{0}$ (and therefore $\langle
r^{2}\rangle$), as expected from analytical theory (equation (missing) 5):
$\Delta E\propto I_{\mathrm{tid}}^{0.98\pm 0.01}$. However, the relation
between $\Delta M$ and $I_{\mathrm{tid}}$ is noticeably non-linear:

$\frac{|\Delta M|}{M_{0}}\propto I_{\mathrm{tid}}^{\gamma}$ (15)

with the power-law index decreasing systematically with $W_{0}$, from
$\gamma=1.99\pm 0.13$ for $W_{0}=2$ to $\gamma=0.88\pm 0.06$ for $W_{0}=8$.
These slopes are very close to those for the mass-energy relation, which means
that the non-linearity of the $\Delta M(I_{\mathrm{tid}})$ relation is
inherited from the $\Delta M(\Delta E)$ relation.

Figure 7: Fractional mass loss vs. tidal heating parameter $I_{\mathrm{tid}}$,
for the most impulsive runs ($x\approx 0.03$). Symbols and colours are the
same as in Fig. 6. Figure 8: Fractional change in the half-mass radius of the
bound stars as a function of the adiabatic parameter, for the perturbation
strength $\Delta E_{\mathrm{imp}}/|E_{0}|=0.1$. Figure 9: Fractional change in
the average density of bound stars within the half-mass radius, $\rho_{\rm
h,bound}$, as a function of adiabatic parameter, for a perturbation strength
$\Delta E_{\mathrm{imp}}/|E_{0}|=0.1$.

###  3.4 Changes in half-mass radius and density

Besides the energy change and mass loss, $r_{\rm h}$ of the cluster is also
modified by the perturbation. Soon after the perturbations, the particles with
the largest energy change will move to larger orbits thereby increasing
$r_{\rm h}$. However, the most energetic particles may be unbound and leave
the cluster, such that the half-mass radius of the remaining bound stars,
$r_{\mathrm{h,bound}}$, could in fact be smaller than the initial half-mass
radius, $r_{\rm h,0}$.

Fig. 8 shows the fractional change for the bound stars, $\Delta
r_{\mathrm{h,bound}}/r_{\rm h,0}$, as a function of $x$. First note that the
change in $r_{\mathrm{h,bound}}$ is larger for the less concentrated clusters
regardless of the perturber. However, the bound system expands if the cluster
is not too concentrated, while for our most concentrated model ($W_{0}=8$) it
contracts – $\Delta r_{\mathrm{h,bound}}$ is negative. This is because for
concentrated clusters, most of the added energy is absorbed by the distant
stars which are already close to being unbound. Removing stars with near-zero
energy results in a change of cluster mass at near-constant total energy, i.e.
$E_{\rm tot}\propto M^{2}/r_{\rm h}$ is approximately constant and therefore
$r_{\rm h}\propto M^{2}$ in the limit that the energy of the remaining bound
stars is the same as the pre-shock energy. In less concentrated clusters, the
mass that is lost also removes some negative binding energy, thereby
increasing the energy of the remaining bound stars which leads to an increase
of the radius (Gieles & Renaud, 2016).

To have the full picture, we combine our measurements of mass and radius
changes to compute the fractional change in the average density of bound stars
within the half-mass radius, $\rho_{\rm h,bound}$. Fig. 9 shows that after
either a disc or PM perturbation, the bound system increases its density for
the case of highest-concentration clusters, while for less concentrated
clusters the density decreases. Whether the density increases or decreases
defines whether the next shock is less or more disruptive, respectively.

As in the previous plots, Fig. 9 also shows clear differences between a disc
perturbation and a PM perturbation outside the impulsive regime. As the
adiabatic parameter $x$ increases, the change in density of the high-
concentration clusters grows in the PM case, while it decreases for the low-
concentration clusters. On the other hand, in the disc case the changes in
density converge to zero, as changes in both mass and radius diminish (Figs. 4
and 8).

Note that the changes in the cluster density imply a change in the value of
the adiabatic parameter $x$. Because $t_{\rm dyn,h}\propto\rho_{h}^{-1/2}$,
the changes in $x$ and $\rho_{h}$ are related by $\Delta
x/x=0.5\Delta\rho_{h}/\rho_{h}$. This means that the fractional changes in
$\rho_{h}$ that we see in Fig. 9 translate to changes in the value of $x$. We
find that the adiabatic parameter oscillates for a few $N$-body times before
returning close to its initial value. In the case of fast shocks, the
amplitude of the oscillations can be up to 25% of the initial value of $x$. In
the case of $x\gtrsim 1$ the variation is negligible. Therefore, the regime of
the encounter (quantified by $x$) remains nearly the same as initial.

##  4 Discussion and interpretation

We have shown that a star cluster responds differently to a PM tidal
perturbation than to a disc perturbation. Moreover, for each perturber and for
a given amplitude of the applied shock (i.e. $I_{\mathrm{tid}}$), the response
depends also on the cluster concentration (i.e. $W_{0}$). Below we discuss
possible reasons for this behaviour.

###  4.1 The dependence of $\Delta E$ on cluster concentration for slow
perturbations

Fig. 1 shows that for slow perturbations with the same $\Delta
E_{\mathrm{imp}}$, $\Delta E\ll\Delta E_{\mathrm{imp}}$ for discs, and $\Delta
E>\Delta E_{\mathrm{imp}}$ for PMs. In previous work, the reduced energy gain
for discs was ascribed to adiabatic damping. We first consider this effect to
see whether it could explain the correlation of $\Delta E$ with $W_{0}$. High
$W_{0}$ clusters may have a larger fraction of stars in the impulsive regime,
reducing the effect of damping. To test this idea, we used the discrete
realisations of King models with different $W_{0}$ from our initial
conditions, and then analogous to equation (missing) 10 we calculated an
estimate of the dynamical time for each particle as $t_{\rm
dyn}(r)=[\pi^{2}r^{3}/4GM(r)]^{1/2}$, where $r$ is the radial position of the
particle and $M(r)$ is the enclosed mass at that radius. The fraction of stars
with the local dynamical time exceeding the shock duration for our slowest
perturbation, $t_{\rm dyn}(r)>\tau\approx 3t_{\rm dyn,h}$ increases
monotonically from 0.7% for $W_{0}=2$ to 21% for $W_{0}=8$. It is consistent
with the trend of the parameter $\epsilon$ in equation (missing) 12 decreasing
with $W_{0}$ from 1.46 to 0.55, indicating that the transition from impulsive
to adiabatic occurs at higher $x$ for larger $W_{0}$. However, this does not
explain fully the qualitative difference between the results for the disc and
PM cases shown in Fig. 1.

We therefore also look at the effect of geometrical distortions during the
perturbations. For this estimate, we first consider the disc shock because of
its simple one-dimensional tidal action, and the two extreme $W_{0}$ clusters
in the suite: $W_{0}=2$ and $W_{0}=8$. We denote their respective quantities
with subscript $2$ and $8$. Consider a slow perturbation with a given $\Delta
E_{\mathrm{imp}}=\frac{1}{6}\,I_{\mathrm{tid}}\,\langle r^{2}\rangle$. Because
$\langle r^{2}\rangle_{8}\simeq 3\langle r^{2}\rangle_{2}$ (Section 2.3) and
$\Delta E_{\mathrm{imp}}$ is the same for both clusters, we have $I_{2}\simeq
3I_{8}$. The velocity change for a particle at $r\simeq r_{\rm t}$ is $\Delta
v\propto\sqrt{I_{\mathrm{tid}}}r_{\rm t}$, and because ${r_{\rm t}}_{,8}\simeq
3{r_{\rm t}}_{,2}$ (Section 2.3), we find that $\Delta v_{8}/\Delta
v_{2}\simeq\sqrt{3}$. The fractional compression along $z$ direction in the
first half of the perturbation is $F\propto\Delta v/r_{\rm t}$, and therefore
the ratio $F_{2}/F_{8}\propto\sqrt{I_{2}/I_{8}}=\sqrt{3}$. This means that,
somewhat counter-intuitively, less concentrated clusters (low $W_{0}$) are
more compressed during the perturbation than concentrated clusters. This could
also explain why the final $\Delta E$ of low-concentration cluster is lower,
because their $\langle z^{2}\rangle$ are reduced more relative to the initial
ones during the second half of the perturbation.

The PM results can be explained with these simple geometrical arguments: half
of the tidal forces are compressive, hence the same arguments of the disc
apply. However, the other half are extensive, and this should lead to a larger
fractional extension for low$-W_{0}$ clusters, resulting in a larger $\Delta
E$ for low concentration clusters. If we apply symmetry arguments, we expect
the effect of compression on $\Delta E$ to cancel the effect of the extension,
predicting that the energy gain for slow PM perturbations is the same as for
fast shocks.

Figure 10: Cumulative gains in kinetic energy $\Delta T/|E_{0}|$ (solid
lines), potential energy $\Delta W/|E_{0}|$ (dashed lines), and total energy
$\Delta E/|E_{0}|$ (dotted lines), computed at the middle of the encounter for
a disc (left) and PM perturbation (right), on a $W_{0}=4$ cluster. The colours
indicate different values of the adiabatic parameter.

###  4.2 On the different response of a star cluster to a disc vs. PM
perturbation

To explore in more detail the cause of the differences between a disc and a PM
tidal perturbation we split the total energy gain into the kinetic ($K$) and
potential ($W$) energy parts. For this exercise we look at the initial
configuration and at the middle of the shock. Each of the clusters is binned
in radial shells such that each shell contains the same number of particles,
in this case we choose 1% of the total of particles per shell. For each shell,
we compute the cumulative kinetic and potential energy of the particles
contained within the shell. Then we subtract the initial $K$ and $W$ from the
values at the middle of the shock. With this procedure we obtain the
cumulative $\Delta K$ and $\Delta W$ across the cluster.

Fig. 10 shows $\Delta T/|E_{0}|$ and $\Delta W/|E_{0}|$ as a function of the
percentage of particles from the cluster centre. Note that the evolution of
kinetic and potential energies depends strongly on the nature of the
perturber. The gain in kinetic energy $\Delta T$ for the PM perturbation
decreases as the duration of the perturbation increases, but has the opposite
behavior for the disc. The gain in potential energy $\Delta W$ for the PM
perturbation increases ($W$ becomes less negative) as the perturbation enters
into the adiabatic regime, while for the disc $W$ becomes more negative.

The reason for the opposing behaviours in both $\Delta T$ and $\Delta W$ lies
in the different directions of the tides for the two types of perturber. The
PM tides have equal extensive and compressive components (Spitzer, 1958),
while they are fully compressive for a disc case. As shown in Fig. 1, the
nature of the tide is not important for the energy gain in the impulsive
regime, but it becomes important when the duration of the perturbation
increases.

Fig. 10 shows that in the adiabatic regime, as $x$ gets closer to 1, the
extensive tide due to the PM has enough time to pull the stars apart,
decreasing the velocity dispersion ($\Delta T<0$) and reducing the binding
energy ($\Delta W>0$). On the other hand, the compressive tide due to the disc
pushes the stars inwards, increasing the velocity dispersion ($\Delta T>0$)
and making the cluster more bound ($\Delta W<0$).

In the formulation of tidal heating by Weinberg (1994b), a star’s energy
changes when it passes through a orbital resonance with the external
gravitational potential. The number and distribution of resonances is
particularly important for a slow perturbation where the gradually evolving
stellar orbit may go through multiple resonances. As described by Murali &
Weinberg (1997), a spherical potential has a different (discrete) spectrum of
resonances than the disc potential (continuous spectrum). This may be another
reason for the different cluster response to the two perturbers.

###  4.3 The effect of the orbit

Another effect that could be important is that in the PM perturbation the
impact parameter and relative velocity are time dependent, while for the disc
we fix it (see Section 2.2). Here we estimate how much the PM orbit deviates
from a straight line. We can do this by realising that for a straight-line
orbit the kinetic energy of the orbit is much larger than its gravitational
energy, i.e. $2GM_{\rm p}/(bV^{2})\ll 1$. With the expression for $M_{\rm p}$
from equation (1) and the definition of $\tau$ and $x$ we then find that the
orbit is approximately a straight line when

$x\lesssim 1.6\sqrt{\frac{0.1}{\Delta E_{\mathrm{imp}}/|E_{0}|}}.$ (16)

Here we used the $\langle r^{2}\rangle$ values for $W_{0}=6$. From this we see
that the 3 slowest encounters in Fig. 1 are affected by this, but the
deviation from a straight-line orbit can not explain the difference between
the behaviour of slow perturbations by a disc and a PM. The dependence on
$\Delta E_{\mathrm{imp}}$ may explain why the weaker (slow) PM perturbations
in Fig. 3 show a larger $\Delta E$, their orbits are more gravitationally
focused and have a smaller impact parameter ($p$) which leads to a larger
energy gain. The maximum velocity ($V_{\rm max}$) is higher than $V$, but this
does not compensate the smaller $p$, because the energy gain depends on
$p^{-4}V_{\rm max}^{-2}$.

###  4.4 The shape of the cluster

In Sec. 3.1 we found that not all our $N$-body models can be described by the
analytic result for the energy gain of GO99 (equation (missing) 11). The
result of GO99 is a fitting function with the theoretically expected power-law
slope at large $x$ ($x^{-3}$) and GO99 did not vary their cluster density
profile to study its effect. Because of this, in this section we contemplate
the possibility that the size and shape of the cluster are changing during the
encounter, which could explain the disagreement between our simulations and
equation (missing) 11, and between the disc and PM perturbation.

To quantify the deformation of the cluster due to the tidal perturbation we
need to measure its shape. Following Zemp et al. (2011), we calculate the
shape tensor, defined by the matrix components:

$S_{jk}=\frac{\sum_{i}m_{i}\,x_{ij}\,x_{ik}}{\sum_{i}m_{i}},$ (17)

where index $i$ runs over $i=1,...,N$ stars, and indices $j$ and $k$ denote
the three Cartesian coordinates for star $i$. The principal axes of the
cluster can be interpreted as the eigenvectors of $S_{jk}$. The lengths of the
primary semi-axes $a\geq b\geq c$ correspond to the square root of the
eigenvalues. We compute the matrix $S_{jk}$ for all stars in our $W_{0}=4$
models under the fast and slow perturbations.

Fig. 11 shows the evolution of the semi-axes of the cluster as a function of
time for the PM case and $\Delta E_{\mathrm{imp}}/|E_{0}|=0.1$. The grey
shaded area indicates the time interval $\pm\tau$ where most of the
perturbation occurs. We see that as the perturbation goes on, the cluster is
slightly stretched along one of its axis, and compressed along the other two.
Notice that the magnitude of the changes in the cluster shape are
significantly larger for slow perturbations than for the fast ones.

Recall that initially the cluster is located at the origin of the coordinate
system while the PM starts moving in the $y$ direction (with a closest
approach of $b$ for fast shocks, see equation (missing) 1), hence, the long
and middle semi-axes lie in the $x-y$ plane and their directions rotate during
the flyby. Fig. 12 illustrates the changing orientation of the semi-axes as a
function of time for one model. For this plot we computed the eigenvectors of
the shape tensor of the cluster for snapshots before, during, and after the
moment of closest approach. In the beginning of the perturbation the cluster
is slightly stretched but still maintains a nearly spherical shape; then, at
the moment of closest approach, the elongation of the cluster is noticeable,
and keeps increasing after the PM has passed. The degree of deformation of the
cluster is quantified with the axis ratio $a/b$, which increases with time.
Note how the semi-major axis rotates due to the gravitational attraction
exerted on the cluster by the moving PM.

Figure 11: Temporal evolution of the cluster semi-axis lengths for a fast PM
shock ($x=0.03$) and a slow PM tidal perturbation ($W_{0}=4,x=0.62$), both
with $\Delta E_{\mathrm{imp}}/|E_{0}|=0.1$. The dashed black line indicates
the time of maximum approach, $t=0$, while the grey shaded area covers the
time interval [$-\tau,+\tau$], where $\tau$ is the duration of the
perturbation. Figure 12: Temporal evolution of the cluster semi-axis direction
for a PM tidal perturbation ($W_{0}=4,x=0.62$), with $\Delta
E_{\mathrm{imp}}/|E_{0}|=0.1$. The cluster is at the origin while the PM moves
along the $y$ direction, which imprints a rotation to the cluster. We show the
same time interval as in Fig. 11, with $t=0$ being the moment of closest
approach. We indicate the angle between the semi-major axis (red arrow) and
the $x$ axis, as well as the ratio of the semi-lengths $a/b$. Figure 13:
Temporal evolution of the cluster semi-axes length for a fast disc shock
($x=0.03$) and a slow disc tidal perturbation ($x=0.62$), both with $\Delta
E_{\mathrm{imp}}/|E_{0}|=0.1$. The dashed black line indicates the time of
maximum approach, $t=0$, while the grey shaded area covers the time interval
[$-\tau,+\tau$], where $\tau$ is the duration of the perturbation.

On the other hand, in the case of a disc perturbation the cluster moves in the
$z$ direction and crosses the disc located in the $x-y$ plane. Fig. 13 shows
that during the perturbation the $x-$axis and $y-$axis of the cluster remain
unchanged, while it is compressed along its $z-$axis. Again, the change in the
size of the cluster is significant for slow perturbations, but negligible for
the fast ones.

The main difference in cluster response to the PM and disc perturbations is in
the sign of the size change. The remaining bound cluster grows in size when
the perturber is a PM, while it is compressed when the perturber is a disc.

###  4.5 Expected change in cluster shape

We note that the fractional change of the semi-axis length depends not only on
the perturbation strength, but also on its duration.

In this section we characterize the expected change of cluster shape as a
function of $x$ and $\Delta E_{\mathrm{imp}}/|E_{0}|$. In order to do this we
use the $W_{0}=4$ subset from the first row of simulations in Table 1
(perturbation strength 0.1), together with the second and third rows, also
with $W_{0}=4$ with perturbation strengths 0.03 and 0.01. Again, we compute
the length of the primary semi-axes of the cluster as the square roots of the
eigenvalues of the shape tensor (equation (missing) 17) at a time $\tau$ after
the impact.

Figure 14: Expected change in cluster shape for a PM shock. Symbols are the
lengths of the semi-axes measured in our sets of simulations. Their values
scale with the duration of the shock as $x^{3/4}$, and with the shock
amplitude as $(\Delta E_{\mathrm{imp}}/|E_{0}|)^{1/2}$. Dotted lines are the
best fits to the data using equation (missing) 18, which combines these two
dependencies.

Fig. 14 shows the fractional change of the length of the semi-axes as a
function of $x$ and different values of $\Delta E_{\mathrm{imp}}$ for a PM
shock. We notice that such fractional changes scale as a power law of
$\tau/t_{\rm dyn,h}$ for a fixed shock amplitude, and as another power law of
$\Delta E_{\mathrm{imp}}/|E_{0}|$ for a fixed duration of the shock. The power
indices are similar for all axes, and average close to 3/4 and 1/2,
respectively. The scaling with the energy change is expected because $\Delta
v\propto\Delta E_{\mathrm{imp}}^{1/2}$ and we analyse all snapshots at the
same time, such that the amount of extension/compression is proportional to
$\Delta v$. The dependence of $x^{3/4}$ appears to work very well, but we do
not have a physical backing for it and it resulted from trying various
functional forms. The combination of these two dependencies allows us to
predict the shape of the cluster after the shock as follows:

$\frac{l}{l_{0}}=1+\mu_{l}\,x^{3/4}\,(\Delta E_{\mathrm{imp}}/|E_{0}|)^{1/2}$
(18)

where $l=\\{a,b,c\\}$ denotes the three semi-axes. We find the values of
constants $\mu_{a}$, $\mu_{b}$, and $\mu_{c}$ by a linear regression fit to
the simulation results. This gives the best-fit values
$\mu_{l}=\\{1.16,-0.28,-0.68\\}$. Fig. 14 shows how well equation (missing) 18
predicts the change in cluster shape for different values of the adiabatic
parameter and shock amplitude.

It is worth noting that for a PM shock the semi-major axis always grows, while
the other two decrease in size. This is because the tides exerted by a PM are
extensive along the line that connects the cluster with the PM ($x$-axis), and
compressive in the $z$ direction. For fast shocks there is no net tidal action
in the $y$ direction (Spitzer, 1958), as seen in the left of Fig. 11. The
slight compression in $y$ seen in the right panel is probably because the
perturbation is slower, breaking the symmetry because the tides are
compressive in the first half and the cluster has time to respond. The nature
of the tides is reflected in the sign of the constants $\mu_{l}$, and well
illustrated by Fig. 14, with $a/a_{0}>1$, and $b/b_{0},c/c_{0}<1$.

Figure 15: Same as in Fig. 14 but for a disc shock.

We follow the same procedure for the disc case, putting together our $W_{0}=4$
models with shock amplitudes 0.1, 0.03, and 0.01. Fig. 15 shows the length of
the cluster semi-axes as a function of shock duration and shock amplitude.
Notice that in all models the semi-lengths $a$ and $b$ change very little
compared with the change in the short axis $c$. Also, the change of all three
semi-lengths is negative, i.e., the cluster is compressed along its three
principal axes after the shock: $a/a_{0},b/b_{0},c/c_{0}<1$. The expected
shape of the cluster after a disc shock can also be expressed by equation
(missing) 18 but with different constants $\mu_{l}=\\{-0.081,-0.093,-1.12\\}$.
These three constants are negative, and reflect a larger compression in the z
direction, as expected for a disc crossing.

Interestingly, the change in cluster shape described above is independent of
the values of $W_{0}$ explored in this work. The reason for this is probably
because our experiments were set up such that for a certain $\Delta
E_{\mathrm{imp}}/|E_{0}|$, stars at a distance $r=\sqrt{\langle r^{2}\rangle}$
from the cluster centre receive the same velocity kick in clusters with
different $W_{0}$. Because we then measure the shape parameter at that radius,
it is expected to be similar across clusters with different $W_{0}$. The
deformation at larger radii for clusters with larger $W_{0}$ is higher. Hence,
either for a PM or disc shock, the expected shape of the cluster follows the
same functional form (equation (missing) 18), where the sign of the constants
$\mu_{l}$ reflects the nature of the tides along each direction.

###  4.6 Comparison with other $N$-body results

A recent study by Webb et al. (2019) investigated the mass loss of a star
cluster experiencing a sequence of two tidal shocks separated by a time
interval. The authors use the same version of nbody6tt code to model clusters
of 50,000 stars with the Plummer density profile and two King models. They
consider fast shocks and low-density clusters, such that the encounters are
mostly in the impulsive regime ($x<0.25$). The clusters are subject to shocks
with only one non-zero component of the tidal tensor $T_{xx}>0$ that varies as
a step function in time, i.e. a one-dimensional extensive shock. We note that
the associated mass density of this tensor is negative and that realistic
perturbers have several non-zero components of the tidal tensor while the
trace of the tensor is smaller than zero (for extended mass distributions) or
equal to zero (for a point mass).

Webb et al. (2019) measure the changes in cluster energy and mass 10 crossing
times (40 Myr) after the shock. They find an almost exactly linear
proportionality between $\Delta M$ and $\Delta E$ but with a higher
normalization than our equation (missing) 13: $|\Delta M|/M_{0}\approx
0.4\,\Delta E/|E_{0}|$. The larger fractional mass loss is likely because of
their choice of the tidal tensor: applying the tidal forces along a single
axis results in a larger velocity increase of individual stars than when the
tides act along multiple axes (for a given energy increase, see the discussion
in Section 3.3). This is why we find a higher fractional mass loss for the
disc case than for the PM case. Because their tides are extensive, particles
are pushed away from the cluster, which makes it easier to become unbound
compared to our disc case, where particles still need to travel through the
cluster. Another difference between our studies is in the cluster density
profile: the Plummer model is different from the King models. However, Gieles
et al. (2006) applied point-mass perturbations to Plummer models and also
found $|\Delta M|/M_{0}\approx 0.22\,\Delta E/|E_{0}|$. We think that the
larger fractional mass loss of Webb et al. is because of their tidal tensor.
In addition, for such extensive tides the cluster may not have reached full
dynamical equilibrium even after 10 crossing times and the value of $\Delta E$
may still be evolving.

More surprisingly, Webb et al. (2019) find222Note that Webb et al. (2019)
performed the fit with $I_{\mathrm{tid}}$ defined as the square root of our
expression, and we have rescaled their slope for consistency with our
definition. a significantly non-linear relation $\Delta M/M_{0}\propto
I_{\mathrm{tid}}^{\gamma}$ with a shallow slope $\gamma\approx 0.6$. In
contrast, our results in equation (missing) 15 point to a steeper slope
$\gamma>1$ for all but the most concentrated model. They attribute the non-
linearity of this relation to a substantial escape time for newly unbound
stars. While the escape time effect plays a role in decreasing mass loss in
tidally limited clusters (e.g., Lee & Ostriker, 1987; Baumgardt, 2001), it is
less likely to affect isolated clusters studied by Webb et al. (2019). We
therefore conclude that the escape time is not the main cause of the non-
linearity of $\Delta M(I_{\mathrm{tid}})$. Note that in our results the
$\Delta M(\Delta E)$ relation is as non-linear as $\Delta
M(I_{\mathrm{tid}})$, while $\Delta E$ is exactly linearly proportional to
$I_{\mathrm{tid}}$ as expected from the tidal theory.

We also note that in this work we considered a spherical perturber as
infinitely compact PM, appropriate for dense molecular clouds or star-forming
regions. More extended perturbers, such as a galactic bulge or a satellite
galaxy, may lead to more significant adiabatic damping of the energy change
(c.f. Gnedin et al., 1999).

##  5 Summary

We studied the response of a star cluster to different tidal perturbations by
using direct $N$-body simulations. We measured the energy gain of the cluster
as it experiences a PM flyby or a disc crossing. We used different density
profiles given by King (1966) models with dimensionless central concentration
parameters $W_{0}=[2,4,6,8]$. Also, we explored the impulsive and adiabatic
regimes by varying the duration of the encounter.

For fast shocks the predictions from the impulse approximation are recovered
for both perturber types and all cluster concentrations. However, in the
adiabatic regime we find important differences.

For slow disc crossings, we find a smaller energy increase than in the
impulsive regime. This has previously been attributed to the effect of
adiabatic damping (Gnedin & Ostriker, 1999). We show that the energy change
depends on the cluster’s concentration (i.e. $W_{0}$), in a way that is
consistent with adiabatic damping: clusters with less stars in the impulsive
regime have a smaller energy increase. However, this $W_{0}$ dependence in the
energy change can also be attributed to a geometrical distortion (i.e.
compression) in the first half of the perturbation, which reduces the
importance of the remaining part of the perturbation because the cluster is
compressed along the $z$ direction. For similar reasons, slow PM perturbations
lead to much larger energy increases, comparable or even larger than what is
expected in the impulsive regime. Because the PM tidal forces have an
extensive component, this result can be explained by geometrical distortion.
Since we see a similar correlation of $\Delta E$ with $W_{0}$, we conclude
that adiabatic damping does play a role, of similar importance to the
geometrical distortion.

We present an accurate parameterization of the expected changes in cluster
shape as a function of the amplitude of the perturbation and the duration of
the shock (equation (missing) 18).

These results are useful for analytical modeling of cluster disruption for a
sub-grid models of cluster evolution in numerical simulations of galaxy
formation. We conclude that the following steps are necessary to correctly
model the impact of tidal shocks. To compute the energy gain $\Delta E$, it is
important first to identify whether the perturbation is compressive (disc-
like) or has an extensive component (PM-like). For a disc-like perturbation,
$\Delta E$ is computed taking into account the cluster density profile as
expressed in equation (missing) 12 and illustrated in Fig. 2. On the other
hand, for a PM-like perturbation the value of $\Delta E$ is nearly constant
for $x<0.7$, then it increases as a function of $x$ and $W_{0}$.

To understand the mass evolution, the most important assumption is the cluster
density profile, which affects the resulting mass loss (for a given
$I_{\mathrm{tid}}$) by more than an order of magnitude. Semi-analytic models
of cluster evolution with tidal perturbations usually assume a fixed value for
the $W_{0}$ parameter (e.g. Gieles et al., 2006; Elmegreen, 2010; Pfeffer et
al., 2018). To increase the predictive power of such models, we need better
understanding of the evolution of the density profile of clusters evolving in
a Galactic tidal field, experiencing repeated shocks, two-body relaxation and
stellar evolution.

## Acknowledgements

We thank Jeremy Webb for helpful comments on the manuscript, and the referee,
Douglas Heggie, for a careful review and insightful suggestions. LAMM
acknowledges support from PAPIIT IA101520 grant. MG and LAMM are acknowledge
support from the European Research Council (ERC StG-335936, CLUSTERS). MG
acknowledges support from the Ministry of Science and Innovation through a
Europa Excelencia grant (EUR2020-112157). OG is supported in part by NSF
through grant 1909063. HL is supported by NASA through the NASA Hubble
Fellowship grant HST-HF2-51438.001-A awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Incorporated, under NASA contract NAS5-26555.

## Data availability

The data that support the findings of this study are available from the
corresponding author, upon reasonable request.

## References

  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations
  * Baumgardt (2001) Baumgardt H., 2001, MNRAS, 325, 1323
  * Choksi & Gnedin (2019) Choksi N., Gnedin O. Y., 2019, MNRAS, 486, 331
  * Choksi et al. (2018) Choksi N., Gnedin O. Y., Li H., 2018, MNRAS, 480, 2343
  * Elmegreen (2010) Elmegreen B. G., 2010, ApJ, 712, L184
  * Elmegreen & Efremov (1997) Elmegreen B. G., Efremov Y. N., 1997, ApJ, 480, 235
  * Gieles & Renaud (2016) Gieles M., Renaud F., 2016, MNRAS, 463, L103
  * Gieles & Zocchi (2015) Gieles M., Zocchi A., 2015, MNRAS, 454, 576
  * Gieles et al. (2006) Gieles M., Portegies Zwart S. F., Baumgardt H., Athanassoula E., Lamers H. J. G. L. M., Sipior M., Leenaarts J., 2006, MNRAS, 371, 793
  * Gnedin (2003) Gnedin O. Y., 2003, ApJ, 582, 141
  * Gnedin & Ostriker (1999) Gnedin O. Y., Ostriker J. P., 1999, ApJ, 513, 626
  * Gnedin et al. (1999) Gnedin O. Y., Hernquist L., Ostriker J. P., 1999, ApJ, 514, 109
  * Hénon (1971) Hénon M. H., 1971, Ap&SS, 14, 151
  * King (1966) King I. R., 1966, AJ, 71, 64
  * Kruijssen (2015) Kruijssen J. M. D., 2015, MNRAS, 454, 1658
  * Lee & Ostriker (1987) Lee H. M., Ostriker J. P., 1987, ApJ, 322, 123
  * Li & Gnedin (2014) Li H., Gnedin O. Y., 2014, ApJ, 796, 10
  * Li et al. (2017) Li H., Gnedin O. Y., Gnedin N. Y., Meng X., Semenov V. A., Kravtsov A. V., 2017, ApJ, 834, 69
  * Li et al. (2018) Li H., Gnedin O. Y., Gnedin N. Y., 2018, ApJ, 861, 107
  * Li et al. (2019) Li H., Vogelsberger M., Marinacci F., Gnedin O. Y., 2019, MNRAS, 487, 364
  * Murali & Weinberg (1997) Murali C., Weinberg M. D., 1997, MNRAS, 288, 749
  * Muratov & Gnedin (2010) Muratov A. L., Gnedin O. Y., 2010, ApJ, 718, 1266
  * Nitadori & Aarseth (2012) Nitadori K., Aarseth S. J., 2012, MNRAS, 424, 545
  * Ostriker et al. (1972) Ostriker J. P., Spitzer L. J., Chevalier R. A., 1972, ApJ, 176, L51
  * Pfeffer et al. (2018) Pfeffer J., Kruijssen J. M. D., Crain R. A., Bastian N., 2018, MNRAS, 475, 4309
  * Plummer (1911) Plummer H. C., 1911, MNRAS, 71, 460
  * Prieto & Gnedin (2008) Prieto J. L., Gnedin O. Y., 2008, ApJ, 689, 919
  * Renaud et al. (2011) Renaud F., Gieles M., Boily C. M., 2011, MNRAS, 418, 759
  * Spitzer (1958) Spitzer L. J., 1958, ApJ, 127, 17
  * Webb et al. (2019) Webb J. J., Reina-Campos M., Kruijssen J. M. D., 2019, MNRAS, 486, 5879
  * Weinberg (1994a) Weinberg M. D., 1994a, AJ, 108, 1398
  * Weinberg (1994b) Weinberg M. D., 1994b, AJ, 108, 1403
  * Weinberg (1994c) Weinberg M. D., 1994c, AJ, 108, 1414
  * Zemp et al. (2011) Zemp M., Gnedin O. Y., Gnedin N. Y., Kravtsov A. V., 2011, ApJS, 197, 30