text:
# Bures Measures over the Spaces of Two and Three-Dimensional Density Matrices

## 1 Introduction

The density matrix  — having origins in early (independent) work of Landau and von Neumann — has proved to be a very useful concept in physics. Researchers have devoted substantial efforts in describing the spaces defined by density matrices , in using them to analyze the separability of quantum systems , in comparing information-theoretic properties of various probability distributions over them , as well as studying the question of parallel transport in this context . There is, of course, a quite straightforward manner in which to parameterize density matrices — simply in terms of the real and complex parts of the entries of these Hermitian matrices. A more indirect approach, but one of particular interest, relies upon Euler angle parameterizations of the special unitary matrices. It has been known for quite some time that one can parameterize density matrices with the help of such angles. To do so, we take a diagonal density matrix, $`\rho `$, which represents our quantum system in a particular basis. We then perform a unitary transformation in Hilbert space that takes $`\rho `$ to an arbitrary basis. This unitary transformation can be assumed to have a determinant equal to one since an overall phase does not affect the physics. Now when we apply the unitary transformation, we operate in the following manner,

$$\rho ^{}=U\rho U^{},$$

where $`USU(n)`$ for an $`n`$-state system. The diagonal matrix $`\rho `$ can be parameterized by its $`n1`$ eigenvalues, and the unitary transformation by $`n^21`$ variables. At this point, however, we have parameterized $`\rho ^{}`$ using too many, that is $`n^2+n2`$ parameters, since the density matrices for the $`n`$-state systems comprise only an $`(n^21)`$-dimensional convex set. A primary objective here will be to demonstrate how such an “over-parameterization” can be avoided, by the elimination (from the unitary transformation) of $`n1`$ of the set of $`n^2+n2`$ parameters.

Below, we will show that the density matrices for the 2-state and 3-state systems can be conveniently and insightfully parameterized in terms of Euler angle coordinates. In doing so, we take full advantage of the group properties associated with the unitary transformation. In particular, we obtain useful (Bures) measures on the spaces of density matrices, which can be used in the integration of functions over spaces and subspaces of quantum systems. Specifically, in Section 2.1 we exhibit the Euler angle-based parameterization of the 2-state quantum systems. In Section 2.2 we present the Euler angle-based parameterization of the 3-state systems. In Section 2.3 we discuss their straightforward generalization to $`n`$-state systems. In Section 3.1 we show how to use this to write down the Bures measures on the space of density matrices for 2-state systems and in Section 3.2 its counterpart for 3-state systems. Finally, we discuss ongoing work which is intended to yield the explicit generalization to 4-state systems. This, we anticipate will be highly useful in analyzing the (fifteen-dimensional) space occupied by pairs of qubits (cf. ).

## 2 Parametrization of Density Matrices

### 2.1 2-state Density Matrices

In this section we review certain properties of density matrices for the 2-state systems and demonstrate how to parameterize the 2-state density matrix in Euler angles. This should make the extension to 3-state systems in the following section more transparent.

Much is known about 2-state density matrices. One of the more obvious properties is that a 2-state system can be spanned by two pure state density matrices of the following form:

$$\rho _1=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right),\rho _2=\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right).$$
(1)
(These can be denoted by antipodal points on a two-sphere, making use of the elegant Riemann sphere representation .) For a mixed state diagonal density matrix, we write the linear combination of these two density matrices as

$$\rho =\rho _ia^i,$$
(2)
where $`a^i=1`$. The $`a^i`$ may be represented in several ways. One common way is to let $`a^1=a`$ and $`a^2=1a`$. What we wish to emphasize in this paper is that they are better represented by $`a^1=\mathrm{cos}^2\theta `$ and $`a^2=\mathrm{sin}^2\theta `$. Why this is true will be more obvious in the next section.

To take $`\rho `$ to an arbitrary configuration (basis) we act with a unitary transformation $`USU(2)`$, in the manner

$$\rho \rho ^{}=U\rho U^1=U\rho U^{}.$$
(3)
One may note that an $`n`$-state density matrix should have $`n^21`$ parameters, whereas this one appears to have 4, that is $`3(\text{from}SU(2))+1(\text{from the diagonal})`$. However, in the Euler angle parameterization of $`SU(2)`$, given by

$$U=e^{i\sigma _3\alpha }e^{i\sigma _2\beta }e^{i\sigma _3\gamma },$$
(4)
where the $`\sigma `$s are the Pauli matrices, we see that the parameter $`\gamma `$ drops out since $`\sigma _3`$ is diagonal and so the matrix exponential $`e^{i\sigma _3\gamma }`$ commutes with $`\rho `$ leaving precisely three parameters, as is appropriate and desired.

An important part of this parameterization is the ranges of the angles. For the 2-state case these are well-known (up to normalization given in eg., ). They are

$$0\alpha \pi ,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}\beta \pi /2,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}\theta \pi /4.$$
(5)
### 2.2 3-state Density Matrices

For 3-state density matrices we can achieve a similar parameterization. The pure states are spanned by

$$\rho _1=\left(\begin{array}{ccc}1\hfill & 0& \hfill 0\\ 0\hfill & 0& \hfill 0\\ 0\hfill & 0& \hfill 0\end{array}\right),\rho _2=\left(\begin{array}{cccc}0\hfill & 0& \hfill 0& \\ 0\hfill & 1& \hfill 0& \\ 0\hfill & 0& \hfill 0& \end{array}\right),\rho _3=\left(\begin{array}{ccc}0\hfill & 0& \hfill 0\\ 0\hfill & 0& \hfill 0\\ 0\hfill & 0& \hfill 1\end{array}\right).$$
(6)
We can then obtain a mixed state by the linear combination $`\rho _ia^i`$, with $`a^i=1`$. Here again we use the squared components of a sphere. In this case, it is a two-sphere, thus for a 3-state system we take a generic mixed state to be

$$\rho =\left(\begin{array}{ccc}\mathrm{cos}^2\theta _1\mathrm{sin}^2\theta _2& 0& 0\\ 0& \mathrm{sin}^2\theta _1\mathrm{sin}^2\theta _2& 0\\ 0& 0& \mathrm{cos}^2\theta _2\end{array}\right).$$
(7)
We then can take this to an arbitrary configuration (basis) by the process

$$\rho \rho ^{}=U\rho U^{},$$
(8)
where $`USU(3)`$. In the Euler angle parameterization of $`SU(3)`$, $`U`$ is given by

$$U=e^{i\lambda _3\alpha }e^{i\lambda _2\beta }e^{i\lambda _3\gamma }e^{i\lambda _5\theta }e^{i\lambda _3a}e^{i\lambda _2b}e^{i\lambda _3c}e^{i\lambda _8\varphi /\sqrt{3}}$$
(9)
where the Gell-Mann matrices have been used;

$$\begin{array}{cccc}\lambda _1=\left(\begin{array}{cccc}0& \hfill 1& 0& \\ 1& \hfill 0& 0& \\ 0& \hfill 0& 0& \end{array}\right),& \hfill \lambda _2=\left(\begin{array}{cccc}0& \hfill i& 0& \\ i& \hfill 0& 0& \\ 0& \hfill 0& 0& \end{array}\right),& \lambda _3=\left(\begin{array}{cccc}1& \hfill 0& 0& \\ 0& \hfill 1& 0& \\ 0& \hfill 0& 0& \end{array}\right),& \\ \lambda _4=\left(\begin{array}{cccc}0& 0\hfill & 1& \\ 0& 0\hfill & 0& \\ 1& 0\hfill & 0& \end{array}\right),& \hfill \lambda _5=\left(\begin{array}{cccc}0& \hfill 0& i& \\ 0& \hfill 0& 0& \\ i& \hfill 0& 0& \end{array}\right),& \lambda _6=\left(\begin{array}{cccc}0& \hfill 0& 0& \\ 0& \hfill 0& 1& \\ 0& \hfill 1& 0& \end{array}\right),& \\ \lambda _7=\left(\begin{array}{cccc}0& \hfill 0& 0& \\ 0& \hfill 0& i& \\ 0& \hfill i& 0& \end{array}\right),& \hfill \lambda _8=\frac{1}{\sqrt{3}}\left(\begin{array}{cccc}1& \hfill 0& 0& \\ 0& \hfill 1& 0& \\ 0& \hfill 0& 2& \end{array}\right).& & \end{array}$$
(10)
Since $`\lambda _3`$, and $`\lambda _8`$ appear on the right in (9), they commute with the diagonal matrix $`\rho `$ and drop out of the parameterization. Thus we are left with $`3^21=8`$ parameters for the 3-state case, as we should expect .

As in the 2-state case, the ranges of the angles in the parameterization are very important. For the case of three states we have

$$0\alpha ,\gamma ,a\pi ,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}\beta ,\theta ,b\pi /2,$$
(11)
$$0\theta _1\pi /4,\text{and}\mathrm{\hspace{0.33em}0}\theta _2\mathrm{cos}^1(1/\sqrt{3}).$$
(12)
Here one should note that using $`a,b`$ and $`1ab`$ for the diagonal elements of the density matrix would give rise to a domain of integration that is non-rectangular. For density matrices of higher dimensional systems this could be extremely akward. Here however, we see that the domain is indeed a rectangular solid with the final angle for a density matrix of $`n`$-dimensions having the range $`0\theta _{n1}\mathrm{cos}^1(1/\sqrt{n})`$.

### 2.3 $`n`$-state Density Matrices

Though there is now no explicit parameterization of $`SU(n)`$ in terms of Euler angles, the generalization is rather obvious. The difficulty in manipulating the one-forms to calculate the volume elements increases as $`n^2`$ so unless some substantial progress is made in getting computers to manipulate trigonometric functions symbolically, it appears that progress will be quite slow. It is anticipated that the density-matrix-parameterization program advanced here will not be achieveable for groups much higher in dimensionality than $`dim(SU(4))=15`$. Perhaps, however, with sufficient interest and effort it could be accomplished for the Lie group $`SU(8)`$, in order to effectively describe the 63-dimensional space of three entangled qubits.

## 3 Bures Measures on the Convex Sets of Density Matrices

### 3.1 2-state Systems

The parameterization of the space of density matrices by Euler angles immediately leads to a particularly simple procedure for integrating over the space. A natural (Bures) measure on the space is given by the product of the (Hall) measure \[15, eqs. (24, 25)\]  on the space of eigenvalues ($`\lambda _1,\mathrm{},\lambda _n`$), that is

$$\text{d}u=\frac{\text{d}\lambda _1\mathrm{}\text{d}\lambda _n}{(\lambda _1\mathrm{}\lambda _n)^{\frac{1}{2}}}\underset{j<k}{}4\frac{(\lambda _j\lambda _k)^2}{\lambda _j+\lambda _k}$$
(13)
and the measure on the the space of unitary matrices, the (modified/truncated) Haar measure. The geometry of the three-dimensional space is obvious from the parameterization. We have a solid section of the two-sphere and the space $`SU(2)/U(1)S^2`$. In other words, our measure is given by

$$DV=du\times dS^2,$$
(14)
where $`DV`$ is the Bures measure on the space of density matrices, $`dS^2`$ is the standard measure on $`S^2`$, ie. $`dS^2=d(G/H)`$, with $`G=SU(2)`$ and $`H=U(1)`$  and $`du`$ is given by (13) with $`n=2`$.

### 3.2 3-state Systems

The measure on the eight-dimensional space of density matrices describing 3-states is a generalization of the previous case. We have

$$DV=du\times d(G/H),$$
(15)
where $`DV`$ is the Bures measure on the space of density matrices, with $`G=SU(3)`$ and $`H=U(1)\times U(1)`$  and $`du`$ is given by (13) with $`n=3`$. The ranges of the Euler angles are given in Section 2.2. This is a new result. In the past the ranges of the angles were not properly specified when using the Haar measure. Additionally, the problem of overparameterization  was not clarified until recently .

In , it has been shown how the Bures measures for the 2- and 3-state systems can be normalized to form probability distributions.

## 4 Applications/Conclusions

The point of this note has been to clarify that: (1) the Euler angle parameterization eliminates any naive over-parameterization () of the density matrix using a unitary group representation; and (2) doing so helps in providing a natural (Bures) measure on the space of density matrices, which incorporates the (truncated) Haar measure on the group manifold of $`SU(n)`$.

So this is not “just another parameterization of the density matrix”, while proving very useful for many calculations. We find that people often do not appreciate the point of view of Feynman, who wrote that “every theoretical physicist who is any good knows six or seven different theoretical representations for exactly the same physics” .

In addition this gives rise to a rather natural (Bures) measure on the space of density matrices, using which one can integrate functions of quantum-mechanical interest (such as the von Neumann entropy, $`\text{Tr}\rho \mathrm{log}\rho `$) over spaces and subspaces of the group with relative ease. Much of the geometry of the space can then be seen in terms of the moduli of the group space , a point apparently not appreciated by a number of readers of .

## 5 Acknowledgements

M.S.B. would like to thank the DOE for partial support of this work under grant number DOE-ER-40757-123 and Southern Methodist University for the hospitality of the Group and for the use of the facilities funded by DE-FG03-95E40908. P.B.S. is grateful to the Institute for Theoretical Physics for computational support.


format: text

text:
# Answer to the comment of Chudnovsky: On the square-root time relaxation in molecular nanomagnets

\[

\]

In Ref. ,we presented a new technique, which we call the hole digging method that can be used to observe the time evolution of molecular states in crystals of molecular clusters. It allows us to measure the statistical distribution of magnetic bias fields that arise from the weak dipole fields of the clusters themselves. A hole can be ’dug’ into the distribution by depleting the available spins at a given applied field. Our method is based on the simple idea that after a rapid field change, the resulting short time relaxation of the magnetization is directly related to the number of molecules which are in resonance at the given applied field.

Prokof’ev and Stamp have suggested that the short time relaxation should follow a $`\sqrt{t}`$relaxation law . However, the hole digging method should work with any short time relaxation law . For the molecular cluster system Fe8 it is an experimental fact that the short time relaxation follows to a good approximation the $`\sqrt{t}`$relaxation law, regardless whether we start from the saturation magnetization , or from an annealed state . This is true as long as the measurements are made below 400 mK, and therefore the relaxation is purely due to quantum tunneling, and not to thermal activation .

An important result of our paper was the observation that the hole line width becomes independent of the initial value of the magnetization for small values of the later. We suggested in Ref.  that this intrinsic hole line width is directly related to the inhomogeneous level broadening due to nuclear spins as predicted by Prokof’ev and Stamp. Since the publication of our article, we have made new measurements on isotopically substituted samples of Fe8. Samples enriched with 57Fe had a larger hole line width, and samples where H is replaced with deuterium had a more narrow hole line width . These measurements confirm our hypotheses, and are in quantitative agreement which numerical simulations which takes into account the altered hyperfine coupling .

For a *saturated* sample, the Prokof’ev–Stamp theory allows us to estimate the tunneling matrix element $`\mathrm{\Delta }`$. Using Eqs. (3), (9) and (12) of , and integration, we find $`\mathrm{\Gamma }_{\mathrm{sqrt}}𝑑\xi =c\frac{\xi _0}{E_D}\frac{\mathrm{\Delta }^2}{\mathrm{}}`$, where $`c`$ is a constant of the order of unity which depends on the sample shape. With $`E_D`$ = 15 mT, $`\xi _0`$ = 0.8 mT, $`c`$ = 1 and $`\mathrm{\Gamma }_{\mathrm{sqrt}}`$ from Fig. 3 in , we find $`\mathrm{\Delta }=0.6\times 10^7`$ K which is close to the result of $`\mathrm{\Delta }=0.5\times 10^7`$ K obtained by using a Landau Zener method .

As for Mn12 in Ref  we did *not* claim that the relaxation follows the $`\sqrt{t}`$ relaxation law. It is well known that the situation in this sample is more complicated due to the fact that there are several coexisting species of Mn12 in any crystal, each with different relaxation times. In Ref  we were able to isolate one faster relaxing species and we stated in that case, that the relaxation can be approximately fit to the $`\sqrt{t}`$ relaxation law, but in fact is better fit to a power law $`t^\alpha `$ with $`0.3<\alpha <0.5`$ (depending on the applied field). We applied the hole digging method to this species, and found evidence for intrinsic line broadening below 0.3K which we suggest comes from nuclear spins in analogy with Fe8. It should be mentioned that a constant external field can only shift the internal fields. We also measured the relaxation of Mn12 at higher temperature (0.04 - 5 K) and found no evidence what so ever for a short time $`\sqrt{t}`$ relaxation.

Finally, we emphasize that the measurements of short time relaxation allows us to study the time evolution of a well-defined initial state, whereas the interpretation of long time relaxation data are far more difficult due to the development of complex intermolecular correlations during the relaxation process, which are not yet well understood .
W. Wernsdorfer1, C. Paulsen2, R. Sessoli3
1Lab. Louis Néel and 2CRTBT
CNRS, BP 166, 38042 Grenoble Cedex 9, France
3Dept. of Chemistry, University of Firenze,
via Maragliano 72, 50144 Firenze, Italy


text:
# Untitled Document

GAMMA-RAY BURSTS AS A COSMIC WINDOW FOR GALAXY EVOLUTION

I.F. Mirabel1,2, D.B. Sanders3, E. Le Floc’h1,4

1CEA/DSM/DAPNIA/SAp. Centre d’Etudes de Saclay. 91191. Gif/Yvette. France

2Intituto de Astronomía y Física del Espacio. Conicet. Argentina

3Institute for Astronomy, University of Hawaii, 2680 Woodlawn Dr., Honolulu, HI 96822

4European Southern Observatory. Casilla 19001, Santiago 19. Chile

Abstract: Present knowledge indicates that gamma-ray bursts are linked with massive stars. They will become invaluable probes of the early universe and galaxy formation. In the future, it will be possible to use gamma-ray bursts for two purposes: 1) to probe the history of massive star formation in the Universe by the rate of occurence of gamma-ray bursts, and 2) for the study of galaxy evolution at all lookback times by determining the nature of the galaxy hosts. Because gamma-rays are not attenuated by intervening dust and gas, the selection of the cosmic sites of massive star formation by this method is less affected by the biases associated with optical-uv surveys (e.g. UV-dropout techniques). Infrared and sub-millimeter follow up studies of the hosts of gamma-ray bursts may: 1) reveal a putative population of reddened ($`RK4`$) galaxies at high redshifts, and 2) detect very massive stars (population III) formed at $`z`$ 5.

Dusty starbursts and AGNs at high redshift

Optical surveys with the Hubble Space Telescope and large ground based telescopes have enabled in recent years the study of luminous UV galaxies out to redshifts z$``$4 (Madau et al. 1996; Steidel et al. 1998). From the early analysis of galaxies in the HDF, Madau et al. (1996) concluded that there seemed to be little evidence for substantial amounts of dust obscuration in high redshift galaxies. It is remarkable how quickly this picture has changed after the recent deep field surveys in the mid-infrared, far-infrared, and submillimiter wavelengths (e.g. Smail et al. 1997; Hughes et al. 1998; Barger et al. 1998; Blain et al. 1999; Elbaz et al. 1999). These new surveys have discovered a strongly evolving population of Luminous Infrared Galaxies (see Sanders & Mirabel 1996, for a review) that could be the progenitors of the present-day population of massive spheroidal galaxies.

The importance of dust at high redshifts has also been pointed out by Fabian (1999), who proposes a large population of dust enshrouded AGNs at $`z>1`$ to account for the X-ray background. This prediction has recently been confirmed by the discovery with Chandra of a large population of X-ray point sources with “extremely faint, or in some cases undetectable” (i.e. $`I>26`$) “ optical counterparts” (Mushotzky et al. 2000). Furthermore, at millimeter and submillimeter wavelengths thermal emission from dust is currently detected from quasars up to redshifts of 5 (e.g. Omont et al. 1996; Carrili et al. 2000). The implied far-infrared luminosities are $`>10^{13}L_{}`$, dust masses $`>10^8M_{}`$, and molecular gas masses of a few $`\times 10^{10}M_{}`$. These major events in galaxy evolution at high redshifts are largely missed by UV/optical surveys, and it is now clear that the use of such surveys alone may lead to a subtantial underestimate of the global star formation rate and certainly to a distorted picture of the early evolution of the most massive galaxies. What may be the fraction of massive star formation in the universe that occured in galaxies so heavily obscured by dust that could not be detected in UV-selected surveys is an issue of current debate (e.g Adelberger & Steidel 2000; Sanders 2000).

Gamma-ray bursts: A new cosmic window for galaxy evolution

Another major timely event in astrophysics has been the dramatic tranformation in our understanding of gamma-ray bursts (GRBs). There is increasing evidence that the most common GRBs, those with durations longer than a few seconds, are at cosmological distances and are associated with sites of massive star formation. The physical properties of the afterglows, their locations at a few kpc from the center of host galaxies (e.g. Djorgovski 1999), and the statistics from the several thousands of GRBs detected so far with BATSE (Fishman et al. 1999), give strong support to the idea that the majority of GRBs are linked to the cataclysmic collapse of massive stars into black holes (Paczyński 1998; MacFayden & Woosley 1999). In this context, GRBs can be used as sign-posts to probe more accurately the history of massive star formation and galaxy evolution.

One of the great advantages of this approach is that GRBs can be detected at all lookback times. The broad band spectral and temporal behavior of the afterglows have confirmed that the bursts are beamed in relativistic jets, with the observer near the jet axis (Mészáros 1999; Kulkarni et al. 1999; Galama et al. 1999; Castro-Tirado et al. 1999). This leads to observed energies that are Doppler boosted by several orders of magnitude (Mirabel & Rodríguez 1999). Therefore, the detection of GRBs depends essentially on the beaming angle rather than on their distance, and they can be observed up to very large redshifts. An example of this was the optical afterglow of GRB990123, which could be seen with a simple pair of binoculars, since it reached an optical brightness of 8.6 mag, despite the fact that it came from z=1.6. High resolution spectroscopy of such optically bright bursts, which are many magnitudes brighter than quasars at the same redshift, could be very valuable to probe with unprecedented sensitivity the Ly$`\alpha `$ forest, and consequently, the chemical evolution of the universe.

The nine GRBs with redshifts determined by optical techniques prior to the end of 1999 are in the redshift range $`z=0.433.47`$, with a mean of $`z=1.3`$ (see also Blain & Natarajan 1999). However, this mean redshift is almost certainly biased toward low values because sources at very early epochs appear rather red and it is difficult to determine their redshifts by optical observations alone. For instance, no \[OIII\] or Ly$`\alpha `$ lines from objects at $`1.3z2.5`$ can be seen in the optical. Furthermore, there are several GRB afterglows that appeared to be extremely red whose redshifts are unknown. Some of these could be related to the dead of a putative class of Population III stars, which are very massive stars brought on by cooling of molecular hydrogen in the dark ages of the very high redshift universe, before galaxies were formed (Lamb & Reichart 1999; Bromm, Coppi & Larson 1999). One of the best cases to explore this possibility is GRB980329, for which Fruchter (1999) hypothesized $`z5`$.

Dust obscuration is very important in studies of recently formed compact objects (e.g. Mirabel et al. 1999). Because GRBs are the last phase of the evolution of the most massive stars (Woosley 1999) which do not live long enough to leave their place of birth, it is expected that a large fraction will be still enshrouded in their placental clouds of molecular gas and dust. The observations of GRB afterglows are consistent with this picture: 1) among the 47 X-ray afterglows with good localizations, about half had no optical counterparts (Greiner 1999), 2) despite persistent optical follow ups, three radio afterglows without optical counterparts have been detected, and 3) some afterglows appear to be extremely reddened (e.g. GRB990705). Dust obscuration is the most likely reason for the non detection of optical transients. Indeed the absence of detectable optical flux accompanying strong X-ray emission in the 1-60 keV energy band of SAX may simply be due to the strong redshift dependence of dust obscuration at optical wavelengths in comparison to the X-rays, i.e. $`\tau _{\mathrm{opt}}`$/$`\tau _{\mathrm{X}\mathrm{ray}}`$ $``$ $`(1+z)^4`$ (Taylor et al. 1998). Therefore, GRBs can be used to detect a putative population of dust-enshrouded forming galaxies at high redshifts.

The future

After the breakthrough produced by the X-ray satellite Beppo-SAX (Costa et al. 1997), it is expected that in the comming years, high energy space missions (SAX, HETE2, CHANDRA, XMM, INTEGRAL and SWIFT) will be providing about 50 GRB X-ray afterglows per year with localizations in the range of 1 arcsec to 1 arcmin. Furthermore, the localizations of the gamma-ray bursts can be improved a posteriori using the Interplanetary Network.

We note that the selection of the initial sample of star-formation sites in this approach is less affected by the problems and biases associated with the optical-UV selections of star formation sites, which are strongly influenced by reddening effects. Indeed the gamma rays are not attenuated by the intervening columns of gas and dust. The bursts thus provide an unbiased sampling of cosmic sites of massive star-formation.

An important caveat is that most X-ray afterglows until now have had error boxes of a few arcmin2, a region large enough to contain from tens up to a few hundred galaxies, therefore it is very difficult without afterglow detections in the optical, infrared or radio wavelengths to identify unambigously which object is the host of the GRB. However, the science programs of the recently launched Chandra and XMM, which carry sophisticated grazing incidence mirrors mounted aboard sophisticated satellites, include ToO observations of X-ray GRB afterglows with subarcminute error boxes. These X-ray telescopes can provide X-ray positions with arcsecond error boxes. A case in point is the recent Chandra observation of GRB000210 on the basis of the sole GRB error box provided by the Beppo-SAX wide field camera (Garcia et al. 2000 in GCN 544).

Studies of the galaxy hosts of GRBs with future instruments (e.g. SIRTF, FIRST, NGST, ALMA, etc) will open a new cosmic window for galaxy formation and evolution.

Conclusion

Gamma-ray bursts can be used as beacons to:

1) study the history of massive star formation in the universe,

2) to detect the putative population III stars at redshifts $``$ 5,

3) to detect a population of dust-enshrouded massive galaxies in formation,

4) to determine the chemical evolution of the universe.

Besides large optical telescopes, this approach will require infrared and sub-millimeter instruments of unprecedented sensitivity.

Acknowledgements: I.F.M. acknowledges support from Conicet/Argentina, and thanks Jacques Paul and Paolo Goldoni for enlightening discussions.

References

Adelberger, K.L. & Steidel, C.C. 2000, ApJ, in press (astro-ph 0001126)

Barger, A.J., et al. 1998, Nature, 394, 241

Blain A.W., et al. 1999, MNRAS, 302, 632

Blain, A.W. & Natarajan, P. 1999, in press (astro-ph/9911468)

Bromm, V., Coppi, P.S. & Larson, R.B. 1999, ApJ, 527, L5

Carrili, C.L., et al. 2000, ApJ, in press (astro-ph 0002386)

Castro-Tirado A., et al. 1999, Science, 283, 2069

Costa E., et al. 1997, Nature, 387, 783

Djorgovski S.G., et al. 1999, AAS, 194, 414

Elbaz D., et al. 1999, A&A, 351, 37

Fabian, A.C., et al. 1999, MNRAS, 308, L39

Fishman G.J. 1998, in The Hot Universe, eds. Koyama et al., IAU S188, 159

Fruchter A.S. 1999, ApJ, 512, L1

Galama T.J., et al. 1999, Nature, 398, 394

Greiner J. 1999, GRB web page http://www.aip.de/~jcg

Hughes, D., et al. 1998, Nature, 394, 241

Kulkarni S.R., et al. 1999, Nature, 398, 389

Lamb, D.Q. & Reichart, D.E. 1999, (astro-ph/9909002)

Madau P., et al. 1996, MNRAS, 283, 1388

MacFayden, A.I. & Woosley, S.E. 1999, ApJ 524, 262

Mészáros P. 1999, Nature, 398, 368

Mirabel I.F. & Rodríguez L.F. 1999, ARAA, 37, 409

Mirabel, I.F., et al. 1999, V Huntsville Conference on Gamma-Ray Bursts, (astro-ph/9912446)

Mushotzky, R.F., Cowie, L.L., Barger, A.J. & Arnaud, K.A. 2000, Nature, in press

Paczyński B. 1998, ApJ, 494, L45

Sanders, D.B. 2000, Galaxy Morphology Conf., ed. Block et al., (astro-ph/9910028)

Sanders D.B. & Mirabel I.F. 1996, ARAA, 34, 749

Smail I., Ivison R.J. & Blain A.W. 1997, ApJ, 490, L5

Steidel C.C., et al. 1999, ApJ, 519, 1

Taylor G.B., et al. 1998, ApJ, 502, L115

Woosley, S.E. 1999, V Huntsville Conference on Gamma-Ray Bursts. (astro-ph/9912484)


text:
# An Ecolab Perspective on the Bedau Evolutionary Statistics

## Introduction

At Alife VI, Mark Bedau proposed some evolutionary statistics (Bedau et al., 1998) as a means of classifying different evolutionary systems. The intent here is to find a general scheme analogous to Wolfram’s (1984) classification scheme of cellular automata. Three statistics are proposed:

The number of species or components in the system

Activity of a species is defined as the population count of that species, the vector $`𝒏`$ in Ecolab terms. Evolutionary activity subtracts from this the neutral or nonadaptive part. This is achieved by running a neutral shadow model, that is identical with the original model, except that natural selection must be “turned off”. Finally, this activity is accumulated over the lifetime of the species, and then averaged over all species.

This corresponds the the number of new species crossing a threshold, divided by the diversity.

Bedau describes four classes of evolutionary behaviour, as in the following table:

| Class | $`D(t)`$ | $`\overline{A}_{\text{cum}}(t)`$ | $`A_{\text{new}}(t)`$ | Description |
| --- | --- | --- | --- | --- |
| 1 | bounded | zero | zero | none |
| 2 | bounded | unbounded | none | unbounded, uncreative |
| 3 | bounded | bounded | positive | bounded, creative |
| 4 | unbounded | positive | positive | unbounded, creative |

Note that in Bedau et al. (1998), only 3 classes are mentioned — class 2 was added later in his presentation at Alife VI. Bedau has applied his statistics to a number of artificial life models, including Echo (Holland, 1995) and Tierra (Ray,1991), none of which exhibit class 4 behaviour. By contrast, the same statistics applied to the fossil record (at least for the Phanerozoic — the period of time since the appearance of multicellular life in the Cambrian) — show a strong class 4 behaviour. Further, Bedau speculates that the global economy and internet traffic are also class 4, particularly as they show strong growth over a significant period of time. Since no artificial life systems to date appear to show class 4 behaviour, the gauntlet has been laid down to discover such a system to work out whether this classification difference is fundamental or not.

Ecolab (Standish, 1994), whilst not an artificial life system, is a model of an evolving ecology that has advantages of mathematical tractability and computational simplicity. It lies in between the extremely simplistic models of (for example) Bak and Sneppen (1993) or Newman (1997) and artificial life models of evolution such as Tierra or Avida. One of its key characteristics is that its dynamics are defined by the ecological interactions between the species, rather than ad hoc exogenous dynamics. The Bedau statistics are well defined for it, so it is interesting to see what class behaviour Ecolab has. Furthermore, an Ecolab-like model is possible for all artificial life systems (valid in a continuum limit). For example, the equations of motion for Tierra are given in Standish(Standish, 1997).

## Ecolab

The Ecolab model (as opposed to the Ecolab simulation system) is based on an evolving Lotka-Volterra ecology. The defining equation is given by:

$$\dot{𝒏}=𝒓𝒏+𝒏𝜷𝒏+𝜸^2𝒏+𝝁(𝒓𝒏)$$
(1)
where $`𝒏`$ is the species density, $`𝒓`$ the effective reproduction rate (difference between the intrinsic birth and death rates in the absence of competition), $`𝜷`$ the matrix of interaction terms between species, $`𝜸`$ the migration rate and $`𝝁`$ the mutation operator. All of these quantities (apart from $`𝜷`$, which is a matrix) are vectors of length $`n_{\mathrm{sp}}`$, the number of species in the ecology. The operator $``$ denotes elementwise multiplication. The mutation operator returns a vector of dimensionality greater than $`n_{\mathrm{sp}}`$, with the first $`n_{\mathrm{sp}}`$ elements set to zero — in effect expanding the dimensionality of the space, a key feature of this system. For a more detailed exposition of the various properties of the model, in particular, the precise form of the mutation operator, the reader is referred to the previous published papers, as well as the Ecolab Technical Report, which are all available from the Ecolab Web Site111http://parallel.hpc.unsw.edu.au/rks/ecolab.html.

For the purposes of this paper, it is worthwhile expounding a little on the properties of the mutation operator. It models point mutations in particular (other mutation types, such as recombination are simply not modeled within Ecolab). Point mutations in genotype space, which satisfy Poisson statistics, give rise random mutations, with locality, in phenotype space. Since the only phenotypic properties of interest to the model are the parameters $`𝒓`$, $`𝜷`$ and $`𝜸`$ the parameters are mutated according to a normal or lognormal distribution (according as the parameters are reals or positive (or negative) respectively), using a sample from the Poisson distribution for the width. The two parameters governing mutation (width of the Poisson distribution, and the rate at which mutations are attempted) are related via a simple proportional factor (called the “species radius (or separation)”) that is kept constant throughout the simulations reported here. Each species has its own mutation rate — given as a vector $`𝝁`$.

Each of these phenotypic parameters are initialised from a uniform distribution. The relevant input parameters for a run are then maximum and minimum values for each of $`𝒓`$, the diagonal of $`𝜷`$, the offdiagonal of $`𝜷`$, $`𝝁`$, $`𝜸`$ and the species radiua $`\rho `$. The complete system may be scaled in the time dimension, fixed by what value is chosen for the timestep. In this case, $`\mathrm{max}_ir_i=0.1`$, so one timestep corresponds to about a 14th of the doubling time of the fastest reproducing organism in the ecology. This is a compromise between continuity of the simulation and computational expense. The ratio $`\frac{\mathrm{max}_ir_i}{\mathrm{max}_i\beta _{ii}}`$ roughly corresponds to the carrying capacity of the ecology. This is chosen to about 100 so that behaviour near the equilibrium is reasonably continuous rather than stochastic. The ratio of offdiagonal to diagonal terms relates to how negative definite $`𝜷`$is. Since mutations tends to drive the matrix away from being negative definite (system stability), the maximum of the offdiagonal terms is chosen to make the initial system marginally unstable. The species radius $`\rho =0.1`$ was chosen empirically to make new species phenotypically distinct from its parent species.

Having fixed the other parameters according to the above criteria, the remaining degrees of freedom are $`𝝁`$ and $`𝜸`$. In this paper, we vary the maximum mutation rate in different simulations, but keep the distribution of migration rates fixed.

One other feature worth noting is that the mutation operator will also randomly add or drop connections between species, according to an exponential distribution. Thus, the mutation operator is in fact highly conservative — with the lognormally mutated parameters capped (in the case of $`𝝁`$ and $`𝜸`$) or restricted by the requirements of boundedness (diagonal components of $`𝜷`$)(Standish, 1998; Ecolab Technical Report).

## Neutral Shadow Model

An important feature for improving the accuracy of the evolutionary statistics is the use of a neutral shadow model. This model should be as similar as possible to the original model, but with all selection turned off. In the case of Ecolab, this is accomplished by running a shadow population density vector $`𝒏`$’, and when $`𝒏`$ is updated, the shadow vector is updated by a random permutation of the updates. Thus each shadow species behaves in the long run like an average species. Activity is also tracked at the same time, with the activity vector being updated by the difference between the population density and the shadow population density, provided that difference is positive.

The new activity statistic $`A_{\text{new}}`$ is computed by summing the number of species that have crossed a threshold. In (Bedau et al., 1998), this threshold is determined by plotting the activity distributions for both the original and the shadow model, and taking the cross-over point as the threshold. This turned out to be 50 individuals, rather than the arbitrary 10 individuals used in other Ecolab studies. In fact the two distributions are nearly equal over the range 10–50, but if an activity is above 50, then it is highly likely to be due to adaptive behaviour.

## Behaviour of Ecolab

Figure 1 shows the Bedau statistics for typical Ecolab runs (panmictic, or spatially independent case), as a function of mutation rate. When the mutation rate is too low, class two behaviour is seen. Diversity remains constant, and activity grows unbounded as the system rapidly sheds unviable organisms and tends to a stable ecology. Conversely, for high levels of mutation, class one behaviour is seen. There is a constant churn of organisms, that do not have any chance to generate activity. For intermediate levels of mutation, an interesting situation arises. Here, the number of mutant organisms that successfully invade the ecosystem roughly balances the number lost through extinction (Standish, 1998). Scale free behaviour is observed in a number of statistics, such as the distribution of species lifetime. These same 3 states of behaviour have been observed in Avida (Adami et al., 1998).

The code used for this simulation is available from the Ecolab web site as version 3.3 of the software. The model including the neutral shadow model is defined in shadow.cc, and a sample experimental script given as bedau.tcl. The only parameters varied are the spatial dimensions and mutation(random,maxval).

The evolutionary statistics were also collected for a spatially dependent Ecolab, however due to some implementation difficulties, run lengths exceeding $`1\times 10^6`$ timesteps have not been achieved prior to this paper’s deadline. Broadly speaking, though, the same behaviour is seen as the panmictic case, although there is a period of diversity growth in the early period prior to settling on a higher level of diversity than the panmictic case.

This can be understood by considering two extremes of spatially dependent Ecolab models, namely zero migration and infinite migration. Infinite migration effectively corresponds to the panmictic case again, whereas zero migration corresponds to a number of cells, independent of each other, each running the panmictic model. So we would expect in the case of zero migration, the diversity (in the long run) should be proportional to the number of cells (or the total area). The in between case of finite nonzero migration should also show an increase in diversity with area, due to partial independence of each cell, but the increase should be sublinear, as migration causes some species to be identified between cells. Island Biogeography (MacArthur and Wilson, 1967) theory postulates that the relationship is $`DA^s`$ for some coefficient $`s`$, which presumably must depend in some fashion on the migration rates, but is generally in the range 0.2–0.35 for most empirical studies.

## May’s Stability Criterion

May (1972) proposed that random Lotka-Volterra webs would be unstable if

$$n_{\mathrm{sp}}<\frac{1}{s^2C}$$
(2)
where $`C`$ is the connectivity, defined as the proportion of nonzero elements in $`𝜷`$, and $`s`$ is the interaction strength, defined as the standard deviation of the offdiagonal terms of $`𝜷`$, divided by the average of the diagonal terms. Cohen and Newman (1985) showed that May’s criteria does not hold for Lotka-Volterra systems in general, only a smaller class related to the models May studied. However, the inverse relationship between species number and connectivity does appear to hold (Pimm, 1982;Cohen and Newman, 1988;Cohen et al., 1990).

Stability is not a relevant property in Ecolab, as really the persistent state (which includes the stable state as a special case) is the attractor. However, the inverse relationship between diversity and connectivity does hold(Standish, 1998), for spatially dependent as well as panmictic cases. Therefore, in order for diversity to show an increasing trend, a corresponding decreasing trend must occur in connectivity. This ought to be true of the biosphere also, given the universality of this relationship.

As mentioned in section Ecolab, the mutation operator is highly conservative with respect to connectivity. It assumes that a new species inherits the same connections as its parent, with random additions or deletions according to a symmetric distribution (just as likely to gain a connection as lose one). This has the effect of preserving the connectivity over time. In order for connectivity to decrease, different dynamics would need to be proposed, for example assuming that the mutant species did not compete with its parent.

One possibility for the cause of this growth in diversity is the mass extinctions, that have occurred a handful of times throughout the Phanerozoic. However, the only reasonable way of modeling this is to remove a random proportion of species from the ecology at a particular time. This operation does not alter the connectivity, as the links lost is exactly balanced by the reduced diversity. When implemented within Ecolab, one gets the characteristic rebound in diversity after the extinction event, however, the rebound is back to about the same diversity level as existed prior to the extinction.

Another possibility that actually would work in the right way is related to the fact that the Phanerozoic era corresponds to the breakup of the Pangaea supercontinent — firstly into Gondwana and Laurasia, then into the six continents we know today. Assuming that there is almost no migration between the continents (thus 6 equal-sized continents would support 6 times the diversity of one continent that size) and that the species-area law within a continent has $`DA^{.3}`$, we would expect that a breakup of a single supercontinent into 6 equal sized pieces should produce $`6^{1.3}=3.5`$ times the diversity of the original supercontinent. This factor accounts for a significant fraction of the diversity growth since the Permian.222In case anyone thinks that this result is an argument in favour of habitat fragmentation for promotion of diversity, this is a question of scale. Over short timescales habitat fragmentation is bad for diversity, as is any major environmental change. Only over evolutionary timescales will the diversity bounce back.(Benton, 1995)

Clearly this is a very rough “back of the envelope” calculation, but it is sufficient to show that continental breakup needs to be allowed for in determining if there is any intrinsic evolutionary processes driving diversity growth.

## Acknowledgements

I wish to thank Mark Bedau for many illuminating discussions, and for assistance in developing the neutral shadow model for Ecolab. I also wish to thank the New South Wales Centre for Parallel Computing for computational resources required for this project.


text:
# Numerical Simulations of Convective Accretion Flows in Three Dimensions

## 1 Introduction

Analytical and numerical calculations have shown that advection-dominated accretion flows (ADAFs) experience strong convection, especially for low values of the viscosity parameter $`\alpha `$ (Narayan & Yi 1994, 1995; Igumenshchev, Chen & Abramowicz 1996; Igumenshchev & Abramowicz 1999, 2000; Stone, Pringle & Begelman 1999). The numerical work reveals that convection significantly modifies the structure of ADAFs (Stone et al. 1999; Narayan, Igumenshchev & Abramowicz 2000, hereafter NIA). The change in the structure is understood to be a consequence of inward transport of angular momentum by convective eddies (NIA; Quataert & Gruzinov 2000, hereafter QG; Ryu & Goodman 1992; Stone & Balbus 1996).

All the numerical simulations done so far have corresponded to axisymmetric two-dimensional flows. NIA and QG emphasized that the inward transport of angular momentum seen in these simulations is a consequence of axisymmetry (cf. Stone & Balbus 1996) and that the results could in principle be qualitatively very different in three dimensions. We present in this paper the first fully three-dimensional numerical simulations of ADAFs. The details of our numerical technique are given in §2 and the results are described in §3. The paper concludes with a discussion in §4.

## 2 Numerical Method

We solve the non-relativistic time-dependent Navier-Stokes equations of an accretion flow in the gravitational potential of a central black hole:

$$\frac{d\rho }{dt}+\rho \stackrel{}{v}=0,$$
$`(2.1)`$
$$\rho \frac{d\stackrel{}{v}}{dt}=P+\rho \mathrm{\Phi }+𝚷,$$
$`(2.2)`$
$$\rho \frac{de}{dt}=P\stackrel{}{v}+Q,$$
$`(2.3)`$
where $`\rho `$ is the density, $`\stackrel{}{v}`$ is the velocity, $`P`$ is the pressure, $`\mathrm{\Phi }=GM/r`$ is the Newtonian potential of the central mass $`M`$, $`e`$ is the specific internal energy, $`𝚷`$ is the viscous stress tensor (we include all components), and $`Q`$ is the viscous dissipation function. Since we assume an ADAF, there is no radiative cooling. We adopt the ideal gas equation of state, $`P=(\gamma 1)\rho e`$, where $`\gamma `$ is the adiabatic index, and write the kinematic shear viscosity coefficient as $`\nu =\alpha c_s^2/\mathrm{\Omega }_K`$, where $`c_s=(P/\rho )^{1/2}`$ is the isothermal sound speed and $`\mathrm{\Omega }_K=(GM/r^3)^{1/2}`$ is the Keplerian angular velocity. Equations (2.1)–(2.3) are independent of the black hole mass $`M`$ when the radius $`r`$ is scaled by the gravitational radius $`r_g=2GM/c^2`$ and time is scaled by $`r_g/c`$.

Each time step of the numerical computation is split into two sub-steps: hydrodynamical and viscous. The hydrodynamical sub-step is calculated by using the PPM algorithm (Colella & Woodward 1984), and the viscous sub-step is solved by using a standard explicit algorithm for the viscous contribution to the equation of motion (2.2) and the energy equation (2.3).

We solve equations (2.1)–(2.3) on a 3D Cartesian grid consisting of a sequence of six nested sub-grids. Each sub-grid $`k`$ is $`26\times 26\times 52`$ in size and is uniform, extending over the range ($`0,R_k^{out}`$), ($`0,R_k^{out}`$), ($`R_k^{out},R_k^{out}`$) in $`x`$, $`y`$, $`z`$, respectively. The innermost grid has $`R_1^{out}=26r_g`$ and each successive $`R_k^{out}`$ is equal to $`2R_{k1}^{out}`$, such that the outermost sub-grid extends to $`R_6^{out}=832r_g`$. Figure 1 illustrates the geometry of the grid. Due to numerical limitations we compute the flow only over a quarter of the full cubic volume and apply azimuthal periodic boundary conditions in the $`(x,z)`$ and $`(y,z)`$ planes. At the outer boundary $`R_6^{out}`$ we use a free outflow boundary condition, and in the innermost sub-grid we assume that matter drains freely inside a sphere of radius $`r3r_g`$. The inner boundary condition mimics the effect of the last stable orbit around a Schwarzschild black hole.

The time-step $`\mathrm{\Delta }t`$ for the numerical integration is chosen in accordance with the Courant condition for the hydrodynamical sub-step. For this $`\mathrm{\Delta }t`$, the viscous integration is numerically stable because we choose sufficiently low values for the viscosity parameter $`\alpha `$. We assume that mass is injected into the computational domain over an equatorial slender torus with radius $`r_{inj}800r_g`$. The matter is injected at a fixed rate $`\dot{M}_{inj}`$ and has rotational velocity (with respect to the $`z`$-axis) equal to $`0.95`$ times the local Keplerian velocity. As a result of viscous interactions, part of the injected matter moves inwards and forms an accretion flow, while the rest leaves the computational domain through the outer boundary. We follow the time evolution of the flow on time scales longer than the viscous time scale at $`r_{inj}`$. This allows the flow to attain steady state in a time-averaged sense. Note that when $`\rho `$ is scaled by $`\dot{M}_{inj}`$, the results are independent of the value of $`\dot{M}_{inj}`$. Therefore, the simulations depend on only two dimensionless parameters: $`\alpha `$, $`\gamma `$.

## 3 Results

We have calculated four models (A, B, C, D) which differ in the choice of $`\alpha `$ and $`\gamma `$ (see Table 1). All the models show complicated time-dependent convective motions, analogous to those described in detail by Igumenshchev et al. (1996), Igumenshchev & Abramowicz (1999, 2000) and Stone et al. (1999). The amplitude of the variability is not large, and the average properties of the flow are well represented by time-averaging. To good accuracy, the time-averaged flow is axisymmetric in all the models.

Figure 2 shows the time-averaged radial profiles in the equatorial plane ($`z=0`$) of $`\rho `$, $`\mathrm{\Omega }`$ and $`c_s`$ for the four models. In the inner regions, $`r10^2r_g`$, the profiles are well approximated by power-laws: $`\mathrm{\Omega }r^{3/2}`$, $`c_sr^{1/2}`$, $`\rho r^\beta `$, where the index $`\beta `$ is $`0.5`$ for Models A and C and $`0.8`$ for Models B and D. (For $`r>2\times 10^2r_g`$, the flows are strongly affected by the outer boundary at $`r8\times 10^2r_g`$ and are not power-law in character.) The scalings of $`\mathrm{\Omega }`$ and $`c_s`$ and that of $`\rho `$ in Models A and C are identical to those predicted by the analytical work of NIA and QG. The $`\beta 0.8`$ scaling of Models B and D is a little surprising, and suggests that $`\beta `$ depends on the value of $`\gamma `$. A similar dependence of $`\beta `$ on $`\gamma `$ was found also in the axisymmetric 2D simulations of Igumenshchev & Abramowicz (2000).

Figure 3 shows a snapshot of the flow in the equatorial plane for Model C. The other three models are qualitatively similar. The most important feature is that the density and velocity fluctuations are almost axisymmetric, as seen by the fact that the eddies are significantly elongated in the azimuthal direction. This suggests that the present 3D results are not likely to be very different from those obtained with 2D axisymmetric simulations. The minimum length-scale of perturbations in the radial $`r`$ and polar $`\theta `$ directions depends on the viscosity parameter $`\alpha `$, being larger for larger values of $`\alpha `$.

The left panel in Fig. 4 shows a snapshot meridional cross-section of Model C in the $`(y,z)`$-plane. The shading shows the distribution of the vorticity vector, $`\stackrel{}{\omega }=\mathrm{rot}\stackrel{}{v}`$. Regions with a positive $`x`$-component of the vorticity $`\omega _x`$ are shown in grey, and regions with negative $`\omega _x`$ are shown in white. It is seen that the flow is dominated by radial streams of inflowing and outflowing material which are sandwiched in the $`\theta `$-direction. The streams are chaotic and move around so that they are washed out when a time average is taken. Models A and B are very similar to Model C, except that they have $`23`$ inflowing streams at a given moment, whereas Model C has usually $`56`$ streams.

Model D is qualitatively different from the other three models, as shown by the right panel in Fig. 4. What is plotted here is the time-averaged flow pattern and vorticity. We see that the model has quasi-stationary structures which survive even after averaging over about 100 periods of the Keplerian rotation at $`r=100r_g`$. The flow pattern is almost symmetrical with respect to the equatorial plane, and shows the presence of a stable and regular arrangement of circulation cells in the meridional plane.

The direction of transport of angular momentum is given by the sign of the $`(r\phi )`$-component of the Reynolds stress tensor, $`v_r^{}v_\phi ^{}`$. NIA showed that this quantity is negative in 2D simulations of convective ADAFs, but noted that it was a consequence of the axisymmetry of the simulations and that real 3D convection might behave differently. In Figure 5 we show the radial dependence of the Reynolds stress averaged over $`\theta `$ for the four 3D models described in this paper. We see that, except for a small boundary region, the stress is negative over the bulk of the flow in all the models. This confirms that the 2D results published in previous papers are physically meaningful. Figure 5 shows that the magnitude of the Reynolds stress is positively correlated with the values of both $`\alpha `$ and $`\gamma `$. Model D has a surprisingly weak stress. This may be connected to the quasi-stationary meridional circulation seen in this model.

## 4 Discussion and Conclusions

We have used numerical simulations to study three-dimensional hydrodynamic accretion flows with inefficient radiative cooling (ADAFs). The simulated flows have nearly Keplerian angular momentum on the outside, relatively weak viscosity, $`\alpha =0.01`$, $`0.03`$, and an adiabatic index $`\gamma `$ equal to either $`5/3`$ or $`4/3`$. We find that the convective motions are nearly axisymmetric (Fig. 3). As a result, the basic properties of the flows are very similar to those found in previous 2D axisymmetric simulations (Igumenshchev et al. 1996; Igumenshchev & Abramowicz 1999, 2000; Stone et al. 1999; NIA).

Our most important result is that convection transports angular momentum inward (Fig. 5). This property of convection was identified by Ryu & Goodman (1992) in analytical work on thin accretion disks. Recently, NIA showed that convection in 2D axisymmetric simulations of ADAFs also moves angular momentum inward. They and QG argued, however, on the basis of previous work by Stone & Balbus (1996), that this result was a consequence of enforcing axisymmetry on the accreting gas. The present work shows that even in 3D, where the flow is not constrained to be axisymmetric, convective motions do transport angular momentum inward. In fact, we find that 2D and 3D simulations show good quantitative agreement in the estimate of the ($`r\phi `$)-component of the Reynolds stress tensor, and in the radial dependences of various fluid variables.

In the models, the turbulent pressure due to convection, $`P_{turb}=\rho v^2/3`$, where $`v^{}`$ is the velocity fluctuation, is $`10^3`$ of the gas pressure $`P`$. Therefore, convection does not have a direct dynamical effect on the flow. However, as argued by NIA and QG, convection can have a very strong indirect effect on the structure of the flow if it causes inward transport of angular momentum. The accretion flow can then achieve a critically balanced state in which the convective transport almost completely cancels the outward transport of angular momentum by viscosity. This leads to a flat radial density profile $`\rho r^{1/2}`$. Since the amount of transport possible with convection is limited, the viscosity parameter $`\alpha `$ must be less than a critical value $`\alpha _{crit}0.1`$ for the new structure to be realized (NIA). The value of $`\alpha `$ in an ADAF is not known at present. Model-fitting to the observed light curve and spectra of soft X-ray transients suggests that $`\alpha 0.20.3`$ (Esin, McClintock & Narayan 1997). Numerical MHD simulations in 3D should provide an independent estimate of $`\alpha `$ (assuming that viscosity in ADAFs is produced by magnetic stresses). If such computations confirm that $`\alpha `$ in ADAFs is large, then convection may turn out to have only limited dynamical effect on these flows. Conversely, if $`\alpha <\alpha _{crit}`$, one may expect ADAFs to resemble the analytical solutions described by NIA and QG.

Our results may be relevant also to the physics of rapidly rotating convective stars. In such stars convection may reduce or even suppress the outward transport of angular momentum by magnetic stresses (see Bisnovatyi-Kogan et al. 1979).

An important property of convective accretion flows is the outward transport of energy by convection. We find that 2D and 3D simulations give similar results; the effective outward luminosity due to convection is $`10^310^2\dot{M}c^2`$ for the values of $`\alpha `$ and $`\gamma `$ considered here, where $`\dot{M}`$ is the net mass accretion rate of the black hole. Very likely the energy will be radiated on the outside of the accretion flow. If so, the results have important implications for the spectra and luminosities of accreting black holes and neutron stars (Igumenshchev & Abramowicz 2000). Convection could also introduce variability on different time-scales (Igumenshchev & Abramowicz 1999).

The simulated flows presented in this paper show some interesting dependence on the adiabatic index $`\gamma `$. Models with $`\gamma =4/3`$ have a steeper density profile (Fig. 2) than models with $`\gamma =5/3`$. In addition, Model D ($`\alpha =0.01`$, $`\gamma =4/3`$) shows a quasi-stationary meridional circulation (Fig. 4), whereas other models have chaotic convective motions. The reason for these differences is unclear.

This study was supported by a grant from the Royal Swedish Academy of Sciences and by grants PHY 9507695 and AST 9820686 from the National Science Foundation.


text:
# III Zw 2, the first superluminal jet in a Seyfert galaxy

## 1 Introduction

The radio properties of quasars with otherwise very similar optical properties can be markedly different. There is a clear dichotomy between radio-loud and radio-quiet quasars in optically selected samples. The radio-loudness is usually characterized by the radio-to-optical flux ratio. In the PG quasar sample, which is probably the best studied quasar sample in the radio and optical (Kellermann et al. kellermann (1989), Boroson & Green boroson (1992)), radio-loud and radio-weak quasars separate cleanly in two distinct populations (e.g. Kellermann et al. 1989).

It is known that radio-loud AGN almost never reside in late type, i.e. spiral galaxies (e.g. Kirhakos et al. kirhakos (1997), Bahcall et al. bahcall (1995)) whereas radio-quiet quasars appear both in spiral and in elliptical host galaxies. Furthermore, all relativistically boosted jets with superluminal motion and typical blazars have been detected in early type galaxies (e.g. Scarpa et al. scarpa (1999)). It is still unclear, why AGN in spiral galaxies, at the same optical luminosity as their elliptical counterparts, should not be able to produce the powerful, relativistic jets seen in radio galaxies.

However, a few sources with intermediate radio-to-optical ratios appear to be neither radio-loud nor radio-quiet. They form a distinct subclass with very similar radio morphological and spectral properties. They all have a compact core at VLA scales and a flat and variable spectrum in common. These properties are very similar to the ones of radio cores in radio-loud quasars, but their low radio-to-optical ratio and their low extended steep-spectrum emission is atypical for radio-loud quasars. Miller et al. (miller (1993)) and Falcke et al. (falcke95 (1995), 1996a &b) have identified a number of these sources, called “radio-intermediate quasars” (RIQs), and suggested that they might be relativistically boosted radio-weak quasars or “radio-weak blazars”. This would imply that most, if not all, radio-quiet quasars also have relativistic jets. In fact, VLBI observations of radio-quiet quasars already have shown high-brightness temperature radio cores and jets (Blundell & Beasly blundell (1998)). A crucial test of the relativistic jet hypothesis is the search for apparent superluminal motion in these sources. A prime candidate for detecting this is the brightest radio source in the RIQ sample, III Zw 2, which we discuss in this paper.

III Zw 2 (PG 0007+106, Mrk 1501, $`z=0.089`$) is one of the most extremely variable radio sources and a very unusual AGN. It was discovered by Zwicky (zwicky (1967)), classified as a Seyfert I galaxy (e.g., Arp arp (1968); Khachikian & Weedman khachikian (1974); Osterbrock osterbrock (1977)), and later also included in the PG quasar sample (Schmidt & Green schmidt (1983)). The host galaxy was classified as a spiral (e.g. Hutchings & Campbell hutchings\_campbell (1983)) and a spiral arm was claimed (Hutchings hutchings (1983)). A disk model was later confirmed by fitting of model isophotes to near-IR images (Taylor et al. taylor (1996)).

The most interesting property of III Zw 2, however, is its extreme variability at radio and other wavelengths with at least 20-fold increases in radio flux density within 4 years (Aller et al. aller (1985)). The source also shows optical (Lloyd lloyd (1984)) and X-ray variability (Kaastra & de Korte kaastra (1988); Pounds pounds (1986)).

III Zw 2 is a core-dominated flat-spectrum AGN with only a faint extended structure (see Unger et al. unger (1987)). The weak extended radio emission and the host galaxy is quite typical for a Seyfert galaxy. Its \[OIII\] luminosity is a mere factor three brighter than that of a bright Seyfert galaxy like Mrk 3 (e.g. Alonso-Herrero et al. alonso (1997)) which explains why it has been classified as either a Seyfert galaxy or a quasar. In this luminosity region a distinction between the two may not be of much significance.

Earlier VLBI observations of the source have only shown a high-brightness temperature core (Falcke et al. 1996b , Kellermann et al. kellermann98 (1998)) and recent Millimeter-VLBI observations by Falcke et al. (falcke99 (1999)) just barely resolved the source into two very compact components.

Based on its average optical-to-radio ratio of $`200`$ and its radio properties, it was suggested that III Zw 2 could be a radio-intermediate quasar and the presence of superluminal motion in this source was predicted. Here we will discuss multi-epoch VLA and millimeter-VLBI observations of this source confirming this hypothesis.

## 2 Observations and Data Reduction

In 1997 we detected the onset of a new major radio outburst in III Zw 2 and we initiated a target of opportunity program to monitor the spectral evolution of the burst with the VLA and its structural evolution with the VLBA.

The VLA observations were made at six frequencies ranging from 1.4 GHz to 43 GHz in A, B, C and D configuration from 1998 September until now. The three epochs discussed in this paper where obtained on 1998 November 04, 1999 March 23 and 1999 July 07 in CnB, D and A configuration respectively. The source 3C48 was used as the primary flux density calibrator, and III Zw 2 was self-calibrated and mapped with the Astronomical Image Processing System (AIPS). Since the 1.4 GHz observation in March 1999 was heavily confused by the nearby sun, we estimated the flux density by interpolation between earlier and later epochs.

We observed III Zw 2 with the VLBA on 1998 February 16, 1998 June 13, 1998 September 22, 1998 December 12, and 1999 July 15 at 43 and 15 GHz. We used a total bandwidth of 64 MHz for a full 8 hour scan, except for the last observation which was a 4 hour scan. We spent three-quarters of the available observing time at 43 GHz and one quarter of it at 15 GHz. For the second epoch, we used the Effelsberg 100 m telescope in combination with the VLBA. We reduced the data using the software packages AIPS and DIFMAP (Shepherd, Pearson, & Taylor shepherd (1994)). Fringes were detected in the III Zw 2 data on all baselines. We calibrated the gains using system temperature information and applied atmospheric opacity corrections. To obtain a reliable total flux estimate, amplitude gains for stations with bad weather conditions were scaled up to match the other antennas. The data were then self-calibrated, first using phase-only and later phase-amplitude self-calibration with solution intervals slowly decreasing down to one minute. A more detailed discussion of the VLA spectra and the 15 GHz VLBA-data will be presented in a forthcoming paper.

## 3 Results

The initial spectrum of III Zw 2 from 1998 May was presented in Falcke et al. (falcke99 (1999)). It was highly inverted at centimeter wavelengths with a spectral index of $`\alpha =+1.9\pm 0.1`$ between 4.8 and 10.5 GHz. The entire outburst spectrum from 1.4 GHz to 666 GHz could basically be fitted by only two homogeneous, synchrotron components which are optically thin at high frequencies and become self-absorbed below 43 GHz. This spectral turnover frequency stayed constant until 1998 November and hence we expected no strong structural change during this time (Fig. 1).

Our first three VLBA epochs were made while the spectral peak was constant at 43 GHz. The core itself is resolved at 43 GHz at all epochs. To represent the extent of the source, the non-zero closure phases at long baselines for the 43 GHz data were fit by two point-like components. A rough estimate of the formal statistical errors of the component separation was obtained by dividing the original beam size by the post-modelfit signal-to-noise ratio (e.g. Fomalont fomalont (1999), Sect. 2.3). The errors were of the order of $`1\mu as`$ for the first four epochs and $`10\mu as`$ for the fifth epoch. Additional to this very small statistical error, there should be a larger systematic error which is difficult to quantify. To minimize this systematic error we used very similar reduction procedures for each epoch.

In accordance to the VLA data, the maps of the first three epochs show no structural change and the separation of the fitted components stayed constant within the statistical errors at $``$76$`\pm 2\mu `$as corresponding to $``$0.11 pc (see Table 1) for an angular size distance of 307.4 Mpc ($`H_0=75`$ km/sec/Mpc, $`q_0=0.5`$ as used in this paper).

After 1998 November, the VLA observations showed a dramatic change in the spectrum. The spectral peak dropped quickly to 15 GHz within a few months (Fig. 1). Since the peak in this source is caused by synchrotron self-absorption (Falcke et al. falcke99 (1999)), the fast change in peak frequency implied a similarly strong morphological change, i.e. a rapid expansion. To roughly estimate the expected expansion speed, we applied a simple equipartition jet model with a $`R\nu _{\mathrm{ssa}}^1`$ dependence (e.g. Blandford & Königl blandford\_koenigl (1979); Falcke & Biermann falcke\_biermann (1995)). For an initial source size $`R=0.11`$ pc and a self-absorption frequency $`\nu _{\mathrm{ssa}}=43`$ GHz in 1998 November, we calculate a source size of 0.32 pc for a self-absorption frequency of 15 GHz in 1999 March. Thus we predicted an apparent expansion speed of 1.9 c after the correction for cosmological time dilatation and asked for further VLBA-observations.

Indeed the fifth epoch of VLBA observations showed a dramatic structural change compared to the earlier epochs (Fig. 2). A model of at least three point-like components is required now. It was not possible to fit the closure phases with a two-component model as in the earlier epochs or to get rid of this third component during self-calibration. The separation of the outer components for all five epochs is plotted in Fig. 3. While the separation at the first three epochs is consistent with an expansion speed of $`0.04c`$, the fifth epoch shows a rapid expansion. The apparent expansion speed between the outer components in the 4th and 5th epoch is $`1.25\pm 0.09c`$.

This value is only a lower limit and increases to 2.66 c if one considers the time range from December until March during which most of the spectral evolution occurred. Applying the standard equation for superluminal motion, $`\beta _{\mathrm{app}}=\frac{\beta \mathrm{sin}\theta }{1\beta \mathrm{cos}\theta }`$ (e.g. Krolik krolik (1999)), to a value of $`\beta _{\mathrm{app}}=2.66`$ constrains the maximal angle between the jet and the line of sight to $`\theta =41^{}`$, since $`\beta <1`$.

## 4 Summary & Discussion

We found a close connection between the spectral and structural evolution of the radio outburst in III Zw 2. While the self-absorption frequency remained constant, we observed no change on VLBA-scales. The excellent consistency of the relative sizes during the first three epochs on the few micro-arcsecond level demonstrates the good quality of our data. The quick drop in peak frequency after 1998 November marked the beginning of a strong structural change. The spectral peak dropped from 43 GHz to 15 GHz within a few months and the VLBA-maps show a rapid expansion with an apparent expansion velocity of 1.25 c. The fact that spectral and spatial evolution are so closely linked also demonstrates that we are dealing with real physical expansion.

With its faint extended radio structure and its spiral host galaxy III Zw 2 is clearly not a radio-loud quasar, but has properties very typical of luminous Seyfert galaxies or radio-quiet quasars. The detection of superluminal motion in this galaxy now clearly shows that we are dealing with a relativistic jet on sub-pc scales. To our knowledge this is the first detection of superluminal motion in a spiral galaxy with a Seyfert nucleus. The maximum aspect angle of $`41^{}`$ from the superluminal motion (see Sec. 3) is in good agreement with the orientation based unified scheme (e.g. Antonucci antonucci (1993)) for AGN where Seyfert I galaxies are seen under intermediate or small aspect angles, so that the nucleus is not obscured by a dusty torus.

For the question of the nature of the radio-loud/radio-quiet dichotomy this means that radio-weak and radio-loud quasars can indeed have central engines that are in many respects very similar. Their optical properties are almost indistinguishable and both types of quasars can produce relativistic jets in their nuclei. The finding of superluminal motion supports the hypothesis of Miller et al. (1993) and Falcke et al. (1996a ) that RIQs are relativistically boosted intrinsically radio-weak AGN.

While these general conclusions seem to be fairly robust, some characteristics of the outburst need to be modeled in more detail. The initial phase of the flux density rise with its millimeter-peaked spectrum and no detectable expansion perhaps has to be explained similar to the physics of Gigahertz-Peaked-Spectrum (GPS) sources, i.e. ultra-compact hotspots pumped up and powered by a jet interacting with the interstellar medium or the torus. The rapid expansion thereafter could have marked the phase where the jet breaks free and starts to propagate relativistically into a lower-density medium. Another explanation of the initial phase could be a fast jet moving through quasi-stationary components as proposed for sources like 4C 39.25 (e.g. Alberdi et al. alberdi (1993)). In any case, this dramatic structural change should go together with a change of the polarization vector. This must be checked by future experiments.

The initial slow expansion has also possible implications for the interpretation of other Seyfert galaxies. For example, in observations of the Seyfert galaxies Mrk 348 and Mrk 231 (Ulvestad et al. ulvestad (1999)) only sub-relativistic expansion was found. Hence, one could raise the question whether the region of relativistic expansion in these two sources is at even smaller scales or whether they just happened to be in a ‘slow’ phase, similar to III Zw 2 early on. If true, it is still possible that Seyfert jets are launched relativistically but are slowed down and disrupted significantly already on the sub-parsec scale. Therefore it would be important to follow the spectral evolution of these sources in more detail.

###### Acknowledgements.

We acknowledge the help of an anonymous referee. This research was supported by DFG grants Fa 358/1-1 and 358/1-2. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.


text:
# EXCITATION OF OSCILLATIONS IN THE MAGNETIC NETWORK ON THE SUN

## 1 INTRODUCTION

It is now well established that the solar photosphere is threaded with strong (kilogauss) magnetic fields in the form of vertical flux tubes located in the network, where they are observed as bright points. Observations have revealed that these network bright points (NBPs) are in a highly dynamical state due to the buffeting action of granules (e.g. Muller 1983; Muller et al. 1994; Berger & Title 1996; Berger et al. 1998; van Ballegooijen et al. 1998).

The jostling of magnetic elements by granules can excite oscillations in flux tubes, which can in principle be an important mechanism for chromospheric and coronal heating, provided that the motions are rapid enough (Choudhuri, Auffret & Priest 1993; Steiner, Knölker & Schüssler 1994; Hasan & Kalkofen 1999, hereafter Paper I). As transverse waves, excited in the photosphere, travel upwards their velocity amplitude increases. In chromospheric layers, the velocity amplitude becomes comparable to the tube speed for kink waves, leading to an efficient coupling with longitudinal waves (Kalkofen 1997). The latter can dissipate by forming shocks (e.g. Zhugzhda, Bromm & Ulmschneider 1995).

In Paper I it was argued that network oscillations can be efficiently excited through the buffeting action on flux tubes by external granules with intermittent motions (of about 2 km s-1) occurring on a time scale of less than about half the cutoff period of kink waves. Typically, a period of 1 min. is consistent with the general result. We assumed an analytic form for the external motions. In the present work, we use observational data obtained at the Swedish Solar Observatory at La Palma during 1995 (Berger et al. 1998) to infer the horizontal footpoint motions of magnetic elements in the photosphere. We model the response of flux tubes to these motions and calculate the energy flux in upward-propagating transverse (kink) waves. We also simulate the effect of hypothetical motions on time scales of a few seconds, which could significantly enhance the energy flux in kink waves. Finally, we examine some implications of our results.

## 2 MEASUREMENT OF BRIGHT POINT MOTION

The data were collected on October 5, 1995 between 10:57 and 12:08 UT, simultaneously in the G band (4305 Å) and the nearby continuum (4686 Å), using the phase-diversity method (see Löfdahl et al. 1998; Berger et al. 1998). After filtering out $`p`$mode oscillations, a time series of 180 images was obtained, with high spatial resolution and with an interval between frames of 23.5 s. Following Berger et al., we subtract the images in the G band and 4686 Å, which clearly highlights the network bright points.

The motion of the bright points was followed using a tracking technique with “corks” as tracers of bright points (van Ballegooijen et al. 1998, hereafter Paper II). The corks, each with a radius of 60 km (the pixel size of the current data set), were advected by an artificial flow field that is proportional to the gradient of the intensity in the magnetic image, and their positions were followed in time (for details see Paper II). In Figure 1 we show the cork positions superposed on a typical frame of a magnetic difference image, where $`x`$ and $`y`$ denote the horizontal coordinates on the Sun. Figures 2a and 2b show the $`x`$ and $`y`$ components of the velocity as a function of time for two representative corks, whose initial locations are indicated by arrows in Figure 1.

## 3 KINK WAVE EXCITATION DUE TO FOOTPOINT MOTION

Let us assume that the motion of the G band bright points can be taken as a proxy for the footpoint motion of flux tubes at the base of the photosphere ($`z=0`$). Furthermore, for simplicity, we assume that flux tubes are isothermal and “thin” compared to the pressure scale height $`H`$. It is convenient to work in terms of the “reduced” transverse displacement $`Q_{}(z,t)=\xi _{}e^{z/4H}`$, where $`\xi _{}`$ is the Lagrangian displacement. Once the footpoint motion is specified, the velocity ($`\dot{Q}_{}`$) at any height $`z`$ and time $`t`$ can be determined by the following expression (Paper I):

$$\dot{Q}_{}(z,t)=\dot{Q}_{}(0,tz/c_\kappa )k_\kappa z_0^{tz/c_\kappa }\dot{Q}_{}(0,t_0)\frac{J_1(\omega _\kappa \zeta _\kappa )}{\zeta _\kappa }dt_0,$$
(1)
where $`J_1`$ is a Bessel function, $`\zeta _\kappa =\sqrt{(tt_0)^2(z/c_\kappa )^2},`$ $`c_\kappa `$ and $`\omega _\kappa `$ are, respectively, the propagation speed and cutoff frequency for kink waves, given by:

$$c_\kappa ^2=\frac{2}{\gamma }\frac{c_s^2}{1+2\beta },$$

$$\omega _\kappa ^2=\frac{g}{8H}\frac{1}{1+2\beta };$$

$`c_s`$ is the sound speed, $`k_\kappa =\omega _\kappa /c_\kappa `$, $`\gamma =5/3`$, $`\beta =8\pi p/B^2`$, $`p(z)`$ is the gas pressure inside the tube and $`B(z)`$ is the magnitude of the vertical component of the magnetic field on the tube axis. The energy flux ($`F_\kappa `$) in a single flux tube is given by (Paper I):

$$F_\kappa =\frac{2p_{\mathrm{e}^,0}}{\beta +1}\dot{Q}_{}\left(\frac{Q_{}}{z}+\frac{1}{4H}Q_{}\right)e^{z/2H},$$
(2)
where $`p_{e,0}`$ denotes the external gas pressure at $`z=0`$.

Figures 3a and 3b show the vertical energy flux in transverse waves versus time at a height $`z=750`$ km for two observed magnetic elements located initially at the coordinates shown above the figure. The left and right panels correspond to footpoint motions with the $`x`$ and $`y`$ components of the velocity, respectively.

We find that the injection of energy into the chromosphere takes place in brief and intermittent bursts, lasting typically 30 s, separated by longer periods (longer than the time scale for radiative losses in the chromosphere) with lower energy flux. The peak energy flux into the chromosphere is as high as $`10^9`$ erg cm-2 s-1 in a single flux tube, although the time-averaged flux is $`10^8`$ erg cm-2 s-1. This scenario for heating the magnetic network would produce strongly intermittent chromospheric emission consisting of brief intense flashes superimposed on generally a very low background. However, in dense network regions the observed chromospheric emission is stably present and exhibits relatively low amplitude, long-period variations (e.g. Damé, Gouttebroze and Malherbe 1984; von Uexküll et al. 1989; Lites, Rutten and Kalkofen 1993). Therefore, the above scenario cannot readily explain the observed persistence of emission from dense network regions. A possible remedy is to consider the effect of high frequency motions.

## 4 SIMULATING HIGH-FREQUENCY MOTIONS

The observations of G band bright points were taken with a time cadence of 23.5 s. However, the possibility of random footpoint motions occurring on a shorter time scale cannot be ruled out. The solar convection has very high Reynolds number and is expected to be highly turbulent, involving a wide range of length and time scales. Indeed, spectroscopic observations of the solar granulation have shown enhanced broadening of spectral lines in intergranular lanes, indicating enhanced turbulence within these region (e.g. Nesis et al. 1992, 1994). Magnetic flux tubes interact with these turbulent flows within the lanes. We speculate that these interactions produce transverse displacements of magnetic field lines on length scales much less than the flux tube diameter ($`100`$ km), i.e., the magnetic field inside a flux tube has fine structure on transverse scales $`\mathrm{}100`$ km. These perturbations likely generate Alfvén waves which propagate upward along the field lines and dissipate their energy higher up in chromosphere.

The velocity at small spatial scale can be estimated as follows. Assuming the turbulence has a Kolmogoroff spectrum, the velocity of eddies on length scale $`\mathrm{}`$ is given by $`v(\mathrm{})v_0(\mathrm{}/\mathrm{}_0)^{1/3}`$, where $`\mathrm{}_0`$ is the outer scale of the turbulence, and $`v_0`$ is the velocity of the largest eddies. We assume that $`\mathrm{}_0`$ is given by the width of an intergranular lane as determined from 3D simulations of the solar granulation (e.g., Stein and Nordlund 1998), and that $`v_0`$ is the convective velocity predicted by such simulations: $`\mathrm{}_0`$ 100 km and $`v_0`$ 2.5 km s-1. This turbulence model predicts that $`v`$ 1 km s-1 on a length scale $`\mathrm{}`$ 10 km, i.e., there is significant power in motions at this scale. The time scale of such motions is $`\mathrm{}/v`$ 10 s, which is significantly shorter than the time scale for granular buffeting.

In the following we simulate the effect of such high-frequency motions on the response of the flux tube (similar to the study of longitudinal wave excitation in flux tubes by Ulmschneider and Musielak 1998). For simplicity we ignore the small transverse scale of the motions and we describe the ”turbulence” in terms of footpoint displacements of the flux tube. We do this by superposing random motions, with a rms amplitude of 1 km s-1 and a zero mean value, on the observationally determined velocity. The random motions vary on a time scale of 2.35 s. Such high-frequency motions (if present) would not have been detected in the present observations. Hence, it is not possible a priori to rule them out. Results from such simulations are shown in Figures 4a and 4b for the $`x`$ component of the velocity as a function of time for two representative corks whose initial locations are indicated above each panel, and in Figures 5a and 5b for the corresponding vertical energy wave flux. We note that the energy flux in upward-propagating transverse waves is much larger and less intermittent than in the case without high-frequency footpoint motions.

The energy flux in the form of high-frequency transverse waves appears adequate to continuously heat the chromosphere. Taking a filling factor of 10% at a height $`z=750`$ km above the photosphere, the predicted flux is approximately $`2\times 10^7`$ erg cm-2 s-1, which is sufficient to balance the observed radiative loss of the chromospheric network ($`10^7`$ erg cm-2 s-1, see Model F in Avrett 1985)

## 5 DISCUSSION

The energy flux in upward propagating transverse waves has been shown for a few representative magnetic elements. We should point out that we carried out similar calculations for several corks, before selecting two as being typical of the large sample that we examined.

We find that the injection of energy into the chromosphere occurs in brief and intermittent bursts (lasting typically 30-60 s), separated by longer periods (longer than the typical time scale for radiative loss) with low energy flux. The peak energy flux into the chromosphere is typically $`10^9`$ erg cm-2 s-1 in a single flux tube. Observationally, such a scenario for heating of the magnetic network on a 1 arc sec scale would lead to high variability in Ca II K emission, contrary to observations. However, the observations of G band bright points used here do not have the time resolution and sensitivity necessary to detect motions on shorter time scales. We simulate such high-frequency motions by adding random velocity fluctuations with a zero mean value and an rms amplitude of 1 km s-1. We conclude that for transverse waves to provide a viable mechanism for sustained chromospheric heating, the main contribution to the heating must come from high-frequency motions (periods 5-50 s).

The proposed high-frequency motions have very small spatial scale (10 - 100 km) and the transverse displacements involved are also very small (for example, a motion with period 50 s and velocity amplitude of 1 km s-1 has a displacement amplitude of only 8 km). Such length scales are well below the angular resolution of present-day solar telescopes. Therefore, it is unlikely that such high-frequency motions will ever be detected in studies of the proper motions of photospheric flux tubes. Spectroscopic methods, such as measurements of spectral line broadening by unresolved turbulent flows, offer greater promise but will still require very high angular resolution.

The observational data used in this study were provided by R.A. Shine of Lockheed-Martin and were obtained at the Swedish Vacuum Solar Telescope. We acknowledge support from the National Science Foundation and NASA. SSH is thankful to the Smithsonian Institution for travel support under SMIN grant no. 9411FR00007. We are grateful to the referee for helpful comments.


text:
# Comment on: rotational properties of trapped bosons

## Abstract

Based on the Hellman-Feynman theorem it is shown that the average square radius of a cloud of interacting bosons in a parabolic well can be derived from their free energy. As an application, the temperature dependence of the moment of inertia of non-interacting bosons in a parabolic trap is determined as a function of the number of bosons. Well below the critical condensation temperature, the Bose-Einstein statistics are found to substantially reduce the moment of inertia of this system, as compared to a gas of “distinguishable” particles in a parabolic well.

Herewith we repost our paper cond-mat/9611090 (1996). It was published in Phys. Rev. A 55, 2453 (March 1997), three years before cond-mat/0003471 (2000) by Schneider and Wallis.


text:
# Quasiperiodic time dependent current in driven superlattices: distorted Poincaré maps and strange attractors

## Abstract

Intriguing routes to chaos have been experimentally observed in semiconductor superlattices driven by an ac field. In this work, a theoretical model of time dependent transport in ac driven superlattices is numerically solved. In agreement with experiments, distorted Poincaré maps in the quasiperiodic regime are found. They indicate the appearance of very complex attractors and routes to chaos as the amplitude of the AC signal increases. Distorted maps are caused by the discrete well-to-well jump motion of a domain wall during spiky high-frequency self-sustained oscillations of the current.

Nonlinear dynamics of weakly coupled semiconductor superlattices (SLs) driven by dc and ac biases has been the research topic of many experimental and theoretical works . The nonlinearity manifests itself in many physical situations as for instance, in their transport properties under finite DC or AC bias. Under appropriate DC voltage bias, the current through the SL displays natural oscillations caused by creation, motion and annihilation of domain walls (DWs) in the SL . The observed oscillation frequencies range from kHz to GHz. Superimposed on a smooth current oscillation, there appear faster current spikes whose frequencies are typically in the GHz regime . Each spike reflects the well-to-well motion of DWs causing the self-oscillation, and therefore their frequency is related to the characteristic tunneling time. Theory and numerical simulations show that DW may be charge monopoles or dipoles , although experimental evidence shows that DW are charge monopoles for self-oscillations observed in the available samples . AC driven SLs display a rich dynamical behavior including quasiperiodic and chaotic oscillations with nontrivial spatial structure .

Many studies fix the frequency of the AC drive as the golden mean number $`(1+\sqrt{5})/21.618`$ times the frequency of the natural oscillations (i.e., the frequency ratio is an irrational number hard to approximate by rational numbers), which is convenient to obtain complex dynamical behavior . In this case, the system presents a rich power spectrum, a complex bifurcation diagram and different routes to chaos including quasiperiodicity, frequency locking, etc. . First return or Poincaré maps are used to analyze unambiguously the underlying attractors . In the quasiperiodic case, Poincaré maps usually consist of smooth loops, whereas they are a set of discrete points in the case of frequency locking. More exotic Poincaré maps resembling distorted double loops in the quasiperiodic case have been experimentally observed in middle of the second plateau of the I–V characteristic of a SL . At the onset of this plateau, Poincaré maps are smooth and not distorted. The origin of distorted maps was not understood at the time of their observation, although disorder and sample imperfections were invoked . Luo et al  showed that a combination of signals with different frequency was needed in order to reproduce experimentally observed distorted double layer Poincaré maps. The origin of this combination was not ascertained in that work.

The aim of the present work is to establish theoretically that high-frequency current spikes of the self-oscillations give rise to these exotic Poincaré maps. In turn, current spikes are due to the well-to-well motion of DW during each period of the self-oscillations. Thus distorted Poincaré maps reflect DW motion in AC driven SLs.

Theoretical studies of nonlinear effects in weakly coupled SL are based in microscopic modelization of the sequential tunneling current between adjacent wells . Early models were discrete drift models with a phenomenological fitting of tunneling current and boundary conditions . More sophisticated models have included the electrostatics at the contacts in a self-consistent framework. Depending on SL configuration and contact doping, we know that undriven self-oscillations are due to recycling of either monopole DWs or dipole waves . In the latter case, a dipole wave consists of two DWs, one corresponding to a charge accumulation layer and another one corresponding to a depletion layer. Between these DWs there is a high-field region encompassing several wells, and the whole dipole wave travels through the SL. Recycling of such dipole waves gives rise to current self-oscillations similar to the well-known Gunn effect in bulk semiconductors .

In this paper we have analyzed the tunneling current through a DC+AC biased SL by means of the microscopic model described in . We have considered a 50-period SL consisting of 13.3 nm GaAs wells and 2.7 nm AlAs barriers, as described in . Doping in the wells and in the contacts are $`N_w=2\times 10^{10}`$ cm-2 and $`N_c=2\times 10^{18}`$ cm-3 respectively. With these doping values, self-oscillations are due to recycling of monopole DWs. We choose not to analyze the original 9/4 sample where the experiments were performed . In this sample, the second plateau of the I–V characteristics ends at a resonant peak due to $`\mathrm{\Gamma }X`$ tunneling. In the simpler SL configuration chosen here, the $`X`$ point does not contribute to the peaks at the first and second plateaus. We analyze the sequential tunneling current with an applied voltage, $`V(t)=V_{AC}(t)+V_{DC}`$, where $`V_{AC}(t)=V_{AC}\mathrm{sin}(2\pi f_{AC}t)`$. The applied AC frequency $`f_{AC}`$ has been set to the golden mean times the natural frequency and is very small compared with typical energy scales of the system. Thus the AC potential modifies adiabatically the potential profile of the SL and our sequential tunneling model holds .

Fig. 1 shows the evolution of the current through the SL, its Fourier spectrum and its Poincaré map for $`V_{DC}=4.2`$ V and $`V_{AC}=19`$ mV. These values correspond to the onset of the second plateau of the I–V characteristic curve. The current trace of Fig. 1(a) is quasiperiodic and does not present observable superimposed high-frequency oscillations (spikes). The natural oscillation near the onset of the plateau is caused by monopole recycling very close to the collector contact. Thus the DW does not move over many wells and the current trace does not present an appreciable number of spikes. In the power spectrum of Fig. 1(b) there are contributions coming from the fundamental frequency $`f_0`$=39 MHz, the frequency of the applied AC field $`f_{AC}`$, the combination of both and their higher harmonics. The Poincaré map depicted in Fig. 1(c) is a smooth loop with a nontrivial double layer structure indicating quasiperiodic oscillations. If we increase the DC bias up to the middle of the second plateau, $`V_{DC}=5.1`$ V, the undriven oscillation is due to recycling of monopole DWs which propagate periodically through part (approximately $`40\%`$) of the structure and disappear at the collector. A similar motion is observed in the AC driven case. The corresponding current trace may show high frequency spikes depending on the chosen initial field profile. A flat initial field profile gives a spiky current trace until the latter settles to the stable oscillation (over which there are no appreciable spikes); see the earlier part of Fig. 2(a). The power spectrum and Poincaré maps corresponding to this case and to the previous one (Figs. 1(b) and (c)) are markedly different. Fig. 2(b) shows that the power spectrum of the case with spikes has a greater harmonic content than that of the case without spikes (Fig. 1), presenting many peaks corresponding to higher harmonics of the fundamental frequencies (i.e., the low fundamental frequency oscillations, $`f_0=33`$ MHz, and the spikes, $`500`$ MHz), the applied frequency and combinations thereof. The corresponding Poincaré map consists of a distorted loop with a double layer structure which shows a strong similarity with previous experiments ; see Fig. 2(c). From these numerical observations, we conclude that distorted Poincaré maps are linked to spiky current traces, even if such traces change to smooth ones after a transient.

The previous conclusion may be reinforced if we change the emitter doping so that undriven self-oscillations are caused by dipole waves, for which the current traces are more spiky . Thus we carried out numerical simulations with a smaller contact doping, $`N_c=2\times 10^{16}`$ cm-3. The frequency of the natural oscillation is now reduced to 4 MHz and the AC intensity is $`2`$ mV. The current at the middle of the first plateau ($`V_{DC}`$=1.5 V) is shown in Fig. 3(a). At finite time, $`I(t)`$ presents dipole-like oscillations and superimposed finite amplitude spikes. The Poincaré map, Fig. 3(c), is much more complicated than in the previous case, showing three well defined distorted loops. Since loops in the Poincaré map are due to combination of strong enough signals of different frequencies , the greater strength of the high-frequency spikes gives rise to the additional loop structure and higher harmonic content (Fig. 3(b)).

As we mentioned above, the double layer structure of Poincaré maps indicates nontrivial attractors for the quasiperiodic case. We have calculated the multifractal dimensions of the attractors, $`D_q`$, for the three cases discussed above; see Fig. 4. In all cases, they correspond to strange attractors  whose $`D_q`$ presents the knee-like structure typical of chaotic attractors with multifractal dimensions . A detailed discussion of the attractor dimensions will be presented elsewhere.

In summary, we have analyzed theoretically the time dependent current through a AC driven SL. We have characterized intriguing Poincaré maps of the quasiperiodic oscillations. The strange attractors which define these Poincaré maps have their physical origin in the complex dynamics of the domain wall. We have shown that distorted loops in the Poincaré maps are related to the presence of high frequency spikes in the current traces. Their frequencies combine with the AC frequency and the low natural frequency to yield a richer power spectrum. Since spikes are associated to extended motion of the DW, we can understand why distorted Poincaré maps are observed at the middle of a I–V plateau and not at its beginning, where the monopole moves over a too small part of the SL. The case of natural self-oscillations due to dipole recycling is different. There the DW of the dipole are generated at the emitter contact and move over the whole SL. In such a case, spikes are more prominent than in the monopole case, and a higher harmonic content and distorted Poincaré maps appear. Our results may help explaining recent experimental evidence showing complex Poincaré maps and intriguing routes to chaos .

Acknowledgments. This work has been supported by the DGES (Spain) grants PB98-0142-C04-01 and PB96-0875, by the European Union TMR contracts ERB FMBX-CT97-0157 and FMRX-CT98-0180 and by the Community of Madrid, project 07N/0026/1998.