text:
# Optical generation of quasisolitons in a single–component gas of neutral fermionic atoms

## I Introduction

Recent pioneering work by DeMarco and Jin Jin  on realization of quantum degeneracy in a trapped gas of fermionic 40K atoms has triggered the interest in properties of cold fermions. As opposed to bosons, fermions cannot occupy the same quantum state, i.e., they obey the Pauli exclusion principle. Since at low temperatures only the $`s`$–wave scattering gets important, the spin–polarized fermions stop to collide at lower temperatures and rethermalization process necessary for evaporation breaks down. This difficulty has been overcome by trapping and cooling the fermions in two hyperfine states, enabling in this way the evaporation by collisions between atoms in a different spin states Jin ; Jin1 . Other groups have achieved quantum degeneracy in a gas of 6Li atoms by sympathetic cooling of a mixture of bosonic (7Li) and fermionic (6Li) atoms reaching the temperatures down to $`0.2T_F`$ Hulet ; Salomon .

Despite of experimental difficulties in achieving the strong quantum degeneracy in fermionic gases, some theoretical work on degenerate Fermi gases in the normal phase at zero and finite temperatures has already been done. The static properties of one– and two–component system were analysed in Refs. Rokhsar , Bruun1  and Roth . Collective excitations in degenerate Fermi gases in the normal phase were discussed within the hydrodynamic approximation in Ref. Bruun  and on the basis of the sum rules in Ref. Vichi . In the absence of $`s`$–wave collisions other forces, like for instance dipole–dipole interaction, start to play role. The fermionic dipoles have been investigated recently, including the analysis of stability conditions Krzysiek . Another aspect of dynamic behavior of a degenerate Fermi gas is studied in Ref. Tomek  where the question about the possibility of generation of solitons and vortices in a normal state of non–interacting Fermi gas is asked. It is shown by solving the many–body Schrödinger equation that the phase imprinting method is capable of generating the solitons and vortices in such a system. An alternative method of generating such structures is discussed in Ref. Damski .

So far, the phase imprinting method has been successfully applied to generate solitons in Bose–Einstein condensates Hannover ; NIST . Although not realized experimentally yet, it offers also a rather rare opportunity of control of the generation of vortices Dobrek ; Andrelczyk . The method consists in passing an off–resonant laser pulse through the appropriately designed absorption plate followed by the impinging it on the atomic gas. The laser pulse is short and strong enough so that the atomic motion is negligible during the pulse and the basic effect of light is a generation of strong optical potential $`V_{opt}(\stackrel{}{r},t)`$ due to the ac–Stark effect. This potential acts for a short time $`\tau _{dur}`$, and essentially imprints a desired phase of atomic wave functions $`\phi (\stackrel{}{r})=_0^{\tau _{dur}}V_{opt}(\stackrel{}{r},t)𝑑t/\mathrm{}`$. After this, the density and the phase try to adopt to each other, which results in generation of structures like solitons and vortices, depending on the pattern written in the absorption plate.

In this paper we explore further the possibility of generating the quasisolitons in a gas of non–interacting Fermi atoms. According to the common knowledge regarding solitons, the nonlinearity of the problem plays essential role, since it cancels the dispersion and allows for a propagation of shape preserving pulses. So, at the first glance it is surprising that solitons can exist in a system of non–interacting particles described by the many–body Schrödinger equation, i.e., the linear equation. However, we show that, due to the Fermi statistics, there exists the regime of parameters where soliton–like solutions can exist in the non–interacting Fermi gas.

The paper is organized as follows. In Sec. II we solve the many–body Schrödinger equation for a system of non–interacting fermions undergoing the phase imprinting, leading to the generation of soliton–like structures. Both zero and finite temperatures are considered. Sec. III offers an analytical approach to the problem of propagation and decay of solitons, based on the Wigner function formalism, for a degenerate as well as a Boltzmann gas. In Sec. IV another approach based on analysis of density matrices, is followed in order to discuss the Thomas–Fermi approximation and to understand the role of dimensionality. Finally, we conclude in Sec. V.

## II Solution of the many–body Schrödinger equation

### II.1 Zero temperature

At zero temperature the many–body wave function of a system of non–interacting particles (fermions) is given by the Slater determinant with the lowest available one–particle orbitals occupied.

$`\mathrm{\Psi }(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)={\displaystyle \frac{1}{\sqrt{N!}}}\left|\begin{array}{ccccc}\phi _1(\stackrel{}{r}_1)\hfill & .\hfill & .\hfill & .\hfill & \phi _1(\stackrel{}{r}_N)\hfill \\ .\hfill & & & & .\hfill \\ .\hfill & & & & .\hfill \\ .\hfill & & & & .\hfill \\ \phi _N(\stackrel{}{r}_1)\hfill & .\hfill & .\hfill & .\hfill & \phi _N(\stackrel{}{r}_N)\hfill \end{array}\right|.`$ (6)
The one–particle orbitals $`(\phi _1(\stackrel{}{r}),\mathrm{},\phi _N(\stackrel{}{r}))`$ are orthogonal and they undergo unitary evolution. Hence, their orthogonality is sustained, both during the phase imprinting and afterwords. Therefore, the one–particle density matrix is at any time given by the following formula

$$\rho _1(\stackrel{}{r}^{},\stackrel{}{r}^{\prime \prime },t)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\phi _i(\stackrel{}{r}^{},t)\phi _i^{}(\stackrel{}{r}^{\prime \prime },t),$$
(7)
and its diagonal part is just the particle density.

Evolution of each single–particle orbital is obtained from the one–particle Schrödinger equation and can be split into two stages. First, the phase imprinting technique is used to disturb the system in a way of writing a desired phase on it. Assuming fast enough phase imprinting, the resulting wave function has the form (this procedure can be also performed numerically within the finite duration of the imprinting process):

$`\phi _k(\stackrel{}{r})\phi _k^{}(\stackrel{}{r})=\phi _k(\stackrel{}{r})\mathrm{exp}(i\varphi (x)),`$
for $`k=1,\mathrm{},N`$, where $`\varphi (x)`$ is the phase imprinted on each atom. Throughout this paper we deal with the generation of quasisolitons, we assume that the imprinted phase depends only on a single coordinate, say ’$`x`$’, and utilize $`\varphi (x)(1+\mathrm{tanh}(x/d_{\mathrm{phase}}))`$. We consider a Fermi gas in a three–dimensional box or harmonic trap characterized by the frequencies $`\omega _x`$, $`\omega _y`$, and $`\omega _z`$. Both are separable in Cartesian coordinates. The one–particle orbitals $`\phi _k(\stackrel{}{r})`$ are taken as a product $`\phi _{n_x}^{(1)}(x)\phi _{n_y}^{(2)}(y)\phi _{n_z}^{(3)}(z)`$ of eigenvectors of the Hamiltonians of one–dimensional boxes or harmonic oscillators.

In the second stage the system undergoes the free evolution in three–dimensional space. However, because the Hamiltonian separates in coordinates ’$`x`$’, ’$`y`$’, and ’$`z`$’, the time propagation of each orbital can be easily reduced to one–dimensional case:

$`\phi _k(x,y,z,t)`$ $`=`$ $`e^{i\mathrm{E}_{\mathrm{n}_\mathrm{y}}\mathrm{t}/\mathrm{}}e^{i\mathrm{E}_{\mathrm{n}_\mathrm{z}}\mathrm{t}/\mathrm{}}`$
$`\times `$ $`\phi _{n_y}^{(2)}(y,0)\phi _{n_z}^{(3)}(z,0)`$
$`\times `$ $`\phi _{n_x}^{(1)}(x,t)`$
where

$`\phi _{n_x}^{(1)}(x,t)`$ $`=`$ $`e^{i\mathrm{H}_\mathrm{x}\mathrm{t}/\mathrm{}}\phi _{n_z}^{(1)}(x,0)e^{i\varphi (x)}`$
and $`\mathrm{H}_\mathrm{x}`$ is the Hamiltonian of one–dimensional box or harmonic oscillator. The diagonal part of one–particle density matrix is given then by the expression:

$`\rho (\stackrel{}{r},t)={\displaystyle \frac{1}{\mathrm{N}}}{\displaystyle \underset{n_x,n_y,n_z}{}}|\phi _{n_y}^{(2)}(y,0)|^2|\phi _{n_z}^{(3)}(z,0)|^2`$
$`\times |\phi _{n_x}^{(1)}(x,t)|^2.`$ (8)
In Fig. 1 we plotted the density profiles of a uniform system of $`500`$ atoms after imprinting the phase designed in the form $`\varphi (x)=\varphi _0(1+\mathrm{tanh}(x/d_{\mathrm{phase}}))/2`$. The characteristic length (the Fermi length) is given by the formula $`\lambda _F=h/p_F=2/NL`$, where $`p_F`$ is the Fermi momentum and $`L`$ is the size of one–dimensional box. For the case of Fig. 1, the Fermi length equals $`0.004L`$ and is $`2.5`$ times smaller than the width characterising the jump of the imprinted phase. As discussed more precisely in Sec. III, the evolution of the density after phase imprinting depends strongly on the ratio $`d_{\mathrm{phase}}/\lambda _F`$ between the width of the phase step and the Fermi length.

It is clear from Fig. 1 (as well as Fig. 2) that two quasisolitons (the bright and the dark one) are generated. They propagate in opposite directions with distinct velocities. We present in Table 1 the values of soliton’s velocities for various widths of the phase steps. The speed of the bright quasisoliton is always higher than the speed of sound whereas dark quasisolitons move slower than the sound. However, increasing the ratio $`d_{\mathrm{phase}}/\lambda _F`$ both velocities are getting closer to the speed of sound ($`p_F/m`$ in one–dimensional space), reaching eventually (i.e., in the limit $`d_{\mathrm{phase}}\lambda _F`$) the value $`p_F/m`$ as shown in Secs. III and IV. In Sec. IV we also analyze the $`d_{\mathrm{phase}}\lambda _F`$ limit by decreasing $`\lambda _F`$ (increasing the number of atoms) and solving the equations of the Thomas–Fermi model.

Fig. 2 shows that going deeper into the regime $`d_{\mathrm{phase}}\lambda _F`$ the soliton structures become less pronounced. This again is explained in terms of Wigner function in Sec. III and based on Thomas–Fermi approach in Sec. IV. The condition $`d_{\mathrm{phase}}\lambda _F`$ means that the height (or depth) of solitons is getting smaller in comparison with unperturbed density. Going in opposite direction (i.e., $`d_{\mathrm{phase}}\lambda _F`$) results, as presented in Fig. 3, in a change of shape of the structures (dispersion) on a time scale of the order of $`\mathrm{}/\epsilon _F`$.

Another aspect of Fig. 1 is connected with the Thomas–Fermi approximation discussed in Sec. IV. Each frame of Fig. 1 shows, in fact, two curves; one of them (solid line) is obtained by solving the many–body Schrödinger equation, and the second (dashed line) is representing the Thomas–Fermi approach. Both curves match well which proves the validity of the Thomas–Fermi approximation in one–dimensional space for large enough number of atoms just as it happens in the case of a harmonic trap (see Fig. 1 in Ref. Tomek , where two– and three–dimensional cases are also discussed).

### II.2 Finite temperatures

One way of including temperature effects on the soliton’s generation and propagation is to allow for populating the single–particle states of higher energy than the Fermi energy. We have realized this by considering a grand canonical ensemble for a system of non–interacting fermions. We generate a number of many–particle configurations according to the Fermi–Dirac statistics, i.e., with the probability of populating a particular one–particle state $`\phi _k`$ given by $`\mathrm{exp}(\beta (\epsilon _k\mu ))+1)^1`$, where $`\epsilon _k`$ is the state energy, $`\beta `$ determines the bath temperature (via $`\beta =1/k_BT`$), and $`\mu `$ is the chemical potential (dependent on temperature in general). More precisely, for each one–particle state (unless the probability of populating this state is lower than $`10^6`$) we generate a random number from the interval $`[0,1]`$. If the number is less than $`\mathrm{exp}(\beta (\epsilon _k\mu ))+1)^1`$, the one–particle state is populated, otherwise it is left empty. Here, $`\epsilon _k`$ is the energy of the state under consideration, the chemical potential

$`\mu =\epsilon _F[1(\pi ^2/3)(k_BT/\epsilon _F)^2],`$ (9)
and $`\epsilon _F`$ equals the Fermi energy corresponding to the average number of atoms. The above formula for the chemical potential applies for the Fermi gas trapped in a harmonic potential for the temperatures up to $`0.55T_F`$ Rokhsar .

In such a way, having given the temperature and the average number of atoms, the many–particle configurations are sampled. Then, for each many–particle configuration the phase imprinting is performed followed by a free evolution. Finally, the density is calculated as an average according to a grand canonical ensemble rules. We have investigated in this way the propagation of quasisolitons for temperatures up to $`0.3T_F`$, already attainable in experiment Jin1 ; Hulet ; Salomon .

In Fig. 4 we plot evolution of the averaged density of the Fermi gas confined in a one–dimensional harmonic trap of frequency $`\omega `$ after phase imprinting of a single step of $`2\pi `$ and the width equal to $`0.5`$ osc. units. ($`\sqrt{\mathrm{}/(m\omega )}`$). Here, the chemical potential is given by the formula $`\mu =\epsilon _F+k_BT\mathrm{log}[1\mathrm{exp}(\epsilon _F/k_BT)]`$, which is a one–dimensional counterpart of formula (9), the temperature equals $`0.2T_F`$, and the averaging procedure is performed over $`5000`$ many–body configurations. Both, bright and dark solitons are still present in the system exhibiting basically the same properties as at zero temperature. For higher temperatures (already $`0.3T_F`$) the soliton structures start to get hidden within the noise originating in the averaging procedure (Fig. 5). However, one has to remember that rather a particular many–body configuration is realized in a single experiment, and not an ensemble of them. Finally, Fig. 6 gives an example that going from one– to three-dimensional space, predictions made based on lower dimensionality survive. Fig. 6 shows density profiles of fermionic gas trapped in a spherically symmetric potential (averaged over $`1000`$ configurations at $`0.1T_F`$ temperature) and can be compared with corresponding figure from Ref. Tomek  at zero temperature since the details of phase imprinting are the same.

## III Wigner function approach

The Wigner function $`w(x,p,t)`$ is defined as

$$w(x,p,t)=𝑑y\widehat{\psi }^{}(x+y/2,t)\widehat{\psi }(xy/2,t)\mathrm{exp}(ipy),$$
(10)
where $`\widehat{\psi }^{}(x,t)`$ and $`\widehat{\psi }(x,t)`$ denote atomic creation and annihilation operators in the Heisenberg picture, respectively. It obeys the equation

$$\left(_t+\frac{p}{m}_x_xV_{ext}(x)_p\right)w(x,p,t)=0,$$
(11)
where $`V_{ext}(x)`$ is the external harmonic trap potential. It immediately follows from the definition, Eq. (10), that the density of particles at time $`t`$ is

$$\rho (x,t)=\frac{dp}{2\pi \mathrm{}}w(x,p,t).$$
(12)
The phase imprinting at $`t=0`$ modifies the Wigner function in the following way:

$`w_0(x,p,t=0)`$
$`w_0(x,p)={\displaystyle 𝑑y\widehat{\psi }^{}(x+y/2,t=0)\widehat{\psi }(xy/2,t=0)}`$
$`\times \mathrm{exp}(ipy/\mathrm{}i\varphi (x+y/2)+i\varphi (xy/2)),`$
where $`\varphi (x)`$ is the imprinted phase. The function $`w_0(x,p)`$ provides the initial condition for the dynamical evolution according to Eq. (11).

Let us consider first the spatially homogeneous case ($`V_{ext}(x)=0`$). Eq. (11) with the initial condition $`w(x,p,t=0)=w_0(x,p)`$ can be easily solved. The solution is

$$w(x,p,t)=w_0(x(p/m)t,p),$$
(14)
and reflects the ballistic motion of particles in the absence of an external potential. By using the fact that in the spatially homogeneous case the correlation function $`\widehat{\psi }^{}(x+y/2,t=0)\widehat{\psi }(xy/2,t=0)=\rho _0(y)`$, i.e., it depends only on ’$`y`$’, the expression for $`w_0(x,p)`$ can be rewritten in the form

$`w_0(x,p)={\displaystyle }dy\rho _0(y)\mathrm{exp}(ipy/\mathrm{}i\varphi (x+y/2)`$
$`+i\varphi (xy/2)).`$ (15)
The above formulae have a general character, and are valid for both degenerate and non-degenerate (Boltzmann) Fermi gas. The difference between these two cases is only in the specific form of the function $`\rho _0(y)`$.

Usually in the experiments, the imprinted phase $`\varphi (x)`$ is a function with transition regions of the size of the order of $`d_{\mathrm{phase}}`$ between the regions with constant values of the phase. Depending on the ratio between $`d_{\mathrm{phase}}`$ and the characteristic width of the function $`\rho _0(y)`$, one has completely different behavior of the solution Eq. (14) and, therefore, the density evolution. In the degenerate regime, the characteristic width of $`\rho _0`$ is given by the Fermi wavelength $`\lambda _F`$, which is inversely proportional to the Fermi momentum, $`\lambda _F=h/p_F`$. In the opposite case of a Boltzmann gas, the width is given by the thermal wavelength $`\lambda _T=h/\sqrt{2mT}`$.

We first analyze the case where $`d_{\mathrm{phase}}\lambda _F,\lambda _T`$. Under this condition, one can expand the exponent in Eq. (15) in powers of $`y`$, and leave only the lowest, linear in $`y`$, contribution:

$`w_0(x,p)`$ $``$ $`{\displaystyle 𝑑y\rho _0(y)\mathrm{exp}(ipy/\mathrm{}i\varphi ^{}(x)y)}`$ (16)
$`=`$ $`\stackrel{~}{\rho }_0(p+\varphi ^{}(x)),`$
where $`\stackrel{~}{\rho }_0(p)`$ is the Fourier transform of the function $`\rho _0(y)`$ and $`\varphi ^{}`$ is the derivative of $`\varphi `$. The density of the particles, Eq. (12), can then be written as

$`\rho (x,t)`$ $`=`$ $`{\displaystyle \frac{dp}{2\pi \mathrm{}}w_0(x(p/m)t,p)}`$ (17)
$``$ $`{\displaystyle \frac{dp}{2\pi \mathrm{}}\stackrel{~}{\rho }_0(p+\mathrm{}\varphi ^{}(x(p/m)t))}.`$
Keeping in mind that $`\varphi ^{}\pi /d_{\mathrm{phase}}`$ and, therefore, is much smaller than the characteristic momenta, $`p_F`$ and $`p_T\sqrt{mT}`$ for the degenerate and Boltzmann gas, respectively, Eq. (17) can be simplified as:

$$\rho (x,t)\rho _0+\frac{dp}{2\pi }\stackrel{~}{\rho }_0^{}(p)\varphi ^{}(x(p/m)t),$$
(18)
where $`\rho _0=(dp/(2\pi \mathrm{}))\stackrel{~}{\rho }_0(p)`$ is the initial constant density and $`\stackrel{~}{\rho }_0^{}`$ is the derivative of $`\stackrel{~}{\rho }_0`$. The function $`\varphi ^{}(x)`$ is only nonzero in the transition regions of the width $`d_{\mathrm{phase}}`$. Therefore, Eq. (18) has a very clear meaning. It represents a superposition of bumps moving with different velocities $`p/m`$ which are distributed with the weights $`\stackrel{~}{\rho }_0^{}(p)`$.

For a degenerate Fermi gas at zero temperature one has $`\stackrel{~}{\rho }_0^{}(p)=\delta (p+p_F)\delta (pp_F)`$, and, hence,

$$\rho (x,t)\rho _0+\frac{1}{2\pi }[\varphi ^{}(x+(p_F/m)t)\varphi ^{}(x(p_F/m)t)].$$
(19)
The above formula describes a ”soliton-like” motion (without broadening and interaction) of bright and dark structures in the density, described by the function $`\varphi ^{}(x)`$, in the opposite directions with the same velocity $`v_F=p_F/m`$. The physical explanation of this phenomenon is very simple. Under the condition $`d_{\mathrm{phase}}\lambda _F`$, the quantum state of the system after the phase imprinting contains excitations, particle-hole pairs, only in the vicinity of the Fermi surface. (The energy width of the excitation distribution is of the order of $`\mathrm{}/d_{\mathrm{phase}}p_F`$.) As a result, they have the same velocity. either $`v_F`$ or $`v_F`$, and, therefore, the corresponding density profile moves in both directions with velocities $`\pm v_F`$, respectively, without any distortion. To be precise, the broadening in this case is of the second order in $`d_{\mathrm{phase}}/\lambda _F`$. To determine the characteristic time scale $`\tau _{F0}`$, one has to keep in mind that the width of the distribution of the excitation momenta after the phase imprinting is of the order of $`\delta p\mathrm{}/d_{\mathrm{phase}}p_F`$. Hence, the time $`\tau _{F0}`$ can be estimated as $`\tau _{F0}(\mathrm{}/T_F)(d_{\mathrm{phase}}/\lambda _F)^2\mathrm{}/T_F`$.

At temperatures $`TT_F`$ (but $`T>T_F\lambda _F/d_{\mathrm{phase}}`$, the condition when the thermal width of the excitation distribution is larger then that after the phase imprinting), the situation is very similar, and Eq. (18) can be written in the form

$`\rho (x,t)`$ $``$ $`\rho _0+{\displaystyle }{\displaystyle \frac{ds}{2\pi }}f_F(s,T)[\varphi ^{}(x+t(p_F+s)/m))`$ (20)
$`\varphi ^{}(xt(p_F+s)/m)],`$
where $`f_F^{}(s,T)=(p_F/4mT)\mathrm{cosh}^2(sp_F/2mT)`$ is related to the derivative of the Fermi-Dirac distribution and describes broadening of the $`\delta `$-like peaks at $`p=\pm p_F`$ in the distribution $`\stackrel{~}{\rho }_0^{}(p)`$ at finite temperatures. As a result, now one has a superposition of bright bumps and dark wells in the density profile, moving with slightly different velocities distributed with two narrow peaks at $`\pm p_F/m`$ with the width of the order of $`mT/p_Fp_FT/T_Fp_F`$. This results in spreading of bumps and wells in the density profile on a time scale $`\tau _{FT}(\mathrm{}/T)(d_{\mathrm{phase}}/\lambda _F)(\mathrm{}/T_F)(T_F/T)(d_{\mathrm{phase}}/\lambda _F)\mathrm{}/T,\mathrm{}/T_F`$ (at $`T=T_F\lambda _F/d_{\mathrm{phase}}`$ one has $`\tau _{FT}=\tau _{F0}`$). Therefore, the time evolution of the gas density after a phase imprinting can be described in terms of a combination of several quasisolitons (bumps) and quasi–darksolitons (wells) which move through each other without any interaction and, for times $`t\tau `$, decay (broaden) very slowly. In Fig. 7 we present the results of numerical calculations of the time evolution of the density and Wigner functions for a spatially homogeneous gas at the temperature $`T=0.01T_F`$ after the $`2\pi `$–step phase imprinting with $`d_{\mathrm{phase}}=160\lambda _F`$. The results clearly demonstrate (see Fig. 7a) that under these conditions the soliton–like structures caused by the phase imprinting behave like real solitons, i.e., they propagate without any noticeable decay (broadening) on a long time scale. The blue regions in Figs. 7b and 7c correspond to negative values of the Wigner function, and, therefore, indicate the quantum origin of the phenomenon.

For high temperatures $`TT_F`$ gas is in the Boltzmann regime, and Eq. (18) takes the form

$$\rho (x,t)\rho _0+\frac{dp}{2\pi }f_B(p,T)\varphi ^{}(x(p/m)t),$$
(21)
where $`f_B(p,t)=\rho _0\sqrt{2\pi /mT}(p/mT)\mathrm{exp}(p^2/2mT)`$ describes the distribution with two broad peaks at $`p=\pm p_T`$ with $`p_T\sqrt{mT}`$, with the width $`\delta pp_T`$. In this case, the decay of the original density profile takes place more rapidly than in the case of a degenerate gas, on the time scale $`\tau _B(\mathrm{}/T)(d_{\mathrm{phase}}/\lambda _T)\tau _F\sqrt{T_F/T}\tau _F`$. The corresponding numerical results for the density evolution at the temperature $`T=5T_F`$ are shown in Fig. 8. Although the parameters of the phase imprinting are exactly the same as for the case of Fig. 7, the temperature (or, in other words, the absence of the sharp Fermi surface) cause a rapid decay of soliton-like structures.

The same consideration is applied to the trapped gas in the harmonic potential. In this case, due to the periodicity of the motion, the essential scale for the time evolution is set by the inverse level spacing: $`t2\pi /\omega `$, where $`\omega `$ is the trap frequency. Therefore, after phase imprinting one would observe a quasisoliton behavior if $`2\pi /\omega \tau _{F0},\tau _{FT}`$, or $`2\pi /\omega \tau _B`$ for the degenerate or Boltzmann regime, respectively. More explicitly, in the degenerate regime one has the following requirement for the temperature $`T`$ and the width $`d_{\mathrm{phase}}`$: $`[\mathrm{}\omega /\mathrm{max}(T,T_F\lambda _F/d_{\mathrm{phase}})](d_{\mathrm{phase}}/\lambda _F)1`$, or $`d_{\mathrm{phase}}/l_0\mathrm{max}(1,T/\mathrm{}\omega \sqrt{N})`$, where $`l_0=\sqrt{\mathrm{}/m\omega }`$ is the oscillator length and $`N`$ is the number of particle. Correspondingly, in the regime of a Boltzmann gas, the requirement reads: $`(\mathrm{}\omega /T)(d_{\mathrm{phase}}/\lambda _T)1`$, or $`d_{\mathrm{phase}}/l_0\sqrt{T/\mathrm{}\omega }`$.

In the regime where $`d_{\mathrm{phase}}\lambda _F,\lambda _T`$, the excitations created during the phase imprinting have a very broad distribution with the width $`\delta p\mathrm{}/d_{\mathrm{phase}}p_F,p_T`$ in a spatially homogeneous gas, and $`\delta \epsilon \mathrm{}\omega (l_0/d_{\mathrm{phase}})^2`$ in a trapped one. In this case one would expect a very rapid decay of ”solitons” on a time scale $`\tau [\mathrm{}/\mathrm{max}(T_F,T)][d_{\mathrm{phase}}/\mathrm{min}(\lambda _F,\lambda _T)]^2\mathrm{}/T_F`$ for both spatially homogeneous and trapped Fermi gases. The corresponding numerical results for the density and Wigner function evolution after the $`2\pi `$–step phase imprinting with $`d_{\mathrm{phase}}=2\lambda _F`$ are shown in Figs. 9 (for the temperature $`T=0.01T_F`$) and 10 (for the temperature $`T=T_F`$). In the quantum degenerate case (Fig. 9), the Wigner function has a well pronounced fringes structure, and the corresponding complicated structure in the density profile decays rather fast. For the non–degenerate gas (Fig. 10), the Wigner function is smooth everywhere except for the narrow region subjected by the phase imprinting. In this case, the density profile has no small-scale structures, as compared to the degenerate case from Fig. 9, but also relaxes very fast to the uniform density distribution.

## IV Density matrix formalism

Another way to discuss the soliton’s formation in a gas of neutral fermionic atoms is to start from a set of equations for reduced density matrices. In particular, the equation of motion for the one–particle density matrix involves the two–particle density matrix $`\rho _2`$ and is given by equ

$`i\mathrm{}{\displaystyle \frac{}{t}}\rho _1(\stackrel{}{r}_1,\stackrel{}{r}_2,t)={\displaystyle \frac{\mathrm{}^2}{2m}}(\stackrel{}{}_1^2\stackrel{}{}_2^2)\rho _1(\stackrel{}{r}_1,\stackrel{}{r}_2,t)`$
$`+{\displaystyle d^3r^{}\left[V(\stackrel{}{r}_1\stackrel{}{r}^{})V(\stackrel{}{r}_2\stackrel{}{r}^{})\right]\rho _2(\stackrel{}{r}_1,\stackrel{}{r}^{};\stackrel{}{r}_2,\stackrel{}{r}^{},t)}`$
$`+\left[V_{ext}(\stackrel{}{r}_1,t)V_{ext}(\stackrel{}{r}_2,t)\right]\rho _1(\stackrel{}{r}_1,\stackrel{}{r}_2,t),`$ (22)
where $`V(\stackrel{}{r}_1\stackrel{}{r}_2)`$ is the two–particle interaction term and $`V_{ext}(\stackrel{}{r},t)`$ is the external potential. In the limit $`\stackrel{}{r}_1\stackrel{}{r}_2`$ Eq. (22) leads to the continuity equation

$`{\displaystyle \frac{\rho (\stackrel{}{r},t)}{t}}+\stackrel{}{}\left[\rho (\stackrel{}{r},t)\stackrel{}{v}(\stackrel{}{r},t)\right]`$ $`=`$ $`0,`$ (23)
where the density and velocity fields are defined as follows:

$`\rho (\stackrel{}{r},t)`$ $`=`$ $`\underset{\stackrel{}{r}_1\stackrel{}{r}_2}{lim}\rho _1(\stackrel{}{r}_1,\stackrel{}{r}_2,t)`$
$`\stackrel{}{v}(\stackrel{}{r},t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2m}}\underset{\stackrel{}{r}_1\stackrel{}{r}_2}{lim}(\stackrel{}{}_1\stackrel{}{}_2)\chi _1(\stackrel{}{r}_1,\stackrel{}{r}_2,t)`$ (24)
and $`\chi _1(\stackrel{}{r}_1,\stackrel{}{r}_2,t)`$ is the phase of the one–particle density matrix.

One can also rewrite Eq. (22) introducing the center–of–mass ($`\stackrel{}{R}=(\stackrel{}{r}_1+\stackrel{}{r}_2)/2`$) and the relative position ($`\stackrel{}{s}=\stackrel{}{r}_1\stackrel{}{r}_2`$) coordinates. By taking the derivative of Eq. (22) with respect to the coordinate $`\stackrel{}{s}`$ the hydrodynamic Euler–type equation of motion is obtained in the limit $`\stackrel{}{s}0`$

$`{\displaystyle \frac{\stackrel{}{v}(\stackrel{}{r},t)}{t}}`$ $`=`$ $`{\displaystyle \frac{\stackrel{}{}\mathrm{T}}{m\rho (\stackrel{}{r},t)}}[\stackrel{}{v}(\stackrel{}{r},t)\stackrel{}{}]\stackrel{}{v}(\stackrel{}{r},t)`$ (25)
$`+`$ $`{\displaystyle \frac{\stackrel{}{F}_{int}(\stackrel{}{r},t)}{m\rho (\stackrel{}{r},t)}}{\displaystyle \frac{\stackrel{}{}V_{ext}(\stackrel{}{r},t)}{m}}.`$
Here, the kinetic–energy stress tensor $`\mathrm{T}`$, whose elements are given by

$`\mathrm{T}_{kl}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{m}}\underset{\stackrel{}{s}0}{lim}{\displaystyle \frac{^{\mathrm{\hspace{0.17em}2}}\sigma _1(\stackrel{}{R},\stackrel{}{s},t)}{s_ks_l}},`$ (26)
depends on $`\sigma _1(\stackrel{}{r}_1,\stackrel{}{r}_2,t)`$, the amplitude of one–particle density matrix and can be calculated based on a local equilibrium assumption which is the substance of Thomas–Fermi approximation. $`\stackrel{}{F}_{int}(\stackrel{}{r},t)`$, on the other hand, represents the force due to fermion–fermion interaction (for spin–polarized atomic fermions the dipole–dipole interaction is a good candidate) and is defined as

$`\stackrel{}{F}_{int}(\stackrel{}{r},t)`$ $`=`$ $`{\displaystyle d^3r^{}\stackrel{}{}_\stackrel{}{r}V(\stackrel{}{r}\stackrel{}{r}^{})\rho _2(\stackrel{}{r},\stackrel{}{r}^{};\stackrel{}{r},\stackrel{}{r}^{},t)}.`$
The one–particle Wigner function within the Thomas–Fermi approximation is given by

$`w(\stackrel{}{r},\stackrel{}{p})`$ $`=`$ $`\eta (\mathrm{}^2[6\pi ^2\rho (\stackrel{}{r})]^{2/3}\stackrel{}{p}^{\mathrm{\hspace{0.33em}2}}),`$
where $`\eta ()`$ is Heaviside’s unit step function and hence the one–particle density matrix is calculated as

$`\rho _1(\stackrel{}{R},\stackrel{}{s})`$ $`=`$ $`{\displaystyle \frac{d^3p}{(2\pi \mathrm{})^3}w(\stackrel{}{R},\stackrel{}{p})e^{i\stackrel{}{p}\stackrel{}{s}/\mathrm{}}}.`$
After straightforward algebra one obtains:

$`\rho _1(\stackrel{}{R},\stackrel{}{s})`$ $`=`$ $`{\displaystyle \frac{2}{(2\pi \mathrm{})^2}}{\displaystyle \frac{1}{s}}\left[{\displaystyle \frac{\mathrm{}p_F}{s}}\mathrm{cos}\left({\displaystyle \frac{p_Fs}{\mathrm{}}}\right)+{\displaystyle \frac{\mathrm{}^2}{s^2}}\mathrm{sin}\left({\displaystyle \frac{p_Fs}{\mathrm{}}}\right)\right].`$
The one–particle density matrix is a real function in this case and depends only on the length $`s=|\stackrel{}{s}|`$ of vector $`\stackrel{}{s}`$. In the limit $`\stackrel{}{s}0`$ the kinetic–energy stress tensor $`\mathrm{T}`$ is getting diagonal with $`(\mathrm{}^2/m)(1/30/\pi ^2)(6\pi ^2\rho )^{5/3}`$ term at each position. Then, assuming the case of non–interacting fermions (spin–polarized atoms at low temperature), the Eqs. (23) and (25) become a closed set of equations for the density and velocity fields:

$`{\displaystyle \frac{\rho }{t}}+\stackrel{}{}\left(\rho \stackrel{}{v}\right)=0`$
$`{\displaystyle \frac{\stackrel{}{v}}{t}}+(\stackrel{}{v}\stackrel{}{})\stackrel{}{v}+\stackrel{}{}\left({\displaystyle \frac{\mathrm{}^2}{2m^2}}(6\pi ^2)^{3/2}\rho ^{2/3}+{\displaystyle \frac{V_{\mathrm{ext}}}{m}}\right)=0,`$
Note that the Eqs. (23) and (25) are general and could be used also for a bosonic system. For example, for the system of bosons described by the many–body wave function $`\mathrm{\Psi }(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)=\phi (\stackrel{}{r}_1)\mathrm{}\phi (\stackrel{}{r}_N)`$ (as in the case of the Bose–Einstein condensate in the mean–field approximation) the one–particle density matrix is given by

$`\rho _1(\stackrel{}{R},\stackrel{}{s})`$ $`=`$ $`\phi (\stackrel{}{R}+{\displaystyle \frac{1}{2}}\stackrel{}{s})\phi ^{}(\stackrel{}{R}{\displaystyle \frac{1}{2}}\stackrel{}{s}).`$
Now, the kinetic–energy stress tensor possesses off–diagonal elements and leads to the so–called “quantum pressure” term in the equation of motion. Assuming contact interaction between the bosonic atoms dominates, the appropriate equations are written as

$`{\displaystyle \frac{\rho }{t}}+\stackrel{}{}\left(\rho \stackrel{}{v}\right)=0`$
$`{\displaystyle \frac{\stackrel{}{v}}{t}}+\stackrel{}{}\left({\displaystyle \frac{4\pi \mathrm{}^2a}{m^2}}\rho +{\displaystyle \frac{\stackrel{}{v}^{\mathrm{\hspace{0.25em}2}}}{2}}+{\displaystyle \frac{V_{\mathrm{ext}}}{m}}{\displaystyle \frac{\mathrm{}^2}{2m^2}}{\displaystyle \frac{\stackrel{}{}^2\sqrt{\rho }}{\sqrt{\rho }}}\right)=0`$ (28)
and are the hydrodynamic Madelung  representation of the Gross–Pitaevskii equation.

It is easy to check that for the spatially homogeneous system, the one–dimensional counterpart of Eqs. (IV) has the following solution

$`\rho (x,t)`$ $`=`$ $`\rho _0+{\displaystyle \frac{1}{2\pi }}[\varphi ^{}(x+(p_F/m)t)\varphi ^{}(x(p_F/m)t)]`$
$`v(x,t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2m}}[\varphi ^{}(x+(p_F/m)t)+\varphi ^{}(x(p_F/m)t)]`$
assuming the slope of the phase being imprinted is much less than the unperturbed density ($`\varphi ^{}\rho _0`$). This assumption is equivalent to the condition $`d_{\mathrm{phase}}\lambda _F`$ considered by us in Sec. III and the solution (IV) coincides with expression (19). It also coincides with the Thomas–Fermi approximation, valid in one–dimensional space when the number of atoms is large enough (Sec. II). Again, the solution (IV) represents two quasisolitons (a bright and a dark one) traveling in opposite directions after imprinting (at $`t=0`$) the phase $`\varphi (x)`$. The speed of quasisolitons is equal to the speed of sound. A similar solution can be found also in two– and three–dimensional case

$`\rho (t)`$ $`=`$ $`\rho _0+\alpha [\varphi ^{}(x+ct)\varphi ^{}(xct)]`$
$`v_x(t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2m}}[\varphi ^{}(x+ct)+\varphi ^{}(xct)],`$ (30)
where $`c`$ equals the speed of sound, $`\alpha \rho _0^{1/2}`$ or $`\alpha \rho _0^{2/3}`$ respectively in two– and three–dimensional space, and the phase being imprinted changes only in ’$`x`$’ direction. Again, the solutions (30) are valid in the limit $`d_{\mathrm{phase}}\lambda _F`$.

It is clear from Fig. 3 that in the regime where $`d_{\mathrm{phase}}\lambda _F`$, the dynamics of fermionic system should be discussed rather in terms of soliton–like structures with dispersion leading to the fast broadening of pulses. Even at zero temperature the phase imprinting excites states with energy above Fermi level, depending on the width of the phase step $`d_{\mathrm{phase}}`$.The deviation of the group velocity of bright and dark wave packets from the zero–temperature sound velocity $`p_F/m`$ maybe qualitatively understood by the following argument. In the spirit of the local density approximation one considers local value of the Fermi momentum determined by the local value of the density. Since the local density is lower at the center of dark quasisoliton its velocity is thus reduced. Of course, the opposite is true for the bright quasisoliton. Using Eq. (IV) one gets local corrections to $`p_F`$ of the form $`\delta p_F=\frac{\mathrm{}}{2}\varphi ^{}(0)`$. This correction agrees with numerical results in the order of magnitude. More precise estimation of the asymmetry of velocities requires, however, taking into account dispersion effects and systematic expansion in $`\lambda _F/d_{\mathrm{phase}}`$ up to higher orders.

There is an equivalent way of calculating the kinetic–energy stress tensor $`\mathrm{T}`$ (26), involving the one–particle Wigner function

$`\mathrm{T}_{kl}`$ $`=`$ $`{\displaystyle \frac{p_kp_l}{m}w(\stackrel{}{r},\stackrel{}{p})d^3p},`$ (31)
which is, however, more suitable for the finite–temperature extension of the Thomas–Fermi approach. In equilibrium, the Wigner function for a uniform system is just the Fermi–Dirac distribution function $`f(\epsilon )=(\mathrm{exp}(\beta (\epsilon \mu ))+1)^1`$ with $`\epsilon =p^2/2m`$ and $`\mu `$ being the particle’s kinetic energy and the chemical potential (constant for a homogeneous case) respectively. The Thomas–Fermi approximation assumes, as usually, that the system is locally uniform and the appropriate equations (in one–dimensional form) are written as

$`{\displaystyle \frac{\rho }{t}}+{\displaystyle \frac{}{x}}\left(\rho v\right)=0`$
$`{\displaystyle \frac{v}{t}}+v{\displaystyle \frac{v}{x}}+{\displaystyle \frac{1}{m\rho }}{\displaystyle \frac{}{x}}{\displaystyle \frac{1}{2\pi \mathrm{}}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{p^2/mdp}{\mathrm{exp}(\beta (\frac{p^2}{2m}\mu ))+1}}=0.`$
Here, the chemical potential is position–dependent and is calculated from the relation

$`\rho ={\displaystyle \frac{1}{2\pi \mathrm{}}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{dp}{\mathrm{exp}(\beta (\frac{p^2}{2m}\mu ))+1}}.`$ (33)
We consider now the case when the temperature is finite but the gas is still strongly degenerate, i.e. $`TT_F`$. By using the Sommerfeld expansion and keeping only the first term one can easily find the expression $`\mu =(2\pi \mathrm{})^2/(8m)\rho ^2`$. It turns out that in the regime where the imprinted phase changes slowly in comparison with unperturbed density, the finite–temperature time–dependent Thomas–Fermi model defined by the Eqs. (IV) and (33) has the following solution

$`\rho =\rho _0+{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d\epsilon }{2\pi }}\left({\displaystyle \frac{f}{\epsilon }}\right)[\varphi ^{}(x+\sqrt{{\displaystyle \frac{2\epsilon }{m}}}t)\varphi ^{}(x\sqrt{{\displaystyle \frac{2\epsilon }{m}}}t)]`$
$`v={\displaystyle \frac{\mathrm{}}{2m}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑\epsilon \left({\displaystyle \frac{f}{\epsilon }}\right)[\varphi ^{}(x+\sqrt{{\displaystyle \frac{2\epsilon }{m}}}t)+\varphi ^{}(x\sqrt{{\displaystyle \frac{2\epsilon }{m}}}t)]`$
For the temperature $`T=0`$ one has $`(f/\epsilon )=\delta (\epsilon \epsilon _F)`$ and the solution (IV) reduces to the one found previously (19, IV). To verify the ansatz (IV) we put it in the set of Eqs. (IV) and keep only linear terms in $`\varphi ^{}`$ (according to $`d_{\mathrm{phase}}\lambda _F`$). The solvability condition has then the form

$`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}d\epsilon ({\displaystyle \frac{f}{\epsilon }})[\sqrt{{\displaystyle \frac{2\epsilon }{m}}}\sqrt{{\displaystyle \frac{2\epsilon _F}{m}}}]\times `$
$`\left[\varphi ^{\prime \prime }\left(x+\sqrt{{\displaystyle \frac{2\epsilon }{m}}}t\right)\pm \varphi ^{\prime \prime }\left(x\sqrt{{\displaystyle \frac{2\epsilon }{m}}}t\right)\right]=0`$ (35)
and is true for low temperatures $`TT_F`$ since then the derivative of the Fermi–Dirac distribution function strongly peaks at $`\epsilon =\epsilon _F`$, whereas the rest of the integrand is very small. The solution (IV) is identical with (20).

### IV.1 Numerical results

We have solved numerically the set of Eqs. (IV). To this end, we applied the inverse Madelung transformation inverse  to Eqs. (IV) and used the Split-Operator technique. Imaginary–time propagation method was chosen to generate the ground state of a nonuniform (trapped) Fermi gas. Real–time propagation, on the other hand, allows to investigate the phase imprinting and developing instabilities after that.

First, by imprinting the phase on a one–dimensional uniform system we again investigate the regime where $`d_{\mathrm{phase}}\lambda _F`$, however, this time by increasing the number of atoms while keeping the same the width of the phase step. Of course, larger number of atoms means smaller Fermi length. Increasing the number of atoms is easily attainable within Thomas–Fermi approach as opposed to the direct solution of the many–body Schrödinger equation. Results regarding the speeds of bright and dark quasisolitons are presented in Table 2. It is clear that reaching the above mentioned regime forces both quasisolitons to move with the same velocity equal to the speed of sound.

## V Conclusions

In this paper we have studied in detail the possibility of generation of soliton–like structures in ultracold Fermi gases using the method of phase imprinting. Generation of such structures is feasible with existing experimental setups. The soliton–like structures are generated because in the process of phase imprinting in the limit in which the characteristic length of the imprint $`d_{\mathrm{phase}}\lambda _F`$ only a narrow band of states with momenta close to $`\pm p_F`$ are excited. The life time of soliton–like structures can be quite long even at not extremely low temperatures allowing for observations. We have presented both numerical results and analytic theory explaining the mechanism of generation and evolution of soliton–like structures both in homogeneous quasi 1D cylinders and in quasi 1D harmonic traps. It will be interesting to generalize the presented theory to the case of two–component Fermi gas in normal phase.

###### Acknowledgements.

We thank C. Salomon, J. Dalibard, L. Santos, A. Sanpera. This paper has been supported by Deutsche Forschungsgemeinschaft (among others in the frame of SFB 407, SP 1116, GK 282, and 436POL), TRN ’Cold Atoms and Molecules’, ESF PESC Programm ’BEC 2000+’. K.R. acknowledges support of the Subsidy by Foundation for Polish Science. Part of the results has been obtained using computers at the Interdisciplinary Center for Mathematical and Computational Modeling at Warsaw University.


format: text

text:
# The Kozai Mechanism and the Evolution of Binary Supermassive Black Holes

## 1 INTRODUCTION

The ubiquity of supermassive black holes (SMBHs) in the nuclei of many galaxies (e.g. Magorrian et al. 1998) suggests that binary and multiple SMBH systems may also be widespread. For example, SMBH binaries should form in galaxy mergers (Begelman, Blandford, & Rees 1980), which is a common process in hierarchical models for galaxy formation and evolution (e.g. Kauffmann & Haehnelt 2000; Menou, Haiman, & Narayanan 2001). SMBH binaries have been invoked to explain periodic wiggles in extragalactic radio jets (e.g. Roos, Kaastra, & Hummel 1993), periodic flares in the BL Lac source OJ 287 (e.g. Lehto & Valtonen 1996), variations in the apparent superluminal transverse velocities and position angles of the 3C 273 radio jet (Romero et al. 2000), and the “core-type” nuclear surface brightness profiles of bright elliptical galaxies (Ebisuzaki, Makino, & Okumura 1991; Quinlan & Hernquist 1997; Faber et al. 1997; Milosavljević & Merritt 2001).

SMBH binary mergers would be powerful sources of gravitational waves (Thorne & Braginsky 1976) which should be easily detectable from planned space-based gravitational wave observatories such as the Laser Interferometer Space Antenna (Bender et al. 1998). However, it is still uncertain whether binary SMBHs can evolve to small enough semimajor axes that gravitational radiation drives them to merge. Energy can be extracted from a wide binary through interactions with surrounding matter. For example, encounters with passing stars will shrink the binary, but this process will eventually be limited by the depletion of stars on sufficiently radial orbits to encounter the ever tightening binary. Detailed work on this aspect of the problem has been done by numerous authors (e.g. Quinlan 1996, Quinlan & Hernquist 1997, Milosavljević & Merritt 2001). A key uncertainty is the amount of stochastic wandering that the binary undergoes at the bottom of the potential well of the galaxy, allowing it to interact with many more stars. Another important process that has received relatively less attention is interaction with surrounding gas, which can also extract energy from the binary and drive it to merge (e.g. Armitage & Natarajan 2002 and references therein).

Large galaxies are typically the product of multiple mergers over the history of the universe. If the characteristic merger time of binary SMBHs is not much less than the characteristic time scale between mergers, then interactions with a third SMBH or a second SMBH binary will likely take place. Strong encounters within these multiple black hole systems may lead to slingshot ejections of SMBHs from the nucleus of the galaxy (Valtonen et al. 1994, Valtonen 1996). On the other hand, the very existence of the additional galaxy merger may introduce more stellar encounters with the original central binary and drive it to merge (Roos 1988). Moreover, repeated encounters between a third black hole and the inner binary can increase the probability that the inner binary attains a large eccentricity, thereby accelerating the energy and angular momentum loss by gravitational radiation (Makino & Ebisuzaki 1994).

A modification of the last scenario is a case where the third black hole has evolved to the point that it has become bound to the SMBH binary, but has not yet come close enough for an unstable three-body interaction to take place. The system then forms a hierarchical triple and can be treated as consisting of an inner binary and an outer binary. In that case, if the mutual inclination angle between the inner and outer binaries is high enough, then the time-averaged tidal gravitational force on the inner binary can induce an oscillation in its eccentricity. This effect is known as the Kozai mechanism (Kozai 1962). Provided all other dynamical influences are negligible (see section 2), an initially very small eccentricity in the inner binary will oscillate through a maximum value given by

$$e_{1,\mathrm{max}}\left(1\frac{5}{3}\mathrm{cos}^2i_0\right)^{1/2},$$
(1)
provided $`|\mathrm{cos}i_0|<(3/5)^{1/2}`$, i.e. the initial mutual inclination angle $`i_0`$ lies between roughly $`39^{}`$ and $`141^{}`$ (Kozai 1962). Given that the gravitational wave merger time is a strong function of eccentricity (approximately proportional to $`[1e_1^2]^{7/2}`$), highly inclined orbits could in principle greatly reduce the merger time. For example, a mutual inclination angle of $`56^{}`$ results in eccentricity oscillations which, when at maximum amplitude, reduce the characteristic merger time by an order of magnitude. Although the true merger time will depend on how long the inner binary spends at high eccentricity, it appears promising that a random distribution of initial inclinations will result in dramatically accelerated mergers in many cases. This is the subject we will explore in the present paper.

While this research was being completed, we learned of recent work by Miller & Hamilton (2002), who investigate a very similar idea: the use of the Kozai mechanism to accelerate binary mergers in triple black hole systems formed in globular clusters.

We begin in section 2 by going through the characteristic time scales which will affect the dynamics of triple supermassive black hole systems. Then in section 3 we summarize the evolution equations of an isolated triple used in our numerical calculations. Section 4 presents the results of those calculations, which help delineate the regions of initial condition parameter space where the Kozai mechanism plays a substantial role in accelerating the merger of the inner binary. We discuss our results and summarize our conclusions in section 5.

## 2 CHARACTERISTIC TIME SCALES

Consider a binary system consisting of two black holes with masses $`m_0`$ and $`m_1`$, with initial semimajor axis $`a_1`$ and eccentricity $`e_1`$. The time it takes for the binary to merge due to gravitational wave emission (Peters 1964) can be written as

$`t_{\mathrm{merge},\mathrm{binary}}`$ $``$ $`2.9\times 10^{12}\mathrm{yr}\left({\displaystyle \frac{m_0}{10^6\mathrm{M}_{}}}\right)^1\left({\displaystyle \frac{m_1}{10^6\mathrm{M}_{}}}\right)^1\left({\displaystyle \frac{m_0+m_1}{2\times 10^6\mathrm{M}_{}}}\right)^1\left({\displaystyle \frac{a_1}{10^2\mathrm{pc}}}\right)^4`$ (2)
$`\times f(e_1)(1e_1^2)^{7/2},`$
where $`f(e_1)`$ is a weak function of the initial eccentricity that is of order unity ($`0.979<f(e_1)<1.81`$ for all $`e_1`$). This gravitational wave merger time is a strong function of both the semimajor axis and the eccentricity. Hence gravitational radiation only becomes important late in the evolution of the binary.

Prior to that time, the evolution is dominated by interactions between the binary and surrounding material, either stars or gas. These interactions are very complex, due to the fact that the surroundings themselves evolve due to their interaction with the binary. (See Milosavljević & Merritt 2001 for a recent discussion of stellar interactions.) As we noted above, it is not yet clear whether these interactions are sufficient to harden the binary to a point where gravitational radiation then causes it to merge.

The scenario we wish to explore in this paper is one where the binary’s semimajor axis evolution has stalled because of insufficient interactions with surrounding material, and where a subsequent galaxy merger then introduces a third black hole in the system. We envisage the third black hole eventually forming a bound, hierarchical triple with the binary. The tidal gravitational torques exerted on the inner binary by the outer black hole can alter the eccentricity of the inner binary, thereby affecting the rate of gravitational wave emission.

General relativity causes periastron precession in the inner binary, with a period

$$P_{\mathrm{GR}}2.3\times 10^6\mathrm{yr}\left(\frac{m_0+m_1}{2\times 10^6\mathrm{M}_{}}\right)^{3/2}\left(\frac{a_1}{10^2\mathrm{pc}}\right)^{5/2}(1e_1^2).$$
(3)
Provided this general relativistic precession is unimportant, and the outer black hole is in a sufficiently inclined orbit around the inner binary, then the eccentricity of the inner binary will oscillate by the Kozai mechanism. The characteristic time scale for these oscillations is given by (e.g. Holman, Touma, & Tremaine 1997)

$$P_e1.3\times 10^5\mathrm{yr}\left(\frac{m_0+m_1}{2\times 10^6\mathrm{M}_{}}\right)^{1/2}\left(\frac{a_1}{10^2\mathrm{pc}}\right)^{3/2}\left(\frac{m_0+m_1}{2m_2}\right)\left(\frac{a_2/a_1}{10}\right)^3(1e_2^2)^{3/2},$$
(4)
where $`m_2`$ is the mass of the third body and $`a_2`$ and $`e_2`$ are the semimajor axis and eccentricity, respectively, of its orbit around the inner binary.

General relativistic precession can stop these eccentricity oscillations by destroying the Kozai resonance (e.g. Holman et al. 1997). A precise criterion for this not to happen is easily derived (see equation \[A6\] in Appendix):111Apart from the dependence on $`e_1`$ and factors of order unity, this is roughly equivalent to $`P_e<P_{\mathrm{GR}}`$. Equation (5) provides a more quantitatively accurate description of the behavior seen in our numerical results.

$$\frac{a_2}{a_1}<34\left(\frac{a_1}{10^2\mathrm{pc}}\right)^{1/3}\left(\frac{m_0+m_1}{2\times 10^6\mathrm{M}_{}}\right)^{1/3}\left(\frac{2m_2}{m_0+m_1}\right)^{1/3}\left(\frac{1e_1^2}{1e_2^2}\right)^{1/2}.$$
(5)
The outer black hole must therefore come quite close to the inner binary for eccentricity oscillations to take place, although the triple is still hierarchical in the sense that the semimajor axis ratio $`a_2/a_1`$ is large.

Interactions with the surrounding stars and gas may cause the outer black hole’s semimajor axis to evolve significantly over the merger time scale of the inner binary, a point to which we shall return in section 5 below. For now, however, we will neglect this fact and consider the evolution of an isolated, hierarchical black hole triple. As we will show, once the outer black hole is sufficiently close that equation (5) is satisfied, the merger time of the inner binary drops substantially for sufficiently high mutual inclinations.

## 3 EQUATIONS OF MOTION OF AN ISOLATED TRIPLE SYSTEM

We assume that the triple is hierarchical, so that the orbit of the outer black hole is much larger in size than the orbit of the inner binary. In this case the triple can be considered to consist of two binaries in approximately Keplerian orbits: the inner binary consisting of black holes with masses $`m_0`$ and $`m_1`$; and the outer binary consisting of the center of mass of the inner binary, with mass $`m_0+m_1`$, and the third black hole with mass $`m_2`$. Let the semimajor axes of the inner and outer binaries be $`a_1`$ and $`a_2`$, respectively, and the eccentricities be $`e_1`$ and $`e_2`$. We also define $`g_1`$ and $`g_2`$ to be the corresponding arguments of periastron.

The magnitude of the angular momenta $`𝐆_1`$ and $`𝐆_2`$ of the inner and outer binaries are given by

$$G_1=m_0m_1\left[\frac{Ga_1(1e_1^2)}{m_0+m_1}\right]^{1/2}$$
(6)
and

$$G_2=(m_0+m_1)m_2\left[\frac{Ga_2(1e_2^2)}{m_0+m_1+m_2}\right]^{1/2}$$
(7)
respectively, where $`G`$ is Newton’s gravitational constant. Let the total orbital angular momentum of the triple be $`𝐇=𝐆_1+𝐆_2`$. In the absence of gravitational radiation, tidal torques, or gravitational interactions with surrounding stars, this vector would be rigorously conserved. Let $`i_1`$ ($`i_2`$) be the inclination of the inner (outer) binary, i.e. the angle between $`𝐆_1`$ ($`𝐆_2`$) and $`𝐇`$. Then

$$H=G_1\mathrm{cos}i_1+G_2\mathrm{cos}i_2$$
(8)
and

$$G_1\mathrm{sin}i_1=G_2\mathrm{sin}i_2.$$
(9)
The mutual inclination angle between the two binaries is $`i=i_1+i_2`$.

Marchal (1990), Krymolowski & Mazeh (1999) and Ford, Kozinsky, & Rasio (2000) have derived the orbit-averaged Hamiltonian of an isolated, Newtonian hierarchical triple of point masses, using secular perturbation theory to octupole order, i.e. to order $`(a_1/a_2)^3`$. We adopt the equations of motion for the orbital elements of Ford et al. here,222The equations of motion of Krymolowski & Mazeh (1999) differ from those of Ford et al. (2000) by terms of order $`(a_1/a_2)^{7/2}`$ resulting from the canonical transformation of the von Zeipel method. We have corrected a sign error in the octupolar terms in the equations of motion of Ford et al. (2000). In equations (22) and (29)-(32) of their paper, all terms involving $`C_3`$ should have the opposite sign. The equations of motion of Krymolowski & Mazeh (1999) also contain the same sign error. but modify them to incorporate two general relativistic effects on the inner binary: the precession of periastron and gravitational radiation. We do this by simply adding orbit-averaged general relativistic correction terms to the Ford et al. (2000) equations of motion for the orbital elements of the inner binary, the same correction terms which would exist if the inner binary were isolated. We must also take into account the fact that $`𝐇`$ is no longer conserved because of the radiative loss $`(dG_1/dt)_{\mathrm{rad}}`$ of the inner binary’s orbital angular momentum. We do this by noting that in the absence of interactions between the inner and outer binaries, gravitational radiation acts to change the magnitude of $`𝐆_1`$, but not its direction. In addition, the vector $`𝐆_2`$ remains unchanged. The resulting equation is

$$\frac{dH}{dt}=\frac{G_1+G_2\mathrm{cos}i}{H}\left(\frac{dG_1}{dt}\right)_{\mathrm{rad}}.$$
(10)
Note that if the black holes are spinning, then the orbital plane of the inner binary can also evolve because of general relativistic spin-orbit coupling. This is a higher order effect than general relativistic precession, and we neglect it entirely, assuming that we are essentially dealing with Schwarzschild black holes.

Thus our evolution equations for the orbital elements are

$$\frac{da_1}{dt}=\frac{64G^3m_0m_1(m_0+m_1)}{5c^5a_1^3(1e_1^2)^{7/2}}\left(1+\frac{73}{24}e_1^2+\frac{37}{96}e_1^4\right),$$
(11)
$`{\displaystyle \frac{dg_1}{dt}}`$ $`=`$ $`6C_2\left\{{\displaystyle \frac{1}{G_1}}[4\theta ^2+(5\mathrm{cos}2g_11)(1e_1^2\theta ^2)]+{\displaystyle \frac{\theta }{G_2}}[2+e_1^2(35\mathrm{cos}2g_1)]\right\}`$ (12)
$`+C_3e_2e_1\left({\displaystyle \frac{1}{G_2}}+{\displaystyle \frac{\theta }{G_1}}\right)\left\{\mathrm{sin}g_1\mathrm{sin}g_2[A+10(3\theta ^21)(1e_1^2)]5\theta B\mathrm{cos}\varphi \right\}`$
$`C_3e_2{\displaystyle \frac{1e_1^2}{e_1G_1}}\left[10\theta (1\theta ^2)(13e_1^2)\mathrm{sin}g_1\mathrm{sin}g_2+\mathrm{cos}\varphi (3A10\theta ^2+2)\right]`$
$`+{\displaystyle \frac{3}{c^2a_1(1e_1^2)}}\left[{\displaystyle \frac{G(m_0+m_1)}{a_1}}\right]^{3/2},`$
$`{\displaystyle \frac{de_1}{dt}}`$ $`=`$ $`30C_2{\displaystyle \frac{e_1(1e_1^2)}{G_1}}(1\theta ^2)\mathrm{sin}2g_1`$ (13)
$`C_3e_2{\displaystyle \frac{1e_1^2}{G_1}}[35\mathrm{cos}\varphi (1\theta ^2)e_1^2\mathrm{sin}2g_110\theta (1e_1^2)(1\theta ^2)\mathrm{cos}g_1\mathrm{sin}g_2`$
$`A(\mathrm{sin}g_1\mathrm{cos}g_2\theta \mathrm{cos}g_1\mathrm{sin}g_2)]`$
$`{\displaystyle \frac{304G^3m_0m_1(m_0+m_1)e_1}{15c^5a_1^4(1e_1^2)^{5/2}}}\left(1+{\displaystyle \frac{121}{304}}e_1^2\right),`$
$$\frac{da_2}{dt}=0,$$
(14)
$`{\displaystyle \frac{dg_2}{dt}}`$ $`=`$ $`3C_2\left\{{\displaystyle \frac{2\theta }{G_1}}[2+e_1^2(35\mathrm{cos}2g_1)]+{\displaystyle \frac{1}{G_2}}[4+6e_1^2+(5\theta ^23)(2+3e_1^25e_1^2\mathrm{cos}2g_1)]\right\}`$ (15)
$`C_3e_1\mathrm{sin}g_1\mathrm{sin}g_2\{{\displaystyle \frac{4e_2^2+1}{e_2G_2}}10\theta (1\theta ^2)(1e_1^2)`$
$`e_2({\displaystyle \frac{1}{G_1}}+{\displaystyle \frac{\theta }{G_2}})[A+10(3\theta ^21)(1e_1^2)]\}`$
$`C_3e_1\mathrm{cos}\varphi \left[5B\theta e_2\left({\displaystyle \frac{1}{G_1}}+{\displaystyle \frac{\theta }{G_2}}\right)+{\displaystyle \frac{4e_2^2+1}{e_2G_2}}A\right],`$
$$\frac{de_2}{dt}=C_3e_1\frac{1e_2^2}{G_2}\left[10\theta (1\theta ^2)(1e_1^2)\mathrm{sin}g_1\mathrm{cos}g_2+A(\mathrm{cos}g_1\mathrm{sin}g_2\theta \mathrm{sin}g_1\mathrm{cos}g_2)\right],$$
(16)
and

$$\frac{dH}{dt}=\frac{32G^3m_0^2m_1^2}{5c^5a_1^3(1e_1^2)^2}\left[\frac{G(m_0+m_1)}{a_1}\right]^{1/2}\left(1+\frac{7}{8}e_1^2\right)\frac{G_1+G_2\theta }{H}.$$
(17)
The quantities $`C_2`$ and $`C_3`$ multiply the quadrupole and octupole perturbation terms, respectively. They are defined by (Ford et al. 2000)

$$C_2=\frac{Gm_0m_1m_2}{16(m_0+m_1)a_2(1e_2^2)^{3/2}}\left(\frac{a_1}{a_2}\right)^2,$$
(18)
and

$$C_3=\frac{15Gm_0m_1m_2(m_0m_1)}{64(m_0+m_1)^2a_2(1e_2^2)^{5/2}}\left(\frac{a_1}{a_2}\right)^3.$$
(19)
Note that the octupole terms vanish if $`m_0=m_1`$. The quantities $`B`$ and $`A`$ in these terms are given by

$$B=2+5e_1^27e_1^2\mathrm{cos}2g_1$$
(20)
and

$$A=4+3e_1^2\frac{5}{2}(1\theta ^2)B.$$
(21)
Finally, the quantity $`\theta `$ is the cosine of the mutual inclination of the binaries,

$$\theta =\mathrm{cos}i=\frac{H^2G_1^2G_2^2}{2G_1G_2},$$
(22)
and $`\varphi `$ is the angle between the periastron directions,

$$\mathrm{cos}\varphi =\mathrm{cos}g_1\mathrm{cos}g_2\theta \mathrm{sin}g_1\mathrm{sin}g_2.$$
(23)
Note that equation (22) determines the evolution of the mutual inclination through the time dependence of the eccentricities, $`a_1`$, and $`H`$. The terms in the equations of motion which remain after setting $`C_2=C_3=0`$ are the general relativistic correction terms. The last term in equation (12) is the general relativistic precession term, and the terms involving the speed of light $`c`$ in equations (11), (13), and (17) are the gravitational radiation terms.

We may immediately deduce some important features of our evolution equations. Equations (13) and (16) imply that if both the inner and outer binaries are circular ($`e_1=e_2=0`$), then they will remain that way. An eccentric outer binary will produce a nonzero eccentricity in an initially circular inner binary, and vice-versa, but only because of the octupole interaction terms.

Note that if we switch off the quadrupole and octupole interaction terms by setting $`C_2=C_3=0`$, then an immediate consequence is $`d\theta /dt=0`$, as may be verified by differentiating equation (22) directly. In other words, general relativistic effects alone do not change the mutual inclination of the inner and outer binaries if there are no interactions between them. This is as one would expect, given our treatment of the gravitational wave angular momentum loss in equation (10).

Our evolution equations also exhibit an important scaling. If all masses and all initial semimajor axes are multiplied by some constant factor, then the merger time and all other time scales change by the same multiplicative factor. This is a direct consequence of the nature of gravity, where mass, length, and time have the same dimensions in units where $`G=c=1`$, and it reduces the parameter space of masses and semimajor axes which need to be explored numerically. Another similarly useful fact is that the equations are unchanged when $`m_0`$ and $`m_1`$ in the inner binary are interchanged and $`g_1`$ or $`g_2`$ is increased by $`180^{}`$. The former corresponds to a relabeling of the inner binary black holes, while the latter represents a spatial inversion of the two orbits.

While our implementation of the general relativistic effects has been heuristic, the general relativistic precession terms for the inner binary can in fact be rigorously justified by orbit-averaging the post-Newtonian Hamiltonian. To octupole order, the resulting Hamiltonian is

$`\overline{}(g_1,g_2,G_1,G_2)`$ $`=`$ $`C_2\left[(2+3e_1^2)(13\theta ^2)15e_1^2(1\theta ^2)\mathrm{cos}2g_1\right]`$ (24)
$`+C_3e_1e_2\left[A\mathrm{cos}\varphi +10\theta (1\theta ^2)(1e_1^2)\mathrm{sin}g_1\mathrm{sin}g_2\right]`$
$`+{\displaystyle \frac{G^2m_0m_1}{c^2a_1^2}}\left[{\displaystyle \frac{15m_1^2+15m_0^2+29m_0m_1}{8(m_0+m_1)}}{\displaystyle \frac{3(m_0+m_1)}{(1e_1^2)^{1/2}}}\right].`$
A rigorous derivation of the radiation terms would be considerably more difficult, and is beyond the scope of this paper.333One of the issues is that the Newtonian portion of the equations of motion were derived by using the invariable plane as the reference plane. Because this plane slowly evolves under the influence of gravitational radiation, it may be that a rigorous set of equations of motion should include terms to reflect this slow evolution. We believe that our equations of motion will capture the overall time evolution of the inner binary, because the early phases of this evolution will be dominated by interactions with the outer black hole, while the late phases will be dominated by gravitational radiation. In both those limits our evolution equations for the orbital elements of the inner binary are rigorously correct. It is possible, however, that interesting effects may occur during the transition between these two phases of evolution, which may not be captured by our equations.

In the absence of gravitational radiation, our equations conserve the orbit averaged Hamiltonian in equation (24). We use this fact to check the accuracy of our numerical integrations.

The reader should bear in mind that our equations are based on orbit averaging and an expansion in the semimajor axis ratio $`a_1/a_2`$, and may therefore be inaccurate when this ratio is not very small. We will show in the next section that the octupole terms usually have fairly small effects on the merger time of the inner binary. Hence our neglect of even higher order terms, and our use of such an expansion in the first place, is probably justified for the primary purpose of this paper.

The major exception to this is the issue of stability. We are investigating whether the presence of a third black hole can cause the inner binary to merge before the third black hole comes close enough to cause an unstable interaction. There is a limit to how small $`a_2/a_1`$ can be for the triple to be stable and for the evolution we calculate to be valid. Completely general stability criteria for hierarchical triples in arbitrary initial configurations do not exist. We adopt as a stability criterion the empirical formula used by Mardling & Aarseth (2001), which may be written as

$$\frac{a_2}{a_1}>\frac{2.8}{1e_2}\left[\left(1+\frac{m_2}{m_0+m_1}\right)\frac{1+e_2}{(1e_2)^{1/2}}\right]^{2/5}.$$
(25)
This criterion was derived for Newtonian coplanar prograde orbits of the inner and outer binaries. Inclined orbits, which are the major focus of this paper, are expected to be more stable, so inequality (25) provides a conservative stability limit. This inequality is also conservative in the sense that it is either close to or above the more complicated stability criterion proposed by Eggleton & Kiseleva (1995), which was empirically verified over a wide range of mass ratios including nearly all those we numerically investigate in the next section. All the calculations we present in this paper are done for systems which are stable according to inequality (25).

The stability criterion (25) may be combined with the constraint of inequality (5) that general relativistic precession not destroy the Kozai resonance to give a lower limit on the inner semimajor axis for which the Kozai mechanism can reduce the inner binary merger time. 444We are grateful to the referee for suggesting this connection between stability and the precession constraint, and for emphasizing to us the importance of discussing the stability of triple systems in general. The result is

$$a_16\times 10^6\mathrm{pc}\left(\frac{m_0+m_1+m_2}{m_0+m_1}\right)^{6/5}\left(\frac{m_0+m_1}{2m_2}\right)\left(\frac{m_0+m_1}{2\times 10^6\mathrm{M}_{}}\right)\frac{(1+e_2)^{27/10}}{(1e_2)^{21/10}(1e_1^2)^{3/2}}.$$
(26)
A stable triple with initial conditions such that the inner semimajor axis $`a_1`$ violates this inequality will not have an accelerated inner binary merger by the Kozai mechanism. However, in that case $`a_1`$ is so small that the binary merger time by equation (2) is

$`t_{\mathrm{merge},\mathrm{binary}}`$ $``$ $`0.3\mathrm{yr}\left({\displaystyle \frac{m_0}{10^6\mathrm{M}_{}}}\right)^1\left({\displaystyle \frac{m_1}{10^6\mathrm{M}_{}}}\right)^1\left({\displaystyle \frac{m_0+m_1}{2\times 10^6\mathrm{M}_{}}}\right)^3\left({\displaystyle \frac{m_0+m_1+m_2}{m_0+m_1}}\right)^{24/5}`$ (27)
$`\left({\displaystyle \frac{m_0+m_1}{2m_2}}\right)^4{\displaystyle \frac{(1+e_2)^{54/5}f(e_1)}{(1e_2)^{42/5}(1e_1^2)^{5/2}}}.`$
At least for roughly equal mass triples and outer eccentricities that are not too high, this implies that the inner binary would already be rapidly merging and needs no help from the Kozai mechanism.

On the other hand, equations (26) and (27) have a rather strong dependence on $`e_2`$, so if the outer black hole’s orbit is very eccentric, then the resulting triple may be too unstable for there to be time for the Kozai mechanism to operate. However, numerical simulations of black hole binary evolution by Milosavljević & Merritt (2001) suggest that this is rather unlikely. They find rather modest initial eccentricities when the binary first forms, and the subsequent evolution of the binary does not go to very high eccentricity. Hence there is probably plenty of parameter space for stable hierarchical triples to form and evolve under the Kozai mechanism.

## 4 NUMERICAL RESULTS

### 4.1 Detailed Evolution of a Triple Consisting of Nearly Equal Mass Black Holes

We have numerically integrated equations (11)-(17) for a wide range of possible black hole masses and initial conditions. Even for fixed black hole masses, the parameter space of possible initial conditions is huge, and it helps to understand which are the most important quantities affecting the merger time. We therefore focus first on the detailed evolution of a particular triple consisting of nearly equal mass black holes. Specifically, we consider an inner binary with a $`2\times 10^6`$ M black hole and a $`10^6`$ M black hole, about which orbits another $`10^6`$ M black hole ten times further out. We have purposely chosen the two inner black holes to have different masses so as to allow for octupole interaction effects in the evolution.

Figure 1 depicts the evolution of inner semimajor axis $`a_1`$, inner eccentricity $`e_1`$, and outer eccentricity $`e_2`$ for a particular choice of initial conditions. The inner binary starts out nearly circular ($`e_1=0.1`$) with a semimajor axis which would give a gravitational wave merger time $`t_{\mathrm{merge},\mathrm{binary}}=9.3\times 10^9`$ yr if the binary were isolated. However, the presence of the outer black hole induces large amplitude eccentricity oscillations in the inner binary through the Kozai mechanism. The time spent at higher eccentricity reduces the gravitational wave merger time of the two inner black holes by roughly a factor thirty for the particular case shown.

In contrast to $`e_1`$, the outer eccentricity $`e_2`$ stays close to its original value throughout the evolution. Equation (16) shows that $`e_2`$ only changes as a result of octupole interactions, which are weaker than the quadrupolar interactions driving the inner binary eccentricity oscillations. The lack of strong evolution in the outer eccentricity is a generic feature of all the numerical calculations we have done, implying that $`e_2`$ is unlikely to be driven to high enough values to destabilize an isolated triple.

As might be expected, the merger time turns out to be rather insensitive to the initial orientation of the roughly circular inner orbit, i.e. the initial value of the inner argument of periastron $`g_1`$. Figure 2 shows the evolution of the same triple as depicted in figure 1, except that the initial value of $`g_1`$ has been shifted by $`90^{}`$. This results in a merger time which is about seven percent longer than that shown in figure 1.

The overall behavior of the eccentricity oscillations depicted in figures 1(b) and 2(b) can be understood by examining the evolution in the $`e_1`$ vs. $`g_1`$ phase space. The initial phase space trajectories for triples with the same total angular momentum as the triples of figures 1 and 2 are depicted in figure 3. Depending on the initial conditions, the evolution in $`g_1`$ is either one of libration about $`g_1=90^{}`$ or $`270^{}`$, or one of circulation. The quadrupolar fixed point of the Kozai resonance lies inside the smallest libration contour at $`\mathrm{cos}g_1=0`$. The only difference between the triples of figures 1 and 2 is that the former starts off circulating while the latter starts off librating. Figure 4 depicts snapshots of the phase space trajectories during the course of the evolution shown in figures 1 and 2. As shown in figure 4(b), gravitational radiation drives an initially librating inner binary over into a circulating inner binary, thereby causing the minimum eccentricity to drop until it crosses the separatrix. This separatrix crossing corresponds to the momentary zero eccentricity spike just before $`10^7`$ years in figure 2(b). Thereafter, the evolution is very similar to the case where the inner binary starts off circulating: the minimum eccentricity rises until gravitational radiation becomes so strong that it starts to circularize the orbit.

Figure 5 shows the overall dependence of the merger time for varying initial inner argument of periastron $`g_1`$, again for an inner orbit which is initially roughly circular: $`e_1=0.1`$. While the overall effect of $`g_1`$ on the merger time is not too great, the exact value of the merger time is extremely sensitive to the initial value of $`g_1`$ for triples close to the libration/circulation separatrix, or for those triples which start out librating ($`g_1`$ initially around $`90^{}`$ or $`270^{}`$) and therefore subsequently cross the separatrix. As shown in figure 3, the region near the separatrix is chaotic, and this chaos arises from the octupole interaction terms. If these terms are neglected, the system then lacks the necessary degrees of freedom to exhibit chaos, and the dependence of the merger time on the initial value of $`g_1`$ is much smoother (the dashed line in figure 5). Note that the separatrix is associated with passing through very small values of $`e_1`$ \[cf. figure 2(b)\]. One might therefore worry that our results suffer from numerical problems associated with the $`1/e_1`$ singularity in the octupole term of equation (12). However, following the suggestion of Ford et al. (2000), we have removed this singularity in our numerical integrations by first transforming the dependent variables from $`(g_1,e_1,g_2,e_2)`$ to $`(e_1\mathrm{cos}g_1,e_1\mathrm{sin}g_1,e_2\mathrm{cos}g_2,e_2\mathrm{sin}g_2)`$. In addition, the evolution of the Hamiltonian, equation (24), is completely smooth through the separatrix crossing, changing only as a result of gravitational radiation losses. We therefore believe the “noise” exhibited in figure 5 is physical, and is caused by chaotic behavior during evolution through the separatrix.

The initial orientation of the inner orbit affects the merger time much more substantially when the inner orbit is more eccentric. Figure 6 depicts the merger time as a function of the initial inner eccentricity $`e_1`$ for two values of the initial inner argument of periastron: $`g_1=0`$ and $`90^{}`$. The former case corresponds to a circulating inner binary which starts out with an eccentricity which is at the minimum in the Kozai oscillation (see figure 3). Hence in this case the Kozai mechanism always speeds up the merger compared to the time $`t_{\mathrm{merge},\mathrm{binary}}`$ it would take an isolated binary with the same initial eccentricity to merge. In contrast, inner binaries with initial $`g_1=90^{}`$ start at the minimum in the Kozai eccentricity oscillation for values of $`e_1`$ which are below the Kozai fixed point, and at the maximum for values of $`e_1`$ above it. We would therefore expect the Kozai mechanism to reduce the merger time in the former case and increase it in the latter. This expectation is confirmed by the behavior shown in figure 6, provided the binary does not start out too near the Kozai fixed point which occurs at $`e_10.8`$ for this set of triples with fixed mutual inclination angles. Near the Kozai fixed point, the merger time is generally increased by the Kozai mechanism, due perhaps to the fact that gravitational radiation drives the eccentricity at the Kozai fixed point down to lower values with time. We therefore conclude that the Kozai mechanism always acts to reduce the merger time of an initially nearly circular binary, but for eccentric inner binaries it will either reduce or increase the merger time, depending on the orientation of the inner binary within the triple system.

The $`g_1=90^{}`$ curve in figure 6 exhibits the chaos we typically find for inner binaries which start out librating, or circulating near the separatrix. This chaos is also present at the lowest values of eccentricity $`e_1`$ in the $`g_1=0`$ curve, as the inner binary then starts out in the chaotic zone shown in figure 3.

We turn now to the effects of the orientation and eccentricity of the outer binary on the merger time of the inner binary. Equations (12)-(13) and (18)-(19) show that increasing the eccentricity $`e_2`$ of the outer binary at fixed semimajor axis $`a_2`$ strengthens both the quadrupolar and octupolar interaction terms, with the latter being enhanced over the former. This is of course physically reasonable as the distance of closest approach of the outer black hole with the inner binary is smaller. On the other hand, the outer argument of periastron $`g_2`$ only affects the evolution of the inner binary through the octupolar interaction terms. We would therefore expect that increasing $`e_2`$ would generally decrease the merger time through the (quadrupolar) Kozai mechanism, but that the initial value of $`g_2`$ could modify this significantly at high eccentricity due to the enhanced octupolar effects. These expectations are confirmed by our numerical calculations shown in figure 7. Note, however, that the effects of $`e_2`$ and $`g_2`$ are not that large, at least for this particular triple. We have therefore chosen to fix their initial values to be $`e_2=0.1`$ and $`g_2=90^{}`$ for all our other numerical calculations.

Perhaps the most important initial condition parameters affecting the inner binary merger time are the semimajor axes $`a_1`$ and $`a_2`$ and the mutual inclination of the inner and outer orbits. Figure 8 shows the dependence of the merger time on the semimajor axis ratio $`a_2/a_1`$ for triples with fixed initial inclination $`i=80^{}`$ (the same as in all previous figures) and various values of the initial inner semimajor axis $`a_1`$. For large initial semimajor axis ratios $`a_2/a_1`$ (to the right of the dashed line in figure 8), general relativistic precession destroys the Kozai resonance and the eccentricity of the inner binary is largely unaffected by the outer black hole. The merger time in this case is then the same as that of an isolated binary, and is given by equation (2). On the other hand, if the outer black hole comes sufficiently close (to the left of the dashed line in figure 8), then the Kozai resonance exists and substantial reduction in the merger time occurs as a result of eccentricity oscillations. The equation for the dashed line itself which separates these two regimes comes from using the binary merger time in equation (2) as a proxy for the inner semimajor axis $`a_1`$, and employing equation (5):

$`t_{\mathrm{merge}}`$ $`=`$ $`1.2\times 10^6\mathrm{yr}\left({\displaystyle \frac{a_2/a_1}{10}}\right)^{12}\left({\displaystyle \frac{2m_2}{m_0+m_1}}\right)^4\left({\displaystyle \frac{m_0}{10^6\mathrm{M}_{}}}\right)^1\left({\displaystyle \frac{m_1}{10^6\mathrm{M}_{}}}\right)^1`$ (28)
$`\left({\displaystyle \frac{m_0+m_1}{2\times 10^6\mathrm{M}_{}}}\right)^3{\displaystyle \frac{(1e_2^2)^6}{(1e_1^2)^{5/2}}}f(e_1).`$
Figure 9 shows how the variation of merger time with initial semimajor axis ratio depends on the initial mutual inclination angle of the triple for initial $`a_1=3.16\times 10^3`$ pc. As expected, significant reduction in the merger time can only occur when the initial inclination is high enough that the Kozai resonance is present: $`|\mathrm{cos}i|(3/5)^{1/2}0.77`$. If this criterion is satisfied, then values of $`a_2/a_1`$ satisfying equation (5), i.e. $`a_2/a_1<17.7`$ for this particular triple, will exhibit accelerated mergers, with higher inclinations showing the fastest merger times. As Miller & Hamilton (2002) point out for the stellar mass case, retrograde triples ($`\mathrm{cos}i<0`$) generally produce faster merger times, at least for nearly equal masses.

### 4.2 Merger Times for Other Combinations of Masses

As we noted above in section 3, our numerical results for the nearly equal mass triples of the previous subsection can be scaled to all triples with the same mass ratio, $`m_0:m_1:m_2=2:1:1`$, provided the initial semimajor axes are scaled by the same factor as the masses. The merger time then scales by exactly the same factor. We have also investigated triples consisting of substantially unequal masses, and our results may also be scaled to other, similar mass ratio triples in the same manner.

Figure 10 depicts the merger time for the same inner binary as shown in figure 8, but with different outer black hole masses. All cases shown start with an initial inner semimajor axis of $`a_1=3.16\times 10^3`$ pc, for which the merger time would be $`9.3\times 10^9`$ yr if the binary were isolated. Because we do not account for gravitational wave losses associated with the outer binary, all the calculations shown in the figure have $`a_2/a_1`$ large enough so that the nominal gravitational wave merger time of the outer binary from equation (2) is at least ten times longer than that of the inner binary.

The $`10^5`$ M curve shows that an outer black hole with a much smaller mass than that of the inner binary does not have a substantial effect on the merger time of that binary. On the other hand, a larger outer black hole mass exerts a stronger tidal perturbation on the inner binary. This reduces the merger time significantly compared to the nearly equal mass case, and also does it at larger semimajor axis ratios. In agreement with equation (5), the value of $`a_2/a_1`$ required for the onset of the Kozai resonance scales with the outer black hole mass as $`m_2^{1/3}`$. Triples in which the outer black hole mass is much larger than that of the inner binary might arise from the merger of a smaller galaxy containing a stalled binary with a larger galaxy containing a bigger black hole.

Figure 11 presents a more detailed look at the dependence of the merger time on the semimajor axes $`a_1`$ and $`a_2`$ for the case of an outer black hole mass of $`m_2=10^8`$ M. With the exception of the larger outer black hole mass, the initial triple parameters are identical to those of figure 8. Note the change in scales on the axes: the larger outer black hole mass greatly reduces the merger time at much larger values of $`a_2/a_1`$. The dashed line comes once again from equation (28), and separates the region on the left where the Kozai resonance exists from that on the right where relativistic periastron precession destroys the resonance.

Figure 12 shows results for a substantially unequal mass inner binary consisting of $`10^8`$ and $`10^6`$ M black holes, in a triple with an outer $`10^8`$ M black hole. We choose initial inner semimajor axes $`a_1`$ to be ten times larger than those of figures 8 and 11 so that $`t_{\mathrm{merge},\mathrm{binary}}`$ is roughly the same. Very substantial reductions in the merger time occur once again, but for semimajor axis ratios $`a_2/a_1`$ that are even smaller than in the roughly equal mass case. Larger outer black hole masses can expand this range of $`a_2/a_1`$ (see equation ), but would also accelerate the gravitational wave merger of the outer binary, increasing the likelihood of an unstable encounter.

We have only shown results here for prograde orbits in triples with substantially unequal masses. In contrast to the nearly equal mass case, retrograde orbits do not produce merger times which are very different from their prograde counterparts. This asymmetry between prograde and retrograde orbits arises from terms in the equations of motion which depend on odd powers of $`\theta `$, which are neglected in the usual Kozai approximation (see Appendix A). These terms are also very small for the unequal mass cases we have calculated here.

## 5 DISCUSSION AND CONCLUSIONS

We have shown that in cases where a bound hierarchical black hole triple forms, substantial reductions in the gravitational wave merger time of the inner binary can take place for reasonably large ranges of initial conditions of the triple. For example, figure 9 shows that for a triple consisting of nearly equal mass black holes, the merger time can be reduced by more than a factor of ten in over fifty percent of all cases, provided the mutual inclination angle is randomly distributed in solid angle. However, we have not addressed the issue of what the distribution of initial conditions is likely to be, as that would entail following the dynamical formation of the triple itself by interactions with the surrounding stars and gas. As we have shown in the previous section, all the initial parameters of the triple affect whether and by how much the merger of the inner binary is altered by the presence of the outer black hole. However, the most important are the mutual inclination angle $`i`$ and the semimajor axes $`a_1`$ and $`a_2`$. There is very little effect on the merger time of the inner binary unless the semimajor axis ratio $`a_2/a_1`$ is small enough that general relativistic precession does not destroy the Kozai resonance (equation ).

An issue related to the uncertainty in initial conditions is the fact that we have not considered the evolution of the semimajor axis of the outer black hole due either to interactions with surrounding material or to gravitational radiation. The latter is not significant for the calculations shown here, as we have always chosen $`a_2/a_1`$ and $`e_2`$ to be such that the gravitational wave evolution time scale of the outer black hole is much longer than the merger time scale of the inner binary. Note that $`e_2`$ evolves purely due to octupolar interactions (equation ), and generally does not undergo the dramatic changes that $`e_1`$ exhibits. Hence we do not expect dramatic reductions in the outer binary merger time scale due to interactions with the inner binary.

However, that still leaves open the issue of whether the inner binary will merge before interactions with surrounding material cause the outer black hole to come sufficiently close that an unstable three body encounter occurs. It would be interesting to follow the evolution of the triple by including the time dependence of $`a_2`$ starting at large enough values that the Kozai resonance is destroyed by general relativistic periastron precession.

In addition to the Kozai mechanism, the presence of a third black hole may affect the merger time scale of the inner binary in other ways. In particular, being bound in a triple forces the inner binary to “wander” through space, albeit in a regular orbit rather than a stochastic trajectory. This may enhance the binary’s interactions with surrounding stars beyond what it would have had if the outer black hole were not present, and therefore the hardening rate might be increased.

Eccentricity oscillations in the inner binary may be induced by other sources of tidal gravitational fields besides that of an orbiting outer black hole, e.g. nearby matter inhomogeneities or an aspherical distribution of surrounding stars and gas. In order to beat general relativistic periastron precession, equation (A6) implies that the quadrupole moment $`\rho _{2m}`$ of these exterior matter distributions must satisfy

$$\rho _{2m}\frac{m_2}{4\pi a_2^3}>2\times 10^6\mathrm{M}_{}\mathrm{pc}^3\left(\frac{m_0+m_1}{2\times 10^6\mathrm{M}_{}}\right)^2\left(\frac{a_1}{10^2\mathrm{pc}}\right)^4(1e_1^2)^{3/2}.$$
(29)
Note that this critical quadrupole moment has a strong scaling with the binary semimajor axis $`a_1`$, and it may be possible to achieve given the actual stellar mass densities observed in galactic nuclei (e.g. Faber et al. 1997).

We thank E. Ford, S. Hughes, M. C. Miller, N. Murray, S. J. Peale, S. Tremaine, and Y. Wu for useful discussions, and the referee for suggestions that greatly improved this paper. This work was supported in part by NASA grants NAG5 7075 and NAG5 3646, and NSF grant AST 9970827. We also acknowledge the hospitality of the Canadian Institute for Theoretical Astrophysics, where part of this work was completed.

## Appendix A APPENDIX: ANALYTIC ESTIMATES NEGLECTING OCTUPOLE TERMS

In this appendix we present analytic calculations which are useful in understanding how general relativistic precession modifies the Kozai mechanism. We neglect the octupole terms in the hierarchical triple evolution equations (11)-(17) throughout this section. The only parameter of the outer binary which then changes is $`g_2`$, and this evolution does not affect the inner binary.

Neglecting gravitational radiation, the evolution of the inner binary is determined solely by equations (12) and (13):

$`{\displaystyle \frac{dg_1}{dt}}`$ $`=`$ $`6C_2\left\{{\displaystyle \frac{1}{G_1}}[4\theta ^2+(5\mathrm{cos}2g_11)(1e_1^2\theta ^2)]+{\displaystyle \frac{\theta }{G_2}}[2+e_1^2(35\mathrm{cos}2g_1)]\right\}`$ (A1)
$`+{\displaystyle \frac{3}{c^2a_1(1e_1^2)}}\left[{\displaystyle \frac{G(m_0+m_1)}{a_1}}\right]^{3/2},`$
and

$$\frac{de_1}{dt}=30C_2\frac{e_1(1e_1^2)}{G_1}(1\theta ^2)\mathrm{sin}2g_1.$$
(A2)
Kozai (1962) has provided an approximate analytic solution to these quadrupole evolution equations when all general relativistic effects are neglected. The approximation neglects the term multiplied by $`\theta /G_2`$ in equation (A1) which is smaller than the first term by $`G_1/G_2(a_1/a_2)^{1/2}`$ for comparable masses. Also, equation (22) implies that

$$2G_1G_2\theta +G_1^2=H^2G_2^2,$$
(A3)
which is a constant when octupole terms and gravitational radiation are neglected. Again taking $`G_2>>G_1`$, this gives an approximate integral of motion

$$\mathrm{\Theta }(1e_1^2)\theta ^2\frac{m_0+m_1}{Ga_1m_0^2m_1^2}\left(\frac{H^2G_2^2}{2G_2}\right)^2.$$
(A4)
For $`\mathrm{\Theta }<3/5`$, Kozai’s solutions exhibit two classes of dynamical behavior in the evolution of the argument of perihelion $`g_1`$: libration about $`g_1=90^{}`$ or $`270^{}`$, and circulation. There then exists a critical inclination angle for which the librating solutions degenerate to a fixed point in the $`e_1`$ vs. $`g_1`$ phase space: $`\mathrm{cos}^2i_{\mathrm{crit}}=3(1e_1^2)/5`$. For $`\mathrm{\Theta }>3/5`$, the resonance vanishes and only circulating solutions exist. In this case the amplitude of eccentricity oscillations is quite small, so the existence of the resonance can be used as a necessary criterion for determining when large amplitude eccentricity oscillations become possible. Because $`\mathrm{\Theta }<3/5`$ implies that $`\mathrm{cos}^2i<3/5/(1e_1^2)`$, large amplitude oscillations are only possible for high inclinations: $`39^{}\mathrm{cos}^1(3/5)^{1/2}<i<180^{}\mathrm{cos}^1(3/5)^{1/2}141^{}`$.

General relativistic precession is itself circulation in $`g_1`$, so general relativity would be expected to reduce the parameter space where librations exist. Equation (A2) implies that the fixed point, if it exists, still occurs at $`g_1=90^{}`$ or $`270^{}`$. If we continue to neglect the term multiplied by $`\theta /G_2`$ in equation (A1), then setting $`dg_1/dt=0`$ gives

$$\mathrm{cos}^2i_{\mathrm{crit}}=\frac{3}{5}(1e_1^2)\frac{4G(m_0+m_1)^2a_2^3(1e_2^2)^{3/2}}{5c^2m_2a_1^4(1e_1^2)^{1/2}}.$$
(A5)
As expected, general relativistic precession decreases $`\mathrm{cos}^2i_{\mathrm{crit}}`$, so that the fixed point can only exist if the inclination angle is pushed higher. This can only be achieved so long as the right hand side of equation (A5) remains positive, implying that

$$\frac{a_2^3}{a_1^3}<\frac{3c^2m_2a_1(1e_1^2)^{3/2}}{4G(m_0+m_1)^2(1e_2^2)^{3/2}},$$
(A6)
which gives equation (5).

Equation (A5) can be rearranged to give an equation for the eccentricity $`e_{1,0}`$ at the fixed point for triples with specified total angular momentum $`H`$, and therefore specified Kozai constant $`\mathrm{\Theta }`$ within the Kozai approximation,

$$1e_{1,0}^2=\left[\frac{5\mathrm{\Theta }}{3}+\frac{4G(m_0+m_1)^2a_2^3(1e_2^2)^{3/2}(1e_{1,0}^2)^{1/2}}{3c^2m_2a_1^4}\right]^{1/2}.$$
(A7)
Solving this equation for $`e_{1,0}`$ requires solving a quartic, but written in this way it is again obvious that general relativistic precession lowers the value of $`e_{1,0}`$ for triples of given angular momentum, thereby increasing the circulating portion of phase space at the expense of libration.