text:
# Correlated 𝚲_𝒄⁺𝚲̄_𝒄⁻ production in 𝒆⁺𝒆⁻ annihilations at √𝒔∼10.5 GeV

## I Introduction

At high energy, quark fragmentation can be calculated using Perturbative QCD (PQCD). The results are often referred to as QCD showers . However, these calculations require, for comparison with experimental results, some model or experimental input to take into account the eventual hadronization of gluons and $`q\overline{q}`$ pairs into actual hadrons of known mass.

Several models have thus been developed, mostly based on a chromodynamic string connecting the initial $`q_0`$ and $`\overline{q}_0`$ , where the subscript indicates that these are the primary quarks formed by $`e^+e^{}`$ annihilation. When stretched because of the initial momentum of $`q_0`$ and $`\overline{q}_0`$, the string breaks by creation of secondary $`q\overline{q}`$ pairs, generating two substrings. The process iterates then for the new substrings until the physical hadrons are produced. This last stage is the critical one where QCD loses its predictive power.

The model most frequently used to simulate the process $`e^+e^{}q_0\overline{q}_0hadrons`$ is the QCD inspired Lund String Model (LSM) , implemented in the Jetset Monte Carlo simulation package . In order to solve the problem of generating hadrons of known mass, the Jetset implementation of LSM, at each step of its iterative process, splits a $`q\overline{q}`$ string into a hadron of known quantum numbers and mass and a remainder string having the leftover quantum numbers and energy-momentum. When, finally, the remainder string has a mass below a suitable limit, this string is made into two known hadrons that, together, carry the leftover quantum numbers.

In this paper we report on correlations between charm baryons in the fragmentation of a $`c\overline{c}`$ system produced in the reaction $`e^+e^{}c\overline{c}`$ at $`\sqrt{s}10`$ GeV, where the annihilation occurs in a largely low-$`Q^2`$, non-perturbative regime. In contrast to inclusive single-particle production, compensation of baryon and charm number is a more subtle aspect of quark fragmentation modeling. One obvious question is whether baryon compensation occurs locally (e.g., small rapidity difference between baryon and antibaryon) or globally (large rapidity difference). Previous studies of $`\mathrm{\Lambda }\overline{\mathrm{\Lambda }}`$ production at $`\sqrt{s}`$=90 GeV found that in events containing both a $`\mathrm{\Lambda }`$ and a $`\overline{\mathrm{\Lambda }}`$, the two particles tended to be produced at very similar rapidities .

Consider the case in which a charm baryon is produced in the first step of fragmentation (e.g. $`e^+e^{}c\overline{c};c\mathrm{\Lambda }_cX`$); it is possible that both baryon and quantum charm numbers be compensated in the opposite hemisphere (e.g. $`e^+e^{}c\overline{c};c\mathrm{\Lambda }_cX,\overline{c}\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}X`$). In the limit that the $`c`$ and $`\overline{c}`$ quarks fragment independently (presumably true at some sufficiently high energy), this type of correlation is not expected. The Jetset implementation of the Lund model, because of the mechanism outlined above, does not produce such a correlation.

In the very reasonable approximation of neglecting string splitting by $`c\overline{c}`$ tunneling from the vacuum, at $`\sqrt{s}10`$ GeV all observed charmed particles must contain a primary quark. We take advantage of this to discriminate between independent vs. correlated fragmentation models. If we assume that primary quarks fragment independently, then the number of times that we find a $`\mathrm{\Lambda }_c^+`$ baryon opposite a $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ antibaryon in an event (i.e., $`cos\theta (\mathrm{\Lambda }_c^+,\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{})<`$0, denoted “$`\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$”), scaled to the total number of observed $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ (denoted “$`\frac{\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}{\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}`$”), should be equal to the number of times that we find a $`\mathrm{\Lambda }_c^+`$ baryon opposite any other anti-charmed hadron $`\overline{H}_c`$, scaled to the total number of observed anti-charmed hadrons ($`\frac{\mathrm{\Lambda }_c^+|\overline{H_c}}{\overline{H_c}}`$). In this analysis, we will check the equality of these ratios; by comparing ratios in this way, we cancel many experimental systematics. Schematically, we are comparing the following event topologies:

$`\begin{array}{cccccccc}& & & \overline{c}& c& & & \\ & & \overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}& & & \mathrm{\Lambda }_c^+& & \end{array}`$

$`\begin{array}{cccccccc}& & & \mathrm{𝑣𝑠}& & & & \end{array}`$

$`\begin{array}{cccccccc}& & & \overline{c}& c& & & \\ & & \overline{H_c}& & & \mathrm{\Lambda }_c^+& & \end{array}`$

Specifically, we compare the rate of $`\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ production to the rate of $`\mathrm{\Lambda }_c^+|\overline{D^0}`$ and $`\mathrm{\Lambda }_c^+|D^{}`$ production. We reconstruct $`\overline{D}^0`$’s and $`D^{}`$’s through the well-measured decay modes $`\overline{D}^0K^+\pi ^{}`$ and $`D^{}K^+\pi ^{}\pi ^{}`$, respectively. $`\mathrm{\Lambda }_c^+`$’s are fully reconstructed in the decay modes $`\mathrm{\Lambda }_c^+pK^{}\pi ^+`$, $`\mathrm{\Lambda }_c^+pK_s^0`$, $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+`$, $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^{}\pi ^+`$, $`\mathrm{\Lambda }_c^+pK_s^0\pi ^+\pi ^{}`$,***Charge conjugation is implicit. and partially reconstructed through $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }`$X.According to JETSET 7.3 simulations, 95% of the $`\mathrm{\Lambda }`$’s opposite an anti-charm tag are charmed baryon daughters (predominantly $`\mathrm{\Lambda }_c`$) after imposing a 1 GeV/c minimum momentum requirement on the $`\mathrm{\Lambda }`$ ($`p_\mathrm{\Lambda }>`$1 GeV/c). $`\mathrm{\Sigma }_c`$’s are also studied; these are reconstructed in $`\mathrm{\Sigma }_c^{++}\mathrm{\Lambda }_c^+\pi ^+`$ and $`\mathrm{\Sigma }_c^0\mathrm{\Lambda }_c^+\pi ^{}`$.

Under the assumption of independent fragmentation, we expect the relative production ratios to satisfy:

$`\frac{\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}{\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}÷\frac{\mathrm{\Lambda }_c^+|\overline{D^0}}{\overline{D^0}}`$ = 1

and,

$`\frac{\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}{\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}÷\frac{\mathrm{\Lambda }_c^+|D^{}}{D^{}}`$ = 1.

## II Apparatus and Event Selection

This analysis was performed using the CLEO II and the upgraded CLEO II.V detectors operating at the Cornell Electron Storage Ring (CESR) at center-of-mass energies $`\sqrt{s}`$ = 10.52–10.58 GeV. For 4.5 fb-1 of the data used for this analysis (“CLEO-II” data), measurements of charged particle momenta were made with three nested coaxial drift chambers consisting of 6, 10, and 51 layers, respectively. In a subsequent upgrade (“CLEO-II.V”), the inner tracking chamber was replaced with a high-precision silicon detector (corresponding to the remaining 9.1 fb-1 of the data used for this analysis). The entire tracking system fills the volume from $`r`$=3 cm to $`r`$=1 m, with $`r`$ the radial coordinate relative to the beam ($`\widehat{z}`$) axis. This system is very efficient ($`ϵ`$98%) for detecting tracks that have transverse momenta ($`p_T`$) relative to the beam axis greater than 200 MeV/c, and that are contained within the good fiducial volume of the drift chamber ($`|\mathrm{cos}\theta |<`$0.94, with $`\theta `$ defined as the polar angle relative to the beam axis). This system achieves a momentum resolution of $`(\delta p/p)^2=(0.0015p)^2+(0.005)^2`$ ($`p`$ is the momentum, measured in GeV/c). Pulse-height measurements in the main drift chamber provide specific ionization resolution of 5.0% (CLEO II.V) or 5.5% (CLEO II) for Bhabha events, giving good $`K/\pi `$ separation for tracks with momenta up to 700 MeV/c and separation of order 2$`\sigma `$ in the relativistic rise region above 2.5 GeV/c. Outside the central tracking chambers are plastic scintillation counters, which are used as a fast element in the trigger system and also provide particle identification information from time-of-flight measurements.

Beyond the time-of-flight system is the electromagnetic calorimeter, consisting of 7800 thallium-doped CsI crystals. The central “barrel” region of the calorimeter covers about 75% of the solid angle and has an energy resolution which is empirically found to follow:

$$\frac{\sigma _\mathrm{E}}{E}(\%)=\frac{0.35}{E^{0.75}}+1.90.1E;$$
(1)
$`E`$ is the shower energy in GeV. This parameterization includes noise effects, and translates to an energy resolution of about 4% at 100 MeV and 1.2% at 5 GeV. Two end-cap regions of the crystal calorimeter extend solid angle coverage to about 95% of $`4\pi `$, although energy resolution is not as good as that of the barrel region. The tracking system, time-of-flight counters, and calorimeter are all contained within a superconducting coil operated at 1.5 Tesla. An iron flux return interspersed with proportional tubes used for muon detection is located immediately outside the coil and in the two end-cap regions.

The event sample used for this measurement is comprised of 9.1 $`\mathrm{fb}^1`$ of data collected at the $`\mathrm{{\rm Y}}`$(4S) resonance and 4.5 $`\mathrm{fb}^1`$ of data collected about 60 MeV below the $`\mathrm{{\rm Y}}`$(4S) resonance. Approximately $`18\times 10^6`$ continuum cc̄ events are included in this sample. Charged track candidates for protons, kaons or pions must pass the following restrictions:

(a) The track must have an impact parameter relative to the estimated event vertex less than 5 mm in a plane perpendicular to the beam axis ($`\widehat{r}\widehat{\varphi }`$) and no more than 5 cm along the beam axis. The estimated event vertex is obtained by averaging the $`e^+e^{}`$ interaction point over a full run.

(b) The track has specific ionization information consistent (at the 99% confidence level) with its assumed particle identity.

(c) The track must have momentum greater than 100 MeV/c.

All reconstructed charmed hadrons must have momentum greater than 2.3 GeV/c to ensure that there is no contamination from $`B`$-meson decays to charm.

### A Production Ratios

We define the single tag yield to be the number of reconstructed events containing one anti-charmed hadron $`\overline{H}_c`$. This yield (Table I) is determined by fitting a double-Gaussian signal over a smooth, low-order polynomial background functional. The signal function is the sum of two Gaussian functions, one narrow and one broad, and is a better representation of the expected signal line shape than a single Gaussian function.

The number of double tags (Table II) is defined as the number of events in which two specific particles are reconstructed in opposite hemispheres (i.e., the opening angle between the two particles must be greater than 90 degrees). The total number of double tags for $`\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$, $`\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}`$, $`\mathrm{\Lambda }_c^+|\overline{D}^0`$ and $`\mathrm{\Lambda }|\overline{D}^0`$ is extracted from two-dimensional invariant mass plots, shown in Figures 1, 2, 3, and 4 respectively.

The total correlated double-tag yield is first determined by performing a sideband subtraction in the two-dimensional invariant mass plot. Consider, for example Fig. 1. We take one-dimensional projections of three slices in the candidate $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ recoil invariant mass “$`M_{recoil}`$” - the $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ signal region: ($`|M_{recoil}2.286|<`$0.03 $`GeV/c^2`$) and the two $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ sideband regions: ($`0.04<|M_{recoil}2.286|<`$0.07 $`GeV/c^2`$). We then subtract the $`\mathrm{\Lambda }_c^+`$ distributions from the $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ sidebands from that of the $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ signal region. Figure 5. shows the signal distribution and the sum of the two sideband distributions. In performing these fits, the double-Gaussian signal shapes are constrained using the parameters determined from fits to the single-tag sample.

We also perform a single fit in two dimensions to extract the signal yields. In this latter fit, a two-dimensional Gaussian signal is used to parametrize the peak region, two single Gaussians are used to fit the “ridges” away from the peak region (corresponding to true charmed hadrons along one axis in association with combinatoric background on the other axis) and a two-dimensional, smooth polynomial is used to parametrize the background. Differences in i) the selection of the sideband/signal regions and the regions projected, as well as ii) the difference between the sideband-subtracted yield compared to the two-dimensional, single-fit yield give a measure of the fitting systematics for each measurement.

We define the “production rate” (Table III) as the percentage of times we find one specific particle in an event opposite a given tag (i.e., the number of double tags, given in Table II divided by the total number of single tags, presented in Table I). Production ratios are formed by taking ratios of production rates (Table IV). Note that, to compare $`\frac{\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}{\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}`$ to $`\frac{\mathrm{\Lambda }_c^+|\overline{D}}{\overline{D}}`$, we must take into account the fact that the latter ratio implicitly includes both sign combinations ($`\frac{\mathrm{\Lambda }_c^+|\overline{D}}{\overline{D}}`$ plus $`\frac{\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}|D}{D}`$), while there is only one unique way to form $`\frac{\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}{\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}`$. To compare these two rates and form production ratios (Table IV), the production rates for $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$-tags presented in the table have therefore been multiplied by this factor of two.

From Table III, we see a clear enhancement in the likelihood of producing a $`\mathrm{\Lambda }_c^+`$ opposite a $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ compared to a $`\overline{D}`$. This is observed in the case where the $`\mathrm{\Lambda }_c^+`$ is fully reconstructed, as well as tagged inclusively by a $`\mathrm{\Lambda }`$. Note that the fractional enhancement is smaller in the case when the $`\mathrm{\Lambda }_c^+`$ is tagged inclusively. This is qualitatively consistent with the expectation that events containing two charmed mesons and two baryons will produce a $`\mathrm{\Lambda }|\overline{D}`$ correlation which will inflate the denominator when we construct the ratio: $`\frac{\mathrm{\Lambda }|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}{\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}}`$ / $`\frac{\mathrm{\Lambda }|\overline{D}}{\overline{D}}`$. Non-cancelling contributions to the numerator, such as $`\mathrm{\Lambda }K\overline{N}D|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$, plus decays of charmed baryons other than $`\mathrm{\Lambda }_c`$ into $`\mathrm{\Lambda }`$, may also be present.

## III Cross Checks

A number of cross checks were conducted in order to verify the accuracy of the double tag signal extraction. These include signal extractions of two-dimensional “wrong sign” plots (expected to have zero signal yield), as well as consistency checks with Monte Carlo and studies of $`D|\overline{D}`$ correlations.

### A Null and Wrong sign signals

Neglecting the doubly Cabibbo suppressed decay $`\overline{D}^0K^{}\pi ^+`$, we expect a null signal yield from a double tag plot of ($`m_{K^{}\pi ^+}|m_{K^{}\pi ^+}`$). Similarly, we expect zero signal yield from a double tag plot of ($`m_{K^{}\pi ^+}m_{K^+\pi ^{}}|`$) in the case where the $`D^0`$ and $`\overline{D}^0`$ candidate are in the same hemisphere. All of these correlations give signal yields consistent with zero as expected (Table V).

### B Monte Carlo studies

The uncertainty in $`\mathrm{\Lambda }_c`$ production characteristics is expected to be somewhat larger than the uncertainty in $`D^0`$ production. We therefore expect the Monte Carlo and data to agree on, e.g., the number of $`D^+|D^{}`$ double tags per $`D^{}`$ ($`\frac{D^+|D^{}}{D^{}}`$), and the number of $`D^0|\overline{D^0}`$ double tags (Figure 6) per $`\overline{D}^0`$ ($`\frac{D^0|\overline{D^0}}{\overline{D}^0}`$), as shown in Table VI.The measured value of 0.0119, for instance, is qualitatively consistent with the expectation that, per $`\overline{D}^0`$ tag, 50% of the time the charm quark will produce a $`D^0`$, which is reconstructed in the $`D^0K^{}\pi ^+`$ mode ($`(D^0K^{}\pi ^+)`$0.04) with approximately 60% efficiency.

Within errors, the agreement is good. We have also further ensured that there is no bias in the signal extraction due to possible peaking of the background in the signal region, by subtracting all the true signal from a Monte Carlo double tag plot and verifying that the measured yield, after subtraction of the true generated particles, is indeed consistent with zero.

### C Check of expected uncorrelated yield

Having obtained values for production rates, we can compare these with the values expected under the assumption of two hemispheres fragmentating independently. In this (hypothetical) limit, the ratio of double-tags to single-tags should be the same as the ratio of single-tags to total charm quarks, since each simply expresses the probability of a charm quark to evolve into a particular final-state particle. Then, using $`(\mathrm{\Lambda }_cX)`$ to designate the branching fraction for $`\mathrm{\Lambda }_c`$ to decay into the final state $`X`$, $``$ as the integrated luminosity of our data set, and $`f_{(c\mathrm{\Lambda }_c)}`$ as the fraction of times a charm quark materializes as a $`\mathrm{\Lambda }_c`$ baryon, we can relate the expected number of single tags and double tags as:

$`N(\mathrm{\Lambda }_c)_{SingleTags}=\sigma (e^+e^{}c\overline{c})f_{(c\mathrm{\Lambda }_c)}(\mathrm{\Lambda }_cX)ϵ_{\mathrm{\Lambda }_cX}(2)`$and, for the probability that, given a reconstructed $`\mathrm{\Lambda }_c`$, we reconstruct an opposite $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$

$`N(\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{})_{DoubleTags}=f_{(c\mathrm{\Lambda }_c)}^{^{}}(\mathrm{\Lambda }_cX)ϵ_{\mathrm{\Lambda }_cX}^{^{}}N(\mathrm{\Lambda }_c)_{SingleTags}(3)`$

In equation (2), for the number of single $`\mathrm{\Lambda }_c^+`$ tags, we use half the total number of observed $`\mathrm{\Lambda }_c^++\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ (70199/2). In equation (3), the term $`f_{(c\mathrm{\Lambda }_c)}^{^{}}`$ indicates that the fraction of times that a charm quark produces a $`\mathrm{\Lambda }_c`$ may be different if the event already contains a $`\overline{\mathrm{\Lambda }}_c`$ produced from the corresponding anticharm quark. This is, of course, exactly the correlation factor we wish to ultimately measure. The prime on $`ϵ_{\mathrm{\Lambda }_cX}^{^{}}`$ indicates that the efficiency for reconstructing a $`\mathrm{\Lambda }_c`$ may be higher for events in which a $`\overline{\mathrm{\Lambda }}_c`$ has already been reconstructed due to geometrical correlations – since the charm and anticharm quarks are back to back in an $`e^+e^{}`$ cc̄ event, reconstruction of one charmed particle ensures that the corresponding antiparticle is in a good acceptance region of the detector. Dividing equation (3) by equation (2) and multiplying each side by $`N(\mathrm{\Lambda }_c)_{SingleTags}`$ we obtain, for the number of double tags:

$`N(\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{})_{DoubleTags}={\displaystyle \frac{(N(\mathrm{\Lambda }_c)_{SingleTags})^2f_{(c\mathrm{\Lambda }_c)}^{^{}}ϵ_{\mathrm{\Lambda }_cX}^{^{}}}{\sigma (e^+e^{}c\overline{c})f_{(c\mathrm{\Lambda }_c)}ϵ_{\mathrm{\Lambda }_cX}}}(4)`$
The difference between $`f_{(c\mathrm{\Lambda }_c)}`$ and $`f_{(c\mathrm{\Lambda }_c)}^{^{}}`$ represents the enhancement in production of the $`\mathrm{\Lambda }_c`$ when a $`\overline{\mathrm{\Lambda }_c}`$ is present. With $`\frac{ϵ_{\mathrm{\Lambda }_cX}^{^{}}}{ϵ_{\mathrm{\Lambda }_cX}}ϵ_{geometry}`$ and $`\frac{f_{(c\mathrm{\Lambda }_c)}^{^{}}}{f_{(c\mathrm{\Lambda }_c)}}f_{correlated}`$ (=1 if independent fragmentation holds), we can estimate our expected double tag to single tag ratio as follows: From the luminosity, we calculate the number of $`e^+e^{}c\overline{c}`$ events using as inputs the $`e^+e^{}q\overline{q}`$ cross-section (3.3 nb), and taking $`\frac{e^+e^{}c\overline{c}}{e^+e^{}q\overline{q}}`$=0.4. Assuming that $`f_{correlated}`$ is equal to unity for $`D|\overline{D}`$ events, we can solve Eqn. (4) for $`ϵ_{geometry}`$ by comparing the total number of $`D^0|\overline{D^0}`$ and $`D^+|D^{}`$ double tags we would expect to find (assuming independent fragmentation) with the actual number of $`D^0|\overline{D^0}`$ and $`D^+|D^{}`$ double tags. This is calculated to be

$$ϵ_{geometry}=\frac{(2099)(0.4)(3.3)(13.6\times 10^6)}{(352294/2)^2}1.21\pm 0.04(statisticalerrorsonly),$$

using the total luminosity $``$=13.6 fb-1, and the single tag and double tag values for our $`D^0`$ and $`D^0|\overline{D}^0`$ samples, respectively. For the $`D^+`$ and $`D^+|D^{}`$ samples, we obtain $`ϵ_{geometry}=1.21\pm 0.09`$. Using a value of $`ϵ_{geometry}1.21`$, a similar calculation can be performed for the expected number of $`\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ double tags, based on the total number of cc̄ events and the single tag yield, giving an expected value of 83.0$`\pm `$2.6 (statistical error only) double tag events. This number is then compared to the measured number of $`\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ double tags to obtain an estimate of $`f_{correlated}=\frac{253\pm 37}{83.0\pm 2.6}=3.04\pm 0.45`$. Despite the roughness of this approach, the observed enhancement is consistent with our calculated production ratios.

## IV Study of the $`\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ signal characteristics

We have compared some of the production characteristics of the observed $`\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ signal with the $`\mathrm{\Lambda }_c^+|\overline{D^0}`$ and $`\mathrm{\Lambda }_c^+|D^{}`$ signals. No obvious difference is found in either the $`\mathrm{\Lambda }_c^+`$ momentum spectrum (or $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$, in the charge conjugate case) or the polar angle distribution of the $`\mathrm{\Lambda }_c^+`$ (see Figs. 7 and 8, respectively) for the $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$-tagged or the $`\overline{D}`$-tagged samples. We also find that the observed $`\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ cross-section is statistically the same for the data taken on the $`\mathrm{{\rm Y}}`$(4S) resonance as data taken just below the $`\mathrm{{\rm Y}}`$(4S). As expected (Figure 9) if the charmed baryons are following the direction of the original charm/anticharm quarks, the observed opening angle between the $`\mathrm{\Lambda }_c^+`$ and $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ peaks at 180. Finally, to investigate the possibility that the $`\mathrm{\Lambda }_c^+|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ signal was associated with the production of a 4-baryon system, we plot the double tag yield in cases where there are well-identified protons found in the same event (Fig. 10). No signal is observed in such a case.

## V $`\mathrm{\Sigma }_c`$ $`(\frac{1}{2})^+`$ production

We have also investigated the possibility of correlated $`\mathrm{\Sigma }_c|\overline{\mathrm{\Sigma }}_c`$ production. Of particular interest here is the possibility of observing an enhanced production rate relative to our observed $`\mathrm{\Lambda }_c|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ correlation. Such an enhancement may be indicative of spin correlations at the first rank in fragmentation. $`\mathrm{\Sigma }_c`$’s are reconstructed through their decay mode into $`\mathrm{\Lambda }_c`$ and a soft pion $`\pi _s`$: $`\mathrm{\Sigma }_c\mathrm{\Lambda }_c\pi _s`$. For the highest efficiency and best resolution, only charged pions $`\pi _s^\pm `$ are used in reconstructing $`\mathrm{\Sigma }_c`$’s. Although statistically limited, a direct search for $`\mathrm{\Sigma }_c|\overline{\mathrm{\Sigma }}_c`$ double tags yields 6 events in the signal region (Figure 11); the extrapolated background under the signal is $`1.8\pm 0.3`$ events.

Using the total number of $`\mathrm{\Sigma }_c`$ single tags observed ($`3804\pm 185`$), the total number of $`\mathrm{\Lambda }_c`$ single tags observed ($`70199\pm 1604`$), and the total number of $`\mathrm{\Lambda }_c|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ double tags observed ($`252.5\pm 37.4`$), we can estimate the number of expected $`\mathrm{\Sigma }_c|\overline{\mathrm{\Sigma }}_c`$ double tags (assuming that the $`\mathrm{\Sigma }_c|\overline{\mathrm{\Sigma }}_c`$ correlated production rate is the same as the correlated $`\mathrm{\Lambda }_c|\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ production rate) as: 252.5$`\times (\frac{3804}{70199})^2`$=$`0.74\pm 0.11`$. We therefore see a statistically very limited indication of $`\mathrm{\Sigma }_c|\overline{\mathrm{\Sigma }}_c`$ correlations, suggesting a spin correlation at the first step of fragmentation. More data are necessary to elucidate this situation.

## VI Systematics

We expect most of the systematic uncertainties associated with measurement of both numerator and denominator of the production rates to cancel in calculation of the production ratios. The sensitivity to the definition of “signal” and “sideband” mass regions in the one- and two-dimensional subtraction was determined by calculating the yields as we vary these parameters. We assign a systematic error of 6% due to “signal parameterization uncertainty” based on this study. We have also performed a second (complete) set of fits using two-dimensional Gaussian parameterizations of the signal and two-dimensional polynomial parameterizations of the background. The two techniques show excellent agreement, and are typically within 5% of each other, for each two-dimensional extraction. When we vary the selection criteria (e.g., particle identification criteria, and minimum momentum requirements) used to define our hadrons ($`\mathrm{\Lambda }_c`$, $`D^0`$, and $`D^+`$), we observe a maximum 7% variation in the production ratios we calculate. As an example of our investigation of the momentum dependence of the $`\mathrm{\Lambda }_c^+\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ production ratio enhancement, Table VII shows the production ratio in the two-dimensional $`(p_{\mathrm{\Lambda }_c},p_{\overline{H}_c})`$ space.

The systematic errors are added, mode-by-mode, to determine the overall systematic error for each calculated production ratio (Table IV).

## VII Summary

The measured values of the production ratios imply that a $`\mathrm{\Lambda }_c`$ is roughly three times more likely to be produced opposite a $`\overline{\mathrm{\Lambda }_c}`$ than opposite either a $`\overline{D}^0`$ or a $`D^{}`$. These results indicate strong evidence in support of correlated production of the $`\mathrm{\Lambda }_c`$. Assuming that 6% of charm quark fragmentation at $`\sqrt{s}`$10 GeV results in a charmed baryon, our result implies that approximately 20% of all $`\mathrm{\Lambda }_c`$’s are produced in association with a $`\overline{\mathrm{\Lambda }}_{c}^{}{}_{}{}^{}`$ in the opposite hemisphere. The effect is not produced by the Jetset simulation but does not contradict the Lund String Model. In fact if the primary $`c\overline{c}`$ were broken at very small proper time by diquark-antidiquark tunneling from the vacuum, producing a (possibly excited) charm baryon system and anti-charm antibaryon system, a correlation of the type observed here would be generated.

## VIII Acknowledgments

We gratefully acknowledge the effort of the CESR staff in providing us with excellent luminosity and running conditions. M. Selen thanks the PFF program of the NSF and the Research Corporation, A.H. Mahmood thanks the Texas Advanced Research Program, F. Blanc thanks the Swiss National Science Foundation, and E. von Toerne thanks the Alexander von Humboldt Stiftung for support. This work was supported by the National Science Foundation, the U.S. Department of Energy, and the Natural Sciences and Engineering Research Council of Canada.


format: text

text:
# Nearly horizon skimming orbits of Kerr black holes

## I Introduction

The space-based gravitational-wave detector LISA  is being designed to make very precise measurements of the characteristics of black hole spacetimes. One source that is particularly well-suited for such measurements is the inspiral of a small (mass $`\mu =110M_{}`$) compact body into a massive ($`M=10^{57}M_{}`$) Kerr black hole. Depending upon the values of $`M`$, $`\mu /M`$, and the hole’s spin $`a`$, such an inspiral will spend several months to several years in LISA’s frequency band , radiating $`10^510^6`$ gravitational-wave cycles. By accurately measuring these cycles, LISA should be able to build a “map”  of the spacetime, testing in detail the predictions of general relativity.

As the community begins developing strategies for analyzing LISA’s datastream, it is important to revisit and carefully analyze the sources one expects to measure. For extreme mass ratio inspirals, this means understanding more deeply the character of Kerr black hole orbits and the nature of gravitational radiation reaction deep in the Kerr metric’s strong field. Of particular interest are features that might simplify data analysis (which is likely to be very difficult, given the many cycles that must be tracked) or that might constitute a strong signature of the spacetime. In this paper, I analyze what might be such a feature — a unique signature of the inspiral of a body through the extreme strong field of rapidly rotating black holes.

The key piece of this analysis was first discussed by Wilkins. Reference  contains a detailed examination of circular orbits of extreme ($`a=M`$)***I use units where $`G=1=c`$. Kerr black holes. (In this context, “circular orbit” means “orbit of constant Boyer-Lindquist coordinate radius”.) One very interesting result given in  is the existence of “horizon skimming” orbits. These are circular orbits of varying inclination angle at the same coordinate radius as the event horizon, $`r=M`$. (Despite being at the same radial coordinate, one can show that these orbits have distinct proper separation, and in particular lie outside the event horizon; see Ref. , particularly Fig. 2.) An extremely interesting feature of the horizon skimming orbits is that as an orbit’s inclination angle $`\iota `$ is increased, its angular momentum component $`L_z`$ likewise increases: $`L_z/\iota >0`$ for the horizon skimming orbits. This property holds over a sizable range of radius, out to $`r1.8M`$. I will call the full set of orbits for which $`L_z/\iota >0`$ the “nearly horizon skimming” (NHS) orbits. This defining property of NHS orbits is opposite to weak-field intuition. For example, in Newtonian theory, orbits at constant radius have $`L_z=|\stackrel{}{L}|\mathrm{cos}\iota `$ (where $`|\stackrel{}{L}|`$ is the same for all orbits at radius $`r`$), which decreases as $`\iota `$ increases.

Intuition from Newtonian theory is highly suspect in the strong field of black holes. However, the true behavior of $`L_z`$ as a function of $`\iota `$ is qualitatively the same as in Newtonian theory over a wide range of orbital radii and spins. For example, when $`a=0.95M`$, $`L_z/\iota <0`$ for all circular orbits, so there are no NHS orbits when $`a=0.95M`$. It turns out that NHS orbits can only exist when $`a0.952412M`$. This is the smallest spin for which stable orbits come close enough to the hole’s event horizon that $`L_z/\iota `$ can switch sign: the property $`L_z/\iota >0`$ arises because, very close to the black hole, all physical processes become “locked” or “frozen” to the hole’s event horizon . In particular, their orbital motion locks to the horizon’s spin. This locking dominates the “Keplerian” tendency of an orbit to move more quickly at smaller radius ($`v_{\mathrm{Kepler}}=\sqrt{M/r}`$), forcing a body to actually slow in the innermost orbits. The locking is particularly strong for the most-bound orbits. As I show below, the least-bound orbits do not strongly lock to the black hole’s spin until they have very nearly reached the innermost orbit. The NHS orbit’s $`L_z/\iota `$ behavior follows from the fact that the most-bound orbit locks to the horizon more readily than the least-bound orbit.

This behavior could have interesting observational consequences. It has been well-understood for some time that the inclination angle of an inspiraling body increases due to gravitational-wave emission . Gravitational waves carry $`L_z`$ away from the orbit, so that $`dL_z/dt<0`$. Since “normal” orbits have $`L_z/\iota <0`$, it follows that $`d\iota /dt>0`$. If $`L_z/\iota `$ switches sign, then $`d\iota /dt`$ will switch sign as well: an inspiraling body will evolve toward an equatorial orbit. If the change in $`\iota `$ is large, it could have a large effect on the gravitational waveform. For example, the spin-orbit modulation of the wave’s amplitude and phase is due to motion in the $`\theta `$ coordinate. This modulation will be reduced as the body’s $`\theta `$ motion reduces.

Since the size of NHS orbits is significant near the “astrophysically maximal” value $`a=0.998M`$ (the value at which a hole’s spin tends to be buffered due to photon capture from thin disk accretion; see Ref. ), astrophysical black holes might spin quickly enough for NHS orbits to play some important role. This motivates a careful analysis to see what role, if any, NHS orbits might play in gravitational-wave sources. I use the code described in Ref.  (which uses the Teukolsky and Sasaki-Nakamura equations  to compute the flux of energy and $`L_z`$ carried away from the orbit by gravitational waves) to study how a small body’s motion evolves as it spirals through the NHS region. By computing the change $`(dr/dt,d\iota /dt)`$ at a large number of points, it is straightforward to construct the inspiral trajectory for a small body . I find that the total change in inclination angle as a body spirals through the NHS region is very small — at most, $`\delta \iota 1^{}2^{}`$. (See Fig. 5; note that the shape of the curves in this figure are independent of the values of $`\mu `$ and $`M`$, although the timescales strongly depend on $`\mu `$ and $`M`$.)

The inspiral code relies on the fact that circular orbits remain circular as they evolve on an adiabatic timescale . It is thus explicitly restricted to the evolution of circular orbits, and cannot say anything about the evolution of eccentric orbits, which are much more realistic as LISA sources . Based on the leading order corrections to $`d\iota /dt`$ seen from a post-Newtonian analysis , I speculate that the change in $`\iota `$ might be much larger when eccentric, inclined orbits evolve through the strong field. Verifying this, however, will require a strong-field radiation reaction formalism that can evolve generic Kerr orbits: orbits that are inclined and eccentric. Such a formalism may need to be based on a local radiation reaction force , although recent discussion suggests it may be possible to evolve generic orbits using gravitational-wave fluxes alone . The possibility that inspiral through the extreme strong field may leave an observationally significant imprint on the system’s gravitational waveform will hopefully motivate future activity and progress on this problem.

Throughout this paper, the quantities $`t`$, $`r`$, $`\theta `$, and $`\varphi `$ refer to the Boyer-Lindquist coordinates. A prime on any quantity denotes $`/r`$. (Note that in Ref.  I erroneously wrote that prime denotes $`d/dr`$.) Section II reviews the properties of circular Kerr orbits, providing formulas that are useful for describing their conserved quantities $`E`$, $`L_z`$, and $`Q`$ in the very strong field. These are used in Sec. III to study the NHS orbits. The NHS orbits are developed and mapped out as functions of spin and radius in Sec. III A. Section III B then re-examines them from the viewpoint of the “zero angular momentum observer”, or ZAMO. The ZAMO makes local measurements of the orbital properties, and can see that their defining behavior $`L_z/\iota >0`$ arises due to the locking of NHS orbits to the black hole’s spin. In Sec. IV, I examine the trajectory of a body that is inspiraling through the NHS region under gravitational-wave emission. Here I show the very small change in $`\iota `$ as the body spirals in, and argue that eccentricity might impact this result greatly. Some concluding discussion is given in Sec. V.

## II Circular orbits of Kerr black holes

Geodesic orbits of a Kerr black hole with mass $`M`$ and spin per unit mass $`a`$ are governed by the following four equations :

$`\mathrm{\Sigma }^2\left({\displaystyle \frac{dr}{d\tau }}\right)^2`$ $`=\left[E(r^2+a^2)aL_z\right]^2\mathrm{\Delta }\left[r^2+(L_zaE)^2+Q\right]R,`$ (2)
$`\mathrm{\Sigma }^2\left({\displaystyle \frac{d\theta }{d\tau }}\right)^2`$ $`=Q\mathrm{cot}^2\theta L_z^2a^2\mathrm{cos}^2\theta (1E^2)\mathrm{\Theta }^2,`$ (3)
$`\mathrm{\Sigma }\left({\displaystyle \frac{d\varphi }{d\tau }}\right)`$ $`=\mathrm{csc}^2\theta L_z+aE\left({\displaystyle \frac{r^2+a^2}{\mathrm{\Delta }}}1\right){\displaystyle \frac{a^2L_z}{\mathrm{\Delta }}},`$ (4)
$`\mathrm{\Sigma }\left({\displaystyle \frac{dt}{d\tau }}\right)`$ $`=E\left[{\displaystyle \frac{(r^2+a^2)^2}{\mathrm{\Delta }}}a^2\mathrm{sin}^2\theta \right]+aL_z\left(1{\displaystyle \frac{r^2+a^2}{\mathrm{\Delta }}}\right).`$ (5)
The quantities $`E`$, $`L_z`$, and $`Q`$ (“energy”, “$`z`$-component of angular momentum”, and “Carter constant”) specify a family of orbits, and are conserved along any orbit of the family. Here, $`\mathrm{\Sigma }=r^2+a^2\mathrm{cos}^2\theta `$ and $`\mathrm{\Delta }=r^22Mr+a^2`$. Eqs. (2) and (3) have been divided by $`\mu ^2`$, and Eqs. (4) and (5) by $`\mu `$ (where $`\mu `$ is the mass of a small body in an orbit); $`E`$, $`L_z`$, and $`Q`$ are thus the specific energy, angular momentum and Carter constant. Also, $`a0`$; prograde and retrograde orbits are distinguished by the orbit’s tilt angle rather than the sign of the hole’s spin.

A circular orbit must satisfy $`R=0`$, $`R^{}=0`$; to be stable, it also must satisfy $`R^{\prime \prime }<0`$. These conditions are met for some set of orbits everywhere outside the innermost stable circular orbit, or ISCO. The ISCO lies at

$`r_{\mathrm{ISCO}}/M`$ $`=`$ $`3+Z_2\left[(3Z_1)(3+Z_1+2Z_2)\right]^{1/2},`$ (6)
$`Z_1`$ $`=`$ $`1+\left[1(a/M)^2\right]^{1/3}\left[(1+a/M)^{1/3}+(1a/M)^{1/3}\right],`$ (7)
$`Z_2`$ $`=`$ $`\left[3(a/M)^2+Z_1^2\right]^{1/2}.`$ (8)
The ISCO varies from $`r=6M`$ for a Schwarzschild black hole to $`r=M`$ for an extreme Kerr hole.

At all $`rr_{\mathrm{ISCO}}`$ there exists a family of circular orbits, each member having a different inclination angle $`\iota `$. We are interested in parameterizing these orbits as functions of $`r`$ and $`\iota `$. Consider first the weak-field limit, $`rr_{\mathrm{ISCO}}`$. Ryan  has provided formulas which, with some manipulation, give $`E`$, $`L_z`$, and $`Q`$ as functions of $`r`$ and $`\iota `$:

$`E`$ $`=`$ $`1{\displaystyle \frac{M}{2r}}2{\displaystyle \frac{a}{M}}\left({\displaystyle \frac{M}{r}}\right)^{3/2}\mathrm{cos}\iota ,`$ (9)
$`\left(L_z^2+Q\right)^{1/2}`$ $`=`$ $`\sqrt{rM}\left[13{\displaystyle \frac{a}{M}}\left({\displaystyle \frac{M}{r}}\right)^{3/2}\mathrm{cos}\iota \right].`$ (10)
\[Note there is a sign error in Eq. (7) of Ref. , as can be seen by taking the zero eccentricity limit of Eq. (6) of Ref. .\] One can then get $`L_z`$ from the definition of the inclination angle:

$$\mathrm{cos}\iota =\frac{L_z}{\sqrt{L_z^2+Q}}.$$
(11)
Two features of these formulas are particularly noteworthy. First, energy monotonically increases as $`\iota `$ increases: $`E/\iota >0`$ for all parameters. This turns out to be true everywhere, not just in the weak field. Second, $`L_z`$ monotically decreases as $`\iota `$ increases provided we don’t abuse the applicability of Eqs. (10) and (11): $`L_z/\iota <0`$ except when $`r/M\left(6a/M\right)^{2/3}`$. This range is not even close to the weak-field, so there is no reason to believe that this result is at all physically relevant. Nonetheless, it foreshadows the behavior of the nearly horizon skimming orbits.

Turn now from the weak field to the strong field. As is conventional , I will use Eq. (11) to define the inclination angleAs discussed in Ref. , this angle doesn’t necessarily accord with intuitive notions of inclination angle. For example, except when $`a=0`$, $`\iota `$ is not the angle at which most observers would see the small body cross the equatorial plane. even in this regime. The most-bound orbitIn this paper, the terms “most-bound”, “least-bound” and “marginally-bound” describe orbits at each radius. This contrasts with the usage in, e.g., Ref.  where these terms refer to properties of all orbits, regardless of radius. (the orbit with the smallest orbital energy) at each radius is the prograde, equatorial orbit. Its constants are given by

$`E^{\mathrm{mb}}`$ $`=`$ $`{\displaystyle \frac{12v^2+qv^3}{\sqrt{13v^2+2qv^3}}},`$ (12)
$`L_z^{\mathrm{mb}}`$ $`=`$ $`rv{\displaystyle \frac{12qv^3+q^2v^4}{\sqrt{13v^2+2qv^3}}},`$ (13)
$`Q^{\mathrm{mb}}`$ $`=`$ $`0,`$ (14)
where $`v\sqrt{M/r}`$ and $`qa/M`$. At fixed radius, the orbital energy increases as the tilt increases from the most-bound orbit at $`\iota =0^{}`$ to the least-bound orbit. The least-bound orbit’s characteristics depend upon $`r`$ and the black hole’s spin. If $`rr_{\mathrm{ret}}`$, where

$$r_{\mathrm{ret}}/M=3+Z_2+\left[(3Z_1)(3+Z_1+2Z_2)\right]^{1/2},$$
(15)
then the least-bound orbit is just the retrograde, equatorial orbit. This orbit has $`Q=0`$ and $`\iota =180^{}`$; expressions for its energy and angular momentum can be found in Ref. . For radii $`r_{\mathrm{ISCO}}rr_{\mathrm{ret}}`$, the least-bound orbit is the marginally-stable orbit: the orbit which satisfies $`R=0`$, $`R^{}=0`$, and $`R^{\prime \prime }=0`$. This orbit has the maximum allowed inclination angle $`\iota _{\mathrm{max}}`$ at that radius. Any orbit tilted at a larger angle is unstable to small perturbations and will quickly plunge into the black hole.

For the rest of this paper, I will focus on the extreme strong field of rapidly rotating black holes. The orbits of interest are well inside the radius of the retrograde orbit. Hence, the energy $`E^{\mathrm{lb}}`$, angular momentum $`L_z^{\mathrm{lb}}`$, and Carter constant $`Q^{\mathrm{lb}}`$ of the least-bound orbit will be determined by numerically solving the equations $`R=0`$, $`R^{}=0`$, $`R^{\prime \prime }=0`$.

To compute the properties of a circular orbit, pick two of its constants — e.g., the orbit’s radius $`r`$ and angular momentum $`L_z`$ — and solve $`R=0=R^{}`$ to find the other two. This yields the following solution for $`E(r,L_z)`$ and $`Q(r,L_z)`$ :

$`E(r,L_z)`$ $`={\displaystyle \frac{a^2L_z^2(rM)+r\mathrm{\Delta }^2}{aL_zM\left(r^2a^2\right)\pm \mathrm{\Delta }\sqrt{r^5(r3M)+a^4r(r+M)+a^2r^2(L_z^22Mr+2r^2)]}}},`$ (16)
$`Q(r,L_z)`$ $`={\displaystyle \frac{\left[(a^2+r^2)E(r,L_z)aL_z\right]^2}{\mathrm{\Delta }}}\left[r^2+a^2E(r,L_z)^22aE(r,L_z)L_z+L_z^2\right].`$ (17)
There is a sign choice in the denominator of Eq. (16). In general, only one choice is physically meaningful at a given value of $`r`$. The argument of the square root in the denominator of Eq. (16) goes to zero at some radius $`r_{\mathrm{branch}}(a)`$; the plus sign corresponds to $`rr_{\mathrm{branch}}(a)`$, and the minus sign to $`rr_{\mathrm{branch}}(a)`$. In Ref. , Eq. (16) was used with the plus sign only since the focus in that paper was on comparatively large radius \[in all cases, $`r_{\mathrm{branch}}(a)`$ is close to $`2M`$\]. In this work, since I will focus on the extreme strong field, both signs are needed.

Rather than using Eq. (16), I will avoid this sign ambiguity by starting with $`r`$ and $`E`$, and re-solving the system $`R=0=R^{}`$ for $`L_z`$ and $`Q`$. Then,

$`L_z(r,E)`$ $`={\displaystyle \frac{EM(r^2a^2)\mathrm{\Delta }\sqrt{r^2(E^21)+rM}}{a(rM)}},`$ (18)
$`Q(r,E)`$ $`={\displaystyle \frac{\left[(a^2+r^2)EaL_z(r,E)\right]^2}{\mathrm{\Delta }}}\left[r^2+a^2E2aEL_z(r,E)+L_z(r,E)^2\right].`$ (19)
(There exists a second solution for $`L_z`$ which has a $`+`$ sign in front of the $`\mathrm{\Delta }`$, but it is not physically meaningful.) Note that Eq. (18) does not behave well as $`a0`$. This is because of a degeneracy in this limit: knowledge of any three of the parameters $`E`$, $`r`$, $`L_z`$, and $`Q`$ suffices to determine the orbit (because of spherical symmetry). Since this paper deals with $`aM`$, this issue is irrelevant here.

Assembling strong-field orbits now reduces to a simple recipe. First, pick an orbital radius. Allow the orbital energy to vary from $`E^{\mathrm{mb}}`$ \[Eq. (13)\] to $`E^{\mathrm{lb}}`$ (found by solving the system $`R=0`$, $`R^{}=0`$, $`R^{\prime \prime }=0`$). For each energy, find $`L_z`$ and $`Q`$ with Eqs. (18) and (19). Parameterize each orbit by its inclination angle $`\iota `$ \[Eq. (11)\]. Repeat at a new radius.

## III Nearly horizon-skimming orbits

### A Overview

Consider for a moment the extreme Kerr limit, $`a=M`$. From Eq. (7), the ISCO is located at $`r=M`$, which is also the coordinate of the event horizon. It is not difficult to show that there exists a set of orbits at this radius, with the parameters

$`2M/\sqrt{3}`$ $``$ $`L_z\sqrt{2}M,`$ (20)
$`E`$ $`=`$ $`L_z/2M,`$ (21)
$`Q`$ $`=`$ $`3L_z^2/4M^2.`$ (22)
These are the horizon skimming orbits. They obey $`L_z/\iota >0`$, similar to the behavior seen when the weak-field results for $`L_z`$ are pushed into the very strong field. This is opposite to the typical behavior, as exemplified by the correct usage of Eqs. (10) and (11).

Focus for now on the most-bound and least-bound orbits. These orbits bound the behavior of all orbits at each radius. As described in Sec. II, it is simple to solve for $`L_z^{\mathrm{lb}}`$ as a function of radius, at least numerically. Doing so for $`a=M`$, we find that there is a region stretching to $`r1.8M`$ in which $`L_z^{\mathrm{lb}}L_z^{\mathrm{mb}}`$. The orbits in this domain have the same dynamical characteristics as Wilkins’ horizon skimming orbits, so we shall call them “nearly horizon skimming” (NHS) orbits.

Figure 1 illustrates the NHS region, and contrasts it with the “usual” behavior of Kerr orbits. It plots $`L_z^{\mathrm{mb}}`$ and $`L_z^{\mathrm{lb}}`$ for black holes with $`a=0.95M`$ (top panel) and $`a=M`$ (bottom panel). For $`a=0.95M`$, the least-bound and most-bound orbits coincide at the ISCO. Moving out in radius, the most-bound orbit’s $`L_z`$ grows and the least-bound orbit’s $`L_z`$ shrinks. This makes sense intuitively, since the inclination angle of the least-bound orbit grows as we move away from the ISCO. (Eventually, it tips over completely to $`\iota =180^{}`$, and becomes the retrograde equatorial orbit.) The lower panel of Fig. 1 shows the behavior when $`a=M`$. We see the horizon skimming orbits at $`r=M`$ and the NHS orbits stretching out to $`r1.8M`$. At that point, $`L_z^{\mathrm{lb}}`$ and $`L_z^{\mathrm{mb}}`$ cross over. Moving further out in radius, they behave in the “normal” way.

There must exist some critical spin value, $`0.95M<a_{\mathrm{NHS}}<M`$, at which NHS orbits first come into existence. The NHS orbits are bounded by two radii, $`r_{\mathrm{ISCO}}`$ and $`r_{\mathrm{cross}}`$ \[where the angular momentum of the least-bound and most-bound orbits cross: $`L_z^{\mathrm{mb}}(r_{\mathrm{cross}})=L_z^{\mathrm{lb}}(r_{\mathrm{cross}})`$\]. Decreasing $`a`$ must decrease the size of the NHS region — $`r_{\mathrm{ISCO}}`$ and $`r_{\mathrm{cross}}`$ approach one another. The spin at which these radii coincide and the NHS region vanishes is $`a_{\mathrm{NHS}}`$.

Recall from Eq. (16) that there are two solutions for orbital energy, only one of which is usually physical. At the crossover point, however, both of these energies must be physical: that with the minus sign (which is larger in magnitude) gives $`E^{\mathrm{lb}}`$, and vice versa. When $`r_{\mathrm{ISCO}}=r_{\mathrm{cross}}`$ the two energies likewise coincide, since the most-bound orbit is also the least-bound orbit at the ISCO. Thus, $`a_{\mathrm{NHS}}`$ is the spin that satisfies

$$E_+[r_{\mathrm{ISCO}}(a_{\mathrm{NHS}}),L_z^{\mathrm{mb}},a_{\mathrm{NHS}}]=E_{}[r_{\mathrm{ISCO}}(a_{\mathrm{NHS}}),L_z^{\mathrm{mb}},a_{\mathrm{NHS}}]$$
(23)
\[where $`E_\pm `$ denotes the two roots given in Eq. (16)\]. The solution to this is $`a_{\mathrm{NHS}}/M=0.952412\mathrm{}`$.

Figure 2 illustrates the change to the NHS region as spin is varied, vanishing altogether when $`a_{\mathrm{NHS}}`$ is reached. The size is still significant near $`a=0.998M`$. This is interesting, since $`a=0.998M`$ is probably the largest spin value that we can encounter in nature . This opens the possibility that NHS orbits may play a role in astrophysical processes.

### B Nearly horizon skimming orbits as seen by a ZAMO

The NHS characteristic $`L_z/\iota >0`$ can be better understood by considering these orbits from the viewpoint of the zero angular momentum observer, or ZAMO . The ZAMO is the observer that corotates with the coordinate system, such that its angular velocity as seen at infinity is

$$\omega _{\mathrm{ZAMO}}=\frac{2aMr}{(r^2+a^2)^2a^2\mathrm{\Delta }\mathrm{sin}\theta ^2}.$$
(24)
If one imagines spacetime to be dragged into a whirlpool-like flow by the black hole’s rotation, then the ZAMO is the observer who simply rides along with the flow. (In accordance with this viewpoint, Wald  calls the ZAMO the “locally non-rotating observer”.)

The ZAMO examines orbits in its local neighborhood. This allows it to interpret the motion of a body in a strong-field orbit with special relativistic formulas. For instance, the ZAMO measures the small body to have velocity $`\stackrel{}{v}`$ and energy $`E_{\mathrm{local}}=(1\stackrel{}{v}\stackrel{}{v})^{1/2}\gamma `$. This local energy is not a constant of the body’s motion. It is related to the conserved energy by the formula

$$E=\alpha E_{\mathrm{local}}+\omega _{\mathrm{ZAMO}}L_z,$$
(25)
where

$$\alpha =\sqrt{\frac{\mathrm{\Sigma }\mathrm{\Delta }}{(r^2+a^2)^2a^2\mathrm{\Delta }\mathrm{sin}\theta ^2}}$$
(26)
is the lapse function. Knowing the energy and the angular momentum of the small body then tells us the body’s speed $`v=\sqrt{\stackrel{}{v}\stackrel{}{v}}`$ as seen by the ZAMO.

Applying Eq. (25) shows that a body’s speed is smallest in the most-bound orbit and highest in the least-bound orbit, varying smoothly between the two. This is due to the dragging of spacetime by the hole’s spin. Let us contrast the prograde and retrograde orbits, $`\iota =0^{}`$ and $`\iota =180^{}`$. The prograde orbit is moving “downstream”: part of the motion needed to keep it in orbit is provided by the dragging of inertial frames. It can orbit with relatively small velocity. The retrograde orbit, by contrast, must “swim upstream”: it must overcome the dragging of inertial frames on top of the motion needed to stay in orbit. It therefore is more energetic than the prograde orbit. The angle $`\iota `$ smoothly varies the orbit between these extremes, so that larger $`\iota `$ corresponds to larger energy (and larger speed).

The velocity of a body in a non-equatorial orbit has components in both the $`\theta `$ and $`\varphi `$ directions. The $`\theta `$ component goes to zero, however, at the orbit’s turning points $`(\theta _{\mathrm{max}},\theta _{\mathrm{min}})`$, when it reverses in $`\theta `$. (Formulas for computing $`\theta _{\mathrm{max}/\mathrm{min}}`$ can be found in Ref. ; they are just the angles at which $`d\theta /d\tau =0`$ \[cf. Eq. (3)\].) At these two points, the velocity is purely along $`\varphi `$, and the body’s motion is fully described by the component $`v_{\widehat{\varphi }}=\stackrel{}{v}\stackrel{}{e}_{\widehat{\varphi }}`$ (where $`\stackrel{}{e}_{\widehat{\varphi }}`$ is the $`\varphi `$-component of the orthonormal basis that the ZAMO uses to make measurements). Evaluating Eq. (25) at $`\theta _{\mathrm{max}/\mathrm{min}}`$ and writing $`E_{\mathrm{local}}=\gamma `$ gives a condition for $`v_{\widehat{\varphi }}`$ at the turning points.

The top panel of Fig. 3 shows $`v_{\widehat{\varphi }}(\theta _{\mathrm{max}/\mathrm{min}})`$ for the most-bound and least-bound orbits at several interesting spins. As expected, bodies move quite a bit faster in the least-bound orbit than in the most-bound orbit. Perhaps more interestingly, $`v_{\widehat{\varphi }}`$ becomes substantially smaller towards the innermost orbits. This is because of the “freezing” of physics near the hole’s event horizon; see Ref.  (particularly Secs. IIC1 and IIIA4) for further discussion. Close to the horizon, a body’s motion locks to the hole’s spin, and it is dragged into rigid corotation. This horizon locking causes $`v_{\widehat{\varphi }}(r)`$ to peak at $`r1.5M`$: the “Keplerian” tendency of a body to move faster as it moves inward dominates at large radii (asymptoting to $`\sqrt{M/r}`$ in the weak field), but is overwhelmed as the body locks onto the horizon close to the hole.

Consider next the orbit’s angular momentum $`L_z`$. The ZAMO sees $`L_z`$ as the product of a radius of gyration $`\varpi `$ and a locally measured azimuthal momentum $`p_{\widehat{\varphi }}`$ :

$$L_z=\varpi p_{\widehat{\varphi }}=\varpi \gamma v_{\widehat{\varphi }}.$$
(27)
The $`\gamma `$ factor causes $`p_{\widehat{\varphi }}`$ to be even more strongly peaked than $`v_{\widehat{\varphi }}`$; see the lower panel of Fig. 3.

The radius of gyration is a purely geometric quantity. It is just the circumference of the ZAMO’s constant $`r`$, constant $`\theta `$ orbit, divided by $`2\pi `$:

$`\varpi `$ $`=`$ $`𝒞(\theta _{\mathrm{max}/\mathrm{min}})/2\pi ={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}\sqrt{g_{\varphi \varphi }(r,\theta _{\mathrm{max}/\mathrm{min}})}𝑑\varphi `$ (28)
$`=`$ $`\sqrt{{\displaystyle \frac{(r^2+a^2)^2a^2\mathrm{\Delta }\mathrm{sin}\theta _{\mathrm{max}/\mathrm{min}}^2}{r^2+a^2\mathrm{cos}\theta _{\mathrm{max}/\mathrm{min}}^2}}}\mathrm{sin}\theta _{\mathrm{max}/\mathrm{min}}.`$ (29)
This function is plotted in Fig. 4. There are no surprises here. In the most-bound orbit, $`\varpi `$ increases monotically with orbital radius, asymptoting to $`r`$ at large radius. The least-bound orbit is more interesting: $`\varpi `$ is smaller (not surprising, since it is a tilted orbit) and is nearly flat as a function of radius, at least over the range of NHS orbits. This near flatness is due to the orbit’s increasing tilt: $`\theta _{\mathrm{max}/\mathrm{min}}`$ of the least-bound orbit changes such that the circumference of the ZAMO’s orbit at $`\theta _{\mathrm{max}/\mathrm{min}}`$ remains nearly constant.

Since $`L_z`$ is just the product of the curves shown in Figs. 3 and 4, any unusual features in the behavior of $`L_z`$ must arise from features in $`p_{\widehat{\varphi }}`$ and $`\varpi `$. Considering these two figures, we immediately see why NHS orbits have $`L_z^{\mathrm{lb}}>L_z^{\mathrm{mb}}`$: the least-bound orbit has so much more linear momentum than the most-bound orbit that it compensates for its smaller radius of gyration. The linear momentum is so much larger, in turn, because the most-bound orbit is strongly locked to the spin of the black hole. The least-bound orbit is also locked for $`r`$ very close to $`r_{\mathrm{ISCO}}`$. However, it does not lock as strongly: being so energetic, the least-bound orbit only locks as the very innermost orbits are approached. Hence, NHS orbits exist because very close orbits are forced to move in rigid corotation with the event horizon.

## IV Application: evolution under gravitational-wave emission

As discussed in the Introduction, binary systems consisting of a small body spiraling into a massive black hole are one of the more anticipated sources of gravitational waves for space-based detectors such as LISA. The NHS region is still rather large for $`a0.998M`$ (cf. Fig. 2), indicating that there might be plenty of time for the properties of NHS orbits to influence the gravitational-wave signal of these sources.

In the extreme mass ratio, radiation reaction should operate adiabatically. In other words, the timescale for gravitational-wave emission to change an orbiting body’s parameters ($`r`$, $`E`$, $`L_z`$, $`Q`$) should be significantly longer than an orbital period: $`T_{\mathrm{orb}}/\tau _{\mathrm{GW}}1`$. (More careful discussion of this point can be found in Ref. .) Because the change in these constants is very slow, the body’s motion is well-approximated as geodesic over small time intervals. It is thus useful to regard the body’s true, radiatively evolving trajectory as motion through a sequence of geodesic orbits. For circular orbits, we regard the body’s inspiral trajectory as the evolution of its radius and inclination angle.

A simple argument shows that a body spiraling through the NHS region should behave rather differently from a body spiraling through “normal” orbits. At a given moment, a circular orbit can be represented as a point on the $`(r,L_z)`$ plane. For example, a point on either the top or the bottom panels of Fig. 1 lying between the most-bound and least-bound orbit curves represents a physically allowed orbit. Radiation reaction will tend to drive this point downward and to the left in this figure: gravitational-wave emission shrinks an object’s radius and reduces $`L_z`$. For a “normal” orbit, as in the top panel of Fig. 1, the body is pushed toward less-bound orbits, increasing $`\iota `$. This behavior was predicted by weak-field calculations  and was recently confirmed in the strong field . The evolution is in the opposite direction inside the NHS region — radiation emission pushes the body toward the sequence of most-bound orbits, so that $`\iota `$ decreases.

Decreasing inclination angle could have important consequences for the gravitational-wave signal. Inclined circular orbits emit waves characterized by harmonics of $`\mathrm{\Omega }_\varphi =2\pi /T_\varphi `$ and $`\mathrm{\Omega }_\theta =2\pi /T_\theta `$ (where $`T_\varphi `$ is the period of a full range of $`\varphi `$ motion, and $`T_\theta `$ is the period of $`\theta `$ motion) . Roughly speaking, the $`\theta `$ motion modulates the gravitational waveform. The waves emitted by a body in an equatorial orbit, on the other hand, depend only on $`\mathrm{\Omega }_\varphi `$ — there is no motion in $`\theta `$, so the amplitudes of $`\mathrm{\Omega }_\theta `$ harmonics vanish. As the orbiting body moves toward the equatorial plane, the importance of the $`\mathrm{\Omega }_\theta `$ harmonics decreases, simplifying the waveform. Observing such simplification would be a clear signal that the small body is evolving toward the equatorial plane.

A rigorous analysis of NHS evolution requires a radiation reaction formalism good deep in the strong field. Such a formalism, specialized to the adiabatic evolution of circular orbits, is given in Ref. . The code described there solves the Teukolsky equation  to calculate the flux of $`E`$ and $`L_z`$ to infinity and down the event horizon. It then applies the proof that circular orbits remain circular (changing only their radius and inclination angle) when they evolve adiabatically under radiation emission  to compute the change in $`Q`$. From this information, it is simple to obtain the change in $`r`$ and $`\iota `$. By computing $`dr/dt`$ and $`d\iota /dt`$ at a large number of points, it is not difficult to compute the inspiral trajectory followed by a body that starts in an inclined, circular orbit. An algorithm and code for computing such trajectories will be presented at a later date .

Figure 5 shows some trajectories in the NHS region for the inspiral of a $`1M_{}`$ compact body into a $`10^6M_{}`$ Kerr black hole. The hole has spin $`a=0.998M`$. The inspiral time for most trajectories is significant. The number of orbits around the spin axis of the black hole by the body varies from $`N_\varphi 73000`$ for the inspiral beginning at $`\iota =0^{}`$ to $`N_\varphi 400`$ for the inspiral beginning at $`\iota =40^{}`$. The number of gravitational-wave cycles is roughly proportional to $`N_\varphi `$, indicating that waves from the NHS region could contribute substantially to a measured signal. The inspiral duration scales with $`M^2/\mu `$, and the number of orbits with $`M/\mu `$. These quantities should be of interesting magnitude for a reasonably wide range of masses, at least for spin $`a=0.998M`$. (As the spin of the black hole is dialed down to $`a=a_{\mathrm{NHS}}`$, the NHS orbits lose significance and inspiral through them becomes irrelevant.) Gravitational waves produced in the NHS region will form a significant part of the data measured from inspirals into rapidly rotating black holes.

Figure 5 also indicates that the signature of NHS inspiral is not very large. Although $`\iota `$ does decrease, as expected, the degree of decrease is quite small. For the inspirals shown, the largest change in $`\iota `$ is for the trajectory that starts at $`\iota =25^{}`$: the final inclination is $`\iota =23.8^{}`$. The change $`\delta \iota =1.2^{}`$ is a rather paltry effect, not at all a robust signature of the strong field.

It seems that the unique signature of orbital evolution through the NHS region is not likely to be of observational interest. Before dismissing the idea entirely, let us speculate for a moment on how eccentricity might affect this conclusion. Using a weak-field radiation reaction force, Ryan has calculated the leading order evolution of the parameters of an inclined, eccentric orbit of a Kerr black hole due to gravitational-wave emission . His results indicate that the eccentricity $`e`$ can significantly magnify the rate at which $`\iota `$ changes:

$$\frac{d\iota }{dt}|_{e0}=\left(1e^2\right)^{7/2}\left(1+\frac{73}{24}e^2+\frac{37}{96}e^4\right)\frac{d\iota }{dt}|_{e=0}.$$
(30)
This relation is derived under the assumption that $`rM`$, so there is no reason to believe it holds in the NHS region. However, it could be indicative of the effect that eccentricity has in the strong field. Notice that as eccentricity varies over the range $`0e0.8`$ (a reasonable range for sources that LISA is likely to measure ), the factor magnifying $`d\iota /dt`$ increases from $`1`$ to about $`100`$. One can speculate that $`d\iota /dt`$ will be magnified by a similar factor (at least to order of magnitude) under strong-field radiation reaction. At this time, there is no way to tell. Hopefully, strong-field radiation reaction programs  will be able to model the evolution of generic Kerr orbits soon. If this speculation proves correct, then the gravitational-waves emitted by inspiral through the NHS region will contain a very strong signature of the Kerr metric strong field.

## V Conclusion

The result of this analysis is somewhat equivocal. I have shown that the “horizon skimming” orbits found by Wilkins are a subgroup of a larger family of orbits, the “nearly horizon skimming” orbits existing around any black hole with spin $`a0.952412M`$. The signature characteristic of these NHS orbits is that, at fixed orbital radius, the $`z`$-component of orbital angular momentum increases as the orbit’s tilt increases, in opposition to the “normal” behavior. This has the consequence that radiation emission, which carries angular momentum away from the orbit, tends to drive the system into an equatorial orbit — the system’s inclination angle decreases rather than increases.

The fact that radiation emission drives these orbits toward the equator is a consequence of the motion of bodies in the strong field of Kerr black holes. As a body spirals in from large radius, its orbital speed initially increases, in accordance with Keplerian intuition that a body in a circular orbit of radius $`r`$ has a speed $`v=\sqrt{M/r}`$. As the body gets very close to the black hole, however, its dynamics become dominated by the nearness of the hole’s event horizon. All physical processes become “locked” to the hole as the horizon is approached. Eventually the motion of orbiting bodies synchronizes with the hole’s spin. The least-bound orbit does not lock as quickly as the the most-bound orbit: close to the horizon, a less-bound orbit moves quite a bit faster than a highly bound orbit, not locking to the hole’s spin until the very innermost orbits are reached. This effect is shown in Fig. 3.

As an inspiraling body’s inclination decreases, the modulation of its gravitational waveform by the $`\theta `$ motion decreases. Measuring this decreasing modulation could be an observational hallmark of the Kerr strong field. Unfortunately, at the moment the most mature formalism for computing the effects of gravitational radiation reaction on the orbits of small bodies can only tackle the evolution of circular, inclined orbits. In this restricted case, the total change in inclination angle as a body evolves through the NHS region is quite small — at most, the orbital inclination changes by about $`12`$ degrees. Although I have not examined the effect of this change on the waveform, it is hard to imagine it will be particularly marked.

It is extremely unlikely that realistic extreme mass ratio inspirals will be circular. When such binaries initially form, their eccentricities are likely to be rather close to $`1`$ , and remain significant by the time that the binary enters the sensitivity band of LISA. Weak field calculations suggest that the rate by which the inclination angle changes might be magnified by a rather large factor. If we take the leading order results at face value \[cf. Eq. (30)\], the change might be so strong that the orbit is driven nearly into the equatorial plane.

Because the NHS region is so deep in the strong field, however, it is a mistake to infer too much from leading order effects. Further progress will require a scheme to calculate radiation reaction in the strong field for generic Kerr orbits — orbits that are inclined and eccentric. This may require local radiation reaction forces , or a clever scheme for extracting the change in the Carter constant from the radiation flux . The fact that observations might be able to detect a strong signature of the Kerr strong field will hopefully motivate future progress.

###### Acknowledgements.

This work grew out of a more general program to compute waveforms and radiation reaction sequences from extreme mass ratio circular orbit inspirals. That work, in turn, has benefitted greatly from conversations with Lior Burko, Curt Cutler, Dan Kennefick, Sam Finn, Lee Lindblom, Sterl Phinney, and Kip Thorne. The package Mathematica was used to aid some of the calculations; all plots were produced using the package sm. This research was supported by NSF Grant AST-9731698 and NASA Grants NAG5-7034 and NAGW-4268 at Caltech, and NSF Grant PHY-9907949 at the ITP.