text:
# Global Optimization by Energy Landscape Paving

## Abstract

We introduce a novel heuristic global optimization method, energy landscape paving (ELP), which combines core ideas from energy surface deformation and tabu search. In appropriate limits, ELP reduces to existing techniques. The approach is very general and flexible and is illustrated here on two protein folding problems. For these examples, the technique gives faster convergence to the global minimum than previous approaches.

Global optimization is one of the key issues in modern science, technology and economy. Typical examples are the problem of optimal transportation routes , finding molecular conformations  or fitting experimental spectra. Consequently, much effort has been spent on designing methods to finding global optima. For this purpose, the system has to be described by an objective function, and optimality is achieved when this function reaches its global minimum. If the objective function is viewed as an ‘energy’ the optimal solution corresponds to the deepest minimum in the energy landscape. For most applications of practical interest, competing interactions and frustration in the system lead to an energy landscape with many local minima separated by high barriers. Since conventional minimization techniques tend to get trapped in whichever local minimum they encounter first it turns out to be extremely difficult to find the global minimum in such cases.

A general characteristic of successful optimization techniques is that they avoid entrapment in local minima and continue to explore the energy landscape for further solutions. For instance, in tabu search  the search is guided away from areas that have already been explored in an effort to cover all important regions of the solution space. The danger with such an approach is that it may result in slow convergence since it does not distinguish between important and less important regions of the landscape.

Entrapment in local minima can also be avoided if the search is performed in a deformed or smoothed energy landscape, for example by lowering diffusion barriers, in stochastic tunneling  or the various generalized ensemble approaches.  In the optimal case the original energy landscape is transformed in a funnel-landscape and convergence toward the global minimum is fast. Although they have been very successful, most of these methods require a considerable amount of fine-tuning or a priori information. Moreover, problems may exist when connecting back to the original landscape since minima on the deformed surface may have been displaced or merged.

Here we introduce a novel approach to the global optimization problem that combines ideas from tabu search and energy landscape deformation. The new method, energy landscape paving (ELP), avoids some of the pitfalls of the other two and has very general applicability. The central idea is to perform low-temperature Monte Carlo (MC) simulations, but with a modified energy expression designed to steer the search away from regions that have already been explored. To be specific, we choose as the statistical weight for a state

$$w(\stackrel{~}{E})=e^{\stackrel{~}{E}/k_BT},$$
(1)
where $`T`$ is a (low) temperature and $`\stackrel{~}{E}`$ the following replacement of the energy $`E`$:

$$E\stackrel{~}{E}=E+f(H(q,t)).$$
(2)
In this expression, $`f(H(q,t))`$ is a function of the histogram $`H(q,t)`$ in a pre-chosen “order parameter” $`q`$. The histogram is updated at each MC step, hence the “time” dependence of $`H(q,t)`$. As a result, the search process keeps track of the number of prior explorations of a particular region in order parameter space and biases against revisiting the same types of states. Rather than using the system states themselves in the histograms an appropriate order parameter is employed. This may be a “natural” quantity for the system under study (such as a spheroidal deformation for a cluster) or the energy itself may be taken as the order parameter.

In a regular low-temperature simulation the probability to escape a local minimum depends only on the height of the surrounding energy barriers. Within ELP the weight of a local minimum state decreases with the time the system stays in that minimum, and consequently the probability to escape the minimum increases. Hence, ELP utilizes the interplay of two factors. Given equal frequencies $`H(q,t)`$, the simulation will favor low energies, thus insuring that no unphysical high-energy conformations are sampled. However, soon the system will run into a local minimum. With time, ELP deforms the energy landscape locally in such way that the local minimum is no longer favored and the system will explore higher energies. It will then either fall in a new local minimum or walk through this high energy region till the corresponding histogram entries all have similar frequencies. At that point the original energy landscape is restored (that is, only shifted by a constant (and irrelevant) factor), and the system again has a bias toward low energies.

ELP bears some similarities to tabu search  in that recently visited regions are not likely to be revisited immediately. Revisitation moves are not completely forbidden, but are given an exponentially lower weight compared to moves that go to regions with a comparable energy that have been explored less. With a short-term memory in the histogram and infinite cost for ‘forbidden’ moves ELP becomes completely equivalent to tabu search. On the other hand, ELP is also akin to an energy deformation approach  in that the additional histogram parameter may be viewed as a (continuously changing) deformation of the energy landscape depending on the frequency with which a particular area (as characterized by its order parameter) has already been explored. Obviously for $`f(H(q,t))=f(H(q))`$ the method reduces to the various generalized-ensemble methods.

We have tested ELP in the context of protein folding, which involves the prediction of the biologically active conformation of a protein solely from the sequence of amino acids. Assuming that this structure is thermodynamically stable, it is reasonable to identify the global-minimum conformation in the free energy at $`T300`$ K with the lowest potential energy conformation and to choose the potential energy of the protein as an objective function. The complexity and importance of the problem make it an ideal target for a test of our new optimization technique.

As with any optimization method, ELP requires the choice of an energy function by which the multitude of protein configurations can be discriminated. Here, we choose the ECEPP/2 force field, a commonly used energy function in protein simulations, as implemented in the computer code SMMP.

In order to test and illustrate ELP, we concentrated on the structure prediction of two molecules. The first system is the pentapeptide Met-enkephalin, which has become a frequently used benchmark model to examine new algorithms. We know from previous work that the ground state of this peptide with the ECEPP/2 force field is given by $`E_0=10.7kcal/mol`$, and that the next higher local minimum has an energy of $`E_1=9.8kcal/mol`$. Hence, we identify any configuration with energy below $`E=9.8kcal/mol`$ as a representative of the ground state.

In the ELP simulations of Met-enkephalin, we used the potential energy itself as an order parameter and thus the deformed energy landscape is generated by $`\stackrel{~}{E}=E+H(E,t)`$, where $`H(E,t)`$ is the histogram in energy at MC sweep $`t`$. We chose a bin size $`E_{bin}=0.25kcal/mol`$ in the histogram, but checked our results also for $`E_{bin}=0.5kcal/mol`$ and $`E_{bin}=1kcal/mol`$ without finding noticeable differences in our results. Setting the temperature to $`T=50K`$, and with $`\beta =1/k_BT`$ we find as a weight for the MC simulation:

$$w(E,t)=e^{\beta (E+H(E,T))}.$$
(3)
The characteristic behavior of our ELP method is exemplified in Fig. 1 which shows for Met-enkephalin the time series of a simulation with 50,000 sweeps. The starting configuration has an energy of $`E_{start}=5.1kcal/mol`$ and was obtained from a random configuration through quenching in initial 100 sweeps. The simulation soon gets trapped in a local minimum of $`E7.8kcal/mol`$ (after only 250 MC sweeps). Through the following MC sweeps entries in the corresponding histogram bin are accumulated and the energy landscape locally deformed, until after about 750 MC sweeps the simulation escapes this local minimum to find a lower local minimum after 2000 MC sweeps. This process is repeated till the simulation finds the global minimum conformation for the first time after 7260 sweeps. Within the 50,000 sweeps of our simulation the ground state region ($`E<9.8kcal/mol`$) was visited 5 times, each visit separated by explorations in the high energy region.

Note that the range of energies covered increases with MC time: ELP starts with filling up the small ‘potholes’ in the energy landscape, but later in the simulation large valleys are also filled up. Hence, our algorithm is self-adjusting: with increasing length of the simulation it becomes possible to overcome higher and higher energy barriers. In that regard, ELP is more efficient than standard techniques such as simulated annealing  where the height of energy barriers that can be overcome shrinks with MC time. In order to test that conjecture we performed 20 independent runs of 50,000 MC sweeps with ELP and compared the results with 20 simulated annealing runs of equal statistics. In the simulated annealing runs the temperature was lowered exponentially in 50,000 sweeps from an initial temperature $`T=1000K`$ to a final temperature $`T=50K`$. In Ref.  this proved to be the optimal annealing schedule for Met-enkephalin. Even with this optimized annealing schedule, the ground state was found only in $`8/20=40\%`$ of the runs (in an average time of 43,000 MC sweeps) and the average value of the lowest energy conformation ($`<E_{min}>=8.5kcal/mol`$) was above our threshold for ground state configurations ($`9.8kcal/mol`$) and subject to large fluctuations (with standard deviation $`\sigma =2.1kcal/mol`$). Better results were obtained in 20 tabu search runs with same statistics. Here the ground state region was found in $`10/20=50\%`$ of the runs and the average value of the lowest energy conformations was $`<E_{min}>=9.5kcal/mol`$. On the other hand, with ELP we found the ground state in each of the 20 runs (in, on average, 12,700 MC sweeps) and the average of lowest energy states $`<E_{min}>=10.3kcal/mol`$ was well below our threshold for ground state configuration, subject to only small fluctuations ($`\sigma =0.3kcal/mol`$).

In order to further evaluate the ELP approach we also studied the much larger villin headpiece subdomain, 36-residue peptide (HP-36) that is with 597 atoms about 8 times larger than Met-enkephalin (75 atoms). HP-36 is one of the smallest peptides that can fold autonomously and was chosen recently for a 1-microsecond molecular dynamics simulation of protein folding. The experimental structure was determined by NMR analyses. Since it is a solvated molecule we also had to take into account the interaction between protein and solvent. We have approximated this contribution to the overall energy by adding a solvent accessible surface term  to the energy function: $`E=E_{Ecepp/2}+_i\sigma _iA_i.`$ Here, the sum goes over all atoms and the $`A_i`$ are the solvent accessible surface areas of the atoms. The parameters $`\sigma _i`$ were chosen from Ref. .

HP-36 allows in a simple way the definition of an order parameter to characterize configurations other than by their energy. This natural order parameter is the number $`n_H`$ of residues in the peptide which are part of an $`\alpha `$helix. Following earlier work  we define a residue as helical if the pair of backbone dihedral angles $`\varphi ,\psi `$ takes a value in the range $`(70\pm 20,37\pm 20)`$. Throughout the search process we tried now to deform the energy landscape by means of a histogram $`H(E,n_H,t)`$ in both helicity and energy: $`\stackrel{~}{E}=E+H(E,n_H,t)`$. Operating again at a temperature $`T=50`$ K, we find as weights for the search algorithm

$$w(E,n_H,t)=e^{\beta (E+H(E,n_H,t))}.$$
(4)
Using this weight we performed simulations with 50,000 MC sweeps (starting from random configurations) keeping track of the lowest energy configuration during the search process.

The structure of HP-36 as obtained from the Protein Data Bank (PDB code 1vii) is shown in Fig. 2. The structure consists of three helices between residues 4-8, 15-18, and 23-32, respectively, which are connected by a loop and a turn. After regularizing this structure with the program FANTOM  we obtained as its energy (ECEPP/2 + solvation term) $`E_{nat}=276`$ kcal/mol. Our new ELP method led after 25,712 MC sweeps to a configuration with lowest energy $`E_{min}=277`$ kcal/mol which we show in Fig. 3. The above structure has a radius of gyration $`R_\gamma =10.1`$ Å indicating that the numerically obtained structure is slightly less compact than the experimental structure ($`R_\gamma =9.6`$Å). It consists of three helices where the first helix stretches from residue 2 to residue 11 and is more elongated than the corresponding one in the native structure (residues 4-8). The second helix consist of residues 13-17 (compared to residue 15-18 in the native structure) and the third helix stretches from residue 23-33 (residues 23-32 in the PDB structure). The structure has 95% of the native helical content, that is 95% of all residues which are part of a helix in the experimental structure are also part of a helix in our structure. We also note that 65% of the native contacts were formed in our structure (two residues $`i`$ and $`j`$ ($`j>i+2`$) are taken to be in contact if their $`C_\alpha `$ atoms are closer than 8.5 Å). Both values are comparable with the results in Ref.  (but required orders of magnitude less computer time) where the optimal structure of a $`1\mu s`$ molecular dynamic folding simulation showed 80% of native helical content and 62 % of native contacts. Similarly comparable were the values of the root-mean-square deviation (RMSD) of both numerically determined conformers to the native structure: $`5.8`$ Å versus $`5.7`$ Å in Ref.  when all backbone atoms where counted.

We conclude that even for large peptides such as HP-36 our novel optimization method is able to find structures that are close to the experimentally determined ones. In passing, we remark that an exploratory simulated annealing run of 100,000 sweeps did not lead to such structures. However, our ELP prediction of the HP-36 structure is limited to an RMSD of $`6`$ Å. This points to a general problem in protein simulations: it is not clear whether the utilized cost function has indeed the biologically active structure of a given protein as its global minimum. In fact, our optimal structure has slightly lower energy than the native one. The problem becomes obvious when solvation effects are neglected. An ELP run of 50,000 sweeps relying only on the ECEPP/2 force field led to a lowest-energy structure with an ECEPP energy of $`E_{GP}=192`$ kcal/mol (found after 25,712 sweeps). That structure, build out of two helices (between residues 2-16 and 23-33) connected by a loop, differs significantly from the regularized PDB-structure with the higher potential energy $`E_{nat}=176`$ kcal/mol. Hence, the native structure of the peptide HP-36 is not the global minimum configuration in ECEPP/2. Only the inclusion of the solvation term led to an essentially correct structure as global minimum configuration.

Summarizing, we have developed a new and general stochastic global optimization method that is easy to implement and combines energy landscape deformation ideas with elements of tabu search. The efficiency of ELP was compared with simulated annealing  and tabu search.  To illustrate the power of our novel approach, we applied it to the structure prediction of HP-36, a 36 residue peptide. For this large peptide an unbiased all-atom simulation using ELP led to a 3D structure very close to the experimentally determined one. In future work we want to extend application of our new approach to other optimization problems.

Acknowledgement
U.H.E. Hansmann acknowledges support by a research grant of the National Science Foundation (CHE-9981874) and by NSF MRI grant #9871133. Part of this article was written while U.H. was visitor at the Department of Physics at the Bielefeld University. He thanks F. Karsch for his hospitality during his stay in Bielefeld.

Figures:

1. Time series of a minimization run of 50,000 sweeps for the pentapeptide Met-enkephalin.
2. NMR derived structure of the 36 residue peptide HP-36 as obtained from the PDB data base (1vii).
3. Lowest-energy structure of HP-36 as obtained with ELP.


format: text

text:
# GRB010921: Localization and Observations by the HETE Satellite

## 1 Introduction

As has long been recognized, accurate locations and rapid follow-up observations in many wavelengths are central to understanding the nature of gamma-ray bursts. For this reason, a strategy evolved in the late 1970’s and early 1980’s (e.g., Woosley et al. 1982) to detect GRBs, not only in gamma-rays, but also in emitted X-rays and optical light. While detection of short transients at the lower energies posed observational challenges, and the strength of the optical signal was unknown, the possibility of arc minute localizations deduced from the X-rays that were known to be present was very appealing. This strategy was implemented in the HETE-1 (High Energy Transient Explorer) satellite, which was unfortunately lost due to a rocket failure on 1996 November 4, and in the highly successful BeppoSAX Mission (Costa et al. 1997).

The HETE-2 satellite (henceforth simply “HETE”), which was successfully launched into equatorial orbit on 9 October 2000, is the first space mission entirely devoted to the study of gamma-ray bursts (GRBs). HETE utilizes a matched suite of low energy X-ray, medium energy X-ray, and gamma-ray detectors mounted on a compact spacecraft. A unique feature of HETE is its capability for localizing GRBs with $``$1-10′ accuracy in real time aboard the spacecraft. GRB locations are transmitted, within seconds to minutes, directly to a dedicated network of telemetry receivers at 13 automated “Burst Alert Stations” (BAS) sited along the satellite ground track (Villasenor et al 2002). The BAS network then re-distributes the GRB locations world-wide to all interested observers via Internet and the GRB Coordinates Network (GCN) in $``$1 s (Vanderspek et al. 2002, Barthelmy et al. 2002). Thus, prompt optical, IR, and radio follow-up identifications can be anticipated for a large fraction of HETE GRBs.

Here we report the localization of the first HETE-discovered 111HETE detects $``$50 GRBs yr-1, of which $``$15 yr-1 are localized by the WXM (Ricker et al. 2002) GRB for which a counterpart has been found. Based on the combined data from HETE and the Interplanetary Network (IPN), a diamond-shaped error box roughly 15 arc minutes on an edge was established for GRB010921. Subsequent searches of this error box by ground-based optical and radio instruments revealed a fading counterpart whose properties are reported elsewhere (Price et al 2001a, 2002; Kulkarni et al 2002; Park et al 2001a, 2002; Djorgovski et al. 2001).

## 2 HETE Mission and Instruments

The HETE spacecraft is a small satellite, measuring roughly a meter high by half a meter in diameter, with a mass of 125 kg. Constructed and launched at a cost less than $`\frac{1}{3}`$ that of a NASA Small Explorer (SMEX), HETE was primarily developed and fabricated in-house at MIT by a small scientific and engineering team, with major hardware and software contributions from international partners in France and Japan (Doty et al. 2002). Contributions to software development were also made by scientific partners in the US, at the Los Alamos National Laboratory, the University of Chicago, and UC-Berkeley. Operation of the HETE satellite and its science instruments, along with a dedicated tracking and data telemetry network, is carried out by the HETE Science Team itself (Crew et al. 2002). Routine activities are largely automated, with oversight and monitoring by members of the Science Team at MIT and at RIKEN.

The gamma-ray burst detection and localization system on HETE consists of three complementary instruments: the French Gamma-ray Telescope (FREGATE), Wide-Field X-Ray Monitor (WXM), and the Soft X-Ray Camera (SXC). The manner in which the three HETE science instruments operate cooperatively is described in Ricker et al. (2002). Since GRB010921 was outside the field-of-view of the SXC, that instrument will not be discussed in this paper.

French Gamma-ray Telescope (FREGATE). FREGATE consists of four, co-aligned cleaved NaI(Tl) scintillators. Its prime objectives are the detection and spectroscopy of GRBs and the monitoring of variable X-ray sources. FREGATE is optimally sensitive in the 6-400 keV energy band, and provides limited ($``$ $`\pi `$ sr) localization information. The FREGATE instrument was developed by Centre d’Etude Spatiale des Rayonnements (CESR; Toulouse, France). Additional details regarding FREGATE are given in Table 4, and in Atteia et al (2002).

Wide-Field X-Ray Monitor (WXM). The Wide-Field X-ray monitor (WXM) consists of two, crossed one-dimensional coded aperture cameras. The cameras are designated WXM-X and WXM-Y, respectively, indicating the spacecraft axis along which each is optimized for GRB localization. Two position sensitive proportional counters (PSPC) are utilized in each camera. The prime objectives of the WXM are the detection and spectroscopy of GRBs and the monitoring of variable X-ray sources in the low energy band. The WXM is optimally sensitive in the 2-25 keV energy band, and typically provides $``$10′ burst localizations. The WXM was developed by RIKEN(Japan) and the Los Alamos National Laboratory. Additional details regarding the WXM are given in Table 4, and in Kawai et al (2002) and Shirasaki et al. (2000).

## 3 Observations

Standard HETE operations consist of a survey for transient gamma-ray and X-ray events in the 2-3 steradians centered on the antisolar point. Observations take place in that half of the HETE orbit where the instruments’ FOV is clear of the Earth. On September 21 at 18950.56 SOD (05:15:50.56) UT, the FREGATE instrument on HETE detected a bright ($`>80\sigma `$) GRB. In this section we describe the localization and properties of GRB010921 as determined from HETE and IPN observations of this burst.

### 3.1 Localization

GRB010921 was recorded as trigger event H1761 by the FREGATE instrument on HETE, and was promptly reported as a GCN Notice at 05:16:08 UT, approximately 17 seconds after the burst. The burst was also detected at a SNR=7.1 in the WXM X detector. The projection of the burst direction was offset 24 degrees from the WXM X-Z plane, resulting in a good ($`\pm `$10′) localization in the X direction and crude limits on Y, thereby localizing the source to a thin, long strip. The projection of the burst direction was offset 39 degrees from the WXM Y-Z plane, outside the fully-coded FOV. Thus, the WXM Y camera provided no useful information regarding the Y location of the burst. Analysis of the Y location required the rapid development and testing by the HETE Team of non-standard ground analysis software that would rely on the limited WXM X camera data, delaying the dissemination of the two-dimensional localization reported in GCN Circular 1096 until approximately 5.1 hours after the GRB. The initial estimate reported in GCN Circular 1096 was that the localization strip extended $``$10o in the long direction, and 20$`\mathrm{}`$ in the short direction (Ricker et al 2001a).

Further refinement of the strip dimensions reduced the WXM error region to a 5.2o by 17$`\mathrm{}`$ box, with the 4 corners at coordinates:
$`\alpha _{J2000}`$ = 22h55m34.7s, $`\delta _{J2000}`$ = 40o25′55″
$`\alpha _{J2000}`$ = 22h54m05.7s, $`\delta _{J2000}`$ = 40o26′42″
$`\alpha _{J2000}`$ = 23h04m30.4s, $`\delta _{J2000}`$ = 45o23′35″
$`\alpha _{J2000}`$ = 23h02m60.0s, $`\delta _{J2000}`$ = 45o24′18″

Both Ulysses and the gamma-ray burst monitor on BeppoSAX observed GRB010921, and an initial IPN error box was circulated 15.2 h after the burst (Hurley et al. 2001). Triangulation has now been performed using the final data, resulting in a Ulysses to HETE-FREGATE annulus centered at $`\alpha _{J2000}`$, $`\delta _{J2000}`$ = 15h29m23.6s, 67o36′14″, with radius 60.082 $`\pm `$ 0.118 o (3 $`\sigma `$). This annulus intersects the 3 $`\sigma `$ WXM error box to give a 310 square arcminute error region whose coordinates are given in Table 1; the diamond-shaped error region is shown in figure 1. This error box is larger than reported in Hurley et al. (2001) due to the fact that the initial WXM error box was quoted for 90% confidence (ie 2$`\sigma `$). Both the Ricker et al. (2001a) 2$`\sigma `$ error region and the 3$`\sigma `$ error region given here are consistent with the position of the optical transient discovered by Price et al. (2001a), with the OT lying 3.6′ from the centerline of the WXM error strip, and 3.1′ from the centerline of the IPN error annulus.

### 3.2 Time History, Peak Flux and Fluence

The time history of GRB010921 observed by the HETE science instruments is shown in figure 2. The six panels in figure 2 illustrate the energy dependence of the burst profile and duration. For FREGATE, all 4 detectors were summed; for the WXM, only the signal from one illuminated anode (out of six total) in the X-camera was utilized, so as to maximize the signal-to-noise ratio, since GRB010921 was quite far (43 degrees) off-axis from the instrument boresight. Quantitative measures of the burst temporal properties (i.e., $`t_{50}`$ and $`t_{90}`$; Paciesas et al. 1999) are given in Table 2. At the highest energies measured by FREGATE (32-400 keV band), the GRB010921 light curve exhibits a high degree of symmetry, with $`t_{50}`$=7.4s and $`t_{90}`$=21.8s. At lower energies, the $`t_{50}`$ duration becomes somewhat larger, increasing by $`30\%`$; on the other hand, $`t_{90}`$ is largely unchanged over the 4-400 keV band. Overall, the light curve for GRB010921 exhibits a hard-to-soft spectral evolution that is common in many GRBs.

The prompt emission by GRB010921 was measured over more than 2 decades in energy by HETE (Table 3). It is notable that GRB010921’s peak energy flux of 2.0 $`10^7`$ erg cm-2 s-1) in the 4-25 keV X-ray band was $``$ $`\frac{1}{3}`$ the peak energy flux in the “traditional” 50-300 keV gamma-ray band. A comparison of the peak photon flux in the two bands is even more striking: $``$4 times as many 4-25 keV photons were emitted as were 50-300 keV photons. (For GRB010921, the smaller SNR in the WXM compared to the SNR in FREGATE arises from the fact that GRB010921 was offset 43o from the boresight of the WXM.)

## 4 Discussion

The temporal characteristics of GRB010921 clearly mark it as belonging to the class of long duration GRBs (Hurley 1992, Kouveliotou et al. 1993). It is also quite “X-ray rich,” in that it exhibits a value of $`L_{x(425keV)}`$/$`L_{\gamma (50300keV)}`$ $``$0.3. However, it is by no means deficient in its 50-300 keV gamma-ray flux: in fact, its peak flux of 3.6 ph cm-2 s-1 in the 50-300 keV range would place it within the top decile of bursts in the 4Br Catalog (Paciesias et al 1999). Thus, GRB010921 does not qualify as an X-ray Flash (XRF; Heise et al. 2001). Based on its apparent redshift of z= 0.450 (Djorgovski et al. 2001), the implied isotropic emitted energy of GRB010921 in the 8-400 keV band is 7.8 1051 ergs (assuming $`\mathrm{\Omega }_M`$=0.3, $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.7, $`H_o`$=65 km s-1 Mpc-1; for beamed emission, the total energy should be multiplied by the as-yet-unknown beaming fraction.) The redshift of GRB010921, the second lowest among established GRBs, makes it a strong candidate for extended searches for a possible associated supernova. Furthermore, GRB010921 was, like all HETE-discovered GRBs, transiting near midnight at the initial epoch of the burst, resulting in its being well-situated for followup studies by ground-based telescopes for several months following the burst. Thus, additional insight into the nature of GRB010921 can be anticipated from long-term optical photometry by large telescopes.

## 5 Conclusions

The localization of GRB010921 by HETE and the IPN has led directly to a successful identification of a low-redshift afterglow counterpart. The prompt initial GCN alert for GRB010921 was distributed worldwide within 17 s of the burst; however, the determination of the HETE WXM localization was delayed by 5.1 hours due to the need to prototype new software to accommodate GRB010921’s extreme off-axis position in one of the WXM cameras. A confirmation, and definitive refinement of the WXM localization, using IPN data available 15.2 hours after the GRB enabled ground-based observers to successfully target GRB010921 during the first night following the burst, and to establish a counterpart well within the HETE-IPN error region. Hopefully, over the coming months the identification of optical counterparts resulting from prompt ($``$100 s), accurate HETE localizations will become routine. With the currently-projected long orbital lifetime ($`>`$10 years) and excellent health of the HETE spacecraft and instruments, we look forward to providing a uniquely valuable service to the worldwide community of GRB observers.

## Acknowledgments

The combined efforts over almost two decades by many talented and dedicated individuals made the HETE mission possible. In particular, the efforts of Bob Dill, Tye Brady, Dave Breslau, Gus Comeyne, Francis Cotin, Jim Francis, Greg Huffman, Frank LaRosa, Fred Miller, Jerry Roberts, and Fuyuki Tokanai were particularly notable.

The HETE mission is supported in the USA by NASA Contract NASW-4690; in Japan, in part by the Ministry of Education, Culture, Sports, Science, and Technology Grant-in-Aid 13440063. KH is grateful for Ulysses support under Contract JPL 958059, and for HETE support under Contract MIT-SC-R-293291. GP acknowledges support by the Italian Space Agency (ASI).


text:
# Effect of disorder on the conductance of a Cu atomic point contact

## I Introduction

Quantum point contacts (QPC’s) are constrictions in a conducting material whose size is of the order of the Fermi wavelength of the electrons. The first quantum point contacts were fabricated in a 2-dimensional electron gas (2DEG) in a semiconductor heterostructure PRL\_vanwees88 ; JPC\_Wharam88 . Their conductance exhibits a series of steps as a function of constriction size. The plateaus in between these steps are positioned at integer multiples of $`G_0=2e^2/h`$ $``$ (12.9 k$`\mathrm{\Omega }`$)-1, the quantum of conductance. This phenomenon is explained as a consequence of the quantization of the electronic energy levels, resulting from the confinement of the electrons in the direction perpendicular to the current. Signs of quantized behavior were also seen in the conductance of quantum point contacts made from metals PRB\_Krans93 ; PRL\_Pascual93 ; PRL\_Olesen94 . When a thin metallic wire is slowly broken, its conductance as a function of contact size shows smooth plateaus, alternated by steps of the order of $`G_0`$. In this case, however, the plateaus are not necessarily horizontal or positioned at integer multiples of $`G_0`$. Moreover, when the two wire ends are pressed together and broken again, the precise shape of the trace differs from trace to trace. In order to obtain fully reproducible results, a histogram is often created by projecting each conductance trace onto the conductance axis and adding the contribution of several thousands of traces. For many metals, the histograms show well-defined peaks which represent preferred conductance values.

The physical origin of these peaks has been debated heavily PRL\_Yanson97 , especially the question whether they are indeed related to the quantization of electronic energy levels or rather to the discrete and abrupt changes in the cross-sectional area of the nanowires, which are caused by the fact that the contact must change its size by at least one, discrete, atom. It was shown both experimentally, by force measurements PRL\_Rubio96 , and theoretically, by numerical simulations PRL\_Sutton93 , that steps are correlated with atomic rearrangements of the contact. Thus, plateaus in the conductance traces correspond to stable atomic configurations. Apart from its geometry, the conductance of such a configuration is determined by the conductance of individual atoms of the metal under investigation. It was shown Nature\_Scheer98  that the conductance of a one-atom point contact is governed by the nature of its valence orbitals. For monovalent s-metals such as Na, Cu or Au, a one-atom contact has a conductance of nearly 1 $`G_0`$. On the other hand, transition metals such as Ni, Pd or Pt, which have d-valence orbitals, have a one-atom conductance of the order of 1.6 $`G_0`$.

In experiments on both semiconductor QPC’s PRL\_vanwees88  and QPC’s made of monovalent s-metals it has been observed that the quantization of the conductance is never perfect: the peaks in the histograms are shifted below quantized values $`mG_0`$ ($`m`$ = 1, 2, 3 $`\mathrm{}`$). In metallic point contacts, such a shift has been observed in gold PRB\_Costa-Kramer97 ; PRB\_Hansen97 , sodium Nature\_Krans95 , silver PRB\_Hansen97  and copper PRB\_Hansen97 . Two possible origins of the shift have been discussed in the literature. It has been attributed to disorder either in the leads to the contact PRL\_vanwees88 , or in the contact itself PRL\_Garcia97 ; PRB\_Brandbyge97 ; PRB\_Burki99 . In all cases, a model was adopted in which the shift is accounted for by an ‘extra’ resistance $`R_s`$ that is placed in series with the true contact resistance $`R_c`$. $`G_c1/R_c`$ is assumed to be an exact integer multiple of $`G_0`$, and it is then straightforward to determine $`R_s`$ from the histograms. Typical experimental values are 400 $`\mathrm{\Omega }`$ for a 2DEG PRL\_vanwees88 , 100 $`\mathrm{\Omega }`$ for sodium Nature\_Krans95 , 100-500 $`\mathrm{\Omega }`$ for gold PRB\_Costa-Kramer97 ; PRB\_Hansen97 , 500 $`\mathrm{\Omega }`$ for silver PRB\_Hansen97  and 700 $`\mathrm{\Omega }`$ for copper PRB\_Hansen97 .

In all theoretical studies into disorder in QPC’s, adding disorder resulted in a downshift and smearing of the peaks in the conductance histogram. In the work of García-Mochales et al. PRL\_GarciaMochales97 , for example, the nanowire was simulated by a tight binding model. The disorder is introduced by random on-site energies, while the overlap between atomic orbitals is described by a hopping term between sites. The conductance was calculated as a function of the Fermi energy, while the geometry was kept constant. A histogram was then constructed, which showed peaks that were indeed shifted. The series resistance (calculated by fitting the best resistance value centering simultaneously all peak maxima at the corresponding integer values of $`G_0`$) were found to range from about 10 $`\mathrm{\Omega }`$ for small disorder to 1000 $`\mathrm{\Omega }`$ for large disorder.

Bürki et al. PRB\_Burki99  used a 2D-free electron model, with a potential term in the Hamiltonian describing randomly distributed impurities. The neck of the constriction in the nanowire was elongated, while keeping the Fermi energy constant. The histogram exhibited downshifted peaks that could be positioned at integer multiples of $`G_0`$ by subtracting a series resistance of 525 $`\mathrm{\Omega }`$.

Thus, while several authors have studied disorder in QPC’s from a theoretical point of view, very little experimental study was devoted to it. Ludoph et al. bas  studied conductance fluctuations as a function of bias voltage in atomic-size contacts for pure metals. A consistent explanation was obtained in terms of scattering of the electrons on defects near the contact and interference of various electron paths with the reflected wave at the contact itself. The values obtained for the elastic scattering length, $`l_e`$, were used to estimate an effective series resistance that would result from this scattering, and this resistance was found to be in agreement with the observed shifts in the conductance peaks near 1 and 3 $`G_0`$ for Cu, Ag and Au. On the other hand, the shift of the second peak was much larger and this led Ludoph et al. to propose that the origin of this peak may be a preferred atomic configuration, rather than conductance quantization. Li et al. PRB\_Li98  studied the influence of the presence of adsorbate molecules such as nitrogen or ethanol on the conductance of a gold nanowire. Plateaus are still observed in individual conductance traces, but their positions scatter more, resulting in a smearing of peaks in the histogram, which depends on the strength of the adsorbate used. This effect is interpreted as a result of the scattering of conduction electrons by the adsorbate molecules. It is noted that the first peak is less affected than the others.

The work presented here is an experimental study of the effect on conductance histograms of various levels of disorder. For this purpose, we have measured the histograms of copper-nickel alloys for various nickel concentrations, and compared them to the histogram of pure copper. The disorder in the point contact is increases with the nickel concentration.

## II Experimental method

Our samples are wires made of random alloys of copper and nickel. Two reasons have motivated the choice for this alloy: Firstly, its phase diagram is particularly simple, copper and nickel being completely miscible in all concentrations. Secondly, the histograms of copper and nickel are very different, as will be shown below. This makes it easy to recognize whether increasing the nickel concentration changes the histogram in a qualitative way, which would make it impossible to see the nickel atoms as impurities in the copper lattice, or whether it induces only quantitative changes, e.g. a shift of the peaks in the copper histogram.

The nickel concentrations used are 1, 10, 25, 35, 45 and 50 atomic %. The required amounts of copper and nickel were melted at 1500 C. The alloys were pulled into wires with a diameter of about 200 $`\mu `$m, which were subsequently annealed for 24 hours at 900C.

Metallic contacts of atomic size were then obtained by using the so called Mechanically Controllable Break Junction (MCBJ) technique. A piece of wire is notched in the middle and glued on a bendable phosphor-bronze substrate by two drops of epoxy. The sample is insulated from the substrate by a layer of kapton foil, and four measurement wires are connected to it with a conducting glue. The sample is first cooled down to 4.2 K and then broken in vacuum ($`10^5`$ mbar at room temperature) by bending the substrate. The bending is achieved by pressing a piezo element against the substrate by means of a differential screw. The high vacuum provided by the low temperature ensures atomically clean contacts after rupture. The two electrodes are then brought back into contact and their separation is controlled using the piezo element. In this way, the displacement between the two electrodes can be controlled with a precision of the order of 0.01 Å. The wire is broken and brought back into contact repeatedly by applying a sawtooth voltage to the piezo element. During the breaking stages current and voltage over the sample are measured. An histogram is build from approximately 3000 individual conductance traces. To ensure reproducibility, for each alloy at least two different samples were measured, while the contact was regularly ‘renewed’ by pressing the electrodes deep into each other.

## III Results

In Fig. 1, the histograms of the pure elements Cu and Ni and those of the copper alloys containing 1, 10, 25, 35, 45 and 50% nickel are presented. These measurements were done using a bias voltage of 10 mV. Topmost, the histogram of pure copper is shown. It is obtained by averaging ten histograms made on four different samples of pure copper. Three peaks can be distinguished, positioned at conductance values $`G`$ = 1.0 $`G_0`$, 1.7 $`G_0`$ and 2.7 $`G_0`$. The first peak has a higher amplitude and is sharper than the others. Close inspection makes clear that this first peak is actually slightly below 1 $`G_0`$. At the bottom of Fig. 1, the histogram of pure nickel is shown. It is very different from the copper histogram: it exhibits a broad first peak centered around 1.6 $`G_0`$, and a second, even broader, structure centered around 3.1 $`G_0`$.

Looking at the histograms of the alloys, it is remarkable that for all nickel concentrations the histograms resemble the histogram of pure copper; in particular, they have a first peak a little below $`1G_0`$. Indeed, for biases lower than 200 mV and for nickel concentrations up to 50%, there is no sign of a ”nickel” peak at 1.6 $`G_0`$. We give a possible interpretation of this effect in the next section. The second peak disappears rapidly when one increases the nickel concentration, due to shifting and broadening. With increasing nickel concentrations, three things can be seen to happen to the other peaks in the histogram: broadening, decreasing in amplitude and downshifting.

The width and height of the first peaks in the histogram are plotted in the right and left panel of Fig. 2, respectively. As can be seen, the amplitude of the peak decreases monotonously with Ni concentration. This can be explained as the result of increasing disorder. The disorder leads to slight variations in the last-plateau conductance of different scans. These variations cause the different plateaus to fall into different bins on the conductance axis, thus causing less high peaks in the histogram. The same effect gives rise to the broadening of the peaks that is observed for Ni concentrations up to 35%. Above this concentration, the width decreases to a minimum at 45%, after which it increases again. Fig. 2 also shows that the results reproduce for different biases. The main feature observed in Fig. 1 is the shift towards lower conductance values of the first and third peaks in the histogram. In Fig. 3, the series resistance corresponding to this shift is plotted as a function of the percentage of nickel in the alloy. Note that we do not calculate the shift with respect to integer multiples of $`G_0`$, but relative to the measured positions of the corresponding peaks in the histogram of pure Cu. Results again reproduce for other biases. The shift increases roughly linearly with nickel concentration, up to 35% nickel. Two surprising features of Fig. 3 must be noted. Firstly, the series resistance is not the same for all the peaks in the histogram. Instead, the resistance corresponding to the shift of the third peak is roughly two or three times higher than the corresponding values for the first peak. The same conclusion can be reached when considering the resistance necessary to bring the first and third peaks of the pure copper histogram to the corresponding integer multiples of $`G_0`$, namely, 60 $`\pm `$ 20 $`\mathrm{\Omega }`$ and 460 $`\pm `$ $`60\mathrm{\Omega }`$ respectively. This is in agreement with earlier findings PRB\_Hansen97 , but contradictory to the widespread practice to bring all conductance peaks to integer multiples of $`G_0`$ by subtracting a constant series resistance. The second surprising feature of Fig. 3 is the decrease of the shift for Ni concentrations higher than 35%. This decrease corresponds to the decrease in width that was already noticed in Fig. 2. Again the minimum is at 45% Ni. We discuss a possible interpretation of this effect in the next section.

## IV Discussion

The main difference between the histograms of copper and nickel, as shown in Fig. 1, is the position of the peaks. In particular, whereas copper shows a sharp first peak nearly exactly at 1 $`G_0`$, nickel has a broad first peak centered around 1.6 $`G_0`$. In both cases, however, it seems reasonable to assume that the first peak in the histogram corresponds to the smallest possible contact: only one atom bridging the gap between the electrodes. This is supported by the fact that the conductance invariably drops sharply to zero at the end of the plateau at $``$ 1 (Cu) or $``$ 1.6 $`G_0`$ (Ni). The difference in conductance between the different kind of atoms can be explained by taking into account their valence orbitals. As already mentioned in the introduction, copper has one spherically symmetric s valence orbital, whereas the valence orbitals of nickel are d-type. It has been shown Nature\_Scheer98  that for monovalent s-type metals there is one, almost completely opened, conductance channel per atom, which gives rise to a 1-atom conductance of nearly 1 $`G_0`$. In contrast, the five valence orbitals of a d-metal give rise to five partially opened channels, with a total conductance in the range of 1.5 - 2.5 $`G_0`$.

We noticed that the histograms for all the alloys resemble the histogram of pure copper, in particular, the first peak is below 1 $`G_0`$, and the third is below 2.7 $`G_0`$. However, if we very crudely model a contact of e.g. Cu50Ni50 as consisting of two pyramids touching each other in the apexes, such that there are two constricting atoms, 25% of all contacts should end in a Cu-Cu constriction, 25% should end in a Ni-Ni constriction, and 50% should have a Cu-Ni constriction. This raises the question why we don’t see the influence of the nickel atoms in the histogram of this alloy, apart from the downshifting and smearing of its peaks. An explanation for this phenomenon might be that the Ni atoms do not take part in the forming of the contact. This might be due to surface segregation, which is known to lead to a very strong copper enrichment of the surface layers in bulk Cu-Ni PRB\_Sakurai86 , due to the lower surface energy of Cu atoms. Calculations for small Cu-Ni clusters PRB\_Tomanek83  show the same phenomenon: all the copper atoms present segregate to the surface of the cluster. All the atoms forming the contact in a MCBJ can be viewed as being at the ‘surface’, to the extent that their coordination number is less than in the bulk. Thus there might be a preference for Cu atoms to form the contact, which would lead to the absence of Ni peaks in the histogram of the alloy. However, this surface segregation takes places only at elevated ($`>`$ 600 K) temperatures, while we performed our measurements at 4 K and at low (10 mV) bias, such that local heating of the contact is probably not sufficient for segregation to happen helko . Moreover, preliminary studies by molecular dynamics simulations jan  don’t show such a segregation between Cu and Ni atoms in the contact. Possibly, however, the repeated contact indentation and breaking provides enough mobility for the atoms to obtain some degree of segregation.

Another possibility to explain the similarity between the histograms of the alloys and that of pure copper is assuming that the conductance of Cu-Ni constriction is the same as that of a Cu-Cu constriction, namely, 1 $`G_0`$. For a contact of Cu50Ni50 in the form of two pyramids, this would lead to at least 75% of all last-atom contacts having a conductance of approximately 1 $`G_0`$. Calculations to verify this are under way.

If the above picture is correct, adding nickel atoms indeed doesn’t change the copper histogram in a qualitative way. Instead, the nickel atoms only act as scatterers for the electrons in the banks to the contact. One would expect the amount of disorder in the banks to increase with the nickel concentration. In Fig. 3, it is shown that the series resistance increases approximately linearly with nickel concentration up to a nickel concentration of 35%. From this we conclude that it is indeed the disorder in the nanowire, in the form of impurities, that causes the downshift of the peaks in the histograms. This is in accordance with the results of the calculations cited above.

However, the sudden decrease of the series resistance when going to 45 % nickel, accompanied by a decrease of the width of the first peak, suggest that the amount of disorder decreases when one increases the nickel concentration to this value. In order to compare with the behavior of bulk copper-nickel, we have measured the bulk resistivity of the alloys at room temperature as well as at $`T`$ = 4.2 K. Results are presented in Fig. 4.

It can be seen that the low-temperature resistivity initially closely follows the RT resistivity and is linear with the nickel concentration. Upon going to 45% nickel, however, the low-temperature resistivity decreases, whereas the room temperature resistivity still increases, although less than would be expected from a linear behavior. The decrease in the shift of the peaks in the copper-nickel histograms therefore probably reflects the decrease is resistivity of bulk copper-nickel. This shows that the origin of the series resistance is probably in the banks to the contact, and is not a (quantum mechanical) effect in the contact itself. A likely explanation for the decrease is the onset of ferromagnetism in the nanowire. The ferromagnetic ordering will decrease the amount of spin-spin scattering events by decreasing the magnetic fluctuations, even though the nickel concentration increases. As a matter of fact, bulk copper-nickel becomes magnetic at 4.2 K for a nickel concentration of approximately 45% handbook . From Fig. 3 it can be seen that the resistance tends to increase again when the nickel concentration is increased to 50%, which may be explained by assuming that once the alloy is ferromagnetic, the amount of disorder in it again increases with increasing nickel concentration. The inset in Fig. 4 shows the bulk mean free paths for the electrons in several of the alloys, calculated from the resistivities using the Drude model. An important difference between the bulk material and the quantum point contacts is the behavior for low nickel concentration. It is well known that the resistivity of metals at low temperatures depends sensitively on the purity of the sample. Thus, the mean free path in pure copper at 4.2 K will be much larger than that in Cu99Ni1. On the other hand, the histograms of these two materials are nearly indistinguishable. This means that in point contacts of both materials, the mean free path is restricted by scattering from the surface and lattice defects related to the geometry of the contact, rather than by impurities. When the concentration of impurities is raised to 10%, they start to have an appreciable effect. Thus, the mean free path in the contact region of a QPC of pure copper must be in between those in bulk Cu99Ni1 and Cu90Ni10, i.e. between 5 and 30 nm. This is consistent with results obtained earlier for quantum point contacts of copper bas .

From Fig. 3 a second important feature can be inferred. It can be seen clearly that the series resistance is not the same for all the peaks in the histogram. The resistance necessary to shift the third peak in the histograms back to its value for pure copper is larger than for the first peak. The difference increases with disorder. This effect is not accounted for by any of the calculations cited above. Possibly this difference in behavior compared to the first peak can be accounted for by considering the fluctuations in Cu/Ni ratio for contacts of this size, or by a non-uniform Ni concentration due to a partial segregation of Ni atoms away from the contact area.

## V Conclusion

In conclusion, we have shown, firstly, that adding nickel impurities to copper leads to an increase of the shift and smearing of the conductance peaks. This is evidence that disorder in the nanowires is a main source for the series resistance, in agreement with several theoretical papers. Up to a nickel concentration of 35%, this shift increases roughly linearly with the nickel concentration. Secondly, when the nickel percentage is increased further, the shift decreases. The latter effect may be due to a ferromagnetic transition leading to a decrease of the resistance of the mesoscopic banks. It was shown that the bulk resistivity of Cu-Ni alloys follows the same behavior as the series resistance does, which is evidence that the series resistance finds its origin in the banks to the contact rather than in the contact itself. Thirdly, we have shown that one single series resistance is not sufficient to shift all the peaks in the histogram back to the respective quantized values. Instead, it is found that the series resistance needed for the third peak is larger than that for the first. Fourthly, our measurements show that adding nickel impurities has remarkably little influence on the global shape of the copper histogram, even if the nickel concentration is raised to 50%. To explain this, we suggest that a Cu-Ni constriction has approximately the same conductance as a Cu-Cu constriction, or that there is a partial segregation of Ni atoms away from the contact.

## VI Acknowledgements

This research has been supported by a Marie Curie fellowship of the European Community under contract number HPMF-CT-1999-00196. We are grateful to I.K. Yanson for fruitful discussions and to R.W.A. Hendrikx for providing the sample materials.


text:
# Gravitational quasinormal modes for Anti-de Sitter black holes

## I Introduction

The spectrum of gravitational waves emitted by a black hole which has been disturbed is dominated, over intermediate timescales, by characteristic frequencies and decay timescales. These can be accurately described in terms of the discrete set of quasinormal modes. The frequencies of these modes can be determined numerically these days to high accuracy up to large order .

A new application of the quasinormal mode spectrum has arisen recently from superstring theory. There is a suggestion, known as the ADS-CFT correspondence, that string theory in anti-de Sitter space (usually with extra internal dimensions) is equivalent to conformal field theory in one less dimension . A black hole, which has a characteristic temperature fixed by the Hawking effect, should correspond to CFT at finite temperature and the quasinormal modes would describe non-equilibrium effects .

The temperature of the black hole is fixed by the surface gravity of the event horizon $`\kappa _1`$. The anti-de Sitter space also has a characteristic length scale $`a`$ which is conventionally refered to as the anti-de Sitter radius. In the non-rotating case, there are two black hole solutions for a given value of the surface gravity , a larger one with positive and smaller one with negative specific heat. Therefore the metric is determined by $`a`$, $`\kappa _1`$ and the choice of the larger or the smaller black hole.

Previous work on the quasinormal mode spectrum in anti-de Sitter space has concentrated on the wave equation for scalar fields. Chan et al.  studied the scalar field problem prior to the interest in the ADS-CFT correspondence. Horowitz et al.  examined the scalar spectrum in four and five dimensions. They found simple scaling behaviour of the spectrum in both limits $`r_1a`$ and $`r_1a`$, the former explained by an approximate symmetry of the metric in that limit.

Gravitational perturbations of the BTZ black hole in three dimensions and the anti-de Sitter black hole in four dimensions have been studied by Cardoso et al. . They make use of the fact that the gravitational perturbation equation is separable in four dimensions . This fact was also used in a preliminary report on our results . The main difference between our work and the work of Cardoso et al. is in the choice of boundary conditions on the axial and polar metric perturbations.

There is a symmetry which relates the axial and polar metric perturbations, whose existence can be explained by the fact that the covariant form of perturbation analysis based on the Weyl tensor is governed by a single equation . We argue that the thermal fluctuations in the energy of the conformal field theory are determined by a single set of modes for which the axial perturbations have dirichlet boundary conditions and the duality between axial and polar perturbations is preserved. We calculate the quasinormal mode frequencies numerically. For fluctuations in other components of the stress energy tensor there is a one parameter family of boundary conditions and the modes depend weakly on this parameter.

We have also investigated the quasinormal mode spectrum for black holes in de Sitter space. An interesting feature of this problem is that the asymptotic behaviour of the modes is under better control than the corresponding problem in flat space . The first few modes where evaluated previously in ref. . Using methods based on continued fractions now gives us many thousands of modes. A Pöshl-Teller approximation for the potential works extremely well in this example .

## II De Sitter Space

The spherically symmetric de Sitter spacetime which contains a black hole has two horizons, one at the event horizon of the black hole and a cosmological horizon further out. In standard coordinates, the metric for a black hole of mass $`M`$ and for a cosmological constant $`\mathrm{\Lambda }`$ can be written

$$ds^2=\frac{\mathrm{\Delta }}{r^2}dt^2+\frac{r^2}{\mathrm{\Delta }}dr^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)$$
(1)
where

$$\mathrm{\Delta }=r^22Mr\frac{1}{/}3\mathrm{\Lambda }r^4.$$
(2)
The event horizon lies at the smallest root $`r_1`$ of $`\mathrm{\Delta }`$, where the surface gravity $`\kappa _1=r^2\mathrm{\Delta }^{}/2`$. We can invert these relations to give

$`M`$ $`=`$ $`{\displaystyle \frac{1}{/}}3r_1(\kappa _1r_1+1)`$ (3)
$`\mathrm{\Lambda }`$ $`=`$ $`r_1^2(12\kappa _1r_1),`$ (4)
showing, incidentally, that $`\kappa _1r_1<1/2`$ in de Sitter space.

The gravitational wave equations were separated originally for charged black holes , then for uncharged  and finally for rotating black holes . The metric perturbations can be classified by their axial or polar symmetry under spatial inversions. The transformation theory relating the two types is described in the next section.

The axial metric perturbations with time dependence $`e^{i\omega t}`$ and angular mode $`l`$ are subject to a one dimensional scattering equation,

$$\frac{d^2\psi }{dr^2}+\left\{\omega ^2V(r)\right\}\psi =0$$
(5)
where

$$V(r)=\frac{\mathrm{\Delta }}{r^4}\left(l(l+1)\frac{6M}{r}\right)$$
(6)
The tortoise coordinate $`r^{}`$, given as usual by $`dr^{}=r^2dr/\mathrm{\Delta }`$, covers the range $`\mathrm{}<r^{}<\mathrm{}`$ as $`r`$ ranges between the two horizons. The potential decays to zero exponentially fast in the $`r^{}`$ coordinate at both horizons. This is an important property which is not present for the black hole in an asymptotically flat spacetime and guarantees that the reflection and transmission coefficients are meromorphic functions of the frequency .

By definition, the quasinormal modes are ingoing waves at the event horizon and outgoing waves at the cosmological horizon,

$$\psi \{\begin{array}{cc}e^{i\omega r^{}}\hfill & \text{as }r^{}\mathrm{}\hfill \\ e^{i\omega r^{}}\hfill & \text{as }r^{}\mathrm{}\hfill \end{array}$$
(7)
The method which we use to determine the complex frequencies follows the continued fraction approach devised by Leaver . The idea is to first scale out the divergent behavior at the cosmological horizon and examine the resulting Fröbenius series. If we begin with a purely divergent solution at the event horizon, the values of $`\omega `$ can be chosen to select the convergent rescaled solution at the cosmological horizon.

The analysis is slightly easier if we use the independent variable $`x=1/r`$. In de Sitter space, the radial equation has three regular singular points $`x_1`$ (the event horizon), $`x_2`$ (the cosmological horizon) and $`x_3`$. The indices of the equation at these points are $`\pm \rho _i=\pm i\omega /(2\kappa _i)`$, where

$`\kappa _1`$ $`=`$ $`2M(x_1x_2)(x_1x_3)`$ (8)
$`\kappa _2`$ $`=`$ $`2M(x_2x_3)(x_2x_1)`$ (9)
$`\kappa _3`$ $`=`$ $`2M(x_3x_1)(x_3x_2)`$ (10)
Note that $`\rho _1+\rho _2+\rho _3=0`$ and

$$e^{i\omega r^{}}=(xx_1)^{\rho _1}(xx_2)^{\rho _2}(xx_3)^{\rho _3}$$
(11)
We define,

$$\psi =e^{i\omega r^{}}u$$
(12)
For a quasinormal mode, $`u(x)`$ must be divergent at $`x_1`$ and finite at $`x_2`$. The Frobenius series is therefore

$$u(x)=(xx_1)^{2\rho _1}\underset{n=1}{}a_n\left(\frac{xx_1}{x_2x_1}\right)^n$$
(13)
where $`a_0=1`$ and the coefficients satisfy a recurrence relation of the form

$$\alpha _na_{n+1}+\beta _na_n+\gamma _na_{n1}=0.$$
(14)
The existence of this recurrence relation with only three terms is very important, because the convergence of the Fröbenius series at $`x=x_2`$ is now dependent on a continued fraction equation,

$$0=\beta _0\frac{\alpha _0\gamma _1}{\beta _1}\frac{\alpha _1\gamma _2}{\beta _2}\mathrm{}$$
(15)
which can be solved very quickly numerically.

Figures 1 and 2 show the real and imaginary parts of the frequency of the $`l=2`$ and $`l=3`$ modes as $`r_1\kappa _1`$ is varied. The quasinormal modes in flat space are recovered at $`r_1\kappa _1=0.5`$.

The simple shape of the potential suggests that a Pöschl-Teller approximation might be effective, as it is for the Schwarzchild black hole . The approximate potential

$$V_{PT}=V_0\mathrm{sech}^2(r^{}/b)$$
(16)
contains two free parameters which are used to fit the height and the second derivative of the potential at the maximum. The location of the maximum has to be found numerically.

The quasinormal modes of the Pöschl-Teller potential can be evaluated analytically,

$$\omega _n=\frac{1}{b}\left\{(n+1/2)i+\sqrt{b^2V_01/4}\right\}.$$
(17)
for $`n=1,2,\mathrm{}`$. These values are compared to the numerical values for $`l=2`$ in figure 1 and $`l=3`$ in figure 2. The approximation works best for the imaginary parts of the frequencies, but the differences between the numerical values and the Pöschl-Teller values decrease as $`l`$ is increased and soon become smaller than the resolution of such graphs.

For large $`l`$ it is possible to obtain the Pöschl-Teller fit to the potential analytically. The maximum of the potential in this case lies close to $`r=3M`$. From the second derivative at $`r=3M`$, we get

$$b^1=\kappa _1(1+\kappa _1r_1)^1(1+\frac{2}{/}3\kappa _1r_1)^{1/2}$$
(18)
and $`b^2V_0=l(l+1)`$. Equations (3) and (4) have been used to express the result in terms of the surface gravity. The analytic expression for the frequencies is therefore

$$\omega _n=\kappa _1\frac{(1+\frac{2}{/}3\kappa _1r_1)^{1/2}}{1+\kappa _1r_1}\left\{(n+1/2)i+(l(l+1)1/4)^{1/2}\right\}.$$
(19)
Figures 1 and 2 show the numerical results and the approximation for small $`l`$. The difference is largest in the flat space limit at the right hand side of the figures. As $`l`$ is increased, the difference beween the approximation and the numerical results soon becomes insignificant.

## III Anti-de Sitter space

The metric for a black hole in anti-de Sitter space takes the same form as for de Sitter space but with $`\mathrm{\Lambda }`$ negative. The cosmological horizon is no longer present. We introduce an anti-de Sitter radius $`a`$ by the definition $`a^2=3/\mathrm{\Lambda }`$. For fixed values of $`\kappa _1`$ and $`a`$ there are two black hole solutions with different radii,

$$\kappa _1r_1=\frac{1}{/}3(\kappa _1a)^2\pm \frac{1}{/}3\left((\kappa _1a)^43(\kappa _1a)^2\right)^{1/2}$$
(20)
The larger solution has positive and the smaller one negative specific heat.

The axial metric perturbations are subject to the same radial equation as for de Sitter space, however in anti-de Sitter space the tortoise coordinate lies in the range $`\mathrm{}<r^{}0`$. The boundary conditions at $`r^{}=0`$ can be dirichlet, neumann or robin depending respectively on whether the field, its derivative or a combination of both vanishes. Previous work has tended to opt for dirichlet boundary conditions , based largely on ADS invariance scalar field theory . However, dirichlet boundary conditions break a duality in the equations between axial and polar metric perturbations . We will introduce a set of self-dual boundary conditions based on the ADS-CFT correspondence.

In the Newman-Penrose formalism, the gravitational perturbations can be described by two Weyl scalars $`\mathrm{\Psi }_0`$ and $`\mathrm{\Psi }_4`$. In flat space, $`\mathrm{\Psi }_0`$ would represent gravitational waves which are incoming at infinity and $`\mathrm{\Psi }_4`$ would represent outgoing gravitational waves. The equations for the Weyl scalars separate to give a generalisation of the Teukolski equations . We put, for example,

$$\mathrm{\Psi }_0=e^{i\omega t}\frac{r^3}{\mathrm{\Delta }}Y(r)S_l(\theta )$$
(21)
where $`S_l(\theta )`$ can be found in Chandrasekhar . The radial equation for $`Y(r)`$ can be related to the radial equation for the metric perturbations (5) if we set

$$Y=V^\pm \psi ^\pm +(W^\pm +2i\omega )\mathrm{\Lambda }_+\psi ^\pm $$
(22)
where

$$\mathrm{\Lambda }_+=\frac{d}{dr^{}}+i\omega .$$
(23)
There are two equivalent formulations corresponding to polar and axial metric perturbations, denoted by the plus and minus signs. The axial form, corresponding to $`\psi ^{}`$, has $`V^{}(r)=V(r)`$ given previously (6).

The axial and polar functions are related by

$$C^{}\psi ^{}=(C^{}+72M^2f)\psi ^\pm 12M\mathrm{\Lambda }_{}\psi ^\pm $$
(24)
with $`C^\pm =\mu ^2(\mu ^2+2)\pm 12iM\omega `$, and

$$f=\frac{\mathrm{\Delta }}{r^3}\frac{1}{\mu ^2r+6M}.$$
(25)
where $`\mu ^2=(l1)(l+2)`$.

The ADS-CFT correspondence conjectures that when the Anti-de Sitter space is conformally compactified, the generating functions of the conformal field theory on the boudary are related to the partition function of the superstring theory on the interior, or

$$\mathrm{exp}\left(\widehat{h}_{\mu \nu }T^{\mu \nu }\right)_{CFT}=Z\left(\widehat{h}_{\mu \nu }\right)$$
(26)
If we calculate the correlation fuctions of a polar quantity, such as the energy, we must therefore set the axial metric perturbation to vanish on the boundary. The transformation theory (24) then implies a Robin boundary condition on the polar metric perturbation, i.e. the normal derivative is determined by the value of the field. The same set of quasinormal modes is obtained from the axial and polar equations and the duality between the two is preserved.

In general, both the axial and polar metric variations are non-vanishing. This is consistent with the conformal theory viewpoint, where the six components of the metric perturbations in three dimensions are subject to three coordinate transformations and the conformal symmetry, leaving two independent functions. If we define the ratio $`\gamma =\psi ^+/\psi ^{}`$ , the two metric variations are related by (24), which gives a robin boundary condition

$$(\xi +\mathrm{\Lambda }_{})\psi ^{}$$
(27)
where

$$\xi =(1\gamma )\frac{C^+}{12M}+\frac{6M}{\mu ^2a^2}$$
(28)
Some special solutions lead to vanishing of the Newmann-Penrose scalars and might be regarded as pure gauge modes from the point of view of the ADS theory. For example,

$$\psi ^{}=e^{i\omega r^{}}\left(1+\frac{6M}{\mu ^2r}\right)$$
(29)
with $`\omega =\mu ^2(\mu ^2+2)i/12M`$, obtained by setting $`\psi ^+=0`$, is a purely axial perturbation with $`Y=0`$. Replacing $`\omega `$ by $`\omega `$ gives a quasinormal mode with frequency

$$\omega =\frac{\mu ^2(\mu ^2+2)}{12M}i$$
(30)
Such modes appear in the numerical solutions discussed below.

The frequencies of the quasinormal modes can be found numerically by adapting the method used for de Sitter space. We first scale out the behaviour of the quasinormal mode at the event horizon,

$$\psi ^{}=e^{i\omega r^{}}u,$$
(31)
where $`u`$ will be regarded as a function of $`x=1/r`$. The radial equation (5) gives the equation

$$pu^{\prime \prime }+p^{}u^{}+2i\omega u^{}+\left\{l(l+1)6Mx\right\}u=0$$
(32)
where $`p=a^2x^2+2Mx^3`$. We take the Fröbenius series which is regular at the event horizon,

$$u=\underset{n=0}{\overset{\mathrm{}}{}}a_n\left(\frac{xx_1}{x_1}\right)^n$$
(33)
with $`a_0=1`$. The coefficients are obtained from a recurrence relation similar to (14). The series can then be substituted into the boundary condition $`u=0`$ at $`x=0`$ to obtain an equation for the frequencies of the quasinormal modes with vanishing axial perturbations.

The frequencies of the first few modes using dirichlet boundary conditions on the axial perturbations have been plotted for $`a\kappa _1=2.0`$ and $`r_1=a`$ in figure 3. For $`l=2`$, one of the frequencies lies on the imaginary axis. (This is also found by Cardoso et al ). We have found a single frequency on the imaginary axis and can trace it for $`\kappa _1r_11.3`$ when $`l=2`$, for $`\kappa _1r_13.4`$ when $`l=3`$ and for $`\kappa _1r_138.0`$ when $`l=10`$.

The dependence of the first three modes on the metric parameters, with $`l=2`$ and $`l=10`$, is displayed in figures 3 and 4. Both the real and imaginary parts of the frequencies become proportional to $`\kappa _1`$ for large $`\kappa _1`$, which agrees with the behaviour for the scalar field quasinormal modes . For $`r_1a`$, the numerical data suggests that imaginary part of the frequency $`\omega _I`$ approaches zero as in the scalar case .

Corresponding results for dirichlet boundary conditions on the polar perturbations ($`\gamma =0`$) are shown in figures 5 and 6. The special mode in this case is the ‘pure gauge’ mode (30) and the analytic solution provides a good check on the accuracy of our numerical technique. The frequencies for these boundary conditions are not drastically different from the dirichlet axial boundary conditions. The two sets of results represent two extreme cases of the one parameter family of boundary conditions (27) given above, and indicate a weak dependence on the parameter $`\gamma `$.

The potential $`V(r)`$ for anti-de Sitter space no longer vanishes as $`r\mathrm{}`$ and it is harder to fit the potential to a Pöschl-Teller form. More importantly, the value of $`r^{}`$ at infinity is zero, and the asymptotic forms of the mode functions for the Pöschl-Teller potential can no longer be used.

## IV conclusion

The quasinormal modes for the gravitational perturbations of a black hole in de Sitter space have much in common with the quasinormal modes for the black hole in asymptotically flat spacetime. The radial scattering problem in de Sitter space is better behaved because the potential falls off exponentially rather than as a power law. This explains why the approximation to the quasinormal mode frequencies, derived in section 2 from a Pöschl-Teller potential, is in such good agreement with the numerical results for even moderate angular order $`l`$.

The corresponding problem in anti-de Sitter space requires careful consideration of the boundary conditions at infinity. The proposal here has been to use the ADS-CFT correspondence and the duality symmetry between axial and polar metric perturbations. For correlation functions of the energy in the conformal field theory, the appropriate boundary conditions are that the axial metric perturbations vanish and the polar perturbations satisfy a robin boundary condition.

In general, two functions have to be specified on the boundary and there is a one parameter family of boundary conditions. The effect on the quasinormal frequencies is quite small.

It would be desirable to analyse thermal widths in three dimensional conformal field theory to look for signs of a correspondence with the quasinormal mode frequencies found here. This could give an impetus for analysing the five dimensional gravitational quasinormal modes, where a full separation of the gravitational perturbation equations is still an unsolved problem.