Implementation of gravitational shocks in twodimensional FokkerPlanck models
^{1} Department of Astronomy & Space ScienceKyung Hee University, Yongin, 446701 Kyungki, Korea
email: sungsoo.kim@khu.ac.kr
^{2} School of Space Research (WCU program), Kyung Hee University, Yongin, 446701 Kyungki, Korea
^{3} Department of Informational Society Studies, Saitama Institute of Technology, Fukaya, 3690293 Saitama, Japan
Received: 21 December 2011
Accepted: 28 February 2012
We derive analytical formulae for drift and dispersion terms of energy and angular momentum (⟨ΔE⟩, ⟨(ΔE)^{2}⟩, ⟨ΔJ⟩, and ⟨(ΔJ)^{2}⟩) as well as their cross term (⟨ΔEΔJ⟩) for stellar systems under an impulsive perturbation. These terms are expressed as functions of E, J, and orbit averages of powers of radius (r) and the product of radius and velocity (rv), and we confirm our formulae with a numerical simulation. Then with another numerical simulation for a timevarying (Gaussian) perturbation, we find that the adiabatic corrections suggeted by Gnedin and Ostriker can be applied not only to the energy changes (drift and dispersion) but also to the angular momentum changes, if the changes are expressed as functions of energy only. The corrections do not describe the changes accurately when the changes are considered as functions of both energy and angular momentum. The deviations between the numerical simulation and analytical expectations are consierable only in the cluster core though, where the effect of perturbation is relatively weaker. These results are to be implemented in twodimensional (E − J) FokkerPlanck models of the evolution of globular clusters.
Key words: methods: numerical / stars: kinematics and dynamics / globular clusters: general
© 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 Nbody simulations and FokkerPlanck (FP) models. While various physical processes such as binary interactions, tidal fields, and gravitational shocks can be more realistically implemented in Nbody simulations, the computational time required for Nbody simulations is considerably longer than for FP models, and is prohibitively long for N ≳ 10^{6}. For these reasons, Nbody 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 onedimensional (1D) FP equation in energy (E) space or a twodimensional (2D) FP equation in energyangular momentum (E − J) space is considered^{1}. Twodimensional 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 halfanalytic, halfempirical formulae (impulsive approximations with adiabatic corrections) for the average energy change ⟨ ΔE ⟩ and energy dispersion ⟨ (ΔE)^{2} ⟩ (⟨ ⟩ denotes orbitaverages) obtained by Gnedin & Ostriker (1999, GO hereafter). For 2D FP models, not only ⟨ ΔE ⟩ and ⟨ (ΔE)^{2} ⟩ , but also ⟨ΔJ⟩, ⟨(ΔJ)^{2}⟩, and ⟨ ΔEΔJ ⟩ from the shock are required.
In the present paper, we will first derive analytical formulae of ⟨ΔE⟩, ⟨(ΔE)^{2}⟩, ⟨ΔJ⟩, ⟨ (ΔJ)^{2}⟩, 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 timevarying perturbation and see if previously suggested adiabatic corrections can describe the energy and angular momentum changes by a timevarying pertubation as functions of both energy and angular momentum as well. We also generalize our analytic formulae for threedimensional 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) lineofnode (LON) in the x–y plane from the xaxis, 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 )^{2} ⟩
For an impulsive perturbation by a thin sheet of mass moving perpendicularly along the zaxis, 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 firstorder (drift) and secondorder (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 )^{2} ⟩
The first and secondorder 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 firstorder term for the angular momentum: (10)Assuming J^{2} ≫ 2J·ΔJ + (ΔJ)^{2}, 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})^{3} 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^{6} ⟩ /35, and ⟨ r^{6} ⟩ /35, respectively. Inserting these and Eq. (24) into Eq. (23) gives (25)The second term in the RHS here corresponds to (Δv_{z})^{4}, 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/J_{c} 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 (J_{c} 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, E_{bind}. Error bars for the simulation results are shown separately from the data points for the sake of clarity. 

Open with DEXTER 
To test the validity of the above derivations and approximations, we performed an Nbody simulation of a single impulsive shock in a spherical stellar system. We integrated the equation of motion for 10^{6} particles in a fixed cluster potential using the fourthorder Hermite integration scheme (Makino & Aarseth 1992). Following GO, we adopted the King model (King 1966) with the structural parameter W_{0} = 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 halfmass 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 J_{c}. For a given J/J_{c} 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)). 

Open with DEXTER 
Fig. 3 Relative firstorder 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/J_{c} 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, E_{bind}. Error bars for the simulation results are shown separately from the data points for the sake of clarity. 

Open with DEXTER 
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. Timevarying 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 timevarying 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 timevarying 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 nonimpulsive 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 zaxis, the squeeze along the zaxis 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 nonimpulsive shock.
Fig. 4 Same as Fig. 1, but for the timevarying (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/J_{c} 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 (J_{c} 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. 

Open with DEXTER 
Fig. 5 Temporal evolutions of relative energy (upper panels) and angular momentum (lower panels) drifts of particles under a timevarying (Gaussian) perturbation given by Eq. (29). Particles in this simulation have the same energy (E/E_{bind} = −0.7) and the same orbit inclination (60deg) from the x–y plane, but three different angular momenta (J/J_{c} = 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). 

Open with DEXTER 
Fig. 6 Same as Fig. 4, but for angular momentum drift and diffusion. In the right panels, crosses are for J/J_{c} 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 (J_{c} 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)). 

Open with DEXTER 
The exact amounts of ⟨ ΔE ⟩ and ⟨ ΔJ ⟩ from the timevarying perturbation are difficult to analytically estimate, though. We have tried several different timevarying 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^{6} 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 timevarying 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 threedimensional 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 FokkerPlanck models with approproiate adiabatic corrections.
5. Summary
We have derived analytical formulae for ⟨ΔE⟩, ⟨(ΔE)^{2}⟩, ⟨ΔJ⟩, ⟨(ΔJ)^{2}⟩, 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 timevarying (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.
EJ_{z} space is considered instead of E − J when the rotation of the system is important (see Einsel & Spurzem 1999, among others).
Acknowledgments
This research was supported by the Kyung Hee University Research Fund in 2007 (KHU20070775). 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
 Baumgardt, H., & Makino, J. 2003, A&A, 340, 227 [NASA ADS] [CrossRef] (In the text)
 Einsel, C., & Spurzem, R. 1999, MNRAS, 302, 81 [NASA ADS] [CrossRef] (In the text)
 Gnedin, O. Y., & Ostriker, J. P. 1999, ApJ, 513, 626 (GO) [NASA ADS] [CrossRef] (In the text)
 Gnedin, O. Y., Hernquist, L., & Ostriker, J. P. 1999a, ApJ, 514, 109 [NASA ADS] [CrossRef] (In the text)
 Gnedin, O. Y., Lee, H. M., & Ostriker, J. P. 1999b, ApJ, 522, 935 [NASA ADS] [CrossRef] (In the text)
 Shin, J., & Kim, S. S. 2007, J. Korean Astron. Soc., 40, 91 [NASA ADS] [CrossRef] (In the text)
 Shin, J., Kim, S. S., & Takahashi, K. 2008, MNRAS, 386, L67 [NASA ADS] (In the text)
 Spitzer, L., Jr. 1987, Dynamical Evolution of Globular Clusters (Princeton: Princeton Univ. Press) (In the text)
 King, I. R. 1966, AJ, 71, 64 [NASA ADS] [CrossRef] (In the text)
 Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141 [NASA ADS] (In the text)
 Ostriker, J. P., Spitzer, L., Jr., & Chevalier, R. A. 1972, ApJ, 176, L51 [NASA ADS] [CrossRef] (In the text)
 Weinberg, M. D. 1994, AJ, 108, 1398 [NASA ADS] [CrossRef] (In the text)
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/J_{c} 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 (J_{c} 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, E_{bind}. Error bars for the simulation results are shown separately from the data points for the sake of clarity. 

Open with DEXTER  
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 J_{c}. For a given J/J_{c} 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)). 

Open with DEXTER  
In the text 
Fig. 3 Relative firstorder 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/J_{c} 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, E_{bind}. Error bars for the simulation results are shown separately from the data points for the sake of clarity. 

Open with DEXTER  
In the text 
Fig. 4 Same as Fig. 1, but for the timevarying (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/J_{c} 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 (J_{c} 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. 

Open with DEXTER  
In the text 
Fig. 5 Temporal evolutions of relative energy (upper panels) and angular momentum (lower panels) drifts of particles under a timevarying (Gaussian) perturbation given by Eq. (29). Particles in this simulation have the same energy (E/E_{bind} = −0.7) and the same orbit inclination (60deg) from the x–y plane, but three different angular momenta (J/J_{c} = 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). 

Open with DEXTER  
In the text 
Fig. 6 Same as Fig. 4, but for angular momentum drift and diffusion. In the right panels, crosses are for J/J_{c} 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 (J_{c} 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)). 

Open with DEXTER  
In the text 