text:
# Quantum-“classical” correspondence in a nonadiabatic-transition system

## Abstract

A nonadiabatic-transition system which exhibits “quantum chaotic” behavior \[Phys. Rev. E 63, 066221 (2001)\] is investigated from quasi-classical aspects. Since such a system does not have a naive classical limit, we take the mapping approach by Stock and Thoss \[Phys. Rev. Lett. 78, 578 (1997)\] to represent the quasi-classical dynamics of the system. We numerically show that there is a sound correspondence between the quantum chaos and classical chaos for the system.

Nonadiabatic transition (NT) is a very fundamental concept in physics and chemistry Nakamura02 ; TAWM00 . In atomic, molecular, and chemical physics literature, NTs occur as a breakdown of the Born-Oppenheimer (BO) approximation, which is essentially an adiabatic approximation to solve quantum systems with many degrees of freedom. It is still a tough problem to analyze the properties of NTs especially for multidimensional systems. One reason preventing us from deeper understanding of NT is the lack of a naive quantum-classical correspondence for NT like tunneling phenomena UT96 ; SI98 . One way to get the classical picture for a NT system is to go back to the original system before the BO approximation: Fuchigami and Someda investigated dynamical properties of H$`{}_{2}{}^{}{}_{}{}^{+}`$ from classical points of view by treating an electron and nuclei as dynamical variables FS03 . Though there is a full quantum study for such a small system to compare with KKF99 , this “purist” way cannot be easily applied to much more “complex” systems.

The mapping method recently advocated by Stock and Thoss ST97  can circumvent this deficiency. (This is reminiscent of the Meyer-Miller method MM79 .) Their method is as follows: After the BO approximation, the discrete electronic degrees of freedom are mapped onto the Schwinger bosons JJSakurai . Since all the degrees of freedom become just bosons, the total system is rather easily treated semiclassically or quasi-classically. Using this method semiclassically, one can obtain, e.g., absorption spectra even for a pyrazin molecule with 24 degrees of freedom TMS00 . One can use it quasi-classically by solving the equations of motion derived from a mapping Hamiltonian. This is a very easy way to simulate NT systems because the additional number of degrees of freedom for electronic parts is rather small. Using the periodic orbit theory Gutzwiller90  or the adiabatic switching method Reinhardt90 , one can obtain even quantum eigenenergies and eigenstates, in principle ZGE86 .

On the other hand, multidimensional NT systems like Jahn-Teller molecules KDC84  are known to show “quantum chaotic” behavior Gutzwiller90 . Fujisaki and Takatsuka investigated this problem deeply employing the two-mode-two-state (TMTS) system which is considered as a system with two vibrational modes and two electronic states FT01b . They calculated the statistical properties of the eigenenergies and eigenfunction for the TMTS system, and found that the system becomes strongly “quantum chaotic” under a certain condition. In addition, they showed that the chaos is not just a reflection of the lower adiabatic system nor that of the diabatic systems. (On the other hand, the chaos of all previous studies is just a reflection of the lower adiabatic systems KDC84 .) This means that conventional classical descriptions do not help to explain the quantum chaotic behavior. Hence this system deserves to be further studied from the “mapping” (extended classical) points of view. Though there are some studies which investigated chaotic properties of this kind of mixed quantum-classical systems BE94 , our focus here is a quantum-“classical” correspondence (if any) for the TMTS system.

The TMTS system FT01b  first introduced by Heller Heller90  is described by the following Hamiltonian:

$$_{\mathrm{TMTS}}=\left(\begin{array}{cc}T_{\mathrm{kin}}+V_A& J\\ J& T_{\mathrm{kin}}+V_B\end{array}\right),$$
(1)
where $`T_{\mathrm{kin}}`$ is the kinetic energy, $`V_i`$ ($`i=A,B`$) is the potential energy for state $`i`$ defined by

$`T_{\mathrm{kin}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(p_x^2+p_y^2),`$ (2)
$`V_A`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\omega _x^2x^2+\omega _y^2y^2)+ϵ_A,`$ (3)
$`V_B`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\omega _x^2\xi ^2+\omega _y^2\eta ^2)+ϵ_B,`$ (4)
with

$`\xi `$ $`=`$ $`(x+2a\mathrm{sin}\theta )\mathrm{cos}2\theta +(y2a\mathrm{cos}\theta )\mathrm{sin}2\theta ,`$ (5)
$`\eta `$ $`=`$ $`(x+2a\mathrm{sin}\theta )\mathrm{sin}2\theta +(y2a\mathrm{cos}\theta )\mathrm{cos}2\theta .`$ (6)
Note that we have just used a harmonic potential for each state. For the geometrical meaning of the parameters, see Fig. 1.

Here the Duschinsky angle $`\theta `$ TLL03  and the nonadiabatic coupling constant $`J`$ are two important parameters for the system. If these are appropriately chosen, the system becomes strongly “quantum chaotic”, i.e., the nearest-neighbor energy-level spacing distribution becomes the Wigner type and the amplitude distribution of the eigenstates becomes Gaussian FT01b .

The mapping Hamiltonian ST97 ; TMS00  for this system is

$`H_{\mathrm{map}}=T_{\mathrm{kin}}+N_AV_A+N_BV_B+J(x_Ax_B+p_Ap_B)`$ (7)
with

$`N_i={\displaystyle \frac{x_i^2+p_i^2}{2}}\gamma (i=A,B).`$ (8)
Note that $`N_A+N_B=1`$. From this relation, one might consider $`N_i`$ ($`i=A,B`$) as a probability but it is not the case. This is because the numerical range for $`N_i`$ is $`\gamma <N_i<1+\gamma `$. The equations of motion for this Hamiltonian are derived as follows:

$`{\displaystyle \frac{d}{dt}}x_A`$ $`=`$ $`V_Ap_A+Jp_B,{\displaystyle \frac{d}{dt}}p_A=V_Ax_AJx_B,`$ (9)
$`{\displaystyle \frac{d}{dt}}x_B`$ $`=`$ $`V_Bp_B+Jp_A,{\displaystyle \frac{d}{dt}}p_B=V_Bx_BJx_A,`$ (10)
$`{\displaystyle \frac{d}{dt}}x`$ $`=`$ $`p_x,{\displaystyle \frac{d}{dt}}p_x={\displaystyle \frac{V_A}{x}}N_A{\displaystyle \frac{V_B}{x}}N_B,`$ (11)
$`{\displaystyle \frac{d}{dt}}y`$ $`=`$ $`p_y,{\displaystyle \frac{d}{dt}}p_y={\displaystyle \frac{V_A}{y}}N_A{\displaystyle \frac{V_B}{y}}N_B.`$ (12)
As a reference system, we take the lower adiabatic system as in FT01b . The lower adiabatic systems is defined by the following Hamiltonian:

$`H_{\mathrm{ad}}^{}=T_{\mathrm{kin}}+V_{\mathrm{ad}}^{},`$ (13)
with

$`V_{\mathrm{ad}}^{}={\displaystyle \frac{V_A+V_B}{2}}\sqrt{\left({\displaystyle \frac{V_AV_B}{2}}\right)^2+J^2}`$ (14)
Since the TMTS system can be quantum chaotic as mentioned above, we investigate the chaotic properties for the mapping system and the lower adiabatic system. As an indicater of chaos, it is very natural to take Lyapunov exponents Ott02 . We focus on a finite-time maximum Lyapunov exponent calculated as

$$\lambda _{\mathrm{max}}(T)\frac{1}{T}\mathrm{log}\frac{\mathrm{\Delta }d(T)}{\mathrm{\Delta }d(0)}$$
(15)
with

$`\mathrm{\Delta }d(t)^2`$ $`=`$ $`\mathrm{\Delta }\stackrel{~}{x}_A(t)^2+\mathrm{\Delta }\stackrel{~}{p}_A(t)^2+\mathrm{\Delta }\stackrel{~}{x}_B(t)^2+\mathrm{\Delta }\stackrel{~}{p}_B(t)^2`$ (16)
$`+`$ $`\mathrm{\Delta }\stackrel{~}{x}(t)^2+\mathrm{\Delta }\stackrel{~}{p}_x(t)^2+\mathrm{\Delta }\stackrel{~}{y}(t)^2+\mathrm{\Delta }\stackrel{~}{p}_y(t)^2`$
where $`\stackrel{~}{x}_A=x_A/\sqrt{2+2\gamma },\stackrel{~}{p}_A=p_A/\sqrt{2+2\gamma },\stackrel{~}{x}_B=x_A/\sqrt{2\gamma },\stackrel{~}{p}_B=p_A/\sqrt{2\gamma },\stackrel{~}{x}=\omega _xx/\sqrt{2E},\stackrel{~}{y}=\omega _yy/\sqrt{2E},\stackrel{~}{p}_x=p_x/\sqrt{2E},\stackrel{~}{p}_y=p_y/\sqrt{2E},`$ and $`E`$ is the total energy for the system. Here tilde variables are introduced for normalization and $`\mathrm{\Delta }`$ means a distance between a trajectory and its auxiliary one. In the following we take the typical value of $`\gamma `$, i.e., $`\gamma =1/2`$ ST97 . For the numerical calculation of Lyapunov exponents, the method by Benettin et al. BGS76  is employed, i.e., we calculate $`\lambda _{\mathrm{max}}(T)`$ for a finite $`T`$, then shorten the distance $`\mathrm{\Delta }d(T)`$ to $`\mathrm{\Delta }d(0)`$, and run the trajectory again and so on. In this study, we took $`T=24`$ because of the numerical divergence of the $`\mathrm{\Delta }d(T)`$ for larger values of $`T`$.

Following the previous studies FT01b , we concentrate on a rather high energy region around $`E=E_0=28`$. (This energy is much higher than that around the crossing seam region.) Varying the Duschinsky angle $`\theta `$ and the nonadiabatic coupling constant $`J`$, we calculate the distribution of the finite-time maximum Lyapunov exponent $`\lambda _{\mathrm{max}}(T)`$. For simplicity, hereafter, we call $`\lambda _{\mathrm{max}}(T)`$ just as a Lyapunov exponent. Since the phase space has a structure especially for the lower adiabatic system, globally averaged Lyapunov exponents are not so useful. Instead, we investigate the properties of the distribution of the Lyapunov exponents, which reflects phase space structure of the system. (Remember that the Berry-Robnik distribution, which reflects the phase space volume of chaotic seas, is useful for mixed systems MHA01 .) To this end, we prepare an initial ensemble of particles which are sampled from a part of the equi-energy potential surface $`V_A(x,y)=E_0`$ with constraints $`p_x=p_y=0`$ and $`y<x\mathrm{tan}\theta +a/\mathrm{cos}\theta `$. (The latter constraint means that we only take points below the crossing seam. See Fig. 1.) We take 40 sample points from this curve and calculate the histogram for the Lyapunov exponents. We believe that these sample points represent a typical situation of the TMTS system because, at least, the characteristic of the lower adiabatic system can be understood from these sampling points FT01b .

First we fix $`\theta =\pi /3`$, and investigate the $`J`$ dependence of the chaotic properties. Let us summarize the corresponding quantum results: the nearest-neighbor (energy-level) spacing distribution is Wigner for $`J=1.5`$ (3rd row of Fig. 3 in FT01b ), whereas it is rather Poisson for $`J=7.5`$ (3rd row of Fig. 4 in FT01b ), and rather mixed for $`J=0.3`$ (3rd row of Fig. 2 in FT01b ). Figure 2 shows the distribution of the Lyapunov exponents for the TMTS system. The distribution for $`J=1.5`$ has a sharp peak around 1, whereas that for $`J=7.5`$ has a rather broad peak around 0.4, and that for $`J=0.3`$ is also a little bit broad. This corresponds to the quantum results, at least, qualitatively. On the other hand, Fig. 3 shows the distribution of the Lyapunov exponents for the lower adiabatic system. As one can see, the values themselves are much smaller than those for the TMTS system and it is difficult to distinguish three distributions. This also corresponds to the quantum mechanical calculation for the lower adiabatic system (Fig. 9 in FT01b ). From this comparison, it is reasonable to conclude that there is a sound quantum-“classical” correspondence between the TMTS system and the mapped system in view of their “chaotic” properties.

Next we investigate the $`\theta `$ dependence of the chaotic properties while fixing $`J=1.5`$. Although there is a strong peak around $`\lambda _{\mathrm{max}}(T)1.3`$ as shown in Fig. 4 (a), we can see that the system with $`\theta =0`$ is not globally chaotic. This is because $`x`$ does not effectively coupled to $`y`$ when $`\theta =0`$, and the motion along $`x`$ axis is regular. (On the other hand, even with $`\theta =0`$, the motion along $`y`$ axis can be chaotic as shown by the peak of the distribution around 1.3.) In such a case, we do not expect that the Wigner type distribution arises in the corresponding quantum system, and this is the case for the TMTS system (1st row of Fig. 3 in FT01b ). On the other hand, for intermediate Duschinsky angles ($`\theta =\pi /6,\pi /3`$), the Lyapunov exponent distributions show that the system is globally chaotic \[Fig. 4 (b),(c)\], and the corresponding quantum system can have the Wiger type distribution, which is also confirmed numerically (2nd and 3rd rows of Fig. 3 in FT01b ). Of course, we have to admit that this correspondence is loosely stated, and there remains a difficult question exactly when the quantum chaos begins in the parameter space. Such an issue must be addressed utilizing semiclassical methods ST97 ; Gutzwiller90  or the adiabatic switching method Reinhardt90 . It is also interesting to investigate phase space structure of this mapping system and to relate it to the nearest-neighbor spacing distribution via the Berry-Robnik distribution MHA01 . We also hope that this study will cast a light on the relation between the statistical reaction theory for NT systems and Lyapunov spectra for them.

In this paper, employing the mapping approach by Stock and Thoss, we investigated a nonadiabatic-transition system which exhibits “quantum chaotic” behavior from quasi-classical aspects. By comparing the statistical properties of the quantum system with the Lyapunov exponent distributions of the mapping system, we found that there is a sound quantum-“classical” correspondence in the system.

The author thanks T. Takami for comments.


format: text

text:
# Absolute frequency measurements of the 𝐷₂ line and fine-structure interval in 39K

## I Introduction

Precise measurements of atomic energy levels continue to play an important role in the development of physics. The energy levels of alkali atoms are particularly important because these atoms can be laser cooled to ultra-low temperatures for subsequent use in precision spectroscopy experiments. For example, precise measurements of the $`D_1`$ line in Cs URH99 , Rb or K BDN03a , in conjunction with photon-recoil shift measurements in an atom interferometer WYC93 , could yield an independent value of the fine-structure constant $`\alpha `$. In addition, a precise value of the frequency of the $`D_2`$ line in Cs URH00  is required for atom-interferometric measurements of the local gravitational acceleration PCC99 .

We have recently developed a technique for measuring the absolute frequencies of optical transitions with sub-MHz precision BDN03 ; BRD03 . The frequency is measured using a ring-cavity resonator whose length is calibrated against a reference laser locked to the $`D_2`$ line of 87Rb. The frequency of the reference laser is known to an accuracy of 10 kHz YSJ96 . We have already used this technique to measure the $`D_1`$ lines in the alkali atoms 39K, 85Rb, and 87Rb with 0.13 ppb uncertainty BDN03a , to facilitate measurements of $`\alpha `$. In this paper, we apply this technique to measure the $`D_2`$ line of 39K with an uncertainty of 0.25 ppb. This represents an improvement of more than two orders of magnitude over tabulated values NIST . Furthermore, in combination with our previous measurement of the $`D_1`$ line, we obtain the fine-structure interval in the $`4P`$ state of 39K with very high precision. Knowledge of fine-structure intervals is useful in the study of atomic collisions and relativistic calculations of atomic energy levels. Indeed, a combination of improved methods of calculation and increasingly precise measurements of atomic fine and hyperfine structure (e.g. in Li WAC03 ) can lead to a precise value of $`\alpha `$ FHH80 .

The $`D_2`$ line in K presents a problem in determining the line center because transitions to individual hyperfine levels of the excited state are not resolved in conventional saturated-absorption spectroscopy. This is because the various hyperfine levels in the $`4P_{3/2}`$ state lie within 30 MHz of each other AIV77 , while the natural linewidth is 6 MHz. We have recently demonstrated a technique to resolve such closely-spaced transitions BAN03 , however, in this work we use computer simulation of the measured spectrum to extract the line center with an accuracy of $``$100 kHz.

## II Experimental details

The experimental schematic is shown in Fig. 1 and has been described extensively in a previous publication BDN03a . The four-mirror ring-cavity is placed inside a vacuum chamber and evacuated to a pressure of $``$10-2 torr to eliminate frequency shifts due to dispersion of air. Lasers 1 and 2 are standard external-cavity diode lasers stabilized using optical feedback from a piezo-mounted grating BRW01 . The reference laser (Laser1) is locked to the $`D_2`$ line of 87Rb using saturated-absorption spectroscopy in a vapor cell, and the cavity is locked to the reference laser. In general, the frequency of the laser to be measured (Laser2) will be offset from the nearest cavity resonance. This offset is compensated by using an acousto-optic modulator (AOM) between Laser2 and the cavity. The AOM is locked to this frequency difference and its frequency is read using a counter. The exact mode number of the cavity is determined by measuring the cavity’s free-spectral range using different lock points of the reference laser, as described in Ref. BDN03 . The mode number and the AOM offset together yield the frequency of Laser2. A second AOM kept between the reference laser and its saturated-absorption spectrometer allows us to vary the reference frequency continuously so that the AOM offset measured for Laser2 is always close to a given value. This eliminates potential systematic errors due to changes in the direction of the beam entering the cavity.

The error signals needed for locking the diode lasers are produced by modulating the injection current at a frequency of 20–50 kHz. The error signal is obtained from the saturated-absorption signal by phase-sensitive detection at the third harmonic of the modulation frequency WAL72 . This produces narrow dispersive signals that are free from effects due to the underlying Doppler profile or intensity fluctuations. The error signal for locking the cavity is only a first-derivative signal since the cavity modes appear against a flat background.

The spectroscopy on 39K is performed in an ultrahigh vacuum glass cell maintained at a pressure below $`10^8`$ torr by an ion pump. K vapor is produced by heating a getter source SAES  with a current of 2.6 A. The ultrahigh vacuum environment is necessary to minimize linewidth broadening arising as a result of background collisions. The getter source also gives us control over the amount of K vapor in the cell, which we optimized to obtain the narrowest linewidth.

Fig. 2 shows a typical saturated-absorption spectrum of the $`D_2`$ line covering the two ground-hyperfine levels. Each peak is actually a convolution of six peaks. However, the individual hyperfine transitions are not resolved because the different hyperfine levels of the excited state lie within 30 MHz of each other. To determine the line center, we therefore performed computer simulations of the spectrum. We first fixed the locations of the 6 peaks (with respect to the line center) according to the known hyperfine shifts AIV77 . We then set the peak amplitudes to correspond to those obtained in saturated-absorption spectroscopy. We further assumed that the linewidth of all the transitions is the same. We calculated the spectrum for a given value of the line center and linewidth.

Fig. 3 shows a close-up of the measured spectrum (open circles) for $`F=2F^{}`$ transitions. The solid curve is a best fit to the spectrum, obtained with a linewidth of 12.3 MHz. In order to check that this linewidth is reasonable, we tuned the diode laser to the $`D_1`$ line of K (at 770 nm), where the individual hyperfine transitions are clearly resolved. The linewidth obtained for those transitions was 14 MHz, close to the fit value for the $`D_2`$ line. To check the linewidth further, we changed the linewidth for the calculated spectrum and held it constant during the fitting. The best-fit calculated spectra for a linewidth of 15 MHz is also shown in Fig. 3. Even with this small change, the lineshape deviates from the measured data. The fit linewidth of 12.3 MHz, however, is larger than the natural linewidth of 6 MHz. The primary causes for this increase are power broadening due to the pump beam and a small angle between the counterpropagating pump and probe beams. The effect of stray magnetic fields and background collisions is negligible.

With the best fit to the spectrum in Fig. 3, we found that its maximum lies 2.60(10) MHz above the line center. Similar modelling for transitions starting from the $`F=1`$ ground level shows that the maximum in the spectrum lies 12.03(10) MHz below the line center.

## III Error analysis

The errors in our frequency measurement technique have been discussed extensively in earlier publications BDN03a ; BDN03 ; BRD03 . We present here a brief overview for the sake of completeness. There are two classes of potential systematic errors that we consider. The first class of errors comes from systematic shifts in the laser frequencies. For the reference laser, changes in the lineshape of the peaks due to optical pumping effects are taken care of by carefully adjusting the pump and probe beam intensities in the saturated-absorption spectrometer (to a ratio of about 3). Shift in its lock point due to peak pulling from neighboring transitions, the underlying Doppler profile, or phase shifts in the feedback loop are minimized by third-harmonic detection for the error signal. Collisional shifts in the Rb vapor cell and the effect of stray magnetic fields are negligible.

The second class of systematic errors is inherent to our technique because we are really comparing the wavelength (and not the frequency) of the two lasers. The most important source of error is dispersion inside the cavity, which is eliminated by using an evacuated cavity. However, there could be wavelength-dependent phase shifts at the dielectric coated mirrors used in the cavity. Such errors can be corrected by repeating the measurement at different cavity lengths. We have shown earlier BDN03a  that this error is negligible when the unknown laser differs from the reference laser by up to 25 nm. In the current work, the wavelength difference was only 13 nm, therefore the measurements were done at a single cavity length of $``$178 mm.

## IV Results

The results of our measurements are listed in Table I. Each value is an average of 50–100 individual measurements, and the quoted error includes our estimate of the systematic error. To check for long-term variations, we repeated these measurements over a period of several weeks. In addition, we used two different lock points of the reference laser, namely the $`F=2F^{}=(2,3)`$ transition and the $`F=1F^{}=(1,2)`$ transition. The frequencies for these transitions differ by 6622.886 MHz. The consistency of our values for the two cases acts as a check on our ability to determine the correct cavity mode number since the cavity free-spectral range is only about 1.7 GHz. Another check on our results is that the measured frequencies for the two transitions should differ by the ground hyperfine splitting in 39K. Our value of 461.569(210) MHz overlaps very well with the accepted value of 461.720 MHz AIV77 .

The values in Table I can be combined to yield the hyperfine-free frequency of the $`D_2`$ line in K

$`4P_{3/2}4S_{1/2}`$ : 391 015 578.040(75) MHz

This value is consistent with the value of 391 015 585(30) MHz from the National Institute of Standards and Technology energy-level tables NIST , but the accuracy is improved by more than two orders of magnitude. The above result can be combined with our earlier measurement of the $`D_1`$ line BDN03a  to yield the fine-structure interval in the $`D`$ line in K

$`4P_{3/2}4P_{1/2}`$ : 1 729 997.132(90) MHz

This value overlaps at the $`1\sigma `$ level with our earlier measurement of the fine-structure interval using a home-built wavemeter BDN04 : 1 729 993 (5) MHz, but again the accuracy is improved considerably.

## V Conclusion

In conclusion, we have measured the $`D_2`$-line frequency and fine-structure interval in $`{}_{}{}^{39}\mathrm{K}`$ with uncertainty less than 100 kHz. The measurement is complicated by the fact that the closely-spaced hyperfine transitions in the line are not resolved in conventional saturated-absorption spectroscopy. We use computer simulation of the measured spectrum to extract the line center. The frequency is measured using a Rb-stabilized ring-cavity resonator that has general applicability to precise frequency measurement of optical transitions BDN03 . The only other technique of measuring optical frequencies with comparable precision is the frequency-comb method URH99 . By comparison, our technique is low-cost and easier to implement. In the future, we plan to measure the $`D`$ lines and fine-structure interval of other alkali atoms such as Li, Na, and Cs. These measurements could play an important role in QED-independent determinations of the fine-structure constant $`\alpha `$.

This work was supported by the Board of Research in Nuclear Sciences (DAE), and the Department of Science and Technology, Government of India.


text:
# Dynamical friction in flattened systems: A numerical test of Binney’s approach

## 1 Introduction

According to current ideas galaxy formation began with small-amplitude Gaussian fluctuations at the early stages of the Universe. In hierarchical cosmological models, these fluctuations decrease with increasing scales, resulting first in the formation of low-mass objects that may merge, building up ever more massive structures. The shape and morphology of these objects are strongly dependent on the cosmological models, as one can conclude from N-body computations, although none of them predict spherical structures. The most successful hierarchical theory is the Cold Dark Matter model (CDM). In this framework, aspherical bound dark matter haloes (DMHs) form as a result of gravitational clustering. The inclusion of gas dynamics in the CDM simulations (Udry & Martinet 1994, Dubinsky 1994) results to a Gaussian distribution of DMH density aspect ratios, $`q_hc/a>0`$, where $`c`$ and $`a`$ are the minor and major axes of an oblate spheroid, of mean $`<q_h>=1/2`$ and dispersion equal to 0.15 (Dubinsky 1994). The degree of asphericity depends on the dark matter nature. For example, numerical computations based on the $`\mathrm{\Lambda }`$CDM predict halo axis-ratios of $`q_h0.7`$ (Bullock 2001), Hot Dark Matter models lead to haloes as round as $`q_h=0.8`$ (Peebles 1993), whereas candidates such as cold molecular gas (Pfenniger, Combes & Martinet 1994) and massive decaying neutrinos (Sciama 1990) may produce halo profiles as flattened as $`q_h=0.2`$.

Observationally, measuring galaxy axis-ratio is complicated: (i) Stellar kinematics: Olling & Merrifield (2000) obtain an axis-ratio of $`q_h0.8`$ for our Galaxy. This method has the disadvantage of having access to information of our Galaxy only on small scales. (ii) The flying gas layer method (Olling 1996, Becquaert, Combes & Viallefond 1997): assuming that the HI emission comes from gas in hydrostatic equilibrium in galactic potentials axis-ratios as low as $`q_h0.3`$ are obtained for the galaxies NGC 891 and 4244, (iii) Warping gas layer: Hofner & Sparke (1994) find axis-ratios of approximately 0.7 for NGC 2903 and of $`q_h0.9`$ for NGC 2841, 3198, 4565 and 4013, (iv) X-ray isophotes: Boute & Canizares (1998) measure values of $`q_h0.5`$ for NGC 3923, 1332 and 720, (v) Polar ring galaxies: Arnaboldi et al. (1993) and Sackett et al. (1994) find an axis-ratio of $`q_h0.3`$ for NGC 4650A, 0.5 for the galaxy A0136-0801 and 0.6 for AM2020-504, (vi) Precessing dusty discs: Steinman-Cameron, Kormendy & Durisen (1992) measure an axis-ratio of 0.9 for the galaxy NGC 4753.

Another method, which we focus on here, is the analysis of satellite dynamics. There are two main different approaches to infer the halo shape from satellites. First, one may attempt to reproduce the observed tidal streams of Milky Way satellites as done, for instance, by Ibata et al. (2001) and Law et al. (2003) who, from measurements of velocity, position and structure of the Sagittarius dwarf galaxy, constrain the initial parameter space and, subsequently, calculate in detail the satellite mass loss. They find that the Milky Way halo potential cannot be more flattened than $`0.8`$, otherwise tidal streams would be too spread out and thick compared to the observations due to orbital precession. The second approach is a statistical study of satellite distribution around spiral galaxies. Holmberg (1969), Zaritsky & González (1999), Prada et al. (2003) and Sales & Lambas (2003) point out that satellites around disc galaxies are found more often aligned with the poles of the host galaxy, the so-called ’Holmberg effect’, whereas Quinn & Goodmann (1986) find in their $`N`$-body study that the disc alone cannot account for the original statistical distribution of Holmberg’s data. A remedy may be sought in the form of an extended non-spherical DMH. An anisotropic velocity (and mass) distribution will cause a satellite’s orbit to align with the axes of the velocity ellipsoid of the host galaxy (Peñarrubia, Kroupa & Boily 2001, hereinafter PKB). In both schemes, a large number of numerical calculations is needed. In the former approach one should integrate several initial orbital parameters to find the best fit to the observed satellite characteristics, whereas in the latter approach an initial satellite sample should statistically reproduce the distributions expected from cosmological models of DMH formation. So far this is prohibitively time-expensive by means of any of the present N-body algorithms.

Several studies of satellite decay have shown that, in spherical systems, Chandrasekhar’s dynamical friction (Chandrasekhar 1943) is accurate enough if the Coulomb logarithm remains as a free parameter to fit to the N-body data (e.g. van den Bosch et al. 1999, Colpi et al. 1999) since it also depends on the code parameters and the number of particles employed. For instance, Prugniel & Combes (1992) and Whade & Donner (1996) find that dynamical friction is artificially increased due to numerical noise if the particle number is small. Semi-analytic methods that include Chandrasekhar’s dynamical friction have been demonstrated to reproduce accurately the overall evolution of satellite galaxies (e.g Velázquez & White 1999, Taylor & Babul 2001) and, therefore, represent a useful tool in order to carry out extensive studies with a large parameter space.
It is, however, still unclear how the inhomogeneity of the system distribution affects the satellite dynamics. Whereas Del Popolo (2003), Del Popolo & Gambera (1998) and Maoz (1993) show that dynamical friction increases the steeper the density profile is, Just & Peñarrubia (2002), following the theoretical scheme proposed by Binney (1977), hereinafter B77, find a negligible effect on the satellite orbit due to the symmetry of the inhomogeneous terms of dynamical friction.

Although Chandrasekhar’s formula, which assumes an isotropic velocity distribution, reproduces accurately dynamical friction in spherical systems, an analytical study of Statler (1991) shows that in the case of Stäckel potentials the velocity anisotropy produce strong effects on the satellite orbit. The N-body computations of PKB confirm that Chandrasekhar’s formula cannot account for the resulting satellite decay and evolution if the halo is flattened. The aim of this paper is to implement a semi-analytic scheme capable of tracking the dynamical evolution of substructures within flattened as well as spherical DMH’s. With this purpose in mind, we implement in our code the analytic expressions of dynamical friction in systems with anisotropic velocity dispersions suggested by Binney (B77) and also used by Statler, which reproduces Chandrasekhar’s for null anisotropy in velocity space. We as well compare the results of using Chandrasekhar’s formula in axi-symmetric systems to determine the effects of the velocity anisotropy on satellite decay.

Section 2 introduces the models. In Section 3 we provide the code and galaxy parameters. We outline the dynamical friction approaches in Section 4 whereas in Section 6 we propose a simple technique to fit a free parameter (in our case the Coulomb logarithm) to the N-body data. In Section 7 we study how flattened DMHs affect satellite decay and calculate the degree of accuracy of Binney’s formula. The paper concludes with Section 8.

## 2 Galaxy and satellite parameters

Our DMH model is that used by PKB to facilitate an inter-comparison of the disc and bulge effects on satellite decay (Peñarrubia 2003). Here we do not add the bulge and disc components in order to distill the dependency of dynamical friction on the velocity anisotropy in a spheroidal DMH.

In order to minimise computational time when constructing flattened DMHs, we apply a highly-efficient technique using multi-pole potential expansions to tailor the local velocity ellipsoid to the required morphology (Boily, Kroupa & Peñarrubia 2001). The MaGalie code scales linearly with particle number and hence we can construct flattened DMHs consisting of $`\begin{array}{c}>\hfill \\ \hfill \end{array}10^6`$ particles or more, in a short computational time.

Following PKB, the flattened DMH is described by a non-singular isothermal profile which, in cylindrical coordinates (for $`m^2[u=0]`$ see eq. 7), can be described as

$`\rho _h[m^2(0)]={\displaystyle \frac{M_h\alpha }{2\pi ^{3/2}r_{\mathrm{cut}}}}{\displaystyle \frac{\mathrm{exp}\left[m^2(0)/r_{\mathrm{cut}}^2\right]}{m^2(0)+\gamma ^2}},`$ (1)
$`\text{with}m^2(0)=R^2+z^2/q_h^2`$
$`M_h`$ being the DMH mass, $`r_{\mathrm{cut}}`$ the cut-off radius and $`\gamma `$ the core radius, and

$`\alpha \{1\sqrt{\pi }\beta \mathrm{exp}(\beta ^2)[1\mathrm{erf}(\beta )]\}^1=`$ (2)
$`1+\sqrt{\pi }\beta +(\pi 2)\beta ^2+O(\beta ^3)`$
where $`\beta =\gamma /r_{\mathrm{cut}}\begin{array}{c}<\hfill \\ \hfill \end{array}1/24`$ in our calculations. For $`\beta =1/24`$ we find $`\alpha 1.0761`$ already and hence thereafter we set $`\alpha =1`$ in our analysis.

We use self-consistent King models (King 1966) to represent our dwarf galaxies. These models fit early-type dwarf galaxies (Binggeli et al. 1984). For a comparison with the work of PKB we adopt $`c=\mathrm{log}_{10}(r_t/r_c)=0.8`$, where $`r_c`$ and $`r_t`$ are the core and tidal radii, respectively.

To construct the models we choose the satellite mass $`M_s`$ and $`r_t`$. The tidal radius is determined by computing the density contrast, $`\rho _s(r_t)/\overline{\rho _g}(r_a)3`$, at the apo-centric distance ($`r_a=55`$ kpc) at $`t=0`$, $`\overline{\rho _g}(r)`$ being the averaged density of the galaxy (same procedure as Velázquez & White 1999). This guarantees that all satellite particles are bound at $`t=0`$.

We employ the system of units of PBK, which refers to the parameters of the Milky Way disc. Adopting $`M_d=R_d=G=1`$ and according to Bahcall, Smith & Soneira (1982), $`M_d=5.6\times 10^{10}M_{}`$ and $`R_d=3.5`$ kpc, the time and velocity units are, respectively, $`1.3\times 10^7`$ yr and $`262\mathrm{kms}^1`$.

The values of the galaxy and the satellite parameters can be found in Table 1. The parent galaxy corresponds to the model G2 of PKB, with the difference that we remove the disc and bulge components here.

## 3 Numerical calculations

### 3.1 Code parameters

The numerical experiments were carried out by using the code Superbox which is a highly efficient particle mesh-algorithm based on a leap-frog scheme (for a detailed description see Fellhauer et al. 2000). Superbox has already been implemented in extensive studies of satellite disruption by Kroupa (1997), Klessen & Kroupa (1998) and PKB. The program calculates the accelerations using a high order NGP (‘nearest grid point’) force calculation scheme based on the second derivatives of the potential. A self-consistent system of several galaxies can be treated by forming sub-grids (3D ’boxes’) which follow the motion of each galaxy. Each sub-grid has three levels of resolution, the two finest levels co-move with the galaxies allowing a high-resolution calculation of the forces acting on the particles, whereas the third one covers the local ’Universe’. The finest levels are centred on the density maximum of the galaxy, which is re-computed at every time-step.

The code parameters are those of PKB. In that paper a detailed description of the system and the grid structure is presented, whereas here we merely give a brief description of the N-body parameters. Our integration time step is $`0.39`$ Myr which is about $`1/40`$th the dynamical time of our satellite. We have three resolution zones, each with $`64^3`$ grid-cells: (i) The inner grid covers out to 3 radial halo scale-lengths, providing a resolution of 350 pc per grid-cell. (ii) The middle grid covers the whole galaxy, with an extension of 24 halo scale-lengths (84 kpc), giving a resolution of 2.8 kpc per grid-cell. The satellite always orbits within this grid except at the very late stages of its evolution, avoiding cross-border effects (see Fellhauer et al. 2000) and (iii) The outermost grid extends to 348 kpc and contains the local universe, at a resolution of 11.6 kpc. As for the satellite grid-structure, the resolutions are 816 pc per grid-cell for the inner grid that extends to 24.48 kpc, 1.2 kpc per grid-cell for the middle grid which extends to 36 kpc, and 11.6 kpc per grid-cell for the outermost grid that covers the local universe.

The selection of grid parameters ensures the conservation of energy and angular momentum for satellites in isolation over times as long as our calculations to better than 1% for all the models.

### 3.2 Orbital parameters

The parameter space of satellite galaxies is extremely large. A complete survey should account for different satellite masses, apo-galacticon distances, orbital eccentricities and inclinations…etc. In this paper, we carry out a set of calculations selecting those parameters that best reflect the effects of the velocity dispersion anisotropy on satellite decay. These parameters are: (i) the initial orbital inclination ($`i`$), defined as the angle between the initial angular momentum vector of a satellite and the axis perpendicular to the axi-symmetry plane (selected as the $`z`$-axis). We expect orbital inclination to decrease in time as predicted by Binney. We note that all the calculations proceed with the same orbital sense, but this is not important since the halo is non-rotating. (ii) The satellite’s initial orbital eccentricity, defined as $`e=(r_ar_p)/(r_a+r_p)`$, where $`r_a,r_p`$ are the apo- and peri-galacticon, respectively.

The parameters of the numerical experiments are listed in Table 2.

## 4 Halo Dynamical Friction

Chandrasekhar’s expression cannot explain some effects observed in N-body calculations of satellite decay within flattened haloes (PKB). Our aim is to check Binney’s approximation (B77) for systems with anisotropic velocity dispersion.

For simplicity, we reproduce here the analytic formulæ employed throughout this study (for a detailed analysis of the friction force see Peñarrubia 2003 and Just & Peñarrubia 2002).
If the distribution function in velocity space is axi-symmetric, the zeroth order specific friction force is (B77)

$`F_i={\displaystyle \frac{2\sqrt{2\pi }\rho _h[m^2(0)]G^2M_s\sqrt{1e_v^2}\mathrm{ln}\mathrm{\Lambda }}{\sigma _R^2\sigma _z}}B_Rv_i,`$ (3)
$`F_z={\displaystyle \frac{2\sqrt{2\pi }\rho _h[m^2(0)]G^2M_s\sqrt{1e_v^2}\mathrm{ln}\mathrm{\Lambda }}{\sigma _R^2\sigma _z}}B_zv_z,`$
where $`i=x,y`$ and $`(\sigma _R,\sigma _z)`$ is the velocity dispersion ellipsoid of a Schwarzschild distribution in cylindrical coordinates with constant ellipticity $`e_v^2=1(\sigma _z/\sigma _R)^2`$, $`\mathrm{ln}\mathrm{\Lambda }`$ is the Coulomb logarithm of the halo and

$$B_R=_0^{\mathrm{}}𝑑q\frac{\mathrm{exp}(\frac{\mathrm{v}_\mathrm{R}^2/2\sigma _\mathrm{R}^2}{1+\mathrm{q}}\frac{\mathrm{v}_\mathrm{z}^2/2\sigma _\mathrm{R}^2}{1\mathrm{e}_\mathrm{v}^2+\mathrm{q}})}{(1+q)^2(1e_v^2+q)^{1/2}},$$

$$B_z=_0^{\mathrm{}}𝑑q\frac{\mathrm{exp}(\frac{\mathrm{v}_\mathrm{R}^2/2\sigma _\mathrm{R}^2}{1+\mathrm{q}}\frac{\mathrm{v}_\mathrm{z}^2/2\sigma _\mathrm{R}^2}{1\mathrm{e}_\mathrm{v}^2+\mathrm{q}})}{(1+q)(1e_v^2+q)^{3/2}},$$

where $`(v_R,v_z)`$ are the coordinates of the satellite velocity in this frame.

As Binney shows, a body with mass $`M_s`$ will suffer a decrease of its orbital inclination whenever $`B_z>B_R`$ (oblate halo). If the orbit is either coplanar or polar, the inclination remains constant since, respectively, the perpendicular and the planar component of $`𝐯`$ is zero. It is straight-forward to show that this expression reproduces Chandrasekhar’s for $`e_v=0`$,

$$𝐅_{\mathrm{ch}}=4\pi GM_s\rho _h[m^2(0)]\mathrm{ln}\mathrm{\Lambda }\left[\mathrm{erf}(X)\frac{2X}{\sqrt{\pi }}e^{X^2}\right]\frac{𝐯_s}{v_s^3},$$
(4)
where $`X=|𝐯_s|/\sqrt{2}\sigma `$.

One important aspect to note is that both expressions of dynamical friction include here an anisotropic halo density, denoted by $`\rho _h[m^2(0)]=\rho _h[R^2+z^2/q_h^2]`$ in cylindrical coordinates. This is the “local approximation” made here but we note that the derivation of eq. 4 assumes a uniform, infinitely extended background medium. In practice, the local approximation implies the only difference between both expressions to be the anisotropic terms of the velocity distribution.

## 5 The semi-analytical code

In order to analyse the accuracy of the analytic expressions of dynamical friction we have constructed a semi-analytic algorithm that solves the equation of motion of a point-mass satellite

$`{\displaystyle \frac{d^2𝐱}{dt^2}}=𝐅_g+𝐅_{\mathrm{df}},`$ (5)
where $`𝐅_g`$ is the specific force from the parent galaxy and $`𝐅_{\mathrm{df}}`$ that due to dynamical friction (eqs. 3 or 4).

If the parent galaxy follows the fixed density profile of eq. (1), the specific force can be written in Cartesian coordinates (Chandrasekhar 1960) as

$`F_{g,i}=2\pi Gx_i{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{du}{(1+u)^2(1+e_h^2+u)^{1/2}}}\rho _h[m^2(u)],`$ (6)
$`F_{g,z}=2\pi Gz{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{du}{(1+u)(1+e_h^2+u)^{3/2}}}\rho _h[m^2(u)],`$
where $`x_i=x,y`$, the ellipticity of the galaxy is $`e_h^21q_h^2`$, and

$$m^2(u)=\frac{R^2}{1+u}+\frac{z^2}{1e_h^2+u}.$$
(7)
The algorithm employed to solve eq. (5) is based on the Bulirsch-Stoer method (for a complete description see Press et al. 1986), which provides high-accurate solutions with minimal computational effort. This method is based on an adaptive step-size scheme, thus, being ideal for systems with non-smooth potentials, as may be the case for satellites on highly eccentric orbits when disc and bulge components are included.

For calculating $`𝐅_{\mathrm{df}}`$ (eqs. 3 and 4) the Coulomb logarithm is treated as a free parameter to be fitted to the numerical orbits. The satellite mass $`M_s(t)`$ is obtained from the N-body data (see Section 7.1) and is treated as an input function. Numerical tests of orbits in a Keplerian potential show that energy and angular momentum are conserved to $`10^8`$ and $`10^9`$ per orbit, respectively, after choosing the accuracy of the Bulirsch-Stoer algorithm to be EPS$`=10^5`$. This value remains fixed for the calculations presented in this paper. An extensive description of our semi-analytic code can be found in Peñarrubia (2003).

## 6 Determining the Coulomb logarithm

The Coulomb logarithm is usually fit to numerical data with the aim of reproducing the overall orbital evolution by means of semi-analytic algorithms (e.g. Fellhauer et al. 2000). This procedure may actually be considered as a “calibration” of the semi-analytic code, which must be done carefully if a detailed inter-comparison between different schemes of dynamical friction is desired. For that reason we present in what follows a method to describe the accuracy of the semi-analytic scheme.

One possible way to quantify similarity of two orbits is via the quantity

$`\chi ^2={\displaystyle \frac{1}{2k}}{\displaystyle \underset{i=1}{\overset{2k}{}}}\left[(𝐫_i𝐫_{i,n})^2+\sigma ^2(r_0)(t_it_{i,n})^2\right],`$ (8)
$`𝐫`$ being the satellite position vector at the peri- and apo-galactica and $`t`$ the time at which the satellite passes by these points. The subindex $`n`$ denotes the numerical values and $`\sigma (r_0)`$ the velocity dispersion of the DMH at the initial galacto-centre distance. The sum is over a given number of orbits $`k`$.

The definition of $`\chi `$ measures the divergence of the numerical and semi-analytical satellite position. The term $`\sigma \mathrm{\Delta }t`$ allows for the possibility that the orbital periods differ during the evolution, which would lead to a secular deviation . By definition, $`\chi `$ is equivalent to the discrepancy between the numerical and semi-analytical position evolution per orbit. The selection of the maximum and minimum galacto-centre distances for comparison permits a direct control over the orbital eccentricity evolution, although the measure of $`\chi `$ can be extended to other points without loss of generality.

The value of $`k`$ depends on the objectives of the study. For instance, if the aim is to find the best calibration for a large period of time, as it may be to reproduce the satellite decay in spiral galaxies, the $`k`$-value must include a number of periods large enough, so that the overall evolution of several satellite orbits can be reproduced with a single value of $`\mathrm{ln}\mathrm{\Lambda }`$. In this paper, we limit our fit to the first satellite periods, namely, $`k=3,4`$, for which the differences between both approaches of dynamical friction can be clearly seen.

In Fig. 1 we plot $`\chi `$ for some of the experiments, concretely, those with inclinations $`60^{}`$, $`45^{}`$ and $`30^{}`$ (rows), with eccentricities 0.3, 0.5 and 0.7 (columns). For each model, the semi-analytic code is employed to generate the satellite orbit using Chandrasekhar’s (dashed and dotted-dashed lines) and Binney’s (full and dotted lines) formula to reproduce dynamical friction. This figure shows that Chandrasekhar’s formula poorly describes the dependence of the satellite orbit with the initial inclination, leading to a wider dispersion of the Coulomb logarithm values (for this range of inclinations, between 30 and 60, ln$`\mathrm{\Lambda }[0.9,2.8]`$). If Binney’s expression is used, the variation of $`\mathrm{ln}\mathrm{\Lambda }`$ is highly reduced (ln$`\mathrm{\Lambda }[2.3,2.5]`$), which shows that this scheme provides a much better description of the effects of anisotropic velocity dispersion on satellite decay, independently of the orbital inclination. The variation range of the Coulomb logarithm becomes larger for a wider inclination spread, but $`\chi `$ barely presents a dependence on the satellite eccentricity if Binney’s formula is applied.

If Chandrasekhar’s approach is used, the Coulomb logarithm that leads to the best fit becomes smaller as the inclination increases. Since dynamical friction is proportional to $`\mathrm{ln}\mathrm{\Lambda }`$, the use of the average value implies an overestimate of the force for low inclinations and vice versa. The final average,

$$(\chi /r_0)=\frac{1}{N_{\mathrm{ln}\mathrm{\Lambda }}}\underset{\mathrm{ln}\mathrm{\Lambda }}{}\chi _i/r_0,$$
(9)
over the $`N_{\mathrm{ln}\mathrm{\Lambda }}`$ numerical experiments of Table 2 is plotted in Fig. 2. This figure shows the large discrepancies produced by Chandrasekhar’s expression if the fit is for a large range of orbital inclinations and eccentricities, as expected. The minima of the curves determine the values of ln$`\mathrm{\Lambda }`$ that lead to the best fits, which we summarise in Table 3. The values of $`\chi _{\mathrm{min}}`$ denote the average error during the first $`k`$ orbits associated with the fit.

For the following analysis the value of the Coulomb logarithm implemented in our semi-analytic code will be found in Table 3. Looking at Fig. 1, we expect that Chandrasekhar’s friction will lead to more accurate fits for low ($`i30^{}`$) than for high inclined orbits after fixing ln$`\mathrm{\Lambda }=2.2`$.

## 7 The velocity anisotropy effects

In this Section, we discuss in more detail various aspects of satellite evolution.

### 7.1 Satellite mass loss

Tidal forces induce satellite mass loss. The satellite mass plays an important role in determining the ultimate fate of its evolution and survival (e.g. Velázquez & White 1999, Klessen & Kroupa 1998, Johnston et al. 1999). The value of $`M_s(t)`$ is treated as an input by the semi-analytical code, and is calculated using the self-consistent Superbox code.

The mass remaining bound to the satellite, $`M_s(t)`$, is determined in the Superbox calculations by computing the potential energy $`\mathrm{\Phi }_i<0`$ of each satellite particle presumed bound to the satellite, and its kinetic energy ($`T_i`$) in the satellite frame. Following PKB, particles with $`E_i=T_i+m_s(\mathrm{\Phi }_i+\mathrm{\Phi }_{\mathrm{ext}})>0`$ are labelled unbound, where $`m_s`$ is the mass of one satellite particle (all have the same mass) and the potentials

$`\mathrm{\Phi }_i={\displaystyle \underset{ij}{}}{\displaystyle \frac{Gm_s}{\sqrt{|r_ir_j|^2+ϵ^2}}},`$ (10)
$`\mathrm{\Phi }_{\mathrm{ext}}=|\mathrm{\Phi }_g(r_s)|,`$
the softening being $`ϵ0.23`$mu $`=0.8`$kpc, which is the resolution of the inner grid focused on the satellite centre-of-density $`r_s`$, and $`\mathrm{\Phi }_g`$ the galaxy potential at this point, neglecting the tidal terms. Thus, all the particles of a satellite are assumed to feel the same external potential, which is a useful and sufficiently accurate approximation, taking into account that most of the bound particles are located very close to this point.

Particles with $`E_i>0`$ are removed and the procedure repeated until only negative energy particles are left. As Johnston et al (1999) show, the energy criterium permits one to distinguish those particles that, though unbound, remain in orbits inside the tidal radius, which will escape from the satellite after some orbital periods.

The mass is calculated each $`\mathrm{\Delta }t=0.312`$ Gyr, so that the semi-analytic code interpolates the value for intermediate points at each time-step. The error is of the order of $`\mathrm{\Delta }M(t)/\mathrm{\Delta }t`$, going linearly with the mass loss. This means that the interpolation might introduce not negligible differences at times where the mass loss is significant (i.e late times of the satellite evolution). In this study we are not concerned in detail with the late phases of evolution.

In the right columns of Fig. 3a, b and c we plot the mass evolution for different orbits. Most of the satellites reach the inner most regions of the parent galaxy with a substantial fraction of their initial mass. A comparison of these curves with those of PKB (where a bulge and a disc component were included) shows that the baryonic subsystems of a galaxy induce a larger mass loss through tidal heating (e.g, Taylor & Babul 2001, Peñarrubia 2003). PKB observe in their numerical experiments that all satellites with $`M_s=0.1M_d`$ and $`r_0=55`$ kpc are destroyed before the remaining bound part of the satellite reaches the central region of the galaxy.

### 7.2 Satellite decay

One of the most important effects of dynamical friction is the monotonic reduction of the orbital angular momentum and energy during the satellite’s evolution that leads to a progressive decrease of the averaged galacto-centre distance. The computations carried out by PKB show a strong dependence of the decay time on the initial inclination that must be compared to analytic estimates.

In Fig. 3a, b and c we plot the radius evolution (left columns) for those models with initial $`e=0.5,0.7`$ and 0.3, respectively. From this figure, we conclude that Binney’s expression clearly produces more accurate results than Chandrasekhar’s for the whole range of orbital inclinations. This result is not surprising due to the small dependence of the Coulomb logarithm on the inclination and eccentricity as shown in Fig. 1. Additionally, the value of $`\mathrm{ln}\mathrm{\Lambda }`$ that produces the best fit for the first few orbits also succeeds in reproducing the decay time of the satellite.

PKB observe that coplanar satellites suffer larger friction than those following polar orbits, leading to survival times 70% shorter. Due to the presence of a disc in their galaxy model, the contribution of the disc anisotropy on the decay differentiation as a function of the inclination cannot be directly measured. The calculations presented here (where the disc and bulge are removed) show survival times that range from 3.7 Gyr (coplanar orbits) to 6 Gyr (polar orbits), using the same orbital parameters and halo flattening as PKB. This implies a decay time difference of around 60% between polar and coplanar satellites, which indicates that the disc contribution might be of the order of 10%. The effects of the disc on the satellite orbit can be found in Peñarrubia (2003) and will be addressed in a following paper.

Depending on the symmetry of the halo distribution, one can observe the following effects:

* Spherical mass distribution and isotropic velocity distribution: Satellites orbiting systems with a spherical distribution function move on orbits that do not depend on their orientation with respect to the symmetry axis.
* Flattened mass distribution and isotropic velocity distribution: By means of the local approximation the value of dynamical friction is determined by the properties of the halo at the satellite’s position, being reproduced by Chandrasekhar’s formula. Our results indicate that the spatial asphericity leads to a strong differentiation of the satellite decay as a function of the orbital inclination.
* Flattened mass distribution and anisotropic velocity distribution: The main influence of the velocity anisotropy on the satellite orbit is the secular reduction of the orbital inclination (see Subsection 7.3) and the reduction of the spatial anisotropy effects, which is equivalent to $`B_R<B_z`$ in Binney’s formula (eq. 3, oblate systems). Taking into account the velocity anisotropy explicitly thus reduces the difference in orbital decay times between satellites on polar and coplanar orbits. As Fig. 3 shows, assuming an isotropic distribution in velocity space ($`B_R=B_z`$) or, equivalently, using Chandrasekhar’s formula to reproduce dynamical friction, leads to an overestimation of dynamical friction for low inclined satellites and to an underestimation for those satellites following highly inclined orbits. This is due to the fact that the anisotropy leads to a reduced effective $`X`$ for high inclined orbits yielding an enhanced friction force and vice versa for low inclination. This competes with the density effect due to the flattening, because highly inclined orbits have lower mean density along the orbit compared to lower inclination.

The evolution of the orbital radius of satellites with initial $`e=0.3,0.7`$ is plotted in Fig. 3b and c. As we can observe, Binney’s approximation reproduces accurately the overall radius evolution independent of the initial eccentricity and orbital inclination.

### 7.3 Evolution of the orbital inclination and eccentricity

Orbits in non-spherical systems have inclinations $`(i\mathrm{arccos}[L_z/L])`$ that do not remain constant but suffer periodical oscillations due to nutation. In addition, the orbital plane precesses at a constant rate. Once the initial conditions are fixed, the amplitude and frequency of nutation and the precession rate remain constant if the friction force is removed from the equations of motion, whereas, if implemented, nutation and precession vary according to the angular momentum and radial distance evolution. Our interest focuses now on the effects induced by the velocity anisotropy on the satellite inclination during the orbit.

Binney (B77) predicts a progressive reduction of $`i`$ due to dynamical friction if the velocity dispersion ellipsoid is axi-symmetric $`(\sigma _R,\sigma _z)`$ and $`\sigma _R>\sigma _z`$. By symmetry, the inclination decrease will not occur if the orbits are either coplanar $`(i=0^{})`$ or polar $`(i=90^{})`$.

The inclination evolution of models with $`e=0.5`$ is plotted in Fig. 6 (left column), where dotted lines denote the numerical data and solid and dashed lines the semi-analytic evolution if dynamical friction is modelled by Binney’s and Chandrasekhar’s formulæ, respectively. This Figure shows the reduction of the mean value of $`i`$ predicted by Binney and observed by PKB in their numerical experiments. After the satellite has sunk to the inner most region of the halo, the inclinations are as low as $`10^{}20^{}`$ independent of their initial value. This large decrease of $`i`$ is well reproduced by Binney’s expression, although the nutation process shows discrepancies with the numerical result, which is connected with the not exact reproduction of the orbit. Despite the accurate description of the overall decay process, this orbital mismatch is also observed when applying Chandrasekhar’s expression in spherical systems (see Just & Peñarrubia 2002). The figure also confirms that the orbital inclination of coplanar and polar satellites remains constant.

If dynamical friction is modelled by Chandrasekhar’s formula, i.e the velocity distribution is assumed to be isotropic, the averaged value of $`i`$ does not change during the orbit, which contradicts the numerical results. Thus, while the difference between the survival times of polar and coplanar orbits is larger for flattened DMHs treated with Chandrasekhar’s dynamical friction formula (Fig. 3), taking into account the DMH anisotropic velocity distribution function explicitly via Binney’s formulae leads to an increase with time of the anisotropy of the satellite distribution due to the kinematical coupling of the satellites to the DMH velocity field. This clearly agrees with the fully self-consistent $`N`$body computations reported here.

In Fig. 6b and 6c (left columns) we plot the comparison for models with $`e=0.7,0.3`$, respectively. The results show barely a dependence on the eccentricity. It is interesting to note that, independently of $`e`$, orbits that are neither coplanar nor polar present large drops of the mean value of $`i`$. After the satellite sinks to the centre, the final orbital inclination lies in between 10-20 for all the models.

We must emphasise the accuracy of Binney’s formula in describing correctly the inclination decrease that satellites suffer in oblate systems. This is crucial for simulating properly satellite motions and for investigating satellite distributions around spiral galaxies.

Like the orbital inclination, the eccentricity is one of the orbital parameters that can be indirectly measured from observations to determine a satellite’s motion around a galaxy. In the right columns of Fig. 6 we show the comparison of the numerical eccentricity evolution with both semi-analytic approaches. The analytic formulæ of dynamical friction lead to a larger eccentricity decrease, which occurs mostly at late-times of the orbit, the so-called orbital circularisation, and becomes stronger for low inclined orbits, those that suffer higher dynamical friction. Fig. 6b and Fig. 6c indicate that the circularisation is more pronounced if the initial orbital eccentricity is higher and decreases for more circular orbits. Both dynamical friction expressions lead to a similar eccentricity evolution for the first few orbital periods, however at late-times the eccentricity exhibits a reduction not present in the numerical calculations that can be as high as $`80\%`$ for low inclined satellites on highly eccentric orbits (Fig. 6b, model H2S100e).

### 7.4 Energy and angular momentum evolution

A flattened system possesses two analytic constants of motion, the energy and the component of the angular momentum perpendicular to the axi-symmetry plane (which we denote as $`L_z`$). The total angular momentum $`L^2=L_R^2+L_z^2`$ is, however, not constant during a satellite orbit (see e.g Binney & Tremaine 1987), but has periodic variations that correspond to a precession and nutation of the orbital plane around the $`z`$-axis.

Since the dynamical friction force has an opposite sense with respect to the satellite velocity, it decreases the angular momentum and the energy which induces a monotonic sinking into the inner regions of the halo potential. The reduction of angular momentum, therefore, implies an increase of the binding energy (in absolute value), since the potential grows with decreasing radius. Due to the small magnitude of dynamical friction compared to the mean field force, we expect an easier comparison between numerical and semi-analytic data by the slow variation of $`L_z`$ and $`E`$ during the orbit.

In Fig. 9 we plot the changes of specific $`E`$ and $`L_z`$ due to dynamical friction for the models with $`e=0.5`$. The results are equivalent to those of the radial evolution. Due to our selection of $`\mathrm{ln}\mathrm{\Lambda }`$ as the average over those Coulomb logarithms that lead to the best fit for each particular model, Chandrasekhar’s formula overestimates dynamical friction for low inclined orbits and underestimates it for highly inclined orbits. For orbits with $`i<30^{}`$, this appears as a stronger reduction of the $`z`$-component of angular momentum and, equivalently, a large increase of the binding energy. The effect is contrary for satellites with $`i>30^{}`$.

This figure illustrates how the kinetic energy of the satellite is lost via friction, being taken-up partially by halo particles and also being deposited in unbound satellite particles. At the end of the simulation the angular momentum has a null value, i.e the satellite remains in the inner most part of the galaxy.

It is interesting to note that the numerical evolution of energy and angular momentum presents small oscillations during the orbit. This behaviour is due to the self-response of the halo to the satellite motion. Since superbox preserves the total energy and angular momentum, the halo also moves around the centre-of-mass of the system. Due to the complexity of the feedback, it cannot be reproduced analytically (the halo centre-of-mass is fixed at the coordinate origin in the semi-analytic code).

## 8 Conclusions

To asses the accuracy of Binney’s equations (B77) and in order to reproduce the decay of satellites in flattened DMHs with a semi-analytical code, we have performed a set of self-consistent numerical experiments for different orbital inclinations and eccentricities of the satellite. The semi-analytic code incorporates dynamical friction either in terms of Chandrasekhar’s expression or as Binney’s formulae. Both treatments include the aspherical density profile by means of the local approximation. This means that the differences on the satellite motion induced by each treatment of dynamical friction comes from the anisotropy in velocity space, which is implemented in the analysis of B77.

The accuracy of Binney’s and Chandrasekhar’s formulæ in comparison with the numerical orbits is determined by the parameter $`\chi ^2=\frac{1}{2k}(\mathrm{\Delta }𝐫^2+\sigma ^2\mathrm{\Delta }t^2)`$ at the peri- and apo-centres for a given number $`k`$ of orbits. If dynamical friction is modelled by Binney’s equation, this quantity shows discrepancies of approximately $`\chi _{\mathrm{min}}=0.009r_0`$ per orbit after averaging over the set of experiments and for the first three orbits, while Chandrasekhar’s formula produces values of around $`\chi =0.05r_0`$ per orbit (Table 3).

We conclude that Binney’s expression faithfully reproduces the process of dynamical friction in anisotropic systems. The fit is as accurate as that employing Chandrasekhar’s formula in isotropic systems (see Just & Peñarrubia 2002).

The comparison of orbits resulting from Chandrasekhar’s and Binney’s expression of dynamical friction gives us the possibility to asses the effects of the DMH velocity anisotropy on satellite dynamics. We have demonstrated that, (i) if the density profile is in both equations $`\rho =\rho (R,z)`$, where $`R,z`$ are the cylindrical coordinates of the satellite position vector, the orbits generated by Chandrasekhar’s formula overestimate the decay time for polar orbits and underestimate it for coplanar ones for an overall best-fit Coulomb logarithm (ln$`\mathrm{\Lambda }=2.2`$). One effect of the velocity anisotropy is then to reduce the interval of decay times as a function of the orbital inclination. Binney’s expression has been shown to reproduce accurately the numerical results independently of the initial eccentricity and with ln$`\mathrm{\Lambda }=2.4`$. (ii) Dynamical friction in systems with anisotropic velocity distribution leads to a marked decrease of the orbital inclination $`i`$ which is well reproduced by Binney’s expression. After the satellite sinks to the inner most region of the galaxy, $`i`$ lies within 10-20, independently of the initial value unless $`i90^\mathrm{o}`$ (polar orbits). (iii) The energy and angular momentum evolution as a function of the orbital inclination confirm the results of (i) and (ii).

The semi-analytic eccentricity evolution, either employing Chandrasekhar’s or Binney’s formula, shows the so-called circularisation process, defined as the progressive reduction of $`e`$ during the orbit. This variation is stronger for increasing friction (as for coplanar orbits or during the late-times of the evolution) and barely takes place in the self-consistent $`N`$body calculations. The small $`N`$body circularisation agrees with the results of van den Bosch et al. (1999). A possible solution may be sought in the position-dependence of the Coulomb logarithm, as proposed by Hashimoto, Funato & Makino (2003). However, despite the improvement in the description of the orbit at early-times, the scheme of Hashimoto et al. overestimates the satellite decay time for all the experiments (see Just & Peñarrubia 2002). We must conclude that the reason for circularisation in the semi-analytic orbits is not yet fully understood.

## 9 Acknowledgements

JP acknowledges support through a SFB 439 grant at the University of Heidelberg. PK acknowledges support through DFG grant KRI635/4. We thank the anonymous referee for his/her advice and useful comments.


text:
# Kerr black holes in horizon-generating form

(28th January 2004)

## Abstract

New coordinates are given which describe non-degenerate Kerr black holes in dual-null foliations based on the outer (or inner) horizons, generalizing the Kruskal form for Schwarzschild black holes. The construction involves an area radius for the transverse surfaces and a generalization of the Regge-Wheeler radial function, both functions of the original radial coordinate only.

It is hard to overstate the importance of the Kerr solution either in the theoretical history of black holes or as a model for actual astrophysical black holes. The solution was discovered Ker  before black holes were understood as such, but it soon became a key ingredient in the classical theory of black holes BCH ; HE ; MTW ; Wal . Uniqueness theorems indicate that this small family of solutions, parametrized by just mass and angular momentum, constitute all stationary black holes in vacuum, and therefore serve as good models for most real black holes for most of the time, from supernova remnants to supermassive black holes mysteriously found in galactic cores. Dynamical black holes have a partial theory being currently developed bhd ; 1st ; in ; cyl ; ql1 ; gwbh ; mg9 ; AK1 ; AK2 ; BF , for which an appropriate understanding of Kerr black holes will also be necessary.

The generic Kerr metric is much more complex than the Schwarzschild metric, to which it reduces in the non-rotating case, and indeed Kerr discovered it almost a half-century after Einstein formulated General Relativity Ein , as compared to Schwarzschild’s few weeks Sch . While it has deservedly been one of the most studied exact solutions, there appears to be a peculiar lack of natural coordinates generating the outer (or inner) horizons, as in the Kruskal extension Kru  of the Schwarzschild solution, which was so crucial for understanding its nature as a black hole. (There is an earlier paper of Fronsdal Fro  which acknowledges priority to unpublished work of Kruskal; the Kruskal paper was actually written by Wheeler Whe ). Specifically, recall that a non-degenerate Kerr black hole has two asymptotically flat external regions and two outer Killing horizons, intersecting in a bifurcation surface or wormhole throat BL ; Car . Kruskal’s coordinates (or rather their sum and difference) describe a dual-null foliation of the space-time, meaning two families of null hypersurfaces intersecting in a two-parameter family of transverse spatial surfaces, such that the horizons are two of the hypersurfaces. Such coordinates are naturally adapted to the horizons and radiation propagation, both into the black hole and out of its vicinity. This Letter reports generalizations to non-degenerate Kerr black holes.

The Kerr solution with mass $`m`$ and angular momentum $`ma`$ in Boyer-Lindquist BL  coordinates $`(t,r,\theta ,\varphi )`$ can be given by the inverse metric

$$g^1=\frac{1}{\mathrm{\Sigma }}\left(\mathrm{\Delta }_r^2+_\theta ^2+\frac{\mathrm{\Delta }a^2\mathrm{sin}^2\theta }{\mathrm{\Delta }\mathrm{sin}^2\theta }_\varphi ^2\frac{4mar}{\mathrm{\Delta }}_\varphi _t\frac{\mathrm{\Pi }}{\mathrm{\Delta }}_t^2\right)$$
(1)
where $`_x`$ are the coordinate derivative vectors, $`_x(dx^{})=\delta _{xx^{}}`$, symmetric tensor products are implied and

$`\mathrm{\Delta }`$ $`=`$ $`r^2+a^22mr`$ (2)
$`\mathrm{\Sigma }`$ $`=`$ $`r^2+a^2\mathrm{cos}^2\theta `$ (3)
$`\mathrm{\Pi }`$ $`=`$ $`(r^2+a^2)^2\mathrm{\Delta }a^2\mathrm{sin}^2\theta .`$ (4)
The inverse metric is crucial in this adventure, since one is seeking null coordinates $`x^\pm `$: $`g^1(dx^\pm ,dx^\pm )=0`$. For $`a<m`$, the Killing horizons are located at $`r=r_\pm `$, where $`r_\pm =m\pm \sqrt{m^2a^2}`$ are the roots of $`\mathrm{\Delta }(r)=0`$. Here the outer horizons $`r=r_+`$ will be studied in detail, with modifications for the inner horizons $`r=r_{}`$ mentioned below. The desired transformation may be given in two steps, firstly to coordinates $`(t^{},r^{},\vartheta ,\phi )`$ defined by

$`t^{}`$ $`=`$ $`ta\mathrm{sin}\theta `$ (5)
$`r^{}`$ $`=`$ $`{\displaystyle \frac{R^2}{\mathrm{\Delta }}𝑑r}`$ (6)
$`\vartheta `$ $`=`$ $`\theta `$ (7)
$`\phi `$ $`=`$ $`\varphi \mathrm{\Omega }(ta\mathrm{sin}\theta )`$ (8)
where $`\mathrm{\Omega }=a/2mr_+`$ is the angular velocity of the horizon and

$$R(r)=\left((r^2+a^2)^2a^2\mathrm{\Delta }\right)^{1/4}=\left(r^4+a^2r^2+2ma^2r\right)^{1/4}.$$
(9)
Secondly one may transform to $`(x^+,x^{},\vartheta ,\phi )`$, where

$$x^\pm =\pm e^{\kappa (r^{}\pm t^{})}$$
(10)
and $`\kappa `$ is determined below to be the surface gravity. The last step constructs dual-null coordinates $`x^\pm `$ by a generalization of Kruskal’s method, having first identified appropriate coordinates $`(t^{},r^{})`$. The colatitude $`\vartheta `$ has been renamed merely for clarity in partial derivatives, and the new longitude $`\phi `$ can be partly understood by noting that it is normal to the horizon-generating Killing vector $`\chi =_t+\mathrm{\Omega }_\varphi `$: $`\chi (d\phi )=\chi (d\vartheta )=0`$. The same is true for $`\varphi \mathrm{\Omega }t`$ alone Cha , but it ceases to be spatial at large radius. The key step is the identification of the new coordinate $`t^{}`$, which absorbs mischievous $`\theta `$-dependent terms to leave the corresponding $`r^{}`$ (such that $`r^{}\pm t^{}`$ are null) miraculously available as a function of $`r`$ alone. One may recognize $`r^{}`$ as a generalization of the Regge-Wheeler RW  tortoise coordinate, named after the ancient paradox of Achilles and the tortoise, as $`r=r_+`$ corresponds to $`r^{}\mathrm{}`$.

The coordinate derivative vectors are related by

$`_t`$ $`=`$ $`\kappa (x^+_+x^{}_{})\mathrm{\Omega }_\phi `$ (11)
$`_r`$ $`=`$ $`{\displaystyle \frac{R^2}{\mathrm{\Delta }}}\kappa (x^+_++x^{}_{})`$ (12)
$`_\theta `$ $`=`$ $`_\vartheta +\mathrm{\Omega }a\mathrm{cos}\vartheta _\phi \kappa a\mathrm{cos}\vartheta (x^+_+x^{}_{})`$ (13)
$`_\varphi `$ $`=`$ $`_\phi `$ (14)
where $`_\pm =_{x^\pm }`$. Here $`r`$ is retained as a function of $`x^+x^{}=e^{2\kappa r^{}}`$, determined implicitly by inverting (6) to give $`r(r^{})`$. By carefully arranged cancellations in $`R^4+\mathrm{\Delta }a^2\mathrm{cos}^2\theta \mathrm{\Pi }=0`$, the inverse metric transforms to

$$g^1=\frac{1}{\mathrm{\Sigma }}\left(\frac{4R^4\kappa ^2x^+x^{}}{\mathrm{\Delta }}_+_{}+2\kappa (x^+_+x^{}_{})(\alpha _\phi a\mathrm{cos}\vartheta _\vartheta )+_\vartheta ^2+2\mathrm{\Omega }a\mathrm{cos}\vartheta _\vartheta _\phi +\left(\frac{1}{\mathrm{sin}^2\vartheta }\beta \right)_\phi ^2\right)$$
(15)
where

$`\alpha (r)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }R^42mar}{\mathrm{\Delta }}}`$ (16)
$`\beta (r)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }^2R^44\mathrm{\Omega }mar+a^2}{\mathrm{\Delta }}}.`$ (17)
Note that $`x^\pm `$ are null by the absence of $`_\pm ^2`$ components, and that $`r=r_+`$ has been mapped to $`x^+x^{}=0`$. Then $`x^\pm =0`$ are the two branches of the outer horizons, with $`x^+=x^{}=0`$ being the bifurcation surface. The Killing vector $`\chi =\kappa (x^+_+x^{}_{})`$ has a saddle point at each angle of the bifurcation surface: it vanishes there, and generates the two branches of the horizons along the principal directions. Having constructed a dual-null foliation based on the horizons in particular coordinates, the remaining coordinate freedom consists of relabelling the null hypersurfaces, $`x^\pm \widehat{x}^\pm (x^\pm )`$, and diffeomorphisms of a transverse surface, i.e. choice of angular coordinates.

The functions $`(\alpha ,\beta ,x^+x^{}/\mathrm{\Delta })`$, each apparently singular at $`\mathrm{\Delta }=0`$, are all actually regular at $`r=r_+`$. Firstly, since $`\mathrm{\Delta }=(rr_+)(rr_{})`$, near $`r=r_+`$ one has $`r^{}(R^2/(rr_{}))|_{r=r_+}\mathrm{ln}(rr_+)`$, so $`e^{2\kappa r^{}}rr_+`$ if

$$\kappa =\frac{rr_{}}{2R^2}|_{r=r_+}=\frac{\sqrt{m^2a^2}}{2mr_+}$$
(18)
which is the usual surface gravity. Then $`x^+x^{}/\mathrm{\Delta }1/(rr_{})`$ is regular at $`r=r_+`$. Secondly one can rearrange

$`\alpha `$ $`=`$ $`\mathrm{\Omega }\left({\displaystyle \frac{(r+r_+)(r^2+a^2)}{rr_{}}}+r_+^2\right)`$ (19)
$`\beta `$ $`=`$ $`\mathrm{\Omega }^2\left({\displaystyle \frac{(r+r_+)^2(rr_+)}{rr_{}}}+2r_+^2+a^2\right)`$ (20)
which are now manifestly regular at $`r=r_+`$. Then $`g^1`$ is regular, except for the polar coordinate breakdown at $`\mathrm{sin}\vartheta =0`$. With this usual caveat, the coordinate chart covers the entire region outside the inner horizons.

For the inner horizons, one can use the same formulae (5)–(17) if one replaces the angular velocity and surface gravity by $`\mathrm{\Omega }=a/2mr_{}`$ and $`\kappa =\sqrt{m^2a^2}/2mr_{}`$ respectively. Thus it is convenient to use a negative surface gravity for the inner horizons, as is anyway natural for dynamical horizons 1st ; in ; cyl ; ql1 ; gwbh ; mg9 . Then the new coordinate chart covers the region inside the outer horizons, up to the discs $`r=0`$, where $`R=0`$. Apart from the ring singularities at $`\mathrm{\Sigma }=0`$, the discs are regular and one can extend through to the internal asymptotically flat regions, which can also be covered by a similar chart. The ring singularities present topological obstructions to any encircling foliation by spheres.

The metric

$$g=\frac{\mathrm{\Sigma }}{\mathrm{\Delta }}dr^2+\mathrm{\Sigma }d\theta ^2+\frac{1}{\mathrm{\Sigma }}\left(\mathrm{\Pi }\mathrm{sin}^2\theta d\varphi ^24mar\mathrm{sin}^2\theta d\varphi dt(\mathrm{\Delta }a^2\mathrm{sin}^2\theta )dt^2\right)$$
(21)
in the new coordinates may be found by inverting the coordinate transformation to give the coordinate differentials

$`dt`$ $`=`$ $`dt^{}+a\mathrm{cos}\vartheta d\vartheta `$ (22)
$`dr`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{R^2}}dr^{}`$ (23)
$`d\theta `$ $`=`$ $`d\vartheta `$ (24)
$`d\varphi `$ $`=`$ $`d\phi +\mathrm{\Omega }dt^{}.`$ (25)
An independent and more systematic method uses the general dual-null decomposition dne

$$g^1=2e^f_+_{}+2e^f(s_+_{}+s_{}_+)2e^fs_+s_{}+h^1$$
(26)
where $`h^1`$ is the inverse metric of the (constant-$`x^\pm `$) transverse surfaces, $`s_\pm `$ are shift vectors and $`f`$ is a normalization function. In this case one finds

$`e^f`$ $`=`$ $`{\displaystyle \frac{2R^4\kappa ^2x^+x^{}}{\mathrm{\Delta }\mathrm{\Sigma }}}`$ (27)
$`s_\pm `$ $`=`$ $`\pm {\displaystyle \frac{\mathrm{\Delta }}{2R^4\kappa x^\pm }}(\alpha _\phi a\mathrm{cos}\vartheta _\vartheta )`$ (28)
$`h^1`$ $`=`$ $`{\displaystyle \frac{1}{R^4\mathrm{\Sigma }}}\left(\mathrm{\Pi }_\vartheta ^2+4ma^2r\mathrm{cos}\vartheta _\vartheta _\phi +\left({\displaystyle \frac{R^4}{\mathrm{sin}^2\vartheta }}a^2r(r+2m)\right)_\phi ^2\right).`$ (29)
The determinant of $`h^1`$ is just $`1/R^4\mathrm{sin}^2\vartheta `$, so it may be easily inverted to yield the transverse metric

$$h=\frac{1}{\mathrm{\Sigma }}\left((r^4+r(r+2m)a^2\mathrm{cos}^2\vartheta )d\vartheta ^24ma^2r\mathrm{cos}\vartheta \mathrm{sin}^2\vartheta d\vartheta d\phi +\mathrm{\Pi }\mathrm{sin}^2\vartheta d\phi ^2\right)$$
(30)
where symmetric tensor products are again implied. Then the metric is given by the general form dne

$$g=h+2h(s^+)dx^++2h(s^{})dx^{}+h(s_+,s_+)(dx^+)^2+h(s_{},s_{})(dx^{})^2+2(h(s_+,s_{})e^f)dx^+dx^{}$$
(31)
by calculating

$`h(s_\pm )`$ $`=`$ $`\pm {\displaystyle \frac{\mathrm{\Delta }}{2\mathrm{\Sigma }\kappa x^\pm }}\left((\alpha +\mathrm{\Omega }a^2\mathrm{cos}^2\vartheta )\mathrm{sin}^2\vartheta d\phi (1+\gamma \mathrm{sin}^2\vartheta )a\mathrm{cos}\vartheta d\vartheta \right)`$ (32)
$`h(s_\pm ,s_\pm )`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Sigma }(2\kappa x^\pm )^2}}\left((\beta +\mathrm{\Omega }^2a^2\mathrm{cos}^2\vartheta )\mathrm{sin}^2\vartheta 1+\mathrm{\Sigma }^2/R^4\right)`$ (33)
$`h(s_+,s_{})`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{4\mathrm{\Sigma }\kappa ^2x^+x^{}}}\left((\beta +\mathrm{\Omega }^2a^2\mathrm{cos}^2\vartheta )\mathrm{sin}^2\vartheta 1+\mathrm{\Sigma }^2/R^4\right)`$ (34)
where

$$\gamma (r)=\mathrm{\Omega }\alpha \beta =\frac{2\mathrm{\Omega }mara^2}{\mathrm{\Delta }}=\frac{2\mathrm{\Omega }ma}{rr_{}}.$$
(35)
The terms involving the shift vectors are quite complex as compared to those in the inverse metric. However, the transverse metric $`h`$ has an unexpectedly nice property: its area form is just $`R^2\mathrm{sin}\vartheta d\vartheta d\phi `$. (All the usual coordinates Ker ; BL ; Car ; Cha  have a transverse area form with more complex $`\theta `$-dependence). Thus the new radial function $`R`$ is revealed as the area radius of the transverse surfaces. This property is useful for evaluating surface integrals, most simply the area $`4\pi R^2`$ of the transverse surfaces. Also, while the coordinates were intended only to cover the outer horizons, the angular coordinates $`(\vartheta ,\phi )`$ become standard spherical polar coordinates on round spheres as $`R\mathrm{}`$: $`h/R^2d\vartheta ^2+\mathrm{sin}^2\vartheta d\phi ^2`$.

One might at first think that this dual-null foliation is unique, since the foliations of each horizon are generated by $`\chi `$ from the bifurcation surface. However, $`\chi `$ vanishes at the bifurcation surface, leaving the possibility to include it in different foliations of the horizons. In fact, another such foliation can be found simply by replacing $`a`$ with $`a`$ in the definitions of $`t^{}`$ (5) and $`\phi `$ (8), which has the effect of changing the sign of $`d\vartheta `$ in the metric, or of $`_\vartheta `$ in the inverse metric. The two sets of transverse surfaces have the same poles, but their equators are displaced relative to each other. The displacements are in opposite directions in the two external regions, so the two foliations are dual under left-right inversion of the space-time. The existence of such dual foliations is due to breaking the explicit symmetry $`(t,\varphi )(t,\varphi )`$ in $`(t^{},\phi )`$. One may choose whether to displace the time coordinate to the future or past in one universe, forcing the opposite choice to be made in the other universe. Despite this peculiarity, it seems likely that this pair of foliations is unique in some sense. For instance, one might conjecture that they are the only dual-null foliations based on the horizons for which the transverse area form is a function of $`r`$ times the area form of a unit sphere.

The new coordinates are so naturally adapted to Kerr horizons and radiation propagation that they are likely to have various applications, while the area radius $`R`$ and the generalized tortoise coordinate $`r^{}`$ may themselves prove useful. Applications might include electromagnetic radiation powered by rotating black holes, gravitational radiation produced by perturbed black holes, quasi-normal modes, Cauchy-horizon instability and black-hole evaporation. It might also help in the search for practical definitions of energy-momentum and angular momentum in the absence of symmetries, instead using extrinsic properties of embedded surfaces. The author’s main interest in this enterprise was to find the Kerr solution in a similar dual-null form to that which has proved useful in understanding dynamical black holes. Its role as a basic test case may now be employed to further develop a general theory of black-hole dynamics.

Supported by research grant “Black holes and gravitational waves” of Ewha Womans University.


text:
# Investigating the Geometrical Structure of Disordered Sphere Packings

## 1 Introduction

For several centuries, the structures generated by tightly packed mono-sized spheres have been studied by scientists as model systems for understanding the emergence of order and crystallisation or -vice versa\- the appearance of disordered and amorphous phases in natural systems . Despite all the efforts, several questions remain unsolved. One of the most intriguing and challenging problems is to understand whether or not the notion of ‘ideal’ or ‘typical’ or ‘common’ disordered packing configurations has any empirical basis . Indeed, disorder does not exclude organization and the positions of the packed spheres are locally highly correlated. Experimentally, it is easy to observe that when spheres with approximately equal sizes are poured into a jar, squeezed into a bag, shaken or mixed, they form structures with densities in a narrow range. No stable packing configurations have been observed with densities below 0.56 and no disordered assemblies of spheres have been observed with densities above 0.645 . Molecular dynamics simulations show that hard sphere systems behave like a gas below $`\rho 0.49`$. Upon compaction, liquid-like behaviour is observed up to $`\rho 0.55`$, where crystallisation starts to occur. If crystallization is avoided, the system undertakes a ‘glass’ transition at $`\rho 0.56`$ and then it can be compacted up to $`\rho 0.645`$ where no further densification can be induced . Empirical and simulated evidence suggest that something very special might happen in the geometry of the packing at densities above $`\rho 0.56`$; a process which must terminate below $`\rho 0.65`$. What makes the understanding of this process particularly challenging is that there are no a priori reasons for the densification process to stop around $`\rho 0.64`$. On the contrary, there are plenty of local configurations which are denser than this limit and any arrangement of spheres in stacked planar hexagonal closed packed layers (the so called Barlow packings ) can reach the density of $`\rho =\pi /\sqrt{18}0.74`$, as achieved in fcc (face-centered cubic) or hcp (hexagonal closed-packed) crystalline packings.

## 2 Experimental Apparatus and Methodology

In this paper we report the analysis on 6 bead pack samples, all in a cylindrical container with an inner diameter of $`D55mm`$ and filled to a height of $`H75mm`$. Specifically:

* Two large samples with 150,000 beads and diameters $`d=1.00\pm 0.05mm`$;
* Four smaller samples with 35,000 beads and diameters $`d=1.59\pm 0.05mm`$.

A packing realization is shown in Fig.1. To our knowledge, the present study is the first empirical investigation of *large* packings of *monosized* spheres. Indeed, there are only two previous investigations which use a similar approach , but the first concerns a single test-sample with $`2,000`$ spheres only, whereas the second investigates polydisperse assemblies with $`15,000`$ spheres.

Sample Preparation: Low density packings are obtained by placing a stick in the middle of the container and pouring the beads into it. Then, the stick is removed leaving a packing density $`\rho 0.58`$ . Higher densities, up to $`\rho 0.63`$, are achieved by gently tapping the sample. The densest sample at $`\rho =0.64`$ is obtained by a combined action of gently tapping and compression from above (with the upper surface left unconfined at the end of the preparation). To reduce boundary effects, the inside of the cylinder have been made rough by randomly gluing spheres on the internal surfaces.

XCT Imaging: A X-ray Computed Tomography apparatus (see Sakellariou et al. in this Issue ) is used to image the samples. The two large samples were analysed by acquiring data sets of $`2000^3`$ voxels with a spatial resolution of $`0.05mm`$; whereas the four small samples were analysed by acquiring data sets of $`1000^3`$ voxels with a spatial resolution of $`60\mu m`$. After segmentation (see Sheppard et al in this issue ) the sample data sets are reduced to three-dimensional binary images, representing two distinct phases, one associated with the spheres and the other with air space.

Sphere Centres: In order to proceed with the analysis of the geometrical and statistical properties the position of all sphere centres are calculated from the binary images. Our approach is to find the sphere centres by moving a reference sphere ($`S`$) throughout the binarised sphere pack ($`P`$) and measuring the local overlap between $`S`$ and $`P`$. This corresponds to a 3-dimensional convolution: $`PS`$. This method is made highly efficient by applying the convolution theorem which allow to transform the convolution into a product in Fourier space: $`[PS]=[P][S]`$, where $``$ represents the (fast)Fourier Transform. The algorithm proceeds in 4 steps: 1) fast Fourier transform of the binary image ($`[P]`$); 2) transform the digitised map of the reference sphere ($`[S]`$, chosen with a diameter about 10 % smaller than $`d`$); 3) perform the direct product between these two; 4) inverse-transform of the product: $`^1[[P][S]]=PS`$. The result is an intensity map of the overlapping between the reference sphere and the bead pack, where the voxels closer to the sphere centres have a higher intensity. A threshold on the intensity map, locates the groups of voxels surrounding the sphere centres. The centre of mass of these grouped voxels is a very good estimation of the sphere centres in the pack.

Central Region: All the analyses reported hereafter have been performed over a central region ($`𝐆`$) at 4 sphere-diameters away from the sample boundaries. Note that spheres outside $`𝐆`$ are considered when computing the neighbouring environment of spheres in $`𝐆`$. The two large samples have about $`N_G80,000`$ spheres in $`𝐆`$, whereas the four smaller have about $`N_G20,000`$. In Table 1 the number of spheres in this region ($`N_G`$) is reported for each sample.

## 3 Study of the Packing Configurations and Discussion

Densities: We calculate the *local-densities* and the *sample-densities*. The local-densities are the fractions between the sphere-volumes and the volumes of the Voronoï cells constructed around the centre of each sphere in the sample (the Voronoï cell is the portion of space closest to a given centre in respect of any other centre). The sample-densities are fractions between the sum over the volumes of the spheres in $`𝐆`$ and the sum over the volumes of the Voronoï cells associated with these spheres. We observe bell-shaped distributions with average densities in the range $`0.586\rho 0.640`$ and standard deviations $`\sigma `$ within 1.5 % (see Table 1). However, the local densities are not homogeneously distributed in the sample. Typically, the densities are relatively smaller than the average in a region close to the cylinders central axis; the density increases going out from the centre then it saturates to rather homogeneous values up to a distance of a few (2-3) sphere diameters from the boundary. Rather inhomogeneous densities are also observed in the vertical direction, but in this case we find different behaviours depending on the sample-preparation. A detailed analysis of these behaviours will be the subject of a future publication.

Local Environment: The identification of touching spheres is, in general, an ill-defined problem from a experimental point of view. Indeed, the result is affected by the precision of the calculated sphere centres and their diameters. At the present we are developing a new technique to overcome this problem. With the present data, we assume that the spheres in contact are located at a radial distance between the diameter $`d`$ and $`d+v`$, where $`v`$ is the voxel-size. Table 1 reports the values of the average number of contacting neighbours ($`N_c`$) computed in $`𝐆`$ at the three different radial distances: $`d`$, $`d+v/2`$ and $`d+v`$. We observe values between $`N_c`$ $`3`$ and $`8`$, and an increasing trend with the packing density.

Packing Efficiency: We compute the number of spheres placed within a certain radial distance from a given sphere. This is a measure of how efficiently local dense agglomerate of spheres are formed and therefore we called it *packing efficiency* $`ϵ`$. It is well known that no more than 12 spheres can be found in contact with one sphere (the ‘Kissing number’ ), but the upper limit for the number of spheres within a given radial distance is, in general, unknown. Here, we empirically investigate the local environment surrounding each sphere starting from the spheres in contact up to the ones at a radial distance of two bead diameters. Fig.2 shows the average, the maximum and the minimum numbers of neighbours within a given radial distance from any sphere in $`𝐆`$. Interestingly, we observe that, in the regions around $`r1.4`$ and $`1.7`$, disordered packings show better average packing efficiency than fcc. Moreover, we observe that in all the samples the upper limits are above or equal to the fcc packing-efficiency in a large range of the radial distances. Let us stress that, these surprisingly high efficiencies in disordered packings could be the inner driving cause for the occurence of such disordered assemblies.

Radial Distribution Function: The radial distribution function $`g(r)`$ (shown in Fig.3) is calculated as the average number of sphere centres, within a radial distance $`r\mathrm{\Delta }/2`$ and $`r+\mathrm{\Delta }/2`$, divided by $`cr^2`$. The constant $`c`$ is fixed by imposing that asymptotically $`g(r)1`$ for $`r\mathrm{}`$. (We have verified that different choices of $`\mathrm{\Delta }`$ within a broad range of $`10^4`$ to $`10^2`$ (in bead diameters units), lead to almost indistinguishable results.) The detail of the close-neighbour region (Fig.3) shows that the two peaks at $`r=\sqrt{3}`$ and $`r=2`$ (in bead diameters) both increase in height with the packing density. This is an indication that there is an increasing number of configurations in the contact network  with edge-sharing in-plane triangles (peak at $`\sqrt{3}`$) and with three (or more) aligned spheres (peak at $`2`$). This might indicate an increasing organisation in the packing structure but, on the other hand, no clear signs of crystallisation were detected (see below).

For all the samples investigated, we found that in the range $`r_0r1.4`$, the radial distributions follow a power law behaviour: $`g(r)(rr_0)^\alpha `$, with exponents $`\alpha `$ between 0.3 and 0.45 and the singularity at $`r_0=1.03`$ (Fig.3). A similar behaviour, but with $`\alpha =0.5`$ and $`r_0=1`$, was reported in  for molecular dynamics simulations. In Fig.3 it is also highlighted the growing trend of the exponent $`\alpha `$ with the density $`\rho `$. To our knowledge, this is the first time that such behaviours (still theoretically unexplained) are observed in experiments.

Shell Analysis: The sphere neighbourhood environment has been also studied topologically by performing a *shell analysis* on the contact network and calculating the topological density ($`T_D`$) . We observe that the number of spheres at a given topological distance ($`j`$ ) from a central one follows quite accurately the quadratic law: $`K_j=3T_Dj^2+c_1j+c_0`$. The best fit values of $`T_D`$’s are reported in Table 1. (This analysis was performed by defining in contact spheres at radial distance within $`rd+v/2`$; the dependence of $`T_D`$ on the choice of the contact-neighbour criteria will be discussed in a forthcoming paper). We note that the topological densities increase with the sample-densities and are consistently smaller than $`10/3=3.33\mathrm{}`$, which is the theoretical lower limit for Barlow packings . A view of the topological shell structure is shown in Fig.1.

Local Order: We estimate the local packing orientation by associating a set of spherical harmonics to the vectors $`\stackrel{}{r}`$ between a sphere and the set of neighbouring spheres ($`|\stackrel{}{r}|1.2`$ \- in beads diameters): $`Y_{l,m}(\theta (\stackrel{}{r}),\varphi (\stackrel{}{r}))`$ (with $`\theta (\stackrel{}{r})`$ and $`\varphi (\stackrel{}{r})`$ the polar and azimuthal angles associated with $`\stackrel{}{r}`$). Following , a set of rotationally invariant quantities are calculated:

$$Q_l=\sqrt{\frac{4\pi }{2l+1}\underset{m=l}{\overset{l}{}}|\underset{i𝐆}{}Y_{l,m}(\theta (\stackrel{}{r}_i),\varphi (\stackrel{}{r}_i))|^2},$$
(1)
We observe that the $`Q_l`$’s (with $`l=4,6`$) calculated over many sub-volumes of the samples tend to vanish as the sub-volumes become larger. This indicates that these packings are macroscopically like *isotropic liquids*  and across the entire system there is no prevailing orientational order. On the other hand, even if *globally* there are no signs of order or organisation, this method is useful to investigate the orientational properties of the *local* configurations. In particular, for each sphere in the central region we measured the coefficients $`Q_4`$ and $`Q_6`$ calculated from the angles between a sphere and all its neighbour up to a radial distance $`r=1.2`$. We first search for signatures of local crystalline orientational order measuring the fraction of local configurations with $`(Q_4,Q_6)`$ in a region close to the fcc $`(0.1909,0.5745)^{(fcc)}`$ and hcp $`(0.0972,0.4848)^{(hcp)}`$ order. We observe that a relatively small fraction of configurations have $`Q_4`$ and $`Q_6`$ within a range $`\pm 0.05`$ of their ideal fcc and hcp values (see ‘fcc’ and ‘hcp’ in Table 1). We also verify that other special configurations like Icosahedral, bcc and Simple Cubic, have significantly lower occurrences. On the other hand, we observe that large percentages (between 27% and 37%) of the local configurations share a common signature in their orientational order with $`0.15Q_40.25`$ and $`0.4Q_60.5`$ (see ‘dis’ in Table 1). To understand whether this is an indication of a ‘typical’ local disordered packing around $`(0.20,0.45)^{(dis)}`$ is the matter of present investigations.

Structural Arrest: We study a quantity which is relevant for system dynamics: the *escape probability* which is the probability that a sphere can moves outward from a given local configuration without re-adjusting the positions of its first neighbours . This quantity is calculated by constructing circles through the centres of the three spheres corresponding to the three faces incident at each vertex in the Voronoï polyhedron. If one of these circles has a radius larger than $`d`$, it implies that the central sphere can pass through that neighbouring configuration and move outward from its local position without displacing the first-neighbours. In other words, the neighbouring ‘cage’ is *open* if at least one radius is lager than $`d`$; vice versa the cage is *closed* when all radii are smaller than $`d`$. The escape probability is defined as the fraction of open cages. We find that all the samples with $`\rho >0.6`$ have zero escape probability in $`𝐆`$. Whereas, the two samples with $`\rho =0.586`$ and $`\rho =0.595`$, have very small fractions of open cages (0.1% and 0.6%, respectively).

Following a theoretical approach recently developed by one of the authors , the packing realization can be considered as an inherent structure and the (thermo)dynamical approach toward this configuration can be reconstructed by virtually decreasing the sphere-diameters reducing the effective density. (Here, the main underlying assumption is that the system dynamics before the structural arrest develops through configurations around the final inherent-structure realization.) We find that, for all the samples with $`\rho >0.6`$, the escape probabilities go to zero at effective densities between 0.61 and 0.63. This strongly suggests that around $`\rho 0.62\pm 0.01`$ an important phase in the system dynamics reaches an end: above this density, local readjustments involving only the displacement of a single sphere are forbidden and the system compaction can proceed only by involving the collective and correlated readjustment of larger set of spheres.

Acknowledgements The authors thank Ajay Limaye for the preparations of Fig.1. Senden gratefully acknowledges the ARC for his Fellowship. This work was partially supported by the ARC discovery project DP0450292.