© ESO, 2012

1. Introduction

Two commonly used methods in studies for dynamical evolution of dense stellar systems such as globular clusters and galactic nuclei are N-body simulations and Fokker-Planck (FP) models. While various physical processes such as binary interactions, tidal fields, and gravitational shocks can be more realistically implemented in N-body simulations, the computational time required for N-body simulations is considerably longer than for FP models, and is prohibitively long for N ≳ 10⁶. For these reasons, N-body simulations are more often chosen for relatively small sets (up to few dozens) of detailed calculations (e.g., Baumgardt & Makino 2003), whereas FP models are preferred for large surveys in statistical studies (e.g., Shin et al. 2008).

The FP model is a statistical method for obtaining the time evolution of a probability density function under the effects of drift and diffusion. For stellar systems, either a one-dimensional (1D) FP equation in energy (E) space or a two-dimensional (2D) FP equation in energy-angular momentum (E − J) space is considered¹. Two-dimensional FP models require much longer computing times and can encounter numerical difficulties more often (see Shin & Kim 2007) than 1D models, but are more realistic since they do not have to assume a velocity isotropy as in 1D models.

Star clusters orbiting in a galaxy can experience gravitational shocks when passing through the galactic disk or bulge (Ostriker et al. 1972). Shocks inject kinetic energy into the cluster and induce a diffusion of stellar energies. Gnedin et al. (1999b) incorporated gravitational shocks into a 1D FP model using half-analytic, half-empirical formulae (impulsive approximations with adiabatic corrections) for the average energy change  ⟨ ΔE ⟩  and energy dispersion  ⟨ (ΔE)² ⟩  (⟨  ⟩  denotes orbit-averages) obtained by Gnedin & Ostriker (1999, GO hereafter). For 2D FP models, not only  ⟨ ΔE ⟩  and  ⟨ (ΔE)² ⟩ , but also  ⟨ΔJ⟩,  ⟨(ΔJ)²⟩, and  ⟨ ΔEΔJ ⟩  from the shock are required.

In the present paper, we will first derive analytical formulae of  ⟨ΔE⟩,  ⟨(ΔE)²⟩,  ⟨ΔJ⟩,  ⟨ (ΔJ)²⟩, and  ⟨ΔEΔJ⟩  as functions of both energy and angular momentum for an impulsive perturbation by a thin sheet of mass, and then perform a numerical simulation to check our derivations. Then we will perform another numerical simulation for a time-varying perturbation and see if previously suggested adiabatic corrections can describe the energy and angular momentum changes by a time-varying pertubation as functions of both energy and angular momentum as well. We also generalize our analytic formulae for three-dimensional shocks such as bulge shocks.

2. Impulsive perturbation

For the orbit averages that will appear throughout the paper, we first define three variables for the location of a star, the inclination of the orbital plane from the x–y plane, i, the angle of the (ascending) line-of-node (LON) in the x–y plane from the x-axis, l, and the phase angle of the star in the orbital plane from the LON, t, such that (1)Note that the unit vector normal to the orbital plane, , has the following relations: (2)When the orbit average is performed for a given set of E and J, i can be treated as a random variable distributed between 0 and π with a weight of sini, and l and t as random variables uniformly distributed between 0 and 2π.

We define one more variable for the velocity, an angle between v and r, p, such that (3)and (4)Note that p is a function of r.

2.1. Energy changes:  ⟨ Δ E  ⟩  and  ⟨ (Δ E )² ⟩ 

For an impulsive perturbation by a thin sheet of mass moving perpendicularly along the z-axis, the change of the vertical velocity of a star at the distance z from the equatorial plane of the cluster can be given by (5)where g_m is the maximum vertical gravitational acceleration by the sheet and V is the perpendicular velocity between the sheet and the cluster (Spitzer 1987). Then, the first-order (drift) and second-order (dispersion) terms for the average energy change of stars are simply given by (6)and (7)With Eqs. (1), (3), and (4), the orbit average of can be evaluated as follows: (8)Here, the averages that do not involve r or v are in fact angle averages, which are subsets of orbit averages. We do not distinguish the two averages throughout this paper. Now, the energy dispersion can be written as (9)

2.2. Angular momentum changes:  ⟨ Δ J  ⟩  and  ⟨ (Δ J )² ⟩ 

The first- and second-order terms for the angular momentum and the cross term between the energy and angular momentum are more complicated to derive than the energy terms since the angular momentum is a vector variable.

Let us start with the first-order term for the angular momentum: (10)Assuming J² ≫ 2J·ΔJ + (ΔJ)², the first term in the right hand side (RHS) can be approximated as (11)using the Taylor expansion. Inserting this into Eq. (10) gives (12)where  ⟨ J·ΔJ ⟩  has vanished by the symmetry.

The first average in the RHS of Eq. (12) can be expressed in terms of spatial variables as (13)using the following relations: (14)The second average in the RHS of Eq. (12) can be written as (15)The three averages here can be written as Now, Eq. (15) simply becomes (19)With Eqs. (13) and (19), one now obtains an angular momentum drift as (20)The dispersion of the angular momentum can be easily obtained from some of the results above. From Eqs. (10), (11), and (19), one finds (21)

2.3. The cross term:  ⟨ Δ E Δ J  ⟩ 

Finally, the cross term between the energy and angular momentum is written as (22)Here, all the terms involving (Δv_z)³ vanished by the symmetry. We use Eq. (2) to obtain (23)Then the first average in the RHS of Eq. (23) now becomes (24)Using Eq. (1), the second and third averages in the RHS of Eq. (23) can be simplified to 2 ⟨ r⁶ ⟩ /35, and  ⟨ r⁶ ⟩ /35, respectively. Inserting these and Eq. (24) into Eq. (23) gives (25)The second term in the RHS here corresponds to (Δv_z)⁴, and we find that this second term is not negligible compared to the first term in the outskirt of the cluster.

2.4. Numerical results

Fig. 1 Relative energy drift (lower panel) and dispersion (upper panel) vs. initial energy of stars after an impulsive perturbation for four different initial angular momentum bins. Crosses are for J/Jc values between 0 and 0.25, triangles for values between 0.25 and 0.5, diamonds for values between 0.5 and 0.75, and squares for values between 0.75 and 1 (Jc is the maximum allowed angular momentum for a given energy). Data points show the energy changes at the end of the simulation. Solid lines are the analytical expectations (Eqs. (6) and (9)). The initial energy in the abscissa is normalized to the binding energy of the cluster, Ebind. Error bars for the simulation results are shown separately from the data points for the sake of clarity.

To test the validity of the above derivations and approximations, we performed an N-body simulation of a single impulsive shock in a spherical stellar system. We integrated the equation of motion for 10⁶ particles in a fixed cluster potential using the fourth-order Hermite integration scheme (Makino & Aarseth 1992). Following GO, we adopted the King model (King 1966) with the structural parameter W₀ = 4 for the potential profile, and the units of G = M = R_c = 1, where G, M, and R_c are the gravitational constant, the total mass of the cluster, and the core radius, respectively. Dynamical timescale at the half-mass radius, t_d,h, of this cluster is  ≈4.5 in the code units, and we chose an integration time step Δt of 0.045.

An impulse of the form (26)was applied to the particles with I_imp = 2g_m/(VΔt) = 1 during one time step at t = 4.5, and the evolution of the cluster was followed until t = 9. The total energy was conserved at the level of ΔE/E ~ 10^-6 during this period.

Figure 1 shows the relative drift and dispersion of the energy as functions of initial energy (normalized to the binding energy of the cluster, E_bind) and angular momentum. Symbols show the values at the end of the simulation, and the solid lines represent the analytical expectations from Eqs. (6) and (9). The figure shows that the energy drift and dispersion from the simulation agree well with the theory within the estimated uncertainties. This agreement was shown by GO as well, but as functions of energy only.

Fig. 2 Relative angular momentum drift (lower panel) and dispersion (upper panel) vs. initial angular momentum of stars after an impulsive perturbation for four different initial energy bins. The initial angular momentum in the abscissa is normalized to Jc. For a given J/Jc bin, stars are binned into four energy bins depending on their initial energies such that the numbers of stars in each energy bin are the same (crosses for the smallest energy bin, then triangles, diamonds, and squares for the largest energy bin). Data points show the angular momentum changes at the end of the simulation. Solid lines are the analytical expectations (Eqs. (20) and (21)).

Fig. 3 Relative first-order cross term between energy and angular momentum vs. initial energy of the stars after an impulsive perturbation for four different initial angular momentum bins. Crosses are for J/Jc values between 0 and 0.25, triangles for values between 0.25 and 0.5, diamonds for values between 0.5 and 0.75, and squares for values between 0.75 and 1. Data points show the energy changes at the end of the simulation. Solid lines are the analytical expectations (Eq. (25)). The initial energy in the abscissa is normalized to the binding energy of the cluster, Ebind. Error bars for the simulation results are shown separately from the data points for the sake of clarity.

The relative drift and dispersion of the angular momentum and the cross term between energy and angular momentum from our simulation are plotted in Figs. 2 and 3. In the former, the abscissa is for the initial angular momentum, instead of energy, and is normalized to the maximum allowed angular momentum for the star’s energy, J_c. These figures show that the simulation result agrees well with our analytical formulae also for angular momentum changes and the cross term (Eqs. (20, 21, 25)).

3. Time-varying perturbation

In the cluster core, rapidly orbiting stars experience a gravitational perturbation over several orbital periods, and the effect of the shocking is not significant because of the conservation of adiabatic invariants. Thus the impulse approximation is not valid and adiabatic corrections are necessary for the cluster core.

Spitzer (1987) found that  ⟨ ΔE ⟩  is exponentially suppressed in the harmonic potential approximation, while Weinberg (1994) asserted that resonances take place in a system with more than one degree of freedom and the effect of shock is stronger than suggested by Spitzer. Using numerical simulations, GO found that for shocks with a characteristic duration, τ, on the order of t_d,h of the cluster, the required adiabatic correction is smaller than that by Weinberg (1994), but larger than that by Spitzer (1987). From their simulations with a fixed cluster potential, GO found that the required adiabatic correction for shocks with τ ~ t_d,h is well described by (27)and (28)where X ≡ ωτ is the adiabatic parameter and ω is the stellar orbital frequency.

We have performed one of the simulations from GO to see if their adiabatic corrections are still valid even in the angular momemtum dimension. Following GO, we applied a time-varying shock of the form (29)which gives the same total energy input as in the impulsive shocking in the previous section. The cluster potential is fixed during the entire shock and τ is set to t_d,h (run C of GO). We found that the simulation for a time-varying perturbation requires higher accuracy than for an impulsive one, so we chose the time step for this simulation to be 0.009, which is five times shorter than that used for the impulsive perturbation.

The relative drift and dispersion of the energy from this simulation are presented in Fig. 4 as a function of energy for all angular momenta (left panels) and for four angular momentum bins (right panels). Also shown are the analytical expectation with the adiabatic corrections by GO (Eq. (28)). The simulation results are reasonably well described by GO’s adiabatic corrections when the energy changes are shown regardless of the angular momentum (left panels). However, significant differences between the simulation and the analytical expectations are seen at low energies when the changes are shown separately for different angular momentum bins. The energy drift even has negative values for orbits with high angular momenta (more circular orbits). This is quite surprising because negative energy drifts were never expected in previous studies, where the energy changes were considered as a function of energy only.

To understand the cause of the negative energy drift, we performed an additional simulation with particles that have the same energy (E/E_bind =  −0.7) and the same orbit inclination (60°) from the x–y plane, but three different angular momenta (R = 0.25, 0.5, and 0.75). The particles in this simulation were set to initially have 36 equally spaced phases in the given orbit and 36 equally spaced right ascension of the apocenter from the line of ascending node (i.e., a total of 1296 particles for each of the three angular momenta), and the Gaussian shock given by Eq. (29) was applied to the particles.

Figure 5 shows that the energy change is indeed a function of angular momentum, with a more negative change for a larger angular momentum. It also shows that the tendency of the energy change follows that of the angular momentum change. We interpret this phenomenon as follows. When a non-impulsive shock is applied, the eccentricy of the orbit of a star can change gradually. Nearly circular orbits (larger J/J_c) will be squeezed in the z direction and become more eccentric (decrease in J/J_c). Highly eccentric orbits (small J/J_c) will have two different responses: when the direction of the orbit elongation is closer to the z-axis, the squeeze along the z-axis will gradually decrease the eccentricity (increase in J/J_c), and when the elongation direction is closer to the x–y plane, the star will spend most of the time near the plane and the change in J/J_c will be small. Then, an increase (decerase) in angular momentum will primarily act on an increase (decrease) in v rather than in r, and thus will increase (decrease) the energy as well. To summarize, orbits with a small J/J_c will have an increase in angular momentum and energy while orbits with a large J/J_c will have a gradual decrease in angular momentum and energy due to a non-impulsive shock.

Fig. 4 Same as Fig. 1, but for the time-varying (Gaussian) perturbation given by Eq. (29). The left panels are for the simulation results and analytical formulae expressed as functions of initial energy only, and the right panels as functions of both initial energy and angular momentum. In the right panels, crosses are for J/Jc values between 0 and 0.25, triangles for values between 0.25 and 0.5, diamonds for values between 0.5 and 0.75, and squares for values between 0.75 and 1 (Jc is the maximum allowed angular momentum for a given energy). Solid lines are the analytical expectations (Eqs. (6) and (9)) with the adiabatic corrections by GO (Eq. (28)). The negative data points are marked with red, thicker symbols.

Fig. 5 Temporal evolutions of relative energy (upper panels) and angular momentum (lower panels) drifts of particles under a time-varying (Gaussian) perturbation given by Eq. (29). Particles in this simulation have the same energy (E/Ebind =  −0.7) and the same orbit inclination (60deg) from the x–y plane, but three different angular momenta (J/Jc = 0.25, 0.5, and 0.75). The particles are set to initially have 36 equally spaced phases in the given orbit and 36 equally spaced right ascension of the apocenter from the line of ascending node (i.e., a total of 1296 particles for each of the three angular momenta).

Fig. 6 Same as Fig. 4, but for angular momentum drift and diffusion. In the right panels, crosses are for J/Jc values between 0 and 0.25, triangles for values between 0.25 and 0.5, diamonds for values between 0.5 and 0.75, and squares for values between 0.75 and 1 (Jc is the maximum allowed angular momentum for a given energy). Solid lines are the analytical expectations (Eqs. (20) and (21)) with the adiabatic corrections by GO (Eq. (28)).

The exact amounts of  ⟨ ΔE ⟩  and  ⟨ ΔJ ⟩  from the time-varying perturbation are difficult to analytically estimate, though. We have tried several different time-varying shocks and found that the dependence of ΔE and ΔJ on angular momentum is quite complicated. Understanding the exact behavior of ΔE and ΔJ as a function of angular momentum requires more detailed studies, and is beyond the scope of this paper.

Figure 6 shows the angular momentum drift and dispersion from our N = 10⁶ simulation with a Gaussian shock along with the GO’s adiabatic corrections. As expected from the discussion above, the angular momentum drifts are negative for the orbits with low J/J_c values, and this is true even for orbits with high energies. As in the case of energy changes, both angular momentum drift and disperion deviate nonnegligibly from the GO’s adiabatic corrections when expressed in terms of both energy and angular momentum (right panels), but they can be reasonably well described by the adiabatic corrections when considered as a function of energy only (left panels).

4. Tidal perturbation by bulges

The impulsive and time-varying perturbations discussed above are specifically derived and calculated for disk shocks. In bulge shocks, velocity changes generally take place in all three dimensions. Here, we present three-dimensional generalizations of the energy and angular momentum changes derived in Sect. 2.

We write the velocity changes as (30)Then, the energy changes are found to be and the angular momentum changes to be where (35)The cross term between the energy and angular momentum up to the fourth order is found to be (36)Gnedin et al. (1999a) presented procedures to obtain the velocity changes (Eq. (30)) of stars in stellar systems that pass a spherically symmetric galaxy (or bulge). Our general formulae above can then be applied to the Fokker-Planck models with approproiate adiabatic corrections.

5. Summary

We have derived analytical formulae for  ⟨ΔE⟩,  ⟨(ΔE)²⟩,  ⟨ΔJ⟩,  ⟨(ΔJ)²⟩, and  ⟨ΔEΔJ⟩  as functions of E, J, and orbit averages of powers of r or rv for an impulsive perturbation, and confirmed the results using a numerical simulation. Then we performed a numerical simulation for a time-varying (Gaussian) perturbation and found that the adiabatic corrections suggested by GO can be applied not only to the energy drift and dispersion but also to the angular momentum drift and disperion, if the energy and angular momentum changes are expressed as functions of energy only. The adiabatic corrections by GO do not describe these changes accurately when the changes are considered as functions of angular momentum as well. The deviations between the numerical simulation and analytical expectations are considerable only in

the cluster core though, where the effect of perturbation is relatively weak.

Acknowledgments

This research was supported by the Kyung Hee University Research Fund in 2007 (KHU-20070775). K.T. thanks Oleg Gnedin for helpful discussion on the analytic derivations of angular momentum changes by gravitational shocks, which initiated our work here.

References

All Figures

Fig. 1 Relative energy drift (lower panel) and dispersion (upper panel) vs. initial energy of stars after an impulsive perturbation for four different initial angular momentum bins. Crosses are for J/Jc values between 0 and 0.25, triangles for values between 0.25 and 0.5, diamonds for values between 0.5 and 0.75, and squares for values between 0.75 and 1 (Jc is the maximum allowed angular momentum for a given energy). Data points show the energy changes at the end of the simulation. Solid lines are the analytical expectations (Eqs. (6) and (9)). The initial energy in the abscissa is normalized to the binding energy of the cluster, Ebind. Error bars for the simulation results are shown separately from the data points for the sake of clarity.

In the text

Fig. 2 Relative angular momentum drift (lower panel) and dispersion (upper panel) vs. initial angular momentum of stars after an impulsive perturbation for four different initial energy bins. The initial angular momentum in the abscissa is normalized to Jc. For a given J/Jc bin, stars are binned into four energy bins depending on their initial energies such that the numbers of stars in each energy bin are the same (crosses for the smallest energy bin, then triangles, diamonds, and squares for the largest energy bin). Data points show the angular momentum changes at the end of the simulation. Solid lines are the analytical expectations (Eqs. (20) and (21)).

In the text

Fig. 3 Relative first-order cross term between energy and angular momentum vs. initial energy of the stars after an impulsive perturbation for four different initial angular momentum bins. Crosses are for J/Jc values between 0 and 0.25, triangles for values between 0.25 and 0.5, diamonds for values between 0.5 and 0.75, and squares for values between 0.75 and 1. Data points show the energy changes at the end of the simulation. Solid lines are the analytical expectations (Eq. (25)). The initial energy in the abscissa is normalized to the binding energy of the cluster, Ebind. Error bars for the simulation results are shown separately from the data points for the sake of clarity.

In the text

Fig. 4 Same as Fig. 1, but for the time-varying (Gaussian) perturbation given by Eq. (29). The left panels are for the simulation results and analytical formulae expressed as functions of initial energy only, and the right panels as functions of both initial energy and angular momentum. In the right panels, crosses are for J/Jc values between 0 and 0.25, triangles for values between 0.25 and 0.5, diamonds for values between 0.5 and 0.75, and squares for values between 0.75 and 1 (Jc is the maximum allowed angular momentum for a given energy). Solid lines are the analytical expectations (Eqs. (6) and (9)) with the adiabatic corrections by GO (Eq. (28)). The negative data points are marked with red, thicker symbols.

In the text

Fig. 5 Temporal evolutions of relative energy (upper panels) and angular momentum (lower panels) drifts of particles under a time-varying (Gaussian) perturbation given by Eq. (29). Particles in this simulation have the same energy (E/Ebind =  −0.7) and the same orbit inclination (60deg) from the x–y plane, but three different angular momenta (J/Jc = 0.25, 0.5, and 0.75). The particles are set to initially have 36 equally spaced phases in the given orbit and 36 equally spaced right ascension of the apocenter from the line of ascending node (i.e., a total of 1296 particles for each of the three angular momenta).

In the text

Fig. 6 Same as Fig. 4, but for angular momentum drift and diffusion. In the right panels, crosses are for J/Jc values between 0 and 0.25, triangles for values between 0.25 and 0.5, diamonds for values between 0.5 and 0.75, and squares for values between 0.75 and 1 (Jc is the maximum allowed angular momentum for a given energy). Solid lines are the analytical expectations (Eqs. (20) and (21)) with the adiabatic corrections by GO (Eq. (28)).

In the text