text:
# Taming the Rugged Landscape: Techniques for the Production, Reordering, and Stabilization of Selected Cluster Inherent Structures

## I Introduction

The minimization/optimization problem is one of the more ubiquitous and challenging in computational science Press et al. (1992). Central to researchers in the physical sciences and engineering, this problem is also of primary importance to social, biological, and economics investigators.

Driven in large measure by such widespread interest, there has been appreciable progress on the minimization problem. Especially notable have been algorithmic advances in the form of annealing and stochastic relaxation approaches Press et al. (1992); Kirkpatrick et al. (1983); Purismia and Scheraga (1986); Li and Scheraga (1987); Faken et al. (1999); Lee and Berne (2000) as well as basin-hopping techniques Wales and Doye (1997); Wales and Scheraga (1999). In both classical and quantum form, these methods offer valuable, complementary alternatives to traditional, gradient or pseudo-gradient approaches Press et al. (1992).

In addition to algorithmic developments relevant to the minimization problem, there have also been notable advances in the tools to classify and analyze the topography of the underlying objective functions. In chemical applications, the principal focus of the remainder of our discussion, the objective function of interest is typically a specified potential or free energy surface. Following Stillinger and Weber Stillinger and Weber (1983, 1984), it is useful to perform an “inherent structure” decomposition of the associated configuration space by employing the minima (local and global) of this surface. These inherent structures, their relative orderings, and their connectivity provide important information concerning the structure, function, and dynamics of the associated physical system. Disconnectivity analysis introduced by Czerminski and Elber Czerminski and Elber (1990), discussed by Becker and Karplus Becker and Karplus (1997) and developed by, among others, Wales, Doye and Miller Miller et al. (1999); Doye et al. (1999a, b) has proved especially valuable with respect to these latter tasks.

As evidenced by the development of classical and quantum annealing methods, there is an important interplay between minimization and the Monte Carlo sampling problem. Both applications, for example, are concerned with overcoming barriers that inhibit the interconversion or isomerization of the various inherent structures of the problem. Consequently, developments in one field contain implications for developments in the other. Advances in rare event sampling methods, such as J-walking Frantz et al. (1992); Tsai and Jordan (1993a) and parallel tempering Marinari and Parisi (1992); Geyer and Thompson (1995); Neirotti et al. (2000); Calvo et al. (2000) methods, thus contain implications for the minimization problem.

The field of atomic and molecular clusters has been and continues to be an important test bed for the development and application of minimization and analysis methods. Utilizing the methods outlined above, researchers have produced a relatively coherent picture of the relationship between the nature of underlying PES and the physical properties of the associated systems. For example, from the studies by Berry et al. Jellinek et al. (1986); Wales and Berry (1990); Berry (1993, 1994), the single component studies of Wales et al. Miller et al. (1999); Doye et al. (1999a, b), the mixed LJ cluster studies of Jordan et al. Munro et al. (2002) as well as the research of others Tsai and Abraham (1978); Frantz (1996, 1997); Serra and Kais (1997), we have begun to understand the nature of systems for which the lowest inherent structures can or cannot be readily located.

In present paper we would like to build upon advances in the minimization problem by effectively turning the logic “upside down”. That is, instead of asking what we have to do in order to locate or sample the global minimum of a specified potential energy surface, we wish instead to ask how we might go about controllably inducing new structures as well as reordering and stabilizing existing structures that are characteristic of higher-lying local minima. Basically, we seek to utilize what we have learned about what it takes to avoid local minima to instead controllably produce them.

In principle, one can envision efforts involving both thermodynamic and kinetic approaches. In the present work we shall focus principally on the thermodynamic issues. Furthermore, we shall limit the discussion in the present work to applications involving clusters. As discussed elsewhere JEL (1999), clusters are of appreciable technological importance, are valuable as prototypes for the study of the properties of extended systems, illuminate issues related to the size-dependence of selected physical properties, and provide valuable test beds for the development and application of emerging computational techniques. This combination of formal, computational and technological interest has produced a vast and growing cluster literature Freeman and Doll (1996); JEL (1999); Knickelbein (1999).

The remainder of the paper is organized as follows: In Section II we outline the computational details of the present study. We discuss the methods we use to determine the inherent structures and transition states of a specified cluster’s potential energy surface. Using these methods, we examine specific results for two prototype systems in Section III. These particular results are designed to demonstrate “proof of principle” with respect to the basic objectives of the present study for selected systems. Finally, in Section IV we summarize our results and speculate about likely future research directions.

## II Computational Details

The present Section describes the computational details of our investigations involving binary clusters of the form $`\mathrm{X}_\mathrm{n}\mathrm{Y}_\mathrm{m}`$. Our overall interest will be to explore the extent to which we can utilize the “adatoms” (i.e. the Y-system) to induce, reorder and stabilize selected inherent structures in the “core” X-system. While one can easily imagine applications involving both more and more complex components, we feel these relatively simple, two-component clusters are a convenient starting point for an initial study of the issues we raise.

We shall assume in what follows that the total potential energy is composed of a pairwise sum of Lennard-Jones interactions. Specifically, we assume that the total potential energy, $`\mathrm{V}_{tot}`$, for an N-particle system is given by

$$V_{tot}=\underset{i<j}{\overset{N}{}}v_{ij}(r_{ij}),$$
(1)
where the pair interaction as a function of the distance between particles i and j, $`r_{ij}`$, is given by

$$v_{ij}(r_{ij})=4ϵ_{ij}[(\frac{\sigma _{ij}}{r_{ij}})^{12}(\frac{\sigma _{ij}}{r_{ij}})^6].$$
(2)
In Eq. (2) the constants $`ϵ_{ij}`$ and $`\sigma _{ij}`$ are the energy and length-scale parameters for the interaction of particles i and j.

For a two-component system, we must specify both the “like” (X-X, Y-Y) as well as the “mixed” (X-Y) interactions. With an eye toward studying trends in the results as opposed to results for particular physical systems, it is convenient to reduce the number of free parameters. To do so, we shall assume in the present study that the “mixed” Lennard-Jones values are obtained from the “like” Lennard-Jones parameters via usual combination rules Hirschfelder et al. (1954)

$$\sigma _{_{XY}}=\frac{1}{2}(\sigma _{_{XX}}+\sigma _{_{YY}})$$
(3)
$$ϵ_{_{XY}}=\sqrt{ϵ_{_{XX}}ϵ_{_{YY}}}.$$
(4)
Furthermore, we note that with the mixed Lennard-Jones parameters specified as in Eqs.(3) and (4), the inherent structure topography of the “reduced” potential energy surface of the binary system (i.e. $`V_{tot}/ϵ_{_{XX}}`$) is a function of only two parameters, ($`\sigma ,ϵ`$), the ratios of the corresponding adatom/core length and energy parameters

$$\sigma =\sigma _{_{YY}}/\sigma _{_{XX}}$$
(5)
$$ϵ=ϵ_{_{YY}}/ϵ_{_{XX}}.$$
(6)
If necessary for a discussion of a specific physical system, the absolute bond lengths, energies, activation energies, etc. can be obtained from the corresponding “reduced” results by a simple rescaling with the appropriate core-system Lennard-Jones parameters.

The computational task in our study is thus one of exploring and characterizing the (reduced) potential energy surface of our binary cluster systems as a function of the number of (core, adatom) particles, (n,m), and for given ($`\sigma ,ϵ`$) ratios. In typical applications the lowest $`\mathrm{N}_{IS}`$ inherent structures and the associated disconnectivity graphs are determined. For the applications reported here, $`\mathrm{N}_{IS}`$ is generally of the order of a few hundred (thousand) or less. Depending on the size of the cluster, inherent structures are found either via conjugate gradient methods starting from randomly chosen initial configurations, or by more systematic surface exploration methods such as those outlined by Wales and co-workers Doye et al. (1999a) and by Jordan et al. Tsai and Jordan (1993b). In all cases, the inherent structures that are located are confirmed to be stable minima via a standard Hessian analysis. To reduce the chance we miss particular local or global minima, we monitor the number of times individual inherent structures are found and demand that each of the NIS inherent structures be located a minimum number of times (at least 10) before we terminate our search. Once we are satisfied we have located the relevant inherent structures, transitions states linking these stable minima are obtained using the eigenvector following methods outlined by Miller and Cerjan Cerjan and Miller (1981) and further developed by Simons et al. Simons et al. (1983); Banerjee et al. (1985); Nichlos et al. (1990), Jordan et al. Tsai and Jordan (1993b) and Wales Wales (1994). Finally, with the requisite inherent structures and barriers in hand, we perform a disconnectivity analysis using methods outlined by Czerminski and Elber Czerminski and Elber (1990), Becker and Karplus Becker and Karplus (1997) and Wales et al. Miller et al. (1999) .

## III Numerical Results

In the present Section, we wish to illustrate the general themes we introduced in Section I . We do so by demonstrating that we can accomplish three basic objectives. Specifically, we show that by adding selected “impurity” atoms to bare “core” systems, we can:

1. induce new “core structures”

2. reorder the energies of existing core inherent structures, and

3. stabilize selected inherent structures by controlling the activation energies that determine their isomerization kinetics.

For purposes of illustration, we shall examine numerical results for a few, simple Lennard-Jones systems involving five and seven core atoms, systems well-known from previous studies to have one and four energetically distinct inherent structures, respectively. The inherent structures and their associated energies for these core systems are illustrated in Fig. 1 and 2.

We first consider mixed clusters of the generic type X5Y2. Here two impurity Y-atoms are added to the parent, five-atom X-core. We have chosen this system because it builds upon the very simple five-atom core, a system that has only a single inherent structure, and because the total system has a total of seven atoms, a magic number for icosahedral growth in homogeneous systems. Using the techniques of Section II, we then determine the lowest several inherent structures for a range of ($`\sigma `$,$`ϵ`$) \[c.f. Eq.(5) and Eq.(6)\]. As can be seen from Fig. 3, the total potential energy \[Eq.(1)\] of the lowest inherent structure for the X5Y2 system shows no appreciable structure as a function of the ($`\sigma `$,$`ϵ`$) parameters.

On the other hand, we see in Fig. 4 that the core potential energy, defined as the potential energy of interaction for only the core X-atoms, of the minimum (total) energy cluster clearly breaks into extended regions, each corresponding to a well-defined core structure. The reader should notice that each region in Fig. 4 contains the same “kind” of core structure but their core energies are slightly different. We have chosen a single “average” core energy value to represent all energies in the corresponding domain for plotting convenience.

The distinct core structures, shown in Fig. 4, have been identified by examining their core energies (Ecore) and their principal moments of inertia. For each structure a triplet of values (Ecore, I2, I3) has been associated, where I2 and I3 are the moments of inertia about the principal axes 2 and 3, respectively. We have defined I2 and I3 in the following way: I2=I$`{}_{2}{}^{}{}_{}{}^{^{}}`$/I$`{}_{1}{}^{}{}_{}{}^{^{}}`$, I3=I$`{}_{3}{}^{}{}_{}{}^{^{}}`$/I$`{}_{1}{}^{}{}_{}{}^{^{}}`$ where I$`{}_{1}{}^{}{}_{}{}^{^{}}`$, I$`{}_{2}{}^{}{}_{}{}^{^{}}`$ and I$`{}_{3}{}^{}{}_{}{}^{^{}}`$ are the principal moments of inertia obtained by diagonalizing the inertia tensor of the system. If the triplet of values has not been sufficient to identify a core structure then we have examined the structure visually.

Selected cluster structures illustrating the core arrangements corresponding to various ($`\sigma `$,$`ϵ`$) values are shown in Fig. 5. We see from Figs. 4 and 5 that the $`\mathrm{X}_5\mathrm{Y}_2`$ cluster exhibits core X-atom structures that include trigonal bipyramidal, planar and square pyramidal core geometries. Of these, only the trigonal bipyramidal form is stable in the parent $`\mathrm{X}_5`$ system. This illustrates that a suitable choice of the ($`\sigma `$,$`ϵ`$) parameters can controllably induce core geometries not present as stable minima in the bare cluster. For example, the square pyramid core structure, seen in Fig. 5.2 as a stable system, corresponds to a transition state in the bare $`\mathrm{X}_5`$ cluster.

Fig. 6 represents the $`\mathrm{X}_5\mathrm{Y}_2`$ cluster at four points in Fig. 4 defined by the ($`\sigma `$,$`ϵ`$) coordinates (0.4,0.5), (0.4,1.0), (0.4,1.5) and (0.4,2.0). Here the pairs of coordinates correspond to (a), (b), (c) and (d) part of Fig. 6, respectively. In other words, we keep value of $`\sigma `$=0.4 fixed, while increasing the value of $`ϵ`$. Each disconnectivity graph shows all inherent structures available to the system for the given ($`\sigma `$,$`ϵ`$) values. The global minimum of each system is labeled by number 1 and contains as a recognizable component the square pyramid core structure (see Fig. 5.2). In Fig. 6.a the square pyramid core structure is connected to two inherent structures, labeled by 2 and 3, by pathways whose energies do not exceed $`13.8`$ (in units of $`ϵ_{_{XX}}`$). Since isomer # 2 contains the same core structure (the square pyramid) as the global minimum the corresponding isomerization thus does not lead to a change in the core structure of the cluster. For present purposes, therefore, the barrier that connects them is not a “relevant” barrier. The relevant barriers are those that connect inherent structures that contain different core structures. The inherent structure # 3 contains as the core structure a (distorted) trigonal bipyramid (see Fig. 5.1). Therefore, the isomerization barrier that connects the inherent structure # 3 with global minimum is the lowest relevant isomerization barrier and its value is $`\mathrm{\Delta }`$E1,3=0.986$`ϵ_{_{XX}}`$. Fig. 6.b and Fig. 6.c show that increasing the value of $`ϵ`$ increases isomerization barriers that connect inherent structure #1 (the square pyramid core structure) with inherent structure #2 (the distorted trigonal bipyramid core structure). Numerically, these barriers are $`\mathrm{\Delta }`$E1,2=1.227$`ϵ_{_{XX}}`$ and $`\mathrm{\Delta }`$E1,2=1.431$`ϵ_{_{XX}}`$, respectively. In Fig. 6.d the square pyramid core structure is connected by an isomerization barrier of $`\mathrm{\Delta }`$E1,2=1.647$`ϵ_{_{XX}}`$ with two (almost degenerate in energy) distorted trigonal bipyramid core structures. As illustrated in Fig. 6 and discussed above, the barriers that determine the isomerization kinetics of these newly induced structures are sensitive to the ($`\sigma `$,$`ϵ`$) values and can thus be at least partially controlled. These two simple results are specific demonstrations of goals 1 and 3 stated above.

As a second illustration, we consider mixed clusters of the type $`\mathrm{X}_7\mathrm{Y}_3`$. This system builds upon a parent, seven-atom, “magic number” system known to exhibit a set of four, energetically distinct inherent structures. The core inherent structures and associated energies for the stable $`\mathrm{X}_7`$ inherent structures are presented in Fig. 2. Figure 7, a $`(\sigma ,ϵ)`$ contour plot of the core-atom potential energies of the lowest total energy $`\mathrm{X}_7\mathrm{Y}_3`$ clusters, again reveals the presence of definite “core-phases”.

As illustrated in Fig. 8, some of these regions correspond to various core structures present in the parent $`\mathrm{X}_7`$ system while others correspond to new structures not seen in the original, single-component cluster. We can see from Figs. 7$``$9

that the impurity Y-atoms provide us with significant control over the relative ordering of the core energies of the parent $`\mathrm{X}_7`$ system. Specifically, by choosing an appropriate range of $`(\sigma ,ϵ)`$ values, we can generate $`\mathrm{X}_7\mathrm{Y}_3`$ clusters in which the lowest (total) energy inherent structure can have core structures that are either pentagonal bipyramid, capped octahedral, or bicapped trigonal bipyramidal in nature. Moreover, since we can manipulate the isomerization barriers in these systems, we can at least partially stabilize clusters that exhibit selected core structures with respect to isomerization. This is illustrated in Fig. 9.

Figures 9.a – 9.d represent the $`\mathrm{X}_7\mathrm{Y}_3`$ cluster at four points in Fig. 7 with $`\mathrm{X}_7\mathrm{Y}_3(\sigma ,ϵ)`$ coordinates (0.4,0.5), (0.4,1.0), (0.4,1.5) and (0.4,2.0), respectively. The number of inherent structures available to the $`\mathrm{X}_7\mathrm{Y}_3`$ cluster varies from more than 800 in Fig. 9.a to 400 in Fig. 9.d. Since we are primarily interested in energetically low-lying inherent structures we show only lowest 70 inherent structures. The global minimum of each system is labeled by number 1 and contains as a recognizable component the core structure shown in Fig. 8.5. We should mention that for a given range of $`\sigma `$ and $`ϵ`$ values, ($`\sigma ,ϵ`$)$``$ \[0.1,2.0\], we have not been able to find a global minimum that would contain as a recognizable component inherent structure # 4 of the parent $`\mathrm{X}_7`$ cluster (see Fig. 2.d). This is the reason why none of the domains in Fig. 7 is labeled by number 4. In Fig. 9.a the global minimum, the core structure # 5 (see Fig. 8.5), is linked to inherent structure # 2 which contains the (distorted) capped octahedron core structure (see Fig. 8.2). The isomerization barrier between them is $`\mathrm{\Delta }`$E1,2=0.494$`ϵ_{_{XX}}`$. Fig. 9.b and Fig. 9.c show that increasing the value of $`ϵ`$ increases isomerization barriers, that connect inherent structure #1 (see core structure in Fig 8.5) with inherent structure #3 and #2 (a distorted capped octahedron core structure), respectively. Numerically, these barriers are $`\mathrm{\Delta }`$E1,3=0.975$`ϵ_{_{XX}}`$ and $`\mathrm{\Delta }`$E1,2=1.136$`ϵ_{_{XX}}`$, respectively. The inherent structure #2 in Fig. 9.b contains the same core structure as the global minimum and, therefore, has not been considered relevant for the isomerization (see explanation above for Fig. 6.a). In Fig. 6.d the core structure # 5 is connected to two, energetically almost degenerate, inherent structures labeled by 2 and 3, by pathways whose energies do not exceed $`40.5`$ (in units of $`ϵ_{_{XX}}`$). Similarly to the case of Fig. 6.a, the isomer # 2 contains the same core structure as the global minimum so the corresponding isomerization does not lead to a change in the core structure of the cluster. The barrier that links them is not a relevant barrier. The inherent structure # 3 contains as the core structure a distorted capped octahedron. Therefore, the isomerization barrier which connects the inherent structure # 3 with global minimum is the lowest relevant isomerization barrier and its value is $`\mathrm{\Delta }`$E1,3=1.304$`ϵ_{_{XX}}`$.

## IV Conclusions

In the present work, we have considered the general task of altering core cluster structures. We are, in effect, attempting to turn the logic of the minimization problem upside down. Rather than seeking the global minimum of complex potential energy surfaces, we are instead attempting to exploit what has been learned about the general minimization problem to controllably alter core cluster structures. Specifically, we are examining the extent to which we can induce new core geometries as well as reorder and stabilize existing, higher-lying, local core structures.

Our approach, in the present discussion, has been thermodynamic in nature. We have utilized selected adatoms to effect our desired core cluster modifications. We have presented results for two simple binary cluster examples, the $`\mathrm{X}_5\mathrm{Y}_2`$ and $`\mathrm{X}_7\mathrm{Y}_3`$ systems, to validate our approach.

We speculate that there are at least two important directions for future theoretical development of the present ideas. One direction will be to explore the use of more complex adsorbates to achieve selected core cluster structures. One could, for example, imagine using “exterior” methods in which encapsulating agents of well-defined geometries were utilized to induce desired core structures. Alternatively, “interior” approaches in which complex objects, perhaps even previously engineered clusters, could be utilized as “seeds” or “templates” to produce a desired structure in the surrounding cluster (either globally or locally). Another important direction will be to explore the extent to which previously engineered cluster structures can be assembled using “cluster assembled materials” methods to produce larger scale, macroscopic structures. If this proves possible, it would seem to offer an important direction in the production of novel materials starting from synthetic precursors whose core structures and properties are highly varied and are under user control.

Acknowledgment

The authors acknowledge support from the National Science Foundation through awards No. CHE-0095053. and CHE-0131114. They would also like to thank Dr. M. Miller for helpful discussions and for his gracious assistance with respect to the preparation of the disconnectivity graphs in the present paper.


format: text

text:
# WMAP data and the curvature of space

(February 28, 2003)

## Abstract

Inter alia, the high precision WMAP data on Cosmic Background Radiation marginally indicate that the universe has positively curved (and hence spherical) spatial sections. In this paper we take this data seriously and consider some of the consequences for the background dynamics. In particular, we show that this implies a limit to the number of $`e`$-foldings that could have taken place in the inflationary epoch; however this limit is consistent with some inflationary models that solve all the usual cosmological problems and are consistent with standard structure formation theory.

The Wilkinson Microwave Anisotropy Probe (WMAP) has recently provided high resolution Cosmic Microwave Background data  of importance for cosmology - indeed they constitute a significant contribution towards the goal of developing precision cosmology. Among the interesting conclusions that have been reached from this data are constraints on present value $`\mathrm{\Omega }_0`$ of the total density parameter of the universe. The WMAP results interestingly indicate that while the universe is close to being flat ($`\mathrm{\Omega }_01`$), a closed universe is marginally preferred: $`\mathrm{\Omega }_0>1`$ . In particular, with a prior on the Hubble constant, one gets that $`\mathrm{\Omega }_0=1.03\pm 0.05`$ at 95% confidence level, while combining WMAP data with type Ia supernovae leads to $`\mathrm{\Omega }_0=1.02\pm 0.04`$ or to $`\mathrm{\Omega }_0=1.02\pm 0.02`$ respectively without and with a prior on the Hubble parameter. The latter may be regarded as the present best estimate of this parameter. Note also that this tendency to point toward a closed universe is not strictly speaking either new or restricted to the WMAP data. For instance with a prior on the nature of the initial conditions, the Hubble parameter and the age of the universe, analysis of the DASI, BOOMERanG, MAXIMA and DMR data  lead to $`\mathrm{\Omega }_0=0.99\pm 0.12`$ but to $`\mathrm{\Omega }_0=1.04\pm 0.05`$ if one takes into account only the DASI, BOOMERanG and CBI data, both at $`1\sigma `$-level. The Archeops balloon experiment  points toward $`\mathrm{\Omega }_0=1.00_{0.02}^{+0.03}`$ with data from the HST and a prior on the Hubble constant but to values in the range $`\mathrm{\Omega }_0=1.16_{0.20}^{+0.24}`$ from combined CBR data alone and to $`\mathrm{\Omega }_0=1.04_{0.04}^{+0.02}`$ from combined CBR and supernova data, all at a 68% confidence level. The improved precision from WMAP is clear.

The confirmation of the existence of Doppler peaks and their respective locations tends to confirm the inflationary paradigm , as does the existence of an almost scale invariant power spectrum. Furthermore a nearly flat universe, argued to be predicted generically by most inflationary models , is usually seen as a further evidence. Since

$$\delta (t)\mathrm{\Omega }(t)1=\frac{K}{a^2H^2}$$
(1)
where $`a(t)`$ is the scale function and $`H(t)=\dot{a}/a`$ the Hubble parameter, the WMAP constraints on the curvature of space imply that the curvature density parameter, $`\mathrm{\Omega }_K(t)K/a^2H^2`$, has been smaller than a few percent from the time the largest observable scale, $`k_{\mathrm{min}}a_0H_0`$, crossed the horizon during inflation up to today. It is then a good approximation to neglect the curvature from the time $`k_{\mathrm{min}}`$ crossed the horizon to today when performing structure formation calculations, mainly because $`\mathrm{\Omega }_{k_{\mathrm{min}}}=\mathrm{\Omega }_{\mathrm{today}}`$ (see Fig. 1). The constraints on the curvature of the universe thus back up the validity of calculations underlying the tests of inflation performed up to now that assume a flat universe and focus on the properties of the perturbations (spectral indices, tensor modes, statistics, etc.).

Nevertheless, when $`K0`$ (which is highly favoured over the case $`K=0`$ by probability considerations: the latter are of measure zero in the space of Friedman-Lemaître-Robertson-Walker universes in most measures), $`\mathrm{\Omega }_K`$ is driven toward 0 exponentially during inflation. The curvature term therefore was necessarily dominant at early enough times, if inflation undergoes a sufficient number of $`e`$-foldings, mainly because the curvature term behaves as $`a^2`$ and will tend to dominate over any kind of matter having an equation of state less stiff than $`1/3`$, and thus over a slow rolling scalar field or a cosmological constant, in the early inflationary era.

We argue that the study of dynamics of the background for such models can give interesting insight into the inflationary era and its dynamics prior to $`k_{\mathrm{min}}`$ horizon exit. Working backward in time from the present, we aim to emphasize the constraints that arise on the inflationary models from the WMAP data. In doing so, we take seriously the indication that a value of $`\mathrm{\Omega }_0`$ slightly larger than 1 is favoured. During a slow-roll inflationary phase, the universe is described by an almost de Sitter spacetime with the $`K=+1`$ scale factor

$$a(t)\mathrm{cosh}(Ht)$$
(2)
instead of the exponential growth obtained in the flat case. It follows that the number of $`e`$-foldings between the onset and the end of inflation cannot be arbitrarily large, when we take present day cosmological parameters into account. The curvature enhances the curvature, $`\ddot{a}`$, of the scale factor in the early stage of inflation so that there exists a turn-around point in the past. It follows  that the number of $`e`$-foldings during inflation is bounded by a maximum possible number of $`e`$-foldings, $`N_{\mathrm{max}}`$, related to the present day density parameter $`\mathrm{\Omega }_0`$.

Assuming that the inflation is driven simply by a cosmological constant, the scale factor of a closed universe in this epoch is given by

$$a(t)=a_{\mathrm{min}}\mathrm{cosh}(\lambda t)$$
(3)
with $`\lambda \sqrt{\mathrm{\Lambda }/3}`$ and $`t=0`$ corresponds to the minimum expansion of the universe. It follows that $`\delta =\mathrm{\Omega }1=1/\mathrm{sinh}^2(\lambda \mathrm{t})`$. The time $`t_{}`$ at which $`k_{\mathrm{min}}=a_0H_0`$ crosses the horizon is thus obtained to be

$$\lambda t_{}=\mathrm{Arcsinh}(1/\sqrt{\delta _0}).$$
(4)
The number of $`e`$-foldings between $`t=0`$ and $`t_{}`$ is simply

$$\mathrm{\Delta }N_{\mathrm{max}}=\mathrm{ln}\frac{a_{k_{\mathrm{min}}}}{a_{\mathrm{min}}}$$
(5)
that is

$$\mathrm{\Delta }N_{\mathrm{max}}=\frac{1}{2}\mathrm{ln}\left(1+\frac{1}{\delta _0}\right).$$
(6)
It typically ranges from 2.3 to 1.5 when $`\delta _0`$ varies from 0.01 to 0.05 (see Fig. 2).

This argument also holds for negatively curved universes, in which case, the scale factor behaves as $`a\mathrm{sinh}\lambda t`$. In that case, there will not be any maximal number of $`e`$-foldings but the universe starts off from a singularity. This will enhance the trans-Planckian problem .

Another way to look at the existence of a maximum number of $`e`$-foldings in the $`K=+1`$ case, is to relate the number of $`e`$-foldings during inflation to the later history of the universe . During inflation, the comoving Hubble length $`aH`$ is decreasing, while it is increasing at all further times. A given comoving scale, $`k`$, crosses the horizon when $`k=a(t_k)H(t_k)a_kH_k`$. The number of $`e`$-foldings, $`N(k)`$, between when this scale crosses the horizon and the end of inflation is obtained to be

$$N(k)=\mathrm{ln}\left(\frac{a_{\mathrm{end}}H_{\mathrm{end}}}{a_kH_k}\right).$$
(7)
The determination of $`N(k)`$ requires that we compute

$$\frac{a_0}{a_{\mathrm{end}}}$$
(8)
which requires a complete history of the evolution of the universe from the end of inflation to today. If one assumes that the universe is matter dominated from the end of (slow roll) inflation to reheating, radiation dominated from reheating to equality and matter dominated up to now 111The results will be only marginally affected by a late-time cosmological constant as indicated by the supernova data., then

$`N(k)`$ $``$ $`62\mathrm{ln}\left({\displaystyle \frac{k}{a_0H_0}}\right)\mathrm{ln}\left({\displaystyle \frac{10^{16}\mathrm{GeV}}{V_k^{1/4}}}\right)`$ (9)
$`+{\displaystyle \frac{1}{4}}\mathrm{ln}{\displaystyle \frac{V_k}{V_{\mathrm{end}}}}{\displaystyle \frac{1}{3}}\mathrm{ln}\left({\displaystyle \frac{V_{\mathrm{end}}^{1/4}}{\rho _{\mathrm{reh}}^{1/4}}}\right)\mathrm{ln}h,`$
$`h`$ being the Hubble parameter in unites of 100 km.s-1/Mpc. For our purpose, it will be more convenient to relate $`N(k)`$ to $``$ which can easily be found to be given by

$$\mathrm{ln}66\frac{1}{3}\mathrm{ln}\left(\frac{V_{\mathrm{end}}^{1/4}}{\rho _{\mathrm{reh}}^{1/4}}\right)+\mathrm{ln}\left(\frac{10^{16}\mathrm{GeV}}{V_{\mathrm{end}}^{1/4}}\right).$$
(10)
It follows that

$`N(k)`$ $``$ $`128\mathrm{ln}\mathrm{ln}\left({\displaystyle \frac{k}{a_0H_0}}\right)`$ (11)
$`2\mathrm{ln}\left({\displaystyle \frac{10^{16}\mathrm{GeV}}{V_k^{1/4}}}\right)\mathrm{ln}h.`$
The cosmic microwave background roughly probe scales from 10 Mpc (i.e. of the order of the thickness of the last scattering surface) to $`10^4`$ Mpc (i.e. of the order of the size of the observable universe) while galaxy surveys can probe the range 1 Mpc to 100 Mpc.

Concerning the characteristic scales involved in the problem, we can assume that $`V_kV_{\mathrm{end}}`$ as long as slow rolling holds. The reheating temperature can be argued to be larger than $`\rho _{\mathrm{reh}}^{1/4}>10^{10}`$ GeV to avoid the gravitino problem  and may be pushed in the extreme case to $`10^3`$ GeV, i.e. just before the electro-weak transition so, that baryogenesis can take place. The amplitude of the cosmological fluctuations (typically of order $`2\times 10^5`$ on Hubble scales) roughly implies that $`V_{\mathrm{end}}^{1/4}`$ is smaller than a few times $`10^{16}`$ GeV and, for the same reason as above, has to be larger than $`10^{10}`$ GeV in the extreme case.

This implies that the number of $`e`$-foldings has approximately to lie between 50 and 70, which is the order of magnitude also required to solve the horizon and flatness problem when $`K=0`$ , but it can be lowered to 25 in the extreme case of thermal inflation , which we are not considering here.

In the following, we restrict our analysis to the case where reheating takes place just at the end of inflation, i.e. $`\rho _{\mathrm{reh}}=V_{\mathrm{end}}`$. It has been shown  that, in that case, the existence of a positive curvature limits the number of $`e`$-foldings to

$$N_{\mathrm{max}}(,\delta _0)=\mathrm{ln}\sqrt{\frac{\alpha }{\delta _0}}$$
(12)
with

$$\alpha =(1+\delta _0\mathrm{\Omega }_{\mathrm{rad0}})+\mathrm{\Omega }_{\mathrm{rad0}}^2$$
(13)
where $`\delta _0\mathrm{\Omega }_01`$ and $`\mathrm{\Omega }_{\mathrm{rad0}}4.17\times 10^5h^2`$ is the radiation density parameter today. Let us emphasize that the derivation of Eq. (12), as well as Eq. (6), does not make any assumption concerning what happened before the start of inflation.

This allows us to estimate the number of $`e`$-foldings allowed prior to the time at which $`k_{\mathrm{min}}`$ crossed the horizon,

$$\mathrm{\Delta }N=N_{\mathrm{max}}[(V_{\mathrm{end}}^{1/4}),\delta _0]N(k_{\mathrm{min}},V_{\mathrm{end}}^{1/4})$$
(14)
where we will let $`\delta _0`$ vary in the range $`0.010.05`$ (see Fig. 2). It can be checked, as expected from Eq. (12) that $`\mathrm{\Delta }N`$ depends slightly on $`V_{\mathrm{end}}^{1/4}`$.

We see that the allowed number of $`e`$-foldings are compatible with the requirements of structure formation. ‘Flatness’ is of course solved to the accuracy represented by $`\delta _0`$. The scenarios sketched here do not lead to as small a value of $`\delta _0`$ as is often supposed - and thereby is compatible with the best-fit WMAP data. However, for the mean value of the curvature, one obtains that

$$\delta _0=0.02\mathrm{\Delta }N_{\mathrm{max}}=1.97.$$
(15)
With standard parameter values this gives the maximum possible number of $`e`$-foldings as about $`N_{\mathrm{max}}=62+1.9764`$ \- compatible with the estimates given above, but not nearly as large as suggested in some models of inflation.

In view of the bounds on $`N_{\mathrm{max}}`$, there are also bounds on the ability of inflation to solve the horizon problem when one considers inflationary dynamics with $`K=+1`$ . Roughly speaking, the horizon problem is only solved during the inflationary era itself if the start of the inflationary era is close enough to the throat at $`t=0`$ in the expansion (see Eq. 3); that is, if the total number $`N`$ of $`e`$-foldings during inflation is close to the total $`N_{\mathrm{max}}`$ allowed, i.e. $`NN_{\mathrm{max}}`$. If inflation starts sufficiently far from that throat, then it does not matter how many $`e`$-foldings take place, they cannot provide causal connectivity on large enough scale to solve the horizon problem; rather that connectivity has to be set up in the pre-inflationary era (and is evidenced by a smooth enough small patch that is then expanded to large scales by inflation). However inflation certainly can start near enough the throat to solve the horizon problem during the inflationary era and still lead to $`\delta _0=0.02`$. Examples of non-slow rolling cases are the emergent and awakening universes  that start form an Einstein static state and are singularity free. The above formula (15) limiting the number of $`e`$-foldins will not apply in these cases.

In this note we have argued that the curvature of space may play an important role in the dynamics of inflation in its early stage. We have argued that a maximum number of 1.5 – 2.5 $`e`$-foldings can take place before the largest scale observed in the universe crossed the horizon during slow-roll inflation. However, despite the fact that the curvature term will be smaller than a few percent during all times ranging from the horizon exit of the largest observable mode to its reentry today, a positive curvature term as small as the one indicated by current data will dramatically change the early phase of inflation and hence lead to the limits on the number of $`e`$-foldings allowed as described above. These restrict the families of inflationary universes compatible with the current data, but do not exclude inflationary models even though the value of $`\delta _0`$ indicated by WMAP is larger than supposed in most inflationary scenarios.

Note a possible way around the argument. In the case that inflation is driven by a scalar field, the existence of a maximal number of $`e`$-foldings may be avoided if it exits the slow-roll regime. It was however shown that the dynamics of closed universes filled with a scalar field is chaotic  and can lead to the existence of singularity free universes  with periodic  or aperiodic  trajectories. In that case, again, a positive curvature affects the dynamics of the early phase of inflation and the slow-roll regime has to be left before $`N_{\mathrm{max}}`$ is reached, which is also at odds with the standard picture.

Let us finally consider further features of spherical universes that may turn to out be useful in understanding the WMAP data. Surprisingly, the WMAP angular correlation function seems to lack signal strength on angular scales larger than 60 degrees . This may indicate a possible discreteness of the initial power spectrum, as expected e.g. from non-standard topology. The case for topology becomes most interesting in the case of spherical spaces . First, the spectrum of the Laplacian in spherical spaces is always discrete . Second, and as emphasized in Ref. , a non trivial topology is most likely to be detectable in the case of spherical spaces and examples of the resulting observational effects are discussed in Ref. .

The same conclusion cannot be reached for negatively curved universe. It is unlikely that we could detect a compact hyperbolic space because, for an observer at a generic point, the topology scale is comparable to the curvature radius or longer. The possibility of a detection would require the observer to sit near a short closed geodesic .

In conclusion, if the indication that $`\mathrm{\Omega }_0`$ differs from unity at a level of a few percent is confirmed then it may imply important results concerning the dynamics of the early stage of our universe. The $`K=+1`$ background de Sitter model (3) differs fundamentally from the scale-free de Sitter model in the $`K=0`$ frame - not least in terms of being geodesically complete. The WMAP observation suggest the former may be the appropriate model to use in investigating early inflationary dynamics.

Acknowledgements: JPU thanks the Department of Mathematics and Applied Mathematics of the University of Cape Town for hospitality. UK acknowledges financial support from the University of Cape Town. We thank Roy Maartens for his comments.


text:
# Features in Evanescent Aharonov-Bohm interferometry.

## I Introduction

The Aharonov-Bohm(AB) interferometer is one of the most exciting topics in mesoscopic physics research, not only because of the fundamental physical insight into questions of quantum non-localitysteinberg  but also for being the basis of many novel phenomena at mesoscopic scales. Some of which are, persistent currents in normal metal ringswebb ; levy , Aharonov-Bohm oscillationsgia , Fano resonanceskobay , and current magnificationcolin\_cm ; deo\_stub . In this work we analyze the local currents and conductance in this interferometer but for evanescent mode propagation. Evanescent means the incident energy of charge carrier’s in the leads attached to the interferometer is less than the potential $`V`$ which characterizes the ring throughout (see FIG. 1, for a schematic representation of our system). Thus when an electron with energy $`E<V`$ impinges on the ring it has to tunnel out of the ring. In certain situations we also consider a quantum well inserted in one of the arms of the interferometer. Only in the well and the connecting leads does the electron propagate with a real wave-vector $`k=\sqrt{E}`$, elsewhere it propagates with a complex wave-vector $`k=i\kappa =i\sqrt{VE}`$. In such situation the contribution to conductance arises simultaneously from two non-classical effects namely, Aharonov-Bohm effect and quantum tunellingswitch . Such situation can arise when the transverse width of the ring is much less than the connecting ideal wires. In this case due to the higher zero point energy arising from transverse confinement, the fundamental sub-band minima in the ring will be at higher energy than the value of few sub-band minima in the ideal connecting wires. Now a situation can arise where propagating modes in the wire have energy less than the minimum propagating sub-band energy of the ring system. Thus the electron propagating in a lower sub-band of the ideal wire feels a barrier to it’s motion (experiences an effective potential barrier $`V`$, arising solely because of the mismatch of the zero-point energies) and tunnels across the system via evanescent mode propagation. For simplicity we have restricted our calculation to single channel case. We are interested only in coherent transport, i.e., absence of dephasing. Using the wave guide method we have solved for all the scattering coefficients and wave amplitudes in the ring, required to calculate the local currents in the arms of the interferometer. For details we refer to Refs.xia ; colin\_cm ; deo\_stub ; coupled . Expressions for the scattering amplitudes can be obtained readily with the use of a simple MAPLE programmaple . However analytical forms of these expressions are too lengthy. Therefore we analyse our results by suitable plots generated from the obtained MAPLE expressions.

Herein we particularly analyze three configurations of the Aharonov-Bohm interferometer and see how the conductance and local currents in the two arms respond to the configurational variations. The first configuration is an AB ring at potential $`V`$ attached to two leads at zero potential. In the second configuration we consider an AB ring at potential $`V`$ throughout apart from the quantum well which is at zero potential, and in the final configuration we consider an AB ring with potential $`V`$ only in the upper arm in addition to the quantum well but the lower arm is set free($`V=0`$). In particular, we show that the first configuration does not support current magnification, while the second configuration supports it. Further the second configuration is marked by absence of Fano resonances as opposed to the third configuration. The magneto conductance shows systematic behavior as a function of flux in the evanescent regime. The detailed analysis of these phenomena is given in the following sections.

## II Evanescent Aharonov-Bohm interferometry

In the first configuration, we compare and contrast the cases of transport in the pure evanescent and pure propagating mode regime, with emphasis on it’s conductance and local currents. In these cases the well in upper arm is absent. In FIG. 2(a), we plot the conductance in dimensionless form which is just the transmission $`G=(2e^2/h)T`$, as a function of the magnetic flux for different Fermi energies below the barrier potential, i.e., in the evanescent regime. The total length of the ring $`l(=lu1+lw+lu2+ld)`$ is taken to be unity and all other physical parameters of the ring are scaled with respect to this (we have set $`\mathrm{}=2m=1`$ throughout), for example energy in dimensionless form is $`El^2E`$, and barrier potential $`Vl^2V`$. As expected the conductance is periodic in flux with period $`\varphi _0`$, the fundamental flux quantum $`\frac{hc}{e}`$. We observe that the small field differential magneto conductance is always negative. In FIG. 2(b), we plot the conductance as a function of the magnetic flux for energies above the barrier potential, i.e., in the propagating regime. We see that the small field differential magneto-conductance (MC) may be positive/negative. From these two figures we observe the systematic behavior of magneto-conductance slopes for pure evanescent wave transport in contrast to the unsystematic nature of magneto-conductance slopes for pure propagating wave transport. The systematics in transport refer to the slope of the magneto-conductance for small fields (i.e., the small field differential magneto-conductance), which for pure evanescent mode transport is always negative irrespective of the position of impurities or system parameters. Unsystematic behavior refers to this slope being very sensitive to system parameters or impurities, i.e., whether it is positive or negative cannot be predicted a prioriswitch .

The reason we observe the systematic nature of magneto-conductance slopes has to do with the harmonics. The harmonics are calculated as followscolin\_de -

$`a_n={\displaystyle \frac{1}{\pi }}{\displaystyle _0^{2\pi }}Tcos(n\varphi )𝑑\varphi `$ (1)
Herein, $`T`$ represents the total transmittance and $`\varphi `$ the enclosed flux. The magnitude of $`nth`$ order harmonic corresponds to the contribution from electronic paths which encircle the flux $`n`$ times. In the case of evanescent mode propagation the first harmonic always dominates the rest, as is seen in FIG. 3(a). Owing to the decaying nature of evanescent modes, the higher harmonics are weighted by a factor $`e^{n\kappa l}`$. Thus higher harmonics, in the evanescent regime which have to encircle the ring many more times naturally decay faster(exponentially), therefore first harmonic completely dominates over the rest, and thus we see systematic behavior in the conductance and hence small field magneto-conductance slope is always negative. This definite feature (zero flux conductance is always greater than the conductance for increasing values of flux upto $`\varphi _0/2`$, i.e., the small field regime) may be utilized in some device application which has been suggested earlierswitch . For propagating waves any harmonic can dominate over the rest (i.e., amplitude of any higher harmonic can dominate the lower one). This is amply clear from FIG. 3(b), and thus we see unsystematic nature of slopes for small field differential magneto-conductance. In the case of propagating modes the nature of the slopes is very sensitive to system parameters and defects if presentswitch .

Still considering transport in the Aharonov-Bohm interferometer without the well, we plot in FIG’s. 4(a) and 4(b), the output current density in dimensionless form)colin\_cm  as well as the current densities in upper $`J_u`$ and lower arms $`J_d`$ as a function of the incident Fermi energy for both evanescent(FIG. 4(a)) and propagating(FIG. 4(b)) modes of transport. In FIG. 4(a), we observe that the currents in the upper and lower arms are individually less than the output current, thus indicating the absence of current magnification. Now what is current magnification? In the AB interferometer, we see that the dimensionless output current density I (=T the transmission probability) around a small Fermi energy interval, flows through the system when voltage bias ($`\mu _1\mu _2=eV`$) is applied. The upper and lower arms of the ring are of different lengths such that current densities $`J_u`$ and $`J_d`$ flow in these with $`J_uJ_d`$, but $`I=J_u+J_d`$, i.e., obeying Kirchoff’s law of current conservation. In particular we see that for some particular ranges of Fermi energy $`J_u`$ or $`J_d`$ can be much larger than $`I`$. Current conservation thus dictates $`J_d`$ or $`J_u`$ to be negativecolin\_cm . This property that current in one of the arms is larger than the transport current is referred to as current magnification effectdeo\_stub . It is an effect without any classical analog in a DC circuitdeo\_stub . The negative current flowing in one arm can be interpreted as a circulating current that flows continually in the ring. The magnitude of the negative current in one of the arms flowing against the direction of applied current is taken to be that of the circulating current. When negative current flows in the upper arm the circulating current direction is taken to be anti-clockwise (or negative) and when it flows in the lower arm the circulating current direction is taken to be clockwise (or positive)colin\_cm . It has been suggested that due to this effect one can observe large orbital magnetic moment of a ring in absence of magnetic field, however, in presence of transport currents (in non-equilibrium state). This effect has been extended to thermal currentsmosk  and to spin currents in the presence of Aharonov-Casher fluxchoi . One can clearly see that current magnification effect is absent in the evanescent regime, as in FIG. 4(a), while being present for the propagating regime, as in FIG. 4(b). In particular from FIG. 4(b), one can see that for Fermi energy range $`\mathrm{80.0.0}<E<100.0`$, $`J_d>I`$ and $`J_u`$ is negative, while $`J_u>I`$ for $`100.0<E<121.0`$ and $`J_d`$ is negative. Preferentially one observes current magnification around the anti-resonances of the ringdeo\_stub . In all the figures we have taken dimensionless parameters, the current densities $`J_u`$ and $`J_d`$ plotted in the figures are in their dimensionless form given by dividing $`J_u`$ and $`J_d`$ by $`\frac{e\mathrm{}k}{m}`$.

## III Evanescent AB interferometry in presence of resonant states

Having established that evanescent mode consideration shows systematics in MC but absence of current magnification while transport in the propagating regime shows opposite behavior, no systematics in MC but current magnification is evident. Now we turn our attention to the problem as we had formulated it in our second configuration, i.e., evanescent mode transport in an Aharonov-Bohm ring with a quantum well in it’s upper arm to simulate resonances. The need to incorporate resonant states arises from the fact that presence of unavoidable defects and other localized resonant states in such small devices can mar switching action and hamper device performanceswitch . To analyze this if a well sustaining a few resonant states can be incorporated in the Aharonov-Bohm interferometer, then device performance can be tested. To this end we justify the inclusion of a quantum well in one of the arms. The parameters defining the quantum well are mentioned in the caption of FIG. 5, in dimensionless form, $`l_w`$ represents width of the quantum well, while $`V`$ defines the potential barrier’s, in the well region the particle is assumed to be free

In FIG. 5, we plot the output current density and the current densities in upper and lower arms as a function of the Fermi energy for asymmetric arm lengths (as in FIG. 1). We encompass three resonant states of the well. The length and height of barrier potential are mentioned in the figure caption. We assume the particle to be free in the well, i.e., in the well, potential is zero. The resonances in FIG. 5 are all of Breit-Wigner type in contrast to transmission resonances for propagating mode transport which can be either of Breit-Wigner or Fano typedeo\_stub . In contrast to evanescent transport without a quantum well, at Fermi energies comparable to barrier potential but still in the evanescent regime, current magnification is seen. The current density($`J_u`$) of the upper arm is more than output current while that of lower arm($`J_d`$) is negative. In this regime a circulating current continually flows in the ring even for evanescent wave transport in the region outside the well. We have also seen separately that if potential well supports many resonant states(second configuration) then also current magnification occurs at Fermi energies comparable to barrier potential and interestingly it is always to the right of the last resonance (E$``$V).

To see the effect of resonant states on the MC we plot conductance for different Fermi energies on both sides of the second resonant peak of FIG. 5, in FIG. 6(a) and FIG. 6(b). In FIG. 6(a) we plot the magneto-conductance as a function of flux for Fermi energies to right of the second resonance and left of third resonance as in FIG. 5. In FIG. 6(a) we see that the small field differential magneto-conductance is still negative and hence systematic behavior can still be seen, while in FIG. 6(b) we plot the magneto-conductance in the energy window in-between first and second resonances, and see that the small field differential magneto-conductance is positive. Thus we can draw some uniform conclusions as to where systematic action is to be seen. In particular nature of slope of magneto-conductance is always negative for Fermi energies to the left of first resonance, and then the nature of this slope (being negative or positive) alternates as we cross each successive resonance. This systematic is possible due to first, a gradual change in the transmission phase by $`\pi `$ across the well while crossing the resonance and second due to the absence of transmission zeros in this system. The transmission zeros can lead to abrupt change in phase by $`\pi `$ of the conductanceimry ; fano\_deo . This behavior is unlike that shown in case of transport via propagating modes as seen in FIG. 2(b), where the systematics in MC cannot even be predicted due to presence of transmission zeros.

## IV Question of Fano resonances

Recently Fano resonances have been observed in a mesoscopic AB ring with a quantum dot embedded in it’s upper armkobay . Fano resonances arise when a discrete set of states are coupled to the continuumkobay . Especially Fano resonances are caused by interference of two alternative pathsimry , a resonant and a non-resonant one. In the geometry considered a resonant path is via a quantum well (having resonant energy levels) and a direct non-resonant path is via the lower arm of the ring (FIG. 1). Transmission resonances are generally of two types Breit-Wigner and Fano. Fano line-shapes are asymmetric as opposed to the symmetric Breit-Wigner line-shapes. In contrast to Breit-Wigner forms, Fano resonances are charecterized by zero-pole pair structure in the complex energy plane of the scattering amplitudeporod ; bagwell .

In FIG. 7(a) we have plotted the conductance (in effect the transmission coefficient in dimensionless form) as function of the Fermi energy. We see resonant transport around the quasi bound states of the well (around $`E=39.0`$). The line shape is that of Breit-Wigner form. We do not observe any Fano line-shapes as long as $`E<V`$. As mentioned earlier Fano line-shapes arise from the interference of two alternative paths, one resonant and another direct. Our results clearly indicate when the transport across the direct path is via evanescent modes, then Fano line-shapes are absent. To recover Fano line-shapes, direct path should be a propagating one. To see this explicitly we set potential in the lower arm to be zero (configuration three). So that electron traverses the lower arm as a propagating wave. For this case we have plotted in FIG. 7(b), the conductance versus Fermi energy for the same physical parameters as in FIG. 7(a). We clearly see the Fano line-shapes, i.e., line-shape is asymmetric and transmission exhibits zero around the same Fermi energy interval. The results presented in FIG.’s 7(a) and 7(b) are in the absence of magnetic field. Increasing the magnetic field which amounts to breaking the existing time reversal symmetry will lead to lifting of zeroes except when the flux piercing the loop is an integral or half integral multiple of the flux quantumbroken .

## V Conclusions

To conclude, when electron traverses the entire circumference of the ring as an evanescent wave, the current magnification effect is absent. In this regime small field differential magneto-conductance is always negative. Inclusion of a well in an otherwise evanescent transport supports current magnification. The magnetoconductance shows systematic behavior even in this regime. The initial slope of the magnetoconductance alternates as we cross successive resonance peaks. The transmission resonances one observes in this set up (resonant states in an evanescent AB interferometer) are of Breit-Wigner type in contrast to the Fano line shapes expected because of the coupling between a resonant and non-resonant path. One recovers Fano resonances if transport through non-resonating path is via propagating modes. This fact emphasizes the need of a direct propagating path for observance of Fano line shapes in transmission.


text:
# Effects of Wandering on the Coalescence of Black Hole Binaries in Galactic Centers

## 1 Introduction

Observations indicate that most galaxies harbor supermassive black holes (BHs) at their centers (e.g., Magorrian et al. 1998). Since galaxy mergers are inevitable, there arises the question of the dynamical fate of the binary black holes (BBHs) that result from these mergers. This problem was first addressed by Begelman, Blandford & Rees (1980) (for recent work, see, e.g., Milosavljević & Merritt 2001, 2002, Yu 2002). Dynamical friction causes the black holes in the original galaxies to sink to the center of the merged galaxy, and then become bound in a binary system. The semi-major axis of the binary, $`a`$, continues to shrink owing to dynamical friction until it reaches the point at which the binary becomes hard: $`a_hGm_2/4\sigma ^2`$, where $`\sigma `$ is the one-dimensional stellar velocity dispersion, $`G`$ is the gravitational constant, and $`m_2`$ ($`<m_1`$) is the mass of the lighter of the two black holes (following Quinlan 1996). After this stage, the dominant mechanism by which the binary shrinks and loses energy is by three-body interactions with stars that pass within a distance $`a`$ of the center of mass of the binary. The binary transfers energy to the stars by ejecting them at a much higher velocity. Once it shrinks to the point at which energy loss to gravitational waves becomes dominant, the binary can quickly coalesce.

Begelman, Blandford & Rees (1980) conjectured that hardening of a tight binary by three-body interactions would continue only until it had interacted with all the stars contained within its loss cone (Frank & Rees 1976); if this was not sufficient to drive it to the gravitational wave emission stage, the binary would stall at a certain separation, and harden very slowly at the two-body relaxation timescale (the timescale at which the loss cone can be repopulated), which turns out for many actual galaxies to be longer than the Hubble time. The lack of observational evidence for the existence of massive BBHs has led to suggestions that BBH merging may be facilitated by other processes, such as the loss of angular momentum to gas (e.g., Begelman, Blandford & Rees 1980, Gould & Rix 2000, Armitage & Natarajan 2002), or triaxiality of the host galaxy (Yu 2002).

It is therefore of interest to inquire whether stellar dynamical processes alone would be sufficient to drive the BBH to the gravitational wave emission stage. Quinlan & Hernquist (1997) pointed out the importance of the wandering of the binary to its hardening. In this paper, we investigate the shrinking of a hard black hole binary in a spherically symmetric galaxy, and ask particularly how Brownian motion of its center of mass can contribute to its continued hardening beyond the stage of loss cone evacuation. In §2, we set up an analytic model that appears to account for the results of N-body simulations described in §3. We comment on the impact of BBH wandering in §4, and derive characteristics of the merger in §5, on the assumption that hardening proceeds at the rate deduced in §2 and §3. We summarize our findings in §6.

## 2 The Model

Consider two black holes of masses $`m_1`$ and $`m_2`$. Let $`\mu _{12}=m_1m_2/(m_1+m_2)`$ be the reduced mass of this BBH system. Also, let $`m_2`$ be the lighter of the two BHs, with the mass ratio being $`q=m_2/m_11`$.

For illustrative purposes, we take the host cluster of stars to be described by a Plummer model of total mass $`M`$ and length parameter $`r_0`$. Thus, the density and potential profiles are given, respectively, by

$$\rho (r)=\frac{3Mr_0^2}{4\pi }\frac{1}{(r^2+r_0^2)^{5/2}},$$
(1)
$$\mathrm{\Phi }(r)=\frac{GM}{(r^2+r_0^2)^{1/2}},$$
(2)
where $`G`$ is the gravitational constant and $`r`$ is the radial distance from the center of the stellar system, which is taken as the origin. In this case, the one-dimensional stellar velocity distribution in the core is given by $`\sigma ^2=GM/6r_0`$.

The phase space distribution function of the stellar system, $`f`$, depends in general both on position $`\stackrel{}{r}`$, and stellar velocity $`\stackrel{}{v_a}`$, and is defined such that $`f(\stackrel{}{r},\stackrel{}{v_a})d^3\stackrel{}{r}d^3\stackrel{}{v_a}`$ is the mass in stars in the phase space volume $`d^3\stackrel{}{r}d^3\stackrel{}{v_a}`$. For the spherically symmetric Plummer model, we assume $`f`$ to be a function of the relative energy per unit mass $``$ only (and independent of specific angular momentum), where $`=\frac{1}{2}v_a^2\mathrm{\Phi }(r)=\mathrm{\Psi }(r)\frac{1}{2}v_a^2`$, $`\mathrm{\Psi }(r)=\mathrm{\Phi }(r)`$ being the relative potential. The distribution function is then (e.g., Binney & Tremaine 1987)

$$f()=\frac{96}{7\sqrt{8}\pi ^3}\frac{Mr_0^2}{(GM)^5}^{7/2}.$$
(3)
Consider an encounter between the BBH and a star of mass $`m_{}`$. Let the impact parameter be $`p`$ and the initial relative velocity (when the star is at infinity) be $`V_0`$. If $`r_{min}`$ is the minimum distance reached by the star to the BBH, and $`V_{max}`$ is the velocity at that point, then the angular momentum is $`l=pV_0=r_{min}V_{max}`$. The reduced mass of the BBH-star system is $`\mu m_{}`$, assuming the BHs to be much more massive than the stars. Then by conservation of energy,

$$E=\frac{1}{2}m_{}V_{max}^2\frac{G(m_1+m_2)m_{}}{r_{min}}=\frac{1}{2}m_{}V_0^2$$
(4)
Eliminating $`V_{max}`$ and rearranging,

$$p^2=r_{min}^2+\frac{2G(m_1+m_2)}{V_0^2}r_{min}.$$
(5)
We set $`r_{min}=a`$, the semi-major axis of the BBH, implying that an encounter occurs (i.e., according to our picture, the star interacts with the BBH and gets thrown out by the slingshot effect, while $`a`$ shrinks) if $`p^2a^2+\frac{2G(m_1+m_2)}{V_0^2}a`$.

The number of stars per unit volume with velocity $`v_a`$ is $`f(v_a)d^3v_a/m_{}`$. Hence, the number of encounters per unit time is

$$\frac{1}{t_{en}}=\frac{f(v_a)}{m_{}}d^3v_a\pi p^2V_0,$$
(6)
where $`V_0=|\stackrel{}{v}\stackrel{}{v}_a|`$, $`\stackrel{}{v}`$ being the velocity of the BBH. In order to take into account the wandering of the BBH’s center of mass, this expression must be averaged over $`\stackrel{}{v}`$ (for a similar calculation, see Binney & Tremaine , §8.5). It was shown in Chatterjee, Hernquist & Loeb (2002a, b) that the distribution of each component of $`\stackrel{}{v}`$ is Gaussian, so that the distribution of the magnitude, $`v`$, is $`W(v)=\frac{2v^2}{\sqrt{2\pi }\sigma _{bh}^3}e^{v^2/2\sigma _{bh}^2}`$, with $`\sigma _{bh}^2=\frac{2GMm_{}}{9r_0(m_1+m_2)}`$. Thus, averaging over $`W`$,

$$\frac{1}{t_{en}}=\frac{2\pi }{\sqrt{2\pi }m_{}\sigma _{bh}^3}v^2e^{v^2/2\sigma ^2}\left(a^2+\frac{2G(m_1+m_2)}{V_0^2}a\right)V_0f(v_a)d^3v_a𝑑v,$$
(7)
or

$$\frac{1}{t_{en}}=\frac{2\pi }{\sqrt{2\pi }m_{}\sigma ^3}v^2e^{v^2/2\sigma ^2}\left(a^2|\stackrel{}{v}\stackrel{}{v_a}|+\frac{2G(m_1+m_2)}{|\stackrel{}{v}\stackrel{}{v_a}|}a\right)A(\psi v_a^2/2)^{7/2}d^3v_a𝑑v,$$
(8)
where $`A=\frac{96}{7\sqrt{8}\pi ^3}\frac{Mr_0^2}{(GM)^5}`$. Note that in the core, $`\psi `$ obtains an approximately constant value, $`\psi GM/r_0`$.

The above integral can be simplified using the approximation $`|\stackrel{}{v}\stackrel{}{v_a}|v_a`$, since the BBH moves slowly relative to the stars when its mass is much larger than that of an individual star. The integration then gives,

$$\frac{1}{t_{en}}=\frac{32\pi ^2A\psi ^{11/2}}{99m_{}}\left(a+\frac{11}{2}\frac{G(m_1+m_2)}{\psi }\right)a.$$
(9)
The BBH’s binding energy is $`E=\frac{Gm_1m_2}{2a}`$. If each star on average extracts energy $`\mathrm{\Delta }E`$, where $`\mathrm{\Delta }E=G\mu _{12}m_{}/a`$, is roughly the potential energy of the BBH-star system at the point of closest approach (giving $`\mathrm{\Delta }E/E=\frac{2m_{}}{m_1+m_2}`$), then

$$\frac{\mathrm{\Delta }E}{\mathrm{\Delta }t}=\frac{Gm_1m_2m_{}}{m_1+m_2}\frac{32\pi ^2A\psi ^{11/2}}{99m_{}}\left(a+\frac{11}{2}\frac{G(m_1+m_2)}{\psi }\right)=\frac{Gm_1m_2}{2a^2}\frac{da}{dt},$$
(10)
and simplifying,

$$\frac{1}{a^2\left(a+\frac{11G(m_1+m_2)}{2\psi }\right)}\frac{da}{dt}=\frac{64\pi ^2A\psi ^{11/2}}{99(m_1+m_2)}.$$
(11)
Thus, there are two limiting regimes for the temporal evolution of $`a`$:

* If $`a11G(m_1+m_2)/2\psi `$, then $`\frac{1}{a^3}\frac{da}{dt}=\frac{C}{m_1+m_2}`$, where $`C`$ is a constant, and we get $`1/at^{1/2}`$.
* If $`a11G(m_1+m_2)/2\psi `$, then $`\frac{d}{dt}(\frac{1}{a})=`$ a constant, and we get $`1/at`$.

For intermediate values, we get a power law, $`1/at^d`$, with $`0.5d1.0`$.

Now, $`11G(m_1+m_2)/2\psi =\frac{11}{2}\frac{m_1+m_2}{M}r_0`$. The BBH becomes hard at $`a_h\frac{3}{2}\frac{m_2}{M}r_0`$, where $`m_2`$ is the lighter of the two BHs. Thus, after the BBH becomes hard, we always tend to the second case above, $`1/at`$.

In our numerical experiments (to be described in the next section), we have found temporal power-law indices for $`1/a`$ which are always bounded between $`0.5`$ (for heavy BBHs) and 1 (for light BBHs). The above simplified theory predicts a power-law index of 1 for hard binaries. However, our experiments showed that heavy BBHs (which contain more than a few percent of the mass of the Plummer galaxy) cause the central density of the galaxy to drop quickly, and so it is not appropriate to assume, as we did above, that the stellar distribution function remains unaffected by the BBH. In particular, the left-hand-side of equation (11) needs to be $`\frac{1}{a^2(a+c_1)}\frac{da}{dt}`$, where $`c_1`$ is now some undetermined term. If $`c_1`$ is not small compared to $`a`$, then it is possible to get $`1/at^d`$, with $`d<1`$. This description is suggestive for why we obtained such limiting behaviors, and it is due to the form of the expression (5) for the impact parameter being a sum of two terms.

For low mass BBHs, the distribution function is not much affected by the hardening of the binary, and so we may proceed to derive the temporal slope of the $`1/a`$ line,

$$s=\frac{d}{dt}\left(\frac{1}{a}\right)=\frac{256\sqrt{2}}{21\pi }\left(\frac{GM}{r_0^5}\right)^{1/2}.$$
(12)
Note that this result is independent of the masses of the two BHs.

We may compare the above value with the hardening rate presented in other work in a slightly different form: $`\frac{d}{dt}\left(\frac{1}{a}\right)=\frac{G\rho _0H}{\sigma }`$, where $`H`$ is a dimensionless hardening rate, $`\rho _0`$ is the central stellar density, and $`\sigma `$ is the one-dimensional stellar velocity dispersion in the core. Comparing with the above equation, we obtain $`H=9.4`$. This is very similar to the hardening rate observed in simulations by Milosavljević & Merritt (2001). It also compares well with the value $`H15`$ (independent of BH mass ratio and orbital eccentricity), which is the hardening rate derived from three-body scattering experiments in the limit of a very hard binary (Quinlan 1996, Hills 1983).

In our simulations, we made the choice $`r_0=3\pi /16`$ and $`GM=1`$, for which one gets $`s=\frac{d}{dt}(\frac{1}{a})=20.6`$. Note, however, the approximations that have been made in the derivation. We have assumed that the distribution function of the stars remains unchanged, which is clearly false: the central density of the stellar system falls, as shown in the next section. For low mass BBHs (say $`m<0.01`$), however, the decline is not as much as for heavy BBHs, and we expect the inverse semi-major axis of the BBH to increase at a constant rate, with $`s`$ given roughly by the expression above.

We have assumed $`\mathrm{\Delta }E/E=2m_{}/(m_1+m_2)`$, a form which is in reasonably good agreement with results of scattering experiments at several mass ratios (see Roos 1988, Hills 1983). It is, of course, an expression for the average energy extracted by a passing star, and need not accurately reflect the actual process of energy extraction for all impact parameters. Given that the precise form of this expression suitable for use in the massive BBH case is uncertain (for example, Hills & Fullerton  find that for equal-mass binaries, the right-hand side of the above expression may be multiplied by 0.7), it is perhaps optimistic to expect this derivation to give the correct slope $`s`$.

## 3 Results from Numerical Experiments

The numerical simulations were carried out using the SCFBDY program, which is described in detail in Quinlan & Hernquist (1997) and Chatterjee, Hernquist & Loeb (2002a). In our code, each BH interacts with the other BH and the stars through the $`1/r^2`$ force, but stars interact with other stars through a mean field expansion of the gravitational potential that is updated self-consistently with time (note that the simulation thus does not treat star-star relaxation accurately \[see Hernquist & Barnes 1990\]; this is justified since the two-body relaxation timescale in real galaxies is much longer than the Hubble time). The particles are moved with individual step-sizes using the fourth-order integrators of Aarseth (1994): NBODY1 for the stars and NBODY1 or NBODY2 for the BHs.

The mean field expansion of the gravitational potential is carried out in a set of basis functions that includes both spherical and non-spherical terms (Hernquist & Ostriker 1992). However, since the galaxy remained nearly spherical during the integration, it made not much difference to the hardening of the BBH whether a few or zero non-spherical terms were included in the potential expansion.

We verified that numerical relaxation was not an important factor in our simulations by forcing the stars to interact only with a fixed background potential, and not with each other, and finding our results to be similar to those obtained using the full SCFBDY code. The insignificance of numerical relaxation in our simulations is further supported by experiments in which we forced the center of mass of the BBH to remain fixed, as described below.

Our results are presented in a system of units in which $`G=M=1`$, and $`r_0=3\pi /16`$, so that the energy of the initial galaxy without the BHs is -1/4 (Heggie & Mathieu 1986). The stars (of equal mass $`m_{}=1/N`$, where $`N`$ is the number of stars) are initially distributed according to the Plummer distribution function, and the black holes are started out on nearly circular orbits at $`\stackrel{}{r}`$ and $`\stackrel{}{r}`$, with $`r=0.3`$. (We tried a range of initial conditions for the BHs, including elliptical orbits; it made no difference to the hardening of the binary.) We have done a range of experiments, with binary mass ranging from $`m_1=m_2=0.00125`$ to $`m_1=m_2=0.01`$, and with particle numbers of up to $`N=400,000`$.

Figure 1 shows the results of such an integration with $`m_1=m_2=0.00125`$, and 400,000 stars. Shown are the separation between the BHs, $`R`$, the inverse semi-major axis of the BBH orbit, $`1/a`$, the eccentricity of the orbit, $`e`$, and the central density, $`\rho `$, of the stellar system averaged over a sphere containing a mass in stars somewhat larger than the mass of the BBH. The BHs become bound at $`t=77.3`$. The BBH becomes hard once $`a`$ shrinks to the value $`a_h=Gm_2/4\sigma ^2=1.104\times 10^3`$, or when $`1/a_h=905`$. As can be seen, the rate of change of $`1/a`$ is nearly constant after the binary has become hard; namely, $`1/at^d`$ with $`d1`$.

Figure 2 shows the results of a similar integration with $`m_1=m_2=0.02`$, and 200,000 stars. The BHs become bound at $`t=4.4`$, and the BBH becomes hard when $`a=a_h=0.018`$, or when $`1/a=1/a_h=56.6`$. The rate of hardening is clearly not constant after this stage. Indeed, this is at the opposite extreme from the previous case, with the hardening being well fitted by $`1/at^d`$ with $`d0.5`$.

For intermediate values of $`m_1`$ and $`m_2`$, the hardening occurs at the rate $`1/at^d`$ with $`0.5\stackrel{<}{}d\stackrel{<}{}1.0`$. We have not found values of $`d`$ significantly outside this range. This is consistent with the explanation given in §2. Notice that for $`m_1=m_2=0.00125`$, the density at the position of the BBH center of mass drops gently over the course of its hardening. This is expected since the BBH expels stars from the center as it hardens. The magnitude of the drop is small, however, and thus to a first approximation, the stellar distribution function remains unaffected. On the other hand, the heavy BBH changes the distribution function to a greater extent, as is clear from the bottom right panel of Figure 2. Thus, in the former case, the theory developed in the previous section holds, and the binary hardens at a constant rate; in the latter case, we get a qualitatively different behavior, for the reasons outlined in §2.

The rate of hardening over time is close to constant in three cases for which we ran simulations: $`m_1=m_2=0.00125`$, $`m_1=m_2=0.0025`$, and $`m_1=m_2=0.005`$. Like Quinlan & Hernquist (1997), we find that the slope of the $`1/a`$ line, given by $`s=\frac{d}{dt}(1/a)`$, varies with the number of stars. As the number of stars is increased, $`s`$ falls systematically, until there are roughly 200,000 – 400,000 stars, when it stabilizes to a particular value, $`s_0`$. This decrease in the slope as the number of stars increases has also been observed by Makino (1997) and Hemsendorf, Sigurdsson & Spurzem (2002), though not by Milosavljević & Merritt (2001). Quinlan & Hernquist (1997) explained it as resulting from two competing factors governing the wandering of the BBH’s center of mass: on the one hand, as the number of particles, $`N`$, rises, the binary wanders less from the center, and thus it reduces the density at the center more, and the hardening rate decreases; on the other hand, the restoring force of the stellar potential goes down as the central stellar density drops, causing an increase in the BBH’s wandering. At a particular value of $`N`$, there is a balance between these two factors, and the hardening rate stops changing for higher values of $`N`$.

Quite apart from the variation of $`s`$ with $`N`$, we found that the limiting value of $`s`$ changes also with the mass of the BBH. We found that for $`m_1=m_2=0.005`$, $`s`$ stabilizes to a value of $`s_08`$; for $`m_1=m_2=0.0025`$, $`s_011`$; and for $`m_1=m_2=0.00125`$, $`s_015`$. Recall that $`s=20.6`$ according to the simplified theory of §2. That value is expected to be obtained in the limit that the stellar distribution function remains unchanged by the influence of the binary. Our numerical experiments indeed reveal that this becomes a better approximation as the binary mass is reduced. However, the question remains whether this value is actually reached for still less massive BBHs.

The orbital eccentricity of the BBH in its hard state varies according to the mass of the binary. For massive binaries, such as for the case $`m_1=m_2=0.02`$ (Figure 2), the BBH orbits are circular to a good approximation. For smaller masses of the BBH, the eccentricity is larger: e.g., for $`m_1=m_2=0.00125`$ (Figure 1), the eccentricity in the hard state is roughly 0.3, and trending very slightly upward. These values are consistent with Quinlan & Hernquist (1997) and Milosavljević & Merritt (2001), and with the results of scattering experiments in the restricted three-body approximation (Quinlan 1996), according to which the BBH eccentricity hardly grows if the initial eccentricity is moderate ($`e_00.3`$), and does not grow by much if it is larger, before the gravitational radiation stage sets in. However, Hemsendorf, Sigurdsson & Spurzem (2002) and Aarseth (2002) have recently found that the eccentricity may grow to values close to 1, resulting in rapid coalescence by the emission of gravitational radiation (see §5). The reason for this discrepancy is not clear, and may have to do with the initial conditions chosen for the simulations.

It is of interest to know the extent of wandering of the BBH. According to the formulae developed in Chatterjee, Hernquist & Loeb (2002 a, b), the mean squared position and velocity of the center of mass of the BBH in the hard stage should, for the case of the Plummer potential, be given, respectively, by

$$r_{cm}^2=\frac{2r_0^2m_{}}{3(m_1+m_2)},$$
(13)
$$v_{cm}^2=\frac{2GMm_{}}{3r_0(m_1+m_2)}.$$
(14)
This assumes that the stellar distribution function is unchanged, which is not the case for a massive BBH; in fact, the wandering should be larger than indicated by these values, since the restoring force on the BBH drops as the central stellar density is diluted by the BBH. However, the above expressions provide a reasonable description for the case of light BBHs. Indeed, we find from simulations that for such cases, the wandering is consistently larger than the values provided by the above formulae, but not by more than a factor of $`2`$ (see Table 1). The reasons for the enhanced wandering could still be the somewhat reduced stellar density at the center of the cluster, as well as the enhanced random motion of a binary compared with that of a point mass, as discussed by Merritt (2001). Thus, the binary is able to interact with a larger region of stellar phase space than calculated according to the simple model of §2.

## 4 The Impact of Wandering on the Hardening of the BBH

Quinlan & Hernquist (1997) reported an experiment in which the center of mass of the binary was held fixed at the origin by the application of an artificial constraint force: they found that hardening halted soon after the binary became hard, following the evacuation of the loss cone. We have confirmed this result in our own experiments (see Figure 2 for an example): the BBH becomes hard when its semimajor axis reaches $`a_h`$, and continues hardening by three-body interactions with the stars until $`a=a_{lc}(a_{lc}<a_h)`$, when the loss cone is empty. These experiments show the importance of wandering in keeping the BBH hardening. (Incidentally, the fact that hardening stops in this case is another indication that numerical relaxation is not causing spurious refilling of the loss cone in our simulations.)

According to traditional theory (e.g., Begelman, Blandford and Rees ), hardening of the binary should stop after the BBH has interacted with all the stars in its loss cone, as long as its wandering is unimportant. However, there is no sign in our experiments that a “hole” develops in the density distribution around the binary as it hardens – instead, the central density drops slowly as the binary hardens (see Figure 1). On the contrary, as explained in the previous section, the wandering of the BBH increases as the central density falls, allowing it to interact with new regions of stellar phase space, and to keep hardening at a rate consistent with $`1/at^d`$. The value of $`d1`$ is realized in the limit where the BBH makes up a small fraction of the host stellar system mass, as is typically the case for galactic bulges.

What if the mean wandering radius of the BBH were smaller than the semi-major axis at which the binary becomes hard ($`a_h`$)? Then the binary would wander “within itself”, so to speak; in this case one may wonder whether it would not encounter fresh regions of phase space, and perhaps hardening would stop (the bottleneck, if there is one, would occur at the stage at which the binary becomes hard, because wandering only becomes more important as the semi-major axis of the binary shrinks). Is there a sign of this happening in the simulations?

Unfortunately, to do this experiment in the regime in which we are most interested (i.e., for BBHs of low mass, so that the hardening rate is constant over time), one would require the particle number to be prohibitively high. However, it is possible to perform this experiment for more massive BBHs. For example, Table 2 lists the experimental root-mean-square value of the binary’s Brownian motion, $`r_{cm}`$, for a number of experiments, with the corresponding value of $`a_h`$; in each of the cases, $`r_{cm}<a_h`$. Hardening does not stop at or near this point, but continues indefinitely beyond this value at the rate $`1/at^d`$, with $`0.5\stackrel{<}{}d\stackrel{<}{}1.0`$. Figure 2 shows the hardening of a BBH with $`m_1=m_2=0.02`$, with 200,000 particles. The rms value of the binary’s Brownian motion is 0.013; the binary becomes hard when $`a=a_h=3m_2a/2M=0.018`$, or when $`1/a_h=56`$.

In each of these cases, however, we find that $`r_{cm}>a_{lc}`$; moreover, the value of $`r_{cm}`$ is very different from the value predicted by equation (13). This indicates that when the binary is allowed to wander, hardening no longer stalls at the semimajor axis $`a_{lc}`$; rather, the amplitude of its wandering increases as the central density of the star cluster declines (see Figure 3 for an explicit example of this), resulting in a lower restoring force holding the binary at the center. This, in turn, allows the binary to keep hardening by interacting with new stars that it would not have encountered if it did not wander. Thus, the “loss cone” is repopulated at a rate higher than the two-body relaxation rate by virtue of the increased wandering of the BBH’s center of mass.

## 5 Implications for the Coalescence of Massive Black Hole Binary Systems

We have evidence, therefore, that the BBH can continue to harden beyond the point at which the loss cone would be evacuated were the binary’s center of mass to be fixed at the origin, and that, moreover, it does so at roughly the rate calculated in §2. If we assume that hardening can continue until the point at which the emission of gravitational radiation becomes the dominant mode of energy loss for the binary, we can calculate when that should occur.

From equation (12), the hardening timescale of the BBH is

$$t_h(a)=\left|\frac{a}{\dot{a}}\right|=\frac{21\pi }{256\sqrt{2}}\left(\frac{r_0^5}{GM}\right)^{1/2}\frac{1}{a}.$$
(15)
The hardening timescale for gravitational radiation is (Peters 1964)

$$t_h(a)=\left|\frac{a}{\dot{a}}\right|=\frac{5}{64}\frac{c^5a^4}{G^3\mu _{12}(m_1+m_2)^2}\frac{(1e^2)^{7/2}}{1+73e^2/24+37e^4/96},$$
(16)
where $`c`$ is the velocity of light, and $`e`$ the eccentricity of the BBH’s orbit. The dependence of the gravitational radiation timescale on eccentricity is weak unless $`e`$ is close to 1. As noted in §3, $`e`$ is very small for heavy BBHs, but rises to larger values ($`0.2`$) for lighter BBHs. In the absence of a precise characterization of the variation of $`e`$ with BBH mass, we take $`e=0`$ in the above expression; thus, our result for the time of onset of gravitational wave emission and coalescence should be regarded as an upper limit, as a higher eccentricity would allow this stage to be reached more quickly.

Equating the two expressions, we have that the semi-major axis at which gravitational radiation losses start to dominate, $`a_{gr}`$, is given by

$$a_{gr}^5=\frac{21\pi }{20\sqrt{2}}\frac{G^{5/2}m_1^3q(1+q)r_0^{5/2}}{c^5M^{1/2}},$$
(17)
where $`q`$ is the mass ratio $`q=m_2/m_11`$.

From equation (12), assuming that $`a_{gr}a_h`$, we can derive the time between the binary’s becoming hard at semi-major axis $`a_h`$, and the time at which $`a=a_{gr}`$,

$$t_{gr}=\frac{21\pi }{256\sqrt{2}}\frac{r_0^{5/2}}{G^{1/2}M^{1/2}}\frac{1}{a_{gr}},$$
(18)
or

$$t_{gr}=\frac{21\pi }{256\sqrt{2}}\left(\frac{20\sqrt{2}}{21\pi }\right)^{1/5}\frac{cr_0^2}{GM^{2/5}m_1^{3/5}[q(1+q)]^{1/5}}.$$
(19)
The dependences on $`M,m_1`$ and $`q`$ are rather weak.

We can also calculate the total number of stars, $`N_{int}`$, the BBH must interact with from the time it becomes hard to the time $`t_{gr}`$. We have $`\mathrm{\Delta }E=G\mu _{12}m_{}/a=\frac{Gm_{}}{a}\frac{m_1q}{1+q}`$. Therefore $`dE=\frac{Gm_{}}{a}\frac{m_1q}{1+q}dN`$, or $`\frac{dE}{dN}=\frac{Gm_{}}{a}\frac{m_1q}{1+q}`$. But $`E=Gm_1m_2/2a`$, and so $`\frac{dE}{dN}=\frac{Gm_1^2q}{2a^2}\frac{da}{dN}`$. Equating the above, we get $`\frac{dN}{da}=\frac{m_1}{m_{}}\frac{1+q}{2a}`$. Integrating, the total number of stars the BBH interacts with is

$$N_{int}=\frac{m_1(1+q)}{2m_{}}\text{ln}\left(\frac{a_h}{a_{gr}}\right),$$
(20)
where $`a_h=3m_2r_0/2M=3m_1qr_0/2M`$. The total mass with which the BBH interacts is then simply

$$M_{int}=\frac{(m_1+m_2)}{2}\text{ln}\left(\frac{a_h}{a_{gr}}\right).$$
(21)
As will be seen from Figure 4c, the total mass in stars the BBH must interact with is roughly of the order of itself. Therefore, according to this theory, the BBH should not change the stellar distribution function much, as long as the mass of the BBH is not significant in comparison with the total mass in stars. This result is consistent with our simulations for small mass BBHs.

Some of the above relations are depicted for physically reasonable stellar systems and BBHs in Figure 4. Instead of presenting results in terms of $`M`$ and $`r_0`$, we use more physically and observationally relevant variables, such as the one-dimensional stellar velocity dispersion, $`\sigma `$, and the central density of the stellar system, $`\rho _0`$. The relations amongst these quantities for the Plummer model are as follows:

$$\rho _0=\frac{3M}{4\pi r_0^3},$$

$$\sigma ^2=\frac{GM}{6r_0}.$$

We then obtain

$$a_{gr}^5=\frac{63\sqrt{3}}{80}\frac{G^2m_1^3q(1+q)\sigma }{c^5\rho _0},$$
(22)
and

$$t_{gr}=\frac{63\sqrt{3}}{1024}\left(\frac{80}{63\sqrt{3}}\right)^{1/5}\frac{c\sigma ^{4/5}}{G^{7/5}m_1^{3/5}[q(1+q)]^{1/5}\rho _0^{4/5}}.$$
(23)
Figures 4a and 4b show that for physically reasonable masses of the BBH and characteristics of the host galaxy, the semi-major axis can shrink to the value $`a_{gr}`$ within a Hubble time, after which coalescence occurs quickly by the emission of gravitational waves. The actual values of $`a_{gr}`$ are very small (Figure 4d).

The above calculations assume that the binary hardens according to $`1/at`$, which is the behavior we see in the simulations when the mass of the BBH is less than about 1% of the total mass in stars. Is this a reliable guide to what happens in actual galaxies?

The ratio of the mass of the BH to the mass of the galaxy is observed to be quite low (on the order of $`10^3`$: see Merritt & Ferrarese 2001, Wandel 1999, Magorrian et al. 1998), and therefore the linear hardening behavior is most likely to apply in the cases of practical interest.

However, it is interesting to ask what would happen if the BBH hardened more slowly than at a constant rate. As an example, consider the case shown in Figure 2, with $`m_1=m_2=0.02`$, and $`N=200,000`$ stars. The hardening behavior is well fit by the expression $`1/a=21.7\times t^{0.5}`$. For purposes of illustration, let us use $`a_{gr_1}`$ to denote the semimajor axis at which gravitational wave radiation becomes dominant if we assume that the BBH hardens at a constant rate, $`1/a20.6\times t`$, and $`a_{gr_2}`$ the semimajor axis at which gravitational wave radiation becomes dominant if we assume that the BBH hardens at the above non-linear rate, $`1/a21.7\times t^{0.5}`$, as is actually observed in the simulations. If we scale our galaxy such that $`r_0=3\pi /16`$ corresponds to 100 pc, and $`M=1`$ corresponds to a total mass in stars in the bulge of $`10^9M_{}`$, then we obtain $`a_{gr_1}/a_{gr_2}0.3`$. Thus, the gravitational radiation stage sets in at a larger semimajor axis if hardening occurs more slowly than linearly. This arises from the different behaviors of the hardening timescales: when $`1/at`$, $`t_h=|a/\dot{a}|1/a`$, whereas when $`1/at^{0.5}`$, $`t_h1/a^2`$. However, the times taken to reach $`a_{gr_1}`$ and $`a_{gr_2}`$ are approximately 1 Gyr and 70 Gyr respectively, since in the latter case the BBH hardens only as $`t^{0.5}`$ rather than $`t`$, as in the former case. Therefore, despite the fact that the gravitational radiation stage is reached at a larger semimajor axis if hardening occurs slower than linearly, the time taken to reach this stage is longer in the former case. Thus, our results in this section hold for light BHs which, as mentioned before, is the case with most galaxies.

## 6 Summary

In summary, our primary results are:

* We have developed a consistent description for the behavior of the semi-major axis of a hard BBH. As long as the BBH mass is small compared to the total mass in stars: $`1/at^d`$, with $`d1`$.
* Our simplified theory gives the value of the proportionality constant of the above linear relation, $`s=20.6`$, independent of the mass of the BBH or the number of stars. We find, however, from numerical experiments that $`s`$ depends on the number of stars – as the number of stars increases, $`s`$ drops, but eventually converges to a particular value for $`N200,000`$ stars. A possible explanation is that this reflects an equilibrium between two competing processes: as the number of stars increases, the wandering of the binary becomes smaller, causing a larger drop in central density; this, however, increases the wandering of the BBH because of the reduction of the central restoring force on it.
* We also find that $`s`$ depends on the mass of the BBH: in our simulations for equal-mass BHs, the convergent value of $`s`$ increases from $`8`$ to $`15`$ as the mass of each BH drops from 0.005 to 0.00125. We interpret this as meaning that the theoretical value 20.6 is reached when no change in the stellar distribution function is caused by the BBH – this will be true only when the mass of the BBH is very small compared to the total mass in stars. If we could run simulations with even lower mass BBHs, we may approach 20.6. However, it should be remembered that we have assumed a somewhat uncertain form for the energy removed from the binary by each stellar interaction. We also note that for calculating the time at which gravitational radiation starts to dominate the energy loss, we do not need to know the precise value of $`s`$, since the fifth root of this number needs to be taken, as shown in §5. There are two difficulties with running simulations with very light BBHs: first, the black holes need to be much more massive than the stars; second, the dynamical friction timescale for light black holes, and consequently the time for them to become hard, is very long.
* For light BBHs, strong disruption of the stellar density does not occur: the central density of the cluster does decline with time but gently. This is where we expect our theory to be applicable. There is no sign that a “hole” is forming at the center, and thus no sign that hardening would stop after the loss cone is evacuated. Wandering allows the BBH to interact with regions of stellar phase space it would not otherwise encounter. According to our experiments, if we keep the BBH center of mass artificially clamped to the center, hardening quickly stops after the BBH has interacted with the stars that could pass within roughly a distance $`a`$ from it.
* For heavy BBHs, there is stronger disruption of the stellar density and distribution function, and strong deviations from the linear behavior of the inverse semi-major axis: $`1/at^d`$, with $`0.5d1`$. This is the regime in which our theory is not applicable, though it gives a qualitative explanation for why $`d`$ should be bracketed by 0.5 and 1.
* We have shown that the hardening of the BBH does not stop even when its typical wandering radius is smaller than the value of the semi-major axis at which it becomes hard, at least in the case when the BBH is fairly massive; however, the wandering radius is always larger than the semi-major axis at which hardening would halt because of loss cone depletion if the binary did not wander. The reason involves the increased wandering of the binary’s center of mass as the central density of the stellar system falls, resulting in a reduction of the restoring force keeping the binary at the center. To verify this in the case of light binaries would require simulations with many millions of particles.
* We have calculated an upper limit on the time required for the BBHs to reach the stage at which energy losses are dominated by gravitational radiation emission. For reasonable characteristics of the stellar system, the gravitational radiation stage is reached within a Hubble time. In order to reach this stage, the BBH must interact with a total mass in stars that is of order its own mass.

Our model uses the spherically symmetric Plummer density profile, which has a flat density structure at the center. However, observations have shown that most elliptical galaxies are characterized by weak or strong cusps (e.g., Faber et al. 1997, Byun et al. 1996, Lauer et al. 1995). Is it then reasonable to use the Plummer model in calculations of BBH coalescence?

The Plummer model is characterized by a simple distribution function, which makes it possible to obtain analytic results for the phenomena discussed in this paper. Here, we are mainly interested in understanding the effect of the binary’s Brownian motion on its hardening, rather than in calculating precisely the time taken for a binary in a real galaxy to coalesce. Despite the fact that our model for galactic bulges is unrealistic, we believe that our results capture generic dynamical processes that take place in more complicated situations for a number of reasons. First, recent simulations have shown that the merger of two galaxies with steep density cusps produce a galaxy with a shallow density cusp (Milosavljević & Merritt 2001), owing to the transfer of energy from the binary to the stars; moreover, the destruction of the steep cusp takes place very rapidly (on a timescale comparable to the local dynamical time), before the binary becomes significantly hard (similar results have been obtained by Nakano and Makino 1999 a, b). Our results apply to the stage when the binary has become hard. Secondly, the weak cusp of the merged galaxy is carried by the BBH as it wanders around; whereas the stellar density and velocity dispersion in our calculations refer to their values in the core of the galaxy but outside the Keplerian rise around the black holes, since most of the encounters that remove energy from the BBH are with stars from outside the cusp. Thirdly, the presence of the density cusp would only increase the supply of stars to the binary, resulting in more efficient hardening; if a binary coalesces within a Hubble time in our model, it would likely coalesce within a shorter time in the presence of a cusp.

We thank G. Quinlan for providing the simulation code. For useful comments, we thank G. Rybicki and the participants of the Massive Black Hole Coalescence Focus Session at Pennsylvania State University, November 2002. AL acknowledges support from the Institute for Advanced Study at Princeton, the John Simon Guggenheim Memorial Fellowship, and NSF grants AST-0071019, AST-0204514.