text:
# Directional tunneling escape from nearly spherical optical resonators

(07/17/2003, Phys. Rev. Lett. 91, 033902 (2003))

## Abstract

We report the surprising observation of directional tunneling escape from nearly spherical fused silica optical resonators, in which most of the phase space is filled with non-chaotic regular trajectories. Experimental and theoretical studies of the dependence of the far-field emission pattern on both the degree of deformation and the excitation condition show that non-perturbative phase space structures in the internal ray dynamics profoundly affect tunneling leakage of the whispering gallery modes.

Electromagnetic fields in a uniform dielectric sphere can be calculated in much the same way as quantum mechanical wave functions in a spherically symmetric potential. The sphere exhibits whispering-gallery (WG) modes, which are long-lived resonances with electromagnetic fields confined near the surface Nussenzveig . For small deviations from the spherical shape, $`ϵ(R_{max}R_{min})/(R_{max}+R_{min})1`$ (where $`R_{max},R_{min}`$ are the maximal and minimal radii), perturbative treatments LaiHM  are routinely employed to infer deformation parameters from the splitting of azimuthal degeneracy of the WG modes Seguin . Certain strongly non-spherical resonators, on the other hand, can be analyzed with methods from quantum chaos such as random-matrix theory or periodic-orbit expansions Stockmann99 . Many generic resonator shapes, however, fall into a transition regime in which none of these known methods apply globally. When entering this regime from the perturbative side, calculations may encounter singularities and undefined limits BerryMVPhysTod . Experimental studies in this regime can thus provide insight into how nature resolves the competition between perturbative and non-perturbative physics, here with the resonator shape as a control parameter.

In this paper we present studies of far-field emission patterns and resonance lifetimes of deformed fused-silica microspheres. The experimental and theoretical studies of the dependence of the far-field emission pattern on both the degree of deformation and the excitation condition show that highly directional tunneling escape can occur in microspheres with small deformation ($`ϵ1\%`$) and with a large size parameter ($`kR=2\pi R/\lambda 785`$, where $`\lambda 800`$nm, $`R100\mu `$m). These results are completely unexpected in the ray optics limit or in the earlier perturbative wave treatment with $`ϵ`$ as a small parameter. The observed emission pattern shows that on one hand, the ray model breaks down in predicting the escape mechanism of WG modes for resonators with small deformation. On the other hand, the nonperturbative phase space structures predicted by the ray model can profoundly affect tunneling leakage of the WG modes. Even at extremely small deformation, nonperturbative phase space structures still persist and are directly responsible for directional tunneling escape. This intricate interplay between ray and wave dynamics provides essential physical insights into properties of weakly deformed WG resonators, especially, for those in which the phase space is filled with non-chaotic regular trajectories. Directional emission patterns have previously been observed in more strongly deformed resonators where a significant fraction of the internal rays shows chaotic motion ChangCSN00 .

Deformed microspheres were formed by melting together two spheres of similar sizes. The individual spheres were fabricated by melting an optical fiber tip with a focused $`\mathrm{CO}_2`$ laser beam. The spheres were brought into contact and heated until surface tension produced a completely convex surface. By carefully controlling the temperature of the glass using the $`\mathrm{CO}_2`$ laser, it was possible to repeatedly reduce the degree of deformation. For reference, we define the elongated axis of the resulting prolate spheroid, which also connects the centers of the two original spheres, to be the $`x`$-axis. The $`z`$-axis is defined by the remaining fiber stem, which breaks the rotational symmetry about the $`x`$-axis, making the deformed microsphere completely non-axisymmetric. The stem was held by a fiber chuck for easy manipulation of the sphere. Images of these deformed microspheres taken from three orthogonal directions have been shown in an earlier study LaceyW01 .

To investigate emission properties and resonance lifetimes of the deformed microspheres, we use frustrated total internal reflection at the surface of a prism with refractive index $`n=1.7`$ to launch individual traveling-wave WG modes near the $`x`$-$`y`$ plane of the spheroid (Fig. 1a inset) newfocus . The initial angle of incidence $`\chi _0`$ with which the WG modes are launched in the microsphere can be controlled by adjusting the angle of incidence $`\psi `$ in the prism.

WG mode spectra were obtained by measuring the far-field emission intensity as a function of the excitation wavelength. Since the stem holding the resonator is a strong leakage pathway, the observed long-lived WG modes must have internal field patterns that do not overlap with the stem region. This is made possible by the fact that the resonator is slightly flattened in the $`z`$ direction (along the stem axis). Experiments as well as ray simulations show that this deformation stabilizes the rays in the vicinity of the $`x`$-$`y`$ plane. Our directionality measurements correspondingly were performed in this plane. By scanning the detector while keeping the resonator fixed, we recorded WG mode spectra at various angles, $`\varphi `$ (from the $`x`$-axis), to construct the far-field emission pattern at a given microsphere deformation $`ϵ`$. The microsphere was then reheated to reduce $`ϵ`$ and the far-field emission pattern was measured again, leading to the evolution shown in Figure 1. The excitation beam is s-polarized, and the detection scheme is polarization-insensitive. The far-field emission patterns are independent of the polarization of the excitation laser beam. Spheroids with similar $`ϵ`$ exhibit qualitatively the same behavior in terms of mode spectra, Q-factors, and emission patterns.

Beginning with a strongly deformed microsphere, the emission was observed to have a strong peak at $`\varphi =45^{}`$ in the far field \[see Fig. 1 (a)\]. Since clockwise traveling waves were excited, as indicated by the white arrow in the sphere image in Fig. 1 (a), this far-field emission direction corresponds to light escaping tangential from the region at $`\varphi =45^{}`$ on the surface of the microsphere, shown by the solid arrow in the sphere image. A much smaller peak ($`5\%`$ of the large peak height) in the far-field emission pattern was observed at $`\varphi =135^{}`$, corresponding to light escaping from the region at $`\varphi =45^{}`$ on the microsphere surface (dotted arrow). As the deformation was reduced, the far-field emission peak at $`\varphi =135^{}`$ grew to about one quarter of the $`\varphi =45^{}`$-peak height in Fig. 1 (b). In Fig. 1 (c), the two peaks reached nearly equal intensity. The bright emission spots at $`\varphi =\pm 45^{}`$ on the resonator surface can be viewed directly with a CCD camera. For reference, we call the pattern in Fig. 1 (a) *asymmetric* and the pattern in Fig. 1 (c) *symmetric* (around $`\varphi =90^{}`$). Figures 1 (e)-(g) also show that as the deformation was reduced, the Q-factor of the relevant WG modes increased by nearly four orders of magnitude.

The measurements discussed thus far were performed with input condition $`\chi _090^{}`$ ($`\mathrm{sin}\chi _01`$). Figure 1 (d) shows the far-field emission pattern from the same microsphere as Fig. 1 (c) but with light launched at $`\mathrm{sin}\chi _00.8`$. The emission pattern is nearly identical to Fig 1 (c) where $`\mathrm{sin}\chi _01`$, although the corresponding Q-factor is nearly two orders of magnitude smaller \[Figs. 1 (g) and (h)\].

The observed emission patterns in Fig. 1 (a)-(d) *disagree* with the intuitive expectation that WG modes in an oval resonator should preferentially emit tangential to the points of highest curvature, into the far-field direction $`\varphi =\pm 90^{}`$. This is what one obtains when modeling our spheroids as triaxial *ellipsoids*, whose internal ray dynamics exhibits no chaos WaalkensWD99 , independently of the axis ratios.

The observed emission pattern of the most deformed microsphere, Fig. 1 (a), is well explained by a ray model. The peak at $`\varphi =45^{}`$ can be attributed to an effect known as dynamical eclipsing ChangCSN00 . In the ray model, a WG mode corresponds to light rays trapped close to the perimeter of the dielectric resonator by total internal reflection, which prevents light escape unless a *critical angle* $`\chi _c\mathrm{arcsin}1/n=43.6^{}`$ is reached, where $`n=1.45`$ is the refractive index of the fused silica. At small but finite $`ϵ`$, any oval can be approximated by the first terms of a multipole expansion, which after proper choice of origin is *quadrupolar*. In this limit, a stable 4-bounce orbit shaped like a diamond forms with its sharp vertices at the highest-curvature points of the resonator with an angle of incidence $`\chi _445^{}`$ \[inset to Fig. 1 (d)\]. This creates “islands of stability” in phase space Stockmann99 : rays launched near the diamond orbit will retain similar reflection points and angles even if they do not close onto themselves after four bounces. At the deformation used in Fig. 1 (a), most of the phase space supporting the WG mode in question is chaotic. However, chaotic rays cannot penetrate the 4-bounce islands; this means chaotic WG rays are prevented from refractively escaping at the points of highest curvature, because $`\chi _4\chi _c`$. Instead, as trajectories flow around the islands in one direction, refractive escape occurs near $`\varphi =45^{}`$ on the resonator surface. From this circulation around the islands, *asymmetric* far-field emission patterns result as a hallmark of refractive escape. This is what we observe in Fig. 1 (a).

This ray mechanism was originally proposed and tested for effectively 2D systems where islands of stability and chaotic regions are mutually exclusive. In our 3D spheroids, the same diamond-shaped stable orbit exists in or near the $`x`$-$`y`$ plane, and refractively escaping rays in its vicinity have lifetimes that correspond to $`10`$ round trips in the resonator, translating to a Q-factor near $`10^4`$. During this time, a given ray can be considered as moving in a cross-sectional plane that may be inclined and slowly precessing around the $`z`$-axis. The variation of cross-sectional shape, $`\mathrm{\Delta }ϵ`$, probed by such rays causes no significant differences in the size of the 4-bounce islands ChangRC99 . Therefore, rays in the 3D resonator behave as in a corresponding fixed 2D resonator during a relatively short time preceding an escape event, and in particular show dynamical eclipsing.

However, *refractive* escape ceases to be the dominant leakage mechanism at the very small deformation used in Figs. 1 (c) and (d), since in this case chaotic diffusion is no longer effective and the stable diamond orbit itself also becomes fully confined by total internal reflection. Those classical rays that do escape refractively, moreover, are not able to produce the *symmetric* emission patterns observed in Figs. 1 (c) and (d) Foot1 . Hence, the emission process in these cases is governed by evanescent escape, i.e. a tunnel coupling between internal and radiation fields.

Tunneling is the only decay mechanism in an ideal sphere where each WG mode in the plane of excitation corresponds to a unique azimuthal angular momentum number $`m`$ with respect to the $`z`$-axis, which is semiclassically related to $`\chi _0`$ by $`\mathrm{sin}\chi _0=m/(nkR)`$, using the fact that the modes of interest remain close to the $`x`$-$`y`$ plane and hence have total angular momentum $`lm`$. At $`kR1`$, the tunneling escape rate is negligible when $`\mathrm{sin}\chi _01`$ but accelerates exponentially as $`\mathrm{sin}\chi _0`$ approaches $`\mathrm{sin}\chi _c`$. In a deformed microsphere, the angle of incidence $`\chi `$ varies as a well-defined function of $`\varphi `$ provided it covers an angular momentum range where chaos is neglibible NockelS97 . A circulating ray with varying $`\chi (\varphi )>\chi _c`$ will then escape with exponentially strong selectivity near the *minima* of $`\chi (\varphi )`$.

The emission locations and directions observed in the experiment indicate that the minima of $`\chi (\varphi )`$ lie near the corners of the unstable rectangular 4-bounce orbit \[inset to Fig. 1 (d)\]. We determined from ray tracing that this only occurs near but above $`\chi _c`$. Thus, the symmetric far-field patterns in Figs. 1 (c) and (d) are a *tunneling probe* for the minima of the *ray-optical* $`\chi (\varphi )`$ in the vicinity of the critical angle $`\chi _c`$. To further corroborate this interpretation, recall that the emission patterns as shown in Figs. 1 (c) and (d) are insensitive to the initial angle of incidence, $`\chi _0`$, as they should be if the detected light originates near $`\mathrm{sin}\chi \mathrm{sin}\chi _c`$. The question remains how the light coupled into the resonator is able to reduce its angular momentum from a high value corresponding to $`\mathrm{sin}\chi _0`$ to the much lower $`\mathrm{sin}\chi _c`$. Our observation shows that such a dynamical process is present in the spheroid, but does not unambiguously reveal its mechanism. A possible mechanism is Arnol’d diffusion, a phase-space transport process that crucially depends on the fact that the resonators are in fact three-dimensional LaceyW01 , with no axial symmetry. Even without resolving this question, we can nevertheless conclude that the phase space region near $`\mathrm{sin}\chi _c`$ is responsible for the observed emission patterns.

Our initial assumption of a 2D model for the emission process, justified by the approximate 2D nature of the 3D ray trajectories for intermediate times as discussed earlier, is therefore consistent with the fact that no long-time properties of the internal ray dynamics are needed to explain the emission directions. With the same justification, we also performed numerical wave calculations by modeling the shape felt by modes in the $`x`$-$`y`$ plane as the cross section of a quadrupolar cylinder. Figures 2 (a) and (b) plot the intensity patterns of two WG modes under traveling-wave excitation NockelS97 ; BraunILNSSVWW00 . The calculation assumes s-polarized light and uses size parameters near $`kR113`$ numericnote .

The emission patterns in Figs. 2 (a) and (b) agree qualitatively with the asymmetric and the symmetric emission patterns shown in Figs. 1 (b) and (d), respectively. To understand the physical mechanism for the emission in Fig. 2 (b), we plot in Fig. 2 (c) the distribution of angular momentum numbers for the WG mode in (b). Figure 2 (c) shows negligible overlap with the window for refractive escape below the critical angle $`\chi _c`$ (corresponding to a *critical angular momentum* $`m_c=kR113`$). This indicates that the emission from this mode is due to *tunneling* escape. By comparison, both refractive and tunneling escape contribute to the emission in Fig. 2 (a), with the refractive escape playing a dominant role.

The peak splitting of width $`\mathrm{\Delta }m4`$ in Fig. 2 (c) is a straightforward consequence of the oval deformation: as the wave circulates around the resonator with varying radius between $`R_{min}`$ and $`R_{max}`$, its angular momentum oscillates but has high probabitity of being near its extrema, given by the extrema of $`mr(\varphi )nkR\mathrm{sin}\chi (\varphi )`$ over the surface. The fact that the minima of $`\chi (\varphi )`$ do *not* occur at $`\varphi =0^{},\mathrm{\hspace{0.17em}180}^{}`$ on the surface is due to the 4-bounce island structure in the vicinity of $`\mathrm{sin}\chi =\mathrm{sin}45^{}`$.

For comparison, Fig. 2 (c) also shows the distribution of angular momentum number of an ultra-high Q WG mode with emission directionality as expected for the ellipse \[see Fig. 2 (d)\]. The mode is confined at a high $`m`$, corresponding to $`\mathrm{sin}\chi _00.91`$. Note that here the tunneling loss is negligible compared with scattering or absorption loss of the material that limits the actual Q-factor of fused silica microspheres to of order $`10^9`$ Gorodetsky . The two modes in Fig. 2 (c) have practically no overlap in angular momentum space. The qualitative difference in the far-field patterns of these two modes further confirms that the mode shown in Fig. 1 (c), which was excited at $`\mathrm{sin}\chi _01`$, is not confined to $`\mathrm{sin}\chi _01`$.

In essence, our experiment exploits the peculiar feature that for fused silica the critical angle lies near the phase-space islands corresponding to the stable 4-bounce orbit. These islands are the dominant structure in the WG region and persist even at small $`ϵ`$ where most of the phase space is filled with non-chaotic, “regular” trajectories. The island formation is a *non-perturbative* consequence of the breakdown of conservation laws WaalkensWD99 , which qualitatively distinguishes the generic quadrupolar shape from an ideal ellipse even though the two shapes differ only to second order in $`ϵ`$. The emission characteristics of a silica spheroid in the seemingly trivial small-$`ϵ`$ regime will *continue* to be strongly affected by ray patterns that wrap around the perimeter in approximately four bounces.

Whether any non-perturbative structure in the ray dynamics can be *resolved* by the wave field, depends on $`kR`$ PrangeNZ99 : a directionality measurement will be able to distinguish the peaks at $`\varphi =90^{}\pm 45^{}`$ if the conjugate angular momentum $`m`$ satisfies the uncertainty relation $`\mathrm{\Delta }\varphi \mathrm{\Delta }m1/2`$ with $`\mathrm{\Delta }\varphi \pi /4`$. This implies $`\mathrm{\Delta }m>2/\pi 1`$, which in our spheroids translates to a fluctuating angle of incidence of $`\mathrm{\Delta }\mathrm{sin}\chi =\mathrm{\Delta }m/(nkR)10^3`$. Even the least deformed resonator ($`ϵ1\%`$) exceeded this estimated resolution threshold significantly.

It is truly remarkable that this intricate interplay between ray and wave dynamics for which chaos plays no dominant role, can be exploited to engineer WG modes that can feature both high-Q and highly directional emission. This makes WG resonators with *small* deformation highly promising for a variety of applications Kippenberg ; Armani , such as microlasers, nonlinear optics, and quantum information processing.

This work was supported in part by NSF under grants No. DMR9733230 and No. DMR0201784.


format: text

text:
# GPCALMA: a Grid Approach to Mammographic Screening

## 1 Introduction

A reduction of breast cancer mortality in asymptomatic women is possible in case of early diagnosis , which is available thanks to screening programs, a periodical mammographic examination performed for $`49`$-$`69`$ years old women. The GPCALMA Collaboration aims at the development of tools that would help in the early diagnosis of breast cancer: Computer Assisted Detection (CAD) would significantly improve the prospects for mammographic screening, by quickly providing reliable information to the radiologists.

A dedicated software to search for massive lesions and microcalcification clusters was developed recently ($`1998`$-$`2001`$): its best results in the search for massive lesions (microcalcification clusters) are $`94\%`$ ($`92\%`$) for sensitivity and $`95\%`$ ($`92\%`$) for specificity. Meanwhile, in view of the huge distributed computing effort required by the CERN/LHC collaborations, several GRID projects were started. It was soon understood that the application of GRID technologies to a database of mammographic images would facilitate a large-scale screening program, providing transparent and real time access to the full data set.

The data collection in a mammographic screening program will intrinsically create a distributed database, involving several sites with different functionality: data collection sites and diagnostic sites, i.e. access points from where radiologists would be able to query/analyze the whole distributed database. The scale is pretty similar to that of LHC projects: taking Italy as an example, a full mammographic screening program would act on a target sample of about $`6.8`$ million women, thus generating $`3.4`$ millions mammographic exams/year. With an average size of $`60MB`$/exam, the amount of raw data would be in the order of $`200TB`$/year: a screening program on the European scale would be a data source comparable to one of the LHC experiments.

GPCALMA was proposed in 2001, with the purpose of developing a G̈RID application,̈ based on technologies similar to those adopted by the CERN/ALICE Collaboration. In each hospital, digital images will be stored in the local database and registered to a common service (Data Catalogue). Data describing the mammograms, also known as metadata, will also be stored in the Data Catalogue and could be used to define an input sample for any kind of epidemiology study. The algorithm for the image analysis will be sent to the remote site where images are stored, rather than moving them to the radiologist’s sites. A preliminary selection of cancer candidates will be quickly performed and only mammograms with cancer probabilities higher than a given threshold would be transferred to the diagnostic sites and interactively analysed by one or more radiologists.

## 2 The GPCALMA CAD Station

The hardware requirements for the GPCALMA CAD Station are very simple: a PC with SCSI bus connected to a planar scanner and to a high resolution monitor. The station can process mammograms directly acquired by the scanner and/or images from file and allows human and/or automatic analysis of the digital mammogram. The software configuration for the use in local mode requires the installation of ROOT and GPCALMA, which can be downloaded either in the form of source code from the respective CVS servers. The functionality is usually accessed through a Graphic User Interface, or, for developers, the ROOT interactive shell. The Graphic User Interface (fig. 1) allows the acquisition of new data, as well as the analysis of existing ones. Three main menus drive the creation of (access to) datasets at the patient and the image level and the execution of CADe algorithms. The images are displayed according to the standard format required by radiologists: for each image, it is possible to insert or modify diagnosis and annotations, manually select the Regions of Interest (ROI) corresponding to the radiologists geometrical indication. An interactive procedure allows zooming, either continously or on a selected region, windowing, gray levels and contrast selection, image inversion, luminosity tuning.

The human analysis produces a diagnosis of the breast lesions in terms of kind, localization on the image, average dimensions and, if present, histological type. The automatic procedure finds the ROI’s on the image with a probability of containing an interesting area larger than a pre-selected threshold value.

## 3 Grid Approach

The amount of data generated by a national or european screening program is so large that they can’t be managed by a single computing centre. In addition, data are generated according to an instrinsically distributed pattern: any hospital participating to the program will collect a small fraction of the data. Still, that amount would be large enough to saturate the available network connections.

The availability of the whole database to a radiologist, regardless of the data distribution, would provide several advantages:

* the CAD algorithms could be trained on a much larger data sample, with an improvement on their performance, in terms of both sensitivity and specificity.
* the CAD algorithms could be used as real time selectors of images with high breast cancer probability (see fig. 2): radiologists would be able to prioritise their work, with a remarkable reduction of the delay between the data acquisition and the human diagnosis (it could be reduced to a few days).
* data associated to the images (i.e., metadata) and stored on the distributed system would be available to select the proper input for epidemiology studies or for the training of young radiologists.

These advantages would be granted by a GRID approach: the configuration of a Virtual Organisation, with common services (Data and Metadata Catalogue, Job Scheduler, Information System) and a number of distributed nodes providing computing and storage resources would allow the implementation of the screening, tele-training and epidemiology use cases. However, with respect to the model applied to High Energy Physics, there are some important differences: the network conditions do not allow the transfer of large amounts of data, the local nodes (hospitals) do not agree on the raw data transfer to other nodes as a standard and, most important, some of the use cases require interactivity.

According to these restrictions, our approach to the implementation of the GPCALMA Grid application is based on two different tools: AliEn  for the management of common services, PROOF  for the interactive analysis of remote data without data transfer.

### 3.1 Data Management

The GPCALMA data model foresees several Data Collection Centres 2 , where mammograms are collected, locally stored and registered in the Data Catalogue. In order to make them available to a radiologist connecting from a Diagnostic Centre, it is mandatory to use a mechanism that identifies the data corresponding to the exam in a site-independent way: they must be selected by means of a set of requirements on the attached metadata and identified through a Logical Name which must be independent of their physical location. AliEn implements these features in its Data Catalogue Services, run by the Server: data are registered making use of a hierarchical namespace for their Logical Names and the system keeps track of their association to the actual name of the physical files. In addition, it is possible to attach metadata to each level of the hierarchical namespace. The Data Catalogue is browsable from the AliEn command line as well as from the Web portal; the C++ Application Program Interface (API) to ROOT is under development. Metadata associated to the images can be classified in several categories: patient and exam identification data, results of the CAD algorithm analysis, radiologist’s diagnosis, histological diagnosis. etc.. Some of these data will be directly stored in the Data Catalogue, but some of them may be stored in dedicated files and registered: the decision will be made after a discussion with the radiologists.

A dedicated AliEn Server for GPCALMA has been configured , in collaboration with the AliEn development team. Fig. 3 shows a screenshot from the WEB Portal.

### 3.2 Remote Data Processing

Tele-diagnosis and tele-training require interactivity in order to be fully exploited, while in the case of screening it would be possible - altough not optimal - to live without. The PROOF Parallel ROOt Facility system  allows to run interactive parallel processes on a distributed cluster of computers. A dedicated cluster of several PCs was configured and the remote analysis of a digitised mammogram without data transfer was recently run. As soon as input selection from the AliEn Data Catalogue will be possible, more complex use cases will be deployed. The basic idea is that, whenever a list of input Logical Names will be selected, that will be split into a number of sub-lists containing all the files stored in a given site and each sub-list will be sent to the corresponding node, where the mammograms will be analysed.

## 4 Present Status and Plans

The project is developing according to the original schedule. The CAD algorithms were rewritten in C++, making use of ROOT, in order to be PROOF-compliant; moreover, the ROOT functionality allowed a significant improvement of the Graphic User Interface, which, thanks to the possibility to manipulate the image and the associated description data, is now considered fully satisfactory by the radiologists involved in the project. The GPCALMA application code is available via CVS server for download and installation; a script to be used for the node configuration is being developed.The AliEn Server, which describes the Virtual Organisation and manages its services, is installed and configured; some AliEn Clients are in use, and they will soon be tested with GPCALMA jobs. The remote analysis of mammograms was successfully accomplished making use of PROOF. Presently, all but one the building blocks required to implement the tele-diagnosis and screening use cases were deployed. As soon as the implementation of the data selection from the ROOT shell through the AliEn C++ API will be available, GPCALMA nodes will be installed in the participating hospitals and connected to the AliEn Server, hosted by INFN. Hopefully, that task will be completed by the end of 2004.

## 5 Acknowledgments

The authors wish to thank the AliEn development team for their support and guidance in the installation and configuration of the GPCALMA server.


text:
# Scaling and renormalization in fault-tolerant quantum computers1footnote 11footnote 1Based on the talk given at the Simons Conference on Quantum and Reversible Computation, Stony Brook NY, May 28-31, 2003

## Acknowledgements. —

The author would like to thank V.E. Korepin for an invitation to speak at the Simons Conference. This work was supported by the Defense Advanced Research Projects Agency and by the U.S. Army Research Office.


text:
# The End of the MACHO Era: Limits on Halo Dark Matter From Stellar Halo Wide Binaries

## 1. Introduction

After formation, wide binaries retain their original orbital parameters except in so far as they are affected by gravitational encounters, or either one or both binary members evolve off the main sequence. These systems can therefore be used to probe inhomogeneities of Galactic potential that may be due to black holes, low luminosity stars, molecular clouds, or other objects, because their low binding energies are easily overcome by gravitational perturbations (Heggie, 1975).

Bahcall, Hut & Tremaine (1985) were the first to apply this principle, using wide binaries in the Galactic disk to investigate disk dark matter in the Solar neighborhood. Weinberg, Shapiro & Wasserman (1987) refined this approach by incorporating several effects that were previously ignored. In particular, they showed that gravitational encounters do not in general induce a sharp cutoff in the binary semi-major axis distribution, but rather generate a break in its power law. Although there have been extensive investigations of disk binary systems, they have not yielded strong conclusions about the disk potential. This is partly due to the relatively small size of the disk binary samples available, partly to the fact that they do not have good sensitivity for separations $`0.1`$ pc, and partly because of the intrinsic complexity of the disk potential. In principle, halo binaries could be used to search for dark matter, but this has never been done, primarily because there were no halo-binary samples adequate to the task.

However, halo dark matter has been investigated by several other techniques. Lacey & Ostriker (1985) suggested a mechanism for disk heating by supermassive black holes and discussed that black holes with $`M>10^6\mathrm{M}_{}`$ could destroy the disk. Null results of a search for “echoes” of gamma-ray bursts induced by gravitational lensing constrained halo dark matter in the mass range $`M10^{6.5}10^{8.5}\mathrm{M}_{}`$ (Nemiroff et al., 1993). Moore (1993) argued that massive black holes $`M>10^3\mathrm{M}_{}`$ could disrupt low-mass globular clusters, which would imply an upper limit $`M<10^3\mathrm{M}_{}`$. However, this argument is somewhat sensitive to assumptions about the initial population of globular clusters. Finally, microlensing experiments by the MACHO collaboration (Alcock et al., 2001), the EROS collaboration (Afonso et al., 2003), and the two collaborations working together (Alcock et al., 1998) found that MAssive Compact Halo Objects (MACHOs) with $`10^{7.5}\mathrm{M}_{}<M<30\mathrm{M}_{}`$ cannot account for the mass of the dark halo. De Rújula, Jetzer & Massó (1992) argued that baryonic dark matter with $`M<10^7\mathrm{M}_{}`$ would have evaporated away in a Galactic time scale.

The publication of the $`\mu >0.\mathrm{}18\mathrm{yr}^1`$ proper motion limited New Luyten Two Tenths (NLTT) catalog (Luyten, 1979; Luyten & Hughes, 1980) has vastly increased the pool of available data on binary systems, but the short color baseline of the photographic photometry in NLTT rendered halo/disk discrimination extremely difficult. However, Gould & Salim (2003) and Salim & Gould (2003) revised NLTT (rNLTT) with improved astrometry and photometry for the 44% of the sky covered by the intersection of the Second Incremental Release of the Two Micron All Sky Survey and the first Palomar Observatory Sky Survey. Chanamé & Gould (2003, hereafter CG) then constructed a homogeneous catalog of binaries from rNLTT and classified each entry as either disk or halo.

In this paper, we use Monte Carlo simulations to evaluate the effect of MACHOs on the halo binary distribution and compare these predictions with the observed halo binary sample of CG by means of a likelihood analysis. We briefly describe our CG halo binary sample in § 2 and make simple analytic estimates of the effects of perturbers in § 3. Detailed justifications of our assumptions and our Monte Carlo algorithm are presented in § 4. We present the results of our numerical simulations in § 5 and the likelihoods of halo dark-matter models in § 6. Finally, we summarize our results in § 7.

## 2. The CG Halo Wide Binaries

Figure 1 shows the halo binaries from the CG catalog upon which the current work is based. This binary distribution function is well approximated by a single power law from $`\mathrm{\Delta }\theta =3.\mathrm{}5`$ to the catalog limit at $`\mathrm{\Delta }\theta =900\mathrm{}`$. However, since the lower threshold of completeness is not precisely established, we consider several different lower limits. There are respectively (90, 68, 59, 47) binaries in the samples with lower limits ($`3.\mathrm{}5`$, $`5.\mathrm{}5`$, $`7.\mathrm{}5`$, $`10.\mathrm{}5`$). For $`\mathrm{\Delta }\theta 3.\mathrm{}5`$, the sample is not complete, and the deviation from a power-law distribution is clearly shown as triangles in Figure 1. However, since the fit is insensitive to the choice of lower limit ($`\mathrm{\Delta }\theta 3.\mathrm{}5`$), we adopt a lower threshold of $`\mathrm{\Delta }\theta =3.\mathrm{}5`$. Since severe incompleteness clearly sets in at $`\mathrm{\Delta }\theta <3.\mathrm{}5`$, one might be concerned about more modest incompleteness just above our adopted boundary $`\mathrm{\Delta }\theta =3.\mathrm{}5`$. However, as we discuss in § 6, incompleteness near $`\mathrm{\Delta }\theta 3.\mathrm{}5`$ would act to lessen the likelihood difference between a halo dark-matter model and the observations. Therefore, our procedure is conservative.

Assuming that the initial binary distribution is characterized by a power law, the slope of which might be different from the final observed one, the observed binary distribution can yield significant constraints on halo dark matter. Since binary systems with large semi-major axes, $`a0.1`$ pc, are more liable to be disturbed by perturbers than those with small $`a`$, the deviation (or lack thereof) from a power-law distribution at wide angular separations provides a test of halo dark-matter models.

## 3. Theoretical Expectations

We begin by considering two extreme regimes; the tidal regime ($`b_{\mathrm{min}}a`$) and the Coulomb regime ($`b_{\mathrm{min}}a`$), where $`a`$ is the binary semi-major axis and,

$`b_{\mathrm{min}}`$ $`=`$ $`\left({\displaystyle \frac{M}{\pi \rho vT}}\right)^{1/2}`$ (1)
$`=`$ $`700\mathrm{AU}\left({\displaystyle \frac{M}{\mathrm{M}_{}}}\right)^{1/2}\left({\displaystyle \frac{\rho }{\rho _H}}\right)^{1/2}\left({\displaystyle \frac{v}{300\mathrm{km}\mathrm{s}^1}}\right)^{1/2}`$
is the typical minimum impact parameter for perturbers of mass $`M`$, density $`\rho `$, and velocity $`v`$ over the duration $`T`$ that the binary is subjected to perturbations. We evaluate equation (1) using the standard local halo density $`\rho _H=0.009\mathrm{M}_{}\mathrm{pc}^3`$ (Bahcall & Soneira, 1980). We adopt $`T=10`$ Gyr. Note that while halo stars (and so halo binaries) are several Gyr older than this, $`T`$ represents the time since the binaries were assembled into the Galactic halo and not the time since their formation.

In the tidal regime, the perturbations are dominated by the single closest encounter, which yields a change of relative velocity,

$$\mathrm{\Delta }\text{v}\frac{2GM}{b^2v}a.$$
(2)
By evaluating $`\mathrm{\Delta }\text{v}`$ at $`b=b_{\mathrm{min}}`$ and equating the result with the internal velocity of the binary system, $`v^2=Gm/a`$, we obtain an estimate of transition separation,

$$a_t=\left(\frac{m}{4\pi ^2G\rho ^2T^2}\right)^{1/3}18,000\mathrm{AU}\left(\frac{\rho }{\rho _H}\right)^{2/3},$$
(3)
where we adopt $`m=1\mathrm{M}_{}`$ as the typical total mass of the binary system. The corresponding transition orbital period is $`P_t2\mathrm{Myr}(\rho /\rho _H)^1`$.

Binary systems whose semi-major axes are less than $`a_t`$ or equivalently whose orbital periods are shorter than $`P_t`$ are insensitive to perturbations. Note that the transition separation in the tidal regime is independent of the mass and velocity of perturbers. This regime approximately applies when $`b_{\mathrm{min}}a_t`$, i.e.,

$$M\left(\frac{m^2v^3}{16\pi G^2\rho T}\right)^{1/3}700\mathrm{M}_{}.$$
(4)
In contrast to the tidal regime, which is dominated by the single closest encounter, the perturbations in the Coulomb regime (named for the Coulomb logarithm, $`\mathrm{ln}\mathrm{\Lambda }\mathrm{ln}(a/b_{\mathrm{min}})`$), are described by continuous weak gravitational encounters,

$$(\mathrm{\Delta }v)^2=2_{b_{\mathrm{min}}}^a2\pi bdb\left(\frac{2GM}{bv}\right)^2\frac{\rho }{M}vT=\frac{16\pi G^2\rho MT}{v}\mathrm{ln}\mathrm{\Lambda },$$
(5)
where the factor 2 in front accounts for the independent perturbations of each component of the binary. The transition separation is then,

$$a_t=\frac{mv}{16\pi G\rho MT\mathrm{ln}\mathrm{\Lambda }}=\frac{3000\mathrm{AU}}{\mathrm{ln}\mathrm{\Lambda }}\left(\frac{M}{10^3\mathrm{M}_{}}\right)^1\left(\frac{\rho }{\rho _H}\right)^1,$$
(6)
and so is roughly inversely proportional to the perturber mass.

## 4. Numerical Simulation

In this section, we summarize our Monte Carlo algorithm, which evaluates the effect of all perturbations using a simple impulse approximation and ignores such effects as large-scale tides and molecular clouds. We include the contribution from the ionized binaries, although as we discuss in § 4.4 and show more fully in the Appendix, disrupted binaries diffuse to separations well beyond the observational range of interests in a time very short compared to $`T`$.

Since previous work, notably that of Weinberg, Shapiro & Wasserman (1987), has focused considerable effort on incorporating the above-mentioned effects, we first justify our decision to ignore them. It is primarily the fact that we are considering halo binaries whereas previous workers were investigating disk binaries that accounts for the difference in importance of these effects.

### 4.1. Impulse Approximation

In the Coulomb regime, the biggest relevant impact parameter is of order $`a`$. Since the perturbers are moving much faster than the binary components, the impulse approximation well describes encounters.

In the tidal regime, the perturbations are dominated by the closest single encounter, whose crossing time is,

$`T_c`$ $`=`$ $`{\displaystyle \frac{2b_{\mathrm{min}}}{v}}`$ (7)
$`=`$ $`20\mathrm{yr}\left({\displaystyle \frac{M}{\mathrm{M}_{}}}\right)^{1/2}\left({\displaystyle \frac{\rho }{\rho _H}}\right)^{1/2}\left({\displaystyle \frac{v}{300\mathrm{km}\mathrm{s}^1}}\right)^{3/2}.`$
Since, for masses $`M10^8M_{}`$, the crossing time is significantly less than the transition orbital period $`P_t2\mathrm{Myr}`$, the impulse approximation is always justified.

### 4.2. Disk and Galactic Tides

When halo binaries pass through the plane of Galactic disk, differential gravitational attractions give rise to tidal effects, also known as disk shocking (Binney & Tremaine, 1987). Since the mean velocity at which halo binaries cross the disk is $`v_z100\mathrm{km}\mathrm{s}^1`$ (Popowski & Gould, 1998a, b), the change in the z-component of velocity of halo binaries is $`\mathrm{\Delta }v_z4\pi G\mathrm{\Sigma }a/v_z`$ where $`\mathrm{\Sigma }=40\mathrm{M}_{}\mathrm{pc}^2`$ is the surface density of the Galactic disk in the Solar neighborhood (Zheng et al. 2001, and references therein). Equating this to the internal velocity of the binary system, $`v^2=Gm/a`$, yields,

$$a_{\mathrm{crit}}=\left(\frac{mv_z^2}{16\pi ^2G\mathrm{\Sigma }^2}\right)^{1/3}2\mathrm{pc},(\mathrm{Disk}\mathrm{Tides}).$$
(8)
For the most favorable (i.e., prograde) orbits, the binary system can be disrupted by Galactic tidal fields if the internal orbital period is longer than a Galactic year, $`P_G230\mathrm{Myr}`$. The critical semi-major axis is therefore,

$$a_{\mathrm{crit}}=\mathrm{AU}\left(\frac{m}{\mathrm{M}_{}}\frac{P_G^2}{\mathrm{yr}^2}\right)^{1/3}2\mathrm{pc},(\mathrm{Galactic}\mathrm{Tides}).$$
(9)
Since these critical values are not much larger than the largest projected separation in our sample, one might be concerned that the effects of disk and Galactic tides are non-negligible. However, the observational data show that the binary distribution is well described by a single power law within the limits of $`3.\mathrm{}5<\mathrm{\Delta }\theta <900\mathrm{}`$, and we are being conservative to ignore both tides because incorporating their effects would only serve to place more stringent constraints on the perturbers. On the other hand, if the binary distribution had shown a power-law break, which would have been indicative of dark matter, it would have been necessary to investigate how much of this signature was actually due to tides.

### 4.3. Giant Molecular Clouds

With typical masses $`10^510^6\mathrm{M}_{}`$ and surface density $`\mathrm{\Sigma }5\mathrm{M}_{}\mathrm{pc}^2`$, corresponding to a mean density $`\rho 6\times 10^4\mathrm{M}_{}\mathrm{pc}^3`$ over typical halo orbits, giant molecular clouds (GMCs) are clearly in the tidal regime. Because of their low density $`\rho 0.07\rho _H`$, they would add very modestly to the effect of halo tidal perturbers. The effect of GMCs is further diminished by the fact that $`b_{\mathrm{min}}`$ (eq. ) is substantially smaller than a typical cloud (and this remains so even if one considers the individual “clumps” inside a GMC). Hence, we ignore GMCs. As with tides, the effect of doing so is conservative: taking account of GMCs would only strengthen any limits we obtain.

### 4.4. Ionized Binaries

Even though binary systems are disrupted by perturbers, the stars of former binary members escape the system with extremely low escape energies due to the interplay of diffusive processes and the tidal boundary rather than disappearing immediately after disruption. As shown by the spectacular images of the globular cluster Palomar 5 suffering on-going tidal disruptions (Odenkirchen et al., 2003), the mass in a tidal tail can be larger than the mass in the parent body. By analogy, ionized binaries could in principle stay well inside the observational range of interest.

However, globular cluster escapees are no longer perturbed and so proceed on deterministic orbits that separate from the cluster at a roughly uniform rate, while the binary members drift away in the Galactic potential and continue to be heated by the omnipresent perturbers after they are well enough separated to be free from each others’ gravitational influence.

In the Appendix, we quantify this argument and analytically show that ionized binaries do not contribute significantly to the observed distribution. We also give a prescription for including them in our simulations, and indeed we find that they contribute negligibly.

### 4.5. Non-Circular Orbits

In our simulations, we assume that the binaries are subjected to perturbations by ambient objects whose density is constant in time. Strictly speaking this would apply only to binaries in circular Galactic orbits. More typical binaries would tend to spend the majority of their time farther from the Galactic center than the Sun, where the perturber density is lower. However, we find by numerically integrating a wide range of orbits, that the average density of perturbers encountered by binaries in non-circular orbit is almost always higher than for those in circular orbits. Hence, our approximation of constant perturber density is conservative.

### 4.6. Monte Carlo Algorithm

To test the consistency of the observed binary distribution with the presence of massive perturbers, we must be able to simulate the effect of these perturbers on arbitrary power-law initial distributions. Our basic method for doing this is Monte Carlo simulations in which the initial binary semi-major axis is drawn uniformly from power-law distribution over the interval, $`1<\mathrm{log}(a/\mathrm{AU})<5.5`$. The eccentricity is drawn randomly from a distribution uniform in $`e^2`$ (Weinberg, Shapiro & Wasserman, 1987), and the phase and orientation of the orbit are assigned randomly. The perturbers are assumed to have an isotropic Maxwellian velocity distribution relative to the binary center of mass, with a one-dimensional dispersion $`\sigma =200\mathrm{km}\mathrm{s}^1`$. This reproduces the true rms velocity, which is a combination of an isothermal-sphere perturber distribution of circular speed $`v_c=\sqrt{2}\sigma =220\mathrm{km}\mathrm{s}^1`$, and the measured dispersions of halo stars111Chiba & Beers (2000) obtained dispersions of halo stars $`(\sigma _\pi ,\sigma _\theta ,\sigma _z)=`$ $`(141,106,95)\mathrm{km}\mathrm{s}^1`$. The corresponding one-dimensional dispersion for perturbers is $`\sigma =194\mathrm{km}\mathrm{s}^1`$, which is therefore insensitive to the precise choice of $`(\sigma _\pi ,\sigma _\theta ,\sigma _z)`$. $`(\sigma _\pi ,\sigma _\theta ,\sigma _z)=`$ $`(170,97,93)\mathrm{km}\mathrm{s}^1`$ (Popowski & Gould, 1998a, b). In principle, the velocity distribution should be treated as anisotropic. However, this is a higher-order effect, which we ignore in the interest of simplicity.

We consider all impact parameters with $`bb_{\mathrm{max}}=\mathrm{max}\{10b_{\mathrm{min}},2a\}`$. In particular, we evaluate the $`100`$ closest impacts in the tidal regime, even though as we discuss in § 3, the single closest encounter dominates. The mass of each binary component is set to be $`0.5\mathrm{M}_{}`$. We evolve 100,000 binary systems in each simulation. To illustrate the dependence on the mass of perturbers, we show their effects on an artificial initial binary distribution that is independent of semi-major axis (see Fig. 2). As expected, for the binary systems that are initially tightly bound, the final distributions are almost the same as the initial ones regardless of the mass of the perturbers. However, at wide separations, the distributions are driven to a new power law, which becomes steeper with increasing mass.

This approach works well for any individual initial power-law distribution, but is too time-consuming to process the very large number of distributions required for the comparison of data and models. In § 5.5, we introduce a scattering-matrix formalism that substantially improves the efficiency of mass-production Monte Carlo simulations.

### 4.7. Fitting Formula

Motivated by the fact that the final distribution is well-approximated by power laws at each extreme, we use a two-line (five-parameter) fitting formula given by,

$$H(x)=\left[f(x)^n+g(x)^n\right]^{1/n},$$
(10)
where $`f(x)`$ and $`g(x)`$ are (two-parameter) straight lines in their argument, $`x=\mathrm{log}(a/\mathrm{AU})`$, each corresponding to the respective asymptotic behaviors of $`H(x)`$. The fifth parameter $`n`$ permits a smooth transition between $`f(x)`$ and $`g(x)`$ in the intermediate region. We calculate the five parameters of equation (10) by minimizing $`\chi ^2`$ for a given data set, and the fitting curves for four different masses of perturber are represented as dashed lines in Figure 2.

## 5. Results

Here, we present our main calculations on the evolution of binary distributions with various initial slopes under the influence of various perturber masses and halo densities, and we evaluate the transition separations. Although the semi-major axis of a binary system is a direct indicator of the binding energy of the system and is the theoretically most tractable quantity, it is not observable. It is the angular separations on the sky that we can directly measure from observations. To compare our results with the data, we calculate physical separations projected on the sky plane and convolve these with an adopted distance distribution to predict the binary distributions as a function of angular separation.

In principle, one could compare models directly to the observed projected separations since CG give individual distance estimates to each binary. However, while the observational selection function is quite simple for angular separations (essentially just a pair of $`\mathrm{\Theta }`$-functions), it is rather complex and would be extremely difficult to model for projected physical separations. Hence, we compare our models to the most directly observed quantity: angular separations.

### 5.1. Semi-Major Axis, $`a`$

We begin the simulation with a flat distribution, $`dN/d\mathrm{log}aconst`$ for calculation of scattering matrices. We then investigate the dependence of the resulting final binary distribution on perturber mass and halo density. Figure 3 shows a sample of our results with the initial power-law distribution, $`dN/d\mathrm{log}aa^{0.567}`$, in agreement with the observations as summarized in § 2. Two models with the same perturber mass but widely different halo densities, $`\rho _H`$ and $`10^1\rho _H`$ are shown. The distributions of semi-major axis are represented by circles.

### 5.2. Projected Physical Separation, $`r_{}`$

At $`T=10`$ Gyr, we calculate the projected physical separation, $`r_{}`$ of each binary taking account of the randomly chosen orientation to the line of sight and the orbital phase. The distributions of $`r_{}`$ are shown by triangles in Figure 3.

The two sets of distributions are very similar to each other as demonstrated by the solid and the dashed lines, which are fitted to the points. This can be understood as follows. The time-averaged physical separation is

$$r(e)=\frac{1}{2\pi }_0^{2\pi }a(1e\mathrm{cos}\psi )^2d\psi =a\left(1+\frac{e^2}{2}\right).$$
(11)
Averaging over the uniform distribution in $`e^2`$ yields $`r=5a/4`$, and projecting onto the sky plane gives,

$$r_{}=\frac{\pi }{4}r=\frac{5\pi }{16}a0.98a.$$
(12)
That is, the triangles are, on average, slightly shifted to left of the circles in Figure 3.

There is, however, a deviation from equation (12) at large projected physical separation. We find that a binary system with a highly eccentric orbit is more likely to be disturbed than one with a circular orbit. Since the highly eccentric systems are selectively destroyed for loosely bound binary systems, the distribution in eccentricity is no longer isotropic, and the actual projected physical separations differ from equation (12). We therefore use numerical calculations of the projected physical separations rather than the analytic estimate of equation (12).

### 5.3. Transition Separation, $`a_t`$

Using the five-parameter fitting formula from equation (10), we calculate the asymptotic slopes for distributions in semi-major axis. These are represented as solid curves in Figure 3. The intersection of these asymptotic slopes characterizes the transition $`a_t`$ from the unperturbed to the perturbed regime. The two arrows in Figure 3 represent the intersections for two different halo densities. As predicted by equation (3), the transition for a denser halo model occurs at smaller separation, $`a_t\rho ^{2/3}`$.

The transition separations for various perturber mass and halo densities are shown in Figure 4, which provides a qualitative understanding of the underlying physics of binary disruption and is in a good accord with the prediction of equation (3) and (6). In particular, in the tidal regime the transition separation is independent of perturber mass, and is a function of halo density $`a_t\rho ^{2/3}`$, although the actual values of $`a_t`$ are a factor of $`2`$ smaller than those predicted by equation (3). At smaller mass, the transition separation grows roughly as $`a_tM^1`$, as predicted by equation (6) for the Coulomb regime.

### 5.4. Angular Separation, $`\mathrm{\Delta }\theta `$

To predict the observed angular distribution function, we must convolve the model’s projected-separation distribution function with an assumed distribution of (inverse) distances to the binaries in the sample.

We derive this distribution from the observed distances of the 90 binaries in the CG sample. To understand our procedure for doing so, consider first the idealized case of binaries drawn from a distance-limited sample of uniform density. Of course, in this case, the distance distribution would be $`dN/d\mathrm{log}RR^3\mathrm{\Theta }(R_{\mathrm{max}}R)`$, where $`\mathrm{\Theta }`$ is a step function, and $`R_{\mathrm{max}}`$ is the distance limit. However, the binaries selected at fixed $`\mathrm{\Delta }\theta `$ would not be distributed as $`R^3`$. This is because the binaries at distances $`R_1`$ and $`R_2`$ would have projected separations $`r_{,1}=R_1\mathrm{\Delta }\theta `$ and $`r_{,2}=R_2\mathrm{\Delta }\theta `$, and these differ in frequency (for fixed log intervals) in the ratio $`(r_{,1}/r_{,2})^{1\alpha }`$, where $`\alpha =1.567`$ is the binary-separation power-law slope. Hence, the sample distribution would be $`dN_{\mathrm{sample}}/d\mathrm{log}RR^{4\alpha }\mathrm{\Theta }(R_{\mathrm{max}}R)`$. Therefore, to obtain an estimate of the underlying distribution of distances that yields the observed sample distribution, we simply adopt the observed distribution, but weight each binary by $`R^{\alpha 1}`$. We estimate $`R`$ for each binary by applying the fourth-order color-magnitude relation of CG to the brighter component, except when that component has no $`J`$-band data, in which case we use the fainter component.

### 5.5. Scattering Matrix, $`𝐒`$

To simultaneously investigate large subsets of initial power-law distributions for each model specified by $`M`$ and $`\rho `$, we use a Monte Carlo simulation (see § 4.6) with an initially flat ($`dN/d\mathrm{log}aconst`$) distribution to construct a scattering matrix, $`𝐒_{ij}(M,\rho )`$, the relative probability that a binary with initial log semi-major axis, $`\mathrm{log}a_j`$ will finally be observed in the angular separation bin, $`\mathrm{log}\mathrm{\Delta }\theta _i`$. Explicitly,

$$𝐒_{ij}(M,\rho )=\underset{k}{}\mathrm{d}\mathrm{log}R\frac{dN}{d\mathrm{log}R}\delta \left[\mathrm{log}R\mathrm{\Delta }\theta _i\mathrm{log}r_{,j}(k;M,\rho )\right],$$
(13)
where $`r_{,j}(k;M,\rho )`$ is the final projected separation of a binary with initial semi-major axis $`a_j`$ in $`k^{th}`$ realization of a Monte Carlo simulation of perturbers with mass $`M`$ and density $`\rho `$. In practice, since the radial profile $`dN/d\mathrm{log}R`$ is estimated directly from the binary data as discussed in § 5.4, we evaluate equation (13) by discretely summing over the $`l=1\mathrm{}90`$ binaries,

$`𝐒_{ij}`$ $`(M,\rho )`$ $`={\displaystyle \underset{k}{}}{\displaystyle \underset{l=1}{\overset{90}{}}}R_l^{\alpha 1}`$
$`\times \mathrm{\Theta }\left({\displaystyle \frac{\delta \mathrm{log}\mathrm{\Delta }\theta }{2}}|\mathrm{log}R_l\mathrm{\Delta }\theta _i\mathrm{log}r_{,j}(k;M,\rho )|\right),`$
where $`\delta \mathrm{log}\mathrm{\Delta }\theta `$ is the logarithmic width of the $`\mathrm{\Delta }\theta `$ bins.

Since the simulations are calculated from a flat distribution, i.e., equal number of binaries per logarithmic semi-major axis bin, the un-normalized probability of finding binaries at $`\mathrm{\Delta }\theta `$ for a given dark-matter model is then given by,

$$P(\mathrm{\Delta }\theta _i;M,\rho ,\alpha )=𝐒_{ij}(M,\rho )w(a_j;\alpha ),$$
(15)
where $`w(a_j;\alpha )=a_j^{1\alpha }`$ is the power-law of the initial binary distribution, $`dN/daa^\alpha `$.

## 6. Dark Matter Limits

For each model, we fix the normalization so that the number expected in the interval $`3.\mathrm{}5<\mathrm{\Delta }\theta <900\mathrm{}`$ is equal to 90, i.e., the total number observed over this interval. We call the normalized resulting function $`P_\mathrm{N}(\mathrm{\Delta }\theta ;M,\rho ,\alpha )`$, evaluate the likelihood of a given model as a sum over the 90 observed binaries,

$$\mathrm{ln}L(M,\rho ,\alpha )=\underset{i=1}{\overset{90}{}}\mathrm{ln}P_\mathrm{N}(\mathrm{\Delta }\theta _i;M,\rho ,\alpha ),$$
(16)
and find the maximum likelihood over various initial slopes $`\alpha `$,

$$\mathrm{ln}L(M,\rho )\underset{\alpha }{\mathrm{max}}\left\{\mathrm{ln}L(M,\rho ,\alpha )\right\}.$$
(17)
Figure 5 shows the normalized resulting functions for various perturber masses with $`\rho =\rho _\mathrm{H}`$. For 1000 $`M_{}`$ perturbers, the likelihood is maximized by an initial distribution that is nearly flat. However, no initial power-law is consistent with the observed distribution. For low mass perturbers, by contrast, the likelihood is maximized by an initial slope that is only slightly shallower than the observed one, and the resulting distribution is consistent with observations.

We compare the resulting likelihoods to $`L(0,0)`$, i.e., the model with no dark-matter perturbers, which is characterized by the best fit slope $`\alpha =1.567`$, and we define confidence levels by

$$\sigma (M,\rho )=\sqrt{2\left[\mathrm{ln}L(0,0)\mathrm{ln}L(M,\rho )\right]}.$$
(18)
As we discussed in § 2, the fact that the binary sample is clearly incomplete for $`\mathrm{\Delta }\theta <3.\mathrm{}5`$ implies that it may also be slightly incomplete for binaries just above this threshold. However, one can see from the form of the normalized resulting functions in Figure 5 that such incompleteness would actually favor models with massive perturbers over the no-perturber model. We verify numerically that this indeed is the sign of the bias. Hence our choice of the $`\mathrm{\Delta }\theta 3.\mathrm{}5`$ threshold is conservative in the face of possible incompleteness.

Figure 6 shows the resulting contour plot excluding various dark-matter models at various confidence ($`\sigma `$) levels. Finally, in Figure 7 we compare our 95% confidence ($`2\sigma `$) limits with those from the EROS (Afonso et al., 2003) and MACHO (Alcock et al., 2001) microlensing collaborations. The halo binaries exclude models that have generally higher MACHO masses than those probed by these microlensing experiments. Moreover, they extend all the way to the Lacey & Ostriker (1985) limit based on stability of the Galactic disk.

All three curves in Figure 7 show limits on models with $`\delta `$-function mass distributions. However any model with broader distribution of mass but total density $`\rho _H`$, would be ruled out provided that the combined limits from the three curves ruled out each mass separately.

Hence, the only remaining window open for MACHOs would be black holes with a mass function strongly peaked at $`M35\mathrm{M}_{}`$. Even this model is more strongly constrained than shown Figure 7 because the MACHO (Alcock et al., 2001) results and our results each separately place weak limits on this model, which when combined comes close to ruling it out.

## 7. Conclusions

In this paper, we have investigated the evolution of halo wide binaries in the presence of MACHOs and estimated upper limits of MACHO density as a function of their assumed mass by comparing our simulations to the sample of halo wide binaries of CG. We exclude MACHOs with masses $`M>43\mathrm{M}_{}`$ at the standard local halo density $`\rho _H`$ at the 95% confidence level.

MACHOs have been a major dark-matter candidate ever since observations first established that this mysterious substance dominates the mass of galaxies. Prodigious efforts over several decades have gradually whittled down the mass range allowed to this dark-matter candidate. However, the window for MACHOs with $`30\mathrm{M}_{}M10^3\mathrm{M}_{}`$ remained completely open while constraints in the range $`10^3\mathrm{M}_{}M10^6\mathrm{M}_{}`$ were somewhat model dependent. Our new limits on MACHOs $`M>43\mathrm{M}_{}`$ all but close this window.

We are grateful to Éric Aubourg and Kim Griest for providing data for Figure 7. We thank John Bahcall, Bohdan Paczyński, and especially Scott Tremaine for valuable comments that significantly improved the paper. A detailed critique by referee Terry Oswalt also greatly improved the paper. This work was supported by grant AST 02-01266 from the NSF.

## Appendix A Ionized Binaries

After disruption of a binary system, the two ionized members remain in similar Galactic orbit, so it is only the separation along the direction of the orbital motion that can keep increasing, while the perpendicular separation oscillates.

In the Coulomb regime, the average post-ionization gain in the relative velocity of binaries along the orbital direction is (see eq.),

$$\sqrt{v_{}^2}=\sqrt{\frac{32\pi G^2\rho M\mathrm{\Delta }t}{3v}\mathrm{ln}\mathrm{\Lambda }},$$
(A1)
where $`\mathrm{\Delta }t`$ is the remaining time to 10 Gyr after disruption. For a diffusive process that is uniform over time $`\mathrm{\Delta }t`$, the root-mean-square separation parallel to the orbital motion is

$$\sqrt{d_{}^2}=\sqrt{v_{}^2\frac{\mathrm{\Delta }t^2}{3}}=2000\mathrm{pc}\left(\frac{\rho }{\rho _H}\right)^{1/2}\left(\frac{M}{30M_{}}\right)^{1/2}\left(\frac{v}{300\mathrm{km}\mathrm{s}^1}\right)^{1/2}\left(\frac{\mathrm{ln}\mathrm{\Lambda }}{5.1}\right)^{1/2}\left(\frac{\mathrm{\Delta }t}{10\mathrm{Gyr}}\right)^{3/2},$$
(A2)
where the Coulomb logarithm is calculated at the scaled quantities.

For a conservative limit of the tidal radius, $`a_t=3`$ pc, the time required for ionized binaries to separate farther than $`a_t`$ is 0.13 Gyr, so that only binaries ionizing within the last $``$ 1% of the age of the Galaxy have a significant chance to be confused with bound systems, even for the lowest mass perturbers that we can effectively probe. In the tidal limit, binaries escape with characteristic velocities of the transition separation $`10^4`$ AU, i.e., 300 $`\mathrm{m}\mathrm{s}^1`$. Hence, they drift one tidal radius in only 10 Myr, so their impact is even smaller. Nevertheless, since there are only a handful of binaries in the widest-separation bins, it is important to make a careful estimate of the contribution from ionizing binaries. We take account of ionized binaries in the simulations as follows. Binaries are considered ionized either when they have positive energy or they have $`a>a_t`$. They are then assigned a relative velocity equal to their escape velocity in the former case, or zero in the latter. A random orbital direction is chosen. The ionized binaries in the Coulomb regime then continue to suffer perturbations to the end of the simulations whose effect we calculate using equation (A1) and (A2). At the end of the simulation, the binary is assigned a transverse separation drawn randomly from a sinusoidal distribution of amplitude,

$$d_{,\mathrm{MAX}}=\frac{v_{}}{\mathrm{\Omega }}=26\mathrm{pc}\left(\frac{v_{}}{\mathrm{km}\mathrm{s}^1}\right),$$
(A3)
where $`v_{}`$ is the final transverse velocity, $`\mathrm{\Omega }=\sqrt{2}v_c/R_0`$ is the epicyclic frequency, and $`R_0`$ is the Galactocentric distance. Finally, we “observe” the ionized binary from a random orientation and record the projected separation. We find that the ionized binaries have no significant effect either on the final distribution or the calculated likelihoods. We ignore ionized binaries in the tidal regime because these escape with much higher initial velocities and because these velocities grow much more rapidly than indicated by equation (A1).


text:
# Fast rotation of an ultra-cold Bose gas

(July 18, 2003)

## Abstract

We study the rotation of a 87Rb Bose-Einstein condensate confined in a magnetic trap to which a focused, off resonant gaussian laser beam is superimposed. The confining potential is well approximated by the sum of a quadratic and a quartic term which allows to increase the rotation frequency of the gas above the trap frequency. In this fast rotation regime we observe a dramatic change in the appearance of the quantum gas. The vortices which were easily detectable for a slower rotation become much less visible, and their surface density is well below the value expected for this rotation frequency domain. We discuss some possible tracks to account for this effect.

The fast rotation of a macroscopic quantum object often involves a dramatic change of its properties. For rotating liquid 4He, superfluidity is expected to disappear for large rotation frequencies Donnelly91 . A similar effect occurs for type-II superconductors placed in a magnetic field, which loose their superconductivity when the field exceeds a critical value Tinkham . For the last few years it has been possible to set gaseous Bose-Einstein condensates into rotation, by nucleating in the gas quantized vortices which arrange themselves in a triangular lattice Matthews99 ; Madison00 ; Ketterle1 ; Hodby01 ; Haljan01 . The theoretical investigation of the properties of these fast rotating gases has led to several possible scenarios. Depending on the confinement of the gas, its dimensionality and the strength of the interactions, theoretical predictions involve the nucleation of multiply charged vortices Fetter01 ; Lundh02 ; Kasamatsu02 ; Kavoulakis03 , melting of the vortex lattice Sinova02 ; Fisher03 , or effects closely connected to Quantum Hall physics Cooper01 ; Paredes01 ; Ho01 ; Regnault03 .

The confinement of cold atomic gases is generally provided by a magnetic trap which creates an isotropic transverse potential $`m\omega _{}^2r_{}^2/2`$ in the plane perpendicular to the rotation axis $`z`$ (we set $`x^2+y^2=r_{}^2`$). As the rotation frequency $`\mathrm{\Omega }`$ varies from $`0`$ to $`\omega _{}`$, the number of quantized vortices increases Raman01 . The limiting case of $`\mathrm{\Omega }=\omega _{}`$, e.g. fast rotation, is singular: the confinement vanishes since the trapping force is compensated by the centrifugal force. In classical terms, the only remaining force on the particle is the Coriolis force. In quantum terms this corresponds to a gauge field $`𝑨(𝒓)=m𝛀\times 𝒓/2`$. This is analog to the physics of charged particles in a uniform magnetic field and one expects to recover for weak interactions an energy spectrum with Landau levels separated by $`2\mathrm{}\omega _{}`$. However the absence of confinement in a pure harmonic trap rotating at $`\mathrm{\Omega }=\omega _{}`$ makes this study experimentally delicate Rosenbusch01 . The Boulder group recently reached $`\mathrm{\Omega }=0.99\omega _{}`$ using evaporative spin-up Boulder03 .

In the present paper, following suggestions in Fetter01 ; Lundh02 ; Kasamatsu02 ; Kavoulakis03  we study the rotation of a Bose-Einstein condensate in a trap whose potential is well approximated by a superposition of a quadratic and a small quartic potential:

$$V(𝒓)\frac{1}{2}m\omega _z^2z^2+\frac{1}{2}m\omega _{}^2r_{}^2+\frac{1}{4}kr_{}^4(k>0).$$
(1)
We can thus explore with no restriction the domain of rotation frequencies $`\mathrm{\Omega }`$ around $`\omega _{}`$. Our results show a strong dependance of the behavior of the gas with $`\mathrm{\Omega }`$. For $`\mathrm{\Omega }/\omega _{}<0.95`$ we observe a regular vortex lattice. This lattice gets disordered when $`\mathrm{\Omega }`$ increases. For $`\mathrm{\Omega }>\omega _{}`$ the number of detectable vortices is dramatically reduced although we have clear evidences that the gas is still ultra-cold and in fast rotation. We conclude this Letter by proposing some possible explanations for this behavior.

Our 87Rb condensate is formed by radio-frequency evaporation in a combined magnetic+laser trap. The pure magnetic trap provides a harmonic confinement along the 3 directions with the frequencies $`\omega _{}^{(0)}/2\pi =75.5`$ Hz and $`\omega _z/2\pi =11.0`$ Hz. We superimpose a blue detuned laser (wavelength 532 nm) propagating along the $`z`$ direction to provide the quartic term in the confinement. The potential created by the laser is

$$U(r_{})=U_0\mathrm{exp}\left(\frac{2r_{}^2}{w^2}\right)U_0\frac{2U_0}{w^2}r_{}^2+\frac{2U_0}{w^4}r_{}^4.$$
(2)
The laser’s waist is $`w=25\mu `$m and its power is 1.2 mW. The first term $`U_0`$ ($`k_B\times 90`$ nK) in the right hand side of (2) is a mere shift of the energy scale. The second term is a correction of the transverse trapping frequency; we get for the combined magnetic-laser trap $`\omega _{}/2\pi =65.6`$ Hz anisotropy . The last term in (2) corresponds to the desired quartic confinement, with $`k=2.6(3)\times 10^{11}`$ J.m-4. The expansion in (2) is valid for $`r_{}w/2`$, which is indeed the case for all data presented here.

We start the experimental sequence with a quasi-pure condensate that we stir using an additional laser beam, also propagating along $`z`$ Madison00 . This laser stirrer creates an anisotropic potential in the $`xy`$ plane which rotates at a frequency $`\mathrm{\Omega }_{\mathrm{stir}}`$. To bring the condensate in rotation at a frequency close to $`\omega _{}`$, we use two stirring phases. First we choose $`\mathrm{\Omega }_{\mathrm{stir}}^{(1)}\omega _{}/\sqrt{2}`$, so that the stirring laser resonantly excites the transverse quadrupole mode $`m=+2`$ of the condensate at rest Madison01 . The duration of this first excitation is 300 ms and we then let the condensate relax for 400 ms in the axisymmetric trap. This procedure sets the condensate in rotation, with a vortex lattice containing typically 15 vortices.

We then apply the laser stirrer for a second time period, now at a rotation frequency $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}`$ close to or even above $`\omega _{}`$. For the condensate with 15 vortices prepared during the first phase, the transverse quadrupole mode is shifted to a higher frequency and can now be resonantly excited by a stirrer rotating at $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}\omega _{}`$ Zambelli98 ; Haljan01 . The second stirring phase lasts for 600 ms. It is followed by a 500 ms period during which we let the condensate equilibrate again in the axisymmetric trap. At this stage the number of trapped atoms is $`3\times 10^5`$.

We find that this double nucleation procedure is a reliable way to produce large vortex arrays and to reach effective rotation frequencies $`\mathrm{\Omega }_{\mathrm{eff}}`$ around or above $`\omega _{}`$. In particular it is much more efficient than the procedure which consists in applying permanently the stirring laser while increasing continuously its frequency. It has one drawback however: the effective rotation frequency $`\mathrm{\Omega }_{\mathrm{eff}}`$ of the gas after equilibration might differ significantly from the rotation frequency $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}`$ that we apply during the second stirring phase. The role of the stirring phase is to inject angular momentum in the system, and the rotation frequency of the atom cloud may subsequently change if the atom distribution (hence the moment of inertia) evolves during the final equilibration phase.

The rotating atom cloud is then probed destructively by switching off the confining potential, letting the cloud expand during $`\tau =18`$ ms and performing absorption imaging. We take simultaneously two images of the atom cloud after time of flight (TOF). One imaging beam is parallel to the rotation axis $`z`$ and the other one is perpendicular to $`z`$ Rosenbusch02 .

Fig. 1a-h shows a series of images taken along the $`z`$ axis and obtained for various $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}`$ around $`\omega _{}`$. When $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}`$ increases, the increasing centrifugal potential weakens the transverse confinement. This leads to an increasing radius in the $`xy`$ plane. As it can be seen from the transverse and longitudinal density profiles of Fig. 2, obtained for $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}/2\pi =66`$ Hz, the gas after TOF has the shape of a flat pancake. The limit of our method for setting the gas in fast rotation is shown on the image (h) at the extreme right of Fig. 1. It corresponds to $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}/2\pi =69`$ Hz, far from the quadrupole resonance of the condensate after the first nucleation. In this case we could not bring the gas in the desired fast rotation regime.

For the pictures (e,f) obtained with $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}/2\pi =66`$ and $`67`$ Hz, the optical thickness of the cloud has a local minimum in the center; this indicates that the confining potential in the rotating frame, $`V_{\mathrm{rot}}(𝒓)=V(𝒓)m\mathrm{\Omega }_{\mathrm{eff}}^2r_{}^2/2`$, has a Mexican hat shape. The effective rotation frequency $`\mathrm{\Omega }_{\mathrm{eff}}`$ thus exceeds $`\omega _{}`$. The striking feature of these images is the small number of visible vortices which seems to conflict with a large value of $`\mathrm{\Omega }_{\mathrm{eff}}`$. The main goal of the remaining part of this Letter is to provide more information on this puzzling regime.

In order to analyze quantitatively the pictures of Fig. 1, we need to model the evolution of the cloud during TOF. For a pure harmonic potential ($`k=0`$), a generalization of the analysis of Castin96  to the case of a rotating condensate shows that the expansion in the $`xy`$ plane is well described by a scaling of the initial distribution by the factor $`(1+\omega _{}^2\tau ^2)^{1/2}`$ Chevy01 . We assume here that this is still approximately the case for a condensate prepared in a trap with a non zero quartic term $`kr_{}^4/4`$.

We have analyzed a series of 60 images such as those of Fig. 1, assuming an initial atomic distribution given by the Thomas-Fermi law:

$$n(𝒓)=\frac{m}{4\pi \mathrm{}^2a}\left(\mu V_{\mathrm{rot}}(𝒓)\right),$$
(3)
where $`a=5.2`$ nm is the scattering length characterizing the atomic interactions and $`\mu `$ the chemical potential. The optical thickness for the imaging beam propagating along $`z`$, proportional to the column atomic density in the $`xy`$ plane, varies as $`(\alpha +\beta r_{}^2+\gamma r_{}^4)^{3/2}`$, where $`\alpha ,\beta ,\gamma `$ can be expressed in terms of the physical parameters of the problem imagetrans . An example of the fit, which takes $`\mu `$ and $`\mathrm{\Omega }_{\mathrm{eff}}`$ as adjustable parameters, is given in Fig. 2a. The agreement is correct, though not as good as for condensates confined in purely harmonic traps. This may be a consequence of the approximate character of the scaling transform that we use to describe the TOF evolution. The resulting values for $`\mathrm{\Omega }_{\mathrm{eff}}`$ as a function of $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}`$ are given in Fig. 3. We find $`\mathrm{\Omega }_{\mathrm{eff}}\mathrm{\Omega }_{\mathrm{stir}}^{(2)}`$ for stirring frequencies $`68`$ Hz. Note that for $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}/2\pi =68`$ Hz (image Fig. 1g) the quality of the fit is comparatively poor due to local inhomogeneities of the atom cloud.

From the value of $`\mu `$ given by the fit, we recover the atom number Fetter01 . For the largest measured rotation frequency, the Thomas-Fermi distribution (3) corresponds to a nearly spherical atom cloud before TOF (diameter 29 $`\mu `$m in the $`xy`$ plane and length 34 $`\mu `$m along $`z`$).

A second determination of the effective rotation frequency $`\mathrm{\Omega }_{\mathrm{eff}}`$ of the condensate is provided by the vortex surface density after TOF. Assuming that the vortex pattern is scaled by the same factor as the condensate density, we deduce the vortex density $`\rho _v`$ before TOF, hence the rotation frequency $`\mathrm{\Omega }_{\mathrm{eff}}=\pi \mathrm{}\rho _v/m`$ Raman01 . For $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}<\omega _{}`$, the value of $`\mathrm{\Omega }_{\mathrm{eff}}`$ deduced in this way and plotted in Fig. 3 is in fair agreement with the one deduced from the fit of the images. On the contrary, for $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}/2\pi =6668`$ Hz the number of distinguishable vortices is much too low to account for the rotation frequencies determined from the fits of the TOF images.

In order to gather more information on the rotational properties of the gas, we now study the two transverse quadrupole modes $`m=\pm 2`$ of the gas. We recall that these two modes have the same frequency for a non rotating gas, due to symmetry. For a rotating gas in the hydrodynamic regime, the frequency difference $`\omega _{+2}\omega _2`$ is proportional to the ratio between the average angular momentum of the gas and its moment of inertia, which is nothing but the desired effective rotation frequency $`\mathrm{\Omega }_{\mathrm{eff}}`$ Zambelli98 . The strength of this approach is that it does not make any assumption on the expansion during TOF.

To study these transverse quadrupole modes, we briefly illuminate the rotating gas using the laser stirrer, now with fixed axes Chevy00 . This excites a superposition of the modes $`m=\pm 2`$ with equal amplitudes. We then let the cloud evolve freely in the trap for an adjustable duration and we perform the TOF analysis. From the time variation of both the ellipticity of the cloud in the $`xy`$ plane and the inclination of its eigenaxes, we deduce $`\omega _{\pm 2}`$, hence the effective rotation frequency $`\mathrm{\Omega }_{\mathrm{eff}}=(\omega _{+2}\omega _2)/2`$. The corresponding results are plotted in Fig. 3. They are in good agreement with the results obtained from the fits of the images: $`\mathrm{\Omega }_{\mathrm{eff}}\mathrm{\Omega }_{\mathrm{stir}}^{(2)}`$. Consequently they conflict with the rotation frequency that one derives from the vortex surface density when $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}/2\pi =6668`$ Hz.

To explain the absence of visible vortices in the regime $`\mathrm{\Omega }_{\mathrm{stir}}^{(2)}/2\pi =6668`$ Hz, one could argue that the gas might be relatively hot, hence described by classical physics. The equilibrium state should then correspond to rigid body rotation. However all experiments shown here are performed in presence of radio-frequency (rf) evaporative cooling. The rf is set 24 kHz above the value which empties the trap and it eliminates all atoms crossing the horizontal plane located at a distance $`x_{\mathrm{ev}}19\mu `$m below the center (one-dimensional evaporation). ¿From evaporative cooling theory Ketterle96 ; Luiten96  we know that the equilibrium temperature $`T`$ of the rotating gas is a fraction of the evaporation threshold $`V_{\mathrm{rot}}(x_{\mathrm{ev}})32`$ nK for $`\mathrm{\Omega }_{\mathrm{eff}}/2\pi =67`$ Hz gaussian . This indicates that $`T`$ is in the range of $`515`$ nK, well below the critical temperature at this density (180 nK for an estimated density $`3\times 10^{13}`$ cm-3) evaporation . This low temperature is also confirmed by the small decay rates of the quadrupole modes ($`20`$ s-1), characteristic of $`TT_c`$.

The estimated temperature $`T`$ is of the order of the splitting between Landau levels $`2\mathrm{}\omega _{}/k_B=6.3`$ nK, so that only the first two or three levels are appreciably populated. For each Landau level, the rf evaporation eliminates states with an angular momentum $`L_z/\mathrm{}>(x_{\mathrm{ev}}/a_{\mathrm{ho}})^2200`$, where $`a_{\mathrm{ho}}=\sqrt{\mathrm{}/(m\omega _{})}`$ is the ground state’s size of the harmonic oscillator with frequency $`\omega _{}`$.

A first possibility to interpret our data consists in assuming that the whole gas shown in Fig. 1e-g is rotating, that it is in the degenerate regime, and that it cannot be described by a single macroscopic wave function since the number of distinguishable vortices is too low to account for the measured rotation frequency. In this point of view, the state of the system is therefore more complex than a single Hartree state. A theoretical analysis along this line, involving the formation of composite bosons, has been very recently proposed in Akkermans03 .

A second possible explanation is that the vortices are still present for $`\mathrm{\Omega }\omega _{}`$, but that the vortex lines are strongly tilted or bent, so that they do not appear as clear density dips in the images of Fig. 1e-g. This could be for example a consequence of the residual anisotropy of the trapping potential in the $`xy`$ plane anisotropy . In an experiment performed at $`\mathrm{\Omega }=0.95\omega _{}`$, the Boulder group has shown that a static anisotropy of $`4`$% reduces considerably the vortex visibility over a time scale of $`0.8`$ s Engels02 . However the residual anisotropy of our trap is $`<1`$% and this effect should be limited, unless the vortex lattice becomes notably more fragile as $`\mathrm{\Omega }`$ approaches $`\omega _{}`$. The distortion of the vortex lattice could also result from a strong increase of the cristallization time for $`\mathrm{\Omega }\omega _{}`$ which would then exceed the 500 ms equilibration time $`\tau _{\mathrm{eq}}`$. To check this effect, we have increased $`\tau _{\mathrm{eq}}`$ up to 2 s without noticing any qualitative change in the atom distribution (for larger $`\tau _{\mathrm{eq}}`$ the gas slows down significantly).

A third possibility is that a small fraction of the gas has stopped rotating during the standard 500 ms equilibration time, and can therefore be at a higher temperature evaporation . This small fraction could decrease the visibility of the vortices present in the condensed part of the cloud. The existence of such a fraction could be tested by a quantitative analysis of the density profiles of Fig. 2ab. Such an analysis should take into account the TOF evolution of an interacting atomic cloud initially confined in a non harmonic potential.

To summarize we have presented in this Letter a direct evidence for a qualitative change in the nature of a degenerate Bose gas when it is rotated around and above the trapping frequency. We plan to complement this study by an analysis of other modes of the rotating gas, such as the transverse breathing mode Stock03 . Several other interesting phenomena have been predicted and remain to be investigated experimentally for this type of potential: existence of several phases involving either a vortex array with a hole, or a giant vortex Fetter01 ; Lundh02 ; Kasamatsu02 ; Fisher03 ; Kavoulakis03 . Another extension of the present work consists in transposing the experimental scheme to a 2D geometry, where the motion along $`z`$ would be frozen. The situation would then be the bosonic analog of the situation leading to the quantum Hall effect, provided the effective rotation frequency of the gas is precisely adjusted to the trap frequency.

###### Acknowledgements.

We thank P. Rosenbusch for participation in earlier stages of this experiment, and S. Stringari, L. Pitaevskii, G. Baym and the ENS group for useful discussions. This work is partially supported by CNRS, Collège de France, Région Ile de France, DAAD, DGA, DRED and EU (CQG network HPRN-CT-2000-00125).