text:
# Equivariant quantization on quotients of simple Lie groups by reductive subgroups of maximal rank

## 1 Introduction

A quantum manifold is obtained from a usual manifold by replacing the original commutative function algebra with a deformed non-commutative algebra. If a manifold $`M`$ is equipped with an action of a Lie group $`G`$, i.e. $`M`$ is a $`G`$-manifold, it is natural to consider deformation quantizations of the function algebra which are equivariant with respect to the action of the group $`G`$ or the corresponding quantum group $`\mathrm{U}_h(𝔤)`$, , where $`𝔤`$ is the Lie algebra of $`G`$.

We consider here homogeneous $`G`$-manifolds of type $`M=G/K`$, where $`G`$ is a simple connected Lie group over $``$ and $`K`$ its reductive subgroup of maximal rank, i.e. $`K`$ contains a maximal torus of $`G`$. The Lie algebra of $`K`$, $`𝔨`$, is generated by a Cartan subalgebra of $`𝔤`$ and the set of contained in $`𝔨`$ root vectors. Such a class of manifolds includes semisimple coadjoint orbits of $`G`$ in $`𝔤^{}`$. For a coadjoint orbit, $`𝔨`$ is a Levi subalgebra generated by simple root vectors.

The notion of $`G`$-equivariant quantization is very natural and supposed to be known. So, we only give the definition of $`\mathrm{U}_h(𝔤)`$-equivariant quantization we admit in the paper.

###### Definition 1.1.

Let $`M`$ be a $`G`$-manifold. A $`\mathrm{U}_h(𝔤)`$-equivariant quantization on $`M`$ is the $`[[h]]`$-module $`(M)[[h]]`$ endowed with a deformed multiplication $`\mu _h`$ such that:

(i) $`\mu _h`$ defines a star-product on $`M`$, i.e. is of the form $`\mu _h=\mu +h\mu _1+h^2\mu _2+\mathrm{}`$, where $`\mu `$ is the usual commutative multiplication on $`(M)`$ and $`\mu _i`$, $`i1`$, are bidifferential operators on $`M`$ vanishing on constants;

(ii) the multiplication $`\mu _h`$ is $`\mathrm{U}_h(𝔤)`$-equivariant, i.e.

$$x\mu _h(ab)=\mu _h\mathrm{\Delta }_h(x)(ab),\text{for any}x\mathrm{U}_h(𝔤),a,bA_h,$$

where $`\mathrm{\Delta }_h`$ denotes the comultiplication in $`\mathrm{U}_h(𝔤)`$. We also suppose that the action of $`\mathrm{U}_h(𝔤)`$ on $`(M)[[h]]`$ is the natural extension of the initial action of $`\mathrm{U}(𝔤)`$ on $`(M)`$ in the following sense. Since $`𝔤`$ is simple, $`\mathrm{U}_h(𝔤)`$ is isomorphic to $`\mathrm{U}(𝔤)[[h]]`$ as an algebra, and the natural action of $`\mathrm{U}(𝔤)[[h]]`$ on $`(M)[[h]]`$ is induced from the action of $`\mathrm{U}(𝔤)`$ on $`(M)`$ by $`[[h]]`$-linearity.

The element $`\nu `$ defined as $`\nu (a,b)=\mu _1(a,b)\mu _1(b,a)`$ for $`a,b(M)`$ is a Poisson bracket on $`M`$. In particular, $`\nu `$ is a bivector field on $`M`$. We call $`\mu _h`$ a quantization corresponding to the Poisson bracket $`\nu `$. One can prove, , that the Poisson bracket of any $`\mathrm{U}_h(𝔤)`$-equivariant quantization has the form

$$\nu =sr_M,$$
(1)
where $`r_M`$ is the bivector field which is the image of the classical r-matrix $`r^2𝔤`$ related to $`\mathrm{U}_h(𝔤)`$ via the action map $`𝔤Vect(M)`$, and $`s`$ is a $`G`$-invariant bivector field on $`M`$ satisfying

$$[s,s]=\phi _M.$$
(2)
By $`\phi _M`$ we denote the $`3`$-vector field induced on $`M`$ by the $`𝔤`$-invariant element $`\phi =[r,r]^3𝔤`$ and by $`[,]`$ the Schouten bracket of two polyvector fields. We call a Poisson bracket of the form (1) admissible.

###### Remark 1.2.

It is easy to see that the admissible Poisson brackets are precisely those which makes $`M`$ a Poisson homogeneous manifold with respect to the Drinfeld–Sklyanin Poisson bracket $`r^rr^l`$ on $`G`$.

The admissible brackets are also related to the dynamical CYBE , .

It is convenient to introduce the following

###### Definition 1.3.

A $`G`$-invariant bivector field $`s`$ on $`M`$ is called *a $`\phi `$-bracket* if it satisfies $`(\text{2})`$.

The first result of the paper is a classification of invariant and admissible Poisson brackets on homogeneous manifolds under consideration. Equation $`(\text{1})`$ reduces the classification of admissible brackets on $`M`$ to a classification of $`\phi `$-brackets on $`M`$.

The second result of the paper is that for any admissible Poisson bracket on $`M`$ there exists the corresponding $`\mathrm{U}_h(𝔤)`$-equivariant quantization. Note that there exists a $`G`$\- invariant connection on $`M`$, since the stabilizer of a point of $`M`$ semisimply acts on $`𝔤`$. Moreover, since $`M`$ is $`G`$-homogeneous, an invariant Poisson bracket on $`M`$ has the same rank at any point of $`M`$. Applying, for example, the well known Fedosov method, we obtain that for any $`G`$-invariant Poisson bracket on $`M`$ the corresponding $`G`$-equivariant quantization always exists.

Among other works relevant to our study, we would like to mention the following. J. Donin and D. Gurevich  proved that on a semisimple coadjoint orbit there exists a quantization of the Sklyanin–Drinfeld Poisson bracket. J. Donin and S. Shnider  constructed one and two parameter equivariant quantization on any symmetric manifold. J. Donin, D. Gurevich, and S. Shnider  solved the problem of one and two parameter equivariant quantization on any semisimple coadjoint orbit, including a complete classification of corresponding Poisson brackets on them.

The paper is organized as follows. In Section 2, we introduce a group $`\mathrm{\Gamma }_𝔨`$ related to a reductive subalgebra $`𝔨𝔤`$ of maximal rank. It is the quotient of the lattice spanned by all roots of $`𝔤`$ by the sublattice spanned by roots of $`𝔨`$. We prove a statement about how the images of roots are lying in $`\mathrm{\Gamma }_𝔨`$.

In Section 3, using this statement we give a classification of invariant Poisson brackets and $`\phi `$-brackets on $`M`$.

In Section 4, we show that those brackets can be quantized.

## 2 Reductive subalgebras of maximal rank and group $`\mathrm{\Gamma }_𝔨`$

Let $`𝔤`$ be a reductive Lie algebra over $``$. Let us fix a Cartan subalgebra $`𝔥𝔤`$ and denote by $`\mathrm{\Omega }𝔥^{}`$ the root system of $`𝔤`$. We are interested in reductive subalgebras, $`𝔨`$, in $`𝔤`$ containing $`𝔥`$, which are called subalgebras of maximal rank. To any subalgebra of this type, it corresponds a subsystem of roots, $`\mathrm{P}`$. Namely, $`\alpha \mathrm{P}`$ if and only if the root vector $`E_\alpha 𝔨`$. Note that $`𝔨`$ is generated by $`𝔥`$ and the root vectors $`E_\alpha `$, $`\alpha \mathrm{P}`$.

###### Example 2.1.

Choose a system of simple roots $`\mathrm{\Pi }\mathrm{\Omega }`$. Fix a subset $`\mathrm{\Sigma }\mathrm{\Pi }`$. Let us put $`\mathrm{P}=(\mathrm{\Sigma })\mathrm{\Omega }`$, where $`(\mathrm{\Sigma })`$ denote the free abelian group (lattice) spanned by $`\mathrm{\Sigma }`$. Let $`𝔨`$ be the Lie subalgebra of $`𝔤`$ generated by the root vectors $`E_\alpha `$ for all $`\alpha \mathrm{P}`$ and by the Cartan subalgebra $`𝔥`$. Such $`𝔨`$ is called a Levi subalgebra of $`𝔤`$.

###### Example 2.2.

Choose a system of simple roots $`\mathrm{\Pi }\mathrm{\Omega }`$. Fix a root $`\alpha \mathrm{\Pi }`$ and a positive integer $`n`$. Consider the set $`\mathrm{P}`$ of all roots in $`\mathrm{\Omega }`$ whose coefficient before $`\alpha `$ in their expansion in basis $`\mathrm{\Pi }`$ is divisible by $`n`$. Let $`𝔨`$ be the Lie subalgebra of $`𝔤`$ generated by the root vectors $`E_\alpha `$ for all $`\alpha \mathrm{P}`$ and by the Cartan subalgebra $`𝔥`$.

The following is a modified version of the well known Dynkin theorem, .

###### Theorem 2.3.

For any reductive subalgebra $`𝔨𝔤`$ of maximal rank, there exists a chain $`𝔤𝔨_1𝔨_2\mathrm{}𝔨_k=𝔨`$ of Lie subalgebras such that the subalgebra $`𝔨_j`$ is obtained from $`𝔨_{j1}`$ by applying one of the procedures described in Examples 2.1 and 2.2.

Let us set $`\mathrm{\Gamma }_𝔨=(\mathrm{\Omega })/(\mathrm{P})`$ and denote by $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ the subgroup of elements of finite order in $`\mathrm{\Gamma }_𝔨`$. The following theorem describes the group $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ and images of roots in it.

###### Proposition 2.4.

The group $`\mathrm{\Gamma }_𝔨`$ is free, i.e. $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)=0`$, if and only if $`𝔨`$ is a Levi subalgebra.

$`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ is cyclic if and only if $`𝔨`$ is obtained as in Theorem 2.3 applying the procedure of Example 2.2 no more than once.

If $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ is cyclic then any its non-zero element is the image of a root.

If $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ is non-cyclic, then the kernel of any character of it contains the image of a root.

###### Proof.

(I) is obvious, since in this case $`\mathrm{P}`$ is generated by a set of roots, $`\mathrm{P}^{}`$, lying in a linear subspace of the vector space spanned by $`\mathrm{\Omega }`$. But it is known that such a $`\mathrm{P}^{}`$ can be included in a set of simple roots.

To prove other points, we use the structure theory of simple Lie algebras. First of all, we come to the following conclusions. Any cyclic summand in $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ is of order $`6`$. If $`𝔤`$ is a classical simple Lie algebra, the group $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ is isomorphic to a direct sum of copies of $`_2=/2`$. For exceptional simple Lie algebras $`𝔤`$, also may be $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)=_3`$, $`_4`$, $`_5`$, $`_2_3`$, $`_2_4`$, or $`_3_3`$. Using this structure of $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$, the proof of the theorem can be completed by a direct computation. ∎

## 3 Classification of invariant Poisson brackets and $`\phi `$-brackets on $`M`$

In this section, we suppose that $`G`$ is a simple connected Lie group and $`K`$ a reductive subgroup of maximal rank. We are going to describe invariant bivector fields, $`s`$, on $`M`$ satisfying the equation

$$[s,s]=\varkappa ^2\phi _M,$$
(3)
where $`\varkappa `$. Note that if $`s`$ satisfies (3) with $`\varkappa ^2=0`$, then it defines an invariant Poisson bracket, while if $`s`$ satisfies (3) with $`\varkappa ^2=1`$ then it defines a $`\phi `$-bracket on $`M`$.

Consider the natural projection $`\pi :GM=G/K`$. It induces the map $`\pi _{}:𝔤T_o`$ where $`T_o`$ is the tangent space to $`M`$ at the point $`o`$ being the image of unity. Since the $`\mathrm{ad}`$-action of $`𝔨`$ on $`𝔤`$ is semisimple, there exists an $`\mathrm{ad}(𝔨)`$-invariant subspace $`𝔪`$ of $`𝔤`$ complementary to $`𝔨`$, and one can identify $`T_o`$ and $`𝔪`$ by means of $`\pi _{}`$. It is easy to see that subspace $`𝔪`$ is uniquely defined and has a basis formed by elements $`E_\gamma `$, $`\gamma \mathrm{\Omega }P`$. Here $`E_\gamma `$ are root vectors satisfying $`(E_\gamma ,E_\gamma )=1`$ for the Killing form $`(,)`$.

It follows from the above that invariant bivector fields on $`M`$ correspond to $`𝔨`$-invariant bivectors of $`^2𝔪`$. It is proven in  and  that any bivector of $`^2𝔪`$ is $`𝔨`$-invariant if and only if it has the form

$$\underset{\alpha \mathrm{\Omega }\mathrm{P}}{}c(\overline{\alpha })E_\alpha E_\alpha ,$$
(4)
where $`\overline{\alpha }`$ is the image of $`\alpha `$ in $`\mathrm{\Gamma }_𝔨`$ and coefficients $`c(\overline{\alpha })`$ satisfy the condition $`c(\overline{\alpha })=c(\overline{\alpha })`$. This means that for invariant bivectors the coefficients before terms $`E_\alpha E_\alpha `$ and $`E_\beta E_\beta `$ are the same if $`\overline{\alpha }=\overline{\beta }`$.

The following result is proved in .

###### Proposition 3.1.

A bivector field $`s`$ presented in the form (4) satisfies (3) if and only if the coefficients $`c(\overline{\beta })`$ obey the following condition: if $`\overline{\alpha }+\overline{\beta }\mathrm{\Gamma }_𝔨`$ then

$$c(\overline{\alpha }+\overline{\beta })(c(\overline{\alpha })+c(\overline{\beta }))=c(\overline{\alpha })c(\overline{\beta })+\varkappa ^2.$$

Using the previous proposition, one can obtain the following description of $`\phi `$-brackets and invariant Poisson brackets on $`M`$, . See also , where this is done for admissible brackets on coadjoint orbits.

###### Proposition 3.2.

A $`\phi `$-bracket on $`M`$ is determined by the following data: a subgroup $`\mathrm{\Psi }\mathrm{\Gamma }_𝔨`$ such that $`\mathrm{\Gamma }_𝔨/\mathrm{\Psi }`$ is free, a group homomorphism $`\chi :\mathrm{\Psi }0`$ such that $`\chi (\overline{\alpha })1`$ for $`\overline{\alpha }`$ being the image of a root, and a linear ordering in the group $`\mathrm{\Gamma }_𝔨/\mathrm{\Psi }`$. The coefficients $`c(\overline{\alpha })`$, $`\overline{\alpha }\mathrm{\Gamma }_𝔨`$, from (4) are given by the formula:

$$c(\overline{\alpha })=\frac{\chi (\overline{\alpha })+1}{\chi (\overline{\alpha })1}\text{for }\overline{\alpha }\mathrm{\Psi }\text{,}$$

$`c(\overline{\alpha })`$ $`=`$ $`1\text{if the projection of }\overline{\alpha }\text{ in }\mathrm{\Gamma }_𝔨/\mathrm{\Psi }\text{ is positive,}`$
$`c(\overline{\alpha })`$ $`=`$ $`1\text{if the projection of }\overline{\alpha }\text{ in }\mathrm{\Gamma }_𝔨/\mathrm{\Psi }\text{ is negative.}`$
###### Proposition 3.3.

An invariant Poisson bracket on $`M`$ is determined by choosing a subgroup $`\mathrm{\Psi }\mathrm{\Gamma }_𝔨`$ containing no torsion elements and a group homomorphism $`\lambda :\mathrm{\Psi }`$ such that $`\lambda (\overline{\alpha })0`$ for $`\alpha `$ being the image of a root. The coefficients $`c(\overline{\alpha })`$ from (4) are given by the formula

$`c(\overline{\alpha })`$ $`=`$ $`{\displaystyle \frac{1}{\lambda (\overline{\alpha })}}\text{for }\overline{\alpha }\mathrm{\Psi }\text{,}`$
$`c(\overline{\alpha })`$ $`=`$ $`0\text{if the projection of }\overline{\alpha }\text{ in }\mathrm{\Gamma }_𝔨/\mathrm{\Psi }\text{ is not zero.}`$
###### Remark 3.4.

Another description of Poisson brackets on quotients of $`G`$ by reductive subgroups of maximal rank is given in .

Combining Proposition 2.4, 3.2, and 3.3, we obtain the following

###### Corollary 3.5.

A $`\phi `$-bracket (and, therefore, an admissible Poisson bracket) on $`M`$ exists if and only if $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ is cyclic.

A nonzero Poisson bracket on $`M`$ exists if and only if $`\mathrm{\Gamma }_𝔨`$ is infinite.

If $`\mathrm{\Gamma }_𝔨`$ is a finite non-cyclic group then there are neither invariant nor admissible Poisson brackets on $`M`$. Therefore, such an $`M`$ does not admit neither $`G`$\- nor $`\mathrm{U}_h(𝔤)`$-equivariant quantization.

###### Proof.

(I) If $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ is non-cyclic then, in virtue of point (IV) of Proposition 2.4, $`\mathrm{Ker}(\chi )`$ must contain the image of a root, what is impossible by Proposition 3.2. On the other hand, if $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ is cyclic, a required homomorphism $`\chi `$ is obviously exists.

(II) If $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ is finite, then $`\mathrm{\Psi }`$ from Proposition 3.3 must be finite too, hence $`\lambda (\mathrm{\Psi })=0`$. It follows that $`\lambda `$ applying to any image of a root gives zero.

(III) follows from (I) and (II). ∎

###### Example 3.6.

The manifold $`M=SO(13)/(SL(2)\times SL(2)\times SL(2)\times SL(2)\times SL(2)\times SL(2))`$, has $`\mathrm{\Gamma }_𝔨=_2\times _2`$. Therefore, it satisfies the point (III) of the Corollary.

## 4 Quantization

###### Theorem 4.1.

Let $`G`$ be a semisimple Lie group, $`K`$ a reductive subgroup of maximal rank. Then, for any admissible Poisson bracket on $`M=G/K`$ there exists a corresponding $`\mathrm{U}_h(𝔤)`$-equivariant quantization.

###### Proof.

As follows from Corollary 3.5, one can suppose that $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)`$ is cyclic. The case when $`\mathrm{Tor}(\mathrm{\Gamma }_𝔨)=0`$, i.e. $`M`$ is a semisimple coadjoint orbit, is considered in , while the case when $`\mathrm{\Gamma }_𝔨`$ is cyclic itself is considered in . The complete proof of the theorem obtains combining methods of the both papers. ∎

###### Remark 4.2.

As was noticed in Introduction, for any invariant Poisson bracket on $`M`$ under consideration there exists a corresponding $`G`$-equivariant quantization.

Question. Describe compatible pairs consisting of invariant and admissible Poisson brackets on $`M`$. Does there exists for any compatible pair a two parameter $`\mathrm{U}_h(𝔤)`$-equivariant quantization? This problem was solved for semisimple coadjoint orbits in .

Acknowledgments.  We thank A. Elashvili, S. Shnider, and A. Stolin for very helpful discussions.


format: text

text:
# Microscopic Structure of a Vortex Line in a Dilute Superfluid Fermi Gas

\[

## Abstract

The microscopic properties of a single vortex in a dilute superfluid Fermi gas at zero temperature are examined within the framework of self-consistent Bogoliubov-de Gennes theory. Using only physical parameters as input, we study the pair potential, the density, the energy, and the current distribution. Comparison of the numerical results with analytical expressions clearly indicates that the energy of the vortex is governed by the zero-temperature BCS coherence length.

\]

The trapping and cooling of dilute Fermi gases is an increasingly active area of research within the field of ultracold atomic gases. Several experimental groups now trap and cool alkali atoms with Fermi statistics reaching temperatures as low as $`0.2T_F`$ with $`T_F`$ being the Fermi temperature of the system . These gases are appealing to study due to the large experimental control of their properties and the microscopically well-understood two-body effective interaction between the atoms. One of the main goals of the present experimental effort is to observe a phase transition to a superfluid state predicted to occur below a certain critical temperature $`T_c`$ . A fascinating prospect of such a superfluid is the formation of quantized vortices. The combination of theoretical and experimental studies of vortices in Bose-Einstein condensates (BECs) have produced a number of beautiful results in recent years . For fermions, the nature of vortices in different systems such as type-II superconductors, superfluid 3He, and neutron stars is a classic problem with a vast literature .

One problem of fundamental importance is to calculate the energy of a vortex for $`T=0`$. This energy is defined as the difference between the energy of the superfluid state with the vortex present and without. Dividing this value by the angular momentum per particle in the vortex state gives the frequency at which the vortex becomes the thermodynamic ground state of a rotating system. It is rather surprising that despite the large amount of work concerned with the structure of a vortex for a fermionic system, there is no clear result regarding this specific problem. This is in contrast with the case of a dilute bosonic superfluid, where the Gross-Pitaevskii equation allows an analytical calculation of the vortex energy for $`T=0`$ if quantum fluctuations are neglected . The equivalent theory relevant for fermions, Ginzburg-Landau theory, is unfortunately only valid for $`|TT_c|/T_c1`$ making an analytical calculation of the energy for $`T=0`$ more complicated. In this letter we present the first ab-initio calculation of the vortex energy based on microscopic theory.

We consider in the following a two-component gas of neutral fermionic atoms with mass $`m`$ in a cylinder of radius $`R`$. We take the number of particles in each spin state $`N_\sigma `$ to be the same as this is the optimum situation for superfluidity . In the dilute limit, the interaction between the atoms in the two spin states $`\sigma =,`$ can be well described by the contact potential $`g\delta (𝐫)`$ where $`g=4\pi \mathrm{}^2a/m`$ and $`a<0`$ is the $`s`$-wave scattering length describing low energy collisions between $``$ atoms and $``$ atoms. In the zero-temperature limit that we are treating here there are no intra-component collisions . Recently, two papers have calculated the $`T=0`$ vortex energy for a Fermi superfluid under these assumptions. Using phenomenological models, the energy of a unit circulation vortex was estimated in Ref.  to be

$$_v\frac{\pi \mathrm{}^2n_\sigma }{2m}\mathrm{ln}\left(D\frac{R}{\xi _{\mathrm{BCS}}}\right),$$
(1)
where $`n_\sigma =\stackrel{~}{k}_F^3/6\pi ^2`$ is the density in each of the two components, and $`\xi _{\mathrm{BCS}}=\mathrm{}^2\stackrel{~}{k}_F/\pi m\mathrm{\Delta }_0`$ is the BCS coherence length with $`\mathrm{\Delta }_0`$ the bulk value of the superfluid gap; BCS theory predicts $`\mathrm{\Delta }_0=8\mathrm{e}^2\stackrel{~}{\mu }\mathrm{exp}(\pi /2\stackrel{~}{k}_F|a|)`$ . The effective Fermi momentum, $`\stackrel{~}{k}_F`$, is defined as $`\mathrm{}^2\stackrel{~}{k}_F^2/2m=\mu gn_\sigma \stackrel{~}{\mu }`$, where $`\mu `$ is the chemical potential, such that it includes the effect of the Hartree mean-field. It was argued in Ref.  that a microscopic calculation for $`T=0`$ would yield $`D`$ to be a *constant* $`𝒪(1)`$ independent of $`\stackrel{~}{k}_F`$ and $`|a|`$ since the characteristic length-scale of a vortex must be expected to be $`𝒪(\xi _{\mathrm{BCS}})`$. The value of $`D`$ depends on the phenomenological model used: If the vortex is modeled as a cylinder of radius $`\xi _{\mathrm{BCS}}`$ containing a normal stationary fluid, surrounded by a rotating superfluid one obtains $`D=1.36`$. We refer to this simple model as the cylinder model. If Ginzburg-Landau theory is applied, we get $`D=1.65`$.

This conclusion was however disputed in the work of Ref. . Here it was argued that the characteristic length scale of the vortex is much smaller than $`\xi _{\mathrm{BCS}}`$ and the energy correspondingly higher. This is because the structure of the vortex core is determined by the lowest lying vortex states. These states are formed out of excitations around the Fermi level with typical wavelengths $`\stackrel{~}{k}_F^1`$, and following the conclusions based on the analytical and numerical solutions of the Bogoliubov-de Gennes (BdG) equations  it was argued that the important length scale of the core region is $`\xi _1=4/\pi \stackrel{~}{k}_F^2|a|\xi _{\mathrm{BCS}}`$ in the dilute regime . Using $`\xi _1`$ as the size of the vortex core, a calculation identical to the one given in Ref.  leads to a vortex energy given by Eq. (1) but with $`D\xi _{\mathrm{BCS}}/\xi _11`$ in the dilute regime. Thus, the energy was predicted to be significantly higher than what was estimated in Ref. . Note that $`D`$ is now not a constant but depends on $`\stackrel{~}{k}_F`$ and $`|a|`$.

It is presently not clear which of the two quite different predictions is correct and thus what the energy of the vortex actually is. In order to settle this question, we now present a *microscopic* calculation of the vortex energy using the assumptions given above. The BdG equations describing the superfluid state read

$`\left[\begin{array}{cc}H_0(𝐫)& \mathrm{\Delta }(𝐫)\\ \mathrm{\Delta }^{}(𝐫)& H_0(𝐫)\end{array}\right]\left[\begin{array}{c}u_\eta (𝐫)\\ v_\eta (𝐫)\end{array}\right]=E_\eta \left[\begin{array}{c}u_\eta (𝐫)\\ v_\eta (𝐫)\end{array}\right],`$ (8)
with $`H_0(𝐫)=\mathrm{}^2^2/2m\mu +U(𝐫)`$, and where $`U(𝐫)=gn_\sigma (𝐫)=g_\eta |v_\eta (𝐫)|^2`$ is the Hartree field and $`\mathrm{\Delta }(𝐫)=\stackrel{~}{g}_\eta u_\eta (𝐫)v_\eta ^{}(𝐫)`$ is the gap function. To avoid having to introduce an arbitrary high energy cut-off in the theory, we have used a zero range pseudo-potential scheme giving rise to a regularized coupling constant $`\stackrel{~}{g}`$, when calculating $`\mathrm{\Delta }(𝐫)`$ . The result is a well-defined theory using only physical parameters as input. Further details of this and the numerical techniques used to solve these equations will be given elsewhere. Once the self-consistent solution is obtained for a given coupling strength and chemical potential, the energy is given by $`E=\widehat{H}\mu \widehat{N}`$, where $`\widehat{H}`$ is the Hamiltonian of the system , and $`\widehat{N}`$ is the number operator. Defining $`E(𝐫)2_\eta |v_\eta (𝐫)|^2E_\eta `$, $`E`$ can be calculated using

$$E=d^3r\left\{E(𝐫)+\frac{1}{\stackrel{~}{g}}\left[|U(𝐫)|^2+|\mathrm{\Delta }(𝐫)|^2\right]\right\}.$$
(9)
We now solve the BdG equations for the gas in a cylinder of height $`L_z`$ and radius $`R\xi _{\mathrm{BCS}}`$. Since $`\xi _{\mathrm{BCS}}`$ is a decreasing function of $`\stackrel{~}{k}_F|a|`$ we use two different cylinder sizes. For $`\stackrel{~}{k}_F|a|0.43`$ $`(\stackrel{~}{k}_F|a|0.43)`$ we take $`R=44.1\mu \mathrm{m}`$ $`(R=12.6\mu \mathrm{m})`$ and $`L_z=12.5\mu \mathrm{m}`$ $`(L_z=5\mu \mathrm{m})`$. We find the lowest energy superfluid state of the system by setting $`\mathrm{\Delta }(𝐫)=\mathrm{\Delta }(\rho ,z)`$ where $`\rho `$ is the perpendicular distance from the axis of symmetry of the cylinder, and $`z`$ is the axial coordinate. For a vortex state, we assume the form $`\mathrm{\Delta }(𝐫)=\mathrm{exp}(i\varphi )\mathrm{\Delta }(\rho ,z)`$ where $`\varphi `$ is the azimuthal angle around the cylinder axis. This corresponds to a vortex line with unit circulation along the cylinder axis. In both cases $`U(𝐫)=U(\rho ,z)`$. In the vortex-free case, the cylindrical symmetry dictates that Cooper pairs form between particles with angular momentum $`\mathrm{}\nu `$ and $`\mathrm{}\nu `$ along the cylinder axis whereas in the vortex case pair constituents have angular momentum $`\mathrm{}\nu `$ and $`\mathrm{}(\nu +1)`$ . Once the two solutions with and without the vortex are obtained, the energy per unit length of the vortex line can be determined as $`_v=(E_vE_0)/L_z`$ where $`E_v`$ denotes the energy of the vortex state and $`E_0`$ the energy of the state without a vortex, both obtained from Eq. (9). Throughout the succeeding analysis we consider a gas of 6Li atoms with $`a=955`$ $`a_0`$. The density in each hyperfine state is chosen to be a few times $`10^{13}\mathrm{cm}^3`$. These parameters are appropriate for on-going experiments . The numerical solution of the BdG equations is obtained within a Bessel function discrete variable representation in $`\rho `$ , and periodic boundary conditions along the vortex axis.

In Fig. 1, we plot the total energy density $`E/V`$ of the superfluid gas, where $`V`$ is the volume of the cylinder, and $`E`$ is given by Eq. (9), as a function of the effective interaction strength $`\stackrel{~}{k}_F|a|`$. The $`\times `$’s and $``$’s are obtained from a self-consistent numerical solution of the BdG equations for two different values of $`R`$, whereas the line is the ground state energy per unit volume of a bulk Fermi superfluid,

$$\frac{E_{bulk}}{V}=\frac{6}{5}n_\sigma \frac{\mathrm{}^2\stackrel{~}{k}_F^2}{2m}\mu 2n_\sigma +gn_\sigma ^2\frac{N(0)\mathrm{\Delta }_0^2}{2}.$$
(10)
Here $`N(0)=m\stackrel{~}{k}_F/2\pi ^2\mathrm{}^2`$ is the density of states at the Fermi level. This expression is obtained by integrating analytically Eq. (9) for a homogeneous gas where the $`u(𝐫)`$’s and $`v(𝐫)`$’s are simple plane wave states. The first three terms in Eq. (10) give to the energy of a homogeneous gas in the normal phase within the Hartree-Fock approximation and the last term is the condensation energy. We see that there is good agreement between our numerical results and the analytical formula. The slight discrepancy is due to boundary effects at the edge of the cylinder, where the density of particles vanishes. Note that we are at the limit of the weak coupling regime $`\stackrel{~}{k}_F|a|1`$, appropriate for dilute gases. For the purposes of comparison with analytical results, however, it is important to calculate properties for the widest possible range of $`\xi _{\mathrm{BCS}}`$, subject to the condition $`\xi _{\mathrm{BCS}}R`$ which ensures that the gap function can heal to its bulk value before becoming suppressed at the cylinder surface.

In Fig. 2, we plot the numerically calculated energy of the vortex $`_v`$ for varying $`\stackrel{~}{k}_F|a|`$. To compare with the analytical predictions, we parameterize $`_v`$ by the variable $`D`$ appearing in Eq. (1). The dashed line corresponds to the prediction $`D=1.36`$  and the solid line to $`D=\xi _{\mathrm{BCS}}/\xi _1`$ . We see that the two predictions for $`D`$ have a completely different dependence on $`\stackrel{~}{k}_F|a|`$. The important conclusion is that the numerical results confirm $`D𝒪(1)`$ being a *constant* independent of $`\stackrel{~}{k}_F|a|`$ in agreement with Ref. . On the other hand, the prediction $`D=\xi _{\mathrm{BCS}}/\xi _1`$ yields a qualitatively incorrect result. We note that the kink in $`_v`$ is due to $`R/\xi _{\mathrm{BCS}}`$ being different for the two cylinder sizes, whereas the spread in $`D`$ at $`\stackrel{~}{k}_F|a|=0.43`$ is indicative of our numerical accuracy. The numerical value $`D2.5`$ is higher than the prediction of phenomenological models in Ref. . This is as expected since these models only can yield the correct order of magnitude of the constant inside the logarithm. Thus, the length scale determining the energy of the vortex is $`\xi _{\mathrm{BCS}}`$ and not $`\xi _1`$.

To examine this in more detail, we plot in Fig. 3 the numerically calculated profile of a vortex for two representative values of $`\stackrel{~}{k}_F|a|`$. Close to the vortex core, only the lowest-energy (bound) states contribute to the order parameter; these give rise to the observed Friedel oscillations, which have a wavelength on the order of $`\stackrel{~}{k}_F^1`$. We see that the length scale defined as $`\xi _1=lim_{\rho 0}[\mathrm{\Delta }(\rho ,z)/\rho \mathrm{\Delta }_{\mathrm{}}]^1`$ giving the slope of $`\mathrm{\Delta }(𝐫)`$ at the vortex core actually is much smaller than $`\xi _{\mathrm{BCS}}`$ as predicted in Ref. . Here, $`\mathrm{\Delta }_{\mathrm{}}`$ is the value of $`|\mathrm{\Delta }(𝐫)|`$ far away from the vortex core with $`\mathrm{\Delta }_{\mathrm{}}\mathrm{\Delta }_0`$ as expected. However, as the distance $`\rho `$ from the vortex core increases, the slope decreases and $`\mathrm{\Delta }(𝐫)`$ reaches the value $`\mathrm{\Delta }_{\mathrm{}}`$ on a length scale $`\xi _{\mathrm{BCS}}`$ and not $`\xi _1`$. To quantify this, we use the cylinder model of the vortex with a vortex radius $`\xi _2=x\xi _{\mathrm{BCS}}`$ to calculate $`_v`$. This yields Eq. (1) but now with $`D=(1.36)^{x^2}/x`$. The equation $`D=2.5`$ then gives $`x=0.42`$. Thus $`\xi _2=0.42\xi _{\mathrm{BCS}}`$ is the length scale determining the energy of the vortex. Again, it should be emphasized that $`x0.42`$ is a *constant* over the large range of $`\xi _{\mathrm{BCS}}`$ used in the calculations thereby verifying that indeed $`\xi _{\mathrm{BCS}}`$ determines the length-scale relevant for the energy as discussed in . The cylinder model of $`\mathrm{\Delta }(\rho )`$ with the correct radius $`\xi _2=0.42\xi _{\mathrm{BCS}}`$ is plotted in Fig. 3.

To examine the superfluid flow associated with the vortex giving rise to the angular momentum, we plot the current density $`𝐣_s(\rho )`$ given by

$$𝐣_s(\rho )=\frac{2\mathrm{}}{mi}\underset{\eta }{}v_\eta ^{}(𝐫)\rho ^1_\varphi v_\eta (𝐫)𝐞_\varphi $$
(11)
in Fig. 4 for $`\stackrel{~}{k}_F|a|=0.59`$. Because the normal component carries no current by construction, the total current density may be written as $`𝐣_s(\rho )=2n_s𝐯_s`$ with the superfluid velocity $`𝐯_s=𝐞_\varphi \mathrm{}/2m\rho `$ and the superfluid (or atom-pair) density $`n_s`$. We plot $`n_s(\rho )`$ defined in this way in Fig. 4. Note that unlike dilute interacting Bose gases at zero temperature, the superfluid density in a dilute Fermi superfluid does not have the same behavior as the order parameter $`\mathrm{\Delta }(𝐫)`$. As expected, $`n_s(\rho )n_\sigma `$ far away from the vortex core and $`n_s(\rho )0`$ for $`\rho 0`$. The dotted line in Fig. 4 gives the analytical result $`j_\varphi (\rho )=\mathrm{}\stackrel{~}{k}_F^3/(6\pi ^2m\rho )`$, which agrees well with the numerics away from the vortex core. As a self-consistency check on the numerics, the angular momentum per particle along $`z`$ is found to be exactly $`\mathrm{}/2`$, corresponding to one unit of angular momentum $`\mathrm{}`$ per Cooper pair as expected.

As shown in Fig. 4, the presence of a vortex in a Fermi system does not lead to any significant change in the particle density $`n_\sigma (\rho )`$ , in contrast with a dilute Bose gas, where the density is minuscule on the vortex line . Direct observation of the vortex core (now commonplace for BECs) is therefore not likely. The quantized currents, and therefore the presence of superfluidity, can be readily detected using at least three approaches, however. One of these is the collective mode spectrum. When no vortex is present, excitations carrying equal and opposite angular momentum along the $`z`$-axis are degenerate in energy. The vortex currents lift this degeneracy since the rotational symmetry is removed. The resultant splitting of the surface modes is proportional to the angular momentum of the gas . This technique has been used to infer the presence of a vortex in a trapped BEC . A second approach was demonstrated in a recent experiment where the precession rate of the scissors oscillation mode was used to measure the quantized angular momentum per particle with great accuracy . A third method is spatially selective Bragg scattering; the superfluid currents modify the Bragg momentum conservation conditions, giving rise to a strongly anisotropic outcoupled atomic beam .

In conclusion, we have accurately determined the energy of an isolated vortex line in a dilute superfluid Fermi gas at zero temperature. The results clearly indicate that the BCS coherence length sets the scale for the vortex energy and therefore the critical frequency for vortex stability.

We acknowledge valuable discussions with Ø. Elgarøy, C. J. Pethick, and B. I. Schneider.


text:
# Parameter identification using the Hilbert transform

## introduction

In the theory of control, it is most common for physical systems to be mathematically modelled using coupled linear systems of ordinary differential equationsLevine (1996). When these equations are transformed using integral transforms, such as those of Laplace or Fourier, then the physical systems are modelled using finite rational polynomials in an auxiliary variable, $`s=j\omega `$:

$$H(s)=\frac{\text{output}(s)}{\text{input(s)}}=\frac{P(s)}{Q(s)}.$$
(1)
The zeros of the polynomial, $`Q(s)`$, are called “poles” and correspond to responses that have finite output for zero input. These are called ”modes.” It is very common for one mode to dominate the response of the whole system. It is also common for this mode to be of a damped oscillatory type, corresponding to a single pair of complex conjugate poles. This can occur whenever the potential energy function of the system has a local minimum Goldstein (1950) In this case, we can approximate a large complicated system, with many poles and zeros, by a simple second-order system with a single pair of complex conjugate poles. This is called the second-order approximation. Many mechanical or electrical systems can be realistically modelled using the second-order approximation. We can write:

$$H(s)\frac{a_2s^2+a_1s+a_0}{s^2+2\alpha s+\omega _{0}^{}{}_{}{}^{2}}.$$
(2)
If we wish to model the behaviour of a real physical system, using an approximate second-order model, then it is necessary for us to estimate the position of the pole pair. This could be done in the frequency domain, by exciting the system with a sinusoidal source and then measuring the magnitude and phase of the response at different frequencies, but this is often not practical. There are situations when the only practical sources are step functions, $`1(t)`$ or impulses $`\delta (t)`$. We can excite the system with steps or impulses and then sample the response in the time domain. The impulse response of a second-order system will generally be of the form:

$$y(t)=Ae^{\alpha t}\mathrm{cos}\left(\omega _dt\right)+Be^{\alpha t}\mathrm{sin}\left(\omega _dt\right)$$
(3)
where $`\omega _{d}^{}{}_{}{}^{2}=\omega _{0}^{}{}_{}{}^{2}\alpha ^2`$.

The problem of plant identification then becomes equivalent to asking: How do we estimate the position of the pair of complex-conjugate poles if the only data at our disposal is a set of time-domain samples of the response of the system to steps or impulses.

As a possible illustration, we could imagine that we strike a bell with a hammer and then record the sound as it gradually decays. We want to estimate the damped frequency of oscillation, $`\omega _d`$, and the damping coefficient, $`\alpha `$, using only the data from our sound recording.

If we knew something about the distribution of the errors of measurement then we could apply the method of maximum likelihood to estimate the parameters , $`\omega _d`$ and $`\alpha `$. If the errors were known to be the result of a very large number of uncorrelated random effects then we could apply the Central Limit Theorem and we could assume that the errors had a Gaussian distribution. The problem of plant identification would reduce to a non-linear least-squares estimation problem Press et al. (1995). The difficulty with this approach is that the resulting equations would be non-linear and would have to be solved iteratively, using a numerical method such as gradient descent. A further weakness of this approach is that it would be an exact solution to an approximation of the real problem. It would be far more reasonable to have a quick but approximate solution to the approximate, second-order, problem. This would tell us most of what we need to know without having to waste a lot of effort.

In this paper, we present a fast, but approximate, algorithm for the estimation of the position of a complex conjugate pair of poles on the $`s`$ plane. We use a discrete approximation to the Hilbert Transform , to reconstruct the complex analytic signal form the sampled real time signal. The analytic signal has a complex-exponential form. We apply very simple statistical techniques to the analytic signal in order to obtain the required parameters. This approach is an alternative to the more conventional, non-linear least squares or Autoregressive, AR, approaches to the problem.

The problem of parameter identification, for a freely vibrating system has been studied by Feldman Feldman (1997) who used the Hilbert Transform to provide information about instantaneous amplitude and phase of a signal. The method that we present here is more simple, and limited, than the approach used by Feldman.

## The analytic signal of the second order response

The response described in Equation 3 is equivalent to

$$y(t)=Ce^{\alpha t}\mathrm{cos}\left(\omega _dt\varphi \right)$$
(4)
where $`C=\sqrt{A^2+B^2}`$, $`\mathrm{cos}(\varphi )=A/\sqrt{A^2+B^2}`$ and $`\mathrm{sin}(\varphi )=B/\sqrt{A^2+B^2}`$. This type of function will apply whenever the input to the system is zero. If the input is a finite sum of step and impulse function then the input will be zero for most of the time. There will be abrupt changes in $`C`$ and $`\varphi `$ but the parameters, $`\alpha `$ and $`\omega _d`$ will be constant as long as the structure of the plant is maintained.

The immediate aim is to reconstruct the analytic signal. The Hilbert transform Bracewell (1965); Hahn (1996); Ersoy (1997) is a standard technique for achieving this. The Hilbert transform of $`u(t)`$ is defined as:

$$v(s)=\frac{1}{\pi }_{\mathrm{}}^+\mathrm{}\frac{u(t)}{st}𝑑t.$$
(5)
It has the important property that it defines the relationship between the real and imaginary parts of a complex analytic function. If we have an analytic function:

$$\mathrm{\Phi }(t)=u(t)+jv(t)$$
(6)
with real and imaginary parts $`u(t)`$ and $`v(t)`$ then the relation for $`v(s)`$ is given by Equation 5. Bedrosian’s theorem tells us that the Hilbert transform of $`a(t)\mathrm{cos}(\omega t)`$ is $`a(t)\mathrm{sin}(\omega t)`$. If we apply Bedrosian’s Hahn (1996)theorem and the shifting property to Equation 4 then we find that the analytic signal is:

$`y_2(t)`$ $`=`$ $`Ce^{\alpha t}\mathrm{cos}\left(\omega _dt\varphi \right)+jCe^{\alpha t}\mathrm{sin}\left(\omega _dt\varphi \right)`$ (7)
$`=`$ $`Ce^{j\varphi }e^{\left(\alpha +j\omega _d\right)t}.`$ (8)
This analytic signal is a pure exponential function and can essentially be “unwrapped” using the $`\mathrm{log}()`$ function. We can write:

$$\mathrm{log}\left(y_2(t)\right)=\mathrm{log}(C)j\varphi +\left(\alpha +j\omega _d\right)t.$$
(9)
If we sample this analytic signal at intervals of $`T_s`$ then we can numerically calculate the slope of Equation 9 to get:

$$\alpha +j\omega _d=\left(\frac{\mathrm{log}\left(y_2\left(t+T_s\right)\right)\mathrm{log}\left(y_2\left(t\right)\right)}{T_s}\right).$$
(10)
This allows us to directly estimate the parameters, $`\alpha `$ and $`\omega _d`$. We note that the use of the difference operation has removed all reference to $`C`$ and to $`\varphi `$. This means that the method is not sensitive to the initial conditions that apply immediately after the shocks that occur then the impulses and steps are fed into the system. Our only requirement is that the input to the system is zero for most of the time.

The Matlab code required to implement this algorithm is very short and simple:

```
% reconstruct the analytic signal, y2
y2 = hilbert( y );
% take the natural logarithm, log(y2)
L = log(y2) ;
% unwrap the phase of the log(y2)
L = real(L) + j*unwrap(imag(L));
% estimate the differences of the log of y2
D = diff(L) ;

```

## A simple statistical technique

Equation 10 suggests that we should be able to precisely estimate the required parameters. There are a few practical problems with the direct application of Equation 10:

* The numerical calculation of the Hilbert transform relies on the Fast Fourier Transform and there are limitations imposed by the finite number of samples. These includes the “Gibbs effect,” due to the finite length of the data set.
* Real samples from a physical process will be subject to noise and errors of measurement. The differencing operation tends to magnify the effect of noise.
* There would be a number of outliers caused by the “shocks” of the steps or impulses. Some measurements will not be reliable.

The authors have found that the median is a very robust measure of location and is less subject to the influence of the outlying values than the arithmetic mean. The Matlab code for this is very simple:

```
% calculate the real and imaginary parts of the differences
d_mag   = real(D) ;
d_phase = imag(D) ;
% calculate the median rates of change
mid_re_slope = median(d_mag) ;
mid_im_slope = median(d_phase) ;

```

## Some results

A second order system was simulated, using known parameters, and the parameters were then estimated using the new algorithm. The reconstructed analytic signal is shown in Figure 1.

This same data can be represented in three dimensions. This is shown in Figure 2.

The logarithmic slopes were estimated and the parameters were calculated.

The quality factor of the simulated system was $`Q=10`$ and so the theoretical value of the damping constant was $`\alpha =0.05`$. The estimated value was $`\alpha _{\text{est}}=0.0472`$. The characteristic frequency of the simulated plant was $`\omega _d=1.0`$ and the estimate was $`\omega _{\text{est}}=0.9948`$.

There was a $`5.5\%`$ error in the estimate of the damping constant, $`\alpha `$ and there was a $`0.4\%`$ error in the estimate of the damped natural angular frequency of oscillation, $`\omega _d`$.

## Summary and limitations

The method does work as an approximation and could estimate $`\alpha `$ to within $`5.5\%`$ and $`\omega _d`$ to within $`0.4\%`$. The error in the estimate of the damping constant, $`\alpha `$ and there was a $`0.4\%`$ error in the estimate of the damped natural angular frequency of oscillation, $`\omega _d`$. There are some notable problems with the method:

* It is sensitive to relative time scale of the constants in question and the width of the sampling window.
* It suffers from boundary effects due to the Gibbs’ phenomenon.
* The method is also numerically unstable for large data sets with large values of time, $`t`$, since $`exp(\alpha t)`$ can cause Matlab to underflow.

Our simulations suggest that this is a fast but approximate technique for the estimation of the position of a single pair of complex conjugate poles but it is sensitive to a number of factors which may limit its practical application.


text:
# Scaling Properties of 1D Anderson Model with Correlated Diagonal Disorder

## I Introduction

It has been known for almost forty years that in standard one-dimensional disordered models all states are localized for any strength of disorder, and that there is no localization transition in such systems.gang ; LGP  Formally this is expressed by the statement that the Lyapunov exponent (LE), $`\gamma `$, defined as

$$\gamma =\underset{L\mathrm{}}{lim}\frac{1}{2L}\mathrm{ln}T,$$
(1)
where $`T`$ is the transmission coefficient through the system of length $`L`$, is always positive.LGP  While it may seem that there is nothing more to be said about one-dimensional models, their localization properties have recently attracted a great deal of attention. This renewed interest has been concentrated in two areas. The first area includes studies of unusual models that demonstrate the presence of extended states. It was noted first in Ref. dimer,  that introducing correlations in the statistical properties of random site energies in the Anderson model, one can create extended states in the one-dimensional Anderson model at certain values of energy. Later on, the authors of Ref. Moura,  showed that regions of extended states exist in 1D systems with long range correlations of the random potential. The authors of Ref. Izrailev,  demonstrated that vanishing LE can be obtained in an arbitrarily chosen spectral region by selecting a special form of the correlation function of the random potential. It should be noted, however, that this result was obtained in Ref. Izrailev,  only in the second order of the weak disorder expansion. Taking into account the next terms of the expansion makes LE positive again.fourth order pt

The second area of research interest in the field of one-dimensional localization is focused upon scaling and statistical properties of finite size systems. The single parameter scaling put forward in Ref. gang,  is a cornerstone of the current approach to the Anderson localization, but its interpretation in the presence of non-self-averaging fluctuations of conductivity was the subject of long debate. Eventually it was understood that the SPS hypothesis in Anderson localization means that the entire distribution function of conductivity (or transmission), and not individual moments, must be parameterized by a single parameter.Anderson ; Shapiro  This distribution in one-dimensional systems has been intensively studiedKramerReview , and it was established that the bulk of the distribution for sufficiently long systems has a log-normal form. It was also understood starting with Ref. Anderson,  that under certain circumstances the two parameters of the distribution, the average value of the finite size LE, $`\stackrel{~}{\gamma }=\mathrm{ln}T/(2L)`$, which coincides with its limiting value, $`\gamma =\stackrel{~}{\gamma }`$, and the variance, $`\sigma ^2=var(\stackrel{~}{\gamma }),`$ are related to each other in a universal way

$$\tau =\frac{\sigma ^2L}{\gamma }=1.$$
(2)
This relation reduces two parameters of the distribution to only one, and provides, therefore, a justification and interpretation for SPS. Conditions, under which Eq. (2) holds, however, have been established only recently, when the authors of Ref. PRLclassics,  showed that the validity of SPS is controlled by a new macroscopic length, $`l_s`$, defined in terms of the integral density of states, $`N(E)`$:

$$l_s^1=\mathrm{sin}(\pi N(E)).$$
(3)
where $`E`$ is energy. It was shown that SPS holds when $`\kappa =(\gamma l_s)^1>1`$, and fails in the opposite case. Recently, the authors of the present paper showed that parameter $`\kappa `$ plays a more important role than just establishing the criterion for SPS. It was found that in the non-SPS region the function $`\tau `$ defined in Eq. (2), and the similar dimensionless combination of the third moment with the length $`L`$ and the Lyapunov exponent $`\gamma `$ are functions of this single parameter.nonSPSlast  In the recent paper Ref. Titov, , higher moments of the distribution of the Lyapunov exponent were calculated analytically for a quantum particle in a potential described by the Gaussian delta-correlated random function. It was found that all moments of the distribution for this model can be expressed as functions of the localization length and the parameter $`E/D^{2/3}`$ , where $`D`$ characterizes the strength of the potential. For the Gaussian white noise potential the parameter $`\kappa `$ is also a function of $`E/D^{2/3},`$PRLclassics  thus in the particular case of the Gaussian white noise potential these two parameters are equivalent.

In the present paper, we combine the two discussed areas and focus upon correlation induced changes in the probability distribution function of the Lyapunov exponent. In order to study this problem, we consider a tight-binding model with the random potential described by a zero-mean dichotomic process. This process is characterized by an exponential correlation function with a correlation radius, which is a free parameter of the model. There is no localization-delocalization transition in this model, and we are interested in effects of the correlation radius on scaling and statistical properties of conductance in the regime of strongly localized states. This model allows for an approximate analytical treatment in the spectral region, where the correlation radius becomes a dominant length. Combining analytical and numerical calculations we show that the presence of this additional length significantly changes the transition between SPS and non-SPS regions, and that in the non-SPS region it results in a new scaling behavior.

## II Quasiclassical approximation for Lyapunov exponent

We consider the tight binding model, described by the equation of motion

$$\psi _{n+1}+\psi _{n1}+(V\zeta _nE)\psi _n=0,$$
(4)
where $`V`$ describes the strength of the random potential. The potential randomly takes one of two values $`\pm V`$ depending upon whether the dichotomic random variable $`\zeta _n`$ is equal to $`1`$ or $`1`$. The assumed value remains constant within regions of random lengths (segments), which are distributed according to the discrete analog of the Poisson distribution. The average length of such regions is

$$l_c=\mathrm{coth}\left(\frac{1}{2r_c}\right),$$
(5)
where $`r_c`$ is the correlation radius defined by

$$\zeta _n\zeta _m=e^{|nm|/r_c}.$$
(6)
The finite size LE is defined through the norm of the transfer matrix:

$$\stackrel{~}{\gamma }(E)=\frac{1}{L}\mathrm{log}T_N\mathrm{}T_1,$$
(7)
where $`T_n`$ are transfer matrices which, in our case, are equal to either $`T_+`$, or $`T_{}`$ depending upon the value of $`\zeta _n`$, where $`T_\pm `$ are

$$T_\pm =\left(\begin{array}{cc}E\pm V& 1\\ 1& 0\end{array}\right).$$
(8)
To consider the spectrum of this system let us note that an addition of a constant potential $`\pm U`$ shifts a conduction band, $`2E2`$, of the homogeneous tight-binding system with no on-site potential by the value of the potential: $`2UE2U`$. Correspondingly, the spectrum of the random system under consideration can be divided into three regions with the qualitatively different spectral and transport properties. One region, $`0<|E|<2V`$, corresponds to the common part of the conduction bands of homogeneous systems with the site potentials equal to either $`+V`$ or $`V`$. In the second spectral region, consisting of energies $`2V<E<2+V`$ and $`2V<E<2+V`$, extended states exist in a homogeneous system with the site potentials equal to $`V`$ or $`+V`$, respectively. Finally the third region, $`2+V<|E|`$ corresponds to the energies, where extended states would not arise in homogeneous systems with either of the two possible values of the potential. Obviously, no states can arise in this region also in the system with the random distribution of the potential.

While all states in this model are localized, different spectral regions still have different transport properties. Transport in the first region is characterized by multiple scattering of propagating modes, and the localization in this region results from the interference of the scattered waves. Transport in the second region can be described as tunneling through states, arising within the segments of the system, where the potential produces potential wells between barriers formed by segments with the opposite value of the potential. If the average length of the “barrier” regions is larger than the penetration length under the barrier, the transport is mostly determined by the under-barrier tunneling, and the interference effects related to phase changes of the wave functions inside potential wells become less important. We are mostly interested in this region, where correlations are expected to result in the most non-trivial effects. One could expect that the transport in this region can be described within the quasi-classical approximation, and we will show below that this is indeed the case.

In the context of our problem, this approximation means neglecting commutators between transfer matrices at different sites. Introducing LE in a homogeneous system with the potential $`U`$ according to

$$\gamma _0(E;U)=\mathrm{Re}\left[\lambda _0(E;U)\right],$$
(9)
where

$$\lambda _0=\mathrm{cosh}^1\left(\frac{EU}{2}\right),$$
(10)
we can write a quasiclassical expression for LE for Eq. (4) in the following form

$$\stackrel{~}{\gamma }(E)\frac{1}{L}\underset{n=1}{\overset{N}{}}\gamma _0(E;V\zeta _n).$$
(11)
Using the fact that the $`\zeta _n=\pm 1`$, Eq. (11) can be rewritten in the form

$$\stackrel{~}{\gamma }=\gamma _a+\frac{\delta }{L}\underset{n=1}{\overset{N}{}}\zeta _n,$$
(12)
where

$`\gamma _a`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\gamma _0(E;V)+\gamma _0(E;V)\right],`$
$`\delta `$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\gamma _0(E;V)\gamma _0(E;V)\right].`$ (13)
The averaged finite-length LE, $`\gamma =\stackrel{~}{\gamma }`$, is determined by $`\gamma _a`$, while its variance $`\sigma ^2`$ depends upon $`\delta `$:

$$\sigma ^2=\gamma ^2\gamma ^2=\frac{\delta ^2}{L}l_c.$$

In the first spectral region this approximation yields vanishing LE, which is quite understandable, because interference effects responsible for the localization in this region are completely neglected.

Eq. (11) can be formally obtained in the following way. The relation between LE and the eigenvalues of the transfer matrix allows us to write

$$2\mathrm{cosh}(\stackrel{~}{\gamma }L)=\mathrm{Tr}\left(T_{N,1}\right).$$
(14)
where the transfer matrix $`T_{N,1}`$ consists of sequences of $`T_\pm `$

$$T_{N,1}=T_N\mathrm{}T_1=\mathrm{}T_+^{n_3}T_{}^{n_2}T_+^{n_1},$$
(15)
Powers $`n_k`$ of transfer matrices $`T_\pm `$ represent the lengths of the regions where the potential remains constant. The matrices $`T_\pm `$ can be diagonalized $`T_\pm ^n=R_\pm S_\pm ^nR_\pm ^1`$, where

$$S_\pm ^n=\left(\begin{array}{cc}e^{\lambda _0(\pm V)n}& 0\\ 0& e^{\lambda _0(\pm V)n}\end{array}\right)$$
(16)
Thus, introducing $`\alpha _{12}=R_+^1R_{}`$ and $`\alpha _{21}=R_{}^1R_+`$, and noting that for long enough samples the effect of matrices $`R_\pm `$ appearing at the end of the structure on LE is negligible, we obtain that LE can be found from Eq. (14) with the transfer matrix replaced by

$$\stackrel{~}{T}_{N,1}=\mathrm{}\alpha _{21}S_+^{n_2}\alpha _{12}S_{}^{n_1}.$$
(17)
The last step in the procedure is to commute, for example, $`S_+`$ and $`\alpha _{12}`$, with the purpose of collecting all $`S_\pm `$ together. If one begins this rearrangement from the right side of Eq. (17), each of such permutations produces two terms. One of them consists of the product of the partially ordered and the remaining non-ordered parts, while the second one contains commutators of $`S_+`$ and $`\alpha _{12}`$. Repeating the procedure leads to the recurrent equation for the part of the transfer matrix containing only products of commuting diagonal matrices $`S_\pm `$

$$T_{N,1}=𝒪_{N,1}+\underset{k\mathrm{even}}{}T_{N,p_{k+1}}\mathrm{\Gamma }_{n_k}𝒪_{p_{k1},1},$$
(18)
where we have used the fact that $`\alpha _{12}\alpha _{21}=1`$, and introduced the commutator $`\mathrm{\Gamma }_n=[S_+^n,\alpha _{12}]`$. Notation $`𝒪_{p_{k1},1}`$ stays for the product of all $`S_\pm `$ between sites $`1`$, and $`p_{k1}`$: $`𝒪_{N,1}=\mathrm{}S_+^{n_2}S_{}^{n_1}`$, where $`p_k=_{l=1}^kn_l`$. Retaining the first term only and using Eq. (14) we arrive at Eq. (11).

Eq. (18) can be used to obtain a representation for the transfer matrix in terms of the polynomial in powers of the parameter $`\delta `$, which enters the expression for commutators $`\mathrm{\Gamma }_n`$, and can be considered as a small parameter in the case of weak disorder. However, the resulting expression for the transfer matrix cannot be interpreted as a weak disorder expansion for LE. The reason for this is the fact that the power of $`\delta `$ is determined by the number of segments with different values of the potential, and the whole polynomial is actually an expansion in terms of the number of “jumps” between $`V`$ and $`V`$. For instance, the term with $`\delta ^2`$ results from configurations with only one such jump, and there are only $`N`$ such configurations. In the limit $`N\mathrm{}`$ they give zero correction to the LE (in addition to the exponentially small probability of such configurations). Thus, we can see that for long enough samples only terms of the order of $`𝒩N/2`$ make a significant contribution to the LE because the number of such terms is proportional to the binomial coefficient $`C_N^𝒩`$ and, therefore they give a linear in $`N`$ contribution to the $`\mathrm{log}T`$. The actual number of configurations contributing to LE depends upon the average length of the segments with constant potential; and the number of jumps, corresponding to the optimal configurations can be estimated as $`N/l_c`$. It is clear, therefore, that corrections to the first term of Eq. (18) decrease with increasing $`l_c`$. We thus expect that the approximation should work when $`l_{loc}/l_c1`$. It actually becomes exact in the limit of infinite $`r_c`$ and $`V.`$Pastur\_Figotin ; LargeCoupling  This can be understood from the following arguments. Deviations from Eq. (12) are caused by the interference of the scattered waves. But in the second spectral region such scattered waves have to tunnel through the barriers, where the magnitude of the scattered waves decreases and the interference is destroyed. Obviously, the effectiveness of the destruction depends on the width of the potential barriers, and the rate of decay of the wave functions in them. These characteristics are determined by $`l_c`$ and the magnitude of the potential respectively.

In the spectral region, $`2V<|E|<2+V`$ only $`\gamma _0(E;V)`$ in Eq. (II) is different from zero, and we have for $`\sigma ^2`$

$$\sigma ^2L=\gamma ^2l_c.$$
(19)
This expression is obviously different from SPS Eq. (2), and is valid when $`\gamma l_c1`$. In this case, scaling described by the function $`\tau (\kappa )`$ is no longer valid, but Eq. (19) suggests that a new type of scaling appears. It can be conveniently described with the help of a new scaling function, $`S`$,

$$S=\frac{\sigma ^2L}{\gamma ^2l_c}1.$$
(20)
The transition between the two scaling regimes is controlled by the parameter $`\gamma l_c`$. If the correlation radius is so small that this parameter remains much less than unity for the entire second spectral region, the regular scaling behavior expressed by $`\tau (\kappa )`$ persists. The behavior expressed by Eq. (20) occurs when $`\gamma l_c`$ is greater than unity for the entire region $`2V<|E|<2+V`$. More detailed analysis of the transition between the two types of behavior requires a numerical approach, which we discuss in the next section of the paper.

In the third spectral region, $`|E|>2+V`$, $`S`$ decreases with energy. This behavior can be approximately described as

$$S\frac{\left(\sqrt{2}\sqrt{|ϵ|1}1\right)^2}{\left(\sqrt{2}\sqrt{|ϵ|1}+1\right)^2}$$
(21)
for $`0<|ϵ|11`$, where $`ϵ=(2|E|)/U`$ and

$$S\left(\frac{1}{ϵ}\right)^2\frac{1}{(\mathrm{log}|U|+\mathrm{log}|ϵ|)^2}$$
(22)
for $`|ϵ|1`$.

## III Scaling properties of LE: numerical simulations

For the numerical analysis, we calculate the LE iteratively using Eq. (7) and the usual technique of renormalization of the resultant vector after every 10 iterations.Kramer  The length $`l_s`$ was calculated according to Eq. (3), where the integrated density of states was obtained using the phase formalism LGP  with consequent averaging over all realizations. To investigate scaling properties we kept the length of chains much larger than the correlation length, $`l_c,`$ and the localization length, $`l_{loc},`$ for all magnitudes of the random potential $`U`$ and values of the energy $`E`$ in the vicinity of the second spectral region. Statistics were collected from $`40000`$ realizations. Figs. 1a and b present results of numerical computations of the Lyapunov exponent. The numerical value of the correlation radius in these figures is such that $`l_c>l_{loc}`$ for the entire $`2V<|E|<2+V`$ spectral region. A comparison of numerical results and the approximation (11) is provided for both average LE and its variance. For energies inside the first region this approximation gives zero, as discussed above. In the second region, one can see that Eq. (11) provides a reasonably good approximation for $`\gamma `$, and even better agreement for $`\sigma ^2`$. In the third spectral region the agreement between Eq. (11) and the numerical calculations becomes perfect.

We also studied numerically the scaling properties of our model. Let us recall that for systems without correlations SPS is controlled by the ratio $`\kappa =l_{loc}/l_s`$ of two macroscopic lengths existing in the system. In the presence of the correlated disorder the situation is more complicated because a new length related to the correlation radius $`r_c`$ should be taken into account. In our model this length is $`l_c`$ the average length of regions with constant potential. If $`l_c`$ is the smallest length in the system, $`l_c<(l_{loc},l_s),`$ the presence of correlations does not affect the scaling properties. Computer simulations show that the scaling parameter $`\tau `$ behaves in this case qualitatively similar to a system with uniformly distributed non-correlated disorder. In Fig. 2 the dependence of $`\tau `$ on $`\kappa `$ is shown for both a tight-binding model with uniformly distributed ($`UV_nU`$) on-site potential (with $`U`$ ranging from $`0.09`$ to $`0.145`$) and the model under consideration. For the latter, potentials ($`V=0.0010.004`$) and correlation radii ($`r_c=2`$ and $`3`$) where chosen such as to have $`l_c`$ the smallest length up to the genuine spectral boundary.

However, when $`l_c>l_{loc}`$, the behavior of $`\sigma ^2`$ changes significantly, and the character of this change depends upon the position of the spectral point, $`E_c,`$ at which $`l_c=l_{loc}(E_c)`$. If $`l_c`$ is so large that $`E_c`$ falls in the first spectral region $`|E|<2V`$, the scaling function $`\tau `$ remains essentially equal to unity up to the very small vicinity of the boundaries of this interval. Inside a very small neighborhood of points $`|E|=2V`$ any scaling behavior disappears, but immediately outside of this neighborhood inside the second spectral interval $`2V<|E|<2+V`$ a new scaling presented by Eq. (20) emerges. This situation is presented in Fig. 3. In the third spectral region the function $`S`$ decreases in agreement with the analytical formulas Eq. (21) and Eq. (22). It can be seen that the scaling in terms of variable $`ϵ=(2|E|)/U`$ predicted by Eq. (21) persists over quite a wide interval of energies.

In the case when the energy $`E_c`$ belongs to the interior of the second spectral region, the scaling properties are determined by the parameter $`\kappa `$. If $`\kappa (E_c)1`$ we return to the situation of the non-correlated disorder and if $`\kappa (E_c)>1`$ the system does not exhibit any scaling behavior in the second spectral region. The destruction of scaling in this case may also be illustrated by the behavior of $`\tau (\kappa )`$ for different values of $`l_c`$. When $`l_c`$ increases from the values corresponding to the non-correlated limit, the function $`\tau (\kappa )`$ starts changing as shown in the insert in Fig. 2 where the dependence of $`\tau `$ is depicted for three different sets of parameters of the random potential $`U`$ and $`r_c`$. It deviates from the universal (SPS) dependence at larger values of $`\kappa `$ with the increasing magnitude of the deviation depending upon both $`r_c`$ and $`U`$. It demonstrates in that way the absence of scaling properties.

## IV Conclusion

We considered statistical characteristics and scaling properties of the Lyapunov exponent in a random medium with correlated disorder, using a tight-binding model with a random potential described by a dichotomic process. A simple theoretical approach for estimation of the statistics of the Lyapunov exponent was suggested. Qualitative analysis based on this model allowed us to estimate the effect of correlations on scaling properties, and express this effect in terms of the ratio of the localization length and the correlation length, $`\theta =l_c/l_{loc}`$. When the correlation radius is much smaller than the localization length ($`\theta 1`$), the system behaves as though correlations are absent and demonstrates single parameter scaling, $`\sigma ^2Ll_{loc}1,`$ whose validity is governed by the ratio $`l_{loc}/l_s`$. In the opposite case, when $`\theta 1`$, the crossover to another scaling, $`\sigma ^2Ll_{loc}^2/l_c1,`$ specific for the model under consideration, takes place at values of energy close to the band gap.

## Acknowledgments

The authors would also like to thank Steve Schwarz for reading and commenting on the manuscript. This work was supported by AFOSR under Contract No. F49620-02-1-0305, PSC-CUNY grants, and partially by NATO Linkage grant No. 978090.


text:
# High-Energy Astrophysics and Cosmology

## 1 Introduction

There are two ways in which high-energy particles have appeared naturally in the Universe: one is via energetic astrophysical sources such as gamma-ray bursters (GRBs) or active galactic nuclei (AGNs), and the other is via the characteristically high particle energies in the very early Universe. In this talk, I illustrate the possible rôles of both these types of sources, and discuss some related open questions in relation to the cosmic-ray energy spectrum shown in Fig. 1. Here are a few examples.

Even relatively ‘low-energy’ parts of this spectrum in the range $`110`$ GeV are directly connected to ultra-high-energy particle physics, via their rôle in producing atmospheric neutrinos , one of our windows on Grand Unified Theories (GUTs). Somewhat higher-energy parts of the spectrum up to energies $`1`$ TeV are relevant to the experimental search for cold dark matter particles. Confirmed observations of GRBs have been limited to the MeV energy range, but there are unconfirmed reports of observations in the GeV and even TeV  energy ranges, and GRBs might even be responsible for the ultra-high-energy cosmic rays (UHECRs). Alternatively, these might be due to some exotic top-down mechanism involving the decays of supermassive particles produced in the very early Universe , or some other extreme astrophysical sources. Either way, they may provide a unique laboratory to look for violations of fundamental principles such as Lorentz invariance .

Before discussing these possible playgrounds for high-energy physics, let us first review the basis for the Big-Bang cosmology that plays an essential rôle in the following sections of this talk.

## 2 Big-Bang Cosmology

According to standard Big-Bang cosmology, the entire visible Universe is expanding homogeneously and isotropically from a very dense and hot initial state. Apart from the present Hubble expansion, the first piece of evidence for the Big Bang was the cosmic microwave background (CMB) radiation, which is thought to have been emitted when the Universe was about 3000 times smaller and hotter than it is today, with age $`3\times 10^5`$ y. The CMB has a dipole deviation from isotropy at the $`10^3`$ level, which this is believed to be due to the Earth’s motion relative to a Machian cosmological frame.

The second piece of evidence for the Big Bang was provided by the abundances of light elements seen in Fig. 2, which are thought to have been established when the Universe was about $`10^8`$ times smaller and hotter than it is today, with age $`1`$ to $`10^2`$ s. This nuclear ‘cooking’ must have occurred when the temperature $`T`$ of the Universe corresponded to characteristic particle energies $`1`$ MeV.

Before this time, when it was $`10^6`$ to $`10^5`$ s old, it is thought that the Universe made a transition from quarks and gluons to hadrons at a temperature $`T100`$ MeV. Previous to that, the electroweak transition when Standard Model particles acquired their masses would have occurred when the Universe was $`10^{12}`$ to $`10^{10}`$ s old, and the temperature $`T100`$ GeV. The CMB is thought to provide a window to an even earlier epoch, via its small-scale fluctuations $`\delta \rho /\rho `$, which show up at the $`10^5`$ level, as seen in Fig. 3. If these are due to quantum fluctuations during an inflationary epoch, $`\delta \rho /\rho (T/m_P)^2`$, telling us that the typical energies of particles in the Universe may once have approached $`10^{16}`$ GeV.

The fact that the CMB fluctuation spectrum is largest for partial waves $`\mathrm{}200`$, as seen in Fig. 3 , suggests that the total energy density of the Universe is close to the critical density: $`\mathrm{\Omega }_{tot}1`$. On the other hand, the Big-Bang nucleosynthesis calculations shown in Fig. 2 indicate that the density of baryons in the Universe today is far less than the critical density: $`\mathrm{\Omega }_B\rho _B/\rho _{crit}\mathrm{few}`$ %, a range supported by observations of the smaller-scale peaks in the CMB fluctuations shown in Fig. 3. What form does the missing energy take?

Observations of large-scale structure suggest that the total matter density $`\mathrm{\Omega }_m0.3`$, a value supported by a combination  of data on CMB fluctuations and high-redshift supernovae . Most of this $`\mathrm{\Omega }_m`$ is thought to be dark non-baryonic matter. What is the nature of this dark matter?

It certainly includes neutrinos, which are now thought to have non-zero masses , but these are probably insufficient to explain most of the dark matter. In any case, people who model the formation of large-scale structures would prefer more massive ‘cold’ dark matter particles that would have been non-relativistic when these structures began to grow. In addition to this cold dark matter, the CMB and other data apparently require about 2/3 of the critical density to be in the form of ‘dark energy’ in the vacuum , but I do not discuss the latter further in this talk. Instead, I concentrate on the dark-matter particles that might have signatures among high-energy cosmic rays.

## 3 Astrophysical Neutrinos

Let us first consider generic features of the density of relic neutrinos or similar neutral weakly-interacting particles. If neutrinos weigh less than $`1`$ MeV, their cosmological relic number density $`n_\nu `$ was fixed when $`T1`$ MeV and is essentially independent of their mass, hence their relic density $`\rho _\nu =m_\nu n_\nu `$ increases linearly with mass, rising above the critical density when $`m_\nu 30`$ eV. The neutrino mass density would be excessive for masses up to $`3`$ GeV, where a Boltzmann factor suppresses the $`\nu `$ number density sufficiently to push $`\mathrm{\Omega }_\nu `$ back below unity. The ‘neutrino’ density would be most suppressed for $`m_\nu m_Z/2`$, when ‘relic’ annihilation is most efficient. For larger ‘neutrino’ masses, the annihilation rate typically falls and the relic density (which is fixed when $`Tm_\nu /25`$) correspondingly rises, reaching the critical density for some ‘neutrino’ mass $`1`$ TeV.

Thus, there are three mass ranges where such a ‘neutrino’ might have a relic density of interest for astrophysics and cosmology: when $`m_\nu 10`$ eV, $`3`$ GeV or $`100`$ GeV to 1 TeV. In the first of these windows, the neutrino would constitute hot dark matter, in the latter two it would be cold dark matter. The middle option is excluded by a combination of experiments at LEP and direct dark-matter searches, and the third option is that exercised by the lightest supersymmetric particle , as we discuss later.

Observations of large-scale structures in the Universe favour the predominance of cold dark matter, as already mentioned, and can be used to set an upper limit on the sum of light neutrino masses :

$$\mathrm{\Sigma }_im_{\nu _i}<\mathrm{\hspace{0.33em}\hspace{0.33em}3}\mathrm{eV},$$
(1)
as seen in Fig. 4. The neutrino oscillation experiments force all three neutrino flavours to be essentially degenerate compared to (1), implying that $`m_{\nu _i}<1`$ eV for each species. This is stronger than the upper limit $`m_{\nu _e}<2.5`$ eV coming from the end-point of Tritium $`\beta `$ decay , and future observations of large-scale structure should improve the sensitivity to $`m_\nu 0.3`$ eV.

Astrophysical sources  have provided us with the first confirmed evidence for neutrino oscillations and (presumably) neutrino masses. The first astronomical image obtained with neutrinos was that of the Sun, via neutrino-electron scattering. Recently, comparative measurements of charged-current and neutral-current reactions by SNO  have established beyond any doubt that solar $`\nu _e`$ oscillate into some combination of $`\nu _\mu `$ and $`\nu _\tau `$, most probably with a relatively large mass-squared difference $`\mathrm{\Delta }m^26\times 10^5`$ eV2 and relatively large, but not maximal, mixing $`\mathrm{sin}^22\theta 0.8`$ (the LMA solution), as seen in Fig. 5 .

As you know, atmospheric neutrinos are produced by cosmic rays with energies mainly in the range $`110`$ GeV. For some time now, it has been established that $`\nu _\mu `$ also oscillate , probably mainly into $`\nu _\tau `$ with $`\mathrm{\Delta }m^22.5\times 10^3`$ eV2 and near-maximal mixing $`\mathrm{sin}^22\theta 1`$, as seen in Fig. 6. The available data on the primary cosmic-ray spectrum and particle production at accelerators enable the atmospheric neutrino flux to be calculated quite reliably, and further improvements will be possible using the data presented here from the AMS  and L3+C  Collaborations.

Long-baseline neutrino oscillation experiments are underway to refine and extend these astrophysical neutrino oscillation results. The K2K experiment already provides some confirmation of atmospheric $`\nu _\mu `$ oscillations, as seen in Fig. 7 , the KamLAND experiment is expected soon to test definitively the LMA solution for solar neutrinos , the MINOS experiment will look explicitly for the oscillatory pattern and have improved sensitivity to $`\nu _\mu \nu _e`$ oscillations , and the OPERA and ICARUS experiments should be able to observe the $`\tau `$ production expected following $`\nu _\mu \nu _\tau `$ oscillations .

## 4 Supersymmetry

As already mentioned, cosmic rays in the range from 100 GeV to 1 TeV may be the best place to look for supersymmetry: why? The primary theoretical motivation for expecting supersymmetry to appear in this energy range is provided by the hierarchy problem : why is the electroweak scale $`m_W`$ so much less than the Planck scale $`m_P10^{19}`$ GeV, which is the only candidate we have for a fundamental mass scale in physics? Equivalently, why is $`G_F1/m_W^2G_N=1/m_P^2`$?

You might just say, why not choose the value of $`m_W`$ and forget about the problem? Life is not as simple as that, because quantum effects in the Standard Model make very large corrections to the electroweak scale:

$$\delta m_W^2𝒪(\frac{\alpha }{\pi })\mathrm{\Lambda }^2,$$
(2)
where $`\mathrm{\Lambda }`$ is a cutoff representing the scale at which the Standard Model must be modified by introducing new physics. A mechanism for cutting the divergence (2) off in a natural way is provided by supersymmetry, which exploits the opposite signs in the quadratically divergent fermionic and bosonic corrections to the electroweak scale:

$$\delta m_W^2𝒪(\frac{\alpha }{\pi })|m_B^2m_F^2|,$$
(3)
which $`m_W^2`$ if $`|m_B^2m_F^2|1`$ TeV2. This argument therefore leads one to expects supersymmetric partners of Standard Model particles to appear at or below the TeV scale .

This argument for low-energy supersymmetry is supported circumstantially by the possibility it offers for unification of the gauge couplings. Such grand unification does not occur in the absence of supersymmetry, but is quite possible if supersymmetric particles weigh about 1 TeV, as seen in Fig. 8 .

In many supersymmetric models, the lightest supersymmetric particle (LSP) $`\chi `$ is stable, and a good candidate for cold dark matter , which provides a third general argument for the TeV mass scale. The relic energy density $`\rho _\chi =m_\chi n_\chi `$, where the relic number density

$$n_\chi \frac{1}{\sigma _{ann}(\chi \chi \mathrm{all})},$$
(4)
where a typical annihilation cross section $`\sigma _{ann}1/m_\chi ^2`$. Thus, the overall relic density increases with mass, and typically becomes too high when $`m_\chi >1`$ TeV.

However, the dark-matter annihilation rate may in exceptional circumstances be enhanced, reducing the supersymmetric relic density for a given mass, and thereby allowing larger relic masses. For example, if the LSP and the next-to-lightest supersymmetric particle $`\stackrel{~}{X}`$ have similar masses, $`n_\chi `$ may be suppressed by coannihilation processes : $`\sigma (\chi \stackrel{~}{X}\mathrm{all})`$, which can be important if $`m_{\stackrel{~}{X}}m_\chi /m_\chi 1/10`$. This coannihilation mechanism can provide an allowed ‘tail’ of parameter space extending out to larger $`m_\chi `$, as seen in Fig. 9. Such a tail may also happen when rapid annihilation through a direct-channel pole is possible, for example if $`m_\chi m_{H,Z,\mathrm{}}/2`$ .

The space of input supersymmetric fermion masses $`m_{1/2}`$ and boson masses $`m_0`$ is illustrated for one particular ratio $`\mathrm{tan}\beta `$ of the Higgs vacuum expectation values in the minimal supersymmetric extension of the Standard Model (MSSM) in Fig. 9 . We see at large values of the ratio $`m_{1/2}/m_0`$ a region excluded because there the relic particle would be charged, a possibility excluded by astrophysics. At small $`m_{1/2}`$ and/or $`m_0`$ we see experimental exclusions from the absences of the supersymmetric partners of the electron $`\stackrel{~}{e}`$ and of the $`W/H^\pm `$, and also of the Higgs boson $`H`$. A dark (green) shaded region is excluded by measurements of $`bs\gamma `$ decay, and a lighter (pink) shading shows regions favoured by the recent measurement of the anomalous magnetic moment of the muon .

Finally, the lightest (turquoise) shading in Fig. 9 picks out the region where the relic LSP density lies within the range favoured by cosmology: $`0.1<\mathrm{\Omega }_\chi h^2<0.3`$. We see that much of this region is disfavoured by the accelerator constraints, particularly the LEP Higgs limit $`m_H>114`$ GeV .

A set of ‘benchmark’ supersymmetric scenarios was recently proposed , that respect all the experimental and cosmological constraints on the minimal supersymmetric extension of the Standard Model (MSSM). As shown in Fig. 10, they indicate the range of options, rather than sample the parameter space in a ‘fair’ manner. The LHC has great possibilities for detecting supersymmetry, principally via events with missing energy and other signatures such as high-energy leptons and/or jets. These prospects for discovering supersymmetry may be compared with those for detecting astrophysical supersmmetric dark matter.

## 5 Search for Supersymmetric Dark Matter

Searches among cosmic rays with energies up to about 1 TeV provide several promising signatures for supersymmetric dark matter particles, that may enable this community to ‘scoop’ the LHC .

One possibility is to look for energetic gamma rays that may be emitted by LSP annihilations in the core of the Milky Way. The benchmark models indicate that these might have typical energies $`10`$ GeV, and detectors such as GLAST with a threshold as low as $`1`$ GeV might have a better chance, as seen in Fig. 11 . The prospects for these searches depend on the degree to which the dark-matter particle density may be enhanced in the core of the Milky Way, which is uncertain by orders of magnitude. In Fig. 11, a middle-of-the-road enhancement by a factor 200 has been assumed.

Another possibility is to look for positrons emitted by LSP annihilations in the halo of the Milky Way. In this case, the benchmark models indicate that energies $`100`$ GeV might be the most interesting, though the signal may be less promising than in the $`\gamma `$ case, as seen in Fig. 12 .

One of the most promising signatures is energetic muons produced in the Earth by energetic neutrinos emitted by LSP annihilations in the centre of the Sun or Earth. In this case, as seen in Fig. 13, all energies up to $`1`$ TeV might be important. According to our calculations , the prospects for detecting relic annihilations in the core of the Sun appear more promising than those in the centre of the Earth. Some upper limits on the energetic solar muon flux have already been produced, most recently by the AMANDA Collaboration , which already begin to exclude some more extreme supersymmetric models.

The most convincing evidence for supersymmetric dark matter might eventually come from direct searches for the scattering of relic particles on nuclei in the laboratory . Here the best chances seem to be offered by spin-independent scattering on relatively heavy nuclei. There has been a claim by the DAMA Collaboration  to have observed an annual modulation effect due to the scattering of dark matter particles, but this interpretation of their data has been largely excluded by the CDMS , EDELWEISS  and UKDMC experiments. In any case, reproducing the DAMA data would have required a scattering cross section much larger than predicted in the simple supersymmetric models studied in  and . As illustrated in Fig. 14, future large cryogenic detectors, such as that proposed by the Heidelberg group, would good chances in many supersymmetric scenarios.

## 6 Space-Time Foam

We know that space-time is essentially flat at large scales, but quantum gravitational effects are expected to cause large fluctuations on small scales in length and time:

$$\mathrm{\Delta }E,\mathrm{\Delta }\chi =𝒪(1)\mathrm{in}\mathrm{\Delta }x,\mathrm{\Delta }t=𝒪(1),$$
(5)
where the energy $`E`$, distance $`x`$ and time $`t`$ are all measured in Planck units $`10^{19}`$ GeV, $`10^{33}`$ cm and $`10^{43}`$ s, respectively, and $`\chi `$ is a generic dimensionless measure of topology. Are there any observable consequences of such quantum gravitational effects?

One suggestion has been that information might be lost across microscopic event horizons associated with such topological fluctuations, causing an apparent modification of microscopic quantum mechanics . Another has been that the gravitational recoil of the vacuum as an energetic particle passes by might modify special relativity, in such a way that the particle’s velocity might be reduced :

$$c(E)=c\times \left(1\frac{E}{M}+\mathrm{}\right),$$
(6)
where $`M`$ is a quantum-gravity mass scale that might $`m_P10^{19}`$ GeV.

The best probes of the possibility (6) may be provided by astrophysical sources of $`\gamma `$ rays: they have large distances $`D`$ and hence light propagation times $`t=D/c`$, and hence their light pulses may exhibit time delays

$$\mathrm{\Delta }t\frac{D}{c^2}(c(E)c).$$
(7)
Distant high-energy sources with short characteristic time scales $`\delta t`$ have the best sensitivity to the quantum-gravity scale :

$$M\frac{ED}{\delta t}.$$
(8)
From this point of view, interesting astrophysical sources include pulsars (relatively short distances, but short time scales), AGNs (moderate redshifts $`z`$ and relatively large time scales, but high energies) and GRBs ($`z1`$, $`\delta t10^2`$ s). There are confirmed observations of GRBs only at relatively low energies, but there have been reports of GeV or even TeV  photons emitted by GRBs.

Some time ago, we published an analysis of all the GRBs whose redshifts had then been measured, using data from the BATSE and OSSE instruments on the Compton Gamma-Ray Observatory (CGRO), finding no significant correlation of time-lag with $`z`$, and inferring that $`M>10^{15}`$ GeV . More recently, we have been making an improved analysis, benefiting from the measurements of more GRB redshifts, and using TTE data from BATSE, which has a finer time resolution, and making a wavelet analysis, whose result is shown in Fig. 15. This time, we find $`M>7.9\times 10^{15}`$ GeV . In the future, GLAST should be able to improve significantly on this sensitivity .

## 7 High-Energy Cosmic Rays and the LHC

Colliding pairs of protons with 14 TeV each, the LHC will be equivalent to a beam of particles with energies $`10^{17}`$ GeV striking a fixed target. The general-purpose detectors, ATLAS  and CMS , are designed to see ‘everything’ emitted with a centre-of-mass rapidity $`|\eta |<3`$. The LHCb experiment  will also be able to measure some diffractive physics. The LHC will also be able to collide beams of heavy ions with energies $`5`$ TeV per nucleon, with the ALICE experiment  as principal user, and could in principle also collide protons or deuterons with heavy ions. Also planned is the TOTEM experiment  to measure the total and elastic $`pp`$ cross sections, as well as some diffractive inelastic events. This experiment will have some capabilities to measure particles with $`|\eta |<3`$ or 5, and the LHC physics community would like to hear from the cosmic-ray community what its needs are . This is the subject of a special workshop at this meeting: we hope that the LHC can contribute to refining models of ultra-high-energy cosmic rays. We should like to know what particles the cosmic-ray community would like to see measured in what rapidity range.

## 8 Ultra-High-Energy Cosmic Rays

As you know, extended-air-shower (EAS) experiments  have reported an apparent excess of events beyond the GZK cutoff, but this is not confirmed by HiRes , the largest fluorescence detector, as seen in Fig. 16. At issue are the absolute and relative energy calibrations of the EAS and fluorescence techniques. The former depends on the models that we hope the LHC will help refine, and the latter rely on normalizations of fluoresecnce lines around 390 nm. A recent experiment seems to find a discrepancy with previous measurements in this region, and the apparent conflict should be clarified .

There are two main classes of models for ultra-high-energy cosmic rays (UHECR): ‘bottom-up’ and ‘top-down’. The former postulate acceleration by astrophysical sources, which should have a minimum size $`R`$:

$$R\frac{100}{Z}\left(\frac{E}{10^{20}\mathrm{eV}}\right)\left(\frac{\mu \mathrm{G}}{B}\right),$$
(9)
where $`Z`$ is the atomic number, $`E`$ is the energy and $`B`$ is the magnetic field bending the UHECR. Alternatively, the maximum energy attainable in any given source is

$$E\mathrm{\hspace{0.33em}10}^{18}Z\left(\frac{R}{\mathrm{Kpc}}\right)\left(\frac{B}{\mu \mathrm{G}}\right)\mathrm{eV}.$$
(10)
Popular candidate sources include GRBs and AGNs, with neutron stars and colliding galaxies also proposed.

The alternative class of ‘top-down’ models postulates the production of UHECR by GUT-scale physics, such as topological defects or the decays of metastable superheavy relic particles  \- which have some chances of being produced with interesting cosmological relic densities . These models predict some anisotropy if the particles have a typical halo distribution , and could also exhibit clustering if the halo is lumpy . A characteristic of ‘top-down’ models is that they predict large numbers of photons among the UHECR, and possibly even supersymmetric particles !

Other suggestions for UHECR include the collisions of ultra-high-energy neutrinos with non-relativistic massive neutrinos to produce $`Z`$ bosons . It has also been suggested  that the type of modification of Lorentz kinematics mentioned earlier might also help distant sources evade the GZK cutoff.

There will be much discussion of these ideas during this meeting. We are all glad that the Auger experiment in Argentina is progressing well . Its high statistics and combination of the EAS and fluorescence techniques should pin down the existence of UHECR and discriminate between rival models. In the longer run, the EUSO project  would be able to provide even higher statistics.

## 9 Conclusions

This talk has provided only a brief review of the interfaces between high-energy physics, astrophysics and cosmology. As we have seen, high-energy cosmic-ray physics may provide much useful information to particle physicists, for example concerning neutrinos and dark matter, and possibly even quantum gravity. Likewise, accelerator experiments can provide useful input for cosmic-ray experiments, for example by testing, calibrating and validating simulation codes.

During this talk, we have discussed both localized astrophysical sources of high-energy particles and global sources provided by the very early Universe. Although lower-energy cosmic rays are believed to originate from local sources within our own galaxy, the origin of the highest-energy cosmic rays is still unknown. Both astrophysical sources and cosmological origins are being actively proposed. Unravelling the origins of the UHECR will surely require an active dialogue between accelerator and cosmic-ray physicists.