text:
# Driven vortices in confined geometry: the Corbino disk

## 1 INTRODUCTION

In the mixed state of type-II superconductors the magnetic field is concentrated in an array of flexible flux bundles that, much like ordinary matter, can form crystalline, liquid and glassy phases. In clean systems the vortex solid melts into a flux liquid via a first order phase transition. If the barriers to vortex line crossing are high, a rapidly cooled vortex liquid can bypass the crystal phase and get trapped in a metastable polymer-like glass phase. The diversity of vortex structures is further increased by pinning from material disorder, which leads to a variety of novel glasses. Disorder-driven glass transitions are continuous, with diverging correlation lengths and universal critical behavior.

Of particular interest is the dynamics of the vortex array in the various phases and in the proximity of a phase transition. In the liquid phase the vortex array flows yielding a linear resistivity. In the presence of large scale spatial inhomogeneities, the liquid flow can be highly nonlocal due to interactions and entanglement. The correlation length controlling the nonlocality of the flow grows with the liquid shear viscosity, which becomes large as the liquid freezes. At a continuous liquid-glass transition this correlation length diverges with a universal critical exponent. In the solid phase the vortex array moves as a single elastic object under uniform drive, provided the shear stresses are not too large. In the presence of strong spatial inhomogeneities, plastic flow occurs for large drives (or even for vanishingly small drives in a glassy solid) and the response is always nonlinear. The dynamical correlation length can be identified with the separation between free dislocations and diverges at a continuous melting transition. Probing spatial velocity correlations can therefore give information on vortex dynamics within a given phase, as well as on the nature of the phase transitions connecting the various phases.

As for ordinary matter, the shear rigidity of the vortex array can be probed by forcing the vortices to flow in confined geometries.  This type of experiments was pioneered by Kes and collaborators to study the shear rigidity of the two-dimensional vortex liquid near freezing in thin films. More recently, patterned irradiation of cuprate superconductors with heavy ions has made it possible to create samples with controlled distributions of damage tracks. We recently showed that an analysis of such experiments that combines an inhomogeneous scaling theory with the hydrodynamics of viscous flux liquids can be used to infer the critical behavior near a continuous glass transition, as well as to distinguish between continuous transitions, such as that to a Bose glass, and nonequilibrium transition to a polymer-like glass driven by interaction and entanglement..

Large scale spatial inhomogeneities can also be introduced in the flow, even in the absence of pinning, by applying a driving force with controlled spatial gradients, as done recently by the Argonne group using the Corbino disk geometry. In this paper we illustrate the analysis of spatially inhomogeneous vortex motion in both the liquid and the solid using the Corbino disk as a prototype of a novel class of experiments exploiting the effect of geometry to study the dynamics of vortex matter.

## 2 LIQUID FLOW IN CHANNELS

In the Corbino disk, with magnetic field along the disk axis ($`z`$ direction), a uniform radial current density of magnitude $`J(r)=I/(2\pi tr)`$ is introduced in the sample by injecting current at the center and removing it at the outer circumference of the disk (inset of Fig. 1). The current drives the vortices to move in circles about the axis. In the flux liquid, the dynamics on scales larger than the intervortex spacing, $`a_0`$, is described by hydrodynamic equations for the flow velocity $`𝐯(𝐫)`$, which determines the local field from flux motion, $`𝐄=n_0\varphi _0\widehat{𝐳}\times 𝐯(𝐫)/c`$, with $`n_0=1/a_0^2`$. For simple geometries like the Corbino disk, where the current is spatially homogeneous in the $`z`$ direction, hydrodynamics reduces to a single equation,

$$\gamma 𝐯+\eta _{}^2𝐯=\frac{1}{c}n_0\varphi _0\widehat{𝐳}\times 𝐉(𝐫),$$
(1)
where $`\gamma (T,H)`$ is the friction, $`\eta (T,H)`$ is the viscosity controlling the viscous drag from interactions and entanglement, and the term on the right hand side is the Lorentz force density driving flux motion. It is instructive to rewrite Eq. (1) as an equation for the local field,

$$\xi ^2_{}^2𝐄+𝐄=\rho _f𝐉,$$
(2)
with $`\xi =\sqrt{\eta /\gamma }`$ the viscous correlation length and $`\rho _f=(n_0\varphi _0/c)^2/\gamma `$ the flux flow resistivity. If the viscous force is negligible, Eq. (2) is simply Ohm’s law and the radial field is $`E_0(r)=(\rho _fI/2\pi t)(1/r)`$.

To probe the viscous drag, it is necessary to force large scale spatial inhomogeneities in the flow. This may be achieved by suitable pinning boundaries. As an example, we imagine selectively irradiating a cylindrical central region and an outer annular region of the disk to obtain the structure sketched in the inset of Fig. 2. Here the vortices in the heavily irradiated central and outer regions (shaded) are in the Bose glass phase, while vortices in the unirradiated (white) annular region are in the flux liquid phase. A radial current drives tangential flow in the resistive flux liquid annulus, which is impeded by the “Bose-glass contacts” at the boundaries. The field profile obtained by solving Eq. (2) with no-slip boundary conditions  is spatially inhomogeneous on length $`\xi `$, as shown in Fig. 1. One can probe this profile and extract $`\xi `$ by placing a string of radial contacts at $`r_n`$, for $`n=1,2,3,\mathrm{}`$, and measuring the voltage $`V_{n+1,n}`$ across each successive pair (inset of Fig. 2). If the viscosity is small ($`\xi <<d`$), the voltage decreases logarithmically as one moves from the inner to the outer contacts, as in a freely flowing uncorrelated liquid, where $`V_{n,+1,n}^0=(\rho _fI/2\pi t)\mathrm{ln}(r_{n+1}/r_n)`$. When $`\xi `$ grows, the onset of rigidity in the liquid becomes apparent (Fig. 2). An elastic vortex solid would rotate as a rigid object under the radial drive, with $`v(r)r`$ and $`V_{n+1,n}^s=(\rho _fI/2\pi tR_2^2)(r_{n+1}^2r_n^2)`$, for $`R_2>>R_1`$. Indeed for $`\xi d`$, $`V_{n+1,n}`$ is no longer monotonic with $`n`$ and it exhibits a solid-like growth with $`n`$ within a boundary layer of width $`\xi `$.

The resistance per unit thickness of the disk, $`\rho _R(T,R_1,R_2)=\mathrm{\Delta }V(R_2,R_1)/(I/2\pi t)`$, is

$$\rho _R(T,R_1,R_2)=\rho _f(T)(\frac{R_1}{\xi },\frac{R_2}{\xi }),$$
(3)
with $`(x_1,x_2)`$ a function determined by the geometry that is obtained from the solution of the hydrodynamic equation. If $`\xi <<R_1,R_2`$, viscous effect are not important and $`\rho _R\rho _f\mathrm{ln}(R_2/R_1)1/\gamma `$. Conversely, when the viscous length exceeds the sample size ($`\xi >>R_1,(R_2R_1)`$), then $`\rho _R\rho _f(R_2^28\pi \xi ^2)1/\eta `$, for $`R_2>>R_1`$. When the effects of geometry are negligible, resistivity measurements directly probe the friction, $`\gamma `$, while measurements in channels narrow compared to $`\xi `$ probe the flux liquid viscosity, $`\eta `$. Experiments in artificial pinning structures can therefore be used to infer both $`\gamma `$ and $`\eta `$. The behavior of these two parameters near a phase transition can in fact be used to identify the transition, as summarized in table 1.

Of particular interest is the case of continuous glass transitions, which are characterized by universal critical behavior. Experiments of the type just described are especially powerful in this case as they can be used to map out the critical behavior. One can probe the Bose glass transition by lightly irradiating the liquid annular region, and then lowering the temperature at constant field from $`T_{BG}^{\mathrm{annulus}}<T<T_{BG}^{\mathrm{contacts}}`$ to $`TT_{BG}^{\mathrm{annulus}}`$. Most physical properties near the transition can be described via a scaling theory in terms of diverging correlation lengths perpendicular and parallel to the field direction, $`\xi _{}(T)|TT_{BG}|^\nu _{}`$ and $`\xi _{}(T)|TT_{BG}|^\nu _{}`$, with $`\nu _{}=2\nu _{}`$, and a diverging correlation time, $`\tau l_{}^z|TT_{BG}|^{z\nu _{}}`$. Scaling can then be used to relate physical quantities to these diverging length and time scales. In particular, the friction coefficient $`\gamma `$ that determines the bulk flux flow resistivity $`\rho _f(T)`$ is predicted to diverge as $`TT_{BG}^+`$ as $`\rho _f|TT_{BG}|^{\nu _{}(z2)}`$.  As shown recently by Marchetti and Nelson, when flux flow in confined geometries is analyzed by combining hydrodynamics with the the Bose glass scaling theory – generalized to the spatially inhomogeneous case – the Bose glass correlation length $`\xi _{}`$ is naturally identified with the viscous length $`\xi `$. It then follows that the liquid shear viscosity diverges at $`T_{BG}`$ as $`\eta |TT_{BG}|^{z\nu _{}}`$. Furthermore, the scaling of the finite-geometry resistivity displayed in Eq. 3 is a general property of the vortex liquid near a continuous glass transition. Hydrodynamics yields the precise form of the scaling function, which depends on the experimental geometry and can be found in  for the Corbino disk.

## 3 PLASTIC FLOW IN DRIVEN SOLIDS

As shown in recent experiments by the Argonne group, the Corbino disk geometry can also be used to study the onset of plastic flow in a driven solid. In this case we consider an unirradiated disk, where the vortex array has a clear melting transition. Below melting the vortex solid moves as a rigid body, with $`v(r)r`$, and the voltage grows as $`r^2`$. The $`1/r`$ dependence of the driving force yields, however, large elastic deformations of the medium, described by the solution of

$$c_{66}^2𝐮=\frac{1}{c}n_0\varphi _0\widehat{𝐳}\times 𝐉(𝐫),$$
(4)
with free boundary conditions ($`c_{66}`$ is the shear modulus of the lattice, assumed incompressible). Elastic deformations yield a finite shear stress, $`\sigma _{r\varphi }(r)=(n_0\varphi _0I/4\pi ct)\left[1+2\mathrm{ln}(R_2/R_1)R_1^2/r^2\right]`$, that can unbind dislocations from bound pairs. Assuming for simplicity that vortices are straight along the field direction, the energy of a pair of dislocations of opposite Burgers vectors, separated by a distance $`x`$, is $`U_0(x)=(c_{66}ta_0^2/\pi )\left[\mathrm{ln}(x/a_0)\mathrm{cos}^2\theta \right]2E_ct`$, with $`\theta `$ the angle between $`𝐱`$ and $`𝐛`$ and $`E_cc_{66}a_0^2`$ the core energy per unit length of an edge dislocation. An applied stress pulls the two dislocations in opposite directions. Ignoring climb, and assuming the spatial variations of the stress field are on scales large compared to $`xa_0`$, the interaction energy is now $`U(x)=U_0(x)a_0x\sigma _{r\varphi }(r)`$, with $`r`$ the radial location of the center of mass of the pair. The applied stress lowers the barrier that confines the bound pair. In the simplest model, unbinding occurs where the location $`x_B`$ of the barrier, defined by $`\left[U/x\right]_{x=x_B}=0`$, becomes comparable to $`a_0`$, or $`x_B=c_{66}a_0/(2\pi \sigma _{r\varphi }(r))a_0`$. By solving for $`r`$, one obtains the critical radius $`R_M(I)`$ where shear-induced dislocation unbinding occurs, yielding slippage of neighboring planes of the vortex lattice, $`R_MR_1\sqrt{2\mathrm{ln}(R_2/R_1)I/(I_0I)}`$, for $`R_2>>R_1`$, with $`I_0=2cc_{66}t/(n_0\varphi _0)`$ a maximum current above which the entire vortex solid shear melts. The melting radius $`R_M`$ increases with current, indicating that, since the stresses are largest near the axis of the disk, “shear-induced melting” occurs first in circular layers close to the axis. This behavior is qualitatively consistent with the observations by the Argonne group. The simple model described here suggests that at high fields the current scale $`I_0`$ is independent of field. A more detailed calculation incorporating field and temperature dependence will be described elsewhere.


format: text

text:
# CU-TP-982hep-th/0007038 Duality and Massless Monopoles11footnote 1Talk delivered at Quarks 2000, Pushkin Russia, May 2000

## 1 Introduction

A quarter century ago ’t Hooft and Polyakov showed  that magnetic monopoles could arise as topologically stable classical solutions in certain spontaneously broken gauge theories. The fact that these included all grand unified theories led to a considerable theoretical effort to study the detailed properties and the astrophysical and cosmological implications of these GUT monopoles, as well as to a number of experimental searches. However, the failure of these searches to detect any monopoles, together with astrophysical arguments that place stringent bounds on the monopole flux at the Earth’s surface, has quite understandably led to a decreased interest from the phenomenological point of view.

Nevertheless, these monopole solutions remain valuable as tools for probing the properties of quantum field theory. In particular, the properties of monopoles in the Bogomol’nyi-Prasad-Sommerfield (BPS) limit  has led to the conjecture  that the electric-magnetic duality of Maxwell’s equations might find a quantum field theoretic generalization that exchanges magnetically charged solitons with electrically charged elementary particles. The prime candidate  for this is $`N=4`$ supersymmetric Yang-Mills theory. In the case of SU(2) broken to U(1), the spectrum of massive particles is invariant under the simultaneous interchange of weak and strong coupling and of electric and magnetic charges. This result generalizes fairly easily to larger gauge groups, as long as the subgroup that remains unbroken is Abelian. However, if the unbroken subgroup is Abelian, the elementary particle sector contains massless particles with electric-type charge (the “gluons” and their superpartners). Duality would then require the existence of massless magnetically charged particles. However, one can easily show that these theories have no massless classical soliton solutions. In this talk, I will argue that, nevertheless, one can find evidence for the required massless monopoles by studying the dynamics of the massive monopoles.

## 2 BPS monopoles in SU(2) and larger groups

Throughout this talk I will work in the BPS limit, with monopoles obeying

$$B_i=D_i\mathrm{\Phi }.$$
(2.1)
While this limit can be obtained by taking a limit of coupling constants, it is most naturally understood in the context of Yang-Mills theory with extended supersymmetry, where Eq. (2.1) is equivalent at the classical level to requiring that the soliton preserve half of the supersymmetry.222Although I have $`N=4`$ supersymmetric Yang-Mills in mind, the only feature of this theory that I will use explicitly is the fact that the Higgs field is in the adjoint representation. Note also that there is no confinement in this theory, even if there is an unbroken non-Abelian subgroup. This ensures that the BPS mass formula

$$M=Q_Mv=𝑑S_i\mathrm{Tr}B_i\mathrm{\Phi }$$
(2.2)
is preserved by quantum corrections.

One important feature of the BPS limit is that the Higgs scalar becomes massless and so can mediate a long-range attractive force between two monopoles. It turns out that this can exactly balance the magnetic repulsion between static monopoles, thus allowing static multimonopole solutions to exist. In fact, for the case of SU(2) broken to U(1), all higher charged BPS solutions can be understood as multimonopole solutions of this sort. Not only does a solution with $`n`$ units of magnetic charge have $`n`$ times the mass of the unit monopole, but index theory methods show  that the number collective coordinates needed to quantize the solution is precisely $`4n`$, corresponding to a three-dimensional position and a U(1) phase for each of the $`n`$ constituent monopoles; excitation of these phase coordinates gives independent dyonic electric charges to the individual monopoles.

In order to have the possibility of symmetry breaking to a non-Abelian subgroup, we must start with a gauge group of rank $`r2`$. The generators of this gauge group can be chosen to be a set of $`r`$ commuting $`H_i`$ that form a basis for the Cartan subalgebra, together with a set of raising and lowering operators associated with the roots $`𝜶`$. The asymptotic adjoint Higgs field in some fixed direction can always be brought into the form

$$\mathrm{\Phi }_0=𝐡𝐇.$$
(2.3)
If the $`r`$-component vector $`𝐡`$ has nonzero inner products with all of the roots, then the gauge symmetry is broken maximally, to U(1)r; if not, then the roots orthogonal to $`𝐡`$ yield the root diagram of an unbroken non-Abelian group $`K`$.

The long range part of the magnetic field must lie in the unbroken part of the gauge group. Hence, in the direction used to define $`\mathrm{\Phi }_0`$, the leading part of asymptotic magnetic field can also be brought into the Cartan subalgebra and written in the form

$$B_i=𝐠𝐇\frac{\widehat{r}_i}{r^2}.$$
(2.4)
In the case of maximal symmetry breaking there are $`r`$ topological charges, one for each of the unbroken U(1) factors. The connection between these and the magnetic charge is most easily seen by choosing a set of $`r`$ simple roots $`𝜷_a`$. These form a basis for the root lattice of the Lie group with the property that any root can be written as a linear combination of simple roots with coefficients all of the same sign. There are many possible choices for the simple roots, but the vector $`𝐡`$ associated with the Higgs field can be used to pick out a unique set satisfying $`𝐡𝜷_a>0`$. The topological quantization condition  then takes the form

$$𝐠=\frac{4\pi }{e}\underset{a}{}n_a\frac{𝜷_a}{𝜷_a^2}$$
(2.5)
where the integers $`n_a`$ are the topological charges. The BPS mass formula, Eq. (2.2), can be written as

$$M=𝐠𝐡=\underset{a}{}n_a\left(\frac{4\pi }{e}𝐡𝜷_a\right)\underset{a}{}n_am_a,$$
(2.6)
while the number of collective coordinates is

$$𝒩=4\underset{a}{}n_a.$$
(2.7)
These results suggest that, just as in the SU(2) case, all higher charged solutions should be viewed as multimonopole configurations. Now, however, there are $`r`$ species of fundamental monopoles, one for each U(1) factor, with the $`a`$th species having mass $`m_a`$, topological charges $`n_b=\delta _{ab}`$, and four degrees of freedom (three for center-of-mass motion and one U(1) phase). These fundamental monopoles can be explicitly obtained by embedding the SU(2) unit monopole in the SU(2) subgroups defined by the various $`𝜷_a`$.

When the symmetry breaking is nonmaximal, to $`K\times \mathrm{U}(1)^{rk}`$, some of the simple roots, which I denote by $`𝜸_i`$, are orthogonal to $`𝐡`$ and form a complete set of simple roots for $`K`$. The remainder, denoted $`\stackrel{~}{𝜷}_a`$, can required, as before, to obey $`𝐡\stackrel{~}{𝜷}_a>0`$. Equation (2.5) is replaced by

$$𝐠=\frac{4\pi }{e}\underset{a}{}\stackrel{~}{n}_a\frac{\stackrel{~}{𝜷}_a}{\stackrel{~}{𝜷}_a^2}+\underset{i}{}q_i\frac{𝜸_j}{𝜸_j^2}.$$
(2.8)
The $`\stackrel{~}{n}_a`$ and $`q_i`$ are integers, with the former being the conserved topological charges.333In contrast with the maximally broken case, the simple roots are not uniquely determined by the requirement that their inner product with $`𝐡`$ be positive. The various allowed sets are related by gauge transformations of the unbroken group. The $`\stackrel{~}{n}_a`$ are invariant under these transformations, but the $`q_i`$ are not.

In general, the corresponding magnetic field has both Abelian and non-Abelian components. In order to avoid certain pathologies  associated with non-Abelian magnetic charges, I will assume for the remainder of this talk that $`𝐠𝜸_i=`$ for all $`i`$, so that the long-range magnetic field is purely Abelian. (There is little loss of generality in this assumption since, given a configuration with nonzero $`𝐠𝜸_i`$, the additional monopoles needed to cancel the non-Abelian part of the total magnetic charge can be placed at an arbitrarily large distance.) With this assumption, Eqs. (2.6) and (2.7) are replaced by

$$M=\underset{a}{}\stackrel{~}{n}_am_a$$
(2.9)
and

$$𝒩=4\underset{a}{}\stackrel{~}{n}_a+4\underset{i}{}q_i.$$
(2.10)
Analogy with the maximally broken case would then suggest that there is one fundamental monopole, with mass $`m_a`$ and four degrees of freedom, associated each of the $`\stackrel{~}{𝜷}_a`$, and one massless monopole, also with four degrees of freedom, associated with each of the $`𝜸_i`$. Indeed, the embedding construction for the massive fundamental monopoles goes through pretty much as before. The massless monopoles, on the other hand, cannot be constructed in this manner: applying the embedding construction to the SU(2) subgroup corresponding to one of the $`𝜸_i`$ simply yields a pure vacuum solution. In fact, it is easy to show that there are no localized classical solutions with zero energy.

## 3 Low energy monopole dynamics

Although we cannot obtain classical solutions corresponding to isolated massless monopoles, one can find evidence of these monopoles in the dynamics of multimonopole systems. The moduli space approximation  is a convenient tool for studying such systems at low energy. The essential idea is to approximate solutions with slowly moving monopoles444These includes dyons with small electric charges, since these correspond to slowly varying U(1) phases. as being motion on the moduli space of static BPS multimonopole solutions. More precisely, let {$`A_i^{\mathrm{BPS}}(𝐫,z)`$, $`\mathrm{\Phi }^{\mathrm{BPS}}(𝐫,z)`$} be a family of gauge-inequivalent BPS solutions parameterized by a set of collective coordinates $`z_a`$. In the moduli space approximation one adopts the Ansatz

$`A_0(𝐫,t)`$ $`=`$ $`0`$ (3.1)
$`A_i(𝐫,t)`$ $`=`$ $`U^1(𝐫,t)A_i^{\mathrm{BPS}}(𝐫,z(t))U(𝐫,t)iU^1(𝐫,t)_iU(𝐫,t)`$ (3.2)
$`\mathrm{\Phi }(𝐫,t)`$ $`=`$ $`U^1(𝐫,t)\mathrm{\Phi }^{\mathrm{BPS}}(𝐫,z(t))U(𝐫,t).`$ (3.3)
With this Ansatz, the time derivatives of the fields are of the form

$`\dot{A}_i`$ $`=`$ $`\dot{z}_j\left[{\displaystyle \frac{A_i}{z_j}}+D_iϵ_j\right]\dot{z}_j\delta _jA_i`$ (3.4)
$`\dot{\mathrm{\Phi }}`$ $`=`$ $`\dot{z}_j\left[{\displaystyle \frac{\mathrm{\Phi }}{z_j}}+[\mathrm{\Phi },ϵ_j]\right]\dot{z}_j\delta _j\mathrm{\Phi }.`$ (3.5)
The terms involving $`ϵ_j`$ arise from differentiating the gauge function $`U(𝐫,t)`$; they are fixed uniquely by Gauss’s law, which turns out to be equivalent to imposing a background gauge condition on $`\dot{A}_i`$ and $`\dot{\mathrm{\Phi }}`$.

Substituting this Ansatz into the Yang-Mills Lagrangian gives

$$L_{MS}=\frac{1}{2}d^3r\mathrm{Tr}\left[\dot{A}_{i}^{}{}_{}{}^{2}+\dot{\mathrm{\Phi }}^2+B_i^2+D_i\mathrm{\Phi }^2\right]=\frac{1}{2}g_{ij}(z)\dot{z}_i\dot{z}_j+M$$
(3.6)
where the static energy $`M`$ is constant on the moduli space and

$$g_{ij}(z)=d^3r\left[\delta _iA_k\delta _jA_k+\delta _i\mathrm{\Phi }\delta _j\mathrm{\Phi }\right].$$
(3.7)
If we interpret $`g_{ij}`$ as a metric on the moduli space, then the solutions to $`L_{MS}`$ are simply geodesic motion on the moduli space. In most cases, it turns out not to be practicable to use Eq. (3.7) directly to determine $`g_{ij}(z)`$. However, by using indirect methods the moduli space metrics for a number of interesting cases have been by more indirect methods. I will make use of some of these results in the next two sections.

## 4 An SO(5) example

The simplest example where one finds evidence of massless monopoles arises with the gauge group SO(5), whose root diagram is shown in Fig. 1. With the Higgs field vector $`𝐡`$ is in the direction shown on the left, the symmetry is broken to U(1)$`\times `$U(1). There are two species of massive fundamental monopoles, corresponding to the simple roots $`𝜷`$ and $`𝜸`$. If instead the Higgs vector is vertical, as shown on the right, the unbroken gauge group is SU(2)$`\times `$U(1), with the SU(2) having roots $`\pm 𝜸`$. The fundamental $`𝜷`$-monopole remains massive, but now $`𝜸`$ corresponds to an elusive massless monopole.

I will focus on solutions with

$$𝐠=\frac{4\pi }{e}\left(\frac{𝜷}{𝜷^2}+\frac{𝜸}{𝜸^2}\right).$$
(4.1)
In the maximally broken case, these are composed of two distinct massive fundamental monopoles, and have an eight-dimensional moduli space whose metric is

> Figure 1: The root diagram of $`SO(5)`$. With the Higgs vector $`𝐡`$ oriented as in (a) the gauge symmetry is broken to U(1)$`\times `$U(1), while with the orientation in (b) the breaking is to SU(2)$`\times `$U(1).

$`ds^2`$ $`=`$ $`Md𝐗_{\mathrm{cm}}^2+{\displaystyle \frac{16\pi ^2}{M}}d\chi _{\mathrm{tot}}^2+\left(\mu +{\displaystyle \frac{k}{r}}\right)\left[dr^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)\right]`$ (4.3)
$`+k^2\left(\mu +{\displaystyle \frac{k}{r}}\right)^1(d\psi +d\mathrm{cos}\theta d\varphi )^2.`$
Here $`𝐗_{\mathrm{cm}}`$ is the center-of-mass position of the two monopoles; $`r`$, $`\theta `$, and $`\varphi `$ are their relative coordinates; and $`\chi _{\mathrm{tot}}`$ and $`\psi `$ are overall and relative U(1) phases. $`M`$ and $`\mu `$ are the total and reduced masses of the monopole pair, and $`k`$ is a numerical constant related to the normalization of roots.

With SU(2)$`\times `$U(1) breaking, Eq. (4.1) corresponds to a combination of one massive and one massless monopole. It turns out to be relatively straightforward to solve the field equations explicitly . The solutions are spherically symmetric and depend on eight parameters: three position coordinates X, four SU(2)$`\times `$U(1) phase angles $`\alpha `$, $`\beta `$, $`\gamma `$, and $`\chi `$. and one last parameter, $`a`$, that can take on any real positive value and whose interpretation will become clear shortly. The fields can be decomposed according their transformations under the unbroken SU(2). The singlet terms are just what would be obtained by embedding the unit SU(2) monopole using the composite root $`2𝜷+𝜸`$; they have no dependence on $`a`$. By themselves, they would describe an object with a massive core surrounded by a simple Coulomb magnetic field. The doublet terms fall exponentially fast outside the monopole core and are relatively uninteresting. Finally, the triplet parts of the gauge and Higgs fields are given by invariant tensors multiplied by the function

$$G(r,a)=\left[\frac{v}{\mathrm{sinh}(evr)}\frac{1}{er}\right]\left[1+\frac{r}{a}\mathrm{coth}(vr/2)\right]^1.$$
(4.4)
For $`r\begin{array}{c}<\hfill \\ \hfill \end{array}a`$, the second factor on the right hand side is approximately unity, so $`G`$, and hence the corresponding $`A_i`$, falls as $`1/r`$. This produces a Coulomb magnetic field in the unbroken SU(2) subgroup. However, for $`r\begin{array}{c}>\hfill \\ \hfill \end{array}a`$, the last factor gives an additional $`1/r`$, so the triplet part of $`A_i`$ falls as $`1/r^2`$, and the non-Abelian component of the magnetic field fall more rapidly.

Hence, we may view the solution as being composed of a massive core of radius $`(ev)^1`$ surrounded by a “cloud” of non-Abelian fields of radius $`a`$. Inside this cloud one finds the magnetic field appropriate to a monopole with both Abelian and non-Abelian magnetic charges; i.e., the charge appropriate to an isolated massive $`𝜷`$-monopole. Outside the cloud, only the Abelian charge is evident. Thus, the massless $`𝜸`$-monopole can be viewed as a shell of radius $`a`$ surrounding the massive monopole. Curiously, the energy of the solution is independent of $`a`$, despite the fact that this parameter is not associated with any symmetry of the system.

We can proceed further and examine the moduli space metric. Using the explicit form of the solutions and Eq. (2.5), one obtains

$$ds^2=Md𝐗^2+\frac{16\pi ^2}{M}d\chi ^2+k\left[\frac{da^2}{a}+a\left(d\alpha ^2+\mathrm{sin}^2\alpha d\beta ^2+(d\gamma +\mathrm{cos}\alpha d\beta )^2\right)\right]$$
(4.5)
where $`k`$ is the same constant as in Eq. (4.3).

Now let us return for a minute to the maximally broken case and imagine approaching the SU(2)$`\times `$U(1) case by “rotating” the Higgs vector until it is vertical. This corresponds to taking the $`𝜸`$-monopole mass, and thus the reduced mass $`\mu `$, to zero. Taking this limit in the metric of Eq. (4.3) gives precisely Eq. (4.5), except for a curious change in notation: The monopole separation $`r`$ becomes the cloud parameter $`a`$, while the spatial angles $`\theta `$ and $`\varphi `$ combine with the relative U(1) phase $`\psi `$ to give the SU(2) Euler angles.

## 5 A more complex system

An example with somewhat more structure is obtained by considering the case of $`\mathrm{SU}(N+2)`$ broken to $`\mathrm{U}(1)\times \mathrm{SU}(N)\times \mathrm{U}(1)`$, with the unbroken $`\mathrm{SU}(N)`$ lying in the middle $`N\times N`$ block in a basis where the eigenvalues of the Higgs field decrease monotonically along the diagonal. A purely Abelian asymptotic magnetic field can be obtained by setting all of the $`\stackrel{~}{n}_j`$ and $`q_i`$ in Eq. (2.8) equal to unity; i.e., by combining one each of the two species of massive and $`N1`$ species of massless monopoles. The moduli space of such solutions is $`4(N+1)`$-dimensional, with the collective coordinates including a position and U(1) phase for each of the massive monopoles, a number of SU($`N`$) orientation angles, and a single cloud parameter, $`b`$, that can take on any value greater than or equal to $`r`$, the separation between the two massive monopoles.

The solutions can be obtained explicitly  by using the Nahm construction . Although their detailed structure is somewhat complex, their behavior well outside the monopole cores is fairly simple. The asymptotic Higgs field can be written as

$$\mathrm{\Phi }_{\mathrm{}}(𝐫)=U^1(𝐫)\mathrm{diag}(v_3,v_2,\mathrm{},v_2,v_1)U(𝐫)$$
(5.1)
where $`U(𝐫)`$ is a gauge transformation whose form will not concern us and the Higgs eigenvalues satisfy $`v_3>v_2>v_1`$. The unbroken SU($`N`$) lies in the middle $`N\times N`$ block corresponding to the repeated eigenvalue $`v_2`$. The form of the magnetic field depends on how the distances $`y_L`$ and $`y_R`$ from the two massive monopoles compare to $`b`$. For $`y_L,y_Rb`$,

$$𝐁(𝐫)=U^1(𝐫)\mathrm{diag}(\frac{\widehat{𝐲}_R}{2y_R^2},\frac{\widehat{𝐲}_R}{2y_R^2},\frac{\widehat{𝐲}_L}{2y_L^2},0,\mathrm{},0,\frac{\widehat{𝐲}_L}{2y_L^2})U(𝐫)+\mathrm{}$$
(5.2)
where the dots represent terms that fall off more rapidly with distance, while for $`y(y_L+y_R)/2b`$

$$𝐁(𝐫)=U^1(𝐫)\mathrm{diag}(\frac{\widehat{𝐲}}{2y^2},0,0,0,\mathrm{},0,\frac{\widehat{𝐲}}{2y^2})U(𝐫)+\mathrm{},$$
(5.3)
Thus, at distances smaller than $`b`$ one sees fields corresponding to both Abelian and non-Abelian magnetic charges, but at larger distances only the Abelian component survives, just as in the SO(5) example.

It is instructive to study the scattering of the two massive monopoles in these solutions. Because these monopoles are associated with orthogonal roots \[or, equivalently, because they correspond to embeddings into commuting $`2\times 2`$ blocks of the SU($`N`$) matrices\], there is no direct interaction between them. Instead, they must interact through the massless monopole cloud. In the low energy limit, we can use the moduli space approximation to study this process . The moduli space metric for the maximally broken case with $`N+1`$ distinct massive fundamental monopoles is known . The SO(5) example of the previous section suggests that this metric is still valid when $`N1`$ of these monopoles become massless, although, just as in that case, the physical interpretations of the collective coordinates change in the massless limit, with the positions of the massless monopoles being transformed into gauge orientations and a single gauge invariant parameter, $`b`$. The various gauge orientation and angular variables are most easily handled by expressing them in terms of the conserved charges and angular momentum. This still leads to a fairly complicated dynamics, but matters simplify considerably if the electric-type charges in the unbroken $`\mathrm{U}(1)\times \mathrm{SU}(N)\times \mathrm{U}(1)`$ all vanish. After separating out the total center-of-mass motion, one is then left with an effective Lagrangian

$$L=\left[\frac{\mu }{2}+\frac{kb}{4(b^2r^2)}\right]\dot{r}^2+\frac{kb}{4(b^2r^2)}\dot{b}^2+\frac{br}{b^2r^2}\dot{b}\dot{r}\frac{bJ^2}{2r^2(k+\mu b)}$$
(5.4)
governing the time development of $`b`$ and $`r`$. Here $`J`$ is the angular momentum, $`\mu `$ is the reduced mass of the massive monopoles and, as in the previous section, $`k`$ is a numerical constant.

Analysis of the Euler-Lagrange equations that follow from Eq. (5.4) shows that $`r`$ and $`b`$ decouple at large times. Thus, at large $`|t|`$ the energy of the system is approximately the sum of a massive monopole kinetic energy $`E_r=\mu \dot{r}^2/2`$ and a cloud energy $`E_b=\dot{b}^2/4b`$, with $`r`$ and $`b`$ varying asymptotically as

$`r`$ $``$ $`v|t|+\mathrm{}`$
$`b`$ $``$ $`ct^2+\mathrm{}`$ (5.5)
with $`v`$ and $`c`$ constants. In a typical scattering process the cloud initially decreases in size while simultaneously the two massive monopoles approach one another. Eventually $`b`$ and $`r`$ reach their minimum values (generally not at the same time) and then begin to to increase again according to Eq. (5.5). In the course of this process, energy is exchanged between the cloud and the massive monopoles, so that the final splitting between $`E_r`$ and $`E_b`$ at $`t\mathrm{}`$ is different from the initial one at $`t\mathrm{}`$.

One surprising feature of Eq. (5.5) is the quadratic growth of $`b`$ with time. It seems likely that this is an artifact of the moduli space approximation, and that radiation of massless gluons will have the effect of reducing this to a linear growth.

## 6 Summary and conclusion

Electric-magnetic duality in gauge theories with the symmetry spontaneously broken to a non-Abelian subgroup requires the existence of massless magnetically charged objects that would be the duals to the electrically charged massless gauge bosons. These massless monopoles cannot be exhibited as isolated classical soliton solutions. However, as I have shown, their existence is manifested through clouds of non-Abelian field surrounding one or more massive monopoles. These clouds are described by a small number of collective coordinates, which correspond to the massless monopole degrees of freedom; for static solutions, the energy is independent of the values of these parameters.

Duality suggests that these cloud parameters should have counterparts in the perturbative sector. It would be instructive to understand more precisely what these are, and to look in perturbative sector scattering for analogies with the low-energy scattering of massive monopoles and massless monopole clouds. Investigations in this direction would certainly lead to deeper insights into the properties of non-Abelian gauge theories.

I would like to express my appreciation to the organizers of this conference. This work was supported in part by the U.S. Department of Energy.


text:
# Discovery of hard X-ray emission from Type II bursts of the Rapid Burster

## 1. Observations

Four Target Of Opportunity (TOO) observations were performed with BeppoSAX (Boella et al. 1997a) on the Rapid Burster (=MXB 1730–335; hereafter RB) during the activity state which started on January 28, 1998 (Fox et al. 1998). These TOOs spanned over one month (from February 18 to March 18) and caught the object in four different snapshots, from the post–maximum decay to the quiescent state. Figure 1, left panel, shows the ASM light curve of the Rossi-XTE satellite with superimposed the times of the four BeppoSAX observations. Here we report on RB data from three of the four instruments mounted on BeppoSAX: LECS (0.1-10 keV; Parmar et al. 1997), MECS (1.5-10 keV; Boella et al. 1997b, and PDS (15-300 keV; Frontera et al. 1997). For the PDS the default rocking collimator law was modified by offsetting the RB by 40$`^{^{}}`$ from the center of the field of view in order to reduce as much as possible the contamination from a nearby variable X–ray source, GX 354-0 (=4U 1728-34), located at about 30$`^{^{}}`$ from the RB. Unfortunately, due to failure of the rocking law setup program, the collimator did not move as requested during TOO2 and TOO3; so, we have only LECS and MECS data for these two observations.

In this paper we report on preliminary results of these observations. Definitive results along with their implications will be the subject of another paper (Masetti et al. 2000). In the following, for the luminosity estimates we will assume that the RB lies at a distance $`d`$ = 8 kpc (Ortolani et al. 1996).

## 2. Spectral analysis and temporal evolution of the RB

During TOO1, the RB was in a strong state of bursting activity. The 2-10 keV light curve obtained with MECS (see Fig. 1, right panel, part a) showed 113 Type II X–ray bursts during 9457 seconds of good observational data. Evidence of Type II bursts was also observed in the 0.1-2 keV data obtained with LECS. We divided the MECS TOO1 data into two subsets: persistent emission (PE; below 5 counts s-1) and bursting emission (BE; above 5 counts s-1).

The MECS PE and BE spectra could be well fit with a photoelectrically absorbed two-component blackbody (2BB); these BB components may originate from the neutron star (NS) surface, a boundary layer between the NS and the inner edge of the accretion disk, or the inner region of the disk itself. The same model was used for the RB by Guerriero et al. (1998) who found values consistent with ours. In Table 1 we report the best-fit parameters along with their 90% confidence errors. The temperature values of the two BB components were slightly higher during the PE than during the BE, while their luminosities were much higher (by a factor 20 for the cooler BB and 60 for the hotter BB) during the BE than during the PE. We also remark that during the BE the hotter BB component was brighter (by more than a factor 3) than the cooler BB, while during the PE they had similar luminosities. This implies that the BE influences more the higher temperature component than the other one. If the BE is due to spasmodic accretion onto the compact object, the higher temperature component should be the one coming from the NS surface.

The source was also visible in the hard X–ray (15-100 keV) energy range. However the statistics of the PDS light curve was much lower and did not allow distinguishing the Type II bursts. In order to construct the BE spectrum we used the time intervals in which the bursts were observed with MECS. Also, we could not derive the correct 15-100 keV flux and spectrum of both BE and PE given the residual source contamination by GX 354-0. Thus, in order to overcome this problem, we used as background level for the 1-100 keV BE spectrum the total count rate level measured during the PE time intervals. The combined LECS+MECS+PDS PE-subtracted bursting spectrum, shown in Fig. 2, was no longer fit with a 2BB model. By adding a power law component we obtained an acceptable fit (see Table 1). The further addition of a Fe K emission line at 6.5 keV slightly improved the fit, with parameter values found for this line in general agreement with the findings by Stella et al. (1988) for Type II bursts.

During TOO2 the RB drastically reduced its bursting activity, and the bursts were concentrated at the beginning of this TOO (Fig. 1, right panel, part b). Also, the emission intensity level decreased. The best–fit model was a photoelectrically absorbed 2BB model (Table 1). No evidence of a Fe emission line was present.

During TOO3 the object further reduced its bursting activity, and no Type II bursts were seen throughout the observation. The best–fit model spectrum was still an absorbed 2BB (Table 1). As in the case of TOO2, no iron emission line at 6.5 keV was found.

The RB was instead no longer visible in MECS/LECS images during TOO4. Stray light from GX354-0 prevented us to get a deep observation of the source. The 3$`\sigma `$ upper limit to the RB X–ray emission was 1.5$`\times `$10-12 erg cm-2 s-1 in the 2-10 keV energy band.

## REFERENCES

Boella G., Butler R.C., Perola G.C. et al., 1997a, A&AS 122, 299
Boella G., Chiappetti L., Conti G. et al., 1997b, A&AS 122, 327
Fox D., Guerriero R., Lewin W.H.G., 1998, ATEL n. 9
Frontera F. Costa E., Dal Fiume D. et al., 1997, A&AS 122, 357
Guerriero R., Fox D., Kommers J. et al., 1999, MNRAS 307, 179
Masetti N., Frontera F., et al., 2000, A&A, submitted
Ortolani S., Bica E., Barbuy B., 1996, A&A 306, 134
Parmar A., Martin D.D.E., Bavdaz M. et al., 1997, A&AS 122, 309
Stella L., Haberl F., Lewin W.H.G. et al., 1988, ApJ 324, 379


text:
# 1 The data on 𝑇₀ obtained by Sastry [] compared with our analytical results.

Comment on ”Liquid Limits: The Glass Transition and Liquid-Gas Spinodal Boundaries of Metastable Liquids”

In a recent paper , Sastry studies by numerical simulations the phase diagram of a simple fragile glass-forming liquid, the binary mixture with Lennard-Jones (LJ) interactions introduced by Kob and Andersen, presenting very interesting and clear results. Among other things, he obtains the position of the liquid-glass transition temperature as function of density $`\rho `$. He studies both the temperature $`T_0(\rho )`$ where the diffusion coefficient (of $`A`$ particles) drops to zero according to the Vogel-Fulcher-Tammann-Hesse law, and the temperature where the extrapolation of data on the configurational entropy $`S_c(\rho ,T)`$ vanishes, $`S_c`$ being evaluated as the difference between the liquid entropy and the harmonic amorphous solid entropy. These two ways give values remarkably close.

We apply to the LJ binary mixture, at various density values, the analytic approach to structural glass thermodynamics recently introduced , which allows to obtain both the liquid-glass transition temperature $`T_K`$ and properties in the glassy phase below $`T_K`$ doing a first principles computation. In one version of this approach, which implements analitically the proposal of inherent structure dominance at low temperature, one needs to compute the liquid free-energy density $`f`$ and two-point correlation functions $`g(r)`$ at $`T\stackrel{>}{}T_K`$.

The theory and different ways for implementing it are discussed in detail elsewhere . We only recall that, within this framework, at $`T_K`$ the system undergoes a second order transition corresponding to a one step replica simmetry breaking. As one could expect intuitively, it is found that in the region $`T\stackrel{>}{}T_K`$, where the phase space can be partitioned in ‘valleys’, the entropy of the liquid corresponds to the entropy of a typical valley $`S_{sol}`$ plus the configurational entropy $`S_c`$, and the Kauzmann temperature is defined by the condition $`S_c(T_K)=0`$.

We use the harmonic approximation scheme, which gives a simple formula for the single valley entropy in the liquid phase:

$$S_{sol}=\frac{3N}{2}\mathrm{ln}\left(2\pi eT\right)\frac{1}{2}\text{T}r\mathrm{ln},$$

where $``$ is the $`3N\times 3N`$ Hessian (i.e. the matrix of second derivative of the potential with respect to the $`Nd`$ particle coordinates) and $``$ means the Boltzmann average. We are therefore approximating $`S_{sol}`$ by only keeping the harmonic vibrational contribution. Further approximations are needed for computing the spectrum. Here we consider fluctuations of the diagonal part of $``$ up to second order , while the whole nondiagonal contribution is summed up by means of a kind of chain approximation.

In the computation of the liquid free energy and of the correlation functions in the high density region we use the Hypernetted Chain Approximation (HNC). Unfortunately this approximation overestimates the gas-liquid transition temperature and a more refined one is necessary at low densities. Liquid quantities at $`\rho =1.2`$ have been computed  within a scheme similar to that introduced by Zerah and Hansen (ZH), which also has the advantage of reducing the error in the approximation (of order $`10\%`$ in the HNC case).

We plot in \[Fig. 1\] the theoretical predictions on $`T_K`$ at different density values, compared with the data on $`T_0`$ by Sastry. Despite of the error in liquid quantities obtained by HNC, there is a quantitative agreement which makes evident that the analytic approach works well and it can be used as a powerful tool for a better understanding of glassy systems.

To conclude, we stress that, at least in this case, the Kauzmann temperature turns out to be directly comparable with the temperature where the diffusion drops to zero, which strongly supports the theoretical expectation of an effective ergodicity breaking.

B. Coluzzi(a), G. Parisi(a) and P. Verrocchio(b)

a) Dipartimento di Fisica, INFN and INFM, Università di Roma La Sapienza, P.le A. Moro 2, 00185 Roma, Italy.

b) Dipartimento di Fisica, Università di Trento, I-3805 Povo, Trento, Italy.

PACS numbers: 64.70.Pf,63.50.+x


text:
# Spectroscopic electric-field measurement in laser-excited plasmas

## Abstract

The electric field in mm-sized clouds of cold Rb+-ions is measured using the Stark effect of Rydberg atoms embedded in the ion clouds. The ion clouds are produced by UV photoionization of laser-cooled rubidium atoms in a magneto-optic trap. In very dense ion clouds, the Rydberg atom excitation spectra reflect the presence of a neutral plasma core and a positively charged ion shell.

Highly excited Rydberg atoms  are very sensitive to static and AC electric and magnetic fields, and can be used for the spectroscopic detection and compensation of such fields. Rydberg atoms have been used to detect microwave and FIR radiation , to measure and to compensate weak static electric fields , and to study QED level shifts and cavity-QED effects involving only few microwave photons . The spectroscopy of Rydberg atoms has been employed to diagnose electric fields in DC and high-frequency discharge plasmas . In this paper, we use Rydberg atom spectroscopy to determine the electric field in millimeter-sized, strongly coupled plasmas that since recently can be generated by the laser excitation of laser-cooled clouds of atoms . These plasmas hold the promise of bearing novel low-temperature recombination  and oscillation  behavior. Further, there is evidence for a Mott transition, by which strongly-coupled plasmas and/or dense gases of Rydberg atoms may form a new state of highly excited, metastable matter . To investigate these and other issues, one could utilize the non-invasive field-measurement method employed in the present paper.

Our experiment has recently been described . In a two-step optical excitation, laser-cooled rubidium atoms confined in a volume $`<1`$mm3 are first excited from the ground state $`5S_{1/2}`$ to the $`5P_{3/2}`$ state, which has an ionization wavelength of $`\lambda _{\mathrm{ion}}=479.1`$nm. A UV laser pulse ($`\lambda _{\mathrm{UV}}=355`$nm, duration$`<`$10ns, pulse energy up to 10mJ) then partially ionizes the $`5P_{3/2}`$-atoms, producing photoelectrons with 0.9eV initial kinetic energy. The photoelectrons escape on a time scale of a few ns, leaving behind a cloud of slow Rb+-ions. About 23ns after the UV pulse and while the lower-step laser is still on, the remaining $`5P_{3/2}`$-atoms that are embedded in the ion plasma are excited by a blue dye laser pulse ($`\lambda 480`$nm, duration$`<`$10ns, bandwidth $`10`$GHz). The ion cloud does not significantly expand before the blue pulse arrives. The blue pulse excites bound Rydberg states, the Stark effect of which we use to measure the plasma electric fields.

We detect the Rydberg atoms using the mechanism explained in the following. The ion plasma represents an electron trap with a depth $`U_0`$ given by the ion number, the cloud size and the profile of the charge distribution . In most of our experiments $`U_0<0.9`$eV, in which case no or only few UV-generated photoelectrons are retained in the electron trap (the ones that are moderated to low enough energy on their way out). However, the plasma electron trap becomes populated with electrons of about $`10`$meV energy due to direct photoionization of $`5P_{3/2}`$-atoms by the ASE with $`\lambda <\lambda _{\mathrm{ion}}`$ contained in the blue laser pulse and due to thermal ionization of a fraction of the bound Rydberg atoms excited by the blue pulse . Within a few $`\mu `$s, $`l`$-changing collisions between the trapped electrons and the Rydberg atoms efficiently promote the latter into long-lived high-$`l`$ states. During the subsequent Coulomb expansion of the plasma, which typically takes a few tens of $`\mu `$s, the depth $`U_0`$ of the plasma electron trap approaches zero. During that time, the electrons trapped in the plasma are gradually released, producing a plasma decay signal in our microchannel-plate (MCP) electron detector, which is located about $`10`$cm from the atomic cloud. A large fraction of the surviving high-$`l`$ Rydberg atoms thermally ionize on a time scale of order $`10`$ms , which exceeds the plasma decay time by about a factor of thousand. The delayed thermal electrons originating in the high-$`l`$ Rydberg atoms are counted in a gate that starts well after the plasma decay, yielding a noise-free signal of the Rydberg atom excitation. Rydberg excitation spectra are obtained by recording the delayed electron counts as a function of the wavelength of the blue laser.

Fig. 1 shows a typical experimental Rydberg excitation spectrum of atoms in an electric-field-bearing plasma. The field manifests itself through a) the appearance of $`p`$\- and hydrogenic lines, labeled $`p`$ and $`h`$ in Fig. 1, b) termination of the discrete Rydberg series at a certain field-dependent value of the principal quantum number $`n`$, and c) continuum lowering. While the effects a) and b) are solely determined by the state of the ion plasma at the time instant of the Rydberg excitation, the feature c) also depends on the history of the Rydberg atoms after their excitation. Since we intend to analyze the plasma at the time instant of the Rydberg atom excitation, we disregard the feature c) as a field measurement tool. We also don’t use b), because we have found that in many of our spectra the laser linewidth terminates the Rydberg series, but not the electric field.

The features of the above type a) develop as follows. As $`n`$ is increased, the electric field first manifests itself in the appearance of $`p`$-lines and of triangular-shaped $`h`$-features (see Fig. 1), caused by electric-field-induced state mixing. Due to the linear Stark effect, the hydrogenic states fan out over an energy range $`3n^2E`$ (atomic units) ( and Fig. 2). As $`n`$ increases further, the $`h`$-features progressively expand and fill in the spectral gaps of originally zero oscillator strength between the discrete non-hydrogenic $`s`$-, $`p`$\- and $`d`$-lines. There are three types of such spectral gaps, which we label $`hs`$-, $`mid`$\- and $`pd`$-gaps. The gaps become filled in with significant oscillator strength in exactly that order, at quite well defined critical wavelengths. Those wavelengths can be readily converted into critical effective quantum numbers $`n_{\mathrm{hs}}`$, $`n_{\mathrm{mid}}`$, and $`n_{\mathrm{pd}}`$. The $`n_\mathrm{i}`$ are robust indicators for the electric field, because they are solely determined by the general spreading behavior of the quasi-continuous $`h`$-features. No absolute line strength or line-strength ratio needs to be determined.

Assuming a constant electric field $`E_0`$, the values of $`n_\mathrm{i}`$ at which the spectral gaps become filled in can be related to $`E_0`$ using the Stark map of Rb. As shown in Fig. 2, the $`hs`$-gap below the hydrogenic levels with principal quantum number $`n`$ disappears when the energy shift of the most red-shifted hydrogenic state, which is approximately $`3n^2E_0/2`$, equals the splitting $`0.13n^3`$ between the unshifted hydrogenic energy and the next-lower $`s`$-state. Further, Fig. 2 shows that the $`mid`$-gaps, i.e. the gaps between the states $`(n+3)s`$ and $`(n+1)d`$ and between the field-free hydrogenic states $`n1`$ and the states $`(n+2)p`$, become entirely covered with hydrogenic states at practically identical electric fields. Finally, the $`pd`$-gaps disappear when neighboring hydrogenic manifolds meet, i.e. when $`3n^2E_0=n^3`$. The relations between $`E_0`$ and the critical quantum numbers $`n_{\mathrm{i0}}`$ are obtained as

$`E_0=`$ $`0.086n_{\mathrm{hs0}}^5`$ (1)
$`E_0=`$ $`0.23n_{\mathrm{mid0}}^5`$ (2)
$`E_0=`$ $`0.33n_{\mathrm{pd0}}^5(\mathrm{atomic}\mathrm{units}).`$ (3)
The zeros in the subscripts indicate that these relations hold for a constant electric field. Atomic species with quantum defects different from Rb would follow similar relations with different numerical factors. Further, the Eqs. 3 share, as expected, the $`n^5`$ dependence of the familiar Inglis-Teller relation .

To model our experimental spectra, we have calculated the photoabsorption spectra $`S(E_0,\lambda )`$ of randomly oriented atoms in constant electric fields $`E_0`$. Since from the $`5P`$ intermediate state we excite three manifolds of the magnetic quantum number $`|m|`$, one spectrum $`S(E_0,\lambda )`$ is a properly weighted sum of an $`m=0`$, an $`m=1`$ and an $`m=2`$ spectrum. The spectra are obtained as outlined in  using a spherical basis set of all states up to $`n=90`$. To account for the laser linewidth, the spectra are convoluted with Gaussians of a few GHz FWHM. The calculated spectra, one of which is shown in Fig. 3 a), display the same qualitative features as our experimental spectra. This applies, in particular, to the triangular shape of the $`h`$-features, which is produced by the unresolved excitation of the states in the hydrogenic manifolds.

An investigation of the calculated photoexcitation spectra reveals that the closing of the $`pd`$-gap appears particularly sudden due to the following peculiarity: for parameters between the $`mid`$-gap and the $`pd`$-gap closing, indicated by a double-arrow in Fig. 2, the blue-shifting hydrogenic levels within the $`pd`$-gap, marked with open circles in Fig. 2, carry much less oscillator strength than all the other nearby hydrogenic levels. As a result, the whole $`pd`$-gap appears virtually oscillator-strength-free up until the $`pd`$-gap closing condition $`3n^2E_0n^3`$ is satisfied. Once the $`pd`$-gap closing condition is satisfied, the oscillator strength of the hydrogenic states inside the $`pd`$-gap increases rather abruptly. This behavior enhances the clarity of the $`pd`$-gap closing.

The electric field experienced by the $`5P_{3/2}`$-atoms is not homogeneous but follows a probability distribution $`P(E)`$. Assuming identical Gaussian density profiles for the ions and the $`5P_{3/2}`$-atoms embedded in between them, the distribution $`P(E)`$ can be calculated by a random placement of ion charges and $`5P_{3/2}`$-atoms using a radial distribution function $`P(r)=4\pi r^2(2\pi \sigma ^2)^{3/2}\mathrm{exp}(r^2/(2\sigma ^2))`$. The resultant field distributions $`P(E)`$ have well defined most probable electric fields $`E_{\mathrm{max}}`$, which practically coincide with the singularities of the distributions produced by the corresponding smoothed charge distributions (compare the solid and dotted curves in Fig. 4). The effect of the microfields of the discrete ion charges largely is to smear out the singularities produced by the smoothed charge distributions, adding fields significantly larger than $`E_{\mathrm{max}}`$ to the field distribution. Deviations of the density profiles of the ions and the $`5P_{3/2}`$-atoms from Gaussians will modify $`P(E)`$, a dominant peak of $`P(E)`$. However, a most probable field $`E_{\mathrm{max}}`$ and a microfield-induced tail of fields $`E>E_{\mathrm{max}}`$ will remain.

We have used simulated electric-field distributions $`P(E)`$ to form averages of the spectra $`S(E_0,\lambda )`$ weighted by $`P(E)`$. As shown in panel b) of Fig. 3, one can still easily obtain critical quantum numbers $`n_{\mathrm{hs}}`$, $`n_{\mathrm{mid}}`$, and $`n_{\mathrm{pd}}`$ at which the respective regions of low oscillator strength disappear. The excitation spectra a) and b) in Fig. 3 are, in fact, quite similar, because the $`P(E)`$ are clearly dominated by fields near the most probable one. The microfield-induced spread of $`P(E)`$ from $`E_{\mathrm{max}}`$ towards higher fields causes the spectral gaps with low oscillator strength to become filled in at values of $`n`$ that are slightly lower than in the case of a homogeneous field $`E_0=E_{\mathrm{max}}`$.

Evaluating a number of spectra as shown in Fig. 3 b), we have determined the relations between the most probable fields $`E_{\mathrm{max}}`$ and the critical effective quantum numbers $`n_{\mathrm{hs}}`$, $`n_{\mathrm{mid}}`$, and $`n_{\mathrm{pd}}`$ for Rydberg atoms excited in Gaussian ion clouds with a FWHM of about 1mm:

$`E_{\mathrm{max}}=`$ $`0.050n_{\mathrm{hs}}^5`$ (4)
$`E_{\mathrm{max}}=`$ $`0.145n_{\mathrm{mid}}^5`$ (5)
$`E_{\mathrm{max}}=`$ $`0.26n_{\mathrm{pd}}^5(\mathrm{atomic}\mathrm{units}),`$ (6)
We have used Eqs. 6 to experimentally determine the values of $`E_{\mathrm{max}}`$ in ion clouds with diameters of $`1`$mm and various ion numbers. The spectra obtained in a few such experiments are displayed in Fig. 5. The critical wavelengths, indicated by arrows, have been determined as shown in Figs. 1 and 3, and have been converted into critical quantum numbers $`n_{\mathrm{hs}}`$, $`n_{\mathrm{mid}}`$ and $`n_{\mathrm{pd}}`$. Using Eq. 6, up to three spectroscopic measurements of $`E_{\mathrm{max}}`$ are obtained for each spectrum. The relative uncertainty of each single electric-field value is $`30\%`$. Within that uncertainty, the different electric-field values obtained for the individual spectra agree well. The electric field values quoted in Fig. 5 are the averages over all field values obtained for the respective spectra.

The most probable field $`E_{\mathrm{max}}`$ in Gaussian ion clouds can be calculated as $`E_{\mathrm{max}}=(3.1\times 10^{10}N_{\mathrm{ion}}\sigma ^2)`$Vm, with the RMS cloud size $`\sigma `$ and ion number $`N_{\mathrm{ion}}`$. We have independently deduced $`N_{\mathrm{ion}}`$ from measurements of the MOT fluorescence. The values of $`\sigma `$ could be determined from cuts through CCD images of the atom cloud fluorescence, yielding values of $`\sigma `$ ranging from 0.3mm to 0.5mm. For large atom clouds, we have found that the measured fields exceed the calculated ones by up to a factor of ten. As the atom cloud size is decreased, the discrepancy drops to about a factor of two. The fields deduced from the spectra are, of course, the ones to be trusted. The most likely explanation for the discrepancy and its dependence on the size of the atomic cloud is that the ions were actually not excited in the whole atomic cloud, but were concentrated in a significantly smaller volume determined by the diameters of the excitation laser pulses.

If the number of ions in a Gaussian cloud exceeds $`7.8\times 10^5\mathrm{mm}^1\sigma `$, the potential depth $`U_0`$ of the ion plasma exceeds $`0.9`$eV, and a fraction of the UV-excited photoelectrons will be retained . This is clearly the case in the lowest two spectra of Fig. 5, where we exceed the critical ion number by factors $`>2.5`$ and $`>10`$, respectively. If the fraction of retained electrons is large, a macroscopically close-to-neutral plasma core is formed, which is surrounded by a positively charged ionic shell. The electric field in the core is close to zero, while the shell carries a large electric field. The distribution $`P(E)`$ will thus become bi-modal, featuring a peak near $`E=0`$ and a high-field peak. Rydberg atoms excited in the shell produce the spectroscopic field signatures discussed in this paper, while atoms excited in the core contribute a low-field excitation spectrum with resolved Rydberg resonances up to large $`n`$. The excitation spectra shown in the two lowest plots of Fig. 5 appear, in fact, as superpositions of a high-field and a low-field spectrum, as expected for bi-modal distributions $`P(E)`$. In the lowest curve, the high-field component of $`P(E)`$ apparently is strongly dominated by the low-field one, making it hard to determine the value of the electric field in the shell. The presence of a strong electric field, can, however, be guessed from the appearance of $`p`$-lines at long wavelengths.

We note that in the case of perfect macroscopic neutrality of the core the microfields in the core would follow a Holtsmark distribution , as displayed in Fig. 4. For our charge carrier densities and laser linewidths, the Holtsmark fields are too low to be detectable.

In this work, we have used spectroscopic tools to probe the electric fields in small, laser-excited plasmas. The technique has a large future potential: using crossed and focused Rydberg atom excitation beams, it will be possible to perform spatially resolved measurements of the plasma electric fields. A variation of the time delay between the UV pulse and the blue probe pulse will allow one to study the plasma expansion via the decay of the electric field. The use of a laser with a narrow bandwidth ($`1`$MHz) will allow one to detect smaller fields, such as the Holtsmark fields, and to make precise measurements of the detailed electric-field distribution $`P(E)`$. Measurements of $`P(E)`$ could reveal inhomogeneities and instabilities in the plasma, which are expected to produce structures in $`P(E)`$ at excessive fields. Spatial order in the charge carrier arrangement, such as Wigner crystals  and the formation of crystalline Rydberg matter  may become detectable via measurements of $`P(E)`$.

Support by NSF is acknowledged. We thank Prof. P. Bucksbaum for inspiring discussions and generous loaning of equipment.


text:
# Direct Simulation of Dense Stellar Systems with GRAPE-6

## 1. Introduction

Direct $`N`$-body simulation of star clusters or other stellar systems has proven itself an extremely powerful tool to study the structure and evolution of stellar systems, since the pioneering work by von Hörner (). The only known way to do experiments on stellar systems is to construct their models in computers, because we cannot do laboratory experiments on stellar systems,

Of course, $`N`$-body simulation is not the only way to construct computer models of star clusters. One could use Monte-Carlo (Henon , for recent development see Giersz ) or direct integration of Fokker-Planck (FP) equations (Cohn , Takahashi , Drukier et al. ). In particular, recent advance in the treatment of the two-dimensional \[$`f(E,J)`$\] FP equation made it possible to use FP code to the study of the evaporation of clusters in the tidal field of the parent galaxy.

The main advantage of these methods over direct $`N`$-body simulation is the calculation cost. While it is still out of reach to perform the direct simulation of, say, a typical globular cluster with $`N=10^6`$, the FP equation is valid in the limit of $`N\mathrm{}`$. Therefore, when the assumption of large-$`N`$ limit is good, methods based on the FP approximation can deliver reliable results for the calculation cost orders of magnitudes smaller than that of direct $`N`$-body calculation. In principle, one can do the $`N`$-body simulation with small $`N`$ and then try to extrapolate that result to larger $`N`$. However, this problem, which is sometimes called “the scaling problem” has proven itself a complex and difficult problem(Heggie et al. ). The main reason of the difficulty is that there are many physical processes of which timescales depend on the number of particles in different ways, and that there is no way to adjust all relevant timescales consistently.

On the other hand, methods based on FP approximation do have their own limitations. For example, it’s pretty difficult to extend the formulation beyond the spherical symmetry. Moreover, they also suffer the problem of the time scaling, only difference is that they are affected from the opposite direction. For example, when we want to include the dynamical effect of the slowly varying potential, FP approximation goes into trouble. Consider the tidal stripping of globular clusters with non-circular orbit. The orbital timescale is of the order of $`10^8`$ years. The orbital timescale of the stars around the tidal boundary of the cluster is also $`10^8`$ years. On the other hand, the half-mass relaxation time is of the order of $`10^9`$ years. Therefore, the timestep to integrate the FP equation is at the largest $`10^8`$ years, and in practice much smaller than that. Clearly, the assumption that the dynamical timescale is smaller than the Fokker-Plank timestep is broken, and FP calculations tends to grossly overestimate the escaper rate.

As beautifully demonstrated by Takahashi and Portegies Zwart (), one can incorporate the escaping timescale into FP calculation. However, in order to do so the model parameter must be adjusted so that the agreement with reference $`N`$-body calculations is achieved. Moreover, this agreement has been so far achieved only with $`N`$-body simulations with simple static “tidal cutoff”, where no tidal field is imposed and stars outside the tidal radius are just removed. Whether or not FP calculations can reproduce the result of $`N`$-body calculations with tidal fields is an open question.

The ultimate solution for the scaling problem is to perform $`N`$-body simulations with sufficiently large number of particles. In the first half of this paper, we describe our effort in that direction, the GRAPE project to develop and use special-purpose computers for $`N`$-body simulations. In the last half of this paper, we’d like to discuss briefly our recent work on $`N`$-body simulations of galactic centers with massive black holes, where again the scaling problem shows up.

## 2. Direct simulation of dense stellar systems

In principle, the direct $`N`$-body simulation of stellar systems is very simple and straightforward. The only thing one has to do is to numerically integrate the following system of equations of motion:

$$\frac{d^2𝐱_i}{dt^2}=\underset{ji}{}Gm_j\frac{𝐱_j𝐱_i}{|𝐱_j𝐱_i|^3},$$
(1)
where $`𝐱_i`$ and $`m_i`$ are the position and mass of the particle with index $`i`$ and $`G`$ is the gravitational constant. The summation is taken over all stars in the system.

In practice, we need complex methods and tricks to integrate the above equation. Since we cannot cover all important issues here, we recommend the reviews by Spurzem () and Aarseth (\[1999a, 1999b\]). Here, we briefly discuss issues directly related to the use of special-purpose hardwares for $`N`$-body problem.

As we wrote above, in principle an $`N`$-body simulation is simple and straightforward. At each timestep, we calculate the forces on all particles in the system, and integrate their orbits using some appropriate integration method. In fact, in Molecular Dynamics simulation, where one solves the $`N`$-body problem of atoms interacting through Coulomb and van der Waals forces, one can rely on this approach. In astrophysics, however, the nature of the gravitational interaction makes such approach impractical.

The problem is that the gravitational force is an attractive force with no characteristic scale length. This fact leads to three complications. The first one is that the inhomogeneity develops as the system evolves. In the case of star clusters, this is known as core collapse or gravitational catastrophe. The central core of the cluster evolves to higher and higher density, while its mass decreases. The timestep has to be small enough to integrate accurately the orbits of particles in the core. Therefore, the timestep decreases as the cluster evolves. The second one is that even when the core density is not very high, random close encounters of two particles can lead to arbitrary short timesteps. The third one is that stable binaries and more complex hierarchical systems are formed through three-body encounters and encounters of larger number of particles (such as binary-binary interactions). These systems have very small orbital timescales. For a quantitative analysis of these issues, see Makino and Hut ().

Fortunately, we can handle these complications by a combination of techniques. The first two are solved by assigning each particle its own time and timestep. This scheme, the individual timestep scheme, was first introduced by Aarseth (), and has been used for four decades.

In the individual timestep scheme, each particle has its own timestep $`\mathrm{\Delta }t_i`$ and maintains its own time $`t_i`$. To integrate the system, one first selects the particle for which the next time ($`t_i+\mathrm{\Delta }t_i`$) is the minimum. Then, one predicts its position at this new time. Positions of all other particles at this time must be predicted also. Then the force on that particle from other particles is calculated. The position and velocity of the particle is then corrected. The new timestep is also calculated. The integration scheme is a variation of Krogh’s scheme (Krogh ) modified for second-order equations. 111Though this scheme is usually called Krogh’s scheme in the field of numerical analysis, to our knowledge the work by Wielen () is the first to discuss the implementation of individual timestep with an arbitrary order integration scheme.

A modification of this individual timestep algorithm is now used to achieve higher efficiency on vector machines (McMillan ) and special-purpose computers (Makino \[1991a\]). This scheme is called the blockstep scheme. In this scheme, the timesteps of particles are forced to integer powers of two. In addition, the new timestep of a particle is chosen so that the time of the particle is an integer multiple of the new timestep. These two criteria make it possible to force multiple particles to share exactly the same time. As a result, the efficiency of vector/parallel hardware is improved significantly. However, it should be noted that the average number of particles which share the same time is not very large, in particular when the central core is small and dense. Therefore, it is necessary that the force calculation procedure can achieve reasonable performance when the number of particles to be integrated in parallel is small. Since the number of particles in the core can be as small as 100 or less, it is necessary that the force calculation procedure can achieve a reasonable speed for that number.

The third problem, the stable binaries and other hierarchical systems, is solved essentially by integrating them as small-$`N`$ subsystem embedded in the cluster. The motion of the center of mass of the subsystem and internal motion of its members are separated, and integrated independently with different timesteps. If the perturbations from the rest of the system is small, we can ignore it. In the case of a binary, this means we can completely skip the numerical integration of the internal motion, since we know the analytic solution. If the perturbation is not negligible but is small, we can apply various techniques such as “Slow KS” (Mikkola and Aarseth ). Here, the hardest problem is not really the calculation cost. The challenge is how we can recognize these subsystems and their internal structures, and how we can decide what method is appropriate for that particular structure.

## 3. GRAPE for the simulation of dense stellar systems

We have developed a series of special-purpose computers for $`N`$-body simulations, which we call GRAPE (GRAvity PipE). Figure 1 gives the basic idea. The host computer, which is usually a general-purpose workstation running UNIX, send the positions and masses of particles to the GRAPE hardware. Then GRAPE hardware calculates the interaction between particles. What GRAPE hardware calculates is the right hand side of equation (1). When we use the Hermite scheme (Makino \[1991b\]), the first time derivative of the force must also be calculated. In this case, velocities of particles are needed. Figure 2 shows the block diagram of the pipeline unit to calculate the force and its time derivative.

In order to combine the individual timestep scheme with GRAPE hardware, one modification of the basic architecture is necessary. Figure 3 shows the change. As described in the previous section, we have to predict the position (and velocity in the case of the Hermite scheme) of all particles to calculate the forces on the particles in the current blockstep. This prediction must also be done on the GRAPE hardware, since otherwize the amount of the calculation the host computer has to do becomes too large.

In the modified architecture shown in figure 3, the particle memory keeps all data necessary for prediction, for all particles in the system. At each blockstep, the host computer writes the position and velocity of particles to be updated, and GRAPE calculates the forces on them and sends the results back to the host. If the number of pipelines is smaller than the number of particles in the block, this step is repeated until the forces on all particles in the block are obtained. Then the host performs the orbit integration using these calculated forces, and updates the data of the integrated particles in the particle memory of the GRAPE hardware.

GRAPE-4 (Makino et al. ) is the first GRAPE hardware to implement this modified architecture. In GRAPE-4, one processor board houses one predictor units and 96 (virtual) force calculation pipelines. The total system in the maximum configuration consisted of 36 boards organized into four clusters, and different boards calculated the forces on the same set of 96 particles. In this way, we met the requirement that the number of forces calculated in parallel is small, even though the number of pipelines is large.

Summation of 9 forces from processor boards in the same cluster is taken care by the communication hardware, and final summation of the forces from four clusters is handled by the host.

GRAPE-4 was completed in 1995, and has been used by many researchers for the study of dense stellar systems.

## 4. GRAPE-6

In 1997, we started the GRAPE-6 project. It’s a five-year project funded by JSPS (Japan Society for the Promotion of Science), and the planned total budget is about 500 M JYE.

The GRAPE-6 is essentially a scaled-up version of GRAPE-4(Makino et al. ), with the peak speed exceeding 100 Tflops. It will consist of around 3000 pipeline chips, each with the peak speed of 40 Gflops. In comparison, GRAPE-4 consists of 1700 pipeline chips, each with 600 Mflops. The increase of a factor of 100 in speed is achieved by integrating six pipelines into one chip (GRAPE-4 chip has one pipeline which needs three cycles to calculate the force from one particle) and using 3–4 times higher clock frequency. The advance of the device technology (from $`1\mu \mathrm{m}`$ to $`0.25\mu \mathrm{m}`$) made these improvements possible. Figure 4 shows the first sample of the processor chip delivered in early 1999. The six pipeline units are visible.

Figures 5 and 6 shows the processor board with 16 processor chips and the prototype four-board system. This four-board system has the theoretical peak speed of 2.1 Tflops, and has achieved the sustained speed of 1.1 Tflops for the simulation of 1 million-body system.

GRAPE-6 will be completed by the year 2001. We plan to make small version of GRAPE-6 (peak speed of around one teraflops) commercially available by that time. We’ve found that the commercial availability of small machines is essential to maximize the scientific outcome from GRAPE hardwares.

Compared to GRAPE-4, GRAPE-6 will give us 100 times more computer power. For simulation of star clusters for the relaxation timescale, this means a factor of five increase in the number of particles we can handle. For short simulations, the increase would be a factor of 10. In the case where we can use tree algorithms, in principle a factor of 100 increase is possible if the host computer has a sufficiently large memory. Table 1 gives a rough idea of what is currently possible with GRAPE-4 and what will be possible soon with GRAPE-6.

## 5. BH binaries in galactic cores

In this section, we briefly describe our resent work on the massive central black holes in the centers of galaxies (Makino and Ebisuzaki , Makino , Nakano and Makino \[1999a, 1999b\]). The problem is what happens when two galaxies, each with a central black hole, merge. The merging of two galaxies is a rather common event. According to the standard scenario of the structure formation in the universe (hierarchical gravitational clustering), galaxies are constructed “bottom up” from smaller objects through merging. On the other hand, almost all galaxies except for some dwarfs seem to have central black holes. Thus, if two such galaxies merge, the merger remnant would have two black holes. These black holes would settle in the center of the merger remnant, and would eventually form a binary. Two questions arise: (a) What is the structure of the central region of the merger galaxy with two black holes? and (b) What will happen to the binary itself? Will two black holes eventually merge?

Makino and Ebisuzaki () studied the first problem, assuming that black holes would eventually merge. They performed the simulation of hierarchical (repeated) merging similar to those by Farouki et al. (), but with central black holes. What Farouki et al. found is that the core radius remains almost unchanged after merging, even though the half-mass radius almost doubles after each merging event. This result is consistent with the theoretical prediction based on the conservation of the central phase space density, but in complete contradiction with the observations of luminous elliptical galaxies. Ground-based observations in 1980s demonstrated that there is strong correlation between the effective radius and the “core” radius of luminous ellipticals (Lauer ). HST observations (Gebhaldt et al. ) have shown that the “cores” are actually all cusps with the profile $`\rho =r^{0.21}`$, where $`\rho `$ is the volume luminosity density.

Figure 7 show the result of simulations with GRAPE-4. The central region of the merger with central black hole looks like the shallow cusps observers found in luminous ellipticals, and the radius of this cusp region is a constant fraction of the half-mass radius. Thus, numerical simulations of merging of galaxies with black holes reproduced the observed characteristic of the central regions of luminous ellipticals rather well.

Nakano and Makino (\[1999a, 1999b\]) investigated why the cusp is formed. First, they tried to understand the formation mechanism better by means of simplified experiments, in which just one black hole was let to sink to the center of a galaxy. Except for the case that the black hole is originally in the center, the cusp with $`\rho =r^{0.5}`$ is universally formed, and the mass within the cusp region is comparable to that of the mass of the central BH.

In the second paper (Nakano and Makino 1999b), they came up with a simple explanation for the formation mechanism. The black hole sink to the center through the dynamical friction, in the dynamical timescale. When the black hole settled at the center of the galaxies, almost no star was strongly bound to the black hole. This is because the black hole sinks to the center relatively slowly. The sinking timescale is not so slow that the adiabatic approximation can be applied to the binding energy of the stars, but sufficiently slow that the binding energy of most of the stars does not change very much. Thus, the distribution function $`f(E)`$ has a sharp cutoff at the energy close to the central potential depth of the initial galaxy model.

If there is a central massive object and $`f(E)`$ has a sharp cutoff, we can show that the central region is a cusp with the slope $`0.5`$ as follows.

The density is obtained by integrating the distribution function in the velocity space as follows:

$$\rho (r)=4\pi _{\psi (r)}^0f(E)\sqrt{2|E\psi (r)|}𝑑E,$$
(2)
where $`\psi (r)`$ is the potential at distance $`r`$ from the center, where the black hole lies. Note that we assumed that the velocity distribution is isotropic.

We assume that $`f(E)=0`$ if $`E<E_0`$. For any $`r`$ with $`\psi (r)<E_0`$, equation (2) can be rewritten as

$$\rho (r)=4\pi _{E_0}^0f(E)\sqrt{2|E\psi (r)|}𝑑E,$$
(3)
which can be expanded as

$$\rho (r)=4\sqrt{2}\pi \sqrt{|\psi (r)|}_{E_0}^0f(E)\left[1\frac{1}{2}\frac{E}{\psi (r)}+O\left(\left[\frac{E}{\psi (r)}\right]^2\right)\right]𝑑E,$$
(4)
Therefore, we have

$$\rho (r)\sqrt{|\psi (r)|}\sqrt{\frac{GM_{\mathrm{BH}}}{r}}.$$
(5)
For the fate of the black hole binary, we performed simulations with very wide range of the number of particles (2,048 to 262,144) to investigate the dependence of the growth timescale of the binding energy of the binary on the number of particles in the galaxy. What we found was intriguing. The timescale seemed to be proportional to $`N^{1/3}`$, while, theoretically, the timescale must be proportional to $`N`$, since the timescale should be limited by the timescale of filling the loss cone through two-body relaxation. We have not yet understood why our numerical experiments did not agree with the theoretical prediction. Simulations with larger $`N`$, which will be possible with GRAPE-6, will give us important clues for the understanding of this problem.

## 6. Discussion

In this paper, we briefly overviewed the computational issues of direct $`N`$-body simulations and the use of special-purpose computers, and we reviewed the plan and development status of GRAPE-6. We believe GRAPE-6 will give us many new important scientific results, much in the same way as GRAPE-4 have done so in the five years since its completion.

Since GRAPE-6 is now close to its completion, it might be worthwhile to speculate on what will come next. The silicon device technology has been advancing as usual, and will do so at least for one more decade or two. If we will start the development of the new system just after the completion of GRAPE-6, we will use the technology safely available by 2002, which is either 0.13 or $`0.09\mu \mathrm{m}`$ technology. This will give us around a factor of six increase in the transistor count. For the clock speed, we expect that the improvement by a factor of three will not be too hard, mainly because the current clock of GRAPE-6 chip is not very high. Thus, we can achieve about a factor of 20 improvement in the speed over that of GRAPE-6 chip, reaching around 0.7 Tflops per chip. With three thousand chips, the peak speed would reach two Petaflops.

The target number of particles for such system would be around 1 million, for star cluster simulations. For that number, however, one might ask whether or not the direct force calculation is really a good approach, even though it guarantees very good cost-performance ratio.

There are two widely used algorithms to reduce the calculation cost. One is the neighbor scheme (Ahmad and Cohen ), and the other is the treecode (Barnes and Hut ). For general-purpose computers, both algorithms give pretty large improvement in the overall calculation speed. The theoretical gain of the neighbor scheme is $`O(N^{1/4})`$ (Makino and Hut ), and the break-even point is as small as $`N=25`$. For the treecode, the gain is $`O(N/\mathrm{log}N)`$. The break-even point depends on the required accuracy, but $`N=10^6`$ is certainly large enough to guarantee a pretty large gain.

However, we (or anybody else) have not yet conceived a feasible way to combine the neighbor scheme and GRAPE, or treecode with individual timestep and GRAPE. In both cases, the problem is that forces on different particles are calculated differently. As a result, $`n`$ pipelines which calculate the forces on $`n`$ different particles need $`n`$ different data from the particle memory, while in the case of the simple direct summation all pipelines can use the same data. In principle, if each pipeline has its own memory, we can reconcile massively parallel GRAPE and these algorithms. With the next-generation system, it will be necessary to integrate the memory and the pipeline processor anyway, since otherwize we cannot achieve required memory bandwidth even for the case of the simple GRAPE architecture. This integrated memory and pipeline might give us the opportunity to implement more efficient force calculation algorithms.

## Acknowledgments

We would like to thank Daiichiro Sugimoto, Toshikazu Ebisuzaki, Makoto Taiji, Tomoyoshi Ito, Toshiyuki Fukushige and many others who were involved in the development of the six generations of GRAPE hardwares, and Yoko Funato, Simon Portegies Zwart, Steve McMillan, Piet Hut and again many others for discussions and collaborations in software development and scientific works on GRAPE hardwares. This work is supported by the Research for the Future Program of Japan Society for the Promotion of Science (JSPS-RFTP97P01102).