text:
# Modelling (001) surfaces of II-VI semiconductors

## Abstract

First, we present a two-dimensional lattice gas model with anisotropic interactions which explains the experimentally observed transition from a dominant $`c(2\times 2)`$ ordering of the CdTe(001) surface to a local $`(2\times 1)`$ arrangement of the Cd atoms as an equilibrium phase transition. Its analysis by means of transfer-matrix and Monte Carlo techniques shows that the small energy difference of the competing reconstructions determines to a large extent the nature of the different phases. Then, this lattice gas is extended to a model of a three-dimensional crystal which qualitatively reproduces many of the characteristic features of CdTe which have been observed during sublimation and atomic layer epitaxy.

Within the last years, potential applications of electronic devices based on II-VI semiconductors have inspired basic research concerning for instance the surface reconstructions of these materials. For CdTe, a fairly complete qualitative phase diagram of $`(001)`$ surfaces has been obtained . An understanding of the interplay of these reconstructions with techniques like molecular beam epitaxy (MBE) or atomic layer epitaxy (ALE) is desirable both for technological applications and for the theory of crystal growth. From the theoretical point of view, the case of CdTe is especially interesting, since it exhibits a phase transition which involves competing reconstructions at typical MBE growth temperatures. While there have been a few models of MBE, especially for GaAs which take into account the effects of surface reconstructions (e.g. ), phase transitions between different reconstructions have not been considered in growth models, so far.

The zinc-blende lattice of CdTe consists of alternating layers of Cd and Te atoms the positions of which lie on regular square lattices. Under vacuum, the CdTe(001) surface is always Cd terminated. The surface is characterized by vacancy structures where the maximum Cd coverage $`\rho _{Cd}`$ is $`1/2`$. At low temperature, one finds a $`c(2\times 2)`$ reconstruction, where the Cd atoms arrange in a checkerboard pattern. Above a critical temperature $`T_c=543K`$ the surface is dominated by a local $`(2\times 1)`$ ordering where the Cd atoms arrange in rows along the $`[110]`$ direction which alternate with rows of vacancies . Electron diffraction experiments have shown that there is a high degree of disorder in the $`(2\times 1)`$ phase . This suggests an interpretation of the reordering as an order-disorder phase transition. An external Cd flux stabilizes the $`c(2\times 2)`$ reconstruction at $`T>T_c`$. If a Te flux is applied, one finds surfaces terminated by Te dimers or trimers .

We first investigate whether the above mentioned transition can be explained within the framework of thermodynamic equilibrium. In a simplifying manner we treat the terminating Cd layer of the surface as a two-dimensional square lattice gas with sites either being occupied by a Cd atom or empty . Chemical bonds and the influence of the underlying crystal are accounted for by effective pairwise anisotropic interactions. An infinite repulsion excludes nearest neighbour (NN) pairs in $`[1\overline{1}0]`$-direction (denoted as $`y`$-direction in the following). This can be justified from considerations like the electron counting rule . In $`[110]`$(x)-direction an attractive interaction $`ϵ_x<0`$ favors the occupation of NN sites. A competing attractive interaction $`ϵ_d<0`$ between diagonal neighbors (NNN) tends to stabilize the $`c(2\times 2)`$ arrangement. The total energy of the system then reads

$`H`$ $`=`$ $`ϵ_dn_{Cd,d}+ϵ_xn_{Cd,x}\mu n_{Cd},`$ (1)
where $`n_{Cd}`$, $`n_{Cd,d}`$ and $`n_{Cd,x}`$ denote the total number of Cd atoms, NNN and NN pairs in x-direction, respectively. The (effective) chemical potential $`\mu `$ controls the mean Cd coverage $`\rho _{Cd}=n_{x,y}`$. By choosing $`ϵ_d=1`$ the energy scale is fixed. In agreement with DFT calculations , $`ϵ_x`$ is chosen such that perfect $`c(2\times 2)`$ and $`(2\times 1)`$ arrangements are nearly degenerate, the groundstate being a $`c(2\times 2)`$ reconstruction.

Besides the coverage $`\rho _{Cd}`$ the system is characterized by the two-point correlations $`C_{Cd}^x=n_{Cd,x}/n_{Cd}`$ and $`C_{Cd}^d=n_{Cd,d}/n_{Cd}`$ which determine the fraction of Cd atoms incorporated in locally $`(2\times 1)`$\- or $`c(2\times 2)`$\- ordered regions.

Through transfer-matrix (TM) calculations and Monte Carlo (MC) simulations at constant $`\rho _{Cd}`$ we obtain an estimate of the phase boundaries in the temperature-coverage plane (see fig. 1). At low temperatures (III), a $`c(2\times 2)`$ ordered phase with $`\rho _{Cd}1/2`$ coexists with a disordered phase of low coverage. At higher temperatures, the system becomes homogeneously disordered (II) or ordered (I) depending on the coverage. For a range of coverages $`\rho _{Cd}`$ the transition from (III) to (II) is accompanied by a sudden decrease of $`C_{Cd}^d`$ and a simultaneous increase of $`C_{Cd}^x`$, meaning that for these values of $`\rho _{Cd}`$ the $`(2\times 1)`$-structure is locally prevailing in the disordered regime (the area left of the dashed line (or squares) in fig. 1). With decreasing $`|ϵ_x|`$ the tricritical point, i.e. the point in the $`\rho _{Cd}`$-$`T`$-diagram where (I),(II) and (III) meet, shifts to smaller coverage and higher temperature.

Even more so does the line which separates $`c(2\times 2)`$ from $`(2\times 1)`$ prevalence. This might explain, why a transition from $`c(2\times 2)`$ to $`(2\times 1)`$ prevalence has not been found in ZnSe, so far. There, the energy difference between perfect $`c(2\times 2)`$ and $`(2\times 1)`$ reconstructions is indeed greater than for CdTe .

In order to investigate the evolution of the morphology of the surface in non-equilibrium processes like growth and sublimation, we extend the planar lattice gas to a model of a three-dimensional crystal. In a previous publication, this has been done for a simple cubic crystal which consists of alternating layers of the two atomic species . In this paper, we present the next step towards a realistic modelling of CdTe by simulating the correct zinc-blende lattice structure. Besides the Cd-Te bond (binding energy $`ϵ_c`$) we consider interactions between atoms on NN and NNN sites in one layer of the crystal. The basic idea here is to model the interactions of Cd atoms on the surface with the anisotropic interactions $`ϵ_x`$, $`ϵ_d`$ and a hardcore repulsion in the $`y`$-direction while there is an isotropic attraction between Te atoms ($`ϵ_{Te,b}`$) and Cd atoms inside the bulk of the crystal ($`ϵ_{Cd,b}`$).

Due to the topology of the zinc-blende lattice the condition of ergodicity can not be fulfilled within the limits of a solid-on-solid (SOS) model, where each atom must be bound to at least two atoms of the opposite species in the layer below. This can be achieved by considering the formation of Te trimers which have been observed experimentally. Additionally, we consider Te atoms which have only one bond to a neighbouring Cd atom. In the following, Te atoms in these configurations will be denoted as Te atoms. Te atoms neutralize the hardcore repulsion between Cd atoms. Otherwise, the deposition of Te on a Cd terminated surface with a $`c(2\times 2)`$ reconstruction would be impossible.

As a simplification, we assume that point defects are not incorporated into the crystal. Then, the bulk can be described in a SOS manner, while the Te atoms are considered seperately. In summary, the Hamiltonian of the model is

$`H`$ $`=`$ $`ϵ_cn_c+ϵ_{Cd,b}n_{Cd,b}+ϵ_{Te,b}+ϵ_hn_{Cd,h}`$ (2)
$`+`$ $`ϵ_xn_{Cd,x}+ϵ_dn_{Cd,d}+ϵ^{}n^{}`$
which is a function of the numbers of Cd-Te bonds ($`n_c`$), Te atoms ($`n^{}`$), NN and NNN pairs of surface Cd atoms ($`n_{Cd,x}`$, $`n_{Cd,d}`$), NN pairs of Cd atoms in the bulk and Te atoms ($`n_{Cd,b}`$, $`n_{Te,b}`$) and surface Cd atoms next to incorporated Cd atoms ($`n_{Cd,h}`$). Pairs of surface Cd atoms which are not bound to Te atoms on neighbour sites in the $`y`$-direction are forbidden.

Besides the deposition of adatoms we consider desorption, diffusion to NN and NNN sites and diffusion across steps. These are Arrhenius-activated processes with rates $`\nu \mathrm{exp}[E/(kT)]`$. Desorption requires an activation energy $`E=\mathrm{\Delta }H`$. We have chosen Metropolis-like energy barriers $`E=\mathrm{max}\{B,B+\mathrm{\Delta }H\}`$ for diffusion processes. There is a barrier $`B=B_0`$ for the diffusion of atoms with two bonds and a barrier $`B^{}`$ for the diffusion of Te atoms.

Unfortunately, the available experimental data are not sufficient for a systematic fit of the model parameters and there are no values of diffusion barriers available from ab initio or semi-empirical calculations. Therefore, the choice of the parameter set is partly based on estimates and physical reasoning. In the following, we present results that have been obtained with the parameter set $`ϵ_d=0.086eV`$, $`ϵ_x=0.168eV`$, $`ϵ_{Cd,b}=ϵ_{Te,b}=0.069eV`$, $`ϵ_h=0`$, $`ϵ_c=0.43eV`$, $`ϵ^{}=0.26eV`$, $`B_0=0.78eV`$, $`B^{}=0.17eV`$ and $`\nu =10^{12}/s`$.

Experiments have revealed a high density of steps on CdTe surfaces with a preferential orientation in the $`(100)`$ direction. This leads to a dominant contribution of step flow sublimation. Therefore, we have simulated sublimation on vicinal surfaces with steps parallel to the $`[100]`$ direction at a distance of $`16`$ lattice constants.

Our model reproduces the observation of Cd terminated surfaces in vacuum. A Te terminated surface decays in a temperature-dependent time intervall $`\tau (T)`$. The corresponding desorption probability $`p(T)1/\tau (T)`$ follows an Arrhenius law with an activation energy of $`1.25\pm 0.02eV`$. After this onset, we obtain a stationary state where Cd and Te evaporate stoichiometrically. Then, the temperature dependence of the sublimation rate follows an Arrhenius law with a higher activation energy of $`1.437\pm 0.003eV`$. Experimentally, activation energies of $`1.54eV`$  and $`0.9eV`$  have been found for stoichiometric step flow sublimation and the decay of a Te terminated surface.

Figure 2 shows $`\rho _{Cd}`$ and the correlations $`C_{Cd}^x`$, $`C_{Cd}^d`$ in the stationary state as functions of $`T`$. In the investigated temperature interval, $`\rho _{Cd}`$ decreases only slightly from $`\rho _{Cd}=0.49`$ at $`T=450K`$ to $`\rho _{Cd}=0.43`$ at $`T=700K`$. At low temperature, the surface is dominated by a $`c(2\times 2)`$ arrangement of the Cd atoms ($`C_{Cd}^d>C_{Cd}^x`$). On the contrary, at high temperature the Cd atoms tend to arrange in the rows of the $`(2\times 1)`$ reconstruction ($`C_{Cd}^d<C_{Cd}^x`$). The crossover between both regimes is at $`T_c=570K`$. This is the counterpart of the phase transition observed in the 2D lattice gas model. However, an investigation of $`C_{Cd}^x`$, $`C_{Cd}^d`$ in the 2D lattice gas at the $`\rho _{Cd}`$ measured in the 3D model shows, that the non-equilibrium conditions of sublimation increase the dominance of $`(2\times 1)`$ over $`c(2\times 2)`$ in the high temperature regime.

ALE provides a scenario where our model can be investigated in a situation present in technical applications. The basic idea is to obtain self-regulated layer-by-layer growth by alternate deposition of pure Cd and Te. Figure 3 shows the evolution of the thickness of the deposited layer and the correlations $`C_{Cd}^x`$, $`C_{Cd}^d`$ at $`T=500K<T_c`$. The initial state was a flat Te terminated surface. In the first phase, a Cd flux of $`5ML/s`$ is applied, in the second phase Te at the same rate. At the end of the first phase, the surface is Cd terminated with a $`c(2\times 2)`$ reconstruction. In the Te phase, Te terminated islands form which, at the end of the cycle, cover $`50\%`$ of the surface. This is due to the fact that only half of the Cd atoms needed for a closed monolayer are present. Note, that the decrease of $`\rho _{Cd}`$ at the beginning of the Te cycle leads to an increase of $`C_{Cd}^x`$, which is consistent with the phase diagram (fig. 1) of the 2D lattice gas. In the following Cd phase, the islands are covered with Cd. Surprisingly, now the reconstruction is $`(2\times 1)`$. This can be understood from the fact that the lattice gas interactions are present only between particles in the same layer. Thus, the island edges impose open boundary conditions to the lattice gas of Cd atoms on the island. In constrast to a $`c(2\times 2)`$ reconstructed domain, a $`(2\times 1)`$ terminated island can reduce the energy of its boundary by elongating in $`x`$-direction. Since the ground state energies of both structures are nearly degenerate, the formation of $`(2\times 1)`$ may thus reduce the surface free energy. In the following Te phase, Cd and Te atoms diffuse into the gaps between the islands, which leads to the formation of a closed, Te covered monolayer again. These observations closely resemble ideas suggested in  to explain experimental observations.

In summary, we have presented a model of CdTe(001) surfaces which reproduces several experimental observations in a semi-quantitative way. In the future, we will extend this model by considering the dimerization of surface Te atoms and the fact that Te is adsorbed and desorbed as Te2. In cooperation with an experimental group, we will try a more quantitative determination of an appropriate parameter set.


format: text

text:
# References

IMSc/2001/07/39

hep-th/0107255

Some Consequences of the Generalised Uncertainty

Principle: Statistical Mechanical, Cosmological,

and Varying Speed of Light

S. Kalyana Rama

Institute of Mathematical Sciences, C. I. T. Campus,

Taramani, CHENNAI 600 113, India.

email: krama@imsc.ernet.in

ABSTRACT

> We study the dynamical consequences of Maggiore’s unique generalised uncertainty principle (GUP). We find that it leads naturally, and generically, to novel consequences. In the high temperature limit, there is a drastic reduction in the degrees of freedom, of the type found, for example, in strings far above the Hagedorn temperature. In view of this, the present GUP may perhaps be taken as the new version of the Heisenberg uncertainty principle, conjectured by Atick and Witten to be responsible for such reduction. Also, the present GUP leads naturally to varying speed of light and modified dispersion relations. They are likely to have novel implications for cosmology and black hole physics, a few of which we discuss qualitatively.

PACS numbers: 11.25.-w, 05.90.+m, 98.80.Cq, 04.70.-s

1. Based on gedanken experiments in string theory  and black holes , the Heisenberg uncertainty principle is found to be modified to

$$\mathrm{\Delta }x\frac{\mathrm{}}{\mathrm{\Delta }p}+(const)\frac{\lambda ^2\mathrm{\Delta }p}{\mathrm{}}$$
(1)
where $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }p`$ denote position and momentum uncertainties, and $`\lambda `$ is a length parameter, given by string length and Planck length in the above contexts. Under certain assumptions, Maggiore has derived in  unique generalised commutation relations (GCRs) which lead to a generalised uncertainty principle (GUP) which, in turn, leads to (1) in a suitable limit. Under a different set of assumptions it is possible to obtain more general commutation relations, e.g. as in , which also lead to (1). However, they are not unique. Hence, in the following, we will consider Maggiore’s generalisation only, although our analysis is applicable to other cases also.

The GCRs are kinematical. The dynamics is determined by specifying a Hamiltonian $`H`$. To illustrate explicitly the dynamical consequences of the GUP, we choose two candidates for $`H`$ and study the statistical mechanics and the particle dynamics of free particle systems obeying the GUP 111 In the case of the GCRs given by , some aspects of the particle dynamics have been studied in . In the case of the GCRs given by , some aspects of the particle dynamics, statistical mechanics, and black hole physics have been studied in , , and  respectively. See also . . We find that the GUP leads naturally to many novel consequences.

In the high temperature limit, we find that there is a drastic reduction in the degrees of freedom (d.o.f) and the corresponding free energy is analogous to that found in certain topological field theories and in strings far above the Hagedorn temperature . In this context, Atick and Witten had indeed conjectured in  that a new version of the Heisenberg uncertainty principle may be responsible for such a reduction in the d.o.f. Here, we see that such a reduction emerges naturally as a consequence of the GUP. In view of this, the present GUP may perhaps be taken as the conjectured new version of the Heisenberg uncertainty principle.

Another consequence of the GUP is the natural emergence of the varying speed of light (VSL) and the modified dispersion relations which, in turn, have non trivial implications for cosmology and black hole physics 222VSL theories were first postulated in  and studied further in . Their implications have been studied in . The implications of modfied dispersion relations for cosmology have been studied in .. For one of the $`H`$’s considered here, the VSL and the free energy together is likely to solve the horizon problem in cosmology, as in . The corresponding VSL is likely to have novel implications for black hole physics also.

The plan of the paper is as follows. In section 2, we present the details of Maggiore’s GCR and GUP, and the two candidate Hamiltonians. In section 3, we study the statistical mechanics and discuss the consequences. In section 4, we study the particle dynamics and discuss the consequences. In section 5, we present a brief summary and close by mentioning a few issues for further study.

2. The uncertainty principle (1) can be thought of as arising from a generalisation of the commutation relations between the position operators $`X_i`$ and the momentum operators $`P_j`$ in $`d`$-dimensional space where $`i,j=1,2,\mathrm{},d`$. Seeking the most general deformed Heisenberg algebra in $`d=3`$, Maggiore has derived in  the generalised commutation relations (GCRs) between $`X_i`$ and $`P_j`$, determined uniquely under the following assumptions: (i) The spatial rotation group and, hence, the commutators $`[J_i,J_j]`$, $`[J_i,X_j]`$, and $`[J_i,P_j]`$ are undeformed. (ii) The translation group and, hence, the commutators $`[P_i,P_j]`$ are undeformed. (iii) The commutators $`[X_i,X_j]`$ and $`[X_i,P_j]`$ depend on a deformation parameter $`\lambda `$, with dimension of length, and reduce to the undeformed ones in the limit $`\lambda 0`$. The GCRs that follow uniquely from these assumptions are given by

$`[X_i,X_j]`$ $`=`$ $`iϵ\mathrm{}^2\lambda ^2ϵ_{ijk}J_k,J_k=iϵ_{klm}P_l{\displaystyle \frac{}{P_m}}`$ (2)
$`[X_i,P_j]`$ $`=`$ $`i\mathrm{}\delta _{ij}f,f=\sqrt{1+{\displaystyle \frac{ϵ\lambda ^2}{\mathrm{}^2}}(P^2+m^2c^2)}`$ (3)
where $`ϵ=\pm 1`$, $`P^2=_iP_i^2`$, and $`\lambda `$ is a length parameter. The GCR (3) then leads to the generalised uncertainty principle (GUP)

$$\mathrm{\Delta }x_i\mathrm{\Delta }p_j\frac{\mathrm{}}{2}\delta _{ij}f.$$
(4)
In the limit $`\lambda ^2(p^2+m^2c^2)\mathrm{}^2`$ and $`\lambda \mathrm{\Delta }p\stackrel{<}{_{}}\mathrm{}`$, where $`p^2`$ is the eigenvalue of $`P^2`$, equation (4) reduces to (1). See  for details.

In the following, we consider $`d`$-dimensional space 333 With the right hand side of (2) written explicitly in terms of $`P_i`$ and $`\frac{}{P_j}`$, equations (2) - (4) are valid in $`d`$-dimensional space also, as can be easily verified.. We set $`\mathrm{}=c=1`$ unless indicated otherwise. The case $`ϵ=1`$ implies a bound $`\lambda ^2(p^2+m^2)<1`$, whose physical significance is not clear. Hence, we will consider the case $`ϵ=1`$ only. Also, we consider non rotating systems only and, therefore, set $`J_k=0`$ in (2). The only nontrivial GCR is then given by (3) with

$$f=\sqrt{1+\lambda ^2(P^2+m^2)},P^2=\underset{i=1}{\overset{d}{}}P_i^2.$$
(5)
We study the consequences of the GCR (3), equivalently of the GUP (4), with $`f`$ given by (5). It is important to note that the commutation relations are kinematical only. The dynamics is determined by the Hamiltonian $`H`$ which, therefore, needs to be specified. In the following, we assume that $`H`$ depends on $`P`$ only through the rotationally invariant combination $`\sqrt{P^2+m^2}`$. Furthermore, for the sake of simplicity, we consider here only the free particle case, for which $`H`$ is independent of $`X`$.

To illustrate the non trivial consequences of the GUP (4), we consider two choices for $`H`$. They are given by

$`H^{}=1`$ $``$ $`H=\sqrt{P^2+m^2}`$ (6)
$`fH^{}=1`$ $``$ $`Sinh\lambda H=\lambda \sqrt{P^2+m^2}`$ (7)
where $`H^{}`$ is the derivative of $`H`$ with respect to $`\sqrt{P^2+m^2}`$ and $`f`$ is given by (5). The Hamiltonian $`H`$ in (6) is, perhaps, the simplest and a natural choice; $`H`$ in (7) is obtained in  from the first Casimir operator and, hence, is also a natural choice from group theoretic point of view.

3. Consider the statistical mechanics of a system of free particles confined in a $`d`$-dimensional volume $`V`$, which obey the GUP (4) with $`f`$ given by (5). The calculations in the microcanonical or canonical ensemble approach are complicated. But they are simple in the grand canonical ensemble approach which we, therefore, use.

In the standard case where $`\lambda =0`$, the one-particle phase space measure is given by $`h^d`$ where $`h`$ is the Planck’s constant. Physically, this is because the phase space is divided into cells of volume $`h^d`$ as a consequence of the Heisenberg uncertainty principle. Various phase space integrals will be of the form $`d^dxd^dph^d()`$ where $`x`$ and $`p`$ denote the eigenvalues of $`X`$ and $`P`$. In the present case, the particles are assumed to obey the GUP (4). Consequently, the phase space must be divided into cells of volume $`h^df^d`$. Then, the one-particle phase space measure is given by $`h^df^d`$, and various phase space integrals are of the form $`d^dxd^dph^df^d()`$.

The function $`f`$ and, in the free particle case, the Hamiltonian $`H`$ depend on $`P`$ only. $`x`$-integration then simply gives a volume factor $`V`$. Writing $`p`$ in terms of energy $`E`$, which is the eigenvalue of $`H`$, the phase space integrals can be written as

$$\frac{d^dxd^dp}{h^df^d}()𝑑Eg(E)().$$
(8)
The measure $`g(E)`$ is the analog of the one-particle density of states. It can be easily calculated and is given by

$$g(E)=\frac{𝒞p^{d2}}{f^dE^{}}\sqrt{p^2+m^2},\mathrm{where}𝒞\frac{\mathrm{\Omega }_{d1}V}{h^d},$$
(9)
$`\mathrm{\Omega }_{d1}`$ is the area of a unit $`(d1)`$-dimensional sphere, and $`E^{}`$ is the derivative of $`E`$ with respect to $`\sqrt{p^2+m^2}`$. In equation (9), $`p`$ is to be expressed in terms of $`E`$.

Consider the grand canonical ensemble. It can be easily verified that the definitions of, and the relations between, various thermodynamical quantities all remain unchanged, with $`g(E)`$ given as in (9) . Thus, we have

$$\beta F=\beta 𝒫V=ln𝒵=\frac{1}{a}_0^{\mathrm{}}𝑑Eg(E)ln(1+ae^{\beta (E\mu )})$$
(10)
where we have used the standard notation: $`\beta =T^1`$ is the inverse temperature, $`F`$ is the free energy, $`𝒫`$ is the pressure, $`𝒵`$ is the grand canonical partition function, and $`\mu `$ is the chemical potential. Also, $`a=1,0`$, or $`+1`$ depending on whether the particles obey, respectively, Bose-Einstein, Maxwell-Boltzmann, or Fermi-Dirac statistics. When $`a=0`$, $`ln𝒵`$ is to be evaluated in the limit $`a0`$. Various thermodynamical quantities can then be calculated using (10): for example, the internal energy $`U=\frac{ln𝒵}{\beta }`$, the particle number $`N=\frac{ln𝒵}{\beta \mu }`$, the entropy $`S=\beta (U+𝒫V\mu N)`$, etc.

It is clear from the above formulae, or from physical arguments, that the effect of $`\lambda `$ will be considerable only when the temperature/energy is of $`𝒪(\lambda ^1)`$ or higher. Therefore, the limit of interest here is the high temperature limit $`\beta \lambda `$. Also, we expect that $`\lambda `$ is extremely small, in particular, $`\lambda m1`$. (For example, $`\lambda string/Planck`$ $`length`$.) Furthermore, in the high temperature limit, the statistics is irrelevent. Therefore, for the sake of simplicity, we set $`m=0`$ and $`a=0`$ in the following. Then $`\mu =0`$ since $`m=0`$. One then obtains $`U,S,N`$ etc. in terms of $`V`$ and $`T`$ using (10) .

To proceed further, $`p`$ and $`f`$ in equation (9) are to be expressed in terms of $`E`$, for which an explicit form of $`H`$ is required. We consider $`H`$ given in (6) and (7), with $`m=0`$. A simple algebra then shows that $`g(E)`$ is given by

$`H^{}=1`$ $``$ $`g(E)={\displaystyle \frac{𝒞E^{d1}}{(1+\lambda ^2E^2)^{\frac{d}{2}}}}`$ (11)
$`fH^{}=1`$ $``$ $`g(E)={\displaystyle \frac{𝒞}{\lambda ^{d1}}}(tanh\lambda E)^{d1}.`$ (12)
The partition function $`𝒵`$ and other quantities can now be evaluated in closed form in terms of special functions, as described in the Appendix.

We expect to obtain the standard results in the limit $`\beta \lambda `$ and to obtain the non trivial features, if any, in the limit $`\beta \lambda `$444In both of these limits, it is easier to evaluate the integral in (10) directly, using equations (11) and (12) and appropriate Taylor expansions, than to work with the special functions.. A simple calculation in the limit $`\beta \lambda `$ shows that, to the leading order in $`\frac{\lambda }{\beta }`$, the thermodynamical quantities are independent of $`\lambda `$, and are indeed the standard ones for a gas of massless free particles in $`d`$-dimensional space obeying the Heisenberg uncertainty principle. They are given by

$$\beta F=\frac{𝒞(d1)!}{\beta ^d},U=\frac{𝒞d!}{\beta ^{d+1}},S=\frac{𝒞(d+1)!}{d\beta ^d}.$$
(13)
Note that the free energy has the behaviour $`|\beta F|T^d`$.

Consider the limit $`\beta \lambda `$. The results now depend on whether $`g(E)`$ is given by equation (11) or (12). In the case where $`g(E)`$ is given by (11), the thermodynamical quantities are given, to the leading order in $`\frac{\beta }{\lambda }`$, by

$$\beta F=const+\frac{𝒞}{\lambda ^d}ln(\lambda T),U=\frac{𝒞}{\beta \lambda ^d},S=const+\frac{𝒞}{\lambda ^d}(1+ln(\lambda T)).$$
(14)
Note that in the high temperature limit $`\beta \lambda `$, the free energy has the behaviour $`|\beta F|(const)+\mathrm{ln}T`$ for any value of $`d`$. This indicates a drastic reduction in the degrees of freedom (d.o.f). Such reduction may be possible, in the context of Heisenberg uncertainty principle, if a continuum field theory is replaced by a lattice theory, with a finite number of bose oscillators at each site ; or, in certain topological theories  with general covariance restored at short distances . Here, we see that such a behaviour emerges naturally as a consequence of the GUP (4), for systems whose Hamiltonian $`H`$ is given by (6).

The thermodynamical relations (14), which, upto polarisation factors, are valid for photons also, may have interesting cosmological consequences. In a recent paper , the authors analyse the thermodynamical behaviour of photons in the framework of non commutative geometry, postulating a model dependent varying speed of light (VSL). Amazingly, although their set up bears no discernible relation to the present one, the equation of state

$$U=\frac{𝒞}{\beta \lambda ^d},\frac{𝒫}{U}=(const)+ln(\lambda T)$$
(15)
in the present case, obtained from (10) and (14), is essentially the same as that obtained in  (equations (38) and (39) of ). It is shown in  that such an equation of state, with an additional ingredient to be mentioned in the next section, solves the horizon problem.

Consider the limit $`\beta \lambda `$, now in the case where $`g(E)`$ is given by (12). The thermodynamical quantities are given, to the leading order in $`\frac{\beta }{\lambda }`$, by

$$\beta F=const+\frac{𝒞}{\lambda ^{d1}\beta },U=\frac{𝒞}{\lambda ^{d1}\beta ^2},S=const+\frac{2𝒞}{\lambda ^{d1}\beta }.$$
(16)
Note that in the high temperature limit $`\beta \lambda `$, the free energy has the behaviour $`|\beta F|(const)+T`$ for any value of $`d`$. This indicates a drastic reduction in the d.o.f. Precisely such a free energy, and hence such a reduction, has been found in  in the case of the strings at temperatures far above the Hagedorn temperature. In this context, Atick and Witten had indeed conjectured in  that a new version of the Heisenberg uncertainty principle may be responsible for such a reduction in the d.o.f. Here, we see that such a drastic reduction in the d.o.f emerges naturally as a consequence of the GUP (4), for systems whose Hamiltonian $`H`$ is given by (7). In view of this, the GUP (4) may perhaps be taken as the conjectured new version of the Heisenberg uncertainty principle.

The physical origin of the reduction in the d.o.f is easy to understand. When the system obeys the GUP (4), the volume of the phase space cells is $`h^df^d`$. It grows at high temperatures/energies and, since the function $`f`$ is given by (5), $`p^d`$. Consequently, the number of available cells is enormously reduced, compared with the standard case. This, essentially, is the origin of the reduction in the d.o.f seen above. The precise amount of reduction depends on the choice of the Hamiltonian $`H`$. For the $`H`$ given by $`f^\alpha H^{}=1`$, where $`\alpha =0,1`$,555For $`\alpha >1`$, $`H`$ becomes bounded, i.e. its eigenvalue $`E(const)`$ as $`p\mathrm{}`$, a behaviour whose significance is not clear. it can be seen from equation (9) that the reduction is such that the resulting d.o.f are equal to that of an $`\alpha `$-dimensional system obeying the Heisenberg uncertainty principle and with an effective volume (effective number of sites in the $`\alpha =0`$ case) $`V_\alpha V\lambda ^{\alpha d}`$. This is precisely the result seen explicitly in equation (14) for $`\alpha =0`$ and in equation (16) for $`\alpha =1`$.

4. Consider the dynamics of particles which obey the GUP (4), equivalently the GCR (3), with $`f`$ given by (5). We first note in passing that the time energy uncertainty relation remains unchanged. The derivation proceeds in the standard way, e.g. as in , with the result that

$$\mathrm{\Delta }t_Q\mathrm{\Delta }E\frac{\mathrm{}}{2}$$
(17)
where $`\mathrm{\Delta }t_Q\mathrm{\Delta }Q\left(\frac{dQ}{dt}\right)^1`$ is the time uncertainty, time being measured by measuring the variation of an observable $`Q`$, which has no explicit time dependence and, hence, obeys $`\frac{dQ}{dt}=\frac{i}{\mathrm{}}[H,Q]`$ where $`H`$ is the Hamiltonian.

We now define the velocity operator $`V_i`$ by (see  also)

$$V_i\frac{dX_i}{dt}=\frac{i}{\mathrm{}}[H,X_i].$$
(18)
In the case of a non rotating system for which $`[X_i,X_j]=0`$, or a free particle system for which $`H`$ is independent of $`X_i`$, one obtains using (3) that

$$V_i=\frac{fH^{}P_i}{\sqrt{P^2+m^2}}$$
(19)
where $`H^{}`$ is the derivative of $`H`$ with respect to $`\sqrt{P^2+m^2}`$. Denoting the eigenvalue of $`V_i`$ by $`v_i`$, the speed $`v`$ of a particle with mass $`m`$ is given by

$$v=\left(\underset{i=1}{\overset{d}{}}v_i^2\right)^{\frac{1}{2}}=\frac{pfE^{}}{\sqrt{p^2+m^2}}.$$
(20)
In Quantum Mechanics, with $`\mathrm{}=1`$, the energy eigenfunction will have a time dependence $`e^{i\omega t}`$ where $`\omega =E`$, the eigenvalue of $`H`$. Its group velocity $`\frac{d\omega }{dk}`$ can then be identified naturally with the velocity $`v`$ in (20) where, now, $`E=\omega `$. This leads to a modified dispersion relation, given by

$$\frac{d\omega }{dk}=\frac{pfE^{}}{\sqrt{p^2+m^2}}$$
(21)
where one sets $`E=\omega `$ on the right hand side. Note that equation (21) can be thought of as defining the wave number $`k`$ in the position space. Indeed, in $`d=1`$, $`k`$ is the eigenvalue of the operator $`K`$ defined by $`f\frac{dK}{dP}=1`$, so that $`[X,K]=i`$. Equation (21) then follows since $`\frac{dH}{dK}=\frac{PH^{}}{\sqrt{P^2+m^2}}\frac{dP}{dK}=\frac{PfH^{}}{\sqrt{P^2+m^2}}`$.

Furthermore, the speed of light, denoted by $`C`$, can be identified naturally with the speed of a particle with mass $`m=0`$. Equation (20) then gives

$$C=fE^{}\mathrm{and}vC.$$
(22)
In equations (20), (21), and (22) $`p`$ and $`f`$ are to be expressed in terms of $`E`$, for which an explicit form of $`H`$ is required. For $`H`$ given by (6) we have, after a simple integration,

$`v`$ $`=`$ $`{\displaystyle \frac{\sqrt{(E^2m^2)(1+\lambda ^2E^2)}}{E}},C=\sqrt{1+\lambda ^2E^2}`$ (23)
$`\lambda ^2\omega ^2`$ $`=`$ $`(1+\lambda ^2m^2)Sinh^2\lambda k+\lambda ^2m^2.`$ (24)
For $`H`$ given by (7) (see  also) we have, after a simple integration,

$`v`$ $`=`$ $`{\displaystyle \frac{\sqrt{Sinh^2\lambda E\lambda ^2m^2}}{Sinh\lambda E}},C=1`$ (25)
$`Cosh\lambda \omega `$ $`=`$ $`\sqrt{1+\lambda ^2m^2}Cosh\lambda k.`$ (26)
We now discuss the physical significance of these results. Consider first the modified dispersion relation (21). For $`\lambda =0`$, it reduces to $`\omega ^2=k^2+m^2`$. For $`\lambda 0`$, the modification is generically non trivial. The exception is when the Hamiltonian $`H`$ is given by (7), in which case the modification is only marginal. These can be seen explicitly by considering the high energy ‘transplanckian’ limit $`\lambda \omega 1`$ of equations (24) and (26).

Recent studies  have suggested that such ‘transplanckian’ modifications of dispersion relation may have observable consequences for the density perturbations that arise during inflation. In these studies, the modified dispersion relations need to be postulated. Here, however, we see that generically the GUP (4) leads naturally to modified dispersion relations. It is clearly of interest to study their consequences for the density fluctuations that arise during inflation.

Consider the speed of light $`C`$. It is clear from equation (22) that, generically, $`C`$ is varying and is a non trivial function of energy $`E`$. The exception is when the Hamiltonian $`H`$ is given by (7), in which case $`C=1`$. Such ‘varying speed of light’ (VSL) theories have been extensively studied , and found to have non trivial implications for cosmology  and black hole physics . In these theories, VSL needs to be postulated. Here, however, we see that generically the GUP (4) leads naturally to VSL. It is clearly of interest to study its consequences, which are likely to be non trivial.

For example, $`C`$ given by (23) increases with energy $`E`$. This is precisely the ingredient, alluded to below equation (15), that is necessary, but postulated, in  to solve the horizon problem. Here, it arises naturally. Moreover, a preliminary analysis shows that the VSL, given by (23), and the photon distribution at energies $`E\lambda ^1`$, deriveable from (10) and (11), both have the right behaviour needed to solve the horizon problem, as in , but within the present framework. A detailed analysis, however, is beyond the scope of the present letter.

Equations (23) are likely to have novel implications for black hole physics also. The horizon of a black hole can naively be thought of as the place where the escape velocity $`=1`$, in units where $`c=1`$. Particles can then escape from, or from even inside, the horizon if their energy $`E`$ is sufficiently high since their speed can then be $`>1`$. This may, therefore, provide a mechanism for the transfer of information from inside the horizon to the outside, a process for which no mechanism is known at present .

Perhaps more correctly, the horizon is to be thought of as the place from where nothing can escape. Then the escape velocity at the horizon must be infinite since $`v`$ and $`C\mathrm{}`$ as $`E\mathrm{}`$. Very likely, therefore, the horizon size must be infinitesimally small or, perhaps equivalently, no black holes can form 666In this context, note that a certain class of VSL theories are shown in  to be equivalent to generalised Brans-Dicke theories. For a certain class of the later theories, it is argued in  that black holes are unlikely to form. For a discussion of black hole physics in VSL theories from another point of view, see .. These implications for the black hole physics of the GUP (4), with $`f`$ given by (5) and the Hamiltonian $`H`$ given by (6), are very interesting but the present arguments are, admittedly, qualitative. Unfortunately, in the absence of a Lorentz and/or general coordinate invariant formulation of the GUP, these issues cannot be addressed rigorously.

5. In summary, we have studied the physical consequences of the GUP (4) that follows from the GCR (3) which is determined uniquely by Maggiore under a set of assumptions. We studied the statistical mechanics and the particle dynamics of systems obeying the GUP (4) and found novel consequences arising in a natural way. For example, the GUP leads naturally to free energies of the form found in certain topological field theories and in strings far above the Hagedorn temperature. It also leads naturally to VSL and to modified dispersion relations. Among other things, these features are likely to solve the horizon problem in cosmology, and may provide novel insights into black hole physics also.

There are numerous issues that require further study. We close by mentioning a few of them. (i) Understanding the physical significance of the bound $`\lambda ^2(p^2+m^2)<1`$ in the $`ϵ=1`$ case. (ii) Finding the physical principle, if any, which selects a given Hamiltonian $`H`$, e.g. the one given by $`f^\alpha H^{}=1`$ for a given value of $`\alpha `$. (iii) Finding a Lorentz and/or general coordinate invariant formulation of the GUP which is crucial, for example, for the study of black hole physics. (iv) Exploring relations, if any, between GUP and string theory, topological field theory, and VSL theories. That such a relation may exist is, perhaps, indicated by equations (16), (14), and (22).

Appendix

The partition function in (10), with $`m=a=0`$, can be obtained in a closed form in terms of special functions both for the cases where $`H^{}=1`$, and where $`fH^{}=1`$. We need to evaluate integrals of the form

$`I_{m,n}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑tt^n(\alpha ^2+t^2)^{\frac{m}{2}}e^{st}\mathrm{for}H^{}=1`$ (27)
$`J_{m,n}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑tt^n(tanht)^me^{st}\mathrm{for}fH^{}=1,`$ (28)
where $`t=\lambda E`$, $`s=\frac{\beta }{\lambda }`$, and $`\alpha ^2=1`$. It is easy to see that

$`I_{2k,n}`$ $`=`$ $`{\displaystyle \frac{(1)^{n+k1}}{(k1)!}}\left({\displaystyle \frac{d}{ds}}\right)^n\left({\displaystyle \frac{d}{d\alpha ^2}}\right)^{k1}I_{2,0}`$
$`I_{2k+1,n}`$ $`=`$ $`{\displaystyle \frac{(1)^{n+k}\mathrm{\hspace{0.33em}2}^{k+1}}{(2k1)!!}}\left({\displaystyle \frac{d}{ds}}\right)^n\left({\displaystyle \frac{d}{d\alpha ^2}}\right)^{k+1}I_{1,0}`$
$`J_{m,n}`$ $`=`$ $`(1)^n\left({\displaystyle \frac{d}{ds}}\right)^nJ_{m,0}.`$
Clearly, $`J_{0,0}=\frac{1}{s}`$. Moreover,

$`I_{2,0}`$ $`=`$ $`{\displaystyle \frac{1}{\alpha }}(ci(\alpha s)sin(\alpha s)si(\alpha s)cos(\alpha s))`$
$`I_{1,0}`$ $`=`$ $`{\displaystyle \frac{\pi \alpha }{2s}}(𝐇_1(\alpha s)𝐘_1(\alpha s))`$
$`J_{1,0}`$ $`=`$ $`\beta ({\displaystyle \frac{s}{2}}){\displaystyle \frac{1}{s}}`$
where $`I_{2,0}`$ is given in equation (3.354.1) of ; $`J_{1,0}`$ in (3.541.7) of ; and $`I_{1,0}`$ is obtained using equation (4.2.27) of  and the properties of the Laplace transform of the derivative of a function. Here, $`ci`$ and $`si`$ are the cosine and sine integrals respectively (see section (8.23) of ), $`𝐇_1`$ is the Struve’s function, $`𝐘_1`$ is the Bessel function of the second kind, and the function $`\beta (x)`$ is related to the derivatives of the gamma function (see section 8.37 of ).

We now obtain a recursion relation for $`J_{m,0}`$. We have, for $`m1`$,

$$T_m(x)_0^x𝑑x^{}(tanhx^{})^m=\frac{(tanhx)^{m1}}{m1}+T_{m2}(x).$$
(29)
Thus, after a partial integration in $`J_{m,0}`$ and using $`T_m(0)=0`$, one obtains

$$J_{m,0}=s_0^{\mathrm{}}𝑑xT_m(x)e^{sx}.$$
(30)
Equations (28), (29), and (30) now lead to the recursion relation

$$J_{m,0}=\frac{s}{m1}J_{m1,0}+J_{m2,0},$$
(31)
using which $`J_{m,0}`$ can be obtained for $`m2`$. Using these formulae, the partition function and other quantities can be evaluated in closed form.


text:
# VLBA Observations and Effelsberg Monitoring of a Flaring Water Maser in Mrk 348

## 1. Introduction

Mrk 348 (NGC 262) is a Seyfert 2 galaxy at a redshift of 0.01503 (Huchra et al. 1999). The galaxy is classified as an S0 with a low inclination, and exhibits a large Hi halo which may have been produced by an interaction with the companion galaxy NGC 266 (e.g. Simkin et al. 1987). VLBA images (Ulvestad et al. 1999) reveal a small-scale double continuum source, the axis of which is aligned with the optical (\[OIII\], Capetti et al. 1996) emission. Astrometry indicates that the optical and VLBI cores are nearly coincident. Ground-based observations (Simpson et al. 1996) show evidence of a dust lane crossing the nucleus and an ionization cone. Attempts to detect the expected obscuring torus at radio wavelengths (e.g. Hi, Gallimore et al 1999; CO, Taniguchi et al. 1990; free-free absorption, Barvainis & Lonsdale 1998) have not been successful.

The compact radio source in Mrk 348 is unique among Seyferts in that it is very bright and extremely variable. The observations presented here were made during a local minimum, when the total continuum flux density at 22 GHz was $``$0.6 Jy.

## 2. Observations

The initial detection of the flaring maser in Mrk 348 using the Effelsberg 100m telescope took place in 2000 March (see Fig. 1, 2nd profile). Re-analysis of previous data on this source (shown in the top profile, Fig. 1) indicates that the maser was also present but only marginally detectable in late 1997. The H2O maser line in Mrk 348 is extremely broad, with a FWHM of $``$130 km s-1, though in many of the monitoring epochs the emission appears to consist of 2 lines which can be fit by Gaussian functions with FWHM of $``$60 km s-1 and $``$100 km s-1 respectively, separated by $``$70 km s-1. There are no detectable high velocity components within 1500 km s-1 on either side of the strong emission line.

Monitoring through June 2000 showed that the maser again decreased to its original level within 2 months. The June profile (bottom profile, Fig. 1) indicates a peak flux of $``$18 mJy. This is consistent with the VLBA results shown in Fig. 2, taken on June 10, 2000. The monitoring observations indicated that the line was too broad (FWZP$`>`$250 km s-1) to fit in a single 16 MHz VLBA IF, so 2 IFs of 16 MHz each were used, overlapped by 5 MHz. Following calibration, the overlapping channels were removed and the 2 IFs were added together to yield a single cube of 174 channels covering 23 MHz. Line profiles of the resulting cube are shown in Fig. 2, superimposed on a continuum map made from 20 line-free channels. The maser emission is clearly seen to lie along the line of sight to the jet, rather than the core which is thought to lie nearer the bright component of the continuum source. No maser emission is seen toward any other region of the radio source, only toward the jet, and here the maser emission is unresolved at our angular resolution of 0.42$`\times `$0.76 mas. This corresponds to a linear size of less than 0.25 pc (assuming H0=75 km s-1Mpc-1). The Gaussian fit to the line shown in the first profile has an amplitude of 16$`\pm `$2 mJy and an integrated flux of 2.4$`\pm `$0.3 Jy/beam/km s-1, indicating that all of the flux measured in the Effelsberg 19 June observation has been recovered. The FWHM is 142$`\pm `$9 km s-1 centered on VLSR=4640$`\pm `$2 km s-1, consistent with the single-dish measurements and redshifted by 131 km s-1 with respect to the systemic velocity. This redshift might indicate that it is the receding jet that impacts the molecular cloud, and that the masing gas is being entrained. A tentative 2 component fit to the H2O profile yields a narrower line at 4683 km s-1 with FWHM $``$60 km s-1 and amplitude $``$8 mJy, and a broader line at 4620 km s-1 with FWHM $``$100 km s-1 and amplitude $``$12 mJy, again consistent with our single-dish measurements.

Resumption of the monitoring program at Effelsberg in December 2001 showed little change in the line flux, then in February another flare appears to have begun. The observations made in early April 2001 (see Fig. 1) show that the FWHM in the line decreased to 65 km s-1, and the centroid increased from 4642 km s-1 to $``$4665 km s-1. Thus it seems likely that only one of the two line components was flaring. Following this peak on 2 April, the line again began to decrease in flux. We are continuing our single-dish monitoring program at the Effelsberg telescope, and will attempt to observe with the VLBA during the next maser flare.

## 3. Conclusions

During early 2000, the H2O emission toward Mrk 348 showed a dramatic intensity increase which coincided with a significant increase in the flux of the nuclear radio continuum source. The unusual line profile leads us to suspect that this source, and possibly NGC 1052, might belong to a class of megamaser galaxies in which the amplified emission is the result of an interaction between the radio jet and an encroaching molecular cloud, rather than occurring in a circumnuclear disk. Analysis of our recent VLBA observations indicates that the emission does indeed arise along the line of sight to the jet in Mrk 348 (see Fig. 2), confirming this prediction. The very high linewidth occurring on such small spatial scales indicates that the H2O emission arises from a shocked region at the interface between the energetic jet material and the molecular gas in the cloud where the jet is boring through. This hypothesis is supported by the spectral evolution of the continuum source (Brunthaler, priv. comm.), which showed an inverted radio spectrum with a peak at 22 GHz, later shifting to lower frequencies. By analogy to IIIZw2 (Brunthaler et al. 2000) this would indicate the formation and evolution of very compact hotspots propagating through a dense medium. In this scenario, the recent high frequency radio continuum flare and the apparent northward movement of the brightest continuum component are attributable to the impact between the jet and the molecular cloud. This is consistent with recent VLBA continuum monitoring observations that show that the continuum source does not appear to be expanding (Ulvestad, priv.comm.), as one would expect from a hotspot or working surface. Rather the jet component toward which the maser appears might be associated with a stationary shock. The very close temporal correlation between the flaring activity in the maser emission and the continuum flare further suggest that the masing region and the continuum hotspot are nearly coincident and may be different manifestations of the same dynamical events.

## References

Barvainis, R. & Lonsdale, C. 1998, AJ, 115, 885

Brunthaler, A., Falcke, H., Bower, G. C., Aller, M. F., Aller, H. D.,Teräsranta, H., Lobanov, A. P., Krichbaum, T. P. & Patnaik, A. R. 2000, A&A, 357, L45

Capetti, A., Axon, D. J., Macchetto, F., Sparks, W. B. & Boksenberg, A. 1996, ApJ, 469, 554

Falcke, H. Henkel, C., Peck, A. B., Hagiwara, Y., Prieto, M. A. & Gallimore, J. F. 2000, A&A, 358, L17

Gallimore, J. F., Baum, S. A., O’Dea C. P., Pedlar, A. & Brinks, E. 1999, ApJ, 524, 684

Huchra, J. P., Vogeley, M. S. & Geller, M. J. 1999, ApJS, 121, 287

Simkin, S. M, Su, H.-J., van Gorkom, J. & Hibbard, J. 1987, Science, 235, 1367

Simpson, C., Mulchaey, J. S., Wilson, A. S., Ward, M. J. & Alonso-Herrero, A. 1996, ApJ, 457, L19

Taniguchi, Y., Kameya, O., Nakai, N. & Kawara, N. 1990, ApJ, 358, 132

Ulvestad, J. S., Wrobel, J. M., Roy, A. L., Wilson, A. S., Falcke, H & Krichbaum, T. P. 1999, ApJ, 517, L81


text:
# Efficient Quantum Computing of Complex Dynamics

## Abstract

We propose a quantum algorithm which uses the number of qubits in an optimal way and efficiently simulates a physical model with rich and complex dynamics described by the quantum sawtooth map. The numerical study of the effect of static imperfections in the quantum computer hardware shows that the main elements of the phase space structures are accurately reproduced up to a time scale which is polynomial in the number of qubits. The errors generated by these imperfections are more dangerous than the errors of random noise in gate operations.

When applied to computation, quantum mechanics opens completely new perspectives: a quantum computer, if constructed, could perform certain computations faster than classical computers, exploiting quantum mechanical features such as entanglement or superposition . Shor discovered a quantum algorithm  which factorizes large integers exponentially faster than any known classical algorithm. It was also shown by Grover  that the search of an item in an unstructured database can be done with a square root speed up over any classical algorithm. However, at present only a few quantum algorithms are known which simulate physical systems with exponential efficiency. Such systems include certain many-body problems , spin lattices , and models of quantum chaos .

In this Letter, we present a quantum algorithm which computes the time evolution of the quantum sawtooth map, exponentially faster than any classical computation. This model has rich and complex dynamics and finds various applications, e.g. for dynamical localization in billiards . Our algorithm, based on the Quantum Fourier Transform (QFT) , simulates the dynamics of a system with $`N`$ levels in $`O((\mathrm{log}_2N)^2)`$ operations per map iteration, while a classical computer, which performs Fast Fourier Transforms (FFT), requires $`O(N\mathrm{log}_2N)`$ operations. A further striking advantage of the algorithm is the optimum utilization of qubits: one needs only $`n_q=\mathrm{log}_2N`$ qubits (without any extra work space). We demonstrate that complex phase space structures can be simulated with less that 10 qubits, while about 40 qubits would allow one to make computations inaccessible to present-day supercomputers. This is particularly important, since experiments with few qubits are being performed at present , in particular the QFT was implemented on a three qubit nuclear magnetic resonance quantum computer . For this reason the investigation of this interesting physical system will be accessible to first quantum computers, operating with few qubits and for which large-scale computations like integer factoring are not possible.

The classical sawtooth map is given by

$$\overline{n}=n+k(\theta \pi ),\overline{\theta }=\theta +T\overline{n},$$
(1)
where $`(n,\theta )`$ are conjugated action-angle variables ($`0\theta <2\pi `$), and the bars denote the variables after one map iteration. Introducing the rescaled variable $`p=Tn`$, one can see that the classical dynamics depends only on the single parameter $`K=kT`$, so that the motion is stable for $`4<K<0`$ and completely chaotic for $`K<4`$ and $`K>0`$. For such a discontinuous map the conditions of the Kolmogorov-Arnold-Moser (KAM) theorem are not satisfied and therefore, for any $`K0`$, the motion is not bounded by KAM tori. The map (1) can be studied on the cylinder ($`p(\mathrm{},+\mathrm{})`$), which can also be closed to form a torus of length $`2\pi L`$, where $`L`$ is an integer. For any $`K>0`$, one has normal diffusion: $`<(\mathrm{\Delta }p)^2>D(K)t`$, where $`t`$ is the discrete time measured in units of map iterations and the average $`<\mathrm{}>`$ is performed over an ensemble of particles with initial momentum $`p_0`$ and random phases $`0\theta <2\pi `$. It is possible to distinguish two different dynamical regimes : for $`K>1`$, the diffusion coefficient is well approximated by the random phase approximation, $`D(K)(\pi ^2/3)K^2`$, while for $`0<K<1`$ diffusion is slowed down, $`D(K)3.3K^{5/2}`$, due to the sticking of trajectories close to broken tori (cantori). For $`4<K<0`$ the motion is stable, the phase space has a complex structure of elliptic islands down to smaller and smaller scales, and we observed anomalous diffusion, $`<(\mathrm{\Delta }p)^2>t^\alpha `$, (for example, $`\alpha =0.57`$ when $`K=0.1`$).

The quantum evolution on one map iteration is described by a unitary operator $`\widehat{U}`$ acting on the wave function $`\psi `$:

$$\overline{\psi }=\widehat{U}\psi =e^{iT\widehat{n}^2/2}e^{ik(\widehat{\theta }\pi )^2/2}\psi ,$$
(2)
where $`\widehat{n}=i/\theta `$ (we set $`\mathrm{}=1`$). The classical limit corresponds to $`k\mathrm{}`$, $`T0`$, and $`K=kT=\text{const}`$. In this quantum model one can observe important physical phenomena like dynamical localization . Indeed, due to quantum interference effects, the chaotic diffusion in momentum is suppressed, in a way similar to Anderson localization in disordered solids. Also in the vicinity of a broken KAM torus, cantori localization takes place, since a cantorus starts to act as a perfect barrier to quantum wave packet evolution, if the flux through it becomes less than $`\mathrm{}`$ .

The most efficient way to simulate the quantum dynamics (2) on a classical computer is based on forward/backward FFT between $`\theta `$ and $`n`$ representations. This is advantageous because the evolution operator $`\widehat{U}`$ is the product of two unitary operators, $`\widehat{U}_k=\mathrm{exp}(ik(\widehat{\theta }\pi )^2/2)`$ (kick) and $`\widehat{U}_T=\mathrm{exp}(iT\widehat{n}^2/2)`$ (free rotation), which are diagonal in the $`\theta `$ and $`n`$ representation, respectively. Therefore, for a system with $`N`$ levels, the one map iteration (2) requires two FFT and two diagonal multiplications and can be performed in $`O(N\mathrm{log}_2(N))`$ operations. The dynamics (2) can be simulated exponentially faster on a quantum computer with $`n_q=\mathrm{log}_2N`$ qubits by means of the following quantum algorithm: (i) the wave function $`|\psi =_{n=0}^{N1}a_n|n`$ (given in the $`n`$ representation) is multiplied by $`\widehat{U}_T`$, so that $`\widehat{U}_T|\psi =_na_n\mathrm{exp}(iTn^2/2)|n`$. This step can be done in $`n_q^2`$ controlled-phase shift gates, as explained in ; (ii) one can get the wave function in the $`\theta `$ representation via the QFT , which requires $`n_q`$ single-qubit (Hadamard) gates and $`n_q(n_q1)/2`$ two-qubit gates (controlled phase-shifts); (iii) the action of $`\widehat{U}_k`$ is diagonal in the angle representation and can be simulated in a way similar to (i) in $`n_q^2`$ two-qubit gates (we note that this is possible thanks to the particular form of $`\widehat{U}_k`$ for the sawtooth map); (iv) we go back to the momentum basis performing backward QFT in $`n_q(n_q+1)/2`$ gates. Therefore the whole algorithm requires $`n_g=3n_q^2+n_q`$ quantum gates per map iteration.

We model the quantum computer hardware as a two-dimensional lattice of qubits (spin halves) with static fluctuations/imperfections in the individual qubit energies and residual short-range inter-qubit couplings. The model is described by the following many-body Hamiltonian (see also ):

$$H_\text{s}=\underset{i}{}(\mathrm{\Delta }_0+\delta _i)\sigma _i^z+\underset{i<j}{}J_{ij}\sigma _i^x\sigma _j^x,$$
(3)
where the $`\sigma _i`$ are the Pauli matrices for the qubit $`i`$, and $`\mathrm{\Delta }_0`$ is the average level spacing for one qubit. The second sum in (3) runs over nearest-neighbor qubit pairs, and $`\delta _i`$, $`J_{ij}`$ are randomly and uniformly distributed in the intervals $`[\delta /2,\delta /2]`$ and $`[J,J]`$, respectively. We study numerically the many-body dynamics of the quantum computer (3) running the quantum algorithm described above. The algorithm is realized by a sequence of instantaneous and perfect one- and two-qubit gates, separated by a time interval $`\tau _g`$, during which the Hamiltonian (3) gives unwanted phase rotations and qubit couplings. A special one-qubit rotation is applied to eliminate the average phase accumulation given by $`\mathrm{\Delta }_0`$. We consider the case of short range inter-qubit interaction, in which two-qubit gates can be performed only between nearest neighbor pairs. As a result, for a square lattice of $`n_q`$ qubits, we make $`O(n_q^{1/2})`$ swap gates to obtain a generic two-qubit gate. Therefore a single iteration of the map (2) requires $`n_g=O(n_q^{5/2})`$ gates.

In this Letter, we study the sawtooth map in the anomalous diffusion regime, with $`K=0.1`$, $`\pi p<\pi `$ (torus geometry). The classical limit is obtained by increasing the number of qubits $`n_q`$, with $`T=2\pi /N`$ ($`k=K/T`$, $`N/2n<N/2`$). We consider as initial state at time $`t=0`$ a momentum eigenstate, $`|\psi (0)=|n_0`$, with $`n_0=[0.38N]`$. Such a state can be prepared in $`O(n_q)`$ one-qubit rotations starting from the ground state $`|0,\mathrm{},0`$. The dynamics of the sawtooth map reveals the complexity of the phase space structure, as shown by the Husimi functions of Fig.1, taken after $`1000`$ map iterations. We note that $`n_q=6`$ qubits are sufficient to observe the quantum localization of the anomalous diffusive propagation through hierarchical integrable islands. At $`n_q=9`$ one can see the appearance of integrable islands, and at $`n_q=16`$ the quantum Husimi function explores the complex hierarchical structure of the classical phase space. The effect of static imperfections (3) for the operability of the quantum computer is shown in Fig.1 (right column). The data are shown for $`J=0`$ (we observed similar structures for $`J=\delta `$). The main features of the wave packet dynamics remain evident even in the presence of significant imperfections, characterized by the dimensionless strength $`ϵ=\delta \tau _g`$. The main manifestation of imperfections is the injection of quantum probability inside integrable islands. This creates characteristic concentric ellipses, which follow classical periodic orbits moving inside integrable islands. These structures become more and more pronounced with the increase of $`n_q`$. Thus quantum errors strongly affect the quantum tunneling inside integrable islands, which in a pure system drops exponentially ($`\mathrm{exp}(CN)`$, $`C=`$const). It is interesting to stress that the effect of quantum errors is qualitatively different from the classical round-off errors, which produce only slow diffusive spreading inside integrable islands (see Fig.1). This difference is related to the fact that spin flips in quantum computation can make direct transfer of probability on a large distance in phase space.

A more quantitative study of the effect of static imperfections is obtained from the fidelity of quantum computation, defined by $`f(t)=|\psi _ϵ(t)|\psi _0(t)|^2`$, where $`\psi _ϵ(t)`$ is the actual quantum wave function in the presence of imperfections and $`\psi _0(t)`$ is the quantum state for a perfect computation. The data of Fig.2 show that $`f(t)`$ drops

It is interesting to compare the effect of static imperfections generated by the Hamiltonian (3) in the regime of perfect gates with the case when the static imperfections are absent but the computer operates with noisy gates. To model the noisy gates we put $`J=0`$ and $`\delta _i`$ fluctuating randomly and independently from one gate to another in the interval $`[\delta /2,\delta /2]`$. In this case the fidelity drops in a qualitatively different way $`f(t)\mathrm{exp}(Bt)`$ (data not shown). This corresponds to the Fermi golden rule regime, where the probability at a given state decays exponentially with time. At each gate operation a probability of order $`ϵ^2`$ is transferred to other states, so that the fidelity time scale is given by $`t_f1/(ϵ^2n_g)`$. This estimate is in good agreement with numerical data shown in Fig.3. We found the same behavior for $`J=\delta =0`$ and noisy gates with unitary rotations on a random angle $`ϵ`$. The dependence on $`ϵ`$ is qualitatively different comparing to the static imperfection case (3) discussed above. We stress that static imperfections give shorter times $`t_f`$ and therefore are more dangerous for quantum computation. We note that the quadratic and linear regimes of fidelity decrease with time were respectively discussed in Refs. .

In summary, we have shown that a realistic quantum computer can simulate with exponential efficiency complex dynamics with rich phase space structures. Our studies demonstrate that the main structures inside the localization domain are rather stable in the presence of static imperfections. We show that the errors from static imperfections are stronger than the errors of noisy gates. Of course it is not possible to extract all exponentially large information hidden in the wave function with $`2^{n_q}`$ states. However, it is possible to have access to coarse grained information. For example from a polynomial number of measurements one can obtain the probability distribution over momentum (or angle) states. This allows one to study the anomalous diffusion in the deep semiclassical regime. Such an information is not accessible for classical computers which cannot simulate more than $`2^{40}`$ quantum states. Recently an efficient algorithm was proposed in Ref. , which allows one to measure the value of the Wigner function at a chosen phase space point. This can also provide important new informations about quantum states in systems with hierarchical phase space structures. We believe that in the near future our results can be experimentally observed in quantum computers operating with a small number of qubits.

This work was supported in part by the EC RTN network contract HPRN-CT-2000-0156 and (for D.L.S.) by the NSA and ARDA under ARO contract No. DAAD19-01-1-0553. Support from the PA INFM “Quantum transport and classical chaos” is gratefully acknowledged.


text:
# Soft Cores in Late-Type Dwarf and LSB Galaxies from H𝛼 Observations

## 1. Introduction

The current scenario for structure formation in the universe is based on the inflationary Cold Dark Matter (CDM) theory. High-resolution N-body simulations predict, for the virialized haloes, density profiles with central cusps with $`\rho r^{1.5}`$ when r $``$ 0 (Moore et al. 1999), even more cuspy than previous simulation results (Navarro, Frenk, and White 1997, hereafter NFW). These profiles are in conflict with the observations: rotation curves of late-type dwarf and LSB galaxies (DM dominated galaxies) seem to rule out singular halo profiles and are in good agreement with density profiles characterized by soft cores, central regions with constant density (de Blok & McGaugh 1997). On the other hand, van den Bosch & Swaters (2001) claim that these HI rotation curves are affected by beam smearing, concluding that the observed rotation curves are consistent with NFW model. A good spatial resolution (easily achievable with H$`\alpha `$ observations) is necessary to explore the galaxy inner regions and rule out beam smearing. We show some preliminary results of our H$`\alpha `$ observations for some dwarf and LSB galaxies. This work assumes $`\mathrm{H}_o`$=75 km s-1 Mpc-1.

## 2. Observations

Observations of late-type dwarf and LSB galaxies were carried out at the 3.6 m Telescopio Nazionale Galileo (TNG) both in photometric and spectroscopic mode using the instrument Dolores (scale = 0.275 Åpxl-1, dispersion = 0.8 Åpxl-1, $`\lambda (6200,7800)`$ Å). In Table 1 we list the properties of some of the observed galaxies.

Our H$`\alpha `$ data (stars in the figures) are plotted together with HI data (from literature, empty squares). As one can see, the HI data suffer of very little beam smearing. The solid lines are the best fits obtained by using the Burkert profile; the dashed lines are NFW profiles and the dotted-dashed lines are Moore profiles (CDM model). The numbers in parenthesis represent the normalized $`\chi ^2`$ for each fit. For each galaxy is stated the spatial resolution for both the H$`\alpha `$ and the HI data. For all the four galaxies, the Burkert profile is the best fit, while the Moore profile does not match data. New observations with high resolution new technology grating (VPHG, Conconi et al. 2001) are being planned.

## 3. Conclusions

Our preliminary analysis of H$`\alpha `$ observations seems to rule out the CDM model in the inner regions of LSB and dwarf galaxies, showing evidence of soft cores.

## References

de Blok, W.J.G., McGaugh, S.S. 1997, MNRAS, 290, 533
Moore,B., Guinn,T., Governato,F., Stadel,J., Lake,G. 1999, MNRAS, 310, 1147
Navarro, J.F., Frenk, C.S., White, S.D.M 1997, ApJ, 490, 493
van den Bosch, F.C., Swaters, R.A. 2001, MNRAS, 325, 1017
Conconi, P. et al. 2001 (SPIE Meeting, San Diego, U.S.A.)


text:
# Model of ejection of matter from nonstationary dense stellar clusters and chaotic motion of gravitating shells.

## 1 Introduction

.

Dynamical processes around supermassive black holes in quasars, blazars and active galactic nuclei (AGN) are characterised by violent phenomena, leading to formation of jets and other outbursts. Jet formation is usually connected with the processes in the magnetized accretion disks (Lovelace, 1976; Bisnovatyi-Kogan & Blinnikov, 1976). Formation of quasispherical outbursts, which are probably observed in quasars with broad absorption lines, could be connected with another mechanism. Here we consider a possibility of a shell outburst from a supermassive black holes (SBH) surrounded by a dense massive stellar cluster, basing on a pure ballistic interaction of gravitating shells oscillating around SBH.

Ballistic ejection may be responsible for appearance of stars and observed SN in the intergalactic space, as well as the existence of stars between galaxies, suggested by many authors.

Investigation of spherical stellar clusters using shell approximation was started by Hénon (1964), and than have been successfully applied for investigation of the stability Hénon (1973), violent relaxation and collapse (Hénon 1964; Gott 1975), leading to formation af a stationary cluster. Investigation of the evolution of spherical stellar cluster with account of different physical processes was done on the base of a shell model in the classical serie of papers of L. Spitzer and his coauthors (Spitzer & Hart 1971a,b; Spitzer & Thuan 1972; Spitzer & Shapiro 1972; Spitzer & Chevalier 1973; Spitzer & Shull 1975a,b; Spitzer & Mathiew 1980).

Numerical calculations of a collapse of stellar clusters in a shell approximation (Yangurazova & Bisnovatyi-Kogan 1984; Bisnovatyi-Kogan & Yangurazova 1987,1988) had shown, that even if all shells are initially gravitationally bound, after a number of intersections some shells obtain sufficient energy to become unbound, and to be thrown to the infinity. In the Newtonian gravity the remnant is formed as a stationary stellar cluster, and in general relativity SBH may be formed as a remnant.

Formation of outbursts due to the ballistic interactions may happen as a result of intersection of shells, oscillating around SBH. In the smooth cluster with or without SBH in the centre stars evaporate, due to stellar encounters mainly with small kinetic energy, and formation of rapidly expelling stars is about 100 times less probable due to predominance of gravitational encounters with small momentum transfer Ambartsumian (1938). If the cluster is strongly aggregated, consisting of few compact pieces, the encounters between these pieces are quite different, collision with large momentum transfer are becoming probable. In this case gravitational interaction between compact pieces may lead to the outburst with large velocity, and when such event happens in the vicinity of SBH the velocity could become a considerable part of the light velocity $`c`$.

Such kind of situation in principle may be produced as a result of galactic collisions with close encounter of nuclei, when one of the nuclei strips the matter from the companion in the form of a collapsing shells. Interaction of such shell with the stellar cluster may lead not only to collapse into SBH, but also to the reverse phenomena: expelling of the shell with a speed, much larger than the average speed of stars in the cluster. The shell is not falling into SBH due to high angular momentum of its stars.

Even more important example of a quasi-spherical mass ejection is a relativistic collapse of a spherical stellar system, which is considered (Zeldovich & Podurets, 1965; Lightman & Shapiro 1978; Ipser 1980; Rasio et al. 1989) as the main mechanism of a formation of supermassive black holes in the galactic centers. Approximation of such collapse by consideration of spherical shells is the simplest approach, which reflects all important features of such collapse (Gott 1975; Bisnovatyi-Kogan & Yangurazova 1987, 1988). Our consideration is related to the motion of stars (shells) which remain outside the newly formed supermassive black hole.

Here we consider a simplified problem of a motion of two massive spherical shells, each consisting of stars with the same specific angular momentum and energies, around SBH. The aim of this paper is to give an elementary treatment of the ballistic ejection, and to estimate the maximal efficiency of this ejection. Note, that in a more complicated case of numerous shell intersection this elementary act is a key process of the energy exchange between stars and of matter ejection. This is also the elementary process leading to the violent relaxation of the cluster Lynden-Bell (1967), studied in the shell approximation by Gott (1975) and Yangurazova & Bisnovatyi-Kogan (1984). Chaotization of the motion of two gravitating intersecting shells appears as a by-prodiuct of our consideration. The appearence of chaos in such a system, which is described by a set of simple algebraic equations, is a rare example. For the oversimplified case with a pure radial motion and reflecting inner boundary the chaotic shell motion was found by Miller and Youngkins (1997). Here we consider realistic shells consisting of stars moving along eliptical trajectories, and their chaotic properties would be examined in more details in another work. We find conditions at which one of two shells is expelling to infinity taking energy from another shell. We find a maximum of the velocity of the outbursting shell as a function of the ratio of its mass $`m`$ to the mass $`M`$ of SBH using only Newtonian theory of the shell’s motion.

However we introduce cutoff fixing the minimal radius of the potentially outbursting shell $`r_m`$ on the level of few $`r_g=2GM/c^2`$. We show, that for equal masses of two shells the expelling velocity reach the value $`v_{max}0.3547v_p`$ at $`m/M=1.0`$ , and $`v_{max}0.3v_p`$ was obtained at the masses of shells $`0.25÷1.5M`$, where the parabolic velocity of a shell in the point of a smallest distance to the black hole $`v_p`$ may be of a considerable part of $`c`$.

In sections 2,3 we describe outburst effect and in sections 4 we present the evidence of the chaos in the system of intersecting shells. The exact solution of these problems in the context of General Relativity will be presented in the subsequent paper.

## 2 Two shells around SBH.

Physically the nature of the ballistic ejection is based on the following four subsequent events. The outer shell is accellerated moving to the center in a strong gravitational field of a central body and inner shell. Somewhere near the inner minimal radius of the trajectory shells intersect. After that the former outer shell is deccellerating moving from the center in the weaker gravitationsl field of only one central body. The second intersection happens somwhere far from the center. That may result in the situation when the total energy (negative) of the initially gravitationally bound outer shell is becoming positive as a result of two subsequent intersections with another shell. The quantitive analysis of this procrss is done in this section.

Equation of motion of a shell with mass $`m`$ and total conserved energy $`E`$ in the field of a central body with mass $`M`$ is

$$E=\frac{mv^2}{2}\frac{Gm(M+m/2)}{r}+\frac{J^2m}{2r^2},$$
(1)
where $`v=dr/dt`$ is the radial velocity of the shell and $`J^2m/2r^2`$ is the total kinetic energy of tangential motions of all particles, the shell make-up from. The constant $`J>0`$ has that interpretation that $`Jm`$ is the sum of the absolute values of the angular momenta of all particles.The term $`m/2`$ in (1) is due to the self-gravity of the shell.

Let’s consider two shells with parameters $`m_1,J_1`$ and $`m_2,J_2`$ moving around SBH with mass $`M`$. Let the shell ”1” be initially outer and the shell ”2” be the inner one. Then equations of motion are:

$$E_{1(0)}=\frac{m_1v_{1(0)}^2}{2}\frac{Gm_1(M+m_1/2+m_2)}{r}+\frac{J_1^2m_1}{2r^2},$$
(2)
$$E_{2(0)}=\frac{m_2v_{2(0)}^2}{2}\frac{Gm_2(M+m_2/2)}{r}+\frac{J_2^2m_2}{2r^2}.$$
(3)
By the index (0) we mark the initial evolution stage before the first intersection of the shells. Assume that both shells are moving to the center and $`v_{1(0)}^2>v_{2(0)}^2`$. Such shells intersect each other at a some radius $`r=a_1`$ and at some time $`t=t_1`$ after wich the shell ”1” becomes inner and shell ”2” – outer. The equations of motion of the shells during this new stage (designate it by the index (1)) are:

$$E_{1(1)}=\frac{m_1v_{1(1)}^2}{2}\frac{Gm_1(M+m_1/2)}{r}+\frac{J_1^2m_1}{2r^2},$$
(4)
$$E_{2(1)}=\frac{m_2v_{2(1)}^2}{2}\frac{Gm_2(M+m_2/2+m_1)}{r}+\frac{J_2^2m_2}{2r^2}.$$
(5)
The crucial point now are the matching conditions at the intersection point $`r=a_1,t=t_1`$ where from one can obtain the initial data to the equations (4), (5) in order to define uniquely the evolution during this new stage. In Newtonian theory these conditions are:

$$E_{1(0)}+E_{2(0)}=E_{1(1)}+E_{2(1)};$$

$$v_{1(0)}(t_1)=v_{1(1)}(t_1);v_{2(0)}(t_1)=v_{2(1)}(t_1),$$
(6)
i.e. the conservation of the total energy of the system and continuity of the velocities through the intersection point. Expressing $`v_{1(0)}(t_1),v_{2(0)}(t_1)`$ from eqs. (2), (3) and $`v_{1(1)}(t_1),v_{2(1)}(t_1)`$ from eqs. (4), (5) and equating them following conditions (6) we get:

$$E_{1(1)}=E_{1(0)}+\frac{Gm_1m_2}{a_1};E_{2(1)}=E_{2(0)}\frac{Gm_1m_2}{a_1}.$$
(7)
At some point $`r=a_2,t=t_2`$ can happen the second intersection after which the shell ”1” again becomes the outer and shell ”2” becomes inner one. It is a simple task to write the equation of motion during this third stage (we designate it by the index (2)) and repeat the same procedure of matching velocities but now at point $`r=a_2,t=t_2`$. The result is:

$$E_{1(2)}=E_{1(1)}\frac{Gm_1m_2}{a_2}=E_{1(0)}+Gm_1m_2\left(\frac{1}{a_1}\frac{1}{a_2}\right),$$
(8)
$$E_{2(2)}=E_{2(1)}+\frac{Gm_1m_2}{a_2}=E_{2(0)}Gm_1m_2\left(\frac{1}{a_1}\frac{1}{a_2}\right).$$
(9)
Let us describe the situation, when one shell is ejected to infinity after intersection of to initially bound shells. We consider a case when $`a_2`$ is larger then $`a_1`$ so, that the second term in (8) has larger absolute value, then the first one, the first shell gains a positive energy and goes to infinity. Both shells have initial negative energies $`E_{1(0)}`$ and $`E_{2(0)}`$, but with small enough absolute values. The first shell takes the energy from the second one, which is becoming more bound with larger absolute value of the negative energy $`E_{2(2)}`$, according to (8).

## 3 Numerical solution.

Let’s illustrate the foregoing scenario by an exact particular example of two shells of equal masses. We choose parameters in the following way:

$$m_1=m_2=m,E_{1(0)}=E_{2(0)}=0,J_1<J_2,$$
(10)
In fact such exact solution represents the first approximation to the more general situation when $`E_{1(0)}`$ and $`E_{2(0)}`$ are nonzero (negative) but small in that sense that both modulus $`|E_{1(0)}|`$ and $`|E_{2(0)}|`$ are much less than $`Gm^2/a_1`$.

We assume that initially both shells are moving towards the center. It follows from eqs. (2), (3) that under condition (10) such shells will intersect inescapably at the point $`r=a_1,t=t_1`$. After the second intersection at $`r=a_2,t=t_2`$ the shell ”1” will be thrown to infinity with expelling velocity $`v_{1exp}`$. It follows from (8) that

$$v_{1exp}=v_{1(2)}|_r\mathrm{}=\sqrt{2Gm\left(\frac{1}{a_1}\frac{1}{a_2}\right)}.$$
(11)
In order to construct a solution with maximal possible $`v_{1exp}`$ we consider a case when initially inner shell ”2” reaches the inner turning point at minimal possible radius $`r=r_m`$ (see introduction), and intersects with the initially outer shell after, during its outward motion.

To find the second intersection point we need to solve equation of motion (4),(5) for shell’s trajectories after the first intersection. First, let us consider the bound shell ”2”. We apply eq. (3) to the point $`r=r_m`$ which gives (see Fig. 1.)

$$J_2^2=Gr_m(2M+m).$$
(12)
Substituting this expression for $`J_2^2`$ into (5) and taking into account (7), (10) we get a solution of the equation (5) in the following parametric form:

$$\begin{array}{c}r(\eta )=b_{2(1)}(1e_{2(1)}\mathrm{cos}\eta ),\\ t(\eta )=\sqrt{\frac{b_{2(1)}^3}{G(M+3m/2)}}(\eta e_{2(1)}\mathrm{sin}\eta )\delta _t.\end{array}$$
(13)
where the constants $`b_{2(1)}`$ and $`e_{2(1)}`$ are:

$$b_{2(1)}=\frac{M+3m/2}{2m}a_1,e_{2(1)}=\sqrt{1\frac{4m(M+m/2)}{(M+3m/2)^2}\frac{r_m}{a_1}},$$
(14)
on the choice of the constant $`\delta _t`$ see below.

Between the first and second intersection there exists the inner turning point of the shell ”1” (now inner shell). We take that additional restriction that the shell ”1” reaches this turning point also at the minimal possible radius $`r=r_m`$. It is easy to show that a wide class of solutions with such restriction really exists. At the point $`r=r_m`$ equation (4) gives:

$$J_1^2=Gr_m(2M+m)+2Gmr_m^2a_1^1.$$
(15)
We substitute this expression for $`J_1^2`$ into eq. (4) and again with conditions (7) and (10) obtain the following parametric form of the solution of this equation:

$$\begin{array}{c}r(\xi )=b_{1(1)}(e_{1(1)}\mathrm{cosh}\xi 1),\\ t(\xi )=\sqrt{\frac{b_{1(1)}^3}{G(M+m/2)}}(e_{1(1)}\mathrm{sinh}\xi \xi ),\end{array}$$
(16)
where the constant $`b_{1(1)}`$ and $`e_{1(1)}`$ are:

$$b_{1(1)}=\frac{M+m/2}{2m}a_1,e_{1(1)}=1+\frac{2m}{M+m/2}\frac{r_m}{a_1},$$
(17)
The first intersection point $`r=a_1,t=t_1`$ corresponds to some parameters $`\eta =\eta _1,\xi =\xi _1`$. At this point the time $`t(\eta _1)`$ calculated from (13) should coincide with time $`t(\xi _1)`$ calculated from (16). This fix the constant $`\delta _t`$:

$$\begin{array}{c}\delta _t=\sqrt{\frac{b_{2(1)}^3}{G(M+3m/2)}}(\eta _1e_{2(1)}\mathrm{sin}\eta _1)\\ \sqrt{\frac{b_{1(1)}^3}{G(M+m/2)}}(e_{1(1)}\mathrm{sinh}\xi _1\xi _1).\end{array}$$
(18)
Also the values of $`r(\eta _1)`$ from (13) and $`r(\xi _1)`$ from (16) should be the same and equal to $`a_1`$. From this follows parameters $`\eta _1`$ and $`\xi _1`$:

$$\mathrm{cosh}\xi _1=\frac{b_{1(1)}+a_1}{b_{1(1)}e_{1(1)}}$$
(19)
$$\mathrm{cos}\eta _1=\frac{b_{2(1)}a_1}{b_{2(1)}e_{2(1)}}$$
(20)
The second intersection point $`r=a_2`$ corresponds to parameters $`\eta =\eta _2,\xi =\xi _2`$. The application of the analogous coincidence conditions to the equations (13) and (16) at this second point gives:

$$a_2=b_{1(1)}(e_{1(1)}\mathrm{cosh}\xi _21)$$
(21)
$$a_2=b_{2(1)}(1e_{2(1)}\mathrm{cos}\eta _2)$$
(22)
$$\begin{array}{c}\sqrt{\frac{b_{1(1)}^3}{G(M+m/2)}}(e_{1(1)}\mathrm{sinh}\xi _2\xi _2e_{1(1)}\mathrm{sinh}\xi _1+\xi _1)=\\ \sqrt{\frac{b_{2(1)}^3}{G(M+3m/2)}}(\eta _2e_{2(1)}\mathrm{sin}\eta _2\eta _1+e_{2(1)}\mathrm{sin}\eta _1).\end{array}$$
(23)
From (14),(17) and (19)-(23) it is easy to see that all these relations determine $`a_2/a_1`$ as a function of $`a_1/r_m`$ and $`m/M`$.

We introduce now the ”parabolic” velocity $`v_p`$ of the any outer shell at the point of its minimal distance to the center $`r=r_m`$ as:

$$v_p=\sqrt{\frac{2G(M+3m/2)}{r_m}}.$$
(24)
Than we have

$$\frac{v_{1exp}}{v_p}=\sqrt{\frac{m}{M+3m/2}\left(\frac{r_m}{a_1}\frac{r_m}{a_2}\right)}.$$
(25)
Consequently this ratio is also a function of $`a_1/r_m`$ and $`m/M`$. We find numerically the value of $`a_1/r_m`$ maximizing the ratio $`v_{1exp}/v_p`$ which value is a function of $`m/M`$ only. Dependencies for $`v_{1exp}/v_p`$, $`a_2/a_1`$, and $`a_1/r_m`$ as functions of the parameter $`m/M`$ for these maximizing solutions are given in Fig. 2-4.

The numerical calculations (illustrated by the Fig. 2) shows that the runaway velocity $`v_{1exp}`$ reaches its maximal possible value at $`m/M1.0`$ and it is $`v_{1exp}0.3547v_p`$. If we consider the shells around SBH then the minimal radius $`r_m`$ of the shell ”1” orbit cannot be less then the two gravitational radiuses $`2r_g=4G(M+3m/2)/c^2`$. In the extreme case when $`r_m=2r_g`$ we have $`v_pc/\sqrt{2}`$ and for the maximal possible runaway velocity in the family of solutions, characterized by (12), (15) we get $`v_{1exp}0.25c`$.

## 4 Chaos in the shell motion.

The first evidence that the motion of two intersecting shells can show chaotic character was given by B.N. Miller & V.P. Youngkins (1997). They investigated the special case when the central body is absent and particles consisting the shells are moving only in radial direction (in our notation $`M=0`$ and $`J_1=J_2=0`$). This situation, however, cannot model astrophysical cluster with massive nuclei, and also the problem of the influence of central Newtonian singularity arise which need some additional care. In any case a study of more physically realistic models with nonzero $`M,J_1`$ and $`J_2`$ from the point of view of possible chaotic behaviour represents essential interest. We report here some results for such more general two shells model which was investigated in the previous sections but again for the shells with equal masses. We consider now only oscillatory regime of motion without any runaway effects.

The shell motion in the Newton gravitational field is completely regular, but at presence of intersections the picture changes qualitatively. The shell intersections result to chaos in their motions. Character of chaos depends, mainly, on mass ratio of a shell and a central body. At a small shell mass the motion of the shells occurs basically in the field of the central body, and after one intersection there is a little change in a trajectory of each shell. Consider a case when angular momentum parameters $`J`$ and energies of shells are close to each other, and as a result of two intersections energy is transferred to the shell ”1”. Since shell parameters changes are not large, the mutual positions and speeds of shells after intersections do not change strongly. Then at the next intersection the part of energy will again be transferred to the shell ”1”. After a large number of intersections, however, the mutual positions changes essentially. Then the energy starts to flow to the shell ”2”, consequently we observe beating in oscillation of shells. In the case when the masses of the intersecting shells are large enough, the energy is transferred very intensively between them. It changes trajectories of shells essentially so that the next intersection does not resemble at all the previous one. In this case a shell’s motion at once gets properties of a randomness.

Shell oscillations with mass ratio $`m/M=0.015`$ are presented in Fig.5. In the case of massive shells the exchange of energy between shells occurs more intensively. As a result we have the obviously expressed chaotic behaviour of shells. Chaotic shell oscillations with $`m/M=0.08`$ are presented in Fig.6. In Fig.6 a full randomness of behaviour of the shells is seen obviously.

Let us note, that when the exchange of energy between shells is sufficiently large the strong randomness appears. Shells may exchange their positions during the motion. The external shell may become inner one and vice versa. One of the shells may give up practically all its energy of the orbital motion, which it may loose at constant angular momentum, i.e. the orbits of its particles become almost circular. It may be seen in Figs.6,7 at those moments when the orbit of the inner shell fills a narrow band by itself.

Let us note that in our calculations the motion of the shells is very sensitive to initial parameters, what gives an additional evidence that intersecting shells represent a really chaotic system. When we change only precision of integrator from $`10^6`$ to $`3\times 10^7`$, after several intersection motion of shells becomes absolutely different from its initial behaviour in Fig.6 at the same initial conditions. If we change the angular momentum parameter of one of shell from 1.40 to 1.41, then in the same way as in previous case, motion of shells becomes absolutely different after several intersection. We can make statement that the behaviour of shells is unstable and is strongly affected by any perturbation in initial parameters of shell motion.

The example of a chaotic behaviour of two intersecting selfgraviting shells, moving in their own gravitational field, without a central mass, is represented in Fig. 7.

Chaotic behaviour of the intersecting selfgravitating shells in the Newton gravity is a rare example of a chaotic dynamical system, described fully analytically by algebraic relations.

## 5 Discussion

We have shown that pure gravitational interaction of two spherical shells and central object may lead to considerable energy exchange and ejection of one of the shell with a subrelativistic speed $`0.25c`$. Motion of stars around SBH may lead to similar effects due to their interaction, and we may expect formation of very rapidly moving stars, which could overcome the galactic gravity and runaway into the intergalactic space. In this way very extended galactic halo may be formed. Ballistic mechanism of energy exchange between gravitating particles may be important in the evolution of spatial distribution of cold dark matter, when gravitational instabilities are developing in masses, much less than galactic ones, and gravitationally bound massive objects, consisting mainly from the dark matter interacting only due to its gravitation, may be formed. Interesting examples of ballistic interaction of satellites with planets and stars, their acceleration and ballistic control are presented in the book of V.V. Beletskyi (1977).

## Acknowledgement

The work of G.S.B.-K. and M.V.B. was partly supported by RFBR grant 99-02-18180, and INTAS grant 00-491. We are grateful to the referee, Prof. Seppo Mikkola for useful comments.

Figure 1. Time dependence of radiuses r of two shells (in units $`r_m`$) on time t (in units $`r_m/v_p,v_p=\sqrt{2G(M+3m/2)/r_m}`$). The shells intersect at points $`a_1=1.1401`$, $`a_2=3.6316`$. Here $`m/M=0.1`$, and after the second intersection the first shell is running away with velocity at infinity $`v_{1exp}=0.23415v_p`$.

Figure 2. Dependence of the maximum runaway velocity $`v_{1exp}`$ in units $`v_p`$ on the ratio of the shell mass to the mass of central body $`m/M`$. The maximum value of $`v_{1exp}=0.3547`$ corresponds to $`m/M=1.0`$.

Figure 3. Dependence of the ratio of intersection radiuses $`a_2/a_1`$ on $`m/M`$ for the solutions with maximum runaway velocities.

Figure 4. The value $`a_1/r_m`$ (the ratio of the radius of the first intersection $`a_1`$ to the minimal radius of the approach to the central body $`r_m`$) as a function of $`m/M`$ for solutions with maximum runaway velocities.

Figure 5. Shell oscillations with mass ratio $`m/M=0.015`$.

Figure 6. The chaotic motion $`r_{1,2}(t)`$ of shells with masses $`m/M=0.08`$.

Figure 7. The chaotic motion $`r_{1,2}(t)`$ of two selfgravitating shells in their own gravitational field without a central mass.


text:
# 900-1200Å interstellar extinction in the Galaxy, LMC, and SMC 1footnote 11footnote 1Based on observations made with the NASA-CNES-CSA Far Ultraviolet Spectroscopic Explorer. FUSE is operated for NASA by the Johns Hopkins University under NASA contract NAS5-3298

## 1 Introduction

The interstellar extinction down to the Lyman limit is of interest for several reasons. Shortward of 1200Å, the extinction rises strongly and differences between galaxy environments appear to increase. This means that corrections to observed (UV restframe) fluxes are potentially very sensitive to the amount and type of interstellar extinction. This is important in understanding the energy budget of galaxies with young stellar populations, and at high redshifts in identifying the full population of young galaxies.

There have been many studies in the IUE/HST range (1200 - 3200Å) where it is known that the extinction curve varies widely between different environments (see e.g. Cardelli, Clayton, and Mathis 1989: CCM; Gordon and Clayton 1998). In particular, the UV extinction laws differ stongly between the Galaxy, LMC, and SMC, in the sense that lower abundance galaxies have higher UV extinction for a given visible reddening, and less pronounced 2175Å ‘bump’. This general rule is thought to indicate that the ISM composition or range of grain sizes changes with continued stellar enrichment. There are, however, variations within the local galaxies that must reflect a range of local ISM conditions or history. It is more difficult to obtain reddening for the much fainter stars in the M31 group galaxies, but work by Bianchi et al (1996) shows that the variations are seen there too.

In the range below 1200Å, observations are rare and the rising extinction makes it difficult to get data on significantly reddened stars. The HUT spectrograph on ASTRO yielded extinction curves in this range, but for only 2 stars (Clayton et al 1996). FUSE is more powerful, has higher resolution, and is building up an archive of spectra suitable for this work. While there are programs planned to observe reddened stars specifically for extinction studies, we present in this paper an initial, but significant derivation of extinction from lightly reddened stars observed for other purposes. This study describes some of the difficulties in deriving extinction curves in this wavelength range, and presents some results from the Galaxy and Magellanic Clouds, where the expected extinction differences are seen and measured.

The sample of stars was chosen to contain pairs matched in spectral type and luminosity with different B-V, assumed to indicate a different amount of extinction. Table 1 shows the stars selected from the FUSE PI team database of OB stars taken for stellar wind and O VI studies. These stars were selected for their spectral class and UV flux, and thus are not highly reddened, except for the Galactic stars, where brighter stars are available. All observations were taken with the large FUSE aperture and the data were processed uniformly with the CALFUSE data pipeline. The data were not obtained through the FUSE archive as most were still proprietary, so the data used may possibly differ from that which were delivered to the archive.

Since the Magellanic Cloud stars are lightly reddened, the precision of the colour photometry is an important part of the error budget of this work. We have used what we regard as the most reliable values from the literature in Table 1, but in all cases we investigated the effects of photometric errors within the range (typically a few hundredths of a magnitude) of values published for the individual stars, or other stars from the same papers. We also rely on spectral types to pick good pairs of stars to compare for differential reddening, and a similar inspection of published spectral types was made. The adopted spectral types are given in Table 1. As the FUSE archive grows, more and better pairs will become available, so this list represents only those that were suitable after the first year of operations.

## 2 Data processing

The FUSE spectra were edited of airglow emission lines, where significant (see Feldman et al 2001), before further processing. Using model spectra of H2 absorption (Tumlinson, private communication), the strongest features were edited out by placing a continuum across them. These absorptions are signifcant at the shorter wavelengths (below 1100Å), and usually in stars with EB-V=0.2 or more. See Hutchings et al (2001) for an illustration of the H2 spectrum. Strong wind lines, present in high luminosity O stars, were identified from the atlas of OB FUSE spectra which the FUSE team is preparing (or as illustrated by Fullerton et al 2000), and also edited out.

The edited spectra were then ratioed with the less reddened star as numerator, and the result was heavily smoothed (5Å or more boxcar). This was done separately for all channels of FUSE data, so that each pair of stars yielded a set of 8 ratioed spectra. These were then plotted in inverse wavelength, normalised to unit EB-V, as is conventional for extinction plots.

This approach usually resulted in 8 non-continuous lines, although all showed the expected rise with decreasing wavelength. The discontinuities arise from the H2 blanketing, which changes strongly with wavelength, flux errors at the edges of channels, FUSE detector ‘worm’ effect, and possibly extraction incompleteness. The small differences of reddening between the stars in our pairs means that these are amplified by large factors (typically 30) in normalising to unit EB-V. Nevertheless, the average of each emsemble of 8 plots appeared to represent credible extinction curves for the FUSE range.

In order to treat each spectrum more uniformly across the whole FUSE range, the spectra were edited of airglow and then all channels resampled to 1A bins and overplotted. H2 models were generated with the same sampling and plotted for a range of absorbing columns. Some stars have had detailed H2 modelling by Tumlinson et al (2001). This process allowed us to a) identify a smooth mean observed spectrum over the whole FUSE range, and b) estimate the continuum level of the star by correlating the strengths of the main H2 absorption features, and interpolating across any strong wind lines.

In this way, each spectrum was modelled as a continuum defined by 12 points across the full FUSE wavelength range, and these were then divided and normalised. The results were generally consistent with the separate channel plots, but without the discontinuities.

The absolute level of the extinction was derived by correcting the numerator spectrum by a factor which amounts to the ratio of unreddened V-band flux, taking into account the observed V magnitudes and B-V colours. This is standard practice and treats the numerator star as having total extinction equal to the difference of extinction between the two stars.

The results for the sample stars are shown in Figures 1 and 2. Because of the large multiplier applied for small extinction differences, we repeated the processing for a set of V and B-V values for the stars covering the range of published values, or their errors. Figures 1 and 2 illustrate these, and show the sensitivity of the derived extinction curves to these quantities.

The figures plot the overall continuum 12-point ‘models’. The results from the individual channel plots are illustrated for the nominal photometric values by the pairs of filled squares. These show the end points of linear fits to the ensembles of 8 individual channel plots. The agreement between these and the plotted continuum models gives a good idea of the reproducibility of our analysis by two different approaches.

## 3 Discussion

The individual extinction curves in Figures 1-2 were assessed for their sensitivity to the uncertainties in continuum fitting, and to the photometric values as described above. There is clearly a wide range of curves present, and some are very dependent on the photometric precision. Figure 3 shows our adopted mean curves for the three galaxies, derived from the individual ranges of extinction curves in Figures 1 and 2, weighted by the colour difference for each pair. These mean curves were calculated for 5 wavelength bins.

The plots also show schematic longer wavelength extinction curves for the Galaxy, LMC, and SMC, that represent those found in the literature. In Figure 3 we indicate the full range of curves for the LMC and SMC from Misselt et al (1999) and Gordon and Clayton (1999), while for simplicity the smaller diagrams just show ‘typical’ curves for each galaxy.

The curves in Figure 3 have ‘formal’ errors per point which are about the separation between the three curves. However, the fact that the curves are nearly parallel and do not cross, add significance to the overall family of curves and their trend with host galaxy abundances. Nevertheless, it is clear from Figures 1 and 2 that more reliable results will require stars with more extinction, more detailed continuum modelling, and more precise photometry.

Clayton et al (1996) published extinction curves for two LMC stars in this range from HUT spectra. These are more heavily reddened than ours, and their extinction curves closely bracket our mean Galactic curve, while lying above the CCM sketched model. However, the general steepening rise is reproduced very closely.

We thank Luciana Bianchi and Alex Fullerton for advice and assistance with the data, spectral types, and H2 modelling.

References

Bianchi L., Clayton G.C., Bohlin R.C., Hutchings J.B., Massey P., 1996, ApJ, 471, 203

Cardelli J.A., Clayton G.C., and Mathis J.S., 1989, ApJ, 345, 245 (CCM)

Clayton G.C. et al, 1996, ApJ, 460, 313

Gordon K.D., and Clayton G.C., 1998, ApJ, 500, 816

Feldman P.D., Sahnow D.S., Kruk J.W., Murphy E.M., Moos H.W., 2001, Journal of Geophysical Research 106, 8119

Fullerton A.W., et al 2000, ApJ, 538, L43

Hutchings J.B., Crampton D., Cowley A.P., Schmidtke P., Fullerton A.W., 2001, AJ, in press

Misselt K.A., Clayton G.C., Gordon K.D., 1999, ApJ, 515, 128

Tumlinson J. et al, 2001, preprint

Captions to figures

1. Extinction from Galactic star pairs as indicated. Schematic extinction curves for SMC, LMC, and CCM are shown in each panel for comparison. The nominal values of $`\mathrm{\Delta }_{BV}`$ are labelled, and the other values are given to show the sensitivity of the extinction curve to the photometric precision. In the case of the B0.5 pair, we also show the effect of a change of V magnitude with the same colour difference. The black squares show the end values for the extinction derived from the separate chennel plots (see text for description). The agreement between the squares and the nominal colour difference plots is an indication of the random errors in continuum fitting.

2. As for Figure 1, for star pairs in the LMC and SMC

3. Mean values from star pairs in the three host galaxies, from Figures 1 and 2. The formal errors on the five values plotted for each galaxy are $`\pm `$1.5 magnitudes, but the general behaviour of the set of curves suggest that their shapes and separation are significant. The star and open square symbols outline the extinction curves from 2 LMC stars by Clayton et al (1996). The longer wavelength schematic plots show the range of extinction curves derived in the SMC and LMC from references given, and the solid line is the CCM curve for Galactic extinction. Note that our results appear to extend the extinction curves for the three galaxies in a straightofrward way. The structure in the FUSE range curves is unlikely to be more than noise.