text:
# 1 Introduction

## 1 Introduction

The phase structure of QCD at high temperature and non-vanishing baryon chemical potential has been subject of intense research in recent years. Heavy-ion collisions at high energies have been and are being explored at SPS/CERN and RHIC/BNL  in search for signals of the Quark-Gluon Plasma (QGP). Large-scale lattice QCD computations at finite temperature have been performed , and first extensions to non-zero baryon chemical potential appear now to be feasible. It has proven possible to trace out the phase boundary line $`T_c(\mu )`$ separating the hadronic phase from the QGP phase for $`N_f=4`$ , $`N_f=2`$  and $`N_f=3`$  flavors of quarks up to quark chemical potentials $`\mu `$ of order $`T_c`$. First numerical results for the QCD Equation of State (EoS), i.e. the pressure $`p(T,\mu )`$ and the quark number density $`n_q(T,\mu )`$ are also available for $`N_f=2+1`$  and $`N_f=2`$ . As well as being of intrinsic theoretical interest, such studies provide conceptual guidance for current heavy ion collision experiments at SPS and RHIC, where the chemical freeze-out occurs at $`\mu _{f.o.}100`$ MeV, (baryon chemical potential $`\mu _B300`$ MeV)  and $`\mu _{f.o.}15`$ MeV, ($`\mu _B45`$ MeV) , respectively.
Systematic perturbative expansions of the QCD equation of state within the framework of thermal field theory show bad convergence even for very large temperatures (several times $`T_c`$) far beyond the region accessible to present experiments . Various techniques, such as dimensional reduction, screened perturbation theory or hard-thermal loop (HTL) perturbation theory show better convergence and good agreement with lattice results for $`T\text{ }>3T_c`$ . Various interpretations of the lattice data have been attempted in terms of physical quantities, most prominently as the EoS of a gas of massive quark and gluon quasiparticles. Their thermally generated masses are based on perturbative calculations carried out in the HTL scheme . This approach has been extended to non-vanishing quark chemical potential and good agreement with finite $`\mu `$ lattice calculations for $`N_f=2+1`$ flavors has been found . More recently, the QGP has also been described in terms of a condensate of $`Z_3`$ Wilson lines  and by more refined quasiparticle models based on the HTL-resummed entropy and (NLO) extensions thereof . These models have found support from resummed perturbation theory  for temperatures $`T\text{ }>3T_c`$. However, they have difficulties explaining the dropping of the thermal gluon screening mass in the vicinity of the phase transition. An improved quasiparticle model  shows the correct temperature dependence of the Debye mass and reproduces lattice thermodynamical quantities such as the pressure, the energy density and the entropy density very well. The main new ingredient of this model is a phenomenological parametrization of (de)confinement.
In the present work, this improved quasiparticle model is extended to finite quark chemical potential $`\mu `$. In section 2, a brief review of the quasiparticle model with confinement is given. The extension of the model to finite quark chemical potential $`\mu `$ is discussed in detail in section 3. Numerical results are presented in chapter 4. The phase boundary line $`T_c(\mu )`$ that separates the hadronic from the QGP-phase is discussed and the quasiparticle model result is compared to recent lattice simulations. Results for the pressure difference from $`\mu =0`$ and the quark number density for various values of the chemical potential $`\mu `$ are also presented and compared to recent lattice simulations. A summary is given in section 5.

## 2 Quasiparticle model with confinement

It is possible to describe the EoS of hot QCD at vanishing quark chemical potential $`\mu `$ to good approximation by the EoS of a gas of quasiparticles with thermally generated masses, incorporating confinement effectively by a temperature-dependent, reduced number of thermodynamically active degrees of freedom. This method is briefly outlined in this section. For a more detailed discussion the reader is referred to ref..
At very high temperatures, spectral functions for gluons or quarks of the form $`\delta (E^2k^2m^2(T))`$ with $`m(T)gT`$ are found in HTL perturbative calculations. Here, $`E`$ is the particle energy, $`k`$ the absolute value of its momentum, $`m(T)`$ its thermally generated mass and $`g`$ the QCD coupling constant. As long as the spectral function at lower temperatures resembles qualitatively this asymptotic form, a quasiparticle description is expected to be applicable. QCD dynamics is then incorporated in the thermal masses of the quark and gluon quasiparticles. These thermal masses are obtained from the self-energies of the corresponding particles, evaluated at thermal momenta $`kT`$:

$`m_q^2=m_{0q}^2+{\displaystyle \frac{N_c^21}{8N_c}}\left(T^2+{\displaystyle \frac{\mu ^2}{\pi ^2}}\right)G^2(T,\mu ),`$ (1)
$`m_g^2=m_{0g}^2+{\displaystyle \frac{1}{6}}\left[\left(N_c+{\displaystyle \frac{N_f}{2}}\right)T^2+{\displaystyle \frac{3}{2\pi ^2}}{\displaystyle \underset{q}{}}\mu _q^2\right]G^2(T,\mu ).`$ (2)
$`N_f`$ is the number of flavors, $`N_c`$ the number of colors. The effective coupling strength $`G`$ is specified as

$$G(T,\mu =0)=\frac{g_0}{\sqrt{11N_c2N_f}}\left([1+\delta ]\frac{T_c}{T}\right)^\beta .$$
(3)
Setting $`g_0=9.4`$, $`\beta =0.1`$, the effective masses as given in equations (1) and (2), approach the HTL result at high temperatures. (A small shift $`\delta =10^6`$ helps fine-tuning at $`TT_c`$). Because of the existence of a heat bath background, new partonic excitations, plasmons (longitudinal gluons) and plasminos (quark-hole excitations) are also present in the plasma. However, their spectral strengths are exponentially suppressed for hard momenta and large temperatures and consequently these states are essentially unpopulated . The functional dependence of $`m_g(T)`$ on $`T`$ is based on the conjecture that the phase transition is second order or weakly first order which suggests an almost power-like behavior $`m(TT_c)^\beta `$ with some critical exponent $`\beta >0`$. It is assumed that the pseudocritical form of the effective coupling constant given in equation (3) also provides the correct approximate expression for the effective quark mass. This is supported by a non-perturbative dispersion relation analysis for a thermal quark interacting with the gluon condensate .
Close to $`T_c`$ the picture of a non-interacting gas is not appropriate because the driving force of the transition, the confinement process, is not taken into account. Below $`T_c`$, the relevant degrees of freedom are pions and other hadrons. Approaching $`T_c`$ from below, deconfinement sets in and the quarks and gluons are liberated, followed by a sudden increase in entropy and energy density. Conversely, when approaching the phase transition from above, the decrease in the thermodynamic quantities is not primarily caused by increasing masses of the quasiparticles, but by the reduction of the number of thermally active degrees of freedom due to the onset of confinement. For example, gluons begin to form heavy clusters (glueballs), so that the gluon density gets reduced as $`T_c`$ is approached from above. This feature can be incorporated in the quasiparticle picture by modifying the number of effective degrees of freedom by a temperature dependent confinement factor $`C(T)`$:

$$C(T,\mu =0)=C_0\left([1+\delta _c]\frac{T_c}{T}\right).$$
(4)
The confinement factor is taken to be universal. The parameters $`C_0`$, $`\delta _c`$ and $`\beta _c`$ are fixed by reproducing the entropy density that results from lattice QCD thermodynamics. Since the results of lattice calculations with dynamical quarks are still dependent on the details of the simulations, $`C_0`$, $`\delta _c`$ and $`\beta _c`$ should be finetuned for different lattice calculations.
For homogeneous systems of large volume $`V`$, the Helmholtz free energy $`F`$ is related to the pressure $`p`$ by $`F(T,V)=p(T)/V`$. In the present framework of a gas of quasiparticles, its explicit expression reads

$$p(T)=\frac{\nu _g}{6\pi ^2}_0^{\mathrm{}}\text{d}kC(T)f_B(E_k^g)\frac{k^4}{E_k^g}+\frac{2N_c}{3\pi ^2}\underset{i=1}{\overset{N_f}{}}_0^{\mathrm{}}\text{d}kC(T)f_D(E_k^i)\frac{k^4}{E_k^i}B(T).$$
(5)
$`\nu _g`$ is the gluon degeneracy factor, $`E_k^g=\sqrt{k^2+m_g^2(T)}`$ is the gluon energy, $`E_k^q=\sqrt{k^2+m_q^2(T)}`$ the quark energy, $`f_B(E_k^g)=(\mathrm{exp}((E_k^g)/T)1)^1`$ the Bose-Einstein distribution function of gluons and $`f_D(E_k^q)=(\mathrm{exp}((E_k^q)/T)+1)^1`$ the Fermi-Dirac distribution function of quarks. The energy density $`ϵ`$ and the entropy density $`s`$ take the form

$$ϵ(T)=\frac{\nu _g}{2\pi ^2}_0^{\mathrm{}}\text{d}kk^2C(T)f_B(E_k^g)E_k^g+\frac{2N_c}{\pi ^2}\underset{i=1}{\overset{N_f}{}}_0^{\mathrm{}}\text{d}kk^2C(T)f_D(E_k^i)E_k^i+B(T).$$
(6)
and

$$s(T)=\frac{\nu _g}{2\pi ^2T}_0^{\mathrm{}}\text{d}kk^2C(T)f_B(E_k^g)\frac{\frac{4}{3}k^2+m_g^2(T)}{E_k^g}+\frac{2N_c}{\pi ^2T}\underset{i=1}{\overset{N_f}{}}_0^{\mathrm{}}\text{d}kk^2C(T)f_D(E_k^i)\frac{\frac{4}{3}k^2+m_q^2(T)}{E_k^i}$$
(7)
The function $`B(T)`$ is introduced to act as a background field. It is necessary in order to maintain thermodynamic consistency: $`p`$, $`ϵ`$ and $`s=p/T`$ have to satisfy the Gibbs-Duhem relation $`ϵ+p=Ts=Tp/T`$. $`B(T)`$ basically compensates the additional $`T`$-derivatives from the temperature-dependent masses in $`p`$ and thus is not an independent quantity. Since $`B(T)`$ adds to the energy density of the quasiparticles, it can be interpreted as the thermal vacuum energy density. The entropy density, as a measure of phase space, is unaffected by $`B(T)`$.

## 3 Finite chemical potential

The quasiparticle model reviewed in the previous section accurately reproduces lattice thermodynamical quantities such as the pressure, the energy density and the entropy density in the temperature range $`T_c<T\text{ }<4T_c`$ at vanishing chemical potential . However, many physical questions, e.g. the structure of quark cores in massive neutron stars, the baryon contrast prior to cosmic confinement or the evolution of the baryon number in the mid-rapidity region of central heavy-ion collisions, require a detailed understanding of the EoS at non-vanishing quark chemical potential. In this section, a thermodynamically self-consistent extension of the quasiparticle model to finite quark chemical potentials is presented. Results for various observables are then computed and compared to finite $`\mu `$ lattice results in the next section.
At vanishing quark chemical potential, it is conjectured from asymptotic freedom that QCD undergoes a phase transition from the hadronic phase to the QGP phase. At extremely high density, cold quark matter is necessarily in the Color-Flavor-Locked (CFL) phase in which quarks of all three colors and all three flavors form cooper pairs. It is expected that this phase is separated from the hadronic phase by the color superconducting 2SC phase. For a review of the QCD phase diagram, the reader is referred to . Our extension of the quasiparticle model provides a straightforward way to map the EoS at finite temperature and vanishing quark chemical potential into the $`T\mu `$ plane without further assumptions. However, since this continuous mapping relies on quark and gluon quasiparticles, it cannot provide information about other possible phases with a different (quasiparticle) structure. It is therefore applicable in a limited range of not too large chemical potentials.
The pressure of an ideal gas of quark and gluon quasiparticles with effective masses depending on temperature and quark chemical potential, is given by

$$p(T,\mu )=\frac{\nu _g}{6\pi ^2}_0^{\mathrm{}}\text{d}kC(T,\mu )f_B(E_k^g)\frac{k^4}{E_k^g}+\frac{N_c}{3\pi ^2}\underset{q=1}{\overset{N_f}{}}_0^{\mathrm{}}𝑑kC(T,\mu )[f_D^+(E_k^q)+f_D^{}(E_k^q)]\frac{k^4}{E_k^q}B(T,\mu ),$$
(8)
with $`f_D^\pm (E_k^q)=(\mathrm{exp}((E_k^q\mu )/T)+1)^1`$. The effective coupling strength $`G(T,\mu )`$, the confinement factor $`C(T,\mu )`$ and the mean field contribution $`B(T,\mu )`$ now also depend on the quark chemical potential $`\mu `$. $`B(T,\mu )`$ is calculated in appendix A. The quark number density (which is related to the baryon number density $`n_B`$ by $`n_q=n_B/3`$) retains the ideal gas form

$$n_q(T,\mu )=\frac{N_c}{\pi ^2}\underset{q=1}{\overset{N_f}{}}_0^{\mathrm{}}𝑑kC(T,\mu )[f_D^+(E_k^q)f_D^{}(E_k^q)]k^2,$$
(9)
but with the confinement factor $`C(T,\mu )`$ included.
In the previous section expressions for the coupling $`G(T,\mu =0)`$ and the confinement factor $`C(T,\mu =0)`$ are given. These expressions can be generalized to finite chemical potential in a thermodynamically self-consistent way using Maxwell relations. Imposing the Maxwell relation between the derivatives of the quark number density and the entropy,

$$\frac{s}{\mu }|_T=\frac{n}{T}|_\mu \underset{i}{}\left(\frac{n}{m_i^2}\frac{m_i^2}{T}\frac{s}{m_i^2}\frac{m_i^2}{\mu }\right)=0\text{and}\left(\frac{n}{C}\frac{C}{T}\frac{s}{C}\frac{C}{\mu }\right)=0,$$
(10)
yields a set of first oder quasilinear partial differential equations for the effective coupling constant $`G^2(T,\mu )`$ and the confinement factor $`C(T,\mu )`$:

$`a_T(T,\mu ;G^2){\displaystyle \frac{G^2}{T}}+a_\mu (T,\mu ;G^2){\displaystyle \frac{G^2}{\mu }}=b(T,\mu ;G^2),`$ (11)
$`c_T(T,\mu ;G^2){\displaystyle \frac{C}{T}}+c_\mu (T,\mu ;G^2){\displaystyle \frac{C}{\mu }}=0.`$ (12)
The coefficients $`a_T`$, $`a_\mu `$, $`b`$, $`c_T`$, $`c_\mu `$ depend on $`T`$, $`\mu `$, $`G^2`$ but not on $`C`$. It can be solved by the method of characteristics (see appendix B). The flow of the effective coupling and the confinement factor is elliptic. In particular, one finds

$$a_T(T,\mu =0)=0,a_\mu (T=0,\mu )=0,c_T(T,\mu =0)=0,c_\mu (T=0,\mu )=0.$$
(13)
Therefore, the characteristics are perpendicular to both the $`T`$ and the $`\mu `$ axis. This guarantees that specifying the coupling constant and the confinement factor on the $`T`$ axis sets up a valid initial condition problem. Plots of the characteristic curves and the confinement factor are shown in figure 1 and 2.

## 4 Comparison with lattice results

Simulations of QCD at finite chemical potential are extremely difficult because the fermion determinant becomes complex. This prohibits Monte Carlo importance sampling, which interprets the measure as a probability factor and thus requires it to be positive. While this problem remains unsolved, there are some approaches which circumvent the sign problem and allow lattice simulations for small chemical potentials $`\mu \text{ }<T_c`$. A review comparing these methods in detail can be found in .

### 4.1 The phase boundary line

In the case of vanishing chemical potential, universal arguments and lattice simulations suggest a phase transition from the hadronic phase to the QGP phase at a critical temperature $`T_c`$. For QCD with three light flavors $`m_um_dm_s5`$ MeV this transition is expected to be first order. For two light flavors $`m_um_d5`$ MeV and an infinitely large $`m_s`$ there is no phase transition, only a smooth crossover . This suggests there is a critical strange mass $`m_s^c`$ at which one finds a second order phase transition. Lattice calculations indicate that $`m_s^c`$ is about half of the physical mass $`m_s`$. At finite quark chemical potential $`\mu `$ and vanishing $`T`$ a first order phase transition is predicted. For the physical $`m_s`$ this implies that there is a first order phase transition for small $`T`$ and large $`\mu `$ which ends at a critical point $`(T^{},\mu ^{})`$. At this point the phase transition is of second order. For large $`T`$ and small $`\mu `$ the two phases are separated by a crossover. We refer to the line $`T_c(\mu )`$ that separates the hadronic phase from the QGP phase as the “phase boundary line”. In the literature  this line is also frequently called the “pseudocritical line”. $`T_c(\mu )`$ has been calculated on the lattice for $`N_f=4`$ , $`N_f=2`$  and $`N_f=3`$  flavors of quarks up to quark chemical potentials $`\mu `$ of order $`T_c`$. In the following we focus on the three-flavor results where the critical line has been calculated with an accuracy up to terms of order $`(\mu /T)^6`$. There, a Wilson gauge action and three degenerate flavors of staggered quarks have been employed, with bare masses in the range $`0.025<am<0.04`$, where $`a`$ denotes the lattice spacing. The finite volume scaling behavior was monitored by using three lattice sizes, $`8^3\times 4`$, $`10^3\times 4`$ and $`12^3\times 4`$.
In our quasiparticle model, the sudden decrease of the pressure, the energy density, the quark number density and the entropy density caused by gluons and quarks getting trapped in glueballs and hadrons when $`T_c`$ is approached from above, is parametrized by the confinement factor $`C(T,\mu )`$. Consequently, it is natural to relate the critical line to the characteristic curve of the confinement factor through $`T_c(\mu )`$, as long as $`\mu `$ is small and the nature of the quasiparticles does not change qualitatively.
In order to calculate the confinement factor at finite chemical potential, we need to specify a valid initial condition, e.g. $`C(T,\mu =0)`$. The functional form of $`C(T,\mu =0)`$ is set by eq.(4). We have employed the following set of parameters, as found in ref.:

|  | $`C_0`$ | $`\delta _c`$ | $`\beta _c`$ |
| --- | --- | --- | --- |
| 3 flavors | 1.03 | 0.02 | 0.2 |

We have checked that the form of the phase boundary line in the quasiparticle model depends only weakly on the exact choice of parameters and a small difference only shows up for values much larger than the range of $`\mu `$ covered by the lattice simulations. The lattice phase boundary line and our result is shown in figure 3.

The quasiparticle result is within the lattice estimate for $`\mu _B\text{ }<2.5T_c`$ and deviates only slightly from the lattice result for larger chemical potentials.

### 4.2 Thermodynamical quantities

There have been lattice calculations of thermodynamical quantities at finite chemical potential for $`N_f=2+1`$  and $`N_f=2`$  flavors of quarks. In the following we focus on results from  where a p4-improved staggered action on a $`16^3\times 4`$ lattice was used. There, the $`N_\tau `$ dependence is known to be small, in contrast to standard staggered fermion actions which show substantially larger cut-off effects. Estimates of the pressure, the quark number density and associated susceptibilities as functions of the quark chemical potential were made via a Taylor series expansion of the thermodynamic grand canonical potential $`\mathrm{\Omega }`$ up to fourth order.
To calculate thermodynamical quantities within the quasiparticle model, we need to fix the parameters of the effective coupling constant and the confinement factor. Our calculations have shown that the results are not sensitive to the detailed choice of parameters for the effective coupling $`G`$. We have therefore used the parameters from ref. in our calculations. In principle, the parameters of the confinement factor can be fixed by comparing our calculations to lattice results at vanishing chemical potential. However, in ref. no $`\mu =0`$ lattice data is given. Since lattice calculations including quarks give slightly different results depending on which action has been used, fitting the parameters in $`C(T,\mu =0)`$ by comparing quasiparticle results to lattice data from a different simulation is not feasible and would lead to large differences. Consequently, we directly used the finite $`\mu `$ lattice results for fitting. Good agreement with the lattice thermodynamical observables was found for the following sets of parameters:

|  | $`C_0`$ | $`\delta _c`$ | $`\beta _c`$ |
| --- | --- | --- | --- |
| Set A | 1.05 | -0.016 | 0.15 |
| Set B | 1.12 | 0.02 | 0.2 |

While set A reproduces the lattice thermodynamical results slightly better, set B is in better agreement with the parameters found in  for $`\mu =0`$ lattice simulations.
The temperature dependence of the normalized pressure difference $`\mathrm{\Delta }p(T,\mu )=(p(T,\mu )p(T,\mu =0))/T^4`$ is shown in figure 4 and that of the normalized quark number density $`n_q(T,\mu )/T^3`$ in figure 5.

Whereas the computation of the quark number density from equation (9) is straightforward, a numerical evaluation of (8) is difficult because of the derivatives of the effective masses and the confinement factor in $`B(T,\mu )`$ (see expressions in appendix A). It turns out that it is simpler to calculate the pressure difference using the following relation:

$$\mathrm{\Delta }p(T,\mu )=\frac{1}{T^4}_0^\mu \text{d}\mu ^{}n_q(T,\mu ^{}).$$
(14)
The lattice pressure difference is well reproduced even for the largest values of the chemical potential. The quark number density is in very good agreement with the lattice data for $`\mu /T_c=0.2`$ and $`0.4`$. For larger values of $`\mu `$ our calculations underestimate the magnitude of the lattice results close to $`T_c`$, but show the same qualitative features.

## 5 Summary

We have presented a quasiparticle description of the QCD EoS at finite temperature and quark chemical potential. Our main modification as compared to previous work is the inclusion of finite quark chemical potential in a thermodynamically consistent way. We have first reviewed our improved quasiparticle model which schematically includes confinement. We have shown how Maxwell relations can be used to construct the effective coupling $`G(T,\mu )`$ and the confinement factor $`C(T,\mu )`$ at finite chemical potential. We then used this model to calculate the phase boundary line $`T_c(\mu )`$ and the normalized pressure difference $`\mathrm{\Delta }p(T,\mu )=(p(T,\mu )p(T,\mu =0))/T^4`$ and the normalized quark number density $`n_q(T,\mu )/T^3`$. We compared our results to recent lattice calculations and found remarkably good agreement even for large quark chemical potentials $`\mu T_c`$.

## Appendix A Calculation of $`B(T,\mu )`$

The “background field” quantity $`B(T,\mu )`$ appearing in eq.(8) can be obtained from the Gibbs-Duhem relation

$$ϵ+p=Ts+\mu n=T\frac{p}{T}+\mu \frac{p}{\mu }.$$
(15)
The left hand side reads:

$$ϵ+p=\frac{N_cN_f}{3\pi ^2}_0^{\mathrm{}}𝑑k[f_D^++f_D^{}]C(T,\mu )k^2\left(\frac{4k^2+3m_q^2}{E_k}\right).$$
(16)
To evaluate the right-hand side, derivatives of $`f_D^\pm (E_k^q)`$ with respect to $`T`$ and $`\mu `$ are rewritten as derivatives with respect to $`k`$. After an integration by parts, the first term on the right-hand side reads

$`T{\displaystyle \frac{p}{T}}`$ $`=`$ $`T{\displaystyle \frac{N_cN_f}{3\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑kf_D^+\left({\displaystyle \frac{C}{T}}{\displaystyle \frac{k^4}{E_k}}C(T,\mu ){\displaystyle \frac{3k^2}{2E_k}}{\displaystyle \frac{m_q^2}{T}}+C(T,\mu )(E_k+\mu )+C(T,\mu ){\displaystyle \frac{k^4}{TE_k}}\right)`$ (17)
$`+`$ $`T{\displaystyle \frac{N_cN_f}{3\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑kf_D^{}\left({\displaystyle \frac{C}{T}}{\displaystyle \frac{k^4}{E_k}}C(T,\mu ){\displaystyle \frac{3k^2}{2E_k}}{\displaystyle \frac{m_q^2}{T}}+C(T,\mu )(E_k\mu )+C(T,\mu ){\displaystyle \frac{k^4}{TE_k}}\right)`$
$``$ $`T{\displaystyle \frac{B(T,\mu )}{T}},`$
and the second term is given by

$`\mu {\displaystyle \frac{p}{\mu }}`$ $`=`$ $`\mu {\displaystyle \frac{N_cN_f}{3\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑kf_D^+\left({\displaystyle \frac{C(T,\mu )}{\mu }}{\displaystyle \frac{k^4}{E_k}}C(T,\mu ){\displaystyle \frac{m_q^2}{\mu }}{\displaystyle \frac{3k^2}{2E_k}}C(T,\mu )3k^2\right)`$ (18)
$`+`$ $`\mu {\displaystyle \frac{N_cN_f}{3\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑kf_D^{}\left({\displaystyle \frac{C(T,\mu )}{\mu }}{\displaystyle \frac{k^4}{E_k}}C(T,\mu ){\displaystyle \frac{m_q^2}{\mu }}{\displaystyle \frac{3k^2}{2E_k}}+C(T,\mu )3k^2\right)\mu {\displaystyle \frac{B(T,\mu )}{\mu }}.`$
Substituting (16), (17) and (18) in the Gibbs-Duhem relation yields a partial differential equation of the type

$$x\frac{f(x,y)}{x}+y\frac{f(x,y)}{y}=(x,y).$$
(19)
It has the general solution

$$f(x,y)=^x𝑑t(t,\frac{y}{x}t)+\left(\frac{y}{x}\right).$$
(20)
Here, $`(y/x)`$ is a solution of the homogeneous equation. Returning to our case, $`(\mu /T)`$ becomes an arbitrary function of the ratio $`\mu /T`$ to be fixed by boundary conditions. For $`\mu 0`$, $`(\mu /T)`$ does not depend on $`T`$ anymore and therefore has to be identified with an integration constant $`B_0`$. Provided that $`(\mu /T)`$ is a continuous function it must be close to $`B_0`$ for small $`\mu /T`$. The first term in a Taylor expansion of $`(\mu /T)`$ vanishes and the series starts only at order $`(\mu /T)^2`$. Therefore we identify $`(\mu /T)`$ with the constant $`B_0`$ for all $`\mu `$ under consideration. Assembling all pieces, the final result reads

$`B(T,\mu )`$ $`=`$ $`B_1(T,\mu )+B_2(T,\mu )+B_0,`$
$`B_1(T,\mu )`$ $`=`$ $`{\displaystyle \frac{N_cN_f}{3\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑k{\displaystyle _{T_c}^T}𝑑\tau \left[f_D^+(E_k^q)+f_D^{}(E_k^q)\right]\left({\displaystyle \frac{C}{\tau }}+{\displaystyle \frac{\mu }{T}}{\displaystyle \frac{C}{\left(\frac{\mu }{T}\tau \right)}}\right){\displaystyle \frac{k^4}{E_k^q}},`$
$`B_2(T,\mu )`$ $`=`$ $`{\displaystyle \frac{N_cN_f}{2\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑k{\displaystyle _{T_c}^T}𝑑\tau C\left[f_D^+(E_k^q)+f_D^{}(E_k^q)\right]\left({\displaystyle \frac{m_q^2}{\tau }}+{\displaystyle \frac{\mu }{T}}{\displaystyle \frac{m_q^2}{\left(\frac{\mu }{T}\tau \right)}}\right){\displaystyle \frac{k^2}{E_k^q}},`$ (21)
where the explicit $`\tau `$-dependence in $`C(\tau ,\mu /T\tau )`$, $`m_q(\tau ,\mu /T\tau )`$ and $`E_k^q(\tau ,\mu /T\tau )`$ has been suppressed for the sake of lucidity.

## Appendix B Method of characteristics

Equations (11) and (12) are a set of coupled quasilinear first order partial differential equations for the effective coupling constant $`G^2(T,\mu )`$ and the confinement factor $`C(T,\mu )`$. Equation (11) does not depend on $`C(T,\mu )`$. Thus we can first solve this equation for $`G^2(T,\mu )`$ and insert the result in equation (12).
The usual method found in textbooks is to reduce a quasilinear partial differential equation of the form

$$a_T(T,\mu ;X)\frac{X}{T}+a_\mu (T,\mu ;X)\frac{X}{\mu }=c(T,\mu ;X)$$
(22)
to a system of coupled ordinary differential equations,

$$\frac{dT(s)}{ds}=a_T,\frac{d\mu (s)}{ds}=a_\mu ,\frac{dX(s)}{ds}=c.$$
(23)
This determines the characteristic curves $`T(s)`$, $`\mu (s)`$, and the evolution of $`X`$ along such a curve, given an initial value. However, this method is not well suited for numerical use which is necessary for non-trivial $`a_T`$, $`a_\mu `$ and $`c`$. Rewriting equation (22) as

$$a_T\left(\frac{dX}{dT}\frac{X}{\mu }\frac{d\mu }{dT}\right)+a_\mu \frac{X}{\mu }=c\frac{X}{\mu }\left(\frac{a_\mu }{a_T}dTd\mu \right)=\frac{c}{a_T}dTdX,$$
(24)
we find the equation $`a_\mu dTa_Td\mu =0`$ for the characteristics and $`cdTa_TdX=0`$ for the evolution of $`X`$. These equations can easily be solved numerically.


format: text

text:
# Local Melting and Drag for a Particle Driven Through a Colloidal Crystal

\[

## Abstract

We numerically investigate a colloidal particle driven through a colloidal crystal as a function of temperature. When the charge of the driven particle is larger or comparable to that of the colloids comprising the crystal, a local melting can occur, characterized by defect generation in the lattice surrounding the driven particle. The generation of the defects is accompanied by an increase in the drag force on the driven particle, as well as large noise fluctuations. We discuss the similarities of these results to the peak effect phenomena observed for vortices in superconductors.

\]

Colloidal crystals are an ideal system for studying a variety of issues that arise in two-dimensional (2D) systems such as melting , defect dynamics , and ordering on 2D  and 1D  periodic substrates. A particular advantage of this system is that, due to the particle size, the individual colloid positions and motions can be directly observed, unlike other systems in which this information is generally inaccessible.

Recently there has been growing interest in controlling colloids individually or in small groups by means of optical techniques such as holographic optical trap arrays . With such methods, individual colloids can be captured and moved, or collections of colloids can be driven through single traps or periodic arrays of traps . Alternative methods for driving individual colloids include moving single magnetic particles through assemblies of non-magnetic colloids . These methods of manipulating and driving individual colloids offer a wealth of new ways to explore the dynamical properties of colloidal crystals and glasses, and can also be used to manipulate particles in other systems such as dusty plasmas

A relatively simple example of manipulating single colloids to probe collective properties of a surrounding colloidal crystal is to drive a single colloid through a crystalline array of other colloids. One question that arises in this system is how the size of the driven particle affects both its motion and the response of the surrounding colloids. A very small driven particle is unlikely to generate enough stress in the surrounding lattice to create defects, and therefore the driven particle motion and the response of the system will be elastic. As the particle size or the temperature is increased, defect generation may become possible and a transition to plastic flow can occur. The defects may remain localized near the driven particle and strongly affect the frictional drag. The drag should also depend on the orientation of the driving direction with respect to the colloidal lattice. For example, in a triangular lattice, easy flow directions should occur along 60 angles. Work on systems of particles flowing over different orientations of rigid substrates has shown that a series of magic angles can arise where the flow locks into particular orientations .

By studying colloids moving through a periodic substrate created by other colloids, it may be possible to gain insight into related phenomena such as the role played by dislocations in friction and depinning. For example, in the case of atomic friction  where one lattice is driven over another, dislocations may nucleate in either of the lattices and can strongly affect the friction coefficient. In the case of driven vortices in superconductors, the onset of dislocations or plasticity at or near melting can give rise to a phenomena known as the peak effect, where the vortices slow considerably or become pinned as a function of external magnetic field or temperature .

In order to investigate the dynamics of a single driven colloidal particle moving through a triangular colloidal lattice we conduct a series of Langevin simulations of the type employed in Refs. . We consider a 2D system of $`N`$ colloidal particles interacting with a repulsive Yukawa potential in the absence of a substrate. We impose periodic boundary conditions in the $`x`$ and $`y`$ directions. The initial configuration of the colloids is a triangular lattice. For a fixed density $`n`$ there is a well defined melting temperature $`T_m`$, measured either by diffusion or by the density of non-sixfold coordinated particles as determined by a Voronoi tessellation. An additional colloid is placed in the system, and a constant driving force $`f_d`$ is applied only to that particle.

The equation of motion for a colloid $`i`$ is

$$\frac{d𝐫_i}{dt}=𝐟_{ij}+𝐟_d+𝐟_T$$
(1)
Here $`𝐟_{ij}=_{ji}^{N_c}_iV(r_{ij})`$ is the interaction force from the other colloids. The colloid-colloid interaction is a Yukawa or screened Coulomb potential given by $`V(r_{ij})=(q_iq_j/|𝐫_i𝐫_j|)\mathrm{exp}(\kappa |𝐫_i𝐫_j|)`$, where $`q_{i(j)}`$ is the particle charge, $`1/\kappa `$ is the screening length, and $`𝐫_{i(j)}`$ is the position of particle $`i(j)`$. All the particles have the same charge $`q`$ except for the driven particle which has charge $`q_d`$ that may differ from $`q`$. The applied driving force $`𝐟_d`$ is finite if colloid $`i`$ is the driven colloid and zero otherwise. The system size is measured in units of the lattice constant $`a_0`$ and we use a screening length of $`\kappa =3/a_0`$. The thermal force $`𝐟_T`$ is a randomly fluctuating force from random kicks. For most of the results in this work we fix the colloid density (thus fixing the melting temperature

$`T_m`$) and vary $`T`$ and $`q_d`$. For results with constant $`f_d`$ we use the same drive for all $`q_d`$. In the absence of any other colloids the driven particle moves at a velocity $`V_0`$.

We first consider the case where a particle with $`q_d/q=3`$ is driven along the zero degree angle with respect to the lattice. Three distinct phases occur for increasing $`T/T_m`$. In Fig. 1(a,c,e) we show the colloids and particle trajectories for different $`T/T_m`$. In Fig. 1(b,d,f) we show the corresponding Voronoi construction or Wigner-Seitz cells, in which a single colloid is located at the center of each polygon. The cells are colored according to the number of nearest neighbors: white for six neighbors, black for five neighbors, and gray for seven neighbors. In Fig. 1(a,b) at $`T/T_m=0.1`$, the flow is elastic and the particle moves along a 1D path while causing small distortions in the surrounding lattice. Fig. 1(b) shows that although a group of defects surrounds the driven particle, which is an extra or interstitial particle in the

triangular lattice, there are no other dislocations induced elsewhere in the lattice. At $`T/T_m=0.3`$, we find a transition to plastic flow where the particle no longer moves strictly in a straight line, and defects form in the surrounding lattice. Additionally, some colloids from the surrounding lattice are pushed in front of the driven particle over distances larger than a lattice constant $`a_0`$. Fig. 1(c,d) illustrates this behavior for $`T/T_m=0.57`$. A localized molten region surrounds the moving particle, and ejects 5-7 defects into the lattice over time. Fig. 1(e,f) shows the system above melting at $`T/T_m=1.07`$, where defects are generated in the sample even in the absence of the driven particle. In this globally molten phase, the driven particle is still able to push other particles in front of it.

We next show the effect of the defect generation on the colloidal velocity in the different phases. In Fig. 2(a) we plot the percentage of sixfold coordinated particles, $`P_6`$, for a system with $`q_d/q=1`$, fixed $`f_d`$, and increasing temperature. We find a transition from elastic to plastic flow at a temperature close to $`T/T_m`$. In Fig. 2(b), we see that the corresponding colloid velocity $`V_x<V_0`$ even for $`T=0`$, indicating that some of the energy of the driven particle is absorbed by vibrations in the surrounding lattice. As $`T`$ increases, $`V_x`$ decreases slightly, culminating in a sharp drop at $`T/T_m=1`$ which coincides with the drop in $`P_6`$. We also find a noticeable decrease in $`V_x`$ at temperatures slightly below $`T_m`$ which is correlated with a decrease in $`P_6`$, indicating that the defects in the colloidal lattice slow the driven particle. As $`T`$ is further increased, $`V_x`$ eventually increases to $`V_0`$ for $`T/T_m>10.`$

We next consider the effect of varying $`q_d`$. In Fig. 2(c) we plot $`P_6`$ for a system with $`q_d/q=0.25`$, and in Fig. 2(d) the corresponding $`V_x`$. Here $`V_x`$ is only slightly lower than $`V_0`$. A small decrease in $`V_x`$ occurs close to melting; however, unlike the case of $`q_d/q=1`$, there is little change in $`V_x`$ across $`T_m`$. In Fig. 2(e) we show $`P_6`$ for $`q_d/q=6`$, and in Fig. 2(f) the corresponding $`V_x`$. $`P_6`$ is more rounded than for $`q_d/q=0.25`$. This is due to the dislocations created by the driven particle, which become more numerous as the melting transition is approached. $`V_x`$ drops dramatically well below the melting temperature, at a much lower $`T`$ than observed for smaller $`q_d`$. The decrease in $`V_x`$ coincides with the creation of dislocations in the lattice near the driven particle. Far from the driven particle, the colloidal lattice remains ordered.

We observe that at the onset of the local dislocation or defect creation, the driven particle often drags an additional colloid in front of it. Once the dragged colloid falls away to one side, another colloid is captured and dragged. The presence of an extra dragged particle produces extra resistance to the motion of the driven particle. In Fig. 2(b), at high temperatures $`TT_m`$, the drag effect is reduced and the particle velocity increases when the thermal fluctuations become so strong that the driven particle can no longer capture other particles. Small driven particles do not have a strong enough colloid-colloid interaction to push any additional particles, so no additional drag effect occurs. The small driven particle of Fig. 2(c) creates little distortion in the surrounding colloidal lattice, so the formation of dislocations occurs only due to $`T`$. In Fig. 2(e) the driven particle is large enough to create a significant amount of distortion in the surrounding lattice, and thus produce dislocations at low $`T`$ well below melting. Thus, the combination of lattice distortions from the driven particle and the temperature produces the onset of dislocations. Dislocations appear at a lower temperature in the presence of a large driven particle than they do for a small or for no particle, when no lattice distortions are introduced. We therefore expect that as $`q_d`$ increases, the local melting transition will shift to lower $`T`$. We have also considered the effect of varying $`f_d`$, and observe the same results whenever the particle velocity is slower than the propagation of disturbances in the surrounding media.

In Fig. 3 we plot the curve separating the crystal phase C, the local melting phase LM, and the global melting phase M. The line separating C and LM is obtained from the decrease of $`V_x`$. Within LM, the moving particle creates dislocations but only locally as seen in Fig. 1(d).

We note that these results have many similarities with a phenomena called the peak effect for vortices in superconductors . An immediate difference between the two systems is that the vortices interact with a background of quenched disorder which is not present in our

system. The peak effect is generally believed to occur due to a transition from an elastic, dislocation free vortex lattice to a disordered glassy state containing many dislocations. The relatively stiff elastic lattice cannot adjust to the disorder and is thus weakly pinned, whereas the soft disordered state can adjust well to the disorder and is strongly pinned. The peak effect is observed by monitoring the vortex velocity in the form of a voltage as a function of temperature. At the peak effect transition the voltage drops abruptly, indicating that the vortices are slowing or being pinned. The peak effect can also be observed by constructing IV curves. There is still controversy over whether the true nature of the peak effect is a thermal melting or a disorder induced transition . Distinguishing between these two scenarios is made difficult by the fact that the entire vortex lattice becomes disordered at the peak effect transition.

In our system there is a transition from an elastic (defect free) flow to a plastic (defected) flow which coincides with a drop in the velocity of the moving particle. The velocity drop occurs even when the defects do not occur in the whole sample but only surround the driven particle. Additionally we find strong $`1/f`$ type velocity noise fluctuations in the disordered flow. In our system, the velocity transition is not caused by thermal fluctuations, but is instead due to the distortions of the surrounding lattice by the driven particle. In the case of the vortices, lattice distortions could arise from the relative motion of vortices trapped at pinning sites. If the distortions in the surrounding lattice are large enough, dislocations nucleate and the overall vortex velocity drops.

To further explore similarities between the peak effect and our system, we slowly increase $`f_d`$ and measure the velocity at fixed $`T`$, in analogy with IV measurements. A

particular characteristic of the peak effect is crossing $`IV`$ curves . Due to thermal activation, the more disordered system has a lower threshold for motion than the ordered system. In contrast, at higher drives when thermal creep is negligible, the overall velocity is lower in the disordered system than in the ordered system. This implies that the velocity-force curves for the two cases must cross at intermediate drives.

In Fig. 4 we plot the velocity vs force curves for $`q_d/q=3`$ at fixed $`T/T_m=0.5`$ in the ordered regime (solid line) and at $`T/T_m=1.1`$ just into the liquid regime (dotted line). For the disordered (melted) system there is no true depinning threshold and we observe a finite velocity almost to zero applied drive. The velocity in the ordered regime for $`f_d<0.2`$ is less than that of the disordered regime and is almost zero, with a clear depinning threshold. By contrast, at higher drives the velocity in the ordered regime is higher than the melted regime. Thus we observe a crossing of the velocity force curves similar the crossing of the IV curves  for the peak effect. We have also examined the velocity in the local melting case and again find a finite depinning threshold; however, here the velocity is always lower than the ordered case, even at high drives. Our results lend support to the idea that the peak effect phenomena in superconductors is not a thermal melting phenomena, but is instead a disorder induced transition.

We have also explored the effect of driving the particle at various angles with respect to the undriven lattice. For certain angles such as $`60^{}`$ the results are unchanged. For incommensurate angles, the plastic flow transition shifts to a lower $`T/T_m`$ since the moving particle generates larger distortions in the surrounding lattice and defect creation is easier.

In conclusion, we have studied a single colloid driven through an ordered colloidal crystal. We find that as a function of temperature and charge of the driven colloid, there is a transition from elastic flow, where no defects are generated in the surrounding lattice, to plastic flow, where defects proliferate. This transition coincides with an increased drag on the driven particle, due to the driven particle trapping and pushing surrounding colloids in front of it. As the charge of the driven particle increases, the elastic-plastic transition shifts to lower temperatures. This system has several similarities to the peak effect phenomena found in vortex matter in superconductors: The onset of defect formation slows the particle, and there is a crossing in the velocity force curves. Our results suggest that the peak effect is due to a disorder induced transition rather than a thermal melting transition. These results should be accessible experimentally by driving single colloids with optical tweezers or dragging a magnetic bead though an ordered colloidal crystal. Another variation on this would be to drive a colloidal crystal past obstacles of varied size. Other systems in which these effects could be observed include driving particles through dusty plasma crystals or ordered foams.

We thank E. Weeks and M.B. Hastings for useful discussions. This work was supported by the U.S. DoE under Contract No. W-7405-ENG-36.


text:
# The Role of Gas in the Merging of Massive Black Holes in Galactic Nuclei. I. Black Hole Merging in a Spherical Gas Cloud.

## 1 Introduction

Massive black holes (MBHs) are believed to be the power sources for active galactic nuclei, turning mass into energy through accretion (Zel’dovich 1964; Salpeter 1964). The dynamical evidence suggests that every galaxy with a significant bulge hosts an MBH at its center (Richstone et al. 1998). Moreover, the data show that the mass of the central black hole is tightly correlated with the velocity dispersion of the host galaxy’s bulge, yielding the so-called $`m_{BH}\sigma _c`$ relation (Ferrarese & Merritt 2000; Gebhardt et al. 2000). The estimated masses of the central MBHs in galaxies range from $`10^6`$ to well above $`10^9`$ M, making them by far the most massive known single objects. However, it remains an open question how these black holes gain their large masses, and why their masses are strongly correlated with the dynamical properties of their host galaxies.

Three possible scenarios have been suggested for massive black hole formation: (1) They might form by the ‘monolithic collapse’ of a massive gas cloud to form a single MBH (Rees 1984). This scenario would require a mechanism to maintain the gas at a high temperature in order to prevent its fragmentation into many low-mass objects. (2) MBHs might form by a series of mergers of smaller black holes, which might be stellar-mass black holes or, perhaps more interestingly, ‘intermediate-mass black holes’ (IMBHs) formed as the remnants of very massive first-generation stars (Bromm, Coppi, & Larson 1999, 2002; Abel, Bryan, & Norman 2002) or by runaway stellar mergers in very dense environments (Portegies Zwart & McMillan 2002.) (3) MBHs might form by the runaway accretional growth of smaller seed black holes in the dense nuclear regions of galaxies. It is also possible that both gas accretion and the merging of smaller black holes into larger ones, possibly associated with a succession of galaxy mergers, might be involved in forming MBHs (Kauffmann & Haehnelt (2000); Haehnelt 2003; Volonteri, Haardt & Madau 2003). In the latter case, it is necessary for black hole mergers to occur on a timescale shorter than that of the galaxy mergers. It remains unclear, however, just how the central MBHs in galaxies are built up.

The possibility of black hole mergers in galactic nuclei was first considered by Begelman, Blandford, & Rees (1980) in a study of the long-term evolution of a black hole binary at the center of a dense stellar system. They showed that the dynamical evolution of such a binary MBH can be divided into three phases: (a) The two black holes sink toward the center of the system because of the gravitational drag or ‘dynamical friction’ created by their interaction with the stellar background (Chandrasekhar 1943; Binney & Tremaine 1987). (b) The orbit of the resulting binary MBH continues to shrink as 3-body interactions with the surrounding stars continue to extract energy from the orbit. (c) If the binary MBH becomes close enough for gravitational radiation to be important, this becomes an efficient mechanism for further angular momentum loss and quickly drives a merger of the two MBHs. The occurrence of a merger is however not ensured because 3-body interactions between the binary MBH and the surrounding stars tend to eject stars from the central region of the galaxy, depleting the phase-space ”loss cone” of stars that can continue to interact with the binary; this may prevent a transition between phases (b) and (c), causing the merger to stall. More recently this problem has been explored numerically by Makino & Ebisuzaki (1996) and Milosavljevic & Merritt (2001, 2003); these authors found the coalescence to stall because of the ejection of stars from the nucleus, although a definitive statement cannot be made because of lack of sufficient numerical resolution.

However, this picture does not include the possible role of gas in driving the evolution of a binary MBH in a galactic nucleus. Observational and theoretical work both indicate that large amounts of gas can be present in the central regions of interacting galaxies, and that this gas can be a dominant component of these regions. Numerical simulations show that in a merger of galaxies containing gas, much of the gas can be driven to the center by gravitational torques that remove angular momentum from the shocked gas (Barnes & Hernquist 1992, 1996; Mihos & Hernquist 1996); as a result, more than 60% of the gas originally present in the merging galaxies can end up in a massive central concentration or cloud with a diameter of several hundred parsecs. Observations of gas-rich interacting galaxies such as the ‘Ultraluminous Infrared Galaxies’ (ULIRGs) confirm that their central regions often contain massive and dense clouds of molecular and atomic gas whose masses are comparable to the total gas content of a large gas-rich galaxy (Sanders & Mirabel 1996). CO data shows that most of the molecular gas is in relatively smooth interclump medium with dense clumps that account for less than half of the total mass (Downes, Solomon & Radford 1993; Downes & Solomon 1998). These massive central gas clouds must have an important influence on the evolution of any central MBH binary that forms following a galaxy merger. Because the gas is strongly dissipative, unlike the stars, it might be expected to remain concentrated near the center and thus to play a continuing role in driving the evolution of a central binary MBH.

So far, the study of the role gas in driving the evolution of a binary MBH has focused on the interaction of the binary with an accretion disk, where the total gaseous mass is typically smaller than the binary mass. The interaction with an accretion disk can be effective in driving the orbital decay for high mass ratios ($`\mathrm{M}_1>>\mathrm{M}_2`$; Armitage & Natarajan 2002; Gould & Rix 2000). For this mechanism to be effective, a disk mass of at least the same order of magnitude as the secondary MBH is needed. For low mass ratios ($`\mathrm{M}_1\mathrm{M}_2`$) this mechanism is less effective because strong tidal and/or resonant forces create a circumbinary gap (Lin & Papaloizou 1979; Goldreich & Tremaine 1980; Artymowicz & Lubow 1994,1996). The case of comparable masses is most relevant for a major galactic merger and the case of our interest.

In this paper, we study numerically the role of a massive gas cloud, like those seen in the central regions of ULIRGs, in driving the evolution of a binary MBH. We follow the evolution of the binary through many orbits and close to the point where gravitational radiation becomes important. As in the stellar case, the evolution can be divided into three regimes. First, each black hole is driven closer to the center by the gravitational drag produced by its interaction with the gaseous background, similar to the stellar case. After a central binary MBH has formed it dominates the dynamics in the region, the orbit of the binary continues to shrink because of its gravitational interaction with a circumbinary gaseous envelope; Unlike the stellar case, this envelope remains concentrated around the binary and is not ejected. Finally, when the holes come close enough together, gravitational radiation becomes important and causes rapid merging of the binary. Here we present results for a relatively simple idealized case in which the gas is assumed to be in a spherical cloud supported by a high virial temperature, so that the gas retains a nearly spherical and relatively smooth distribution. Our model studies the effect of the relatively smooth medium that accounts for most of the molecular mass in the center of a typical ULIRG (Downes, Solomon & Radford 1993; Downes & Solomon 1998). Subsequent papers will present results of calculations in progress for cases in which the gas is cooler and has a more flattened or clumpy spatial distribution with rotation.

We start our analysis validating the use of the code in §2 by performing several tests. In §3 we perform a numerical study of the differences between the dynamical friction in a gaseous background and the dynamical friction in the stellar case. In §4 we study the evolution of a binary MBH in a gaseous sphere. Finally in §5 we apply our results to the inner regions of a ULIRG, showing that the merging of two massive black holes in such a system can occur in a time as short as $`10^7`$ years.

## 2 Description of the Code

We use the Smooth Particle Hydrodynamics (SPH) code called GADGET (Springel, Yoshida & White 2001) to study the dynamical evolution of a black hole binary in a gaseous medium. GADGET has been successfully tested in a broad range of problems such as galaxy interactions and collisions, cosmological studies of large scale structure and high resolution simulations of the formation of galaxy clusters (Springel, Yoshida & White 2001). The code has many special features such as shear-reduced artificial viscosity, efficient cell-opening criteria, individual time-steps of arbitrary size for all particles, good scalability and an efficient integrator in the linear regime of gravitational clustering. All of these features make GADGET suitable for our work.

The SPH code allows one to study both collisional and collisionless systems. Stars and MBHs are modeled as a self-gravitating collisionless fluid governed by the Vlasov Equation (Binney & Tremaine 1987, eq 4-13c)

$$\frac{df}{dt}=\frac{f}{t}+𝐯\frac{f}{𝐫}\frac{\mathrm{\Phi }}{𝐫}\frac{f}{𝐯}=0,$$
(1)
where the gravitational potential $`\mathrm{\Phi }`$ satisfies Poisson’s equation

$$^2\mathrm{\Phi }(𝐫,t)=4\pi Gf(𝐫,𝐯,t)𝑑𝐯,$$
(2)
and $`f(𝐫,𝐯,t)d𝐯`$ is the phase space mass density. It is very difficult to solve this set of coupled differential equations directly. Instead we use the N-body method, a montecarlo approach where the fluid is represented by N particles. The accuracy in this method depends on the number of particles used.

The gaseous medium is modeled as an ideal, viscid gas, governed by the standard laws of fluid dynamics: the equation of continuity

$$\frac{d\rho }{dt}+\rho 𝐯=0,$$
(3)
and the momentum equation

$$\frac{d𝐯}{dt}=\frac{P}{\rho }+𝐚^{\mathrm{𝐯𝐢𝐬𝐜}}\mathrm{\Phi },$$
(4)
where the time derivatives are Lagrangian $`(d/dt=/t+𝐯)`$. To close the set of equations, we use a polytropic equation of state

$$p=K\rho ^\gamma ,$$
(5)
with polytropic index $`\gamma `$ equal to 1 (isothermal) or 5/3 (adiabatic). We solve the equations of fluid dynamics by the Lagrangian technique called SPH (see Monaghan 1992 for a review). In SPH the fluid elements constituting the system are sampled and represented by particles. By design, the SPH technique is free of geometrical restrictions and the resolution is automatically increased in the denser regions, because particles are concentrated in those regions.

In order to validate the proper use of GADGET, we did several additional tests to those performed by the authors of the code (Springel, Yoshida & White 2001); e.g. the pulsation of an adiabatic sphere, the Evrard collapse (Evrard 1988), and the collapse of a primordial rotating cloud with cooling, all of which yielded with satisfactory results. A full description of the implementation of these tests is given in Appendix A.

## 3 Dynamical Friction in a Gaseous Medium

Dynamical friction is the gravitational drag force that a massive body suffers while moving through a collisionless background medium. It has been widely studied in the literature since the concept was introduced by Chandrasekhar in 1943. He computed the deceleration that a massive body suffers as it moves through an infinite homogeneous medium of smaller particles that exchange momentum with the body, ignoring the self-gravity of the medium. Because Chandrasekhar assumed an infinite homogeneous background medium, the formula diverges. An heuristic device was introduced to remove the divergence, the maximum impact parameter $`(\mathrm{b}_{\mathrm{max}})`$, which is related to the size of the system. The combination of this parameter with another that can be interpreted dimensionally as a minimum impact parameter appears in Chandrasekhar’s formula as a Coulomb logarithm (ln($`\mathrm{b}_{\mathrm{max}}`$/$`\mathrm{b}_{\mathrm{min}}`$) = ln$`\mathrm{\Lambda }`$).

Despite the heuristic treatment of the Coulomb logarithm, numerical simulations have shown that Chandrasekhar’s formula provides a remarkably accurate description of the drag that a massive body orbiting in a stellar system experiences(Lin & Tremaine 1983; Bontekoe & van Albada 1987; Cora et al. 1997). The procedure used to compare Chandrasekhar’s formula with simulations is done by evaluating the formula using local quantities at the instantaneous position of the body, leaving ln$`\mathrm{\Lambda }`$ as a parameter that can be adjusted to achieve the best fit.

Less effort has been devoted to the study of the drag force experienced by a gravitating body moving through a gaseous background medium. Ostriker (1999) did an analytical study of the gravitational drag that a massive body suffers while traveling through a infinite uniform gaseous medium. The drag force arises from the gravitational attraction between the perturber and the density wake behind the body. The force is computed using time-dependent linear perturbation theory. Ostriker’s main result was that the gaseous drag is enhanced for supersonic motions (Mach numbers of the gravitating body in the range of 1–2.5), compared to Chandrasekhar’s formula. Additionally, Ostriker found a non-vanishing gravitational drag for subsonic motions, although less efficient than the stellar case. Previous work (Rephaeli & Salpeter 1980) based on a steady state theory had found that the gravitational drag vanishes for subsonic motions due to the front-back symmetry of the perturbed density distribution.

Several groups have numerically computed the gravitational drag force that a massive body suffers while traveling through a infinite uniform gaseous medium, as a by-product of a study of the Bondi-Hoyle accretion problem (see Shima et al. 1985, Ruffert 1996, and references therein), and they found results consistent with Chandrasekhar’s formula. Sanchez-Salcedo & Brandenburg (2001) extensively numerically tested the validity of Ostriker’s formula, in a study of orbital decay analogous to the problem of orbital decay in spherical stellar system where Chandrasekhar’s formula is a good approximation (Lin & Tremaine 1983; Bontekoe & van Albada 1987; Cora et al. 1997). They found that Ostriker’s formula is able to explain successfully the evolution of the perturber when the motion is highly supersonic, although this theory fails if the motion is only slightly supersonic. They also found that the dynamical friction time-scale is a factor of two smaller in the case of a gaseous sphere than in the corresponding stellar system.

In this paper we model the evolution of a binary MBH in a gaseous sphere. This is basically a gravitational drag problem, as long as the evolution of the binary is determined by the gravitational interaction of each black hole with its local surrounding medium. During the later stages, when the binary MBH dominates the gravitational potential around the orbit, the evolution of the binary is determined by the interaction of a binary MBH with a circumbinary gaseous envelope; we call this regime a ”hard binary” in analogy to the stellar dynamics problem (Quinlan 1996, Milosvljevic & Merritt 2001). For both stellar and gaseous backgrounds, the dynamical friction approximation is no longer valid when the MBH’s are at distances where $`M_{\mathrm{binary}}M_{(\mathrm{star},\mathrm{gas})}(r<r_{\mathrm{binary}})`$.

### 3.1 Model Setup

In our model we assume that the background consists of a self-gravitating gas sphere in hydrostatic equilibrium, within which a massive particle is initially orbiting in a circular orbit. The equilibrium density profile is the nonsingular isothermal sphere (NSIS)

$$\rho (r)=\rho _ce^{\mathrm{\Phi }(r)/\mathrm{\Phi }_c},$$
(6)
where $`\mathrm{\Phi }(r)/\mathrm{\Phi }_c`$ is the gravitational potential obtained by solving the isothermal Lane-Emden equation numerically. A Bonnor-Ebert stability analysis was carried out to guarantee that the NSIS is truncated at a radius where the equilibrium is stable. The SPH particles are distributed according the density profile given by Eq 6, by employing a standard comparison rejection method (Press et al. 1992). We relax the system using the technique developed by Lucy (1977) in order to reduce the numerical noise. The units used in this paper are such $`M_{\mathrm{total}}^{\mathrm{gas}}`$ = 2, $`C_\mathrm{S}`$ = 2.5 and the sound crossing time is 1.

We check that our calculations converge as we increase the resolution of the code, running the same problem with different numbers of SPH particles (25000, 50000, 100000, 200000, and 400000). On each run we introduce a massive collisionless particle, the MBH, in the same initial circular orbit and follow its evolution through several orbits. The mass of the MBH is 1% of the total gas mass.

Figure 1 shows the evolution of the radial distance of the MBH from the center of mass of the gaseous system. The curves are the results for the different resolution levels, the upper line has 400,000 particles and lower lines indicate the results for runs with decreasing number of particles. The solution rapidly converges as we increase the number of particles, and has an adequate level of accuracy with 100,000 particles or more.

The interaction between the massive body and the gas produces a density enhancement in the gas sphere as shown in Fig. 2. The colors represent the density enhancement $`[\rho (r,t)/\rho (r,0)]`$ in the plane of the orbit (z=0) for $`\mathrm{t}=1`$. The dashed black curve is the orbit of the MBH between $`\mathrm{t}=0`$ and $`\mathrm{t}=1`$.

The initial circular velocity of the massive body is that predicted for an isothermal sphere: $`v_{\mathrm{BH}}=\sqrt{2}C_\mathrm{S}`$. The supersonic motion produces a spiral shock that propagates outwards with a mildly supersonic speed $`(1.1C_\mathrm{S})`$. A qualitative explanation for the spiral form of the shock is that the MBH accelerates the gas by the reaction to the gravitational drag force. The direction of the resulting impulse is tangential to the orbit and it is assumed that the disturbance propagates in a straight line with a speed of $`1.1C_\mathrm{S}`$. The solid black line shows the locus of the propagated disturbances, showing a good agreement with the spiral density enhancement. This feature is qualitatively the same as that found by Sanchez-Salcedo & Brandenburg (2001) in a similar analysis (see their Fig. 1). However, we have calculated that 90% of the gravitational drag comes from the region within a radius of 0.3 measured from the MBH; therefore the details of the spiral shock do not affect the motion of the MBH.

### 3.2 Comparison of the Dynamical Friction in a Stellar vs a Gaseous Medium

We want to compare the dynamical friction in a gaseous background with that in a stellar medium. The dynamical friction can be expressed as the product of two factors: a factor with the dimensions of force, which is the same in the stellar and gaseous cases, and a dimensionless factor, $`f^{(\mathrm{star},\mathrm{gas})}`$, which depends on the nature of the background medium:

$$F_{\mathrm{DF}}^{(star,gas)}=4\pi \rho \left(\frac{GM_{\mathrm{BH}}}{C_\mathrm{S}}\right)^2\times f^{(\mathrm{star},\mathrm{gas})}\left(\right)$$
(7)
where $``$=$`v_{\mathrm{BH}}`$/$`C_\mathrm{S}`$. In this context, $`C_\mathrm{S}`$ is the sound speed for the gaseous background or the velocity dispersion for the stellar background. In the stellar case, the expression $`f^{(\mathrm{star})}`$ is the following (Chadrasekhar 1943)

$$f^{(\mathrm{star})}=\mathrm{ln}\left(\frac{r_{\mathrm{max}}}{r_{\mathrm{min}}}\right)^2\left[\mathrm{erf}\left(\frac{}{\sqrt{2}}\right)\sqrt{\frac{2}{\pi }}e^{\frac{^2}{2}}\right];$$
(8)
similarly, Ostriker (1999) computed $`f^{(\mathrm{gas})}`$ for both supersonic and subsonic perturbers

$$f_{\mathrm{subsonic}}^{(\mathrm{gas})}=\frac{1}{2^2}\mathrm{ln}\left(\frac{1+}{1}\right)\frac{1}{},$$
(9)
$$f_{\mathrm{supersonic}}^{(\mathrm{gas})}=\frac{1}{^2}\left[0.5\mathrm{ln}\left(\frac{+1}{1}\right)\mathrm{ln}\left(\frac{r_{\mathrm{max}}}{r_{\mathrm{min}}}\right)\right].$$
(10)
In order to numerically evaluate the differences, we set up models with exactly the same dimensional parameters $`(M_{\mathrm{BH}},\rho andv_{\mathrm{BH}}=C_\mathrm{S})`$, for the gaseous and stellar mediums. To do this we start with self-gravitating isothermal spheres that are in hydrostatic equilibrium (for a gaseous system) and in dynamical equilibrium (for a stellar system; §4.4.3b of Binney & Tremaine 1987). The goal is to use these models to compare directly the gravitational drag in a collisionless medium with that in a collisional one.

To investigate the difference between the gaseous and stellar cases, we follow the evolution of a single MBH (with 1% of the total gaseous/stellar mass) in the same initial circular orbit, in both the gaseous and stellar equilibrium configurations. The velocity of the MBH is initially supersonic $`(v_{\mathrm{BH}}1.4C_\mathrm{S})`$ and remains barely supersonic through most of the simulation. We initialize both runs with the same value for $`4\pi \rho \frac{(GM_{\mathrm{BH}})^2}{C_\mathrm{S}^2}`$, thus, any difference in the evolution reflects purely the effects of the term $`f^{(\mathrm{star},\mathrm{gas})}`$. Each of these runs has a resolution of 100,000 stars/gas particles.

Fig. 3 shows the evolution of the radial distance of the MBH from the center of mass of the system, represented by a dashed line in the stellar case and by the continuous line in the gaseous medium. The distance is plotted in code units and time in units of the initial orbital period. The distance of the MBH from the center changes from 1.4 at t=0 to 0.6 at t=1.5 $`\mathrm{t}_{\mathrm{orbital}}`$ in the gaseous case. However, in the stellar case the average distance is not 0.6 until t=2.3 $`\mathrm{t}_{\mathrm{orbital}}`$. This suggests that the drag timescale is approximately a factor of 1.5 shorter for the gas configuration. The velocity of the BH remains supersonic through most of the run, thus we are measuring $`f^{(\mathrm{star},\mathrm{gas})}`$ close to its peak value (see Fig 4b).

An attempt to explain the behavior of the MBH (for both mediums, in terms of Eqs 8-10) is shown in Fig. 4a. The dashed and continuous black lines represent the radial distance of the MBH in the stellar and gaseous medium. The green line is the fit of the Chandrasekhar formula (Eq 8) to the evolution of the MBH in the stellar case, with $`\mathrm{ln}(\mathrm{b}_{\mathrm{max}}/\mathrm{b}_{\mathrm{min}})`$=$`\mathrm{ln}\mathrm{\Lambda }`$=3.1 as the best fit value for the coulomb logarithm. This value is consistent with results found by other authors for orbiting objects with masses of the order of 1% of the total mass in stars (Lin & Tremaine 1983, Cora et al. 1997).

We try to reproduce the gaseous simulation using Ostriker’s formula (Eqs 9 & 10) for $`\mathrm{ln}\mathrm{\Lambda }`$=3.1, and the result is shown by the blue line in Fig. 4a. It clearly overestimates our prediction for the gaseous drag by a factor of 1.5 (or 2.3 relative to the stellar case). This result is not surprising, because Sanchez-Salcedo & Brandenburg (2001) had already found that Ostriker’s formula overestimates in predicting the gaseous gravitational drag for a body moving at transonic velocities ($`v_{\mathrm{BH}}=\sqrt{2}C_\mathrm{S}`$ at the start of our run). Fig. 4b shows the dimensionless factor of the dynamical friction force as a function of $`v_{\mathrm{BH}}/C_\mathrm{S}`$, as used in different fits in Fig. 4a. In particular Ostriker’s formula is represented by the blue line in Fig. 4b.

Due to the failure of Ostriker’s formula, we consider a parametric function for $`f^{(\mathrm{gas})}()`$ as a better approximation. This function has the same form of the Eq 8 but with two choices for $`\mathrm{ln}\mathrm{\Lambda }`$, with values of 4.7 for $``$ $``$ 0.8 and 1.5 $``$ $`<`$ 0.8. The red curve in Fig. 4b is a plot of this function and is qualitatively like Ostriker’s formula but less peaked. The result obtained with this parametric function is considerably better, showing good agreement with our SPH run as seen in Fig. 4a.

## 4 Evolution of a binary MBH

In §3 we looked at the effect of the gravitational drag on the evolution of a single MBH in a gaseous sphere. In this section, we study the effect of gravitational drag on the evolution of a binary MBH in a gaseous sphere, with the aim of determining the timescales for hardening of the binary and possible merging. To investigate this problem, we follow the evolution of two massive bodies of equal mass in the same stable gas configuration as in §3. As in the single MBH case, the gas mass is much greater than the mass in the binary and the gas mass within the binary produces the acceleration necessary to maintain the circular orbit. As long as the condition $`M_{\mathrm{gas}}(r<r_{\mathrm{binary}})M_{\mathrm{binary}}`$ is satisfied, the evolution of the binary will be determined by the gravitational interaction of each black hole with its surrounding medium, as in the problem studied in $`\mathrm{\S }`$ 3.3. We use a resolution of 100,000 SPH particles to model the gas sphere and run the code for various MBH masses.

Fig. 5 shows the response of the gaseous medium to the presence of the MBH binary. The colors in the figure represent the density enhancement $`(\rho (r,t)/\rho (r,0))`$ in the plane of the orbit of the binary. A couple of spiral shocks propagate outwards, one caused by each MBH. The black curves show the individual orbits of the MBHs between t=0 and t=1. This figure shows the simulation for the case of a MBH with 1% of the total mass of the gas sphere. These density enhancements appear similar to that formed by a single MBH with the same mass (see Fig. 2).

In our study we consider three cases with different black hole masses of 1%, 3% and 5% of the total gas mass. Figure 6 shows the evolution of the binary separation in these three cases over several orbits. The different curves represent the binary separation for different binary to gas mass ratios: $`M_{\mathrm{binary}}=0.02M_{\mathrm{gas}}`$ (red), $`M_{\mathrm{binary}}=0.06M_{\mathrm{gas}}`$ (blue), $`M_{\mathrm{binary}}=0.1M_{\mathrm{gas}}`$ (green). The binary evolution has a strong dependence on the MBH masses, the more massive cases suffering a stronger deceleration. The time required for the binary separation to decrease by a factor of 2 decreases from $`\mathrm{t}1.25`$ for $`M_{\mathrm{binary}}=0.02M_{\mathrm{gas}}`$ to $`\mathrm{t}0.5`$ for $`M_{\mathrm{binary}}=0.1M_{\mathrm{gas}}`$ (times given in units of the initial orbital period).

Fig. 6 shows clearly the fast decay of the binary separation. For all mass ratios, the binary MBH reaches the resolution limit, the gravitational softening length, in only 2 to 4 initial orbital periods, depending on the binary mass. These timescales are shorter than in the stellar case because, as seen in §3.3, the timescales are a factor 1.5 faster due to the collisional nature of the background medium.

At the final separation, gravitational radiation is still an inefficient mechanism for angular momentum loss, but these results are limited by the gravitational softening length. Simulation of later stages of evolution therefore requires higher numerical resolution, as is described below.

### 4.1 Further Evolution Using Higher Resolution Simulation

To continue the evolution of the binary it is necessary to reduce the softening length, something that is extremely expensive computationally. We choose to reduce it from $`ϵ_{\mathrm{soft}}`$=0.01 to 0.001 in order to follow the evolution of the binary separation by one more order of magnitude. We also need to increase the resolution to resolve the region inside r = 0.1 where most of the dynamics occurs. We use the particle splitting procedure described by Bromm (2000), Kitsionas (2000), and in Appendix B to achieve this goal. The procedure is applied to the case where each MBH has 1% of the total gas mass in gas. We split the SPH particles at a time of 3 initial orbital periods when the binary enters the r $`<`$ 0.1 region; for each parent particle we introduce $`N_{\mathrm{split}}`$=20 child particles and therefore we increase the number of particles to 2,000,000.

Fig. 7 shows the evolution of the binary, from the beginning of the low-resolution simulation to the end of the high-resolution simulation. The red curve is the same calculation shown in Fig. 6 but on a logarithmic scale, and the green curve shows the MBH’s separation in the high-resolution calculation. The green curve shows an almost linear trend in logarithmic separation that differs considerably from the behavior shown by the low-resolution simulation. At the end of this high resolution simulation, the binary separation has decreased by almost three orders of magnitude.

The overall evolution of the binary MBH can be divided into two physical regimes: i)Up to t=3 $`\mathrm{t}_{\mathrm{orbital}}`$ where $`M_{\mathrm{bin}}M_{\mathrm{gas}}(r<r_{\mathrm{bin}})`$, the angular momentum loss of the binary is governed by the gravitational drag exerted by the gaseous medium. ii) After t=4 $`\mathrm{t}_{\mathrm{orbital}}`$ the binary satisfies the condition $`(M_{\mathrm{bin}}>>M_{\mathrm{gas}}(r<r_{\mathrm{bin}}))`$, therefore the binary completely dominates the gravitational potential in its vicinity and the response of the medium to that gravitational field is an ellipsoidal density distribution as seen in Fig. 8. The axis of the ellipsoid is not coincident with the binary axis but lags behind it, and this offset produces a torque on the binary that is now responsible for the angular momentum loss. The evolution of the system is approximately self similar because the size of the ellipsoid depends on the binary separation and this is responsible for the nearly exponential decay after t=4 $`\mathrm{t}_{\mathrm{orbital}}`$ (Fig. 7). In the next section we study this torque in more detail. From t=3 $`\mathrm{t}_{\mathrm{orbital}}`$ to 4 $`\mathrm{t}_{\mathrm{orbital}}`$ is a transition period between these two physical regimes.

### 4.2 Ellipsoidal Torque

Our numerical results for the final stage of evolution can be understood in terms of a simple analytical model, consisting of a binary system embedded in a uniform ellipsoidal density distribution. The major axis of the ellipsoid lags behind that of the binary by a constant angle, while the axis of rotation of the ellipsoid coincides with that of the binary. The gravitational potential of an uniform ellipsoid is given by:

$$\mathrm{\Phi }(x,y,z)=\pi G\rho (\alpha _0x^2+\beta _0y^2+\gamma _0z^2+\chi _0),$$
(11)
where $`\alpha _0`$, $`\beta _0`$, $`\gamma _0`$, $`\chi _0`$ are constants given in Lamb (1879) and depend on the ratio between the principal axes of the ellipsoid.

The treatment of the problem is an extension of that given in Goldstein (1950) for the two body problem and is analogous to the development by Boss (1984) in the context of star formation. The binary system can be treated as an equivalent one body problem, subject to an external gravitational potential (Eq 11). The Lagrangian of the system, in cylindrical coordinates, has the form

$$L=\frac{1}{2}\mu \left[\left(\frac{dr}{dt}\right)^2+r^2\left(\frac{d\varphi }{dt}\right)^2+\left(\frac{dz}{dt}\right)^2\right]+\frac{Gm^2}{r}\frac{\pi }{2}G\rho m(\alpha _0r^2\mathrm{cos}^2(\mathrm{\Delta }\varphi )+\beta _0r^2\mathrm{sin}^2(\mathrm{\Delta }\varphi )+\gamma _0z^2+\chi _0),$$
(12)
where $`\mathrm{\Delta }\varphi `$ is the angle between the major axis of the ellipsoid and the axis between the binary members. As expected from the self similar behavior, in our simulation (§4.1) we found that $`\mathrm{\Delta }\varphi `$ is approximately constant. In this simple model we assume that the angle $`\mathrm{\Delta }\varphi `$ is constant and this produces an explicit dependence of the Lagrangian on the coordinate $`\varphi `$, therefore its canonical momentum is not conserved. The Euler-Lagrange equation for the coordinate $`\varphi `$ is

$$\frac{d}{dt}\left(r^2\frac{d\varphi }{dt}\right)=2\pi G(\beta _0\alpha _0)\mathrm{cos}(\mathrm{\Delta }\varphi )\mathrm{sin}(\mathrm{\Delta }\varphi )\rho r^2.$$
(13)
Throughout our simulation, the mean density of the ellipsoid in the vicinity of the binary continues to increase as the binary separation decreases, owing to the gravitational attraction of the black holes. On the right side of equation $`(\text{13})`$ the term that evolves with time is $`\rho r^2`$, which in our simulation varies roughly as $`\mathrm{t}^2`$, as is seen in Fig. 9. The product $`\rho r^2`$ decreases only weakly with time because $`\rho `$ increases strongly as r decreases, increasing almost as $`r^2`$. Also in this later stage, the motions of the MBHs are almost circular, i.e.: $`v=\left((\frac{dr}{dt})^2+(r\frac{d\varphi }{dt})^2\right)^{1/2}r\frac{d\varphi }{dt}`$ and the circular velocity is $`v=\left(\frac{GM_{12}}{r}\right)^{1/2}`$. Under those conditions Eq $`(\text{13})`$ has the following solution

$$r(t)=\left[r_0^{1/2}+\mathrm{\hspace{0.17em}\hspace{0.17em}2}\pi (\beta _0\alpha _0)\mathrm{cos}(\mathrm{\Delta }\varphi )\mathrm{sin}(\mathrm{\Delta }\varphi )\rho _0r_0^2t_0^2\sqrt{\frac{G}{M_{12}}}\left(\frac{1}{t}\frac{1}{t_0}\right)\right]^2,$$
(14)
where $`r_0`$ & $`\rho _0`$ are respectively the distance and mean density at the time $`\mathrm{t}_0=4t_{\mathrm{orbital}}`$ and $`M_{12}=M_{\mathrm{BH1}}+M_{\mathrm{BH2}}`$. We apply this model to a spheroid with principal axes in the ratios 2:1 (i.e. a=2b=2c). In such a configuration $`\beta _0`$=0.825 and $`\alpha _0`$=0.35. The mean angle between the binary and the major axis $`(\mathrm{\Delta }\varphi )`$ is $`22.5^o`$ and the initial mean density of the ellipsoid $`(\rho _0)`$ is 10. The black crosses in Fig. 10 represents the separation given by Eq $`(\text{14})`$, and by comparing with the SPH run (green line), we see that this simple model describes successfully the evolution of the binary MBH with 10% accuracy. Thus, we can interpret the SPH result shown in Fig. $`(\text{7})`$ as the angular momentum loss due to the torque produced by the ellipsoid.

## 5 Application to a Galaxy Center

We are interested in studying the evolution of a binary MBH after a recent merger of gas-rich galaxies. During the merger gravitational torques remove angular momentum from shocked galactic gas causing most of the gas to fall into the center (Barnes & Hernquist 1991, 1996). Those nuclear inflows fuel starburst activity manifested by anomalous colors (Larson & Tinsley 1978) and enormous infrared luminosities (Sanders & Mirabel 1996) characteristic of merging galaxies. Numerical simulations show that typically 60% of the gas from the progenitor galaxies ends up in a nuclear cloud with dimensions of several hundred pc (Barnes & Hernquist 1991, 1996). On these scales, the galactic gas dominates the gravitational potential. The best examples of such configurations are the Ultra Luminous Infrared Galaxies (ULIRG), which typically have $`10^{10}M_{}`$ of atomic and molecular gas in their inner regions (r $``$ few 100pc). Our model can be scaled to represent the conditions in the inner regions of an ULIRG, the time conversion of the initial orbital period into physical units being

$$\mathrm{t}_{\mathrm{orbital}}=\left(\frac{L}{100\mathrm{p}\mathrm{c}}\right)^{3/2}\left(\frac{2.510^{10}\mathrm{M}_{}}{M_{\mathrm{gas}}}\right)^{1/2}\mathrm{\hspace{0.17em}\hspace{0.17em}7.4}10^5\mathrm{yr}$$
(15)
Thus, if we convert the simulation with $`M_{\mathrm{BH}}`$=0.01 $`M_{\mathrm{gas}}`$ into physical units, a binary with $`M_{\mathrm{BH}}=10^8M_{}`$ will reduce its separation from 200 pc to 0.1 pc in only 5.5 $`10^6`$ yrs.

We don’t find any sign of ejection of the surrounding gas, as happens with the stars in the later stages in the evolution of a binary MBH in a stellar system (Begelman, Blandford & Rees 1980; Makino & Ebisuzaki 1996; Milosavljevic & Merritt 2001,2003). Therefore we expect that the torque mechanism that we have described will continue to operate until additional physical effects become relevant (e.g. radiation pressure). In any case, at the end of our simulation, the black hole binary is entering the regime where gravitational radiation is efficient (see §5.1), driving the final coalescence of the MBH. Therefore it is reasonable to expect that the binary MBH will coalesce in a merger of gas-rich galaxies.

### 5.1 Gravitational Radiation

The final phase in the evolution of a binary MBH occurs when eventually the MBHs become close enough to allow gravitational radiation to become an efficient mechanism for angular momentum loss.

The gravitational radiation timescale (Peters 1964) is

$$t_{\mathrm{gr}}=\left(\frac{\dot{a}}{a}\right)^1=\frac{5}{64}\frac{c^5a^4}{G^3M_{12}^3}F(e)$$
(16)
where $`M_{12}=M_{\mathrm{BH1}}+M_{\mathrm{BH2}}`$, $`a`$ is the binary separation and $`F(e)`$ is:

$$F(e)=\frac{(1e^2)^{7/2}}{1+\frac{73}{24}e^2+\frac{37}{96}e^4}$$
(17)
$`F(e)`$ contains the eccentricity dependence which is a weak for small e; e.g. $`F`$(0)=1 and $`F`$(0.5)$``$0.205.

At the end of our simulations we reach a final separation of 0.1 pc. At this distance the gravitational radiation timescale is $`\mathrm{t}_{\mathrm{gr}}1.210^8`$yr. We feel confident that we can extrapolate our results to smaller separations, because the ellipsoidal torque model describes successfully the evolution of the binary MBH. If we extrapolate our results to a separation of 0.05pc, we predict that the binary MBH will merge in only $`10^7`$yr. This timescale is smaller than the time that star formation takes to consume most of the gas in such a dense cloud (a few $`10^8`$ yr), justifiying our assumption that the binary interacts with a gaseous background instead of with a stellar one.

### 5.2 Another Application: Formation of Close Binary Systems

An important problem in star formation theory is that of the formation of close binary stars. Orbital decay of forming binary systems due to gravitational drag may help to explain the distribution of binary separations, since all of the proposed formation mechanisms (fission, fragmentation, capture, disintegration of a larger stellar system) predict an insufficient number of close binary systems (see Bonnell 2001 for a review). Bate, Bonnell & Bromm (2002) show that gravitational interaction of a binary with its surrounding medium is a possible mechanism to form close binary systems. Our work has enough resolution to study such systems in more detail than Bate, Bonnell & Bromm (2002).

We scale our model to represent the conditions in star forming regions. Thus, the time conversion of the initial orbital period into physical units is

$$\mathrm{t}_{\mathrm{orbital}}=\left(\frac{\mathrm{M}_{\mathrm{gas}}}{10\mathrm{M}_{\mathrm{}}}\right)^{1/2}\left(\frac{\mathrm{D}}{0.01\mathrm{pc}}\right)^{3/2}\mathrm{\hspace{0.17em}3.6}\times 10^4\mathrm{yr},$$
(18)
where $`M_{gas}`$ is the total mass of the gas and D is the distance unit.

We apply our result to the case of low mass star formation, assuming a total gas mass of 10$`M_{\mathrm{}}`$ and an initial separation between the protostars of 0.02pc. This gives an initial orbital period of $`3.6\times 10^4`$ yr. Thus we predict that a binary with 0.1$`M_{\mathrm{}}`$ protostars will reduce its separation from 0.02 pc to 2 AU in 7 initial orbital periods (Fig. 7), or only $`2.6\times 10^5`$ years.

Since the gravitational drag timescales are predicted to be comparable with those of the star formation process, gravitational drag can play an important role in the formation of stars and can determine the final separations of binaries. Therefore, the orbital evolution of binary systems due to gravitational drag is a promising mechanism for the formation of close binary stars.

## 6 Conclusions

In this paper we study the role of gas in driving the evolution of a binary MBH, and we follow its evolution through many orbits and close to the point where gravitational radiation becomes important. We present results for a relatively idealized case in which the gas is assumed to be in a nearly spherical and relatively smooth distribution.

There are important differences in the long term evolution of a binary MBH in a gaseous medium, compared to a stellar background. From N-body simulations, it is expected that the merging of an MBH binary eventually stalls in a stellar background. However, our simulations suggest that the MBH binary will finally merge if it’s embedded in a gaseous medium.

In the early evolution of a binary MBH, the separation diminishes due to the gravitational drag exerted by the background gas. This decay is faster than in a stellar background by a factor $``$1.5 for a given density profile. In the later stages, when the binary MBH dominates the gravitational potential in its vicinity, the medium responds forming an ellipsoidal configuration. The axis of the ellipsoid lags behind the binary axis, and this offset produces a torque on the binary that is now responsible for the continuing loss of angular momentum. The evolution of the system is approximately self-similar because the size of the ellipsoid shrinks with the binary separation. This torque is able to reduce the binary separation to distances where gravitational radiation is efficient. We don’t find any sign of ejection of the surrounding gas, as happens with stars in the later stages in the evolution of a binary MBH in a stellar system. Moreover, the gas density $`\rho `$ in the vicinity of the binary increases strongly as the separation r decreases, increasing almost as $`r^2`$.

Assuming typical parameters from observations of ULIRG, the overall binary evolution has a merger timescale of $`10^7`$ yrs. This timescale is fast enough to support the assumption that the binary is in a mostly gaseous background medium, because the timescale for gas depletion in a starburst region is typically several times $`10^7`$ yrs. Galaxies typically merge in $`10^8`$ yrs, and therefore these results imply that in a merger of gas-rich galaxies the BHs will coalesce soon after the galaxies merge. The final coalescence has crucial implications for possible scenarios of massive black hole evolution and growth. In particular, this result supports scenarios where hierarchical build-up of massive black holes plays an important role. The final coalescence of the black holes lead to gravitational radiation emission that would be detectable up to high redshift by LISA.

In subsequent papers, we will present results of calculations in progress for cases in which the gas is cooler, and has a more flattened and clumpy spatial distribution. Preliminary results show again, as in the case of smooth gas, a large decrease of the binary separation in a few initial orbital periods, and thus support the more general applicability of the results presented here for an idealized case.

## 7 Acknowledgments

A. E. thanks Fundación Andes for fellowship support. We are grateful to Volker Springel for making available to us a version of GADGET and we thank Priya Natarajan for valuable comments. Some of the simulations were performed at the Universidad de Chile computer cluster, funded by FONDAP grant 11980002.

## Appendix A Code Testing

In order to validate the proper use of GADGET, we have done additional tests to those performed by the code authors (Springel, Yoshida & White 2001).

### A.1 Effects of Numerical Diffusion

One of the disadvantages of using multidimensional Lagrangian schemes is that they suffer from numerical diffusion, introduced by the appearance of vortices that requires remapping the tangled and twisted code grid. The SPH technique can be viewed as an intrinsic remapped Lagrangian grid which implies a sizable numerical diffusion.

To test the effects of numerical diffusion in the GADGET code, we simulate the adiabatic oscillation of a polytrope following Steinmetz & Muller (1993) and investigate the damping due to numerical diffusion. We first need to establish the hydrostatic equilibrium configuration, which is also a test of the interaction between the gravitational and the hydrodynamic elements of the code. We chose a density profile given by the solution to the Lane-Emden equation for n=3/2, and let it relax under the polytropic equation of state

$$p=K\rho ^{5/3}.$$
(A1)
Fig. 11 compares the resulting density distribution from the SPH code using 5000 particles with the Lane-Emden solution (Kippenham and Weigert 1991). The dots follow the curve throughout the cloud radius to within a 5% of spread in density.

In order to test the damping of an adiabatic oscillation due to numerical diffusion, we imposed (on the equilibrium model) a spherically homologous contraction by a factor 0.8 in radius in order to produce strong pulsations of the polytropic sphere. Fig. 12 shows the evolution of the central density as a function of time, and the oscillation frequency has the predicted value $`\nu _{\mathrm{osc}}1/(4t_{\mathrm{sound}})`$.

The pulsation amplitude decreases in time due to the numerical diffusion of the code. After 5 pulsation periods it has decreased by a factor 2.5. This diffusion depends on the artificial viscosity used, in this case (shear-reduced artificial viscosity) the results are similar to that found by Bate (1995). Other prescriptions for artificial viscosity produce less damping in the oscillation, but those are not suitable for our problem because they produce a higher loss of angular momentum due to the extra shear viscosity.

### A.2 Shock Modeling

The supersonic motion of MBHs in a gaseous medium will produce strong shocks, and therefore it is crucial to know if the code is able to resolve shocks. We test the ability of GADGET to handle shocks by following the adiabatic collapse of an initially isothermal gas sphere, which is a common test case for SPH codes (Evrard 1988, Hernquist & Katz 1989, Steinmetz & Muller 1993, Thacker et al. 2000). The initial (non equilibrium) state consists of an adiabatic gas sphere with total mass M, radius R and density profile

$$\rho (r)=\frac{M}{2\pi R^2}\frac{1}{r}$$
(A2)
The gas is initially isothermal, with a specific internal energy of

$$u=0.05\frac{GM}{R}$$
(A3)
The thermal energy is not enough to provide pressure support, therefore the gas sphere begins to collapse. As the collapse proceeds the gas is heated until the core temperature rises enough to produce a central bounce and a shock wave is formed and propagates outwards (Evrard 1988). After the shock has passed through most of the gas, the gas sphere should reach virial equilibrium. The evolution of the kinetic, thermal, gravitational and total energies have been widely discussed in the literature (Evrard 1988, Hernquist & Katz 1989, Steinmetz & Muller 1993, Thacker et al. 2000) and the results for the GADGET code are published in Springel, Yoshida & White (2001). We focus our attention on the outward propagating shock and test how well is handled.

Fig. 13 shows the density of the collapsing system obtained for a simulation with 30,000 particles at t=0.8. Initially the gas sphere is at rest. We obtain the coordinates from a distorted regular grid that reproduces the density profile and we give our results in units of M = R = G = 1. The shock propagating outward is clearly visible between r=0.15 and r=0.25 in the figure. The results agree very well with the results obtained by a finite difference method (Thomas 1987; red solid line). The results from Fig. 13 are also in very good agreement with those found by Hernquist & Katz (1989), Steinmetz & Muller (1993) and even with the results obtained by the best artificial viscosity implementation of Thacker et al. (2000).

### A.3 Collapse of rotating sphere

The collapse of a rotating sphere has become a standard test for galaxy formation codes (Navarro & White 1993). One relevant application of our work is to study the formation of MBH in protogalaxies. The spherical ’protogalaxy’ is composed of dark matter providing 90% of the mass and baryons providing the remaining 10%. Initially the cloud has a density profile of $`\rho (r)r^1`$ . The velocities are chosen so that the sphere is in solid rotation with a spin parameter

$$\lambda =\frac{LE^{1/2}}{GM^{5/2}}\mathrm{\hspace{0.17em}0.1}$$
(A4)
the initial radius of the cloud is $`100\mathrm{Kpc}`$, and the total mass (gas and dark matter) is $`10^{12}\mathrm{M}_{}.`$ The dynamical time-scale of the system is then $`\tau _{\mathrm{dyn}}3.7\times 10^8`$yr. The gas is initially at a uniform temperature of $`T=10^3K`$ that is well below the equilibrium virial temperature of the system ($`T_{\mathrm{virial}}10^6`$K). We use the same units as Navarro & White (1993): G=1, $`[\mathrm{mass}]=10^{10}\mathrm{M}_{}`$, $`[\mathrm{time}]=\mathrm{4.71\hspace{0.17em}10}^6\mathrm{yr}`$, $`[\mathrm{distance}]=1\mathrm{Kpc}.`$

As in the case of Evrard collapse (A.2), the gas evolution is dominated by a shock wave that propagates outward and transforms most of the kinetic energy into heat. Fig. 14(a,b) shows the evolution of the gas and the dark matter components, at three given times: t=32, 160, 320. The dark matter is virialized soon after collapse forming a tight core (Fig. 14a). The gas shows a distribution similar to the dominant dark matter except for being slightly more spherical. The results are consistent with those of Navarro & White (1993).

We repeat the simulation including radiative cooling; we use a cooling function for primordial gas interpolated from Sutherland & Dopita (1993). When radiative cooling is included the evolution of the system changes considerably because the cooling timescale in the denser regions becomes much smaller than the dynamical timescale of the gas; This produces radiation losses that dissipate the thermal energy gained through shocks. With no pressure support, the gas collapses until it is centrifugally supported forming a flat disk structure, as is clearly seen in Fig. 14c.

## Appendix B Particle Splitting

Particle splitting (Bromm 2000, Kitsionas 2000, Bromm & Loeb 2003) is a procedure designed to increase the resolution of a SPH code in some regions where it is required. In this procedure, all particles in the region of interest are replaced by set of smaller particles, increasing the resolution of the simulations. We refer to the smaller particles as ’fine’ particles and the original ones as ’coarse’ particles.

In our simulations we need to increase the resolution in the later stages, to be able to resolve the region inside r = 0.1, where most of the relevant dynamics occur. We stopped the simulation just before the MBHs arrive to the inner region (r $`<`$ 0.1). We then resampled the fluid by splitting each SPH particle from the original simulation (called ’parent’ particle and denoted by p) into $`N_{\mathrm{split}}`$ particles (called ’child’ particles and denoted by c).

The child particles are randomly distributed according to the SPH smoothing kernel W$`(r_\mathrm{c}r_\mathrm{p};h_\mathrm{p})`$ where $`h_\mathrm{p}`$ is the smoothing length of the parent particle (see Springel, Yoshida & White 2001 for the kernel implementation used by GADGET). The velocity is directly inherited from the parent particle, i.e. $`\stackrel{}{v_\mathrm{c}}=\stackrel{}{v_\mathrm{p}}`$ and each child particle is assigned a mass of $`m_\mathrm{c}`$ = $`m_\mathrm{p}/N_{\mathrm{split}}.`$

We apply the resampling procedure to the whole sphere instead of only the inner region. The reason is that in the boundary between the coarse and fine regions the density is overestimated for the coarse particles and underestimated for the fine ones. The net result is an inter-penetration of the coarse particles into the r $`<`$ 0.1 region, when they become close enough to the MBH to produce a numerical overestimate of the drag force. Kitsionas (2000) tried to solve this artifact by changing the SPH formalism from a fixed total number of neighbors to a fixed total mass of the neighbors. This reduced the overestimate of the fluid properties for the coarse particles but was not able to prevent the penetration of coarse particles into the fine region. For that reason we choose to resample the whole fluid although it is numerically more expensive.