text:
# A survey of spinning test particle orbits in Kerr spacetime

## I INTRODUCTION

The present paper continues the investigation initiated in Hartl (2003), which considered the dynamics of spinning test particles (i.e., compact astrophysical objects) orbiting a spinning black hole (Kerr spacetime). The primary objective is to determine whether the orbits of spinning compact objects spiraling into much more massive black holes are chaotic or not. Evidence for chaotic orbits in relativistic binaries has been presented in Suzuki and Maeda (1997); Levin (2000, ); Cornish (2001); Cornish and Levin (2002, ) (although only Suzuki and Maeda (1997) was directly concerned with the extreme mass ratio systems we consider here). An astrophysical example of the systems we consider is a solar-mass black hole or neutron star orbiting a supermassive black hole in a galactic nucleus. Such systems are potentially important sources of gravitational radiation in frequency bands detectable by space-based gravitational wave detectors such as the proposed Laser Interferometer Space Antenna (LISA) mission. The methods of data analysis for signals from such detectors typically rely on matched filters, which use a discrete mesh in parameter space to achieve a given signal-to-noise ratio. The presence of chaos would cause the number of templates needed to fit a given signal to grow exponentially with decreasing mesh size, potentially overwhelming the computers performing the analysis.

Many studies have used the Papapetrou equations to investigate the dynamics of spinning test particles in the background spacetime of a central black hole (including most recently Tanaka et al. (1996); Suzuki and Maeda (1997); Hartl (2003); Semerák (1999); Suzuki and Maeda (2000); see Semerák (1999) and Suzuki and Maeda (2000) for a more thorough literature review). We found in Hartl (2003) that the Papapetrou equations formally admit chaotic solutions in Kerr spacetime, but the chaos seemed to disappear for physically realistic spins. This conclusion was tentative, since we investigated only a few values of the many parameters appearing in the equations. In the present study, we undertake a thorough search of parameter space in order to determine the prevalence of chaotic orbits, both for dynamically interesting (but physically unrealistic) orbits with high spin parameter $`S`$ and for astrophysically relevant systems with smaller spin.

As in Hartl (2003), we use Lyapunov exponents to measure the divergence of nearby trajectories and to provide an estimate for the timescale of the divergence. For chaotic systems, initial conditions separated by a small distance $`ϵ_0`$ diverge exponentially, with the separation growing as $`ϵ(\tau )=ϵ_0e^{\lambda \tau }`$, where $`\lambda `$ is the Lyapunov exponent. For nonchaotic orbits, $`\lambda =0`$, but when $`\lambda `$ is nonzero it provides the $`e`$-folding timescale $`\tau _\lambda =1/\lambda `$ for chaotic behavior to manifest itself. In this paper we use two rigorous methods for determining the maximum Lyapunov exponent $`\lambda _{\mathrm{max}}`$ described in Hartl (2003). Further details appear below in Sec. IV and in Sec. III of Hartl (2003).

Our parameter space search makes use of a convenient orbital parameterization technique, which allows us to specify desired values of pericenter $`r_p`$ and orbital inclination $`\iota `$ (Sec. III). We vary these orbital parameters along with spin magnitude, eccentricity, Kerr spin angular momentum $`a`$, and spin inclination, and calculate the largest Lyapunov exponent for each set of initial conditions. A complete description of our methods for searching parameter space appears in Sec. V. Though not exhaustive, the resulting survey of parameter space gives a thorough picture of chaos in the Papapetrou equations with a Kerr background.

We set $`G=c=1`$, use sign conventions as in MTW (Misner et al., 1973), and use Boyer-Lindquist coordinates $`(r,\theta ,\varphi )`$ for Kerr spacetime. We use vector arrows for 4-vectors and boldface for Euclidean vectors. The symbol $`\mathrm{log}`$ denotes the natural logarithm.

## II Spinning test particles

We use the Papapetrou equations of motion (Papapetrou, 1951) to model a spinning test particle. As reformulated by Dixon (Dixon, 1970), these equations describe the motion of a “pole-dipole particle,” which neglects effects smaller than those due to the mass monopole and spin dipole (thus neglecting the tidal coupling, which is a mass quadrupole effect). The Papapetrou equations describe a particle with negligible mass compared to the masses responsible for generating the background spacetime. In our case, we write $`\mu `$ for the test particle’s mass and $`M`$ for the mass of the central Kerr black hole; the test particle limit then requires that $`\mu M`$.

### II.1 Equations of motion

Dixon’s formulation for the equations of motion uses an antisymmetric spin tensor $`S^{\mu \nu }`$ to represent the spin angular momentum of the particle. The covariant derivative of the 4-momentum $`p^\mu `$ deviates from the nonspinning (geodesic) case by a term representing a coupling of the spin to the background spacetime curvature. If we write $`v^\mu `$ for the 4-velocity, the full equations of motion appear as follows:

$`{\displaystyle \frac{dx^\mu }{d\tau }}`$ $`=`$ $`v^\mu `$
$`_\stackrel{}{v}p^\mu `$ $`=`$ $`\frac{1}{2}R_{\nu \alpha \beta }^\mu v^\nu S^{\alpha \beta }`$ (1)
$`_\stackrel{}{v}S^{\mu \nu }`$ $`=`$ $`2p^{[\mu }v^{\nu ]}.`$
Here $`R_{\nu \alpha \beta }^\mu `$ is the Riemann curvature tensor of the spacetime, which in our case corresponds to the Riemann tensor for the Kerr metric. Note that in the case of vanishing spin (all components of $`S^{\mu \nu }`$ equal to zero) or flat spacetime (all components of $`R_{\nu \alpha \beta }^\mu `$ equal to zero), we recover the geodesic equation $`_\stackrel{}{v}p^\mu =0`$.

For numerical and conceptual purposes, we use a reformulation of the equations of motion in terms of the momentum 1-form $`p_\mu `$ and spin 1-form $`S_\mu `$. The spin 1-form can be derived from the 4-momentum and spin tensor using

$$S_\mu =\frac{1}{2}ϵ_{\mu \nu \alpha \beta }u^\nu S^{\alpha \beta },$$
(2)
where $`u^\nu =p^\nu /\mu `$. The spin 1-form is easier to visualize than the spin tensor, and the Papapetrou equations are better-behaved numerically in the $`S0`$ limit when expressed in terms of the spin 1-form. Details of this formulation appear in Hartl (2003) and the appendix.

### II.2 General constraints

The Papapetrou system of equations (II.1) is highly constrained, even in arbitrary spacetimes. The 4-momentum and spin satisfy normalization conditions that are conserved by the equations of motion:

$$p^\mu p_\mu =\mu ^2$$
(3)
and

$$\frac{1}{2}S^{\mu \nu }S_{\mu \nu }=S^2.$$
(4)
A further condition is required to identify a unique worldline for the center of mass:

$$p_\mu S^{\mu \nu }=0.$$
(5)
A more detailed discussion of this “spin supplementary condition” appears in Hartl (2003).

### II.3 Kerr constraints

The Papapetrou equations share an important property with geodesic motion, namely, to each symmetry of the background spacetime (typically represented by a Killing vector) there corresponds a constant of the motion. Kerr spacetime has two such symmetries, which provide two constraints in addition to those described in the previous section. (The Killing tensor that leads to the Carter constant in the geodesic does not correspond to a conserved quantity when the spin is nonzero; see Sec. III.2.3 below.)

For each Killing vector $`\stackrel{}{\xi }`$ there corresponds a constant given by the standard expression for geodesics plus a contribution due to a coupling with the spin tensor:

$$C_\xi =\xi ^\mu p_\mu \frac{1}{2}\xi _{\mu ;\nu }S^{\mu \nu }.$$
(6)
Kerr spacetime’s two Killing vectors $`\stackrel{}{\xi }^t=/t`$ and $`\stackrel{}{\xi }^\varphi =/\varphi `$ then give energy and $`z`$ angular momentum conservation:

$$E=p_t\frac{1}{2}g_{t\mu ,\nu }S^{\mu \nu }$$
(7)
and

$$J_z=p_\varphi \frac{1}{2}g_{\varphi \mu ,\nu }S^{\mu \nu }.$$
(8)
### II.4 The spin parameter $`S`$

The spin magnitude $`S`$ that appears in Eq. (4) quantifies the size of the spin and thus plays a crucial role in determining the behavior of spinning test particle systems. As in Hartl (2003), we measure distances and times in terms of $`M`$ and momenta in terms of $`\mu `$. In these units, $`S`$ is measured in units of $`\mu M`$. When measured in traditional geometrized units, the spin of a black hole can be as large as its mass squared, i.e., $`S_{\mathrm{geom},\mathrm{max}}=\mu ^2`$. In this case the spin parameter goes like $`S\mu ^2/(\mu M)=\mu /M`$, which is necessarily small for a test particle system. This result applies to all astrophysically relevant systems, as shown in Hartl (2003). Thus, we have the important physical constraint that *the spin parameter $`S`$ must satisfy $`S1`$ for all physically realistic systems*. A thorough discussion of various possibilities (including neutron stars and white dwarfs) in Hartl (2003) showed that realistic values of $`S`$ for LISA sources fall in the range $`10^4`$$`10^7`$.

The smallness of the spin’s effect is not dependent on our choice to measure $`S`$ in terms of $`\mu M`$. If instead we measure $`S`$ in terms of $`\mu ^2`$, then the equation of motion for $`S^{\mu \nu }`$ becomes (rewriting $`_\stackrel{}{v}`$ as $`D/d\tau `$)

$$\frac{D(S^{\mu \nu }/\mu ^2)}{d(\tau /M)}=2\frac{1}{\mu }p^{[\mu }v^{\nu ]},$$
(9)
where we have factored out from each variable its corresponding scale factor. This gives

$$\frac{DS^{\mu \nu }}{d\tau }=\frac{\mu }{M}\left(2p^{[\mu }v^{\nu ]}\right),$$
(10)
so that in these units the spin’s effective magnitude is suppressed by a factor of $`\mu /M`$. Measuring $`S`$ in terms of $`\mu M`$ effectively absorbs the unavoidable smallness of the spin effect into the spin parameter itself.

## III Parameterizing initial conditions

The many constraints on the equations of motion lead to considerable complexity in parameterizing the initial conditions. We summarize here the primary parameterization method described in Hartl (2003), and then discuss in detail a method for parameterizing initial conditions using the geometrical properties of the orbit.

### III.1 Energy and angular momentum parameterization

It is easy to parameterize the Papapetrou equations directly in terms of the momentum and spin components, but more physically relevant parameterizations are more difficult. The parameterization method discussed in Hartl (2003) solves for the initial conditions using the integrals of the motion $`E`$ and $`J_z`$. This method also serves as the foundation for a more sophisticated parameterization in terms of orbital parameters (Sec. III.2).

The energy and angular momentum parameterization method proceeds as follows. If we think of the system in terms of the spin 1-form $`S_\mu `$, we have twelve variables, four each for position, 4-momentum, and spin. Since Kerr spacetime is static and axially symmetric, without loss of generality we may set the initial time and axial angle to zero: $`\tau =\varphi =0`$. (We use the proper time $`\tau `$ as our time parameter, as discussed in Sec. IV below.) We then provide the initial values for $`r`$, $`\theta `$, and $`p_r`$, together with the constants $`E`$, $`J_z`$, and $`S`$. Finally, we provide two components of the spin vector in an orthonormal basis. It is easiest to specify the radial and $`\theta `$ components $`S^{\widehat{r}}`$ and $`S^{\widehat{\theta }}`$; since the $`r`$-$`\theta `$ part of the Kerr metric is diagonal, we may then easily derive the 1-form components $`S_r`$ and $`S_\theta `$.

Having specified values for seven of the variables, we must now use the five equations that relate them. Since we measure the particle’s momentum in terms of its rest mass $`\mu `$, the momentum normalization relation is

$$p^\mu p_\mu =1,$$
(11)
which we use to eliminate $`p_\theta `$. We then solve the spin normalization condition

$$S^\mu S_\mu =S^2$$
(12)
to eliminate $`S_\varphi `$. The spin-momentum orthogonality constraint and the two integrals of the motion then give three equations in the three remaining unknowns $`p_t`$, $`p_\varphi `$, and $`S_t`$:

$`0`$ $`=`$ $`p_\mu S_\nu g^{\mu \nu }`$ (13)
$`E`$ $`=`$ $`p_t+\frac{1}{2}g_{t\mu ,\nu }S^{\mu \nu }`$ (14)
$`J_z`$ $`=`$ $`p_\varphi \frac{1}{2}g_{\varphi \mu ,\nu }S^{\mu \nu }`$ (15)
(The first of these equations is equivalent to the condition $`p_\mu S^{\mu \nu }=0`$, as discussed in the appendix.)

We solve Eqs. (13)–(15) with a Newton-Raphson root finder, which works robustly given good initial guesses. The terms $`\frac{1}{2}g_{t\mu ,\nu }S^{\mu \nu }`$ and $`\frac{1}{2}g_{\varphi \mu ,\nu }S^{\mu \nu }`$ are typically small, even in the physically unrealistic case of $`S1`$, so $`E`$ and $`J_z`$ are good initial guesses for their corresponding momenta. We typically use $`S_\theta `$ as the initial guess for $`S_t`$. In all cases, we verify *a posteriori* that all constraints are satisfied to within machine precision.

### III.2 Orbital geometry parameterization

While the method described above is sufficient for determining a set of valid initial conditions, parameterizing orbits by energy and angular momentum is not particularly intuitive. It is much more natural to think in terms of the orbital geometry, so we prefer to parameterize by the pericenter $`r_p`$, the eccentricity $`e`$, and the orbital inclination angle $`\iota `$. Such parameters only have precise meaning for geodesic orbits, but are nevertheless still useful for the case of spinning test particles.

#### III.2.1 Kerr geodesics

The first step in parameterizing solutions of the Papapetrou-Dixon equations using orbital parameters is to solve the geodesic case. The traditional method for specifying a geodesic in terms of conserved quantities uses the energy $`E`$, the $`z`$ angular momentum $`L_z`$, and the Carter constant $`Q`$ Misner et al. (1973). In order to use the orbital parameters, we adopt a mapping from $`(r_p,e,\iota )`$ to $`(E,L_z,Q)`$ based on unpublished notes supplied by Teviet Creighton and Scott Hughes (and implemented in Hughes’s Kerr geodesic integrator Hughes ).

In order to use the more intuitive orbital parameters, we must determine the set $`(E,L_z,Q)`$ given the set $`(r_p,e,\iota )`$. We obtain two of the necessary equations by noting that the radial velocity $`dr/d\tau `$ vanishes at pericenter and apocenter, since these radii correspond to turning points in the radial motion. The equation for the time-evolution of the Boyer-Lindquist radius $`r`$ is Misner et al. (1973):

$$\mathrm{\Sigma }^2\left(\frac{dr}{d\tau }\right)^2=R(r),$$
(16)
where

$$R(r)=[E(r^2+a^2)aL_z]^2\mathrm{\Delta }[r^2+(L_zaE)^2+Q],$$
(17)
and we use the standard auxiliary variables

$$\mathrm{\Sigma }=r^2+a^2\mathrm{cos}^2\theta $$
(18)
and

$$\mathrm{\Delta }=r^22Mr+a^2.$$
(19)
The quantity $`a`$ is the Kerr spin parameter $`J/M`$, i.e., the central black hole’s spin angular momentum per unit mass, which is dimensionless in our normalized units. From Eq. (17) we see that $`dr/d\tau =0`$ implies that $`R(r)=0`$, so we obtain one equation at each turning point:111In this paper we never consider exactly circular orbits, but we note that our prescription fails in this case: the conditions (20) and (21) are identical when $`e=0`$, since $`r_a=r_p`$. For exactly circular orbits one must use the additional condition $`R^{}(r_p)=0`$, where $`R^{}=dR/dr`$.

$$R(r_p)=0$$
(20)
and

$$R(r_a)=0,$$
(21)
where the apocenter is defined by

$$r_a=\left(\frac{1+e}{1e}\right)r_p.$$
(22)
The final equation required to complete the mapping is Ryan (1996)

$$Q=L_z^2\mathrm{tan}^2\iota .$$
(23)
The value of $`\iota `$ resulting from this definition agrees closely with the maximum value of $`|\pi /2\theta |`$ for a numerically determined solution to the equations of motion, i.e., it faithfully captures a geometric property of the orbit.

Eqs. (20)–(23) give three equations in three unknowns, which are easy to solve using a nonlinear root finder as long as good initial guesses for the energy, angular momentum, and Carter constant can be found. The approach we adopt uses the degenerate cases of circular equatorial orbits to provide the raw material for analytical guesses. The energies for prograde and retrograde circular orbits in the equatorial plane are

$$E^{\mathrm{pro}}(r)=\frac{12v^2+\stackrel{~}{a}v^3}{\sqrt{13v^2+2\stackrel{~}{a}v^3}}$$
(24)
and

$$E^{\mathrm{ret}}(r)=\frac{12v^2\stackrel{~}{a}v^3}{\sqrt{13v^22\stackrel{~}{a}v^3}},$$
(25)
where we write $`v\sqrt{M/r}`$ and $`\stackrel{~}{a}=a/M`$ for notational simplicity. The initial guess for the energy is then an average of these energies, weighted using the inclination angle,222We adopt the convention that $`a0`$, so that $`\iota `$ indicates whether the orbit is prograde ($`0\iota <\frac{\pi }{2}`$) or retrograde ($`\frac{\pi }{2}<\iota \pi `$). with “radius” given by the semimajor axis of an ellipse with pericenter $`r_p`$ and eccentricity $`e`$:

$$E^{\mathrm{guess}}=\frac{1}{2}\left[\alpha _+E^{\mathrm{pro}}(r_{\mathrm{semi}})+\alpha _{}E^{\mathrm{ret}}(r_{\mathrm{semi}})\right],$$
(26)
where

$$\alpha _\pm =1\pm \mathrm{cos}\iota ,$$
(27)
and

$$r_{\mathrm{semi}}=\frac{r_p}{1e}.$$
(28)
It is possible (though rare) for Eq. (26) to yield an energy guess greater than 1; in this case, we simply set $`E^{\mathrm{guess}}=1`$.

Once we have a guess for the energy, we can find the corresponding guess for the angular momentum. Using the value from Eq. (26) and the expression for the angular momentum for a circular equatorial orbit gives

$$L_z^{\mathrm{guess}}=\mathrm{cos}\iota \sqrt{\frac{1e^2}{2(1E^{\mathrm{guess}})}}$$
(29)
as an initial guess for the angular momentum. Finally, the guess for the Carter constant is

$$Q^{\mathrm{guess}}=\left(L_z^{\mathrm{guess}}\right)^2\mathrm{tan}^2\iota .$$
(30)
Plugging Eqs. (26), (29), and (30) into the nonlinear root finder yields the actual values of $`E`$, $`L_z`$, and $`Q`$ to within machine precision in fewer than 10 iterations.

One caveat about our parameterization method is worth mentioning: some values of $`(r_p,e,\iota )`$ correspond to unstable Kerr orbits, and in this case the method returns a set $`(E,L_z,Q)`$ corresponding to an orbit with a different pericenter from the one requested. We can illustrate this behavior by factoring Eq. (17), which is a quartic function in $`r`$:

$$R(r)=(1E^2)(r_1r)(rr_2)(rr_3)(rr_4).$$
(31)
The roots are ordered so that $`r_1r_2r_3r_4`$. Bound motion occurs for $`r_1rr_2`$, which implies that $`r_1=r_a`$ and $`r_2=r_p`$, but this only works for stable orbits. In the case that the orbit *requested* is unstable, the set $`(r_p,e,\iota )`$ returned by the algorithm instead corresponds to a nearby, stable orbit. In this case, the numerically calculated roots of $`R(r)`$ satisfy $`r_3=r_{p,\mathrm{requested}}`$, but the parameterization method returns $`r_2`$ as the pericenter, i.e., it returns a nearby stable orbit with pericenter $`r_p=r_2>r_{p,\mathrm{requested}}`$. As a result, the *actual* pericenter is larger than the value requested.

We use routines from Hughes  to identify the boundary between stable and unstable orbits (so that the latter may be excluded), but the code has a few minor bugs, and the identification procedure is not infallible. As a result, some such orbits appear in the results below (Sec. V) and can be identified by having pericenters different from those requested. It is essential to understand, however, that the orbits returned by the parameterization algorithm are never unstable, and represent perfectly valid (stable) solutions to the equations of motion.

#### III.2.2 From Kerr geodesics to Papapetrou initial conditions

Once we have the set $`(E,L_z,Q)`$ for the geodesic case, we can use the well-known properties of the Kerr metric to solve for all the parameters necessary for the method described in Sec. III.1 above. The first step is to force the constants of the motion to agree by setting $`E_{\mathrm{Papapetrou}}=E_{\mathrm{Kerr}}`$ and $`J_z=L_z`$. Next, we set $`\theta _0=\pi /2`$ (corresponding to the equatorial plane), since this is the only angle guaranteed to be shared by all Kerr geodesic orbits.333This restricts our sample to Papapetrou orbits that cross the equatorial plane. This almost certainly represents the vast majority of valid Papapetrou solutions, but there remains the intriguing possibility of spinning particle solutions that orbit permanently above or below the equatorial plane. We leave an examination of this possibility to future investigators. Finally, we must choose an initial value $`r_0`$ for the Boyer-Lindquist radius, which in the case of $`\theta _0=\pi /2`$ coincides with the crossing of the equatorial plane. One possibility is simply to use the average value of the pericenter and apocenter,

$$r_0=\frac{1}{2}(r_p+r_a).$$
(32)
This radius is guaranteed to lie between pericenter and apocenter, and the prescription in Eq. (32) works fine for the mildly eccentric orbits considered in this paper, but a more robust method must take into account that highly eccentric orbits should in general cross the equatorial plane near pericenter. (Requiring a plane-crossing far from pericenter would force the inclination angle to be small, which is a constraint we do not wish to impose.) This suggests choosing an initial value of $`r`$ close to pericenter. One flexible method is to choose

$$r_0=r_p(1+\alpha e),$$
(33)
where $`\alpha `$ is a number of order unity. This reduces correctly to $`r_p`$ in the $`e=0`$ (circular) limit, and selects an equatorial plane crossing near pericenter in the $`e1`$ limit. If we set $`\alpha =2`$ in Eq. (33), then this method agrees exactly with Eq. (32) in the cases $`e=0`$ and $`e=0.5`$, and differs from Eq. (32) by less than 15% when $`e=0.6`$. Most of the results in this paper use Eq. (32), but Eq. (33) is preferable in general.

Once the initial $`\theta `$ and $`r`$ are known, we can determine all components of the Kerr 4-momentum $`p_K^\mu `$, but three of the Papapetrou momenta are determined by the constraints (Sec. III.1). We are therefore free to force only one of the four components of the Papapetrou momentum $`p_P^\mu `$ to be the same as its Kerr counterpart. For the bulk of this paper, we set the radial components equal ($`p_K^r=p_P^r`$), since it is the radial momentum that is most closely tied to the stability and boundedness of the orbits. This choice results in Papapetrou orbits with pericenters and eccentricities fairly close to the Kerr geodesic values, but with much higher orbital inclination angles (Fig. 7). The alternate choice of $`p_K^\theta =p_P^\theta `$ results in Papapetrou orbits with inclinations similar to their Kerr counterparts, but with very different pericenters (Fig. 8). See Sec. V for further discussion.

We determine a value for the radial Papapetrou momentum $`p_P^r`$ using the Kerr geodesic parameters $`E_K`$, $`L_z`$, and $`Q`$ by applying Eq. (17) and the equation

$$\rho ^2p^r=\mu \sqrt{R(r)},$$
(34)
where $`\rho =r^2+a^2\mathrm{cos}^2\theta `$ Misner et al. (1973). We then convert to $`p_r`$ using $`p_r=p^rg_{rr}`$. Proceeding exactly as in Sec. III.1, we specify two of the spin components and eliminate two variables using 4-momentum and spin normalization, and then solve numerically for $`p_t`$, $`p_\varphi `$, and $`S_t`$. The result is a set of initial conditions for the Papapetrou equations with the same energy and angular momentum as a Kerr geodesic with the desired values of $`r_p`$, $`e`$, and $`\iota `$.

It sometimes happens that the Papapetrou initial conditions derived in this manner specify an orbit that is unstable against plunge into the black hole. Since there is no “effective potential” for a generic Papapetrou orbit as there is for Kerr geodesics, there is no way *a priori* to detect this instability. Plunge orbits are detected at runtime by testing for radial coordinates less than the horizon radius.444In practice, the most common runtime error is actually a numerical underflow in the integration stepsize as the particle approaches the horizon. Orbits that plunge are removed from consideration since by definition they cannot be chaotic.

#### III.2.3 Empirical orbital parameters

In making the transition from geodesics to solutions of the Papapetrou equations, we are able to enforce the conditions $`E_{\mathrm{Papapetrou}}=E_{\mathrm{Kerr}}`$ and $`J_z=L_z`$, but this is no guarantee that the Papapetrou orbit has the corresponding orbital parameters $`r_p`$ and $`\iota `$: the spin-coupling term $`\frac{1}{2}\xi _{\mu ;\nu }S^{\mu \nu }`$ \[cf. Eq. (6)\] has a potentially large effect on the empirical values of the pericenter and orbital inclination. This effect is most pronounced when we consider high spin parameter values, i.e., $`S1`$. In these dynamically interesting cases, the *empirical* pericenter and inclination will differ in general from the values requested in the parameterization.

The empirical Papapetrou pericenter $`r_{p,P}`$ is easy to find: we simply integrate the orbit with a small stepsize for a large number of periods, and then record the minimum radius achieved. In practice, this works robustly, reproducing almost exactly the requested Kerr value of $`r_{p,K}`$ in the limit $`S1`$. The only exception involves values of $`r_{p,K}`$ that correspond to requested unstable orbits, as discussed in Sec. III.2.1 above. Each of these orbits has an empirical pericenter larger than the pericenter requested: the requested pericenter corresponds to the root $`r_3`$ in Eq. (31), but the empirical pericenter returned by the algorithm corresponds to the larger root $`r_2`$.

Having found the empirical pericenter for an orbit, we must next find its empirical orbital inclination angle $`\iota _P`$. In order to reproduce the definition from Eq. (23), we need to find an analogue of the Carter constant $`Q`$ for spinning test particles. Kerr spacetime has a Killing tensor $`K_{\mu \nu }`$ Wald (1984) that satisfies

$$K_{\mu (\nu ;\alpha )}=0,$$
(35)
which gives rise to an extra conserved quantity in the case of geodesic motion:

$$K=K_{\mu \nu }p^\mu p^\nu .$$
(36)
This quantity is related to the traditional Carter constant $`Q`$ by Misner et al. (1973)

$$K=Q+(L_zaE)^2.$$
(37)
When the test particle has nonzero spin, the quantity defined by Eq. (37) is no longer constant, but there is an analogous expression that is conserved to linear order. Adapting a result from Tanaka et al. (1996), we can write this approximately conserved quantity as

$$C=K_{\mu \nu }p^\mu p^\nu 2p^\mu S^{\rho \sigma }(f_\sigma ^\nu f_{\mu \rho \nu }f_\mu ^\nu f_{\rho \sigma \nu }),$$
(38)
where

$$f_{\mu \nu }=2a\mathrm{cos}\theta e_{[\mu }^1e_{\nu ]}^0+2re_{[\mu }^2e_{\nu ]}^3,$$
(39)
$$f_{\mu \nu \sigma }=6\left(\frac{a\mathrm{sin}\theta }{\sqrt{\mathrm{\Sigma }}}e_{[\mu }^0e_\nu ^1e_{\sigma ]}^2+\sqrt{\frac{\mathrm{\Delta }}{\mathrm{\Sigma }}}e_{[\mu }^1e_\nu ^2e_{\sigma ]}^3\right),$$
(40)
and the $`\{e_\mu ^a\}`$ are the standard orthonormal tetrad for the Kerr metric. The effective Carter “constant” for spinning particles is then

$$Q_{\mathrm{eff}}=C(J_zaE)^2,$$
(41)
where we use the full angular momentum $`J_z`$ (which includes the contribution from the spin) in place of $`L_z`$. The quantity $`Q_{\mathrm{eff}}`$ is nearly but not exactly not constant, so in order to define an empirical inclination angle we find the maximum effective $`Q`$ over an orbit, and then define $`\iota _P`$ by

$$Q_{\mathrm{eff},\mathrm{max}}=J_z^2\mathrm{tan}^2\iota _P.$$
(42)
As in the geodesic case, the value of $`\iota `$ obtained from Eq. (42) agrees well with the maximum value of $`|\pi /2\theta |`$ over an orbit.555This is true only when we force the Kerr and Papapetrou values of $`p^r`$ to agree (Sec. III.2.2), which is the case for all orbits considered except for the initial conditions shown in Fig. 8. In that case we revert to the simpler method of finding the maximum value of $`|\pi /2\theta |`$ over several orbits. When $`S=0`$, Eq. (42) reduces to the definition of the orbital inclination for geodesics, Eq. (23).

## IV Lyapunov exponents

### IV.1 The principal exponent

The Lyapunov exponents for a chaotic dynamical system quantify the chaos and give insight into its dynamics (revealing, for example, whether it is Hamiltonian or dissipative). For a dynamical system of $`n`$ degrees of freedom, in general there are $`n`$ Lyapunov exponents, which describe the time-evolution of an infinitesimal ball centered on an initial condition. This initial ball evolves into an ellipsoid under the action of the Jacobian matrix of the system, and the Lyapunov exponents are related to the average stretching of the ellipsoid’s principal axes. We described in Hartl (2003) a general method for calculating all $`n`$ of the system’s exponents, but in the present study we are interested only in the presence or absence of chaos in the Papapetrou system, so we need only calculate the principal Lyapunov exponent, i.e., the exponent corresponding to the direction of greatest stretching. A nonzero principal exponent indicates the presence of chaos.

When studying a differentiable dynamical system, we typically introduce a set of variables $`𝐲=\{y_i\}`$ to represent the system’s phase space, together with the autonomous set of differential equations

$$\frac{d𝐲}{d\tau }=𝐟(𝐲)$$
(43)
which determine the dynamics. Associated with each initial condition is a solution (or *flow*). The principal Lyapunov exponent quantifies the *local* divergence of nearby initial conditions, so any fully rigorous method necessarily involves the local behavior of the system, i.e., its derivative. For a multidimensional system, this derivative map is given by the Jacobian matrix,

$$(\mathrm{𝐃𝐟})_{ij}=\frac{f_i}{x^j}.$$
(44)
It is the Jacobian map that determines the time-evolution of infinitesimally separated trajectories. If we consider a point $`𝐲`$ on the flow and a nearby point $`𝐲+\delta 𝐲`$, then we have

$$𝐟(𝐲+\delta 𝐲)=\frac{d}{d\tau }(𝐲+\delta 𝐲)=\frac{d𝐲}{d\tau }+\frac{d(\delta 𝐲)}{d\tau },$$
(45)
so that the separation $`\delta 𝐲`$ satisfies

$$\frac{d(\delta 𝐲)}{d\tau }=𝐟(𝐲+\delta 𝐲)\frac{d𝐲}{d\tau }=𝐟(𝐲+\delta 𝐲)𝐟(𝐲).$$
(46)
Since

$$𝐟(𝐲+\delta 𝐲)𝐟(𝐲)=\mathrm{𝐃𝐟}\delta 𝐲+𝒪(\delta 𝐲^2),$$
(47)
we can write the time-evolution of the deviation vector as

$$\frac{d(\delta 𝐲)}{d\tau }=\mathrm{𝐃𝐟}\delta 𝐲+𝒪(\delta 𝐲^2).$$
(48)
If we identify the “infinitesimal” deviation $`\delta 𝐲`$ with an element $`𝝃`$ of the tangent space at $`𝐲`$, we effectively take the limit as $`\delta 𝐲0`$; the equation of motion for $`𝝃`$ is then

$$\frac{d𝝃}{d\tau }=\mathrm{𝐃𝐟}𝝃.$$
(49)
This equation describes the time-evolution of the separation between nearby initial conditions in a rigorous way. We will refer to $`𝝃`$ as a tangent vector, since formally is an element of the tangent space to the phase space.

If the Jacobian $`\mathrm{𝐃𝐟}`$ in Eq. (49) were some constant matrix $`𝐀`$, the solution would be the matrix exponential

$$𝝃(\tau )=\mathrm{exp}(𝐀\tau )𝝃_0.$$
(50)
For long times, the solution would be dominated by the largest eigenvector of $`𝐀`$, and would grow like

$$𝝃(\tau )e^{\lambda _{\mathrm{max}}\tau },$$
(51)
where $`\lambda _{\mathrm{max}}`$ is the largest eigenvalue.666The most common choice for the norm $``$ is the Euclidean norm, but we use a slightly different norm in the case of the Papapetrou equations (Sec. IV.2 below). (Here we have used $`𝝃_0=1`$.) Turning things around, if we measured the time-evolution of $`𝝃`$, we could find an approximation for the exponent using

$$\lambda _{\mathrm{max}}=\frac{\mathrm{log}𝝃(\tau )}{\tau }.$$
(52)
The general case follows by considering Jacobian matrices that are time-dependent. In this case, we are unable to define a unique principal exponent valid for all times, but there still is a unique *average* exponent that describes the average stretching of the principal eigenvector. Our method is to track the evolution of a tangent vector as it evolves into the principal axis of the ellipsoid. If we use $`r_e(\tau )=𝝃(\tau )`$ to denote the length of the longest principal ellipsoid axis, the principal Lyapunov exponent is then defined by

$$\lambda _{\mathrm{max}}=\underset{\tau \mathrm{}}{lim}\frac{\mathrm{log}[r_e(\tau )]}{\tau }.$$
(53)
We use an infinite time limit in this formal definition, but of course in practice a numerical approach relies on a finite cutoff to obtain a numerical approximation.

The *Jacobian method* for determining the largest Lyapunov exponent involves solving Eqs. (43) and (49) as a coupled system of differential equations in order to follow the time-evolution of a ball of initial conditions. One possibility is then to use Eq. (53) to estimate the system’s largest exponent. A related technique, which provides more accurate exponents, is to sample $`\mathrm{log}[r_e(\tau )]`$ at regular time intervals, and then perform a least-squares fit on the simulation data. Since $`\mathrm{log}[r_e(\tau )]=\lambda \tau `$, the slope of the resulting line then gives an estimate for the Lyapunov exponent. This is the method we implement in practice.

A less rigorous but still useful technique, which we call the *deviation vector* approach, involves solving only Eq. (43), but for two initial conditions: $`𝐲_0`$ and $`𝐲_0+\delta 𝐲_0`$. If the solutions to these initial conditions a time $`\tau `$ later are $`𝐲`$ and $`𝐲^{}`$, respectively, then the approximate principal Lyapunov exponent is

$$\lambda _{\mathrm{max}}=\frac{\mathrm{log}\delta 𝐲/\delta 𝐲_0}{\tau },$$
(54)
where $`\delta 𝐲𝐲^{}𝐲`$. This approach has a serious drawback: no matter how small the initial size of the deviation, eventually the method saturates as the difference between $`𝐲`$ and $`𝐲^{}`$ grows so large that it no longer samples the *local* difference between two trajectories.777It is possible to re-scale the deviation once it reaches a certain size, but this method is error-prone since it can depend sensitively on the precise method of rescaling. The constrained nature of the Papapetrou equations also presents difficulties for rescaling, as discussed in Sec. IV.2. On the other hand, because we need only solve Eq. (43) and not Eq. (49), the deviation vector method is significantly faster than the Jacobian method (by a factor of approximately 5 for the system considered here). We therefore adopt the deviation vector method as our principal tool for broad surveys of parameter space. The method for handling the saturation problem is discussed in Sec. IV.3.

A comparison of the Jacobian and deviation vector methods appears in Fig. 1. It is apparent that the two methods agree closely until the deviation vector approach reaches the saturation limit.

### IV.2 The Papapetrou case

The discussion in the previous section was of a general nature, applicable to virtually any dynamical system described by differential equations. Here we describe some of the details needed to apply the general methods to the Papapetrou equations. In particular, we must discuss further the key ideas of constrained deviation vectors and phase space metric.

For the Papapetrou system, the phase space vector $`𝐲`$ has 12 components:

$$𝐲=(t,r,\mu ,\varphi ,p_t,p_r,p_\theta ,p_\varphi ,S_t,S_r,S_\mu ,S_\varphi ).$$
(55)
The tangent vector $`𝝃`$ has one component for each variable. The set of equations Eq. (43) is simply the Papapetrou equations written out in full:

$`\dot{x}^\mu `$ $`=`$ $`v^\mu `$
$`\dot{p}_\mu `$ $`=`$ $`R_{\mu \nu }^{\alpha \beta }v^\nu p_\alpha S_\beta +\mathrm{\Gamma }_{\beta \mu }^\alpha p_\alpha v^\beta `$
$`\dot{S}_\mu `$ $`=`$ $`p_\mu \left(R_\beta ^{\alpha \gamma \delta }S_\alpha v^\beta p_\gamma S_\delta \right)+\mathrm{\Gamma }_{\beta \mu }^\alpha S_\alpha v^\beta ,`$
where we have used the formulation in terms of the spin 1-form described in the appendix. The second necessary equation, Eq. (49), requires the Jacobian matrix,

$$\left(\begin{array}{ccc}\hfill \frac{\dot{x}^\mu }{x^\nu }& \frac{\dot{x}^\mu }{p_\nu }& \frac{\dot{x}^\mu }{S_\nu }\hfill \\ \hfill \frac{\dot{p}_\mu }{x^\nu }& \frac{\dot{p}_\mu }{p_\nu }& \frac{\dot{p}_\mu }{S_\nu }\hfill \\ \hfill \frac{\dot{S}_\mu }{x^\nu }& \frac{\dot{S}_\mu }{p_\nu }& \frac{\dot{S}_\mu }{S_\nu }\hfill \end{array}\right),$$
(56)
whose explicit form appears in Hartl (2003).

It is important to mention that the tangent vector $`𝝃`$—or, equivalently, the deviation vector $`\delta 𝐲`$—cannot have completely arbitrary components. On the contrary, the deviation must be chosen carefully, in order to ensure that, given a point $`𝐲`$ that satisfies the constraints from Sec. II.2, the deviated vector $`𝐲+\delta 𝐲`$ also satisfies the constraints. Otherwise, the relation $`𝐟(𝐲+\delta 𝐲)𝐟(𝐲)=\mathrm{𝐃𝐟}\delta 𝐲+𝒪(\delta 𝐲^2)`$ is not satisfied. (This is the principal complication in implementing a rescaled version of the deviation vector method: the rescaled vector would violate the constraint.) In practice, we are able to find a valid deviation vector by applying the same techniques used to satisfy the constraints in the first place (Sec. III); for details, see Hartl (2003).

One final detail is the notion of metric: implicit in the definition of the Lyapunov exponent, Eq. (53), is a metric on phase space used to calculate the norm of the tangent vector $`𝝃`$. We adopt a metric introduced in Karas and Vokrouhlický (1992), which involves projecting the deviation vector onto the spacelike hypersurface defined by the zero angular momentum observers (ZAMOs). The projection is effected using the projection tensor $`P_\nu ^\mu =\delta _\nu ^\mu +U^\mu U_\nu `$, where $`U^\mu `$ is the ZAMO 4-velocity. The spatial and momentum variables are then projected according to $`x^\mu \stackrel{~}{x}^i=P_\mu ^ix^\mu `$, $`p_\mu \stackrel{~}{p}_i=P_i^\mu p_\mu `$, and $`S_\mu \stackrel{~}{S}_i=P_i^\mu S_\mu `$. After the projection, we calculate the Euclidean norm in the 3-dimensional hypersurface. We note that while this prescription is convenient, and reduces correctly in the nonrelativistic limit, the magnitudes of the Lyapunov exponents are similar for several other possible choices of metric Hartl (2003).

### IV.3 Chaos detector

Since we are concerned with calculating Lyapunov exponents for a large number of parameter values, we use the (unrescaled) deviation vector method because of its speed. As mentioned in Sec. IV.1, this method has the property of saturation (as illustrated in Fig. 1), which is ordinarily a problem, but here we use it to our advantage as a sensitive detector of chaos.

Our method for determining whether a particular initial condition is chaotic is to consider a nearby initial condition separated by a small vector $`\delta 𝐲_0`$ (with norm $`\delta 𝐲_0`$ typically of order $`10^7`$ or $`10^8`$) and then integrate until the system reaches 90% of the saturation limit, defined as a separation $`\delta 𝐲`$ with unit norm. If we write $`r_e=\delta 𝐲/\delta 𝐲_0`$, the approximate Lyapunov exponent satisfies

$$\mathrm{log}[r_e(\tau )]=\lambda \tau ,$$
(57)
so that $`\lambda `$ is the slope of the line $`\mathrm{log}[r_e(\tau )]`$ vs. $`\tau `$. Saturation occurs when $`r_e=1/\delta 𝐲_0`$, so that the integration ends when $`\mathrm{log}[r_e(\tau )]=\mathrm{log}(0.9\delta 𝐲_0)`$. We record the value of $`\mathrm{log}[r_e(\tau )]`$ at regular time intervals (typically every time $`T=100M`$ in our case), and upon reaching 90% saturation perform a least-squares fit on the simulation data to determine the exponent.

We note that the cutoff value of $`0.9`$ is somewhat arbitrary and is the result of numerical experimentation. When using the (unrescaled) deviation vector approach, most of the chaotic systems saturate—that is, plots of $`\mathrm{log}[r_e(\tau )]`$ vs. $`\tau `$ flatten out (Fig. 1)—when the separation is of order unity, corresponding to radial separations of order $`M`$, velocity separations of order $`1`$, and angular separations of order $`1\mathrm{radian}`$. The 90% prescription ends the integration before the growth flattens out, so that the numerical estimate for the exponent is still reasonably accurate.

For the large maps of parameter space, we typically integrate up to $`\tau _{\mathrm{final}}=10^5M`$ or saturation, whichever comes first. We choose this maximum time mainly for practical reasons: it is the longest integration possible in a reasonable amount of time. \[We integrate as deeply as $`10^7M`$ (Sec. V.7) for individual orbits, but such long integrations are impractical for more than a handful of parameter combinations.\] Dramatically longer integrations are also not particularly useful, since the timescale for gravitational radiation reaction is on the order of $`\tau _{\mathrm{GW}}=M^2/\mu `$ Hughes (2001), which for the most relevant LISA sources is $`10^4`$$`10^6M`$ (i.e., $`\mu 10^6`$$`10^4M`$). Searching for chaos in such systems on a timescale longer than $`10^7M`$ is probably pointless, since the radiation reaction would dominate the dynamics in this case. Finally, it appears that chaos, when present in the Papapetrou system, manifests itself on relatively short timescales ($`10^2`$$`10^4M`$), or else not at all. The onset of chaos is marked by a qualitative change from power law growth (which appears as logarithmic growth on our plots of $`\mathrm{log}[r_e(\tau )]`$ vs. $`\tau `$) to exponential growth (which is linear on the same plots). An example of two similar initial conditions giving rise to qualitatively different dynamical behavior appears in Fig. 2.

One important refinement to the technique described above is to require several 90% saturation points *in a row* before declaring the orbit to be chaotic. This is necessary because some nonchaotic orbits have very high amplitudes on plots of $`\mathrm{log}[r_e(\tau )]`$ vs. $`\tau `$ (Fig. 3). This phenomenon occurs mainly for orbits with many periods in the deeply relativistic zone near the horizon. Such orbits may reach “90% saturation” briefly as part of their oscillation, but quickly return to separations far below our chaotic cutoff. We therefore adopt the criterion of *three 90% saturation points in a row* (with a time $`T=100M`$ sampling interval) as a robust practical test for chaos. See Sec. V.8 for further discussion.

Our confidence of this method’s robustness derives from comparing the method above to the Jacobian method for the same initial conditions. Since the Jacobian method does not saturate, the agreement of the two methods indicates that our procedure provides an accurate detector for chaos (as in Fig. 12 below).

### IV.4 Implementation notes

We integrate the Papapetrou equations on a computer using Bulirsch-Stoer and Runge-Kutta integrators implemented in the C programming language, as described in Hartl (2003). The derivatives and Jacobian matrix are extensively hand-optimized for speed. We monitor errors using constraints and conserved quantities, with a global error goal of $`10^{13}`$. The errors are at the $`10^{11}`$ level or better for highly chaotic orbits after $`10^5M`$. Orbits with low spin or high pericenter are even more accurate, often achieving the error goal of $`10^{13}`$.

The many plots in Sec. V are typically generated using driver programs written in the Perl programming language, which in turns calls the C integrator repeatedly. This general paradigm—using an interpreted language such as Perl to call optimized routines in a compiled language such as C—is one we warmly recommend.

## V Results

We present here the results of parameter variation in the Papapetrou system of equations. We represent the effects compactly using several different kinds of plots, most of which involve plotting inclination vs. pericenter, with other parameters held fixed. We refer to these as $`r_p`$-$`\iota `$ maps. As discussed in Sec. III, starting with the Kerr values $`r_{p,K}`$, $`\iota _K`$, and $`e_K`$, we find the corresponding energy $`E_K`$, angular momentum $`L_{z,K}`$, and Carter constant $`Q_K`$, and then force the Papapetrou values to satisfy $`E_P=E_K`$, $`J_{z,P}=L_{z,K}`$ and $`p_P^r=p_K^r`$. Each $`r_p`$-$`\iota `$ map has two components: part (a), shown on the left, uses the Kerr parameters $`\iota _K`$ and $`r_{p,K}`$ requested by the parameterization (Sec. III.2.2), while part (b), shown on the right, always shows the empirical Papapetrou values $`\iota _P`$ and $`r_{p,P}`$ in the sense of Sec. III.2.3.

One important feature of $`r_p`$-$`\iota `$ parameter space apparent in the empirical plots is the prevalence of large empirical inclination angles for all values of the Kerr inclination $`\iota _K`$. Fig. 7 shows the mapping for $`S=1`$ between the requested Kerr parameters \[Fig. 7(a)\] and the empirical Papapetrou parameters \[Fig. 7(b)\] using a shading scheme (which appears as a more informative color scale in electronic versions of this paper). We see that even the orbits at the bottom of Fig. 7(a) get mapped to high empirical inclination angles; $`\iota _K=1^{}`$ orbits are mapped to inclinations of order $`\iota _P=40^{}`$.

This compression of parameter space is the result of our choice to force the Kerr and Papapetrou value of the radial momentum to agree (Sec. III.2.2). The price we pay for this choice is that the $`\theta `$ component of the Papapetrou momentum—which is constrained to satisfy the equations in Sec. II.2—is relatively large. This high $`p_P^\theta `$ flings even orbits with low requested values of $`\iota _K`$ to high inclinations (Fig. 7). It is possible to find low-inclination Papapetrou orbits by requiring that the Kerr and Papapetrou values of $`p^\theta `$ agree, but at the cost of forcing $`p^r`$ to be very different—again a result of the constraints. The resulting parameter space (Fig. 8) is not nearly as compressed in inclination angle, but the Papapetrou pericenters are compressed and shifted down, and many requested values of $`r_p,K`$ are lost as plunge orbits. In addition, the empirical values of the eccentricity are typically not close to the value requested (reaching, e.g., $`e_P=0.75`$ for $`e_K=0.5`$). Because of the deficiencies of this alternate parameterization method, we choose the fixed $`p^r`$ plots are our primary investigative tool in this paper.

### V.1 Varying pericenter and orbital inclination

In this section we show $`r_p`$-$`\iota `$ maps for orbits with fixed eccentricity $`e=0.5`$ and spin parameters $`S^{\widehat{r}}=S^{\widehat{z}}=0.2S`$, where $`S`$ is the total spin.888These “fixed” spin parameters give rise to a variety of spin inclination angles in the fiducial (ZAMO) rest frame; see Sec. V.6 below. We indicate the strength of the chaos at each point with a grayscale rectangle, with the darkest colors representing the strongest chaos (and with white indicating no chaos). An example of such a plot appears in Fig. 9. (Fig. 10 shows a similar plot for the alternate parameterization method that forces $`p_K^\theta =p_P^\theta `$. The maximum exponents in the two cases are virtually identical.) The most important general results evident from the plots are two-fold. First, the largest exponents occur for orbits with pericenters deep in the relativistic zone near the horizon. Second, the prevalence of chaotic orbits is a decreasing function of spin parameter $`S`$, with virtually all chaos gone by the time $`S=0.1`$.

An example of the value of the empirical $`r_{p,P}`$-$`\iota _P`$ plots appears in Fig. 26, which shows orbits of particles with spin $`S=0.5`$. The appearance of a strongly chaotic point in Fig. 26 seems perplexing, given that it is surrounded but many points with much smaller exponents. As is evident from the empirical plot, this point of strong chaos (which is, in fact, the largest Lyapunov exponent for any of the plots) maps to a compressed area of initial conditions with small empirical pericenters \[Fig. 26(b)\].

From an astrophysical standpoint, the most interesting parameter to vary is the spin $`S`$, since only $`S1`$ orbits are physically realistic (Sec. II.4). From Figs. 2529, which involve varying $`S`$ from 0.9 down to $`10^4`$, we see that both the prevalence and strength of chaos decrease significantly as $`S`$ is decreased. The Lyapunov exponent ranges as high as $`10^2M^1`$ when $`S=0.5`$ (Fig. 26), but the chaos is weak when $`S=0.2`$ (Fig. 27) and is gone below $`S=0.1`$ (Figs. 28 and 29).

### V.2 Varying eccentricity

The choice of $`e=0.5`$ in the previous section is partially motivated by likely gravitational wave sources for LISA LIS , e.g., a neutron star or small black hole in an eccentric orbit around a supermassive black hole. In this section, we consider a second series of eccentric orbits at fixed $`e=0.6`$ and varying spin parameter. We also investigate the case of a near-circular orbit ($`e=0.01`$) more appropriate for the circularized gravitational wave sources important for ground-based detectors such as LIGO.

The $`r_p`$-$`\iota `$ plots for $`e=0.6`$ follow the same pattern as those with $`e=0.5`$. Chaos is strongest for orbits with small pericenters and values of $`S`$ of order unity (Figs. 3033). There is a single orbit at $`S=0.1`$ that appears to be chaotic (Fig. 32), but other than this one exception there is apparently no chaos below $`S=0.1`$. A close examination of the single $`S=0.1`$ chaotic orbit appears in Figs. 1112, which shows that the chaos is real.

The effect of chaos is smaller for the near-circular ($`e=0.01`$) orbits considered (Figs. 3437), with typical Lyapunov exponents of order $`2\times 10^3M^1`$ when $`S=1`$ (Fig 34). Moreover, we find only four points with nonzero exponents for $`S=0.5`$ (Fig. 35), in contrast to the more eccentric orbits, which have strong chaos for $`S=0.5`$. By the time $`S=0.1`$ (Fig. 36), all chaos has completely disappeared for the near-circular orbits.

### V.3 Varying Kerr parameter $`a`$

Here we investigate the effect of the Kerr parameter $`a`$ on the strength and prevalence of chaos. Although Suzuki & Maeda found in Suzuki and Maeda (1997) that there is chaos even in Schwarzschild spacetime, to which Kerr spacetime reduces when $`a=0`$, the chaotic orbits found in Suzuki and Maeda (1997) are exceptional orbits lying on the edge of a generalized effective potential. We found evidence in Hartl (2003) that such chaotic orbits are rare.

The conclusion that chaotic orbits become less prevalent as $`a0`$ is supported by a more thorough examination, as illustrated in Figs. 3843. All of these orbits have eccentricity $`e=0.5`$ and spin $`S=1`$, and $`a`$ varies from $`0.9`$ down to 0. Even for low values of $`a`$, the parameter space is strongly affected by the spin, with plots of the empirical pericenter and orbital inclination showing significant distortion. Nevertheless, we see unambiguously that the chaotic orbits prevalent when $`a=1`$ are greatly suppressed as $`a`$ decreases, with no chaotic orbits at all below $`a=0.2`$. This appears to be a result of the increase of the marginally stable radius $`r_{\mathrm{ms}}`$ as $`a0`$. When $`a=0`$, the minimum stable radius is at $`r_{\mathrm{ms}}=4`$ in units of the central black hole’s mass, fully $`2M`$ away from the horizon at $`r_H=2M`$. As discussed in Hartl (2003), the extra (spherical) symmetry of the Schwarzschild metric leads to an additional integral of the motion, which also has a suppressive effect on chaos.

### V.4 Spin cutoffs for chaos

In this section we provide *spin cutoff values* for chaos, i.e., the minimum spin values for which chaos exists. For a given initial condition defined in terms of fixed orbital parameters (as described in Sec. III.2), we vary the spin parameter $`S`$ and find the maximum value for which chaos occurs. The smaller this cutoff value, the stronger the chaos: nonchaotic orbits have a cutoff value of $`S=1`$, i.e., they are not chaotic even in the extreme $`S=1`$ limit; conversely, a cutoff value of $`S=10^5`$ would indicate chaos for the (physically realistic) value $`S=10^5`$, but not for any smaller values.999Implicit in this scheme is an assumption of monotonicity, i.e., monotonically decreasing chaos as $`S`$ decreases. While not strictly true (as discussed in Hartl (2003)), this assumption is still broadly applicable, and exceptions are rare.

Fig. 13 is an example of a spin cutoff map. The procedure for producing such a map is similar to the method used to make the Lyapunov $`r_p`$-$`\iota `$ maps: we consider a grid of points in $`r_p`$-$`\iota `$ space, and for each point we find an approximate value for the spin cutoff. We begin by finding out if the system is chaotic for $`S=1`$, using the Lyapunov map as a start. If the orbit for $`S=1`$ is chaotic, we halve the spin value and calculate the Lyapunov exponent for $`S=0.5`$. If the system is still chaotic, we consider $`S=0.25`$; otherwise, we consider $`S=0.75`$. The procedure continues until the difference between chaotic and nonchaotic spin values is smaller than some threshold, which we choose to be $`0.05`$. (This value is chosen to achieve reasonable accuracy while still completing the calculations in a tolerable amount of time.)

Figs. 1315 show the spin cutoff values for several parameter combinations corresponding to Lyapunov maps from Sec. V.1 above. The plots are color-coded with grayscale so that the most chaotic points—those with the *smallest* spin cutoff values—appear darkest. The points surrounded by white are not chaotic even for $`S=1`$. The cutoffs for the darkest points depend on the parameter values: for Fig. 13 ($`a=1`$ and $`e=0.5`$), the same as Fig. 9 above), we find $`S_{\mathrm{cutoff}}=0.18457`$; Fig. 14 ($`a=1`$ and $`e=0.01`$, the same as Fig. 34 above) has $`S_{\mathrm{cutoff}}=0.65625`$; and Fig. 15 ($`a=0.5`$ and $`e=0.5`$, the same as Fig. 40) has $`S_{\mathrm{cutoff}}=0.28125`$. We should not take the digits of precision seriously, since these values are only accurate to within $`0.05`$, but *in all cases the spin cutoff values are significantly above the physically realistic range of $`S10^4`$$`10^7`$*.

### V.5 Retrograde orbits

We have considered a wide variety of orbits—varying eccentricity and Kerr parameter for different pericenter, orbital inclination, and spin parameters—but so far all orbits have satisfied $`0<\iota <\pi /2`$, i.e., they have all been prograde orbits, moving in the same direction as the central black hole’s spin. We investigate now the case of retrograde orbits, and show that they are poor candidates for chaos.

It is evident from looking at an $`r_p`$-$`\iota `$ plot of a retrograde orbit (Fig. 16) that the pericenters are much larger than their prograde counterparts. For the $`S=1`$ particle illustrated in Fig. 16, the minimum empirical pericenter is larger than $`6M`$, in contrast to prograde orbits, which get below $`1.5M`$. Furthermore, although it is clear from Fig. 16(b) that the parameter space is severely distorted, the chaos is extremely weak. The largest Lyapunov exponent, even in this extreme $`S=1`$ case, is $`\lambda _{\mathrm{max}}=3.5\times 10^4`$, two orders of magnitude smaller than the prograde case. Unsurprisingly, all chaos disappears when $`S1`$. The smallest value of $`S_{\mathrm{cutoff}}`$ is $`0.65265`$ for the parameter values shown in Fig. 16; we find no evidence of chaos below this value of $`S`$.

### V.6 Varying spin inclination

So far we have always used the same values for the two spin components passed to the parameterization procedure (scaled by the total spin $`S`$): $`S^{\widehat{r}}=S^{\widehat{z}}=0.2S`$. We consider now the effect of varying these parameters, and also discuss the corresponding initial spin inclination angles in a fiducial rest frame.

We begin with an initial condition that is chaotic for $`S=1`$ but is otherwise arbitrary. We then vary $`S^{\widehat{r}}`$ and $`S^{\widehat{z}}`$ and calculate the Lyapunov exponent for each configuration. The result for $`a=1`$, $`e_K=0.5`$, $`r_{p,K}=2.3`$ and $`\iota _K=20^{}`$ appears in Fig. 17. When $`S^{\widehat{r}}=S^{\widehat{z}}=0.2`$, the parameter values chosen correspond to a point from Fig. 9. There are not valid initial conditions for all choices of parameter; in particular, negative values of $`S^{\widehat{z}}`$ are often unstable or are unable to satisfy the spin constraints. Nevertheless, there is a large variety of parameters that do give rise to valid orbits, and many of them are chaotic. The strongest chaos exists for orbits that have small values of $`S^{\widehat{z}}`$, but this appears to be mainly because such orbits are able to achieve small empirical pericenters. As Fig. 18 clearly shows, the Lyapunov exponent generally decreases as pericenter increases, with no chaotic orbits above $`r_{p,P}=2.5M`$.

While we are forced by the constraints to parameterize the equations of motion in terms of spin components, it is more convenient to visualize the spin in a fiducial rest frame that is hypersurface-orthogonal to the particle’s trajectory. In Kerr spacetime, the hypersurface-orthogonal observers are the zero angular momentum observers (ZAMO), which is the same frame we use when calculating the Lyapunov exponents. By projecting the components of the spin vector $`S^\mu `$ into the ZAMO frame in the same way as we do for the projected norm (Sec. IV.2), we can find the *local* value of the spin inclination angle $`\theta _{\mathrm{local}}`$ (i.e., the angle between the spin axis and the axis of the central black hole). The results appear in Table 1. It is clear that our variation of spin components samples a large variety of initial spin inclination angles.

### V.7 Deep integrations

Since we adopted a maximum time of $`10^5M`$ for the Lyapunov integrations, it is reasonable to ask whether chaos might manifest itself on a longer timescale. This is certainly possible, but it appears that most initial conditions are either chaotic on a timescale of order $`10^2`$$`10^4M`$ or are not chaotic at all, as discussed in Sec. IV.3. An example appears in Fig. 2, where one initial condition is unambiguously chaotic, while a second located close by is not chaotic, even on a much longer timescale.

To convince ourselves that slow chaos is not lurking in the apparently nonchaotic regions, we performed a few longer-time integrations. In particular, we calculated the Lyapunov exponents using $`\tau _{\mathrm{final}}=10^7M`$ for all the innermost ($`r_p=1.32`$) orbits from Fig. 28 ($`S=0.1`$) and Fig. 29 ($`S=10^4`$), which are strongly chaotic when $`S=1`$ but are apparently without chaos below $`S=0.1`$. The largest exponent occurs for $`\iota =28.5^{}`$, which is therefore the worst-case scenario. Plots of the Lyapunov exponents vs. time appear in Fig. 19 ($`S=0.1`$) and Fig. 20 ($`S=10^4`$); their magnitudes are on the order of $`\lambda _{\mathrm{max}}=3\times 10^7M^1`$, corresponding to an $`e`$-folding timescale of approximately $`3.3\times 10^6M`$. For comparison, we show $`\lambda `$ vs. $`\tau `$ for a chaotic initial condition in Fig. 21; it is clear that the Lyapunov exponent asymptotes to a nonzero value in much less than $`10^5M`$, even for weak chaos. (The initial condition in Fig. 21 is the $`S=0.1`$ orbit illustrated in Figs. 11 and 12, whose Lyapunov exponent is actually quite small compared to analogous $`S=1`$ orbits.) The long integrations thus provide strong evidence that the disappearance of nearly all chaotic orbits below $`S=0.1`$ is a real effect.

### V.8 $`S=0`$ and chaos mimics

Since the case of $`S=0`$ corresponds exactly to geodesic orbits in Kerr spacetime, such systems cannot be chaotic—the return of the Carter constant $`Q`$ and the loss of the spin degrees of freedom make the system fully integrable. Nevertheless, even some geodesic orbits can have a large separation of nearby initial conditions, which can appear to be chaotic. These chaos mimics typically spend many orbital periods whirling around deep in the strongly relativistic zone near the horizon, only occasionally zooming out to higher radii. These so-called *zoom-whirl* orbits may provide significant challenges to detection despite their formal integrability.

An example of how much divergence an $`S=0`$ orbit can experience appears in Fig. 3. A picture of the corresponding orbit (visualized in Boyer-Lindquist coordinates embedded in ordinary space) appears in Fig. 22, which makes clear the large number of low-radius $`\varphi `$-periods characteristic of zoom-whirl orbits. A second example of a chaos mimic appears in Figs. 24 and 23. This orbit, in contrast to the previous one, does not have a particularly small pericenter, but its high inclination angle and zoom-whirliness allow it to mimic chaotic orbits.

The chaos mimics can exhibit large growth of the initial deviation vector, approximately a factor of $`10^6`$$`10^8`$, after $`10^5M`$. The principal means for detecting them is by requiring several high-separation points in a row (on a time-$`T100M`$ basis) as mentioned in Sec. IV.3; the mimics have oscillations with high amplitudes due to their zoom-whirliness, but they do not represent true saturations of the separation vector.

## VI Conclusions

The Papapetrou equations, which model a spinning test particle, exhibit chaotic solutions in Kerr spacetime for a wide range of parameters. In terms of the mass $`M`$ of the central (Kerr) black hole, the largest Lyapunov exponents are of order $`\lambda _{\mathrm{max}}=0.01M^1`$, which represents an exponential divergence of trajectories on a timescale of $`\tau _\lambda =1/\lambda =100M`$. Furthermore, there are many chaotic orbits with exponents in the range $`10^3`$$`10^4M^1`$. Despite the large number of chaotic orbits, we find that values of $`\lambda `$ corresponding to unambiguous chaos occur exclusively when the spin parameter $`S`$ is not small compared to unity. In particular, we find virtually no chaos for spin values below $`S=0.1`$, and no evidence of any chaos for spins below $`S=10^4`$.

The strongest determinant of chaotic behavior, apart from the spin parameter, is the pericenter of the orbit in question. The most highly chaotic orbits are those that reach pericenters near the horizon of the black hole. This is due to the high spacetime curvature in these regions, which maximizes the size of the coupling of the spin to the Riemann curvature tensor \[Eq. (II.1\]. When the Kerr parameter $`a`$ is small, so that the Kerr metric differs only slightly from the Schwarzschild metric, the pericenters are much higher than in the extreme Kerr ($`a=1`$) case. Chaos in the Papapetrou system is therefore weak when $`a`$ is small. The prevalence of chaos is also dependent on orbital eccentricity. Near-circular ($`e=0.01`$) orbits have many fewer regions of chaotic orbits than those with higher eccentricities ($`e=0.5`$ or $`e=0.6`$). This seems due primarily to the lower pericenters accessible to high-eccentricity orbits.

The dependence of the Lyapunov exponents on $`S`$ is our most important result: in all cases considered, physically realistic values of $`S`$ (satisfying $`S1`$) are not chaotic. We have shown conclusively that the Papapetrou equations admit many solutions that are formally chaotic, but without exception such chaotic solutions occur only for relatively large values of $`S`$. Below the upper limit for physically realistic spins ($`S10^4`$), we find no evidence of chaotic solutions. As a practical matter, this means that chaos will not manifest itself in the gravitational radiation from extreme mass-ratio binary inspirals.

###### Acknowledgements.

I would like to thank Scott Hughes and Teviet Creighton for providing notes on parameterizing Kerr geodesics in terms of orbital parameters. I also thank Sterl Phinney for his careful reading of the manuscript and perceptive comments. This work was supported in part by NASA grant NAG5-10707.

*

## Appendix A Spin vector formulation

We summarize here the formulation of the Papapetrou equations in terms of the spin 1-form $`S_\mu `$ (often referred to loosely as the “spin vector”), as mentioned in Sec. II.1. In this paper we consider a spinning particle of rest mass $`\mu `$ orbiting a central Kerr black hole of mass $`M`$, and it is convenient to measure all times and lengths in terms of $`M`$ and all momenta in terms of $`\mu `$. In these normalized units, the equations of motion in terms of the spin 1-form are

$`{\displaystyle \frac{dx^\mu }{d\tau }}`$ $`=`$ $`v^\mu `$
$`_\stackrel{}{v}p_\mu `$ $`=`$ $`R_{\mu \nu }^{\alpha \beta }v^\nu p_\alpha S_\beta `$ (58)
$`_\stackrel{}{v}S_\mu `$ $`=`$ $`p_\mu \left(R_\beta ^{\alpha \gamma \delta }S_\alpha v^\beta p_\gamma S_\delta \right)`$
where

$$R_\beta ^{\alpha \mu \nu }=\frac{1}{2}R_{\beta \rho \sigma }^\alpha ϵ^{\rho \sigma \mu \nu }.$$
(59)
The supplementary condition Eq. (5) allows for an explicit solution for the 4-velocity $`v^\mu `$ in terms of $`p^\mu `$:

$$v^\mu =N(p^\mu +w^\mu ),$$
(60)
where

$$w^\mu ={}_{}{}^{}R_{}^{\mu \alpha \beta \gamma }S_\alpha p_\beta S_\gamma $$
(61)
and

$${}_{}{}^{}R_{}^{\alpha \beta \mu \nu }=\frac{1}{2}R^{\alpha \beta \rho \sigma }ϵ_{\rho \sigma }^{\mu \nu }.$$
(62)
The normalization constant $`N`$ is fixed by the constraint $`v_\mu v^\mu =1`$.

The spin 1-form satisfies two orthogonality constraints:

$$p^\mu S_\mu =0$$
(63)
and

$$v^\mu S_\mu =0.$$
(64)
These two constraints are equivalent as long as $`v^\mu `$ is given by Eq. (60): $`0=v^\mu S_\mu p^\mu S_\mu +w^\mu S_\mu =p^\mu S_\mu `$, since by definition of $`w^\mu `$ \[Eq. (61)\] the second term involves the contraction of a symmetric tensor with an antisymmetric tensor and therefore vanishes. We enforce Eq. (63) in our parameterization scheme, and we use Eq. (60) in the equations of motion, so Eq. (64) is then automatically satisfied.


format: text

text:
# Subgraphs in random networks

## I Introduction

Many natural systems are described as networks of interacting components (Strogatz 2001 -Maslov 1998 ). Random networks have been studied as models of these complex systems. The classic model for a random network is the Erdős model (Erdos1959 -Bollobas1985 ), in which each of the possible edges in the network exists with probability $`p`$. There exists an analytical solution to many of the properties of Erdős networks, such as the diameter, clustering coefficient, component size distributions, and subgraph distributions (Erdos1959 -Bollobas1985 ). The average number of appearances $`G`$ of a subgraph with n nodes and g edges in a directed network of $`N`$ nodes is

$`G=\lambda \left({\displaystyle \genfrac{}{}{0pt}{}{N}{n}}\right)p^g(1p)^{n(n1)g}\lambda N^n\left({\displaystyle \frac{K}{N}}\right)^g`$
$`GN^{ng}`$ (1)
assuming a fixed mean connectivity $`K=pN`$. $`\lambda `$ is a term of order 1 which stems from the symmetry of each subgraph. Erdős networks have been extensively used as models for analyzing real networks. An excellent example is the work of Davis, Holland and Leinhardt on subgraphs in social networks(Holland -Wasserman ).

Erdős networks exhibit a Poissonian degree distribution: the distribution of the number of edges per node is $`P(k)=k^ke^k/k!`$. Nodes with a number of edges much higher than the mean are exponentially rare. Many naturally occurring networks, on the other hand, obey a long-tailed degree sequence, often described as a power law, $`P(k)k^\gamma `$, with $`\gamma `$ often between 2 and 3 (Barabasi1999 -Sole 2001 ). These networks, termed scale-free networks, are characterized by the existence of nodes with high degree, termed hubs (Fig1). The existence of hubs dramatically influences the properties of these networks. Some of the global properties of random networks with arbitrary degree distribution, and specifically scale-free networks, have been calculated. These include sizes of connected components (Newman 2001 ,Aiello 2001 ,Molloy 1998 ), distances (Chung\_diameter ), percolation thresholds (Cohen 2002 -Newman 1999 percolation ) and clustering coefficients (Newman book ,dorogovtsev 2002a ,Ravasz 2002 ).

There is much current interest in the local structure of networks (Watts 1998 ,Shen or 2002 -Maslov 1998 ,Newman book -Eckmann2 ,Maslov B -Newman 1998 ). Recently subgraph structure was analyzed in biological and technological networks(Shen or 2002 ,Milo 2002 ). It was found that these natural or designed networks contain network motifs, subgraphs that occur much more often than in an ensemble of randomized networks with the same degree sequence. In biological networks, the network motifs were suggested to be elementary building blocks which carry out key information processing functions (Shen or 2002 ,Milo 2002 ). In these studies, random networks generation and the enumeration of their subgraphs were performed numerically. To complement this numerical work, it would be important to theoretically characterize the subgraph distribution of random networks. Here we present approximate formulas for the average number of subgraphs in an ensemble of random networks with an arbitrary degree sequence. In the random ensemble each node has a specified indegree, outdegree, and mutual degree. These formulas give a very good approximation for random networks which allow for multiple edges between nodes (more than one edge in a given direction), as in the well-studied configuration model (Bollobas1985 ,Molloy 1998 -Chung\_diameter ). We also show that they provide a reasonable approximation for networks where multiple edges are not allowed, which represent more realistically many naturally occurring networks. We apply these formulas to arrive at new scaling laws for networks with a scale-free degree distribution. We find that each subgraph has its own scaling exponent, influenced by its topology. Considering the power exponent of the degree distribution, $`\gamma `$, as a control parameter, we show that the random networks exhibit transitions between 3 regimes. In each regime the subgraph number of appearances follows a different scaling law. We find that certain subgraphs appear much more frequently than in Erdős networks.

## II Number of Subgraphs, approximate solution

The following approximation assumes sparse networks ($`KN`$). The network degree sequence is given by the outdegree $`\{K_i\}_{i=1}^N`$ (the number of edges outgoing from each node), indegree $`\{R_i\}_{i=1}^N`$ (the number of incoming edges at each node), and mutual degree $`\{M_i\}_{i=1}^N`$ (the number of mutual edges at each node). Mutual edges are cases where there is a pair of edges in both directions between two nodes. This property has been studied in social networks (Holland -Wasserman ) and in the world wide web (Eckmann ). We begin by computing the probability of obtaining an n-node subgraph with $`g_a`$ single edges, $`g_m`$ mutual edges, subgraph outdegree sequence $`\{k_j\}_{j=1}^n`$, subgraph indegree sequence $`\{r_j\}_{j=1}^n`$ and subgraph mutual degree sequence $`\{m_j\}_{j=1}^n`$ in a given set of nodes. Consider the example of (Fig 2). The probability of obtaining a directed edge from node 1 to node 2 is approximately

$`P(edge1)={\displaystyle \frac{K_1R_2}{NK}}`$ (2)
assuming $`K_1R_2NK`$ (see Appendix A). The probability of obtaining a second edge from node 1 to node 3 is :

$`P(edge2|edge1)={\displaystyle \frac{(K_11)R_3}{NK}}`$ (3)
This reasoning applies to all the subgraph edges. The mean number of appearances of a subgraph is found by taking the average of the resulting expression with respect to all choices of $`n`$ distinct nodes $`\{\sigma _1\mathrm{}\sigma _n\}`$, and multiplying by the number of possible choices of n nodes out of N:

$`G={\displaystyle \frac{aN^{ng_ag_m}}{K^{g_a}M^{g_m}}}{\displaystyle \underset{j=1}{\overset{n}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{K_{\sigma _j}}{k_j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{R_{\sigma _j}}{r_j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{M_{\sigma _j}}{m_j}}\right)_{\{\sigma \}}`$ (4)
Where $`K`$ is the average outdegree (equals the average indegree $`R`$), and $`M`$ is the average mutual edge degree. The symmetry factor $`a`$ is $`a_0^1_{j=1}^nk_j!r_j!m_j!`$ , where $`a_0`$ is the number of different permutations of the nodes that give an isomorphic subgraph.

The average (4) reduces to a product of moments of different orders of the indegree, outdegree and mutual degree distributions:

$$G=\frac{aN^{ng_ag_m}}{K^{g_a}M^{g_m}}\underset{j=1}{\overset{n}{}}\left(\genfrac{}{}{0pt}{}{K_i}{k_j}\right)\left(\genfrac{}{}{0pt}{}{R_i}{r_j}\right)\left(\genfrac{}{}{0pt}{}{M_i}{m_j}\right)_i$$
(5)
where the fact that each node should participate in the summation of only one term j introduces higher order corrections which we neglect. For example, subgraph id102 (Table I), has n=3 nodes, $`g_a=2`$ single edges and $`g_m=1`$ mutual edge. The subgraph degree sequences are $`k_j=\{1,1,0\}`$, $`r_j=\{0,1,1\}`$, and $`m_j=\{1,0,1\}`$. Using (5) we find :

$$G=id102=\frac{KMRMRK}{K^2M}$$
(6)
The approximation (Eq. 5) is exact in the case of Erdős networks. In Erdős networks, both indegree and outdegree are Poisson distributed and independent, and Eq. (5) reduces to Eq. (1) .
For non-sparse networks, a more accurate approximation takes into account the probabilities of a non-existent edge between two nodes (see Appendix B).
We tested the equations on random networks taken with the degree sequence of real world networks - transcription interactions in the yeast S. cerevisiae(Milo 2002 ), synaptic connections between neurons in C.elegans(White ) and world-wide-web hyperlinks between web pages in a single domain(Barabasi1999 ). When multiple edges in the same direction are allowed, as in the configuration model, the equations (5) are within a few percent of the numerical simulation results (Table I). We have also simulated random networks in which only one edge was allowed in each direction between any two nodes. As can be seen in Table I, the equations (5) are still within a few percent of the numerical simulation results for most subgraphs. There are some discrepancies (most notably a factor of almost 4 for subgraph id38 in the randomized world wide web network). In addition, we find good agreement between our approximation and numerical enumeration of subgraphs in simulated random networks with scale-free outdegree (Fig. 3).

## III Scale-free Networks

Scale-free networks have degree distributions that follow $`P(k)k^\gamma `$ at large $`k`$ (Barabasi1999 -Sole 2001 ). We consider directed networks in which the outgoing edge degree is scale-free, while the incoming edge degree distribution is Poissonian. Our results can be easily extended to scale-free indegree. For simplicity we choose the following form for the outgoing degree distribution for a network with N nodes (this function was used in Newman 2001  to fit world-wide web data) :

$`P(k)={\displaystyle \frac{\gamma 1}{k_0^{1\gamma }}}(k+k_0)^\gamma k<N`$ (7)
The mean connectivity $`K`$ is determined by $`k_0`$.

The hub is the node with the maximal number of outgoing edges, T. The hub size distribution (Fig. 4) is :

$`P(T)=NP(k=T)[P(kT)]^{N1}=`$
$`{\displaystyle \frac{N(\gamma 1)}{k_0}}(T/k_0)^\gamma \left(1(T/k_0)^{\gamma +1}\right)^{N1}`$ (8)
assuming $`Tk_0`$. For $`2<\gamma <3`$, the mean hub scales as:

$`T={\displaystyle _1^{N1}}TP(T)𝑑TN^{\frac{1}{\gamma 1}}`$ (9)
where the mean is over an ensemble of random networks with the same $`\gamma `$ and mean connectivity (see also dorogovtsev\_mesoscopic ,Cohen\_resilience  for an alternative method of deriving this result).

At $`\gamma 2`$, there is a condensation effect Bianconi 2001 , where a finite fraction of the nodes have outdegree $`1`$, and the mean hub size becomes proportional to N. Using (5), and assuming a compact distribution for the number of mutual edges, we find that the subgraph distribution is dominated by the hubs, and that the dominant term is that of the subgraph node with maximal outdegree, $`s`$. The number of appearances of each subgraph can be shown to scale as :

$`GaN^{ng1}K^{gs}{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{K_i}{s}}\right)N^\alpha `$ (10)
where $`g=g_a+2g_m`$ is the total number of edges in the subgraph Footnote1 . We derive the scaling exponent $`\alpha `$ in the following section.

## IV Transitions at different $`\gamma `$

The subgraph numbers scale as

$$GN^\alpha $$
(11)
We find three different regimes, in each of which the scaling exponent $`\alpha `$ behaves differently. Taking an ensemble average by integrating the largest term in Eq. (10) over the hub distribution (8) we get:

$$GN^{ng1}_1^{N1}T^sP(T)𝑑T$$
(12)
For $`\gamma 2`$ the network is in a condensed regime, where the hub $`T=O(N)`$. In this regime :

$`GN^{ng+s1}`$ (13)
For $`2<\gamma <\gamma _c`$ substituting Eq. 8 in Eq. 12 yields :

$`GN^{ng+s\gamma +1}`$ (14)
In this regime, the tail of P(T) is the dominant contribution to the integral. Finally at values above a critical $`\gamma `$, another transition occurs, where $`\alpha `$ equals the scaling exponent in Erdős networks, $`\alpha =ng=\alpha _{erdos}`$. The critical $`\gamma `$ is $`\gamma _c`$:

$`\gamma _c=s+1`$ (15)
In this regime, the hubs no longer contribute significantly to the subgraph distribution. In summary, $`GN^\alpha `$, where $`\alpha `$ is :

$`\alpha =\{\begin{array}{cc}ng+s1& \gamma 2\\ ng+s\gamma +1& 2<\gamma <s+1\\ ng& \gamma s+1\end{array}`$ (16)
Table 2 shows the expected scaling exponent for the 13 connected directed 3-node subgraphs, as well as for several 4-node subgraphs. The scaling laws agree very well with numerical results (Fig. 5). The three regimes of scaling are clearly seen. Note that the topology of each subgraph effects its scaling, through the subgraph maximal outdegree, $`s`$. These results can be easily extended to the case of scale-free indegree and non-directed networks. For loops of any size in non-directed networks the critical $`\gamma `$ is $`\gamma _c=3`$. At $`\gamma >3`$, loop numbers scale as $`N^0`$. This is consistent with bianconi\_hloops , which showed logarithmic corrections for the number of loops in Barabasi-Albert scale-free networks, which have $`\gamma =3`$.

## V Discussion

To summarize, we have presented an approximate solution for the average number of directed connected subgraphs in an ensemble of random networks with arbitrary degree sequence. We have presented scaling formulas for the number of subgraphs in scale-free random networks, and showed that the subgraph numbers can be very different from those in Erdős random networks. Whereas in Erdős random networks the scaling exponent is strictly determined by the number of nodes and edges of the subgraph, in scale-free random networks, the exact topology of the subgraph determines the scaling exponent. We showed that the scaling exponent $`\alpha `$ exhibits three different scaling laws in three regimes, depending on the control parameter $`\gamma `$ (the power of the degree distribution). In the common case of scale-free networks with $`\gamma `$ between 2 and 3, there are many more subgraphs which contain a node connected to more than one other node than in the corresponding Erdős networks with the same mean connectivity. For example, the feed-forward loop, (id38 in Table I) is much more common for $`\gamma <3`$. At $`\gamma =2.5`$, the number of feed-forward-loops scales as $`N^{0.5}`$, as opposed to $`N^0`$ in Erdős networks. On the other hand, subgraphs such as the 3-node cycle (id98 in Table I) have the same scaling, $`N^0`$, as in Erdős networks.

This study adds to our understanding of the random network models to which real-world networks are compared. It highlights the importance of using random networks that preserve the single and mutual degree sequence of the real network. Our approach may be readily extended to networks with multiple colors of edges. The present results may be useful for enumerating subgraphs in very large random networks which are beyond the reach of current numerical algorithms.

## Appendix A Edge Probabilities

Here we give a more detailed derivation for the edge probabilities used in Eq. (2,3). Without loss of generality we treat a network with no mutual edges. We denote by $`E=NK`$ the total number of edges. We begin by calculating the probability that no edge connects a source node with $`K`$ outgoing edges and a target node with $`R`$ incoming edges. This happens when all $`K`$ edges connect to a set of nodes $`\{\sigma _i\}_{i=1}^k`$ which does not contain the target node:

$$p(noedge|\{\sigma _i\})=\underset{k=0}{\overset{K1}{}}\left(1\frac{R}{ER^{}_{i=1}^kR_{\sigma _i}}\right)$$
(17)
where $`R^{}`$ is the indegree of the source node (we do not allow self edges). The probability of having no edge is obtained by summing over all possible sets $`\{\sigma _i\}_{i=1}^k`$ :

$$p(noedge)=\frac{1}{K!\left(\genfrac{}{}{0pt}{}{N2}{K}\right)}\underset{\{\sigma \}}{}\underset{k=0}{\overset{K1}{}}\left(1\frac{R}{ER^{}_{i=1}^kR_{\sigma _i}}\right)$$
(18)
Assuming $`\mathrm{max}_{i=1}^kR_{\sigma _i}E`$, and taking the complement as the probability of an edge existing, we obtain:

$$p(edge)=1(1\frac{R}{NK})^K=1e^{KR/NK}\frac{KR}{NK}$$
(19)
where our last approximation assumes $`KRNK`$. Intuitively, this result can be understood as $`K`$ attempts for the source node to connect to the target node with a probability of $`R/NK`$ at each attempt. $`R/NK`$ is the probability of an arbitrary edge connecting into the target node. Pairs of nodes in which $`KR`$ is of the order of $`NK`$ will contribute multiple edges in the same direction in the approximation, leading to over-estimation of subgraph numbers in the simulated networks where multiple edges are not allowed (Table I).

## Appendix B Non-sparse networks

In calculating the number of appearances of subgraphs in non-sparse networks, a more accurate approximation takes into account the probabilities of a non-existent edge between two nodes. For such subgraphs, in addition to the specified subgraph, Eq. (5) counts a set of subgraphs, with the null edges replaced by single or mutual edges. The corrections for the 3-node subgraphs are :

$`id6=id6id38id108`$
$`id12=id12id38id102`$
$`id14=id14id46id102id110`$
$`id36=id36id38id46`$
$`id74=id74id102id108id110`$
$`id78=id78id110id238`$
(20)
where $`G`$ represents the values obtained from Eq. (5), and $`G`$ is the corrected value. Generally, for larger subgraphs the corrections made will be of an inclusion-exclusion type.

## Appendix C Subgraph Enumeration

In numerically enumerating the subgraphs we combine a dynamic programming method(Milo 2002 ), which is applied generally for n-node subgraphs with $`n4`$, and a more rapid calculation, based on adjacency matrix operations, used for 3-node subgraphs. The method generalizes the results of (Harari ). Here we give formulas for the thirteen 3-node connected directed subgraphs based on the adjacency matrix. The network adjacency matrix is denoted by $`M`$, where $`M_{ij}=1`$ if a directed edge exists from node i to node j. We begin by dividing the network into a network containing only antisymmetric arrows, whose adjacency matrix will be denoted by $`A`$, and a network containing only mutual arrows, whose symmetric adjacency matrix will be denoted as $`S`$ .

$`M=A+S`$ (21)
We denote by $`AB`$ the matrix multiplication of matrices $`A`$ and $`B`$, and by $`AB`$ the dot multiplication. $`\stackrel{~}{A}`$ is the logical inverse of matrix $`A`$, where the 0 elements of $`A`$ are the 1 of $`\stackrel{~}{A}`$ and vice-versa. $`A^{}`$ is the transpose matrix of $`A`$. A summation denotes summation of all the matrix indices. The matrix formulas for the 13 directed connected 3-node subgraphs are given in Table 3. For example id38 has two nodes which are connected by a path of 2 edges and a path of one edge. $`A_{}^{2}{}_{ij}{}^{}`$ is the number of length 2 paths between node i and node j. Dot-multiplication with matrix $`A`$ and summation of the terms of the resultant matrix gives the correct count. In some of the subgraphs a correction is made for the terms on the diagonal (id6,id36,id78).

###### Acknowledgements.

We thank S. Maslov, R. Cohen, A. Mayo, A. Natan, M. Itzkovitz and all members of our lab for valuable discussions. We acknowledge support from the Israel Science Foundation, the Human Frontier Science Program, and the Minerva Foundation.