text:
# Cluster Growth in two- and three-dimensional Granular Gases

## I Introduction

Granular media are interesting multi-particle systems with a rich phenomenology Herrmann et al. (1998); Cafiero et al. (2000); Pöschel and Luding (2001); Vermeer et al. (2001); Kishino (2001). They can form a hybrid state between a fluid and a solid, where the behavior is controlled by the balance between energy input and energy dissipation. Energy input leads to a reduction of the density due to more collisions and increasing pressure, so that the material can flow. In the absence of energy input, e. g., in freely cooling systems, granular materials become denser, i. e. they locally “solidify” due to dissipation. Because of mass conservation, the local densification is accompanied by a density decrease in other parts of the system, giving rise to complex patterns and structures, with an interesting time evolution. However, theoretical approaches are non-classical and appear often extremely difficult, so there is still active research directed towards the better understanding of granular media.

The subject of this paper is the pattern formation via clustering in a dissipative, freely cooling granular gas Goldhirsch and Zanetti (1993); Goldhirsch et al. (1993); McNamara and Young (1996); Luding and Herrmann (1999); Nie et al. (2002); Das and Puri (2003a, b). The basic idea of clustering is that in an initially homogeneous freely cooling granular gas, fluctuations in density, velocity, and temperature cause a position dependent energy loss. Due to strong locally inhomogeneous dissipation, pressure and energy drop rapidly and material moves from “hot” to “cold” regions, leading to even stronger dissipation and thus causing the density instability with ever growing clusters until eventually, clusters reach system size (see Fig. 1).

We investigate this clustering instability with respect to different dissipation rates and different system sizes. Even though many of our findings are empirical, we attempt to reduce the complexity of the evolution of the system by use of simple scaling laws.

In section II, we explain in detail the simulation approach. The results of our numerical experiments are discussed in section III. Finally, in section IV a short summary and a discussion are given.

## II Simulation Details

A granular gas can be idealized as an ensemble of hard spheres in which the energy loss that accompanies the collision of macroscopic particles is modeled with a single coefficient of restitution $`r`$. In the simplest case the particles are identical in size and mass and there are no inter-particle forces between collisions.

Details about initial and boundary conditions are given in section II.1. The microscopic dynamics of the motion and the collision of the particles is discussed in section II.2 and the simulation method is explained in section II.3. Section II.4 deals with the inelastic collapse, a problematic artefact of the hard sphere model with dissipation.

### II.1 Initial and Boundary Conditions

The simulation volume consists of a box with equal side length and periodic boundary conditions in two or three dimensions.

An initial state with random particle positions and velocities is prepared in the following way: The particles first sit on a regular lattice and have a Maxwellian velocity distribution with a total momentum of zero. Then the simulation is started without dissipation and runs for about $`10^2`$ collisions per particle. This state is now used as initial configuration for the dissipative simulations.

### II.2 Microscopic Dynamics

Between collisions no forces act upon the particles and they move at constant velocity.

The particles are idealized as hard spheres. This means that collisions take infinitesimal time and involve only two particles. Conservation of momentum leads to the collision rule Goldhirsch and Zanetti (1993); Luding (1998)

$`𝐯_{}^{}{}_{1/2}{}^{}`$ $`=𝐯_{1/2}{\displaystyle \frac{1+r}{2}}\left(\widehat{𝐤}(𝐯_1𝐯_2)\right)\widehat{𝐤},`$ (1)
where a prime indicates the velocities $`𝐯`$ after the collision, $`\widehat{𝐤}`$ is a unit vector pointing along the line of centers from particle 1 to particle 2, and $`r`$ is the coefficient of restitution. The relative tangential velocity does not change during a collision, the relative normal velocity changes its sign and is reduced by a factor $`1r`$. Energy dissipation is proportional to $`\lambda =1r^2`$, so that the elastic limit $`r=1`$ implies $`\lambda =0`$, i. e. no dissipation, while $`r<1`$ implies $`\lambda >0`$.

### II.3 Event-Driven Molecular Dynamics

The simulation of hard spheres can be handled efficiently with event-driven molecular dynamics Lubachevsky (1991); Miller and Luding (2004). The collisions are the events which have to be treated by the algorithm. Between these collisions the particles move on trivial trajectories and so the algorithm can easily compute the point of time $`t_{12}`$ of the next collision of two particles 1 and 2 as

$$\begin{array}{c}t_{12}=t_0𝐫_{12}𝐯_{12}/v_{12}^2\hfill \\ \hfill +\sqrt{(𝐫_{12}𝐯_{12})^2\left(r_{12}^2d^2\right)v_{12}^2}/v_{12}^2,\end{array}$$
(2)
where $`𝐯_{12}=𝐯_2(t_0)𝐯_1(t_0)`$ and $`𝐫_{12}=𝐫_2(t_0)𝐫_1(t_0)`$ are the relative velocities and positions of the particles at time $`t_0`$, and $`d`$ is the diameter of a particle.

The algorithm processes the events one after another. After a collision the positions and velocities of the two involved particles are updated; the state of all other particles remains unchanged. For the two colliding particles, new events are calculated and the next future event is stored in the event priority queue for both particles. The next event is obtained from the priority queue, the new positions and velocities after the collision for the collision partners are updated, and so on. Neighborhood search is enhanced with standard linked cell methods Allen and Tildesley (1987), where the cell change of a particle is treated as a new event type. The details of the algorithm can be found in Lubachevsky (1991, 1992); Miller and Luding (2004).

### II.4 Avoiding the Inelastic Collapse with the TC Model

Our model makes use of hard spheres with an infinitely stiff interaction potential. This is an idealization of real physical particles and can lead to the dramatic consequence of inelastic collapse: An infinite number of collisions occurs in finite time. This singularity is unphysical, of course, and a major drawback for numerical simulations, too. But it has been shown that one can circumvent this artefact of the model in the following way Luding and McNamara (1998): If two consecutive collisions of a particle happen within a small time $`t_c`$, dissipation is switched off for the second collision. This timestep $`t_c`$ corresponds to the duration of the contact of physical particles.

There exist other deterministic and random models which prevent inelastic collapse. Even though many of them lack a solid theoretical background and physical motivation, their details should be insignificant for the physical evolution of the system anyway, since only a small negligible fraction of the particles in the system is involved in the inelastic collapse. For an extensive discussion see Luding and McNamara (1998).

## III Numerical Experiments

The simulation is started from a homogeneous system, prepared as described in section II.1. Depending on the dissipation $`\lambda `$, the density $`\rho `$, and the number $`N`$ of particles, the system remains in the homogeneous cooling state (HCS) for some time, until clustering starts and the system becomes inhomogeneous.

We discuss the evolution of the kinetic energy and the collision frequency of the system in section III.1. In section III.2 we investigate the clustering by means of appropriate measures of the cluster size distribution. In section III.3 and III.4 we encounter interesting parallels with percolation theory and discuss several critical exponents.

### III.1 Kinetic Energy and Collision Frequency

Dissipative collisions lead to a decay of the kinetic energy and the collision frequency (see Fig. 2 and Fig. 3). From these figures three different regimes can be clearly distinguished. They are:

(1) The homogenous cooling state (HCS) at the beginning, when no clusters have formed yet, is well understood Haff (1983); Luding et al. (1998). The decay of the kinetic energy $`E`$ is governed by the equation

$`E(\tau )`$ $`={\displaystyle \frac{E(0)}{(1+\tau )^2}},`$ (3)
with the scaled time $`\tau =\frac{\lambda }{2D}\frac{t}{t_E}`$. $`D`$ is the dimension of the system and $`t_E`$ is the initial Enskog collision time $`t_E=(\sqrt{\pi }d)/(2\sqrt{D}v\rho g_d(\rho ))`$. $`v`$ is the mean velocity of a particle, $`d`$ its diameter, $`\rho `$ is the volume fraction, and $`g_d`$ is the contact probability. In 2D $`g_d(\rho )=(17\rho /16)/(1\rho )^2`$ and in 3D $`g_d(\rho )=(1\rho /2)/(1\rho )^3`$.

The evolution of the collision frequency per particle $`f_c(\tau )`$ with time is given by

$`f_c(\tau )`$ $`=t_E^1(0)\sqrt{{\displaystyle \frac{E(\tau )}{E(0)}}}`$ (4)
and it is the natural time scale controlling the evolution of the system in the HCS.

(2) When clusters start to grow, the decays deviate from these laws. In the cluster growth regime, the decay of the kinetic energy slows down. Furthermore, the collision frequency fluctuates a lot and can even increase. Note that the collision frequency is no longer a natural time scale of the system, since it is mainly determined by cluster-cluster collisions, where it increases strongly. The deviation from Eq. (3) occurs earlier and is more dramatic for larger dissipation $`\lambda `$, i. e. smaller $`r`$. However, the cluster growth regime is characterized by an energy decay around $`E\tau ^1`$, independently of $`r`$, cf. Nie et al. (2002). (In 3D the exponent might be slightly larger than unity; a best fit yields $`E\tau ^{1.1\pm 0.1}`$.) In contrast, the evolution of the collision frequency with time depends on $`r`$ during cluster growth. This regime is more distinct for large dissipation $`\lambda `$ and large system sizes $`N^{1/D}`$ and can e. g. barely be seen for $`r=0.98`$ in Fig. 3.

(3) Finally, when the largest cluster in the population has reached system size in the /em saturation regime, the cooling resembles the homogeneous cooling state in so far that $`E(\tau )\tau ^2`$ and $`f_c(\tau )\tau ^1`$, even if the latter still shows large fluctuations. In Fig. 2 this regime is not clearly visible, because the simulation times are not long enough for this system.

More details of the 2D situation were discussed in McNamara and Young (1996); Luding and McNamara (1998) and references therein. Further studies on the three-dimensional systems are in progress.

### III.2 Cluster Growth

The energy loss of the particles first leads to a reduced separation velocity after collision and eventually to the formation of clusters. But the definition of a cluster suffers from the fact that it takes an infinite number of collisions for the particles to stay in permanent contact with each other Luding and Herrmann (1999). So we use the following (geometrical) definition: Two particles belong to the same cluster, if their distance is smaller than $`s=0.1`$ particle diameters. The choice of $`s`$ is arbitrary and shifts the results; the qualitative behavior of the quantities discussed below does not depend on $`s`$, as long as it is neither too small nor too large 111 For a detailed study of different $`s`$ in 2D see Luding and Herrmann (1999); in 3D, the $`s`$ dependence of some quantities seems to be stronger than for others..

The moments $`M_k`$ of the cluster size distribution are defined as

$`M_k`$ $`:={\displaystyle \frac{1}{n_c}}{\displaystyle \underset{i}{}}i^kn_i,`$ (5)
where $`n_c`$ denotes the total number of clusters and $`n_i`$ the number of clusters of size $`i`$.

In many cases there are a lot of small clusters and one large cluster of size $`N_x`$, which contains a macroscopic fraction $`m_x:=N_x/N`$ of the total number $`N`$ of particles. Therefore, we also define reduced moments $`M_k^{}`$, which do not include the largest cluster.

In Fig. 4 and Fig. 5 the growth of the clusters can be seen on the basis of the first and second moments. After several collisions, particles start to cluster and the moments of the cluster size distribution grow until they reach their “saturation” value 222Saturation means that there is a transition to a regime with a much slower change of the moments. The dynamics is not finished, however, since the large cluster still can collect particles or break into pieces.. A numerical analysis reveals that the increase in $`M_1`$ and $`M_2`$ is mainly due to one large cluster which grows until it reaches its maximum size. In this final state this cluster contains a macroscopic fraction $`m_x`$ of the particles (see Fig. 6).

The onset time of cluster growth and also the final size $`m_x`$ of the large cluster depend strongly on the restitution coefficient $`r`$ (see Fig. 4 and Fig. 5). At low dissipation rates, for a long time nothing interesting happens and finally small and strongly fluctuating clusters appear. High dissipation leads to almost immediate cluster growth and a very large cluster at last. On the other hand, the system size $`N`$ does not seem to affect the behavior of the system provided that $`N`$ is not too small (see Fig. 6).

It has been shown McNamara and Young (1996); Sela and Goldhirsch (1995); Goldhirsch and Zanetti (1993), that for infinite system size the system is always unstable to the formation of clusters, whereas for a finite system size the density $`\rho `$ and the dissipation $`\lambda =1r^2`$ must not be too small. With our choice $`N10^5`$ and $`\rho =0.25`$ we expect no cluster formation below the critical dissipation $`\lambda _c=10^4`$ in our 2D system and $`\lambda _c=10^2`$ in our 3D system McNamara and Young (1996). This is in good agreement with our results for the 3D system, where cluster formation still can be seen with dissipation $`\lambda =0.04`$, but not for $`\lambda =0.004`$ (see Fig. 5).

As can be seen in Fig. 7 (left) the onset $`\tau _c`$ of the growth of the large cluster in 2D does not depend on the restitution coefficient $`r`$ explicitly. (Implicitly, $`\tau _c`$ is dependent on $`r`$ via rescaled time, of course.) In contrast $`\tau _c`$ does strongly depend on $`r`$ in a 3D system (see Fig. 7 (right)).

An empirical formula, which gives the dependency rather accurately over orders of magnitude in $`\lambda `$, is

$`\tau _c(\lambda )`$ $`=\tau _1+c{\displaystyle \frac{(1\lambda )^2}{\lambda \lambda _c}},`$ (6)
where $`\tau _1=0.24\pm 0.02`$ and $`c=36\pm 5`$ are fit parameters (for $`s=0.1`$).

In the limit $`\lambda 1`$, one has $`\tau _c\tau _1`$, where $`\tau _1`$ is proportional to the rescaled collision frequency and corresponds to $`t=2D\tau _1t_E/\lambda 1.5t_E`$. This constant stems from the fact that the particles need at least one collision to cluster, however, it also depends on $`s`$. In the limit $`\lambda \lambda _c`$ the time $`\tau _c`$ diverges: the particles never start to cluster 333 Since we focused on strong dissipation in this study, rather than the limit $`\lambda \lambda _c`$, we can not draw a conclusion on the functional behavior from our limited data..

of magnitude in Rescaling time in Fig. 7 (right) according to Eq. (6) results in Fig. 8 (right). There, the onset of the clustering happens approximately at the same rescaled time, whereas these times differ by 4 orders of magnitude in Fig. 7 (right). The fluctuations in $`\tau _c`$ are not systematically dependent on $`r`$.

Another difference between the 2D and 3D system in Fig. 7 is the fact that the cluster growth in 2D starts very smoothly, whereas the beginning of the clustering is very sharp in 3D.

### III.3 Critical exponents in 3D

As the moments of the cluster size distribution are clearly dominated by the large cluster, we will now study the reduced second moment $`M_2^{}`$, which does not include the large cluster. Fig. 9 shows that $`M_2^{}`$ is small most of the time, except for a peak at $`\tau _c`$. This gives us a clean definition of $`\tau _c`$. But what is more interesting, parallels to percolation theory arise. In percolation theory Stauffer and Aharony (1994) one studies the scaling behavior of certain quantities depending on the occupation probability $`p`$ around the percolation threshold $`p_c`$. One result for the reduced second moment $`M_2^{}(p)`$ is, e. g.,

$`M_2^{}(p)\left|pp_c\right|^\gamma ,`$ (7)
where $`\gamma `$ is a universal critical exponent. In order to transfer this result to cluster growth, we replace $`p`$ with $`\tau `$ and $`p_c`$ with $`\tau _c`$:

$`M_2^{}(\tau )\left|\tau \tau _c\right|^\gamma `$ (8)
Indeed we find in Fig. 9 the same power law behaviour. All data for different system sizes collapse on the same master-curve and even the exponent $`\gamma =1.8\pm 0.1`$ is in agreement with the 3D percolation problem result $`\gamma =1.80`$ after the beginning of clustering.

Now, if we study the amplitude of this peak, percolation theory tells us that the maximum sits at

$`p_{max}p_c(1aL^{1/\nu }),`$ (9)
where $`L`$ is the system size, $`a`$ is a constant, and $`\nu `$ is another critical exponent. Thus we expect:

$`M_2^{}(p_{max})L^{\gamma /\nu }.`$ (10)
As the system size $`LN^{1/3}`$, the maximum value of the peak in Fig. 9 scales as

$`M_{2,max}^{}(N)N^{\gamma /3\nu }.`$ (11)
Here we can also verify the power law behavior (see Fig. 10). Our simulations yield the result $`\gamma /3\nu =0.77\pm 0.09`$, which leads to $`\nu =0.78\pm 0.14`$. Within the rather large margin of errors this result is close to the 3D percolation problem $`\nu =0.88`$, too.

A third exponent $`q`$ 444The common notation for this exponent is $`\tau `$. In order to avoid name conflicts we call it $`q`$. is given by the cluster size population at time $`\tau _c`$:

$`n_i(\tau _c)i^q,`$ (12)
where $`n_i`$ is the number of clusters of size $`i`$.

Fig. 11 yields $`q=2.2\pm 0.2`$. This result is in agreement with the prediction of percolation theory $`q=2.2`$, too.

Is this similarity of clustering in 3D and percolation coincidence or is there a deeper reason?

### III.4 Consequences

Percolation theory makes many universal predictions for randomly disordered systems. The particles in a granular gas form such a disordered system. There have been other attempts to apply percolation theory to granular systems Tobochnik (1999), and also more advanced phase ordering models have been used to parallel the clustering dynamics Das and Puri (2003a, b). Our results indicate that such attemps can be justified, even though the percolation problem is purely static, whereas the clustering of granular gases also involves the dynamics, i. e. momentum and energy.

The only thing that might seem strange at first sight is the derivation of Eq. (8). Why can we replace the purely static quantity $`p`$ with time $`\tau `$?

In order to introduce the occupation probability $`p`$ in the clustering problem, we define it as

$`p(\tau ):={\displaystyle \frac{M_1(\tau )M_1(\tau ^{})}{M_1(\tau )}}`$ (13)
where $`\tau ^{}\tau _c`$ is a time shortly after the first few collisions. Now, we examine $`p(\tau )`$ around $`\tau _c`$.

Fig. 12 (right) shows that in 3D, $`p`$ is almost proportional to $`\tau `$. If we make use of this linear relation and insert it in Eq. (7), we arrive at the postulated formula Eq. (8).

In contrast, in our 2D system the growth of the large cluster happens at a rather large time $`\tau _c`$. (Note, that $`\tau _c`$ is not very clearly defined in the 2D case, because the onset of the cluster growth is very smooth.) In this late stage, the relation between $`p`$ and $`\tau `$ is highly non-linear, see Fig. 12 (left). Thus we could not find universal critical exponents for $`\tau `$ there. However, it might be possible to find a linear relation for a different set of parameters, e. g. for a different volume fraction, too. Whether a 2D system and a 3D system with the same parameters (like, e. g., volume fraction) are in fact comparable is another open issue.

## IV Summary and Discussion

The evolution of freely cooling granular systems can be divided into three regimes. First, the system is in the homogeneous cooling state (HCS). Then, in the cluster-growth regime, clusters begin to develop and grow. Finally, in the saturation regime, the clusters merge to practically one large cluster, which grows until it reaches system size. Besides the macroscopic fraction of particles in the large cluster, there are still many small clusters with interesting statistics.

In the HCS, the decay of the kinetic energy $`E`$ and the collision frequency $`f_c`$ can be described by the simple analytical expressions $`E(\tau )(1+\tau )^2`$ and $`f_c(\tau )(1+\tau )^1`$. The collision frequency is the natural time scale here, mainly determined by the density and the dissipation rate of the system. For strong dissipation, at short range, particles already stay closer together after only one collision, so that the moments of the size distribution change rather early due to a (short-range) change in the radial pair-distribution (data not shown).

In the cluster-growth regime, the collision frequency shows large fluctuations because of cluster-cluster collisions and cannot be predicted during cluster growth, because it changes erratically and discontinuously. The energy decay is characterized by $`E\tau ^1`$, where the accuracy of the exponent is limited in 3D due to the comparatively short duration of the cluster growth regime. However, this regime shows interesting differences between two and three dimensions. In 3D, cluster growth can be described by a power law behavior with the same critical exponents as in percolation theory. The onset of this cluster growth $`\tau _c`$ is very sharp and does strongly depend on dissipation $`\lambda `$. In contrast, we could not find a similar behavior in 2D, where the onset of the cluster growth is very smooth and depends on dissipation $`\lambda `$ only implicitly.

When cluster growth has reached a dynamic equilibrium in the saturation regime, the system is dominated by one large cluster which contains a macroscopic fraction of the system. Note that this regime still has an interesting dynamics and smaller clusters interact with the large cluster. Neither the small clusters nor the large one are static and the latter loses or eats up particles or smaller clusters. Kinetic energy and collision frequency still fluctuate, but are governed by the equations $`E(\tau )\tau ^2`$ and $`f_c(\tau )\tau ^1`$. This means the evolution in time is similar to the homogeneous cooling state.

A very interesting observation is the similarity of clustering in 3D and percolation with respect to the critical exponents. Even for the dynamics of cluster growth we have found the same exponents as for the occupation probability in percolation theory. We have provided a tentative explanation which is based on the linear relation of the occupation probability $`p`$ and time $`\tau `$ around the onset $`\tau _c`$ of cluster growth in 3D. In 2D we did not find a power law behavior because the growth of the large cluster happens when the linear relation between $`p`$ and $`\tau `$ has disappeared — at least for the given set of parameters.

However, there still remain some open questions: Why does $`\tau _c`$ depend on $`r`$ in 3D, but not in 2D? Why does the growth of the large cluster proceed smoothly in 2D and rather sharply in 3D? And why are all these differences present while the energy decay is proportional to $`\tau ^1`$ in both cases? Detailed studies of these questions and also of the cluster size probability distribution are in progress as well as a more systematic study of the cluster definition in 3D.

###### Acknowledgements.

This research was supported by the DFG, the SFB 382, the SFB 404, and LU 450/9-1. We thank Hans Herrmann and Sean McNamara for helpful discussions.


format: text

text:
# Shape and surface anisotropy effects on the hysteresis of ferrimagnetic nanoparticles

## Abstract

We present the results of Monte Carlo simulations of a model of a single maghemite ferrimagnetic nanoparticle with the aim to clarify the role played by the increased anisotropy at the surface and by the shape (spherical or elliptical) of the particle on the magnetization processes at low temperatures. The formation of hedgehog-like structures for high enough surface anisotropy is responsible for a change in the reversal mechanism of the particles.

url\]http://www.ffn.ub.es/oscar and

The progressive size reduction of magnetic particles with technological application has brought renewed interest on the issue of surface and finite-size effects in their magnetic properties. In particular, it is an old-known fact  that, when entering the nanometric range, the reduction of lattice symmetry at the surface of the particle results in an enhancement of magnetic surface anisotropy with respect to bulk values which is not clearly understood at present, although a number of computer simulations have been recently performed in order to clarify this issue . Here we present the results of Monte Carlo simulations of a model of a ferrimagnetic nanoparticle with the spinel lattice structure of $`\gamma `$Fe2O3 (maghemite). We consider an extension of our previous model in  to the Heisenberg case, including the term $`H_{\mathrm{anis}}/k_B=_{i=1}^N\left(K_C(S_i^z)^2+K_S(\stackrel{}{S}_i\widehat{n}_i)^2\right)`$ to account for finite anisotropy. Here surface spins are allowed to have radial surface anisotropy $`K_S`$ distinct from the uniaxial anisotropy $`K_C`$ of those in the core.

By using the standard Metropolis algorithm for continuous spins , we have computed the spin configurations of ellipsoidal particles attained after a cooling from a high temperature phase. We have found that, for high enough $`K_S/K_C`$ ratios, these configurations change from an AF state in which spins in each sublattice are almost aligned along the core easy-axis to a hedgehog-like structure with spins pointing along radial direction and AF order among sublattices (not shown here due to space limitations ). These structures are not observed when uniaxial surface anisotropy is considered and therefore can be attributed to the increased surface anisotropy.

Zero temperature hysteresis loops for a spherical particle with diameter $`D=3a`$ and an ellipsoidal particle with the same diameter and major axis $`L=8a`$ with $`K_C=1`$ and $`K_S=10,50`$ are compared in Figs. 1 and 2. In these figures, we also show the sum of projections surface and core magnetizations into the respective anisotropy directions $`M_n=_{i=1}^N\left|\stackrel{}{S}_i\widehat{n}_i\right|`$.

First, we notice that as particles become more elongated, although the coercive field and remanence magnetization remain almost constant, the hysteresis loops become more elongated ressembling those observed experimentally. Also the reversal of the core (dashed line) becomes less uniform, a fact that can be attributed to the increasing proportion of spins with anisotropy direction far from the field direction when increasing $`L`$. As $`K_S`$ increases from $`10`$ to $`50`$ the coercive field $`h_c`$ increases and, more importantly, there is a change in the magnetization reversal mechanism both for spherical and elliptical particles. Whereas in the first case the particle core reverses in a quasi-uniform form with the spins pointing mostly along the z-axis ($`M_n1`$) during the process except near $`h_c`$ and with the surface spins following the core reversal (the departure ). However, at higher $`K_S`$ \[(b) panels\], surface spins remain close to the local radial direction ($`M_n1`$) during all the reversal process, driving the core spins out of their local easy axis and making their reversal non-uniform ($`0.5<M_n<1`$) due to the appearance of the hedgehog-like structures mentioned above.

We acknowledge CESCA and CEPBA under coordination of C4 for the computer facilities. This work has been supported by SEEUID through project MAT2000-0858 and CIRIT under project 2001SGR00066.


text:
# 1 Fermion masses in MeV.

Mass and mixing of fermions from a rotating mass matrix111Updated version of a talk given by the first author at the International Symposium on Frontiers of Science - In Celebration of the 80th Birthday of C N Yang, June 2002, Beijing.

TSOU Sheung Tsun

and

Jonathan S. PALMER

Mathematical Institute, Oxford University

24–29 St. Giles’, Oxford OX1 3LB

United Kingdom

tsou @ maths.ox.ac.uk palmerj @ maths.ox.ac.uk

Abstract

We summarize the results showing that all existing data on mixing between up and down fermion states (i.e. CKM matrix and neutrino oscillations) and on the hierarchical quark and lepton mass ratios between generations are consistent with the two phenomena being both consequences of a mass matrix rotating in generation space with changing energy scale. Thus the rotation of the mass matrix can be traced over some 14 orders of magnitude in energy from the mass scale of the $`t`$-quark at 175 GeV to below that of the atmospheric neutrino at 0.05 eV.

By a rotating mass matrix we mean that the mass matrix of a fermion species, for example any one of $`U`$ the up quarks ($`t,c,u`$), $`D`$ the down quarks ($`b,s,c`$), $`L`$ the charged leptons ($`\tau ,\mu ,e`$) and $`N`$ the neutrinos ($`\nu _\tau ,\nu _\mu ,\nu _e`$), has not only scale-dependent eigenvalues but also scale-dependent eigenvectors. This last property means that the eigenvectors change orientation (or rotate) in generation space.

This is not just a hypothetical situation, because already in the Standard Model the renormalization group equation for the $`U`$ mass matrix has a term

$$16\pi ^2\frac{dU}{dt}=\frac{3}{2}DD^{}U+\mathrm{},$$
(1)
where the factor $`D^{}U=V`$ is the CKM quark mixing matrix. Now this matrix is experimentally found to be non-diagonal, meaning that even if $`U`$ is diagonal at a certain scale, it will become nondiagonal as a result of running to a different scale, at which when it is re-diagonalized the eigenvector will have rotated from its previous orientation. However, given that the off-diagonal elements of the CKM matrix are small, this effect can usually be neglected in practice in the Standard Model. This may very well not be the case in extensions of the Standard Model, where other forces may contribute to the mass matrix rotation to an observable extent.

In this note we shall turn the argument around and seek experimental evidence of a rotating mass matrix, in the manner we shall now detail.

There are two well established experimental facts which the Standard Model does not explain, namely fermion mass hierarchy and fermion mass mixing. The masses of the up and down quarks, and the charged fermions differ by more than 1–2 orders of magnitude as one goes from one generation to the next , as can be seen in Table 1.

The absolute values of the various CKM matrix elements have been measured to a good accuracy :

$$\left(\begin{array}{ccc}V_{ud}\hfill & V_{us}\hfill & V_{ub}\hfill \\ V_{cd}\hfill & V_{cs}\hfill & V_{cb}\hfill \\ V_{td}\hfill & V_{ts}\hfill & V_{tb}\hfill \end{array}\right)=\left(\begin{array}{ccc}0.97420.9757\hfill & 0.2190.226\hfill & 0.0020.005\hfill \\ 0.2190.225\hfill & 0.97340.9749\hfill & 0.0370.043\hfill \\ 0.0040.014\hfill & 0.0350.043\hfill & 0.99900.9993\hfill \end{array}\right)$$

The absolute values of the leptonic MNS mixing matrix, on the other hand, are not so well determined at present, but from neutrino oscillation experiments  we know that the mixing is in general considerably larger than the corresponding CKM mixing. The experimental estimates are:

$$\left(\begin{array}{ccc}U_{e1}\hfill & U_{e2}\hfill & U_{e3}\hfill \\ U_{\mu 1}\hfill & U_{\mu 2}\hfill & U_{\mu 3}\hfill \\ U_{\tau 1}\hfill & U_{\tau 2}\hfill & U_{\tau 3}\hfill \end{array}\right)=\left(\begin{array}{ccc}& 0.40.7& 0.00.15\\ & & 0.560.83\\ & & \end{array}\right)$$

Since the Standard Model gives no explanation for these two very remarkable experimental facts, one is tempted to seek other origins for them, and if one could find a single mechanism which could give rise to both, it would certainly provide a very interesting and promising avenue of study .

We now make the following rotating hypothesis:

> The rotating mass matrix gives rise to all observed mixing and lower generation masses.

We note as a side remark that this hypothesis is made and implemented in the Dualized Standard Model , from which we obtain good agreement with experiment222Specifically, with 3 real parameters fitted to data, it gives correctly to within present experimental bounds the following measured quantities: the mass ratios $`m_c/m_t`$, $`m_s/m_b`$, $`m_\mu /m_\tau `$, all 9 elements $`|V_{rs}|`$ of the CKM matrix, plus the 2 elements $`|U_{\mu 3}|`$ and $`|U_{e3}|`$ of the MNS matrix. It gives further by interpolation sensible though inaccurate estimates for the following: the mass ratios $`m_u/m_t`$, $`m_d/m_b`$, $`m_e/m_\tau `$ and the solar neutrino angle $`U_{e2}`$. Moreover, numerous detailed predictions have been made, and tested against data, in flavour-violation effects over a wide area comprsising meson mass differences, rare hadron decays, $`e^+e^{}`$ collisions, and muon-electron conversion in nuclei, with further predictions on effects as far apart in energy as neutrinoless double-beta decays in nuclei and cosmic ray air showers beyond the GZK cut-off of $`10^{20}`$ eV at the extreme end of the present observable energy range..

So we start with a mass matrix which is of rank 1, and which remains factorized on running, so that at any scale we have only one non-zero eigenvalue $`m_T`$:

$$m=m_T\left(\begin{array}{c}\xi \\ \eta \\ \zeta \end{array}\right)(\xi \eta \zeta )=m_T\left(\begin{array}{ccc}\xi ^2& \xi \eta & \xi \zeta \\ \xi \eta & \eta ^2& \eta \zeta \\ \xi \zeta & \eta \zeta & \zeta ^2\end{array}\right),$$
(2)
Thus the whole mass matrix can be encapsulated by one single normalized rotating vector $`v=(\xi ,\eta ,\zeta )`$ (without loss of generality we may assume $`\xi \eta \zeta `$). Note that the eigenvalue $`m_T`$ depends on the fermion species $`T=U,D,L,N`$.

Care has to be taken in defining physical masses and states when the mass matrix runs with scale. In our case when the mass matrix is given by a single vector we can adopt the following procedure  which guarantees that:

* lepton state vectors are always orthogonal,
* mixing matrix is always unitary,
* mass hierarchy is automatic.

First we run $`m`$ to a scale $`\mu `$ such that $`\mu =m_T(\mu )`$; this value we can reasonably call the mass of the highest generation. To fix ideas, let us concentrate on the $`U`$ type quarks $`t,c,u`$. The corresponding eigenvector $`v_t`$ is then the state vector of the $`t`$ quark. Having fixed this, we now know that the $`c`$ and $`u`$ lie in the 2-dimensional subspace $`V`$ orthogonal to $`v_t`$. As we go down in scale the eigenvector $`v`$ rotates so that $`V`$ is no longer the null eigenspace, and the projection of $`m`$ onto $`V`$ is a $`2\times 2`$ matrix, again of rank 1. We now repeat the procedure for this submatrix and determine both the mass and the state of the $`c`$ quark (Figure 1).

Once the $`c`$ state is determined, we know the $`u`$ state as well, as being the third vector of the orthonormal triad $`(t,c,u)`$. The $`u`$ mass is similarly determined.

The above procedure can be repeated for the other three types of fermions333Neutrinos will need further special treatment which we shall omit here. Instead see reference .. The direction cosines of the two triads $`(t,c,u)`$ and $`(b,s,d)`$, in other words, the 9 inner products, will give us the CKM matrix, Figure 2.

Similarly for the MNS mixing matrix of the leptons.

The above procedure takes a much simpler form when we have only two generations, and we shall in fact study this case first, as a planar approximation for the full three generation case, to better illustrate the methodology and the results .

Consider the up quarks first. When applied to this simple case, the procedure detailed above for defining masses and state vectors of flavour states gives the state vector $`v_t`$ of $`t`$ as the single massive eigenstate $`v`$ of the $`U`$-quark mass matrix at the scale $`\mu =m_t`$, and the vector $`v_c`$ as a vector orthogonal to $`v_t`$, as depicted in Figure 3.

But this should not be interpreted to mean that $`c`$ has zero mass, because $`m_c`$ is given by the nonzero eigenvalue of the mass submatrix, now only $`1\times 1`$, and is hence just the expectation value of $`m`$ in the state vector $`v_c`$, namely

$$m_c=m_t\mathrm{sin}^2\theta _{tc}0,$$
(3)
where $`\theta _{tc}`$ is the rotation angle between the scales $`\mu =m_t`$ and $`\mu =m_c`$. Similarly, although the mass matrices of the $`U`$ and $`D`$ quarks are aligned in orientation at all scales, one sees from

Figure 4 that by virtue of the rotation of the vector $`v`$ from the scale $`\mu =m_t`$ to the scale $`\mu =m_b`$ where the state vectors $`v_t`$ and $`v_b`$ are respectively defined, the two state vectors will not point in the same direction, having rotated by an angle $`\theta _{tb}`$ between the two scales:

$$V_{tb}=𝐯_t.𝐯_b=\mathrm{cos}\theta _{tb}1.$$
(4)
Other fermion species are similar.

Using (3) and (4), and Figures 3 and 4, and similar ones for the other fermions, we obtain the following relations:

$$V_{tb}=\mathrm{cos}\theta _{tb},|V_{ts}|=|V_{cb}|=\mathrm{sin}\theta _{tb},$$
(5)
from the mixing, and from the masses

$$m_c/m_t=\mathrm{sin}^2\theta _{tc},m_s/m_b=\mathrm{sin}^2\theta _{bs},m_\mu /m_\tau =\mathrm{sin}^2\theta _{\tau \mu }.$$
(6)
Inputting from data  on mixing

$$|V_{tb}|=0.99900.9993|V_{ts}|=0.0350.043,|V_{cb}|=0.0370.043,$$
(7)
and on masses for up quarks

$$m_t=174.3\pm 5.1\mathrm{GeV},m_c=1.151.35\mathrm{GeV},$$
(8)
for down quarks

$$m_b=4.04.4\mathrm{GeV},m_s=75170\mathrm{MeV},$$
(9)
and for charged leptons

$$m_\tau =1.777\mathrm{GeV},m_\mu =105.66\mathrm{MeV};$$
(10)
we obtain compatible values for

$$\theta _{tb}=0.03740.0447,0.03500.0430,0.03700.0430$$
(11)
and

$$\theta _{tc}=0.08010.0894,\theta _{bs}=0.13090.2076,\theta _{\tau \mu }=0.2463.$$
(12)
These values suggest a smooth curve, which will be the case if these angles are swept out by a rotating vector $`v`$. This is clearly borne out by Figure 5,

where the dotted curve represents the best-fit to data using MINUIT:

$$\theta =\mathrm{exp}(2.2670.509\mathrm{ln}\mu )0.0075$$
(13)
$`\mu `$ in GeV, with an excellent $`\chi ^2`$ of 0.21 per degree of freedom. The solid curve shown is an earlier calculation  to 1-loop in the Dualized Standard Model (see Footnote 1). The two curves are almost indistinguishable, and hence fit the data points equally well.

The positive result from the two-generation case is encouraging, but below the scale of roughly the $`s`$ quark mass, nonplanar effects begin to be appreciable, as can be estimated by the square of the Cabibbo angle which gives about 4 %. Therefore to go further down the scale it is necessary to study the full three-generation case . Writing

$$v(\mu )=(\xi (\mu ),\eta (\mu ),\zeta (\mu )),$$
(14)
we shall present the extracted angles by plotting $`\eta (\mu ),\zeta (\mu )`$ against scale $`\mu `$.

First we fix the $`U`$ triad to be

$$U:(1,0,0),(0,1,0),(0,0,1).$$
(15)
Then the $`D`$ triad is just given by the elements of the CKM matrix:

$$D:(V_{tb},V_{cb},V_{ub}),(V_{ts},V_{cs},V_{us}),(V_{td},V_{cd},V_{ud}).$$
(16)
We thus obtain the points $`t`$ and $`b`$ on the plot (Figure 6).

Next we consider $`v(\mu =m_c)`$. Since the mass $`m_c`$ comes from ‘leakage’ from $`m_t`$, we can write

$`v(\mu =m_c)`$ $`=`$ $`\mathrm{cos}\theta _{tc}v_t+\mathrm{sin}\theta _{tc}v_c`$ (17)
$`=`$ $`\sqrt{1m_c/m_t}v_t+\sqrt{m_c/m_t}v_c,`$
which gives the point $`c`$ in Figure 6. Note that the existence of another disjoint branch of the square root does not affect the question of interest to us here, namely, whether the allowed region is consistent with data lying on a smooth rotation curve, so long as the first branch already does, and can therefore be ignored.

Similar considerations give us a line as the allowed region for the $`s`$ quark. We shall leave open the question of the light quarks $`u`$ and $`d`$, since they are tightly bound and both the definition and the value of their masses are quite uncertain.

In principle the leptons are not linked to the quarks, but the planar approximation above (Figure 5) and the DSM calculations  both suggest that they can be put on the same trajectory as the quarks. In fact, by putting the $`\tau `$ point by interpolation between the $`b`$ and $`c`$ points, we can easily continue the rotation curve with the leptons, as is again evident from Figure 6. In Figure 7 we show the same curve with its three projections onto the three planes: $`\mu \eta `$, $`\mu \zeta `$, and $`\eta \zeta `$.

We have not given in any detail on how the actual angles are extracted from the quark and lepton data, only the theory of it. Full details can be found in .

In conclusion we can say that, using the rotation hypothesis, and all the fermion mass and mixing data (apart from $`m_u`$ and $`m_d`$), we exhibit a smooth curve in 3-d space, passing through all the allowed regions, and spanning some 14 orders of magnitude in energy. Moreover, this curve coincides with the DSM curve  where we expect it to do, that is, where the 1-loop approximation is good.

Most of the work reported here was done in collaboration with José Bordes and Chan Hong-Mo.


text:
# General Relativity and Gravitation — Millisecond oscillators in accreting neutron stars and black holes —

## 1 Millisecond variability in LMXBs

Frequencies as large as 1.2 kHz in the X-ray flux of LMXBs have been detected in 1996 by the first large satellite detector to have the requisite time resolution, the Rossi X-ray Timing Explorer. For a review of observations see van der Klis 2000, and for that of the previously discovered lower frequency variability, van der Klis 1989.

Millisecond variability of X-ray emission from accreting neutron stars and black holes has been expected at least from the time of Shvartsman’s 1971 publication. Bath (1973) pointed out that in cataclysmic variables a natural timescale of variation is given by the Keplerian frequency and that clumps orbiting in the accretion disk could modulate the flux. For stellar-mass black holes and for neutron stars such clumps would induce millisecond variability of the X-ray flux (Sunyaev 1973). After quasi-periodic signals (QPOs) had been observed at frequencies corresponding to Keplerian motion in the inner parts of the accretion disk in a strongly magnetized accreting neutron star, the X-ray pulsar EXO 2030+375 (Angelini et al. 1989), Kluźniak, Michelson and Wagoner (1990) pointed out that scaling to the case of a moderately rotating compact neutron star of negligible magnetic field would yield a QPO at the frequency of orbital motion in the marginally stable orbit, $`f_K=2.2\mathrm{kHz}(1+0.75j)M_{}/M`$ (with $`j`$ the dimensionless angular momentum of the star), and that this frequency, when observed, could be used to obtain a direct estimate of the stellar mass $`M`$.

Another line of reasoning originates with Okazaki, Kato and Fukue (1987), who pointed out that the radial epicyclic frequency in the space time metric around a black hole (or a neutron star) has a maximum, and goes to zero in the marginally stable orbit, which allows certain modes of vibration to be trapped in the inner region of the accretion disk. Diskoseismology, the theory of vibrations of accretion disks, has been developed in detail in the following years and has been recently reviewed by Wagoner (1999) and Kato (2001). In particular, it predicts the appearance of certain stable frequencies of quasi-periodic motion in black-hole accretion disks. The fundamental g-mode and c-mode, especially, have frequencies which are a function of the black-hole mass and the Kerr parameter $`a`$, alone, and which agree well with those observed in the micro-quasars.

The kHz QPOs observed in accreting neutron stars are more difficult to understand, as the observed frequencies vary considerably on a timescale of minutes or hours. Further, typically, two high frequency peaks in the power spectra are observed, where only one had at first been reported in black holes. It has been thought, initially, that the separation between the two peaks corresponds to the spin frequency of the neutron star but, in fact, this difference frequency is not constant either (e.g., van der Klis et al. 1997), whereas the spin frequency can vary only on the accretion time-scale, which is millions of years in LMXBs (even in the brightest sources the accretion rate is less that $`10^7M_{}`$/y, and the neutron star mass is on the order of $`1M_{}`$).

In spite of these differences, the power spectra of both the neutron star sources and the black hole ones are quite similar, the disparity in the characteristic frequencies being easy to explain by the higher mass of the black holes. Also, in the neutron star sources, the maximum frequency observed varies remarkably little with the accretion rate, or the type of the source. All this points to an accretion disk origin of the phenomenon.

## 2 The accretion disk in LMXBs

There is little doubt that mass transfer from the binary companion onto a neutron star or a black hole is responsible for the observed X-ray emission in low-mass X-ray binaries (Shklovsky 1967, Zel’dovich and Guseynov 1966). Although the central region around the compact object is too small to be resolved with current instruments (being only several tens of kilometers across, observed from a distance of kiloparsecs) it is generally agreed that the accreting fluid forms a disk around the central compact object (Prendergast and Burbidge 1968). Whether this is a geometrically thin, and hence rotation-supported, accretion disk (Shakura and Sunyaev 1973), or a sub-Keplerian, geometrically thick disk (Jaroszyński et al. 1980), is a matter of some dispute. For an artist’s impression of the general appearance of the system see, e.g., the cover of the book by Kato, Fukue and Mineshige (1998).

There are two analogous types of systems in which the disk can be resolved, and in these the disk is geometrically thin. These are cataclysmic variables (binaries in which the accreting object is a non-magnetic white dwarf), and the celebrated extragalactic quasars. The correspondence between the physical processes in quasars and the Galactic X-ray novae, dubbed micro-quasars, has already been remarked upon in the literature (e.g., Mirabel and Rodriguez 1998), the main difference between the real and the “micro” quasars being in the mass of the black hole: millions or billions of solar masses as opposed to about ten $`M_{}`$.

Water-vapour masers have been observed in some of the quasars, and radio observations allowed the mapping of the geometrical distribution of the emitters, as well as of their velocity distribution. The results leave no doubt—in NGC 4258 the masers are present in a geometrically thin and warped disk, and they orbit the central source with Keplerian velocities (Miyoshi et al. 1995). Perhaps the accretion disks in LMXBs are also Keplerian, although most of what we have to say here will apply equally well to some accretion tori.

## 3 Two variable kHz frequencies: non-linear resonance?

If a common framework can be found for understanding the rapid variability of the X-ray flux in both the black hole and the neutron star sources in LMXBs, the behavior of QPOs in the latter may offer the key to unlocking the secrets of the former.

The most characteristic property of kHz QPOs in neutron stars is that there are typically two peaks in the power spectrum, located at frequencies which vary in time, but only within a limited, characteristic range (Fig. 2). A single variable frequency can be identified already in the Newtonian regime. It is also easy to find two characteristic constant frequencies. But we would like to ask what may be the origin of two, characteristic, but variable frequencies, in nearly Keplerian motion around a neutron star.

In Newtonian gravity of a spherically symmetric point mass, there is no preferred frequency (Fig. 3). Further, there is degeneracy between the orbital frequency in circular and not-quite-circular motion—this is why bound orbits are closed curves (Kepler’s first law). So if neutron star exteriors were described by this theory, one would expect temporal fluctuations on all time-scales in the accretion disk, corresponding to the (huge) range of radii covered by the accreting fluid, but no characteristic frequencies.

However, characteristic frequencies can be induced in a Keplerian distribution externally—the gaps in the rings of Saturn are a well-known example. The radius of a neutron star introduces a characteristic length-scale to the problem, and a corresponding frequency, the Keplerian frequency at the surface of the star. This frequency is fixed for a given star, of course, and together with the spin frequency (which, in LMXBS, is practically constant on a human time-scale) allows two fixed frequencies to be identified in the Newtonian problem of accretion onto a neutron star.

If the accretion disk is terminated by the magnetosphere of the neutron star, its inner radius and the corresponding orbital frequency can vary somewhat with the accretion rate. This would lead to the appearance of a single characteristic, and variable, frequency, in addition to the constant stellar spin frequency. As an example take the X-ray pulsar EXO 2030+375, already mentioned in Section 1: the power spectrum of its X-ray flux exhibits prominent harmonics of the rotational frequency of the neutron star, in addition to the $`20`$ mHz QPO with frequency varying in step with the luminosity (Angelini et al. 1989).

In general relativity orbital degeneracy is broken, bound orbits are no longer closed curves (“precession of the perihelion”), the frequency of epicyclic motion induced by a radial perturbation differs from orbital frequency. There is a preferred length-scale, the gravitational radius $`M(G/c^2)`$, and corresponding characteristic frequencies: the radial epicyclic frequency goes to zero at the innermost (marginally) bound circular orbit, and has a maximum at a slightly larger orbit (Fig. 4). However, these are fixed frequencies.

In any sufficiently complicated system, resonances appear at characteristic frequencies. Of particular interest are resonances in non-linear systems, these occur not only when the driving frequency is close to the eigenfrequency of the oscillator, but also when the two frequencies are (nearly) in a rational ratio, $`m:n`$, with the strongest resonances occurring at the lowest values of the integers (Landau and Lifshitz 1976). The appearance of two, somewhat variable, frequencies in the rather complex accreting flow in LMXBs is strongly suggestive of a non-linear resonance.

## 4 Epicyclic resonance?

We note that there are two types of fluid perturbations in nearly circular three-dimensional flow which correspond to the two epicyclic frequencies familiar from study of test-particle orbits—the radial epicyclic frequency of planar motion, and the meridional epicyclic frequency of slightly non-planar motion. For a spherically symmetric gravitating body, the meridional epicyclic frequency is equal to the orbital frequency of the unperturbed circular orbit. For co-rotating orbits in the Kerr metric the meridional epicyclic frequency is always lower than the orbital frequency and higher than the radial epicyclic frequency (e.g., Wagoner 1999), and a similar relation obtains for rapidly rotating neutron stars. Thus, at characteristic radii, the ratio of the two epicyclic frequencies is rational, allowing the possibility of a non-linear resonance between the two types of fluid perturbations. We have suggested that this is the origin of the two kHz QPOs in neutron stars (Kluźniak and Abramowicz 2001a).

What would be the implications of this suggestion for the QPOs observed in black hole systems? When we first argued for the presence of a non-linear epicyclic resonance in LMXBs, only a single high-frequency QPO had been observed in black hole systems. The most obvious prediction following from our hypothesis was that a second QPO should be present also in black holes, and that the two frequencies should be in a rational ratio.

We are greatly encouraged that a second QPO has indeed been discovered in GRO J1655-40 (which reportedly has a $`6M_{}`$ black hole) and that the two simultaneously observed frequencies have been reported to be 450 Hz and 300 Hz (Strohmayer 2001a), clearly in a 3:2 ratio (Abramowicz and Kluźniak 2001). Another system, GRO J 1550-564 with a $`10M_{}`$ black hole, in one of its states also shows two QPOs at the same frequency ratio, and the frequencies of 270 Hz and 180 Hz can be obtained by mass scaling from those of 1655-40, implying that the Kerr parameter of the two holes may be nearly identical (Remillard et al. 2002). Fig. 5 (reproduced from Remillard et al. 2002) illustrates how these ideas may be used to constrain the spin of the black hole, when the mass is known and the observed frequencies are identified with linear combinations of the orbital frequency and the radial epicyclic frequency.

## 5 Parametric epicyclic resonance?

We have seen that the frequencies of these first twin QPOs to be discovered in black hole systems are in a 3:2 ratio. Why not 1:2 or 1:3 ? There is a particular type of non-linear resonance in which large-amplitude motions are easily excited by a small perturbation of the right frequency. This is parametric resonance, in which the eigenfrequency of an oscillator is varied periodically, i.e., when the eigenfrequency of a system is of the form $`\omega =\omega _0(1+h\mathrm{cos}\omega _1t)`$, with $`h<<1`$. The displacement of a harmonic oscillator with this eigenfrequency satisfies the Mathieu equation. Resonance occurs when

$$\omega _1/\omega _0=2/n,$$

with $`n`$ integer; the amplitude grows intially exponentially to a saturation value determined by non-linearities, the fastest growing modes correspond to the lowest values of the integer $`n`$ (Landau and Lifshitz 1976). For resonance to occur in the presence of dissipation, the size $`h`$ of the perturbation has to exceed a certain threshold value.

One can write down the equations of disk perturbation in a form in which the the vertical and radial components of motion are coupled by a “pressure” term proportional to the speed of sound squared. Typically, in a linear analysis, the term can be neglected. But when a state of motion arises in which this coupling is modulated at a frequency obeying the resonant condition, this term becomes the most important one. Specifically, there are states of disk motion in which the dominant frequency is close to the radial epicyclic frequency $`\omega _r`$, and other states of “vertical” motion, in which the frequency is close to the meridional epicyclic frequency $`\omega _z`$. As we have seen, $`\omega _r<\omega _z`$ for black holes, always. Identifying $`\omega _r`$ with $`\omega _1`$ in the Mathieu equation, and $`\omega _z`$ with $`\omega _0`$, we see that for parametric resonance to occur in the Kerr metric, both

$$\omega _1<\omega _0$$

and

$$\omega _1/\omega _0=2/n$$

must hold. The lowest value of $`n`$ for which both eqs. are satisfied is $`n=3`$. Hence, if 180 Hz is the value of the radial epicyclic frequency, vertical motions at a frequency of 270 Hz are easily excited in the disk.

Incidentally, it is often assumed that accretion disks in these systems have a hot corona. The observed energy dependence of the QPOs in GRO J1655-40 (Fig. 6) is easy to square with the notion that motions perpendicular to the disk plane, i.e., the ones most likely to disturb the corona, occur predominantly at the higher frequency.

## 6 Frequencies in the non-linear coupled system

What is the expected harmonic content of the epicyclic oscillators in resonance? First let us examine the non-linearities of the effective potential. If the unperturbed disk is symmetric with respect to reflections in the equatorial plane, in agreement with the symmetry of the black hole, the anharmonicity of the meridional oscillator is of third order in the equations of motion. In contrast, the equation of motion of the radial oscillator has second order anharmonic terms (the potential has cubic terms). This implies that subharmonic content is possible in the motion of these oscillators (Landau and Lifshitz 1976). The meridional oscillator may have a subharmonic at $`\omega =(1/3)\omega _0`$ and the radial oscillator a subharmonic at $`\omega =(1/2)\omega _1`$. If the two oscillators are in a 2:3 resonance, $`\omega _1=(2/3)\omega _0`$, and the two subharmonics coincide!

It is noteworthy, that Remillard et al. (2002) report a third QPO frequency in one of the states of XTE J1550-564. In addition to 270 Hz and 180 Hz, a weaker oscillation at 90 Hz is present (Fig. 5). The third frequency coincides with the predicted subharmonic.

## 7 1:2:3:5

If the two oscillators are in parametric resonance, as we have argued above, the frequency of 90 Hz in XTE J1550-564 is predicted for a more direct reason. The Mathieu equation does not admit strictly harmonic solutions. Instead, if the leading term in the solution is harmonic at frequency $`\omega _0`$, additional terms must be included at “combination” frequencies, $`\omega _0+\omega _1`$ and $`\omega _0\omega _1`$. For a 2:3 resonance these are at $`(1/3)\omega _0`$ and $`(5/3)\omega _0`$.

It is interesting to note that of these two combination frequencies one has been detected in XTE J1550-564, while the other may have been detected in GRS 1915+105. Why this is the case may be difficult to understand. Experience with asteroseismology shows that it is by far easier to compute the spectrum of allowed frequencies than to predict which of the frequencies will actually show up in the source at any given time. Primarily this is because energy flow in the system is very sensitive to the detailed, and poorly understood, coupling of modes in the non-linear system. These difficulties are exacerbated in accretion disk theory, which can offer only educated guesses as to the physical state of the accreting fluid.

The state of affairs is this. The theory of parametric resonance predicts a sequence of initial primes for the frequencies, 1:2:3:5, with 2 and 3 the most prominent. The 2:3 have been detected in two sources, GRO 1655-40 (300 Hz : 450 Hz) and XTE J1550-564 (180 Hz : 270 Hz). The “1” (90 Hz) also in XTE J1550-564, 1:2:3 really (Fig. 7). In a third system, 1915+105, the frequencies are 41.5 Hz and 69.2 Hz (Strohmayer 2001b), i.e., they are in a 3:5 ratio (Kluźniak and Abramowicz 2002).

## 8 Neutron star systems

The observed ratios of QPO frequencies in black hole systems agree very well with the theory of non-linear resonance, and parametric resonance between epicyclic oscillators in particular. We have been led to this theory by a search for the origin of two variable QPO frequencies in neutron star systems, and we view the subsequent discovery in black hole systems of double frequencies in rational ratios as confirmation of this theory. However, in neutron stars, the ratio of the two frequencies is not fixed (e.g., Fig. 2). Is the kHz QPO phenomenon in neutron stars and in black holes really a manifestation of the same physics? If so, then a 2:3 ratio should also be discernible in LMXBs containing neutron stars.

Let us begin by examining anecdotal evidence. In the van der Klis 2000 review, only two examples are shown of twin peaks in the power spectra. We have reproduced the first one, showing the kHz QPOs of Sco X-1, as Fig. 1. The second one, showing the twin QPOs of the neutron star source 4U 1608-52 (Méndez et al. 1998), is reproduced here as Fig. 8. Note that in both cases, the frequencies of the QPOs are close to 600 Hz and 900 Hz, i.e., close to a 2:3 ratio. This could be coincidental, or it could reflect the fact that certain frequencies are more likely to appear than others. Abramowicz et al. (2002b) argue that the latter is the case, at least for Sco X-1.

For Sco X-1, the two frequencies seem to be scattered over a wide range (Fig. 2), but the distribution of points along the line of correlation is not random (Abramowicz et al. 2002b). One way to demonstrate this is to plot a histogram of the frequency ratios (Fig. 9). It is striking that the observations cluster about the value 2/3. We take this to be a very strong indication that a resonance is involved in the simultaneous appearance of the two frequencies. A corotation resonance has been considered (Kato 2002), but it seems to be damped (Kato 2003a). Kato 2003b has shown that a parametric resonance may excite both an axisymmetric and a non-axisymmetric disk (one-armed) modes whose frequencies are in a $`1/\sqrt{2}0.71`$ ratio for the Schwarzschild metric.

Abramowicz et al. (2003) have investigated the parametric epicyclic resonance for nearly circular and nearly geodesic motion. Choosing an arbitrary, but fairly general form of coupling between the radial and vertical components of motion, they have demonstrated the presence of parametric resonance, at frequencies somewhat depending on the strength of coupling, which is taken to simulate pressure effects in the disk. In resonance, the ratio of the frequencies is accurately 2:3, but sizable amplitudes of motion are excited also slightly off-resonance. Quantifying with a single adjustable parameter how far off resonance the system is, the authors have reproduced the slope of the frequency correlation found in Sco X-1 (Fig. 10).


text:
# Exact, explicit and entire solutions to a nontrivial finite-difference equation and their applications

## I Introduction

Explicit, exact and entire solutions to a finite-difference equation

$$u(z\delta )/2u(z+\delta )/2\delta ^2u(z)/z=\lambda ^{(\delta )}u(z)$$
(1)
under the condition $`u(0)=0`$ have been found. Alternatively, the functions $`u(z)`$ can be considered as solutions of the differential equation of infinite order

$$\left(\mathrm{cosh}(\delta d/dz)\delta ^2/z\right)u(z)=\lambda ^{(\delta )}u(z).$$
(2)
The solutions become uniquely defined if we require that they coincide with the solutions of the limiting equation

$$u^{\prime \prime }(z)/2u(z)/z=[(\lambda ^{(\delta )}1)/\delta ^2]u(z)$$
(3)
in the limit $`\delta 0`$. On the positive real axis this equation is the radial Schrödinger equation for the ($`l=0`$)-states of the hydrogen atom.

By using the power series ansatz introduced in \[M. Aunola, J. Math Phys 44, 1913 (2003); math-ph/0207038\] we look for solutions of

$$\underset{m=0}{\overset{\mathrm{}}{}}\frac{\delta ^{2m}}{(2m)!}\frac{d^{2m}u}{dz^{2m}}\frac{\delta ^2u}{z}=\lambda ^{(\delta )}u(z)$$
(4)
and obtain the leading terms of the solution and eigenvalues. By reinterpreting the different parts as analytical functions we find the general form of the solution as

$$u_n^{(\delta )}(z)=\left(\underset{k=1}{\overset{n}{}}L_k^{(n)}\alpha _k^{(n,\delta )}z^k\right)\mathrm{exp}(z\mathrm{arsinh}(\delta /n)/\delta ),$$
(5)
where $`\mathrm{arsinh}`$ is the inverse of the hyperbolic sine, $`L_k^{(n)}`$ are the coefficients of the associated Laguerre polynomial of the first kind and the corresponding eigenvalues are given by

$$\lambda _n^{(\delta )}=\sqrt{1+\delta ^2/n^2}.$$
(6)
The power series expansions converge for $`|\delta |<n`$, but they also apply to all finite, real and positive values of $`\delta `$.

Simultaneously we diagonalise the infinite matrices $`V`$ with elements

$$V_{jj}=\delta /j,V_{j,j+1}=V_{j+1,j}=1/2,j=1,2,\mathrm{}.$$
(7)
The $`j^{\mathrm{th}}`$ component of eigenvector corresponds to $`u_j:=u(\delta j)`$. In the matrix representation the condition $`u_0=0`$ has been explicitly inserted. We can also construct isospectral matrices with simple eigenvalues, i.e. $`\{1/n^{2k}\}_{n=1}^{\mathrm{}}`$ and, subsequently, infinite matrices with desired eigenvalues. In addition, the results of \[V. Berezin, Phys. Rev. D 55, 2139 (1997)\], concerning quantum black holes are correctly reproduced. The discretised, associated Laguerre polynomials corresponding to $`L_n^1(x)=dL_n(x)/dx`$ are defined also.

For further details, see the source file Notebook/Exactsolu.nb containing all explicit formulae.

Acknowledgement: Discussions with Prof. J. Timonen, Prof. O. Civitarese and Dr. L. Kahanpää are gratefully acknowledged.


text:
# 1 Introduction

## 1 Introduction

The lightest supersymmetric particle (LSP) in the standard R-parity conserving supersymmetric model is the leading particle physics candidate for the dark matter (DM) of the universe . The most popular supersymmetry (SUSY) breaking model is the minimal supergravity (SUGRA) model having universal scalar, gaugino masses and trilinear couplings at the GUT scale. Over most of the parameter space of this model the LSP is dominantly a bino ($`\stackrel{~}{\chi }_1^0\stackrel{~}{B}`$) which does not couple to W or Z-boson. Hence they can only pair-annihilate via the exchange of superparticles like squarks or sleptons, $`\stackrel{~}{B}\stackrel{~}{B}\stackrel{\stackrel{~}{q}(\stackrel{~}{l})}{}q\overline{q}(l^+l^{})`$. The experimental limits on these particle masses, $`m_{\stackrel{~}{q}}\begin{array}{c}>\hfill \\ \hfill \end{array}200\mathrm{GeV}`$ and $`m_{\stackrel{~}{l}}\begin{array}{c}>\hfill \\ \hfill \end{array}100\mathrm{GeV}`$ , imply a rather slow rate of pair annihilation. Consequently, the model predicts an over-abundance of the DM relic density over most of the parameter space . This has led to several recent works, extending the SUSY DM investigations to nonminimal SUGRA models \[References-References\]. While many of them explore models with nonuniversal scalar masses, we shall concentrate here on nonuniversal gaugino mass models. In particular, we shall focus on a model leading generically to a Higgsino-like LSP. Because of its unsuppressed coupling to W and Z bosons the $`\stackrel{~}{H}\stackrel{~}{H}W^+W^{},ZZ`$ annihilation rates via s-channel Z-boson and t-channel Higgsino exchanges are large. Besides, there is a near degeneracy of the lighter neutralinos and lightest chargino masses in this case,

$`m_{\stackrel{~}{\chi }_1^0}m_{\stackrel{~}{\chi }_2^0}m_{\stackrel{~}{\chi }_1^\pm }|\mu |,`$ (1)
where $`\mu `$ is the supersymmetric Higgsino mass parameter. This leads to large coannihilation cross sections . Consequently, the Higgsino DM density falls below the cosmologically favoured range ,

$`0.05<\mathrm{\Omega }_mh^2<0.2,`$ (2)
where the lower limit corresponds to the galactic density of DM ($`\mathrm{\Omega }_m0.1`$) from rotation curves. Thus the Higgsino DM can only be a subdominant component of the galactic DM density. However, its large coupling to Z implies a large rate of capture inside the Sun. Hence, the model predicts a sizable indirect detection rate of Higgsino DM via high energy neutrinos coming from their pair annihilation in the solar core. This is much larger than the minimal SUGRA model prediction and should be detectable at the future neutrino telescopes, as shown below.

## 2 Non-Universal Gaugino Mass Model

SUGRA model with nonuniversal gaugino masses at the GUT scale have been discussed in many earlier works . We shall only quote the main results here, focusing on the SU(5) GUT. In this model the gauge kinetic function depends on a nonsinglet chiral superfield $`\mathrm{\Phi }`$, whose auxiliary $`F`$-component acquires a large vacuum expectation value (vev). Then the gaugino masses come from the following dimension five term in the Lagrangian:

$`L={\displaystyle \frac{<F_\mathrm{\Phi }>_{ij}}{M_{Planck}}}\lambda _i\lambda _j`$ (3)
where $`\lambda _{1,2,3}`$ are the $`U(1)`$, $`SU(2)`$ and $`SU(3)`$ gaugino fields i.e. the bino $`\stackrel{~}{B}`$, the wino $`\stackrel{~}{W}`$ and the gluino $`\stackrel{~}{g}`$ respectively. Since the gauginos belong to the adjoint representation of SU(5), $`\mathrm{\Phi }`$ and $`F_\mathrm{\Phi }`$ can belong to any of the irreducible representations appearing in their symmetric product, i.e.

$`(24\times 24)_{symm}=1+24+75+200`$ (4)
The minimal SUGRA model assumes $`\mathrm{\Phi }`$ to be a singlet, which implies equal gaugino masses at the GUT scale. On the other hand if $`\mathrm{\Phi }`$ belongs to one of the nonsinglet representations of $`SU(5)`$, then these gaugino masses are unequal but related to one another via the representation invariants. Thus the three gaugino masses at the GUT scale in a given representation $`n`$ are determined in terms of a single SUSY breaking mass parameter $`M_{1/2}`$ by

$`M_{1,2,3}^{G,n}=C_{1,2,3}^nM_{1/2}`$ (5)
where $`C_{1,2,3}^1=(1,1,1)`$, $`C_{1,2,3}^{24}=(1,3,2)`$, $`C_{1,2,3}^{75}=(5,3,1)`$ and $`C_{1,2,3}^{200}=(10,2,1)`$. The resulting ratios of $`M_i^G`$’s for each $`n`$ are listed in Table 1.

Of course in general the gauge kinetic function can involve several chiral superfields belonging to different representations of $`SU(5)`$ which gives us the freedom to vary mass ratios continuously. We shall explore such a possibility in a future work. But let us concentrate here on the representations 1, 24, 75 and 200 individually. While the singlet representation corresponds to universal gaugino masses, each of the nonsinglet representations corresponds to definite mass ratios and is therefore as predictive as the former.

These nonuniversal gaugino mass models are known to be consistent with the observed universality of the gauge couplings at the GUT scale \[ReferencesReferences\]

$`\alpha _3^G=\alpha _2^G=\alpha _1^G=\alpha ^G(1/25)`$ (6)
Since the gaugino masses evolve like the gauge couplings at one loop level of the renormalisation group equations (RGE), the three gaugino masses at the electroweak (EW) scale are proportional to the corresponding gauge couplings, i.e.

$`M_1`$ $`=`$ $`(\alpha _1/\alpha _G)M_1^G(25/60)C_1^nM_{1/2}`$
$`M_2`$ $`=`$ $`(\alpha _2/\alpha _G)M_2^G(25/30)C_2^nM_{1/2}`$
$`M_3`$ $`=`$ $`(\alpha _3/\alpha _G)M_3^G(25/9)C_3^nM_{1/2}`$ (7)
For simplicity we shall assume a universal SUSY breaking scalar mass $`m_0`$ at the GUT scale. Then the corresponding scalar masses at the EW scale are given by the renormalisation group evolution formulae . A very important SUSY breaking mass parameter at this scale is $`m_{H_2}`$ , as it appears in the EW symmetry breaking condition,

$`\mu ^2+M_Z^2/2={\displaystyle \frac{m_{H_1}^2m_{H_2}^2\mathrm{tan}^2\beta }{\mathrm{tan}^2\beta 1}}m_{H_2}^2,`$ (8)
where the last equality holds for the $`\mathrm{tan}\beta \begin{array}{c}>\hfill \\ \hfill \end{array}`$ 5 region, which is favoured by the Higgs mass limit from LEP . Expressing $`m_{H_2}^2`$ at the right hand side in terms of the GUT scale mass parameters at a representative value of $`\mathrm{tan}\beta =10`$ gives

$`\mu ^2+{\displaystyle \frac{1}{2}}M_Z^2`$ $``$ $`0.1m_0^2+2.1M_{3}^{G}{}_{}{}^{2}0.22M_{2}^{G}{}_{}{}^{2}0.006M_{1}^{G}{}_{}{}^{2}+0.006M_1^GM_2^G+`$ (9)
$`0.19M_2^GM_3^G+0.03M_1^GM_3^G,`$
neglecting the contribution from the trilinear coupling term at the GUT scale. Moreover, the coefficients vary rather mildly over the moderate $`\mathrm{tan}\beta `$ region. Although we shall use exact numerical solutions to the two-loop RGE, two points are worth noting from this simple equation.

Firstly, eq.(9) gives a measure of fine-tuning from the required degree of cancellation between the dominant terms $`\mu ^2`$ and $`M_{3}^{G}{}_{}{}^{2}`$ to give the right EW scale $`M_Z^2/2`$. The LEP limit on the lightest chargino mass, $`m_{\stackrel{~}{\chi }_1^\pm }>100\mathrm{GeV}`$  implies

$`|\mu |,|M_2|>100\mathrm{GeV}`$ (10)
while eq.(9) implies

$`\mu ^2+M_Z^2/22.1M_{1/2}^2(\mathrm{for}\mathrm{mSUGRA})\mathrm{and},`$
$`\mu ^2+M_Z^2/21.4M_{1/2}^2(\mathrm{for}\mathrm{n}=200\mathrm{model})`$ (11)
Thus for the universal gaugino mass case of mSUGRA eqs.(7), (10) and (11) imply fine-tuning at least at the level of $`10`$ . On the other hand one sees from these equations that the fine-tuning problem is significantly alleviated in the nonuniversal models with $`n=75`$ and 200, corresponding to $`C_2^n=3`$ and 2 respectively .

Secondly, the universal gaugino mass model corresponds to $`|\mu |>|M_2|>|M_1|`$, which implies that the lighter chargino and neutralinos are dominantly gauginos with hierarchical masses -i.e. $`m_{\chi _1^0}M_1`$ and $`m_{\chi _2^0,\chi _1^\pm }M_22M_1`$. There is however a narrow strip of very large $`m_0`$ region where the first term of eq.(9) pushes down $`|\mu |`$ towards the LEP limit as given in eq.(10). This is the so called focus point region , where the lighter chargino and neutralinos are mixed Higgsino-gaugino states. But, over the bulk of the parameter space the LSP is dominantly a $`\stackrel{~}{B}`$, which leads to an over-abundance of the DM relic density, as discussed earlier. One expects from Table 1 a similar result for the $`n=24`$ model. In fact it predicts a larger hierarchy between the $`\chi _1^0`$ and $`\chi _2^0(\chi _1^\pm )`$ masses, $`M_1`$ and $`M_2`$. Consequently the LSP is completely dominated by $`\stackrel{~}{B}`$ and can be relatively light. The SUSY DM phenomenology for this case has been discussed recently in Ref. . In contrast we see from eqs.(7), (11) and Table 1 that the $`n=200`$ model predict the opposite hierarchy $`|\mu |<|M_2|<|M_1|`$, while the $`n=75`$ model has $`|\mu |<|M_1|<|M_2|`$. Thus the lighter chargino and neutralino states are Higgsino dominated and roughly degenerate (eq.(1)) over the bulk of the parameter space for both $`n=75`$ and 200 models. This leads to a DM relic density somewhat below the cosmologically favoured range (eq.(2)). We should mention here that unlike the case of mSUGRA, there is no possibility of having stau coannihilations in the $`n=75`$ and the $`n=200`$ models. This is related to staus being significantly heavier than $`m_{\stackrel{~}{\chi }_1^0}`$ in these scenarios for all values of $`\mathrm{tan}\beta `$. This in fact originates from the gauge sector running of the slepton RGEs due to the specific gaugino mass nonuniversalities.

It is for the above reason that the Higgsino DM phenomenology of the nonuniversal SUGRA models, corresponding to $`n=75`$ and 200 have not been explored in detail so far. We feel that this is important for two reasons: (i) Even though the Higgsino may be a subdominant component of the galactic DM density, its large coupling to $`Z`$-boson implies large rate of capture inside the Sun. Hence the model predicts a sizable indirect detection rate via high energy neutrinos coming from their pair annihilation inside the sun. Even after rescaling by the low DM density factor the indirect detection rate comes out to be larger than the minimal SUGRA predictions . (ii) It is possible that the thermal relic density of the Higgsino DM is enhanced by either a modification of the freeze-out temperature of the standard cosmological model due to a quintessence field as suggested in Ref. or by nonthermal production mechanisms of the type suggested in Ref.. In that case it can be the dominant component of the galactic DM density. Therefore we have computed the indirect and direct detection rates both with and without rescaling.

## 3 Higgsino Dark Matter in the n=200 SUGRA Model

We shall concentrate on the nonuniversal SUGRA model corresponding to $`n=200`$ because it can generate radiative electroweak symmetry breaking (EWSB) over a much wider range of parameters compared to the $`n=75`$ case, the latter being restricted to have small $`\mathrm{tan}\beta `$ solutions only. Fig. 1 shows the allowed regions in the $`m_0M_{1/2}(=M_3^G)`$ plane for $`\mathrm{tan}\beta =`$5, 10, 30 and 50. The area marked I at the top is disallowed because $`\mu ^2`$ falls below the LEP limit (eq.10) and then becomes negative, signalling the absence of EWSB. The area marked II at the bottom is disallowed because the Higgs potential becomes unstable at the GUT scale. We have chosen $`\mu >0`$ since the $`\mu <0`$ branch is strongly disfavoured by the $`bs+\gamma `$ branching ratio, along with the muon anomalous magnetic moment ($`a_\mu `$) constraint. This figure shows that the $`bs+\gamma `$ constraint to be rather mild for $`\mu >0`$. The lower limit from $`a_\mu `$ (not shown) is even milder.

We see from the contours of $`\mu `$ in Fig. 1 that one generally has $`\mu <500`$ GeV in the model as anticipated. It also shows the gaugino component of the LSP,

$`Z_g=N_{11}^2+N_{12}^2,`$ (12)
where

$`\stackrel{~}{\chi }_1^0=N_{11}\stackrel{~}{B}+N_{12}\stackrel{~}{W}+N_{13}\stackrel{~}{H}_1+N_{14}\stackrel{~}{H}_2.`$ (13)
We see that the bulk of the allowed parameter space corresponds to $`Z_g<10\%`$, which means that the LSP is dominated by the Higgsino component to more that $`90\%`$. Finally Fig. 1 shows the contours of neutralino relic density $`\mathrm{\Omega }_\chi h^2`$ which was computed using the micrOMEGAs of Ref. . It is seen to generally lie below the lower limit of the cosmologically desirable range of eq.(2) by a factor of 2 to 4. This is due to the rapid pair annihilation processes $`\stackrel{~}{H}\stackrel{~}{H}W^+W^{},ZZ`$ via s-channel Z-boson and t-channel Higgsino exchanges, as mentioned earlier. Besides, the near degeneracy of the lighter chargino and neutralino masses (eq.1) leads to large coannihilations. In view of the mass degeneracy it is important to include radiative corrections to the $`\stackrel{~}{\chi }_{1,2}^0`$ and $`\stackrel{~}{\chi }_1^\pm `$ masses in the Higgsino LSP scenario . We have included this using the code of Manuel Drees. But, it does not enhance the neutralino relic density significantly. It should be noted here that the dominant gaugino component of the LSP (eq.13) comes from $`\stackrel{~}{W}`$ instead of $`\stackrel{~}{B}`$, in view of the inverted mass hierarchy $`\mu <M_2<M_1`$ in this model. Since, the winos have very similar annihilation mechanisms like the Higgsinos there is no increase in the relic density in the mixed Higgsino-gaugino region ($`Z_g>10\%`$).

## 4 Indirect Detection Rates

Since the Z-boson couples only to the Higgsino component of neutralino, the $`Z\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_1^0`$ coupling is proportional to $`N_{13,14}^2`$ . Moreover, the spin dependent force from Z-exchange is known to dominate the $`\stackrel{~}{\chi }_1^0`$ interaction rate with the solar matter, which is predominantly Hydrogen . Hence the solar capture rate of the Higgsino DM is predicted to be enormously larger than the bino DM of the minimal SUGRA model. This implies in turn an enormously higher rate of pair annihilation, $`\stackrel{~}{H}\stackrel{~}{H}W^+W^{}(ZZ)`$, since the capture and annihilation rates balance one another at equilibrium. The high energy neutrinos coming from W(Z) decay are expected to be detected at the large area neutrino telescopes like the IceCubes  and the ANTARES  via their charged current interaction ($`\nu \mu `$). The resulting muons constitute the so called indirect DM detection signal. Fig. 2 shows the indirect signal rate contours over the full parameter space which was computed by using DARKSUSY of Ref. . The signal contours are shown both with and without rescaling by a factor $`\xi =\mathrm{\Omega }_\chi h^2/0.05`$ . The denominator corresponds to $`\mathrm{\Omega }_m=0.1`$ which is the galactic DM relic density assumed in this computation. Even with rescaling one expects muon flux $`\varphi _\mu `$ of 5-100 $`\mathrm{events}/\mathrm{km}^2/\mathrm{year}`$ over practically the full parameter space of the model. In contrast the minimal SUGRA model predicts a $`\varphi _\mu <1\mathrm{event}/\mathrm{km}^2/\mathrm{year}`$ over the entire parameter space, except for a very narrow strip at the $`|\mu |=100\mathrm{GeV}`$ boundary corresponding to a larger Higgsino content in the LSP . The proposed large area neutrino telescopes like IceCube and ANTARES are expected to cover a detection area of 1$`\mathrm{km}^2`$. The irreducible background for these experiments, coming from the high energy neutrinos produced by the cosmic ray interaction with the solar corona is estimated to be $`\varphi _\mu 5\mathrm{events}/\mathrm{km}^2/\mathrm{year}`$ . Therefore these experiments can probe a signal of $`\varphi _\mu \begin{array}{c}>\hfill \\ \hfill \end{array}5\mathrm{events}/\mathrm{km}^2/\mathrm{year}`$, as expected over practically the full parameter space of this nonuniversal SUGRA model. It may be mentioned here that in the presence of some enhancement mechanism for the Higgsino DM relic density , it can become the dominant component of the cold dark matter. This will enhance the indirect detection rate further, as indicated by the contours without rescaling.

## 5 Direct Detection Rates

For the sake of completeness we have computed the $`\stackrel{~}{\chi }_1^0p`$ elastic scattering cross-sections in this model, which determine the signal rate in direct detection experiments. Both the spin-dependent and the spin-independent cross-sections have been computed using the DARKSUSY code .

The spin-dependent cross section is known to be dominated by Z-exchange. Therefore the spin-dependent cross-sections are much larger here compared to the minimal SUGRA model. Fig. 3 gives scatter plots of spin-dependent cross section against the LSP mass, both with and without rescaling for $`\mathrm{tan}\beta =`$10 and 50. Even the rescaled cross-section is 1-2 orders of magnitude larger than the minimal SUGRA predictions . Unfortunately, the direct detection experiments are not sensitive to the spin-dependent cross-section as they are based on heavy nuclei. For example the UKDMC detector is only sensitive to a spin-dependent cross section $`\begin{array}{c}>\hfill \\ \hfill \end{array}0.5\mathrm{pb}`$ , which is much above any SUSY model prediction.

The spin -independent (scalar) cross-section is dominated by Higgs exchanges. Since the Higgs couplings to a $`\stackrel{~}{\chi }^0\stackrel{~}{\chi }^0`$ pair is proportional to the product of their Higgsino and gaugino components , they are suppressed for both Higgsino and gaugino dominated DM. Fig. 4 gives scatter plots of the scalar cross-section with and without rescaling for $`\mathrm{tan}\beta =`$10 and 50. The upper range of the scatter-plots correspond to the mixed Higgsino-gaugino region ($`Z_g>10\%`$) of Fig. 1, as expected. The cross sections without the rescaling factors are moderately larger than the minimal SUGRA prediction . But the rescaled cross-sections are similar in size to the latter. Fig. 4 also shows that a significant part of the unrescaled cross-section lies above the discovery limits of the future CDMS  and GENIUS  experiments; but the rescaled cross-sections generally lie below these limits except for the mixed Higgsino-gaugino region. The DAMA  and present CDMS limits are also shown. In other words the upcoming experiments can detect the Higgsino DM if it is the dominant component of the galactic DM, but not if it is only a subdominant component of the latter.

## 6 Summary

We have investigated the dark matter phenomenology of a SUGRA model with nonuniversal gaugino masses. Its gauge kinetic function is a function of a nonsinglet chiral superfield, belonging to the 200-plet representation of SU(5). It is as predictive as the minimal SUGRA model and has less fine-tuning problem than the latter. It predicts a dominantly Higgsino LSP over the practically entire parameter space. The resulting thermal relic density of the Higgsino dark matter lies moderately below the cosmologically favoured range of eq.(2). Thus the Higgsino can only be a subleading component of the cold dark matter in the standard cosmological scenario. On the other hand its unsuppressed coupling to the Z-boson implies an enhanced rate of capture by the Sun. Consequently the predicted rate of indirect detection via high energy neutrinos coming from its pair-annihilation inside the Sun is much larger than the minimal SUGRA even after rescaling by the low density factor. This signal can be detected by the proposed $`\mathrm{km}^2`$ neutrino telescopes like IceCube and ANTARES over practically the full parameter space of the model. For the direct detections the predicted rate after rescaling is rougly similar to the minimal SUGRA prediction.

We thank Manuel Drees for discussion and for the use of his radiative correction code for chargino and neutralinos.


text:
# Comment on “Relating atomic-scale electronic phenomena to wave-like quasiparticle states in superconducting Bi2Sr2CaCu2O8+δ”

\[

\]

In a recent article , McElroy $`etal.`$ argued that the peaks observed by scanning tunneling measurements in the reciprocal space of the spatial modulation of local density-of-states (LDOS) in superconducting Bi2Sr2CaCu2O8 can be understood as due to quasiparticle interference between regions of high DOS on the Fermi surface. They proposed an “octet” model and have managed to describe many essential features of the peaks. In their paper, they pointed out that the location of the peaks possess mirror symmetry about the lattice axes and diagonals. In this Comment, we point out that there should exist a further symmetry, the $`umklapp`$ symmetry in a periodic lattice. As a consequence, some extra peaks which were not mentioned by the authors, should exist within their model and may have escaped from being detected. Moreover, these peaks can also disrupt the distinguishability of individual peaks.

In Fig. 1, we have illustratively shown some direct and $`umklapp`$ scattering wavevectors in the octet model at a bias voltage. The location of the wavevectors has now additional reflection symmetries about $`q_x,q_y=\pm \pi ,\pm 3\pi ,\pm 5\pi ,etc`$. There are two kinds of nonequivalent $`umklapp`$ generated peaks. Focusing on the quadrant $`0<q_x<\pi `$ and $`0<q_y<\pi `$ (peaks in other quadrants are obtained through symmetry operation on this), an extra peak denoted by $`q_5^{}`$ is seen on the lattice axis when the bias voltage exceeds a threshold. Peak $`q_5`$ was reported in Ref. 1 but $`q_5^{}`$ was not. Since $`umklapp`$ symmetry is dictated by the lattice and estabished by say, angle-resolved photoemission spectroscopy, $`q_5^{}`$ should exist correspondingly. More disturbingly, there also exists an $`umklapp`$ partner of $`q_4`$ (denoted by $`q_4^{}`$) which may touch $`q_{2,6}`$, or vice versa speaking. We draw the attention of the readers to the fact that it was mentioned in Ref.  that $`q_4`$ had been hard to detect, and we propose that this may be due to the close neighboring of two peaks.

In conclusion, there should exist $`umklapp`$ peaks to complete the quasiparticle interference picture, and those peaks may already have disrupted current data analysis. The crowding in of the $`umklapp`$ peaks into the quadrant can make the peaks hard to discern above some bias voltage. Raising the temperature at a fixed bias voltage will also be detrimental to the peaks, since the superconducting gap is closed and as a result the umklapp peaks crowd in. This could be part of the reason of the observed degradation of those peaks at higher biases mentioned in Ref. .